| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // 23 //
>> 24 // $Id: G4ComptonScattering.cc,v 1.18 2004/03/10 16:48:45 vnivanch Exp $
>> 25 // GEANT4 tag $Name: geant4-06-01 $
>> 26 //
27 // 27 //
28 //------------ G4ComptonScattering physics pro 28 //------------ G4ComptonScattering physics process -----------------------------
29 // by Michel Maire, April 19 29 // by Michel Maire, April 1996
30 // 30 //
31 // Modified: Michel Maire and Vladimir Ivanche << 31 // 28-05-96, DoIt() small change in ElecDirection, by M.Maire
32 // << 32 // 10-06-96, simplification in ComputeMicroscopicCrossSection(), by M.Maire
>> 33 // 21-06-96, SetCuts implementation, M.Maire
>> 34 // 13-09-96, small changes in DoIt for better efficiency. Thanks to P.Urban
>> 35 // 06-01-97, crossection table + meanfreepath table, M.Maire
>> 36 // 05-03-97, new Physics scheme, M.Maire
>> 37 // 28-03-97, protection in BuildPhysicsTable, M.Maire
>> 38 // 07-04-98, remove 'tracking cut' of the scattered gamma, MMa
>> 39 // 04-06-98, in DoIt, secondary production condition:
>> 40 // range>std::min(threshold,safety)
>> 41 // 13-08-98, new methods SetBining() PrintInfo()
>> 42 // 15-12-98, cross section=0 below 10 keV
>> 43 // 28-05-01, V.Ivanchenko minor changes to provide ANSI -wall compilation
>> 44 // 13-07-01, DoIt: suppression of production cut for the electron (mma)
>> 45 // 03-08-01, new methods Store/Retrieve PhysicsTable (mma)
>> 46 // 06-08-01, BuildThePhysicsTable() called from constructor (mma)
>> 47 // 17-09-01, migration of Materials to pure STL (mma)
>> 48 // 20-09-01, DoIt: fminimalEnergy = 1*eV (mma)
>> 49 // 01-10-01, come back to BuildPhysicsTable(const G4ParticleDefinition&)
>> 50 // 17-04-02, LowestEnergyLimit = 1*keV
33 // ------------------------------------------- 51 // -----------------------------------------------------------------------------
34 52
35 #include "G4ComptonScattering.hh" 53 #include "G4ComptonScattering.hh"
36 #include "G4SystemOfUnits.hh" << 54 #include "G4UnitsTable.hh"
37 #include "G4KleinNishinaModel.hh" <<
38 #include "G4KleinNishinaCompton.hh" <<
39 #include "G4Electron.hh" <<
40 #include "G4EmParameters.hh" <<
41 55
42 //....oooOO0OOooo........oooOO0OOooo........oo 56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
43 57
44 using namespace std; << 58 // constructor
>> 59
>> 60 G4ComptonScattering::G4ComptonScattering(const G4String& processName,
>> 61 G4ProcessType type):G4VDiscreteProcess (processName, type),
>> 62 theCrossSectionTable(NULL),
>> 63 theMeanFreePathTable(NULL),
>> 64 LowestEnergyLimit ( 1*keV),
>> 65 HighestEnergyLimit(100*GeV),
>> 66 NumbBinTable(80),
>> 67 fminimalEnergy(1*eV)
>> 68 {}
45 69
46 G4ComptonScattering::G4ComptonScattering(const << 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47 G4ProcessType type):G4VEmProcess (processNam << 71
>> 72 // destructor
>> 73
>> 74 G4ComptonScattering::~G4ComptonScattering()
48 { 75 {
49 SetStartFromNullFlag(true); << 76 if (theCrossSectionTable) {
50 SetBuildTableFlag(true); << 77 theCrossSectionTable->clearAndDestroy();
51 SetSecondaryParticle(G4Electron::Electron()) << 78 delete theCrossSectionTable;
52 SetProcessSubType(fComptonScattering); << 79 }
53 SetMinKinEnergyPrim(1*CLHEP::MeV); << 80
54 SetSplineFlag(true); << 81 if (theMeanFreePathTable) {
>> 82 theMeanFreePathTable->clearAndDestroy();
>> 83 delete theMeanFreePathTable;
>> 84 }
55 } 85 }
56 86
57 //....oooOO0OOooo........oooOO0OOooo........oo 87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 88
>> 89 void G4ComptonScattering::SetPhysicsTableBining(
>> 90 G4double lowE, G4double highE, G4int nBins)
>> 91 {
>> 92 LowestEnergyLimit = lowE; HighestEnergyLimit = highE; NumbBinTable = nBins;
>> 93 }
>> 94
>> 95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58 96
59 G4ComptonScattering::~G4ComptonScattering() = << 97 void G4ComptonScattering::BuildPhysicsTable(const G4ParticleDefinition&)
>> 98 // Build cross section and mean free path tables
>> 99 {
>> 100 G4double LowEdgeEnergy, Value;
>> 101 G4PhysicsLogVector* ptrVector;
>> 102
>> 103 // Build cross section per atom tables for the Compton Scattering process
>> 104
>> 105 if (theCrossSectionTable) {
>> 106 theCrossSectionTable->clearAndDestroy(); delete theCrossSectionTable;}
60 107
61 //....oooOO0OOooo........oooOO0OOooo........oo << 108 theCrossSectionTable = new G4PhysicsTable(G4Element::GetNumberOfElements());
>> 109 const G4ElementTable* theElementTable = G4Element::GetElementTable();
>> 110 G4double AtomicNumber;
>> 111 size_t J;
>> 112
>> 113 for ( J=0 ; J < G4Element::GetNumberOfElements(); J++ )
>> 114 {
>> 115 //create physics vector then fill it ....
>> 116 ptrVector = new G4PhysicsLogVector(LowestEnergyLimit,HighestEnergyLimit,
>> 117 NumbBinTable );
>> 118 AtomicNumber = (*theElementTable)[J]->GetZ();
>> 119
>> 120 for ( G4int i = 0 ; i < NumbBinTable ; i++ )
>> 121 {
>> 122 LowEdgeEnergy = ptrVector->GetLowEdgeEnergy(i);
>> 123 Value = ComputeCrossSectionPerAtom(LowEdgeEnergy, AtomicNumber);
>> 124 ptrVector->PutValue(i,Value);
>> 125 }
>> 126
>> 127 theCrossSectionTable->insertAt( J , ptrVector ) ;
>> 128
>> 129 }
>> 130
>> 131 // Build mean free path table for the Compton Scattering process
>> 132
>> 133 if (theMeanFreePathTable) {
>> 134 theMeanFreePathTable->clearAndDestroy(); delete theMeanFreePathTable;}
>> 135
>> 136 theMeanFreePathTable= new G4PhysicsTable(G4Material::GetNumberOfMaterials());
>> 137 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
>> 138 G4Material* material;
>> 139
>> 140 for ( J=0 ; J < G4Material::GetNumberOfMaterials(); J++ )
>> 141 {
>> 142 //create physics vector then fill it ....
>> 143 ptrVector = new G4PhysicsLogVector(LowestEnergyLimit,HighestEnergyLimit,
>> 144 NumbBinTable ) ;
>> 145 material = (*theMaterialTable)[J];
>> 146
>> 147 for ( G4int i = 0 ; i < NumbBinTable ; i++ )
>> 148 {
>> 149 LowEdgeEnergy = ptrVector->GetLowEdgeEnergy( i ) ;
>> 150 Value = ComputeMeanFreePath( LowEdgeEnergy, material);
>> 151 ptrVector->PutValue( i , Value ) ;
>> 152 }
>> 153
>> 154 theMeanFreePathTable->insertAt( J , ptrVector ) ;
>> 155 }
>> 156
>> 157 PrintInfoDefinition();
>> 158
>> 159 }
>> 160
>> 161 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
62 162
63 G4bool G4ComptonScattering::IsApplicable(const << 163 G4double G4ComptonScattering::ComputeCrossSectionPerAtom
>> 164 (G4double GammaEnergy, G4double Z)
>> 165
>> 166 // Calculates the cross section per atom in GEANT4 internal units.
>> 167 // A parametrized formula from L. Urban is used to estimate
>> 168 // the total cross section.
>> 169 // It gives a good description of the data from 10 keV to 100/Z GeV.
>> 170
>> 171 {
>> 172 G4double CrossSection = 0.0 ;
>> 173 if ( Z < 1. ) return CrossSection;
>> 174 if ( GammaEnergy < 1.*keV ) return CrossSection;
>> 175 if ( GammaEnergy > (100.*GeV/Z) ) return CrossSection;
>> 176
>> 177 static const G4double a = 20.0 , b = 230.0 , c = 440.0;
>> 178
>> 179 static const G4double
>> 180 d1= 2.7965e-1*barn, d2=-1.8300e-1*barn, d3= 6.7527 *barn, d4=-1.9798e+1*barn,
>> 181 e1= 1.9756e-5*barn, e2=-1.0205e-2*barn, e3=-7.3913e-2*barn, e4= 2.7079e-2*barn,
>> 182 f1=-3.9178e-7*barn, f2= 6.8241e-5*barn, f3= 6.0480e-5*barn, f4= 3.0274e-4*barn;
>> 183
>> 184 G4double p1Z = Z*(d1 + e1*Z + f1*Z*Z), p2Z = Z*(d2 + e2*Z + f2*Z*Z),
>> 185 p3Z = Z*(d3 + e3*Z + f3*Z*Z), p4Z = Z*(d4 + e4*Z + f4*Z*Z);
>> 186
>> 187 G4double X = GammaEnergy / electron_mass_c2 ;
>> 188
>> 189 return CrossSection = p1Z*log(1.+2*X)/X
>> 190 + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X);
>> 191 }
>> 192
>> 193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 194
>> 195 G4VParticleChange* G4ComptonScattering::PostStepDoIt(const G4Track& aTrack,
>> 196 const G4Step& aStep)
>> 197 //
>> 198 // The scattered gamma energy is sampled according to Klein - Nishina formula.
>> 199 // The random number techniques of Butcher & Messel are used
>> 200 // (Nuc Phys 20(1960),15).
>> 201 // GEANT4 internal units
>> 202 //
>> 203 // Note : Effects due to binding of atomic electrons are negliged.
>> 204
64 { 205 {
65 return (&p == G4Gamma::Gamma()); << 206 aParticleChange.Initialize(aTrack);
>> 207
>> 208 const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
>> 209 G4double GammaEnergy0 = aDynamicGamma->GetKineticEnergy();
>> 210 G4double E0_m = GammaEnergy0 / electron_mass_c2 ;
>> 211
>> 212 G4ParticleMomentum GammaDirection0 = aDynamicGamma->GetMomentumDirection();
>> 213
>> 214 //
>> 215 // sample the energy rate of the scattered gamma
>> 216 //
>> 217
>> 218 G4double epsilon, epsilonsq, onecost, sint2, greject ;
>> 219
>> 220 G4double epsilon0 = 1./(1. + 2*E0_m) , epsilon0sq = epsilon0*epsilon0;
>> 221 G4double alpha1 = - log(epsilon0) , alpha2 = 0.5*(1.- epsilon0sq);
>> 222
>> 223 do {
>> 224 if ( alpha1/(alpha1+alpha2) > G4UniformRand() )
>> 225 { epsilon = exp(-alpha1*G4UniformRand()); // epsilon0**r
>> 226 epsilonsq = epsilon*epsilon; }
>> 227 else {
>> 228 epsilonsq = epsilon0sq + (1.- epsilon0sq)*G4UniformRand();
>> 229 epsilon = sqrt(epsilonsq);
>> 230 };
>> 231 onecost = (1.- epsilon)/(epsilon*E0_m);
>> 232 sint2 = onecost*(2.-onecost);
>> 233 greject = 1. - epsilon*sint2/(1.+ epsilonsq);
>> 234 } while (greject < G4UniformRand());
>> 235
>> 236 //
>> 237 // scattered gamma angles. ( Z - axis along the parent gamma)
>> 238 //
>> 239
>> 240 G4double cosTeta = 1. - onecost , sinTeta = sqrt (sint2);
>> 241 G4double Phi = twopi * G4UniformRand();
>> 242 G4double dirx = sinTeta*cos(Phi), diry = sinTeta*sin(Phi), dirz = cosTeta;
>> 243
>> 244 //
>> 245 // update G4VParticleChange for the scattered gamma
>> 246 //
>> 247
>> 248 G4ThreeVector GammaDirection1 ( dirx,diry,dirz );
>> 249 GammaDirection1.rotateUz(GammaDirection0);
>> 250 aParticleChange.SetMomentumChange( GammaDirection1 );
>> 251 G4double GammaEnergy1 = epsilon*GammaEnergy0;
>> 252 G4double localEnergyDeposit = 0.;
>> 253
>> 254 if (GammaEnergy1 > fminimalEnergy)
>> 255 {
>> 256 aParticleChange.SetEnergyChange( GammaEnergy1 );
>> 257 }
>> 258 else
>> 259 {
>> 260 localEnergyDeposit += GammaEnergy1;
>> 261 aParticleChange.SetEnergyChange(0.) ;
>> 262 aParticleChange.SetStatusChange(fStopAndKill);
>> 263 }
>> 264
>> 265 //
>> 266 // kinematic of the scattered electron
>> 267 //
>> 268
>> 269 G4double ElecKineEnergy = GammaEnergy0 - GammaEnergy1;
>> 270
>> 271 if (ElecKineEnergy > fminimalEnergy)
>> 272 {
>> 273 G4double ElecMomentum = sqrt(ElecKineEnergy*
>> 274 (ElecKineEnergy+2.*electron_mass_c2));
>> 275 G4ThreeVector ElecDirection (
>> 276 (GammaEnergy0*GammaDirection0 - GammaEnergy1*GammaDirection1)
>> 277 *(1./ElecMomentum) );
>> 278
>> 279 // create G4DynamicParticle object for the electron.
>> 280 G4DynamicParticle* aElectron= new G4DynamicParticle(
>> 281 G4Electron::Electron(),ElecDirection,ElecKineEnergy);
>> 282
>> 283 aParticleChange.SetNumberOfSecondaries(1);
>> 284 aParticleChange.AddSecondary( aElectron );
>> 285 }
>> 286 else
>> 287 {
>> 288 aParticleChange.SetNumberOfSecondaries(0);
>> 289 localEnergyDeposit += ElecKineEnergy;
>> 290 }
>> 291 aParticleChange.SetLocalEnergyDeposit (localEnergyDeposit);
>> 292
>> 293 // Reset NbOfInteractionLengthLeft and return aParticleChange
>> 294 return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep);
66 } 295 }
67 296
68 //....oooOO0OOooo........oooOO0OOooo........oo << 297 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69 298
70 void G4ComptonScattering::InitialiseProcess(co << 299 G4bool G4ComptonScattering::StorePhysicsTable(G4ParticleDefinition* particle,
>> 300 const G4String& directory,
>> 301 G4bool ascii)
71 { 302 {
72 if(!isInitialised) { << 303 G4String filename;
73 isInitialised = true; << 304
74 if(nullptr == EmModel(0)) { SetEmModel(new << 305 // store cross section table
75 G4EmParameters* param = G4EmParameters::In << 306 filename = GetPhysicsTableFileName(particle,directory,"CrossSection",ascii);
76 EmModel(0)->SetLowEnergyLimit(param->MinKi << 307 if ( !theCrossSectionTable->StorePhysicsTable(filename, ascii) ){
77 EmModel(0)->SetHighEnergyLimit(param->MaxK << 308 G4cout << " FAIL theCrossSectionTable->StorePhysicsTable in " << filename
78 AddEmModel(1, EmModel(0)); << 309 << G4endl;
79 } << 310 return false;
>> 311 }
>> 312
>> 313 // store mean free path table
>> 314 filename = GetPhysicsTableFileName(particle,directory,"MeanFreePath",ascii);
>> 315 if ( !theMeanFreePathTable->StorePhysicsTable(filename, ascii) ){
>> 316 G4cout << " FAIL theMeanFreePathTable->StorePhysicsTable in " << filename
>> 317 << G4endl;
>> 318 return false;
>> 319 }
>> 320
>> 321 G4cout << GetProcessName() << " for " << particle->GetParticleName()
>> 322 << ": Success to store the PhysicsTables in "
>> 323 << directory << G4endl;
>> 324 return true;
80 } 325 }
81 326
82 //....oooOO0OOooo........oooOO0OOooo........oo 327 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
83 328
84 void G4ComptonScattering::ProcessDescription(s << 329 G4bool G4ComptonScattering::RetrievePhysicsTable(G4ParticleDefinition* particle,
>> 330 const G4String& directory,
>> 331 G4bool ascii)
85 { 332 {
86 out << " Compton scattering"; << 333 // delete theCrossSectionTable and theMeanFreePathTable
87 G4VEmProcess::ProcessDescription(out); << 334 if (theCrossSectionTable != 0) {
>> 335 theCrossSectionTable->clearAndDestroy();
>> 336 delete theCrossSectionTable;
>> 337 }
>> 338 if (theMeanFreePathTable != 0) {
>> 339 theMeanFreePathTable->clearAndDestroy();
>> 340 delete theMeanFreePathTable;
>> 341 }
>> 342
>> 343 G4String filename;
>> 344
>> 345 // retreive cross section table
>> 346 filename = GetPhysicsTableFileName(particle,directory,"CrossSection",ascii);
>> 347 theCrossSectionTable = new G4PhysicsTable(G4Element::GetNumberOfElements());
>> 348 if ( !theCrossSectionTable->RetrievePhysicsTable(filename, ascii) ){
>> 349 G4cout << " FAIL theCrossSectionTable->RetrievePhysicsTable in " << filename
>> 350 << G4endl;
>> 351 return false;
>> 352 }
>> 353
>> 354 // retreive mean free path table
>> 355 filename = GetPhysicsTableFileName(particle,directory,"MeanFreePath",ascii);
>> 356 theMeanFreePathTable = new G4PhysicsTable(G4Material::GetNumberOfMaterials());
>> 357 if ( !theMeanFreePathTable->RetrievePhysicsTable(filename, ascii) ){
>> 358 G4cout << " FAIL theMeanFreePathTable->RetrievePhysicsTable in " << filename
>> 359 << G4endl;
>> 360 return false;
>> 361 }
>> 362
>> 363 G4cout << GetProcessName() << " for " << particle->GetParticleName()
>> 364 << ": Success to retrieve the PhysicsTables from "
>> 365 << directory << G4endl;
>> 366 return true;
88 } 367 }
>> 368
>> 369 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
89 370
90 //....oooOO0OOooo........oooOO0OOooo........oo << 371 void G4ComptonScattering::PrintInfoDefinition()
>> 372 {
>> 373 G4String comments = "Total cross sections from a parametrisation. ";
>> 374 comments += "Good description from 10 KeV to (100/Z) GeV. \n";
>> 375 comments += " Scattered gamma energy according Klein-Nishina.";
>> 376
>> 377 G4cout << G4endl << GetProcessName() << ": " << comments
>> 378 << "\n PhysicsTables from "
>> 379 << G4BestUnit(LowestEnergyLimit,"Energy")
>> 380 << " to " << G4BestUnit(HighestEnergyLimit,"Energy")
>> 381 << " in " << NumbBinTable << " bins. \n";
>> 382 }
>> 383
>> 384 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
91 385