ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4BetheHeitlerModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4BetheHeitlerModel.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: G4BetheHeitlerModel
33 //
34 // Author: Vladimir Ivanchenko on base of Michel Maire code
35 //
36 // Creation date: 15.03.2005
37 //
38 // Modifications by Vladimir Ivanchenko, Michel Maire, Mihaly Novak
39 //
40 // Class Description:
41 //
42 // -------------------------------------------------------------------
43 //
44 
45 #include "G4BetheHeitlerModel.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "G4Electron.hh"
49 #include "G4Positron.hh"
50 #include "G4Gamma.hh"
51 #include "Randomize.hh"
53 #include "G4Pow.hh"
54 #include "G4Exp.hh"
55 #include "G4ModifiedTsai.hh"
56 
58 std::vector<G4BetheHeitlerModel::ElementData*> G4BetheHeitlerModel::gElementData;
59 
61  const G4String& nam)
62 : G4VEmModel(nam),
63  fG4Calc(G4Pow::GetInstance()), fTheGamma(G4Gamma::Gamma()),
64  fTheElectron(G4Electron::Electron()), fThePositron(G4Positron::Positron()),
65  fParticleChange(nullptr)
66 {
68 }
69 
71 {
72  if (IsMaster()) {
73  // clear ElementData container
74  for (size_t iz = 0; iz < gElementData.size(); ++iz) {
75  if (gElementData[iz]) delete gElementData[iz];
76  }
77  gElementData.clear();
78  }
79 }
80 
82  const G4DataVector& cuts)
83 {
84  if (IsMaster()) {
86  }
88  if (IsMaster()) {
90  }
91 }
92 
94  G4VEmModel* masterModel)
95 {
97 }
98 
99 // Calculates the microscopic cross section in GEANT4 internal units.
100 // A parametrized formula from L. Urban is used to estimate
101 // the total cross section.
102 // It gives a good description of the data from 1.5 MeV to 100 GeV.
103 // below 1.5 MeV: sigma=sigma(1.5MeV)*(GammaEnergy-2electronmass)
104 // *(GammaEnergy-2electronmass)
105 G4double
107  G4double gammaEnergy, G4double Z,
109 {
110  G4double xSection = 0.0 ;
111  // short versions
112  static const G4double kMC2 = CLHEP::electron_mass_c2;
113  // zero cross section below the kinematical limit: Eg<2mc^2
114  if (Z < 0.9 || gammaEnergy <= 2.0*kMC2) { return xSection; }
115  //
116  static const G4double gammaEnergyLimit = 1.5*CLHEP::MeV;
117  // set coefficients a, b c
118  static const G4double a0 = 8.7842e+2*CLHEP::microbarn;
119  static const G4double a1 = -1.9625e+3*CLHEP::microbarn;
120  static const G4double a2 = 1.2949e+3*CLHEP::microbarn;
121  static const G4double a3 = -2.0028e+2*CLHEP::microbarn;
122  static const G4double a4 = 1.2575e+1*CLHEP::microbarn;
123  static const G4double a5 = -2.8333e-1*CLHEP::microbarn;
124 
125  static const G4double b0 = -1.0342e+1*CLHEP::microbarn;
126  static const G4double b1 = 1.7692e+1*CLHEP::microbarn;
127  static const G4double b2 = -8.2381 *CLHEP::microbarn;
128  static const G4double b3 = 1.3063 *CLHEP::microbarn;
129  static const G4double b4 = -9.0815e-2*CLHEP::microbarn;
130  static const G4double b5 = 2.3586e-3*CLHEP::microbarn;
131 
132  static const G4double c0 = -4.5263e+2*CLHEP::microbarn;
133  static const G4double c1 = 1.1161e+3*CLHEP::microbarn;
134  static const G4double c2 = -8.6749e+2*CLHEP::microbarn;
135  static const G4double c3 = 2.1773e+2*CLHEP::microbarn;
136  static const G4double c4 = -2.0467e+1*CLHEP::microbarn;
137  static const G4double c5 = 6.5372e-1*CLHEP::microbarn;
138  // check low energy limit of the approximation (1.5 MeV)
139  G4double gammaEnergyOrg = gammaEnergy;
140  if (gammaEnergy < gammaEnergyLimit) { gammaEnergy = gammaEnergyLimit; }
141  // compute gamma energy variables
142  const G4double x = G4Log(gammaEnergy/kMC2);
143  const G4double x2 = x *x;
144  const G4double x3 = x2*x;
145  const G4double x4 = x3*x;
146  const G4double x5 = x4*x;
147  //
148  const G4double F1 = a0 + a1*x + a2*x2 + a3*x3 + a4*x4 + a5*x5;
149  const G4double F2 = b0 + b1*x + b2*x2 + b3*x3 + b4*x4 + b5*x5;
150  const G4double F3 = c0 + c1*x + c2*x2 + c3*x3 + c4*x4 + c5*x5;
151  // compute the approximated cross section
152  xSection = (Z + 1.)*(F1*Z + F2*Z*Z + F3);
153  // check if we are below the limit of the approximation and apply correction
154  if (gammaEnergyOrg < gammaEnergyLimit) {
155  const G4double dum = (gammaEnergyOrg-2.*kMC2)/(gammaEnergyLimit-2.*kMC2);
156  xSection *= dum*dum;
157  }
158  // make sure that the cross section is never negative
159  xSection = std::max(xSection, 0.);
160  return xSection;
161 }
162 
163 // The secondaries e+e- energies are sampled using the Bethe - Heitler
164 // cross sections with Coulomb correction.
165 // A modified version of the random number techniques of Butcher & Messel
166 // is used (Nuc Phys 20(1960),15).
167 //
168 // GEANT4 internal units.
169 //
170 // Note 1 : Effects due to the breakdown of the Born approximation at
171 // low energy are ignored.
172 // Note 2 : The differential cross section implicitly takes account of
173 // pair creation in both nuclear and atomic electron fields.
174 // However triplet prodution is not generated.
175 void G4BetheHeitlerModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
176  const G4MaterialCutsCouple* couple,
177  const G4DynamicParticle* aDynamicGamma,
179 {
180  // set some constant values
181  const G4double gammaEnergy = aDynamicGamma->GetKineticEnergy();
182  const G4double eps0 = CLHEP::electron_mass_c2/gammaEnergy;
183  //
184  // check kinematical limit: gamma energy(Eg) must be at least 2 e- rest mass
185  if (eps0 > 0.5) { return; }
186  //
187  // select target element of the material (probs. are based on partial x-secs)
188  const G4Element* anElement = SelectTargetAtom(couple, fTheGamma, gammaEnergy,
189  aDynamicGamma->GetLogKineticEnergy());
190 
191  //
192  // get the random engine
193  CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
194  //
195  // 'eps' is the total energy transferred to one of the e-/e+ pair in initial
196  // gamma energy units Eg. Since the corresponding DCS is symmetric on eps=0.5,
197  // the kinematical limits for eps0=mc^2/Eg <= eps <= 0.5
198  // 1. 'eps' is sampled uniformly on the [eps0, 0.5] inteval if Eg<Egsmall
199  // 2. otherwise, on the [eps_min, 0.5] interval according to the DCS (case 2.)
200  G4double eps;
201  // case 1.
202  static const G4double Egsmall = 2.*CLHEP::MeV;
203  if (gammaEnergy < Egsmall) {
204  eps = eps0 + (0.5-eps0)*rndmEngine->flat();
205  } else {
206  // case 2.
207  // get the Coulomb factor for the target element (Z) and gamma energy (Eg)
208  // F(Z) = 8*ln(Z)/3 if Eg <= 50 [MeV] => no Coulomb correction
209  // F(Z) = 8*ln(Z)/3 + 8*fc(Z) if Eg > 50 [MeV] => fc(Z) is the Coulomb cor.
210  //
211  // The screening variable 'delta(eps)' = 136*Z^{-1/3}*eps0/[eps(1-eps)]
212  // Due to the Coulomb correction, the DCS can go below zero even at
213  // kinematicaly allowed eps > eps0 values. In order to exclude this eps
214  // range with negative DCS, the minimum eps value will be set to eps_min =
215  // max[eps0, epsp] with epsp is the solution of SF(delta(epsp)) - F(Z)/2 = 0
216  // with SF being the screening function (SF1=SF2 at high value of delta).
217  // The solution is epsp = 0.5 - 0.5*sqrt[ 1 - 4*136*Z^{-1/3}eps0/deltap]
218  // with deltap = Exp[(42.038-F(Z))/8.29]-0.958. So the limits are:
219  // - when eps=eps_max = 0.5 => delta_min = 136*Z^{-1/3}*eps0/4
220  // - epsp = 0.5 - 0.5*sqrt[ 1 - delta_min/deltap]
221  // - and eps_min = max[eps0, epsp]
222  static const G4double midEnergy = 50.*CLHEP::MeV;
223  const G4int iZet = std::min(gMaxZet, anElement->GetZasInt());
224  const G4double deltaFactor = 136.*eps0/anElement->GetIonisation()->GetZ3();
225  G4double deltaMax = gElementData[iZet]->fDeltaMaxLow;
226  G4double FZ = 8.*anElement->GetIonisation()->GetlogZ3();
227  if (gammaEnergy > midEnergy) {
228  FZ += 8.*(anElement->GetfCoulomb());
229  deltaMax = gElementData[iZet]->fDeltaMaxHigh;
230  }
231  const G4double deltaMin = 4.*deltaFactor;
232  //
233  // compute the limits of eps
234  const G4double epsp = 0.5 - 0.5*std::sqrt(1. - deltaMin/deltaMax) ;
235  const G4double epsMin = std::max(eps0,epsp);
236  const G4double epsRange = 0.5 - epsMin;
237  //
238  // sample the energy rate (eps) of the created electron (or positron)
239  G4double F10, F20;
240  ScreenFunction12(deltaMin, F10, F20);
241  F10 -= FZ;
242  F20 -= FZ;
243  const G4double NormF1 = std::max(F10 * epsRange * epsRange, 0.);
244  const G4double NormF2 = std::max(1.5 * F20 , 0.);
245  const G4double NormCond = NormF1/(NormF1 + NormF2);
246  // we will need 3 uniform random number for each trial of sampling
247  G4double rndmv[3];
248  G4double greject = 0.;
249  do {
250  rndmEngine->flatArray(3, rndmv);
251  if (NormCond > rndmv[0]) {
252  eps = 0.5 - epsRange * fG4Calc->A13(rndmv[1]);
253  const G4double delta = deltaFactor/(eps*(1.-eps));
254  greject = (ScreenFunction1(delta)-FZ)/F10;
255  } else {
256  eps = epsMin + epsRange*rndmv[1];
257  const G4double delta = deltaFactor/(eps*(1.-eps));
258  greject = (ScreenFunction2(delta)-FZ)/F20;
259  }
260  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
261  } while (greject < rndmv[2]);
262  } // end of eps sampling
263  //
264  // select charges randomly
265  G4double eTotEnergy, pTotEnergy;
266  if (rndmEngine->flat() > 0.5) {
267  eTotEnergy = (1.-eps)*gammaEnergy;
268  pTotEnergy = eps*gammaEnergy;
269  } else {
270  pTotEnergy = (1.-eps)*gammaEnergy;
271  eTotEnergy = eps*gammaEnergy;
272  }
273  //
274  // sample pair kinematics
275  const G4double eKinEnergy = std::max(0.,eTotEnergy - CLHEP::electron_mass_c2);
276  const G4double pKinEnergy = std::max(0.,pTotEnergy - CLHEP::electron_mass_c2);
277  //
278  G4ThreeVector eDirection, pDirection;
279  //
281  eKinEnergy, pKinEnergy,
282  eDirection, pDirection);
283  // create G4DynamicParticle object for the particle1
284  G4DynamicParticle* aParticle1= new G4DynamicParticle(
285  fTheElectron,eDirection,eKinEnergy);
286  // create G4DynamicParticle object for the particle2
287  G4DynamicParticle* aParticle2= new G4DynamicParticle(
288  fThePositron,pDirection,pKinEnergy);
289  // Fill output vector
290  fvect->push_back(aParticle1);
291  fvect->push_back(aParticle2);
292  // kill incident photon
295 }
296 
297 // should be called only by the master and at initialisation
299 {
300  G4int size = gElementData.size();
301  if (size < gMaxZet+1) {
302  gElementData.resize(gMaxZet+1, nullptr);
303  }
304  // create for all elements that are in the detector
305  const G4ElementTable* elemTable = G4Element::GetElementTable();
306  size_t numElems = (*elemTable).size();
307  for (size_t ie = 0; ie < numElems; ++ie) {
308  const G4Element* elem = (*elemTable)[ie];
309  const G4int iz = std::min(gMaxZet, elem->GetZasInt());
310  if (!gElementData[iz]) { // create it if doesn't exist yet
311  G4double FZLow = 8.*elem->GetIonisation()->GetlogZ3();
312  G4double FZHigh = FZLow + 8.*elem->GetfCoulomb();
313  ElementData* elD = new ElementData();
314  elD->fDeltaMaxLow = G4Exp((42.038 - FZLow )/8.29) - 0.958;
315  elD->fDeltaMaxHigh = G4Exp((42.038 - FZHigh)/8.29) - 0.958;
316  gElementData[iz] = elD;
317  }
318  }
319 }
320