#include <G4ParticleHPMadlandNixSpectrum.hh>

Inheritance diagram for G4ParticleHPMadlandNixSpectrum:

Public Member Functions
	G4ParticleHPMadlandNixSpectrum ()

G4double	GetFractionalProbability (G4double anEnergy)

void	Init (std::istream &aDataFile)

G4double	Sample (G4double anEnergy)

	~G4ParticleHPMadlandNixSpectrum ()

Private Member Functions
G4double	E1 (G4double aValue)

G4double	FissionIntegral (G4double tm, G4double anEnergy)

G4double	Gamma05 (G4double aValue)

G4double	Gamma15 (G4double aValue)

G4double	Gamma25 (G4double aValue)

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

G4double	Madland (G4double aSecEnergy, G4double tm)

Private Attributes
G4double	expm1

G4double	theAvarageKineticPerNucleonForHeavyFragments

G4double	theAvarageKineticPerNucleonForLightFragments

G4ParticleHPVector	theFractionalProb

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.);
    theAvarageKineticPerNucleonForLightFragments = 0.0;
    theAvarageKineticPerNucleonForHeavyFragments = 0.0;
  }

References expm1, G4Exp(), theAvarageKineticPerNucleonForHeavyFragments, and theAvarageKineticPerNucleonForLightFragments.

◆ ~G4ParticleHPMadlandNixSpectrum()

G4ParticleHPMadlandNixSpectrum::~G4ParticleHPMadlandNixSpectrum ( )

inline

Definition at line 60 of file G4ParticleHPMadlandNixSpectrum.hh.

61 {

62 }

Member Function Documentation

◆ E1()

G4double G4ParticleHPMadlandNixSpectrum::E1 ( G4double aValue )

inlineprivate

Definition at line 119 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;
  }

References G4Log().

Referenced by GIntegral(), and Madland().

◆ FissionIntegral()

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

inlineprivate

Definition at line 85 of file G4ParticleHPMadlandNixSpectrum.hh.

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

References GIntegral(), theAvarageKineticPerNucleonForHeavyFragments, theAvarageKineticPerNucleonForLightFragments, and G4InuclParticleNames::tm.

Referenced by Sample().

◆ Gamma05()

G4double G4ParticleHPMadlandNixSpectrum::Gamma05 ( G4double aValue )

inlineprivate

Definition at line 93 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;
  }

References G4Exp(), and CLHEP::pi.

Referenced by Gamma15().

◆ Gamma15()

G4double G4ParticleHPMadlandNixSpectrum::Gamma15 ( G4double aValue )

inlineprivate

Definition at line 104 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;
  }

References G4Exp(), and Gamma05().

Referenced by Gamma25(), GIntegral(), and Madland().

◆ Gamma25()

G4double G4ParticleHPMadlandNixSpectrum::Gamma25 ( G4double aValue )

inlineprivate

Definition at line 112 of file G4ParticleHPMadlandNixSpectrum.hh.

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

References G4Exp(), Gamma15(), G4Pow::GetInstance(), and G4Pow::powA().

Referenced by GIntegral().

◆ GetFractionalProbability()

G4double G4ParticleHPMadlandNixSpectrum::GetFractionalProbability ( G4double anEnergy )

inlinevirtual

Implements G4VParticleHPEDis.

Definition at line 74 of file G4ParticleHPMadlandNixSpectrum.hh.

  {
    return theFractionalProb.GetY(anEnergy);
  }

References G4ParticleHPVector::GetY(), and theFractionalProb.

◆ GIntegral()

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

private

Definition at line 111 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;
  }

References A, alpha, alpha2, B(), anonymous_namespace{G4PionRadiativeDecayChannel.cc}::beta, E1(), eV, G4Exp(), Gamma15(), Gamma25(), G4Pow::GetInstance(), G4Pow::powA(), and G4InuclParticleNames::tm.

Referenced by FissionIntegral().

◆ Init()

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

inlinevirtual

Implements G4VParticleHPEDis.

Definition at line 64 of file G4ParticleHPMadlandNixSpectrum.hh.

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

References CLHEP::eV, G4ParticleHPVector::Init(), theAvarageKineticPerNucleonForHeavyFragments, theAvarageKineticPerNucleonForLightFragments, theFractionalProb, and theMaxTemp.

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

References E1(), G4INCL::KinematicsUtils::energy(), eV, Gamma15(), G4Pow::GetInstance(), G4Pow::powA(), theAvarageKineticPerNucleonForHeavyFragments, theAvarageKineticPerNucleonForLightFragments, and G4InuclParticleNames::tm.

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

References FissionIntegral(), G4cout, G4endl, G4UniformRand, G4ParticleHPVector::GetY(), MeV, theMaxTemp, and G4InuclParticleNames::tm.

Field Documentation

◆ expm1

G4double G4ParticleHPMadlandNixSpectrum::expm1

private

Definition at line 147 of file G4ParticleHPMadlandNixSpectrum.hh.

Referenced by G4ParticleHPMadlandNixSpectrum().

◆ theAvarageKineticPerNucleonForHeavyFragments

G4double G4ParticleHPMadlandNixSpectrum::theAvarageKineticPerNucleonForHeavyFragments

private

Definition at line 154 of file G4ParticleHPMadlandNixSpectrum.hh.

Referenced by FissionIntegral(), G4ParticleHPMadlandNixSpectrum(), Init(), and Madland().

◆ theAvarageKineticPerNucleonForLightFragments

G4double G4ParticleHPMadlandNixSpectrum::theAvarageKineticPerNucleonForLightFragments

private

Definition at line 153 of file G4ParticleHPMadlandNixSpectrum.hh.

Referenced by FissionIntegral(), G4ParticleHPMadlandNixSpectrum(), Init(), and Madland().

◆ theFractionalProb

G4ParticleHPVector G4ParticleHPMadlandNixSpectrum::theFractionalProb

private

Definition at line 151 of file G4ParticleHPMadlandNixSpectrum.hh.

Referenced by GetFractionalProbability(), and Init().

◆ theMaxTemp

G4ParticleHPVector G4ParticleHPMadlandNixSpectrum::theMaxTemp

private

Definition at line 156 of file G4ParticleHPMadlandNixSpectrum.hh.

Referenced by Init(), and Sample().

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

source/processes/hadronic/models/particle_hp/include/G4ParticleHPMadlandNixSpectrum.hh
source/processes/hadronic/models/particle_hp/src/G4ParticleHPMadlandNixSpectrum.cc

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ G4ParticleHPMadlandNixSpectrum()

◆ ~G4ParticleHPMadlandNixSpectrum()

Member Function Documentation

◆ E1()

◆ FissionIntegral()

◆ Gamma05()

◆ Gamma15()

◆ Gamma25()

◆ GetFractionalProbability()

◆ GIntegral()

◆ Init()

◆ Madland()

◆ Sample()

Field Documentation

◆ expm1

◆ theAvarageKineticPerNucleonForHeavyFragments

◆ theAvarageKineticPerNucleonForLightFragments

◆ theFractionalProb

◆ theMaxTemp