G4eIonisationSpectrum Class Reference

#include <G4eIonisationSpectrum.hh>

Inheritance diagram for G4eIonisationSpectrum:

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 ()
	G4VEnergySpectrum (const G4VEnergySpectrum &)=delete
G4VEnergySpectrum &	operator= (const G4VEnergySpectrum &right)=delete

Detailed Description

Definition at line 63 of file G4eIonisationSpectrum.hh.

Constructor & Destructor Documentation

G4eIonisationSpectrum()

G4eIonisationSpectrum::G4eIonisationSpectrum

(

)

Definition at line 61 of file G4eIonisationSpectrum.cc.

                                            :G4VEnergySpectrum(),
  lowestE(0.1*eV),
  factor(1.3),
  iMax(24),
  verbose(0)
{
  theParam = new G4eIonisationParameters();
}

~G4eIonisationSpectrum()

G4eIonisationSpectrum::~G4eIonisationSpectrum

(

)

Definition at line 71 of file G4eIonisationSpectrum.cc.

{
  delete theParam;
}

Member Function Documentation

AverageEnergy()

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

virtual

Implements G4VEnergySpectrum.

Definition at line 172 of file G4eIonisationSpectrum.cc.

{
  // Please comment what AverageEnergy does and what are the three
  // functions mentioned below
  // Describe the algorithms used
 
  G4double eMax = MaxEnergyOfSecondaries(e);
  G4double t0 = std::max(tMin, lowestE);
  G4double tm = std::min(tMax, eMax);
  if(t0 >= tm) return 0.0;
 
  G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
                           Shell(Z, shell)->BindingEnergy();
 
  if(e <= bindingEnergy) return 0.0;
 
  G4double energy = e + bindingEnergy;
 
  G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
  G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
 
  if(verbose > 1) {
    G4cout << "G4eIonisationSpectrum::AverageEnergy: Z= " << Z
           << "; shell= " << shell
           << "; E(keV)= " << e/keV
           << "; bindingE(keV)= " << bindingEnergy/keV
           << "; x1= " << x1
           << "; x2= " << x2
           << G4endl;
  }
 
  G4DataVector p;
 
  // Access parameters
  for (G4int i=0; i<iMax; i++)
  {
    G4double x = theParam->Parameter(Z, shell, i, e);
    if(i<4) x /= energy;
    p.push_back(x);
  }
 
  if(p[3] > 0.5) p[3] = 0.5;
 
  G4double gLocal2 = energy/electron_mass_c2 + 1.;
  p.push_back((2.0*gLocal2 - 1.0)/(gLocal2*gLocal2));
 
 
  //Add protection against division by zero: actually in Function()
  //parameter p[3] appears in the denominator. Bug report 1042
  if (p[3] > 0)
    p[iMax-1] = Function(p[3], p);
  else
    {
      G4cout << "WARNING: G4eIonisationSpectrum::AverageEnergy "
         << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = "
         << Z << ". Please check and/or update it " << G4endl;
    }
 
  G4double val = AverageValue(x1, x2, p);
  G4double x0  = (lowestE + bindingEnergy)/energy;
  G4double nor = IntSpectrum(x0, 0.5, p);
  val *= energy;
 
  if(verbose > 1) {
    G4cout << "tcut(MeV)= " << tMin/MeV
           << "; tMax(MeV)= " << tMax/MeV
           << "; x0= " << x0
           << "; x1= " << x1
           << "; x2= " << x2
           << "; val= " << val
           << "; nor= " << nor
           << "; sum= " << p[0]
           << "; a= " << p[1]
           << "; b= " << p[2]
           << "; c= " << p[3]
           << G4endl;
  }
 
  p.clear();
 
  if(nor > 0.0) val /= nor;
  else          val  = 0.0;
 
  return val;
}

Excitation()

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

inlinevirtual

Implements G4VEnergySpectrum.

Definition at line 128 of file G4eIonisationSpectrum.hh.

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

MaxEnergyOfSecondaries()

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

virtual

Implements G4VEnergySpectrum.

