ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PairProductionRelModel.hh
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PairProductionRelModel.hh
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 header file
30 //
31 //
32 // File name: G4PairProductionRelModel
33 //
34 // Author: Andreas Schaelicke
35 //
36 // Creation date: 02.04.2009
37 //
38 // Modifications:
39 // 28-05-18 New version with improved screening function approximation, improved
40 // LPM function approximation, efficiency, documentation and cleanup.
41 // Corrected call to selecting target atom in the final state sampling.
42 // (M. Novak)
43 //
44 // Class Description:
45 //
46 // Implementation of gamma convertion to e+e- in the field of a nucleus
47 // relativistic approximation
48 //
49 
50 // -------------------------------------------------------------------
51 //
52 
53 #ifndef G4PairProductionRelModel_h
54 #define G4PairProductionRelModel_h 1
55 
56 #include <CLHEP/Units/PhysicalConstants.h>
57 
58 #include "G4VEmModel.hh"
59 #include "G4Log.hh"
60 #include "G4Exp.hh"
61 #include "G4Pow.hh"
62 
63 #include <vector>
64 
65 class G4ParticleChangeForGamma;
66 
67 class G4PairProductionRelModel : public G4VEmModel
68 {
69 
70 public:
71 
72  explicit G4PairProductionRelModel(const G4ParticleDefinition* p = nullptr,
73  const G4String& nam = "BetheHeitlerLPM");
74 
75  virtual ~G4PairProductionRelModel();
76 
77  virtual void Initialise(const G4ParticleDefinition*,
78  const G4DataVector&) override;
79 
80  virtual void InitialiseLocal(const G4ParticleDefinition*,
81  G4VEmModel* masterModel) override;
82 
83  virtual G4double ComputeCrossSectionPerAtom(
84  const G4ParticleDefinition*,
85  G4double kinEnergy,
86  G4double Z,
87  G4double A=0.,
88  G4double cut=0.,
89  G4double emax=DBL_MAX) override;
90 
91  virtual void SampleSecondaries(std::vector<G4DynamicParticle*>*,
92  const G4MaterialCutsCouple*,
93  const G4DynamicParticle*,
94  G4double tmin,
95  G4double maxEnergy) override;
96 
97  virtual void SetupForMaterial(const G4ParticleDefinition*,
98  const G4Material*,G4double) override;
99 
100  inline void SetLPMflag(G4bool val) { fIsUseLPMCorrection = val; }
101  inline G4bool LPMflag() const { return fIsUseLPMCorrection; }
102 
103 protected:
104 
105  // for evaluating screening related functions
106  inline void ComputePhi12(const G4double delta, G4double &phi1, G4double &phi2);
107  inline G4double ScreenFunction1(const G4double delta);
108  inline G4double ScreenFunction2(const G4double delta);
109  inline void ScreenFunction12(const G4double delta, G4double &f1, G4double &f2);
110  // helper methods for cross-section computation under different approximations
111  G4double ComputeParametrizedXSectionPerAtom(G4double gammaEnergy, G4double Z);
112  G4double ComputeXSectionPerAtom(G4double gammaEnergy, G4double Z);
113  G4double ComputeDXSectionPerAtom(G4double eplusEnergy, G4double gammaEnergy,
114  G4double Z);
115  G4double ComputeRelDXSectionPerAtom(G4double eplusEnergy, G4double gammaEnergy,
116  G4double Z);
117 
118 
119 private:
120 
121  // for creating some data structure per Z with often used comp. intensive data
122  void InitialiseElementData();
123  struct ElementData {
124  G4double fLogZ13;
125  G4double fCoulomb;
126  G4double fLradEl;
127  G4double fDeltaFactor;
128  G4double fDeltaMaxLow;
129  G4double fDeltaMaxHigh;
130  G4double fEtaValue;
131  G4double fLPMVarS1Cond;
132  G4double fLPMILVarS1Cond;
133  };
134  // for precomputing comp. intensive parts of LPM suppression functions and
135  // using them at run-time
136  void InitLPMFunctions();
137  void ComputeLPMGsPhis(G4double &funcGS, G4double &funcPhiS,
138  const G4double varShat);
139  void GetLPMFunctions(G4double &lpmGs, G4double &lpmPhis, const G4double sval);
140  void ComputeLPMfunctions(G4double &fXiS, G4double &fGS, G4double &fPhiS,
141  const G4double eps, const G4double egamma,
142  const G4int izet);
143  struct LPMFuncs {
144  LPMFuncs() : fIsInitialized(false), fISDelta(100.), fSLimit(2.) {}
145  G4bool fIsInitialized;
146  G4double fISDelta;
147  G4double fSLimit;
148  std::vector<G4double> fLPMFuncG;
149  std::vector<G4double> fLPMFuncPhi;
150  };
151 
152 
153 private:
154 
155  // hide assignment operator
156  G4PairProductionRelModel & operator=
157  (const G4PairProductionRelModel &right) = delete;
158  G4PairProductionRelModel(const G4PairProductionRelModel&) = delete;
159 
160 protected:
161  static const G4int gMaxZet;
162  //
163  static const G4double gLPMconstant;
164  //
165  static const G4double gXGL[8];
166  static const G4double gWGL[8];
167  static const G4double gFelLowZet[8];
168  static const G4double gFinelLowZet[8];
169  //
170  static const G4double gXSecFactor;
171  static const G4double gEgLPMActivation;
172  //
173  static std::vector<ElementData*> gElementData;
174  static LPMFuncs gLPMFuncs;
175  //
176  G4bool fIsUseLPMCorrection;
177  G4bool fIsUseCompleteScreening;
178  //
179  G4double fLPMEnergy;
180  //
181  G4double fParametrizedXSectionThreshold;
182  G4double fCoulombCorrectionThreshold;
183  //
184  G4Pow* fG4Calc;
185  G4ParticleDefinition* fTheGamma;
186  G4ParticleDefinition* fTheElectron;
187  G4ParticleDefinition* fThePositron;
188  G4ParticleChangeForGamma* fParticleChange;
189 };
190 
191 
192 
193 //
194 // Bethe screening functions for the elastic (coherent) scattering:
195 // Bethe's phi1, phi2 coherent screening functions were computed numerically
196 // by using (the universal) atomic form factors computed based on the Thomas-
197 // Fermi model of the atom (using numerical solution of the Thomas-Fermi
198 // screening function instead of Moliere's analytical approximation). The
199 // numerical results can be well approximated (better than Butcher & Messel
200 // especially near the delta=1 limit) by:
201 // ## if delta <= 1.4
202 // phi1(delta) = 20.806 - delta*(3.190 - 0.5710*delta)
203 // phi2(delta) = 20.234 - delta*(2.126 - 0.0903*delta)
204 // ## if delta > 1.4
205 // phi1(delta) = phi2(delta) = 21.0190 - 4.145*ln(delta + 0.958)
206 // with delta = 136mc^2kZ^{-1/3}/[E(Eg-E)] = 136Z^{-1/3}eps0/[eps(1-eps)] where
207 // Eg is the initial photon energy, E is the total energy transferred to one of
208 // the e-/e+ pair, eps0 = mc^2/Eg and eps = E/Eg.
209 
210 inline
211 void G4PairProductionRelModel::ComputePhi12(const G4double delta, G4double &phi1,
212  G4double &phi2)
213 {
214  if (delta > 1.4) {
215  phi1 = 21.0190 - 4.145*G4Log(delta + 0.958);
216  phi2 = phi1;
217  } else {
218  phi1 = 20.806 - delta*(3.190 - 0.5710*delta);
219  phi2 = 20.234 - delta*(2.126 - 0.0903*delta);
220  }
221 }
222 
223 
224 // Compute the value of the screening function 3*PHI1(delta) - PHI2(delta):
225 inline G4double G4PairProductionRelModel::ScreenFunction1(const G4double delta)
226 {
227  return (delta > 1.4) ? 42.038 - 8.29*G4Log(delta + 0.958)
228  : 42.184 - delta*(7.444 - 1.623*delta);
229 }
230 
231 
232 // Compute the value of the screening function 1.5*PHI1(delta) +0.5*PHI2(delta):
233 inline G4double G4PairProductionRelModel::ScreenFunction2(const G4double delta)
234 {
235  return (delta > 1.4) ? 42.038 - 8.29*G4Log(delta + 0.958)
236  : 41.326 - delta*(5.848 - 0.902*delta);
237 }
238 
239 
240 // Same as ScreenFunction1 and ScreenFunction2 but computes them at once
241 inline void G4PairProductionRelModel::ScreenFunction12(const G4double delta,
242  G4double &f1, G4double &f2)
243 {
244  if (delta > 1.4) {
245  f1 = 42.038 - 8.29*G4Log(delta + 0.958);
246  f2 = f1;
247  } else {
248  f1 = 42.184 - delta*(7.444 - 1.623*delta);
249  f2 = 41.326 - delta*(5.848 - 0.902*delta);
250  }
251 }
252 
253 #endif