G4PreCompoundEmission Class Reference

#include <G4PreCompoundEmission.hh>

Collaboration diagram for G4PreCompoundEmission:

Public Member Functions
	G4PreCompoundEmission ()
	~G4PreCompoundEmission ()
void	SetDefaultModel ()
void	SetHETCModel ()
G4ReactionProduct *	PerformEmission (G4Fragment &aFragment)
G4double	GetTotalProbability (const G4Fragment &aFragment)
void	SetOPTxs (G4int)
void	UseSICB (G4bool)

Private Member Functions
void	AngularDistribution (G4VPreCompoundFragment *theFragment, const G4Fragment &aFragment, G4double KineticEnergy)
G4double	rho (G4int p, G4int h, G4double gg, G4double E, G4double Ef) const
	G4PreCompoundEmission (const G4PreCompoundEmission &right)
const G4PreCompoundEmission &	operator= (const G4PreCompoundEmission &right)
G4bool	operator== (const G4PreCompoundEmission &right) const
G4bool	operator!= (const G4PreCompoundEmission &right) const

Private Attributes
G4Pow *	g4pow
G4double	fLevelDensity
G4double	fFermiEnergy
G4PreCompoundFragmentVector *	theFragmentsVector
G4VPreCompoundEmissionFactory *	theFragmentsFactory
G4ThreeVector	theFinalMomentum

Detailed Description

Definition at line 50 of file G4PreCompoundEmission.hh.

Constructor & Destructor Documentation

G4PreCompoundEmission() [1/2]

G4PreCompoundEmission::G4PreCompoundEmission

(

)

Definition at line 57 of file G4PreCompoundEmission.cc.

{
  theFragmentsFactory = new G4PreCompoundEmissionFactory();
  theFragmentsVector = 
    new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
  g4pow = G4Pow::GetInstance();
  G4PreCompoundParameters param;
  fLevelDensity = param.GetLevelDensity();
  fFermiEnergy  = param.GetFermiEnergy();
}

Here is the call graph for this function:

~G4PreCompoundEmission()

G4PreCompoundEmission::~G4PreCompoundEmission

(

)

Definition at line 68 of file G4PreCompoundEmission.cc.

{
  delete theFragmentsFactory; 
  delete theFragmentsVector; 
}

G4PreCompoundEmission() [2/2]

G4PreCompoundEmission::G4PreCompoundEmission ( const G4PreCompoundEmission &  right )

private

Member Function Documentation

AngularDistribution()

void G4PreCompoundEmission::AngularDistribution ( G4VPreCompoundFragment *  theFragment, const G4Fragment &  aFragment, G4double  KineticEnergy  )

private

Definition at line 174 of file G4PreCompoundEmission.cc.

{
  G4int p = aFragment.GetNumberOfParticles();
  G4int h = aFragment.GetNumberOfHoles();
  G4double U = aFragment.GetExcitationEnergy();

  // Emission particle separation energy
  G4double Bemission = thePreFragment->GetBindingEnergy();
    
  G4double gg = (6.0/pi2)*aFragment.GetA_asInt()*fLevelDensity;
    
  // Average exciton energy relative to bottom of nuclear well
  G4double Eav = 2*p*(p+1)/((p+h)*gg);
    
  // Excitation energy relative to the Fermi Level
  G4double Uf = std::max(U - (p - h)*fFermiEnergy , 0.0);
  //  G4double Uf = U - KineticEnergyOfEmittedFragment - Bemission;

  G4double w_num = rho(p+1, h, gg, Uf, fFermiEnergy);
  G4double w_den = rho(p,   h, gg, Uf, fFermiEnergy);
  if (w_num > 0.0 && w_den > 0.0)
    {
      Eav *= (w_num/w_den);
      Eav += - Uf/(p+h) + fFermiEnergy;
    }
  else 
    {
      Eav = fFermiEnergy;
    }
  
  // VI + JMQ 19/01/2010 update computation of the parameter an
  //
  G4double an = 0.0;
  G4double Eeff = ekin + Bemission + fFermiEnergy;
  if(ekin > DBL_MIN && Eeff > DBL_MIN) {

    G4double zeta = std::max(1.0,9.3/std::sqrt(ekin/MeV));
  
    // This should be the projectile energy. If I would know which is 
    // the projectile (proton, neutron) I could remove the binding energy. 
    // But, what happens if INC precedes precompound? This approximation
    // seems to work well enough
    G4double ProjEnergy = aFragment.GetExcitationEnergy();

    an = 3*std::sqrt((ProjEnergy+fFermiEnergy)*Eeff)/(zeta*Eav);

    G4int ne = aFragment.GetNumberOfExcitons() - 1;
    if ( ne > 1 ) { an /= (G4double)ne; }
            
    // protection of exponent
    if ( an > 10. ) { an = 10.; }
  }

  // sample cosine of theta and not theta as in old versions  
  G4double random = G4UniformRand();
  G4double cost;
 
  if(an < 0.1) { cost = 1. - 2*random; }
  else {
    G4double exp2an = G4Exp(-2*an);
    cost = 1. + G4Log(1-random*(1-exp2an))/an;
    if(cost > 1.) { cost = 1.; }
    else if(cost < -1.) {cost = -1.; }
  }  

  G4double phi = CLHEP::twopi*G4UniformRand();
  
  // Calculate the momentum magnitude of emitted fragment   
  G4double pmag = 
    std::sqrt(ekin*(ekin + 2.0*thePreFragment->GetNuclearMass()));
  
  G4double sint = std::sqrt((1.0-cost)*(1.0+cost));

  theFinalMomentum.set(pmag*std::cos(phi)*sint,pmag*std::sin(phi)*sint,
               pmag*cost);

  // theta is the angle wrt the incident direction
  G4ThreeVector theIncidentDirection = aFragment.GetMomentum().vect().unit();
  theFinalMomentum.rotateUz(theIncidentDirection);
}

