#include <G4EvaporationProbability.hh>

Inheritance diagram for G4EvaporationProbability:

Collaboration diagram for G4EvaporationProbability:

Public Member Functions
	G4EvaporationProbability (G4int anA, G4int aZ, G4double aGamma, G4VCoulombBarrier *)

virtual	~G4EvaporationProbability ()

G4int	GetZ (void) const

G4int	GetA (void) const

G4double	EmissionProbability (const G4Fragment &fragment, G4double maxKineticEnergy)

G4double	TotalProbability (const G4Fragment &fragment, G4double minKineticEnergy, G4double maxKineticEnergy)

G4double	ProbabilityDistributionFunction (G4double K)

G4double	SampleKineticEnergy (G4double minKineticEnergy, G4double maxKineticEnergy)

Public Member Functions inherited from G4VEmissionProbability
	G4VEmissionProbability ()

virtual	~G4VEmissionProbability ()

void	SetOPTxs (G4int opt)

void	UseSICB (G4bool use)

Protected Member Functions
virtual G4double	CalcAlphaParam (const G4Fragment &fragment)=0

virtual G4double	CalcBetaParam (const G4Fragment &fragment)=0

Private Member Functions
G4double	IntegrateEmissionProbability (G4double low, G4double up)

G4double	CrossSection (G4double K)

	G4EvaporationProbability (const G4EvaporationProbability &right)

const G4EvaporationProbability &	operator= (const G4EvaporationProbability &right)

G4bool	operator== (const G4EvaporationProbability &right) const

G4bool	operator!= (const G4EvaporationProbability &right) const

Private Attributes
G4int	theA

G4int	theZ

G4int	fragA

G4int	fragZ

G4int	resA

G4int	resZ

G4int	index

G4int	nbins

G4double	resA13

G4double	muu

G4double	partMass

G4double	resMass

G4double	fragMass

G4double	U

G4double	delta0

G4double	delta1

G4double	a0

G4double	Gamma

G4double	probability [11]

Additional Inherited Members
Protected Attributes inherited from G4VEmissionProbability
G4int	OPTxs

G4bool	useSICB

G4Pow *	fG4pow

G4PairingCorrection *	fPairCorr

G4EvaporationLevelDensityParameter *	theEvapLDPptr

Detailed Description

Definition at line 42 of file G4EvaporationProbability.hh.

Constructor & Destructor Documentation

◆ G4EvaporationProbability() [1/2]

G4EvaporationProbability::G4EvaporationProbability	(	G4int	anA,
		G4int	aZ,
		G4double	aGamma,
		G4VCoulombBarrier *
	)

