Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4LindhardSorensenIonModel.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/G4LindhardSorensenIonModel.cc (Version 11.3.0) and /processes/electromagnetic/standard/src/G4LindhardSorensenIonModel.cc (Version 9.4.p2)


  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 // GEANT4 Class header file                       
 29 //                                                
 30 //                                                
 31 // File name:     G4LindhardSorensenIonModel      
 32 //                                                
 33 // Author:        Alexander Bagulya & Vladimir    
 34 //                                                
 35 // Creation date: 16.04.2018                      
 36 //                                                
 37 //                                                
 38 // -------------------------------------------    
 39 //                                                
 40                                                   
 41 //....oooOO0OOooo........oooOO0OOooo........oo    
 42 //....oooOO0OOooo........oooOO0OOooo........oo    
 43                                                   
 44 #include "G4LindhardSorensenIonModel.hh"          
 45 #include "Randomize.hh"                           
 46 #include "G4PhysicalConstants.hh"                 
 47 #include "G4SystemOfUnits.hh"                     
 48 #include "G4Electron.hh"                          
 49 #include "G4LossTableManager.hh"                  
 50 #include "G4EmCorrections.hh"                     
 51 #include "G4ParticleChangeForLoss.hh"             
 52 #include "G4Log.hh"                               
 53 #include "G4DeltaAngle.hh"                        
 54 #include "G4LindhardSorensenData.hh"              
 55 #include "G4BraggModel.hh"                        
 56 #include "G4BetheBlochModel.hh"                   
 57 #include "G4IonICRU73Data.hh"                     
 58 #include "G4AutoLock.hh"                          
 59                                                   
 60 //....oooOO0OOooo........oooOO0OOooo........oo    
 61                                                   
 62 using namespace std;                              
 63                                                   
 64 G4LindhardSorensenData* G4LindhardSorensenIonM    
 65 G4IonICRU73Data* G4LindhardSorensenIonModel::f    
 66                                                   
 67 namespace                                         
 68 {                                                 
 69   G4Mutex ionXSMutex = G4MUTEX_INITIALIZER;       
 70 }                                                 
 71                                                   
 72 G4LindhardSorensenIonModel::G4LindhardSorensen    
 73                                                   
 74   : G4VEmModel(nam),                              
 75     twoln10(2.0*G4Log(10.0))                      
 76 {                                                 
 77   theElectron = G4Electron::Electron();           
 78   corr = G4LossTableManager::Instance()->EmCor    
 79   nist = G4NistManager::Instance();               
 80   fBraggModel = new G4BraggModel();               
 81   fBBModel = new G4BetheBlochModel();             
 82   fElimit = 2.0*CLHEP::MeV;                       
 83 }                                                 
 84                                                   
 85 //....oooOO0OOooo........oooOO0OOooo........oo    
 86                                                   
 87 G4LindhardSorensenIonModel::~G4LindhardSorense    
 88   if(isFirst) {                                   
 89     delete lsdata;                                
 90     delete fIonData;                              
 91     lsdata = nullptr;                             
 92     fIonData = nullptr;                           
 93   }                                               
 94 }                                                 
 95                                                   
 96 //....oooOO0OOooo........oooOO0OOooo........oo    
 97                                                   
 98 void G4LindhardSorensenIonModel::Initialise(co    
 99                                             co    
100 {                                                 
101   fBraggModel->Initialise(p, ptr);                
102   fBBModel->Initialise(p, ptr);                   
103   SetParticle(p);                                 
104                                                   
105   // always false before the run                  
106   SetDeexcitationFlag(false);                     
107                                                   
108   if(nullptr == fParticleChange) {                
109     fParticleChange = GetParticleChangeForLoss    
110     if(UseAngularGeneratorFlag() && nullptr ==    
111       SetAngularDistribution(new G4DeltaAngle(    
112     }                                             
113   }                                               
114   if(nullptr == lsdata) {                         
115     G4AutoLock l(&ionXSMutex);                    
116     if(nullptr == lsdata) {                       
117       isFirst = true;                             
118       lsdata = new G4LindhardSorensenData();      
119       fIonData = new G4IonICRU73Data();           
120     }                                             
121     l.unlock();                                   
122   }                                               
123   if(isFirst) {                                   
124     fIonData->Initialise();                       
125   }                                               
126 }                                                 
127                                                   
128 //....oooOO0OOooo........oooOO0OOooo........oo    
129                                                   
130 G4double                                          
131 G4LindhardSorensenIonModel::GetChargeSquareRat    
132                                                   
133                                                   
134 {                                                 
135   chargeSquare = corr->EffectiveChargeSquareRa    
136   return chargeSquare;                            
137 }                                                 
138                                                   
139 //....oooOO0OOooo........oooOO0OOooo........oo    
140                                                   
141 G4double                                          
142 G4LindhardSorensenIonModel::GetParticleCharge(    
143                                                   
144                                                   
145 {                                                 
146   return corr->GetParticleCharge(p,mat,kinEner    
147 }                                                 
148                                                   
149 //....oooOO0OOooo........oooOO0OOooo........oo    
150                                                   
151 void G4LindhardSorensenIonModel::SetupParamete    
152 {                                                 
153   mass = particle->GetPDGMass();                  
154   spin = particle->GetPDGSpin();                  
155   charge = particle->GetPDGCharge()*inveplus;     
156   Zin = G4lrint(std::abs(charge));                
157   chargeSquare = charge*charge;                   
158   eRatio = CLHEP::electron_mass_c2/mass;          
159   pRatio = CLHEP::proton_mass_c2/mass;            
160   const G4double aMag =                           
161     1./(0.5*CLHEP::eplus*CLHEP::hbar_Planck*CL    
162   G4double magmom = particle->GetPDGMagneticMo    
163   magMoment2 = magmom*magmom - 1.0;               
164   G4double x = 0.8426*CLHEP::GeV;                 
165   if(spin == 0.0 && mass < GeV) { x = 0.736*CL    
166   else if (Zin > 1) { x /= nist->GetA27(Zin);     
167                                                   
168   formfact = 2.0*CLHEP::electron_mass_c2/(x*x)    
169   tlimit = 2.0/formfact;                          
170 }                                                 
171                                                   
172 //....oooOO0OOooo........oooOO0OOooo........oo    
173                                                   
174 G4double G4LindhardSorensenIonModel::MinEnergy    
175                                          const    
176 {                                                 
177   return couple->GetMaterial()->GetIonisation(    
178 }                                                 
179                                                   
180 //....oooOO0OOooo........oooOO0OOooo........oo    
181                                                   
182 G4double                                          
183 G4LindhardSorensenIonModel::ComputeCrossSectio    
184                                                   
185                                                   
186                                                   
187                                                   
188 {                                                 
189   // take into account formfactor                 
190   G4double tmax = MaxSecondaryEnergy(p, kinEne    
191   G4double emax = std::min(tmax, maxKinEnergy)    
192   G4double escaled = kinEnergy*pRatio;            
193   G4double cross = (escaled <= fElimit)           
194     ? fBraggModel->ComputeCrossSectionPerElect    
195     : fBBModel->ComputeCrossSectionPerElectron    
196   // G4cout << "LS: e= " << kinEnergy << " tmi    
197   //        << " tmax= " << maxEnergy << " cro    
198   return cross;                                   
199 }                                                 
200                                                   
201 //....oooOO0OOooo........oooOO0OOooo........oo    
202                                                   
203 G4double G4LindhardSorensenIonModel::ComputeCr    
204                                            con    
205                                                   
206                                                   
207                                                   
208                                                   
209 {                                                 
210   return Z*ComputeCrossSectionPerElectron(p,ki    
211 }                                                 
212                                                   
213 //....oooOO0OOooo........oooOO0OOooo........oo    
214                                                   
215 G4double G4LindhardSorensenIonModel::CrossSect    
216                                            con    
217                                            con    
218                                                   
219                                                   
220                                                   
221 {                                                 
222   return material->GetElectronDensity()           
223     *ComputeCrossSectionPerElectron(p,kineticE    
224 }                                                 
225                                                   
226 //....oooOO0OOooo........oooOO0OOooo........oo    
227                                                   
228 G4double                                          
229 G4LindhardSorensenIonModel::ComputeDEDXPerVolu    
230                                                   
231                                                   
232                                                   
233 {                                                 
234   // formfactor is taken into account in Corre    
235   G4double tmax      = MaxSecondaryEnergy(p, k    
236   G4double cutEnergy = std::min(std::min(cut,t    
237                                                   
238   G4double escaled = kinEnergy*pRatio;            
239   G4double dedx = (escaled <= fElimit)            
240     ? fBraggModel->ComputeDEDXPerVolume(mat, p    
241     : fBBModel->ComputeDEDXPerVolume(mat, p, k    
242                                                   
243   //G4cout << "E(MeV)=" << kinEnergy/MeV << "     
244   //  << "  " << mat->GetName() << " Ecut(MeV)    
245   return dedx;                                    
246 }                                                 
247                                                   
248 //....oooOO0OOooo........oooOO0OOooo........oo    
249                                                   
250 void G4LindhardSorensenIonModel::CorrectionsAl    
251                                  const G4Mater    
252                                  const G4Dynam    
253                                  const G4doubl    
254                                        G4doubl    
255 {                                                 
256   // no correction at the last step               
257   const G4double preKinEnergy = dp->GetKinetic    
258   if(eloss >= preKinEnergy) { return; }           
259                                                   
260   const G4ParticleDefinition* p = dp->GetDefin    
261   SetParticle(p);                                 
262   const G4Material* mat = couple->GetMaterial(    
263   const G4double eDensity = mat->GetElectronDe    
264   const G4double e = std::max(preKinEnergy - e    
265   const G4double tmax = MaxSecondaryEnergy(p,     
266   const G4double escaled = e*pRatio;              
267   const G4double tau = e/mass;                    
268   const G4double q2 = corr->EffectiveChargeSqu    
269   const G4int Z = p->GetAtomicNumber();           
270                                                   
271   G4double res = 0.0;                             
272   if(escaled <= fElimit) {                        
273     // data from ICRU73 or ICRU90                 
274     if(Z > 2 && Z <= 80) {                        
275       res = fIonData->GetDEDX(mat, Z, escaled,    
276       /*                                          
277   G4cout << "  GetDEDX for Z=" << Z << " in "     
278   << " Escaled=" << escaled << " E="              
279   << e << " dEdx=" << res << G4endl;              
280       */                                          
281     }                                             
282     if(res > 0.0) {                               
283       auto pcuts = couple->GetProductionCuts()    
284       G4double cut = (nullptr == pcuts) ? tmax    
285       if(cut < tmax) {                            
286   const G4double x = cut/tmax;                    
287   res += (G4Log(x)*(tau + 1.)*(tau + 1.)/(tau     
288     *q2*CLHEP::twopi_mc2_rcl2*eDensity;           
289       }                                           
290       res *= length;                              
291     } else {                                      
292       // simplified correction                    
293       res = eloss*q2/chargeSquare;                
294     }                                             
295   } else {                                        
296     // Lindhard-Sorensen model                    
297     const G4double gam = tau + 1.0;               
298     const G4double beta2 = tau * (tau+2.0)/(ga    
299     G4double deltaL0 = 2.0*corr->BarkasCorrect    
300     G4double deltaL = lsdata->GetDeltaL(Zin, g    
301                                                   
302     res = eloss +                                 
303       CLHEP::twopi_mc2_rcl2*q2*eDensity*(delta    
304     /*                                            
305     G4cout << "  E(GeV)=" << preKinEnergy/GeV     
306      << " L= " << eloss*beta2/(twopi_mc2_rcl2*    
307      << " dL0= " << deltaL0                       
308      << " dL= " << deltaL << " dE(MeV)=" << re    
309     */                                            
310   }                                               
311   if(res > preKinEnergy || 2*res < eloss) { re    
312   /*                                              
313   G4cout << "  G4LindhardSorensenIonModel::Cor    
314          << preKinEnergy/GeV << " eloss(MeV)="    
315    << " res(MeV)=" << res << G4endl;              
316   */                                              
317   eloss = res;                                    
318 }                                                 
319                                                   
320 //....oooOO0OOooo........oooOO0OOooo........oo    
321                                                   
322 void G4LindhardSorensenIonModel::SampleSeconda    
323             std::vector<G4DynamicParticle*>* v    
324                                           cons    
325                                           cons    
326                                           G4do    
327                                           G4do    
328 {                                                 
329   G4double kineticEnergy = dp->GetKineticEnerg    
330   // take into account formfactor                 
331   const G4double tmax = MaxSecondaryEnergy(dp-    
332   const G4double minKinEnergy = std::min(cut,     
333   const G4double maxKinEnergy = std::min(maxEn    
334   if(minKinEnergy >= maxKinEnergy) { return; }    
335                                                   
336   //G4cout << "G4LindhardSorensenIonModel::Sam    
337   // << minKinEnergy << " Emax= " << maxKinEne    
338                                                   
339   G4double totEnergy = kineticEnergy + mass;      
340   G4double etot2 = totEnergy*totEnergy;           
341   G4double beta2 = kineticEnergy*(kineticEnerg    
342                                                   
343   G4double deltaKinEnergy, f;                     
344   G4double f1 = 0.0;                              
345   G4double fmax = 1.0;                            
346   if( 0.0 < spin ) { fmax += 0.5*maxKinEnergy*    
347                                                   
348   CLHEP::HepRandomEngine* rndmEngineMod = G4Ra    
349   G4double rndm[2];                               
350                                                   
351   // sampling without nuclear size effect         
352   do {                                            
353     rndmEngineMod->flatArray(2, rndm);            
354     deltaKinEnergy = minKinEnergy*maxKinEnergy    
355                     /(minKinEnergy*(1.0 - rndm    
356                                                   
357     f = 1.0 - beta2*deltaKinEnergy/tmax;          
358     if( 0.0 < spin ) {                            
359       f1 = 0.5*deltaKinEnergy*deltaKinEnergy/e    
360       f += f1;                                    
361     }                                             
362                                                   
363     // Loop checking, 03-Aug-2015, Vladimir Iv    
364   } while( fmax*rndm[1] > f);                     
365                                                   
366   // projectile formfactor - suppresion of hig    
367   // delta-electron production at high energy     
368                                                   
369   G4double x = formfact*deltaKinEnergy;           
370   if(x > 1.e-6) {                                 
371                                                   
372     G4double x1 = 1.0 + x;                        
373     G4double grej  = 1.0/(x1*x1);                 
374     if( 0.0 < spin ) {                            
375       G4double x2 = 0.5*electron_mass_c2*delta    
376       grej *= (1.0 + magMoment2*(x2 - f1/f)/(1    
377     }                                             
378     if(grej > 1.1) {                              
379       G4cout << "### G4LindhardSorensenIonMode    
380              << "  " << dp->GetDefinition()->G    
381              << " Ekin(MeV)= " <<  kineticEner    
382              << " delEkin(MeV)= " << deltaKinE    
383              << G4endl;                           
384     }                                             
385     if(rndmEngineMod->flat() > grej) { return;    
386   }                                               
387                                                   
388   G4ThreeVector deltaDirection;                   
389                                                   
390   if(UseAngularGeneratorFlag()) {                 
391                                                   
392     const G4Material* mat =  couple->GetMateri    
393     G4int Z = SelectRandomAtomNumber(mat);        
394                                                   
395     deltaDirection =                              
396       GetAngularDistribution()->SampleDirectio    
397                                                   
398   } else {                                        
399                                                   
400     G4double deltaMomentum =                      
401       std::sqrt(deltaKinEnergy * (deltaKinEner    
402     G4double cost = deltaKinEnergy * (totEnerg    
403       (deltaMomentum * dp->GetTotalMomentum())    
404     cost = std::min(cost, 1.0);                   
405     G4double sint = std::sqrt((1.0 - cost)*(1.    
406                                                   
407     G4double phi = CLHEP::twopi*rndmEngineMod-    
408                                                   
409     deltaDirection.set(sint*std::cos(phi),sint    
410     deltaDirection.rotateUz(dp->GetMomentumDir    
411   }                                               
412   /*                                              
413     G4cout << "### G4LindhardSorensenIonModel     
414            << dp->GetDefinition()->GetParticle    
415            << " Ekin(MeV)= " <<  kineticEnergy    
416            << " delEkin(MeV)= " << deltaKinEne    
417            << " tmin(MeV)= " << minKinEnergy      
418            << " tmax(MeV)= " << maxKinEnergy      
419            << " dir= " << dp->GetMomentumDirec    
420            << " dirDelta= " << deltaDirection     
421            << G4endl;                             
422   */                                              
423   // create G4DynamicParticle object for delta    
424   auto delta = new G4DynamicParticle(theElectr    
425                                                   
426   vdp->push_back(delta);                          
427                                                   
428   // Change kinematics of primary particle        
429   kineticEnergy -= deltaKinEnergy;                
430   G4ThreeVector finalP = dp->GetMomentum() - d    
431   finalP = finalP.unit();                         
432                                                   
433   fParticleChange->SetProposedKineticEnergy(ki    
434   fParticleChange->SetProposedMomentumDirectio    
435 }                                                 
436                                                   
437 //....oooOO0OOooo........oooOO0OOooo........oo    
438                                                   
439 G4double                                          
440 G4LindhardSorensenIonModel::MaxSecondaryEnergy    
441                                                   
442 {                                                 
443   // here particle type is checked for any met    
444   SetParticle(pd);                                
445   G4double tau  = kinEnergy/mass;                 
446   return 2.0*CLHEP::electron_mass_c2*tau*(tau     
447     (1. + 2.0*(tau + 1.)*eRatio + eRatio*eRati    
448 }                                                 
449                                                   
450 //....oooOO0OOooo........oooOO0OOooo........oo    
451