Geant4  10.03
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: G4eIonisationSpectrum.cc 95950 2016-03-03 10:42:48Z gcosmo $
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 #include "G4Exp.hh"
61 
62 G4eIonisationSpectrum::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 
72 G4eIonisationSpectrum::~G4eIonisationSpectrum()
73 {
74  delete theParam;
75 }
76 
77 
78 G4double 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 
173 G4double 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 
265 G4double 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= G4Exp(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 
433 G4double 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= G4Exp(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 
517 G4double 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= G4Exp(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 
596 void G4eIonisationSpectrum::PrintData() const
597 {
598  theParam->PrintData();
599 }
600 
601 G4double G4eIonisationSpectrum::MaxEnergyOfSecondaries(G4double kineticEnergy,
602  G4int, // Z = 0,
603  const G4ParticleDefinition* ) const
604 {
605  return 0.5 * kineticEnergy;
606 }
G4Exp.hh
G4eIonisationSpectrum::SampleEnergy
G4double SampleEnergy(G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell, const G4ParticleDefinition *pd=0) const
Definition: G4eIonisationSpectrum.cc:265
G4eIonisationParameters::Parameter
G4double Parameter(G4int Z, G4int shellIndex, G4int parameterIndex, G4double e) const
Definition: G4eIonisationParameters.cc:92
G4AtomicShell.hh
G4int
int G4int
Definition: G4Types.hh:78
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
eV
static constexpr double eV
Definition: G4SIunits.hh:215
Randomize.hh
G4PhysicalConstants.hh
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4AtomicTransitionManager.hh
G4eIonisationSpectrum::Probability
G4double Probability(G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell, const G4ParticleDefinition *pd=0) const
Definition: G4eIonisationSpectrum.cc:78
G4eIonisationSpectrum::theParam
G4eIonisationParameters * theParam
Definition: G4eIonisationSpectrum.hh:111
G4eIonisationSpectrum::AverageValue
G4double AverageValue(G4double xMin, G4double xMax, const G4DataVector &p) const
Definition: G4eIonisationSpectrum.cc:517
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition: G4INCLKinematicsUtils.cc:157
G4eIonisationSpectrum::MaxEnergyOfSecondaries
G4double MaxEnergyOfSecondaries(G4double kineticEnergy, G4int Z=0, const G4ParticleDefinition *pd=0) const
Definition: G4eIonisationSpectrum.cc:601
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:117
G4endl
#define G4endl
Definition: G4ios.hh:61
G4eIonisationSpectrum::Function
G4double Function(G4double x, const G4DataVector &p) const
Definition: G4eIonisationSpectrum.hh:119
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
G4eIonisationSpectrum::AverageEnergy
G4double AverageEnergy(G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell, const G4ParticleDefinition *pd=0) const
Definition: G4eIonisationSpectrum.cc:173
G4eIonisationSpectrum.hh
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
G4DataVector.hh
G4InuclSpecialFunctions::bindingEnergy
G4double bindingEnergy(G4int A, G4int Z)
Definition: bindingEnergy.cc:37
G4AtomicTransitionManager::Instance
static G4AtomicTransitionManager * Instance()
Definition: G4AtomicTransitionManager.cc:45
keV
static constexpr double keV
Definition: G4SIunits.hh:216
G4eIonisationSpectrum::IntSpectrum
G4double IntSpectrum(G4double xMin, G4double xMax, const G4DataVector &p) const
Definition: G4eIonisationSpectrum.cc:433