ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MollerBhabhaModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4MollerBhabhaModel.cc
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(p) { SetParticle(p); }
86  fParticleChange = nullptr;
87 }
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90 
91 G4MollerBhabhaModel::~G4MollerBhabhaModel()
92 {}
93 
94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
95 
96 G4double G4MollerBhabhaModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
97  G4double kinEnergy)
98 {
99  G4double tmax = kinEnergy;
100  if(isElectron) { tmax *= 0.5; }
101  return tmax;
102 }
103 
104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
105 
106 void G4MollerBhabhaModel::Initialise(const G4ParticleDefinition* p,
107  const G4DataVector&)
108 {
109  if(!particle) { SetParticle(p); }
110 
111  if(isInitialised) { return; }
112 
113  isInitialised = true;
114  fParticleChange = GetParticleChangeForLoss();
115  if(UseAngularGeneratorFlag() && !GetAngularDistribution()) {
116  SetAngularDistribution(new G4DeltaAngle());
117  }
118 }
119 
120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
121 
122 G4double
123 G4MollerBhabhaModel::ComputeCrossSectionPerElectron(const G4ParticleDefinition* p,
124  G4double kineticEnergy,
125  G4double cutEnergy,
126  G4double maxEnergy)
127 {
128  if(!particle) { SetParticle(p); }
129 
130  G4double cross = 0.0;
131  G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
132  tmax = std::min(maxEnergy, tmax);
133  //G4cout << "E= " << kineticEnergy << " cut= " << cutEnergy
134  // << " Emax= " << tmax << G4endl;
135  if(cutEnergy < tmax) {
136 
137  G4double xmin = cutEnergy/kineticEnergy;
138  G4double xmax = tmax/kineticEnergy;
139  G4double tau = kineticEnergy/electron_mass_c2;
140  G4double gam = tau + 1.0;
141  G4double gamma2= gam*gam;
142  G4double beta2 = tau*(tau + 2)/gamma2;
143 
144  //Moller (e-e-) scattering
145  if (isElectron) {
146 
147  G4double gg = (2.0*gam - 1.0)/gamma2;
148  cross = ((xmax - xmin)*(1.0 - gg + 1.0/(xmin*xmax)
149  + 1.0/((1.0-xmin)*(1.0 - xmax)))
150  - gg*G4Log( xmax*(1.0 - xmin)/(xmin*(1.0 - xmax)) ) ) / beta2;
151 
152  //Bhabha (e+e-) scattering
153  } else {
154 
155  G4double y = 1.0/(1.0 + gam);
156  G4double y2 = y*y;
157  G4double y12 = 1.0 - 2.0*y;
158  G4double b1 = 2.0 - y2;
159  G4double b2 = y12*(3.0 + y2);
160  G4double y122= y12*y12;
161  G4double b4 = y122*y12;
162  G4double b3 = b4 + y122;
163 
164  cross = (xmax - xmin)*(1.0/(beta2*xmin*xmax) + b2
165  - 0.5*b3*(xmin + xmax)
166  + b4*(xmin*xmin + xmin*xmax + xmax*xmax)/3.0)
167  - b1*G4Log(xmax/xmin);
168  }
169 
170  cross *= twopi_mc2_rcl2/kineticEnergy;
171  }
172  return cross;
173 }
174 
175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
176 
177 G4double G4MollerBhabhaModel::ComputeCrossSectionPerAtom(
178  const G4ParticleDefinition* p,
179  G4double kineticEnergy,
180  G4double Z, G4double,
181  G4double cutEnergy,
182  G4double maxEnergy)
183 {
184  return Z*ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
185 }
186 
187 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
188 
189 G4double G4MollerBhabhaModel::CrossSectionPerVolume(
190  const G4Material* material,
191  const G4ParticleDefinition* p,
192  G4double kinEnergy,
193  G4double cutEnergy,
194  G4double maxEnergy)
195 {
196  G4double eDensity = material->GetElectronDensity();
197  return eDensity*ComputeCrossSectionPerElectron(p,kinEnergy,cutEnergy,maxEnergy);
198 }
199 
200 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
201 
202 G4double G4MollerBhabhaModel::ComputeDEDXPerVolume(
203  const G4Material* material,
204  const G4ParticleDefinition* p,
205  G4double kineticEnergy,
206  G4double cut)
207 {
208  if(nullptr == particle) { SetParticle(p); }
209  // calculate the dE/dx due to the ionization by Seltzer-Berger formula
210  // checl low-energy limit
211  G4double electronDensity = material->GetElectronDensity();
212 
213  G4double Zeff = material->GetIonisation()->GetZeffective();
214  G4double th = 0.25*sqrt(Zeff)*keV;
215  G4double tkin = std::max(kineticEnergy, th);
216 
217  G4double tau = tkin/electron_mass_c2;
218  G4double gam = tau + 1.0;
219  G4double gamma2= gam*gam;
220  G4double bg2 = tau*(tau + 2);
221  G4double beta2 = bg2/gamma2;
222 
223  G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
224  eexc /= electron_mass_c2;
225  G4double eexc2 = eexc*eexc;
226 
227  G4double d = std::min(cut, MaxSecondaryEnergy(p, tkin))/electron_mass_c2;
228  G4double dedx;
229 
230  // electron
231  if (isElectron) {
232 
233  dedx = G4Log(2.0*(tau + 2.0)/eexc2) - 1.0 - beta2
234  + G4Log((tau-d)*d) + tau/(tau-d)
235  + (0.5*d*d + (2.0*tau + 1.)*G4Log(1. - d/tau))/gamma2;
236 
237  //positron
238  } else {
239 
240  G4double d2 = d*d*0.5;
241  G4double d3 = d2*d/1.5;
242  G4double d4 = d3*d*0.75;
243  G4double y = 1.0/(1.0 + gam);
244  dedx = G4Log(2.0*(tau + 2.0)/eexc2) + G4Log(tau*d)
245  - beta2*(tau + 2.0*d - y*(3.0*d2
246  + y*(d - d3 + y*(d2 - tau*d3 + d4))))/tau;
247  }
248 
249  //density correction
250  G4double x = G4Log(bg2)/twoln10;
251  dedx -= material->GetIonisation()->DensityCorrection(x);
252 
253  // now you can compute the total ionization loss
254  dedx *= twopi_mc2_rcl2*electronDensity/beta2;
255  if (dedx < 0.0) { dedx = 0.0; }
256 
257  // lowenergy extrapolation
258 
259  if (kineticEnergy < th) {
260  x = kineticEnergy/th;
261  if(x > 0.25) { dedx /= sqrt(x); }
262  else { dedx *= 1.4*sqrt(x)/(0.1 + x); }
263  }
264  return dedx;
265 }
266 
267 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
268 
269 void
270 G4MollerBhabhaModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
271  const G4MaterialCutsCouple* couple,
272  const G4DynamicParticle* dp,
273  G4double cutEnergy,
274  G4double maxEnergy)
275 {
276  G4double kineticEnergy = dp->GetKineticEnergy();
277  //G4cout << "G4MollerBhabhaModel::SampleSecondaries: E= " << kineticEnergy
278  // << " in " << couple->GetMaterial()->GetName() << G4endl;
279  G4double tmax;
280  G4double tmin = cutEnergy;
281  if(isElectron) {
282  tmax = 0.5*kineticEnergy;
283  } else {
284  tmax = kineticEnergy;
285  }
286  if(maxEnergy < tmax) { tmax = maxEnergy; }
287  if(tmin >= tmax) { return; }
288 
289  G4double energy = kineticEnergy + electron_mass_c2;
290  G4double xmin = tmin/kineticEnergy;
291  G4double xmax = tmax/kineticEnergy;
292  G4double gam = energy/electron_mass_c2;
293  G4double gamma2 = gam*gam;
294  G4double beta2 = 1.0 - 1.0/gamma2;
295  G4double x, z, grej;
296  CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
297  G4double rndm[2];
298 
299  //Moller (e-e-) scattering
300  if (isElectron) {
301 
302  G4double gg = (2.0*gam - 1.0)/gamma2;
303  G4double y = 1.0 - xmax;
304  grej = 1.0 - gg*xmax + xmax*xmax*(1.0 - gg + (1.0 - gg*y)/(y*y));
305 
306  do {
307  rndmEngine->flatArray(2, rndm);
308  x = xmin*xmax/(xmin*(1.0 - rndm[0]) + xmax*rndm[0]);
309  y = 1.0 - x;
310  z = 1.0 - gg*x + x*x*(1.0 - gg + (1.0 - gg*y)/(y*y));
311  /*
312  if(z > grej) {
313  G4cout << "G4MollerBhabhaModel::SampleSecondary Warning! "
314  << "Majorant " << grej << " < "
315  << z << " for x= " << x
316  << " e-e- scattering"
317  << G4endl;
318  }
319  */
320  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
321  } while(grej * rndm[1] > z);
322 
323  //Bhabha (e+e-) scattering
324  } else {
325 
326  G4double y = 1.0/(1.0 + gam);
327  G4double y2 = y*y;
328  G4double y12 = 1.0 - 2.0*y;
329  G4double b1 = 2.0 - y2;
330  G4double b2 = y12*(3.0 + y2);
331  G4double y122= y12*y12;
332  G4double b4 = y122*y12;
333  G4double b3 = b4 + y122;
334 
335  y = xmax*xmax;
336  grej = 1.0 + (y*y*b4 - xmin*xmin*xmin*b3 + y*b2 - xmin*b1)*beta2;
337  do {
338  rndmEngine->flatArray(2, rndm);
339  x = xmin*xmax/(xmin*(1.0 - rndm[0]) + xmax*rndm[0]);
340  y = x*x;
341  z = 1.0 + (y*y*b4 - x*y*b3 + y*b2 - x*b1)*beta2;
342  /*
343  if(z > grej) {
344  G4cout << "G4MollerBhabhaModel::SampleSecondary Warning! "
345  << "Majorant " << grej << " < "
346  << z << " for x= " << x
347  << " e+e- scattering"
348  << G4endl;
349  }
350  */
351  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
352  } while(grej * rndm[1] > z);
353  }
354 
355  G4double deltaKinEnergy = x * kineticEnergy;
356 
357  G4ThreeVector deltaDirection;
358 
359  if(UseAngularGeneratorFlag()) {
360  const G4Material* mat = couple->GetMaterial();
361  G4int Z = SelectRandomAtomNumber(mat);
362 
363  deltaDirection =
364  GetAngularDistribution()->SampleDirection(dp, deltaKinEnergy, Z, mat);
365 
366  } else {
367 
368  G4double deltaMomentum =
369  sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
370  G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
371  (deltaMomentum * dp->GetTotalMomentum());
372  if(cost > 1.0) { cost = 1.0; }
373  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
374 
375  G4double phi = twopi * rndmEngine->flat() ;
376 
377  deltaDirection.set(sint*cos(phi),sint*sin(phi), cost) ;
378  deltaDirection.rotateUz(dp->GetMomentumDirection());
379  }
380 
381  // create G4DynamicParticle object for delta ray
382  G4DynamicParticle* delta =
383  new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
384  vdp->push_back(delta);
385 
386  // primary change
387  kineticEnergy -= deltaKinEnergy;
388  G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum();
389  finalP = finalP.unit();
390 
391  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
392  fParticleChange->SetProposedMomentumDirection(finalP);
393 }
394 
395 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......