G4PreCompoundEmission Class Reference

#include <G4PreCompoundEmission.hh>

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

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

Private Attributes
G4double	fFermiEnergy
G4int	fModelID
G4NuclearLevelData *	fNuclData
G4bool	fUseAngularGenerator
G4Pow *	g4calc
G4ThreeVector	theFinalMomentum
G4VPreCompoundEmissionFactory *	theFragmentsFactory
G4PreCompoundFragmentVector *	theFragmentsVector

Detailed Description

Definition at line 50 of file G4PreCompoundEmission.hh.

Constructor & Destructor Documentation

G4PreCompoundEmission() [1/2]

G4PreCompoundEmission::G4PreCompoundEmission

(

)

Definition at line 60 of file G4PreCompoundEmission.cc.

{
  theFragmentsFactory = new G4PreCompoundEmissionFactory();
  theFragmentsVector = 
    new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
  g4calc = G4Pow::GetInstance();
  fNuclData = G4NuclearLevelData::GetInstance();
  G4DeexPrecoParameters* param = fNuclData->GetParameters();
  fFermiEnergy  = param->GetFermiEnergy();
  fUseAngularGenerator = param->UseAngularGen();
  fModelID = G4PhysicsModelCatalog::GetModelID("model_PRECO");
}

References fFermiEnergy, fModelID, fNuclData, fUseAngularGenerator, g4calc, G4DeexPrecoParameters::GetFermiEnergy(), G4VPreCompoundEmissionFactory::GetFragmentVector(), G4Pow::GetInstance(), G4NuclearLevelData::GetInstance(), G4PhysicsModelCatalog::GetModelID(), G4NuclearLevelData::GetParameters(), theFragmentsFactory, theFragmentsVector, and G4DeexPrecoParameters::UseAngularGen().

~G4PreCompoundEmission()

G4PreCompoundEmission::~G4PreCompoundEmission

(

)

Definition at line 73 of file G4PreCompoundEmission.cc.

{
  delete theFragmentsFactory; 
  delete theFragmentsVector; 
}

References theFragmentsFactory, and 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 180 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)*fNuclData->GetLevelDensity(aFragment.GetZ_asInt(),
                                                     aFragment.GetA_asInt(),U);
    
  // 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/CLHEP::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 /= static_cast<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);
}

References G4InuclParticleNames::an, DBL_MIN, fFermiEnergy, fNuclData, G4Exp(), G4Log(), G4UniformRand, G4Fragment::GetA_asInt(), G4VPreCompoundFragment::GetBindingEnergy(), G4Fragment::GetExcitationEnergy(), G4NuclearLevelData::GetLevelDensity(), G4Fragment::GetMomentum(), G4VPreCompoundFragment::GetNuclearMass(), G4Fragment::GetNumberOfExcitons(), G4Fragment::GetNumberOfHoles(), G4Fragment::GetNumberOfParticles(), G4Fragment::GetZ_asInt(), G4INCL::Math::max(), CLHEP::MeV, pi2, rho(), CLHEP::Hep3Vector::rotateUz(), CLHEP::Hep3Vector::set(), theFinalMomentum, CLHEP::twopi, CLHEP::Hep3Vector::unit(), and CLHEP::HepLorentzVector::vect().

Referenced by PerformEmission().

GetTotalProbability()

G4double G4PreCompoundEmission::GetTotalProbability ( const G4Fragment &  aFragment )

inline

Definition at line 105 of file G4PreCompoundEmission.hh.

{
  return theFragmentsVector->CalculateProbabilities(aFragment);
}

References G4PreCompoundFragmentVector::CalculateProbabilities(), and theFragmentsVector.

Referenced by G4PreCompoundModel::DeExcite().

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 104 of file G4PreCompoundEmission.cc.

{
  // Choose a Fragment for emission
  G4VPreCompoundFragment * thePreFragment = 
    theFragmentsVector->ChooseFragment();
  if (thePreFragment == nullptr)
    {
      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);
  kinEnergy = std::max(kinEnergy, 0.0);
  
  // Calculate the fragment momentum (three vector)
  if(fUseAngularGenerator) {
    AngularDistribution(thePreFragment,aFragment,kinEnergy);
  } else {
    G4double pmag = 
      std::sqrt(kinEnergy*(kinEnergy + 2.0*thePreFragment->GetNuclearMass()));
    theFinalMomentum = pmag*G4RandomDirection();
  }
 
  // 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();
 
  // Set the creator model ID
  aFragment.SetCreatorModelID(fModelID);
  if (MyRP != nullptr) MyRP->SetCreatorModelID(fModelID);
 
  return MyRP;
}

