| Geant4 Cross Reference |
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
38 /** \file G4INCLClusterDecay.cc
39 * \brief Static class for carrying out cluster decays
40 *
41 * \date 6th July 2011
42 * \author Davide Mancusi
43 */
44
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
59 /// \brief Carries out two-body decays
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
136 /// \brief Carries out three-body decays
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 /** \brief Disassembles unbound nuclei using a simple phase-space model
239 *
240 * The decay products are assumed to uniformly populate the momentum space
241 * (the "phase-space" naming is a long-standing but misleading convention).
242 * The generation of the final state is done by rejection, using the
243 * Raubold-Lynch method. Parts of our implementation were "inspired" by
244 * ROOT's TGenPhaseSpace class, which in turn is a translation of CERNLIB's
245 * historical GENBOD routine [CERN report 68-15 (1968)]. The ROOT
246 * implementation is documented at the following URL:
247 *
248 * http://root.cern.ch/root/html/TGenPhaseSpace.html#TGenPhaseSpace
249 */
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
400 /** \brief Disassembles unbound nuclei using the phase-space model
401 *
402 * This implementation uses the Kopylov algorithm, defined in namespace
403 * PhaseSpaceGenerator.
404 */
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
504 /** \brief Recursively decay clusters
505 *
506 * \param c cluster that should decay
507 * \param decayProducts decay products are appended to the end of this list
508 */
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
577 /** \brief Table for cluster decays
578 *
579 * Definition of "Stable": halflife > 1 ms
580 *
581 * These table includes decay data for clusters that INCL presently does
582 * not produce. It can't hurt.
583 *
584 * Unphysical nuclides (A<Z) are marked as stable, but should never be
585 * produced by INCL. If you find them in the output, something is fishy.
586 */
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
661