G4ParticleHPMadlandNixSpectrum Class Reference

#include <G4ParticleHPMadlandNixSpectrum.hh>

Inheritance diagram for G4ParticleHPMadlandNixSpectrum:

Collaboration diagram for G4ParticleHPMadlandNixSpectrum:

Public Member Functions
	G4ParticleHPMadlandNixSpectrum ()
	~G4ParticleHPMadlandNixSpectrum ()
void	Init (std::istream &aDataFile)
G4double	GetFractionalProbability (G4double anEnergy)
G4double	Sample (G4double anEnergy)
Public Member Functions inherited from G4VParticleHPEDis
	G4VParticleHPEDis ()
virtual	~G4VParticleHPEDis ()

Private Member Functions
G4double	Madland (G4double aSecEnergy, G4double tm)
G4double	FissionIntegral (G4double tm, G4double anEnergy)
G4double	GIntegral (G4double tm, G4double anEnergy, G4double aMean)
G4double	Gamma05 (G4double aValue)
G4double	Gamma15 (G4double aValue)
G4double	Gamma25 (G4double aValue)
G4double	E1 (G4double aValue)

Private Attributes
G4double	expm1
G4ParticleHPVector	theFractionalProb
G4double	theAvarageKineticPerNucleonForLightFragments
G4double	theAvarageKineticPerNucleonForHeavyFragments
G4ParticleHPVector	theMaxTemp

Detailed Description

Definition at line 51 of file G4ParticleHPMadlandNixSpectrum.hh.

Constructor & Destructor Documentation

G4ParticleHPMadlandNixSpectrum()

G4ParticleHPMadlandNixSpectrum::G4ParticleHPMadlandNixSpectrum ( )

inline

Definition at line 54 of file G4ParticleHPMadlandNixSpectrum.hh.

  {
    expm1 = G4Exp(-1.);
  }

Here is the call graph for this function:

~G4ParticleHPMadlandNixSpectrum()

G4ParticleHPMadlandNixSpectrum::~G4ParticleHPMadlandNixSpectrum ( )

inline

Definition at line 58 of file G4ParticleHPMadlandNixSpectrum.hh.

  {
  }

Member Function Documentation

E1()

G4double G4ParticleHPMadlandNixSpectrum::E1 ( G4double  aValue )

inlineprivate

Definition at line 117 of file G4ParticleHPMadlandNixSpectrum.hh.

  {
  // good only for rather low aValue @@@ replace by the corresponding NAG function for the
  // exponential integral. (<5 seems ok.
    G4double gamma = 0.577216;
    G4double precision = 0.000001;
    G4double result =-gamma - G4Log(aValue);
    G4double term = -aValue;
    //110527TKDB  Unnessary codes, Detected by gcc4.6 compiler 
    //G4double last;
    G4int count = 1;
    result -= term;
    for(;;)
    {
      count++;
      //110527TKDB  Unnessary codes, Detected by gcc4.6 compiler 
      //last = result;
      term = -term*aValue*(count-1)/(count*count);
      result -=term;
      if(std::fabs(term)/std::fabs(result)<precision) break;
    }
//    NagError *fail; @
//    result = nag_exp_integral(aValue, fail); @
    return result;
  }

Here is the call graph for this function:

Here is the caller graph for this function:

FissionIntegral()

G4double G4ParticleHPMadlandNixSpectrum::FissionIntegral ( G4double  tm, G4double  anEnergy  )

inlineprivate

Definition at line 83 of file G4ParticleHPMadlandNixSpectrum.hh.

  {
    return 0.5*(  GIntegral(tm, anEnergy, theAvarageKineticPerNucleonForLightFragments) 
                       +GIntegral(tm, anEnergy, theAvarageKineticPerNucleonForHeavyFragments) );
  }

Here is the call graph for this function:

Here is the caller graph for this function:

Gamma05()

G4double G4ParticleHPMadlandNixSpectrum::Gamma05 ( G4double  aValue )

inlineprivate

Definition at line 91 of file G4ParticleHPMadlandNixSpectrum.hh.

  {
    G4double result;
    // gamma(1.2,x*X) = std::sqrt(CLHEP::pi)*Erf(x)
    G4double x = std::sqrt(aValue);
    G4double t = 1./(1+0.47047*x);
    result = 1- (0.3480242*t - 0.0958798*t*t + 0.7478556*t*t*t)*G4Exp(-aValue); // @ check
    result *= std::sqrt(CLHEP::pi);
    return result;
  }

Here is the call graph for this function:

Here is the caller graph for this function:

Gamma15()

G4double G4ParticleHPMadlandNixSpectrum::Gamma15 ( G4double  aValue )

inlineprivate

Definition at line 102 of file G4ParticleHPMadlandNixSpectrum.hh.

  {
    G4double result;
    // gamma(a+1, x) = a*gamma(a,x)-x**a*std::exp(-x)
    result = 0.5*Gamma05(aValue) - std::sqrt(aValue)*G4Exp(-aValue); // @ check
    return result;
  }

Here is the call graph for this function:

