G4eIonisationSpectrum Class Reference

#include <G4eIonisationSpectrum.hh>

Inheritance diagram for G4eIonisationSpectrum:

Collaboration 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 ()

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

Private Attributes
G4eIonisationParameters *	theParam
G4double	lowestE
G4double	factor
G4int	iMax
G4int	verbose

Detailed Description

Definition at line 64 of file G4eIonisationSpectrum.hh.

Constructor & Destructor Documentation

G4eIonisationSpectrum() [1/2]

G4eIonisationSpectrum::G4eIonisationSpectrum

(

)

Definition at line 62 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 72 of file G4eIonisationSpectrum.cc.

{
  delete theParam;
}

G4eIonisationSpectrum() [2/2]

G4eIonisationSpectrum::G4eIonisationSpectrum ( const G4eIonisationSpectrum &  )

private

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 173 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;
}

Here is the call graph for this function:

AverageValue()

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

private

Definition at line 517 of file G4eIonisationSpectrum.cc.

{
  G4double sum = 0.0;
  if(xMin >= xMax) return sum;

  G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2;

  // Integral over interpolation aria
  if(xMin < p[3]) {

    x1 = p[1];
    y1 = p[4];

    G4double dx = (p[2] - p[1]) / 3.0;
    G4double dx1= G4Exp(std::log(p[3]/p[2]) / 16.0);

    for (size_t i=0; i<19; i++) {

      if (i < 3) {
        x2 = x1 + dx;
      } else if(18 == i) {
        x2 = p[3];
      } else {
        x2 = x1*dx1;
      }

      y2 = p[5 + i];

      if (xMax <= x1) {
        break;
      } else if (xMin < x2) {

        xs1 = x1;
        xs2 = x2;
        ys1 = y1;
        ys2 = y2;

        if (x2 > x1) {
          if (xMin > x1) {
            xs1 = xMin;
            ys1 += (xs1 - x1)*(y2 - y1)/(x2 - x1);
      }
          if (xMax < x2) {
            xs2 = xMax;
            ys2 += (xs2 - x2)*(y1 - y2)/(x1 - x2);
      }
          if (xs2 > xs1) {
            sum += std::log(xs2/xs1)*(ys1*xs2 - ys2*xs1)/(xs2 - xs1)
                +  ys2 - ys1;
      }
    }
      }
      x1 = x2;
      y1 = y2;

    }
  }

  // Integral over aria with parametrised formula

  x1 = std::max(xMin, p[3]);
  if(x1 >= xMax) return sum;
  x2 = xMax;

  xs1 = 1./x1;
  xs2 = 1./x2;

  sum  += std::log(x2/x1)*(1.0 - p[0])
        + 0.5*(1. - p[iMax])*(x2*x2 - x1*x1)
        + 1./(1. - x2) - 1./(1. - x1)
        + (1. + p[iMax])*std::log((1. - x2)/(1. - x1))
        + 0.5*p[0]*(xs1 - xs2);

  return sum;
}

Here is the call graph for this function:

Here is the caller graph for this function:

Excitation()

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

inlinevirtual

Implements G4VEnergySpectrum.

Definition at line 129 of file G4eIonisationSpectrum.hh.

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

Here is the call graph for this function:

Function()

G4double G4eIonisationSpectrum::Function ( G4double  x, const G4DataVector &  p  ) const

inlineprivate

Definition at line 119 of file G4eIonisationSpectrum.hh.

{
  G4double f = 1.0 - p[0] - p[iMax]*x 
             + x*x*(1.0 - p[iMax] + (1.0/(1.0 - x) - p[iMax])/(1.0 - x))
             + 0.5*p[0]/x;

  return f;
} 

Here is the caller graph for this function:

IntSpectrum()

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

private

Definition at line 433 of file G4eIonisationSpectrum.cc.

