34#define INCLXX_IN_GEANT4_MODE 1
55 namespace ClusterDecay {
66 switch(theDecayMode) {
74 decayParticle =
new Cluster(2,4,0,
false);
80 INCL_ERROR(
"Unrecognized cluster-decay mode in two-body decay: " << theDecayMode <<
'\n'
91#ifdef INCLXX_IN_GEANT4_MODE
92 if ((c->
getZ() == 1) && (c->
getA() == 2) && (c->
getS() == -1)) {
120 if(motherMass-daughterMass-decayMass>0.)
129 decayParticle->
boost(velocity);
133 decayProducts->push_back(decayParticle);
144 switch(theDecayMode) {
154 INCL_ERROR(
"Unrecognized cluster-decay mode in three-body decay: " << theDecayMode <<
'\n'
170 const G4int decayZ1 = decayParticle1->
getZ();
171 const G4int decayA1 = decayParticle1->
getA();
172 const G4int decayS1 = decayParticle1->
getS();
173 const G4int decayZ2 = decayParticle2->
getZ();
174 const G4int decayA2 = decayParticle2->
getA();
175 const G4int decayS2 = decayParticle2->
getS();
176 const G4int decayZ = decayZ1 + decayZ2;
177 const G4int decayA = decayA1 + decayA2;
178 const G4int decayS = decayS1 + decayS2;
179 const G4int daughterZ = c->
getZ() - decayZ;
180 const G4int daughterA = c->
getA() - decayA;
181 const G4int daughterS = c->
getS() - decayS;
187 G4double qValue = motherMass - daughterMass - decayMass1 - decayMass2;
195 const G4double decayMass = decayMass1 + decayMass2 + qValueB;
211 const ThreeVector decayBoostVector = momentumA/std::sqrt(decayMass*decayMass + momentumA.
mag2());
224 decayParticle1->
boost(decayBoostVector);
225 decayParticle2->
boost(decayBoostVector);
228 decayParticle1->
boost(velocity);
229 decayParticle2->
boost(velocity);
233 decayProducts->push_back(decayParticle1);
234 decayProducts->push_back(decayParticle2);
237#ifdef INCL_DO_NOT_COMPILE
259 switch(theDecayMode) {
269 INCL_ERROR(
"Unrecognized cluster-decay mode in phase-space decay: " << theDecayMode <<
'\n'
276 G4int finalDaughterZ, finalDaughterA;
282 finalDaughterZ -= theZStep;
295 const G4int nSplits = theA-finalDaughterA;
298 G4double availableEnergy = c->
getMass() - finalDaughterMass - nSplits*ejectileMass;
300 if(availableEnergy<0.)
305 G4double eMax = finalDaughterMass + availableEnergy;
306 G4double eMin = finalDaughterMass - ejectileMass;
307 for(
G4int iSplit=0; iSplit<nSplits; ++iSplit) {
308 eMax += ejectileMass;
309 eMin += ejectileMass;
315 std::vector<G4double> theCMMomenta;
316 std::vector<G4double> invariantMasses;
331 if(nTries++>maxTries) {
332 INCL_WARN(
"Phase-space decay exceeded the maximum number of rejections (" << maxTries
334 <<
", availableEnergy=" << availableEnergy
335 <<
", nSplits=" << nSplits
341 std::vector<G4double> randomNumbers;
342 for(
G4int iSplit=0; iSplit<nSplits-1; ++iSplit)
344 std::sort(randomNumbers.begin(), randomNumbers.end());
347 invariantMasses.clear();
348 invariantMasses.push_back(finalDaughterMass);
349 for(
G4int iSplit=0; iSplit<nSplits-1; ++iSplit)
350 invariantMasses.push_back(finalDaughterMass + (iSplit+1)*ejectileMass + randomNumbers.at(iSplit)*availableEnergy);
351 invariantMasses.push_back(c->
getMass());
354 theCMMomenta.clear();
355 for(
G4int iSplit=0; iSplit<nSplits; ++iSplit) {
356 G4double motherMass = invariantMasses.at(nSplits-iSplit);
357 const G4double daughterMass = invariantMasses.at(nSplits-iSplit-1);
360 if(motherMass-daughterMass-ejectileMass>0.)
362 theCMMomenta.push_back(pCM);
367 for(
G4int iSplit=0; iSplit<nSplits; ++iSplit) {
370#if !defined(NDEBUG) && !defined(INCLXX_IN_GEANT4_MODE)
375 c->
setMass(invariantMasses.at(nSplits-iSplit-1));
377 Particle *ejectile =
new Particle(theEjectileType, mom,
pos);
378 ejectile->setRealMass();
381 ThreeVector momentum;
385 ejectile->setMomentum(-momentum);
386 ejectile->adjustEnergyFromMomentum();
390 ejectile->boost(velocity);
393 decayProducts->push_back(ejectile);
415 switch(theDecayMode) {
426 if(theZ==0) theEjectileType =
Neutron;
427 else theEjectileType =
Proton;
430 INCL_ERROR(
"Unrecognized cluster-decay mode in phase-space decay: " << theDecayMode <<
'\n'
437 G4int finalDaughterZ, finalDaughterA, finalDaughterL;
444 finalDaughterZ -= theZStep;
464 const G4int nLambda = theL-finalDaughterL;
465 const G4int nSplits = theA-finalDaughterA-nLambda;
474 c->
setA(finalDaughterA);
475 c->
setZ(finalDaughterZ);
476 c->
setS((-1)*finalDaughterL);
480 products.push_back(c);
482 for(
G4int j=0; j<nLambda; ++j) {
485 products.push_back(ejectile);
487 for(
G4int i=0; i<nSplits; ++i) {
490 products.push_back(ejectile);
494 products.
boost(boostVector);
497 ParticleList::iterator productsIter = products.begin();
498 std::advance(productsIter, 1);
499 decayProducts->insert(decayProducts->end(), productsIter, products.end());
520 switch(theDecayMode) {
522 INCL_ERROR(
"Unrecognized cluster-decay mode: " << theDecayMode <<
'\n'
557 INCL_DEBUG(
"Cluster is outside the decay-mode table." << c->
print() <<
'\n');
559 INCL_DEBUG(
"Z==A, will decompose it in free protons." <<
'\n');
562 INCL_DEBUG(
"Z==0, will decompose it in free neutrons." <<
'\n');
590 {
StableCluster,
StableCluster,
NeutronDecay,
NeutronUnbound,
NeutronUnbound,
NeutronUnbound,
NeutronUnbound,
NeutronUnbound,
NeutronUnbound,
NeutronUnbound,
NeutronUnbound,
NeutronUnbound,
NeutronUnbound},
591 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
NeutronDecay,
TwoNeutronDecay,
NeutronDecay,
TwoNeutronDecay,
NeutronUnbound,
NeutronUnbound,
NeutronUnbound,
NeutronUnbound,
NeutronUnbound},
592 {
StableCluster,
StableCluster,
ProtonDecay,
StableCluster,
StableCluster,
NeutronDecay,
StableCluster,
NeutronDecay,
StableCluster,
NeutronDecay,
TwoNeutronDecay,
NeutronUnbound,
NeutronUnbound},
593 {
StableCluster,
StableCluster,
StableCluster,
ProtonUnbound,
ProtonDecay,
ProtonDecay,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
NeutronDecay,
StableCluster,
NeutronDecay},
594 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
ProtonUnbound,
ProtonDecay,
TwoProtonDecay,
StableCluster,
AlphaDecay,
StableCluster,
StableCluster,
StableCluster,
StableCluster},
595 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
ProtonUnbound,
TwoProtonDecay,
ProtonDecay,
StableCluster,
ProtonDecay,
StableCluster,
StableCluster,
StableCluster},
596 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
ProtonUnbound,
ProtonUnbound,
TwoProtonDecay,
StableCluster,
StableCluster,
StableCluster,
StableCluster},
597 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
ProtonUnbound,
ProtonUnbound,
ProtonUnbound,
ProtonDecay,
ProtonDecay,
StableCluster},
598 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
ProtonUnbound,
ProtonUnbound,
ProtonUnbound,
ProtonUnbound,
ProtonDecay}
602 {
StableCluster,
StableCluster,
NeutronDecay,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound},
603 {
StableCluster,
StableCluster,
LambdaDecay,
StableCluster,
StableCluster,
NeutronDecay,
StableCluster,
StableCluster,
StableCluster,
NeutronDecay,
NeutronUnbound,
NeutronUnbound,
NeutronUnbound},
604 {
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
NeutronDecay,
StableCluster,
NeutronUnbound},
605 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
ProtonDecay,
ProtonDecay,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster},
606 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
ProtonUnbound,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster},
607 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
ProtonUnbound,
ProtonDecay,
StableCluster,
StableCluster,
StableCluster,
StableCluster},
608 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
ProtonUnbound,
StableCluster,
StableCluster,
StableCluster,
StableCluster},
609 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
ProtonUnbound,
ProtonDecay,
ProtonDecay,
ProtonDecay},
610 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
ProtonUnbound,
ProtonUnbound,
ProtonUnbound}
614 {
StableCluster,
StableCluster,
LambdaDecay,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound},
615 {
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
StableCluster,
StableCluster,
NeutronDecay,
StableCluster,
StableCluster,
StableCluster,
NeutronDecay,
NeutronUnbound,
NeutronUnbound},
616 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster},
617 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
ProtonDecay,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster},
618 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
ProtonUnbound,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster},
619 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
ProtonUnbound,
StableCluster,
StableCluster,
StableCluster,
StableCluster},
620 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
ProtonUnbound,
StableCluster,
StableCluster,
StableCluster},
621 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
ProtonUnbound,
ProtonDecay,
ProtonDecay},
622 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
ProtonUnbound,
ProtonUnbound}
626 {
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound,
LambdaUnbound},
627 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster},
628 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster},
629 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
ProtonDecay,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster},
630 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
ProtonUnbound,
StableCluster,
StableCluster,
StableCluster,
StableCluster},
631 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
ProtonUnbound,
StableCluster,
StableCluster,
StableCluster},
632 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
ProtonUnbound,
StableCluster,
StableCluster},
633 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
ProtonUnbound,
ProtonDecay},
634 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
LambdaUnbound,
ProtonUnbound}
648 else if(c->
getS()==-1)
655 return decayProducts;
G4double S(G4double temp)
static const G4double pos
Static class for carrying out cluster decays.
static constexpr double L
void boost(const ThreeVector &aBoostVector)
Boost the cluster with the indicated velocity.
void setExcitationEnergy(const G4double e)
Set the excitation energy of the cluster.
void setS(const G4int S)
Set the strangess number of the cluster.
G4double getExcitationEnergy() const
Get the excitation energy of the cluster.
void setA(const G4int A)
Set the mass number of the cluster.
void setZ(const G4int Z)
Set the charge number of the cluster.
std::string print() const
void boost(const ThreeVector &b) const
ThreeVector boostVector() const
G4int getS() const
Returns the strangeness number.
G4double getEnergy() const
void setMass(G4double mass)
G4int getZ() const
Returns the charge number.
const G4INCL::ThreeVector & getPosition() const
std::vector< G4int > getBiasCollisionVector() const
Get the vector list of biased vertices on the particle path.
void setEmissionTime(G4double t)
G4double adjustEnergyFromMomentum()
Recompute the energy to match the momentum.
void setNumberOfDecays(G4int n)
Set the number of decays undergone by the particle.
const G4INCL::ThreeVector & getMomentum() const
virtual void setMomentum(const G4INCL::ThreeVector &momentum)
G4double getEmissionTime()
virtual void makeParticipant()
void setType(ParticleType t)
G4double getMass() const
Get the cached particle mass.
void boost(const ThreeVector &aBoostVector)
virtual void setPosition(const G4INCL::ThreeVector &position)
void setRealMass()
Set the mass of the Particle to its real mass.
G4int getA() const
Returns the baryon number.
void threeBodyDecay(Cluster *const c, ClusterDecayType theDecayMode, ParticleList *decayProducts)
Carries out three-body decays.
void phaseSpaceDecay(Cluster *const c, ClusterDecayType theDecayMode, ParticleList *decayProducts)
Disassembles unbound nuclei using the phase-space model.
void twoBodyDecay(Cluster *const c, ClusterDecayType theDecayMode, ParticleList *decayProducts)
Carries out two-body decays.
void recursiveDecay(Cluster *const c, ParticleList *decayProducts)
Recursively decay clusters.
G4ThreadLocal ClusterDecayType clusterDecayMode[ParticleTable::clusterTableSSize][ParticleTable::clusterTableZSize][ParticleTable::clusterTableASize]
Table for cluster decays.
ParticleList decay(Cluster *const c)
Carries out a cluster decay.
G4bool isStable(Cluster const *const c)
True if the cluster is stable.
G4double momentumInCM(Particle const *const p1, Particle const *const p2)
gives the momentum in the CM frame of two particles.
const G4int clusterTableZSize
G4double getRealMass(const G4INCL::ParticleType t)
Get particle mass (in MeV/c^2)
const G4int clusterTableSSize
const G4int clusterTableASize
void generate(const G4double sqrtS, ParticleList &particles)
Generate an event in the CM system.
ThreeVector normVector(G4double norm=1.)
ParticleList::const_iterator ParticleIter