ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLClusterDecay.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLClusterDecay.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35 
36 #include "globals.hh"
37 
45 #include "G4INCLClusterDecay.hh"
46 #include "G4INCLParticleTable.hh"
47 #include "G4INCLKinematicsUtils.hh"
48 #include "G4INCLRandom.hh"
49 #include "G4INCLPhaseSpaceGenerator.hh"
50 // #include <cassert>
51 #include <algorithm>
52 
53 namespace G4INCL {
54 
55  namespace ClusterDecay {
56 
57  namespace {
58 
60  void twoBodyDecay(Cluster * const c, ClusterDecayType theDecayMode, ParticleList *decayProducts) {
61  Particle *decayParticle = 0;
62  const ThreeVector mom(0.0, 0.0, 0.0);
63  const ThreeVector pos = c->getPosition();
64 
65  // Create the emitted particle
66  switch(theDecayMode) {
67  case ProtonDecay:
68  decayParticle = new Particle(Proton, mom, pos);
69  break;
70  case NeutronDecay:
71  decayParticle = new Particle(Neutron, mom, pos);
72  break;
73  case AlphaDecay:
74  decayParticle = new Cluster(2,4,0,false);
75  break;
76  case LambdaDecay:
77  decayParticle = new Particle(Lambda, mom, pos);
78  break;
79  default:
80  INCL_ERROR("Unrecognized cluster-decay mode in two-body decay: " << theDecayMode << '\n'
81  << c->print());
82  return;
83  }
84  decayParticle->makeParticipant();
85  decayParticle->setNumberOfDecays(1);
86  decayParticle->setPosition(c->getPosition());
87  decayParticle->setEmissionTime(c->getEmissionTime());
88  decayParticle->setRealMass();
89 
90  // Save some variables of the mother cluster
91 #ifdef INCLXX_IN_GEANT4_MODE
92  if ((c->getZ() == 1) && (c->getA() == 2) && (c->getS() == -1)) { // no Mass for A=2,Z=1,S=-1 in Geant4
93  c->setMass(2053.952);
94  if (c->getEnergy() < 2053.952) // Energy can be lower than the sum of p and Lambda masses (2053.952)...
95  c->setMomentum(c->getMomentum() * 0.) ;
96  else
97  c->setMomentum(c->getMomentum() / (std::sqrt(c->getMomentum().mag2())/std::sqrt(c->getMomentum().mag2() - 2053.952*2053.952))) ;
98  }
99 #endif
100  G4double motherMass = c->getMass();
101  const ThreeVector velocity = -c->boostVector();
102 
103  // Characteristics of the daughter particle
104  const G4int daughterZ = c->getZ() - decayParticle->getZ();
105  const G4int daughterA = c->getA() - decayParticle->getA();
106  const G4int daughterS = c->getS() - decayParticle->getS();
107  const G4double daughterMass = ParticleTable::getRealMass(daughterA,daughterZ,daughterS);
108 
109  // The mother cluster becomes the daughter
110  c->setZ(daughterZ);
111  c->setA(daughterA);
112  c->setS(daughterS);
113  c->setMass(daughterMass);
114  c->setExcitationEnergy(0.);
115 
116  // Decay kinematics in the mother rest frame
117  const G4double decayMass = decayParticle->getMass();
118 // assert(motherMass-daughterMass-decayMass>-1.e-5); // Q-value should be >0
119  G4double pCM = 0.;
120  if(motherMass-daughterMass-decayMass>0.)
121  pCM = KinematicsUtils::momentumInCM(motherMass, daughterMass, decayMass);
122  const ThreeVector momentum = Random::normVector(pCM);
123  c->setMomentum(momentum);
124  c->adjustEnergyFromMomentum();
125  decayParticle->setMomentum(-momentum);
126  decayParticle->adjustEnergyFromMomentum();
127 
128  // Boost to the lab frame
129  decayParticle->boost(velocity);
130  c->boost(velocity);
131 
132  // Add the decay particle to the list of decay products
133  decayProducts->push_back(decayParticle);
134  }
135 
137  void threeBodyDecay(Cluster * const c, ClusterDecayType theDecayMode, ParticleList *decayProducts) {
138  Particle *decayParticle1 = 0;
139  Particle *decayParticle2 = 0;
140  const ThreeVector mom(0.0, 0.0, 0.0);
141  const ThreeVector pos = c->getPosition();
142 
143  // Create the emitted particles
144  switch(theDecayMode) {
145  case TwoProtonDecay:
146  decayParticle1 = new Particle(Proton, mom, pos);
147  decayParticle2 = new Particle(Proton, mom, pos);
148  break;
149  case TwoNeutronDecay:
150  decayParticle1 = new Particle(Neutron, mom, pos);
151  decayParticle2 = new Particle(Neutron, mom, pos);
152  break;
153  default:
154  INCL_ERROR("Unrecognized cluster-decay mode in three-body decay: " << theDecayMode << '\n'
155  << c->print());
156  return;
157  }
158  decayParticle1->makeParticipant();
159  decayParticle2->makeParticipant();
160  decayParticle1->setNumberOfDecays(1);
161  decayParticle2->setNumberOfDecays(1);
162  decayParticle1->setRealMass();
163  decayParticle2->setRealMass();
164 
165  // Save some variables of the mother cluster
166  const G4double motherMass = c->getMass();
167  const ThreeVector velocity = -c->boostVector();
168 
169  // Masses and charges of the daughter particle and of the decay products
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;
182  const G4double decayMass1 = decayParticle1->getMass();
183  const G4double decayMass2 = decayParticle2->getMass();
184  const G4double daughterMass = ParticleTable::getRealMass(daughterA,daughterZ,daughterS);
185 
186  // Q-values
187  G4double qValue = motherMass - daughterMass - decayMass1 - decayMass2;
188 // assert(qValue > -1e-5); // Q-value should be >0
189  if(qValue<0.)
190  qValue=0.;
191  const G4double qValueB = qValue * Random::shoot();
192 
193  // The decay particles behave as if they had more mass until the second
194  // decay
195  const G4double decayMass = decayMass1 + decayMass2 + qValueB;
196 
197  /* Stage A: mother --> daughter + (decay1+decay2) */
198 
199  // The mother cluster becomes the daughter
200  c->setZ(daughterZ);
201  c->setA(daughterA);
202  c->setS(daughterS);
203  c->setMass(daughterMass);
204  c->setExcitationEnergy(0.);
205 
206  // Decay kinematics in the mother rest frame
207  const G4double pCMA = KinematicsUtils::momentumInCM(motherMass, daughterMass, decayMass);
208  const ThreeVector momentumA = Random::normVector(pCMA);
209  c->setMomentum(momentumA);
210  c->adjustEnergyFromMomentum();
211  const ThreeVector decayBoostVector = momentumA/std::sqrt(decayMass*decayMass + momentumA.mag2());
212 
213  /* Stage B: (decay1+decay2) --> decay1 + decay2 */
214 
215  // Decay kinematics in the (decay1+decay2) rest frame
216  const G4double pCMB = KinematicsUtils::momentumInCM(decayMass, decayMass1, decayMass2);
217  const ThreeVector momentumB = Random::normVector(pCMB);
218  decayParticle1->setMomentum(momentumB);
219  decayParticle2->setMomentum(-momentumB);
220  decayParticle1->adjustEnergyFromMomentum();
221  decayParticle2->adjustEnergyFromMomentum();
222 
223  // Boost decay1 and decay2 to the Stage-A decay frame
224  decayParticle1->boost(decayBoostVector);
225  decayParticle2->boost(decayBoostVector);
226 
227  // Boost all particles to the lab frame
228  decayParticle1->boost(velocity);
229  decayParticle2->boost(velocity);
230  c->boost(velocity);
231 
232  // Add the decay particles to the list of decay products
233  decayProducts->push_back(decayParticle1);
234  decayProducts->push_back(decayParticle2);
235  }
236 
237 #ifdef INCL_DO_NOT_COMPILE
238 
250  void phaseSpaceDecayLegacy(Cluster * const c, ClusterDecayType theDecayMode, ParticleList *decayProducts) {
251  const G4int theA = c->getA();
252  const G4int theZ = c->getZ();
253 // assert(c->getS() == 0);
254  const ThreeVector mom(0.0, 0.0, 0.0);
255  const ThreeVector pos = c->getPosition();
256 
257  G4int theZStep;
258  ParticleType theEjectileType;
259  switch(theDecayMode) {
260  case ProtonUnbound:
261  theZStep = 1;
262  theEjectileType = Proton;
263  break;
264  case NeutronUnbound:
265  theZStep = 0;
266  theEjectileType = Neutron;
267  break;
268  default:
269  INCL_ERROR("Unrecognized cluster-decay mode in phase-space decay: " << theDecayMode << '\n'
270  << c->print());
271  return;
272  }
273 
274  // Find the daughter cluster (first cluster which is not
275  // proton/neutron-unbound, in the sense of the table)
276  G4int finalDaughterZ, finalDaughterA;
277  if(theZ<ParticleTable::clusterTableZSize && theA<ParticleTable::clusterTableASize) {
278  finalDaughterZ=theZ;
279  finalDaughterA=theA;
280  while(clusterDecayMode[0][finalDaughterZ][finalDaughterA]==theDecayMode) { /* Loop checking, 10.07.2015, D.Mancusi */
281  finalDaughterA--;
282  finalDaughterZ -= theZStep;
283  }
284  } else {
285  finalDaughterA = 1;
286  if(theDecayMode==ProtonUnbound)
287  finalDaughterZ = 1;
288  else
289  finalDaughterZ = 0;
290  }
291 // assert(finalDaughterZ<=theZ && finalDaughterA<theA && finalDaughterA>0 && finalDaughterZ>=0);
292  const G4double finalDaughterMass = ParticleTable::getRealMass(finalDaughterA, finalDaughterZ);
293 
294  // Compute the available decay energy
295  const G4int nSplits = theA-finalDaughterA;
296  const G4double ejectileMass = ParticleTable::getRealMass(1, theZStep);
297  // c->getMass() can possibly contain some excitation energy, too
298  G4double availableEnergy = c->getMass() - finalDaughterMass - nSplits*ejectileMass;
299 // assert(availableEnergy>-1.e-5);
300  if(availableEnergy<0.)
301  availableEnergy=0.;
302 
303  // Compute an estimate of the maximum event weight
304  G4double maximumWeight = 1.;
305  G4double eMax = finalDaughterMass + availableEnergy;
306  G4double eMin = finalDaughterMass - ejectileMass;
307  for(G4int iSplit=0; iSplit<nSplits; ++iSplit) {
308  eMax += ejectileMass;
309  eMin += ejectileMass;
310  maximumWeight *= KinematicsUtils::momentumInCM(eMax, eMin, ejectileMass);
311  }
312 
313  // Sample decays until the weight cutoff is satisfied
314  G4double weight;
315  std::vector<G4double> theCMMomenta;
316  std::vector<G4double> invariantMasses;
317  G4int nTries=0;
318  /* Maximum number of trials dependent on nSplits. 50 trials seems to be
319  * sufficient for small nSplits. For nSplits>=5, maximumWeight is a gross
320  * overestimate of the actual maximum weight, leading to unreasonably high
321  * rejection rates. For these cases, we set nSplits=1000, although the sane
322  * thing to do would be to improve the importance sampling (maybe by
323  * parametrising maximumWeight?).
324  */
325  G4int maxTries;
326  if(nSplits<5)
327  maxTries=50;
328  else
329  maxTries=1000;
330  do {
331  if(nTries++>maxTries) {
332  INCL_WARN("Phase-space decay exceeded the maximum number of rejections (" << maxTries
333  << "). Z=" << theZ << ", A=" << theA << ", E*=" << c->getExcitationEnergy()
334  << ", availableEnergy=" << availableEnergy
335  << ", nSplits=" << nSplits
336  << '\n');
337  break;
338  }
339 
340  // Construct a sorted vector of random numbers
341  std::vector<G4double> randomNumbers;
342  for(G4int iSplit=0; iSplit<nSplits-1; ++iSplit)
343  randomNumbers.push_back(Random::shoot0());
344  std::sort(randomNumbers.begin(), randomNumbers.end());
345 
346  // Divide the available decay energy in the right number of steps
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());
352 
353  weight = 1.;
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);
358 // assert(motherMass-daughterMass-ejectileMass>-1.e-5);
359  G4double pCM = 0.;
360  if(motherMass-daughterMass-ejectileMass>0.)
361  pCM = KinematicsUtils::momentumInCM(motherMass, daughterMass, ejectileMass);
362  theCMMomenta.push_back(pCM);
363  weight *= pCM;
364  }
365  } while(maximumWeight*Random::shoot()>weight); /* Loop checking, 10.07.2015, D.Mancusi */
366 
367  for(G4int iSplit=0; iSplit<nSplits; ++iSplit) {
368  ThreeVector const velocity = -c->boostVector();
369 
370 #if !defined(NDEBUG) && !defined(INCLXX_IN_GEANT4_MODE)
371  const G4double motherMass = c->getMass();
372 #endif
373  c->setA(c->getA() - 1);
374  c->setZ(c->getZ() - theZStep);
375  c->setMass(invariantMasses.at(nSplits-iSplit-1));
376 
377  Particle *ejectile = new Particle(theEjectileType, mom, pos);
378  ejectile->setRealMass();
379 
380 // assert(motherMass-c->getMass()-ejectileMass>-1.e-5);
381  ThreeVector momentum;
382  momentum = Random::normVector(theCMMomenta.at(iSplit));
383  c->setMomentum(momentum);
384  c->adjustEnergyFromMomentum();
385  ejectile->setMomentum(-momentum);
386  ejectile->adjustEnergyFromMomentum();
387 
388  // Boost to the lab frame
389  c->boost(velocity);
390  ejectile->boost(velocity);
391 
392  // Add the decay particle to the list of decay products
393  decayProducts->push_back(ejectile);
394  }
395 // assert(std::abs(c->getRealMass()-c->getMass())<1.e-3);
396  c->setExcitationEnergy(0.);
397  }
398 #endif // INCL_DO_NOT_COMPILE
399 
405  void phaseSpaceDecay(Cluster * const c, ClusterDecayType theDecayMode, ParticleList *decayProducts) {
406  const G4int theA = c->getA();
407  const G4int theZ = c->getZ();
408  const G4int theL = (-1)*(c->getS());
409  const ThreeVector mom(0.0, 0.0, 0.0);
410  const ThreeVector pos = c->getPosition();
411 
412  G4int theZStep;
413 
414  ParticleType theEjectileType;
415  switch(theDecayMode) {
416  case ProtonUnbound:
417  theZStep = 1;
418  theEjectileType = Proton;
419  break;
420  case NeutronUnbound:
421  theZStep = 0;
422  theEjectileType = Neutron;
423  break;
424  case LambdaUnbound: // Will always completly decay. Append only if theA == 0 and/or theZ == 0
425  theZStep = -99;
426  if(theZ==0) theEjectileType = Neutron;
427  else theEjectileType = Proton;
428  break;
429  default:
430  INCL_ERROR("Unrecognized cluster-decay mode in phase-space decay: " << theDecayMode << '\n'
431  << c->print());
432  return;
433  }
434 
435  // Find the daughter cluster (first cluster which is not
436  // proton/neutron-unbound, in the sense of the table)
437  G4int finalDaughterZ, finalDaughterA, finalDaughterL;
438  if(theZ<ParticleTable::clusterTableZSize && theA<ParticleTable::clusterTableASize && theZStep != -99) {
439  finalDaughterZ=theZ;
440  finalDaughterA=theA;
441  finalDaughterL=theL;
442  while(finalDaughterA>0 && clusterDecayMode[finalDaughterL][finalDaughterZ][finalDaughterA]!=StableCluster) { /* Loop modified, 15.01.18, J. Hirtz */
443  finalDaughterA--;
444  finalDaughterZ -= theZStep;
445  }
446  } else {
447  finalDaughterA = 1;
448  if(theDecayMode==ProtonUnbound){
449  finalDaughterZ = 1;
450  finalDaughterL = 0;
451  }
452  else if(theDecayMode==NeutronUnbound){
453  finalDaughterZ = 0;
454  finalDaughterL = 0;
455  }
456  else {
457  finalDaughterZ = 0;
458  finalDaughterL = 1;
459  }
460  }
461 // assert(finalDaughterZ<=theZ && finalDaughterA<theA && finalDaughterA>0 && finalDaughterZ>=0 && finalDaughterL>=0);
462 
463  // Compute the available decay energy
464  const G4int nLambda = theL-finalDaughterL;
465  const G4int nSplits = theA-finalDaughterA-nLambda;
466  // c->getMass() can possibly contain some excitation energy, too
467  const G4double availableEnergy = c->getMass();
468 
469  // Save the boost vector for the cluster
470  const ThreeVector boostVector = - c->boostVector();
471 
472  // Create a list of decay products
473  ParticleList products;
474  c->setA(finalDaughterA);
475  c->setZ(finalDaughterZ);
476  c->setS((-1)*finalDaughterL);
477  c->setRealMass();
478  c->setMomentum(ThreeVector());
479  c->adjustEnergyFromMomentum();
480  products.push_back(c);
481 
482  for(G4int j=0; j<nLambda; ++j) {
483  Particle *ejectile = new Particle(Lambda, mom, pos);
484  ejectile->setRealMass();
485  products.push_back(ejectile);
486  }
487  for(G4int i=0; i<nSplits; ++i) {
488  Particle *ejectile = new Particle(theEjectileType, mom, pos);
489  ejectile->setRealMass();
490  products.push_back(ejectile);
491  }
492 
493  PhaseSpaceGenerator::generate(availableEnergy, products);
494  products.boost(boostVector);
495 
496  // Copy decay products in the output list (but skip the residue)
497  ParticleList::iterator productsIter = products.begin();
498  std::advance(productsIter, 1);
499  decayProducts->insert(decayProducts->end(), productsIter, products.end());
500 
501  c->setExcitationEnergy(0.);
502  }
503 
509  void recursiveDecay(Cluster * const c, ParticleList *decayProducts) {
510  const G4int Z = c->getZ();
511  const G4int A = c->getA();
512  const G4int S = c->getS();
513 // assert(c->getExcitationEnergy()>-1.e-5);
514  if(c->getExcitationEnergy()<0.)
515  c->setExcitationEnergy(0.);
516 
517  if(Z<ParticleTable::clusterTableZSize && A<ParticleTable::clusterTableASize && (S*(-1))<ParticleTable::clusterTableSSize) {
518  ClusterDecayType theDecayMode = clusterDecayMode[(S*(-1))][Z][A];
519 
520  switch(theDecayMode) {
521  default:
522  INCL_ERROR("Unrecognized cluster-decay mode: " << theDecayMode << '\n'
523  << c->print());
524  return;
525  break;
526  case StableCluster:
527  // For stable clusters, just return
528  return;
529  break;
530  case ProtonDecay:
531  case NeutronDecay:
532  case LambdaDecay:
533  case AlphaDecay:
534  // Two-body decays
535  twoBodyDecay(c, theDecayMode, decayProducts);
536  break;
537  case TwoProtonDecay:
538  case TwoNeutronDecay:
539  // Three-body decays
540  threeBodyDecay(c, theDecayMode, decayProducts);
541  break;
542  case ProtonUnbound:
543  case NeutronUnbound:
544  case LambdaUnbound:
545  // Phase-space decays
546  phaseSpaceDecay(c, theDecayMode, decayProducts);
547  break;
548  }
549 
550  // Calls itself recursively in case the produced remnant is still unstable.
551  // Sneaky, isn't it.
552  recursiveDecay(c,decayProducts);
553 
554  } else {
555  // The cluster is too large for our decay-mode table. Decompose it only
556  // if Z==0 || Z==A.
557  INCL_DEBUG("Cluster is outside the decay-mode table." << c->print() << '\n');
558  if(Z==A) {
559  INCL_DEBUG("Z==A, will decompose it in free protons." << '\n');
560  phaseSpaceDecay(c, ProtonUnbound, decayProducts);
561  } else if(Z==0) {
562  INCL_DEBUG("Z==0, will decompose it in free neutrons." << '\n');
563  phaseSpaceDecay(c, NeutronUnbound, decayProducts);
564  }
565  }
566  }
567 
568  } // namespace
569 
570  G4bool isStable(Cluster const * const c) {
571  const G4int Z = c->getZ();
572  const G4int A = c->getA();
573  const G4int L = ((-1)*c->getS());
574  return (clusterDecayMode[L][Z][A]==StableCluster);
575  }
576 
587  G4ThreadLocal ClusterDecayType clusterDecayMode[ParticleTable::clusterTableSSize][ParticleTable::clusterTableZSize][ParticleTable::clusterTableASize] =
588  {{/* S = 0 */
589  /* A = 0 1 2 3 4 5 6 7 8 9 10 11 12 */
590  /* Z = 0 */ {StableCluster, StableCluster, NeutronDecay, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound},
591  /* Z = 1 */ {StableCluster, StableCluster, StableCluster, StableCluster, NeutronDecay, TwoNeutronDecay, NeutronDecay, TwoNeutronDecay, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound},
592  /* Z = 2 */ {StableCluster, StableCluster, ProtonDecay, StableCluster, StableCluster, NeutronDecay, StableCluster, NeutronDecay, StableCluster, NeutronDecay, TwoNeutronDecay, NeutronUnbound, NeutronUnbound},
593  /* Z = 3 */ {StableCluster, StableCluster, StableCluster, ProtonUnbound, ProtonDecay, ProtonDecay, StableCluster, StableCluster, StableCluster, StableCluster, NeutronDecay, StableCluster, NeutronDecay},
594  /* Z = 4 */ {StableCluster, StableCluster, StableCluster, StableCluster, ProtonUnbound, ProtonDecay, TwoProtonDecay, StableCluster, AlphaDecay, StableCluster, StableCluster, StableCluster, StableCluster},
595  /* Z = 5 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, ProtonUnbound, TwoProtonDecay, ProtonDecay, StableCluster, ProtonDecay, StableCluster, StableCluster, StableCluster},
596  /* Z = 6 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, ProtonUnbound, ProtonUnbound, TwoProtonDecay, StableCluster, StableCluster, StableCluster, StableCluster},
597  /* Z = 7 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, ProtonUnbound, ProtonUnbound, ProtonUnbound, ProtonDecay, ProtonDecay, StableCluster},
598  /* Z = 8 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, ProtonUnbound, ProtonUnbound, ProtonUnbound, ProtonUnbound, ProtonDecay}
599  },
600  { /* S = -1 */
601  /* A = 0 1 2 3 4 5 6 7 8 9 10 11 12 */
602  /* Z = 0 */ {StableCluster, StableCluster, NeutronDecay, LambdaUnbound, LambdaUnbound, LambdaUnbound, LambdaUnbound, LambdaUnbound, LambdaUnbound, LambdaUnbound, LambdaUnbound, LambdaUnbound, LambdaUnbound},
603  /* Z = 1 */ {StableCluster, StableCluster, LambdaDecay, StableCluster, StableCluster, NeutronDecay, StableCluster, StableCluster, StableCluster, NeutronDecay, NeutronUnbound,NeutronUnbound,NeutronUnbound},
604  /* Z = 2 */ {StableCluster, StableCluster, StableCluster, LambdaUnbound, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, NeutronDecay, StableCluster, NeutronUnbound},
605  /* Z = 3 */ {StableCluster, StableCluster, StableCluster, StableCluster, LambdaUnbound, ProtonDecay, ProtonDecay, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster},
606  /* Z = 4 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, LambdaUnbound, ProtonUnbound, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster},
607  /* Z = 5 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, LambdaUnbound, ProtonUnbound, ProtonDecay, StableCluster, StableCluster, StableCluster, StableCluster},
608  /* Z = 6 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, LambdaUnbound, ProtonUnbound, StableCluster, StableCluster, StableCluster, StableCluster},
609  /* Z = 7 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, LambdaUnbound, ProtonUnbound, ProtonDecay, ProtonDecay, ProtonDecay},
610  /* Z = 8 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, LambdaUnbound, ProtonUnbound, ProtonUnbound, ProtonUnbound}
611  },
612  { /* S = -2 */
613  /* A = 0 1 2 3 4 5 6 7 8 9 10 11 12 */
614  /* Z = 0 */ {StableCluster, StableCluster, LambdaDecay, LambdaUnbound, LambdaUnbound, LambdaUnbound, LambdaUnbound, LambdaUnbound, LambdaUnbound, LambdaUnbound, LambdaUnbound, LambdaUnbound, LambdaUnbound},
615  /* Z = 1 */ {StableCluster, StableCluster, StableCluster, LambdaUnbound, StableCluster, StableCluster, NeutronDecay, StableCluster, StableCluster, StableCluster, NeutronDecay,NeutronUnbound,NeutronUnbound},
616  /* Z = 2 */ {StableCluster, StableCluster, StableCluster, StableCluster, LambdaUnbound, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster},
617  /* Z = 3 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, LambdaUnbound, ProtonDecay, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster},
618  /* Z = 4 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, LambdaUnbound, ProtonUnbound, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster},
619  /* Z = 5 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, LambdaUnbound, ProtonUnbound, StableCluster, StableCluster, StableCluster, StableCluster},
620  /* Z = 6 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, LambdaUnbound, ProtonUnbound, StableCluster, StableCluster, StableCluster},
621  /* Z = 7 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, LambdaUnbound, ProtonUnbound, ProtonDecay, ProtonDecay},
622  /* Z = 8 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, LambdaUnbound, ProtonUnbound, ProtonUnbound}
623  },
624  { /* S = -3 */
625  /* A = 0 1 2 3 4 5 6 7 8 9 10 11 12 */
626  /* Z = 0 */ {StableCluster, StableCluster, StableCluster, LambdaUnbound, LambdaUnbound, LambdaUnbound, LambdaUnbound, LambdaUnbound, LambdaUnbound, LambdaUnbound, LambdaUnbound, LambdaUnbound, LambdaUnbound},
627  /* Z = 1 */ {StableCluster, StableCluster, StableCluster, StableCluster, LambdaUnbound, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster},
628  /* Z = 2 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, LambdaUnbound, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster},
629  /* Z = 3 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, LambdaUnbound, ProtonDecay, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster},
630  /* Z = 4 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, LambdaUnbound, ProtonUnbound, StableCluster, StableCluster, StableCluster, StableCluster},
631  /* Z = 5 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, LambdaUnbound, ProtonUnbound, StableCluster, StableCluster, StableCluster},
632  /* Z = 6 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, LambdaUnbound, ProtonUnbound, StableCluster, StableCluster},
633  /* Z = 7 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, LambdaUnbound, ProtonUnbound, ProtonDecay},
634  /* Z = 8 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, LambdaUnbound, ProtonUnbound}
635  }};
636 
637  ParticleList decay(Cluster * const c) {
638  ParticleList decayProducts;
639  recursiveDecay(c, &decayProducts);
640 
641  for(ParticleIter i = decayProducts.begin(), e =decayProducts.end(); i!=e; i++) (*i)->setBiasCollisionVector(c->getBiasCollisionVector());
642 
643  // Correctly update the particle type
644  if(c->getA()==1) {
645 // assert(c->getZ()==1 || c->getZ()==0);
646  if(c->getZ()==1)
647  c->setType(Proton);
648  else if(c->getS()==-1)
649  c->setType(Lambda);
650  else
651  c->setType(Neutron);
652  c->setRealMass();
653  }
654 
655  return decayProducts;
656  }
657 
658  } // namespace ClusterDecay
659 } // namespace G4INCL
660