{
  // Please comment what IntSpectrum does
  G4double sum = 0.0;
  if(xMin >= xMax) return sum;

  G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2, q;

  // Integral over interpolation aria
  if(xMin < p[3]) {

    x1 = p[1];
    y1 = p[4];

    G4double dx = (p[2] - p[1]) / 3.0;
    G4double dx1= G4Exp(std::log(p[3]/p[2]) / 16.0);

    for (size_t i=0; i<19; i++) {

      q = 0.0;
      if (i < 3) {
        x2 = x1 + dx;
      } else if(18 == i) {
        x2 = p[3];
      } else {
        x2 = x1*dx1;
      }

      y2 = p[5 + i];

      if (xMax <= x1) {
        break;
      } else if (xMin < x2) {

        xs1 = x1;
        xs2 = x2;
        ys1 = y1;
        ys2 = y2;

        if (x2 > x1) {
          if (xMin > x1) {
            xs1 = xMin;
            ys1 += (xs1 - x1)*(y2 - y1)/(x2 - x1);
      }
          if (xMax < x2) {
            xs2 = xMax;
            ys2 += (xs2 - x2)*(y1 - y2)/(x1 - x2);
      }
          if (xs2 > xs1) {
            q = (ys1*xs2 - ys2*xs1)/(xs1*xs2)
              +  std::log(xs2/xs1)*(ys2 - ys1)/(xs2 - xs1);
            sum += q;
            if(p.size() == 26) G4cout << "i= " << i << "  q= " << q << " sum= " << sum << G4endl;
      }
    }
      }
      x1 = x2;
      y1 = y2;
    }
  }

  // Integral over aria with parametrised formula

  x1 = std::max(xMin, p[3]);
  if(x1 >= xMax) return sum;
  x2 = xMax;

  xs1 = 1./x1;
  xs2 = 1./x2;
  q = (xs1 - xs2)*(1.0 - p[0])
       - p[iMax]*std::log(x2/x1)
       + (1. - p[iMax])*(x2 - x1)
       + 1./(1. - x2) - 1./(1. - x1)
       + p[iMax]*std::log((1. - x2)/(1. - x1))
       + 0.25*p[0]*(xs1*xs1 - xs2*xs2);
  sum += q;
  if(p.size() == 26) G4cout << "param...  q= " << q <<  " sum= " << sum << G4endl;

  return sum;
}

Here is the call graph for this function:

Here is the caller graph for this function:

MaxEnergyOfSecondaries()

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

virtual

Implements G4VEnergySpectrum.

Definition at line 601 of file G4eIonisationSpectrum.cc.

{
  return 0.5 * kineticEnergy;
}

Here is the caller graph for this function:

operator=()

G4eIonisationSpectrum& G4eIonisationSpectrum::operator= ( const G4eIonisationSpectrum &  right )

private

PrintData()

void G4eIonisationSpectrum::PrintData ( void  ) const

virtual

Implements G4VEnergySpectrum.

Definition at line 596 of file G4eIonisationSpectrum.cc.

{
  theParam->PrintData();
}

Here is the call graph for this function:

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 78 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;
}

Here is the call graph for this function:

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 265 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;
}

Here is the call graph for this function:

Member Data Documentation

factor

G4double G4eIonisationSpectrum::factor

private

Definition at line 113 of file G4eIonisationSpectrum.hh.

iMax

G4int G4eIonisationSpectrum::iMax

private

Definition at line 114 of file G4eIonisationSpectrum.hh.

lowestE

G4double G4eIonisationSpectrum::lowestE

private

Definition at line 112 of file G4eIonisationSpectrum.hh.

theParam

G4eIonisationParameters* G4eIonisationSpectrum::theParam

private

Definition at line 111 of file G4eIonisationSpectrum.hh.

verbose

G4int G4eIonisationSpectrum::verbose

private

Definition at line 115 of file G4eIonisationSpectrum.hh.

---

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

• G4eIonisationSpectrum.hh
• G4eIonisationSpectrum.cc