| 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 // 26 //
27 // ------------------------------------------- 27 // -------------------------------------------------------------------
28 // 28 //
29 // GEANT4 Class file 29 // GEANT4 Class file
30 // 30 //
31 // 31 //
32 // File name: G4eBremsstrahlungRelModel 32 // File name: G4eBremsstrahlungRelModel
33 // 33 //
34 // Author: Andreas Schaelicke 34 // Author: Andreas Schaelicke
35 // 35 //
36 // Creation date: 12.08.2008 36 // Creation date: 12.08.2008
37 // 37 //
38 // Modifications: 38 // Modifications:
39 // 39 //
40 // 13.11.08 add SetLPMflag and SetLPMconsta 40 // 13.11.08 add SetLPMflag and SetLPMconstant methods
41 // 13.11.08 change default LPMconstant valu 41 // 13.11.08 change default LPMconstant value
42 // 13.10.10 add angular distributon interfa 42 // 13.10.10 add angular distributon interface (VI)
43 // 31.05.16 change LPMconstant such that it 43 // 31.05.16 change LPMconstant such that it gives suppression variable 's'
44 // that consistent to Migdal's one 44 // that consistent to Migdal's one; fix a small bug in 'logTS1'
45 // computation; better agreement w 45 // computation; better agreement with exp.(M.Novak)
46 // 15.07.18 improved LPM suppression functi 46 // 15.07.18 improved LPM suppression function approximation (no artificial
47 // steps), code cleanup and optimi 47 // steps), code cleanup and optimizations,more implementation and
48 // model related comments, consist 48 // model related comments, consistent variable naming (M.Novak)
49 // 49 //
50 // Main References: 50 // Main References:
51 // Y.-S.Tsai, Rev. Mod. Phys. 46 (1974) 815; 51 // Y.-S.Tsai, Rev. Mod. Phys. 46 (1974) 815; Rev. Mod. Phys. 49 (1977) 421.
52 // S.Klein, Rev. Mod. Phys. 71 (1999) 1501. 52 // S.Klein, Rev. Mod. Phys. 71 (1999) 1501.
53 // T.Stanev et.al., Phys. Rev. D25 (1982) 129 53 // T.Stanev et.al., Phys. Rev. D25 (1982) 1291.
54 // M.L.Ter-Mikaelian, High-energy Electromagn 54 // M.L.Ter-Mikaelian, High-energy Electromagnetic Processes in Condensed Media,
55 // Wiley, 1972. 55 // Wiley, 1972.
56 // 56 //
57 // ------------------------------------------- 57 // -------------------------------------------------------------------
58 // 58 //
59 59
60 #include "G4eBremsstrahlungRelModel.hh" 60 #include "G4eBremsstrahlungRelModel.hh"
61 #include "G4PhysicalConstants.hh" 61 #include "G4PhysicalConstants.hh"
62 #include "G4SystemOfUnits.hh" 62 #include "G4SystemOfUnits.hh"
63 #include "G4Electron.hh" 63 #include "G4Electron.hh"
64 #include "G4Gamma.hh" 64 #include "G4Gamma.hh"
65 #include "Randomize.hh" 65 #include "Randomize.hh"
66 #include "G4Material.hh" 66 #include "G4Material.hh"
67 #include "G4Element.hh" 67 #include "G4Element.hh"
68 #include "G4ElementVector.hh" 68 #include "G4ElementVector.hh"
69 #include "G4ParticleChangeForLoss.hh" 69 #include "G4ParticleChangeForLoss.hh"
70 #include "G4ModifiedTsai.hh" 70 #include "G4ModifiedTsai.hh"
71 #include "G4Exp.hh" << 71 //#include "G4DipBustGenerator.hh"
72 #include "G4Log.hh" <<
73 #include "G4Pow.hh" <<
74 #include "G4EmParameters.hh" <<
75 #include "G4AutoLock.hh" <<
76 #include <thread> <<
77 72
78 const G4int G4eBremsstrahlungRelModel::gMaxZet 73 const G4int G4eBremsstrahlungRelModel::gMaxZet = 120;
79 74
80 // constant DCS factor: 16\alpha r_0^2/3 75 // constant DCS factor: 16\alpha r_0^2/3
81 const G4double G4eBremsstrahlungRelModel::gBre 76 const G4double G4eBremsstrahlungRelModel::gBremFactor
82 = 16. * CLHEP::fine_structure_const * CLHEP: 77 = 16. * CLHEP::fine_structure_const * CLHEP::classic_electr_radius
83 * CLHEP::classic_electr_radius/3.; 78 * CLHEP::classic_electr_radius/3.;
84 79
85 // Migdal's constant: 4\pi r_0*electron_reduce 80 // Migdal's constant: 4\pi r_0*electron_reduced_compton_wavelength^2
86 const G4double G4eBremsstrahlungRelModel::gMig 81 const G4double G4eBremsstrahlungRelModel::gMigdalConstant
87 = 4. * CLHEP::pi * CLHEP::classic_electr_rad 82 = 4. * CLHEP::pi * CLHEP::classic_electr_radius
88 * CLHEP::electron_Compton_length * CLHEP:: 83 * CLHEP::electron_Compton_length * CLHEP::electron_Compton_length;
89 84
90 // LPM constant: \alpha(mc^2)^2/(4\pi*\hbar c) 85 // LPM constant: \alpha(mc^2)^2/(4\pi*\hbar c)
91 const G4double G4eBremsstrahlungRelModel::gLPM 86 const G4double G4eBremsstrahlungRelModel::gLPMconstant
92 = CLHEP::fine_structure_const * CLHEP::elect 87 = CLHEP::fine_structure_const * CLHEP::electron_mass_c2
93 * CLHEP::electron_mass_c2 / (4. * CLHEP::p 88 * CLHEP::electron_mass_c2 / (4. * CLHEP::pi * CLHEP::hbarc);
94 89
95 // abscissas and weights of an 8 point Gauss-L 90 // abscissas and weights of an 8 point Gauss-Legendre quadrature
96 // for numerical integration on [0,1] 91 // for numerical integration on [0,1]
97 const G4double G4eBremsstrahlungRelModel::gXGL 92 const G4double G4eBremsstrahlungRelModel::gXGL[] = {
98 1.98550718e-02, 1.01666761e-01, 2.37233795e- 93 1.98550718e-02, 1.01666761e-01, 2.37233795e-01, 4.08282679e-01,
99 5.91717321e-01, 7.62766205e-01, 8.98333239e- 94 5.91717321e-01, 7.62766205e-01, 8.98333239e-01, 9.80144928e-01
100 }; 95 };
101 const G4double G4eBremsstrahlungRelModel::gWGL 96 const G4double G4eBremsstrahlungRelModel::gWGL[] = {
102 5.06142681e-02, 1.11190517e-01, 1.56853323e- 97 5.06142681e-02, 1.11190517e-01, 1.56853323e-01, 1.81341892e-01,
103 1.81341892e-01, 1.56853323e-01, 1.11190517e- 98 1.81341892e-01, 1.56853323e-01, 1.11190517e-01, 5.06142681e-02
104 }; 99 };
105 100
106 // elastic and inelatic radiation logarithms f 101 // elastic and inelatic radiation logarithms for light elements (where the
107 // Thomas-Fermi model doesn't work): computed 102 // Thomas-Fermi model doesn't work): computed by using Dirac-Fock model of atom.
108 const G4double G4eBremsstrahlungRelModel::gFel 103 const G4double G4eBremsstrahlungRelModel::gFelLowZet [] = {
109 0.0, 5.3104, 4.7935, 4.7402, 4.7112, 4.6694, 104 0.0, 5.3104, 4.7935, 4.7402, 4.7112, 4.6694, 4.6134, 4.5520
110 }; 105 };
111 const G4double G4eBremsstrahlungRelModel::gFin 106 const G4double G4eBremsstrahlungRelModel::gFinelLowZet[] = {
112 0.0, 5.9173, 5.6125, 5.5377, 5.4728, 5.4174, 107 0.0, 5.9173, 5.6125, 5.5377, 5.4728, 5.4174, 5.3688, 5.3236
113 }; 108 };
114 109
115 // LPM supression functions evaluated at initi 110 // LPM supression functions evaluated at initialisation time
116 std::shared_ptr<G4eBremsstrahlungRelModel::LPM << 111 G4eBremsstrahlungRelModel::LPMFuncs G4eBremsstrahlungRelModel::gLPMFuncs;
117 { <<
118 // We have to use shared pointer for the LPM <<
119 // by the G4eBremsstrahlungRelModel used in <<
120 // model is owned (well deleted) by (at leas <<
121 // a G4SeltzerBergerModel which is owned by <<
122 // which owned by a G4ThreadLocalSingleton<G <<
123 // which is a static global and thus deleted <<
124 // is deleted. <<
125 static auto _instance = std::make_shared<G4e <<
126 return _instance; <<
127 } <<
128 112
129 // special data structure per element i.e. per 113 // special data structure per element i.e. per Z
130 std::shared_ptr<std::vector<G4eBremsstrahlungR << 114 std::vector<G4eBremsstrahlungRelModel::ElementData*> G4eBremsstrahlungRelModel::gElementData;
131 { <<
132 // Same code comment as for gLPMFuncs. <<
133 static auto _instance = std::make_shared<std <<
134 return _instance; <<
135 } <<
136 <<
137 static std::once_flag applyOnce; <<
138 <<
139 namespace <<
140 { <<
141 G4Mutex theBremRelMutex = G4MUTEX_INITIALIZE <<
142 } <<
143 115
144 G4eBremsstrahlungRelModel::G4eBremsstrahlungRe 116 G4eBremsstrahlungRelModel::G4eBremsstrahlungRelModel(const G4ParticleDefinition* p,
145 117 const G4String& nam)
146 : G4VEmModel(nam), fLPMFuncs(gLPMFuncs()), fEl << 118 : G4VEmModel(nam), fIsElectron(true), fIsScatOffElectron(false),
>> 119 fIsLPMActive(false), fPrimaryParticle(nullptr), fIsUseCompleteScreening(false)
147 { 120 {
>> 121 fCurrentIZ = 0;
>> 122 //
>> 123 fPrimaryParticleMass = 0.;
>> 124 fPrimaryKinEnergy = 0.;
>> 125 fPrimaryTotalEnergy = 0.;
>> 126 fDensityFactor = 0.;
>> 127 fDensityCorr = 0.;
>> 128 fNucTerm = 0.;
>> 129 fSumTerm = 0.;
>> 130 //
>> 131 fPrimaryParticle = nullptr;
148 fGammaParticle = G4Gamma::Gamma(); 132 fGammaParticle = G4Gamma::Gamma();
>> 133 fParticleChange = nullptr;
149 // 134 //
150 fLowestKinEnergy = 1.0*CLHEP::MeV; << 135 fLowestKinEnergy = 1.0*MeV;
151 SetLowEnergyLimit(fLowestKinEnergy); 136 SetLowEnergyLimit(fLowestKinEnergy);
152 // 137 //
153 fLPMEnergyThreshold = 1.e+39; 138 fLPMEnergyThreshold = 1.e+39;
154 fLPMEnergy = 0.; 139 fLPMEnergy = 0.;
>> 140
>> 141 SetLPMFlag(true);
>> 142 //
155 SetAngularDistribution(new G4ModifiedTsai()) 143 SetAngularDistribution(new G4ModifiedTsai());
>> 144 //SetAngularDistribution(new G4DipBustGenerator());
156 // 145 //
157 if (nullptr != p) { << 146 if (p) {
158 SetParticle(p); 147 SetParticle(p);
159 } 148 }
160 } 149 }
161 150
162 G4eBremsstrahlungRelModel::~G4eBremsstrahlungR 151 G4eBremsstrahlungRelModel::~G4eBremsstrahlungRelModel()
163 { 152 {
164 if (fIsInitializer) { << 153 if (IsMaster()) {
165 // clear ElementData container 154 // clear ElementData container
166 for (auto const & ptr : *fElementData) { d << 155 for (size_t iz = 0; iz < gElementData.size(); ++iz) {
167 fElementData->clear(); << 156 if (gElementData[iz]) {
>> 157 delete gElementData[iz];
>> 158 }
>> 159 }
>> 160 gElementData.clear();
168 // clear LPMFunctions (if any) 161 // clear LPMFunctions (if any)
169 if (fLPMFuncs->fIsInitialized) { << 162 if (LPMFlag()) {
170 fLPMFuncs->fLPMFuncG.clear(); << 163 gLPMFuncs.fLPMFuncG.clear();
171 fLPMFuncs->fLPMFuncPhi.clear(); << 164 gLPMFuncs.fLPMFuncPhi.clear();
172 fLPMFuncs->fIsInitialized = false; << 165 gLPMFuncs.fIsInitialized = false;
173 } 166 }
174 } 167 }
175 } 168 }
176 169
177 void G4eBremsstrahlungRelModel::Initialise(con 170 void G4eBremsstrahlungRelModel::Initialise(const G4ParticleDefinition* p,
178 con 171 const G4DataVector& cuts)
179 { 172 {
180 // parameters in each thread << 173 if (p) {
181 if (fPrimaryParticle != p) { <<
182 SetParticle(p); 174 SetParticle(p);
183 } 175 }
184 fUseLPM = G4EmParameters::Instance()->LPM(); <<
185 fCurrentIZ = 0; 176 fCurrentIZ = 0;
186 << 177 // init element data and precompute LPM functions (only if lpmflag is true)
187 // init static element data and precompute L <<
188 std::call_once(applyOnce, [this]() { fIsInit <<
189 <<
190 // for all treads and derived classes <<
191 if (fIsInitializer || fElementData->empty()) <<
192 G4AutoLock l(&theBremRelMutex); <<
193 if (fElementData->empty()) { <<
194 fElementData->resize(gMaxZet+1, nullptr) <<
195 } <<
196 InitialiseElementData(); <<
197 InitLPMFunctions(); <<
198 l.unlock(); <<
199 } <<
200 <<
201 // element selectors are initialized in the <<
202 if (IsMaster()) { 178 if (IsMaster()) {
203 InitialiseElementSelectors(p, cuts); << 179 InitialiseElementData();
204 } << 180 if (LPMFlag()) { InitLPMFunctions(); }
205 // initialisation in all threads << 181 if (LowEnergyLimit() < HighEnergyLimit()) {
206 if (nullptr == fParticleChange) { << 182 InitialiseElementSelectors(p, cuts);
207 fParticleChange = GetParticleChangeForLoss << 183 }
208 } 184 }
>> 185 if (!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
209 if (GetTripletModel()) { 186 if (GetTripletModel()) {
210 GetTripletModel()->Initialise(p, cuts); 187 GetTripletModel()->Initialise(p, cuts);
211 fIsScatOffElectron = true; 188 fIsScatOffElectron = true;
212 } 189 }
213 } 190 }
214 191
215 void G4eBremsstrahlungRelModel::InitialiseLoca 192 void G4eBremsstrahlungRelModel::InitialiseLocal(const G4ParticleDefinition*,
216 193 G4VEmModel* masterModel)
217 { 194 {
218 SetElementSelectors(masterModel->GetElementS << 195 if (LowEnergyLimit() < HighEnergyLimit()) {
>> 196 SetElementSelectors(masterModel->GetElementSelectors());
>> 197 }
219 } 198 }
220 199
221 void G4eBremsstrahlungRelModel::SetParticle(co 200 void G4eBremsstrahlungRelModel::SetParticle(const G4ParticleDefinition* p)
222 { 201 {
223 fPrimaryParticle = p; 202 fPrimaryParticle = p;
224 fPrimaryParticleMass = p->GetPDGMass(); 203 fPrimaryParticleMass = p->GetPDGMass();
225 fIsElectron = (p==G4Electron::Elect 204 fIsElectron = (p==G4Electron::Electron());
226 } 205 }
227 206
228 // Sets kinematical variables like E_kin, E_t 207 // Sets kinematical variables like E_kin, E_t and some material dependent
229 // variables like LPM energy and characteristi 208 // variables like LPM energy and characteristic photon energy k_p (more exactly
230 // k_p^2) for the Ter-Mikaelian suppression ef 209 // k_p^2) for the Ter-Mikaelian suppression effect.
231 void G4eBremsstrahlungRelModel::SetupForMateri 210 void G4eBremsstrahlungRelModel::SetupForMaterial(const G4ParticleDefinition*,
232 211 const G4Material* mat,
233 212 G4double kineticEnergy)
234 { 213 {
235 fDensityFactor = gMigdalConstant*mat->GetEle 214 fDensityFactor = gMigdalConstant*mat->GetElectronDensity();
236 fLPMEnergy = gLPMconstant*mat->GetRadlen 215 fLPMEnergy = gLPMconstant*mat->GetRadlen();
237 // threshold for LPM effect (i.e. below whic 216 // threshold for LPM effect (i.e. below which LPM hidden by density effect)
238 if (fUseLPM) { << 217 if (LPMFlag()) {
239 fLPMEnergyThreshold = std::sqrt(fDensityFa 218 fLPMEnergyThreshold = std::sqrt(fDensityFactor)*fLPMEnergy;
240 } else { 219 } else {
241 fLPMEnergyThreshold = 1.e+39; // i.e. do 220 fLPMEnergyThreshold = 1.e+39; // i.e. do not use LPM effect
242 } 221 }
243 // calculate threshold for density effect: k 222 // calculate threshold for density effect: k_p = sqrt(fDensityCorr)
244 fPrimaryKinEnergy = kineticEnergy; 223 fPrimaryKinEnergy = kineticEnergy;
245 fPrimaryTotalEnergy = kineticEnergy+fPrimary 224 fPrimaryTotalEnergy = kineticEnergy+fPrimaryParticleMass;
246 fDensityCorr = fDensityFactor*fPrimar 225 fDensityCorr = fDensityFactor*fPrimaryTotalEnergy*fPrimaryTotalEnergy;
247 // set activation flag for LPM effects in th 226 // set activation flag for LPM effects in the DCS
248 fIsLPMActive = (fPrimaryTotalEnergy>fLPMEner << 227 fIsLPMActive = (fPrimaryTotalEnergy>fLPMEnergyThreshold);
249 } 228 }
250 229
251 // minimum primary (e-/e+) energy at which dis 230 // minimum primary (e-/e+) energy at which discrete interaction is possible
252 G4double G4eBremsstrahlungRelModel::MinPrimary 231 G4double G4eBremsstrahlungRelModel::MinPrimaryEnergy(const G4Material*,
253 232 const G4ParticleDefinition*,
254 233 G4double cut)
255 { 234 {
256 return std::max(fLowestKinEnergy, cut); 235 return std::max(fLowestKinEnergy, cut);
257 } 236 }
258 237
259 // Computes the restricted dE/dx as the approp 238 // Computes the restricted dE/dx as the appropriate weight of the individual
260 // element contributions that are computed by 239 // element contributions that are computed by numerically integrating the DCS.
261 G4double 240 G4double
262 G4eBremsstrahlungRelModel::ComputeDEDXPerVolum 241 G4eBremsstrahlungRelModel::ComputeDEDXPerVolume(const G4Material* material,
263 242 const G4ParticleDefinition* p,
264 243 G4double kineticEnergy,
265 244 G4double cutEnergy)
266 { 245 {
267 G4double dedx = 0.0; 246 G4double dedx = 0.0;
268 if (nullptr == fPrimaryParticle) { << 247 if (!fPrimaryParticle) {
269 SetParticle(p); 248 SetParticle(p);
270 } 249 }
271 if (kineticEnergy < LowEnergyLimit()) { 250 if (kineticEnergy < LowEnergyLimit()) {
272 return dedx; 251 return dedx;
273 } 252 }
274 // maximum value of the dE/dx integral (the 253 // maximum value of the dE/dx integral (the minimum is 0 of course)
275 G4double tmax = std::min(cutEnergy, kineticE 254 G4double tmax = std::min(cutEnergy, kineticEnergy);
276 if (tmax == 0.0) { 255 if (tmax == 0.0) {
277 return dedx; 256 return dedx;
278 } 257 }
279 // sets kinematical and material related var 258 // sets kinematical and material related variables
280 SetupForMaterial(fPrimaryParticle, material, 259 SetupForMaterial(fPrimaryParticle, material,kineticEnergy);
281 // get element compositions of the material 260 // get element compositions of the material
282 const G4ElementVector* theElemVector = mater 261 const G4ElementVector* theElemVector = material->GetElementVector();
283 const G4double* theAtomNumDensVector = mater 262 const G4double* theAtomNumDensVector = material->GetAtomicNumDensityVector();
284 const std::size_t numberOfElements = theElem << 263 const size_t numberOfElements = theElemVector->size();
285 // loop over the elements of the material an 264 // loop over the elements of the material and compute their contributions to
286 // the restricted dE/dx by numerical integra 265 // the restricted dE/dx by numerical integration of the dependent part of DCS
287 for (std::size_t ie = 0; ie < numberOfElemen << 266 for (size_t ie = 0; ie < numberOfElements; ++ie) {
288 G4VEmModel::SetCurrentElement((*theElemVec 267 G4VEmModel::SetCurrentElement((*theElemVector)[ie]);
289 G4int zet = (*theElemVector)[ie]->GetZasIn << 268 //SetCurrentElement((*theElementVector)[i]->GetZasInt());
290 fCurrentIZ = std::min(zet, gMaxZet); << 269 const G4double zet = (*theElemVector)[ie]->GetZ();
291 dedx += (zet*zet)*theAtomNumD << 270 fCurrentIZ = std::min(G4lrint(zet), gMaxZet);
>> 271 dedx += theAtomNumDensVector[ie]*zet*zet*ComputeBremLoss(tmax);
292 } 272 }
293 // apply the constant factor C/Z = 16\alpha 273 // apply the constant factor C/Z = 16\alpha r_0^2/3
294 dedx *= gBremFactor; 274 dedx *= gBremFactor;
295 return std::max(dedx,0.); 275 return std::max(dedx,0.);
296 } 276 }
297 277
298 // Computes the integral part of the restricte 278 // Computes the integral part of the restricted dE/dx contribution from a given
299 // element (Z) by numerically integrating the 279 // element (Z) by numerically integrating the k dependent part of the DCS between
300 // k_min=0 and k_max = tmax = min[gamma-cut, e 280 // k_min=0 and k_max = tmax = min[gamma-cut, electron-kinetic-eenrgy].
301 // The numerical integration is done by dividi 281 // The numerical integration is done by dividing the integration range into 'n'
302 // subintervals and an 8 pint GL integral (on 282 // subintervals and an 8 pint GL integral (on [0,1]) is performed on each sub-
303 // inteval by tranforming k to alpha=k/E_t (E_ 283 // inteval by tranforming k to alpha=k/E_t (E_t is the total energy of the e-)
304 // and each sub-interavl is transformed to [0, 284 // and each sub-interavl is transformed to [0,1]. So the integrastion is done
305 // in xi(alpha) = xi(k) = [k/E_t-alpha_i]/delt 285 // in xi(alpha) = xi(k) = [k/E_t-alpha_i]/delta where alpha_i=(i-1)*delta for
306 // the i = 1,2,..,n-th sub-interval so xi(k) i 286 // the i = 1,2,..,n-th sub-interval so xi(k) in [0,1] on each sub-intevals.
307 // This transformation from 'k' to 'xi(k)' res 287 // This transformation from 'k' to 'xi(k)' results in a multiplicative factor
308 // of E_t*delta at each step. 288 // of E_t*delta at each step.
309 // The restricted dE/dx = N int_{0}^{k_max} k* 289 // The restricted dE/dx = N int_{0}^{k_max} k*ds/dk dk. There are 2 DCS model
310 // one with LPM and one without LPM effects (s 290 // one with LPM and one without LPM effects (see them below). In both case not
311 // the ds/dk(Z,k) but ds/dk(Z,k)*[F*k/C] is co 291 // the ds/dk(Z,k) but ds/dk(Z,k)*[F*k/C] is computed since:
312 // (i) what we need here is ds/dk*k and not 292 // (i) what we need here is ds/dk*k and not k so this multiplication is done
313 // (ii) the Ter-Mikaelian suppression i.e. F 293 // (ii) the Ter-Mikaelian suppression i.e. F related factor is done here
314 // (iii) the constant factor C (includes Z^2 294 // (iii) the constant factor C (includes Z^2 as well)is accounted in the caller
315 G4double G4eBremsstrahlungRelModel::ComputeBre 295 G4double G4eBremsstrahlungRelModel::ComputeBremLoss(G4double tmax)
316 { 296 {
317 // number of intervals and integration step 297 // number of intervals and integration step
318 const G4double alphaMax = tmax/fPrimaryTotal 298 const G4double alphaMax = tmax/fPrimaryTotalEnergy;
319 const G4int nSub = (G4int)(20*alphaMa 299 const G4int nSub = (G4int)(20*alphaMax)+3;
320 const G4double delta = alphaMax/((G4doubl 300 const G4double delta = alphaMax/((G4double)nSub);
321 // set minimum value of the first sub-inteva 301 // set minimum value of the first sub-inteval
322 G4double alpha_i = 0.0; 302 G4double alpha_i = 0.0;
323 G4double dedxInteg = 0.0; 303 G4double dedxInteg = 0.0;
324 for (G4int l = 0; l < nSub; ++l) { 304 for (G4int l = 0; l < nSub; ++l) {
325 for (G4int igl = 0; igl < 8; ++igl) { 305 for (G4int igl = 0; igl < 8; ++igl) {
326 // compute the emitted photon energy k 306 // compute the emitted photon energy k
327 const G4double k = (alpha_i+gXGL[igl]* 307 const G4double k = (alpha_i+gXGL[igl]*delta)*fPrimaryTotalEnergy;
328 // compute the DCS value at k (without t 308 // compute the DCS value at k (without the constant, the 1/k, 1/F factors)
329 const G4double dcs = fIsLPMActive 309 const G4double dcs = fIsLPMActive
330 ? ComputeRelDXSectio 310 ? ComputeRelDXSectionPerAtom(k) // DCS WITHOUT LPM
331 : ComputeDXSectionPe 311 : ComputeDXSectionPerAtom(k); // DCS WITH LPM
332 // account Ter-Mikaelian suppression: ti 312 // account Ter-Mikaelian suppression: times 1/F with F = 1+(k_p/k)^2
333 dedxInteg += gWGL[igl]*dcs/(1.0+fDensity 313 dedxInteg += gWGL[igl]*dcs/(1.0+fDensityCorr/(k*k));
334 } 314 }
335 // update sub-interval minimum value 315 // update sub-interval minimum value
336 alpha_i += delta; 316 alpha_i += delta;
337 } 317 }
338 // apply corrections due to variable transfo 318 // apply corrections due to variable transformation i.e. E_t*delta
339 dedxInteg *= delta*fPrimaryTotalEnergy; 319 dedxInteg *= delta*fPrimaryTotalEnergy;
340 return std::max(dedxInteg,0.); 320 return std::max(dedxInteg,0.);
341 } 321 }
342 322
343 // Computes restrected atomic cross section by 323 // Computes restrected atomic cross section by numerically integrating the
344 // DCS between the proper kinematical limits a 324 // DCS between the proper kinematical limits accounting the gamma production cut
345 G4double G4eBremsstrahlungRelModel::ComputeCro 325 G4double G4eBremsstrahlungRelModel::ComputeCrossSectionPerAtom(
346 326 const G4ParticleDefinition* p,
347 327 G4double kineticEnergy,
348 328 G4double Z,
349 329 G4double,
350 330 G4double cut,
351 331 G4double maxEnergy)
352 { 332 {
353 G4double crossSection = 0.0; 333 G4double crossSection = 0.0;
354 if (nullptr == fPrimaryParticle) { << 334 if (!fPrimaryParticle) {
355 SetParticle(p); 335 SetParticle(p);
356 } 336 }
357 if (kineticEnergy < LowEnergyLimit()) { 337 if (kineticEnergy < LowEnergyLimit()) {
358 return crossSection; 338 return crossSection;
359 } 339 }
360 // min/max kinetic energy limits of the DCS 340 // min/max kinetic energy limits of the DCS integration:
361 const G4double tmin = std::min(cut, kineticE 341 const G4double tmin = std::min(cut, kineticEnergy);
362 const G4double tmax = std::min(maxEnergy, ki 342 const G4double tmax = std::min(maxEnergy, kineticEnergy);
363 // zero restricted x-section if e- kinetic e 343 // zero restricted x-section if e- kinetic energy is below gamma cut
364 if (tmin >= tmax) { 344 if (tmin >= tmax) {
365 return crossSection; 345 return crossSection;
366 } 346 }
367 fCurrentIZ = std::min(G4lrint(Z), gMaxZet); << 347 fCurrentIZ = std::min(G4lrint(Z), gMaxZet);
368 // integrate numerically (dependent part of) 348 // integrate numerically (dependent part of) the DCS between the kin. limits:
369 // a. integrate between tmin and kineticEner 349 // a. integrate between tmin and kineticEnergy of the e-
370 crossSection = ComputeXSectionPerAtom(tmin); 350 crossSection = ComputeXSectionPerAtom(tmin);
371 // allow partial integration: only if maxEne 351 // allow partial integration: only if maxEnergy < kineticEnergy
372 // b. integrate between tmax and kineticEner 352 // b. integrate between tmax and kineticEnergy (tmax=maxEnergy in this case)
373 // (so the result in this case is the integr 353 // (so the result in this case is the integral of DCS between tmin and
374 // maxEnergy) 354 // maxEnergy)
375 if (tmax < kineticEnergy) { 355 if (tmax < kineticEnergy) {
376 crossSection -= ComputeXSectionPerAtom(tma 356 crossSection -= ComputeXSectionPerAtom(tmax);
377 } 357 }
378 // multiply with the constant factors: 16\al 358 // multiply with the constant factors: 16\alpha r_0^2/3 Z^2
379 crossSection *= Z*Z*gBremFactor; 359 crossSection *= Z*Z*gBremFactor;
380 return std::max(crossSection, 0.); 360 return std::max(crossSection, 0.);
381 } 361 }
382 362
383 // Numerical integral of the (k dependent part 363 // Numerical integral of the (k dependent part of) DCS between k_min=tmin and
384 // k_max = E_k (where E_k is the kinetic energ 364 // k_max = E_k (where E_k is the kinetic energy of the e- and tmin is the
385 // minimum of energy of the emitted photon). 365 // minimum of energy of the emitted photon). The integration is done in the
386 // transformed alpha(k) = ln(k/E_t) variable ( 366 // transformed alpha(k) = ln(k/E_t) variable (with E_t being the total energy of
387 // the primary e-). The integration range is d 367 // the primary e-). The integration range is divided into n sub-intervals with
388 // delta = [ln(k_min/E_t)-ln(k_max/E_t)]/n wid 368 // delta = [ln(k_min/E_t)-ln(k_max/E_t)]/n width each. An 8 point GL integral
389 // on [0,1] is applied on each sub-inteval so 369 // on [0,1] is applied on each sub-inteval so alpha is transformed to
390 // xi(alpha) = xi(k) = [ln(k/E_t)-alpha_i]/del 370 // xi(alpha) = xi(k) = [ln(k/E_t)-alpha_i]/delta where alpha_i = ln(k_min/E_t) +
391 // (i-1)*delta for the i = 1,2,..,n-th sub-int 371 // (i-1)*delta for the i = 1,2,..,n-th sub-interval and xi(k) in [0,1] on each
392 // sub-intevals. From the transformed xi, k(xi 372 // sub-intevals. From the transformed xi, k(xi) = E_t exp[xi*delta+alpha_i].
393 // Since the integration is done in variable x 373 // Since the integration is done in variable xi instead of k this
394 // transformation results in a multiplicative 374 // transformation results in a multiplicative factor of k*delta at each step.
395 // However, DCS differential in k is ~1/k so t 375 // However, DCS differential in k is ~1/k so the multiplicative factor is simple
396 // becomes delta and the 1/k factor is dropped 376 // becomes delta and the 1/k factor is dropped from the DCS computation.
397 // NOTE: 377 // NOTE:
398 // - LPM suppression is accounted above thre 378 // - LPM suppression is accounted above threshold e- energy (corresponidng
399 // flag is set in SetUpForMaterial() => 2 379 // flag is set in SetUpForMaterial() => 2 DCS with/without LPM
400 // - Ter-Mikaelian suppression is always acc 380 // - Ter-Mikaelian suppression is always accounted
401 G4double G4eBremsstrahlungRelModel::ComputeXSe 381 G4double G4eBremsstrahlungRelModel::ComputeXSectionPerAtom(G4double tmin)
402 { 382 {
403 G4double xSection = 0.0; 383 G4double xSection = 0.0;
404 const G4double alphaMin = G4Log(tmin/fPrimar 384 const G4double alphaMin = G4Log(tmin/fPrimaryTotalEnergy);
405 const G4double alphaMax = G4Log(fPrimaryKinE << 385 const G4double alphaMax = G4Log(fPrimaryKinEnergy/fPrimaryTotalEnergy);
406 const G4int nSub = std::max((G4int)(0.45*alp << 386 const G4int nSub = (G4int)(0.45*(alphaMax-alphaMin))+4;
407 const G4double delta = alphaMax/((G4double)n << 387 const G4double delta = (alphaMax-alphaMin)/((G4double)nSub);
408 // set minimum value of the first sub-inteva 388 // set minimum value of the first sub-inteval
409 G4double alpha_i = alphaMin; << 389 G4double alpha_i = alphaMin;
410 for (G4int l = 0; l < nSub; ++l) { 390 for (G4int l = 0; l < nSub; ++l) {
411 for (G4int igl = 0; igl < 8; ++igl) { 391 for (G4int igl = 0; igl < 8; ++igl) {
412 // compute the emitted photon energy k 392 // compute the emitted photon energy k
413 const G4double k = G4Exp(alpha_i+gXGL[ 393 const G4double k = G4Exp(alpha_i+gXGL[igl]*delta)*fPrimaryTotalEnergy;
414 // compute the DCS value at k (without t 394 // compute the DCS value at k (without the constant, the 1/k, 1/F factors)
415 const G4double dcs = fIsLPMActive 395 const G4double dcs = fIsLPMActive
416 ? ComputeRelDXSectio 396 ? ComputeRelDXSectionPerAtom(k) // DCS WITHOUT LPM
417 : ComputeDXSectionPe 397 : ComputeDXSectionPerAtom(k); // DCS WITH LPM
418 // account Ter-Mikaelian suppression: ti 398 // account Ter-Mikaelian suppression: times 1/F with F = 1+(k_p/k)^2
419 xSection += gWGL[igl]*dcs/(1.0+fDensityC 399 xSection += gWGL[igl]*dcs/(1.0+fDensityCorr/(k*k));
420 } 400 }
421 // update sub-interval minimum value 401 // update sub-interval minimum value
422 alpha_i += delta; 402 alpha_i += delta;
423 } 403 }
424 // apply corrections due to variable transfo 404 // apply corrections due to variable transformation
425 xSection *= delta; 405 xSection *= delta;
426 // final check 406 // final check
427 return std::max(xSection, 0.); 407 return std::max(xSection, 0.);
428 } 408 }
429 409
430 // DCS WITH LPM EFFECT: complete screening apr 410 // DCS WITH LPM EFFECT: complete screening aprx. and includes LPM suppression
431 // ds/dk(Z,k) = C/[F*k]*{ Xi(s*F)*[y^2*G/4 +(1 411 // ds/dk(Z,k) = C/[F*k]*{ Xi(s*F)*[y^2*G/4 +(1-y+y^2/3)Phi]*[L_el-f_c+L_inel/Z]
432 // +(1-y)*[1+1/Z]/12} 412 // +(1-y)*[1+1/Z]/12} with C = 16\alpha r_0^2/3 Z^2 and
433 // Xi(s),G(s), Phi(s) are LPM suppression func 413 // Xi(s),G(s), Phi(s) are LPM suppression functions:
434 // 414 //
435 // LPM SUPPRESSION: The 's' is the suppression 415 // LPM SUPPRESSION: The 's' is the suppression variable and F = F(k,k_p) =
436 // 1+(k_p/k)^2 with k_p = hbar*w_p*E/(m*c^2) i 416 // 1+(k_p/k)^2 with k_p = hbar*w_p*E/(m*c^2) is a material (e- density)
437 // dependent constant. F accounts the Ter-Mika 417 // dependent constant. F accounts the Ter-Mikaelian suppression with a smooth
438 // transition in the emitted photon energy. Al 418 // transition in the emitted photon energy. Also, the LPM suppression functions
439 // goes to 0 when s goes to 0 and goes to 1 wh 419 // goes to 0 when s goes to 0 and goes to 1 when s is increasing (=1 at s=~2)
440 // So evaluating the LPM suppression functions 420 // So evaluating the LPM suppression functions at 'sF' instead of 's' ensures a
441 // smooth transition depending on the emitted 421 // smooth transition depending on the emitted photon energy 'k': LPM effect is
442 // smoothly turned off i.e. Xi(sF)=G(sF)=Phi(s 422 // smoothly turned off i.e. Xi(sF)=G(sF)=Phi(sF)=1 when k << k_p because F >> 1
443 // and sF ~ s when k >> k_p since F ~ 1 in tha 423 // and sF ~ s when k >> k_p since F ~ 1 in that case.
444 // HERE, ds/dk(Z,k)*[F*k/C] is computed since: 424 // HERE, ds/dk(Z,k)*[F*k/C] is computed since:
445 // (i) DCS ~ 1/k factor will disappear due 425 // (i) DCS ~ 1/k factor will disappear due to the variable transformation
446 // v(k)=ln(k/E_t) -> dk/dv=E_t*e^v=k -> 426 // v(k)=ln(k/E_t) -> dk/dv=E_t*e^v=k -> ds/dv= ds/dk*dk/dv=ds/dk*k so it
447 // would cnacell out the 1/k factor => 427 // would cnacell out the 1/k factor => 1/k don't included here
448 // (ii) the constant factor C and Z don't de 428 // (ii) the constant factor C and Z don't depend on 'k' => not included here
449 // (iii) the 1/F(k) factor is accounted in th 429 // (iii) the 1/F(k) factor is accounted in the callers: explicitly (cross sec-
450 // tion computation) or implicitly thro 430 // tion computation) or implicitly through further variable transformaton
451 // (in the final state sampling algorit 431 // (in the final state sampling algorithm)
452 // COMPLETE SCREENING: see more at the DCS wit 432 // COMPLETE SCREENING: see more at the DCS without LPM effect below.
453 G4double 433 G4double
454 G4eBremsstrahlungRelModel::ComputeRelDXSection 434 G4eBremsstrahlungRelModel::ComputeRelDXSectionPerAtom(G4double gammaEnergy)
455 { 435 {
456 G4double dxsec = 0.0; 436 G4double dxsec = 0.0;
457 if (gammaEnergy < 0.) { 437 if (gammaEnergy < 0.) {
458 return dxsec; 438 return dxsec;
459 } 439 }
460 const G4double y = gammaEnergy/fPrimaryT 440 const G4double y = gammaEnergy/fPrimaryTotalEnergy;
461 const G4double onemy = 1.-y; 441 const G4double onemy = 1.-y;
462 const G4double dum0 = 0.25*y*y; 442 const G4double dum0 = 0.25*y*y;
463 // evaluate LPM functions (combined with the 443 // evaluate LPM functions (combined with the Ter-Mikaelian effect)
464 G4double funcGS, funcPhiS, funcXiS; 444 G4double funcGS, funcPhiS, funcXiS;
465 ComputeLPMfunctions(funcXiS, funcGS, funcPhi 445 ComputeLPMfunctions(funcXiS, funcGS, funcPhiS, gammaEnergy);
466 const ElementData* elDat = (*fElementData)[f << 446 const ElementData* elDat = gElementData[fCurrentIZ];
467 const G4double term1 = funcXiS*(dum0*fun 447 const G4double term1 = funcXiS*(dum0*funcGS+(onemy+2.0*dum0)*funcPhiS);
468 dxsec = term1*elDat->fZFactor1+onemy*elDat-> 448 dxsec = term1*elDat->fZFactor1+onemy*elDat->fZFactor2;
469 // 449 //
470 if (fIsScatOffElectron) { 450 if (fIsScatOffElectron) {
471 fSumTerm = dxsec; 451 fSumTerm = dxsec;
472 fNucTerm = term1*elDat->fZFactor11 + onemy 452 fNucTerm = term1*elDat->fZFactor11 + onemy/12.;
473 } 453 }
474 return std::max(dxsec,0.0); 454 return std::max(dxsec,0.0);
475 } 455 }
476 456
477 // DCS WITHOUT LPM EFFECT: DCS with sceening ( 457 // DCS WITHOUT LPM EFFECT: DCS with sceening (Z>5) and Coulomb cor. no LPM
478 // ds/dk(Z,k)=C/[F*k]*{(1-y+3*y^2/4)*[(0.25*ph 458 // ds/dk(Z,k)=C/[F*k]*{(1-y+3*y^2/4)*[(0.25*phi1(g)-ln(Z)/3-f_c)+(0.25*psi1(e)
479 // -2*ln(Z)/3)/Z]+ (1-y)*[(phi1(g)-phi2(g))+(p 459 // -2*ln(Z)/3)/Z]+ (1-y)*[(phi1(g)-phi2(g))+(psi1(e)-psi2(e))/Z]/8}
480 // where f_c(Z) is the Coulomb correction fact 460 // where f_c(Z) is the Coulomb correction factor and phi1(g),phi2(g) and psi1(e),
481 // psi2(e) are coherent and incoherent screeni 461 // psi2(e) are coherent and incoherent screening functions. In the Thomas-Fermi
482 // model of the atom, the screening functions 462 // model of the atom, the screening functions will have a form that do not
483 // depend on Z (not explicitly). These numeric 463 // depend on Z (not explicitly). These numerical screening functions can be
484 // approximated as Tsai Eqs. [3.38-3.41] with 464 // approximated as Tsai Eqs. [3.38-3.41] with the variables g=gamma and
485 // e=epsilon given by Tsai Eqs. [3.30 and 3.31 465 // e=epsilon given by Tsai Eqs. [3.30 and 3.31] (see more details at the method
486 // ComputeScreeningFunctions()). Note, that in 466 // ComputeScreeningFunctions()). Note, that in case of complete screening i.e.
487 // g = e = 0 => 0.25*phi1(0)-ln(Z)/3 = ln(184. 467 // g = e = 0 => 0.25*phi1(0)-ln(Z)/3 = ln(184.149/Z^(1/3)) = L_el and
488 // 0.25*psi1(0)-2*ln(Z)/3=ln(1193.923/Z^(2/3)) 468 // 0.25*psi1(0)-2*ln(Z)/3=ln(1193.923/Z^(2/3))=L_inel and phi1(0)-phi2(0) =
489 // psi1(0)-psi2(0) = 2/3 so the DCS in complet 469 // psi1(0)-psi2(0) = 2/3 so the DCS in complete screening =>
490 // COMPLETE SCREENING: 470 // COMPLETE SCREENING:
491 // ds/dk(Z,k)=C/k*{(1-y+3*y^2/4)*[L_el-f_c+L_i 471 // ds/dk(Z,k)=C/k*{(1-y+3*y^2/4)*[L_el-f_c+L_inel/Z] + (1-y)*[1+1/Z]/12} that is
492 // used in case of DCS with LPM above (if all 472 // used in case of DCS with LPM above (if all the suprression functions are
493 // absent i.e. their value = 1). 473 // absent i.e. their value = 1).
494 // Since the Thomas-Fermi model of the atom is 474 // Since the Thomas-Fermi model of the atom is not accurate at low Z, the DCS in
495 // complete screening is used here at low Z(<5 475 // complete screening is used here at low Z(<5) with L_el(Z), L_inel(Z) values
496 // computed by using the Dirac-Fock model of t 476 // computed by using the Dirac-Fock model of the atom.
497 // NOTE: that the Ter-Mikaelian suppression is 477 // NOTE: that the Ter-Mikaelian suppression is accounted in the DCS through the
498 // 1/F factor but it is included in the caller 478 // 1/F factor but it is included in the caller and not considered here.
499 // HERE, ds/dk(Z,k)*[F*k/C] is computed exactl 479 // HERE, ds/dk(Z,k)*[F*k/C] is computed exactly like in the DCS with LPM case.
500 G4double 480 G4double
501 G4eBremsstrahlungRelModel::ComputeDXSectionPer 481 G4eBremsstrahlungRelModel::ComputeDXSectionPerAtom(G4double gammaEnergy)
502 { 482 {
503 G4double dxsec = 0.0; 483 G4double dxsec = 0.0;
504 if (gammaEnergy < 0.) { 484 if (gammaEnergy < 0.) {
505 return dxsec; 485 return dxsec;
506 } 486 }
507 const G4double y = gammaEnergy/fPrim 487 const G4double y = gammaEnergy/fPrimaryTotalEnergy;
508 const G4double onemy = 1.-y; 488 const G4double onemy = 1.-y;
509 const G4double dum0 = onemy+0.75*y*y; 489 const G4double dum0 = onemy+0.75*y*y;
510 const ElementData* elDat = (*fElementData)[f << 490 const ElementData* elDat = gElementData[fCurrentIZ];
511 // use complete screening and L_el, L_inel f 491 // use complete screening and L_el, L_inel from Dirac-Fock model instead of TF
512 if (fCurrentIZ < 5 || fIsUseCompleteScreenin 492 if (fCurrentIZ < 5 || fIsUseCompleteScreening) {
513 dxsec = dum0*elDat->fZFactor1; 493 dxsec = dum0*elDat->fZFactor1;
514 dxsec += onemy*elDat->fZFactor2; 494 dxsec += onemy*elDat->fZFactor2;
515 if (fIsScatOffElectron) { 495 if (fIsScatOffElectron) {
516 fSumTerm = dxsec; 496 fSumTerm = dxsec;
517 fNucTerm = dum0*elDat->fZFactor11+onemy/ 497 fNucTerm = dum0*elDat->fZFactor11+onemy/12.;
518 } 498 }
519 } else { 499 } else {
520 // use Tsai's analytical approx. (Tsai Eqs 500 // use Tsai's analytical approx. (Tsai Eqs. [3.38-3.41]) to the 'universal'
521 // numerical screening functions computed 501 // numerical screening functions computed by using the TF model of the atom
522 const G4double invZ = 1./(G4double)fCur 502 const G4double invZ = 1./(G4double)fCurrentIZ;
523 const G4double Fz = elDat->fFz; 503 const G4double Fz = elDat->fFz;
524 const G4double logZ = elDat->fLogZ; 504 const G4double logZ = elDat->fLogZ;
525 const G4double dum1 = y/(fPrimaryTotalE 505 const G4double dum1 = y/(fPrimaryTotalEnergy-gammaEnergy);
526 const G4double gamma = dum1*elDat->fGamm 506 const G4double gamma = dum1*elDat->fGammaFactor;
527 const G4double epsilon = dum1*elDat->fEpsi 507 const G4double epsilon = dum1*elDat->fEpsilonFactor;
528 // evaluate the screening functions 508 // evaluate the screening functions
529 G4double phi1, phi1m2, psi1, psi1m2; 509 G4double phi1, phi1m2, psi1, psi1m2;
530 ComputeScreeningFunctions(phi1, phi1m2, ps 510 ComputeScreeningFunctions(phi1, phi1m2, psi1, psi1m2, gamma, epsilon);
531 dxsec = dum0*((0.25*phi1-Fz) + (0.25*psi1 511 dxsec = dum0*((0.25*phi1-Fz) + (0.25*psi1-2.*logZ/3.)*invZ);
532 dxsec += 0.125*onemy*(phi1m2 + psi1m2*invZ 512 dxsec += 0.125*onemy*(phi1m2 + psi1m2*invZ);
533 if (fIsScatOffElectron) { 513 if (fIsScatOffElectron) {
534 fSumTerm = dxsec; 514 fSumTerm = dxsec;
535 fNucTerm = dum0*(0.25*phi1-Fz) + 0.125*o 515 fNucTerm = dum0*(0.25*phi1-Fz) + 0.125*onemy*phi1m2;
536 } 516 }
537 } 517 }
538 return std::max(dxsec,0.0); 518 return std::max(dxsec,0.0);
539 } 519 }
540 520
541 // Coherent and incoherent screening function 521 // Coherent and incoherent screening function approximations (see Tsai
542 // Eqs.[3.38-3.41]). Tsai's analytical approxi 522 // Eqs.[3.38-3.41]). Tsai's analytical approximations to the numerical screening
543 // functions computed by using the Thomas-Ferm 523 // functions computed by using the Thomas-Fermi model of atom (Moliere's appro-
544 // ximation to the numerical TF screening func 524 // ximation to the numerical TF screening function). In the TF-model, these
545 // screening functions can be expressed in a ' 525 // screening functions can be expressed in a 'universal' i.e. Z (directly) inde-
546 // pendent variable (see Tsai Eqs. Eqs. [3.30 526 // pendent variable (see Tsai Eqs. Eqs. [3.30 and 3.31]).
547 void G4eBremsstrahlungRelModel::ComputeScreeni 527 void G4eBremsstrahlungRelModel::ComputeScreeningFunctions(G4double& phi1,
548 528 G4double& phi1m2,
549 529 G4double& psi1,
550 530 G4double& psi1m2,
551 531 const G4double gam,
552 532 const G4double eps)
553 { 533 {
554 const G4double gam2 = gam*gam; 534 const G4double gam2 = gam*gam;
555 phi1 = 16.863-2.0*G4Log(1.0+0.311877*gam2) 535 phi1 = 16.863-2.0*G4Log(1.0+0.311877*gam2)+2.4*G4Exp(-0.9*gam)
556 +1.6*G4Exp(-1.5*gam); 536 +1.6*G4Exp(-1.5*gam);
557 phi1m2 = 2.0/(3.0+19.5*gam+18.0*gam2); // 537 phi1m2 = 2.0/(3.0+19.5*gam+18.0*gam2); // phi1-phi2
558 const G4double eps2 = eps*eps; 538 const G4double eps2 = eps*eps;
559 psi1 = 24.34-2.0*G4Log(1.0+13.111641*eps2) 539 psi1 = 24.34-2.0*G4Log(1.0+13.111641*eps2)+2.8*G4Exp(-8.0*eps)
560 +1.2*G4Exp(-29.2*eps); 540 +1.2*G4Exp(-29.2*eps);
561 psi1m2 = 2.0/(3.0+120.0*eps+1200.0*eps2); // 541 psi1m2 = 2.0/(3.0+120.0*eps+1200.0*eps2); //psi1-psi2
562 } 542 }
563 543
564 void 544 void
565 G4eBremsstrahlungRelModel::SampleSecondaries(s 545 G4eBremsstrahlungRelModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
566 c 546 const G4MaterialCutsCouple* couple,
567 c 547 const G4DynamicParticle* dp,
568 G 548 G4double cutEnergy,
569 G 549 G4double maxEnergy)
570 { 550 {
571 const G4double kineticEnergy = dp->GetKin 551 const G4double kineticEnergy = dp->GetKineticEnergy();
>> 552 // const G4double logKineticEnergy = dp->GetLogKineticEnergy();
572 if (kineticEnergy < LowEnergyLimit()) { 553 if (kineticEnergy < LowEnergyLimit()) {
573 return; 554 return;
574 } 555 }
575 // min, max kinetic energy limits 556 // min, max kinetic energy limits
576 const G4double tmin = std::min(cutEnergy, ki 557 const G4double tmin = std::min(cutEnergy, kineticEnergy);
577 const G4double tmax = std::min(maxEnergy, ki 558 const G4double tmax = std::min(maxEnergy, kineticEnergy);
578 if (tmin >= tmax) { 559 if (tmin >= tmax) {
579 return; 560 return;
580 } 561 }
581 // 562 //
582 SetupForMaterial(fPrimaryParticle, couple->G 563 SetupForMaterial(fPrimaryParticle, couple->GetMaterial(), kineticEnergy);
583 const G4Element* elm = SelectTargetAtom(coup 564 const G4Element* elm = SelectTargetAtom(couple,fPrimaryParticle,kineticEnergy,
584 dp-> 565 dp->GetLogKineticEnergy(),tmin,tmax);
585 // 566 //
586 fCurrentIZ = elm->GetZasInt(); 567 fCurrentIZ = elm->GetZasInt();
587 const ElementData* elDat = (*fElementData)[f << 568 const ElementData* elDat = gElementData[fCurrentIZ];
588 const G4double funcMax = elDat->fZFactor1+el 569 const G4double funcMax = elDat->fZFactor1+elDat->fZFactor2;
589 // get the random engine 570 // get the random engine
590 G4double rndm[2]; 571 G4double rndm[2];
591 CLHEP::HepRandomEngine* rndmEngine = G4Rando 572 CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
592 // min max of the transformed variable: x(k) 573 // min max of the transformed variable: x(k) = ln(k^2+k_p^2) that is in [ln(k_c^2+k_p^2), ln(E_k^2+k_p^2)]
593 const G4double xmin = G4Log(tmin*tmin+fDen 574 const G4double xmin = G4Log(tmin*tmin+fDensityCorr);
594 const G4double xrange = G4Log(tmax*tmax+fDen 575 const G4double xrange = G4Log(tmax*tmax+fDensityCorr)-xmin;
595 G4double gammaEnergy, funcVal; 576 G4double gammaEnergy, funcVal;
596 do { 577 do {
597 rndmEngine->flatArray(2, rndm); 578 rndmEngine->flatArray(2, rndm);
598 gammaEnergy = std::sqrt(std::max(G4Exp(xmi 579 gammaEnergy = std::sqrt(std::max(G4Exp(xmin+rndm[0]*xrange)-fDensityCorr, 0.0));
599 funcVal = fIsLPMActive 580 funcVal = fIsLPMActive
600 ? ComputeRelDXSectionPerAtom( 581 ? ComputeRelDXSectionPerAtom(gammaEnergy)
601 : ComputeDXSectionPerAtom(gam 582 : ComputeDXSectionPerAtom(gammaEnergy);
602 // cross-check of proper function maximum 583 // cross-check of proper function maximum in the rejection
603 // if (funcVal > funcMax) { 584 // if (funcVal > funcMax) {
604 // G4cout << "### G4eBremsstrahlungRelMod 585 // G4cout << "### G4eBremsstrahlungRelModel Warning: Majoranta exceeded! "
605 // << funcVal << " > " << funcMax 586 // << funcVal << " > " << funcMax
606 // << " Egamma(MeV)= " << gammaEnergy 587 // << " Egamma(MeV)= " << gammaEnergy
607 // << " Ee(MeV)= " << kineticEnergy 588 // << " Ee(MeV)= " << kineticEnergy
608 // << " " << GetName() 589 // << " " << GetName()
609 // << G4endl; 590 // << G4endl;
610 // } 591 // }
611 // Loop checking, 03-Aug-2015, Vladimir Iv 592 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
612 } while (funcVal < funcMax*rndm[1]); 593 } while (funcVal < funcMax*rndm[1]);
613 // 594 //
614 // scattering off nucleus or off e- by tripl 595 // scattering off nucleus or off e- by triplet model
615 if (fIsScatOffElectron && rndmEngine->flat() 596 if (fIsScatOffElectron && rndmEngine->flat()*fSumTerm>fNucTerm) {
616 GetTripletModel()->SampleSecondaries(vdp, 597 GetTripletModel()->SampleSecondaries(vdp, couple, dp, cutEnergy, maxEnergy);
617 return; 598 return;
618 } 599 }
619 // 600 //
620 // angles of the emitted gamma. ( Z - axis a 601 // angles of the emitted gamma. ( Z - axis along the parent particle)
621 // use general interface 602 // use general interface
622 G4ThreeVector gamDir = 603 G4ThreeVector gamDir =
623 GetAngularDistribution()->SampleDirection( 604 GetAngularDistribution()->SampleDirection(dp,fPrimaryTotalEnergy-gammaEnergy,
624 605 fCurrentIZ, couple->GetMaterial());
625 // create G4DynamicParticle object for the G 606 // create G4DynamicParticle object for the Gamma
626 auto gamma = new G4DynamicParticle(fGammaPar << 607 G4DynamicParticle* gamma = new G4DynamicParticle(fGammaParticle, gamDir,
>> 608 gammaEnergy);
627 vdp->push_back(gamma); 609 vdp->push_back(gamma);
628 // compute post-interaction kinematics of pr 610 // compute post-interaction kinematics of primary e-/e+ based on
629 // energy-momentum conservation 611 // energy-momentum conservation
630 const G4double totMomentum = std::sqrt(kinet 612 const G4double totMomentum = std::sqrt(kineticEnergy*(
631 fPrimaryTotalEn 613 fPrimaryTotalEnergy + CLHEP::electron_mass_c2));
632 G4ThreeVector dir = 614 G4ThreeVector dir =
633 (totMomentum*dp->GetMomentumDirec 615 (totMomentum*dp->GetMomentumDirection()-gammaEnergy*gamDir).unit();
634 const G4double finalE = kineticEnergy-gamm 616 const G4double finalE = kineticEnergy-gammaEnergy;
635 // if secondary gamma energy is higher than 617 // if secondary gamma energy is higher than threshold(very high by default)
636 // then stop tracking the primary particle a 618 // then stop tracking the primary particle and create new secondary e-/e+
637 // instead of the primary one 619 // instead of the primary one
638 if (gammaEnergy > SecondaryThreshold()) { 620 if (gammaEnergy > SecondaryThreshold()) {
639 fParticleChange->ProposeTrackStatus(fStopA 621 fParticleChange->ProposeTrackStatus(fStopAndKill);
640 fParticleChange->SetProposedKineticEnergy( 622 fParticleChange->SetProposedKineticEnergy(0.0);
641 auto el = new G4DynamicParticle( << 623 G4DynamicParticle* el = new G4DynamicParticle(
642 const_cast<G4ParticleDefinition* 624 const_cast<G4ParticleDefinition*>(fPrimaryParticle), dir, finalE);
643 vdp->push_back(el); 625 vdp->push_back(el);
644 } else { // continue tracking the primary e- 626 } else { // continue tracking the primary e-/e+ otherwise
645 fParticleChange->SetProposedMomentumDirect 627 fParticleChange->SetProposedMomentumDirection(dir);
646 fParticleChange->SetProposedKineticEnergy( 628 fParticleChange->SetProposedKineticEnergy(finalE);
647 } 629 }
648 } 630 }
649 631
650 void G4eBremsstrahlungRelModel::InitialiseElem 632 void G4eBremsstrahlungRelModel::InitialiseElementData()
651 { 633 {
>> 634 const G4int size = gElementData.size();
>> 635 if (size < gMaxZet+1) {
>> 636 gElementData.resize(gMaxZet+1, nullptr);
>> 637 }
652 // create for all elements that are in the d 638 // create for all elements that are in the detector
653 auto elemTable = G4Element::GetElementTable( << 639 const G4ElementTable* elemTable = G4Element::GetElementTable();
654 for (auto const & elem : *elemTable) { << 640 size_t numElems = (*elemTable).size();
655 const G4double zet = elem->GetZ(); << 641 for (size_t ielem=0; ielem<numElems; ++ielem) {
656 const G4int izet = std::min(elem->GetZasIn << 642 const G4Element* elem = (*elemTable)[ielem];
657 if (nullptr == (*fElementData)[izet]) { << 643 const G4double zet = elem->GetZ();
658 auto elemData = new ElementData(); << 644 const G4int izet = std::min(G4lrint(zet),gMaxZet);
>> 645 if (!gElementData[izet]) {
>> 646 ElementData *elemData = new ElementData();
659 const G4double fc = elem->GetfCoulomb(); 647 const G4double fc = elem->GetfCoulomb();
660 G4double Fel = 1.; 648 G4double Fel = 1.;
661 G4double Finel = 1.; 649 G4double Finel = 1.;
662 elemData->fLogZ = G4Log(zet); 650 elemData->fLogZ = G4Log(zet);
663 elemData->fFz = elemData->fLogZ/3.+f 651 elemData->fFz = elemData->fLogZ/3.+fc;
664 if (izet < 5) { 652 if (izet < 5) {
665 Fel = gFelLowZet[izet]; 653 Fel = gFelLowZet[izet];
666 Finel = gFinelLowZet[izet]; 654 Finel = gFinelLowZet[izet];
667 } else { 655 } else {
668 Fel = G4Log(184.15) - elemData->f 656 Fel = G4Log(184.15) - elemData->fLogZ/3.;
669 Finel = G4Log(1194) - 2.*elemData->f 657 Finel = G4Log(1194) - 2.*elemData->fLogZ/3.;
670 } 658 }
671 const G4double z13 = G4Pow::GetInstance( << 659 const G4double z23 = std::pow(zet,2./3.);
672 const G4double z23 = z13*z13; << 660 const G4double z13 = std::pow(zet,1./3.);
673 elemData->fZFactor1 = (Fel-fc)+Fine 661 elemData->fZFactor1 = (Fel-fc)+Finel/zet;
674 elemData->fZFactor11 = (Fel-fc); // 662 elemData->fZFactor11 = (Fel-fc); // used only for the triplet
675 elemData->fZFactor2 = (1.+1./zet)/1 663 elemData->fZFactor2 = (1.+1./zet)/12.;
676 elemData->fVarS1 = z23/(184.15*1 664 elemData->fVarS1 = z23/(184.15*184.15);
677 elemData->fILVarS1Cond = 1./(G4Log(std 665 elemData->fILVarS1Cond = 1./(G4Log(std::sqrt(2.0)*elemData->fVarS1));
678 elemData->fILVarS1 = 1./G4Log(elem 666 elemData->fILVarS1 = 1./G4Log(elemData->fVarS1);
679 elemData->fGammaFactor = 100.0*electro 667 elemData->fGammaFactor = 100.0*electron_mass_c2/z13;
680 elemData->fEpsilonFactor = 100.0*electro 668 elemData->fEpsilonFactor = 100.0*electron_mass_c2/z23;
681 (*fElementData)[izet] = elemData; << 669 gElementData[izet] = elemData;
682 } 670 }
683 } 671 }
684 } 672 }
685 673
686 void G4eBremsstrahlungRelModel::ComputeLPMfunc 674 void G4eBremsstrahlungRelModel::ComputeLPMfunctions(G4double& funcXiS,
687 675 G4double& funcGS,
688 676 G4double& funcPhiS,
689 677 const G4double egamma)
690 { 678 {
691 static const G4double sqrt2 = std::sqrt(2.); 679 static const G4double sqrt2 = std::sqrt(2.);
692 const G4double redegamma = egamma/fPrimar 680 const G4double redegamma = egamma/fPrimaryTotalEnergy;
693 const G4double varSprime = std::sqrt(0.12 681 const G4double varSprime = std::sqrt(0.125*redegamma*fLPMEnergy/
694 ((1.0-redegamm 682 ((1.0-redegamma)*fPrimaryTotalEnergy));
695 const ElementData* elDat = (*fElementData << 683 const ElementData* elDat = gElementData[fCurrentIZ];
696 const G4double varS1 = elDat->fVarS1; 684 const G4double varS1 = elDat->fVarS1;
697 const G4double condition = sqrt2*varS1; 685 const G4double condition = sqrt2*varS1;
698 G4double funcXiSprime = 2.0; 686 G4double funcXiSprime = 2.0;
699 if (varSprime > 1.0) { 687 if (varSprime > 1.0) {
700 funcXiSprime = 1.0; 688 funcXiSprime = 1.0;
701 } else if (varSprime > condition) { 689 } else if (varSprime > condition) {
702 const G4double ilVarS1Cond = elDat->fILVar 690 const G4double ilVarS1Cond = elDat->fILVarS1Cond;
703 const G4double funcHSprime = G4Log(varSpri 691 const G4double funcHSprime = G4Log(varSprime)*ilVarS1Cond;
704 funcXiSprime = 1.0 + funcHSprime - 0.08*(1 692 funcXiSprime = 1.0 + funcHSprime - 0.08*(1.0-funcHSprime)*funcHSprime
705 *(2.0-fu 693 *(2.0-funcHSprime)*ilVarS1Cond;
706 } 694 }
707 const G4double varS = varSprime/std::sqrt 695 const G4double varS = varSprime/std::sqrt(funcXiSprime);
708 // - include dielectric suppression effect i 696 // - include dielectric suppression effect into s according to Migdal
709 const G4double varShat = varS*(1.0+fDensityC 697 const G4double varShat = varS*(1.0+fDensityCorr/(egamma*egamma));
710 funcXiS = 2.0; 698 funcXiS = 2.0;
711 if (varShat > 1.0) { 699 if (varShat > 1.0) {
712 funcXiS = 1.0; 700 funcXiS = 1.0;
713 } else if (varShat > varS1) { 701 } else if (varShat > varS1) {
714 funcXiS = 1.0+G4Log(varShat)*elDat->fILVar 702 funcXiS = 1.0+G4Log(varShat)*elDat->fILVarS1;
715 } 703 }
716 GetLPMFunctions(funcGS, funcPhiS, varShat); 704 GetLPMFunctions(funcGS, funcPhiS, varShat);
717 //ComputeLPMGsPhis(funcGS, funcPhiS, varShat 705 //ComputeLPMGsPhis(funcGS, funcPhiS, varShat);
718 // 706 //
719 //MAKE SURE SUPPRESSION IS SMALLER THAN 1: d 707 //MAKE SURE SUPPRESSION IS SMALLER THAN 1: due to Migdal's approximation on xi
720 if (funcXiS*funcPhiS > 1. || varShat > 0.57) 708 if (funcXiS*funcPhiS > 1. || varShat > 0.57) {
721 funcXiS=1./funcPhiS; 709 funcXiS=1./funcPhiS;
722 } 710 }
723 } 711 }
724 712
725 void G4eBremsstrahlungRelModel::ComputeLPMGsPh 713 void G4eBremsstrahlungRelModel::ComputeLPMGsPhis(G4double& funcGS,
726 714 G4double& funcPhiS,
727 715 const G4double varShat)
728 { 716 {
729 if (varShat < 0.01) { 717 if (varShat < 0.01) {
730 funcPhiS = 6.0*varShat*(1.0-CLHEP::pi*varS 718 funcPhiS = 6.0*varShat*(1.0-CLHEP::pi*varShat);
731 funcGS = 12.0*varShat-2.0*funcPhiS; 719 funcGS = 12.0*varShat-2.0*funcPhiS;
732 } else { 720 } else {
733 const G4double varShat2 = varShat*varShat; 721 const G4double varShat2 = varShat*varShat;
734 const G4double varShat3 = varShat*varShat2 722 const G4double varShat3 = varShat*varShat2;
735 const G4double varShat4 = varShat2*varShat 723 const G4double varShat4 = varShat2*varShat2;
736 // use Stanev approximation: for \psi(s) a 724 // use Stanev approximation: for \psi(s) and compute G(s)
737 if (varShat < 0.415827) { 725 if (varShat < 0.415827) {
738 funcPhiS = 1.0-G4Exp(-6.0*varShat*(1.0+v 726 funcPhiS = 1.0-G4Exp(-6.0*varShat*(1.0+varShat*(3.0-CLHEP::pi))
739 + varShat3/(0.623+0.796*varSha 727 + varShat3/(0.623+0.796*varShat+0.658*varShat2));
740 // 1-\exp \left\{-4s-\frac{8s^2}{1+3.936 728 // 1-\exp \left\{-4s-\frac{8s^2}{1+3.936s+4.97s^2-0.05s^3+7.5s^4} \right\}
741 const G4double funcPsiS = 1.0 - G4Exp(-4 729 const G4double funcPsiS = 1.0 - G4Exp(-4.0*varShat
742 - 8.0*varShat2/ 730 - 8.0*varShat2/(1.0+3.936*varShat+4.97*varShat2
743 - 0.05*varShat3 731 - 0.05*varShat3 + 7.5*varShat4));
744 // G(s) = 3 \psi(s) - 2 \phi(s) 732 // G(s) = 3 \psi(s) - 2 \phi(s)
745 funcGS = 3.0*funcPsiS - 2.0*funcPhiS; 733 funcGS = 3.0*funcPsiS - 2.0*funcPhiS;
746 } else if (varShat<1.55) { 734 } else if (varShat<1.55) {
747 funcPhiS = 1.0-G4Exp(-6.0*varShat*(1.0+v 735 funcPhiS = 1.0-G4Exp(-6.0*varShat*(1.0+varShat*(3.0-CLHEP::pi))
748 + varShat3/(0.623+0.796*varSha 736 + varShat3/(0.623+0.796*varShat+0.658*varShat2));
749 const G4double dum0 = -0.160723 737 const G4double dum0 = -0.160723 + 3.755030*varShat
750 -1.798138*varShat 738 -1.798138*varShat2 + 0.672827*varShat3
751 -0.120772*varShat 739 -0.120772*varShat4;
752 funcGS = std::tanh(dum0); 740 funcGS = std::tanh(dum0);
753 } else { 741 } else {
754 funcPhiS = 1.0-0.011905/varShat4; 742 funcPhiS = 1.0-0.011905/varShat4;
755 if (varShat<1.9156) { 743 if (varShat<1.9156) {
756 const G4double dum0 = -0.160723 744 const G4double dum0 = -0.160723 + 3.755030*varShat
757 -1.798138*varSha 745 -1.798138*varShat2 + 0.672827*varShat3
758 -0.120772*varSha 746 -0.120772*varShat4;
759 funcGS = std::tanh(dum0); 747 funcGS = std::tanh(dum0);
760 } else { 748 } else {
761 funcGS = 1.0-0.023065/varShat4; 749 funcGS = 1.0-0.023065/varShat4;
762 } 750 }
763 } 751 }
764 } 752 }
765 } 753 }
766 754
767 // s goes up to 2 with ds = 0.01 to be the def 755 // s goes up to 2 with ds = 0.01 to be the default bining
768 void G4eBremsstrahlungRelModel::InitLPMFunctio 756 void G4eBremsstrahlungRelModel::InitLPMFunctions()
769 { 757 {
770 if (!fLPMFuncs->fIsInitialized) { << 758 if (!gLPMFuncs.fIsInitialized) {
771 const G4int num = fLPMFuncs->fSLimit*fLPMF << 759 const G4int num = gLPMFuncs.fSLimit*gLPMFuncs.fISDelta+1;
772 fLPMFuncs->fLPMFuncG.resize(num); << 760 gLPMFuncs.fLPMFuncG.resize(num);
773 fLPMFuncs->fLPMFuncPhi.resize(num); << 761 gLPMFuncs.fLPMFuncPhi.resize(num);
774 for (G4int i = 0; i < num; ++i) { 762 for (G4int i = 0; i < num; ++i) {
775 const G4double sval=i/fLPMFuncs->fISDelt << 763 const G4double sval=i/gLPMFuncs.fISDelta;
776 ComputeLPMGsPhis(fLPMFuncs->fLPMFuncG[i] << 764 ComputeLPMGsPhis(gLPMFuncs.fLPMFuncG[i],gLPMFuncs.fLPMFuncPhi[i],sval);
777 } 765 }
778 fLPMFuncs->fIsInitialized = true; << 766 gLPMFuncs.fIsInitialized = true;
779 } 767 }
780 } 768 }
781 769
782 void G4eBremsstrahlungRelModel::GetLPMFunction 770 void G4eBremsstrahlungRelModel::GetLPMFunctions(G4double& lpmGs,
783 771 G4double& lpmPhis,
784 772 const G4double sval)
785 { 773 {
786 if (sval < fLPMFuncs->fSLimit) { << 774 if (sval < gLPMFuncs.fSLimit) {
787 G4double val = sval*fLPMFuncs->fISDelt << 775 G4double val = sval*gLPMFuncs.fISDelta;
788 const G4int ilow = (G4int)val; 776 const G4int ilow = (G4int)val;
789 val -= ilow; 777 val -= ilow;
790 lpmGs = (fLPMFuncs->fLPMFuncG[ilow+1]-fL << 778 lpmGs = (gLPMFuncs.fLPMFuncG[ilow+1]-gLPMFuncs.fLPMFuncG[ilow])*val
791 + fLPMFuncs->fLPMFuncG[ilow]; << 779 + gLPMFuncs.fLPMFuncG[ilow];
792 lpmPhis = (fLPMFuncs->fLPMFuncPhi[ilow+1]- << 780 lpmPhis = (gLPMFuncs.fLPMFuncPhi[ilow+1]-gLPMFuncs.fLPMFuncPhi[ilow])*val
793 + fLPMFuncs->fLPMFuncPhi[ilow]; << 781 + gLPMFuncs.fLPMFuncPhi[ilow];
794 } else { 782 } else {
795 G4double ss = sval*sval; 783 G4double ss = sval*sval;
796 ss *= ss; 784 ss *= ss;
797 lpmPhis = 1.0-0.01190476/ss; 785 lpmPhis = 1.0-0.01190476/ss;
798 lpmGs = 1.0-0.0230655/ss; 786 lpmGs = 1.0-0.0230655/ss;
799 } 787 }
800 } 788 }
801 789
802 790