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