Geant4  10.01
G4MollerBhabhaModel.cc
Go to the documentation of this file.
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 // $Id: G4MollerBhabhaModel.cc 84397 2014-10-15 07:19:14Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4MollerBhabhaModel
34 //
35 // Author: Vladimir Ivanchenko on base of Laszlo Urban code
36 //
37 // Creation date: 03.01.2002
38 //
39 // Modifications:
40 //
41 // 13-11-02 Minor fix - use normalised direction (V.Ivanchenko)
42 // 04-12-02 Change G4DynamicParticle constructor in PostStepDoIt (V.Ivanchenko)
43 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
44 // 27-01-03 Make models region aware (V.Ivanchenko)
45 // 13-02-03 Add name (V.Ivanchenko)
46 // 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
47 // 25-07-05 Add protection in calculation of recoil direction for the case
48 // of complete energy transfer from e+ to e- (V.Ivanchenko)
49 // 06-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
50 // 15-05-06 Fix MinEnergyCut (V.Ivanchenko)
51 //
52 //
53 // Class Description:
54 //
55 // Implementation of energy loss and delta-electron production by e+/e-
56 //
57 // -------------------------------------------------------------------
58 //
59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 
62 #include "G4MollerBhabhaModel.hh"
63 #include "G4PhysicalConstants.hh"
64 #include "G4SystemOfUnits.hh"
65 #include "G4Electron.hh"
66 #include "G4Positron.hh"
67 #include "Randomize.hh"
69 #include "G4Log.hh"
70 #include "G4DeltaAngle.hh"
71 
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
73 
74 using namespace std;
75 
77  const G4String& nam)
78  : G4VEmModel(nam),
79  particle(0),
80  isElectron(true),
81  twoln10(2.0*G4Log(10.0)),
82  lowLimit(0.02*keV),
83  isInitialised(false)
84 {
86  if(p) { SetParticle(p); }
87  fParticleChange = 0;
88 }
89 
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
91 
93 {}
94 
95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
96 
98  G4double kinEnergy)
99 {
100  G4double tmax = kinEnergy;
101  if(isElectron) { tmax *= 0.5; }
102  return tmax;
103 }
104 
105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
106 
108  const G4DataVector&)
109 {
110  if(!particle) { SetParticle(p); }
111 
112  if(isInitialised) { return; }
113 
114  isInitialised = true;
118  }
119 }
120 
121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
122 
123 G4double
125  G4double kineticEnergy,
126  G4double cutEnergy,
127  G4double maxEnergy)
128 {
129  if(!particle) { SetParticle(p); }
130 
131  G4double cross = 0.0;
132  G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
133  tmax = std::min(maxEnergy, tmax);
134  //G4cout << "E= " << kineticEnergy << " cut= " << cutEnergy
135  // << " Emax= " << tmax << G4endl;
136  if(cutEnergy < tmax) {
137 
138  G4double xmin = cutEnergy/kineticEnergy;
139  G4double xmax = tmax/kineticEnergy;
140  G4double tau = kineticEnergy/electron_mass_c2;
141  G4double gam = tau + 1.0;
142  G4double gamma2= gam*gam;
143  G4double beta2 = tau*(tau + 2)/gamma2;
144 
145  //Moller (e-e-) scattering
146  if (isElectron) {
147 
148  G4double gg = (2.0*gam - 1.0)/gamma2;
149  cross = ((xmax - xmin)*(1.0 - gg + 1.0/(xmin*xmax)
150  + 1.0/((1.0-xmin)*(1.0 - xmax)))
151  - gg*G4Log( xmax*(1.0 - xmin)/(xmin*(1.0 - xmax)) ) ) / beta2;
152 
153  //Bhabha (e+e-) scattering
154  } else {
155 
156  G4double y = 1.0/(1.0 + gam);
157  G4double y2 = y*y;
158  G4double y12 = 1.0 - 2.0*y;
159  G4double b1 = 2.0 - y2;
160  G4double b2 = y12*(3.0 + y2);
161  G4double y122= y12*y12;
162  G4double b4 = y122*y12;
163  G4double b3 = b4 + y122;
164 
165  cross = (xmax - xmin)*(1.0/(beta2*xmin*xmax) + b2
166  - 0.5*b3*(xmin + xmax)
167  + b4*(xmin*xmin + xmin*xmax + xmax*xmax)/3.0)
168  - b1*G4Log(xmax/xmin);
169  }
170 
171  cross *= twopi_mc2_rcl2/kineticEnergy;
172  }
173  return cross;
174 }
175 
176 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
177 
179  const G4ParticleDefinition* p,
180  G4double kineticEnergy,
181  G4double Z, G4double,
182  G4double cutEnergy,
183  G4double maxEnergy)
184 {
185  return Z*ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
186 }
187 
188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
189 
191  const G4Material* material,
192  const G4ParticleDefinition* p,
193  G4double kinEnergy,
194  G4double cutEnergy,
195  G4double maxEnergy)
196 {
197  G4double eDensity = material->GetElectronDensity();
198  return eDensity*ComputeCrossSectionPerElectron(p,kinEnergy,cutEnergy,maxEnergy);
199  /*
200  G4double Zeff = eDensity/material->GetTotNbOfAtomsPerVolume();
201  G4double th = 0.25*sqrt(Zeff)*keV;
202  G4double cut;
203  if(isElectron) { cut = std::max(th*0.5, cutEnergy); }
204  else { cut = std::max(th, cutEnergy); }
205  G4double res = 0.0;
206  // below this threshold no bremsstrahlung
207  if (kinEnergy > th) {
208  res = eDensity*ComputeCrossSectionPerElectron(p,kinEnergy,cut,maxEnergy);
209  }
210  return res;
211  */
212 }
213 
214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
215 
217  const G4Material* material,
218  const G4ParticleDefinition* p,
219  G4double kineticEnergy,
220  G4double cut)
221 {
222  if(!particle) { SetParticle(p); }
223  // calculate the dE/dx due to the ionization by Seltzer-Berger formula
224  // checl low-energy limit
225  G4double electronDensity = material->GetElectronDensity();
226 
227  G4double Zeff = electronDensity/material->GetTotNbOfAtomsPerVolume();
228  G4double th = 0.25*sqrt(Zeff)*keV;
229  // G4double cut;
230  // if(isElectron) { cut = std::max(th*0.5, cutEnergy); }
231  // else { cut = std::max(th, cutEnergy); }
232  G4double tkin = kineticEnergy;
233  if (kineticEnergy < th) { tkin = th; }
234 
235  G4double tau = tkin/electron_mass_c2;
236  G4double gam = tau + 1.0;
237  G4double gamma2= gam*gam;
238  G4double bg2 = tau*(tau + 2);
239  G4double beta2 = bg2/gamma2;
240 
241  G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
242  eexc /= electron_mass_c2;
243  G4double eexc2 = eexc*eexc;
244 
245  G4double d = std::min(cut, MaxSecondaryEnergy(p, tkin))/electron_mass_c2;
246  G4double dedx;
247 
248  // electron
249  if (isElectron) {
250 
251  dedx = G4Log(2.0*(tau + 2.0)/eexc2) - 1.0 - beta2
252  + G4Log((tau-d)*d) + tau/(tau-d)
253  + (0.5*d*d + (2.0*tau + 1.)*G4Log(1. - d/tau))/gamma2;
254 
255  //positron
256  } else {
257 
258  G4double d2 = d*d*0.5;
259  G4double d3 = d2*d/1.5;
260  G4double d4 = d3*d*0.75;
261  G4double y = 1.0/(1.0 + gam);
262  dedx = G4Log(2.0*(tau + 2.0)/eexc2) + G4Log(tau*d)
263  - beta2*(tau + 2.0*d - y*(3.0*d2
264  + y*(d - d3 + y*(d2 - tau*d3 + d4))))/tau;
265  }
266 
267  //density correction
268  G4double x = G4Log(bg2)/twoln10;
269  dedx -= material->GetIonisation()->DensityCorrection(x);
270 
271  // now you can compute the total ionization loss
272  dedx *= twopi_mc2_rcl2*electronDensity/beta2;
273  if (dedx < 0.0) { dedx = 0.0; }
274 
275  // lowenergy extrapolation
276 
277  if (kineticEnergy < th) {
278  x = kineticEnergy/th;
279  if(x > 0.25) { dedx /= sqrt(x); }
280  else { dedx *= 1.4*sqrt(x)/(0.1 + x); }
281  }
282  return dedx;
283 }
284 
285 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
286 
287 void
288 G4MollerBhabhaModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
289  const G4MaterialCutsCouple* couple,
290  const G4DynamicParticle* dp,
291  G4double cutEnergy,
292  G4double maxEnergy)
293 {
294  G4double kineticEnergy = dp->GetKineticEnergy();
295  //const G4Material* mat = couple->GetMaterial();
296  //G4double Zeff = mat->GetElectronDensity()/mat->GetTotNbOfAtomsPerVolume();
297  // G4double th = 0.25*sqrt(Zeff)*keV;
298  G4double tmax;
299  G4double tmin = cutEnergy;
300  if(isElectron) {
301  tmax = 0.5*kineticEnergy;
302  } else {
303  tmax = kineticEnergy;
304  }
305  if(maxEnergy < tmax) { tmax = maxEnergy; }
306  if(tmin >= tmax) { return; }
307 
308  G4double energy = kineticEnergy + electron_mass_c2;
309  G4double xmin = tmin/kineticEnergy;
310  G4double xmax = tmax/kineticEnergy;
311  G4double gam = energy/electron_mass_c2;
312  G4double gamma2 = gam*gam;
313  G4double beta2 = 1.0 - 1.0/gamma2;
314  G4double x, z, q, grej;
315 
316  //Moller (e-e-) scattering
317  if (isElectron) {
318 
319  G4double gg = (2.0*gam - 1.0)/gamma2;
320  G4double y = 1.0 - xmax;
321  grej = 1.0 - gg*xmax + xmax*xmax*(1.0 - gg + (1.0 - gg*y)/(y*y));
322 
323  do {
324  q = G4UniformRand();
325  x = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
326  y = 1.0 - x;
327  z = 1.0 - gg*x + x*x*(1.0 - gg + (1.0 - gg*y)/(y*y));
328  /*
329  if(z > grej) {
330  G4cout << "G4MollerBhabhaModel::SampleSecondary Warning! "
331  << "Majorant " << grej << " < "
332  << z << " for x= " << x
333  << " e-e- scattering"
334  << G4endl;
335  }
336  */
337  } while(grej * G4UniformRand() > z);
338 
339  //Bhabha (e+e-) scattering
340  } else {
341 
342  G4double y = 1.0/(1.0 + gam);
343  G4double y2 = y*y;
344  G4double y12 = 1.0 - 2.0*y;
345  G4double b1 = 2.0 - y2;
346  G4double b2 = y12*(3.0 + y2);
347  G4double y122= y12*y12;
348  G4double b4 = y122*y12;
349  G4double b3 = b4 + y122;
350 
351  y = xmax*xmax;
352  grej = 1.0 + (y*y*b4 - xmin*xmin*xmin*b3 + y*b2 - xmin*b1)*beta2;
353  do {
354  q = G4UniformRand();
355  x = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
356  y = x*x;
357  z = 1.0 + (y*y*b4 - x*y*b3 + y*b2 - x*b1)*beta2;
358  /*
359  if(z > grej) {
360  G4cout << "G4MollerBhabhaModel::SampleSecondary Warning! "
361  << "Majorant " << grej << " < "
362  << z << " for x= " << x
363  << " e+e- scattering"
364  << G4endl;
365  }
366  */
367  } while(grej * G4UniformRand() > z);
368  }
369 
370  G4double deltaKinEnergy = x * kineticEnergy;
371 
372  G4ThreeVector deltaDirection;
373 
375  const G4Material* mat = couple->GetMaterial();
376  G4int Z = SelectRandomAtomNumber(mat);
377 
378  deltaDirection =
379  GetAngularDistribution()->SampleDirection(dp, deltaKinEnergy, Z, mat);
380 
381  } else {
382 
383  G4double deltaMomentum =
384  sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
385  G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
386  (deltaMomentum * dp->GetTotalMomentum());
387  if(cost > 1.0) { cost = 1.0; }
388  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
389 
390  G4double phi = twopi * G4UniformRand() ;
391 
392  deltaDirection.set(sint*cos(phi),sint*sin(phi), cost) ;
393  deltaDirection.rotateUz(dp->GetMomentumDirection());
394  }
395 
396  // create G4DynamicParticle object for delta ray
397  G4DynamicParticle* delta =
398  new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
399  vdp->push_back(delta);
400 
401  // primary change
402  kineticEnergy -= deltaKinEnergy;
403  G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum();
404  finalP = finalP.unit();
405 
408 }
409 
410 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
static const G4double d3
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
G4IonisParamMat * GetIonisation() const
Definition: G4Material.hh:224
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:133
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4bool isElectron(G4int ityp)
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double z
Definition: TRTMaterials.hh:39
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kineticEnergy, G4double Z, G4double A, G4double cutEnergy, G4double maxEnergy)
void SetParticle(const G4ParticleDefinition *p)
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:598
G4MollerBhabhaModel(const G4ParticleDefinition *p=0, const G4String &nam="MollerBhabha")
int G4int
Definition: G4Types.hh:78
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy)
G4double GetTotalMomentum() const
#define G4UniformRand()
Definition: Randomize.hh:95
G4ParticleChangeForLoss * fParticleChange
virtual G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, const G4Material *)=0
static const G4double b3
G4double GetElectronDensity() const
Definition: G4Material.hh:215
G4bool UseAngularGeneratorFlag() const
Definition: G4VEmModel.hh:678
const G4ParticleDefinition * particle
const G4ThreeVector & GetMomentumDirection() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
static const G4double b2
void SetProposedMomentumDirection(const G4ThreeVector &dir)
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:207
G4ParticleDefinition * theElectron
virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy)
G4double DensityCorrection(G4double x)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4double energy(const ThreeVector &p, const G4double m)
virtual G4double ComputeCrossSectionPerElectron(const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:605
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static const G4double b1
G4double GetMeanExcitationEnergy() const
static G4Electron * Electron()
Definition: G4Electron.cc:94
static const G4double twoln10
static const double keV
Definition: G4SIunits.hh:195
double G4double
Definition: G4Types.hh:76
static const G4double d2
static const G4double b4
static const G4double d4
G4ThreeVector GetMomentum() const
const G4Material * GetMaterial() const
G4int SelectRandomAtomNumber(const G4Material *)
Definition: G4VEmModel.hh:545