Geant4-11
Public Member Functions | Private Member Functions | Private Attributes
G4eIonisationSpectrum Class Reference

#include <G4eIonisationSpectrum.hh>

Inheritance diagram for G4eIonisationSpectrum:
G4VEnergySpectrum

Public Member Functions

G4double AverageEnergy (G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell, const G4ParticleDefinition *pd=0) const
 
G4double Excitation (G4int Z, G4double e) const
 
 G4eIonisationSpectrum ()
 
G4double MaxEnergyOfSecondaries (G4double kineticEnergy, G4int Z=0, const G4ParticleDefinition *pd=0) const
 
void PrintData () const
 
G4double Probability (G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell, const G4ParticleDefinition *pd=0) const
 
G4double SampleEnergy (G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell, const G4ParticleDefinition *pd=0) const
 
 ~G4eIonisationSpectrum ()
 

Private Member Functions

G4double AverageValue (G4double xMin, G4double xMax, const G4DataVector &p) const
 
G4double Function (G4double x, const G4DataVector &p) const
 
 G4eIonisationSpectrum (const G4eIonisationSpectrum &)
 
G4double IntSpectrum (G4double xMin, G4double xMax, const G4DataVector &p) const
 
G4eIonisationSpectrumoperator= (const G4eIonisationSpectrum &right)
 

Private Attributes

G4double factor
 
G4int iMax
 
G4double lowestE
 
G4eIonisationParameterstheParam
 
G4int verbose
 

Detailed Description

Definition at line 63 of file G4eIonisationSpectrum.hh.

Constructor & Destructor Documentation

◆ G4eIonisationSpectrum() [1/2]

G4eIonisationSpectrum::G4eIonisationSpectrum ( )

◆ ~G4eIonisationSpectrum()

G4eIonisationSpectrum::~G4eIonisationSpectrum ( )

Definition at line 71 of file G4eIonisationSpectrum.cc.

72{
73 delete theParam;
74}

References theParam.

◆ G4eIonisationSpectrum() [2/2]

G4eIonisationSpectrum::G4eIonisationSpectrum ( const G4eIonisationSpectrum )
private

Member Function Documentation

◆ AverageEnergy()

G4double G4eIonisationSpectrum::AverageEnergy ( G4int  Z,
G4double  tMin,
G4double  tMax,
G4double  kineticEnergy,
G4int  shell,
const G4ParticleDefinition pd = 0 
) const
virtual

Implements G4VEnergySpectrum.

Definition at line 172 of file G4eIonisationSpectrum.cc.

