| 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,v 1.9 2010-10-26 10:35:22 vnivanch Exp $
>> 27 // GEANT4 tag $Name: geant4-09-04-patch-02 $
26 // 28 //
27 // ------------------------------------------- 29 // -------------------------------------------------------------------
28 // 30 //
29 // GEANT4 Class header file 31 // GEANT4 Class header file
30 // 32 //
31 // 33 //
32 // File name: G4PairProductionRelModel 34 // File name: G4PairProductionRelModel
33 // 35 //
34 // Author: Andreas Schaelicke 36 // Author: Andreas Schaelicke
35 // 37 //
36 // Creation date: 02.04.2009 38 // Creation date: 02.04.2009
37 // 39 //
38 // Modifications: 40 // 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 // 41 //
44 // Class Description: 42 // Class Description:
45 // 43 //
46 // Implementation of gamma conversion to e+e- << 44 // Implementation of gamma convertion to e+e- in the field of a nucleus
47 // relativistic approximation 45 // relativistic approximation
48 // 46 //
49 47
50 // ------------------------------------------- 48 // -------------------------------------------------------------------
51 // 49 //
52 50
53 #ifndef G4PairProductionRelModel_h 51 #ifndef G4PairProductionRelModel_h
54 #define G4PairProductionRelModel_h 1 52 #define G4PairProductionRelModel_h 1
55 53
56 #include <CLHEP/Units/PhysicalConstants.h> <<
57 <<
58 #include "G4VEmModel.hh" 54 #include "G4VEmModel.hh"
59 #include "G4Log.hh" << 55 #include "G4PhysicsTable.hh"
60 #include <vector> << 56 #include "G4NistManager.hh"
61 57
62 class G4ParticleChangeForGamma; 58 class G4ParticleChangeForGamma;
63 class G4Pow; <<
64 59
65 class G4PairProductionRelModel : public G4VEmM 60 class G4PairProductionRelModel : public G4VEmModel
66 { 61 {
67 62
68 public: 63 public:
69 64
70 explicit G4PairProductionRelModel(const G4Pa << 65 G4PairProductionRelModel(const G4ParticleDefinition* p = 0,
71 const G4St << 66 const G4String& nam = "BetheHeitlerLPM");
72 67
73 ~G4PairProductionRelModel() override; << 68 virtual ~G4PairProductionRelModel();
74 69
75 void Initialise(const G4ParticleDefinition*, << 70 virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&);
76 71
77 void InitialiseLocal(const G4ParticleDefinit << 72 virtual G4double ComputeCrossSectionPerAtom(
78 G4VEmModel* masterModel) override; << 73 const G4ParticleDefinition*,
>> 74 G4double kinEnergy,
>> 75 G4double Z,
>> 76 G4double A=0.,
>> 77 G4double cut=0.,
>> 78 G4double emax=DBL_MAX);
>> 79
>> 80 virtual void SampleSecondaries(std::vector<G4DynamicParticle*>*,
>> 81 const G4MaterialCutsCouple*,
>> 82 const G4DynamicParticle*,
>> 83 G4double tmin,
>> 84 G4double maxEnergy);
>> 85
>> 86 virtual void SetupForMaterial(const G4ParticleDefinition*,
>> 87 const G4Material*,G4double);
>> 88
>> 89 // * fast inline functions *
>> 90 inline void SetCurrentElement(G4double /*Z*/);
>> 91
>> 92 // set / get methods
>> 93 inline void SetLPMconstant(G4double val);
>> 94 inline G4double LPMconstant() const;
79 95
80 G4double ComputeCrossSectionPerAtom(const G4 << 96 inline void SetLPMflag(G4bool);
81 G4double kinEnergy, << 97 inline G4bool LPMflag() const;
82 G4double Z, <<
83 G4double A=0., <<
84 G4double cut=0., <<
85 G4double emax=DBL_MAX) override; <<
86 <<
87 void SampleSecondaries(std::vector<G4Dynamic <<
88 const G4MaterialCutsCouple*, <<
89 const G4DynamicParticle*, <<
90 G4double tmin, <<
91 G4double maxEnergy) override; <<
92 <<
93 void SetupForMaterial(const G4ParticleDefini <<
94 const G4Material*,G4double) override; <<
95 <<
96 inline void SetLPMflag(G4bool val) { fIsUs <<
97 inline G4bool LPMflag() const { retur <<
98 <<
99 G4PairProductionRelModel & operator= <<
100 (const G4PairProductionRelModel &right) = de <<
101 G4PairProductionRelModel(const G4PairProduc <<
102 98
103 protected: 99 protected:
>> 100 // screening functions
>> 101 inline G4double Phi1(G4double delta) const;
>> 102 inline G4double Phi2(G4double delta) const;
>> 103 inline G4double DeltaMax() const;
>> 104 inline G4double DeltaMin(G4double) const;
>> 105 // lpm functions
>> 106 void CalcLPMFunctions(G4double k, G4double eplus);
>> 107
>> 108 // obsolete
>> 109 G4double ScreenFunction1(G4double ScreenVariable);
>> 110 G4double ScreenFunction2(G4double ScreenVariable);
>> 111
>> 112 G4double ComputeXSectionPerAtom(G4double totalEnergy, G4double Z);
>> 113
>> 114 G4double ComputeDXSectionPerAtom(G4double eplusEnergy, G4double totalEnergy, G4double Z);
>> 115 G4double ComputeRelDXSectionPerAtom(G4double eplusEnergy, G4double totalEnergy, G4double Z);
>> 116
>> 117 // hide assignment operator
>> 118 G4PairProductionRelModel & operator=(const G4PairProductionRelModel &right);
>> 119 G4PairProductionRelModel(const G4PairProductionRelModel&);
>> 120
>> 121 G4NistManager* nist;
>> 122
>> 123 G4ParticleDefinition* theGamma;
>> 124 G4ParticleDefinition* theElectron;
>> 125 G4ParticleDefinition* thePositron;
>> 126 G4ParticleChangeForGamma* fParticleChange;
>> 127
>> 128 G4double fLPMconstant;
>> 129 G4bool fLPMflag;
>> 130
>> 131 // cash
>> 132 G4double z13, z23, lnZ;
>> 133 G4double Fel, Finel, fCoulomb;
>> 134 G4double currentZ;
>> 135
>> 136 // LPM effect
>> 137 G4double lpmEnergy;
>> 138 G4double xiLPM, phiLPM, gLPM;
>> 139
>> 140 // consts
>> 141 G4bool use_completescreening;
>> 142
>> 143 static const G4double xgi[8], wgi[8];
>> 144 static const G4double Fel_light[5];
>> 145 static const G4double Finel_light[5];
>> 146 static const G4double facFel;
>> 147 static const G4double facFinel;
104 148
105 // for evaluating screening related function << 149 static const G4double preS1, logTwo;
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 150
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 }; 151 };
184 // << 152
185 // Bethe screening functions for the elastic ( << 153
186 // Bethe's phi1, phi2 coherent screening funct << 154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
187 // by using (the universal) atomic form factor << 155
188 // Fermi model of the atom (using numerical so << 156 inline
189 // screening function instead of Moliere's ana << 157 void G4PairProductionRelModel::SetLPMconstant(G4double val)
190 // numerical results can be well approximated << 158 {
191 // especially near the delta=1 limit) by: << 159 fLPMconstant = val;
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 } 160 }
213 161
214 // Compute the value of the screening function << 162 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
215 inline G4double G4PairProductionRelModel::Scre << 163
>> 164 inline
>> 165 G4double G4PairProductionRelModel::LPMconstant() const
216 { 166 {
217 return (delta > 1.4) ? 42.038 - 8.29*G4Log(d << 167 return fLPMconstant;
218 : 42.184 - delta*(7.444 <<
219 } 168 }
220 169
221 // Compute the value of the screening function << 170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
222 inline G4double G4PairProductionRelModel::Scre << 171
>> 172 inline
>> 173 void G4PairProductionRelModel::SetLPMflag(G4bool val)
>> 174 {
>> 175 fLPMflag = val;
>> 176 }
>> 177
>> 178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 179
>> 180 inline
>> 181 G4bool G4PairProductionRelModel::LPMflag() const
223 { 182 {
224 return (delta > 1.4) ? 42.038 - 8.29*G4Log(d << 183 return fLPMflag;
225 : 41.326 - delta*(5.848 <<
226 } 184 }
227 185
228 // Same as ScreenFunction1 and ScreenFunction2 << 186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
229 inline void G4PairProductionRelModel::ScreenFu << 187
230 << 188 inline void G4PairProductionRelModel::SetCurrentElement(G4double Z)
231 { 189 {
232 if (delta > 1.4) { << 190 if(Z != currentZ) {
233 f1 = 42.038 - 8.29*G4Log(delta + 0.958); << 191 currentZ = Z;
234 f2 = f1; << 192
235 } else { << 193 G4int iz = G4int(Z);
236 f1 = 42.184 - delta*(7.444 - 1.623*delta); << 194 z13 = nist->GetZ13(iz);
237 f2 = 41.326 - delta*(5.848 - 0.902*delta); << 195 z23 = z13*z13;
>> 196 lnZ = nist->GetLOGZ(iz);
>> 197
>> 198 if (iz <= 4) {
>> 199 Fel = Fel_light[iz];
>> 200 Finel = Finel_light[iz] ;
>> 201 }
>> 202 else {
>> 203 Fel = facFel - lnZ/3. ;
>> 204 Finel = facFinel - 2.*lnZ/3. ;
>> 205 }
>> 206 fCoulomb=GetCurrentElement()->GetfCoulomb();
238 } 207 }
239 } 208 }
>> 209
>> 210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 211
>> 212 inline G4double G4PairProductionRelModel::Phi1(G4double delta) const
>> 213 {
>> 214 G4double screenVal;
>> 215
>> 216 if (delta > 1.)
>> 217 screenVal = 21.12 - 4.184*std::log(delta+0.952);
>> 218 else
>> 219 screenVal = 20.868 - delta*(3.242 - 0.625*delta);
>> 220
>> 221 return screenVal;
>> 222 }
>> 223
>> 224 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 225
>> 226 inline G4double G4PairProductionRelModel::Phi2(G4double delta) const
>> 227 {
>> 228 G4double screenVal;
>> 229
>> 230 if (delta > 1.)
>> 231 screenVal = 21.12 - 4.184*std::log(delta+0.952);
>> 232 else
>> 233 screenVal = 20.209 - delta*(1.930 + 0.086*delta);
>> 234
>> 235 return screenVal;
>> 236 }
>> 237
>> 238 inline G4double G4PairProductionRelModel::ScreenFunction1(G4double ScreenVariable)
>> 239
>> 240 // compute the value of the screening function 3*PHI1 - PHI2
>> 241
>> 242 {
>> 243 G4double screenVal;
>> 244
>> 245 if (ScreenVariable > 1.)
>> 246 screenVal = 42.24 - 8.368*std::log(ScreenVariable+0.952);
>> 247 else
>> 248 screenVal = 42.392 - ScreenVariable*(7.796 - 1.961*ScreenVariable);
>> 249
>> 250 return screenVal;
>> 251 }
>> 252
>> 253 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 254
>> 255 inline G4double G4PairProductionRelModel::ScreenFunction2(G4double ScreenVariable)
>> 256
>> 257 // compute the value of the screening function 1.5*PHI1 + 0.5*PHI2
>> 258
>> 259 {
>> 260 G4double screenVal;
>> 261
>> 262 if (ScreenVariable > 1.)
>> 263 screenVal = 42.24 - 8.368*std::log(ScreenVariable+0.952);
>> 264 else
>> 265 screenVal = 41.405 - ScreenVariable*(5.828 - 0.8945*ScreenVariable);
>> 266
>> 267 return screenVal;
>> 268 }
>> 269
>> 270 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 271
>> 272
>> 273 inline G4double G4PairProductionRelModel::DeltaMax() const
>> 274 {
>> 275 // k > 50 MeV
>> 276 G4double FZ = 8.*(lnZ/3. + fCoulomb);
>> 277 return std::exp( (42.24-FZ)/8.368 ) + 0.952;
>> 278 }
>> 279
>> 280 inline G4double G4PairProductionRelModel::DeltaMin(G4double k) const
>> 281 {
>> 282 return 4.*136./z13*(electron_mass_c2/k);
>> 283 }
>> 284
>> 285
240 286
241 #endif 287 #endif
242 288