G4WilsonAblationModel Class Reference

#include <G4WilsonAblationModel.hh>

Inheritance diagram for G4WilsonAblationModel:

Collaboration diagram for G4WilsonAblationModel:

Public Types
typedef std::vector< G4ParticleDefinition * >	VectorOfFragmentTypes

Public Member Functions
	G4WilsonAblationModel ()
virtual	~G4WilsonAblationModel ()
G4FragmentVector *	BreakItUp (const G4Fragment &theNucleus)
void	SetProduceSecondaries (G4bool)
G4bool	GetProduceSecondaries ()
void	SetVerboseLevel (G4int)
G4int	GetVerboseLevel ()
Public Member Functions inherited from G4VEvaporation
	G4VEvaporation ()
virtual	~G4VEvaporation ()
virtual void	BreakFragment (G4FragmentVector *, G4Fragment *theNucleus)
virtual void	InitialiseChannels ()
virtual void	SetPhotonEvaporation (G4VEvaporationChannel *ptr)
G4VEvaporationChannel *	GetPhotonEvaporation ()
G4VEvaporationChannel *	GetFissionChannel ()
void	SetOPTxs (G4int opt)
void	UseSICB (G4bool use)

Private Member Functions
void	SelectSecondariesByEvaporation (G4Fragment *)
void	SelectSecondariesByDefault (G4ThreeVector)
void	PrintWelcomeMessage ()

Private Attributes
G4bool	produceSecondaries
G4int	verboseLevel
G4double	B
G4int	nFragTypes
G4ParticleDefinition *	fragType [6]
G4FragmentVector *	fragmentVector
VectorOfFragmentTypes	evapType
G4double	fSig [200]

Additional Inherited Members
Protected Member Functions inherited from G4VEvaporation
void	CleanChannels ()
Protected Attributes inherited from G4VEvaporation
G4VEvaporationChannel *	thePhotonEvaporation
G4int	OPTxs
G4bool	useSICB
std::vector< G4VEvaporationChannel * > *	theChannels
G4VEvaporationFactory *	theChannelFactory

Detailed Description

Definition at line 84 of file G4WilsonAblationModel.hh.

Member Typedef Documentation

VectorOfFragmentTypes

typedef std::vector<G4ParticleDefinition*> G4WilsonAblationModel::VectorOfFragmentTypes

Definition at line 90 of file G4WilsonAblationModel.hh.

Constructor & Destructor Documentation

G4WilsonAblationModel()

G4WilsonAblationModel::G4WilsonAblationModel

(

)

Definition at line 113 of file G4WilsonAblationModel.cc.

{
//
//
// Send message to stdout to advise that the G4Abrasion model is being used.
//
  PrintWelcomeMessage();
//
//
// Set the default verbose level to 0 - no output.
//
  verboseLevel = 0;  
//
//
// Set the binding energy per nucleon .... did I mention that this is a crude
// model for nuclear de-excitation?
//
  B = 10.0 * MeV;
//
//
// It is possuble to switch off secondary particle production (other than the
// final nuclear fragment).  The default is on.
//
  produceSecondaries = true;
//
//
// Now we need to define the decay modes.  We're using the G4Evaporation model
// to help determine the kinematics of the decay.
//
  nFragTypes  = 6;
  fragType[0] = G4Alpha::Alpha();
  fragType[1] = G4He3::He3();
  fragType[2] = G4Triton::Triton();
  fragType[3] = G4Deuteron::Deuteron();
  fragType[4] = G4Proton::Proton();
  fragType[5] = G4Neutron::Neutron();
  for(G4int i=0; i<200; ++i) { fSig[i] = 0.0; }
//
//
// Set verboseLevel default to no output.
//
  verboseLevel = 0;
  theChannelFactory = new G4EvaporationFactory(new G4PhotonEvaporation());
  theChannels = theChannelFactory->GetChannel();
//
//
// Set defaults for evaporation classes.  These can be overridden by user
// "set" methods.
//
  OPTxs   = 3;
  useSICB = false;
  fragmentVector = 0;
}

Here is the call graph for this function:

~G4WilsonAblationModel()

G4WilsonAblationModel::~G4WilsonAblationModel ( )

virtual

Definition at line 168 of file G4WilsonAblationModel.cc.

{}

Here is the call graph for this function:

Member Function Documentation

BreakItUp()

G4FragmentVector * G4WilsonAblationModel::BreakItUp ( const G4Fragment &  theNucleus )

virtual

Implements G4VEvaporation.

Definition at line 174 of file G4WilsonAblationModel.cc.