178{
179 // Please comment what AverageEnergy does and what are the three
180 // functions mentioned below
181 // Describe the algorithms used
182
184 G4double t0 = std::max(tMin, lowestE);
185 G4double tm = std::min(tMax, eMax);
186 if(t0 >= tm) return 0.0;
187
189 Shell(Z, shell)->BindingEnergy();
190
191 if(e <= bindingEnergy) return 0.0;
192
194
195 G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
197
198 if(verbose > 1) {
199 G4cout << "G4eIonisationSpectrum::AverageEnergy: Z= " << Z
200 << "; shell= " << shell
201 << "; E(keV)= " << e/keV
202 << "; bindingE(keV)= " << bindingEnergy/keV
203 << "; x1= " << x1
204 << "; x2= " << x2
205 << G4endl;
206 }
207
208 G4DataVector p;
209
210 // Access parameters
211 for (G4int i=0; i<iMax; i++)
212 {
213 G4double x = theParam->Parameter(Z, shell, i, e);
214 if(i<4) x /= energy;
215 p.push_back(x);
216 }
217
218 if(p[3] > 0.5) p[3] = 0.5;
219
220 G4double gLocal2 = energy/electron_mass_c2 + 1.;
221 p.push_back((2.0*gLocal2 - 1.0)/(gLocal2*gLocal2));
222
223
224 //Add protection against division by zero: actually in Function()
225 //parameter p[3] appears in the denominator. Bug report 1042
226 if (p[3] > 0)
227 p[iMax-1] = Function(p[3], p);
228 else
229 {
230 G4cout << "WARNING: G4eIonisationSpectrum::AverageEnergy "
231 << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = "
232 << Z << ". Please check and/or update it " << G4endl;
233 }
234
235 G4double val = AverageValue(x1, x2, p);
237 G4double nor = IntSpectrum(x0, 0.5, p);
238 val *= energy;
239
240 if(verbose > 1) {
241 G4cout << "tcut(MeV)= " << tMin/MeV
242 << "; tMax(MeV)= " << tMax/MeV
243 << "; x0= " << x0
244 << "; x1= " << x1
245 << "; x2= " << x2
246 << "; val= " << val
247 << "; nor= " << nor
248 << "; sum= " << p[0]
249 << "; a= " << p[1]
250 << "; b= " << p[2]
251 << "; c= " << p[3]
252 << G4endl;
253 }
254
255 p.clear();
256
257 if(nor > 0.0) val /= nor;
258 else val = 0.0;
259
260 return val;
261}
static constexpr double keV
Definition: G4SIunits.hh:202
static constexpr double MeV
Definition: G4SIunits.hh:200
double G4double
Definition: G4Types.hh:83
int G4int
Definition: G4Types.hh:85
const G4int Z[17]
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
static G4AtomicTransitionManager * Instance()
G4double Parameter(G4int Z, G4int shellIndex, G4int parameterIndex, G4double e) const
G4double IntSpectrum(G4double xMin, G4double xMax, const G4DataVector &p) const
G4double MaxEnergyOfSecondaries(G4double kineticEnergy, G4int Z=0, const G4ParticleDefinition *pd=0) const
G4double Function(G4double x, const G4DataVector &p) const
G4double AverageValue(G4double xMin, G4double xMax, const G4DataVector &p) const
G4double energy(const ThreeVector &p, const G4double m)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double bindingEnergy(G4int A, G4int Z)
float electron_mass_c2
Definition: hepunit.py:273

References AverageValue(), G4InuclSpecialFunctions::bindingEnergy(), source.hepunit::electron_mass_c2, G4INCL::KinematicsUtils::energy(), Function(), G4cout, G4endl, iMax, G4AtomicTransitionManager::Instance(), IntSpectrum(), keV, lowestE, G4INCL::Math::max(), MaxEnergyOfSecondaries(), MeV, G4INCL::Math::min(), G4eIonisationParameters::Parameter(), theParam, G4InuclParticleNames::tm, verbose, and Z.

◆ AverageValue()

G4double G4eIonisationSpectrum::AverageValue ( G4double  xMin,
G4double  xMax,
const G4DataVector p 
) const
private

Definition at line 516 of file G4eIonisationSpectrum.cc.

519{
520 G4double sum = 0.0;
521 if(xMin >= xMax) return sum;
522
523 G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2;
524
525 // Integral over interpolation aria
526 if(xMin < p[3]) {
527
528 x1 = p[1];
529 y1 = p[4];
530
531 G4double dx = (p[2] - p[1]) / 3.0;
532 G4double dx1= G4Exp(std::log(p[3]/p[2]) / 16.0);
533
534 for (size_t i=0; i<19; i++) {
535
536 if (i < 3) {
537 x2 = x1 + dx;
538 } else if(18 == i) {
539 x2 = p[3];
540 } else {
541 x2 = x1*dx1;
542 }
543
544 y2 = p[5 + i];
545
546 if (xMax <= x1) {
547 break;
548 } else if (xMin < x2) {
549
550 xs1 = x1;
551 xs2 = x2;
552 ys1 = y1;
553 ys2 = y2;
554
555 if (x2 > x1) {
556 if (xMin > x1) {
557 xs1 = xMin;
558 ys1 += (xs1 - x1)*(y2 - y1)/(x2 - x1);
559 }
560 if (xMax < x2) {
561 xs2 = xMax;
562 ys2 += (xs2 - x2)*(y1 - y2)/(x1 - x2);
563 }
564 if (xs2 > xs1) {
565 sum += std::log(xs2/xs1)*(ys1*xs2 - ys2*xs1)/(xs2 - xs1)
566 + ys2 - ys1;
567 }
568 }
569 }
570 x1 = x2;
571 y1 = y2;
572
573 }
574 }
575
576 // Integral over aria with parametrised formula
577
578 x1 = std::max(xMin, p[3]);
579 if(x1 >= xMax) return sum;
580 x2 = xMax;
581
582 xs1 = 1./x1;
583 xs2 = 1./x2;
584
585 sum += std::log(x2/x1)*(1.0 - p[0])
586 + 0.5*(1. - p[iMax])*(x2*x2 - x1*x1)
587 + 1./(1. - x2) - 1./(1. - x1)
588 + (1. + p[iMax])*std::log((1. - x2)/(1. - x1))
589 + 0.5*p[0]*(xs1 - xs2);
590
591 return sum;
592}
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:179

