Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4eIonisationSpectrum.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// -------------------------------------------------------------------
28//
29// GEANT4 Class file
30//
31//
32// File name: G4eIonisationSpectrum
33//
34// Author: V.Ivanchenko (Vladimir.Ivanchenko@cern.ch)
35//
36// Creation date: 29 September 2001
37//
38// Modifications:
39// 10.10.2001 MGP Revision to improve code quality and
40// consistency with design
41// 02.11.2001 VI Optimize sampling of energy
42// 29.11.2001 VI New parametrisation
43// 19.04.2002 VI Add protection in case of energy below binding
44// 30.05.2002 VI Update to 24-parameters data
45// 11.07.2002 VI Fix in integration over spectrum
46// 23.03.2009 LP Added protection against division by zero (for
47// faulty database files), Bug Report 1042
48//
49// -------------------------------------------------------------------
50//
51
52#include "G4eIonisationSpectrum.hh"
53#include "G4AtomicTransitionManager.hh"
54#include "G4AtomicShell.hh"
55#include "G4DataVector.hh"
56#include "Randomize.hh"
57#include "G4PhysicalConstants.hh"
58#include "G4SystemOfUnits.hh"
59#include "G4Exp.hh"
60
61G4eIonisationSpectrum::G4eIonisationSpectrum():G4VEnergySpectrum(),
62 lowestE(0.1*eV),
63 factor(1.3),
64 iMax(24),
65 verbose(0)
66{
67 theParam = new G4eIonisationParameters();
68}
69
70
71G4eIonisationSpectrum::~G4eIonisationSpectrum()
72{
73 delete theParam;
74}
75
76
77G4double G4eIonisationSpectrum::Probability(G4int Z,
78 G4double tMin,
79 G4double tMax,
80 G4double e,
81 G4int shell,
82 const G4ParticleDefinition* ) const
83{
84 // Please comment what Probability does and what are the three
85 // functions mentioned below
86 // Describe the algorithms used
87
88 G4double eMax = MaxEnergyOfSecondaries(e);
89 G4double t0 = std::max(tMin, lowestE);
90 G4double tm = std::min(tMax, eMax);
91 if(t0 >= tm) return 0.0;
92
93 G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
94 Shell(Z, shell)->BindingEnergy();
95
96 if(e <= bindingEnergy) return 0.0;
97
98 G4double energy = e + bindingEnergy;
99
100 G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
101 G4double x2 = std::min(0.5,(tm + 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);
144 G4double x0 = (lowestE + bindingEnergy)/energy;
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}
170
171
172G4double G4eIonisationSpectrum::AverageEnergy(G4int Z,
173 G4double tMin,
174 G4double tMax,
175 G4double e,
176 G4int shell,
177 const G4ParticleDefinition* ) const
178{
179 // Please comment what AverageEnergy does and what are the three
180 // functions mentioned below
181 // Describe the algorithms used
182
183 G4double eMax = MaxEnergyOfSecondaries(e);
184 G4double t0 = std::max(tMin, lowestE);
185 G4double tm = std::min(tMax, eMax);
186 if(t0 >= tm) return 0.0;
187
188 G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
189 Shell(Z, shell)->BindingEnergy();
190
191 if(e <= bindingEnergy) return 0.0;
192
193 G4double energy = e + bindingEnergy;
194
195 G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
196 G4double x2 = std::min(0.5,(tm + 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);
236 G4double x0 = (lowestE + bindingEnergy)/energy;
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}
262
263
264G4double G4eIonisationSpectrum::SampleEnergy(G4int Z,
265 G4double tMin,
266 G4double tMax,
267 G4double e,
268 G4int shell,
269 const G4ParticleDefinition* ) const
270{
271 // Please comment what SampleEnergy does
272 G4double tDelta = 0.0;
273 G4double t0 = std::max(tMin, lowestE);
274 G4double tm = std::min(tMax, MaxEnergyOfSecondaries(e));
275 if(t0 > tm) return tDelta;
276
277 G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
278 Shell(Z, shell)->BindingEnergy();
279
280 if(e <= bindingEnergy) return 0.0;
281
282 G4double energy = e + bindingEnergy;
283
284 G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
285 G4double x2 = std::min(0.5,(tm + 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}
430
431
432G4double G4eIonisationSpectrum::IntSpectrum(G4double xMin,
433 G4double xMax,
434 const G4DataVector& p) const
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}
514
515
516G4double G4eIonisationSpectrum::AverageValue(G4double xMin,
517 G4double xMax,
518 const G4DataVector& p) const
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}
593
594
595void G4eIonisationSpectrum::PrintData() const
596{
597 theParam->PrintData();
598}
599
600G4double G4eIonisationSpectrum::MaxEnergyOfSecondaries(G4double kineticEnergy,
601 G4int, // Z = 0,
602 const G4ParticleDefinition* ) const
603{
604 return 0.5 * kineticEnergy;
605}
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
G4double
double G4double
Definition: G4Types.hh:83
G4int
int G4int
Definition: G4Types.hh:85
Z
const G4int Z[17]
Definition: G4WaterStopping.cc:51
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
G4AtomicTransitionManager::Instance
static G4AtomicTransitionManager * Instance()
Definition: G4AtomicTransitionManager.cc:48
G4eIonisationParameters::Parameter
G4double Parameter(G4int Z, G4int shellIndex, G4int parameterIndex, G4double e) const
Definition: G4eIonisationParameters.cc:90
G4eIonisationSpectrum::AverageEnergy
G4double AverageEnergy(G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell, const G4ParticleDefinition *pd=0) const
Definition: G4eIonisationSpectrum.cc:172
G4eIonisationSpectrum::SampleEnergy
G4double SampleEnergy(G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell, const G4ParticleDefinition *pd=0) const
Definition: G4eIonisationSpectrum.cc:264
G4eIonisationSpectrum::MaxEnergyOfSecondaries
G4double MaxEnergyOfSecondaries(G4double kineticEnergy, G4int Z=0, const G4ParticleDefinition *pd=0) const
Definition: G4eIonisationSpectrum.cc:600
G4eIonisationSpectrum::Probability
G4double Probability(G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell, const G4ParticleDefinition *pd=0) const
Definition: G4eIonisationSpectrum.cc:77