Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/adjoint/src/G4AdjointhIonisationModel.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/adjoint/src/G4AdjointhIonisationModel.cc (Version 11.3.0) and /processes/electromagnetic/adjoint/src/G4AdjointhIonisationModel.cc (Version 4.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 #include "G4AdjointhIonisationModel.hh"           
 28                                                   
 29 #include "G4AdjointCSManager.hh"                  
 30 #include "G4AdjointElectron.hh"                   
 31 #include "G4AdjointProton.hh"                     
 32 #include "G4BetheBlochModel.hh"                   
 33 #include "G4BraggModel.hh"                        
 34 #include "G4NistManager.hh"                       
 35 #include "G4ParticleChange.hh"                    
 36 #include "G4PhysicalConstants.hh"                 
 37 #include "G4Proton.hh"                            
 38 #include "G4SystemOfUnits.hh"                     
 39 #include "G4TrackStatus.hh"                       
 40                                                   
 41 //////////////////////////////////////////////    
 42 G4AdjointhIonisationModel::G4AdjointhIonisatio    
 43   : G4VEmAdjointModel("Adjoint_hIonisation")      
 44 {                                                 
 45   fUseMatrix               = true;                
 46   fUseMatrixPerElement     = true;                
 47   fApplyCutInRange         = true;                
 48   fOneMatrixForAllElements = true;                
 49   fSecondPartSameType      = false;               
 50                                                   
 51   // The direct EM Model is taken as BetheBloc    
 52   // computation of the differential cross sec    
 53   // The Bragg model could be used as an alter    
 54   // differential cross section                   
 55                                                   
 56   fDirectModel              = new G4BetheBloch    
 57   fBraggDirectEMModel       = new G4BraggModel    
 58   fAdjEquivDirectSecondPart = G4AdjointElectro    
 59   fDirectPrimaryPart        = pDef;               
 60                                                   
 61   if(pDef == G4Proton::Proton())                  
 62   {                                               
 63     fAdjEquivDirectPrimPart = G4AdjointProton:    
 64   }                                               
 65                                                   
 66   DefineProjectileProperty();                     
 67 }                                                 
 68                                                   
 69 //////////////////////////////////////////////    
 70 G4AdjointhIonisationModel::~G4AdjointhIonisati    
 71                                                   
 72 //////////////////////////////////////////////    
 73 void G4AdjointhIonisationModel::SampleSecondar    
 74   const G4Track& aTrack, G4bool isScatProjToPr    
 75   G4ParticleChange* fParticleChange)              
 76 {                                                 
 77   if(!fUseMatrix)                                 
 78     return RapidSampleSecondaries(aTrack, isSc    
 79                                                   
 80   const G4DynamicParticle* theAdjointPrimary =    
 81                                                   
 82   // Elastic inverse scattering                   
 83   G4double adjointPrimKinEnergy = theAdjointPr    
 84   G4double adjointPrimP         = theAdjointPr    
 85                                                   
 86   if(adjointPrimKinEnergy > GetHighEnergyLimit    
 87   {                                               
 88     return;                                       
 89   }                                               
 90                                                   
 91   // Sample secondary energy                      
 92   G4double projectileKinEnergy =                  
 93     SampleAdjSecEnergyFromCSMatrix(adjointPrim    
 94   CorrectPostStepWeight(fParticleChange, aTrac    
 95                         adjointPrimKinEnergy,     
 96                         isScatProjToProj);        
 97   // Caution!!! this weight correction should     
 98                                                   
 99   // Kinematic:                                   
100   // we consider a two body elastic scattering    
101   // the projectile knock on an e- at rest and    
102   G4double projectileM0          = fAdjEquivDi    
103   G4double projectileTotalEnergy = projectileM    
104   G4double projectileP2 =                         
105     projectileTotalEnergy * projectileTotalEne    
106                                                   
107   // Companion                                    
108   G4double companionM0 = fAdjEquivDirectPrimPa    
109   if(isScatProjToProj)                            
110   {                                               
111     companionM0 = fAdjEquivDirectSecondPart->G    
112   }                                               
113   G4double companionTotalEnergy =                 
114     companionM0 + projectileKinEnergy - adjoin    
115   G4double companionP2 =                          
116     companionTotalEnergy * companionTotalEnerg    
117                                                   
118   // Projectile momentum                          
119   G4double P_parallel =                           
120     (adjointPrimP * adjointPrimP + projectileP    
121     (2. * adjointPrimP);                          
122   G4double P_perp = std::sqrt(projectileP2 - P    
123   G4ThreeVector dir_parallel = theAdjointPrima    
124   G4double phi               = G4UniformRand()    
125   G4ThreeVector projectileMomentum =              
126     G4ThreeVector(P_perp * std::cos(phi), P_pe    
127   projectileMomentum.rotateUz(dir_parallel);      
128                                                   
129   if(!isScatProjToProj)                           
130   {  // kill the primary and add a secondary      
131     fParticleChange->ProposeTrackStatus(fStopA    
132     fParticleChange->AddSecondary(                
133       new G4DynamicParticle(fAdjEquivDirectPri    
134   }                                               
135   else                                            
136   {                                               
137     fParticleChange->ProposeEnergy(projectileK    
138     fParticleChange->ProposeMomentumDirection(    
139   }                                               
140 }                                                 
141                                                   
142 //////////////////////////////////////////////    
143 void G4AdjointhIonisationModel::RapidSampleSec    
144   const G4Track& aTrack, G4bool isScatProjToPr    
145   G4ParticleChange* fParticleChange)              
146 {                                                 
147   const G4DynamicParticle* theAdjointPrimary =    
148   DefineCurrentMaterial(aTrack.GetMaterialCuts    
149                                                   
150   G4double adjointPrimKinEnergy = theAdjointPr    
151   G4double adjointPrimP         = theAdjointPr    
152                                                   
153   if(adjointPrimKinEnergy > GetHighEnergyLimit    
154   {                                               
155     return;                                       
156   }                                               
157                                                   
158   G4double projectileKinEnergy = 0.;              
159   G4double eEnergy             = 0.;              
160   G4double newCS =                                
161     fCurrentMaterial->GetElectronDensity() * t    
162   if(!isScatProjToProj)                           
163   {  // 1/E^2 distribution                        
164                                                   
165     eEnergy       = adjointPrimKinEnergy;         
166     G4double Emax = GetSecondAdjEnergyMaxForPr    
167     G4double Emin = GetSecondAdjEnergyMinForPr    
168     if(Emin >= Emax)                              
169       return;                                     
170     G4double a = 1. / Emax;                       
171     G4double b = 1. / Emin;                       
172     newCS      = newCS * (b - a) / eEnergy;       
173                                                   
174     projectileKinEnergy = 1. / (b - (b - a) *     
175   }                                               
176   else                                            
177   {                                               
178     G4double Emax =                               
179       GetSecondAdjEnergyMaxForScatProjToProj(a    
180     G4double Emin =                               
181       GetSecondAdjEnergyMinForScatProjToProj(a    
182     if(Emin >= Emax)                              
183       return;                                     
184     G4double diff1 = Emin - adjointPrimKinEner    
185     G4double diff2 = Emax - adjointPrimKinEner    
186                                                   
187     G4double t1    = adjointPrimKinEnergy * (1    
188     G4double t2    = adjointPrimKinEnergy * (1    
189     G4double t3    = 2. * std::log(Emax / Emin    
190     G4double sum_t = t1 + t2 + t3;                
191     newCS      = newCS * sum_t / adjointPrimKi    
192     G4double t = G4UniformRand() * sum_t;         
193     if(t <= t1)                                   
194     {                                             
195       G4double q          = G4UniformRand() *     
196       projectileKinEnergy = adjointPrimKinEner    
197     }                                             
198     else if(t <= t2)                              
199     {                                             
200       G4double q          = G4UniformRand() *     
201       projectileKinEnergy = 1. / (1. / Emin -     
202     }                                             
203     else                                          
204     {                                             
205       projectileKinEnergy = Emin * std::pow(Em    
206     }                                             
207     eEnergy = projectileKinEnergy - adjointPri    
208   }                                               
209                                                   
210   G4double diffCS_perAtom_Used = twopi_mc2_rcl    
211                                  projectileKin    
212                                  eEnergy / eEn    
213                                                   
214   // Weight correction                            
215   // First w_corr is set to the ratio between     
216   G4double w_corr =                               
217     G4AdjointCSManager::GetAdjointCSManager()-    
218                                                   
219   w_corr *= newCS / fLastCS;                      
220   // Then another correction is needed due to     
221   // differential CS has been used rather than    
222   // direct model. Here we consider the true d    
223   // numerical differentiation over Tcut of th    
224                                                   
225   G4double diffCS =                               
226     DiffCrossSectionPerAtomPrimToSecond(projec    
227   w_corr *= diffCS / diffCS_perAtom_Used;         
228                                                   
229   if (isScatProjToProj && fTcutSecond>0.005) w    
230                                                   
231   G4double new_weight = aTrack.GetWeight() * w    
232   fParticleChange->SetParentWeightByProcess(fa    
233   fParticleChange->SetSecondaryWeightByProcess    
234   fParticleChange->ProposeParentWeight(new_wei    
235                                                   
236   // Kinematic:                                   
237   // we consider a two body elastic scattering    
238   // the projectile knocks on an e- at rest an    
239   G4double projectileM0          = fAdjEquivDi    
240   G4double projectileTotalEnergy = projectileM    
241   G4double projectileP2 =                         
242     projectileTotalEnergy * projectileTotalEne    
243                                                   
244   // Companion                                    
245   G4double companionM0 = fAdjEquivDirectPrimPa    
246   if(isScatProjToProj)                            
247   {                                               
248     companionM0 = fAdjEquivDirectSecondPart->G    
249   }                                               
250   G4double companionTotalEnergy =                 
251     companionM0 + projectileKinEnergy - adjoin    
252   G4double companionP2 =                          
253     companionTotalEnergy * companionTotalEnerg    
254                                                   
255   // Projectile momentum                          
256   G4double P_parallel =                           
257     (adjointPrimP * adjointPrimP + projectileP    
258     (2. * adjointPrimP);                          
259   G4double P_perp = std::sqrt(projectileP2 - P    
260   G4ThreeVector dir_parallel = theAdjointPrima    
261   G4double phi               = G4UniformRand()    
262   G4ThreeVector projectileMomentum =              
263     G4ThreeVector(P_perp * std::cos(phi), P_pe    
264   projectileMomentum.rotateUz(dir_parallel);      
265                                                   
266   if(!isScatProjToProj)                           
267   {  // kill the primary and add a secondary      
268     fParticleChange->ProposeTrackStatus(fStopA    
269     fParticleChange->AddSecondary(                
270       new G4DynamicParticle(fAdjEquivDirectPri    
271   }                                               
272   else                                            
273   {                                               
274     fParticleChange->ProposeEnergy(projectileK    
275     fParticleChange->ProposeMomentumDirection(    
276   }                                               
277 }                                                 
278                                                   
279 //////////////////////////////////////////////    
280 G4double G4AdjointhIonisationModel::DiffCrossS    
281   G4double kinEnergyProj, G4double kinEnergyPr    
282 {  // Probably here the Bragg Model should be     
283    // kinEnergyProj/nuc < 2 MeV                   
284                                                   
285   G4double dSigmadEprod = 0.;                     
286   G4double Emax_proj    = GetSecondAdjEnergyMa    
287   G4double Emin_proj    = GetSecondAdjEnergyMi    
288                                                   
289   // the produced particle should have a kinet    
290   // projectile                                   
291   if(kinEnergyProj > Emin_proj && kinEnergyPro    
292   {                                               
293     G4double Tmax = kinEnergyProj;                
294     G4double E1   = kinEnergyProd;                
295     //1.0006 factor seems to give the best dif    
296     G4double E2   = kinEnergyProd *1.0006;        
297     G4double sigma1, sigma2;                      
298     if(kinEnergyProj > 2. * MeV)                  
299     {                                             
300       sigma1 = fDirectModel->ComputeCrossSecti    
301         fDirectPrimaryPart, kinEnergyProj, Z,     
302       sigma2 = fDirectModel->ComputeCrossSecti    
303         fDirectPrimaryPart, kinEnergyProj, Z,     
304     }                                             
305     else                                          
306     {                                             
307       sigma1 = fBraggDirectEMModel->ComputeCro    
308         fDirectPrimaryPart, kinEnergyProj, Z,     
309       sigma2 = fBraggDirectEMModel->ComputeCro    
310         fDirectPrimaryPart, kinEnergyProj, Z,     
311     }                                             
312                                                   
313     dSigmadEprod = (sigma1 - sigma2) / (E2 - E    
314     if(dSigmadEprod > 1.)                         
315     {                                             
316       G4cout << "sigma1 " << kinEnergyProj / M    
317              << '\t' << sigma1 << G4endl;         
318       G4cout << "sigma2 " << kinEnergyProj / M    
319              << '\t' << sigma2 << G4endl;         
320       G4cout << "dsigma " << kinEnergyProj / M    
321              << '\t' << dSigmadEprod << G4endl    
322     }                                             
323                                                   
324     // correction of differential cross sectio    
325     // the suppression of particle at secondar    
326     // Bloch Model. This correction consists o    
327     // function used to test the rejection of     
328     // from G4BetheBlochModel::SampleSecondari    
329     G4double deltaKinEnergy = kinEnergyProd;      
330                                                   
331     // projectile formfactor - suppression of     
332     // delta-electron production at high energ    
333     G4double x = fFormFact * deltaKinEnergy;      
334     if(x > 1.e-6)                                 
335     {                                             
336       G4double totEnergy = kinEnergyProj + fMa    
337       G4double etot2     = totEnergy * totEner    
338       G4double beta2 = kinEnergyProj * (kinEne    
339       G4double f     = 1.0 - beta2 * deltaKinE    
340       G4double f1    = 0.0;                       
341       if(0.5 == fSpin)                            
342       {                                           
343         f1 = 0.5 * deltaKinEnergy * deltaKinEn    
344         f += f1;                                  
345       }                                           
346       G4double x1 = 1.0 + x;                      
347       G4double gg = 1.0 / (x1 * x1);              
348       if(0.5 == fSpin)                            
349       {                                           
350         G4double x2 = 0.5 * electron_mass_c2 *    
351         gg *= (1.0 + fMagMoment2 * (x2 - f1 /     
352       }                                           
353       if(gg > 1.0)                                
354       {                                           
355         G4cout << "### G4BetheBlochModel in Ad    
356                << G4endl;                         
357         gg = 1.;                                  
358       }                                           
359       dSigmadEprod *= gg;                         
360     }                                             
361   }                                               
362                                                   
363   return dSigmadEprod;                            
364 }                                                 
365                                                   
366 //////////////////////////////////////////////    
367 void G4AdjointhIonisationModel::DefineProjecti    
368 {                                                 
369   // Slightly modified code taken from G4Bethe    
370   G4String pname = fDirectPrimaryPart->GetPart    
371                                                   
372   fMass           = fDirectPrimaryPart->GetPDG    
373   fSpin           = fDirectPrimaryPart->GetPDG    
374   fMassRatio      = electron_mass_c2 / fMass;     
375   fOnePlusRatio2  = (1. + fMassRatio) * (1. +     
376   fOneMinusRatio2 = (1. - fMassRatio) * (1. -     
377   G4double magmom = fDirectPrimaryPart->GetPDG    
378                     (0.5 * eplus * hbar_Planck    
379   fMagMoment2 = magmom * magmom - 1.0;            
380   fFormFact   = 0.0;                              
381   if(fDirectPrimaryPart->GetLeptonNumber() ==     
382   {                                               
383     G4double x = 0.8426 * GeV;                    
384     if(fSpin == 0.0 && fMass < GeV)               
385     {                                             
386       x = 0.736 * GeV;                            
387     }                                             
388     else if(fMass > GeV)                          
389     {                                             
390       x /= G4NistManager::Instance()->GetZ13(f    
391     }                                             
392     fFormFact = 2.0 * electron_mass_c2 / (x *     
393   }                                               
394 }                                                 
395                                                   
396 //////////////////////////////////////////////    
397 G4double G4AdjointhIonisationModel::AdjointCro    
398   const G4MaterialCutsCouple* aCouple, G4doubl    
399   G4bool isScatProjToProj)                        
400 {                                                 
401   if(fUseMatrix)                                  
402     return G4VEmAdjointModel::AdjointCrossSect    
403                                                   
404   DefineCurrentMaterial(aCouple);                 
405                                                   
406   G4double Cross =                                
407     fCurrentMaterial->GetElectronDensity() * t    
408                                                   
409   if(!isScatProjToProj)                           
410   {                                               
411     G4double Emax_proj = GetSecondAdjEnergyMax    
412     G4double Emin_proj = GetSecondAdjEnergyMin    
413     if(Emax_proj > Emin_proj && primEnergy > f    
414     {                                             
415       Cross *= (1. / Emin_proj - 1. / Emax_pro    
416     }                                             
417     else                                          
418       Cross = 0.;                                 
419   }                                               
420   else                                            
421   {                                               
422     G4double Emax_proj = GetSecondAdjEnergyMax    
423     G4double Emin_proj =                          
424       GetSecondAdjEnergyMinForScatProjToProj(p    
425     G4double diff1 = Emin_proj - primEnergy;      
426     G4double diff2 = Emax_proj - primEnergy;      
427     G4double t1 =                                 
428       (1. / diff1 + 1. / Emin_proj - 1. / diff    
429     G4double t2 =                                 
430       2. * std::log(Emax_proj / Emin_proj) / p    
431     Cross *= (t1 + t2);                           
432   }                                               
433   fLastCS = Cross;                                
434   return Cross;                                   
435 }                                                 
436                                                   
437 //////////////////////////////////////////////    
438 G4double G4AdjointhIonisationModel::GetSecondA    
439   G4double primAdjEnergy)                         
440 {                                                 
441   G4double Tmax = primAdjEnergy * fOnePlusRati    
442                   (fOneMinusRatio2 - 2. * fMas    
443   return Tmax;                                    
444 }                                                 
445                                                   
446 //////////////////////////////////////////////    
447 G4double G4AdjointhIonisationModel::GetSecondA    
448   G4double primAdjEnergy, G4double tcut)          
449 {                                                 
450   return primAdjEnergy + tcut;                    
451 }                                                 
452                                                   
453 //////////////////////////////////////////////    
454 G4double G4AdjointhIonisationModel::GetSecondA    
455 {                                                 
456   return GetHighEnergyLimit();                    
457 }                                                 
458                                                   
459 //////////////////////////////////////////////    
460 G4double G4AdjointhIonisationModel::GetSecondA    
461   G4double primAdjEnergy)                         
462 {                                                 
463   G4double Tmin =                                 
464     (2. * primAdjEnergy - 4. * fMass +            
465      std::sqrt(4. * primAdjEnergy * primAdjEne    
466                8. * primAdjEnergy * fMass * (1    
467     4.;                                           
468   return Tmin;                                    
469 }                                                 
470