G4WilsonAbrasionModel Class Reference

#include <G4WilsonAbrasionModel.hh>

Inheritance diagram for G4WilsonAbrasionModel:

Collaboration diagram for G4WilsonAbrasionModel:

Public Member Functions
	G4WilsonAbrasionModel (G4bool useAblation1=false)
	G4WilsonAbrasionModel (G4ExcitationHandler *)
	~G4WilsonAbrasionModel ()
	G4WilsonAbrasionModel (const G4WilsonAbrasionModel &right)
const G4WilsonAbrasionModel &	operator= (G4WilsonAbrasionModel &right)
virtual G4HadFinalState *	ApplyYourself (const G4HadProjectile &, G4Nucleus &)
void	SetVerboseLevel (G4int)
void	SetUseAblation (G4bool)
G4bool	GetUseAblation ()
void	SetConserveMomentum (G4bool)
G4bool	GetConserveMomentum ()
void	SetExcitationHandler (G4ExcitationHandler *)
G4ExcitationHandler *	GetExcitationHandler ()
virtual void	ModelDescription (std::ostream &) const
Public Member Functions inherited from G4HadronicInteraction
	G4HadronicInteraction (const G4String &modelName="HadronicModel")
virtual	~G4HadronicInteraction ()
virtual G4double	SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
virtual G4bool	IsApplicable (const G4HadProjectile &, G4Nucleus &)
G4double	GetMinEnergy () const
G4double	GetMinEnergy (const G4Material *aMaterial, const G4Element *anElement) const
void	SetMinEnergy (G4double anEnergy)
void	SetMinEnergy (G4double anEnergy, const G4Element *anElement)
void	SetMinEnergy (G4double anEnergy, const G4Material *aMaterial)
G4double	GetMaxEnergy () const
G4double	GetMaxEnergy (const G4Material *aMaterial, const G4Element *anElement) const
void	SetMaxEnergy (const G4double anEnergy)
void	SetMaxEnergy (G4double anEnergy, const G4Element *anElement)
void	SetMaxEnergy (G4double anEnergy, const G4Material *aMaterial)
const G4HadronicInteraction *	GetMyPointer () const
virtual G4int	GetVerboseLevel () const
const G4String &	GetModelName () const
void	DeActivateFor (const G4Material *aMaterial)
void	ActivateFor (const G4Material *aMaterial)
void	DeActivateFor (const G4Element *anElement)
void	ActivateFor (const G4Element *anElement)
G4bool	IsBlocked (const G4Material *aMaterial) const
G4bool	IsBlocked (const G4Element *anElement) const
void	SetRecoilEnergyThreshold (G4double val)
G4double	GetRecoilEnergyThreshold () const
G4bool	operator== (const G4HadronicInteraction &right) const
G4bool	operator!= (const G4HadronicInteraction &right) const
virtual const std::pair< G4double, G4double >	GetFatalEnergyCheckLevels () const
virtual std::pair< G4double, G4double >	GetEnergyMomentumCheckLevels () const
void	SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)
virtual void	BuildPhysicsTable (const G4ParticleDefinition &)

Private Member Functions
void	PrintWelcomeMessage ()
G4Fragment *	GetAbradedNucleons (G4int, G4double, G4double, G4double)
G4double	GetNucleonInducedExcitation (G4double, G4double, G4double)
void	SetConserveEnergy (G4bool)
G4bool	GetConserveEnergy ()

Private Attributes
G4double	r0sq
G4double	npK
G4bool	useAblation
G4WilsonAblationModel *	theAblation
G4ExcitationHandler *	theExcitationHandler
G4ExcitationHandler *	theExcitationHandlerx
G4bool	conserveEnergy
G4bool	conserveMomentum
G4double	B
G4double	third
G4double	fradius

Additional Inherited Members
Protected Member Functions inherited from G4HadronicInteraction
void	SetModelName (const G4String &nam)
G4bool	IsBlocked () const
void	Block ()
Protected Attributes inherited from G4HadronicInteraction
G4HadFinalState	theParticleChange
G4int	verboseLevel
G4double	theMinEnergy
G4double	theMaxEnergy
G4bool	isBlocked

Detailed Description

Definition at line 77 of file G4WilsonAbrasionModel.hh.

Constructor & Destructor Documentation

G4WilsonAbrasionModel() [1/3]

G4WilsonAbrasionModel::G4WilsonAbrasionModel

(

G4bool

useAblation1 = false

)

