62 G4cout <<
" **************************************************** " <<
G4endl;
63 G4cout <<
" * The RPG model is currently under development and * " <<
G4endl;
65 G4cout <<
" **************************************************** " <<
G4endl;
78 for( i=2; i<=np; i++ )npf +=
G4Log((
double)i);
79 for( i=2; i<=nneg; i++ )nmf +=
G4Log((
double)i);
80 for( i=2; i<=nz; i++ )nzf +=
G4Log((
double)i);
82 r =
std::min( expxu,
std::max( expxl, -(np-nneg+nz+b)*(np-nneg+nz+b)/(2*c*c*n*n)-npf-nmf-nzf ) );
91 if (j <= 1)
return result;
92 for (
G4int i = 2; i <= j; ++i) result *= i;
114 leadParticle = currentParticle;
121 leadParticle = targetParticle;
132 if( np+nneg+nz == 0 )
return;
135 for( i=0; i<np; ++i )
142 for( i=np; i<np+nneg; ++i )
149 for( i=np+nneg; i<np+nneg+nz; ++i )
166 const G4int numSec = 60;
182 n = 3.62567 + aleab*(0.665843 + aleab*(0.336514 + aleab*(0.117712 + 0.0136912*aleab)));
189 for(
G4int i=iBegin; i<=numSec; ++i )
191 temp =
pi*i/(2.0*n*
n);
195 if( test >= 1.0e-10 )anpn += temp*test;
211 G4bool& incidentHasChanged,
230 incidentHasChanged, originalTarget,
231 targetParticle, targetHasChanged,
232 targetNucleus, currentParticle,
234 false, leadingStrangeParticle);
240 leadingStrangeParticle );
245 G4bool finishedGenXPt =
false;
246 G4bool annihilation =
false;
248 currentParticle.
GetMass() == 0.0 && targetParticle.
GetMass() == 0.0 )
257 if( ek > 1.0*
GeV )ekcor = 1./(ek/
GeV);
259 ek = 2*tarmas + ek*(1.+ekcor/atomicWeight);
274 G4double p = std::sqrt( std::abs(et*et-amas*amas) );
281 if( ekOrg <= 0.0001 )
290 const G4double twsup[] = { 1.0, 0.7, 0.5, 0.3, 0.2, 0.1 };
297 std::vector<G4ReactionProduct> savevec;
298 for (
G4int i = 0; i < vecLen; i++) savevec.push_back(*vec[i]);
307 if( annihilation || vecLen > 5 ||
315 || rand2 > twsup[vecLen]) ) )
319 incidentHasChanged, originalTarget,
320 targetParticle, targetHasChanged,
321 targetNucleus, currentParticle,
323 leadFlag, leadingStrangeParticle);
325 if (finishedGenXPt)
return;
327 G4bool finishedTwoClu =
false;
330 for (
G4int i = 0; i < vecLen; i++)
delete vec[i];
338 if (!finishedGenXPt && annihilation) {
339 currentParticle = saveCurrent;
340 targetParticle = saveTarget;
341 for (
G4int i = 0; i < vecLen; i++)
delete vec[i];
344 for (
G4int i = 0; i <
G4int(savevec.size()); i++) {
364 incidentHasChanged, originalTarget,
365 targetParticle, targetHasChanged,
366 targetNucleus, currentParticle,
368 leadFlag, leadingStrangeParticle);
377 if (finishedTwoClu)
return;
380 incidentHasChanged, originalTarget,
381 targetParticle, targetHasChanged,
382 targetNucleus, currentParticle,
384 false, leadingStrangeParticle);
408 G4bool& incidentHasChanged )
418 incidentHasChanged =
true;
427 incidentHasChanged =
true;
440 for (i = 0; i < vecLen; ++i) {
444 vec[i]->SetDefinitionAndUpdateE(aKaonZL);
446 vec[i]->SetDefinitionAndUpdateE(aKaonZS);
451 if (incidentHasChanged) {
468 if (std::fabs(aE)<.1*
eV) aE=.1*
eV;
472 if (targetParticle.
GetMass() > 0.0)
481 dir = momentum/momentum.mag();
490 for (i = 0; i < vecLen; ++i) {
491 G4double secKE = vec[i]->GetKineticEnergy();
497 dir = momentum/momentum.mag();
506 std::pair<G4int, G4double>
512 for (
G4int i = 1; i < 30; i++) {
519 return std::pair<G4int, G4double>(index, fraction);
528 for (i = 0; i <
G4int(sigma.size()); i++) sum += sigma[i];
534 for (i = 0; i <
G4int(sigma.size()); i++) {
535 partialSum += sigma[i];
536 if (fsum < partialSum) {
562 for (
G4int i = 0; i < vecLen; i++) {
563 secDef = vec[i]->GetDefinition();
571 if (chargeSum != Q) {
575 if (baryonSum != B) {
579 if (strangenessSum != S) {
587 for (
G4int i = 0; i < vecLen; i++) {
588 secDef = vec[i]->GetDefinition();
598 0.0, 0.01, 0.013, 0.018, 0.024, 0.032, 0.042, 0.056, 0.075, 0.1,
599 0.13, 0.18, 0.24, 0.32, 0.42, 0.56, 0.75, 1.0, 1.3, 1.8,
600 2.4, 3.2, 4.2, 5.6, 7.5, 10.0, 13.0, 18.0, 24.0, 32.0 };
static G4KaonZero * KaonZero()
void SetElement(G4int anIndex, Type *anElement)
G4bool MarkLeadingStrangeParticle(const G4ReactionProduct ¤tParticle, const G4ReactionProduct &targetParticle, G4ReactionProduct &leadParticle)
void SetMomentum(const G4ThreeVector &momentum)
G4double GetTotalMomentum() const
G4bool ReactionStage(const G4HadProjectile *, G4ReactionProduct &, G4bool &, const G4DynamicParticle *, G4ReactionProduct &, G4bool &, const G4Nucleus &, G4ReactionProduct &, G4FastVector< G4ReactionProduct, 256 > &, G4int &, G4bool, G4ReactionProduct &)
CLHEP::Hep3Vector G4ThreeVector
void SetUpChange(G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen, G4ReactionProduct ¤tParticle, G4ReactionProduct &targetParticle, G4bool &incidentHasChanged)
void SetKineticEnergy(const G4double en)
void SetMomentum(const G4double x, const G4double y, const G4double z)
static G4OmegaMinus * OmegaMinus()
static G4KaonZeroLong * KaonZeroLong()
void SetSide(const G4int sid)
G4int GetAntiQuarkContent(G4int flavor) const
void CalculateMomenta(G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen, const G4HadProjectile *originalIncident, const G4DynamicParticle *originalTarget, G4ReactionProduct &modifiedOriginal, G4Nucleus &targetNucleus, G4ReactionProduct ¤tParticle, G4ReactionProduct &targetParticle, G4bool &incidentHasChanged, G4bool &targetHasChanged, G4bool quasiElastic)
G4int GetPDGEncoding() const
G4int sampleFlat(std::vector< G4double > sigma) const
G4ParticleDefinition * GetDefinition() const
double B(double temperature)
void Initialize(G4int items)
const G4String & GetParticleName() const
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
static G4SigmaZero * SigmaZero()
static G4KaonMinus * KaonMinus()
void SetStatusChange(G4HadFinalStateStatus aS)
G4double Pmltpc(G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)
static G4XiZero * XiZero()
const G4ParticleDefinition * GetDefinition() const
G4GLOB_DLL std::ostream G4cout
const G4ParticleDefinition * GetDefinition() const
G4ParticleDefinition * particleDef[18]
static G4KaonZeroShort * KaonZeroShort()
static const G4double energyScale[30]
static G4AntiProton * AntiProton()
G4double GetKineticEnergy() const
static G4XiMinus * XiMinus()
G4int GetQuarkContent(G4int flavor) const
static G4Proton * Proton()
static G4PionPlus * PionPlus()
void SetDefinitionAndUpdateE(const G4ParticleDefinition *aParticleDefinition)
static G4Neutron * Neutron()
const G4LorentzVector & Get4Momentum() const
static G4PionZero * PionZero()
static G4SigmaMinus * SigmaMinus()
G4double GetKineticEnergy() const
G4double G4Log(G4double x)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
void SetEnergyChange(G4double anEnergy)
G4double AnnihilationEvaporationEffects(G4double kineticEnergy, G4double ekOrg)
static G4AntiKaonZero * AntiKaonZero()
G4double GetPDGMass() const
G4RPGTwoCluster twoCluster
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
G4bool ReactionStage(const G4HadProjectile *, G4ReactionProduct &, G4bool &, const G4DynamicParticle *, G4ReactionProduct &, G4bool &, const G4Nucleus &, G4ReactionProduct &, G4FastVector< G4ReactionProduct, 256 > &, G4int &, G4bool, G4ReactionProduct &)
G4RPGInelastic(const G4String &modelName="RPGInelastic")
G4double Cinema(G4double kineticEnergy)
static G4PionMinus * PionMinus()
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4ThreeVector GetMomentum() const
G4HadFinalState theParticleChange
void CheckQnums(G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen, G4ReactionProduct ¤tParticle, G4ReactionProduct &targetParticle, G4double Q, G4double B, G4double S)
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
static G4SigmaPlus * SigmaPlus()
void GetNormalizationConstant(const G4double availableEnergy, G4double &n, G4double &anpn)
static G4Lambda * Lambda()
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
static G4KaonPlus * KaonPlus()
G4double GetPDGCharge() const
void SetUpPions(const G4int np, const G4int nm, const G4int nz, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
G4RPGFragmentation fragmentation
std::pair< G4int, G4double > interpolateEnergy(G4double ke) const
void SetMomentumChange(const G4ThreeVector &aV)
void Report(std::ostream &aS)
static G4AntiNeutron * AntiNeutron()
G4int GetBaryonNumber() const
G4bool ReactionStage(const G4HadProjectile *, G4ReactionProduct &, G4bool &, const G4DynamicParticle *, G4ReactionProduct &, G4bool &, const G4Nucleus &, G4ReactionProduct &, G4FastVector< G4ReactionProduct, 256 > &, G4int &, G4bool, G4ReactionProduct &)