Definition at line 600 of file G4eIonisationSpectrum.cc.

{
  return 0.5 * kineticEnergy;
}

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

PrintData()

void G4eIonisationSpectrum::PrintData ( ) const

virtual

Implements G4VEnergySpectrum.

Definition at line 595 of file G4eIonisationSpectrum.cc.

{
  theParam->PrintData();
}

Probability()

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

virtual

Implements G4VEnergySpectrum.

Definition at line 77 of file G4eIonisationSpectrum.cc.

{
  // Please comment what Probability does and what are the three
  // functions mentioned below
  // Describe the algorithms used
 
  G4double eMax = MaxEnergyOfSecondaries(e);
  G4double t0 = std::max(tMin, lowestE);
  G4double tm = std::min(tMax, eMax);
  if(t0 >= tm) return 0.0;
 
  G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
                           Shell(Z, shell)->BindingEnergy();
 
  if(e <= bindingEnergy) return 0.0;
 
  G4double energy = e + bindingEnergy;
 
  G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
  G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
 
  if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0.0)) {
    G4cout << "G4eIonisationSpectrum::Probability: Z= " << Z
           << "; shell= " << shell
           << "; E(keV)= " << e/keV
           << "; Eb(keV)= " << bindingEnergy/keV
           << "; x1= " << x1
           << "; x2= " << x2
           << G4endl;
 
  }
 
  G4DataVector p;
 
  // Access parameters
  for (G4int i=0; i<iMax; i++)
  {
    G4double x = theParam->Parameter(Z, shell, i, e);
    if(i<4) x /= energy;
    p.push_back(x);
  }
 
  if(p[3] > 0.5) p[3] = 0.5;
 
  G4double gLocal = energy/electron_mass_c2 + 1.;
  p.push_back((2.0*gLocal - 1.0)/(gLocal*gLocal));
 
  //Add protection against division by zero: actually in Function()
  //parameter p[3] appears in the denominator. Bug report 1042
  if (p[3] > 0)
    p[iMax-1] = Function(p[3], p);
  else
    {
      G4cout << "WARNING: G4eIonisationSpectrum::Probability "
         << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = "
         << Z << ". Please check and/or update it " << G4endl;
    }
 
  if(e >= 1. && e <= 0. && Z == 4) p.push_back(0.0);
 
 
  G4double val = IntSpectrum(x1, x2, p);
  G4double x0  = (lowestE + bindingEnergy)/energy;
  G4double nor = IntSpectrum(x0, 0.5, p);
 
  if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0.0)) {
    G4cout << "tcut= " << tMin
           << "; tMax= " << tMax
           << "; x0= " << x0
           << "; x1= " << x1
           << "; x2= " << x2
           << "; val= " << val
           << "; nor= " << nor
           << "; sum= " << p[0]
           << "; a= " << p[1]
           << "; b= " << p[2]
           << "; c= " << p[3]
           << G4endl;
    if(shell == 1) G4cout << "============" << G4endl;
  }
 
  p.clear();
 
  if(nor > 0.0) val /= nor;
  else          val  = 0.0;
 
  return val;
}

SampleEnergy()

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

virtual

Implements G4VEnergySpectrum.

Definition at line 264 of file G4eIonisationSpectrum.cc.

