| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 // $Id: G4PairProductionRelModel.hh 106628 2017-10-17 06:25:38Z gcosmo $
26 // 27 //
27 // ------------------------------------------- 28 // -------------------------------------------------------------------
28 // 29 //
29 // GEANT4 Class header file 30 // GEANT4 Class header file
30 // 31 //
31 // 32 //
32 // File name: G4PairProductionRelModel 33 // File name: G4PairProductionRelModel
33 // 34 //
34 // Author: Andreas Schaelicke 35 // Author: Andreas Schaelicke
35 // 36 //
36 // Creation date: 02.04.2009 37 // Creation date: 02.04.2009
37 // 38 //
38 // Modifications: 39 // Modifications:
39 // 28-05-18 New version with improved screenin <<
40 // LPM function approximation, effici <<
41 // Corrected call to selecting target <<
42 // (M. Novak) <<
43 // 40 //
44 // Class Description: 41 // Class Description:
45 // 42 //
46 // Implementation of gamma conversion to e+e- << 43 // Implementation of gamma convertion to e+e- in the field of a nucleus
47 // relativistic approximation 44 // relativistic approximation
48 // 45 //
49 46
50 // ------------------------------------------- 47 // -------------------------------------------------------------------
51 // 48 //
52 49
53 #ifndef G4PairProductionRelModel_h 50 #ifndef G4PairProductionRelModel_h
54 #define G4PairProductionRelModel_h 1 51 #define G4PairProductionRelModel_h 1
55 52
56 #include <CLHEP/Units/PhysicalConstants.h> 53 #include <CLHEP/Units/PhysicalConstants.h>
57 54
58 #include "G4VEmModel.hh" 55 #include "G4VEmModel.hh"
>> 56 #include "G4PhysicsTable.hh"
59 #include "G4Log.hh" 57 #include "G4Log.hh"
60 #include <vector> << 58 #include "G4Exp.hh"
>> 59 #include "G4Pow.hh"
61 60
62 class G4ParticleChangeForGamma; 61 class G4ParticleChangeForGamma;
63 class G4Pow; <<
64 62
65 class G4PairProductionRelModel : public G4VEmM 63 class G4PairProductionRelModel : public G4VEmModel
66 { 64 {
67 65
68 public: 66 public:
69 67
70 explicit G4PairProductionRelModel(const G4Pa 68 explicit G4PairProductionRelModel(const G4ParticleDefinition* p = nullptr,
71 const G4St << 69 const G4String& nam = "BetheHeitlerLPM");
72 70
73 ~G4PairProductionRelModel() override; << 71 virtual ~G4PairProductionRelModel();
74 72
75 void Initialise(const G4ParticleDefinition*, << 73 virtual void Initialise(const G4ParticleDefinition*,
>> 74 const G4DataVector&) override;
76 75
77 void InitialiseLocal(const G4ParticleDefinit << 76 virtual void InitialiseLocal(const G4ParticleDefinition*,
78 G4VEmModel* masterModel) override; << 77 G4VEmModel* masterModel) override;
79 78
80 G4double ComputeCrossSectionPerAtom(const G4 << 79 virtual G4double ComputeCrossSectionPerAtom(
81 G4double kinEnergy, << 80 const G4ParticleDefinition*,
82 G4double Z, << 81 G4double kinEnergy,
83 G4double A=0., << 82 G4double Z,
84 G4double cut=0., << 83 G4double A=0.,
85 G4double emax=DBL_MAX) override; << 84 G4double cut=0.,
86 << 85 G4double emax=DBL_MAX) override;
87 void SampleSecondaries(std::vector<G4Dynamic << 86
88 const G4MaterialCutsCouple*, << 87 virtual void SampleSecondaries(std::vector<G4DynamicParticle*>*,
89 const G4DynamicParticle*, << 88 const G4MaterialCutsCouple*,
90 G4double tmin, << 89 const G4DynamicParticle*,
91 G4double maxEnergy) override; << 90 G4double tmin,
>> 91 G4double maxEnergy) override;
>> 92
>> 93 virtual void SetupForMaterial(const G4ParticleDefinition*,
>> 94 const G4Material*,G4double) override;
>> 95
>> 96 // * fast inline functions *
>> 97 inline void SetCurrentElement(G4double Z);
>> 98
>> 99 // set / get methods
>> 100 inline void SetLPMconstant(G4double val);
>> 101 inline G4double LPMconstant() const;
92 102
93 void SetupForMaterial(const G4ParticleDefini << 103 inline void SetLPMflag(G4bool);
94 const G4Material*,G4double) override; << 104 inline G4bool LPMflag() const;
95 105
96 inline void SetLPMflag(G4bool val) { fIsUs << 106 protected:
97 inline G4bool LPMflag() const { retur << 107
>> 108 // screening functions
>> 109 inline G4double Phi1(G4double delta) const;
>> 110 inline G4double Phi2(G4double delta) const;
>> 111 inline G4double ScreenFunction1(G4double ScreenVariable);
>> 112 inline G4double ScreenFunction2(G4double ScreenVariable);
>> 113 inline G4double DeltaMax() const;
>> 114 inline G4double DeltaMin(G4double) const;
>> 115
>> 116 // lpm functions
>> 117 void CalcLPMFunctions(G4double k, G4double eplus);
>> 118
>> 119 G4double ComputeXSectionPerAtom(G4double totalEnergy, G4double Z);
98 120
>> 121 G4double ComputeDXSectionPerAtom(G4double eplusEnergy, G4double totalEnergy, G4double Z);
>> 122 G4double ComputeRelDXSectionPerAtom(G4double eplusEnergy, G4double totalEnergy, G4double Z);
>> 123
>> 124 // hide assignment operator
99 G4PairProductionRelModel & operator= 125 G4PairProductionRelModel & operator=
100 (const G4PairProductionRelModel &right) = de 126 (const G4PairProductionRelModel &right) = delete;
101 G4PairProductionRelModel(const G4PairProduc 127 G4PairProductionRelModel(const G4PairProductionRelModel&) = delete;
102 128
103 protected: << 129 G4Pow* g4calc;
>> 130 G4ParticleDefinition* theGamma;
>> 131 G4ParticleDefinition* theElectron;
>> 132 G4ParticleDefinition* thePositron;
>> 133 G4ParticleChangeForGamma* fParticleChange;
>> 134
>> 135 G4double fLPMconstant;
>> 136 G4bool fLPMflag;
>> 137
>> 138 // cash
>> 139 G4double z13, z23, lnZ;
>> 140 G4double Fel, Finel, fCoulomb;
>> 141 G4double currentZ;
>> 142
>> 143 // LPM effect
>> 144 G4double lpmEnergy;
>> 145 G4double xiLPM, phiLPM, gLPM;
>> 146
>> 147 // consts
>> 148 G4bool use_completescreening;
>> 149
>> 150 static const G4double xgi[8], wgi[8];
>> 151 static const G4double Fel_light[5];
>> 152 static const G4double Finel_light[5];
>> 153 static const G4double facFel;
>> 154 static const G4double facFinel;
104 155
105 // for evaluating screening related function << 156 static const G4double preS1, logTwo, xsfactor, Egsmall, Eghigh;
106 inline void ComputePhi12(const G4double delt <<
107 G4double &phi1, G4double &phi2); <<
108 inline G4double ScreenFunction1(const G4doub <<
109 inline G4double ScreenFunction2(const G4doub <<
110 inline void ScreenFunction12(const G4double <<
111 G4double &f1, G4double &f2); <<
112 // helper methods for cross-section computat <<
113 G4double ComputeParametrizedXSectionPerAtom( <<
114 G4double ComputeXSectionPerAtom(G4double gam <<
115 G4double ComputeDXSectionPerAtom(G4double ep <<
116 G4double Z) <<
117 G4double ComputeRelDXSectionPerAtom(G4double <<
118 G4double gammaEnergy, G4double Z <<
119 <<
120 private: <<
121 <<
122 // for creating some data structure per Z <<
123 void InitialiseElementData(); <<
124 struct ElementData { <<
125 G4double fLogZ13; <<
126 G4double fCoulomb; <<
127 G4double fLradEl; <<
128 G4double fDeltaFactor; <<
129 G4double fDeltaMaxLow; <<
130 G4double fDeltaMaxHigh; <<
131 G4double fEtaValue; <<
132 G4double fLPMVarS1Cond; <<
133 G4double fLPMILVarS1Cond; <<
134 }; <<
135 // for precomputing comp. intensive parts of <<
136 // using them at run-time <<
137 void InitLPMFunctions(); <<
138 void ComputeLPMGsPhis(G4double &funcGS, G4do <<
139 const G4double varShat <<
140 void GetLPMFunctions(G4double &lpmGs, G4doub <<
141 void ComputeLPMfunctions(G4double &fXiS, G4d <<
142 const G4double eps, <<
143 const G4int izet); <<
144 struct LPMFuncs { <<
145 LPMFuncs() : fIsInitialized(false), fISDel <<
146 G4bool fIsInitialized; <<
147 G4double fISDelta; <<
148 G4double fSLimit; <<
149 std::vector<G4double> fLPMFuncG; <<
150 std::vector<G4double> fLPMFuncPhi; <<
151 }; <<
152 157
153 protected: <<
154 static const G4int gMaxZet; <<
155 // <<
156 static const G4double gLPMconsta <<
157 // <<
158 static const G4double gXGL[8]; <<
159 static const G4double gWGL[8]; <<
160 static const G4double gFelLowZet <<
161 static const G4double gFinelLowZ <<
162 // <<
163 static const G4double gXSecFacto <<
164 static const G4double gEgLPMActi <<
165 // <<
166 static std::vector<ElementData*> gElementDa <<
167 static LPMFuncs gLPMFuncs; <<
168 // <<
169 G4bool isFirstInstance{false}; <<
170 G4bool fIsUseLPMC <<
171 G4bool fIsUseComp <<
172 // <<
173 G4double fLPMEnergy <<
174 // <<
175 G4double fParametri <<
176 G4double fCoulombCo <<
177 // <<
178 G4Pow* fG4Calc; <<
179 G4ParticleDefinition* fTheGamma; <<
180 G4ParticleDefinition* fTheElectr <<
181 G4ParticleDefinition* fThePositr <<
182 G4ParticleChangeForGamma* fParticleC <<
183 }; 158 };
184 // << 159
185 // Bethe screening functions for the elastic ( << 160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
186 // Bethe's phi1, phi2 coherent screening funct << 161
187 // by using (the universal) atomic form factor << 162 inline
188 // Fermi model of the atom (using numerical so << 163 void G4PairProductionRelModel::SetLPMconstant(G4double val)
189 // screening function instead of Moliere's ana << 164 {
190 // numerical results can be well approximated << 165 fLPMconstant = val;
191 // especially near the delta=1 limit) by: <<
192 // ## if delta <= 1.4 <<
193 // phi1(delta) = 20.806 - delta*(3.190 - 0.57 <<
194 // phi2(delta) = 20.234 - delta*(2.126 - 0.09 <<
195 // ## if delta > 1.4 <<
196 // phi1(delta) = phi2(delta) = 21.0190 - 4.14 <<
197 // with delta = 136mc^2kZ^{-1/3}/[E(Eg-E)] = 1 <<
198 // Eg is the initial photon energy, E is the t <<
199 // the e-/e+ pair, eps0 = mc^2/Eg and eps = E/ <<
200 <<
201 inline void G4PairProductionRelModel::ComputeP <<
202 G4double &phi1, <<
203 G4double &phi2) <<
204 { <<
205 if (delta > 1.4) { <<
206 phi1 = 21.0190 - 4.145*G4Log(delta + 0.9 <<
207 phi2 = phi1; <<
208 } else { <<
209 phi1 = 20.806 - delta*(3.190 - 0.5710*de <<
210 phi2 = 20.234 - delta*(2.126 - 0.0903*de <<
211 } <<
212 } 166 }
213 167
214 // Compute the value of the screening function << 168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
215 inline G4double G4PairProductionRelModel::Scre << 169
>> 170 inline
>> 171 G4double G4PairProductionRelModel::LPMconstant() const
216 { 172 {
217 return (delta > 1.4) ? 42.038 - 8.29*G4Log(d << 173 return fLPMconstant;
218 : 42.184 - delta*(7.444 <<
219 } 174 }
220 175
221 // Compute the value of the screening function << 176 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
222 inline G4double G4PairProductionRelModel::Scre << 177
>> 178 inline
>> 179 void G4PairProductionRelModel::SetLPMflag(G4bool val)
223 { 180 {
224 return (delta > 1.4) ? 42.038 - 8.29*G4Log(d << 181 fLPMflag = val;
225 : 41.326 - delta*(5.848 <<
226 } 182 }
227 183
228 // Same as ScreenFunction1 and ScreenFunction2 << 184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
229 inline void G4PairProductionRelModel::ScreenFu << 185
230 << 186 inline
>> 187 G4bool G4PairProductionRelModel::LPMflag() const
231 { 188 {
232 if (delta > 1.4) { << 189 return fLPMflag;
233 f1 = 42.038 - 8.29*G4Log(delta + 0.958); << 190 }
234 f2 = f1; << 191
235 } else { << 192 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
236 f1 = 42.184 - delta*(7.444 - 1.623*delta); << 193
237 f2 = 41.326 - delta*(5.848 - 0.902*delta); << 194 inline void G4PairProductionRelModel::SetCurrentElement(G4double Z)
>> 195 {
>> 196 if(Z != currentZ) {
>> 197 currentZ = Z;
>> 198
>> 199 G4int iz = G4lrint(Z);
>> 200 z13 = g4calc->Z13(iz);
>> 201 z23 = z13*z13;
>> 202 lnZ = g4calc->logZ(iz);
>> 203
>> 204 if (iz <= 4) {
>> 205 Fel = Fel_light[iz];
>> 206 Finel = Finel_light[iz] ;
>> 207 }
>> 208 else {
>> 209 Fel = facFel - lnZ/3. ;
>> 210 Finel = facFinel - 2.*lnZ/3. ;
>> 211 }
>> 212 fCoulomb=GetCurrentElement()->GetfCoulomb();
238 } 213 }
>> 214 }
>> 215
>> 216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 217
>> 218 inline G4double G4PairProductionRelModel::Phi1(G4double delta) const
>> 219 {
>> 220 return (delta > 1.)
>> 221 ? 21.12 - 4.184*G4Log(delta+0.952)
>> 222 : 20.868 - delta*(3.242 - 0.625*delta);
>> 223 }
>> 224
>> 225 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 226
>> 227 inline G4double G4PairProductionRelModel::Phi2(G4double delta) const
>> 228 {
>> 229 return (delta > 1.)
>> 230 ? 21.12 - 4.184*G4Log(delta+0.952)
>> 231 : 20.209 - delta*(1.930 + 0.086*delta);
>> 232 }
>> 233
>> 234 inline G4double G4PairProductionRelModel::ScreenFunction1(G4double ScreenVariable)
>> 235 // compute the value of the screening function 3*PHI1 - PHI2
>> 236 {
>> 237 return (ScreenVariable > 1.)
>> 238 ? 42.24 - 8.368*G4Log(ScreenVariable+0.952)
>> 239 : 42.392 - ScreenVariable*(7.796 - 1.961*ScreenVariable);
>> 240 }
>> 241
>> 242 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 243
>> 244 inline G4double G4PairProductionRelModel::ScreenFunction2(G4double ScreenVariable)
>> 245 // compute the value of the screening function 1.5*PHI1 + 0.5*PHI2
>> 246 {
>> 247 return (ScreenVariable > 1.)
>> 248 ? 42.24 - 8.368*G4Log(ScreenVariable+0.952)
>> 249 : 41.405 - ScreenVariable*(5.828 - 0.8945*ScreenVariable);
>> 250 }
>> 251
>> 252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 253
>> 254 inline G4double G4PairProductionRelModel::DeltaMax() const
>> 255 {
>> 256 // k > 50 MeV
>> 257 G4double FZ = 8.*(lnZ/3. + fCoulomb);
>> 258 return G4Exp( (42.24-FZ)/8.368 ) + 0.952;
>> 259 }
>> 260
>> 261 inline G4double G4PairProductionRelModel::DeltaMin(G4double k) const
>> 262 {
>> 263 return 544.*CLHEP::electron_mass_c2/(z13*k);
239 } 264 }
240 265
241 #endif 266 #endif
242 267