References G4Exp(), iMax, and G4INCL::Math::max().

Referenced by AverageEnergy().

◆ Excitation()

G4double G4eIonisationSpectrum::Excitation ( G4int  Z,
G4double  e 
) const
inlinevirtual

Implements G4VEnergySpectrum.

Definition at line 128 of file G4eIonisationSpectrum.hh.

129{
130 return theParam->Excitation(Z, e);
131}
G4double Excitation(G4int Z, G4double e) const

References G4eIonisationParameters::Excitation(), theParam, and Z.

◆ Function()

G4double G4eIonisationSpectrum::Function ( G4double  x,
const G4DataVector p 
) const
inlineprivate

Definition at line 118 of file G4eIonisationSpectrum.hh.

120{
121 G4double f = 1.0 - p[0] - p[iMax]*x
122 + x*x*(1.0 - p[iMax] + (1.0/(1.0 - x) - p[iMax])/(1.0 - x))
123 + 0.5*p[0]/x;
124
125 return f;
126}

References iMax.

Referenced by AverageEnergy(), Probability(), and SampleEnergy().

◆ IntSpectrum()

G4double G4eIonisationSpectrum::IntSpectrum ( G4double  xMin,
G4double  xMax,
const G4DataVector p 
) const
private

Definition at line 432 of file G4eIonisationSpectrum.cc.

435{
436 // Please comment what IntSpectrum does
437 G4double sum = 0.0;
438 if(xMin >= xMax) return sum;
439
440 G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2, q;
441
442 // Integral over interpolation aria
443 if(xMin < p[3]) {
444
445 x1 = p[1];
446 y1 = p[4];
447
448 G4double dx = (p[2] - p[1]) / 3.0;
449 G4double dx1= G4Exp(std::log(p[3]/p[2]) / 16.0);
450
451 for (size_t i=0; i<19; i++) {
452
453 q = 0.0;
454 if (i < 3) {
455 x2 = x1 + dx;
456 } else if(18 == i) {
457 x2 = p[3];
458 } else {
459 x2 = x1*dx1;
460 }
461
462 y2 = p[5 + i];
463
464 if (xMax <= x1) {
465 break;
466 } else if (xMin < x2) {
467
468 xs1 = x1;
469 xs2 = x2;
470 ys1 = y1;
471 ys2 = y2;
472
473 if (x2 > x1) {
474 if (xMin > x1) {
475 xs1 = xMin;
476 ys1 += (xs1 - x1)*(y2 - y1)/(x2 - x1);
477 }
478 if (xMax < x2) {
479 xs2 = xMax;
480 ys2 += (xs2 - x2)*(y1 - y2)/(x1 - x2);
481 }
482 if (xs2 > xs1) {
483 q = (ys1*xs2 - ys2*xs1)/(xs1*xs2)
484 + std::log(xs2/xs1)*(ys2 - ys1)/(xs2 - xs1);
485 sum += q;
486 if(p.size() == 26) G4cout << "i= " << i << " q= " << q << " sum= " << sum << G4endl;
487 }
488 }
489 }
490 x1 = x2;
491 y1 = y2;
492 }
493 }
494
495 // Integral over aria with parametrised formula
496
497 x1 = std::max(xMin, p[3]);
498 if(x1 >= xMax) return sum;
499 x2 = xMax;
500
501 xs1 = 1./x1;
502 xs2 = 1./x2;
503 q = (xs1 - xs2)*(1.0 - p[0])
504 - p[iMax]*std::log(x2/x1)
505 + (1. - p[iMax])*(x2 - x1)
506 + 1./(1. - x2) - 1./(1. - x1)
507 + p[iMax]*std::log((1. - x2)/(1. - x1))
508 + 0.25*p[0]*(xs1*xs1 - xs2*xs2);
509 sum += q;
510 if(p.size() == 26) G4cout << "param... q= " << q << " sum= " << sum << G4endl;
511
512 return sum;
513}

