Geant4 11.4.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//
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
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) {
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 G4double val = IntSpectrum(x1, x2, p);
141 G4double x0 = (lowestE + bindingEnergy)/energy;
142 G4double nor = IntSpectrum(x0, 0.5, p);
143
144 if (verbose > 1) {
145 G4cout << "tcut= " << tMin
146 << "; tMax= " << tMax
147 << "; x0= " << x0
148 << "; x1= " << x1
149 << "; x2= " << x2
150 << "; val= " << val
151 << "; nor= " << nor
152 << "; sum= " << p[0]
153 << "; a= " << p[1]
154 << "; b= " << p[2]
155 << "; c= " << p[3]
156 << G4endl;
157 if(shell == 1) G4cout << "============" << G4endl;
158 }
159
160 p.clear();
161
162 if(nor > 0.0) val /= nor;
163 else val = 0.0;
164
165 return val;
166}
167
168
169G4double G4eIonisationSpectrum::AverageEnergy(G4int Z,
170 G4double tMin,
171 G4double tMax,
172 G4double e,
173 G4int shell,
174 const G4ParticleDefinition* ) const
175{
176 // Please comment what AverageEnergy does and what are the three
177 // functions mentioned below
178 // Describe the algorithms used
179
180 G4double eMax = MaxEnergyOfSecondaries(e);
181 G4double t0 = std::max(tMin, lowestE);
182 G4double tm = std::min(tMax, eMax);
183 if(t0 >= tm) return 0.0;
184
185 G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
186 Shell(Z, shell)->BindingEnergy();
187
188 if(e <= bindingEnergy) return 0.0;
189
190 G4double energy = e + bindingEnergy;
191
192 G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
193 G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
194
195 if(verbose > 1) {
196 G4cout << "G4eIonisationSpectrum::AverageEnergy: Z= " << Z
197 << "; shell= " << shell
198 << "; E(keV)= " << e/keV
199 << "; bindingE(keV)= " << bindingEnergy/keV
200 << "; x1= " << x1
201 << "; x2= " << x2
202 << G4endl;
203 }
204
205 G4DataVector p;
206
207 // Access parameters
208 for (G4int i=0; i<iMax; i++)
209 {
210 G4double x = theParam->Parameter(Z, shell, i, e);
211 if(i<4) x /= energy;
212 p.push_back(x);
213 }
214
215 if(p[3] > 0.5) p[3] = 0.5;
216
217 G4double gLocal2 = energy/electron_mass_c2 + 1.;
218 p.push_back((2.0*gLocal2 - 1.0)/(gLocal2*gLocal2));
219
220
221 //Add protection against division by zero: actually in Function()
222 //parameter p[3] appears in the denominator. Bug report 1042
223 if (p[3] > 0)
224 p[iMax-1] = Function(p[3], p);
225 else
226 {
227 G4cout << "WARNING: G4eIonisationSpectrum::AverageEnergy "
228 << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = "
229 << Z << ". Please check and/or update it " << G4endl;
230 }
231
232 G4double val = AverageValue(x1, x2, p);
233 G4double x0 = (lowestE + bindingEnergy)/energy;
234 G4double nor = IntSpectrum(x0, 0.5, p);
235 val *= energy;
236
237 if(verbose > 1) {
238 G4cout << "tcut(MeV)= " << tMin/MeV
239 << "; tMax(MeV)= " << tMax/MeV
240 << "; x0= " << x0
241 << "; x1= " << x1
242 << "; x2= " << x2
243 << "; val= " << val
244 << "; nor= " << nor
245 << "; sum= " << p[0]
246 << "; a= " << p[1]
247 << "; b= " << p[2]
248 << "; c= " << p[3]
249 << G4endl;
250 }
251
252 p.clear();
253
254 if(nor > 0.0) val /= nor;
255 else val = 0.0;
256
257 return val;
258}
259
260
261G4double G4eIonisationSpectrum::SampleEnergy(G4int Z,
262 G4double tMin,
263 G4double tMax,
264 G4double e,
265 G4int shell,
266 const G4ParticleDefinition* ) const
267{
268 // Please comment what SampleEnergy does
269 G4double tDelta = 0.0;
270 G4double t0 = std::max(tMin, lowestE);
271 G4double tm = std::min(tMax, MaxEnergyOfSecondaries(e));
272 if(t0 > tm) return tDelta;
273
274 G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
275 Shell(Z, shell)->BindingEnergy();
276
277 if(e <= bindingEnergy) return 0.0;
278
279 G4double energy = e + bindingEnergy;
280
281 G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
282 G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
283 if(x1 >= x2) return tDelta;
284
285 if(verbose > 1) {
286 G4cout << "G4eIonisationSpectrum::SampleEnergy: Z= " << Z
287 << "; shell= " << shell
288 << "; E(keV)= " << e/keV
289 << G4endl;
290 }
291
292 // Access parameters
293 G4DataVector p;
294
295 // Access parameters
296 for (G4int i=0; i<iMax; i++)
297 {
298 G4double x = theParam->Parameter(Z, shell, i, e);
299 if(i<4) x /= energy;
300 p.push_back(x);
301 }
302
303 if(p[3] > 0.5) p[3] = 0.5;
304
305 G4double gLocal3 = energy/electron_mass_c2 + 1.;
306 p.push_back((2.0*gLocal3 - 1.0)/(gLocal3*gLocal3));
307
308
309 //Add protection against division by zero: actually in Function()
310 //parameter p[3] appears in the denominator. Bug report 1042
311 if (p[3] > 0)
312 p[iMax-1] = Function(p[3], p);
313 else
314 {
315 G4cout << "WARNING: G4eIonisationSpectrum::SampleSpectrum "
316 << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = "
317 << Z << ". Please check and/or update it " << G4endl;
318 }
319
320 G4double aria1 = 0.0;
321 G4double a1 = std::max(x1,p[1]);
322 G4double a2 = std::min(x2,p[3]);
323 if(a1 < a2) aria1 = IntSpectrum(a1, a2, p);
324 G4double aria2 = 0.0;
325 G4double a3 = std::max(x1,p[3]);
326 G4double a4 = x2;
327 if(a3 < a4) aria2 = IntSpectrum(a3, a4, p);
328
329 G4double aria = (aria1 + aria2)*G4UniformRand();
330 G4double amaj, fun, q, x, z1, z2, dx, dx1;
331
332 //======= First aria to sample =====
333
334 if(aria <= aria1) {
335
336 amaj = p[4];
337 for (G4int j=5; j<iMax; j++) {
338 if(p[j] > amaj) amaj = p[j];
339 }
340
341 a1 = 1./a1;
342 a2 = 1./a2;
343
344 G4int i;
345 do {
346
347 x = 1./(a2 + G4UniformRand()*(a1 - a2));
348 z1 = p[1];
349 z2 = p[3];
350 dx = (p[2] - p[1]) / 3.0;
351 dx1= G4Exp(std::log(p[3]/p[2]) / 16.0);
352 for (i=4; i<iMax-1; i++) {
353
354 if (i < 7) {
355 z2 = z1 + dx;
356 } else if(iMax-2 == i) {
357 z2 = p[3];
358 break;
359 } else {
360 z2 = z1*dx1;
361 }
362 if(x >= z1 && x <= z2) break;
363 z1 = z2;
364 }
365 fun = p[i] + (x - z1) * (p[i+1] - p[i])/(z2 - z1);
366
367 if(fun > amaj) {
368 G4cout << "WARNING in G4eIonisationSpectrum::SampleEnergy:"
369 << " Majoranta " << amaj
370 << " < " << fun
371 << " in the first aria at x= " << x
372 << G4endl;
373 }
374
375 q = amaj*G4UniformRand();
376
377 } while (q >= fun);
378
379 //======= Second aria to sample =====
380
381 } else {
382
383 amaj = std::max(p[iMax-1], Function(0.5, p)) * factor;
384 a1 = 1./a3;
385 a2 = 1./a4;
386
387 do {
388
389 x = 1./(a2 + G4UniformRand()*(a1 - a2));
390 fun = Function(x, p);
391
392 if(fun > amaj) {
393 G4cout << "WARNING in G4eIonisationSpectrum::SampleEnergy:"
394 << " Majoranta " << amaj
395 << " < " << fun
396 << " in the second aria at x= " << x
397 << G4endl;
398 }
399
400 q = amaj*G4UniformRand();
401
402 } while (q >= fun);
403
404 }
405
406 p.clear();
407
408 tDelta = x*energy - bindingEnergy;
409
410 if(verbose > 1) {
411 G4cout << "tcut(MeV)= " << tMin/MeV
412 << "; tMax(MeV)= " << tMax/MeV
413 << "; x1= " << x1
414 << "; x2= " << x2
415 << "; a1= " << a1
416 << "; a2= " << a2
417 << "; x= " << x
418 << "; be= " << bindingEnergy
419 << "; e= " << e
420 << "; tDelta= " << tDelta
421 << G4endl;
422 }
423
424
425 return tDelta;
426}
427
428
429G4double G4eIonisationSpectrum::IntSpectrum(G4double xMin,
430 G4double xMax,
431 const G4DataVector& p) const
432{
433 // Please comment what IntSpectrum does
434 G4double sum = 0.0;
435 if(xMin >= xMax) return sum;
436
437 G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2, q;
438
439 // Integral over interpolation aria
440 if(xMin < p[3]) {
441
442 x1 = p[1];
443 y1 = p[4];
444
445 G4double dx = (p[2] - p[1]) / 3.0;
446 G4double dx1= G4Exp(std::log(p[3]/p[2]) / 16.0);
447
448 for (size_t i=0; i<19; i++) {
449
450 q = 0.0;
451 if (i < 3) {
452 x2 = x1 + dx;
453 } else if(18 == i) {
454 x2 = p[3];
455 } else {
456 x2 = x1*dx1;
457 }
458
459 y2 = p[5 + i];
460
461 if (xMax <= x1) {
462 break;
463 } else if (xMin < x2) {
464
465 xs1 = x1;
466 xs2 = x2;
467 ys1 = y1;
468 ys2 = y2;
469
470 if (x2 > x1) {
471 if (xMin > x1) {
472 xs1 = xMin;
473 ys1 += (xs1 - x1)*(y2 - y1)/(x2 - x1);
474 }
475 if (xMax < x2) {
476 xs2 = xMax;
477 ys2 += (xs2 - x2)*(y1 - y2)/(x1 - x2);
478 }
479 if (xs2 > xs1) {
480 q = (ys1*xs2 - ys2*xs1)/(xs1*xs2)
481 + std::log(xs2/xs1)*(ys2 - ys1)/(xs2 - xs1);
482 sum += q;
483 if(p.size() == 26) G4cout << "i= " << i << " q= " << q << " sum= " << sum << G4endl;
484 }
485 }
486 }
487 x1 = x2;
488 y1 = y2;
489 }
490 }
491
492 // Integral over aria with parametrised formula
493
494 x1 = std::max(xMin, p[3]);
495 if(x1 >= xMax) return sum;
496 x2 = xMax;
497
498 xs1 = 1./x1;
499 xs2 = 1./x2;
500 q = (xs1 - xs2)*(1.0 - p[0])
501 - p[iMax]*std::log(x2/x1)
502 + (1. - p[iMax])*(x2 - x1)
503 + 1./(1. - x2) - 1./(1. - x1)
504 + p[iMax]*std::log((1. - x2)/(1. - x1))
505 + 0.25*p[0]*(xs1*xs1 - xs2*xs2);
506 sum += q;
507 if(p.size() == 26) G4cout << "param... q= " << q << " sum= " << sum << G4endl;
508
509 return sum;
510}
511
512
513G4double G4eIonisationSpectrum::AverageValue(G4double xMin,
514 G4double xMax,
515 const G4DataVector& p) const
516{
517 G4double sum = 0.0;
518 if(xMin >= xMax) return sum;
519
520 G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2;
521
522 // Integral over interpolation aria
523 if(xMin < p[3]) {
524
525 x1 = p[1];
526 y1 = p[4];
527
528 G4double dx = (p[2] - p[1]) / 3.0;
529 G4double dx1= G4Exp(std::log(p[3]/p[2]) / 16.0);
530
531 for (size_t i=0; i<19; i++) {
532
533 if (i < 3) {
534 x2 = x1 + dx;
535 } else if(18 == i) {
536 x2 = p[3];
537 } else {
538 x2 = x1*dx1;
539 }
540
541 y2 = p[5 + i];
542
543 if (xMax <= x1) {
544 break;
545 } else if (xMin < x2) {
546
547 xs1 = x1;
548 xs2 = x2;
549 ys1 = y1;
550 ys2 = y2;
551
552 if (x2 > x1) {
553 if (xMin > x1) {
554 xs1 = xMin;
555 ys1 += (xs1 - x1)*(y2 - y1)/(x2 - x1);
556 }
557 if (xMax < x2) {
558 xs2 = xMax;
559 ys2 += (xs2 - x2)*(y1 - y2)/(x1 - x2);
560 }
561 if (xs2 > xs1) {
562 sum += std::log(xs2/xs1)*(ys1*xs2 - ys2*xs1)/(xs2 - xs1)
563 + ys2 - ys1;
564 }
565 }
566 }
567 x1 = x2;
568 y1 = y2;
569
570 }
571 }
572
573 // Integral over aria with parametrised formula
574
575 x1 = std::max(xMin, p[3]);
576 if(x1 >= xMax) return sum;
577 x2 = xMax;
578
579 xs1 = 1./x1;
580 xs2 = 1./x2;
581
582 sum += std::log(x2/x1)*(1.0 - p[0])
583 + 0.5*(1. - p[iMax])*(x2*x2 - x1*x1)
584 + 1./(1. - x2) - 1./(1. - x1)
585 + (1. + p[iMax])*std::log((1. - x2)/(1. - x1))
586 + 0.5*p[0]*(xs1 - xs2);
587
588 return sum;
589}
590
591
592void G4eIonisationSpectrum::PrintData() const
593{
594 theParam->PrintData();
595}
596
597G4double G4eIonisationSpectrum::MaxEnergyOfSecondaries(G4double kineticEnergy,
598 G4int, // Z = 0,
599 const G4ParticleDefinition* ) const
600{
601 return 0.5 * kineticEnergy;
602}
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition G4Exp.hh:132
G4double
double G4double
Definition G4Types.hh:83
G4int
int G4int
Definition G4Types.hh:85
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition Randomize.hh:52
G4AtomicTransitionManager::Instance
static G4AtomicTransitionManager * Instance()
Definition G4AtomicTransitionManager.cc:49
G4eIonisationSpectrum::AverageEnergy
G4double AverageEnergy(G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell, const G4ParticleDefinition *pd=0) const
Definition G4eIonisationSpectrum.cc:169
G4eIonisationSpectrum::SampleEnergy
G4double SampleEnergy(G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell, const G4ParticleDefinition *pd=0) const
Definition G4eIonisationSpectrum.cc:261
G4eIonisationSpectrum::MaxEnergyOfSecondaries
G4double MaxEnergyOfSecondaries(G4double kineticEnergy, G4int Z=0, const G4ParticleDefinition *pd=0) const
Definition G4eIonisationSpectrum.cc:597
G4eIonisationSpectrum::Probability
G4double Probability(G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell, const G4ParticleDefinition *pd=0) const
Definition G4eIonisationSpectrum.cc:77