G4RPGInelastic Class Reference

#include <G4RPGInelastic.hh>

Inheritance diagram for G4RPGInelastic:

Public Member Functions
	G4RPGInelastic (const G4String &modelName="RPGInelastic")
virtual	~G4RPGInelastic ()
Public Member Functions inherited from G4HadronicInteraction
	G4HadronicInteraction (const G4String &modelName="HadronicModel")
virtual	~G4HadronicInteraction ()
virtual G4HadFinalState *	ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)=0
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
virtual void	SetVerboseLevel (G4int value)
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	ModelDescription (std::ostream &outFile) const

Protected Types
enum	{   pi0, pip, pim, kp,   km, k0, k0b, pro,   neu, lam, sp, s0,   sm, xi0, xim, om,   ap, an }

Protected Member Functions
G4double	Pmltpc (G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)
G4int	Factorial (G4int n)
G4bool	MarkLeadingStrangeParticle (const G4ReactionProduct &currentParticle, const G4ReactionProduct &targetParticle, G4ReactionProduct &leadParticle)
void	SetUpPions (const G4int np, const G4int nm, const G4int nz, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
void	GetNormalizationConstant (const G4double availableEnergy, G4double &n, G4double &anpn)
void	CalculateMomenta (G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen, const G4HadProjectile *originalIncident, const G4DynamicParticle *originalTarget, G4ReactionProduct &modifiedOriginal, G4Nucleus &targetNucleus, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4bool &incidentHasChanged, G4bool &targetHasChanged, G4bool quasiElastic)
void	SetUpChange (G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4bool &incidentHasChanged)
std::pair< G4int, G4double >	interpolateEnergy (G4double ke) const
G4int	sampleFlat (std::vector< G4double > sigma) const
void	CheckQnums (G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4double Q, G4double B, G4double S)
Protected Member Functions inherited from G4HadronicInteraction
void	SetModelName (const G4String &nam)
G4bool	IsBlocked () const
void	Block ()

Protected Attributes
G4RPGFragmentation	fragmentation
G4RPGTwoCluster	twoCluster
G4RPGPionSuppression	pionSuppression
G4RPGStrangeProduction	strangeProduction
G4RPGTwoBody	twoBody
G4ParticleDefinition *	particleDef [18]
Protected Attributes inherited from G4HadronicInteraction
G4HadFinalState	theParticleChange
G4int	verboseLevel
G4double	theMinEnergy
G4double	theMaxEnergy
G4bool	isBlocked

Detailed Description

Definition at line 53 of file G4RPGInelastic.hh.

Member Enumeration Documentation

anonymous enum

protected

Enumerator
pi0
pip
pim
kp
km
k0
k0b
pro
neu
lam
sp
s0
sm
xi0
xim
om
ap
an

Definition at line 121 of file G4RPGInelastic.hh.

        {pi0, pip, pim, kp, km, k0, k0b, pro, neu, 
         lam, sp, s0, sm, xi0, xim, om, ap, an};

Constructor & Destructor Documentation

G4RPGInelastic::G4RPGInelastic

(

const G4String &

modelName = "RPGInelastic"

)

Definition at line 38 of file G4RPGInelastic.cc.

References G4AntiKaonZero::AntiKaonZero(), G4AntiNeutron::AntiNeutron(), G4AntiProton::AntiProton(), G4cout, G4endl, G4KaonMinus::KaonMinus(), G4KaonPlus::KaonPlus(), G4KaonZero::KaonZero(), G4Lambda::Lambda(), G4Neutron::Neutron(), G4OmegaMinus::OmegaMinus(), particleDef, G4PionMinus::PionMinus(), G4PionPlus::PionPlus(), G4PionZero::PionZero(), G4Proton::Proton(), G4SigmaMinus::SigmaMinus(), G4SigmaPlus::SigmaPlus(), G4SigmaZero::SigmaZero(), G4XiMinus::XiMinus(), and G4XiZero::XiZero().

  : G4HadronicInteraction(modelName)
{
  cache = 0.0;
  particleDef[0] = G4PionZero::PionZero();
  particleDef[1] = G4PionPlus::PionPlus();
  particleDef[2] = G4PionMinus::PionMinus();
  particleDef[3] = G4KaonPlus::KaonPlus();
  particleDef[4] = G4KaonMinus::KaonMinus();
  particleDef[5] = G4KaonZero::KaonZero();
  particleDef[6] = G4AntiKaonZero::AntiKaonZero();
  particleDef[7] = G4Proton::Proton();
  particleDef[8] = G4Neutron::Neutron();
  particleDef[9] = G4Lambda::Lambda();
  particleDef[10] = G4SigmaPlus::SigmaPlus();
  particleDef[11] = G4SigmaZero::SigmaZero();
  particleDef[12] = G4SigmaMinus::SigmaMinus();
  particleDef[13] = G4XiZero::XiZero();
  particleDef[14] = G4XiMinus::XiMinus();
  particleDef[15] = G4OmegaMinus::OmegaMinus();
  particleDef[16] = G4AntiProton::AntiProton();
  particleDef[17] = G4AntiNeutron::AntiNeutron();

  G4cout << " **************************************************** " << G4endl; 
  G4cout << " * The RPG model is currently under development and * " << G4endl; 
  G4cout << " * should not be used.                              * " << G4endl; 
  G4cout << " **************************************************** " << G4endl; 
}

virtual G4RPGInelastic::~G4RPGInelastic ( )

inlinevirtual

Definition at line 59 of file G4RPGInelastic.hh.

   { }

Member Function Documentation

void G4RPGInelastic::CalculateMomenta ( G4FastVector< G4ReactionProduct, 256 > &  vec, G4int &  vecLen, const G4HadProjectile *  originalIncident, const G4DynamicParticle *  originalTarget, G4ReactionProduct &  modifiedOriginal, G4Nucleus &  targetNucleus, G4ReactionProduct &  currentParticle, G4ReactionProduct &  targetParticle, G4bool &  incidentHasChanged, G4bool &  targetHasChanged, G4bool  quasiElastic  )

protected

Definition at line 202 of file G4RPGInelastic.cc.

References G4Nucleus::AnnihilationEvaporationEffects(), G4Nucleus::Cinema(), fragmentation, G4cout, G4UniformRand, G4HadProjectile::Get4Momentum(), G4Nucleus::GetA_asInt(), G4HadProjectile::GetDefinition(), G4DynamicParticle::GetDefinition(), G4HadProjectile::GetKineticEnergy(), G4ReactionProduct::GetKineticEnergy(), G4ReactionProduct::GetMass(), G4ReactionProduct::GetMomentum(), G4ParticleDefinition::GetPDGEncoding(), G4ParticleDefinition::GetPDGMass(), G4ReactionProduct::GetTotalMomentum(), python.hepunit::GeV, G4FastVector< Type, N >::Initialize(), G4KaonMinus::KaonMinus(), G4KaonPlus::KaonPlus(), G4KaonZeroLong::KaonZeroLong(), G4KaonZeroShort::KaonZeroShort(), CLHEP::Hep3Vector::mag(), MarkLeadingStrangeParticle(), G4INCL::Math::max(), G4InuclParticleNames::pp, G4RPGTwoCluster::ReactionStage(), G4RPGTwoBody::ReactionStage(), G4RPGFragmentation::ReactionStage(), G4HadReentrentException::Report(), G4FastVector< Type, N >::SetElement(), G4ReactionProduct::SetKineticEnergy(), G4ReactionProduct::SetMomentum(), twoBody, twoCluster, and CLHEP::HepLorentzVector::vect().

Referenced by G4RPGPiMinusInelastic::ApplyYourself(), G4RPGPiPlusInelastic::ApplyYourself(), G4RPGKMinusInelastic::ApplyYourself(), G4RPGProtonInelastic::ApplyYourself(), G4RPGNeutronInelastic::ApplyYourself(), G4RPGXiMinusInelastic::ApplyYourself(), G4RPGAntiLambdaInelastic::ApplyYourself(), G4RPGAntiProtonInelastic::ApplyYourself(), G4RPGSigmaPlusInelastic::ApplyYourself(), G4RPGAntiSigmaPlusInelastic::ApplyYourself(), G4RPGOmegaMinusInelastic::ApplyYourself(), G4RPGAntiOmegaMinusInelastic::ApplyYourself(), G4RPGAntiXiMinusInelastic::ApplyYourself(), G4RPGAntiSigmaMinusInelastic::ApplyYourself(), G4RPGAntiNeutronInelastic::ApplyYourself(), G4RPGKZeroInelastic::ApplyYourself(), G4RPGLambdaInelastic::ApplyYourself(), G4RPGSigmaMinusInelastic::ApplyYourself(), G4RPGKPlusInelastic::ApplyYourself(), G4RPGAntiXiZeroInelastic::ApplyYourself(), G4RPGXiZeroInelastic::ApplyYourself(), and G4RPGAntiKZeroInelastic::ApplyYourself().

{
  cache = 0;
  what = originalIncident->Get4Momentum().vect();

  G4ReactionProduct leadingStrangeParticle;

  //  strangeProduction.ReactionStage(originalIncident, modifiedOriginal,
  //                                  incidentHasChanged, originalTarget,
  //                                  targetParticle, targetHasChanged,
  //                                  targetNucleus, currentParticle,
  //                                  vec, vecLen,
  //                  false, leadingStrangeParticle);

  if( quasiElastic )
  {
    twoBody.ReactionStage(originalIncident, modifiedOriginal,
                          incidentHasChanged, originalTarget,
                          targetParticle, targetHasChanged,
                          targetNucleus, currentParticle,
                          vec, vecLen,
                          false, leadingStrangeParticle);
    return;
  }

  G4bool leadFlag = MarkLeadingStrangeParticle(currentParticle,
                                               targetParticle,
                                               leadingStrangeParticle );
  //
  // Note: the number of secondaries can be reduced in GenerateXandPt 
  // and TwoCluster
  //
  G4bool finishedGenXPt = false;
  G4bool annihilation = false;
  if( originalIncident->GetDefinition()->GetPDGEncoding() < 0 &&
      currentParticle.GetMass() == 0.0 && targetParticle.GetMass() == 0.0 )
  {
    // original was an anti-particle and annihilation has taken place
    annihilation = true;
    G4double ekcor = 1.0;
    G4double ek = originalIncident->GetKineticEnergy();
    G4double ekOrg = ek;
      
    const G4double tarmas = originalTarget->GetDefinition()->GetPDGMass();
    if( ek > 1.0*GeV )ekcor = 1./(ek/GeV);
    const G4double atomicWeight = targetNucleus.GetA_asInt();
    ek = 2*tarmas + ek*(1.+ekcor/atomicWeight);
    G4double tkin = targetNucleus.Cinema(ek);
    ek += tkin;
    ekOrg += tkin;
    //      modifiedOriginal.SetKineticEnergy( ekOrg );
    //
    // evaporation --  re-calculate black track energies
    //                 this was Done already just before the cascade
    //
    tkin = targetNucleus.AnnihilationEvaporationEffects(ek, ekOrg);
    ekOrg -= tkin;
    ekOrg = std::max( 0.0001*GeV, ekOrg );
    modifiedOriginal.SetKineticEnergy( ekOrg );
    G4double amas = originalIncident->GetDefinition()->GetPDGMass();
    G4double et = ekOrg + amas;
    G4double p = std::sqrt( std::abs(et*et-amas*amas) );
    G4double pp = modifiedOriginal.GetMomentum().mag();
    if( pp > 0.0 )
    {
      G4ThreeVector momentum = modifiedOriginal.GetMomentum();
      modifiedOriginal.SetMomentum( momentum * (p/pp) );
    }
    if( ekOrg <= 0.0001 )
    {
      modifiedOriginal.SetKineticEnergy( 0.0 );
      modifiedOriginal.SetMomentum( 0.0, 0.0, 0.0 );
    }
  }

  // twsup gives percentage of time two-cluster model is called

  const G4double twsup[] = { 1.0, 0.7, 0.5, 0.3, 0.2, 0.1 };
  G4double rand1 = G4UniformRand();
  G4double rand2 = G4UniformRand();

  // Cache current, target, and secondaries
  G4ReactionProduct saveCurrent = currentParticle;
  G4ReactionProduct saveTarget = targetParticle;
  std::vector<G4ReactionProduct> savevec;
  for (G4int i = 0; i < vecLen; i++) savevec.push_back(*vec[i]);

  // Call fragmentation code if
  //   1) there is annihilation, or
  //   2) there are more than 5 secondaries, or
  //   3) incident KE is > 1 GeV AND
  //        ( incident is a kaon AND rand < 0.5 OR twsup )
  //

  if( annihilation || vecLen > 5 ||
      ( modifiedOriginal.GetKineticEnergy()/GeV >= 1.0 &&

      (((originalIncident->GetDefinition() == G4KaonPlus::KaonPlus() ||
         originalIncident->GetDefinition() == G4KaonMinus::KaonMinus() ||
         originalIncident->GetDefinition() == G4KaonZeroLong::KaonZeroLong() ||
         originalIncident->GetDefinition() == G4KaonZeroShort::KaonZeroShort()) &&
      rand1 < 0.5) 
       || rand2 > twsup[vecLen]) ) )

    finishedGenXPt =
      fragmentation.ReactionStage(originalIncident, modifiedOriginal,
                                  incidentHasChanged, originalTarget,
                                  targetParticle, targetHasChanged,
                                  targetNucleus, currentParticle,
                                  vec, vecLen,
                                  leadFlag, leadingStrangeParticle);

  if (finishedGenXPt) return;

  G4bool finishedTwoClu = false;

  if (modifiedOriginal.GetTotalMomentum() < 1.0) {
    for (G4int i = 0; i < vecLen; i++) delete vec[i];
    vecLen = 0;

  } else {
    // Occaisionally, GenerateXandPt will fail in the annihilation channel.
    // Restore current, target and secondaries to pre-GenerateXandPt state
    // before trying annihilation in TwoCluster

    if (!finishedGenXPt && annihilation) {
      currentParticle = saveCurrent;
      targetParticle = saveTarget;
      for (G4int i = 0; i < vecLen; i++) delete vec[i];
      vecLen = 0;
      vec.Initialize( 0 );
      for (G4int i = 0; i < G4int(savevec.size()); i++) {
        G4ReactionProduct* p = new G4ReactionProduct;
        *p = savevec[i];
        vec.SetElement( vecLen++, p );
      }
    }

    // Big violations of energy conservation in this method - don't use
    // 
    //    pionSuppression.ReactionStage(originalIncident, modifiedOriginal,
    //                                  incidentHasChanged, originalTarget,
    //                                  targetParticle, targetHasChanged,
    //                                  targetNucleus, currentParticle,
    //                                  vec, vecLen,
    //                                  false, leadingStrangeParticle);

    try
    {
      finishedTwoClu = 
        twoCluster.ReactionStage(originalIncident, modifiedOriginal,
                                 incidentHasChanged, originalTarget,
                                 targetParticle, targetHasChanged,
                                 targetNucleus, currentParticle,
                                 vec, vecLen,
                                 leadFlag, leadingStrangeParticle);
    }
     catch(G4HadReentrentException aC)
    {
       aC.Report(G4cout);
       throw G4HadReentrentException(__FILE__, __LINE__, "Failing to calculate momenta");
    }
  }

  if (finishedTwoClu) return;

  twoBody.ReactionStage(originalIncident, modifiedOriginal,
                        incidentHasChanged, originalTarget,
                        targetParticle, targetHasChanged,
                        targetNucleus, currentParticle,
                        vec, vecLen,
                        false, leadingStrangeParticle);
}

void G4RPGInelastic::CheckQnums ( G4FastVector< G4ReactionProduct, 256 > &  vec, G4int &  vecLen, G4ReactionProduct &  currentParticle, G4ReactionProduct &  targetParticle, G4double  Q, G4double  B, G4double  S  )

protected

Definition at line 545 of file G4RPGInelastic.cc.

References G4cout, G4endl, G4ParticleDefinition::GetAntiQuarkContent(), G4ParticleDefinition::GetBaryonNumber(), G4ReactionProduct::GetDefinition(), G4ParticleDefinition::GetParticleName(), G4ParticleDefinition::GetPDGCharge(), and G4ParticleDefinition::GetQuarkContent().

{
  G4ParticleDefinition* projDef = currentParticle.GetDefinition();
  G4ParticleDefinition* targDef = targetParticle.GetDefinition();
  G4double chargeSum = projDef->GetPDGCharge() + targDef->GetPDGCharge();
  G4double baryonSum = projDef->GetBaryonNumber() + targDef->GetBaryonNumber();
  G4double strangenessSum = projDef->GetQuarkContent(3) - 
                            projDef->GetAntiQuarkContent(3) + 
                            targDef->GetQuarkContent(3) -
                            targDef->GetAntiQuarkContent(3);

  G4ParticleDefinition* secDef = 0;
  for (G4int i = 0; i < vecLen; i++) {
    secDef = vec[i]->GetDefinition();
    chargeSum += secDef->GetPDGCharge();
    baryonSum += secDef->GetBaryonNumber();
    strangenessSum += secDef->GetQuarkContent(3) 
                    - secDef->GetAntiQuarkContent(3);
  }

  G4bool OK = true;
  if (chargeSum != Q) {
    G4cout << " Charge not conserved " << G4endl;
    OK = false;
  }
  if (baryonSum != B) {
    G4cout << " Baryon number not conserved " << G4endl;
    OK = false;
  }
  if (strangenessSum != S) {
    G4cout << " Strangeness not conserved " << G4endl;
    OK = false;
  } 

  if (!OK) {
    G4cout << " projectile: " << projDef->GetParticleName() 
           << "  target: " << targDef->GetParticleName() << G4endl;
    for (G4int i = 0; i < vecLen; i++) {
      secDef = vec[i]->GetDefinition();
      G4cout << secDef->GetParticleName() << " " ;
    }
    G4cout << G4endl;
  }

}

G4int G4RPGInelastic::Factorial ( G4int  n )

protected

Definition at line 86 of file G4RPGInelastic.cc.

References G4INCL::Math::min().

{
  G4int j = std::min(n,10);
  G4int result = 1;
  if (j <= 1) return result;
  for (G4int i = 2; i <= j; ++i) result *= i;
  return result;
}

void G4RPGInelastic::GetNormalizationConstant ( const G4double  availableEnergy, G4double &  n, G4double &  anpn  )

protected

Definition at line 158 of file G4RPGInelastic.cc.

References energy(), python.hepunit::GeV, G4INCL::Math::max(), G4INCL::Math::min(), n, python.hepunit::pi, and mcscore::test().

  {
    const G4double expxu =  82.;          // upper bound for arg. of exp
    const G4double expxl = -expxu;        // lower bound for arg. of exp
    const G4int numSec = 60;
    //
    // the only difference between the calculation for annihilation channels
    // and normal is the starting value, iBegin, for the loop below
    //
    G4int iBegin = 1;
    G4double en = energy;
    if( energy < 0.0 )
    {
      iBegin = 2;
      en *= -1.0;
    }
    //
    // number of total particles vs. centre of mass Energy - 2*proton mass
    //
    G4double aleab = std::log(en/GeV);
    n = 3.62567 + aleab*(0.665843 + aleab*(0.336514 + aleab*(0.117712 + 0.0136912*aleab)));
    n -= 2.0;
    //
    // normalization constant for kno-distribution
    //
    anpn = 0.0;
    G4double test, temp;
    for( G4int i=iBegin; i<=numSec; ++i )
    {
      temp = pi*i/(2.0*n*n);
      test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(i*i)/(n*n) ) ) );
      if( temp < 1.0 )
      {
        if( test >= 1.0e-10 )anpn += temp*test;
      }
      else
        anpn += temp*test;
    }
  }

std::pair< G4int, G4double > G4RPGInelastic::interpolateEnergy ( G4double  ke ) const

protected

Definition at line 506 of file G4RPGInelastic.cc.

Referenced by G4RPGNucleonInelastic::GetFSPartTypesForT0(), G4RPGNucleonInelastic::GetFSPartTypesForT1(), G4RPGPionInelastic::GetFSPartTypesForT12(), G4RPGPionInelastic::GetFSPartTypesForT32(), G4RPGNucleonInelastic::GetMultiplicityT0(), G4RPGNucleonInelastic::GetMultiplicityT1(), G4RPGPionInelastic::GetMultiplicityT12(), and G4RPGPionInelastic::GetMultiplicityT32().

{
  G4int index = 29;
  G4double fraction = 0.0;

  for (G4int i = 1; i < 30; i++) {
    if (e < energyScale[i]) {
      index = i-1;
      fraction = (e - energyScale[index]) / (energyScale[i] - energyScale[index]);
      break;
    }
  }
  return std::pair<G4int, G4double>(index, fraction);
}

G4bool G4RPGInelastic::MarkLeadingStrangeParticle ( const G4ReactionProduct &  currentParticle, const G4ReactionProduct &  targetParticle, G4ReactionProduct &  leadParticle  )

protected

Definition at line 96 of file G4RPGInelastic.cc.

References G4ReactionProduct::GetDefinition(), G4ReactionProduct::GetMass(), G4ParticleDefinition::GetPDGMass(), G4KaonPlus::KaonPlus(), G4Neutron::Neutron(), and G4Proton::Proton().

Referenced by CalculateMomenta().

{
  // The following was in GenerateXandPt and TwoCluster.
  // Add a parameter to the GenerateXandPt function telling it about the 
  // strange particle.
  //
  // Assumes that the original particle was a strange particle
  //
  G4bool lead = false;
  if( (currentParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
      (currentParticle.GetDefinition() != G4Proton::Proton()) &&
      (currentParticle.GetDefinition() != G4Neutron::Neutron()) )
  {
    lead = true;
    leadParticle = currentParticle;              //  set lead to the incident particle
  }
  else if( (targetParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
           (targetParticle.GetDefinition() != G4Proton::Proton()) &&
           (targetParticle.GetDefinition() != G4Neutron::Neutron()) )
  {
    lead = true;
    leadParticle = targetParticle;              //   set lead to the target particle
  }
  return lead;
}

G4double G4RPGInelastic::Pmltpc ( G4int  np, G4int  nm, G4int  nz, G4int  n, G4double  b, G4double  c  )

protected

Definition at line 68 of file G4RPGInelastic.cc.

References G4INCL::Math::max(), and G4INCL::Math::min().

{
  const G4double expxu =  82.;           // upper bound for arg. of exp
  const G4double expxl = -expxu;         // lower bound for arg. of exp
  G4double npf = 0.0;
  G4double nmf = 0.0;
  G4double nzf = 0.0;
  G4int i;
  for( i=2; i<=np; i++ )npf += std::log((double)i);
  for( i=2; i<=nneg; i++ )nmf += std::log((double)i);
  for( i=2; i<=nz; i++ )nzf += std::log((double)i);
  G4double r;
  r = std::min( expxu, std::max( expxl, -(np-nneg+nz+b)*(np-nneg+nz+b)/(2*c*c*n*n)-npf-nmf-nzf ) );
  return std::exp(r);
}

G4int G4RPGInelastic::sampleFlat ( std::vector< G4double >  sigma ) const

protected

Definition at line 523 of file G4RPGInelastic.cc.

References G4UniformRand.

Referenced by G4RPGNucleonInelastic::GetFSPartTypesForT0(), G4RPGNucleonInelastic::GetFSPartTypesForT1(), G4RPGPionInelastic::GetFSPartTypesForT12(), G4RPGPionInelastic::GetFSPartTypesForT32(), G4RPGNucleonInelastic::GetMultiplicityT0(), G4RPGNucleonInelastic::GetMultiplicityT1(), G4RPGPionInelastic::GetMultiplicityT12(), and G4RPGPionInelastic::GetMultiplicityT32().

{
  G4int i;
  G4double sum(0.);
  for (i = 0; i < G4int(sigma.size()); i++) sum += sigma[i];

  G4double fsum = sum*G4UniformRand();
  G4double partialSum = 0.0;
  G4int channel = 0;

  for (i = 0; i < G4int(sigma.size()); i++) {
    partialSum += sigma[i];
    if (fsum < partialSum) {
      channel = i;
      break;
    }
  }

  return channel;
}

void G4RPGInelastic::SetUpChange ( G4FastVector< G4ReactionProduct, 256 > &  vec, G4int &  vecLen, G4ReactionProduct &  currentParticle, G4ReactionProduct &  targetParticle, G4bool &  incidentHasChanged  )

protected

Definition at line 403 of file G4RPGInelastic.cc.

References G4HadFinalState::AddSecondary(), G4HadFinalState::Clear(), DBL_MIN, python.hepunit::eV, G4UniformRand, G4ReactionProduct::GetDefinition(), G4ReactionProduct::GetKineticEnergy(), G4ReactionProduct::GetMass(), G4ReactionProduct::GetMomentum(), k0, k0b, G4KaonZeroLong::KaonZeroLong(), G4KaonZeroShort::KaonZeroShort(), CLHEP::Hep3Vector::mag(), python.hepunit::MeV, particleDef, CLHEP::Hep3Vector::rotate(), G4DynamicParticle::SetDefinition(), G4ReactionProduct::SetDefinitionAndUpdateE(), G4HadFinalState::SetEnergyChange(), G4DynamicParticle::SetMomentum(), G4HadFinalState::SetMomentumChange(), G4HadFinalState::SetStatusChange(), stopAndKill, G4HadronicInteraction::theParticleChange, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

Referenced by G4RPGPiMinusInelastic::ApplyYourself(), G4RPGPiPlusInelastic::ApplyYourself(), G4RPGKMinusInelastic::ApplyYourself(), G4RPGProtonInelastic::ApplyYourself(), G4RPGNeutronInelastic::ApplyYourself(), G4RPGXiMinusInelastic::ApplyYourself(), G4RPGAntiLambdaInelastic::ApplyYourself(), G4RPGAntiProtonInelastic::ApplyYourself(), G4RPGSigmaPlusInelastic::ApplyYourself(), G4RPGAntiSigmaPlusInelastic::ApplyYourself(), G4RPGOmegaMinusInelastic::ApplyYourself(), G4RPGAntiOmegaMinusInelastic::ApplyYourself(), G4RPGAntiXiMinusInelastic::ApplyYourself(), G4RPGAntiSigmaMinusInelastic::ApplyYourself(), G4RPGAntiNeutronInelastic::ApplyYourself(), G4RPGKZeroInelastic::ApplyYourself(), G4RPGLambdaInelastic::ApplyYourself(), G4RPGSigmaMinusInelastic::ApplyYourself(), G4RPGKPlusInelastic::ApplyYourself(), G4RPGAntiXiZeroInelastic::ApplyYourself(), G4RPGXiZeroInelastic::ApplyYourself(), and G4RPGAntiKZeroInelastic::ApplyYourself().

{
  theParticleChange.Clear();
  G4ParticleDefinition* aKaonZL = G4KaonZeroLong::KaonZeroLong();
  G4ParticleDefinition* aKaonZS = G4KaonZeroShort::KaonZeroShort();
  G4int i;

  if (currentParticle.GetDefinition() == particleDef[k0]) {
    if (G4UniformRand() < 0.5) {
      currentParticle.SetDefinitionAndUpdateE(aKaonZL);
      incidentHasChanged = true;
    } else {
      currentParticle.SetDefinitionAndUpdateE(aKaonZS);
    }
  } else if (currentParticle.GetDefinition() == particleDef[k0b]) {
    if (G4UniformRand() < 0.5) {
      currentParticle.SetDefinitionAndUpdateE(aKaonZL);
    } else {
      currentParticle.SetDefinitionAndUpdateE(aKaonZS);
      incidentHasChanged = true;
    }
  }

  if (targetParticle.GetDefinition() == particleDef[k0] || 
      targetParticle.GetDefinition() == particleDef[k0b] ) {
    if (G4UniformRand() < 0.5) {
      targetParticle.SetDefinitionAndUpdateE(aKaonZL);
    } else {
      targetParticle.SetDefinitionAndUpdateE(aKaonZS);
    }
  }

  for (i = 0; i < vecLen; ++i) {
    if (vec[i]->GetDefinition() == particleDef[k0] ||
        vec[i]->GetDefinition() == particleDef[k0b] ) {
      if (G4UniformRand() < 0.5) {
        vec[i]->SetDefinitionAndUpdateE(aKaonZL);
      } else {
        vec[i]->SetDefinitionAndUpdateE(aKaonZS);
      }
    }
  }

  if (incidentHasChanged) {
    G4DynamicParticle* p0 = new G4DynamicParticle;
    p0->SetDefinition(currentParticle.GetDefinition() );
    p0->SetMomentum(currentParticle.GetMomentum() );
    theParticleChange.AddSecondary( p0 );
    theParticleChange.SetStatusChange( stopAndKill );
    theParticleChange.SetEnergyChange( 0.0 );

  } else {
    G4double p = currentParticle.GetMomentum().mag()/MeV;
    G4ThreeVector mom = currentParticle.GetMomentum();
    if (p > DBL_MIN)
      theParticleChange.SetMomentumChange(mom.x()/p, mom.y()/p, mom.z()/p );
    else
      theParticleChange.SetMomentumChange(0.0, 0.0, 1.0);

    G4double aE = currentParticle.GetKineticEnergy();
    if (std::fabs(aE)<.1*eV) aE=.1*eV;
    theParticleChange.SetEnergyChange( aE );
  }

  if (targetParticle.GetMass() > 0.0)  // Tgt particle can be eliminated in TwoBody
  {
    G4ThreeVector momentum = targetParticle.GetMomentum();
    momentum = momentum.rotate(cache, what);
    G4double targKE = targetParticle.GetKineticEnergy();
    G4ThreeVector dir(0.0, 0.0, 1.0);
    if (targKE < DBL_MIN)
      targKE = DBL_MIN;
    else
      dir = momentum/momentum.mag();

    G4DynamicParticle* p1 = 
      new G4DynamicParticle(targetParticle.GetDefinition(), dir, targKE);

    theParticleChange.AddSecondary( p1 );
  }

  G4DynamicParticle* p;
  for (i = 0; i < vecLen; ++i) {
    G4double secKE = vec[i]->GetKineticEnergy();
    G4ThreeVector momentum = vec[i]->GetMomentum();
    G4ThreeVector dir(0.0, 0.0, 1.0);
    if (secKE < DBL_MIN)
      secKE = DBL_MIN;
    else
      dir = momentum/momentum.mag();

    p = new G4DynamicParticle(vec[i]->GetDefinition(), dir, secKE);
    theParticleChange.AddSecondary( p );
    delete vec[i];
  }
}

void G4RPGInelastic::SetUpPions ( const G4int  np, const G4int  nm, const G4int  nz, G4FastVector< G4ReactionProduct, 256 > &  vec, G4int &  vecLen  )

protected

Definition at line 126 of file G4RPGInelastic.cc.

References G4UniformRand, G4PionMinus::PionMinus(), G4PionPlus::PionPlus(), G4PionZero::PionZero(), G4ReactionProduct::SetDefinition(), G4FastVector< Type, N >::SetElement(), and G4ReactionProduct::SetSide().

 {
   if( np+nneg+nz == 0 )return;
   G4int i;
   G4ReactionProduct *p;
   for( i=0; i<np; ++i )
   {
     p = new G4ReactionProduct;
     p->SetDefinition( G4PionPlus::PionPlus() );
     (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
     vec.SetElement( vecLen++, p );
   }
   for( i=np; i<np+nneg; ++i )
   {
     p = new G4ReactionProduct;
     p->SetDefinition( G4PionMinus::PionMinus() );
     (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
     vec.SetElement( vecLen++, p );
   }
   for( i=np+nneg; i<np+nneg+nz; ++i )
   {
     p = new G4ReactionProduct;
     p->SetDefinition( G4PionZero::PionZero() );
     (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
     vec.SetElement( vecLen++, p );
   }
 }

Field Documentation

G4RPGFragmentation G4RPGInelastic::fragmentation

protected

Definition at line 101 of file G4RPGInelastic.hh.

Referenced by CalculateMomenta().

G4ParticleDefinition* G4RPGInelastic::particleDef[18]

protected

Definition at line 126 of file G4RPGInelastic.hh.

Referenced by G4RPGInelastic(), and SetUpChange().

G4RPGPionSuppression G4RPGInelastic::pionSuppression

protected

Definition at line 105 of file G4RPGInelastic.hh.

G4RPGStrangeProduction G4RPGInelastic::strangeProduction

protected

Definition at line 107 of file G4RPGInelastic.hh.

G4RPGTwoBody G4RPGInelastic::twoBody

protected

Definition at line 109 of file G4RPGInelastic.hh.

Referenced by CalculateMomenta().

G4RPGTwoCluster G4RPGInelastic::twoCluster

protected

Definition at line 103 of file G4RPGInelastic.hh.

Referenced by CalculateMomenta().

---

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

• G4RPGInelastic.hh
• G4RPGInelastic.cc