Definition at line 56 of file G4EvaporationProbability.cc.

   : theA(anA),
     theZ(aZ),
     Gamma(aGamma)
 {
   resZ = resA = fragA = fragZ = nbins = 0;
   resA13 = muu = resMass = fragMass = U = delta0 = delta1 = a0 = 0.0;
   partMass = G4NucleiProperties::GetNuclearMass(theA, theZ);
   index = 0;
   if(1 == theZ) { index = theA; }
   else { index = theA + 1; }
   for(G4int i=0; i<11; ++i) { probability[i] = 0.0; }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ ~G4EvaporationProbability()

G4EvaporationProbability::~G4EvaporationProbability ( )

virtual

Definition at line 72 of file G4EvaporationProbability.cc.

73 {}

◆ G4EvaporationProbability() [2/2]

G4EvaporationProbability::G4EvaporationProbability ( const G4EvaporationProbability & right )

private

Member Function Documentation

◆ CalcAlphaParam()

virtual G4double G4EvaporationProbability::CalcAlphaParam ( const G4Fragment & fragment )

protectedpure virtual

Implemented in G4DeuteronEvaporationProbability, G4NeutronEvaporationProbability, G4TritonEvaporationProbability, G4AlphaEvaporationProbability, G4He3EvaporationProbability, and G4ProtonEvaporationProbability.

Here is the caller graph for this function:

◆ CalcBetaParam()

virtual G4double G4EvaporationProbability::CalcBetaParam ( const G4Fragment & fragment )

protectedpure virtual

Implemented in G4DeuteronEvaporationProbability, G4NeutronEvaporationProbability, G4TritonEvaporationProbability, G4AlphaEvaporationProbability, G4He3EvaporationProbability, and G4ProtonEvaporationProbability.

Here is the caller graph for this function:

◆ CrossSection()

G4double G4EvaporationProbability::CrossSection ( G4double K )

private

Definition at line 193 of file G4EvaporationProbability.cc.

 {
   G4double res;
   if(OPTxs <= 2) { 
     res = G4ChatterjeeCrossSection::ComputeCrossSection(K, resA13, muu,
                             index, theZ, 
                             resZ, resA); 
   } else { 
     res = G4KalbachCrossSection::ComputeCrossSection(K, resA13, muu,
                              index, theZ, theA, 
                              resZ, resA);
   }
   return res;
 }  

Here is the call graph for this function:

Here is the caller graph for this function:

◆ EmissionProbability()

G4double G4EvaporationProbability::EmissionProbability	(	const G4Fragment &	fragment,
		G4double	maxKineticEnergy
	)

virtual

Implements G4VEmissionProbability.

Definition at line 76 of file G4EvaporationProbability.cc.

 {
   return 0.0;
 }

Here is the caller graph for this function:

◆ GetA()

G4int G4EvaporationProbability::GetA ( void ) const

inline

Definition at line 53 of file G4EvaporationProbability.hh.

53 { return theA; }

G4EvaporationProbability::theA

G4int theA

Definition: G4EvaporationProbability.hh:88

Here is the call graph for this function:

Here is the caller graph for this function:

◆ GetZ()

G4int G4EvaporationProbability::GetZ ( void ) const

inline

Definition at line 51 of file G4EvaporationProbability.hh.

51 { return theZ; }

G4EvaporationProbability::theZ

G4int theZ

Definition: G4EvaporationProbability.hh:89

Here is the caller graph for this function:

◆ IntegrateEmissionProbability()

G4double G4EvaporationProbability::IntegrateEmissionProbability	(	G4double	low,
		G4double	up
	)

private

Definition at line 146 of file G4EvaporationProbability.cc.

 {
   G4double sum = 0.0;
   G4double del = up - low;
   nbins = std::min(10, G4int(del*invmev));
   nbins = std::max(nbins, 2);
   del /= G4double(nbins);
   G4double T = low - del*0.5;
   for(G4int i=1; i<=nbins; ++i) {
     T += del;
     sum +=  ProbabilityDistributionFunction(T);
     probability[i] = sum;
   }
   sum *= del;
   return sum;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ operator!=()

G4bool G4EvaporationProbability::operator!= ( const G4EvaporationProbability & right ) const

private

Here is the caller graph for this function:

◆ operator=()

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

private

Here is the caller graph for this function:

◆ operator==()

G4bool G4EvaporationProbability::operator== ( const G4EvaporationProbability & right ) const

private

Here is the caller graph for this function:

◆ ProbabilityDistributionFunction()

G4double G4EvaporationProbability::ProbabilityDistributionFunction ( G4double K )

Definition at line 163 of file G4EvaporationProbability.cc.

 { 
   //G4cout << "### G4EvaporationProbability::ProbabilityDistributionFunction" 
   //  << G4endl;
 
   G4double E0 = U - delta0;
   G4double E1 = U + fragMass - partMass - resMass - delta1 - K;
   /*
   G4cout << "PDF: FragZ= " << fragment.GetZ_asInt() << " FragA= " 
      << fragment.GetA_asInt()
      << " Z= " << theZ << "  A= " << theA 
      << " K= " << K << " E0= " << E0 << " E1= " << E1 << G4endl;
   */
   if(E1 < 0.0) { return 0.0; }
 
   G4double a1 = theEvapLDPptr->LevelDensityParameter(resA, resZ, E1);
 
   //JMQ 14/02/09 BUG fixed: hbarc should be in the denominator instead 
   //             of hbar_Planck; without 1/hbar_Panck remains as a width
 
   static const G4double pcoeff = millibarn/((pi*hbarc)*(pi*hbarc)); 
 
   // Fixed numerical problem
   G4double Prob = pcoeff*Gamma*partMass
     *G4Exp(2*(std::sqrt(a1*E1) - std::sqrt(a0*E0)))*K*CrossSection(K);
 
   return Prob;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ SampleKineticEnergy()

G4double G4EvaporationProbability::SampleKineticEnergy	(	G4double	minKineticEnergy,
		G4double	maxKineticEnergy
	)

Definition at line 209 of file G4EvaporationProbability.cc.

 {
   if(maxKinEnergy <= minKinEnergy) { return 0.0; }
 
   G4double T = 0.0;
   if (OPTxs==0) {
     // JMQ:
     // It uses Dostrovsky's approximation for the inverse reaction cross
     // in the probability for fragment emission
     // MaximalKineticEnergy energy in the original version (V.Lara) was 
     // calculated at the Coulomb barrier.
     
     G4double Rb = 4.0*a0*maxKinEnergy;
     G4double RbSqrt = std::sqrt(Rb);
     G4double PEX1 = 0.0;
     if (RbSqrt < 160.0) PEX1 = G4Exp(-RbSqrt);
     G4double Rk = 0.0;
     G4double FRk = 0.0;
     do {
       G4double RandNumber = G4UniformRand();
       Rk = 1.0 + (1./RbSqrt)*G4Log(RandNumber + (1.0-RandNumber)*PEX1);
       G4double Q1 = 1.0;
       G4double Q2 = 1.0;
       if (theZ == 0) { // for emitted neutron
         G4double Beta = (2.12/(resA*resA) - 0.05)*MeV/(0.76 + 2.2/resA13);
         Q1 = 1.0 + Beta/maxKinEnergy;
         Q2 = Q1*std::sqrt(Q1);
       } 
       
       FRk = ssqr3 * Rk * (Q1 - Rk*Rk)/Q2;
       
       // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
     } while (FRk < G4UniformRand());
     
     T = maxKinEnergy * (1.0-Rk*Rk) + minKinEnergy;
 
   } else { 
 
     G4double p =  probability[nbins]*G4UniformRand();
     G4int i = 0;
     for(; i<nbins; ++i) { if(p <= probability[i+1]) { break; } }
     G4double delta = (maxKinEnergy - minKinEnergy)/G4double(nbins);
     T = minKinEnergy + delta*i 
       + delta*(p - probability[i])/(probability[i+1] - probability[i]);
 
   }
   //G4cout << " T= " << T << G4endl;
   return T;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ TotalProbability()

G4double G4EvaporationProbability::TotalProbability	(	const G4Fragment &	fragment,
		G4double	minKineticEnergy,
		G4double	maxKineticEnergy
	)

Definition at line 82 of file G4EvaporationProbability.cc.

 {
   if (maxEnergy <= minEnergy) { return 0.0; }
 
   fragA = fragment.GetA_asInt();
   fragZ = fragment.GetZ_asInt();
   resA  = fragA - theA;
   resZ  = fragZ - theZ;
 
   U = fragment.GetExcitationEnergy();
   delta0 = std::max(0.0, fPairCorr->GetPairingCorrection(fragA,fragZ));
   if(U < delta0) { return 0.0; }
 
   fragMass = fragment.GetGroundStateMass();
   delta1 = std::max(0.0, fPairCorr->GetPairingCorrection(resA,resZ));
   resMass = G4NucleiProperties::GetNuclearMass(resA, resZ);
   resA13 = fG4pow->Z13(resA);
    
   G4double Width = 0.0;
   if (OPTxs==0) {
 
     G4double SystemEntropy = 2.0*std::sqrt(
       theEvapLDPptr->LevelDensityParameter(fragA,fragZ,U)*(U-delta0));
                                   
     static const G4double RN2 = 
       2.25*fermi*fermi/(twopi* hbar_Planck*hbar_Planck);
 
     G4double Alpha = CalcAlphaParam(fragment);
     G4double Beta = CalcBetaParam(fragment);
     
     a0 = theEvapLDPptr->LevelDensityParameter(resA,resZ,maxEnergy);
     G4double GlobalFactor = Gamma*Alpha*partMass*RN2*resA13*resA13/(a0*a0);
     
     G4double maxea = maxEnergy*a0;
     G4double Term1 = Beta*a0 - 1.5 + maxea;
     G4double Term2 = (2.0*Beta*a0-3.0)*std::sqrt(maxea) + 2*maxea;
     
     G4double ExpTerm1 = 0.0;
     if (SystemEntropy <= explim) { ExpTerm1 = G4Exp(-SystemEntropy); }
     
     G4double ExpTerm2 = 2.*std::sqrt(maxea) - SystemEntropy;
     ExpTerm2 = std::min(ExpTerm2, explim);
     ExpTerm2 = G4Exp(ExpTerm2);
     
     Width = GlobalFactor*(Term1*ExpTerm1 + Term2*ExpTerm2);
              
   } else {
 
     a0 = theEvapLDPptr->LevelDensityParameter(fragA,fragZ,U - delta0);
     // compute power once
     if(OPTxs <= 2) { 
       muu =  G4ChatterjeeCrossSection::ComputePowerParameter(resA, index);
     } else {
       muu = G4KalbachCrossSection::ComputePowerParameter(resA, index);
     }
     // if Coulomb barrier cutoff is superimposed for all cross sections 
     // then the limit is the Coulomb Barrier
     Width = IntegrateEmissionProbability(minEnergy, maxEnergy);
   }
   return Width;
 }