Here is the call graph for this function:

Here is the caller graph for this function:

GetTotalProbability()

G4double G4PreCompoundEmission::GetTotalProbability ( const G4Fragment &  aFragment )

inline

Definition at line 101 of file G4PreCompoundEmission.hh.

{
  return theFragmentsVector->CalculateProbabilities(aFragment);
}

Here is the call graph for this function:

Here is the caller graph for this function:

operator!=()

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

private

operator=()

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

private

operator==()

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

private

PerformEmission()

G4ReactionProduct * G4PreCompoundEmission::PerformEmission

(

G4Fragment &

aFragment

)

Definition at line 99 of file G4PreCompoundEmission.cc.

{
  // Choose a Fragment for emission
  G4VPreCompoundFragment * thePreFragment = 
    theFragmentsVector->ChooseFragment();
  if (thePreFragment == 0)
    {
      G4cout << "G4PreCompoundEmission::PerformEmission : "
         << "I couldn't choose a fragment\n"
         << "while trying to de-excite\n" 
         << aFragment << G4endl;
      throw G4HadronicException(__FILE__, __LINE__, "");
    }

  //G4cout << "Chosen fragment: " << G4endl;
  //G4cout << *thePreFragment << G4endl;

  // Kinetic Energy of emitted fragment
  G4double kinEnergy = thePreFragment->SampleKineticEnergy(aFragment);
  //  if(kinEnergy < MeV) {
  //  G4cout << "Chosen fragment: " << G4endl;
  //  G4cout << *thePreFragment << G4endl;
  //  G4cout << "Ekin= " << kinEnergy << G4endl;
    // }
  if(kinEnergy < 0.0) { kinEnergy = 0.0; }
  
  // Calculate the fragment momentum (three vector)
  AngularDistribution(thePreFragment,aFragment,kinEnergy);
  
  // Mass of emittef fragment
  G4double EmittedMass = thePreFragment->GetNuclearMass();
  // Now we can calculate the four momentum 
  // both options are valid and give the same result but 2nd one is faster
  G4LorentzVector Emitted4Momentum(theFinalMomentum,EmittedMass + kinEnergy);
    
  // Perform Lorentz boost
  G4LorentzVector Rest4Momentum = aFragment.GetMomentum();
  Emitted4Momentum.boost(Rest4Momentum.boostVector());  

  // Set emitted fragment momentum
  thePreFragment->SetMomentum(Emitted4Momentum);    

  // NOW THE RESIDUAL NUCLEUS
  // ------------------------

  Rest4Momentum -= Emitted4Momentum;
    
  // Update nucleus parameters:
  // --------------------------

  // Z and A
  aFragment.SetZandA_asInt(thePreFragment->GetRestZ(),
               thePreFragment->GetRestA());
    
  // Number of excitons
  aFragment.SetNumberOfParticles(aFragment.GetNumberOfParticles()-
                 thePreFragment->GetA());
  // Number of charges
  aFragment.SetNumberOfCharged(aFragment.GetNumberOfCharged()-
                   thePreFragment->GetZ());
    
  // Update nucleus momentum 
  // A check on consistence of Z, A, and mass will be performed
  aFragment.SetMomentum(Rest4Momentum);
    
  // Create a G4ReactionProduct 
  G4ReactionProduct * MyRP = thePreFragment->GetReactionProduct();

  //  if(kinEnergy < MeV) {
  //  G4cout << "G4PreCompoundEmission::Fragment emitted" << G4endl;
  //  G4cout << *thePreFragment << G4endl;
    // }
  return MyRP;
}

