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