Geant4 Cross Reference

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


** Warning: Cannot open xref database.

  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:   G4MollerBhabhaModel               
 33 //                                                
 34 // Author:        Vladimir Ivanchenko on base     
 35 //                                                
 36 // Creation date: 03.01.2002                      
 37 //                                                
 38 // Modifications:                                 
 39 //                                                
 40 // 13-11-02 Minor fix - use normalised directi    
 41 // 04-12-02 Change G4DynamicParticle construct    
 42 // 23-12-02 Change interface in order to move     
 43 // 27-01-03 Make models region aware (V.Ivanch    
 44 // 13-02-03 Add name (V.Ivanchenko)               
 45 // 08-04-05 Major optimisation of internal int    
 46 // 25-07-05 Add protection in calculation of r    
 47 //          of complete energy transfer from e    
 48 // 06-02-06 ComputeCrossSectionPerElectron, Co    
 49 // 15-05-06 Fix MinEnergyCut (V.Ivanchenko)       
 50 //                                                
 51 //                                                
 52 // Class Description:                             
 53 //                                                
 54 // Implementation of energy loss and delta-ele    
 55 //                                                
 56 // -------------------------------------------    
 57 //                                                
 58 //....oooOO0OOooo........oooOO0OOooo........oo    
 59 //....oooOO0OOooo........oooOO0OOooo........oo    
 60                                                   
 61 #include "G4MollerBhabhaModel.hh"                 
 62 #include "G4PhysicalConstants.hh"                 
 63 #include "G4SystemOfUnits.hh"                     
 64 #include "G4Electron.hh"                          
 65 #include "G4Positron.hh"                          
 66 #include "Randomize.hh"                           
 67 #include "G4ParticleChangeForLoss.hh"             
 68 #include "G4Log.hh"                               
 69 #include "G4DeltaAngle.hh"                        
 70                                                   
 71 //....oooOO0OOooo........oooOO0OOooo........oo    
 72                                                   
 73 using namespace std;                              
 74                                                   
 75 G4MollerBhabhaModel::G4MollerBhabhaModel(const    
 76                                          const    
 77   : G4VEmModel(nam),                              
 78     particle(nullptr),                            
 79     isElectron(true),                             
 80     twoln10(2.0*G4Log(10.0)),                     
 81     lowLimit(0.02*keV),                           
 82     isInitialised(false)                          
 83 {                                                 
 84   theElectron = G4Electron::Electron();           
 85   if(nullptr != p) { SetParticle(p); }            
 86   fParticleChange = nullptr;                      
 87 }                                                 
 88                                                   
 89 //....oooOO0OOooo........oooOO0OOooo........oo    
 90                                                   
 91 G4MollerBhabhaModel::~G4MollerBhabhaModel() =     
 92                                                   
 93 //....oooOO0OOooo........oooOO0OOooo........oo    
 94                                                   
 95 G4double G4MollerBhabhaModel::MaxSecondaryEner    
 96                                                   
 97 {                                                 
 98   G4double tmax = kinEnergy;                      
 99   if(isElectron) { tmax *= 0.5; }                 
100   return tmax;                                    
101 }                                                 
102                                                   
103 //....oooOO0OOooo........oooOO0OOooo........oo    
104                                                   
105 void G4MollerBhabhaModel::Initialise(const G4P    
106                                      const G4D    
107 {                                                 
108   if(p != particle) { SetParticle(p); }           
109                                                   
110   if(isInitialised) { return; }                   
111                                                   
112   isInitialised = true;                           
113   fParticleChange = GetParticleChangeForLoss()    
114   if(UseAngularGeneratorFlag() && !GetAngularD    
115     SetAngularDistribution(new G4DeltaAngle())    
116   }                                               
117 }                                                 
118                                                   
119 //....oooOO0OOooo........oooOO0OOooo........oo    
120                                                   
121 G4double G4MollerBhabhaModel::ComputeCrossSect    
122          const G4ParticleDefinition* p, G4doub    
123    G4double cutEnergy, G4double maxEnergy)        
124 {                                                 
125   if(p != particle) { SetParticle(p); }           
126                                                   
127   G4double cross = 0.0;                           
128   G4double tmax = MaxSecondaryEnergy(p, kineti    
129   tmax = std::min(maxEnergy, tmax);               
130   //G4cout << "E= " << kineticEnergy << " cut=    
131   //         << " Emax= " << tmax << G4endl;      
132   if(cutEnergy < tmax) {                          
133                                                   
134     G4double xmin  = cutEnergy/kineticEnergy;     
135     G4double xmax  = tmax/kineticEnergy;          
136     G4double tau   = kineticEnergy/electron_ma    
137     G4double gam   = tau + 1.0;                   
138     G4double gamma2= gam*gam;                     
139     G4double beta2 = tau*(tau + 2)/gamma2;        
140                                                   
141     //Moller (e-e-) scattering                    
142     if (isElectron) {                             
143                                                   
144       G4double gg = (2.0*gam - 1.0)/gamma2;       
145       cross = ((xmax - xmin)*(1.0 - gg + 1.0/(    
146                               + 1.0/((1.0-xmin    
147             - gg*G4Log( xmax*(1.0 - xmin)/(xmi    
148                                                   
149     //Bhabha (e+e-) scattering                    
150     } else {                                      
151                                                   
152       G4double y   = 1.0/(1.0 + gam);             
153       G4double y2  = y*y;                         
154       G4double y12 = 1.0 - 2.0*y;                 
155       G4double b1  = 2.0 - y2;                    
156       G4double b2  = y12*(3.0 + y2);              
157       G4double y122= y12*y12;                     
158       G4double b4  = y122*y12;                    
159       G4double b3  = b4 + y122;                   
160                                                   
161       cross = (xmax - xmin)*(1.0/(beta2*xmin*x    
162             - 0.5*b3*(xmin + xmax)                
163             + b4*(xmin*xmin + xmin*xmax + xmax    
164             - b1*G4Log(xmax/xmin);                
165     }                                             
166                                                   
167     cross *= twopi_mc2_rcl2/kineticEnergy;        
168   }                                               
169   return cross;                                   
170 }                                                 
171                                                   
172 //....oooOO0OOooo........oooOO0OOooo........oo    
173                                                   
174 G4double G4MollerBhabhaModel::ComputeCrossSect    
175                                            con    
176                                                   
177                                                   
178                                                   
179                                                   
180 {                                                 
181   return Z*ComputeCrossSectionPerElectron(p,ki    
182 }                                                 
183                                                   
184 //....oooOO0OOooo........oooOO0OOooo........oo    
185                                                   
186 G4double G4MollerBhabhaModel::CrossSectionPerV    
187                                            con    
188                                            con    
189                                                   
190                                                   
191                                                   
192 {                                                 
193   G4double eDensity = material->GetElectronDen    
194   return eDensity*ComputeCrossSectionPerElectr    
195 }                                                 
196                                                   
197 //....oooOO0OOooo........oooOO0OOooo........oo    
198                                                   
199 G4double G4MollerBhabhaModel::ComputeDEDXPerVo    
200                                           cons    
201                                           cons    
202                                                   
203                                                   
204 {                                                 
205   if(p != particle) { SetParticle(p); }           
206   // calculate the dE/dx due to the ionization    
207   // checl low-energy limit                       
208   G4double electronDensity = material->GetElec    
209                                                   
210   G4double Zeff  = material->GetIonisation()->    
211   G4double th    = 0.25*sqrt(Zeff)*keV;           
212   G4double tkin = std::max(kineticEnergy, th);    
213                                                   
214   G4double tau   = tkin/electron_mass_c2;         
215   G4double gam   = tau + 1.0;                     
216   G4double gamma2= gam*gam;                       
217   G4double bg2   = tau*(tau + 2);                 
218   G4double beta2 = bg2/gamma2;                    
219                                                   
220   G4double eexc  = material->GetIonisation()->    
221   eexc          /= electron_mass_c2;              
222   G4double eexc2 = eexc*eexc;                     
223                                                   
224   G4double d = std::min(cut, MaxSecondaryEnerg    
225   G4double dedx;                                  
226                                                   
227   // electron                                     
228   if (isElectron) {                               
229                                                   
230     dedx = G4Log(2.0*(tau + 2.0)/eexc2) - 1.0     
231          + G4Log((tau-d)*d) + tau/(tau-d)         
232          + (0.5*d*d + (2.0*tau + 1.)*G4Log(1.     
233                                                   
234   //positron                                      
235   } else {                                        
236                                                   
237     G4double d2 = d*d*0.5;                        
238     G4double d3 = d2*d/1.5;                       
239     G4double d4 = d3*d*0.75;                      
240     G4double y  = 1.0/(1.0 + gam);                
241     dedx = G4Log(2.0*(tau + 2.0)/eexc2) + G4Lo    
242          - beta2*(tau + 2.0*d - y*(3.0*d2         
243          + y*(d - d3 + y*(d2 - tau*d3 + d4))))    
244   }                                               
245                                                   
246   //density correction                            
247   G4double x = G4Log(bg2)/twoln10;                
248   dedx -= material->GetIonisation()->DensityCo    
249                                                   
250   // now you can compute the total ionization     
251   dedx *= twopi_mc2_rcl2*electronDensity/beta2    
252   if (dedx < 0.0) { dedx = 0.0; }                 
253                                                   
254   // lowenergy extrapolation                      
255                                                   
256   if (kineticEnergy < th) {                       
257     x = kineticEnergy/th;                         
258     if(x > 0.25) { dedx /= sqrt(x); }             
259     else { dedx *= 1.4*sqrt(x)/(0.1 + x); }       
260   }                                               
261   return dedx;                                    
262 }                                                 
263                                                   
264 //....oooOO0OOooo........oooOO0OOooo........oo    
265                                                   
266 void                                              
267 G4MollerBhabhaModel::SampleSecondaries(std::ve    
268                                        const G    
269                                        const G    
270                                        G4doubl    
271                                        G4doubl    
272 {                                                 
273   G4double kineticEnergy = dp->GetKineticEnerg    
274   //G4cout << "G4MollerBhabhaModel::SampleSeco    
275   //       << " in " << couple->GetMaterial()-    
276   G4double tmax;                                  
277   G4double tmin = cutEnergy;                      
278   if(isElectron) {                                
279     tmax = 0.5*kineticEnergy;                     
280   } else {                                        
281     tmax = kineticEnergy;                         
282   }                                               
283   if(maxEnergy < tmax) { tmax = maxEnergy; }      
284   if(tmin >= tmax) { return; }                    
285                                                   
286   G4double energy = kineticEnergy + electron_m    
287   G4double xmin   = tmin/kineticEnergy;           
288   G4double xmax   = tmax/kineticEnergy;           
289   G4double gam    = energy/electron_mass_c2;      
290   G4double gamma2 = gam*gam;                      
291   G4double beta2  = 1.0 - 1.0/gamma2;             
292   G4double x, z, grej;                            
293   CLHEP::HepRandomEngine* rndmEngine = G4Rando    
294   G4double rndm[2];                               
295                                                   
296   //Moller (e-e-) scattering                      
297   if (isElectron) {                               
298                                                   
299     G4double gg = (2.0*gam - 1.0)/gamma2;         
300     G4double y = 1.0 - xmax;                      
301     grej = 1.0 - gg*xmax + xmax*xmax*(1.0 - gg    
302                                                   
303     do {                                          
304       rndmEngine->flatArray(2, rndm);             
305       x = xmin*xmax/(xmin*(1.0 - rndm[0]) + xm    
306       y = 1.0 - x;                                
307       z = 1.0 - gg*x + x*x*(1.0 - gg + (1.0 -     
308       /*                                          
309       if(z > grej) {                              
310         G4cout << "G4MollerBhabhaModel::Sample    
311                << "Majorant " << grej << " < "    
312                << z << " for x= " << x            
313                << " e-e- scattering"              
314                << G4endl;                         
315       }                                           
316       */                                          
317       // Loop checking, 03-Aug-2015, Vladimir     
318     } while(grej * rndm[1] > z);                  
319                                                   
320   //Bhabha (e+e-) scattering                      
321   } else {                                        
322                                                   
323     G4double y   = 1.0/(1.0 + gam);               
324     G4double y2  = y*y;                           
325     G4double y12 = 1.0 - 2.0*y;                   
326     G4double b1  = 2.0 - y2;                      
327     G4double b2  = y12*(3.0 + y2);                
328     G4double y122= y12*y12;                       
329     G4double b4  = y122*y12;                      
330     G4double b3  = b4 + y122;                     
331                                                   
332     y    = xmax*xmax;                             
333     grej = 1.0 + (y*y*b4 - xmin*xmin*xmin*b3 +    
334     do {                                          
335       rndmEngine->flatArray(2, rndm);             
336       x = xmin*xmax/(xmin*(1.0 - rndm[0]) + xm    
337       y = x*x;                                    
338       z = 1.0 + (y*y*b4 - x*y*b3 + y*b2 - x*b1    
339       /*                                          
340       if(z > grej) {                              
341         G4cout << "G4MollerBhabhaModel::Sample    
342                << "Majorant " << grej << " < "    
343                << z << " for x= " << x            
344                << " e+e- scattering"              
345                << G4endl;                         
346       }                                           
347       */                                          
348       // Loop checking, 03-Aug-2015, Vladimir     
349     } while(grej * rndm[1] > z);                  
350   }                                               
351                                                   
352   G4double deltaKinEnergy = x * kineticEnergy;    
353                                                   
354   G4ThreeVector deltaDirection;                   
355                                                   
356   if(UseAngularGeneratorFlag()) {                 
357     const G4Material* mat =  couple->GetMateri    
358     G4int Z = SelectRandomAtomNumber(mat);        
359                                                   
360     deltaDirection =                              
361       GetAngularDistribution()->SampleDirectio    
362                                                   
363   } else {                                        
364                                                   
365     G4double deltaMomentum =                      
366       sqrt(deltaKinEnergy * (deltaKinEnergy +     
367     G4double cost = deltaKinEnergy * (energy +    
368       (deltaMomentum * dp->GetTotalMomentum())    
369     if(cost > 1.0) { cost = 1.0; }                
370     G4double sint = sqrt((1.0 - cost)*(1.0 + c    
371                                                   
372     G4double phi = twopi * rndmEngine->flat()     
373                                                   
374     deltaDirection.set(sint*cos(phi),sint*sin(    
375     deltaDirection.rotateUz(dp->GetMomentumDir    
376   }                                               
377                                                   
378   // create G4DynamicParticle object for delta    
379   auto delta = new G4DynamicParticle(theElectr    
380   vdp->push_back(delta);                          
381                                                   
382   // primary change                               
383   kineticEnergy -= deltaKinEnergy;                
384   G4ThreeVector finalP = dp->GetMomentum() - d    
385   finalP               = finalP.unit();           
386                                                   
387   fParticleChange->SetProposedKineticEnergy(ki    
388   fParticleChange->SetProposedMomentumDirectio    
389 }                                                 
390                                                   
391 //....oooOO0OOooo........oooOO0OOooo........oo    
392