Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
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// $Id$
27//
28// -------------------------------------------------------------------
29//
30// GEANT4 Class file
31//
32//
33// File name: G4eIonisationSpectrum
34//
35// Author: V.Ivanchenko (Vladimir.Ivanchenko@cern.ch)
36//
37// Creation date: 29 September 2001
38//
39// Modifications:
40// 10.10.2001 MGP Revision to improve code quality and
41// consistency with design
42// 02.11.2001 VI Optimize sampling of energy
43// 29.11.2001 VI New parametrisation
44// 19.04.2002 VI Add protection in case of energy below binding
45// 30.05.2002 VI Update to 24-parameters data
46// 11.07.2002 VI Fix in integration over spectrum
47// 23.03.2009 LP Added protection against division by zero (for
48// faulty database files), Bug Report 1042
49//
50// -------------------------------------------------------------------
51//
52
53#include "G4eIonisationSpectrum.hh"
54#include "G4AtomicTransitionManager.hh"
55#include "G4AtomicShell.hh"
56#include "G4DataVector.hh"
57#include "Randomize.hh"
58#include "G4PhysicalConstants.hh"
59#include "G4SystemOfUnits.hh"
60
61
62G4eIonisationSpectrum::G4eIonisationSpectrum():G4VEnergySpectrum(),
63 lowestE(0.1*eV),
64 factor(1.3),
65 iMax(24),
66 verbose(0)
67{
68 theParam = new G4eIonisationParameters();
69}
70
71
72G4eIonisationSpectrum::~G4eIonisationSpectrum()
73{
74 delete theParam;
75}
76
77
78G4double G4eIonisationSpectrum::Probability(G4int Z,
79 G4double tMin,
80 G4double tMax,
81 G4double e,
82 G4int shell,
83 const G4ParticleDefinition* ) const
84{
85 // Please comment what Probability does and what are the three
86 // functions mentioned below
87 // Describe the algorithms used
88
89 G4double eMax = MaxEnergyOfSecondaries(e);
90 G4double t0 = std::max(tMin, lowestE);
91 G4double tm = std::min(tMax, eMax);
92 if(t0 >= tm) return 0.0;
93
94 G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
95 Shell(Z, shell)->BindingEnergy();
96
97 if(e <= bindingEnergy) return 0.0;
98
99 G4double energy = e + bindingEnergy;
100
101 G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
102 G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
103
104 if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0.0)) {
105 G4cout << "G4eIonisationSpectrum::Probability: Z= " << Z
106 << "; shell= " << shell
107 << "; E(keV)= " << e/keV
108 << "; Eb(keV)= " << bindingEnergy/keV
109 << "; x1= " << x1
110 << "; x2= " << x2
111 << G4endl;
112
113 }
114
115 G4DataVector p;
116
117 // Access parameters
118 for (G4int i=0; i<iMax; i++)
119 {
120 G4double x = theParam->Parameter(Z, shell, i, e);
121 if(i<4) x /= energy;
122 p.push_back(x);
123 }
124
125 if(p[3] > 0.5) p[3] = 0.5;
126
127 G4double gLocal = energy/electron_mass_c2 + 1.;
128 p.push_back((2.0*gLocal - 1.0)/(gLocal*gLocal));
129
130 //Add protection against division by zero: actually in Function()
131 //parameter p[3] appears in the denominator. Bug report 1042
132 if (p[3] > 0)
133 p[iMax-1] = Function(p[3], p);
134 else
135 {
136 G4cout << "WARNING: G4eIonisationSpectrum::Probability "
137 << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = "
138 << Z << ". Please check and/or update it " << G4endl;
139 }
140
141 if(e >= 1. && e <= 0. && Z == 4) p.push_back(0.0);
142
143
144 G4double val = IntSpectrum(x1, x2, p);
145 G4double x0 = (lowestE + bindingEnergy)/energy;
146 G4double nor = IntSpectrum(x0, 0.5, p);
147
148 if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0.0)) {
149 G4cout << "tcut= " << tMin
150 << "; tMax= " << tMax
151 << "; x0= " << x0
152 << "; x1= " << x1
153 << "; x2= " << x2
154 << "; val= " << val
155 << "; nor= " << nor
156 << "; sum= " << p[0]
157 << "; a= " << p[1]
158 << "; b= " << p[2]
159 << "; c= " << p[3]
160 << G4endl;
161 if(shell == 1) G4cout << "============" << G4endl;
162 }
163
164 p.clear();
165
166 if(nor > 0.0) val /= nor;
167 else val = 0.0;
168
169 return val;
170}
171
172
173G4double G4eIonisationSpectrum::AverageEnergy(G4int Z,
174 G4double tMin,
175 G4double tMax,
176 G4double e,
177 G4int shell,
178 const G4ParticleDefinition* ) const
179{
180 // Please comment what AverageEnergy does and what are the three
181 // functions mentioned below
182 // Describe the algorithms used
183
184 G4double eMax = MaxEnergyOfSecondaries(e);
185 G4double t0 = std::max(tMin, lowestE);
186 G4double tm = std::min(tMax, eMax);
187 if(t0 >= tm) return 0.0;
188
189 G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
190 Shell(Z, shell)->BindingEnergy();
191
192 if(e <= bindingEnergy) return 0.0;
193
194 G4double energy = e + bindingEnergy;
195
196 G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
197 G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
198
199 if(verbose > 1) {
200 G4cout << "G4eIonisationSpectrum::AverageEnergy: Z= " << Z
201 << "; shell= " << shell
202 << "; E(keV)= " << e/keV
203 << "; bindingE(keV)= " << bindingEnergy/keV
204 << "; x1= " << x1
205 << "; x2= " << x2
206 << G4endl;
207 }
208
209 G4DataVector p;
210
211 // Access parameters
212 for (G4int i=0; i<iMax; i++)
213 {
214 G4double x = theParam->Parameter(Z, shell, i, e);
215 if(i<4) x /= energy;
216 p.push_back(x);
217 }
218
219 if(p[3] > 0.5) p[3] = 0.5;
220
221 G4double gLocal2 = energy/electron_mass_c2 + 1.;
222 p.push_back((2.0*gLocal2 - 1.0)/(gLocal2*gLocal2));
223
224
225 //Add protection against division by zero: actually in Function()
226 //parameter p[3] appears in the denominator. Bug report 1042
227 if (p[3] > 0)
228 p[iMax-1] = Function(p[3], p);
229 else
230 {
231 G4cout << "WARNING: G4eIonisationSpectrum::AverageEnergy "
232 << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = "
233 << Z << ". Please check and/or update it " << G4endl;
234 }
235
236 G4double val = AverageValue(x1, x2, p);
237 G4double x0 = (lowestE + bindingEnergy)/energy;
238 G4double nor = IntSpectrum(x0, 0.5, p);
239 val *= energy;
240
241 if(verbose > 1) {
242 G4cout << "tcut(MeV)= " << tMin/MeV
243 << "; tMax(MeV)= " << tMax/MeV
244 << "; x0= " << x0
245 << "; x1= " << x1
246 << "; x2= " << x2
247 << "; val= " << val
248 << "; nor= " << nor
249 << "; sum= " << p[0]
250 << "; a= " << p[1]
251 << "; b= " << p[2]
252 << "; c= " << p[3]
253 << G4endl;
254 }
255
256 p.clear();
257
258 if(nor > 0.0) val /= nor;
259 else val = 0.0;
260
261 return val;
262}
263
264
265G4double G4eIonisationSpectrum::SampleEnergy(G4int Z,
266 G4double tMin,
267 G4double tMax,
268 G4double e,
269 G4int shell,
270 const G4ParticleDefinition* ) const
271{
272 // Please comment what SampleEnergy does
273 G4double tDelta = 0.0;
274 G4double t0 = std::max(tMin, lowestE);
275 G4double tm = std::min(tMax, MaxEnergyOfSecondaries(e));
276 if(t0 > tm) return tDelta;
277
278 G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
279 Shell(Z, shell)->BindingEnergy();
280
281 if(e <= bindingEnergy) return 0.0;
282
283 G4double energy = e + bindingEnergy;
284
285 G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
286 G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
287 if(x1 >= x2) return tDelta;
288
289 if(verbose > 1) {
290 G4cout << "G4eIonisationSpectrum::SampleEnergy: Z= " << Z
291 << "; shell= " << shell
292 << "; E(keV)= " << e/keV
293 << G4endl;
294 }
295
296 // Access parameters
297 G4DataVector p;
298
299 // Access parameters
300 for (G4int i=0; i<iMax; i++)
301 {
302 G4double x = theParam->Parameter(Z, shell, i, e);
303 if(i<4) x /= energy;
304 p.push_back(x);
305 }
306
307 if(p[3] > 0.5) p[3] = 0.5;
308
309 G4double gLocal3 = energy/electron_mass_c2 + 1.;
310 p.push_back((2.0*gLocal3 - 1.0)/(gLocal3*gLocal3));
311
312
313 //Add protection against division by zero: actually in Function()
314 //parameter p[3] appears in the denominator. Bug report 1042
315 if (p[3] > 0)
316 p[iMax-1] = Function(p[3], p);
317 else
318 {
319 G4cout << "WARNING: G4eIonisationSpectrum::SampleSpectrum "
320 << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = "
321 << Z << ". Please check and/or update it " << G4endl;
322 }
323
324 G4double aria1 = 0.0;
325 G4double a1 = std::max(x1,p[1]);
326 G4double a2 = std::min(x2,p[3]);
327 if(a1 < a2) aria1 = IntSpectrum(a1, a2, p);
328 G4double aria2 = 0.0;
329 G4double a3 = std::max(x1,p[3]);
330 G4double a4 = x2;
331 if(a3 < a4) aria2 = IntSpectrum(a3, a4, p);
332
333 G4double aria = (aria1 + aria2)*G4UniformRand();
334 G4double amaj, fun, q, x, z1, z2, dx, dx1;
335
336 //======= First aria to sample =====
337
338 if(aria <= aria1) {
339
340 amaj = p[4];
341 for (G4int j=5; j<iMax; j++) {
342 if(p[j] > amaj) amaj = p[j];
343 }
344
345 a1 = 1./a1;
346 a2 = 1./a2;
347
348 G4int i;
349 do {
350
351 x = 1./(a2 + G4UniformRand()*(a1 - a2));
352 z1 = p[1];
353 z2 = p[3];
354 dx = (p[2] - p[1]) / 3.0;
355 dx1= std::exp(std::log(p[3]/p[2]) / 16.0);
356 for (i=4; i<iMax-1; i++) {
357
358 if (i < 7) {
359 z2 = z1 + dx;
360 } else if(iMax-2 == i) {
361 z2 = p[3];
362 break;
363 } else {
364 z2 = z1*dx1;
365 }
366 if(x >= z1 && x <= z2) break;
367 z1 = z2;
368 }
369 fun = p[i] + (x - z1) * (p[i+1] - p[i])/(z2 - z1);
370
371 if(fun > amaj) {
372 G4cout << "WARNING in G4eIonisationSpectrum::SampleEnergy:"
373 << " Majoranta " << amaj
374 << " < " << fun
375 << " in the first aria at x= " << x
376 << G4endl;
377 }
378
379 q = amaj*G4UniformRand();
380
381 } while (q >= fun);
382
383 //======= Second aria to sample =====
384
385 } else {
386
387 amaj = std::max(p[iMax-1], Function(0.5, p)) * factor;
388 a1 = 1./a3;
389 a2 = 1./a4;
390
391 do {
392
393 x = 1./(a2 + G4UniformRand()*(a1 - a2));
394 fun = Function(x, p);
395
396 if(fun > amaj) {
397 G4cout << "WARNING in G4eIonisationSpectrum::SampleEnergy:"
398 << " Majoranta " << amaj
399 << " < " << fun
400 << " in the second aria at x= " << x
401 << G4endl;
402 }
403
404 q = amaj*G4UniformRand();
405
406 } while (q >= fun);
407
408 }
409
410 p.clear();
411
412 tDelta = x*energy - bindingEnergy;
413
414 if(verbose > 1) {
415 G4cout << "tcut(MeV)= " << tMin/MeV
416 << "; tMax(MeV)= " << tMax/MeV
417 << "; x1= " << x1
418 << "; x2= " << x2
419 << "; a1= " << a1
420 << "; a2= " << a2
421 << "; x= " << x
422 << "; be= " << bindingEnergy
423 << "; e= " << e
424 << "; tDelta= " << tDelta
425 << G4endl;
426 }
427
428
429 return tDelta;
430}
431
432
433G4double G4eIonisationSpectrum::IntSpectrum(G4double xMin,
434 G4double xMax,
435 const G4DataVector& p) const
436{
437 // Please comment what IntSpectrum does
438 G4double sum = 0.0;
439 if(xMin >= xMax) return sum;
440
441 G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2, q;
442
443 // Integral over interpolation aria
444 if(xMin < p[3]) {
445
446 x1 = p[1];
447 y1 = p[4];
448
449 G4double dx = (p[2] - p[1]) / 3.0;
450 G4double dx1= std::exp(std::log(p[3]/p[2]) / 16.0);
451
452 for (size_t i=0; i<19; i++) {
453
454 q = 0.0;
455 if (i < 3) {
456 x2 = x1 + dx;
457 } else if(18 == i) {
458 x2 = p[3];
459 } else {
460 x2 = x1*dx1;
461 }
462
463 y2 = p[5 + i];
464
465 if (xMax <= x1) {
466 break;
467 } else if (xMin < x2) {
468
469 xs1 = x1;
470 xs2 = x2;
471 ys1 = y1;
472 ys2 = y2;
473
474 if (x2 > x1) {
475 if (xMin > x1) {
476 xs1 = xMin;
477 ys1 += (xs1 - x1)*(y2 - y1)/(x2 - x1);
478 }
479 if (xMax < x2) {
480 xs2 = xMax;
481 ys2 += (xs2 - x2)*(y1 - y2)/(x1 - x2);
482 }
483 if (xs2 > xs1) {
484 q = (ys1*xs2 - ys2*xs1)/(xs1*xs2)
485 + std::log(xs2/xs1)*(ys2 - ys1)/(xs2 - xs1);
486 sum += q;
487 if(p.size() == 26) G4cout << "i= " << i << " q= " << q << " sum= " << sum << G4endl;
488 }
489 }
490 }
491 x1 = x2;
492 y1 = y2;
493 }
494 }
495
496 // Integral over aria with parametrised formula
497
498 x1 = std::max(xMin, p[3]);
499 if(x1 >= xMax) return sum;
500 x2 = xMax;
501
502 xs1 = 1./x1;
503 xs2 = 1./x2;
504 q = (xs1 - xs2)*(1.0 - p[0])
505 - p[iMax]*std::log(x2/x1)
506 + (1. - p[iMax])*(x2 - x1)
507 + 1./(1. - x2) - 1./(1. - x1)
508 + p[iMax]*std::log((1. - x2)/(1. - x1))
509 + 0.25*p[0]*(xs1*xs1 - xs2*xs2);
510 sum += q;
511 if(p.size() == 26) G4cout << "param... q= " << q << " sum= " << sum << G4endl;
512
513 return sum;
514}
515
516
517G4double G4eIonisationSpectrum::AverageValue(G4double xMin,
518 G4double xMax,
519 const G4DataVector& p) const
520{
521 G4double sum = 0.0;
522 if(xMin >= xMax) return sum;
523
524 G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2;
525
526 // Integral over interpolation aria
527 if(xMin < p[3]) {
528
529 x1 = p[1];
530 y1 = p[4];
531
532 G4double dx = (p[2] - p[1]) / 3.0;
533 G4double dx1= std::exp(std::log(p[3]/p[2]) / 16.0);
534
535 for (size_t i=0; i<19; i++) {
536
537 if (i < 3) {
538 x2 = x1 + dx;
539 } else if(18 == i) {
540 x2 = p[3];
541 } else {
542 x2 = x1*dx1;
543 }
544
545 y2 = p[5 + i];
546
547 if (xMax <= x1) {
548 break;
549 } else if (xMin < x2) {
550
551 xs1 = x1;
552 xs2 = x2;
553 ys1 = y1;
554 ys2 = y2;
555
556 if (x2 > x1) {
557 if (xMin > x1) {
558 xs1 = xMin;
559 ys1 += (xs1 - x1)*(y2 - y1)/(x2 - x1);
560 }
561 if (xMax < x2) {
562 xs2 = xMax;
563 ys2 += (xs2 - x2)*(y1 - y2)/(x1 - x2);
564 }
565 if (xs2 > xs1) {
566 sum += std::log(xs2/xs1)*(ys1*xs2 - ys2*xs1)/(xs2 - xs1)
567 + ys2 - ys1;
568 }
569 }
570 }
571 x1 = x2;
572 y1 = y2;
573
574 }
575 }
576
577 // Integral over aria with parametrised formula
578
579 x1 = std::max(xMin, p[3]);
580 if(x1 >= xMax) return sum;
581 x2 = xMax;
582
583 xs1 = 1./x1;
584 xs2 = 1./x2;
585
586 sum += std::log(x2/x1)*(1.0 - p[0])
587 + 0.5*(1. - p[iMax])*(x2*x2 - x1*x1)
588 + 1./(1. - x2) - 1./(1. - x1)
589 + (1. + p[iMax])*std::log((1. - x2)/(1. - x1))
590 + 0.5*p[0]*(xs1 - xs2);
591
592 return sum;
593}
594
595
596void G4eIonisationSpectrum::PrintData() const
597{
598 theParam->PrintData();
599}
600
601G4double G4eIonisationSpectrum::MaxEnergyOfSecondaries(G4double kineticEnergy,
602 G4int, // Z = 0,
603 const G4ParticleDefinition* ) const
604{
605 return 0.5 * kineticEnergy;
606}
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
G4AtomicTransitionManager::Instance
static G4AtomicTransitionManager * Instance()
Definition: G4AtomicTransitionManager.cc:163
G4eIonisationParameters::Parameter
G4double Parameter(G4int Z, G4int shellIndex, G4int parameterIndex, G4double e) const
Definition: G4eIonisationParameters.cc:92
G4eIonisationSpectrum::AverageEnergy
G4double AverageEnergy(G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell, const G4ParticleDefinition *pd=0) const
Definition: G4eIonisationSpectrum.cc:173
G4eIonisationSpectrum::SampleEnergy
G4double SampleEnergy(G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell, const G4ParticleDefinition *pd=0) const
Definition: G4eIonisationSpectrum.cc:265
G4eIonisationSpectrum::MaxEnergyOfSecondaries
G4double MaxEnergyOfSecondaries(G4double kineticEnergy, G4int Z=0, const G4ParticleDefinition *pd=0) const
Definition: G4eIonisationSpectrum.cc:601
G4eIonisationSpectrum::Probability
G4double Probability(G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell, const G4ParticleDefinition *pd=0) const
Definition: G4eIonisationSpectrum.cc:78