Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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 ([email protected])
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