ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PairProductionRelModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PairProductionRelModel.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: G4PairProductionRelModel
33 //
34 // Author: Andreas Schaelicke
35 //
36 // Creation date: 02.04.2009
37 //
38 // Modifications:
39 // 20.03.17 Change LPMconstant such that it gives suppression variable 's'
40 // that consistent to Migdal's one; fix a small bug in 'logTS1'
41 // computation; suppression is consistent now with the one in the
42 // brem. model (F.Hariri)
43 // 28-05-18 New version with improved screening function approximation, improved
44 // LPM function approximation, efficiency, documentation and cleanup.
45 // Corrected call to selecting target atom in the final state sampling.
46 // (M. Novak)
47 //
48 // Class Description:
49 //
50 // Main References:
51 // J.W.Motz et.al., Rev. Mod. Phys. 41 (1969) 581.
52 // S.Klein, Rev. Mod. Phys. 71 (1999) 1501.
53 // T.Stanev et.al., Phys. Rev. D25 (1982) 1291.
54 // M.L.Ter-Mikaelian, High-energy Electromagnetic Processes in Condensed Media,
55 // Wiley, 1972.
56 //
57 // -------------------------------------------------------------------
58 
60 #include "G4PhysicalConstants.hh"
61 #include "G4SystemOfUnits.hh"
62 #include "G4Gamma.hh"
63 #include "G4Electron.hh"
64 #include "G4Positron.hh"
66 #include "G4LossTableManager.hh"
67 #include "G4ModifiedTsai.hh"
68 
69 
71 
72 // LPM constant: \alpha(mc^2)^2/(4\pi*\hbar c)
75  /(4.*CLHEP::pi*CLHEP::hbarc);
76 
77 // abscissas and weights of an 8 point Gauss-Legendre quadrature
78 // for numerical integration on [0,1]
80  1.98550718e-02, 1.01666761e-01, 2.37233795e-01, 4.08282679e-01,
81  5.91717321e-01, 7.62766205e-01, 8.98333239e-01, 9.80144928e-01
82 };
84  5.06142681e-02, 1.11190517e-01, 1.56853323e-01, 1.81341892e-01,
85  1.81341892e-01, 1.56853323e-01, 1.11190517e-01, 5.06142681e-02
86 };
87 
88 // elastic and inelatic radiation logarithms for light elements (where the
89 // Thomas-Fermi model doesn't work): computed by using Dirac-Fock model of atom.
91  0.0, 5.3104, 4.7935, 4.7402, 4.7112, 4.6694, 4.6134, 4.5520
92 };
94  0.0, 5.9173, 5.6125, 5.5377, 5.4728, 5.4174, 5.3688, 5.3236
95 };
96 
97 // constant cross section factor
101 
102 // gamma energy limit above which LPM suppression will be applied (if the
103 // fIsUseLPMCorrection flag is true)
105 
106 // special data structure per element i.e. per Z
107 std::vector<G4PairProductionRelModel::ElementData*> G4PairProductionRelModel::gElementData;
108 
109 // LPM supression functions evaluated at initialisation time
111 
112 // CTR
114  const G4String& nam)
115  : G4VEmModel(nam), fIsUseLPMCorrection(true), fIsUseCompleteScreening(false),
116  fLPMEnergy(0.), fG4Calc(G4Pow::GetInstance()), fTheGamma(G4Gamma::Gamma()),
117  fTheElectron(G4Electron::Electron()), fThePositron(G4Positron::Positron()),
118  fParticleChange(nullptr)
119 {
120  // gamma energy below which the parametrized atomic x-section is used (80 GeV)
122  // gamma energy below the Coulomb correction is turned off (50 MeV)
124  // set angular generator used in the final state kinematics computation
126 }
127 
128 // DTR
130 {
131  if (IsMaster()) {
132  // clear ElementData container
133  for (size_t iz = 0; iz < gElementData.size(); ++iz) {
134  if (gElementData[iz]) delete gElementData[iz];
135  }
136  gElementData.clear();
137  // clear LPMFunctions (if any)
138  if (fIsUseLPMCorrection) {
139  gLPMFuncs.fLPMFuncG.clear();
140  gLPMFuncs.fLPMFuncPhi.clear();
141  gLPMFuncs.fIsInitialized = false;
142  }
143  }
144 }
145 
147  const G4DataVector& cuts)
148 {
149  if (IsMaster()) {
150  // init element data and LPM funcs
151  if (IsMaster()) {
153  if (fIsUseLPMCorrection) {
155  }
156  }
157  }
159  if(IsMaster() && LowEnergyLimit() < HighEnergyLimit()) {
161  }
162 }
163 
165  G4VEmModel* masterModel)
166 {
167  if(LowEnergyLimit() < HighEnergyLimit()) {
169  }
170 }
171 
173  G4double Z)
174 {
175  G4double xSection = 0.0;
176  // check if LPM suppression needs to be used
177  const G4bool isLPM = (fIsUseLPMCorrection && gammaEnergy>gEgLPMActivation);
178  // determine the kinematical limits (taken into account the correction due to
179  // the way in which the Coulomb correction is applied i.e. avoid negative DCS)
180  const G4int iz = std::min(gMaxZet, G4lrint(Z));
181  const G4double eps0 = CLHEP::electron_mass_c2/gammaEnergy;
182  // Coulomb correction is always included in the DCS even below 50 MeV (note:
183  // that this DCS is only used to get the integrated x-section)
184  const G4double dmax = gElementData[iz]->fDeltaMaxHigh;
185  const G4double dmin = 4.*eps0*gElementData[iz]->fDeltaFactor;
186  const G4double eps1 = 0.5 - 0.5*std::sqrt(1.-dmin/dmax);
187  const G4double epsMin = std::max(eps0, eps1);
188  const G4double epsMax = 0.5; // DCS is symmetric around eps=0.5
189  // let Et be the total energy transferred to the e- or to the e+
190  // the [Et-min, Et-max] interval will be divided into i=1,2,..,n subintervals
191  // with width of dInterv = (Et-max - Et-min)/n and numerical integration will
192  // be done in each sub-inteval using the xi = (Et - Et_i-min)/dInterv variable
193  // that is in [0,1]. The 8-point GL q. is used for the integration on [0,1].
194  const G4int numSub = 2;
195  const G4double dInterv= (epsMax - epsMin)*gammaEnergy/G4double(numSub);
196  G4double minEti = epsMin*gammaEnergy; // Et-min i.e. Et_0-min
197  for (G4int i = 0; i < numSub; ++i) {
198  for (G4int ngl = 0; ngl < 8; ++ngl) {
199  const G4double Et = (minEti + gXGL[ngl]*dInterv);
200  const G4double xs = isLPM ? ComputeRelDXSectionPerAtom(Et, gammaEnergy, Z)
201  : ComputeDXSectionPerAtom(Et, gammaEnergy, Z);
202  xSection += gWGL[ngl]*xs;
203  }
204  // update minimum Et of the sub-inteval
205  minEti += dInterv;
206  }
207  // apply corrections of variable transformation and half interval integration
208  xSection = std::max(2.*xSection*dInterv, 0.);
209  return xSection;
210 }
211 
212 // DCS WITHOUT LPM SUPPRESSION
213 // Computes DCS value for a given target element (Z), initial gamma energy (Eg),
214 // total energy transferred to one of the e-/e+ pair(Et) WITHOUT LPM suppression
215 // The constant factor 4 \alpha r_0^2 Z (Z +\eta(Z)) is not included here and
216 // the returned value will be differential in total energy transfer instead of
217 // the eps=Et/Eg. The computed part of the DCS
218 // NORMAL CASE: DEFAULT STTING (i.e. fIsUseCompleteScreening = FALSE)
219 // ds/deps(Et,Eg,Z) = ds/deps(eps,Z) = (eps^2+(1-eps)^2)*[phi1(d)/4-ln(Z)/3-fc]
220 // + 2*eps(1-eps)*[phi2(d)/4-ln(Z)/3-fc]/3 where the universal (in the TF model)
221 // screening variable d=d(eps)=136Z^(-1/3)eps0/[eps*(1-eps)] with eps0=mc^2/Eg.
222 // COMPLETE SCREENING (when d(eps) approx-equal-to 0) : NEED TO BE SET BY USER
223 // ds/deps(Et,Eg,Z) = ds/deps(eps,Z) = (eps^2+(1-eps)^2+eps*(1-eps)/3)*[Lel-fc]
224 // -eps(1-eps)/9 where Lel=phi1(0)/4-ln(Z)/3 is the elastic(coherent) radiation
225 // logarithm, fc is the Coulomb correction and the relation phi2(0)/4-ln(Z)/3 =
226 // phi1(0)/4-1/6-ln(Z)/3 = Lel-1/6 (due to phi2(0)=phi1(0)-2/3) was used.
228  G4double gammaEnergy,
229  G4double Z)
230 {
231  G4double xSection = 0.;
232  const G4int iz = std::min(gMaxZet, G4lrint(Z));
233  const G4double eps = pEnergy/gammaEnergy;
234  const G4double epsm = 1.-eps;
235  const G4double dum = eps*epsm;
237  // complete screening:
238  const G4double Lel = gElementData[iz]->fLradEl;
239  const G4double fc = gElementData[iz]->fCoulomb;
240  xSection = (eps*eps + epsm*epsm + 2.*dum/3.)*(Lel-fc) - dum/9.;
241  } else {
242  // normal case:
243  const G4double eps0 = CLHEP::electron_mass_c2/gammaEnergy;
244  const G4double fc = gElementData[iz]->fCoulomb;
245  const G4double lnZ13 = gElementData[iz]->fLogZ13;
246  const G4double delta = gElementData[iz]->fDeltaFactor*eps0/dum;
247  G4double phi1, phi2;
248  ComputePhi12(delta, phi1, phi2);
249  xSection = (eps*eps + epsm*epsm)*(0.25*phi1-lnZ13-fc)
250  + 2.*dum*(0.25*phi2-lnZ13-fc)/3.;
251  }
252  // non-const. part of the DCS differential in total energy transfer not in eps
253  // ds/dEt=ds/deps deps/dEt with deps/dEt=1/Eg
254  return std::max(xSection, 0.0)/gammaEnergy;
255 }
256 
257 // DCS WITH POSSIBLE LPM SUPPRESSION
258 // Computes DCS value for a given target element (Z), initial gamma energy (Eg),
259 // total energy transferred to one of the e-/e+ pair(Et) WITH LPM suppression.
260 // For a given Z, the LPM suppression will depend on the material through the
261 // LMP-Energy. This will determine the suppression variable s and the LPM sup-
262 // pression functions xi(s), fi(s) and G(s).
263 // The constant factor 4 \alpha r_0^2 Z (Z +\eta(Z)) is not included here and
264 // the returned value will be differential in total energy transfer instead of
265 // the eps=Et/Eg. The computed part of the DCS
266 // NORMAL CASE: DEFAULT STTING (i.e. fIsUseCompleteScreening = FALSE)
267 // ds/deps(Et,Eg,Z)=ds/deps(eps,Z) = xi(s)*{ (eps^2+(1-eps)^2)*[2fi(s)/3+G(s)/3]
268 // *[phi1(d)/4-ln(Z)/3-fc] + 2*eps(1-eps)*G(s)*[phi2(d)/4-ln(Z)/3-fc]/3 } where
269 // the universal (in the TF model) screening variable d=d(eps)=136Z^(-1/3)eps0
270 // /[eps*(1-eps)] with eps0=mc^2/Eg.
271 // COMPLETE SCREENING (when d(eps) approx-equal-to 0) : NEED TO BE SET BY USER
272 // ds/deps(Et,Eg,Z) = ds/deps(eps,Z) = xi(s)*{ [Lel-fc]*[ (eps^2+(1-eps)^2+eps
273 // *(1-eps)/3)*2fi(s)/3 + G(s)/3] - eps(1-eps)*G(s)/9 }
274 // Note, that when the LPM suppression is absent i.e. xi(s)=fi(s)=G(s)=1, both
275 // the normal and the complete screening DCS give back the NO-LMP case above.
277  G4double gammaEnergy,
278  G4double Z)
279 {
280  G4double xSection = 0.;
281  const G4int iz = std::min(gMaxZet, G4lrint(Z));
282  const G4double eps = pEnergy/gammaEnergy;
283  const G4double epsm = 1.-eps;
284  const G4double dum = eps*epsm;
285  // evaluate LPM suppression functions
286  G4double fXiS, fGS, fPhiS;
287  ComputeLPMfunctions(fXiS, fGS, fPhiS, eps, gammaEnergy, iz);
289  // complete screening:
290  const G4double Lel = gElementData[iz]->fLradEl;
291  const G4double fc = gElementData[iz]->fCoulomb;
292  xSection = (Lel-fc)*((eps*eps+epsm*epsm)*2.*fPhiS + fGS)/3. - dum*fGS/9.;
293  } else {
294  // normal case:
295  const G4double eps0 = CLHEP::electron_mass_c2/gammaEnergy;
296  const G4double fc = gElementData[iz]->fCoulomb;
297  const G4double lnZ13 = gElementData[iz]->fLogZ13;
298  const G4double delta = gElementData[iz]->fDeltaFactor*eps0/dum;
299  G4double phi1, phi2;
300  ComputePhi12(delta, phi1, phi2);
301  xSection = (eps*eps + epsm*epsm)*(2.*fPhiS+fGS)*(0.25*phi1-lnZ13-fc)/3.
302  + 2.*dum*fGS*(0.25*phi2-lnZ13-fc)/3.;
303  }
304  // non-const. part of the DCS differential in total energy transfer not in eps
305  // ds/dEt=ds/deps deps/dEt with deps/dEt=1/Eg
306  return std::max(fXiS*xSection, 0.0)/gammaEnergy;
307 }
308 
309 G4double
311  G4double gammaEnergy, G4double Z, G4double, G4double, G4double)
312 {
313  G4double crossSection = 0.0 ;
314  // check kinematical limit
315  if ( gammaEnergy <= 2.0*electron_mass_c2 ) { return crossSection; }
316  // compute the atomic cross section either by using x-section parametrization
317  // or by numerically integrationg the DCS (with or without LPM)
318  if ( gammaEnergy < fParametrizedXSectionThreshold) {
319  // using the parametrized cross sections (max up to 80 GeV)
320  crossSection = ComputeParametrizedXSectionPerAtom(gammaEnergy, Z);
321  } else {
322  // by numerical integration of the DCS:
323  // Computes the cross section with or without LPM suppression depending on
324  // settings (by default with if the gamma energy is above a given threshold)
325  // and using or not using complete sreening approximation (by default not).
326  // Only the dependent part is computed in the numerical integration of the DCS
327  // i.e. the result must be multiplied here with 4 \alpha r_0^2 Z(Z+\eta(Z))
328  crossSection = ComputeXSectionPerAtom(gammaEnergy, Z);
329  // apply the constant factors:
330  // - eta(Z) is a correction to account interaction in the field of e-
331  // - gXSecFactor = 4 \alpha r_0^2
332  const G4int iz = std::min(gMaxZet, G4lrint(Z));
333  const G4double eta = gElementData[iz]->fEtaValue;
334  crossSection *= gXSecFactor*Z*(Z+eta);
335  }
336  // final protection
337  return std::max(crossSection, 0.);
338 }
339 
341  const G4Material* mat, G4double)
342 {
344 }
345 
346 void
347 G4PairProductionRelModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
348  const G4MaterialCutsCouple* couple,
349  const G4DynamicParticle* aDynamicGamma,
350  G4double,
351  G4double)
352 // The secondaries e+e- energies are sampled using the Bethe - Heitler
353 // cross sections with Coulomb correction.
354 // A modified version of the random number techniques of Butcher & Messel
355 // is used (Nuc Phys 20(1960),15).
356 //
357 // GEANT4 internal units.
358 //
359 // Note 1 : Effects due to the breakdown of the Born approximation at
360 // low energy are ignored.
361 // Note 2 : The differential cross section implicitly takes account of
362 // pair creation in both nuclear and atomic electron fields.
363 // However triplet prodution is not generated.
364 {
365  const G4Material* mat = couple->GetMaterial();
366  const G4double gammaEnergy = aDynamicGamma->GetKineticEnergy();
367  const G4double eps0 = CLHEP::electron_mass_c2/gammaEnergy ;
368  //
369  // check kinematical limit: gamma energy(Eg) must be at least 2 e- rest mass
370  // (but the model should be used at higher energies above 100 MeV)
371  if (eps0 > 0.5) { return; }
372  //
373  // select target atom of the material
374  const G4Element* anElement = SelectTargetAtom(couple, fTheGamma, gammaEnergy,
375  aDynamicGamma->GetLogKineticEnergy());
376  CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
377  //
378  // 'eps' is the total energy transferred to one of the e-/e+ pair in initial
379  // gamma energy units Eg. Since the corresponding DCS is symmetric on eps=0.5,
380  // the kinematical limits for eps0=mc^2/Eg <= eps <= 0.5
381  // 1. 'eps' is sampled uniformly on the [eps0, 0.5] inteval if Eg<Egsmall
382  // 2. otherwise, on the [eps_min, 0.5] interval according to the DCS (case 2.)
383  G4double eps;
384  // case 1.
385  static const G4double Egsmall = 2.*CLHEP::MeV;
386  if (gammaEnergy < Egsmall) {
387  eps = eps0 + (0.5-eps0)*rndmEngine->flat();
388  } else {
389  // case 2.
390  // get the Coulomb factor for the target element (Z) and gamma energy (Eg)
391  // F(Z) = 8*ln(Z)/3 if Eg <= 50 [MeV] => no Coulomb correction
392  // F(Z) = 8*ln(Z)/3 + 8*fc(Z) if Eg > 50 [MeV] => fc(Z) is the Coulomb cor.
393  //
394  // The screening variable 'delta(eps)' = 136*Z^{-1/3}*eps0/[eps(1-eps)]
395  // Due to the Coulomb correction, the DCS can go below zero even at
396  // kinematicaly allowed eps > eps0 values. In order to exclude this eps
397  // range with negative DCS, the minimum eps value will be set to eps_min =
398  // max[eps0, epsp] with epsp is the solution of SF(delta(epsp)) - F(Z)/2 = 0
399  // with SF being the screening function (SF1=SF2 at high value of delta).
400  // The solution is epsp = 0.5 - 0.5*sqrt[ 1 - 4*136*Z^{-1/3}eps0/deltap]
401  // with deltap = Exp[(42.038-F(Z))/8.29]-0.958. So the limits are:
402  // - when eps=eps_max = 0.5 => delta_min = 136*Z^{-1/3}*eps0/4
403  // - epsp = 0.5 - 0.5*sqrt[ 1 - delta_min/deltap]
404  // - and eps_min = max[eps0, epsp]
405  const G4int iZet = std::min(gMaxZet, anElement->GetZasInt());
406  const G4double deltaFactor = gElementData[iZet]->fDeltaFactor*eps0;
407  const G4double deltaMin = 4.*deltaFactor;
408  G4double deltaMax = gElementData[iZet]->fDeltaMaxLow;
409  G4double FZ = 8.*gElementData[iZet]->fLogZ13;
410  if ( gammaEnergy > fCoulombCorrectionThreshold ) { // Eg > 50 MeV ?
411  FZ += 8.*gElementData[iZet]->fCoulomb;
412  deltaMax = gElementData[iZet]->fDeltaMaxHigh;
413  }
414  // compute the limits of eps
415  const G4double epsp = 0.5 - 0.5*std::sqrt(1. - deltaMin/deltaMax) ;
416  const G4double epsMin = std::max(eps0,epsp);
417  const G4double epsRange = 0.5 - epsMin;
418  //
419  // sample the energy rate (eps) of the created electron (or positron)
420  G4double F10, F20;
421  ScreenFunction12(deltaMin, F10, F20);
422  F10 -= FZ;
423  F20 -= FZ;
424  const G4double NormF1 = std::max(F10 * epsRange * epsRange, 0.);
425  const G4double NormF2 = std::max(1.5 * F20 , 0.);
426  const G4double NormCond = NormF1/(NormF1 + NormF2);
427  // check if LPM correction is active
428  const G4bool isLPM = (fIsUseLPMCorrection && gammaEnergy>gEgLPMActivation);
430  // we will need 3 uniform random number for each trial of sampling
431  G4double rndmv[3];
432  G4double greject = 0.;
433  do {
434  rndmEngine->flatArray(3, rndmv);
435  if (NormCond > rndmv[0]) {
436  eps = 0.5 - epsRange * fG4Calc->A13(rndmv[1]);
437  const G4double delta = deltaFactor/(eps*(1.-eps));
438  if (isLPM) {
439  G4double lpmXiS, lpmGS, lpmPhiS, phi1, phi2;
440  ComputePhi12(delta, phi1, phi2);
441  ComputeLPMfunctions(lpmXiS, lpmGS, lpmPhiS, eps, gammaEnergy, iZet);
442  greject = lpmXiS*((2.*lpmPhiS+lpmGS)*phi1-lpmGS*phi2-lpmPhiS*FZ)/F10;
443  } else {
444  greject = (ScreenFunction1(delta)-FZ)/F10;
445  }
446  } else {
447  eps = epsMin + epsRange*rndmv[1];
448  const G4double delta = deltaFactor/(eps*(1.-eps));
449  if (isLPM) {
450  G4double lpmXiS, lpmGS, lpmPhiS, phi1, phi2;
451  ComputePhi12(delta, phi1, phi2);
452  ComputeLPMfunctions(lpmXiS, lpmGS, lpmPhiS, eps, gammaEnergy, iZet);
453  greject = lpmXiS*( (lpmPhiS+0.5*lpmGS)*phi1 + 0.5*lpmGS*phi2
454  -0.5*(lpmGS+lpmPhiS)*FZ )/F20;
455  } else {
456  greject = (ScreenFunction2(delta)-FZ)/F20;
457  }
458  }
459  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
460  } while (greject < rndmv[2]);
461  // end of eps sampling
462  }
463  //
464  // select charges randomly
465  G4double eTotEnergy, pTotEnergy;
466  if (rndmEngine->flat() > 0.5) {
467  eTotEnergy = (1.-eps)*gammaEnergy;
468  pTotEnergy = eps*gammaEnergy;
469  } else {
470  pTotEnergy = (1.-eps)*gammaEnergy;
471  eTotEnergy = eps*gammaEnergy;
472  }
473  //
474  // sample pair kinematics
475  //
476  const G4double eKinEnergy = std::max(0.,eTotEnergy - CLHEP::electron_mass_c2);
477  const G4double pKinEnergy = std::max(0.,pTotEnergy - CLHEP::electron_mass_c2);
478  //
479  G4ThreeVector eDirection, pDirection;
480  //
482  eKinEnergy, pKinEnergy, eDirection, pDirection);
483  // create G4DynamicParticle object for the particle1
484  G4DynamicParticle* aParticle1= new G4DynamicParticle(
485  fTheElectron,eDirection,eKinEnergy);
486 
487  // create G4DynamicParticle object for the particle2
488  G4DynamicParticle* aParticle2= new G4DynamicParticle(
489  fThePositron,pDirection,pKinEnergy);
490  // Fill output vector
491  fvect->push_back(aParticle1);
492  fvect->push_back(aParticle2);
493  // kill incident photon
496 }
497 
498 // should be called only by the master and at initialisation
500 {
501  G4int size = gElementData.size();
502  if (size < gMaxZet+1) {
503  gElementData.resize(gMaxZet+1, nullptr);
504  }
505  // create for all elements that are in the detector
506  const G4ElementTable* elemTable = G4Element::GetElementTable();
507  size_t numElems = (*elemTable).size();
508  for (size_t ie = 0; ie < numElems; ++ie) {
509  const G4Element* elem = (*elemTable)[ie];
510  const G4int iz = std::min(gMaxZet, elem->GetZasInt());
511  if (!gElementData[iz]) { // create it if doesn't exist yet
512  const G4double logZ13 = elem->GetIonisation()->GetlogZ3();
513  const G4double Z13 = elem->GetIonisation()->GetZ3();
514  const G4double fc = elem->GetfCoulomb();
515  const G4double FZLow = 8.*logZ13;
516  const G4double FZHigh = 8.*(logZ13 + fc);
517  G4double Fel;
518  G4double Finel;
519  if (iz<5) { // use data from Dirac-Fock atomic model
520  Fel = gFelLowZet[iz];
521  Finel = gFinelLowZet[iz];
522  } else { // use the results of the Thomas-Fermi-Moliere model
523  Fel = G4Log(184.) - logZ13;
524  Finel = G4Log(1194.) - 2.*logZ13;
525  }
526  ElementData* elD = new ElementData();
527  elD->fLogZ13 = logZ13;
528  elD->fCoulomb = fc;
529  elD->fLradEl = Fel;
530  elD->fDeltaFactor = 136./Z13;
531  elD->fDeltaMaxLow = G4Exp((42.038 - FZLow)/8.29) - 0.958;
532  elD->fDeltaMaxHigh = G4Exp((42.038 - FZHigh)/8.29) - 0.958;
533  elD->fEtaValue = Finel/(Fel-fc);
534  elD->fLPMVarS1Cond = std::sqrt(2.)*Z13*Z13/(184.*184.);
535  elD->fLPMILVarS1Cond = 1./G4Log(elD->fLPMVarS1Cond);
536  gElementData[iz] = elD;
537  }
538  }
539 }
540 
541 // s goes up to 2 with ds = 0.01 be default
543  if (!gLPMFuncs.fIsInitialized) {
545  gLPMFuncs.fLPMFuncG.resize(num);
546  gLPMFuncs.fLPMFuncPhi.resize(num);
547  for (G4int i=0; i<num; ++i) {
548  const G4double sval = i/gLPMFuncs.fISDelta;
550  }
551  gLPMFuncs.fIsInitialized = true;
552  }
553 }
554 
555 // used only at initialisation time
556 void G4PairProductionRelModel::ComputeLPMGsPhis(G4double &funcGS, G4double &funcPhiS, const G4double varShat) {
557  if (varShat < 0.01) {
558  funcPhiS = 6.0*varShat*(1.0-CLHEP::pi*varShat);
559  funcGS = 12.0*varShat-2.0*funcPhiS;
560  } else {
561  const G4double varShat2 = varShat*varShat;
562  const G4double varShat3 = varShat*varShat2;
563  const G4double varShat4 = varShat2*varShat2;
564  if (varShat < 0.415827397755) { // Stanev ap.: for \psi(s) and compute G(s)
565  funcPhiS = 1.0-G4Exp( -6.0*varShat*(1.0+varShat*(3.0-CLHEP::pi))
566  + varShat3/(0.623+0.796*varShat+0.658*varShat2));
567  // 1-\exp \left\{-4s-\frac{8s^2}{1+3.936s+4.97s^2-0.05s^3+7.5s^4} \right\}
568  const G4double funcPsiS = 1.0-G4Exp( -4.0*varShat - 8.0*varShat2/(1.0
569  + 3.936*varShat+4.97*varShat2-0.05*varShat3+7.5*varShat4));
570  // G(s) = 3 \psi(s) - 2 \phi(s)
571  funcGS = 3.0*funcPsiS - 2.0*funcPhiS;
572  } else if (varShat < 1.55) {
573  funcPhiS = 1.0-G4Exp( -6.0*varShat*(1.0+varShat*(3.0-CLHEP::pi))
574  + varShat3/(0.623+0.796*varShat+0.658*varShat2));
575  const G4double dum0 = -0.16072300849123999+3.7550300067531581*varShat
576  -1.7981383069010097 *varShat2
577  +0.67282686077812381*varShat3
578  -0.1207722909879257 *varShat4;
579  funcGS = std::tanh(dum0);
580  } else {
581  funcPhiS = 1.0-0.01190476/varShat4;
582  if (varShat < 1.9156) {
583  const G4double dum0 = -0.16072300849123999+3.7550300067531581*varShat
584  -1.7981383069010097 *varShat2
585  +0.67282686077812381*varShat3
586  -0.1207722909879257 *varShat4;
587  funcGS = std::tanh(dum0);
588  } else {
589  funcGS = 1.0-0.0230655/varShat4;
590  }
591  }
592  }
593 }
594 
595 // used at run-time to get some pre-computed LPM function values
597  G4double &lpmPhis,
598  const G4double sval) {
599  if (sval < gLPMFuncs.fSLimit) {
600  G4double val = sval*gLPMFuncs.fISDelta;
601  const G4int ilow = (G4int)val;
602  val -= ilow;
603  lpmGs = (gLPMFuncs.fLPMFuncG[ilow+1]-gLPMFuncs.fLPMFuncG[ilow])*val
604  + gLPMFuncs.fLPMFuncG[ilow];
605  lpmPhis = (gLPMFuncs.fLPMFuncPhi[ilow+1]-gLPMFuncs.fLPMFuncPhi[ilow])*val
606  + gLPMFuncs.fLPMFuncPhi[ilow];
607  } else {
608  G4double ss = sval*sval;
609  ss *= ss;
610  lpmPhis = 1.0-0.01190476/ss;
611  lpmGs = 1.0-0.0230655/ss;
612  }
613 }
614 
616  G4double &funcGS, G4double &funcPhiS, const G4double eps,
617  const G4double egamma, const G4int izet)
618 {
619  // 1. y = E_+/E_{\gamma} with E_+ being the total energy transfered
620  // to one of the e-/e+ pair
621  // s' = \sqrt{ \frac{1}{8} \frac{1}{y(1-y)} \frac{E^{KL}_{LPM}}{E_{\gamma}} }
622  const G4double varSprime = std::sqrt(0.125*fLPMEnergy/(eps*egamma*(1.0-eps)));
623  const G4double condition = gElementData[izet]->fLPMVarS1Cond;
624  funcXiS = 2.0;
625  if (varSprime > 1.0) {
626  funcXiS = 1.0;
627  } else if (varSprime > condition) {
628  const G4double dum = gElementData[izet]->fLPMILVarS1Cond;
629  const G4double funcHSprime = G4Log(varSprime)*dum;
630  funcXiS = 1.0 + funcHSprime
631  - 0.08*(1.0-funcHSprime)*funcHSprime*(2.0-funcHSprime)*dum;
632  }
633  // 2. s=\frac{s'}{\sqrt{\xi(s')}}
634  const G4double varShat = varSprime / std::sqrt(funcXiS);
635  GetLPMFunctions(funcGS, funcPhiS, varShat);
636  // MAKE SURE SUPPRESSION IS SMALLER THAN 1: due to Migdal's approximation on xi
637  if (funcXiS * funcPhiS > 1. || varShat > 0.57) {
638  funcXiS = 1. / funcPhiS;
639  }
640 }
641 
642 // Calculates the microscopic cross section in GEANT4 internal units. Same as in
643 // G4BetheHeitlerModel and should be used below 80 GeV since it start to deverge
644 // from the cross section data above 80-90 GeV:
645 // Parametrized formula (L. Urban) is used to estimate the atomic cross sections
646 // given numerically in the table of [Hubbell, J. H., Heinz Albert Gimm, and I.
647 // Overbo: "Pair, Triplet, and Total Atomic Cross Sections (and Mass Attenuation
648 // Coefficients) for 1 MeV‐100 GeV Photons in Elements Z= 1 to 100." Journal of
649 // physical and chemical reference data 9.4 (1980): 1023-1148.]
650 //
651 // The formula gives a good approximation of the data from 1.5 MeV to 100 GeV.
652 // below 1.5 MeV: sigma=sigma(1.5MeV)*(GammaEnergy-2electronmass)
653 // *(GammaEnergy-2electronmass)
654 G4double
656  G4double Z)
657 {
658  G4double xSection = 0.0 ;
659  // short versions
660  static const G4double kMC2 = CLHEP::electron_mass_c2;
661  // zero cross section below the kinematical limit: Eg<2mc^2
662  if (Z < 0.9 || gammaE <= 2.0*kMC2) { return xSection; }
663  //
664  static const G4double gammaEnergyLimit = 1.5*CLHEP::MeV;
665  // set coefficients a, b c
666  static const G4double a0 = 8.7842e+2*CLHEP::microbarn;
667  static const G4double a1 = -1.9625e+3*CLHEP::microbarn;
668  static const G4double a2 = 1.2949e+3*CLHEP::microbarn;
669  static const G4double a3 = -2.0028e+2*CLHEP::microbarn;
670  static const G4double a4 = 1.2575e+1*CLHEP::microbarn;
671  static const G4double a5 = -2.8333e-1*CLHEP::microbarn;
672 
673  static const G4double b0 = -1.0342e+1*CLHEP::microbarn;
674  static const G4double b1 = 1.7692e+1*CLHEP::microbarn;
675  static const G4double b2 = -8.2381 *CLHEP::microbarn;
676  static const G4double b3 = 1.3063 *CLHEP::microbarn;
677  static const G4double b4 = -9.0815e-2*CLHEP::microbarn;
678  static const G4double b5 = 2.3586e-3*CLHEP::microbarn;
679 
680  static const G4double c0 = -4.5263e+2*CLHEP::microbarn;
681  static const G4double c1 = 1.1161e+3*CLHEP::microbarn;
682  static const G4double c2 = -8.6749e+2*CLHEP::microbarn;
683  static const G4double c3 = 2.1773e+2*CLHEP::microbarn;
684  static const G4double c4 = -2.0467e+1*CLHEP::microbarn;
685  static const G4double c5 = 6.5372e-1*CLHEP::microbarn;
686  // check low energy limit of the approximation (1.5 MeV)
687  G4double gammaEnergyOrg = gammaE;
688  if (gammaE < gammaEnergyLimit) { gammaE = gammaEnergyLimit; }
689  // compute gamma energy variables
690  const G4double x = G4Log(gammaE/kMC2);
691  const G4double x2 = x *x;
692  const G4double x3 = x2*x;
693  const G4double x4 = x3*x;
694  const G4double x5 = x4*x;
695  //
696  const G4double F1 = a0 + a1*x + a2*x2 + a3*x3 + a4*x4 + a5*x5;
697  const G4double F2 = b0 + b1*x + b2*x2 + b3*x3 + b4*x4 + b5*x5;
698  const G4double F3 = c0 + c1*x + c2*x2 + c3*x3 + c4*x4 + c5*x5;
699  // compute the approximated cross section
700  xSection = (Z + 1.)*(F1*Z + F2*Z*Z + F3);
701  // check if we are below the limit of the approximation and apply correction
702  if (gammaEnergyOrg < gammaEnergyLimit) {
703  const G4double dum = (gammaEnergyOrg-2.*kMC2)/(gammaEnergyLimit-2.*kMC2);
704  xSection *= dum*dum;
705  }
706  return xSection;
707 }
708 
709