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