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 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 //
 27 // -------------------------------------------------------------------
 28 //
 29 // GEANT4 Class file
 30 //
 31 //
 32 // File name:   G4MollerBhabhaModel
 33 //
 34 // Author:        Vladimir Ivanchenko on base of Laszlo Urban code
 35 //
 36 // Creation date: 03.01.2002
 37 //
 38 // Modifications:
 39 //
 40 // 13-11-02 Minor fix - use normalised direction (V.Ivanchenko)
 41 // 04-12-02 Change G4DynamicParticle constructor in PostStepDoIt (V.Ivanchenko)
 42 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
 43 // 27-01-03 Make models region aware (V.Ivanchenko)
 44 // 13-02-03 Add name (V.Ivanchenko)
 45 // 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
 46 // 25-07-05 Add protection in calculation of recoil direction for the case 
 47 //          of complete energy transfer from e+ to e- (V.Ivanchenko)
 48 // 06-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
 49 // 15-05-06 Fix MinEnergyCut (V.Ivanchenko)
 50 //
 51 //
 52 // Class Description:
 53 //
 54 // Implementation of energy loss and delta-electron production by e+/e-
 55 //
 56 // -------------------------------------------------------------------
 57 //
 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 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........oooOO0OOooo........oooOO0OOooo......
 72 
 73 using namespace std;
 74 
 75 G4MollerBhabhaModel::G4MollerBhabhaModel(const G4ParticleDefinition* p,
 76                                          const G4String& nam)
 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........oooOO0OOooo........oooOO0OOooo......
 90 
 91 G4MollerBhabhaModel::~G4MollerBhabhaModel() = default;
 92 
 93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 94 
 95 G4double G4MollerBhabhaModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
 96                                                  G4double kinEnergy) 
 97 {
 98   G4double tmax = kinEnergy;
 99   if(isElectron) { tmax *= 0.5; }
100   return tmax;
101 }
102 
103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
104 
105 void G4MollerBhabhaModel::Initialise(const G4ParticleDefinition* p,
106                                      const G4DataVector&)
107 {
108   if(p != particle) { SetParticle(p); }
109 
110   if(isInitialised) { return; }
111 
112   isInitialised = true;
113   fParticleChange = GetParticleChangeForLoss();
114   if(UseAngularGeneratorFlag() && !GetAngularDistribution()) {
115     SetAngularDistribution(new G4DeltaAngle());
116   }
117 }
118 
119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
120 
121 G4double G4MollerBhabhaModel::ComputeCrossSectionPerElectron(
122          const G4ParticleDefinition* p, G4double kineticEnergy,
123    G4double cutEnergy, G4double maxEnergy)
124 {
125   if(p != particle) { SetParticle(p); }
126 
127   G4double cross = 0.0;
128   G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
129   tmax = std::min(maxEnergy, tmax);
130   //G4cout << "E= " << kineticEnergy << " cut= " << cutEnergy 
131   //         << " Emax= " << tmax << G4endl;
132   if(cutEnergy < tmax) {
133 
134     G4double xmin  = cutEnergy/kineticEnergy;
135     G4double xmax  = tmax/kineticEnergy;
136     G4double tau   = kineticEnergy/electron_mass_c2;
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/(xmin*xmax)
146                               + 1.0/((1.0-xmin)*(1.0 - xmax)))
147             - gg*G4Log( xmax*(1.0 - xmin)/(xmin*(1.0 - xmax)) ) ) / beta2;
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*xmax) + b2
162             - 0.5*b3*(xmin + xmax)
163             + b4*(xmin*xmin + xmin*xmax + xmax*xmax)/3.0)
164             - b1*G4Log(xmax/xmin);
165     }
166 
167     cross *= twopi_mc2_rcl2/kineticEnergy;
168   }
169   return cross;
170 }
171 
172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
173 
174 G4double G4MollerBhabhaModel::ComputeCrossSectionPerAtom(
175                                            const G4ParticleDefinition* p,
176                                                  G4double kineticEnergy,
177                                                  G4double Z, G4double,
178                                                  G4double cutEnergy,
179                                                  G4double maxEnergy)
180 {
181   return Z*ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
182 }
183 
184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
185 
186 G4double G4MollerBhabhaModel::CrossSectionPerVolume(
187                                            const G4Material* material,
188                                            const G4ParticleDefinition* p,
189                                                  G4double kinEnergy,
190                                                  G4double cutEnergy,
191                                                  G4double maxEnergy)
192 {
193   G4double eDensity = material->GetElectronDensity();
194   return eDensity*ComputeCrossSectionPerElectron(p,kinEnergy,cutEnergy,maxEnergy);
195 }
196 
197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
198 
199 G4double G4MollerBhabhaModel::ComputeDEDXPerVolume(
200                                           const G4Material* material,
201                                           const G4ParticleDefinition* p,
202                                                 G4double kineticEnergy,
203                                                 G4double cut)
204 {
205   if(p != particle) { SetParticle(p); }
206   // calculate the dE/dx due to the ionization by Seltzer-Berger formula
207   // checl low-energy limit
208   G4double electronDensity = material->GetElectronDensity();
209   
210   G4double Zeff  = material->GetIonisation()->GetZeffective();
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()->GetMeanExcitationEnergy();
221   eexc          /= electron_mass_c2;
222   G4double eexc2 = eexc*eexc; 
223 
224   G4double d = std::min(cut, MaxSecondaryEnergy(p, tkin))/electron_mass_c2;
225   G4double dedx;
226 
227   // electron
228   if (isElectron) {
229 
230     dedx = G4Log(2.0*(tau + 2.0)/eexc2) - 1.0 - beta2
231          + G4Log((tau-d)*d) + tau/(tau-d)
232          + (0.5*d*d + (2.0*tau + 1.)*G4Log(1. - d/tau))/gamma2;
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) + G4Log(tau*d)
242          - beta2*(tau + 2.0*d - y*(3.0*d2 
243          + y*(d - d3 + y*(d2 - tau*d3 + d4))))/tau;
244   } 
245 
246   //density correction 
247   G4double x = G4Log(bg2)/twoln10;
248   dedx -= material->GetIonisation()->DensityCorrection(x); 
249 
250   // now you can compute the total ionization loss
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........oooOO0OOooo........oooOO0OOooo......
265 
266 void 
267 G4MollerBhabhaModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
268                                        const G4MaterialCutsCouple* couple,
269                                        const G4DynamicParticle* dp,
270                                        G4double cutEnergy,
271                                        G4double maxEnergy)
272 {
273   G4double kineticEnergy = dp->GetKineticEnergy();
274   //G4cout << "G4MollerBhabhaModel::SampleSecondaries: E= " << kineticEnergy
275   //       << " in " << couple->GetMaterial()->GetName() << G4endl;
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_mass_c2;
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 = G4Random::getTheEngine();
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 + (1.0 - gg*y)/(y*y));
302 
303     do {
304       rndmEngine->flatArray(2, rndm);
305       x = xmin*xmax/(xmin*(1.0 - rndm[0]) + xmax*rndm[0]);
306       y = 1.0 - x;
307       z = 1.0 - gg*x + x*x*(1.0 - gg + (1.0 - gg*y)/(y*y));
308       /*
309       if(z > grej) {
310         G4cout << "G4MollerBhabhaModel::SampleSecondary Warning! "
311                << "Majorant " << grej << " < "
312                << z << " for x= " << x
313                << " e-e- scattering"
314                << G4endl;
315       }
316       */
317       // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
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 + y*b2 - xmin*b1)*beta2; 
334     do {
335       rndmEngine->flatArray(2, rndm);
336       x = xmin*xmax/(xmin*(1.0 - rndm[0]) + xmax*rndm[0]);
337       y = x*x;
338       z = 1.0 + (y*y*b4 - x*y*b3 + y*b2 - x*b1)*beta2; 
339       /*
340       if(z > grej) {
341         G4cout << "G4MollerBhabhaModel::SampleSecondary Warning! "
342                << "Majorant " << grej << " < "
343                << z << " for x= " << x
344                << " e+e- scattering"
345                << G4endl;
346       }
347       */
348       // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
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->GetMaterial();
358     G4int Z = SelectRandomAtomNumber(mat);
359 
360     deltaDirection = 
361       GetAngularDistribution()->SampleDirection(dp, deltaKinEnergy, Z, mat);
362 
363   } else {
364  
365     G4double deltaMomentum =
366       sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
367     G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
368       (deltaMomentum * dp->GetTotalMomentum());
369     if(cost > 1.0) { cost = 1.0; }
370     G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
371 
372     G4double phi = twopi * rndmEngine->flat() ;
373 
374     deltaDirection.set(sint*cos(phi),sint*sin(phi), cost) ;
375     deltaDirection.rotateUz(dp->GetMomentumDirection());
376   }  
377 
378   // create G4DynamicParticle object for delta ray
379   auto delta = new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
380   vdp->push_back(delta);
381 
382   // primary change
383   kineticEnergy -= deltaKinEnergy;
384   G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum();
385   finalP               = finalP.unit();
386 
387   fParticleChange->SetProposedKineticEnergy(kineticEnergy);
388   fParticleChange->SetProposedMomentumDirection(finalP);
389 }
390 
391 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
392