G4WilsonAblationModel Class Reference

#include <G4WilsonAblationModel.hh>

Inheritance diagram for G4WilsonAblationModel:

Public Types
typedef std::vector< G4ParticleDefinition * >	VectorOfFragmentTypes
Public Member Functions
	G4WilsonAblationModel ()
	~G4WilsonAblationModel ()
G4FragmentVector *	BreakItUp (const G4Fragment &theNucleus)
void	SetProduceSecondaries (G4bool)
G4bool	GetProduceSecondaries ()
void	SetVerboseLevel (G4int)
G4int	GetVerboseLevel ()

---

Detailed Description

Definition at line 83 of file G4WilsonAblationModel.hh.

---

Member Typedef Documentation

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

Definition at line 89 of file G4WilsonAblationModel.hh.

---

Constructor & Destructor Documentation

G4WilsonAblationModel::G4WilsonAblationModel

(

)

Definition at line 105 of file G4WilsonAblationModel.cc.

References G4Alpha::Alpha(), G4Deuteron::Deuteron(), G4VEvaporationFactory::GetChannel(), G4He3::He3(), G4Neutron::Neutron(), G4VEvaporation::OPTxs, G4Proton::Proton(), G4Triton::Triton(), and G4VEvaporation::useSICB.

{
//
//
// 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();
//
//
// Set verboseLevel default to no output.
//
  verboseLevel = 0;
  theChannelFactory = new G4EvaporationFactory();
  theChannels = theChannelFactory->GetChannel();
//
//
// Set defaults for evaporation classes.  These can be overridden by user
// "set" methods.
//
  OPTxs   = 3;
  useSICB = false;
  fragmentVector = 0;
}

G4WilsonAblationModel::~G4WilsonAblationModel

(

)

Definition at line 159 of file G4WilsonAblationModel.cc.

{
  if (theChannels != 0) theChannels = 0;
  if (theChannelFactory != 0) delete theChannelFactory;
}

---

Member Function Documentation

G4FragmentVector * G4WilsonAblationModel::BreakItUp

(

const G4Fragment &

theNucleus

)

[virtual]

Implements G4VEvaporation.

Definition at line 167 of file G4WilsonAblationModel.cc.

References G4cout, G4endl, G4UniformRand, G4Fragment::GetA_asInt(), G4ParticleDefinition::GetBaryonNumber(), G4Fragment::GetExcitationEnergy(), G4Fragment::GetGroundStateMass(), G4ParticleTable::GetIonTable(), G4Fragment::GetMomentum(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGCharge(), G4Fragment::GetZ_asInt(), CLHEP::detail::n, G4Neutron::Neutron(), and G4Proton::Proton().

{
//
//
// 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)
  {
    G4double AFd = (G4double) AF;
    G4double R = 11.8 / std::pow(AFd, 0.45);
    G4int minZ = Z - DAabl;
    if (minZ <= 0) minZ = 1;
//
//
// Here we define an integral probability distribution based on the Rudstam
// equation assuming a constant AF.
//    
    G4double sig[100];
    G4double sum = 0.0;
    for (G4int ii=minZ; ii<= Z; ii++)
    {
      sum   += std::exp(-R*std::pow(std::abs(ii - 0.486*AFd + 3.8E-04*AFd*AFd),1.5));
      sig[ii] = sum;
    }
//
//
// Now sample that distribution to determine a value for ZF.
//
    G4double xi  = G4UniformRand();
    G4int iz     = minZ;
    G4bool found = false;
    while (iz <= Z && !found)
    {
      found = (xi <= sig[iz]/sum);
      if (!found) iz++;
    }
    if (iz > Z)
      ZF = Z;
    else
      ZF = iz;
  }
  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;    
}

G4bool G4WilsonAblationModel::GetProduceSecondaries

(

)

[inline]

Definition at line 122 of file G4WilsonAblationModel.hh.

  {return produceSecondaries;}

G4int G4WilsonAblationModel::GetVerboseLevel

(

)

[inline]

Definition at line 130 of file G4WilsonAblationModel.hh.

  {return verboseLevel;}

void G4WilsonAblationModel::SetProduceSecondaries

(

G4bool

)

[inline]

Definition at line 118 of file G4WilsonAblationModel.hh.

  {produceSecondaries = produceSecondaries1;}

void G4WilsonAblationModel::SetVerboseLevel

(

G4int

)

[inline]

Definition at line 126 of file G4WilsonAblationModel.hh.

Referenced by G4WilsonAbrasionModel::G4WilsonAbrasionModel(), G4WilsonAbrasionModel::SetUseAblation(), and G4WilsonAbrasionModel::SetVerboseLevel().

  {verboseLevel = verboseLevel1;}

---

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

• G4WilsonAblationModel.hh
• G4WilsonAblationModel.cc

---