ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4eIonisationSpectrum.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4eIonisationSpectrum.cc
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 
61 G4eIonisationSpectrum::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 
71 G4eIonisationSpectrum::~G4eIonisationSpectrum()
72 {
73  delete theParam;
74 }
75 
76 
77 G4double 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 
172 G4double 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 
264 G4double 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 
432 G4double 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 
516 G4double 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 
595 void G4eIonisationSpectrum::PrintData() const
596 {
597  theParam->PrintData();
598 }
599 
600 G4double G4eIonisationSpectrum::MaxEnergyOfSecondaries(G4double kineticEnergy,
601  G4int, // Z = 0,
602  const G4ParticleDefinition* ) const
603 {
604  return 0.5 * kineticEnergy;
605 }