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