G4eIonisationSpectrum Class Reference

#include <G4eIonisationSpectrum.hh>

Inheritance diagram for G4eIonisationSpectrum:

G4VEnergySpectrum

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

Detailed Description

Definition at line 64 of file G4eIonisationSpectrum.hh.

Constructor & Destructor Documentation

G4eIonisationSpectrum::G4eIonisationSpectrum (  ) 

Definition at line 62 of file G4eIonisationSpectrum.cc.

00062                                             :G4VEnergySpectrum(),
00063   lowestE(0.1*eV),
00064   factor(1.3),
00065   iMax(24),
00066   verbose(0)
00067 {
00068   theParam = new G4eIonisationParameters();
00069 }

G4eIonisationSpectrum::~G4eIonisationSpectrum (  ) 

Definition at line 72 of file G4eIonisationSpectrum.cc.

00073 {
00074   delete theParam;
00075 }

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.

References G4InuclSpecialFunctions::bindingEnergy(), G4cout, G4endl, G4AtomicTransitionManager::Instance(), MaxEnergyOfSecondaries(), and G4eIonisationParameters::Parameter().

00179 {
00180   // Please comment what AverageEnergy does and what are the three 
00181   // functions mentioned below
00182   // Describe the algorithms used
00183 
00184   G4double eMax = MaxEnergyOfSecondaries(e);
00185   G4double t0 = std::max(tMin, lowestE);
00186   G4double tm = std::min(tMax, eMax);
00187   if(t0 >= tm) return 0.0;
00188 
00189   G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
00190                            Shell(Z, shell)->BindingEnergy();
00191 
00192   if(e <= bindingEnergy) return 0.0;
00193 
00194   G4double energy = e + bindingEnergy;
00195 
00196   G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
00197   G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
00198 
00199   if(verbose > 1) {
00200     G4cout << "G4eIonisationSpectrum::AverageEnergy: Z= " << Z
00201            << "; shell= " << shell
00202            << "; E(keV)= " << e/keV
00203            << "; bindingE(keV)= " << bindingEnergy/keV
00204            << "; x1= " << x1 
00205            << "; x2= " << x2 
00206            << G4endl;
00207   }
00208 
00209   G4DataVector p;
00210 
00211   // Access parameters
00212   for (G4int i=0; i<iMax; i++) 
00213   {
00214     G4double x = theParam->Parameter(Z, shell, i, e);
00215     if(i<4) x /= energy; 
00216     p.push_back(x);
00217   }
00218 
00219   if(p[3] > 0.5) p[3] = 0.5;
00220 
00221   G4double gLocal2 = energy/electron_mass_c2 + 1.;
00222   p.push_back((2.0*gLocal2 - 1.0)/(gLocal2*gLocal2));
00223 
00224   
00225   //Add protection against division by zero: actually in Function() 
00226   //parameter p[3] appears in the denominator. Bug report 1042
00227   if (p[3] > 0)
00228     p[iMax-1] = Function(p[3], p);
00229   else
00230     {
00231       G4cout << "WARNING: G4eIonisationSpectrum::AverageEnergy "
00232              << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = " 
00233              << Z << ". Please check and/or update it " << G4endl;      
00234     }
00235     
00236   G4double val = AverageValue(x1, x2, p);
00237   G4double x0  = (lowestE + bindingEnergy)/energy;
00238   G4double nor = IntSpectrum(x0, 0.5, p);
00239   val *= energy;
00240 
00241   if(verbose > 1) {
00242     G4cout << "tcut(MeV)= " << tMin/MeV 
00243            << "; tMax(MeV)= " << tMax/MeV 
00244            << "; x0= " << x0 
00245            << "; x1= " << x1 
00246            << "; x2= " << x2 
00247            << "; val= " << val 
00248            << "; nor= " << nor 
00249            << "; sum= " << p[0] 
00250            << "; a= " << p[1] 
00251            << "; b= " << p[2] 
00252            << "; c= " << p[3] 
00253            << G4endl;
00254   }
00255 
00256   p.clear();
00257 
00258   if(nor > 0.0) val /= nor;
00259   else          val  = 0.0;
00260 
00261   return val; 
00262 }

G4double G4eIonisationSpectrum::Excitation ( G4int  Z,
G4double  e 
) const [inline, virtual]

Implements G4VEnergySpectrum.

Definition at line 129 of file G4eIonisationSpectrum.hh.

References G4eIonisationParameters::Excitation().

00130 {
00131   return theParam->Excitation(Z, e);
00132 }

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

Implements G4VEnergySpectrum.

Definition at line 601 of file G4eIonisationSpectrum.cc.

Referenced by AverageEnergy(), Probability(), and SampleEnergy().

00604 {
00605   return 0.5 * kineticEnergy;
00606 }

void G4eIonisationSpectrum::PrintData (  )  const [virtual]

Implements G4VEnergySpectrum.

Definition at line 596 of file G4eIonisationSpectrum.cc.

References G4eIonisationParameters::PrintData().

00597 {
00598   theParam->PrintData();
00599 }

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.

References G4InuclSpecialFunctions::bindingEnergy(), G4cout, G4endl, G4AtomicTransitionManager::Instance(), MaxEnergyOfSecondaries(), and G4eIonisationParameters::Parameter().

00084 {
00085   // Please comment what Probability does and what are the three 
00086   // functions mentioned below
00087   // Describe the algorithms used
00088 
00089   G4double eMax = MaxEnergyOfSecondaries(e);
00090   G4double t0 = std::max(tMin, lowestE);
00091   G4double tm = std::min(tMax, eMax);
00092   if(t0 >= tm) return 0.0;
00093 
00094   G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
00095                            Shell(Z, shell)->BindingEnergy();
00096 
00097   if(e <= bindingEnergy) return 0.0;
00098 
00099   G4double energy = e + bindingEnergy;
00100 
00101   G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
00102   G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
00103 
00104   if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0.0)) {
00105     G4cout << "G4eIonisationSpectrum::Probability: Z= " << Z
00106            << "; shell= " << shell
00107            << "; E(keV)= " << e/keV
00108            << "; Eb(keV)= " << bindingEnergy/keV
00109            << "; x1= " << x1 
00110            << "; x2= " << x2 
00111            << G4endl;
00112      
00113   }
00114 
00115   G4DataVector p;
00116 
00117   // Access parameters
00118   for (G4int i=0; i<iMax; i++) 
00119   {
00120     G4double x = theParam->Parameter(Z, shell, i, e);
00121     if(i<4) x /= energy; 
00122     p.push_back(x); 
00123   }
00124 
00125   if(p[3] > 0.5) p[3] = 0.5;
00126   
00127   G4double gLocal = energy/electron_mass_c2 + 1.;
00128   p.push_back((2.0*gLocal - 1.0)/(gLocal*gLocal));
00129   
00130   //Add protection against division by zero: actually in Function() 
00131   //parameter p[3] appears in the denominator. Bug report 1042
00132   if (p[3] > 0)
00133     p[iMax-1] = Function(p[3], p);
00134   else
00135     {
00136       G4cout << "WARNING: G4eIonisationSpectrum::Probability "
00137              << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = " 
00138              << Z << ". Please check and/or update it " << G4endl;      
00139     }
00140 
00141   if(e >= 1. && e <= 0. && Z == 4) p.push_back(0.0);
00142 
00143   
00144   G4double val = IntSpectrum(x1, x2, p);
00145   G4double x0  = (lowestE + bindingEnergy)/energy;
00146   G4double nor = IntSpectrum(x0, 0.5, p);
00147   
00148   if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0.0)) {
00149     G4cout << "tcut= " << tMin 
00150            << "; tMax= " << tMax 
00151            << "; x0= " << x0 
00152            << "; x1= " << x1 
00153            << "; x2= " << x2 
00154            << "; val= " << val 
00155            << "; nor= " << nor 
00156            << "; sum= " << p[0] 
00157            << "; a= " << p[1] 
00158            << "; b= " << p[2] 
00159            << "; c= " << p[3] 
00160            << G4endl;
00161     if(shell == 1) G4cout << "============" << G4endl; 
00162   }
00163 
00164   p.clear();
00165 
00166   if(nor > 0.0) val /= nor;
00167   else          val  = 0.0;
00168 
00169   return val; 
00170 }

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.

References G4InuclSpecialFunctions::bindingEnergy(), G4cout, G4endl, G4UniformRand, G4AtomicTransitionManager::Instance(), MaxEnergyOfSecondaries(), and G4eIonisationParameters::Parameter().