Here is the caller graph for this function:

Gamma25()

G4double G4ParticleHPMadlandNixSpectrum::Gamma25 ( G4double  aValue )

inlineprivate

Definition at line 110 of file G4ParticleHPMadlandNixSpectrum.hh.

  {
    G4double result;
    result = 1.5*Gamma15(aValue) - G4Pow::GetInstance()->powA(aValue,1.5)*G4Exp(aValue); // @ check
    return result;
  }

Here is the call graph for this function:

Here is the caller graph for this function:

GetFractionalProbability()

G4double G4ParticleHPMadlandNixSpectrum::GetFractionalProbability ( G4double  anEnergy )

inlinevirtual

Implements G4VParticleHPEDis.

Definition at line 72 of file G4ParticleHPMadlandNixSpectrum.hh.

  {
    return theFractionalProb.GetY(anEnergy);
  }

Here is the call graph for this function:

GIntegral()

G4double G4ParticleHPMadlandNixSpectrum::GIntegral ( G4double  tm, G4double  anEnergy, G4double  aMean  )

private

Definition at line 112 of file G4ParticleHPMadlandNixSpectrum.cc.

  {
    G4Pow* Pow=G4Pow::GetInstance();
    if(aMean<1*eV) return 0;
    G4double b = anEnergy/eV;
    G4double sb = std::sqrt(b);
    G4double EF = aMean/eV;
    
    G4double alpha = std::sqrt(tm); 
    G4double beta = std::sqrt(EF);
    G4double A = EF/tm;
    G4double B =  (sb+beta)*(sb+beta)/tm;
    G4double Ap = A;
    G4double Bp = (sb-beta)*(sb-beta)/tm;
    
    G4double result;
    G4double alpha2 = alpha*alpha;
    G4double alphabeta = alpha*beta;
    if(b<EF)
    {
      result =
      (
        (0.4*alpha2*Pow->powA(B,2.5) - 0.5*alphabeta*B*B)*E1(B) -  
        (0.4*alpha2*Pow->powA(A,2.5) - 0.5*alphabeta*A*A)*E1(A) 
      )
       -
      (
        (0.4*alpha2*Pow->powA(Bp,2.5) + 0.5*alphabeta*Bp*Bp)*E1(Bp) -  
        (0.4*alpha2*Pow->powA(Ap,2.5) + 0.5*alphabeta*Ap*Ap)*E1(Ap) 
      )
      +
      (
        (alpha2*B-2*alphabeta*std::sqrt(B))*Gamma15(B)  - 
        (alpha2*A-2*alphabeta*std::sqrt(A))*Gamma15(A) 
      )
      -
      (
        (alpha2*Bp-2*alphabeta*std::sqrt(Bp))*Gamma15(Bp) -
        (alpha2*Ap-2*alphabeta*std::sqrt(Ap))*Gamma15(Ap)
      )
      - 0.6*alpha2*(Gamma25(B) - Gamma25(A) - Gamma25(Bp) + Gamma25(Ap))
      - 1.5*alphabeta*(G4Exp(-B)*(1+B) - G4Exp(-A)*(1+A) + G4Exp(-Bp)*(1+Bp) + G4Exp(-Ap)*(1+Ap)) ;
    }
    else
    {
      result =
      (
        (0.4*alpha2*Pow->powA(B,2.5) - 0.5*alphabeta*B*B)*E1(B) -  
        (0.4*alpha2*Pow->powA(A,2.5) - 0.5*alphabeta*A*A)*E1(A) 
      );
       result -=
      (
        (0.4*alpha2*Pow->powA(Bp,2.5) + 0.5*alphabeta*Bp*Bp)*E1(Bp) -  
        (0.4*alpha2*Pow->powA(Ap,2.5) + 0.5*alphabeta*Ap*Ap)*E1(Ap) 
      );
      result +=
      (
        (alpha2*B-2*alphabeta*std::sqrt(B))*Gamma15(B)  - 
        (alpha2*A-2*alphabeta*std::sqrt(A))*Gamma15(A) 
      );
      result -=
      (
        (alpha2*Bp+2*alphabeta*std::sqrt(Bp))*Gamma15(Bp) -
        (alpha2*Ap+2*alphabeta*std::sqrt(Ap))*Gamma15(Ap)
      );
      result -= 0.6*alpha2*(Gamma25(B) - Gamma25(A) - Gamma25(Bp) + Gamma25(Ap));
      result -= 1.5*alphabeta*(G4Exp(-B)*(1+B) - G4Exp(-A)*(1+A) + G4Exp(-Bp)*(1+Bp) + G4Exp(-Ap)*(1+Ap) - 2.) ;
    }
    result = result / (3.*std::sqrt(tm*EF));
    return result;
  }

Here is the call graph for this function:

Here is the caller graph for this function:

Init()

void G4ParticleHPMadlandNixSpectrum::Init ( std::istream &  aDataFile )

inlinevirtual

Implements G4VParticleHPEDis.

Definition at line 62 of file G4ParticleHPMadlandNixSpectrum.hh.

  {
    theFractionalProb.Init(aDataFile);
    aDataFile>> theAvarageKineticPerNucleonForLightFragments;
    theAvarageKineticPerNucleonForLightFragments*=CLHEP::eV;
    aDataFile>> theAvarageKineticPerNucleonForHeavyFragments;
    theAvarageKineticPerNucleonForHeavyFragments*=CLHEP::eV;
    theMaxTemp.Init(aDataFile);
  }

Here is the call graph for this function:

Madland()

G4double G4ParticleHPMadlandNixSpectrum::Madland ( G4double  aSecEnergy, G4double  tm  )

private

Definition at line 34 of file G4ParticleHPMadlandNixSpectrum.cc.

  {
    G4Pow* Pow=G4Pow::GetInstance();
    G4double result;
    G4double energy = aSecEnergy/eV;
    G4double EF;
    
    EF = theAvarageKineticPerNucleonForLightFragments/eV;
    G4double lightU1 = std::sqrt(energy)-std::sqrt(EF);
    lightU1 *= lightU1/tm;
    G4double lightU2 = std::sqrt(energy)+std::sqrt(EF);
    lightU2 *= lightU2/tm;
    G4double lightTerm=0;
    if(theAvarageKineticPerNucleonForLightFragments>1*eV)
    {
      lightTerm  = Pow->powA(lightU2, 1.5)*E1(lightU2);
      lightTerm -= Pow->powA(lightU1, 1.5)*E1(lightU1);
      lightTerm += Gamma15(lightU2)-Gamma15(lightU1);
      lightTerm /= 3.*std::sqrt(tm*EF);
    }
    
    EF = theAvarageKineticPerNucleonForHeavyFragments/eV;
    G4double heavyU1 = std::sqrt(energy)-std::sqrt(EF);
    heavyU1 *= heavyU1/tm;
    G4double heavyU2 = std::sqrt(energy)+std::sqrt(EF);
    heavyU2 *= heavyU2/tm;
    G4double heavyTerm=0    ;
    if(theAvarageKineticPerNucleonForHeavyFragments> 1*eV)
    {
      heavyTerm  = Pow->powA(heavyU2, 1.5)*E1(heavyU2);
      heavyTerm -= Pow->powA(heavyU1, 1.5)*E1(heavyU1);
      heavyTerm += Gamma15(heavyU2)-Gamma15(heavyU1);
      heavyTerm /= 3.*std::sqrt(tm*EF);
    }
    
    result = 0.5*(lightTerm+heavyTerm);
    
    return result;
  }

Here is the call graph for this function:

Here is the caller graph for this function:

Sample()

G4double G4ParticleHPMadlandNixSpectrum::Sample ( G4double  anEnergy )

virtual

Implements G4VParticleHPEDis.

Definition at line 74 of file G4ParticleHPMadlandNixSpectrum.cc.

  {
    G4double tm = theMaxTemp.GetY(anEnergy);
    G4double last=0, buff, current = 100*MeV;
    G4double precision = 0.001;
    G4double newValue = 0., oldValue=0.;
    G4double random = G4UniformRand();
    
    G4int icounter=0;
    G4int icounter_max=1024;
    do
    {
      icounter++;
      if ( icounter > icounter_max ) {
     G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
         break;
      }
      oldValue = newValue;
      newValue = FissionIntegral(tm, current);
      if(newValue < random) 
      {
        buff = current;
    current+=std::abs(current-last)/2.;
    last = buff;
        if(current>190*MeV) throw G4HadronicException(__FILE__, __LINE__, "Madland-Nix Spectrum has not converged in sampling");
      }
      else
      {
        buff = current;
        current-=std::abs(current-last)/2.;
    last = buff;
      }
    }
    while (std::abs(oldValue-newValue)>precision*newValue); // Loop checking, 11.05.2015, T. Koi
    return current;
  }

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

expm1

G4double G4ParticleHPMadlandNixSpectrum::expm1

private

Definition at line 145 of file G4ParticleHPMadlandNixSpectrum.hh.

theAvarageKineticPerNucleonForHeavyFragments

G4double G4ParticleHPMadlandNixSpectrum::theAvarageKineticPerNucleonForHeavyFragments

private

Definition at line 152 of file G4ParticleHPMadlandNixSpectrum.hh.

theAvarageKineticPerNucleonForLightFragments

G4double G4ParticleHPMadlandNixSpectrum::theAvarageKineticPerNucleonForLightFragments

private

Definition at line 151 of file G4ParticleHPMadlandNixSpectrum.hh.

theFractionalProb

G4ParticleHPVector G4ParticleHPMadlandNixSpectrum::theFractionalProb

private

Definition at line 149 of file G4ParticleHPMadlandNixSpectrum.hh.

theMaxTemp

G4ParticleHPVector G4ParticleHPMadlandNixSpectrum::theMaxTemp

private

Definition at line 154 of file G4ParticleHPMadlandNixSpectrum.hh.

---

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

• G4ParticleHPMadlandNixSpectrum.hh
• G4ParticleHPMadlandNixSpectrum.cc