| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // -------------------------------------------
27 //
28 // GEANT4 Class file
29 //
30 //
31 // File name: G4SeltzerBergerModel
32 //
33 // Author: Vladimir Ivanchenko use inhe
34 // base class implementing ultr
35 // model
36 //
37 // Creation date: 04.10.2011
38 //
39 // Modifications:
40 //
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 // -------------------------------------------
50 //
51
52 #include "G4SeltzerBergerModel.hh"
53 #include "G4PhysicalConstants.hh"
54 #include "G4SystemOfUnits.hh"
55 #include "Randomize.hh"
56 #include "G4Material.hh"
57 #include "G4Element.hh"
58 #include "G4ElementVector.hh"
59 #include "G4ParticleChangeForLoss.hh"
60 #include "G4SBBremTable.hh"
61 #include "G4ModifiedTsai.hh"
62
63 #include "G4EmParameters.hh"
64 #include "G4ProductionCutsTable.hh"
65 #include "G4Gamma.hh"
66 #include "G4Electron.hh"
67
68 #include "G4Physics2DVector.hh"
69 #include "G4Exp.hh"
70 #include "G4Log.hh"
71 #include "G4AutoLock.hh"
72
73 #include "G4ios.hh"
74
75 #include <fstream>
76 #include <iomanip>
77 #include <sstream>
78 #include <thread>
79
80 G4double G4SeltzerBergerModel::gYLimitData[] =
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
98 namespace
99 {
100 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 }
112
113 G4SeltzerBergerModel::G4SeltzerBergerModel(con
114 con
115 : G4VEmModel(nam),
116 fGammaParticle(G4Gamma::Gamma()),
117 fLowestKinEnergy(1.0*CLHEP::keV)
118 {
119 SetLowEnergyLimit(fLowestKinEnergy);
120 SetAngularDistribution(new G4ModifiedTsai())
121 if (fPrimaryParticle != p) { SetParticle(p);
122 }
123
124 G4SeltzerBergerModel::~G4SeltzerBergerModel()
125 {
126 // delete SB-DCS data per Z
127 if (isInitializer) {
128 for (std::size_t iz = 0; iz < gMaxZet; ++i
129 if (gSBDCSData[iz]) {
130 delete gSBDCSData[iz];
131 gSBDCSData[iz] = nullptr;
132 }
133 }
134 if (gSBSamplingTable) {
135 delete gSBSamplingTable;
136 gSBSamplingTable = nullptr;
137 }
138 }
139 }
140
141 void G4SeltzerBergerModel::Initialise(const G4
142 const G4
143 {
144 // parameters in each thread
145 if (fPrimaryParticle != p) {
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
165 // init sampling tables if it was requeste
166 if (fIsUseSamplingTables) {
167 if (nullptr == gSBSamplingTable) {
168 gSBSamplingTable = new G4SBBremTable()
169 }
170 gSBSamplingTable->Initialize(std::max(fL
171 HighEnergyL
172 }
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 }
188 }
189
190 void G4SeltzerBergerModel::InitialiseLocal(con
191
192 {
193 SetElementSelectors(masterModel->GetElementS
194 }
195
196 void G4SeltzerBergerModel::SetParticle(const G
197 {
198 fPrimaryParticle = p;
199 fIsElectron = (p == G4Electron::Electron());
200 }
201
202 void G4SeltzerBergerModel::ReadData(G4int Z) {
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 {
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 }
421
422 G4double G4SeltzerBergerModel::ComputeDXSectio
423 {
424 G4double dxsec = 0.0;
425 if (gammaEnergy < 0.0 || fPrimaryKinEnergy <
426 return dxsec;
427 }
428 // reduced photon energy
429 const G4double x = gammaEnergy/fPrimaryKinEn
430 // l-kinetic energy of the e-/e+
431 const G4double y = G4Log(fPrimaryKinEnergy/C
432 // make sure that the Z-related SB-DCS are l
433 // NOTE: fCurrentIZ should have been set bef
434 fCurrentIZ = std::max(std::min(fCurrentIZ, g
435 if (nullptr == gSBDCSData[fCurrentIZ]) {
436 G4AutoLock l(&theSBMutex);
437 ReadData(fCurrentIZ);
438 l.unlock();
439 }
440 // NOTE: SetupForMaterial should have been c
441 const G4double pt2 = fPrimaryKinEnergy*(fPri
442 const G4double invb2 = fPrimaryTotalEnergy*f
443 G4double val = gSBDCSData[fCurrentIZ]->Value
444 dxsec = val*invb2*CLHEP::millibarn/gBremFact
445 // e+ correction
446 if (!fIsElectron) {
447 const G4double invbeta1 = std::sqrt(invb2)
448 const G4double e2 = fPrimaryKinEnergy - ga
449 if (e2 > 0.0) {
450 const G4double invbeta2 =
451 (e2 + CLHEP::electron_mass_c2)/std::sqrt(e2*
452 const G4double dum0 = kAlpha*fCurrentIZ*
453 if (dum0 < gExpNumLimit) {
454 dxsec = 0.0;
455 } else {
456 dxsec *= G4Exp(dum0);
457 }
458 } else {
459 dxsec = 0.0;
460 }
461 }
462 return dxsec;
463 }
464
465 void
466 G4SeltzerBergerModel::SampleSecondaries(std::v
467 const
468 const
469 G4dou
470 G4dou
471 {
472 const G4double kinEnergy = dp->GetKineticEne
473 const G4double logKinEnergy = dp->GetLogKine
474 const G4double tmin = std::min(cutEnergy, ki
475 const G4double tmax = std::min(maxEnergy, ki
476 if (tmin >= tmax) {
477 return;
478 }
479 // set local variables and select target ele
480 SetupForMaterial(fPrimaryParticle, couple->G
481 const G4Element* elm = SelectTargetAtom(coup
482 logK
483 fCurrentIZ = std::max(std::min(elm->GetZasIn
484 //
485 const G4double totMomentum = std::sqrt(kinEn
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 /*
515 G4cout << "### G4SBModel: v= "
516 << " Eg(MeV)= " << gammaEnergy
517 << " Ee(MeV)= " << kineticEnergy
518 << " DirE " << direction << " DirG "
519 << G4endl;
520 */
521 // if secondary gamma energy is higher than
522 // then stop tracking the primary particle a
523 // instead of the primary
524 if (gammaEnergy > SecondaryThreshold()) {
525 fParticleChange->ProposeTrackStatus(fStopA
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
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
549 const G4double x0 = tmin/kinEnergy;
550 G4double vmax;
551 if (nullptr == gSBDCSData[fCurrentIZ]) {
552 ReadData(fCurrentIZ);
553 }
554 vmax = gSBDCSData[fCurrentIZ]->Value(x0, y,
555 //
556 static const G4double kEPeakLim = 300.*CLHEP
557 static const G4double kELowLim = 20.*CLHEP
558 // majoranta corrected for e-
559 if (fIsElectron && x0 < 0.97 &&
560 ((kinEnergy>kEPeakLim) || (kinEnergy<kEL
561 G4double ylim = std::min(gYLimitData[fCurr
562 1.1*gSBDCSData[fCurre
563 vmax = std::max(vmax, ylim);
564 }
565 if (x0 < 0.05) {
566 vmax *= 1.2;
567 }
568 //G4cout<<"y= "<<y<<" xmin= "<<xmin<<" xmax=
569 //<<" vmax= "<<vmax<<G4endl;
570 static const G4int kNCountMax = 100;
571 CLHEP::HepRandomEngine* rndmEngine = G4Rando
572 G4double rndm[2];
573 G4double gammaEnergy, v;
574 for (G4int nn = 0; nn < kNCountMax; ++nn) {
575 rndmEngine->flatArray(2, rndm);
576 gammaEnergy =
577 std::sqrt(std::max(G4Exp(xmin + rndm[0]*
578 v = gSBDCSData[fCurrentIZ]->Value(gammaEne
579 // e+ correction
580 if (!fIsElectron) {
581 const G4double e1 = kinEnergy - tmin;
582 const G4double invbeta1 =
583 (e1 + CLHEP::electron_mass_c2)/std::sqrt(e1*
584 const G4double e2 = kinEnergy-gamm
585 const G4double invbeta2 =
586 (e2 + CLHEP::electron_mass_c2)/std::sqrt(e2*
587 const G4double dum0 = kAlpha*fCurren
588 if (dum0 < gExpNumLimit) {
589 v = 0.0;
590 } else {
591 v *= G4Exp(dum0);
592 }
593 }
594 if (v > 1.05*vmax && fNumWarnings < 11) {
595 ++fNumWarnings;
596 G4ExceptionDescription ed;
597 ed << "### G4SeltzerBergerModel Warning:
598 << v << " > " << vmax << " by " << v/
599 << " Niter= " << nn
600 << " Egamma(MeV)= " << gammaEnergy
601 << " Ee(MeV)= " << kinEnergy
602 << " Z= " << fCurrentIZ << " " << fP
603 //
604 if (10 == fNumWarnings) {
605 ed << "\n ### G4SeltzerBergerModel War
606 }
607 G4Exception("G4SeltzerBergerModel::Sampl
608 JustWarning, ed,"");
609 }
610 if (v >= vmax*rndm[1]) {
611 break;
612 }
613 }
614 return gammaEnergy;
615 }
616