References G4cout, G4endl, G4Exp(), iMax, and G4INCL::Math::max().

Referenced by AverageEnergy(), Probability(), and SampleEnergy().

◆ MaxEnergyOfSecondaries()

G4double G4eIonisationSpectrum::MaxEnergyOfSecondaries ( G4double  kineticEnergy,
G4int  Z = 0,
const G4ParticleDefinition pd = 0 
) const
virtual

Implements G4VEnergySpectrum.

Definition at line 600 of file G4eIonisationSpectrum.cc.

603{
604 return 0.5 * kineticEnergy;
605}

Referenced by AverageEnergy(), Probability(), and SampleEnergy().

◆ operator=()

G4eIonisationSpectrum & G4eIonisationSpectrum::operator= ( const G4eIonisationSpectrum right)
private

◆ PrintData()

void G4eIonisationSpectrum::PrintData ( void  ) const
virtual

Implements G4VEnergySpectrum.

Definition at line 595 of file G4eIonisationSpectrum.cc.

596{
598}

References G4eIonisationParameters::PrintData(), and theParam.

◆ Probability()

G4double G4eIonisationSpectrum::Probability ( G4int  Z,
G4double  tMin,
G4double  tMax,
G4double  kineticEnergy,
G4int  shell,
const G4ParticleDefinition pd = 0 
) const
virtual

Implements G4VEnergySpectrum.

Definition at line 77 of file G4eIonisationSpectrum.cc.

83{
84 // Please comment what Probability does and what are the three
85 // functions mentioned below
86 // Describe the algorithms used
87
89 G4double t0 = std::max(tMin, lowestE);
90 G4double tm = std::min(tMax, eMax);
91 if(t0 >= tm) return 0.0;
92
94 Shell(Z, shell)->BindingEnergy();
95
96 if(e <= bindingEnergy) return 0.0;
97
99
100 G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
102
103 if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0.0)) {
104 G4cout << "G4eIonisationSpectrum::Probability: Z= " << Z
105 << "; shell= " << shell
106 << "; E(keV)= " << e/keV
107 << "; Eb(keV)= " << bindingEnergy/keV
108 << "; x1= " << x1
109 << "; x2= " << x2
110 << G4endl;
111
112 }
113
114 G4DataVector p;
115
116 // Access parameters
117 for (G4int i=0; i<iMax; i++)
118 {
119 G4double x = theParam->Parameter(Z, shell, i, e);
120 if(i<4) x /= energy;
121 p.push_back(x);
122 }
123
124 if(p[3] > 0.5) p[3] = 0.5;
125
126 G4double gLocal = energy/electron_mass_c2 + 1.;
127 p.push_back((2.0*gLocal - 1.0)/(gLocal*gLocal));
128
129 //Add protection against division by zero: actually in Function()
130 //parameter p[3] appears in the denominator. Bug report 1042
131 if (p[3] > 0)
132 p[iMax-1] = Function(p[3], p);
133 else
134 {
135 G4cout << "WARNING: G4eIonisationSpectrum::Probability "
136 << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = "
137 << Z << ". Please check and/or update it " << G4endl;
138 }
139
140 if(e >= 1. && e <= 0. && Z == 4) p.push_back(0.0);
141
142
143 G4double val = IntSpectrum(x1, x2, p);
145 G4double nor = IntSpectrum(x0, 0.5, p);
146
147 if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0.0)) {
148 G4cout << "tcut= " << tMin
149 << "; tMax= " << tMax
150 << "; x0= " << x0
151 << "; x1= " << x1
152 << "; x2= " << x2
153 << "; val= " << val
154 << "; nor= " << nor
155 << "; sum= " << p[0]
156 << "; a= " << p[1]
157 << "; b= " << p[2]
158 << "; c= " << p[3]
159 << G4endl;
160 if(shell == 1) G4cout << "============" << G4endl;
161 }
162
163 p.clear();
164
165 if(nor > 0.0) val /= nor;
166 else val = 0.0;
167
168 return val;
169}

