| 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 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4eBremsstrahlungRelModel
33 //
34 // Author: Andreas Schaelicke
35 //
36 // Creation date: 12.08.2008
37 //
38 // Modifications:
39 //
40 // 13.11.08 add SetLPMflag and SetLPMconsta
41 // 13.11.08 change default LPMconstant valu
42 // 13.10.10 add angular distributon interfa
43 // 31.05.16 change LPMconstant such that it
44 // that consistent to Migdal's one
45 // computation; better agreement w
46 // 15.07.18 improved LPM suppression functi
47 // steps), code cleanup and optimi
48 // model related comments, consist
49 //
50 // Main References:
51 // Y.-S.Tsai, Rev. Mod. Phys. 46 (1974) 815;
52 // S.Klein, Rev. Mod. Phys. 71 (1999) 1501.
53 // T.Stanev et.al., Phys. Rev. D25 (1982) 129
54 // M.L.Ter-Mikaelian, High-energy Electromagn
55 // Wiley, 1972.
56 //
57 // -------------------------------------------
58 //
59
60 #include "G4eBremsstrahlungRelModel.hh"
61 #include "G4PhysicalConstants.hh"
62 #include "G4SystemOfUnits.hh"
63 #include "G4Electron.hh"
64 #include "G4Gamma.hh"
65 #include "Randomize.hh"
66 #include "G4Material.hh"
67 #include "G4Element.hh"
68 #include "G4ElementVector.hh"
69 #include "G4ParticleChangeForLoss.hh"
70 #include "G4ModifiedTsai.hh"
71 #include "G4Exp.hh"
72 #include "G4Log.hh"
73 #include "G4Pow.hh"
74 #include "G4EmParameters.hh"
75 #include "G4AutoLock.hh"
76 #include <thread>
77
78 const G4int G4eBremsstrahlungRelModel::gMaxZet
79
80 // constant DCS factor: 16\alpha r_0^2/3
81 const G4double G4eBremsstrahlungRelModel::gBre
82 = 16. * CLHEP::fine_structure_const * CLHEP:
83 * CLHEP::classic_electr_radius/3.;
84
85 // Migdal's constant: 4\pi r_0*electron_reduce
86 const G4double G4eBremsstrahlungRelModel::gMig
87 = 4. * CLHEP::pi * CLHEP::classic_electr_rad
88 * CLHEP::electron_Compton_length * CLHEP::
89
90 // LPM constant: \alpha(mc^2)^2/(4\pi*\hbar c)
91 const G4double G4eBremsstrahlungRelModel::gLPM
92 = CLHEP::fine_structure_const * CLHEP::elect
93 * CLHEP::electron_mass_c2 / (4. * CLHEP::p
94
95 // abscissas and weights of an 8 point Gauss-L
96 // for numerical integration on [0,1]
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 {
148 fGammaParticle = G4Gamma::Gamma();
149 //
150 fLowestKinEnergy = 1.0*CLHEP::MeV;
151 SetLowEnergyLimit(fLowestKinEnergy);
152 //
153 fLPMEnergyThreshold = 1.e+39;
154 fLPMEnergy = 0.;
155 SetAngularDistribution(new G4ModifiedTsai())
156 //
157 if (nullptr != p) {
158 SetParticle(p);
159 }
160 }
161
162 G4eBremsstrahlungRelModel::~G4eBremsstrahlungR
163 {
164 if (fIsInitializer) {
165 // clear ElementData container
166 for (auto const & ptr : *fElementData) { d
167 fElementData->clear();
168 // clear LPMFunctions (if any)
169 if (fLPMFuncs->fIsInitialized) {
170 fLPMFuncs->fLPMFuncG.clear();
171 fLPMFuncs->fLPMFuncPhi.clear();
172 fLPMFuncs->fIsInitialized = false;
173 }
174 }
175 }
176
177 void G4eBremsstrahlungRelModel::Initialise(con
178 con
179 {
180 // parameters in each thread
181 if (fPrimaryParticle != p) {
182 SetParticle(p);
183 }
184 fUseLPM = G4EmParameters::Instance()->LPM();
185 fCurrentIZ = 0;
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 }
213 }
214
215 void G4eBremsstrahlungRelModel::InitialiseLoca
216
217 {
218 SetElementSelectors(masterModel->GetElementS
219 }
220
221 void G4eBremsstrahlungRelModel::SetParticle(co
222 {
223 fPrimaryParticle = p;
224 fPrimaryParticleMass = p->GetPDGMass();
225 fIsElectron = (p==G4Electron::Elect
226 }
227
228 // Sets kinematical variables like E_kin, E_t
229 // variables like LPM energy and characteristi
230 // k_p^2) for the Ter-Mikaelian suppression ef
231 void G4eBremsstrahlungRelModel::SetupForMateri
232
233
234 {
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 {
267 G4double dedx = 0.0;
268 if (nullptr == fPrimaryParticle) {
269 SetParticle(p);
270 }
271 if (kineticEnergy < LowEnergyLimit()) {
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 zet = (*theElemVector)[ie]->GetZasIn
290 fCurrentIZ = std::min(zet, gMaxZet);
291 dedx += (zet*zet)*theAtomNumD
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 // 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 }
335 // update sub-interval minimum value
336 alpha_i += delta;
337 }
338 // apply corrections due to variable transfo
339 dedxInteg *= delta*fPrimaryTotalEnergy;
340 return std::max(dedxInteg,0.);
341 }
342
343 // Computes restrected atomic cross section by
344 // DCS between the proper kinematical limits a
345 G4double G4eBremsstrahlungRelModel::ComputeCro
346
347
348
349
350
351
352 {
353 G4double crossSection = 0.0;
354 if (nullptr == fPrimaryParticle) {
355 SetParticle(p);
356 }
357 if (kineticEnergy < LowEnergyLimit()) {
358 return crossSection;
359 }
360 // min/max kinetic energy limits of the DCS
361 const G4double tmin = std::min(cut, kineticE
362 const G4double tmax = std::min(maxEnergy, ki
363 // zero restricted x-section if e- kinetic e
364 if (tmin >= tmax) {
365 return crossSection;
366 }
367 fCurrentIZ = std::min(G4lrint(Z), gMaxZet);
368 // integrate numerically (dependent part of)
369 // a. integrate between tmin and kineticEner
370 crossSection = ComputeXSectionPerAtom(tmin);
371 // allow partial integration: only if maxEne
372 // b. integrate between tmax and kineticEner
373 // (so the result in this case is the integr
374 // maxEnergy)
375 if (tmax < kineticEnergy) {
376 crossSection -= ComputeXSectionPerAtom(tma
377 }
378 // multiply with the constant factors: 16\al
379 crossSection *= Z*Z*gBremFactor;
380 return std::max(crossSection, 0.);
381 }
382
383 // Numerical integral of the (k dependent part
384 // k_max = E_k (where E_k is the kinetic energ
385 // minimum of energy of the emitted photon).
386 // transformed alpha(k) = ln(k/E_t) variable (
387 // the primary e-). The integration range is d
388 // delta = [ln(k_min/E_t)-ln(k_max/E_t)]/n wid
389 // on [0,1] is applied on each sub-inteval so
390 // xi(alpha) = xi(k) = [ln(k/E_t)-alpha_i]/del
391 // (i-1)*delta for the i = 1,2,..,n-th sub-int
392 // sub-intevals. From the transformed xi, k(xi
393 // Since the integration is done in variable x
394 // transformation results in a multiplicative
395 // However, DCS differential in k is ~1/k so t
396 // becomes delta and the 1/k factor is dropped
397 // NOTE:
398 // - LPM suppression is accounted above thre
399 // flag is set in SetUpForMaterial() => 2
400 // - Ter-Mikaelian suppression is always acc
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 }
421 // update sub-interval minimum value
422 alpha_i += delta;
423 }
424 // apply corrections due to variable transfo
425 xSection *= delta;
426 // final check
427 return std::max(xSection, 0.);
428 }
429
430 // DCS WITH LPM EFFECT: complete screening apr
431 // ds/dk(Z,k) = C/[F*k]*{ Xi(s*F)*[y^2*G/4 +(1
432 // +(1-y)*[1+1/Z]/12}
433 // Xi(s),G(s), Phi(s) are LPM suppression func
434 //
435 // LPM SUPPRESSION: The 's' is the suppression
436 // 1+(k_p/k)^2 with k_p = hbar*w_p*E/(m*c^2) i
437 // dependent constant. F accounts the Ter-Mika
438 // transition in the emitted photon energy. Al
439 // goes to 0 when s goes to 0 and goes to 1 wh
440 // So evaluating the LPM suppression functions
441 // smooth transition depending on the emitted
442 // smoothly turned off i.e. Xi(sF)=G(sF)=Phi(s
443 // and sF ~ s when k >> k_p since F ~ 1 in tha
444 // HERE, ds/dk(Z,k)*[F*k/C] is computed since:
445 // (i) DCS ~ 1/k factor will disappear due
446 // v(k)=ln(k/E_t) -> dk/dv=E_t*e^v=k ->
447 // would cnacell out the 1/k factor =>
448 // (ii) the constant factor C and Z don't de
449 // (iii) the 1/F(k) factor is accounted in th
450 // tion computation) or implicitly thro
451 // (in the final state sampling algorit
452 // COMPLETE SCREENING: see more at the DCS wit
453 G4double
454 G4eBremsstrahlungRelModel::ComputeRelDXSection
455 {
456 G4double dxsec = 0.0;
457 if (gammaEnergy < 0.) {
458 return dxsec;
459 }
460 const G4double y = gammaEnergy/fPrimaryT
461 const G4double onemy = 1.-y;
462 const G4double dum0 = 0.25*y*y;
463 // evaluate LPM functions (combined with the
464 G4double funcGS, funcPhiS, funcXiS;
465 ComputeLPMfunctions(funcXiS, funcGS, funcPhi
466 const ElementData* elDat = (*fElementData)[f
467 const G4double term1 = funcXiS*(dum0*fun
468 dxsec = term1*elDat->fZFactor1+onemy*elDat->
469 //
470 if (fIsScatOffElectron) {
471 fSumTerm = dxsec;
472 fNucTerm = term1*elDat->fZFactor11 + onemy
473 }
474 return std::max(dxsec,0.0);
475 }
476
477 // DCS WITHOUT LPM EFFECT: DCS with sceening (
478 // ds/dk(Z,k)=C/[F*k]*{(1-y+3*y^2/4)*[(0.25*ph
479 // -2*ln(Z)/3)/Z]+ (1-y)*[(phi1(g)-phi2(g))+(p
480 // where f_c(Z) is the Coulomb correction fact
481 // psi2(e) are coherent and incoherent screeni
482 // model of the atom, the screening functions
483 // depend on Z (not explicitly). These numeric
484 // approximated as Tsai Eqs. [3.38-3.41] with
485 // e=epsilon given by Tsai Eqs. [3.30 and 3.31
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 }
519 } else {
520 // use Tsai's analytical approx. (Tsai Eqs
521 // numerical screening functions computed
522 const G4double invZ = 1./(G4double)fCur
523 const G4double Fz = elDat->fFz;
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 }
537 }
538 return std::max(dxsec,0.0);
539 }
540
541 // Coherent and incoherent screening function
542 // Eqs.[3.38-3.41]). Tsai's analytical approxi
543 // functions computed by using the Thomas-Ferm
544 // ximation to the numerical TF screening func
545 // screening functions can be expressed in a '
546 // pendent variable (see Tsai Eqs. Eqs. [3.30
547 void G4eBremsstrahlungRelModel::ComputeScreeni
548
549
550
551
552
553 {
554 const G4double gam2 = gam*gam;
555 phi1 = 16.863-2.0*G4Log(1.0+0.311877*gam2)
556 +1.6*G4Exp(-1.5*gam);
557 phi1m2 = 2.0/(3.0+19.5*gam+18.0*gam2); //
558 const G4double eps2 = eps*eps;
559 psi1 = 24.34-2.0*G4Log(1.0+13.111641*eps2)
560 +1.2*G4Exp(-29.2*eps);
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 }
581 //
582 SetupForMaterial(fPrimaryParticle, couple->G
583 const G4Element* elm = SelectTargetAtom(coup
584 dp->
585 //
586 fCurrentIZ = elm->GetZasInt();
587 const ElementData* elDat = (*fElementData)[f
588 const G4double funcMax = elDat->fZFactor1+el
589 // get the random engine
590 G4double rndm[2];
591 CLHEP::HepRandomEngine* rndmEngine = G4Rando
592 // min max of the transformed variable: x(k)
593 const G4double xmin = G4Log(tmin*tmin+fDen
594 const G4double xrange = G4Log(tmax*tmax+fDen
595 G4double gammaEnergy, funcVal;
596 do {
597 rndmEngine->flatArray(2, rndm);
598 gammaEnergy = std::sqrt(std::max(G4Exp(xmi
599 funcVal = fIsLPMActive
600 ? ComputeRelDXSectionPerAtom(
601 : ComputeDXSectionPerAtom(gam
602 // cross-check of proper function maximum
603 // if (funcVal > funcMax) {
604 // G4cout << "### G4eBremsstrahlungRelMod
605 // << funcVal << " > " << funcMax
606 // << " Egamma(MeV)= " << gammaEnergy
607 // << " Ee(MeV)= " << kineticEnergy
608 // << " " << GetName()
609 // << G4endl;
610 // }
611 // Loop checking, 03-Aug-2015, Vladimir Iv
612 } while (funcVal < funcMax*rndm[1]);
613 //
614 // scattering off nucleus or off e- by tripl
615 if (fIsScatOffElectron && rndmEngine->flat()
616 GetTripletModel()->SampleSecondaries(vdp,
617 return;
618 }
619 //
620 // angles of the emitted gamma. ( Z - axis a
621 // use general interface
622 G4ThreeVector gamDir =
623 GetAngularDistribution()->SampleDirection(
624
625 // create G4DynamicParticle object for the G
626 auto gamma = new G4DynamicParticle(fGammaPar
627 vdp->push_back(gamma);
628 // compute post-interaction kinematics of pr
629 // energy-momentum conservation
630 const G4double totMomentum = std::sqrt(kinet
631 fPrimaryTotalEn
632 G4ThreeVector dir =
633 (totMomentum*dp->GetMomentumDirec
634 const G4double finalE = kineticEnergy-gamm
635 // if secondary gamma energy is higher than
636 // then stop tracking the primary particle a
637 // instead of the primary one
638 if (gammaEnergy > SecondaryThreshold()) {
639 fParticleChange->ProposeTrackStatus(fStopA
640 fParticleChange->SetProposedKineticEnergy(
641 auto el = new G4DynamicParticle(
642 const_cast<G4ParticleDefinition*
643 vdp->push_back(el);
644 } else { // continue tracking the primary e-
645 fParticleChange->SetProposedMomentumDirect
646 fParticleChange->SetProposedKineticEnergy(
647 }
648 }
649
650 void G4eBremsstrahlungRelModel::InitialiseElem
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 {
733 const G4double varShat2 = varShat*varShat;
734 const G4double varShat3 = varShat*varShat2
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 }
765 }
766
767 // s goes up to 2 with ds = 0.01 to be the def
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
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
802