Geant4  10.01.p02
G4RDeIonisationSpectrum.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 // GEANT4 tag $Name: geant4-09-01-ref-00 $
28 //
29 // -------------------------------------------------------------------
30 //
31 // GEANT4 Class file
32 //
33 //
34 // File name: G4RDeIonisationSpectrum
35 //
36 // Author: V.Ivanchenko (Vladimir.Ivanchenko@cern.ch)
37 //
38 // Creation date: 29 September 2001
39 //
40 // Modifications:
41 // 10.10.2001 MGP Revision to improve code quality and
42 // consistency with design
43 // 02.11.2001 VI Optimize sampling of energy
44 // 29.11.2001 VI New parametrisation
45 // 19.04.2002 VI Add protection in case of energy below binding
46 // 30.05.2002 VI Update to 24-parameters data
47 // 11.07.2002 VI Fix in integration over spectrum
48 //
49 // -------------------------------------------------------------------
50 //
51 
52 #include "G4RDeIonisationSpectrum.hh"
53 #include "G4RDAtomicTransitionManager.hh"
54 #include "G4RDAtomicShell.hh"
55 #include "G4PhysicalConstants.hh"
56 #include "G4SystemOfUnits.hh"
57 #include "G4DataVector.hh"
58 #include "Randomize.hh"
59 
60 
61 G4RDeIonisationSpectrum::G4RDeIonisationSpectrum():G4RDVEnergySpectrum(),
62  lowestE(0.1*eV),
63  factor(1.3),
64  iMax(24),
65  verbose(0)
66 {
67  theParam = new G4RDeIonisationParameters();
68 }
69 
70 
71 G4RDeIonisationSpectrum::~G4RDeIonisationSpectrum()
72 {
73  delete theParam;
74 }
75 
76 
77 G4double G4RDeIonisationSpectrum::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 = (G4RDAtomicTransitionManager::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 << "G4RDeIonisationSpectrum::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 g = energy/electron_mass_c2 + 1.;
127  p.push_back((2.0*g - 1.0)/(g*g));
128 
129  p[iMax-1] = Function(p[3], p);
130 
131  if(e >= 1. && e <= 0. && Z == 4) p.push_back(0.0);
132 
133 
134  G4double val = IntSpectrum(x1, x2, p);
135  G4double x0 = (lowestE + bindingEnergy)/energy;
136  G4double nor = IntSpectrum(x0, 0.5, p);
137 
138  if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0.0)) {
139  G4cout << "tcut= " << tMin
140  << "; tMax= " << tMax
141  << "; x0= " << x0
142  << "; x1= " << x1
143  << "; x2= " << x2
144  << "; val= " << val
145  << "; nor= " << nor
146  << "; sum= " << p[0]
147  << "; a= " << p[1]
148  << "; b= " << p[2]
149  << "; c= " << p[3]
150  << G4endl;
151  if(shell == 1) G4cout << "============" << G4endl;
152  }
153 
154  p.clear();
155 
156  if(nor > 0.0) val /= nor;
157  else val = 0.0;
158 
159  return val;
160 }
161 
162 
163 G4double G4RDeIonisationSpectrum::AverageEnergy(G4int Z,
164  G4double tMin,
165  G4double tMax,
166  G4double e,
167  G4int shell,
168  const G4ParticleDefinition* ) const
169 {
170  // Please comment what AverageEnergy does and what are the three
171  // functions mentioned below
172  // Describe the algorithms used
173 
174  G4double eMax = MaxEnergyOfSecondaries(e);
175  G4double t0 = std::max(tMin, lowestE);
176  G4double tm = std::min(tMax, eMax);
177  if(t0 >= tm) return 0.0;
178 
179  G4double bindingEnergy = (G4RDAtomicTransitionManager::Instance())->
180  Shell(Z, shell)->BindingEnergy();
181 
182  if(e <= bindingEnergy) return 0.0;
183 
184  G4double energy = e + bindingEnergy;
185 
186  G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
187  G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
188 
189  if(verbose > 1) {
190  G4cout << "G4RDeIonisationSpectrum::AverageEnergy: Z= " << Z
191  << "; shell= " << shell
192  << "; E(keV)= " << e/keV
193  << "; bindingE(keV)= " << bindingEnergy/keV
194  << "; x1= " << x1
195  << "; x2= " << x2
196  << G4endl;
197  }
198 
199  G4DataVector p;
200 
201  // Access parameters
202  for (G4int i=0; i<iMax; i++)
203  {
204  G4double x = theParam->Parameter(Z, shell, i, e);
205  if(i<4) x /= energy;
206  p.push_back(x);
207  }
208 
209  if(p[3] > 0.5) p[3] = 0.5;
210 
211  G4double g = energy/electron_mass_c2 + 1.;
212  p.push_back((2.0*g - 1.0)/(g*g));
213 
214  p[iMax-1] = Function(p[3], p);
215 
216  G4double val = AverageValue(x1, x2, p);
217  G4double x0 = (lowestE + bindingEnergy)/energy;
218  G4double nor = IntSpectrum(x0, 0.5, p);
219  val *= energy;
220 
221  if(verbose > 1) {
222  G4cout << "tcut(MeV)= " << tMin/MeV
223  << "; tMax(MeV)= " << tMax/MeV
224  << "; x0= " << x0
225  << "; x1= " << x1
226  << "; x2= " << x2
227  << "; val= " << val
228  << "; nor= " << nor
229  << "; sum= " << p[0]
230  << "; a= " << p[1]
231  << "; b= " << p[2]
232  << "; c= " << p[3]
233  << G4endl;
234  }
235 
236  p.clear();
237 
238  if(nor > 0.0) val /= nor;
239  else val = 0.0;
240 
241  return val;
242 }
243 
244 
245 G4double G4RDeIonisationSpectrum::SampleEnergy(G4int Z,
246  G4double tMin,
247  G4double tMax,
248  G4double e,
249  G4int shell,
250  const G4ParticleDefinition* ) const
251 {
252  // Please comment what SampleEnergy does
253  G4double tDelta = 0.0;
254  G4double t0 = std::max(tMin, lowestE);
255  G4double tm = std::min(tMax, MaxEnergyOfSecondaries(e));
256  if(t0 > tm) return tDelta;
257 
258  G4double bindingEnergy = (G4RDAtomicTransitionManager::Instance())->
259  Shell(Z, shell)->BindingEnergy();
260 
261  if(e <= bindingEnergy) return 0.0;
262 
263  G4double energy = e + bindingEnergy;
264 
265  G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
266  G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
267  if(x1 >= x2) return tDelta;
268 
269  if(verbose > 1) {
270  G4cout << "G4RDeIonisationSpectrum::SampleEnergy: Z= " << Z
271  << "; shell= " << shell
272  << "; E(keV)= " << e/keV
273  << G4endl;
274  }
275 
276  // Access parameters
277  G4DataVector p;
278 
279  // Access parameters
280  for (G4int i=0; i<iMax; i++)
281  {
282  G4double x = theParam->Parameter(Z, shell, i, e);
283  if(i<4) x /= energy;
284  p.push_back(x);
285  }
286 
287  if(p[3] > 0.5) p[3] = 0.5;
288 
289  G4double g = energy/electron_mass_c2 + 1.;
290  p.push_back((2.0*g - 1.0)/(g*g));
291 
292  p[iMax-1] = Function(p[3], p);
293 
294  G4double aria1 = 0.0;
295  G4double a1 = std::max(x1,p[1]);
296  G4double a2 = std::min(x2,p[3]);
297  if(a1 < a2) aria1 = IntSpectrum(a1, a2, p);
298  G4double aria2 = 0.0;
299  G4double a3 = std::max(x1,p[3]);
300  G4double a4 = x2;
301  if(a3 < a4) aria2 = IntSpectrum(a3, a4, p);
302 
303  G4double aria = (aria1 + aria2)*G4UniformRand();
304  G4double amaj, fun, q, x, z1, z2, dx, dx1;
305 
306  //======= First aria to sample =====
307 
308  if(aria <= aria1) {
309 
310  amaj = p[4];
311  for (G4int j=5; j<iMax; j++) {
312  if(p[j] > amaj) amaj = p[j];
313  }
314 
315  a1 = 1./a1;
316  a2 = 1./a2;
317 
318  G4int i;
319  do {
320 
321  x = 1./(a2 + G4UniformRand()*(a1 - a2));
322  z1 = p[1];
323  z2 = p[3];
324  dx = (p[2] - p[1]) / 3.0;
325  dx1= std::exp(std::log(p[3]/p[2]) / 16.0);
326  for (i=4; i<iMax-1; i++) {
327 
328  if (i < 7) {
329  z2 = z1 + dx;
330  } else if(iMax-2 == i) {
331  z2 = p[3];
332  break;
333  } else {
334  z2 = z1*dx1;
335  }
336  if(x >= z1 && x <= z2) break;
337  z1 = z2;
338  }
339  fun = p[i] + (x - z1) * (p[i+1] - p[i])/(z2 - z1);
340 
341  if(fun > amaj) {
342  G4cout << "WARNING in G4RDeIonisationSpectrum::SampleEnergy:"
343  << " Majoranta " << amaj
344  << " < " << fun
345  << " in the first aria at x= " << x
346  << G4endl;
347  }
348 
349  q = amaj*G4UniformRand();
350 
351  } while (q >= fun);
352 
353  //======= Second aria to sample =====
354 
355  } else {
356 
357  amaj = std::max(p[iMax-1], Function(0.5, p)) * factor;
358  a1 = 1./a3;
359  a2 = 1./a4;
360 
361  do {
362 
363  x = 1./(a2 + G4UniformRand()*(a1 - a2));
364  fun = Function(x, p);
365 
366  if(fun > amaj) {
367  G4cout << "WARNING in G4RDeIonisationSpectrum::SampleEnergy:"
368  << " Majoranta " << amaj
369  << " < " << fun
370  << " in the second aria at x= " << x
371  << G4endl;
372  }
373 
374  q = amaj*G4UniformRand();
375 
376  } while (q >= fun);
377 
378  }
379 
380  p.clear();
381 
382  tDelta = x*energy - bindingEnergy;
383 
384  if(verbose > 1) {
385  G4cout << "tcut(MeV)= " << tMin/MeV
386  << "; tMax(MeV)= " << tMax/MeV
387  << "; x1= " << x1
388  << "; x2= " << x2
389  << "; a1= " << a1
390  << "; a2= " << a2
391  << "; x= " << x
392  << "; be= " << bindingEnergy
393  << "; e= " << e
394  << "; tDelta= " << tDelta
395  << G4endl;
396  }
397 
398 
399  return tDelta;
400 }
401 
402 
403 G4double G4RDeIonisationSpectrum::IntSpectrum(G4double xMin,
404  G4double xMax,
405  const G4DataVector& p) const
406 {
407  // Please comment what IntSpectrum does
408  G4double sum = 0.0;
409  if(xMin >= xMax) return sum;
410 
411  G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2, q;
412 
413  // Integral over interpolation aria
414  if(xMin < p[3]) {
415 
416  x1 = p[1];
417  y1 = p[4];
418 
419  G4double dx = (p[2] - p[1]) / 3.0;
420  G4double dx1= std::exp(std::log(p[3]/p[2]) / 16.0);
421 
422  for (size_t i=0; i<19; i++) {
423 
424  q = 0.0;
425  if (i < 3) {
426  x2 = x1 + dx;
427  } else if(18 == i) {
428  x2 = p[3];
429  } else {
430  x2 = x1*dx1;
431  }
432 
433  y2 = p[5 + i];
434 
435  if (xMax <= x1) {
436  break;
437  } else if (xMin < x2) {
438 
439  xs1 = x1;
440  xs2 = x2;
441  ys1 = y1;
442  ys2 = y2;
443 
444  if (x2 > x1) {
445  if (xMin > x1) {
446  xs1 = xMin;
447  ys1 += (xs1 - x1)*(y2 - y1)/(x2 - x1);
448  }
449  if (xMax < x2) {
450  xs2 = xMax;
451  ys2 += (xs2 - x2)*(y1 - y2)/(x1 - x2);
452  }
453  if (xs2 > xs1) {
454  q = (ys1*xs2 - ys2*xs1)/(xs1*xs2)
455  + std::log(xs2/xs1)*(ys2 - ys1)/(xs2 - xs1);
456  sum += q;
457  if(p.size() == 26) G4cout << "i= " << i << " q= " << q << " sum= " << sum << G4endl;
458  }
459  }
460  }
461  x1 = x2;
462  y1 = y2;
463  }
464  }
465 
466  // Integral over aria with parametrised formula
467 
468  x1 = std::max(xMin, p[3]);
469  if(x1 >= xMax) return sum;
470  x2 = xMax;
471 
472  xs1 = 1./x1;
473  xs2 = 1./x2;
474  q = (xs1 - xs2)*(1.0 - p[0])
475  - p[iMax]*std::log(x2/x1)
476  + (1. - p[iMax])*(x2 - x1)
477  + 1./(1. - x2) - 1./(1. - x1)
478  + p[iMax]*std::log((1. - x2)/(1. - x1))
479  + 0.25*p[0]*(xs1*xs1 - xs2*xs2);
480  sum += q;
481  if(p.size() == 26) G4cout << "param... q= " << q << " sum= " << sum << G4endl;
482 
483  return sum;
484 }
485 
486 
487 G4double G4RDeIonisationSpectrum::AverageValue(G4double xMin,
488  G4double xMax,
489  const G4DataVector& p) const
490 {
491  G4double sum = 0.0;
492  if(xMin >= xMax) return sum;
493 
494  G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2;
495 
496  // Integral over interpolation aria
497  if(xMin < p[3]) {
498 
499  x1 = p[1];
500  y1 = p[4];
501 
502  G4double dx = (p[2] - p[1]) / 3.0;
503  G4double dx1= std::exp(std::log(p[3]/p[2]) / 16.0);
504 
505  for (size_t i=0; i<19; i++) {
506 
507  if (i < 3) {
508  x2 = x1 + dx;
509  } else if(18 == i) {
510  x2 = p[3];
511  } else {
512  x2 = x1*dx1;
513  }
514 
515  y2 = p[5 + i];
516 
517  if (xMax <= x1) {
518  break;
519  } else if (xMin < x2) {
520 
521  xs1 = x1;
522  xs2 = x2;
523  ys1 = y1;
524  ys2 = y2;
525 
526  if (x2 > x1) {
527  if (xMin > x1) {
528  xs1 = xMin;
529  ys1 += (xs1 - x1)*(y2 - y1)/(x2 - x1);
530  }
531  if (xMax < x2) {
532  xs2 = xMax;
533  ys2 += (xs2 - x2)*(y1 - y2)/(x1 - x2);
534  }
535  if (xs2 > xs1) {
536  sum += std::log(xs2/xs1)*(ys1*xs2 - ys2*xs1)/(xs2 - xs1)
537  + ys2 - ys1;
538  }
539  }
540  }
541  x1 = x2;
542  y1 = y2;
543 
544  }
545  }
546 
547  // Integral over aria with parametrised formula
548 
549  x1 = std::max(xMin, p[3]);
550  if(x1 >= xMax) return sum;
551  x2 = xMax;
552 
553  xs1 = 1./x1;
554  xs2 = 1./x2;
555 
556  sum += std::log(x2/x1)*(1.0 - p[0])
557  + 0.5*(1. - p[iMax])*(x2*x2 - x1*x1)
558  + 1./(1. - x2) - 1./(1. - x1)
559  + (1. + p[iMax])*std::log((1. - x2)/(1. - x1))
560  + 0.5*p[0]*(xs1 - xs2);
561 
562  return sum;
563 }
564 
565 
566 void G4RDeIonisationSpectrum::PrintData() const
567 {
568  theParam->PrintData();
569 }
570 
571 G4double G4RDeIonisationSpectrum::MaxEnergyOfSecondaries(G4double kineticEnergy,
572  G4int, // Z = 0,
573  const G4ParticleDefinition* ) const
574 {
575  return 0.5 * kineticEnergy;
576 }
MeV
static const double MeV
Definition: G4SIunits.hh:193
a1
static const G4double a1
Definition: G4BetheHeitlerModel.cc:72
a4
static const G4double a4
Definition: G4BetheHeitlerModel.cc:73
G4RDeIonisationSpectrum.hh
G4int
int G4int
Definition: G4Types.hh:78
G4RDeIonisationSpectrum::IntSpectrum
G4double IntSpectrum(G4double xMin, G4double xMax, const G4DataVector &p) const
Definition: G4RDeIonisationSpectrum.cc:403
G4RDAtomicTransitionManager::Instance
static G4RDAtomicTransitionManager * Instance()
Definition: G4RDAtomicTransitionManager.cc:163
G4RDeIonisationSpectrum::SampleEnergy
G4double SampleEnergy(G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell, const G4ParticleDefinition *pd=0) const
Definition: G4RDeIonisationSpectrum.cc:245
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:93
G4cout
G4GLOB_DLL std::ostream G4cout
Randomize.hh
a3
static const G4double a3
Definition: G4BetheHeitlerModel.cc:73
G4PhysicalConstants.hh
eV
static const double eV
Definition: G4SIunits.hh:194
factor
static const G4double factor
Definition: G4ContinuumGammaTransition.cc:63
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
g
static const double g
Definition: G4SIunits.hh:162
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:117
G4RDeIonisationSpectrum::MaxEnergyOfSecondaries
G4double MaxEnergyOfSecondaries(G4double kineticEnergy, G4int Z=0, const G4ParticleDefinition *pd=0) const
Definition: G4RDeIonisationSpectrum.cc:571
G4endl
#define G4endl
Definition: G4ios.hh:61
G4RDeIonisationSpectrum::theParam
G4RDeIonisationParameters * theParam
Definition: G4RDeIonisationSpectrum.hh:112
G4RDeIonisationParameters::Parameter
G4double Parameter(G4int Z, G4int shellIndex, G4int parameterIndex, G4double e) const
Definition: G4RDeIonisationParameters.cc:91
G4RDeIonisationSpectrum::AverageValue
G4double AverageValue(G4double xMin, G4double xMax, const G4DataVector &p) const
Definition: G4RDeIonisationSpectrum.cc:487
keV
static const double keV
Definition: G4SIunits.hh:195
G4RDAtomicShell.hh
G4double
double G4double
Definition: G4Types.hh:76
G4RDeIonisationSpectrum::Probability
G4double Probability(G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell, const G4ParticleDefinition *pd=0) const
Definition: G4RDeIonisationSpectrum.cc:77
G4SystemOfUnits.hh
G4DataVector.hh
G4InuclSpecialFunctions::bindingEnergy
G4double bindingEnergy(G4int A, G4int Z)
Definition: bindingEnergy.cc:37
G4RDeIonisationSpectrum::Function
G4double Function(G4double x, const G4DataVector &p) const
Definition: G4RDeIonisationSpectrum.hh:120
a2
static const G4double a2
Definition: G4BetheHeitlerModel.cc:72
G4RDAtomicTransitionManager.hh
G4RDeIonisationSpectrum::AverageEnergy
G4double AverageEnergy(G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell, const G4ParticleDefinition *pd=0) const
Definition: G4RDeIonisationSpectrum.cc:163