Geant4  10.02.p03
Public Member Functions | Private Member Functions | Private Attributes | List of all members
G4RDeIonisationSpectrum Class Reference

#include <G4RDeIonisationSpectrum.hh>

Inheritance diagram for G4RDeIonisationSpectrum:
Collaboration diagram for G4RDeIonisationSpectrum:

Public Member Functions

 G4RDeIonisationSpectrum ()
 
 ~G4RDeIonisationSpectrum ()
 
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 G4RDVEnergySpectrum
 G4RDVEnergySpectrum ()
 
virtual ~G4RDVEnergySpectrum ()
 

Private Member Functions

G4double IntSpectrum (G4double xMin, G4double xMax, const G4DataVector &p) const
 
G4double AverageValue (G4double xMin, G4double xMax, const G4DataVector &p) const
 
G4double Function (G4double x, const G4DataVector &p) const
 
 G4RDeIonisationSpectrum (const G4RDeIonisationSpectrum &)
 
G4RDeIonisationSpectrum & operator= (const G4RDeIonisationSpectrum &right)
 

Private Attributes

G4RDeIonisationParameters * theParam
 
G4double lowestE
 
G4double factor
 
G4int iMax
 
G4int verbose
 

Detailed Description

Definition at line 65 of file G4RDeIonisationSpectrum.hh.

Constructor & Destructor Documentation

◆ G4RDeIonisationSpectrum() [1/2]

G4RDeIonisationSpectrum::G4RDeIonisationSpectrum ( )

Definition at line 61 of file G4RDeIonisationSpectrum.cc.

61  :G4RDVEnergySpectrum(),
62  lowestE(0.1*eV),
63  factor(1.3),
64  iMax(24),
65  verbose(0)
66 {
67  theParam = new G4RDeIonisationParameters();
68 }
eV
static const double eV
Definition: G4SIunits.hh:212
G4RDeIonisationSpectrum::theParam
G4RDeIonisationParameters * theParam
Definition: G4RDeIonisationSpectrum.hh:112

◆ ~G4RDeIonisationSpectrum()

G4RDeIonisationSpectrum::~G4RDeIonisationSpectrum ( )

Definition at line 71 of file G4RDeIonisationSpectrum.cc.

72 {
73  delete theParam;
74 }
G4RDeIonisationSpectrum::theParam
G4RDeIonisationParameters * theParam
Definition: G4RDeIonisationSpectrum.hh:112

◆ G4RDeIonisationSpectrum() [2/2]

G4RDeIonisationSpectrum::G4RDeIonisationSpectrum ( const G4RDeIonisationSpectrum &  )
private

Member Function Documentation

◆ AverageEnergy()

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

Implements G4RDVEnergySpectrum.

Definition at line 163 of file G4RDeIonisationSpectrum.cc.

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 }
G4RDeIonisationSpectrum::MaxEnergyOfSecondaries
G4double MaxEnergyOfSecondaries(G4double kineticEnergy, G4int Z=0, const G4ParticleDefinition *pd=0) const
Definition: G4RDeIonisationSpectrum.cc:571
MeV
static const double MeV
Definition: G4SIunits.hh:211
x2
Double_t x2[nxs]
Definition: Hadr00/scripts/Style.C:19
G4RDeIonisationSpectrum::IntSpectrum
G4double IntSpectrum(G4double xMin, G4double xMax, const G4DataVector &p) const
Definition: G4RDeIonisationSpectrum.cc:403
G4int
int G4int
Definition: G4Types.hh:78
G4RDAtomicTransitionManager::Instance
static G4RDAtomicTransitionManager * Instance()
Definition: G4RDAtomicTransitionManager.cc:163
G4RDeIonisationParameters::Parameter
G4double Parameter(G4int Z, G4int shellIndex, G4int parameterIndex, G4double e) const
Definition: G4RDeIonisationParameters.cc:91
g
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5206
G4cout
G4GLOB_DLL std::ostream G4cout
energy
double energy
Definition: plottest35.C:25
Z
Float_t Z
Definition: advanced/microbeam/plot.C:39
x1
Double_t x1[nxs]
Definition: Hadr00/scripts/Style.C:18
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4RDeIonisationSpectrum::Function
G4double Function(G4double x, const G4DataVector &p) const
Definition: G4RDeIonisationSpectrum.hh:120
G4RDeIonisationSpectrum::AverageValue
G4double AverageValue(G4double xMin, G4double xMax, const G4DataVector &p) const
Definition: G4RDeIonisationSpectrum.cc:487
G4endl
#define G4endl
Definition: G4ios.hh:61
G4RDeIonisationSpectrum::theParam
G4RDeIonisationParameters * theParam
Definition: G4RDeIonisationSpectrum.hh:112
keV
static const double keV
Definition: G4SIunits.hh:213
G4double
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:

◆ AverageValue()

G4double G4RDeIonisationSpectrum::AverageValue ( G4double  xMin,
G4double  xMax,
const G4DataVector &  p 
) const
private

Definition at line 487 of file G4RDeIonisationSpectrum.cc.

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 }
y2
Double_t y2[nxs]
Definition: Hadr00/scripts/Style.C:21
y1
Double_t y1[nxs]
Definition: Hadr00/scripts/Style.C:20
x2
Double_t x2[nxs]
Definition: Hadr00/scripts/Style.C:19
x1
Double_t x1[nxs]
Definition: Hadr00/scripts/Style.C:18
G4double
double G4double
Definition: G4Types.hh:76
Here is the caller graph for this function:

◆ Excitation()

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

Implements G4RDVEnergySpectrum.

Definition at line 130 of file G4RDeIonisationSpectrum.hh.

131 {
132  return theParam->Excitation(Z, e);
133 }
G4RDeIonisationParameters::Excitation
G4double Excitation(G4int Z, G4double e) const
Definition: G4RDeIonisationParameters.cc:124
Z
Float_t Z
Definition: advanced/microbeam/plot.C:39
G4RDeIonisationSpectrum::theParam
G4RDeIonisationParameters * theParam
Definition: G4RDeIonisationSpectrum.hh:112
Here is the call graph for this function:

◆ Function()

G4double G4RDeIonisationSpectrum::Function ( G4double  x,
const G4DataVector &  p 
) const
inlineprivate

Definition at line 120 of file G4RDeIonisationSpectrum.hh.

122 {
123  G4double f = 1.0 - p[0] - p[iMax]*x
124  + x*x*(1.0 - p[iMax] + (1.0/(1.0 - x) - p[iMax])/(1.0 - x))
125  + 0.5*p[0]/x;
126 
127  return f;
128 }
G4double
double G4double
Definition: G4Types.hh:76
Here is the caller graph for this function:

◆ IntSpectrum()

G4double G4RDeIonisationSpectrum::IntSpectrum ( G4double  xMin,
G4double  xMax,
const G4DataVector &  p 
) const
private

Definition at line 403 of file G4RDeIonisationSpectrum.cc.

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 }
y2
Double_t y2[nxs]
Definition: Hadr00/scripts/Style.C:21
y1
Double_t y1[nxs]
Definition: Hadr00/scripts/Style.C:20
x2
Double_t x2[nxs]
Definition: Hadr00/scripts/Style.C:19
G4cout
G4GLOB_DLL std::ostream G4cout
x1
Double_t x1[nxs]
Definition: Hadr00/scripts/Style.C:18
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
Here is the caller graph for this function:

◆ MaxEnergyOfSecondaries()

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

Implements G4RDVEnergySpectrum.

Definition at line 571 of file G4RDeIonisationSpectrum.cc.

574 {
575  return 0.5 * kineticEnergy;
576 }
Here is the caller graph for this function:

◆ operator=()

G4RDeIonisationSpectrum& G4RDeIonisationSpectrum::operator= ( const G4RDeIonisationSpectrum &  right)
private

◆ PrintData()

void G4RDeIonisationSpectrum::PrintData ( void  ) const
virtual

Implements G4RDVEnergySpectrum.

Definition at line 566 of file G4RDeIonisationSpectrum.cc.

567 {
568  theParam->PrintData();
569 }
G4RDeIonisationSpectrum::theParam
G4RDeIonisationParameters * theParam
Definition: G4RDeIonisationSpectrum.hh:112
Here is the call graph for this function:

◆ Probability()

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

Implements G4RDVEnergySpectrum.