Definition at line 117 of file G4WilsonAbrasionModel.cc.

 :G4HadronicInteraction("G4WilsonAbrasion")
{
  // 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;
  useAblation  = useAblation1;

  // No de-excitation handler has been supplied - define the default handler.

  theExcitationHandler  = new G4ExcitationHandler;
  theExcitationHandlerx = new G4ExcitationHandler;
  if (useAblation)
  {
    theAblation = new G4WilsonAblationModel;
    theAblation->SetVerboseLevel(verboseLevel);
    theExcitationHandler->SetEvaporation(theAblation);
    theExcitationHandlerx->SetEvaporation(theAblation);
  }
  else
  {
    theAblation                      = NULL;
    G4Evaporation * theEvaporation   = new G4Evaporation;
    G4FermiBreakUp * theFermiBreakUp = new G4FermiBreakUp;
    G4StatMF * theMF                 = new G4StatMF;
    theExcitationHandler->SetEvaporation(theEvaporation);
    theExcitationHandler->SetFermiModel(theFermiBreakUp);
    theExcitationHandler->SetMultiFragmentation(theMF);
    theExcitationHandler->SetMaxAandZForFermiBreakUp(12, 6);
    theExcitationHandler->SetMinEForMultiFrag(5.0*MeV);

    theEvaporation  = new G4Evaporation;
    theFermiBreakUp = new G4FermiBreakUp;
    theExcitationHandlerx->SetEvaporation(theEvaporation);
    theExcitationHandlerx->SetFermiModel(theFermiBreakUp);
    theExcitationHandlerx->SetMaxAandZForFermiBreakUp(12, 6);
  }

  // Set the minimum and maximum range for the model (despite nomanclature,
  // this is in energy per nucleon number).  

  SetMinEnergy(70.0*MeV);
  SetMaxEnergy(10.1*GeV);
  isBlocked = false;

  // npK, when mutiplied by the nuclear Fermi momentum, determines the range of
  // momentum over which the secondary nucleon momentum is sampled.

  r0sq = 0.0;
  npK = 5.0;
  B = 10.0 * MeV;
  third = 1.0 / 3.0;
  fradius = 0.99;
  conserveEnergy = false;
  conserveMomentum = true;
}

Here is the call graph for this function:

G4WilsonAbrasionModel() [2/3]

G4WilsonAbrasionModel::G4WilsonAbrasionModel

(

G4ExcitationHandler *

aExcitationHandler

)

Definition at line 188 of file G4WilsonAbrasionModel.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;

  theAblation = NULL;   //A.R. 26-Jul-2012 Coverity fix.
  useAblation = false;  //A.R. 14-Aug-2012 Coverity fix.
                      
//
// The user is able to provide the excitation handler as well as an argument
// which is provided in this instantiation is used to determine
// whether the spectators of the interaction are free following the abrasion.
//
  theExcitationHandler             = aExcitationHandler;
  theExcitationHandlerx            = new G4ExcitationHandler;
  G4Evaporation * theEvaporation   = new G4Evaporation;
  G4FermiBreakUp * theFermiBreakUp = new G4FermiBreakUp;
  theExcitationHandlerx->SetEvaporation(theEvaporation);
  theExcitationHandlerx->SetFermiModel(theFermiBreakUp);
  theExcitationHandlerx->SetMaxAandZForFermiBreakUp(12, 6);
//
//
// Set the minimum and maximum range for the model (despite nomanclature, this
// is in energy per nucleon number).  
//
  SetMinEnergy(70.0*MeV);
  SetMaxEnergy(10.1*GeV);
  isBlocked = false;
//
//
// npK, when mutiplied by the nuclear Fermi momentum, determines the range of
// momentum over which the secondary nucleon momentum is sampled.
//
  r0sq             = 0.0;  //A.R. 14-Aug-2012 Coverity fix. 
  npK              = 5.0;
  B                = 10.0 * MeV;
  third            = 1.0 / 3.0;
  fradius          = 0.99;
  conserveEnergy   = false;
  conserveMomentum = true;
}

Here is the call graph for this function:

~G4WilsonAbrasionModel()

G4WilsonAbrasionModel::~G4WilsonAbrasionModel

(

)

Definition at line 236 of file G4WilsonAbrasionModel.cc.

{
//
//
// The destructor doesn't have to do a great deal!
//
  if (theExcitationHandler)  delete theExcitationHandler;
  if (theExcitationHandlerx) delete theExcitationHandlerx;
  if (useAblation)           delete theAblation;
//  delete theExcitationHandler;
//  delete theExcitationHandlerx;
}

G4WilsonAbrasionModel() [3/3]

G4WilsonAbrasionModel::G4WilsonAbrasionModel

(

const G4WilsonAbrasionModel &

right

)

Member Function Documentation

ApplyYourself()

G4HadFinalState * G4WilsonAbrasionModel::ApplyYourself ( const G4HadProjectile &  theTrack, G4Nucleus &  theTarget  )

virtual

Implements G4HadronicInteraction.

Definition at line 250 of file G4WilsonAbrasionModel.cc.