00271 {
00272   // Please comment what SampleEnergy does
00273   G4double tDelta = 0.0;
00274   G4double t0 = std::max(tMin, lowestE);
00275   G4double tm = std::min(tMax, MaxEnergyOfSecondaries(e));
00276   if(t0 > tm) return tDelta;
00277 
00278   G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
00279                            Shell(Z, shell)->BindingEnergy();
00280 
00281   if(e <= bindingEnergy) return 0.0;
00282 
00283   G4double energy = e + bindingEnergy;
00284 
00285   G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
00286   G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
00287   if(x1 >= x2) return tDelta;
00288 
00289   if(verbose > 1) {
00290     G4cout << "G4eIonisationSpectrum::SampleEnergy: Z= " << Z
00291            << "; shell= " << shell
00292            << "; E(keV)= " << e/keV
00293            << G4endl;
00294   }
00295 
00296   // Access parameters
00297   G4DataVector p;
00298 
00299   // Access parameters
00300   for (G4int i=0; i<iMax; i++) 
00301   {
00302     G4double x = theParam->Parameter(Z, shell, i, e);
00303     if(i<4) x /= energy; 
00304     p.push_back(x);
00305   }
00306 
00307   if(p[3] > 0.5) p[3] = 0.5;
00308 
00309   G4double gLocal3 = energy/electron_mass_c2 + 1.;
00310   p.push_back((2.0*gLocal3 - 1.0)/(gLocal3*gLocal3));
00311 
00312   
00313   //Add protection against division by zero: actually in Function() 
00314   //parameter p[3] appears in the denominator. Bug report 1042
00315   if (p[3] > 0)
00316     p[iMax-1] = Function(p[3], p);
00317   else
00318     {
00319       G4cout << "WARNING: G4eIonisationSpectrum::SampleSpectrum "
00320              << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = " 
00321              << Z << ". Please check and/or update it " << G4endl;      
00322     }
00323 
00324   G4double aria1 = 0.0;
00325   G4double a1 = std::max(x1,p[1]);
00326   G4double a2 = std::min(x2,p[3]);
00327   if(a1 < a2) aria1 = IntSpectrum(a1, a2, p);
00328   G4double aria2 = 0.0;
00329   G4double a3 = std::max(x1,p[3]);
00330   G4double a4 = x2;
00331   if(a3 < a4) aria2 = IntSpectrum(a3, a4, p);
00332 
00333   G4double aria = (aria1 + aria2)*G4UniformRand();
00334   G4double amaj, fun, q, x, z1, z2, dx, dx1;
00335 
00336   //======= First aria to sample =====
00337 
00338   if(aria <= aria1) { 
00339 
00340     amaj = p[4];
00341     for (G4int j=5; j<iMax; j++) {
00342       if(p[j] > amaj) amaj = p[j];
00343     }
00344 
00345     a1 = 1./a1;
00346     a2 = 1./a2;
00347 
00348     G4int i;
00349     do {
00350 
00351       x = 1./(a2 + G4UniformRand()*(a1 - a2));
00352       z1 = p[1];
00353       z2 = p[3];
00354       dx = (p[2] - p[1]) / 3.0;
00355       dx1= std::exp(std::log(p[3]/p[2]) / 16.0);
00356       for (i=4; i<iMax-1; i++) {
00357 
00358         if (i < 7) {
00359           z2 = z1 + dx;
00360         } else if(iMax-2 == i) {
00361           z2 = p[3];
00362           break;
00363         } else {
00364           z2 = z1*dx1;
00365         }
00366         if(x >= z1 && x <= z2) break;
00367         z1 = z2;
00368       }
00369       fun = p[i] + (x - z1) * (p[i+1] - p[i])/(z2 - z1);
00370 
00371       if(fun > amaj) {
00372           G4cout << "WARNING in G4eIonisationSpectrum::SampleEnergy:" 
00373                  << " Majoranta " << amaj 
00374                  << " < " << fun
00375                  << " in the first aria at x= " << x
00376                  << G4endl;
00377       }
00378 
00379       q = amaj*G4UniformRand();
00380 
00381     } while (q >= fun);
00382 
00383   //======= Second aria to sample =====
00384 
00385   } else {
00386 
00387     amaj = std::max(p[iMax-1], Function(0.5, p)) * factor;
00388     a1 = 1./a3;
00389     a2 = 1./a4;
00390 
00391     do {
00392 
00393       x = 1./(a2 + G4UniformRand()*(a1 - a2));
00394       fun = Function(x, p);
00395 
00396       if(fun > amaj) {
00397           G4cout << "WARNING in G4eIonisationSpectrum::SampleEnergy:" 
00398                  << " Majoranta " << amaj 
00399                  << " < " << fun
00400                  << " in the second aria at x= " << x
00401                  << G4endl;
00402       }
00403 
00404       q = amaj*G4UniformRand();
00405 
00406     } while (q >= fun);
00407 
00408   }
00409 
00410   p.clear();
00411 
00412   tDelta = x*energy - bindingEnergy;
00413 
00414   if(verbose > 1) {
00415     G4cout << "tcut(MeV)= " << tMin/MeV 
00416            << "; tMax(MeV)= " << tMax/MeV 
00417            << "; x1= " << x1 
00418            << "; x2= " << x2 
00419            << "; a1= " << a1 
00420            << "; a2= " << a2 
00421            << "; x= " << x 
00422            << "; be= " << bindingEnergy 
00423            << "; e= " << e 
00424            << "; tDelta= " << tDelta 
00425            << G4endl;
00426   }
00427 
00428 
00429   return tDelta; 
00430 }

The documentation for this class was generated from the following files:
Generated on Mon May 27 17:51:50 2013 for Geant4 by  doxygen 1.4.7