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