61 G4cout <<
" **************************************************** " <<
G4endl;
62 G4cout <<
" * The RPG model is currently under development and * " <<
G4endl;
64 G4cout <<
" **************************************************** " <<
G4endl;
77 for( i=2; i<=np; i++ )npf += std::log((
double)i);
78 for( i=2; i<=nneg; i++ )nmf += std::log((
double)i);
79 for( i=2; i<=nz; i++ )nzf += std::log((
double)i);
81 r =
std::min( expxu,
std::max( expxl, -(np-nneg+nz+b)*(np-nneg+nz+b)/(2*c*c*n*n)-npf-nmf-nzf ) );
90 if (j <= 1)
return result;
91 for (
G4int i = 2; i <= j; ++i) result *= i;
113 leadParticle = currentParticle;
120 leadParticle = targetParticle;
131 if( np+nneg+nz == 0 )
return;
134 for( i=0; i<np; ++i )
141 for( i=np; i<np+nneg; ++i )
148 for( i=np+nneg; i<np+nneg+nz; ++i )
165 const G4int numSec = 60;
181 n = 3.62567 + aleab*(0.665843 + aleab*(0.336514 + aleab*(0.117712 + 0.0136912*aleab)));
188 for(
G4int i=iBegin; i<=numSec; ++i )
190 temp =
pi*i/(2.0*n*
n);
194 if( test >= 1.0e-10 )anpn += temp*test;
210 G4bool& incidentHasChanged,
229 incidentHasChanged, originalTarget,
230 targetParticle, targetHasChanged,
231 targetNucleus, currentParticle,
233 false, leadingStrangeParticle);
239 leadingStrangeParticle );
244 G4bool finishedGenXPt =
false;
245 G4bool annihilation =
false;
247 currentParticle.
GetMass() == 0.0 && targetParticle.
GetMass() == 0.0 )
256 if( ek > 1.0*
GeV )ekcor = 1./(ek/
GeV);
258 ek = 2*tarmas + ek*(1.+ekcor/atomicWeight);
273 G4double p = std::sqrt( std::abs(et*et-amas*amas) );
280 if( ekOrg <= 0.0001 )
289 const G4double twsup[] = { 1.0, 0.7, 0.5, 0.3, 0.2, 0.1 };
296 std::vector<G4ReactionProduct> savevec;
297 for (
G4int i = 0; i < vecLen; i++) savevec.push_back(*vec[i]);
306 if( annihilation || vecLen > 5 ||
314 || rand2 > twsup[vecLen]) ) )
318 incidentHasChanged, originalTarget,
319 targetParticle, targetHasChanged,
320 targetNucleus, currentParticle,
322 leadFlag, leadingStrangeParticle);
324 if (finishedGenXPt)
return;
326 G4bool finishedTwoClu =
false;
329 for (
G4int i = 0; i < vecLen; i++)
delete vec[i];
337 if (!finishedGenXPt && annihilation) {
338 currentParticle = saveCurrent;
339 targetParticle = saveTarget;
340 for (
G4int i = 0; i < vecLen; i++)
delete vec[i];
343 for (
G4int i = 0; i <
G4int(savevec.size()); i++) {
363 incidentHasChanged, originalTarget,
364 targetParticle, targetHasChanged,
365 targetNucleus, currentParticle,
367 leadFlag, leadingStrangeParticle);
376 if (finishedTwoClu)
return;
379 incidentHasChanged, originalTarget,
380 targetParticle, targetHasChanged,
381 targetNucleus, currentParticle,
383 false, leadingStrangeParticle);
407 G4bool& incidentHasChanged )
417 incidentHasChanged =
true;
426 incidentHasChanged =
true;
439 for (i = 0; i < vecLen; ++i) {
443 vec[i]->SetDefinitionAndUpdateE(aKaonZL);
445 vec[i]->SetDefinitionAndUpdateE(aKaonZS);
450 if (incidentHasChanged) {
467 if (std::fabs(aE)<.1*
eV) aE=.1*
eV;
471 if (targetParticle.
GetMass() > 0.0)
480 dir = momentum/momentum.mag();
489 for (i = 0; i < vecLen; ++i) {
490 G4double secKE = vec[i]->GetKineticEnergy();
496 dir = momentum/momentum.mag();
505 std::pair<G4int, G4double>
511 for (
G4int i = 1; i < 30; i++) {
518 return std::pair<G4int, G4double>(index, fraction);
527 for (i = 0; i <
G4int(sigma.size()); i++) sum += sigma[i];
533 for (i = 0; i <
G4int(sigma.size()); i++) {
534 partialSum += sigma[i];
535 if (fsum < partialSum) {
561 for (
G4int i = 0; i < vecLen; i++) {
562 secDef = vec[i]->GetDefinition();
570 if (chargeSum != Q) {
574 if (baryonSum != B) {
578 if (strangenessSum != S) {
586 for (
G4int i = 0; i < vecLen; i++) {
587 secDef = vec[i]->GetDefinition();
597 0.0, 0.01, 0.013, 0.018, 0.024, 0.032, 0.042, 0.056, 0.075, 0.1,
598 0.13, 0.18, 0.24, 0.32, 0.42, 0.56, 0.75, 1.0, 1.3, 1.8,
599 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
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
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 &)