{
//
//
// Initilise the pointer to the G4FragmentVector used to return the information
// about the breakup.
//
  fragmentVector = new G4FragmentVector;
  fragmentVector->clear();
//
//
// Get the A, Z and excitation of the nucleus.
//
  G4int A     = theNucleus.GetA_asInt();
  G4int Z     = theNucleus.GetZ_asInt();
  G4double ex = theNucleus.GetExcitationEnergy();
  if (verboseLevel >= 2)
  {
    G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
           <<"oooooooooooooooooooooooooooooooooooooooo"
           <<G4endl;
    G4cout.precision(6);
    G4cout <<"IN G4WilsonAblationModel" <<G4endl;
    G4cout <<"Initial prefragment A=" <<A
           <<", Z=" <<Z
           <<", excitation energy = " <<ex/MeV <<" MeV"
           <<G4endl; 
  }
//
//
// Check that there is a nucleus to speak of.  It's possible there isn't one
// or its just a proton or neutron.  In either case, the excitation energy
// (from the Lorentz vector) is not used.
//
  if (A == 0)
  {
    if (verboseLevel >= 2)
    {
      G4cout <<"No nucleus to decay" <<G4endl;
      G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
             <<"oooooooooooooooooooooooooooooooooooooooo"
             <<G4endl;
    }
    return fragmentVector;
  }
  else if (A == 1)
  {
    G4LorentzVector lorentzVector = theNucleus.GetMomentum();
    lorentzVector.setE(lorentzVector.e()-ex+10.0*eV);
    if (Z == 0)
    {
      G4Fragment *fragment = new G4Fragment(lorentzVector,G4Neutron::Neutron());
      fragmentVector->push_back(fragment);
    }
    else
    {
      G4Fragment *fragment = new G4Fragment(lorentzVector,G4Proton::Proton());
      fragmentVector->push_back(fragment);
    }
    if (verboseLevel >= 2)
    {
      G4cout <<"Final fragment is in fact only a nucleon) :" <<G4endl;
      G4cout <<(*fragmentVector)[0] <<G4endl;
      G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
             <<"oooooooooooooooooooooooooooooooooooooooo"
             <<G4endl;
    }
    return fragmentVector;
  }
//
//
// Then the number of nucleons ablated (either as nucleons or light nuclear
// fragments) is based on a simple argument for the binding energy per nucleon.
//
  G4int DAabl = (G4int) (ex / B);
  if (DAabl > A) DAabl = A;
// The following lines are no longer accurate given we now treat the final fragment
//  if (verboseLevel >= 2)
//    G4cout <<"Number of nucleons ejected = " <<DAabl <<G4endl;

//
//
// Determine the nuclear fragment from the ablation process by sampling the
// Rudstam equation.
//
  G4int AF = A - DAabl;
  G4int ZF = 0;
  
  if (AF > 0)
  {
    G4Pow* g4pow = G4Pow::GetInstance(); 
    G4double AFd = (G4double) AF;
    G4double R = 11.8 / g4pow->powZ(AF, 0.45);
    G4int minZ = std::max(1, Z - DAabl);
//
//
// Here we define an integral probability distribution based on the Rudstam
// equation assuming a constant AF.
//    
    G4int zmax = std::min(199, Z);
    G4double sum = 0.0;
    for (ZF=minZ; ZF<=zmax; ++ZF)
    {
      sum += G4Exp(-R*g4pow->powA(std::abs(ZF - 0.486*AFd + 3.8E-04*AFd*AFd),1.5));
      fSig[ZF] = sum;
    }
//
//
// Now sample that distribution to determine a value for ZF.
//
    sum *= G4UniformRand();
    for (ZF=minZ; ZF<=zmax; ++ZF) {
      if(sum <= fSig[ZF]) { break; }
    } 
  }
  G4int DZabl = Z - ZF;
//
//
// Now determine the nucleons or nuclei which have bee ablated.  The preference
// is for the production of alphas, then other nuclei in order of decreasing
// binding energy. The energies assigned to the products of the decay are
// provisional for the moment (the 10eV is just to avoid errors with negative
// excitation energies due to rounding).
//
  G4double totalEpost = 0.0;
  evapType.clear();
  for (G4int ift=0; ift<nFragTypes; ift++)
  {
    G4ParticleDefinition *type = fragType[ift];
    G4double n  = std::floor((G4double) DAabl / type->GetBaryonNumber() + 1.0E-10);
    G4double n1 = 1.0E+10;
    if (fragType[ift]->GetPDGCharge() > 0.0)
      n1 = std::floor((G4double) DZabl / type->GetPDGCharge() + 1.0E-10);
    if (n > n1) n = n1;
    if (n > 0.0)
    {
      G4double mass = type->GetPDGMass();
      for (G4int j=0; j<(G4int) n; j++)
      {
        totalEpost += mass;
        evapType.push_back(type);
      }
      DAabl -= (G4int) (n * type->GetBaryonNumber() + 1.0E-10);
      DZabl -= (G4int) (n * type->GetPDGCharge() + 1.0E-10);
    }
  }
//
//
// Determine the properties of the final nuclear fragment.  Note that if
// the final fragment is predicted to have a nucleon number of zero, then
// really it's the particle last in the vector evapType which becomes the
// final fragment.  Therefore delete this from the vector if this is the
// case.
//
  G4double massFinalFrag = 0.0;
  if (AF > 0)
    massFinalFrag = G4ParticleTable::GetParticleTable()->GetIonTable()->
      GetIonMass(ZF,AF);
  else
  {
    G4ParticleDefinition *type = evapType[evapType.size()-1];
    AF                         = type->GetBaryonNumber();
    ZF                         = (G4int) (type->GetPDGCharge() + 1.0E-10);
    evapType.erase(evapType.end()-1);
  }
  totalEpost   += massFinalFrag;
//
//
// Provide verbose output on the nuclear fragment if requested.
//
  if (verboseLevel >= 2)
  {
    G4cout <<"Final fragment      A=" <<AF
           <<", Z=" <<ZF
           <<G4endl;
    for (G4int ift=0; ift<nFragTypes; ift++)
    {
      G4ParticleDefinition *type = fragType[ift];
      G4int n                    = std::count(evapType.begin(),evapType.end(),type);
      if (n > 0) 
        G4cout <<"Particle type: " <<std::setw(10) <<type->GetParticleName()
               <<", number of particles emitted = " <<n <<G4endl;
    }
  }
//
// Add the total energy from the fragment.  Note that the fragment is assumed
// to be de-excited and does not undergo photo-evaporation .... I did mention
// this is a bit of a crude model?
//
  G4double massPreFrag      = theNucleus.GetGroundStateMass();
  G4double totalEpre        = massPreFrag + ex;
  G4double excess           = totalEpre - totalEpost;
//  G4Fragment *resultNucleus(theNucleus);
  G4Fragment *resultNucleus = new G4Fragment(A, Z, theNucleus.GetMomentum());
  G4ThreeVector boost(0.0,0.0,0.0);
  G4int nEvap = 0;
  if (produceSecondaries && evapType.size()>0)
  {
    if (excess > 0.0)
    {
      SelectSecondariesByEvaporation (resultNucleus);
      nEvap = fragmentVector->size();
      boost = resultNucleus->GetMomentum().findBoostToCM();
      if (evapType.size() > 0)
        SelectSecondariesByDefault (boost);
    }
    else
      SelectSecondariesByDefault(G4ThreeVector(0.0,0.0,0.0));
  }

  if (AF > 0)
  {
    G4double mass = G4ParticleTable::GetParticleTable()->GetIonTable()->
      GetIonMass(ZF,AF);
    G4double e    = mass + 10.0*eV;
    G4double p    = std::sqrt(e*e-mass*mass);
    G4ThreeVector direction(0.0,0.0,1.0);
    G4LorentzVector lorentzVector = G4LorentzVector(direction*p, e);
    lorentzVector.boost(-boost);
    G4Fragment* frag = new G4Fragment(AF, ZF, lorentzVector);
    fragmentVector->push_back(frag);
  }
  delete resultNucleus;
//
//
// Provide verbose output on the ablation products if requested.
//
  if (verboseLevel >= 2)
  {
    if (nEvap > 0)
    {
      G4cout <<"----------------------" <<G4endl;
      G4cout <<"Evaporated particles :" <<G4endl;
      G4cout <<"----------------------" <<G4endl;
    }
    G4int ie = 0;
    G4FragmentVector::iterator iter;
    for (iter = fragmentVector->begin(); iter != fragmentVector->end(); iter++)
    {
      if (ie == nEvap)
      {
//        G4cout <<*iter  <<G4endl;
        G4cout <<"---------------------------------" <<G4endl;
        G4cout <<"Particles from default emission :" <<G4endl;
        G4cout <<"---------------------------------" <<G4endl;
      }
      G4cout <<*iter <<G4endl;
    }
    G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
           <<"oooooooooooooooooooooooooooooooooooooooo"
           <<G4endl;
  }

  return fragmentVector;    
}