{
//
//
// The secondaries will be returned in G4HadFinalState &theParticleChange -
// initialise this.  The original track will always be discontinued and
// secondaries followed.
//
  theParticleChange.Clear();
  theParticleChange.SetStatusChange(stopAndKill);
//
//
// Get relevant information about the projectile and target (A, Z, energy/nuc,
// momentum, etc).
//
  const G4ParticleDefinition *definitionP = theTrack.GetDefinition();
  const G4double AP  = definitionP->GetBaryonNumber();
  const G4double ZP  = definitionP->GetPDGCharge();
  G4LorentzVector pP = theTrack.Get4Momentum();
  G4double E         = theTrack.GetKineticEnergy()/AP;
  G4double AT        = theTarget.GetA_asInt();
  G4double ZT        = theTarget.GetZ_asInt();
  G4double TotalEPre = theTrack.GetTotalEnergy() +
    theTarget.AtomicMass(AT, ZT) + theTarget.GetEnergyDeposit();
  G4double TotalEPost = 0.0;
//
//
// Determine the radii of the projectile and target nuclei.
//
  G4WilsonRadius aR;
  G4double rP   = aR.GetWilsonRadius(AP);
  G4double rT   = aR.GetWilsonRadius(AT);
  G4double rPsq = rP * rP;
  G4double rTsq = rT * rT;
  if (verboseLevel >= 2)
  {
    G4cout <<"########################################"
           <<"########################################"
           <<G4endl;
    G4cout.precision(6);
    G4cout <<"IN G4WilsonAbrasionModel" <<G4endl;
    G4cout <<"Initial projectile A=" <<AP 
           <<", Z=" <<ZP
           <<", radius = " <<rP/fermi <<" fm"
           <<G4endl; 
    G4cout <<"Initial target     A=" <<AT
           <<", Z=" <<ZT
           <<", radius = " <<rT/fermi <<" fm"
           <<G4endl;
    G4cout <<"Projectile momentum and Energy/nuc = " <<pP <<" ," <<E <<G4endl;
  }
//
//
// The following variables are used to determine the impact parameter in the
// near-field (i.e. taking into consideration the electrostatic repulsion).
//
  G4double rm   = ZP * ZT * elm_coupling / (E * AP);
  G4double r    = 0.0;
  G4double rsq  = 0.0;
//
//
// Initialise some of the variables which wll be used to calculate the chord-
// length for nucleons in the projectile and target, and hence calculate the
// number of abraded nucleons and the excitation energy.
//
  G4NuclearAbrasionGeometry *theAbrasionGeometry = NULL;
  G4double CT   = 0.0;
  G4double F    = 0.0;
  G4int Dabr    = 0;
//
//
// The following loop is performed until the number of nucleons which are
// abraded by the process is >1, i.e. an interaction MUST occur.
//
  G4bool skipInteraction = false;  // It will be set true if the two nuclei fail to collide 
  const G4int maxNumberOfLoops = 1000;
  G4int loopCounter = -1;
  while (Dabr == 0 && ++loopCounter < maxNumberOfLoops)  /* Loop checking, 07.08.2015, A.Ribon */
  {
// Added by MHM 20050119 to fix leaking memory on second pass through this loop
    if (theAbrasionGeometry)
    {
      delete theAbrasionGeometry;
      theAbrasionGeometry = NULL;
    }
//
//
// Sample the impact parameter.  For the moment, this class takes account of
// electrostatic effects on the impact parameter, but (like HZETRN AND NUCFRG2)
// does not make any correction for the effects of nuclear-nuclear repulsion.
//
    G4double rPT   = rP + rT;
    G4double rPTsq = rPT * rPT;
//
//
// This is a "catch" to make sure we don't go into an infinite loop because the
// energy is too low to overcome nuclear repulsion.  PRT 20091023.  If the
// value of rm < fradius * rPT then we're unlikely to sample a small enough
// impact parameter (energy of incident particle is too low).
//
    if (rm >= fradius * rPT) {
      skipInteraction = true;
    }
//
//
// Now sample impact parameter until the criterion is met that projectile
// and target overlap, but repulsion is taken into consideration.
//
    G4int evtcnt   = 0;
    r              = 1.1 * rPT;
    while (r > rPT && ++evtcnt < 1000)  /* Loop checking, 07.08.2015, A.Ribon */
    {
      G4double bsq = rPTsq * G4UniformRand();
      r            = (rm + std::sqrt(rm*rm + 4.0*bsq)) / 2.0;
    }
//
//
// We've tried to sample this 1000 times, but failed.  
//
    if (evtcnt >= 1000) {
      skipInteraction = true;
    }

    rsq = r * r;
//
//
// Now determine the chord-length through the target nucleus.
//
    if (rT > rP)
    {
      G4double x = (rPsq + rsq - rTsq) / 2.0 / r;
      if (x > 0.0) CT = 2.0 * std::sqrt(rTsq - x*x);
      else         CT = 2.0 * std::sqrt(rTsq - rsq);
    }
    else
    {
      G4double x = (rTsq + rsq - rPsq) / 2.0 / r;
      if (x > 0.0) CT = 2.0 * std::sqrt(rTsq - x*x);
      else         CT = 2.0 * rT;
    }
//
//
// Determine the number of abraded nucleons.  Note that the mean number of
// abraded nucleons is used to sample the Poisson distribution.  The Poisson
// distribution is sampled only ten times with the current impact parameter,
// and if it fails after this to find a case for which the number of abraded
// nucleons >1, the impact parameter is re-sampled.
//
    theAbrasionGeometry = new G4NuclearAbrasionGeometry(AP,AT,r);
    F                   = theAbrasionGeometry->F();
    G4double lambda     = 16.6*fermi / G4Pow::GetInstance()->powA(E/MeV,0.26);
    G4double Mabr       = F * AP * (1.0 - G4Exp(-CT/lambda));
    G4long n            = 0;
    for (G4int i = 0; i<10; i++)
    {
      n = G4Poisson(Mabr);
      if (n > 0)
      {
        if (n>AP) Dabr = (G4int) AP;
        else      Dabr = (G4int) n;
        break;
      }
    }
  }  // End of while loop

  if ( loopCounter >= maxNumberOfLoops || skipInteraction ) {
    // Assume nuclei do not collide and return unchanged primary. 
    theParticleChange.SetStatusChange(isAlive);
    theParticleChange.SetEnergyChange(theTrack.GetKineticEnergy());
    theParticleChange.SetMomentumChange(theTrack.Get4Momentum().vect().unit());
    if (verboseLevel >= 2) {
      G4cout <<"Particle energy too low to overcome repulsion." <<G4endl;
      G4cout <<"Event rejected and original track maintained" <<G4endl;
      G4cout <<"########################################"
             <<"########################################"
             <<G4endl;
    }
    return &theParticleChange;
  }

  if (verboseLevel >= 2)
  {
    G4cout <<G4endl;
    G4cout <<"Impact parameter    = " <<r/fermi <<" fm" <<G4endl;
    G4cout <<"# Abraded nucleons  = " <<Dabr <<G4endl;
  }
//
//
// The number of abraded nucleons must be no greater than the number of
// nucleons in either the projectile or the target.  If AP - Dabr < 2 or 
// AT - Dabr < 2 then either we have only a nucleon left behind in the
// projectile/target or we've tried to abrade too many nucleons - and Dabr
// should be limited.
//
  if (AP - (G4double) Dabr < 2.0) Dabr = (G4int) AP;
  if (AT - (G4double) Dabr < 2.0) Dabr = (G4int) AT;
//
//
// Determine the abraded secondary nucleons from the projectile.  *fragmentP
// is a pointer to the prefragment from the projectile and nSecP is the number
// of nucleons in theParticleChange which have been abraded.  The total energy
// from these is determined.
//
  G4ThreeVector boost   = pP.findBoostToCM();
  G4Fragment *fragmentP = GetAbradedNucleons (Dabr, AP, ZP, rP); 
  G4int nSecP           = theParticleChange.GetNumberOfSecondaries();
  G4int i               = 0;
  for (i=0; i<nSecP; i++)
  {
    TotalEPost += theParticleChange.GetSecondary(i)->
      GetParticle()->GetTotalEnergy();
  }
//
//
// Determine the number of spectators in the interaction region for the
// projectile.
//
  G4int DspcP = (G4int) (AP*F) - Dabr;
  if (DspcP <= 0)           DspcP = 0;
  else if (DspcP > AP-Dabr) DspcP = ((G4int) AP) - Dabr;
//
//
// Determine excitation energy associated with excess surface area of the
// projectile (EsP) and the excitation due to scattering of nucleons which are
// retained within the projectile (ExP).  Add the total energy from the excited
// nucleus to the total energy of the secondaries.
//
  G4bool excitationAbsorbedByProjectile = false;
  if (fragmentP != NULL)
  {
    G4double EsP = theAbrasionGeometry->GetExcitationEnergyOfProjectile();
    G4double ExP = 0.0;
    if (Dabr < AT)
      excitationAbsorbedByProjectile = G4UniformRand() < 0.5;
    if (excitationAbsorbedByProjectile)
      ExP = GetNucleonInducedExcitation(rP, rT, r);
    G4double xP = EsP + ExP;
    if (xP > B*(AP-Dabr)) xP = B*(AP-Dabr);
    G4LorentzVector lorentzVector = fragmentP->GetMomentum();
    lorentzVector.setE(lorentzVector.e()+xP);
    fragmentP->SetMomentum(lorentzVector);
    TotalEPost += lorentzVector.e();
  }
  G4double EMassP = TotalEPost;
//
//
// Determine the abraded secondary nucleons from the target.  Note that it's
// assumed that the same number of nucleons are abraded from the target as for
// the projectile, and obviously no boost is applied to the products. *fragmentT
// is a pointer to the prefragment from the target and nSec is the total number
// of nucleons in theParticleChange which have been abraded.  The total energy
// from these is determined.
//
  G4Fragment *fragmentT = GetAbradedNucleons (Dabr, AT, ZT, rT); 
  G4int nSec = theParticleChange.GetNumberOfSecondaries();
  for (i=nSecP; i<nSec; i++)
  {
    TotalEPost += theParticleChange.GetSecondary(i)->
      GetParticle()->GetTotalEnergy();
  }
//
//
// Determine the number of spectators in the interaction region for the
// target.
//
  G4int DspcT = (G4int) (AT*F) - Dabr;
  if (DspcT <= 0)           DspcT = 0;
  else if (DspcT > AP-Dabr) DspcT = ((G4int) AT) - Dabr;
//
//
// Determine excitation energy associated with excess surface area of the
// target (EsT) and the excitation due to scattering of nucleons which are
// retained within the target (ExT).  Add the total energy from the excited
// nucleus to the total energy of the secondaries.
//
  if (fragmentT != NULL)
  {
    G4double EsT = theAbrasionGeometry->GetExcitationEnergyOfTarget();
    G4double ExT = 0.0;
    if (!excitationAbsorbedByProjectile)
      ExT = GetNucleonInducedExcitation(rT, rP, r);
    G4double xT = EsT + ExT;
    if (xT > B*(AT-Dabr)) xT = B*(AT-Dabr);
    G4LorentzVector lorentzVector = fragmentT->GetMomentum();
    lorentzVector.setE(lorentzVector.e()+xT);
    fragmentT->SetMomentum(lorentzVector);
    TotalEPost += lorentzVector.e();
  }
//
//
// Now determine the difference between the pre and post interaction
// energy - this will be used to determine the Lorentz boost if conservation
// of energy is to be imposed/attempted.
//
  G4double deltaE = TotalEPre - TotalEPost;
  if (deltaE > 0.0 && conserveEnergy)
  {
    G4double beta = std::sqrt(1.0 - EMassP*EMassP/G4Pow::GetInstance()->powN(deltaE+EMassP,2));
    boost = boost / boost.mag() * beta;
  }
//
//
// Now boost the secondaries from the projectile.
//
  G4ThreeVector pBalance = pP.vect();
  for (i=0; i<nSecP; i++)
  {
    G4DynamicParticle *dynamicP = theParticleChange.GetSecondary(i)->
      GetParticle();
    G4LorentzVector lorentzVector = dynamicP->Get4Momentum();
    lorentzVector.boost(-boost);
    dynamicP->Set4Momentum(lorentzVector);
    pBalance -= lorentzVector.vect();
  }
//
//
// Set the boost for the projectile prefragment.  This is now based on the
// conservation of momentum.  However, if the user selected momentum of the
// prefragment is not to be conserved this simply boosted to the velocity of the
// original projectile times the ratio of the unexcited to the excited mass
// of the prefragment (the excitation increases the effective mass of the
// prefragment, and therefore modifying the boost is an attempt to prevent
// the momentum of the prefragment being excessive).
//
  if (fragmentP != NULL)
  {
    G4LorentzVector lorentzVector = fragmentP->GetMomentum();
    G4double fragmentM            = lorentzVector.m();
    if (conserveMomentum)
      fragmentP->SetMomentum
        (G4LorentzVector(pBalance,std::sqrt(pBalance.mag2()+fragmentM*fragmentM+1.0*eV*eV)));
    else
    {
      G4double fragmentGroundStateM = fragmentP->GetGroundStateMass();
      fragmentP->SetMomentum(lorentzVector.boost(-boost * fragmentGroundStateM/fragmentM));
    }
  }
//
//
// Output information to user if verbose information requested.
//
  if (verboseLevel >= 2)
  {
    G4cout <<G4endl;
    G4cout <<"-----------------------------------" <<G4endl;
    G4cout <<"Secondary nucleons from projectile:" <<G4endl;
    G4cout <<"-----------------------------------" <<G4endl;
    G4cout.precision(7);
    for (i=0; i<nSecP; i++)
    {
      G4cout <<"Particle # " <<i <<G4endl;
      theParticleChange.GetSecondary(i)->GetParticle()->DumpInfo();
      G4DynamicParticle *dyn = theParticleChange.GetSecondary(i)->GetParticle();
      G4cout <<"New nucleon (P) " <<dyn->GetDefinition()->GetParticleName()
             <<" : "              <<dyn->Get4Momentum()
             <<G4endl;
    }
    G4cout <<"---------------------------" <<G4endl;
    G4cout <<"The projectile prefragment:" <<G4endl;
    G4cout <<"---------------------------" <<G4endl;
    if (fragmentP != NULL)
      G4cout <<*fragmentP <<G4endl;
    else
      G4cout <<"(No residual prefragment)" <<G4endl;
    G4cout <<G4endl;
    G4cout <<"-------------------------------" <<G4endl;
    G4cout <<"Secondary nucleons from target:" <<G4endl;
    G4cout <<"-------------------------------" <<G4endl;
    G4cout.precision(7);
    for (i=nSecP; i<nSec; i++)
    {
      G4cout <<"Particle # " <<i <<G4endl;
      theParticleChange.GetSecondary(i)->GetParticle()->DumpInfo();
      G4DynamicParticle *dyn = theParticleChange.GetSecondary(i)->GetParticle();
      G4cout <<"New nucleon (T) " <<dyn->GetDefinition()->GetParticleName()
             <<" : "              <<dyn->Get4Momentum()
             <<G4endl;
    }
    G4cout <<"-----------------------" <<G4endl;
    G4cout <<"The target prefragment:" <<G4endl;
    G4cout <<"-----------------------" <<G4endl;
    if (fragmentT != NULL)
      G4cout <<*fragmentT <<G4endl;
    else
      G4cout <<"(No residual prefragment)" <<G4endl;
  }
//
//
// Now we can decay the nuclear fragments if present.  The secondaries are
// collected and boosted as well.  This is performed first for the projectile...
//
  if (fragmentP !=NULL)
  {
    G4ReactionProductVector *products = NULL;
    if (fragmentP->GetZ_asInt() != fragmentP->GetA_asInt())
      products = theExcitationHandler->BreakItUp(*fragmentP);
    else
      products = theExcitationHandlerx->BreakItUp(*fragmentP);      
    delete fragmentP;
    fragmentP = NULL;
  
    G4ReactionProductVector::iterator iter;
    for (iter = products->begin(); iter != products->end(); ++iter)
    {
      G4DynamicParticle *secondary =
        new G4DynamicParticle((*iter)->GetDefinition(),
        (*iter)->GetTotalEnergy(), (*iter)->GetMomentum());
      theParticleChange.AddSecondary (secondary);  // Added MHM 20050118
      G4String particleName = (*iter)->GetDefinition()->GetParticleName();
      delete (*iter); // get rid of leftover particle def!  // Added MHM 20050118
      if (verboseLevel >= 2 && particleName.find("[",0) < particleName.size())
      {
        G4cout <<"------------------------" <<G4endl;
        G4cout <<"The projectile fragment:" <<G4endl;
        G4cout <<"------------------------" <<G4endl;
        G4cout <<" fragmentP = " <<particleName
               <<" Energy    = " <<secondary->GetKineticEnergy()
               <<G4endl;
      }
    }
    delete products;  // Added MHM 20050118
  }
//
//
// Now decay the target nucleus - no boost is applied since in this
// approximation it is assumed that there is negligible momentum transfer from
// the projectile.
//
  if (fragmentT != NULL)
  {
    G4ReactionProductVector *products = NULL;
    if (fragmentT->GetZ_asInt() != fragmentT->GetA_asInt())
      products = theExcitationHandler->BreakItUp(*fragmentT);
    else
      products = theExcitationHandlerx->BreakItUp(*fragmentT);      
    delete fragmentT;
    fragmentT = NULL;
  
    G4ReactionProductVector::iterator iter;
    for (iter = products->begin(); iter != products->end(); ++iter)
    {
      G4DynamicParticle *secondary =
        new G4DynamicParticle((*iter)->GetDefinition(),
        (*iter)->GetTotalEnergy(), (*iter)->GetMomentum());
      theParticleChange.AddSecondary (secondary);
      G4String particleName = (*iter)->GetDefinition()->GetParticleName();
      delete (*iter); // get rid of leftover particle def!  // Added MHM 20050118
      if (verboseLevel >= 2 && particleName.find("[",0) < particleName.size())
      {
        G4cout <<"--------------------" <<G4endl;
        G4cout <<"The target fragment:" <<G4endl;
        G4cout <<"--------------------" <<G4endl;
        G4cout <<" fragmentT = " <<particleName
               <<" Energy    = " <<secondary->GetKineticEnergy()
               <<G4endl;
      }
    }
    delete products;  // Added MHM 20050118
  }

  if (verboseLevel >= 2)
     G4cout <<"########################################"
            <<"########################################"
            <<G4endl;
  
  delete theAbrasionGeometry;
  
  return &theParticleChange;
}

Here is the call graph for this function:

GetAbradedNucleons()

G4Fragment * G4WilsonAbrasionModel::GetAbradedNucleons ( G4int  Dabr, G4double  A, G4double  Z, G4double  r  )

private

Definition at line 722 of file G4WilsonAbrasionModel.cc.

{
//
//
// Initialise variables.  tau is the Fermi radius of the nucleus.  The variables
// p..., C... and gamma are used to help sample the secondary nucleon
// spectrum.
//
  
  G4double pK   = hbarc * G4Pow::GetInstance()->A13(9.0 * pi / 4.0 * A) / (1.29 * r);
  if (A <= 24.0) pK *= -0.229*G4Pow::GetInstance()->A13(A) + 1.62;
  G4double pKsq = pK * pK;
  G4double p1sq = 2.0/5.0 * pKsq;
  G4double p2sq = 6.0/5.0 * pKsq;
  G4double p3sq = 500.0 * 500.0;
  G4double C1   = 1.0;
  G4double C2   = 0.03;
  G4double C3   = 0.0002;
  G4double gamma = 90.0 * MeV;
  G4double maxn = C1 + C2 + C3;
//
//
// initialise the number of secondary nucleons abraded to zero, and initially set
// the type of nucleon abraded to proton ... just for now.
//  
  G4double Aabr                     = 0.0;
  G4double Zabr                     = 0.0; 
  G4ParticleDefinition *typeNucleon = G4Proton::ProtonDefinition();
  G4DynamicParticle *dynamicNucleon = NULL;
  G4ParticleMomentum pabr(0.0, 0.0, 0.0);
//
//
// Now go through each abraded nucleon and sample type, spectrum and angle.
//
  G4bool isForLoopExitAnticipated = false;
  for (G4int i=0; i<Dabr; i++)
  {
//
//
// Sample the nucleon momentum distribution by simple rejection techniques.  We
// reject values of p == 0.0 since this causes bad behaviour in the sinh term.
//
    G4double p   = 0.0;
    G4bool found = false;
    const G4int maxNumberOfLoops = 100000;
    G4int loopCounter = -1;
    while (!found && ++loopCounter < maxNumberOfLoops)  /* Loop checking, 07.08.2015, A.Ribon */
    {
      while (p <= 0.0) p = npK * pK * G4UniformRand();  /* Loop checking, 07.08.2015, A.Ribon */
      G4double psq = p * p;
      found = maxn * G4UniformRand() < C1*G4Exp(-psq/p1sq/2.0) +
        C2*G4Exp(-psq/p2sq/2.0) + C3*G4Exp(-psq/p3sq/2.0) + p/gamma/(0.5*(G4Exp(p/gamma)-G4Exp(-p/gamma)));
    }
    if ( loopCounter >= maxNumberOfLoops )
    {
      isForLoopExitAnticipated = true;
      break;
    }
//
//
// Determine the type of particle abraded.  Can only be proton or neutron,
// and the probability is determine to be proportional to the ratio as found
// in the nucleus at each stage.
//
    G4double prob = (Z-Zabr)/(A-Aabr);
    if (G4UniformRand()<prob)
    {
      Zabr++;
      typeNucleon = G4Proton::ProtonDefinition();
    }
    else
      typeNucleon = G4Neutron::NeutronDefinition();
    Aabr++;
//
//
// The angular distribution of the secondary nucleons is approximated to an
// isotropic distribution in the rest frame of the nucleus (this will be Lorentz
// boosted later.
//
    G4double costheta = 2.*G4UniformRand()-1.0;
    G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
    G4double phi      = 2.0*pi*G4UniformRand()*rad;
    G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
    G4double nucleonMass = typeNucleon->GetPDGMass();
    G4double E           = std::sqrt(p*p + nucleonMass*nucleonMass)-nucleonMass;
    dynamicNucleon = new G4DynamicParticle(typeNucleon,direction,E);
    theParticleChange.AddSecondary (dynamicNucleon);
    pabr += p*direction;
  }
//
//
// Next determine the details of the nuclear prefragment .. that is if there
// is one or more protons in the residue.  (Note that the 1 eV in the total
// energy is a safety factor to avoid any possibility of negative rest mass
// energy.)
//
  G4Fragment *fragment = NULL;
  if ( ! isForLoopExitAnticipated && Z-Zabr>=1.0 )
  {
    G4double ionMass = G4ParticleTable::GetParticleTable()->GetIonTable()->
                       GetIonMass(G4lrint(Z-Zabr),G4lrint(A-Aabr));
    G4double E       = std::sqrt(pabr.mag2() + ionMass*ionMass);
    G4LorentzVector lorentzVector = G4LorentzVector(-pabr, E + 1.0*eV);
    fragment =
      new G4Fragment((G4int) (A-Aabr), (G4int) (Z-Zabr), lorentzVector);
  }

  return fragment;
}

Here is the call graph for this function:

Here is the caller graph for this function:

GetConserveEnergy()

G4bool G4WilsonAbrasionModel::GetConserveEnergy ( )

inlineprivate

Definition at line 133 of file G4WilsonAbrasionModel.hh.

  {return conserveEnergy;}

GetConserveMomentum()

G4bool G4WilsonAbrasionModel::GetConserveMomentum ( )

inline

Definition at line 139 of file G4WilsonAbrasionModel.hh.

  {return conserveMomentum;}

GetExcitationHandler()

G4ExcitationHandler * G4WilsonAbrasionModel::GetExcitationHandler ( )

inline

Definition at line 124 of file G4WilsonAbrasionModel.hh.

  {return theExcitationHandler;}

GetNucleonInducedExcitation()

G4double G4WilsonAbrasionModel::GetNucleonInducedExcitation ( G4double  rP, G4double  rT, G4double  r  )

private

Definition at line 835 of file G4WilsonAbrasionModel.cc.

{
//
//
// Initialise variables.
//
  G4double Cl   = 0.0;
  G4double rPsq = rP * rP;
  G4double rTsq = rT * rT;
  G4double rsq  = r * r;
//
//
// Depending upon the impact parameter, a different form of the chord length is
// is used.
//  
  if (r > rT) Cl = 2.0*std::sqrt(rPsq + 2.0*r*rT - rsq - rTsq);
  else        Cl = 2.0*rP;
//
//
// The next lines have been changed to include a "catch" to make sure if the
// projectile and target are too close, Ct is set to twice rP or twice rT.
// Otherwise the standard Wilson algorithm should work fine.
// PRT 20091023.
//
  G4double Ct = 0.0;
  if      (rT > rP && rsq < rTsq - rPsq) Ct = 2.0 * rP;
  else if (rP > rT && rsq < rPsq - rTsq) Ct = 2.0 * rT;
  else {
    G4double bP = (rPsq+rsq-rTsq)/2.0/r;
    G4double x  = rPsq - bP*bP;
    if (x < 0.0) {
      G4cerr <<"########################################"
             <<"########################################"
             <<G4endl;
      G4cerr <<"ERROR IN G4WilsonAbrasionModel::GetNucleonInducedExcitation"
             <<G4endl;
      G4cerr <<"rPsq - bP*bP < 0.0 and cannot be square-rooted" <<G4endl;
      G4cerr <<"Set to zero instead" <<G4endl;
      G4cerr <<"########################################"
             <<"########################################"
             <<G4endl;
    }
    Ct = 2.0*std::sqrt(x);
  }
  
  G4double Ex = 13.0 * Cl / fermi;
  if (Ct > 1.5*fermi)
    Ex += 13.0 * Cl / fermi /3.0 * (Ct/fermi - 1.5);

  return Ex;
}

Here is the caller graph for this function:

GetUseAblation()

G4bool G4WilsonAbrasionModel::GetUseAblation ( )

inline

Definition at line 127 of file G4WilsonAbrasionModel.hh.

  {return useAblation;}

ModelDescription()

void G4WilsonAbrasionModel::ModelDescription ( std::ostream &  outFile ) const

virtual

Reimplemented from G4HadronicInteraction.

Definition at line 176 of file G4WilsonAbrasionModel.cc.

{
  outFile << "G4WilsonAbrasionModel is a macroscopic treatment of\n"
          << "nucleus-nucleus collisions using simple geometric arguments.\n"
          << "The smaller projectile nucleus gouges out a part of the larger\n"
          << "target nucleus, leaving a residual nucleus and a fireball\n"
          << "region where the projectile and target intersect.  The fireball"
          << "is then treated as a highly excited nuclear fragment.  This\n"
          << "model is based on the NUCFRG2 model and is valid for all\n"
          << "projectile energies between 70 MeV/n and 10.1 GeV/n. \n";
}

operator=()

const G4WilsonAbrasionModel& G4WilsonAbrasionModel::operator=

(

G4WilsonAbrasionModel &

right

)

PrintWelcomeMessage()

void G4WilsonAbrasionModel::PrintWelcomeMessage ( )

private

Definition at line 927 of file G4WilsonAbrasionModel.cc.

{
  G4cout <<G4endl;
  G4cout <<" *****************************************************************"
         <<G4endl;
  G4cout <<" Nuclear abrasion model for nuclear-nuclear interactions activated"
         <<G4endl;
  G4cout <<" (Written by QinetiQ Ltd for the European Space Agency)"
         <<G4endl;
  G4cout <<" *****************************************************************"
         <<G4endl;
  G4cout << G4endl;

  return;
}

Here is the caller graph for this function:

SetConserveEnergy()

void G4WilsonAbrasionModel::SetConserveEnergy ( G4bool  conserveEnergy1 )

inlineprivate

Definition at line 130 of file G4WilsonAbrasionModel.hh.

  {conserveEnergy = conserveEnergy1;}

SetConserveMomentum()

void G4WilsonAbrasionModel::SetConserveMomentum ( G4bool  conserveMomentum1 )

inline

Definition at line 136 of file G4WilsonAbrasionModel.hh.

  {conserveMomentum = conserveMomentum1;}

SetExcitationHandler()

void G4WilsonAbrasionModel::SetExcitationHandler ( G4ExcitationHandler *  aExcitationHandler )

inline

Definition at line 121 of file G4WilsonAbrasionModel.hh.

  {theExcitationHandler = aExcitationHandler;}