References G4InuclSpecialFunctions::bindingEnergy(), source.hepunit::electron_mass_c2, G4INCL::KinematicsUtils::energy(), Function(), G4cout, G4endl, iMax, G4AtomicTransitionManager::Instance(), IntSpectrum(), keV, lowestE, G4INCL::Math::max(), MaxEnergyOfSecondaries(), G4INCL::Math::min(), G4eIonisationParameters::Parameter(), theParam, G4InuclParticleNames::tm, verbose, and Z.

◆ SampleEnergy()

G4double G4eIonisationSpectrum::SampleEnergy ( G4int  Z,
G4double  tMin,
G4double  tMax,
G4double  kineticEnergy,
G4int  shell,
const G4ParticleDefinition pd = 0 
) const
virtual

Implements G4VEnergySpectrum.

Definition at line 264 of file G4eIonisationSpectrum.cc.

270{
271 // Please comment what SampleEnergy does
272 G4double tDelta = 0.0;
273 G4double t0 = std::max(tMin, lowestE);
275 if(t0 > tm) return tDelta;
276
278 Shell(Z, shell)->BindingEnergy();
279
280 if(e <= bindingEnergy) return 0.0;
281
283
284 G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
286 if(x1 >= x2) return tDelta;
287
288 if(verbose > 1) {
289 G4cout << "G4eIonisationSpectrum::SampleEnergy: Z= " << Z
290 << "; shell= " << shell
291 << "; E(keV)= " << e/keV
292 << G4endl;
293 }
294
295 // Access parameters
296 G4DataVector p;
297
298 // Access parameters
299 for (G4int i=0; i<iMax; i++)
300 {
301 G4double x = theParam->Parameter(Z, shell, i, e);
302 if(i<4) x /= energy;
303 p.push_back(x);
304 }
305
306 if(p[3] > 0.5) p[3] = 0.5;
307
308 G4double gLocal3 = energy/electron_mass_c2 + 1.;
309 p.push_back((2.0*gLocal3 - 1.0)/(gLocal3*gLocal3));
310
311
312 //Add protection against division by zero: actually in Function()
313 //parameter p[3] appears in the denominator. Bug report 1042
314 if (p[3] > 0)
315 p[iMax-1] = Function(p[3], p);
316 else
317 {
318 G4cout << "WARNING: G4eIonisationSpectrum::SampleSpectrum "
319 << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = "
320 << Z << ". Please check and/or update it " << G4endl;
321 }
322
323 G4double aria1 = 0.0;
324 G4double a1 = std::max(x1,p[1]);
325 G4double a2 = std::min(x2,p[3]);
326 if(a1 < a2) aria1 = IntSpectrum(a1, a2, p);
327 G4double aria2 = 0.0;
328 G4double a3 = std::max(x1,p[3]);
329 G4double a4 = x2;
330 if(a3 < a4) aria2 = IntSpectrum(a3, a4, p);
331
332 G4double aria = (aria1 + aria2)*G4UniformRand();
333 G4double amaj, fun, q, x, z1, z2, dx, dx1;
334
335 //======= First aria to sample =====
336
337 if(aria <= aria1) {
338
339 amaj = p[4];
340 for (G4int j=5; j<iMax; j++) {
341 if(p[j] > amaj) amaj = p[j];
342 }
343
344 a1 = 1./a1;
345 a2 = 1./a2;
346
347 G4int i;
348 do {
349
350 x = 1./(a2 + G4UniformRand()*(a1 - a2));
351 z1 = p[1];
352 z2 = p[3];
353 dx = (p[2] - p[1]) / 3.0;
354 dx1= G4Exp(std::log(p[3]/p[2]) / 16.0);
355 for (i=4; i<iMax-1; i++) {
356
357 if (i < 7) {
358 z2 = z1 + dx;
359 } else if(iMax-2 == i) {
360 z2 = p[3];
361 break;
362 } else {
363 z2 = z1*dx1;
364 }
365 if(x >= z1 && x <= z2) break;
366 z1 = z2;
367 }
368 fun = p[i] + (x - z1) * (p[i+1] - p[i])/(z2 - z1);
369
370 if(fun > amaj) {
371 G4cout << "WARNING in G4eIonisationSpectrum::SampleEnergy:"
372 << " Majoranta " << amaj
373 << " < " << fun
374 << " in the first aria at x= " << x
375 << G4endl;
376 }
377
378 q = amaj*G4UniformRand();
379
380 } while (q >= fun);
381
382 //======= Second aria to sample =====
383
384 } else {
385
386 amaj = std::max(p[iMax-1], Function(0.5, p)) * factor;
387 a1 = 1./a3;
388 a2 = 1./a4;
389
390 do {
391
392 x = 1./(a2 + G4UniformRand()*(a1 - a2));
393 fun = Function(x, p);
394
395 if(fun > amaj) {
396 G4cout << "WARNING in G4eIonisationSpectrum::SampleEnergy:"
397 << " Majoranta " << amaj
398 << " < " << fun
399 << " in the second aria at x= " << x
400 << G4endl;
401 }
402
403 q = amaj*G4UniformRand();
404
405 } while (q >= fun);
406
407 }
408
409 p.clear();
410
411 tDelta = x*energy - bindingEnergy;
412
413 if(verbose > 1) {
414 G4cout << "tcut(MeV)= " << tMin/MeV
415 << "; tMax(MeV)= " << tMax/MeV
416 << "; x1= " << x1
417 << "; x2= " << x2
418 << "; a1= " << a1
419 << "; a2= " << a2
420 << "; x= " << x
421 << "; be= " << bindingEnergy
422 << "; e= " << e
423 << "; tDelta= " << tDelta
424 << G4endl;
425 }
426
427
428 return tDelta;
429}
#define G4UniformRand()
Definition: Randomize.hh:52

