G4RPGInelastic Class Reference

#include <G4RPGInelastic.hh>

Inheritance diagram for G4RPGInelastic:

G4HadronicInteraction G4RPGAntiKZeroInelastic G4RPGAntiLambdaInelastic G4RPGAntiNeutronInelastic G4RPGAntiOmegaMinusInelastic G4RPGAntiProtonInelastic G4RPGAntiSigmaMinusInelastic G4RPGAntiSigmaPlusInelastic G4RPGAntiXiMinusInelastic G4RPGAntiXiZeroInelastic G4RPGKLongInelastic G4RPGKMinusInelastic G4RPGKPlusInelastic G4RPGKShortInelastic G4RPGKZeroInelastic G4RPGLambdaInelastic G4RPGNucleonInelastic G4RPGOmegaMinusInelastic G4RPGPionInelastic G4RPGSigmaMinusInelastic G4RPGSigmaPlusInelastic G4RPGXiMinusInelastic G4RPGXiZeroInelastic

Public Member Functions

 G4RPGInelastic (const G4String &modelName="RPGInelastic")
virtual ~G4RPGInelastic ()

Protected Types

 pi0
 pip
 pim
 kp
 km
 k0
 k0b
 pro
 neu
 lam
 sp
 s0
 sm
 xi0
 xim
 om
 ap
 an
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, G4doubleinterpolateEnergy (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 Attributes

G4RPGFragmentation fragmentation
G4RPGTwoCluster twoCluster
G4RPGPionSuppression pionSuppression
G4RPGStrangeProduction strangeProduction
G4RPGTwoBody twoBody
G4ParticleDefinitionparticleDef [18]

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.

00121         {pi0, pip, pim, kp, km, k0, k0b, pro, neu, 
00122          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().

00039   : G4HadronicInteraction(modelName)
00040 {
00041   cache = 0.0;
00042   particleDef[0] = G4PionZero::PionZero();
00043   particleDef[1] = G4PionPlus::PionPlus();
00044   particleDef[2] = G4PionMinus::PionMinus();
00045   particleDef[3] = G4KaonPlus::KaonPlus();
00046   particleDef[4] = G4KaonMinus::KaonMinus();
00047   particleDef[5] = G4KaonZero::KaonZero();
00048   particleDef[6] = G4AntiKaonZero::AntiKaonZero();
00049   particleDef[7] = G4Proton::Proton();
00050   particleDef[8] = G4Neutron::Neutron();
00051   particleDef[9] = G4Lambda::Lambda();
00052   particleDef[10] = G4SigmaPlus::SigmaPlus();
00053   particleDef[11] = G4SigmaZero::SigmaZero();
00054   particleDef[12] = G4SigmaMinus::SigmaMinus();
00055   particleDef[13] = G4XiZero::XiZero();
00056   particleDef[14] = G4XiMinus::XiMinus();
00057   particleDef[15] = G4OmegaMinus::OmegaMinus();
00058   particleDef[16] = G4AntiProton::AntiProton();
00059   particleDef[17] = G4AntiNeutron::AntiNeutron();
00060 
00061   G4cout << " **************************************************** " << G4endl; 
00062   G4cout << " * The RPG model is currently under development and * " << G4endl; 
00063   G4cout << " * should not be used.                              * " << G4endl; 
00064   G4cout << " **************************************************** " << G4endl; 
00065 }

virtual G4RPGInelastic::~G4RPGInelastic (  )  [inline, virtual]

Definition at line 59 of file G4RPGInelastic.hh.

00060    { }


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(), G4DynamicParticle::GetDefinition(), G4HadProjectile::GetDefinition(), G4ReactionProduct::GetKineticEnergy(), G4HadProjectile::GetKineticEnergy(), G4ReactionProduct::GetMass(), G4ReactionProduct::GetMomentum(), G4ParticleDefinition::GetPDGEncoding(), G4ParticleDefinition::GetPDGMass(), G4ReactionProduct::GetTotalMomentum(), G4FastVector< Type, N >::Initialize(), G4KaonMinus::KaonMinus(), G4KaonPlus::KaonPlus(), G4KaonZeroLong::KaonZeroLong(), G4KaonZeroShort::KaonZeroShort(), MarkLeadingStrangeParticle(), G4InuclParticleNames::pp, G4RPGTwoCluster::ReactionStage(), G4RPGFragmentation::ReactionStage(), G4RPGTwoBody::ReactionStage(), G4HadReentrentException::Report(), G4FastVector< Type, N >::SetElement(), G4ReactionProduct::SetKineticEnergy(), G4ReactionProduct::SetMomentum(), twoBody, and twoCluster.

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

00213 {
00214   cache = 0;
00215   what = originalIncident->Get4Momentum().vect();
00216 
00217   G4ReactionProduct leadingStrangeParticle;
00218 
00219   //  strangeProduction.ReactionStage(originalIncident, modifiedOriginal,
00220   //                                  incidentHasChanged, originalTarget,
00221   //                                  targetParticle, targetHasChanged,
00222   //                                  targetNucleus, currentParticle,
00223   //                                  vec, vecLen,
00224   //                              false, leadingStrangeParticle);
00225 
00226   if( quasiElastic )
00227   {
00228     twoBody.ReactionStage(originalIncident, modifiedOriginal,
00229                           incidentHasChanged, originalTarget,
00230                           targetParticle, targetHasChanged,
00231                           targetNucleus, currentParticle,
00232                           vec, vecLen,
00233                           false, leadingStrangeParticle);
00234     return;
00235   }
00236 
00237   G4bool leadFlag = MarkLeadingStrangeParticle(currentParticle,
00238                                                targetParticle,
00239                                                leadingStrangeParticle );
00240   //
00241   // Note: the number of secondaries can be reduced in GenerateXandPt 
00242   // and TwoCluster
00243   //
00244   G4bool finishedGenXPt = false;
00245   G4bool annihilation = false;
00246   if( originalIncident->GetDefinition()->GetPDGEncoding() < 0 &&
00247       currentParticle.GetMass() == 0.0 && targetParticle.GetMass() == 0.0 )
00248   {
00249     // original was an anti-particle and annihilation has taken place
00250     annihilation = true;
00251     G4double ekcor = 1.0;
00252     G4double ek = originalIncident->GetKineticEnergy();
00253     G4double ekOrg = ek;
00254       
00255     const G4double tarmas = originalTarget->GetDefinition()->GetPDGMass();
00256     if( ek > 1.0*GeV )ekcor = 1./(ek/GeV);
00257     const G4double atomicWeight = targetNucleus.GetA_asInt();
00258     ek = 2*tarmas + ek*(1.+ekcor/atomicWeight);
00259     G4double tkin = targetNucleus.Cinema(ek);
00260     ek += tkin;
00261     ekOrg += tkin;
00262     //      modifiedOriginal.SetKineticEnergy( ekOrg );
00263     //
00264     // evaporation --  re-calculate black track energies
00265     //                 this was Done already just before the cascade
00266     //
00267     tkin = targetNucleus.AnnihilationEvaporationEffects(ek, ekOrg);
00268     ekOrg -= tkin;
00269     ekOrg = std::max( 0.0001*GeV, ekOrg );
00270     modifiedOriginal.SetKineticEnergy( ekOrg );
00271     G4double amas = originalIncident->GetDefinition()->GetPDGMass();
00272     G4double et = ekOrg + amas;
00273     G4double p = std::sqrt( std::abs(et*et-amas*amas) );
00274     G4double pp = modifiedOriginal.GetMomentum().mag();
00275     if( pp > 0.0 )
00276     {
00277       G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00278       modifiedOriginal.SetMomentum( momentum * (p/pp) );
00279     }
00280     if( ekOrg <= 0.0001 )
00281     {
00282       modifiedOriginal.SetKineticEnergy( 0.0 );
00283       modifiedOriginal.SetMomentum( 0.0, 0.0, 0.0 );
00284     }
00285   }
00286 
00287   // twsup gives percentage of time two-cluster model is called
00288 
00289   const G4double twsup[] = { 1.0, 0.7, 0.5, 0.3, 0.2, 0.1 };
00290   G4double rand1 = G4UniformRand();
00291   G4double rand2 = G4UniformRand();
00292 
00293   // Cache current, target, and secondaries
00294   G4ReactionProduct saveCurrent = currentParticle;
00295   G4ReactionProduct saveTarget = targetParticle;
00296   std::vector<G4ReactionProduct> savevec;
00297   for (G4int i = 0; i < vecLen; i++) savevec.push_back(*vec[i]);
00298 
00299   // Call fragmentation code if
00300   //   1) there is annihilation, or
00301   //   2) there are more than 5 secondaries, or
00302   //   3) incident KE is > 1 GeV AND
00303   //        ( incident is a kaon AND rand < 0.5 OR twsup )
00304   //
00305 
00306   if( annihilation || vecLen > 5 ||
00307       ( modifiedOriginal.GetKineticEnergy()/GeV >= 1.0 &&
00308 
00309       (((originalIncident->GetDefinition() == G4KaonPlus::KaonPlus() ||
00310          originalIncident->GetDefinition() == G4KaonMinus::KaonMinus() ||
00311          originalIncident->GetDefinition() == G4KaonZeroLong::KaonZeroLong() ||
00312          originalIncident->GetDefinition() == G4KaonZeroShort::KaonZeroShort()) &&
00313           rand1 < 0.5) 
00314        || rand2 > twsup[vecLen]) ) )
00315 
00316     finishedGenXPt =
00317       fragmentation.ReactionStage(originalIncident, modifiedOriginal,
00318                                   incidentHasChanged, originalTarget,
00319                                   targetParticle, targetHasChanged,
00320                                   targetNucleus, currentParticle,
00321                                   vec, vecLen,
00322                                   leadFlag, leadingStrangeParticle);
00323 
00324   if (finishedGenXPt) return;
00325 
00326   G4bool finishedTwoClu = false;
00327 
00328   if (modifiedOriginal.GetTotalMomentum() < 1.0) {
00329     for (G4int i = 0; i < vecLen; i++) delete vec[i];
00330     vecLen = 0;
00331 
00332   } else {
00333     // Occaisionally, GenerateXandPt will fail in the annihilation channel.
00334     // Restore current, target and secondaries to pre-GenerateXandPt state
00335     // before trying annihilation in TwoCluster
00336 
00337     if (!finishedGenXPt && annihilation) {
00338       currentParticle = saveCurrent;
00339       targetParticle = saveTarget;
00340       for (G4int i = 0; i < vecLen; i++) delete vec[i];
00341       vecLen = 0;
00342       vec.Initialize( 0 );
00343       for (G4int i = 0; i < G4int(savevec.size()); i++) {
00344         G4ReactionProduct* p = new G4ReactionProduct;
00345         *p = savevec[i];
00346         vec.SetElement( vecLen++, p );
00347       }
00348     }
00349 
00350     // Big violations of energy conservation in this method - don't use
00351     // 
00352     //    pionSuppression.ReactionStage(originalIncident, modifiedOriginal,
00353     //                                  incidentHasChanged, originalTarget,
00354     //                                  targetParticle, targetHasChanged,
00355     //                                  targetNucleus, currentParticle,
00356     //                                  vec, vecLen,
00357     //                                  false, leadingStrangeParticle);
00358 
00359     try
00360     {
00361       finishedTwoClu = 
00362         twoCluster.ReactionStage(originalIncident, modifiedOriginal,
00363                                  incidentHasChanged, originalTarget,
00364                                  targetParticle, targetHasChanged,
00365                                  targetNucleus, currentParticle,
00366                                  vec, vecLen,
00367                                  leadFlag, leadingStrangeParticle);
00368     }
00369      catch(G4HadReentrentException aC)
00370     {
00371        aC.Report(G4cout);
00372        throw G4HadReentrentException(__FILE__, __LINE__, "Failing to calculate momenta");
00373     }
00374   }
00375 
00376   if (finishedTwoClu) return;
00377 
00378   twoBody.ReactionStage(originalIncident, modifiedOriginal,
00379                         incidentHasChanged, originalTarget,
00380                         targetParticle, targetHasChanged,
00381                         targetNucleus, currentParticle,
00382                         vec, vecLen,
00383                         false, leadingStrangeParticle);
00384 }

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().

00550 {
00551   G4ParticleDefinition* projDef = currentParticle.GetDefinition();
00552   G4ParticleDefinition* targDef = targetParticle.GetDefinition();
00553   G4double chargeSum = projDef->GetPDGCharge() + targDef->GetPDGCharge();
00554   G4double baryonSum = projDef->GetBaryonNumber() + targDef->GetBaryonNumber();
00555   G4double strangenessSum = projDef->GetQuarkContent(3) - 
00556                             projDef->GetAntiQuarkContent(3) + 
00557                             targDef->GetQuarkContent(3) -
00558                             targDef->GetAntiQuarkContent(3);
00559 
00560   G4ParticleDefinition* secDef = 0;
00561   for (G4int i = 0; i < vecLen; i++) {
00562     secDef = vec[i]->GetDefinition();
00563     chargeSum += secDef->GetPDGCharge();
00564     baryonSum += secDef->GetBaryonNumber();
00565     strangenessSum += secDef->GetQuarkContent(3) 
00566                     - secDef->GetAntiQuarkContent(3);
00567   }
00568 
00569   G4bool OK = true;
00570   if (chargeSum != Q) {
00571     G4cout << " Charge not conserved " << G4endl;
00572     OK = false;
00573   }
00574   if (baryonSum != B) {
00575     G4cout << " Baryon number not conserved " << G4endl;
00576     OK = false;
00577   }
00578   if (strangenessSum != S) {
00579     G4cout << " Strangeness not conserved " << G4endl;
00580     OK = false;
00581   } 
00582 
00583   if (!OK) {
00584     G4cout << " projectile: " << projDef->GetParticleName() 
00585            << "  target: " << targDef->GetParticleName() << G4endl;
00586     for (G4int i = 0; i < vecLen; i++) {
00587       secDef = vec[i]->GetDefinition();
00588       G4cout << secDef->GetParticleName() << " " ;
00589     }
00590     G4cout << G4endl;
00591   }
00592 
00593 }

G4int G4RPGInelastic::Factorial ( G4int  n  )  [protected]

Definition at line 86 of file G4RPGInelastic.cc.

00087 {
00088   G4int j = std::min(n,10);
00089   G4int result = 1;
00090   if (j <= 1) return result;
00091   for (G4int i = 2; i <= j; ++i) result *= i;
00092   return result;
00093 }

void G4RPGInelastic::GetNormalizationConstant ( const G4double  availableEnergy,
G4double n,
G4double anpn 
) [protected]

Definition at line 158 of file G4RPGInelastic.cc.

References G4INCL::Math::pi.

00162   {
00163     const G4double expxu =  82.;          // upper bound for arg. of exp
00164     const G4double expxl = -expxu;        // lower bound for arg. of exp
00165     const G4int numSec = 60;
00166     //
00167     // the only difference between the calculation for annihilation channels
00168     // and normal is the starting value, iBegin, for the loop below
00169     //
00170     G4int iBegin = 1;
00171     G4double en = energy;
00172     if( energy < 0.0 )
00173     {
00174       iBegin = 2;
00175       en *= -1.0;
00176     }
00177     //
00178     // number of total particles vs. centre of mass Energy - 2*proton mass
00179     //
00180     G4double aleab = std::log(en/GeV);
00181     n = 3.62567 + aleab*(0.665843 + aleab*(0.336514 + aleab*(0.117712 + 0.0136912*aleab)));
00182     n -= 2.0;
00183     //
00184     // normalization constant for kno-distribution
00185     //
00186     anpn = 0.0;
00187     G4double test, temp;
00188     for( G4int i=iBegin; i<=numSec; ++i )
00189     {
00190       temp = pi*i/(2.0*n*n);
00191       test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(i*i)/(n*n) ) ) );
00192       if( temp < 1.0 )
00193       {
00194         if( test >= 1.0e-10 )anpn += temp*test;
00195       }
00196       else
00197         anpn += temp*test;
00198     }
00199   }

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().

00507 {
00508   G4int index = 29;
00509   G4double fraction = 0.0;
00510 
00511   for (G4int i = 1; i < 30; i++) {
00512     if (e < energyScale[i]) {
00513       index = i-1;
00514       fraction = (e - energyScale[index]) / (energyScale[i] - energyScale[index]);
00515       break;
00516     }
00517   }
00518   return std::pair<G4int, G4double>(index, fraction);
00519 }

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().

00100 {
00101   // The following was in GenerateXandPt and TwoCluster.
00102   // Add a parameter to the GenerateXandPt function telling it about the 
00103   // strange particle.
00104   //
00105   // Assumes that the original particle was a strange particle
00106   //
00107   G4bool lead = false;
00108   if( (currentParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
00109       (currentParticle.GetDefinition() != G4Proton::Proton()) &&
00110       (currentParticle.GetDefinition() != G4Neutron::Neutron()) )
00111   {
00112     lead = true;
00113     leadParticle = currentParticle;              //  set lead to the incident particle
00114   }
00115   else if( (targetParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
00116            (targetParticle.GetDefinition() != G4Proton::Proton()) &&
00117            (targetParticle.GetDefinition() != G4Neutron::Neutron()) )
00118   {
00119     lead = true;
00120     leadParticle = targetParticle;              //   set lead to the target particle
00121   }
00122   return lead;
00123 }

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

Definition at line 68 of file G4RPGInelastic.cc.

00070 {
00071   const G4double expxu =  82.;           // upper bound for arg. of exp
00072   const G4double expxl = -expxu;         // lower bound for arg. of exp
00073   G4double npf = 0.0;
00074   G4double nmf = 0.0;
00075   G4double nzf = 0.0;
00076   G4int i;
00077   for( i=2; i<=np; i++ )npf += std::log((double)i);
00078   for( i=2; i<=nneg; i++ )nmf += std::log((double)i);
00079   for( i=2; i<=nz; i++ )nzf += std::log((double)i);
00080   G4double r;
00081   r = std::min( expxu, std::max( expxl, -(np-nneg+nz+b)*(np-nneg+nz+b)/(2*c*c*n*n)-npf-nmf-nzf ) );
00082   return std::exp(r);
00083 }

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().

00524 {
00525   G4int i;
00526   G4double sum(0.);
00527   for (i = 0; i < G4int(sigma.size()); i++) sum += sigma[i];
00528 
00529   G4double fsum = sum*G4UniformRand();
00530   G4double partialSum = 0.0;
00531   G4int channel = 0;
00532 
00533   for (i = 0; i < G4int(sigma.size()); i++) {
00534     partialSum += sigma[i];
00535     if (fsum < partialSum) {
00536       channel = i;
00537       break;
00538     }
00539   }
00540 
00541   return channel;
00542 }

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, G4UniformRand, G4ReactionProduct::GetDefinition(), G4ReactionProduct::GetKineticEnergy(), G4ReactionProduct::GetMass(), G4ReactionProduct::GetMomentum(), k0, k0b, G4KaonZeroLong::KaonZeroLong(), G4KaonZeroShort::KaonZeroShort(), particleDef, G4DynamicParticle::SetDefinition(), G4ReactionProduct::SetDefinitionAndUpdateE(), G4HadFinalState::SetEnergyChange(), G4DynamicParticle::SetMomentum(), G4HadFinalState::SetMomentumChange(), G4HadFinalState::SetStatusChange(), stopAndKill, and G4HadronicInteraction::theParticleChange.

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

00408 {
00409   theParticleChange.Clear();
00410   G4ParticleDefinition* aKaonZL = G4KaonZeroLong::KaonZeroLong();
00411   G4ParticleDefinition* aKaonZS = G4KaonZeroShort::KaonZeroShort();
00412   G4int i;
00413 
00414   if (currentParticle.GetDefinition() == particleDef[k0]) {
00415     if (G4UniformRand() < 0.5) {
00416       currentParticle.SetDefinitionAndUpdateE(aKaonZL);
00417       incidentHasChanged = true;
00418     } else {
00419       currentParticle.SetDefinitionAndUpdateE(aKaonZS);
00420     }
00421   } else if (currentParticle.GetDefinition() == particleDef[k0b]) {
00422     if (G4UniformRand() < 0.5) {
00423       currentParticle.SetDefinitionAndUpdateE(aKaonZL);
00424     } else {
00425       currentParticle.SetDefinitionAndUpdateE(aKaonZS);
00426       incidentHasChanged = true;
00427     }
00428   }
00429 
00430   if (targetParticle.GetDefinition() == particleDef[k0] || 
00431       targetParticle.GetDefinition() == particleDef[k0b] ) {
00432     if (G4UniformRand() < 0.5) {
00433       targetParticle.SetDefinitionAndUpdateE(aKaonZL);
00434     } else {
00435       targetParticle.SetDefinitionAndUpdateE(aKaonZS);
00436     }
00437   }
00438 
00439   for (i = 0; i < vecLen; ++i) {
00440     if (vec[i]->GetDefinition() == particleDef[k0] ||
00441         vec[i]->GetDefinition() == particleDef[k0b] ) {
00442       if (G4UniformRand() < 0.5) {
00443         vec[i]->SetDefinitionAndUpdateE(aKaonZL);
00444       } else {
00445         vec[i]->SetDefinitionAndUpdateE(aKaonZS);
00446       }
00447     }
00448   }
00449 
00450   if (incidentHasChanged) {
00451     G4DynamicParticle* p0 = new G4DynamicParticle;
00452     p0->SetDefinition(currentParticle.GetDefinition() );
00453     p0->SetMomentum(currentParticle.GetMomentum() );
00454     theParticleChange.AddSecondary( p0 );
00455     theParticleChange.SetStatusChange( stopAndKill );
00456     theParticleChange.SetEnergyChange( 0.0 );
00457 
00458   } else {
00459     G4double p = currentParticle.GetMomentum().mag()/MeV;
00460     G4ThreeVector mom = currentParticle.GetMomentum();
00461     if (p > DBL_MIN)
00462       theParticleChange.SetMomentumChange(mom.x()/p, mom.y()/p, mom.z()/p );
00463     else
00464       theParticleChange.SetMomentumChange(0.0, 0.0, 1.0);
00465 
00466     G4double aE = currentParticle.GetKineticEnergy();
00467     if (std::fabs(aE)<.1*eV) aE=.1*eV;
00468     theParticleChange.SetEnergyChange( aE );
00469   }
00470 
00471   if (targetParticle.GetMass() > 0.0)  // Tgt particle can be eliminated in TwoBody
00472   {
00473     G4ThreeVector momentum = targetParticle.GetMomentum();
00474     momentum = momentum.rotate(cache, what);
00475     G4double targKE = targetParticle.GetKineticEnergy();
00476     G4ThreeVector dir(0.0, 0.0, 1.0);
00477     if (targKE < DBL_MIN)
00478       targKE = DBL_MIN;
00479     else
00480       dir = momentum/momentum.mag();
00481 
00482     G4DynamicParticle* p1 = 
00483       new G4DynamicParticle(targetParticle.GetDefinition(), dir, targKE);
00484 
00485     theParticleChange.AddSecondary( p1 );
00486   }
00487 
00488   G4DynamicParticle* p;
00489   for (i = 0; i < vecLen; ++i) {
00490     G4double secKE = vec[i]->GetKineticEnergy();
00491     G4ThreeVector momentum = vec[i]->GetMomentum();
00492     G4ThreeVector dir(0.0, 0.0, 1.0);
00493     if (secKE < DBL_MIN)
00494       secKE = DBL_MIN;
00495     else
00496       dir = momentum/momentum.mag();
00497 
00498     p = new G4DynamicParticle(vec[i]->GetDefinition(), dir, secKE);
00499     theParticleChange.AddSecondary( p );
00500     delete vec[i];
00501   }
00502 }

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().

00130  {
00131    if( np+nneg+nz == 0 )return;
00132    G4int i;
00133    G4ReactionProduct *p;
00134    for( i=0; i<np; ++i )
00135    {
00136      p = new G4ReactionProduct;
00137      p->SetDefinition( G4PionPlus::PionPlus() );
00138      (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
00139      vec.SetElement( vecLen++, p );
00140    }
00141    for( i=np; i<np+nneg; ++i )
00142    {
00143      p = new G4ReactionProduct;
00144      p->SetDefinition( G4PionMinus::PionMinus() );
00145      (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
00146      vec.SetElement( vecLen++, p );
00147    }
00148    for( i=np+nneg; i<np+nneg+nz; ++i )
00149    {
00150      p = new G4ReactionProduct;
00151      p->SetDefinition( G4PionZero::PionZero() );
00152      (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
00153      vec.SetElement( vecLen++, p );
00154    }
00155  }


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:
Generated on Mon May 27 17:53:18 2013 for Geant4 by  doxygen 1.4.7