Here is the call graph for this function:

Here is the caller graph for this function:

GetProduceSecondaries()

G4bool G4WilsonAblationModel::GetProduceSecondaries ( )

inline

Definition at line 120 of file G4WilsonAblationModel.hh.

  {return produceSecondaries;}

GetVerboseLevel()

G4int G4WilsonAblationModel::GetVerboseLevel ( )

inline

Definition at line 128 of file G4WilsonAblationModel.hh.

  {return verboseLevel;}

PrintWelcomeMessage()

void G4WilsonAblationModel::PrintWelcomeMessage ( )

private

Definition at line 586 of file G4WilsonAblationModel.cc.

{
  G4cout <<G4endl;
  G4cout <<" *****************************************************************"
         <<G4endl;
  G4cout <<" Nuclear ablation model for nuclear-nuclear interactions activated"
         <<G4endl;
  G4cout <<" (Written by QinetiQ Ltd for the European Space Agency)"
         <<G4endl;
  G4cout <<" !!! WARNING: This model is not well validation and should not be used for accurate simulation !!!"
         <<G4endl;
  G4cout <<" *****************************************************************"
         <<G4endl;
  G4cout << G4endl;

  return;
}

Here is the caller graph for this function:

SelectSecondariesByDefault()

void G4WilsonAblationModel::SelectSecondariesByDefault ( G4ThreeVector  boost )