References G4InuclSpecialFunctions::bindingEnergy(), source.hepunit::electron_mass_c2, G4INCL::KinematicsUtils::energy(), factor, Function(), G4cout, G4endl, G4Exp(), G4UniformRand, iMax, G4AtomicTransitionManager::Instance(), IntSpectrum(), keV, lowestE, G4INCL::Math::max(), MaxEnergyOfSecondaries(), MeV, G4INCL::Math::min(), G4eIonisationParameters::Parameter(), theParam, G4InuclParticleNames::tm, verbose, and Z.

Field Documentation

◆ factor

G4double G4eIonisationSpectrum::factor
private

Definition at line 112 of file G4eIonisationSpectrum.hh.

Referenced by SampleEnergy().

◆ iMax

G4int G4eIonisationSpectrum::iMax
private

◆ lowestE

G4double G4eIonisationSpectrum::lowestE
private

Definition at line 111 of file G4eIonisationSpectrum.hh.

Referenced by AverageEnergy(), Probability(), and SampleEnergy().

◆ theParam

G4eIonisationParameters* G4eIonisationSpectrum::theParam
private

◆ verbose

G4int G4eIonisationSpectrum::verbose
private

Definition at line 114 of file G4eIonisationSpectrum.hh.

Referenced by AverageEnergy(), Probability(), and SampleEnergy().


The documentation for this class was generated from the following files: