| 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: G4SeltzerBergerModel.cc 98737 2016-08-09 12:51:38Z gcosmo $
>> 27 //
26 // ------------------------------------------- 28 // -------------------------------------------------------------------
27 // 29 //
28 // GEANT4 Class file 30 // GEANT4 Class file
29 // 31 //
30 // 32 //
31 // File name: G4SeltzerBergerModel 33 // File name: G4SeltzerBergerModel
32 // 34 //
33 // Author: Vladimir Ivanchenko use inhe 35 // Author: Vladimir Ivanchenko use inheritance from Andreas Schaelicke
34 // base class implementing ultr 36 // base class implementing ultra relativistic bremsstrahlung
35 // model << 37 // model
36 // 38 //
37 // Creation date: 04.10.2011 39 // Creation date: 04.10.2011
38 // 40 //
39 // Modifications: 41 // Modifications:
40 // 42 //
41 // 24.07.2018 Introduced possibility to use sa <<
42 // emitted photon energy (instead o <<
43 // Seltzer-Berger scalled DCS for b <<
44 // Using these sampling tables opti <<
45 // state generation than the origin <<
46 // extra memory (+ ~6MB in the case <<
47 // (M Novak) <<
48 // <<
49 // ------------------------------------------- 43 // -------------------------------------------------------------------
50 // 44 //
>> 45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51 47
52 #include "G4SeltzerBergerModel.hh" 48 #include "G4SeltzerBergerModel.hh"
53 #include "G4PhysicalConstants.hh" 49 #include "G4PhysicalConstants.hh"
54 #include "G4SystemOfUnits.hh" 50 #include "G4SystemOfUnits.hh"
>> 51 #include "G4Electron.hh"
>> 52 #include "G4Positron.hh"
>> 53 #include "G4Gamma.hh"
55 #include "Randomize.hh" 54 #include "Randomize.hh"
56 #include "G4Material.hh" 55 #include "G4Material.hh"
57 #include "G4Element.hh" 56 #include "G4Element.hh"
58 #include "G4ElementVector.hh" 57 #include "G4ElementVector.hh"
>> 58 #include "G4ProductionCutsTable.hh"
59 #include "G4ParticleChangeForLoss.hh" 59 #include "G4ParticleChangeForLoss.hh"
60 #include "G4SBBremTable.hh" <<
61 #include "G4ModifiedTsai.hh" 60 #include "G4ModifiedTsai.hh"
62 61
63 #include "G4EmParameters.hh" <<
64 #include "G4ProductionCutsTable.hh" <<
65 #include "G4Gamma.hh" <<
66 #include "G4Electron.hh" <<
67 <<
68 #include "G4Physics2DVector.hh" 62 #include "G4Physics2DVector.hh"
69 #include "G4Exp.hh" 63 #include "G4Exp.hh"
70 #include "G4Log.hh" 64 #include "G4Log.hh"
71 #include "G4AutoLock.hh" <<
72 65
73 #include "G4ios.hh" 66 #include "G4ios.hh"
74 <<
75 #include <fstream> 67 #include <fstream>
76 #include <iomanip> 68 #include <iomanip>
77 #include <sstream> <<
78 #include <thread> <<
79 69
80 G4double G4SeltzerBergerModel::gYLimitData[] = << 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
81 G4Physics2DVector* G4SeltzerBergerModel::gSBDC <<
82 G4SBBremTable* G4SeltzerBergerModel::gSBSampli <<
83 <<
84 // constant DCS factor: 16\alpha r_0^2/3 <<
85 const G4double G4SeltzerBergerModel::gBremFact <<
86 = 16. * CLHEP::fine_structure_const * CLHEP: <<
87 * CLHEP::classic_electr_radius/3.; <<
88 <<
89 // Migdal's constant: 4\pi r_0*electron_reduce <<
90 const G4double G4SeltzerBergerModel::gMigdalCo <<
91 = 4. * CLHEP::pi * CLHEP::classic_electr_rad <<
92 * CLHEP::electron_Compton_length * CLHEP:: <<
93 <<
94 static constexpr G4double twoMass = 2* CLHEP:: <<
95 static constexpr G4double kAlpha = CLHEP::twop <<
96 static std::once_flag applyOnce; <<
97 71
98 namespace << 72 using namespace std;
99 { <<
100 G4Mutex theSBMutex = G4MUTEX_INITIALIZER; <<
101 73
102 // for numerical integration on [0,1] << 74 G4Physics2DVector* G4SeltzerBergerModel::dataSB[] = {nullptr};
103 const G4double gXGL[8] = { << 75 G4double G4SeltzerBergerModel::ylimit[] = {0.0};
104 1.98550718e-02, 1.01666761e-01, 2.37233795 << 76 G4double G4SeltzerBergerModel::expnumlim = -12.;
105 5.91717321e-01, 7.62766205e-01, 8.98333239 <<
106 }; <<
107 const G4double gWGL[8] = { <<
108 5.06142681e-02, 1.11190517e-01, 1.56853323 <<
109 1.81341892e-01, 1.56853323e-01, 1.11190517 <<
110 }; <<
111 } <<
112 77
113 G4SeltzerBergerModel::G4SeltzerBergerModel(con 78 G4SeltzerBergerModel::G4SeltzerBergerModel(const G4ParticleDefinition* p,
114 con 79 const G4String& nam)
115 : G4VEmModel(nam), << 80 : G4eBremsstrahlungRelModel(p,nam),useBicubicInterpolation(false)
116 fGammaParticle(G4Gamma::Gamma()), <<
117 fLowestKinEnergy(1.0*CLHEP::keV) <<
118 { 81 {
119 SetLowEnergyLimit(fLowestKinEnergy); << 82 SetLowestKinEnergy(1.0*keV);
120 SetAngularDistribution(new G4ModifiedTsai()) << 83 SetLowEnergyLimit(LowestKinEnergy());
121 if (fPrimaryParticle != p) { SetParticle(p); << 84 SetLPMFlag(false);
>> 85 nwarn = 0;
>> 86 idx = idy = 0;
122 } 87 }
123 88
>> 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 90
124 G4SeltzerBergerModel::~G4SeltzerBergerModel() 91 G4SeltzerBergerModel::~G4SeltzerBergerModel()
125 { 92 {
126 // delete SB-DCS data per Z << 93 if(IsMaster()) {
127 if (isInitializer) { << 94 for(size_t i=0; i<101; ++i) {
128 for (std::size_t iz = 0; iz < gMaxZet; ++i << 95 if(dataSB[i]) {
129 if (gSBDCSData[iz]) { << 96 delete dataSB[i];
130 delete gSBDCSData[iz]; << 97 dataSB[i] = nullptr;
131 gSBDCSData[iz] = nullptr; << 98 }
132 } <<
133 } <<
134 if (gSBSamplingTable) { <<
135 delete gSBSamplingTable; <<
136 gSBSamplingTable = nullptr; <<
137 } 99 }
138 } 100 }
139 } 101 }
140 102
>> 103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 104
141 void G4SeltzerBergerModel::Initialise(const G4 105 void G4SeltzerBergerModel::Initialise(const G4ParticleDefinition* p,
142 const G4 106 const G4DataVector& cuts)
143 { 107 {
144 // parameters in each thread << 108 // Access to elements
145 if (fPrimaryParticle != p) { << 109 if(IsMaster()) {
146 SetParticle(p); <<
147 } <<
148 fIsUseSamplingTables = G4EmParameters::Insta <<
149 fCurrentIZ = 0; <<
150 <<
151 // initialise static tables for the Seltzer- <<
152 std::call_once(applyOnce, [this]() { isIniti <<
153 <<
154 if (isInitializer) { <<
155 G4AutoLock l(&theSBMutex); <<
156 <<
157 // initialisation per element is done only <<
158 auto elemTable = G4Element::GetElementTabl <<
159 for (auto const & elm : *elemTable) { <<
160 G4int Z = std::max(1,std::min(elm->GetZa <<
161 // load SB-DCS data for this atomic numb <<
162 if (gSBDCSData[Z] == nullptr) ReadData(Z <<
163 } <<
164 110
165 // init sampling tables if it was requeste << 111 // check environment variable
166 if (fIsUseSamplingTables) { << 112 // Build the complete string identifying the file with the data set
167 if (nullptr == gSBSamplingTable) { << 113 char* path = getenv("G4LEDATA");
168 gSBSamplingTable = new G4SBBremTable() << 114
>> 115 const G4ElementTable* theElmTable = G4Element::GetElementTable();
>> 116 size_t numOfElm = G4Element::GetNumberOfElements();
>> 117 if(numOfElm > 0) {
>> 118 for(size_t i=0; i<numOfElm; ++i) {
>> 119 G4int Z = G4lrint(((*theElmTable)[i])->GetZ());
>> 120 if(Z < 1) { Z = 1; }
>> 121 else if(Z > 100) { Z = 100; }
>> 122 //G4cout << "Z= " << Z << G4endl;
>> 123 // Initialisation
>> 124 if(nullptr == dataSB[Z]) { ReadData(Z, path); }
169 } 125 }
170 gSBSamplingTable->Initialize(std::max(fL <<
171 HighEnergyL <<
172 } 126 }
173 l.unlock(); <<
174 } <<
175 // element selectors are initialized in the <<
176 if (IsMaster()) { <<
177 InitialiseElementSelectors(p, cuts); <<
178 } <<
179 // initialisation in all threads <<
180 if (nullptr == fParticleChange) { <<
181 fParticleChange = GetParticleChangeForLoss <<
182 } <<
183 auto trmodel = GetTripletModel(); <<
184 if (nullptr != trmodel) { <<
185 trmodel->Initialise(p, cuts); <<
186 fIsScatOffElectron = true; <<
187 } 127 }
188 } <<
189 128
190 void G4SeltzerBergerModel::InitialiseLocal(con << 129 G4eBremsstrahlungRelModel::Initialise(p, cuts);
191 <<
192 { <<
193 SetElementSelectors(masterModel->GetElementS <<
194 } 130 }
195 131
196 void G4SeltzerBergerModel::SetParticle(const G << 132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
197 { <<
198 fPrimaryParticle = p; <<
199 fIsElectron = (p == G4Electron::Electron()); <<
200 } <<
201 133
202 void G4SeltzerBergerModel::ReadData(G4int Z) { << 134 G4String G4SeltzerBergerModel::DirectoryPath() const
203 // return if it has been already loaded <<
204 if (gSBDCSData[Z] != nullptr) return; <<
205 <<
206 if (gSBDCSData[Z] == nullptr) { <<
207 std::ostringstream ost; <<
208 ost << G4EmParameters::Instance()->GetDirL <<
209 std::ifstream fin(ost.str().c_str()); <<
210 if (!fin.is_open()) { <<
211 G4ExceptionDescription ed; <<
212 ed << "Bremsstrahlung data file <" << os <<
213 << "> is not opened!"; <<
214 G4Exception("G4SeltzerBergerModel::ReadD <<
215 ed,"G4LEDATA version should be G4EMLOW6. <<
216 return; <<
217 } <<
218 //G4cout << "G4SeltzerBergerModel read fro <<
219 // << ">" << G4endl; <<
220 auto v = new G4Physics2DVector(); <<
221 if (v->Retrieve(fin)) { <<
222 v->SetBicubicInterpolation(fIsUseBicubic <<
223 static const G4double emaxlog = 4*G4Log( <<
224 gYLimitData[Z] = v->Value(0.97, emaxlog, <<
225 gSBDCSData[Z] = v; <<
226 } else { <<
227 G4ExceptionDescription ed; <<
228 ed << "Bremsstrahlung data file <" << os <<
229 << "> is not retrieved!"; <<
230 G4Exception("G4SeltzerBergerModel::ReadD <<
231 ed,"G4LEDATA version should be G4EMLOW6. <<
232 delete v; <<
233 } <<
234 } <<
235 } <<
236 <<
237 // minimum primary (e-/e+) energy at which dis <<
238 G4double G4SeltzerBergerModel::MinPrimaryEnerg <<
239 <<
240 <<
241 { 135 {
242 return std::max(fLowestKinEnergy, cut); << 136 return "/brem_SB/br";
243 } 137 }
244 138
245 // Sets kinematical variables like E_kin, E_t << 139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
246 // for characteristic photon energy k_p (more <<
247 // k_p^2) for the Ter-Mikaelian suppression ef <<
248 void G4SeltzerBergerModel::SetupForMaterial(co <<
249 co <<
250 G4double <<
251 { <<
252 fDensityFactor = gMigdalConstant*mat->GetEle <<
253 // calculate threshold for density effect: k <<
254 fPrimaryKinEnergy = kinEnergy; <<
255 fPrimaryTotalEnergy = kinEnergy + CLHEP::ele <<
256 fDensityCorr = fDensityFactor*fPrimaryTotalE <<
257 } <<
258 140
259 // Computes the restricted dE/dx as the approp << 141 void G4SeltzerBergerModel::ReadData(G4int Z, const char* path)
260 // element contributions that are computed by <<
261 G4double <<
262 G4SeltzerBergerModel::ComputeDEDXPerVolume(con <<
263 con <<
264 G4d <<
265 G4d <<
266 { 142 {
267 G4double dedx = 0.0; << 143 // G4cout << "ReadData Z= " << Z << G4endl;
268 if (nullptr == fPrimaryParticle) { << 144 // G4cout << "Status for Z= " << dataSB[Z] << G4endl;
269 SetParticle(p); << 145 //if(path) { G4cout << path << G4endl; }
270 } << 146 if(dataSB[Z]) { return; }
271 if (kineticEnergy <= fLowestKinEnergy) { << 147 const char* datadir = path;
272 return dedx; <<
273 } <<
274 // maximum value of the dE/dx integral (the <<
275 G4double tmax = std::min(cutEnergy, kineticE <<
276 if (tmax == 0.0) { <<
277 return dedx; <<
278 } <<
279 // sets kinematical and material related var <<
280 SetupForMaterial(fPrimaryParticle, material, <<
281 // get element compositions of the material <<
282 const G4ElementVector* theElemVector = mater <<
283 const G4double* theAtomNumDensVector = mater <<
284 const std::size_t numberOfElements = theElem <<
285 // loop over the elements of the material an <<
286 // the restricted dE/dx by numerical integra <<
287 for (std::size_t ie = 0; ie < numberOfElemen <<
288 G4VEmModel::SetCurrentElement((*theElemVec <<
289 G4int Z = (*theElemVector)[ie]->GetZasInt( <<
290 fCurrentIZ = std::min(Z, gMaxZet); <<
291 dedx += (Z*Z)*theAtomNumDensVector[ie]*Com <<
292 } <<
293 // apply the constant factor C/Z = 16\alpha <<
294 dedx *= gBremFactor; <<
295 return std::max(dedx, 0.); <<
296 } <<
297 148
298 // Computes the integral part of the restricte << 149 if(!datadir) {
299 // element (Z) by numerically integrating the << 150 datadir = getenv("G4LEDATA");
300 // k_min=0 and k_max = tmax = min[gamma-cut, e << 151 if(!datadir) {
301 // The numerical integration is done by dividi << 152 G4Exception("G4SeltzerBergerModel::ReadData()","em0006",FatalException,
302 // subintervals and an 8 pint GL integral (on << 153 "Environment variable G4LEDATA not defined");
303 // inteval by tranforming k to alpha=k/E_t (E_ << 154 return;
304 // and each sub-interavl is transformed to [0, <<
305 // in xi(alpha) = xi(k) = [k/E_t-alpha_i]/delt <<
306 // the i = 1,2,..,n-th sub-interval so xi(k) i <<
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 // the ds/dk(Z,k) but ds/dk(Z,k)*[F*k/C] is co <<
311 // (i) what we need here is ds/dk*k and not <<
312 // (ii) the Ter-Mikaelian suppression i.e. F <<
313 // (iii) the constant factor C (includes Z^2 <<
314 G4double G4SeltzerBergerModel::ComputeBremLoss <<
315 { <<
316 // number of intervals and integration step <<
317 const G4double alphaMax = tmax/fPrimaryTotal <<
318 const G4int nSub = (G4int)(20*alphaMax)+3; <<
319 const G4double delta = alphaMax/((G4double)n <<
320 // set minimum value of the first sub-inteva <<
321 G4double alpha_i = 0.0; <<
322 G4double dedxInteg = 0.0; <<
323 for (G4int l = 0; l < nSub; ++l) { <<
324 for (G4int igl = 0; igl < 8; ++igl) { <<
325 // compute the emitted photon energy k <<
326 const G4double k = (alpha_i+gXGL[igl]* <<
327 // compute the DCS value at k (without t <<
328 const G4double dcs = ComputeDXSectionPer <<
329 // account Ter-Mikaelian suppression: ti <<
330 dedxInteg += gWGL[igl]*dcs/(1.0+fDensity <<
331 } 155 }
332 // update sub-interval minimum value <<
333 alpha_i += delta; <<
334 } <<
335 // apply corrections due to variable transfo <<
336 dedxInteg *= delta*fPrimaryTotalEnergy; <<
337 return std::max(dedxInteg,0.); <<
338 } <<
339 <<
340 // Computes restrected atomic cross section by <<
341 // DCS between the proper kinematical limits a <<
342 G4double <<
343 G4SeltzerBergerModel::ComputeCrossSectionPerAt <<
344 G4double kineticEnergy, <<
345 G4double Z, <<
346 G4double, <<
347 G4double cut, <<
348 G4double maxEnergy) <<
349 { <<
350 G4double crossSection = 0.0; <<
351 if (nullptr == fPrimaryParticle) { <<
352 SetParticle(p); <<
353 } <<
354 if (kineticEnergy <= fLowestKinEnergy) { <<
355 return crossSection; <<
356 } <<
357 // min/max kinetic energy limits of the DCS <<
358 const G4double tmin = std::min(cut, kineticE <<
359 const G4double tmax = std::min(maxEnergy, ki <<
360 // zero restricted x-section if e- kinetic e <<
361 if (tmin >= tmax) { <<
362 return crossSection; <<
363 } 156 }
364 fCurrentIZ = std::min(G4lrint(Z), gMaxZet); << 157 std::ostringstream ost;
365 // integrate numerically (dependent part of) << 158 ost << datadir << DirectoryPath() << Z;
366 // a. integrate between tmin and kineticEner << 159 std::ifstream fin(ost.str().c_str());
367 crossSection = ComputeXSectionPerAtom(tmin); << 160 if( !fin.is_open()) {
368 // allow partial integration: only if maxEne << 161 G4ExceptionDescription ed;
369 // b. integrate between tmax and kineticEner << 162 ed << "Bremsstrahlung data file <" << ost.str().c_str()
370 // (so the result in this case is the integr << 163 << "> is not opened!";
371 // maxEnergy) << 164 G4Exception("G4SeltzerBergerModel::ReadData()","em0003",FatalException,
372 if (tmax < kineticEnergy) { << 165 ed,"G4LEDATA version should be G4EMLOW6.23 or later.");
373 crossSection -= ComputeXSectionPerAtom(tma << 166 return;
>> 167 }
>> 168 //G4cout << "G4SeltzerBergerModel read from <" << ost.str().c_str()
>> 169 // << ">" << G4endl;
>> 170 G4Physics2DVector* v = new G4Physics2DVector();
>> 171 if(v->Retrieve(fin)) {
>> 172 if(useBicubicInterpolation) { v->SetBicubicInterpolation(true); }
>> 173 dataSB[Z] = v;
>> 174 static const G4double emaxlog = 4*G4Log(10.);
>> 175 ylimit[Z] = v->Value(0.97, emaxlog, idx, idy);
>> 176 } else {
>> 177 G4ExceptionDescription ed;
>> 178 ed << "Bremsstrahlung data file <" << ost.str().c_str()
>> 179 << "> is not retrieved!";
>> 180 G4Exception("G4SeltzerBergerModel::ReadData()","em0005",FatalException,
>> 181 ed,"G4LEDATA version should be G4EMLOW6.23 or later.");
>> 182 delete v;
374 } 183 }
375 // multiply with the constant factors: 16\al << 184 // G4cout << dataSB[Z] << G4endl;
376 crossSection *= Z*Z*gBremFactor; <<
377 return std::max(crossSection, 0.); <<
378 } 185 }
379 186
380 // Numerical integral of the (k dependent part << 187 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
381 // k_max = E_k (where E_k is the kinetic energ <<
382 // minimum of energy of the emitted photon). <<
383 // transformed alpha(k) = ln(k/E_t) variable ( <<
384 // the primary e-). The integration range is d <<
385 // delta = [ln(k_min/E_t)-ln(k_max/E_t)]/n wid <<
386 // on [0,1] is applied on each sub-inteval so <<
387 // xi(alpha) = xi(k) = [ln(k/E_t)-alpha_i]/del <<
388 // (i-1)*delta for the i = 1,2,..,n-th sub-int <<
389 // sub-intevals. From the transformed xi, k(xi <<
390 // Since the integration is done in variable x <<
391 // transformation results in a multiplicative <<
392 // However, DCS differential in k is ~1/k so t <<
393 // becomes delta and the 1/k factor is dropped <<
394 // Ter-Mikaelian suppression is always account <<
395 G4double G4SeltzerBergerModel::ComputeXSection <<
396 { <<
397 G4double xSection = 0.0; <<
398 const G4double alphaMin = G4Log(tmin/fPrimar <<
399 const G4double alphaMax = G4Log(fPrimaryKinE <<
400 const G4int nSub = (G4int)(0.45*(alphaMax <<
401 const G4double delta = (alphaMax-alphaMin)/( <<
402 // set minimum value of the first sub-inteva <<
403 G4double alpha_i = alphaMin; <<
404 for (G4int l = 0; l < nSub; ++l) { <<
405 for (G4int igl = 0; igl < 8; ++igl) { <<
406 // compute the emitted photon energy k <<
407 const G4double k = G4Exp(alpha_i+gXGL[ig <<
408 // compute the DCS value at k (without t <<
409 const G4double dcs = ComputeDXSectionPer <<
410 // account Ter-Mikaelian suppression: ti <<
411 xSection += gWGL[igl]*dcs/(1.0+fDensityC <<
412 } <<
413 // update sub-interval minimum value <<
414 alpha_i += delta; <<
415 } <<
416 // apply corrections due to variable transfo <<
417 xSection *= delta; <<
418 // final check <<
419 return std::max(xSection, 0.); <<
420 } <<
421 188
422 G4double G4SeltzerBergerModel::ComputeDXSectio 189 G4double G4SeltzerBergerModel::ComputeDXSectionPerAtom(G4double gammaEnergy)
423 { 190 {
424 G4double dxsec = 0.0; << 191
425 if (gammaEnergy < 0.0 || fPrimaryKinEnergy < << 192 if(gammaEnergy < 0.0 || kinEnergy <= 0.0) { return 0.0; }
426 return dxsec; << 193 G4double x = gammaEnergy/kinEnergy;
427 } << 194 G4double y = G4Log(kinEnergy/MeV);
428 // reduced photon energy << 195 G4int Z = G4lrint(currentZ);
429 const G4double x = gammaEnergy/fPrimaryKinEn << 196
430 // l-kinetic energy of the e-/e+ << 197 //G4cout << "G4SeltzerBergerModel::ComputeDXSectionPerAtom Z= " << Z
431 const G4double y = G4Log(fPrimaryKinEnergy/C << 198 // << " x= " << x << " y= " << y << " " << dataSB[Z] << G4endl;
432 // make sure that the Z-related SB-DCS are l << 199 if(nullptr == dataSB[Z]) { InitialiseForElement(0, Z); }
433 // NOTE: fCurrentIZ should have been set bef << 200 /*
434 fCurrentIZ = std::max(std::min(fCurrentIZ, g << 201 G4ExceptionDescription ed;
435 if (nullptr == gSBDCSData[fCurrentIZ]) { << 202 ed << "Bremsstrahlung data for Z= " << Z
436 G4AutoLock l(&theSBMutex); << 203 << " are not initialized!";
437 ReadData(fCurrentIZ); << 204 G4Exception("G4SeltzerBergerModel::ComputeDXSectionPerAtom()","em0005",
438 l.unlock(); << 205 FatalException, ed,
>> 206 "G4LEDATA version should be G4EMLOW6.23 or later.");
439 } 207 }
440 // NOTE: SetupForMaterial should have been c << 208 */
441 const G4double pt2 = fPrimaryKinEnergy*(fPri << 209 G4double invb2 =
442 const G4double invb2 = fPrimaryTotalEnergy*f << 210 totalEnergy*totalEnergy/(kinEnergy*(kinEnergy + 2*particleMass));
443 G4double val = gSBDCSData[fCurrentIZ]->Value << 211 G4double cross = dataSB[Z]->Value(x,y,idx,idy)*invb2*millibarn/bremFactor;
444 dxsec = val*invb2*CLHEP::millibarn/gBremFact << 212
445 // e+ correction << 213 if(!isElectron) {
446 if (!fIsElectron) { << 214 G4double invbeta1 = sqrt(invb2);
447 const G4double invbeta1 = std::sqrt(invb2) << 215 G4double e2 = kinEnergy - gammaEnergy;
448 const G4double e2 = fPrimaryKinEnergy - ga << 216 if(e2 > 0.0) {
449 if (e2 > 0.0) { << 217 G4double invbeta2 = (e2 + particleMass)/sqrt(e2*(e2 + 2*particleMass));
450 const G4double invbeta2 = << 218 static const G4double alpha = CLHEP::twopi*CLHEP::fine_structure_const;
451 (e2 + CLHEP::electron_mass_c2)/std::sqrt(e2* << 219 G4double xxx = alpha*currentZ*(invbeta1 - invbeta2);
452 const G4double dum0 = kAlpha*fCurrentIZ* << 220 if(xxx < expnumlim) { cross = 0.0; }
453 if (dum0 < gExpNumLimit) { << 221 else { cross *= G4Exp(xxx); }
454 dxsec = 0.0; <<
455 } else { <<
456 dxsec *= G4Exp(dum0); <<
457 } <<
458 } else { 222 } else {
459 dxsec = 0.0; << 223 cross = 0.0;
460 } 224 }
461 } 225 }
462 return dxsec; << 226
>> 227 return cross;
463 } 228 }
464 229
465 void << 230 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
466 G4SeltzerBergerModel::SampleSecondaries(std::v << 231
467 const << 232 void
468 const << 233 G4SeltzerBergerModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
469 G4dou << 234 const G4MaterialCutsCouple* couple,
470 G4dou << 235 const G4DynamicParticle* dp,
>> 236 G4double cutEnergy,
>> 237 G4double maxEnergy)
471 { 238 {
472 const G4double kinEnergy = dp->GetKineticEne << 239 G4double kineticEnergy = dp->GetKineticEnergy();
473 const G4double logKinEnergy = dp->GetLogKine << 240 G4double cut = std::min(cutEnergy, kineticEnergy);
474 const G4double tmin = std::min(cutEnergy, ki << 241 G4double emax = std::min(maxEnergy, kineticEnergy);
475 const G4double tmax = std::min(maxEnergy, ki << 242 if(cut >= emax) { return; }
476 if (tmin >= tmax) { << 243
477 return; << 244 SetupForMaterial(particle, couple->GetMaterial(), kineticEnergy);
478 } << 245
479 // set local variables and select target ele << 246 const G4Element* elm =
480 SetupForMaterial(fPrimaryParticle, couple->G << 247 SelectRandomAtom(couple,particle,kineticEnergy,cut,emax);
481 const G4Element* elm = SelectTargetAtom(coup << 248 SetCurrentElement(elm->GetZasInt());
482 logK << 249
483 fCurrentIZ = std::max(std::min(elm->GetZasIn << 250 totalEnergy = kineticEnergy + particleMass;
484 // << 251 densityCorr = densityFactor*totalEnergy*totalEnergy;
485 const G4double totMomentum = std::sqrt(kinEn << 252 G4double totMomentum = sqrt(kineticEnergy*(totalEnergy + electron_mass_c2));
486 /* <<
487 G4cout << "G4SeltzerBergerModel::SampleSecon <<
488 << kinEnergy/MeV <<
489 << " Z= " << fCurrentIZ << " cut(MeV) <<
490 << " emax(MeV)= " << tmax/MeV << " co <<
491 */ <<
492 // sample emitted photon energy either by re <<
493 const G4double gammaEnergy = !fIsUseSampling <<
494 ? SampleEnergyTransfer(kinEnergy, logK <<
495 : gSBSamplingTable->SampleEnergyTransf <<
496 fDensityCorr, fCurrentIZ, <<
497 // should never happen under normal conditio <<
498 if (gammaEnergy <= 0.) { <<
499 return; <<
500 } <<
501 // <<
502 // angles of the emitted gamma. ( Z - axis a <<
503 // general interface <<
504 G4ThreeVector gamDir = GetAngularDistributio <<
505 fPrimaryTotalEnergy-gammaEnergy, fCurrentIZ <<
506 // create G4DynamicParticle object for the e <<
507 auto gamma = new G4DynamicParticle(fGammaPar <<
508 vdp->push_back(gamma); <<
509 // <<
510 // compute post-interaction kinematics of th <<
511 G4ThreeVector dir = <<
512 (totMomentum*dp->GetMomentumDirection() - <<
513 const G4double finalE = kinEnergy - gammaEne <<
514 /* 253 /*
515 G4cout << "### G4SBModel: v= " << 254 G4cout << "G4SeltzerBergerModel::SampleSecondaries E(MeV)= "
516 << " Eg(MeV)= " << gammaEnergy << 255 << kineticEnergy/MeV
517 << " Ee(MeV)= " << kineticEnergy << 256 << " Z= " << Z << " cut(MeV)= " << cut/MeV
518 << " DirE " << direction << " DirG " << 257 << " emax(MeV)= " << emax/MeV << " corr= " << densityCorr << G4endl;
519 << G4endl; <<
520 */ 258 */
521 // if secondary gamma energy is higher than << 259 G4double xmin = G4Log(cut*cut + densityCorr);
522 // then stop tracking the primary particle a << 260 G4double xmax = G4Log(emax*emax + densityCorr);
523 // instead of the primary << 261 G4double y = G4Log(kineticEnergy/MeV);
524 if (gammaEnergy > SecondaryThreshold()) { << 262
525 fParticleChange->ProposeTrackStatus(fStopA << 263 G4double gammaEnergy, v;
526 fParticleChange->SetProposedKineticEnergy( <<
527 auto el = new G4DynamicParticle( <<
528 const_cast<G4ParticleDefinition* <<
529 vdp->push_back(el); <<
530 } else { // continue tracking the primary e- <<
531 fParticleChange->SetProposedMomentumDirect <<
532 fParticleChange->SetProposedKineticEnergy( <<
533 } <<
534 } <<
535 264
536 // sample emitted photon energy by usign rejec <<
537 G4double <<
538 G4SeltzerBergerModel::SampleEnergyTransfer(con <<
539 con <<
540 con <<
541 con <<
542 { <<
543 // min max of the transformed variable: x(k) <<
544 // [ln(k_c^2+k_p^2), ln(E_k^2+k_p^2)] <<
545 const G4double xmin = G4Log(tmin*tmin+fDen <<
546 const G4double xrange = G4Log(tmax*tmax+fDen <<
547 const G4double y = logKinEnergy; <<
548 // majoranta 265 // majoranta
549 const G4double x0 = tmin/kinEnergy; << 266 G4double x0 = cut/kineticEnergy;
550 G4double vmax; 267 G4double vmax;
551 if (nullptr == gSBDCSData[fCurrentIZ]) { << 268 if(currentZ <= 92) {
552 ReadData(fCurrentIZ); << 269 vmax = dataSB[currentZ]->Value(x0, y, idx, idy)*1.02;
>> 270 } else {
>> 271 idx = idy = 0;
>> 272 vmax = dataSB[currentZ]->Value(x0, y, idx, idy)*1.2;
553 } 273 }
554 vmax = gSBDCSData[fCurrentIZ]->Value(x0, y, << 274
555 // << 275 static const G4double epeaklimit= 300*CLHEP::MeV;
556 static const G4double kEPeakLim = 300.*CLHEP << 276 static const G4double elowlimit = 20*CLHEP::keV;
557 static const G4double kELowLim = 20.*CLHEP << 277
558 // majoranta corrected for e- 278 // majoranta corrected for e-
559 if (fIsElectron && x0 < 0.97 && << 279 if(isElectron && x0 < 0.97 &&
560 ((kinEnergy>kEPeakLim) || (kinEnergy<kEL << 280 ((kineticEnergy > epeaklimit) || (kineticEnergy < elowlimit))) {
561 G4double ylim = std::min(gYLimitData[fCurr << 281 G4double ylim = std::min(ylimit[currentZ],1.1*dataSB[currentZ]->Value(0.97,y,idx,idy));
562 1.1*gSBDCSData[fCurre << 282 if(ylim > vmax) { vmax = ylim; }
563 vmax = std::max(vmax, ylim); <<
564 } <<
565 if (x0 < 0.05) { <<
566 vmax *= 1.2; <<
567 } 283 }
>> 284 if(x0 < 0.05) { vmax *= 1.2; }
>> 285
568 //G4cout<<"y= "<<y<<" xmin= "<<xmin<<" xmax= 286 //G4cout<<"y= "<<y<<" xmin= "<<xmin<<" xmax= "<<xmax
569 //<<" vmax= "<<vmax<<G4endl; 287 //<<" vmax= "<<vmax<<G4endl;
570 static const G4int kNCountMax = 100; << 288 static const G4int ncountmax = 100;
571 CLHEP::HepRandomEngine* rndmEngine = G4Rando 289 CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
572 G4double rndm[2]; 290 G4double rndm[2];
573 G4double gammaEnergy, v; << 291
574 for (G4int nn = 0; nn < kNCountMax; ++nn) { << 292 for(G4int nn=0; nn<ncountmax; ++nn) {
575 rndmEngine->flatArray(2, rndm); 293 rndmEngine->flatArray(2, rndm);
576 gammaEnergy = << 294 G4double x = G4Exp(xmin + rndm[0]*(xmax - xmin)) - densityCorr;
577 std::sqrt(std::max(G4Exp(xmin + rndm[0]* << 295 if(x < 0.0) { x = 0.0; }
578 v = gSBDCSData[fCurrentIZ]->Value(gammaEne << 296 gammaEnergy = sqrt(x);
579 // e+ correction << 297 G4double x1 = gammaEnergy/kineticEnergy;
580 if (!fIsElectron) { << 298 v = dataSB[currentZ]->Value(x1, y, idx, idy);
581 const G4double e1 = kinEnergy - tmin; << 299
582 const G4double invbeta1 = << 300 // correction for positrons
583 (e1 + CLHEP::electron_mass_c2)/std::sqrt(e1* << 301 if(!isElectron) {
584 const G4double e2 = kinEnergy-gamm << 302 G4double e1 = kineticEnergy - cut;
585 const G4double invbeta2 = << 303 G4double invbeta1 = (e1 + particleMass)/sqrt(e1*(e1 + 2*particleMass));
586 (e2 + CLHEP::electron_mass_c2)/std::sqrt(e2* << 304 G4double e2 = kineticEnergy - gammaEnergy;
587 const G4double dum0 = kAlpha*fCurren << 305 G4double invbeta2 = (e2 + particleMass)/sqrt(e2*(e2 + 2*particleMass));
588 if (dum0 < gExpNumLimit) { << 306 G4double xxx = twopi*fine_structure_const*currentZ*(invbeta1 - invbeta2);
589 v = 0.0; << 307
590 } else { << 308 if(xxx < expnumlim) { v = 0.0; }
591 v *= G4Exp(dum0); << 309 else { v *= G4Exp(xxx); }
592 } << 310 }
593 } << 311
594 if (v > 1.05*vmax && fNumWarnings < 11) { << 312 if (v > 1.05*vmax && nwarn < 5) {
595 ++fNumWarnings; << 313 ++nwarn;
596 G4ExceptionDescription ed; 314 G4ExceptionDescription ed;
597 ed << "### G4SeltzerBergerModel Warning: 315 ed << "### G4SeltzerBergerModel Warning: Majoranta exceeded! "
598 << v << " > " << vmax << " by " << v/ 316 << v << " > " << vmax << " by " << v/vmax
599 << " Niter= " << nn << 317 << " Niter= " << nn
600 << " Egamma(MeV)= " << gammaEnergy 318 << " Egamma(MeV)= " << gammaEnergy
601 << " Ee(MeV)= " << kinEnergy << 319 << " Ee(MeV)= " << kineticEnergy
602 << " Z= " << fCurrentIZ << " " << fP << 320 << " Z= " << currentZ << " " << particle->GetParticleName();
603 // << 321
604 if (10 == fNumWarnings) { << 322 if ( 20 == nwarn ) {
605 ed << "\n ### G4SeltzerBergerModel War 323 ed << "\n ### G4SeltzerBergerModel Warnings stopped";
606 } 324 }
607 G4Exception("G4SeltzerBergerModel::Sampl 325 G4Exception("G4SeltzerBergerModel::SampleScattering","em0044",
608 JustWarning, ed,""); 326 JustWarning, ed,"");
>> 327
609 } 328 }
610 if (v >= vmax*rndm[1]) { << 329 if(v >= vmax*rndm[1]) { break; }
611 break; << 330 }
612 } << 331
>> 332 //
>> 333 // angles of the emitted gamma. ( Z - axis along the parent particle)
>> 334 // use general interface
>> 335 //
>> 336
>> 337 G4ThreeVector gammaDirection =
>> 338 GetAngularDistribution()->SampleDirection(dp, totalEnergy-gammaEnergy,
>> 339 currentZ, couple->GetMaterial());
>> 340
>> 341 // create G4DynamicParticle object for the Gamma
>> 342 G4DynamicParticle* gamma =
>> 343 new G4DynamicParticle(theGamma,gammaDirection,gammaEnergy);
>> 344 vdp->push_back(gamma);
>> 345
>> 346 G4ThreeVector direction = (totMomentum*dp->GetMomentumDirection()
>> 347 - gammaEnergy*gammaDirection).unit();
>> 348
>> 349 /*
>> 350 G4cout << "### G4SBModel: v= "
>> 351 << " Eg(MeV)= " << gammaEnergy
>> 352 << " Ee(MeV)= " << kineticEnergy
>> 353 << " DirE " << direction << " DirG " << gammaDirection
>> 354 << G4endl;
>> 355 */
>> 356 // energy of primary
>> 357 G4double finalE = kineticEnergy - gammaEnergy;
>> 358
>> 359 // stop tracking and create new secondary instead of primary
>> 360 if(gammaEnergy > SecondaryThreshold()) {
>> 361 fParticleChange->ProposeTrackStatus(fStopAndKill);
>> 362 fParticleChange->SetProposedKineticEnergy(0.0);
>> 363 G4DynamicParticle* el =
>> 364 new G4DynamicParticle(const_cast<G4ParticleDefinition*>(particle),
>> 365 direction, finalE);
>> 366 vdp->push_back(el);
>> 367
>> 368 // continue tracking
>> 369 } else {
>> 370 fParticleChange->SetProposedMomentumDirection(direction);
>> 371 fParticleChange->SetProposedKineticEnergy(finalE);
613 } 372 }
614 return gammaEnergy; <<
615 } 373 }
>> 374
>> 375 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 376
>> 377 #include "G4AutoLock.hh"
>> 378 namespace { G4Mutex SeltzerBergerModelMutex = G4MUTEX_INITIALIZER; }
>> 379 void G4SeltzerBergerModel::InitialiseForElement(const G4ParticleDefinition*,
>> 380 G4int Z)
>> 381 {
>> 382 G4AutoLock l(&SeltzerBergerModelMutex);
>> 383 // G4cout << "G4SeltzerBergerModel::InitialiseForElement Z= " << Z << G4endl;
>> 384 if(nullptr == dataSB[Z]) { ReadData(Z); }
>> 385 }
>> 386
>> 387 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 388
>> 389
616 390