Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4KleinNishinaModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/standard/src/G4KleinNishinaModel.cc (Version 11.3.0) and /processes/electromagnetic/standard/src/G4KleinNishinaModel.cc (Version 7.0.p1)


  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 24 // *******************************************    
 25 //                                                
 26 //                                                
 27 // -------------------------------------------    
 28 //                                                
 29 // GEANT4 Class file                              
 30 //                                                
 31 //                                                
 32 // File name:     G4KleinNishinaModel             
 33 //                                                
 34 // Author:        Vladimir Ivanchenko on base     
 35 //                                                
 36 // Creation date: 13.06.2010                      
 37 //                                                
 38 // Modifications:                                 
 39 //                                                
 40 // Class Description:                             
 41 //                                                
 42 // -------------------------------------------    
 43 //                                                
 44 //....oooOO0OOooo........oooOO0OOooo........oo    
 45 //....oooOO0OOooo........oooOO0OOooo........oo    
 46                                                   
 47 #include "G4KleinNishinaModel.hh"                 
 48 #include "G4PhysicalConstants.hh"                 
 49 #include "G4SystemOfUnits.hh"                     
 50 #include "G4Electron.hh"                          
 51 #include "G4Gamma.hh"                             
 52 #include "Randomize.hh"                           
 53 #include "G4RandomDirection.hh"                   
 54 #include "G4DataVector.hh"                        
 55 #include "G4ParticleChangeForGamma.hh"            
 56 #include "G4VAtomDeexcitation.hh"                 
 57 #include "G4AtomicShells.hh"                      
 58 #include "G4LossTableManager.hh"                  
 59 #include "G4Log.hh"                               
 60 #include "G4Exp.hh"                               
 61                                                   
 62 //....oooOO0OOooo........oooOO0OOooo........oo    
 63                                                   
 64 using namespace std;                              
 65                                                   
 66 G4KleinNishinaModel::G4KleinNishinaModel(const    
 67   : G4VEmModel(nam),                              
 68     lv1(0.,0.,0.,0.),                             
 69     lv2(0.,0.,0.,0.),                             
 70     bst(0.,0.,0.)                                 
 71 {                                                 
 72   theGamma = G4Gamma::Gamma();                    
 73   theElectron = G4Electron::Electron();           
 74   lowestSecondaryEnergy = 10*eV;                  
 75   limitFactor       = 4;                          
 76   fProbabilities.resize(9,0.0);                   
 77   SetDeexcitationFlag(true);                      
 78   fParticleChange = nullptr;                      
 79   fAtomDeexcitation = nullptr;                    
 80 }                                                 
 81                                                   
 82 //....oooOO0OOooo........oooOO0OOooo........oo    
 83                                                   
 84 G4KleinNishinaModel::~G4KleinNishinaModel() =     
 85                                                   
 86 //....oooOO0OOooo........oooOO0OOooo........oo    
 87                                                   
 88 void G4KleinNishinaModel::Initialise(const G4P    
 89                                      const G4D    
 90 {                                                 
 91   fAtomDeexcitation = G4LossTableManager::Inst    
 92   if(IsMaster()) { InitialiseElementSelectors(    
 93   if(nullptr == fParticleChange) {                
 94     fParticleChange = GetParticleChangeForGamm    
 95   }                                               
 96 }                                                 
 97                                                   
 98 //....oooOO0OOooo........oooOO0OOooo........oo    
 99                                                   
100 void G4KleinNishinaModel::InitialiseLocal(cons    
101                                           G4VE    
102 {                                                 
103   SetElementSelectors(masterModel->GetElementS    
104 }                                                 
105                                                   
106 //....oooOO0OOooo........oooOO0OOooo........oo    
107                                                   
108 G4double                                          
109 G4KleinNishinaModel::ComputeCrossSectionPerAto    
110                                                   
111                                                   
112                                                   
113 {                                                 
114   G4double xSection = 0.0 ;                       
115   if (gammaEnergy <= LowEnergyLimit()) { retur    
116                                                   
117   static const G4double a = 20.0 , b = 230.0 ,    
118                                                   
119 static const G4double                             
120   d1= 2.7965e-1*CLHEP::barn, d2=-1.8300e-1*CLH    
121   d3= 6.7527   *CLHEP::barn, d4=-1.9798e+1*CLH    
122   e1= 1.9756e-5*CLHEP::barn, e2=-1.0205e-2*CLH    
123   e3=-7.3913e-2*CLHEP::barn, e4= 2.7079e-2*CLH    
124   f1=-3.9178e-7*CLHEP::barn, f2= 6.8241e-5*CLH    
125   f3= 6.0480e-5*CLHEP::barn, f4= 3.0274e-4*CLH    
126                                                   
127   G4double p1Z = Z*(d1 + e1*Z + f1*Z*Z), p2Z =    
128            p3Z = Z*(d3 + e3*Z + f3*Z*Z), p4Z =    
129                                                   
130   G4double T0  = 15.0*keV;                        
131   if (Z < 1.5) { T0 = 40.0*keV; }                 
132                                                   
133   G4double X   = max(gammaEnergy, T0) / electr    
134   xSection = p1Z*G4Log(1.+2.*X)/X                 
135                + (p2Z + p3Z*X + p4Z*X*X)/(1. +    
136                                                   
137   //  modification for low energy. (special ca    
138   static const G4double dT0 = keV;                
139   if (gammaEnergy < T0) {                         
140     X = (T0+dT0) / electron_mass_c2 ;             
141     G4double sigma = p1Z*G4Log(1.+2*X)/X          
142                     + (p2Z + p3Z*X + p4Z*X*X)/    
143     G4double   c1 = -T0*(sigma-xSection)/(xSec    
144     G4double   c2 = 0.150;                        
145     if (Z > 1.5) { c2 = 0.375-0.0556*G4Log(Z);    
146     G4double    y = G4Log(gammaEnergy/T0);        
147     xSection *= G4Exp(-y*(c1+c2*y));              
148   }                                               
149                                                   
150   if(xSection < 0.0) { xSection = 0.0; }          
151   //  G4cout << "e= " << GammaEnergy << " Z= "    
152   //  << " cross= " << xSection << G4endl;        
153   return xSection;                                
154 }                                                 
155                                                   
156 //....oooOO0OOooo........oooOO0OOooo........oo    
157                                                   
158 void G4KleinNishinaModel::SampleSecondaries(      
159                              std::vector<G4Dyn    
160                              const G4MaterialC    
161                              const G4DynamicPa    
162                              G4double,            
163                              G4double)            
164 {                                                 
165   // primary gamma                                
166   G4double energy = aDynamicGamma->GetKineticE    
167                                                   
168   // do nothing below the threshold               
169   if(energy <= LowEnergyLimit()) { return; }      
170                                                   
171   G4ThreeVector direction = aDynamicGamma->Get    
172                                                   
173   // select atom                                  
174   const G4Element* elm = SelectRandomAtom(coup    
175                                                   
176   // select shell first                           
177   G4int nShells = elm->GetNbOfAtomicShells();     
178   if(nShells > (G4int)fProbabilities.size()) {    
179   G4double totprob = 0.0;                         
180   G4int i;                                        
181   for(i=0; i<nShells; ++i) {                      
182     //G4double bindingEnergy = elm->GetAtomicS    
183     totprob += elm->GetNbOfShellElectrons(i);     
184     //totprob += elm->GetNbOfShellElectrons(i)    
185     fProbabilities[i] = totprob;                  
186   }                                               
187                                                   
188   // Loop on sampling                             
189   static const G4int nlooplim = 1000;             
190   G4int nloop = 0;                                
191                                                   
192   G4double bindingEnergy, ePotEnergy, eKinEner    
193   G4double gamEnergy0, gamEnergy1;                
194                                                   
195   CLHEP::HepRandomEngine* rndmEngineMod = G4Ra    
196   G4double rndm[4];                               
197                                                   
198   do {                                            
199     ++nloop;                                      
200                                                   
201     // 4 random numbers to select e-              
202     rndmEngineMod->flatArray(4, rndm);            
203     G4double xprob = totprob*rndm[0];             
204                                                   
205     // select shell                               
206     for(i=0; i<nShells; ++i) { if(xprob <= fPr    
207                                                   
208     bindingEnergy = elm->GetAtomicShell(i);       
209     lv1.set(0.0,0.0,energy,energy);               
210     /*                                            
211     G4cout << "nShells= " << nShells << " i= "    
212        << " Egamma= " << energy << " Ebind= "     
213        << G4endl;                                 
214     */                                            
215     // for rest frame of the electron             
216     G4double x = -G4Log(rndm[1]);                 
217     eKinEnergy = bindingEnergy*x;                 
218     ePotEnergy = bindingEnergy*(1.0 + x);         
219                                                   
220     // for rest frame of the electron             
221     G4double eTotMomentum = sqrt(eKinEnergy*(e    
222     G4double phi = rndm[2]*twopi;                 
223     G4double costet = 2*rndm[3] - 1;              
224     G4double sintet = sqrt((1 - costet)*(1 + c    
225     lv2.set(eTotMomentum*sintet*cos(phi),eTotM    
226             eTotMomentum*costet,eKinEnergy + e    
227     bst = lv2.boostVector();                      
228     lv1.boost(-bst);                              
229                                                   
230     gamEnergy0 = lv1.e();                         
231                                                   
232     // In the rest frame of the electron          
233     // The scattered gamma energy is sampled a    
234     // The random number techniques of Butcher    
235     // (Nuc Phys 20(1960),15).                    
236     G4double E0_m = gamEnergy0/electron_mass_c    
237                                                   
238     //G4cout << "Nloop= "<< nloop << " Ecm(keV    
239     //                                            
240     // sample the energy rate of the scattered    
241     //                                            
242                                                   
243     G4double epsilon, epsilonsq, onecost, sint    
244                                                   
245     G4double eps0       = 1./(1 + 2*E0_m);        
246     G4double epsilon0sq = eps0*eps0;              
247     G4double alpha1     = - G4Log(eps0);          
248     G4double alpha2     = alpha1 + 0.5*(1 - ep    
249                                                   
250     do {                                          
251       ++nloop;                                    
252       // false interaction if too many iterati    
253       if(nloop > nlooplim) { return; }            
254                                                   
255       // 3 random numbers to sample scattering    
256       rndmEngineMod->flatArray(3, rndm);          
257                                                   
258       if ( alpha1 > alpha2*rndm[0] ) {            
259         epsilon   = G4Exp(-alpha1*rndm[1]);       
260         epsilonsq = epsilon*epsilon;              
261                                                   
262       } else {                                    
263         epsilonsq = epsilon0sq + (1.- epsilon0    
264         epsilon   = sqrt(epsilonsq);              
265       }                                           
266                                                   
267       onecost = (1.- epsilon)/(epsilon*E0_m);     
268       sint2   = onecost*(2.-onecost);             
269       greject = 1. - epsilon*sint2/(1.+ epsilo    
270                                                   
271       // Loop checking, 03-Aug-2015, Vladimir     
272     } while (greject < rndm[2]);                  
273     gamEnergy1 = epsilon*gamEnergy0;              
274                                                   
275     // before scattering total 4-momentum in e    
276     lv2.set(0.0,0.0,0.0,electron_mass_c2);        
277     lv2 += lv1;                                   
278                                                   
279     //                                            
280     // scattered gamma angles. ( Z - axis alon    
281     //                                            
282     if(sint2 < 0.0) { sint2 = 0.0; }              
283     costet = 1. - onecost;                        
284     sintet = sqrt(sint2);                         
285     phi  = twopi * rndmEngineMod->flat();         
286                                                   
287     // e- recoil                                  
288     //                                            
289     // in  rest frame of the electron             
290     G4ThreeVector gamDir = lv1.vect().unit();     
291     G4ThreeVector v = G4ThreeVector(sintet*cos    
292     v.rotateUz(gamDir);                           
293     lv1.set(gamEnergy1*v.x(),gamEnergy1*v.y(),    
294     lv2 -= lv1;                                   
295     //G4cout<<"Egam(keV)= " << lv1.e()/keV        
296     //          <<" Ee(keV)= " << (lv2.e()-ele    
297     lv2.boost(bst);                               
298     eKinEnergy = lv2.e() - electron_mass_c2 -     
299     //G4cout << "Nloop= " << nloop << " eKinEn    
300                                                   
301     // Loop checking, 03-Aug-2015, Vladimir Iv    
302   } while ( eKinEnergy < 0.0 );                   
303                                                   
304   //                                              
305   // update G4VParticleChange for the scattere    
306   //                                              
307                                                   
308   lv1.boost(bst);                                 
309   gamEnergy1 = lv1.e();                           
310   if(gamEnergy1 > lowestSecondaryEnergy) {        
311     G4ThreeVector gamDirection1 = lv1.vect().u    
312     gamDirection1.rotateUz(direction);            
313     fParticleChange->ProposeMomentumDirection(    
314   } else {                                        
315     fParticleChange->ProposeTrackStatus(fStopA    
316     gamEnergy1 = 0.0;                             
317   }                                               
318   fParticleChange->SetProposedKineticEnergy(ga    
319                                                   
320   //                                              
321   // kinematic of the scattered electron          
322   //                                              
323                                                   
324   if(eKinEnergy > lowestSecondaryEnergy) {        
325     G4ThreeVector eDirection = lv2.vect().unit    
326     eDirection.rotateUz(direction);               
327     auto dp = new G4DynamicParticle(theElectro    
328     fvect->push_back(dp);                         
329   } else { eKinEnergy = 0.0; }                    
330                                                   
331   G4double edep = energy - gamEnergy1 - eKinEn    
332   G4double esec = 0.0;                            
333                                                   
334   // sample deexcitation                          
335   //                                              
336   if(nullptr != fAtomDeexcitation) {              
337     G4int index = couple->GetIndex();             
338     if(fAtomDeexcitation->CheckDeexcitationAct    
339       G4int Z = elm->GetZasInt();                 
340       auto as = (G4AtomicShellEnumerator)(i);     
341       const G4AtomicShell* shell = fAtomDeexci    
342       G4int nbefore = (G4int)fvect->size();       
343       fAtomDeexcitation->GenerateParticles(fve    
344       G4int nafter = (G4int)fvect->size();        
345       //G4cout << "N1= " << nbefore << "  N2=     
346       for (G4int j=nbefore; j<nafter; ++j) {      
347         G4double e = ((*fvect)[j])->GetKinetic    
348         if(esec + e > edep) {                     
349           // correct energy in order to have e    
350           e = edep - esec;                        
351           ((*fvect)[j])->SetKineticEnergy(e);     
352           esec += e;                              
353           /*                                      
354             G4cout << "### G4KleinNishinaModel    
355                    << " Esec(eV)= " << esec/eV    
356                    << " E["<< j << "](eV)= " <    
357                    << " N= " << nafter            
358                    << " Z= " << Z << " shell=     
359                    << "  Ebind(keV)= " << bind    
360                    << "  Eshell(keV)= " << she    
361                    << G4endl;                     
362           */                                      
363           // delete the rest of secondaries (s    
364           for (G4int jj=nafter-1; jj>j; --jj)     
365             delete (*fvect)[jj];                  
366             fvect->pop_back();                    
367           }                                       
368           break;                                  
369         }                                         
370         esec += e;                                
371       }                                           
372       edep -= esec;                               
373     }                                             
374   }                                               
375   if(std::abs(energy - gamEnergy1 - eKinEnergy    
376     G4cout << "### G4KleinNishinaModel dE(eV)=    
377            << (energy - gamEnergy1 - eKinEnerg    
378            << " shell= " << i                     
379            << "  E(keV)= " << energy/keV          
380            << "  Ebind(keV)= " << bindingEnerg    
381            << "  Eg(keV)= " << gamEnergy1/keV     
382            << "  Ee(keV)= " << eKinEnergy/keV     
383            << "  Esec(keV)= " << esec/keV         
384            << "  Edep(keV)= " << edep/keV         
385            << G4endl;                             
386   }                                               
387   // energy balance                               
388   if(edep > 0.0) {                                
389     fParticleChange->ProposeLocalEnergyDeposit    
390   }                                               
391 }                                                 
392                                                   
393 //....oooOO0OOooo........oooOO0OOooo........oo    
394                                                   
395