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