Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | List of all members
G4eIonisationSpectrum Class Reference

#include <G4eIonisationSpectrum.hh>

Inheritance diagram for G4eIonisationSpectrum:
Inheritance graph
[legend]
Collaboration diagram for G4eIonisationSpectrum:
Collaboration graph
[legend]

Public Member Functions

 G4eIonisationSpectrum ()
 
 ~G4eIonisationSpectrum ()
 
G4double Probability (G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell, const G4ParticleDefinition *pd=0) const
 
G4double AverageEnergy (G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell, const G4ParticleDefinition *pd=0) const
 
G4double SampleEnergy (G4int Z, G4double tMin, G4double tMax, G4double kineticEnergy, G4int shell, const G4ParticleDefinition *pd=0) const
 
G4double MaxEnergyOfSecondaries (G4double kineticEnergy, G4int Z=0, const G4ParticleDefinition *pd=0) const
 
G4double Excitation (G4int Z, G4double e) const
 
void PrintData () const
 
- Public Member Functions inherited from G4VEnergySpectrum
 G4VEnergySpectrum ()
 
virtual ~G4VEnergySpectrum ()
 

Detailed Description

Definition at line 64 of file G4eIonisationSpectrum.hh.

Constructor & Destructor Documentation

G4eIonisationSpectrum::G4eIonisationSpectrum ( )

Definition at line 62 of file G4eIonisationSpectrum.cc.

62  :G4VEnergySpectrum(),
63  lowestE(0.1*eV),
64  factor(1.3),
65  iMax(24),
66  verbose(0)
67 {
68  theParam = new G4eIonisationParameters();
69 }
eV
static constexpr double eV
Definition: G4SIunits.hh:215
G4eIonisationSpectrum::~G4eIonisationSpectrum ( )

Definition at line 72 of file G4eIonisationSpectrum.cc.

73 {
74  delete theParam;
75 }

Member Function Documentation

G4double G4eIonisationSpectrum::AverageEnergy ( G4int  Z,
G4double  tMin,
G4double  tMax,
G4double  kineticEnergy,
G4int  shell,
const G4ParticleDefinition *  pd = 0 
) const
virtual

Implements G4VEnergySpectrum.

Definition at line 173 of file G4eIonisationSpectrum.cc.

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 }
p
const char * p
Definition: xmltok.h:285
G4eIonisationParameters::Parameter
G4double Parameter(G4int Z, G4int shellIndex, G4int parameterIndex, G4double e) const
Definition: G4eIonisationParameters.cc:92
test.x
tuple x
Definition: test.py:50
G4int
int G4int
Definition: G4Types.hh:78
G4cout
G4GLOB_DLL std::ostream G4cout
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
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
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
G4double
double G4double
Definition: G4Types.hh:76
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
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262

Here is the call graph for this function:

G4double G4eIonisationSpectrum::Excitation ( G4int  Z,
G4double  e 
) const
inlinevirtual

Implements G4VEnergySpectrum.

Definition at line 129 of file G4eIonisationSpectrum.hh.

130 {
131  return theParam->Excitation(Z, e);
132 }
G4eIonisationParameters::Excitation
G4double Excitation(G4int Z, G4double e) const
Definition: G4eIonisationParameters.cc:125
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262

Here is the call graph for this function:

G4double G4eIonisationSpectrum::MaxEnergyOfSecondaries ( G4double  kineticEnergy,
G4int  Z = 0,
const G4ParticleDefinition *  pd = 0 
) const
virtual

Implements G4VEnergySpectrum.

Definition at line 601 of file G4eIonisationSpectrum.cc.

604 {
605  return 0.5 * kineticEnergy;
606 }

Here is the caller graph for this function:

void G4eIonisationSpectrum::PrintData ( void  ) const
virtual

Implements G4VEnergySpectrum.

Definition at line 596 of file G4eIonisationSpectrum.cc.

597 {
598  theParam->PrintData();
599 }

Here is the call graph for this function:

G4double G4eIonisationSpectrum::Probability ( G4int  Z,
G4double  tMin,
G4double  tMax,
G4double  kineticEnergy,
G4int  shell,
const G4ParticleDefinition *  pd = 0 
) const
virtual

Implements G4VEnergySpectrum.

Definition at line 78 of file G4eIonisationSpectrum.cc.

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 }
p
const char * p
Definition: xmltok.h:285
G4eIonisationParameters::Parameter
G4double Parameter(G4int Z, G4int shellIndex, G4int parameterIndex, G4double e) const
Definition: G4eIonisationParameters.cc:92
test.x
tuple x
Definition: test.py:50
G4int
int G4int
Definition: G4Types.hh:78
G4cout
G4GLOB_DLL std::ostream G4cout
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
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
G4double
double G4double
Definition: G4Types.hh:76
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
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262

Here is the call graph for this function:

G4double G4eIonisationSpectrum::SampleEnergy ( G4int  Z,
G4double  tMin,
G4double  tMax,
G4double  kineticEnergy,
G4int  shell,
const G4ParticleDefinition *  pd = 0 
) const
virtual

Implements G4VEnergySpectrum.

Definition at line 265 of file G4eIonisationSpectrum.cc.

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 }
p
const char * p
Definition: xmltok.h:285
G4eIonisationParameters::Parameter
G4double Parameter(G4int Z, G4int shellIndex, G4int parameterIndex, G4double e) const
Definition: G4eIonisationParameters.cc:92
test.x
tuple x
Definition: test.py:50
G4int
int G4int
Definition: G4Types.hh:78
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
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
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
G4double
double G4double
Definition: G4Types.hh:76
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
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262

Here is the call graph for this function:

The documentation for this class was generated from the following files: