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 ) );
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.0
e-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 G4PionMinus * PionMinus()
void Initialize(G4int items)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
static G4Lambda * Lambda()
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4bool ReactionStage(const G4HadProjectile *, G4ReactionProduct &, G4bool &, const G4DynamicParticle *, G4ReactionProduct &, G4bool &, const G4Nucleus &, G4ReactionProduct &, G4FastVector< G4ReactionProduct, 256 > &, G4int &, G4bool, G4ReactionProduct &)
void SetUpChange(G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen, G4ReactionProduct ¤tParticle, G4ReactionProduct &targetParticle, G4bool &incidentHasChanged)
G4bool ReactionStage(const G4HadProjectile *, G4ReactionProduct &, G4bool &, const G4DynamicParticle *, G4ReactionProduct &, G4bool &, const G4Nucleus &, G4ReactionProduct &, G4FastVector< G4ReactionProduct, 256 > &, G4int &, G4bool, G4ReactionProduct &)
static G4KaonZero * KaonZero()
static constexpr double MeV
static G4AntiNeutron * AntiNeutron()
G4double Cinema(G4double kineticEnergy)
void SetMomentumChange(const G4ThreeVector &aV)
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
void SetMomentum(const G4double x, const G4double y, const G4double z)
static G4XiZero * XiZero()
static G4AntiProton * AntiProton()
static G4OmegaMinus * OmegaMinus()
void SetDefinitionAndUpdateE(const G4ParticleDefinition *aParticleDefinition)
const G4String & GetParticleName() const
static G4PionPlus * PionPlus()
void GetNormalizationConstant(const G4double availableEnergy, G4double &n, G4double &anpn)
G4double GetPDGCharge() const
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
const G4ParticleDefinition * GetDefinition() const
static G4SigmaZero * SigmaZero()
static G4Proton * Proton()
static G4XiMinus * XiMinus()
static G4KaonMinus * KaonMinus()
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)
G4double GetPDGMass() const
static G4KaonPlus * KaonPlus()
static G4SigmaMinus * SigmaMinus()
G4double G4Log(G4double x)
void CheckQnums(G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen, G4ReactionProduct ¤tParticle, G4ReactionProduct &targetParticle, G4double Q, G4double B, G4double S)
static G4AntiKaonZero * AntiKaonZero()
static G4SigmaPlus * SigmaPlus()
void SetEnergyChange(G4double anEnergy)
G4double Pmltpc(G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)
static G4KaonZeroLong * KaonZeroLong()
void SetKineticEnergy(const G4double en)
G4double GetKineticEnergy() const
G4ParticleDefinition * GetDefinition() const
G4bool MarkLeadingStrangeParticle(const G4ReactionProduct ¤tParticle, const G4ReactionProduct &targetParticle, G4ReactionProduct &leadParticle)
G4int GetBaryonNumber() const
G4RPGInelastic(const G4String &modelName="RPGInelastic")
G4RPGTwoCluster twoCluster
G4int GetQuarkContent(G4int flavor) const
G4double GetTotalMomentum() const
void SetUpPions(const G4int np, const G4int nm, const G4int nz, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
static constexpr double eV
G4ThreeVector GetMomentum() const
void SetSide(const G4int sid)
Hep3Vector & rotate(double, const Hep3Vector &)
const G4LorentzVector & Get4Momentum() const
G4double G4ParticleHPJENDLHEData::G4double result
void Report(std::ostream &aS)
G4double AnnihilationEvaporationEffects(G4double kineticEnergy, G4double ekOrg)
G4double GetKineticEnergy() const
G4ParticleDefinition * particleDef[18]
void SetElement(G4int anIndex, Type *anElement)
static G4KaonZeroShort * KaonZeroShort()
const G4ParticleDefinition * GetDefinition() const
G4RPGFragmentation fragmentation
std::pair< G4int, G4double > interpolateEnergy(G4double ke) const
static const G4double energyScale[30]
G4int sampleFlat(std::vector< G4double > sigma) const
static G4Neutron * Neutron()
void SetMomentum(const G4ThreeVector &momentum)
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
G4int GetPDGEncoding() const
G4GLOB_DLL std::ostream G4cout
static G4PionZero * PionZero()
G4bool ReactionStage(const G4HadProjectile *, G4ReactionProduct &, G4bool &, const G4DynamicParticle *, G4ReactionProduct &, G4bool &, const G4Nucleus &, G4ReactionProduct &, G4FastVector< G4ReactionProduct, 256 > &, G4int &, G4bool, G4ReactionProduct &)
static constexpr double pi
double B(double temperature)
G4int GetAntiQuarkContent(G4int flavor) const
static constexpr double GeV
G4HadFinalState theParticleChange
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void SetStatusChange(G4HadFinalStateStatus aS)