private

Definition at line 558 of file G4WilsonAblationModel.cc.

{
  for (unsigned i=0; i<evapType.size(); i++)
  {
    G4ParticleDefinition *type = evapType[i];
    G4double mass              = type->GetPDGMass();
    G4double e                 = mass + 10.0*eV;
    G4double p                 = std::sqrt(e*e-mass*mass);
    G4double costheta          = 2.0*G4UniformRand() - 1.0;
    G4double sintheta          = std::sqrt((1.0 - costheta)*(1.0 + costheta));
    G4double phi               = twopi * G4UniformRand() * rad;
    G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
    G4LorentzVector lorentzVector = G4LorentzVector(direction*p, e);
    lorentzVector.boost(-boost);
// Possibility that the following line is not correctly carrying over A and Z
// from particle definition.  Force values.  PRT 03/12/2009.
//    G4Fragment *fragment          = 
//      new G4Fragment(lorentzVector, type);
    G4int A = type->GetBaryonNumber();
    G4int Z = (G4int) (type->GetPDGCharge() + 1.0E-10);
    G4Fragment *fragment          = 
      new G4Fragment(A, Z, lorentzVector);

    fragmentVector->push_back(fragment);
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

SelectSecondariesByEvaporation()

void G4WilsonAblationModel::SelectSecondariesByEvaporation ( G4Fragment *  intermediateNucleus )

private

Definition at line 432 of file G4WilsonAblationModel.cc.

{
  G4Fragment theResidualNucleus = *intermediateNucleus;
  G4bool evaporate = true;
  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
  while (evaporate && evapType.size() != 0)
  {
//
//
// Here's the cheaky bit.  We're hijacking the G4Evaporation model, in order to
// more accurately sample to kinematics, but the species of the nuclear
// fragments will be the ones of our choosing as above.
//
    std::vector <G4VEvaporationChannel*>  theChannels1;
    theChannels1.clear();
    std::vector <G4VEvaporationChannel*>::iterator i;
    VectorOfFragmentTypes::iterator iter;
    std::vector <VectorOfFragmentTypes::iterator> iters;
    iters.clear();
    iter = std::find(evapType.begin(), evapType.end(), G4Alpha::Alpha());
    if (iter != evapType.end())
    {
      theChannels1.push_back(new G4AlphaEvaporationChannel);
      i = theChannels1.end() - 1;
      (*i)->SetOPTxs(OPTxs);
      (*i)->UseSICB(useSICB);
//      (*i)->Initialize(theResidualNucleus);
      iters.push_back(iter);
    }
    iter = std::find(evapType.begin(), evapType.end(), G4He3::He3());
    if (iter != evapType.end())
    {
      theChannels1.push_back(new G4He3EvaporationChannel);
      i = theChannels1.end() - 1;
      (*i)->SetOPTxs(OPTxs);
      (*i)->UseSICB(useSICB);
//      (*i)->Initialize(theResidualNucleus);
      iters.push_back(iter);
    }
    iter = std::find(evapType.begin(), evapType.end(), G4Triton::Triton());
    if (iter != evapType.end())
    {
      theChannels1.push_back(new G4TritonEvaporationChannel);
      i = theChannels1.end() - 1;
      (*i)->SetOPTxs(OPTxs);
      (*i)->UseSICB(useSICB);
//      (*i)->Initialize(theResidualNucleus);
      iters.push_back(iter);
    }
    iter = std::find(evapType.begin(), evapType.end(), G4Deuteron::Deuteron());
    if (iter != evapType.end())
    {
      theChannels1.push_back(new G4DeuteronEvaporationChannel);
      i = theChannels1.end() - 1;
      (*i)->SetOPTxs(OPTxs);
      (*i)->UseSICB(useSICB);
//      (*i)->Initialize(theResidualNucleus);
      iters.push_back(iter);
    }
    iter = std::find(evapType.begin(), evapType.end(), G4Proton::Proton());
    if (iter != evapType.end())
    {
      theChannels1.push_back(new G4ProtonEvaporationChannel);
      i = theChannels1.end() - 1;
      (*i)->SetOPTxs(OPTxs);
      (*i)->UseSICB(useSICB);
//      (*i)->Initialize(theResidualNucleus);
      iters.push_back(iter);
    }
    iter = std::find(evapType.begin(), evapType.end(), G4Neutron::Neutron());
    if (iter != evapType.end())
    {
      theChannels1.push_back(new G4NeutronEvaporationChannel);
      i = theChannels1.end() - 1;
      (*i)->SetOPTxs(OPTxs);
      (*i)->UseSICB(useSICB);
//      (*i)->Initialize(theResidualNucleus);
      iters.push_back(iter);
    }
    G4int nChannels = theChannels1.size();

    G4double totalProb = 0.0;
    G4int ich = 0;
    G4double probEvapType[6] = {0.0};
    std::vector<G4VEvaporationChannel*>::iterator iterEv;
    for (iterEv=theChannels1.begin(); iterEv!=theChannels1.end(); iterEv++) {
      totalProb += (*iterEv)->GetEmissionProbability(intermediateNucleus);
      probEvapType[ich] = totalProb;
      ++ich;
    }
    if (totalProb > 0.0) {
//
//
// The emission probability for at least one of the evaporation channels is
// positive, therefore work out which one should be selected and decay
// the nucleus.
//
      G4double xi = totalProb*G4UniformRand();
      G4int ii     = 0;
      for (ii=0; ii<nChannels; ii++) {
        if (xi < probEvapType[ii]) { break; }
      }
      if (ii >= nChannels) { ii = nChannels - 1; }
      G4FragmentVector *evaporationResult = theChannels1[ii]->
        BreakUp(*intermediateNucleus);
      fragmentVector->push_back((*evaporationResult)[0]);
      *intermediateNucleus = *(*evaporationResult)[1];
      //delete evaporationResult->back();
      delete evaporationResult;
      //evapType.erase(iters[ii]);
    }
    else
    {
//
//
// Probability for further evaporation is nil so have to escape from this
// routine and set the energies of the secondaries to 10eV.
//
      evaporate = false;
    }
  }
  
  return;
}

Here is the call graph for this function:

Here is the caller graph for this function:

SetProduceSecondaries()

void G4WilsonAblationModel::SetProduceSecondaries ( G4bool  produceSecondaries1 )

inline

Definition at line 116 of file G4WilsonAblationModel.hh.

  {produceSecondaries = produceSecondaries1;}

SetVerboseLevel()

void G4WilsonAblationModel::SetVerboseLevel ( G4int  verboseLevel1 )

inline

Definition at line 124 of file G4WilsonAblationModel.hh.

  {verboseLevel = verboseLevel1;}

Here is the caller graph for this function:

Member Data Documentation

B

G4double G4WilsonAblationModel::B

private

Definition at line 106 of file G4WilsonAblationModel.hh.

evapType

VectorOfFragmentTypes G4WilsonAblationModel::evapType

private

Definition at line 110 of file G4WilsonAblationModel.hh.

fragmentVector

G4FragmentVector* G4WilsonAblationModel::fragmentVector

private

Definition at line 109 of file G4WilsonAblationModel.hh.

fragType

G4ParticleDefinition* G4WilsonAblationModel::fragType[6]

private

Definition at line 108 of file G4WilsonAblationModel.hh.

fSig

G4double G4WilsonAblationModel::fSig[200]

private

Definition at line 111 of file G4WilsonAblationModel.hh.

nFragTypes

G4int G4WilsonAblationModel::nFragTypes

private

Definition at line 107 of file G4WilsonAblationModel.hh.

produceSecondaries

G4bool G4WilsonAblationModel::produceSecondaries

private

Definition at line 104 of file G4WilsonAblationModel.hh.

verboseLevel

G4int G4WilsonAblationModel::verboseLevel

private

Definition at line 105 of file G4WilsonAblationModel.hh.

---

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

• G4WilsonAblationModel.hh
• G4WilsonAblationModel.cc