SetUseAblation()

void G4WilsonAbrasionModel::SetUseAblation

(

G4bool

useAblation1

)

Definition at line 888 of file G4WilsonAbrasionModel.cc.

{
  if (useAblation != useAblation1)
  {
    useAblation = useAblation1;
    delete theExcitationHandler;
    delete theExcitationHandlerx;
    theExcitationHandler  = new G4ExcitationHandler;
    theExcitationHandlerx = new G4ExcitationHandler;
    if (useAblation)
    {
      theAblation = new G4WilsonAblationModel;
      theAblation->SetVerboseLevel(verboseLevel);
      theExcitationHandler->SetEvaporation(theAblation);
      theExcitationHandlerx->SetEvaporation(theAblation);
    }
    else
    {
      theAblation                      = NULL;
      G4Evaporation * theEvaporation   = new G4Evaporation;
      G4FermiBreakUp * theFermiBreakUp = new G4FermiBreakUp;
      G4StatMF * theMF                 = new G4StatMF;
      theExcitationHandler->SetEvaporation(theEvaporation);
      theExcitationHandler->SetFermiModel(theFermiBreakUp);
      theExcitationHandler->SetMultiFragmentation(theMF);
      theExcitationHandler->SetMaxAandZForFermiBreakUp(12, 6);
      theExcitationHandler->SetMinEForMultiFrag(5.0*MeV);

      theEvaporation  = new G4Evaporation;
      theFermiBreakUp = new G4FermiBreakUp;
      theExcitationHandlerx->SetEvaporation(theEvaporation);
      theExcitationHandlerx->SetFermiModel(theFermiBreakUp);
      theExcitationHandlerx->SetMaxAandZForFermiBreakUp(12, 6);
    }
  }
  return; 
}

Here is the call graph for this function:

SetVerboseLevel()

void G4WilsonAbrasionModel::SetVerboseLevel ( G4int  verboseLevel1 )

inlinevirtual

Reimplemented from G4HadronicInteraction.

Definition at line 142 of file G4WilsonAbrasionModel.hh.

{
  verboseLevel = verboseLevel1;
  if (useAblation) theAblation->SetVerboseLevel(verboseLevel);
}

Here is the call graph for this function:

Member Data Documentation

B

G4double G4WilsonAbrasionModel::B

private

Definition at line 115 of file G4WilsonAbrasionModel.hh.

conserveEnergy

G4bool G4WilsonAbrasionModel::conserveEnergy

private

Definition at line 113 of file G4WilsonAbrasionModel.hh.

conserveMomentum

G4bool G4WilsonAbrasionModel::conserveMomentum

private

Definition at line 114 of file G4WilsonAbrasionModel.hh.

fradius

G4double G4WilsonAbrasionModel::fradius

private

Definition at line 117 of file G4WilsonAbrasionModel.hh.

npK

G4double G4WilsonAbrasionModel::npK

private

Definition at line 108 of file G4WilsonAbrasionModel.hh.

r0sq

G4double G4WilsonAbrasionModel::r0sq

private

Definition at line 107 of file G4WilsonAbrasionModel.hh.

theAblation

G4WilsonAblationModel* G4WilsonAbrasionModel::theAblation

private

Definition at line 110 of file G4WilsonAbrasionModel.hh.

theExcitationHandler

G4ExcitationHandler* G4WilsonAbrasionModel::theExcitationHandler

private

Definition at line 111 of file G4WilsonAbrasionModel.hh.

theExcitationHandlerx

G4ExcitationHandler* G4WilsonAbrasionModel::theExcitationHandlerx

private

Definition at line 112 of file G4WilsonAbrasionModel.hh.

third

G4double G4WilsonAbrasionModel::third

private

Definition at line 116 of file G4WilsonAbrasionModel.hh.

useAblation

G4bool G4WilsonAbrasionModel::useAblation

private

Definition at line 109 of file G4WilsonAbrasionModel.hh.

---

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

• G4WilsonAbrasionModel.hh
• G4WilsonAbrasionModel.cc