Here is the call graph for this function:

Here is the caller graph for this function:

rho()

G4double G4PreCompoundEmission::rho ( G4int  p, G4int  h, G4double  gg, G4double  E, G4double  Ef  ) const

private

Definition at line 258 of file G4PreCompoundEmission.cc.

{   
  // 25.02.2010 V.Ivanchenko added more protections
  G4double Aph   = (p*p + h*h + p - 3.0*h)/(4.0*gg);
  //  G4double alpha = (p*p + h*h)/(2.0*gg);
  
  if ( E - Aph < 0.0) { return 0.0; }
  
  G4double logConst =  (p+h)*G4Log(gg) 
    - g4pow->logfactorial(p+h-1) - g4pow->logfactorial(p) 
    - g4pow->logfactorial(h);

  // initialise values using j=0

  G4double t1=1;
  G4double t2=1;
  G4double logt3 = (p+h-1) * G4Log(E-Aph) + logConst;
  const G4double logmax = 200.;
  if(logt3 > logmax) { logt3 = logmax; }
  G4double tot = G4Exp( logt3 );

  // and now sum rest of terms
  // 25.02.2010 V.Ivanchenko change while to for loop and cleanup 
  G4double Eeff = E - Aph; 
  for(G4int j=1; j<=h; ++j) 
    {
      Eeff -= Ef;
      if(Eeff < 0.0) { break; }
      t1 *= -1.;
      t2 *= (G4double)(h+1-j)/(G4double)j;
      logt3 = (p+h-1) * G4Log( Eeff) + logConst;
      if(logt3 > logmax) { logt3 = logmax; }
      tot += t1*t2*G4Exp(logt3);
    }
        
  return tot;
}

Here is the call graph for this function:

Here is the caller graph for this function:

SetDefaultModel()

void G4PreCompoundEmission::SetDefaultModel

(

)

Definition at line 74 of file G4PreCompoundEmission.cc.

{
  if (theFragmentsFactory) { delete theFragmentsFactory; }
  theFragmentsFactory = new G4PreCompoundEmissionFactory();
  if (theFragmentsVector) {
    theFragmentsVector->SetVector(theFragmentsFactory->GetFragmentVector());
  } else {
    theFragmentsVector = 
      new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

SetHETCModel()

void G4PreCompoundEmission::SetHETCModel

(

)

Definition at line 86 of file G4PreCompoundEmission.cc.

{
  if (theFragmentsFactory) delete theFragmentsFactory;
  theFragmentsFactory = new G4HETCEmissionFactory();
  if (theFragmentsVector) {
    theFragmentsVector->SetVector(theFragmentsFactory->GetFragmentVector());
  } else {
    theFragmentsVector = 
      new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

SetOPTxs()

void G4PreCompoundEmission::SetOPTxs ( G4int  opt )

inline

Definition at line 106 of file G4PreCompoundEmission.hh.

{
  theFragmentsVector->SetOPTxs(opt);
}

Here is the call graph for this function:

Here is the caller graph for this function:

UseSICB()

void G4PreCompoundEmission::UseSICB ( G4bool  use )

inline

Definition at line 111 of file G4PreCompoundEmission.hh.

{
  theFragmentsVector->UseSICB(use);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

fFermiEnergy

G4double G4PreCompoundEmission::fFermiEnergy

private

Definition at line 90 of file G4PreCompoundEmission.hh.

fLevelDensity

G4double G4PreCompoundEmission::fLevelDensity

private

Definition at line 89 of file G4PreCompoundEmission.hh.

g4pow

G4Pow* G4PreCompoundEmission::g4pow

private

Definition at line 88 of file G4PreCompoundEmission.hh.

theFinalMomentum

G4ThreeVector G4PreCompoundEmission::theFinalMomentum

private

Definition at line 97 of file G4PreCompoundEmission.hh.

theFragmentsFactory

G4VPreCompoundEmissionFactory* G4PreCompoundEmission::theFragmentsFactory

private

Definition at line 94 of file G4PreCompoundEmission.hh.

theFragmentsVector

G4PreCompoundFragmentVector* G4PreCompoundEmission::theFragmentsVector

private

Definition at line 93 of file G4PreCompoundEmission.hh.

---

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

• G4PreCompoundEmission.hh
• G4PreCompoundEmission.cc