References AngularDistribution(), CLHEP::HepLorentzVector::boost(), CLHEP::HepLorentzVector::boostVector(), G4PreCompoundFragmentVector::ChooseFragment(), fModelID, fUseAngularGenerator, G4cout, G4endl, G4RandomDirection(), G4VPreCompoundFragment::GetA(), G4Fragment::GetMomentum(), G4VPreCompoundFragment::GetNuclearMass(), G4Fragment::GetNumberOfCharged(), G4Fragment::GetNumberOfParticles(), G4VPreCompoundFragment::GetReactionProduct(), G4VPreCompoundFragment::GetRestA(), G4VPreCompoundFragment::GetRestZ(), G4VPreCompoundFragment::GetZ(), G4INCL::Math::max(), G4VPreCompoundFragment::SampleKineticEnergy(), G4ReactionProduct::SetCreatorModelID(), G4Fragment::SetCreatorModelID(), G4VPreCompoundFragment::SetMomentum(), G4Fragment::SetMomentum(), G4Fragment::SetNumberOfCharged(), G4Fragment::SetNumberOfParticles(), G4Fragment::SetZandA_asInt(), theFinalMomentum, and theFragmentsVector.

Referenced by G4PreCompoundModel::DeExcite().

rho()

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

private

Definition at line 265 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);
  
  if ( E - Aph < 0.0) { return 0.0; }
  
  G4double logConst =  (p+h)*G4Log(gg) 
    - g4calc->logfactorial(p+h-1) - g4calc->logfactorial(p) 
    - g4calc->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 *= static_cast<G4double>(h+1-j)/static_cast<G4double>(j);
      logt3 = (p+h-1) * G4Log( Eeff) + logConst;
      if(logt3 > logmax) { logt3 = logmax; }
      tot += t1*t2*G4Exp(logt3);
    }
        
  return tot;
}

References g4calc, G4Exp(), G4Log(), and G4Pow::logfactorial().

Referenced by AngularDistribution().

SetDefaultModel()

void G4PreCompoundEmission::SetDefaultModel

(

)

Definition at line 79 of file G4PreCompoundEmission.cc.

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

References G4VPreCompoundEmissionFactory::GetFragmentVector(), G4PreCompoundFragmentVector::SetVector(), theFragmentsFactory, and theFragmentsVector.

SetHETCModel()

void G4PreCompoundEmission::SetHETCModel

(

)

Definition at line 91 of file G4PreCompoundEmission.cc.

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

References G4VPreCompoundEmissionFactory::GetFragmentVector(), G4PreCompoundFragmentVector::SetVector(), theFragmentsFactory, and theFragmentsVector.

Referenced by G4PreCompoundModel::InitialiseModel().

SetOPTxs()

void G4PreCompoundEmission::SetOPTxs ( G4int  opt )

inline

Definition at line 110 of file G4PreCompoundEmission.hh.

{
  theFragmentsVector->SetOPTxs(opt);
}

References G4PreCompoundFragmentVector::SetOPTxs(), and theFragmentsVector.

Referenced by G4PreCompoundModel::InitialiseModel().

UseSICB()

void G4PreCompoundEmission::UseSICB ( G4bool  use )

inline

Definition at line 115 of file G4PreCompoundEmission.hh.

{
  theFragmentsVector->UseSICB(use);
}

References theFragmentsVector, and G4PreCompoundFragmentVector::UseSICB().

Field Documentation

fFermiEnergy

G4double G4PreCompoundEmission::fFermiEnergy

private

Definition at line 91 of file G4PreCompoundEmission.hh.

Referenced by AngularDistribution(), and G4PreCompoundEmission().

fModelID

G4int G4PreCompoundEmission::fModelID

private

Definition at line 101 of file G4PreCompoundEmission.hh.

Referenced by G4PreCompoundEmission(), and PerformEmission().

fNuclData

G4NuclearLevelData* G4PreCompoundEmission::fNuclData

private

Definition at line 89 of file G4PreCompoundEmission.hh.

Referenced by AngularDistribution(), and G4PreCompoundEmission().

fUseAngularGenerator

G4bool G4PreCompoundEmission::fUseAngularGenerator

private

Definition at line 99 of file G4PreCompoundEmission.hh.

Referenced by G4PreCompoundEmission(), and PerformEmission().

g4calc

G4Pow* G4PreCompoundEmission::g4calc

private

Definition at line 88 of file G4PreCompoundEmission.hh.

Referenced by G4PreCompoundEmission(), and rho().

theFinalMomentum

G4ThreeVector G4PreCompoundEmission::theFinalMomentum

private

Definition at line 98 of file G4PreCompoundEmission.hh.

Referenced by AngularDistribution(), and PerformEmission().

theFragmentsFactory

G4VPreCompoundEmissionFactory* G4PreCompoundEmission::theFragmentsFactory

private

Definition at line 95 of file G4PreCompoundEmission.hh.

Referenced by G4PreCompoundEmission(), SetDefaultModel(), SetHETCModel(), and ~G4PreCompoundEmission().

theFragmentsVector

G4PreCompoundFragmentVector* G4PreCompoundEmission::theFragmentsVector

private

Definition at line 94 of file G4PreCompoundEmission.hh.

Referenced by G4PreCompoundEmission(), GetTotalProbability(), PerformEmission(), SetDefaultModel(), SetHETCModel(), SetOPTxs(), UseSICB(), and ~G4PreCompoundEmission().

---

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

• source/processes/hadronic/models/pre_equilibrium/exciton_model/include/G4PreCompoundEmission.hh
• source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundEmission.cc