Definition at line 77 of file G4RDeIonisationSpectrum.cc.

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 }
G4RDeIonisationSpectrum::MaxEnergyOfSecondaries
G4double MaxEnergyOfSecondaries(G4double kineticEnergy, G4int Z=0, const G4ParticleDefinition *pd=0) const
Definition: G4RDeIonisationSpectrum.cc:571
x2
Double_t x2[nxs]
Definition: Hadr00/scripts/Style.C:19
G4RDeIonisationSpectrum::IntSpectrum
G4double IntSpectrum(G4double xMin, G4double xMax, const G4DataVector &p) const
Definition: G4RDeIonisationSpectrum.cc:403
G4int
int G4int
Definition: G4Types.hh:78
G4RDAtomicTransitionManager::Instance
static G4RDAtomicTransitionManager * Instance()
Definition: G4RDAtomicTransitionManager.cc:163
G4RDeIonisationParameters::Parameter
G4double Parameter(G4int Z, G4int shellIndex, G4int parameterIndex, G4double e) const
Definition: G4RDeIonisationParameters.cc:91
g
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5206
G4cout
G4GLOB_DLL std::ostream G4cout
energy
double energy
Definition: plottest35.C:25
Z
Float_t Z
Definition: advanced/microbeam/plot.C:39
x1
Double_t x1[nxs]
Definition: Hadr00/scripts/Style.C:18
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4RDeIonisationSpectrum::Function
G4double Function(G4double x, const G4DataVector &p) const
Definition: G4RDeIonisationSpectrum.hh:120
G4endl
#define G4endl
Definition: G4ios.hh:61
G4RDeIonisationSpectrum::theParam
G4RDeIonisationParameters * theParam
Definition: G4RDeIonisationSpectrum.hh:112
keV
static const double keV
Definition: G4SIunits.hh:213
G4double
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:

◆ SampleEnergy()

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

Implements G4RDVEnergySpectrum.

Definition at line 245 of file G4RDeIonisationSpectrum.cc.

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 }
G4RDeIonisationSpectrum::MaxEnergyOfSecondaries
G4double MaxEnergyOfSecondaries(G4double kineticEnergy, G4int Z=0, const G4ParticleDefinition *pd=0) const
Definition: G4RDeIonisationSpectrum.cc:571
MeV
static const double MeV
Definition: G4SIunits.hh:211
x2
Double_t x2[nxs]
Definition: Hadr00/scripts/Style.C:19
a1
static const G4double a1
Definition: G4BetheHeitlerModel.cc:72
G4RDeIonisationSpectrum::IntSpectrum
G4double IntSpectrum(G4double xMin, G4double xMax, const G4DataVector &p) const
Definition: G4RDeIonisationSpectrum.cc:403
a4
static const G4double a4
Definition: G4BetheHeitlerModel.cc:73
G4int
int G4int
Definition: G4Types.hh:78
G4RDAtomicTransitionManager::Instance
static G4RDAtomicTransitionManager * Instance()
Definition: G4RDAtomicTransitionManager.cc:163
G4RDeIonisationParameters::Parameter
G4double Parameter(G4int Z, G4int shellIndex, G4int parameterIndex, G4double e) const
Definition: G4RDeIonisationParameters.cc:91
g
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5206
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
energy
double energy
Definition: plottest35.C:25
Z
Float_t Z
Definition: advanced/microbeam/plot.C:39
x1
Double_t x1[nxs]
Definition: Hadr00/scripts/Style.C:18
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
a3
static const G4double a3
Definition: G4BetheHeitlerModel.cc:73
G4RDeIonisationSpectrum::Function
G4double Function(G4double x, const G4DataVector &p) const
Definition: G4RDeIonisationSpectrum.hh:120
G4endl
#define G4endl
Definition: G4ios.hh:61
G4RDeIonisationSpectrum::theParam
G4RDeIonisationParameters * theParam
Definition: G4RDeIonisationSpectrum.hh:112
keV
static const double keV
Definition: G4SIunits.hh:213
G4double
double G4double
Definition: G4Types.hh:76
a2
static const G4double a2
Definition: G4BetheHeitlerModel.cc:72
Here is the call graph for this function:

Member Data Documentation

◆ factor

G4double G4RDeIonisationSpectrum::factor
private

Definition at line 114 of file G4RDeIonisationSpectrum.hh.

◆ iMax

G4int G4RDeIonisationSpectrum::iMax
private

Definition at line 115 of file G4RDeIonisationSpectrum.hh.

◆ lowestE

G4double G4RDeIonisationSpectrum::lowestE
private

Definition at line 113 of file G4RDeIonisationSpectrum.hh.

◆ theParam

G4RDeIonisationParameters* G4RDeIonisationSpectrum::theParam
private

Definition at line 112 of file G4RDeIonisationSpectrum.hh.

◆ verbose

G4int G4RDeIonisationSpectrum::verbose
private

Definition at line 116 of file G4RDeIonisationSpectrum.hh.

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