{
  // Please comment what SampleEnergy does
  G4double tDelta = 0.0;
  G4double t0 = std::max(tMin, lowestE);
  G4double tm = std::min(tMax, MaxEnergyOfSecondaries(e));
  if(t0 > tm) return tDelta;
 
  G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
                           Shell(Z, shell)->BindingEnergy();
 
  if(e <= bindingEnergy) return 0.0;
 
  G4double energy = e + bindingEnergy;
 
  G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
  G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
  if(x1 >= x2) return tDelta;
 
  if(verbose > 1) {
    G4cout << "G4eIonisationSpectrum::SampleEnergy: Z= " << Z
           << "; shell= " << shell
           << "; E(keV)= " << e/keV
           << G4endl;
  }
 
  // Access parameters
  G4DataVector p;
 
  // Access parameters
  for (G4int i=0; i<iMax; i++)
  {
    G4double x = theParam->Parameter(Z, shell, i, e);
    if(i<4) x /= energy;
    p.push_back(x);
  }
 
  if(p[3] > 0.5) p[3] = 0.5;
 
  G4double gLocal3 = energy/electron_mass_c2 + 1.;
  p.push_back((2.0*gLocal3 - 1.0)/(gLocal3*gLocal3));
 
 
  //Add protection against division by zero: actually in Function()
  //parameter p[3] appears in the denominator. Bug report 1042
  if (p[3] > 0)
    p[iMax-1] = Function(p[3], p);
  else
    {
      G4cout << "WARNING: G4eIonisationSpectrum::SampleSpectrum "
         << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = "
         << Z << ". Please check and/or update it " << G4endl;
    }
 
  G4double aria1 = 0.0;
  G4double a1 = std::max(x1,p[1]);
  G4double a2 = std::min(x2,p[3]);
  if(a1 < a2) aria1 = IntSpectrum(a1, a2, p);
  G4double aria2 = 0.0;
  G4double a3 = std::max(x1,p[3]);
  G4double a4 = x2;
  if(a3 < a4) aria2 = IntSpectrum(a3, a4, p);
 
  G4double aria = (aria1 + aria2)*G4UniformRand();
  G4double amaj, fun, q, x, z1, z2, dx, dx1;
 
  //======= First aria to sample =====
 
  if(aria <= aria1) {
 
    amaj = p[4];
    for (G4int j=5; j<iMax; j++) {
      if(p[j] > amaj) amaj = p[j];
    }
 
    a1 = 1./a1;
    a2 = 1./a2;
 
    G4int i;
    do {
 
      x = 1./(a2 + G4UniformRand()*(a1 - a2));
      z1 = p[1];
      z2 = p[3];
      dx = (p[2] - p[1]) / 3.0;
      dx1= G4Exp(std::log(p[3]/p[2]) / 16.0);
      for (i=4; i<iMax-1; i++) {
 
        if (i < 7) {
          z2 = z1 + dx;
        } else if(iMax-2 == i) {
          z2 = p[3];
          break;
        } else {
          z2 = z1*dx1;
        }
        if(x >= z1 && x <= z2) break;
        z1 = z2;
      }
      fun = p[i] + (x - z1) * (p[i+1] - p[i])/(z2 - z1);
 
      if(fun > amaj) {
          G4cout << "WARNING in G4eIonisationSpectrum::SampleEnergy:"
                 << " Majoranta " << amaj
                 << " < " << fun
                 << " in the first aria at x= " << x
                 << G4endl;
      }
 
      q = amaj*G4UniformRand();
 
    } while (q >= fun);
 
  //======= Second aria to sample =====
 
  } else {
 
    amaj = std::max(p[iMax-1], Function(0.5, p)) * factor;
    a1 = 1./a3;
    a2 = 1./a4;
 
    do {
 
      x = 1./(a2 + G4UniformRand()*(a1 - a2));
      fun = Function(x, p);
 
      if(fun > amaj) {
          G4cout << "WARNING in G4eIonisationSpectrum::SampleEnergy:"
                 << " Majoranta " << amaj
                 << " < " << fun
                 << " in the second aria at x= " << x
                 << G4endl;
      }
 
      q = amaj*G4UniformRand();
 
    } while (q >= fun);
 
  }
 
  p.clear();
 
  tDelta = x*energy - bindingEnergy;
 
  if(verbose > 1) {
    G4cout << "tcut(MeV)= " << tMin/MeV
           << "; tMax(MeV)= " << tMax/MeV
           << "; x1= " << x1
           << "; x2= " << x2
           << "; a1= " << a1
           << "; a2= " << a2
           << "; x= " << x
           << "; be= " << bindingEnergy
           << "; e= " << e
           << "; tDelta= " << tDelta
           << G4endl;
  }
 
 
  return tDelta;
}

---

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

• G4eIonisationSpectrum.hh
• G4eIonisationSpectrum.cc