33 #define INCLXX_IN_GEANT4_MODE 1
54 namespace ClusterDecay {
60 Particle *decayParticle = 0;
61 const ThreeVector mom(0.0, 0.0, 0.0);
62 const ThreeVector pos = c->getPosition();
65 switch(theDecayMode) {
67 decayParticle =
new Particle(
Proton, mom, pos);
70 decayParticle =
new Particle(
Neutron, mom, pos);
73 decayParticle =
new Cluster(2,4,
false);
76 INCL_ERROR(
"Unrecognized cluster-decay mode in two-body decay: " << theDecayMode << std::endl
80 decayParticle->makeParticipant();
81 decayParticle->setNumberOfDecays(1);
82 decayParticle->setPosition(c->getPosition());
83 decayParticle->setEmissionTime(c->getEmissionTime());
84 decayParticle->setRealMass();
88 const ThreeVector velocity = -c->boostVector();
91 const G4int daughterZ = c->getZ() - decayParticle->getZ();
92 const G4int daughterA = c->getA() - decayParticle->getA();
98 c->setMass(daughterMass);
99 c->setExcitationEnergy(0.);
102 const G4double decayMass = decayParticle->getMass();
105 if(motherMass-daughterMass-decayMass>0.)
108 c->setMomentum(momentum);
109 c->adjustEnergyFromMomentum();
110 decayParticle->setMomentum(-momentum);
111 decayParticle->adjustEnergyFromMomentum();
114 decayParticle->boost(velocity);
118 decayProducts->push_back(decayParticle);
123 Particle *decayParticle1 = 0;
124 Particle *decayParticle2 = 0;
125 const ThreeVector mom(0.0, 0.0, 0.0);
126 const ThreeVector pos = c->getPosition();
129 switch(theDecayMode) {
131 decayParticle1 =
new Particle(
Proton, mom, pos);
132 decayParticle2 =
new Particle(
Proton, mom, pos);
135 decayParticle1 =
new Particle(
Neutron, mom, pos);
136 decayParticle2 =
new Particle(
Neutron, mom, pos);
139 INCL_ERROR(
"Unrecognized cluster-decay mode in three-body decay: " << theDecayMode << std::endl
143 decayParticle1->makeParticipant();
144 decayParticle2->makeParticipant();
145 decayParticle1->setNumberOfDecays(1);
146 decayParticle2->setNumberOfDecays(1);
147 decayParticle1->setRealMass();
148 decayParticle2->setRealMass();
151 const G4double motherMass = c->getMass();
152 const ThreeVector velocity = -c->boostVector();
155 const G4int decayZ1 = decayParticle1->getZ();
156 const G4int decayA1 = decayParticle1->getA();
157 const G4int decayZ2 = decayParticle2->getZ();
158 const G4int decayA2 = decayParticle2->getA();
159 const G4int decayZ = decayZ1 + decayZ2;
160 const G4int decayA = decayA1 + decayA2;
161 const G4int daughterZ = c->getZ() - decayZ;
162 const G4int daughterA = c->getA() - decayA;
163 const G4double decayMass1 = decayParticle1->getMass();
164 const G4double decayMass2 = decayParticle2->getMass();
168 G4double qValue = motherMass - daughterMass - decayMass1 - decayMass2;
176 const G4double decayMass = decayMass1 + decayMass2 + qValueB;
183 c->setMass(daughterMass);
184 c->setExcitationEnergy(0.);
189 c->setMomentum(momentumA);
190 c->adjustEnergyFromMomentum();
191 const ThreeVector decayBoostVector = momentumA/std::sqrt(decayMass*decayMass + momentumA.mag2());
198 decayParticle1->setMomentum(momentumB);
199 decayParticle2->setMomentum(-momentumB);
200 decayParticle1->adjustEnergyFromMomentum();
201 decayParticle2->adjustEnergyFromMomentum();
204 decayParticle1->boost(decayBoostVector);
205 decayParticle2->boost(decayBoostVector);
208 decayParticle1->boost(velocity);
209 decayParticle2->boost(velocity);
213 decayProducts->push_back(decayParticle1);
214 decayProducts->push_back(decayParticle2);
217 #ifndef INCLXX_IN_GEANT4_MODE
231 const G4int theA = c->getA();
232 const G4int theZ = c->getZ();
233 const ThreeVector mom(0.0, 0.0, 0.0);
234 const ThreeVector pos = c->getPosition();
238 switch(theDecayMode) {
248 INCL_ERROR(
"Unrecognized cluster-decay mode in phase-space decay: " << theDecayMode << std::endl
255 G4int finalDaughterZ, finalDaughterA;
261 finalDaughterZ -= theZStep;
274 const G4int nSplits = theA-finalDaughterA;
277 G4double availableEnergy = c->getMass() - finalDaughterMass - nSplits*ejectileMass;
279 if(availableEnergy<0.)
284 G4double eMax = finalDaughterMass + availableEnergy;
285 G4double eMin = finalDaughterMass - ejectileMass;
286 for(
G4int iSplit=0; iSplit<nSplits; ++iSplit) {
287 eMax += ejectileMass;
288 eMin += ejectileMass;
294 std::vector<G4double> theCMMomenta;
295 std::vector<G4double> invariantMasses;
310 if(nTries++>maxTries) {
311 INCL_WARN(
"Phase-space decay exceeded the maximum number of rejections (" << maxTries
312 <<
"). Z=" << theZ <<
", A=" << theA <<
", E*=" << c->getExcitationEnergy()
313 <<
", availableEnergy=" << availableEnergy
314 <<
", nSplits=" << nSplits
320 std::vector<G4double> randomNumbers;
321 for(
G4int iSplit=0; iSplit<nSplits-1; ++iSplit)
323 std::sort(randomNumbers.begin(), randomNumbers.end());
326 invariantMasses.clear();
327 invariantMasses.push_back(finalDaughterMass);
328 for(
G4int iSplit=0; iSplit<nSplits-1; ++iSplit)
329 invariantMasses.push_back(finalDaughterMass + (iSplit+1)*ejectileMass + randomNumbers.at(iSplit)*availableEnergy);
330 invariantMasses.push_back(c->getMass());
333 theCMMomenta.clear();
334 for(
G4int iSplit=0; iSplit<nSplits; ++iSplit) {
335 G4double motherMass = invariantMasses.at(nSplits-iSplit);
336 const G4double daughterMass = invariantMasses.at(nSplits-iSplit-1);
339 if(motherMass-daughterMass-ejectileMass>0.)
341 theCMMomenta.push_back(pCM);
346 for(
G4int iSplit=0; iSplit<nSplits; ++iSplit) {
347 ThreeVector
const velocity = -c->boostVector();
349 #if !defined(NDEBUG) && !defined(INCLXX_IN_GEANT4_MODE)
350 const G4double motherMass = c->getMass();
352 c->setA(c->getA() - 1);
353 c->setZ(c->getZ() - theZStep);
354 c->setMass(invariantMasses.at(nSplits-iSplit-1));
356 Particle *ejectile =
new Particle(theEjectileType, mom, pos);
357 ejectile->setRealMass();
360 ThreeVector momentum;
362 c->setMomentum(momentum);
363 c->adjustEnergyFromMomentum();
364 ejectile->setMomentum(-momentum);
365 ejectile->adjustEnergyFromMomentum();
369 ejectile->boost(velocity);
372 decayProducts->push_back(ejectile);
375 c->setExcitationEnergy(0.);
377 #endif // INCLXX_IN_GEANT4_MODE
385 const G4int theA = c->getA();
387 const ThreeVector mom(0.0, 0.0, 0.0);
388 const ThreeVector pos = c->getPosition();
392 switch(theDecayMode) {
402 INCL_ERROR(
"Unrecognized cluster-decay mode in phase-space decay: " << theDecayMode << std::endl
409 G4int finalDaughterZ, finalDaughterA;
413 while(finalDaughterA>0 &&
clusterDecayMode[finalDaughterZ][finalDaughterA]==theDecayMode) {
415 finalDaughterZ -= theZStep;
427 const G4int nSplits = theA-finalDaughterA;
429 const G4double availableEnergy = c->getMass();
432 const ThreeVector boostVector = - c->boostVector();
436 c->setA(finalDaughterA);
437 c->setZ(finalDaughterZ);
439 c->setMomentum(ThreeVector());
440 c->adjustEnergyFromMomentum();
441 products.push_back(c);
442 for(
G4int i=0; i<nSplits; ++i) {
443 Particle *ejectile =
new Particle(theEjectileType, mom, pos);
444 ejectile->setRealMass();
445 products.push_back(ejectile);
451 ParticleList::iterator productsIter = products.begin();
452 std::advance(productsIter, 1);
453 decayProducts->insert(decayProducts->end(), productsIter, products.end());
455 c->setExcitationEnergy(0.);
463 void recursiveDecay(Cluster *
const c,
ParticleList *decayProducts) {
464 const G4int Z = c->getZ();
465 const G4int A = c->getA();
467 if(c->getExcitationEnergy()<0.)
468 c->setExcitationEnergy(0.);
473 switch(theDecayMode) {
475 INCL_ERROR(
"Unrecognized cluster-decay mode: " << theDecayMode << std::endl
487 twoBodyDecay(c, theDecayMode, decayProducts);
492 threeBodyDecay(c, theDecayMode, decayProducts);
497 phaseSpaceDecay(c, theDecayMode, decayProducts);
503 recursiveDecay(c,decayProducts);
508 INCL_DEBUG(
"Cluster is outside the decay-mode table." << c->print() << std::endl);
510 INCL_DEBUG(
"Z==A, will decompose it in free protons." << std::endl);
513 INCL_DEBUG(
"Z==0, will decompose it in free neutrons." << std::endl);
540 {
StableCluster,
StableCluster,
NeutronDecay,
NeutronUnbound,
NeutronUnbound,
NeutronUnbound,
NeutronUnbound,
NeutronUnbound,
NeutronUnbound,
NeutronUnbound,
NeutronUnbound,
NeutronUnbound, NeutronUnbound},
541 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
NeutronDecay,
TwoNeutronDecay,
NeutronDecay,
TwoNeutronDecay,
NeutronUnbound,
NeutronUnbound,
NeutronUnbound,
NeutronUnbound, NeutronUnbound},
542 {
StableCluster,
StableCluster,
ProtonDecay,
StableCluster,
StableCluster,
NeutronDecay,
StableCluster,
NeutronDecay,
StableCluster,
NeutronDecay,
TwoNeutronDecay,
NeutronUnbound, NeutronUnbound},
543 {
StableCluster,
StableCluster,
StableCluster,
ProtonUnbound,
ProtonDecay,
ProtonDecay,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
NeutronDecay,
StableCluster, NeutronDecay},
544 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
ProtonUnbound,
ProtonDecay,
TwoProtonDecay,
StableCluster,
AlphaDecay,
StableCluster,
StableCluster,
StableCluster, StableCluster},
545 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
ProtonUnbound,
TwoProtonDecay,
ProtonDecay,
StableCluster,
ProtonDecay,
StableCluster,
StableCluster, StableCluster},
546 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
ProtonUnbound,
ProtonUnbound,
TwoProtonDecay,
StableCluster,
StableCluster,
StableCluster, StableCluster},
547 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
ProtonUnbound,
ProtonUnbound,
ProtonUnbound,
ProtonDecay,
ProtonDecay, StableCluster},
548 {
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
StableCluster,
ProtonUnbound,
ProtonUnbound,
ProtonUnbound,
ProtonUnbound,
ProtonDecay}
553 recursiveDecay(c, &decayProducts);
565 return decayProducts;
G4int getA() const
Returns the baryon number.
const G4int clusterTableASize
Static class for carrying out cluster decays.
G4ThreadLocal ClusterDecayType clusterDecayMode[ParticleTable::clusterTableZSize][ParticleTable::clusterTableASize]
Table for cluster decays.
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.
ThreeVector normVector(G4double norm=1.)
G4double getRealMass(const G4INCL::ParticleType t)
Get particle mass (in MeV/c^2)
UnorderedVector< Particle * > ParticleList
ParticleList decay(Cluster *const c)
Carries out a cluster decay.
void decay(G4double initialMass, const ThreeVector &theBoostVector, ParticleList &particles)
Generate decay momenta according to a uniform phase-space model.
G4int getZ() const
Returns the charge number.
void setRealMass()
Set the mass of the Particle to its real mass.
void setType(ParticleType t)
const G4int clusterTableZSize