| 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 /* 38 /*
39 * G4INCLNucleus.cc 39 * G4INCLNucleus.cc
40 * 40 *
41 * \date Jun 5, 2009 41 * \date Jun 5, 2009
42 * \author Pekka Kaitaniemi 42 * \author Pekka Kaitaniemi
43 */ 43 */
44 44
45 #ifndef G4INCLNucleus_hh 45 #ifndef G4INCLNucleus_hh
46 #define G4INCLNucleus_hh 1 46 #define G4INCLNucleus_hh 1
47 47
48 #include "G4INCLGlobals.hh" 48 #include "G4INCLGlobals.hh"
49 #include "G4INCLLogger.hh" 49 #include "G4INCLLogger.hh"
50 #include "G4INCLParticle.hh" 50 #include "G4INCLParticle.hh"
51 #include "G4INCLIAvatar.hh" 51 #include "G4INCLIAvatar.hh"
52 #include "G4INCLNucleus.hh" 52 #include "G4INCLNucleus.hh"
53 #include "G4INCLKinematicsUtils.hh" 53 #include "G4INCLKinematicsUtils.hh"
54 #include "G4INCLDecayAvatar.hh" 54 #include "G4INCLDecayAvatar.hh"
55 #include "G4INCLStrangeAbsorbtionChannel.hh" 55 #include "G4INCLStrangeAbsorbtionChannel.hh"
56 #include "G4INCLCluster.hh" 56 #include "G4INCLCluster.hh"
57 #include "G4INCLClusterDecay.hh" 57 #include "G4INCLClusterDecay.hh"
58 #include "G4INCLDeJongSpin.hh" 58 #include "G4INCLDeJongSpin.hh"
59 #include "G4INCLParticleSpecies.hh" 59 #include "G4INCLParticleSpecies.hh"
60 #include "G4INCLParticleTable.hh" 60 #include "G4INCLParticleTable.hh"
61 #include <iterator> 61 #include <iterator>
62 #include <cstdlib> 62 #include <cstdlib>
63 #include <sstream> 63 #include <sstream>
64 // #include <cassert> 64 // #include <cassert>
65 #include "G4INCLPbarAtrestEntryChannel.hh" 65 #include "G4INCLPbarAtrestEntryChannel.hh"
66 #include "G4INCLBinaryCollisionAvatar.hh" 66 #include "G4INCLBinaryCollisionAvatar.hh"
67 67
68 namespace G4INCL { 68 namespace G4INCL {
69 69
70 Nucleus::Nucleus(G4int mass, G4int charge, G 70 Nucleus::Nucleus(G4int mass, G4int charge, G4int strangess, Config const * const conf, const G4double universeRadius,
71 AnnihilationType AType) //D 71 AnnihilationType AType) //D
72 : Cluster(charge,mass,strangess,true), 72 : Cluster(charge,mass,strangess,true),
73 theInitialZ(charge), theInitialA(mass), t 73 theInitialZ(charge), theInitialA(mass), theInitialS(strangess),
74 theNpInitial(0), theNnInitial(0), 74 theNpInitial(0), theNnInitial(0),
75 theNpionplusInitial(0), theNpionminusInit 75 theNpionplusInitial(0), theNpionminusInitial(0),
76 theNkaonplusInitial(0), theNkaonminusInit 76 theNkaonplusInitial(0), theNkaonminusInitial(0),
77 theNantiprotonInitial(0), 77 theNantiprotonInitial(0),
78 initialInternalEnergy(0.), 78 initialInternalEnergy(0.),
79 incomingAngularMomentum(0.,0.,0.), incomi 79 incomingAngularMomentum(0.,0.,0.), incomingMomentum(0.,0.,0.),
80 initialCenterOfMass(0.,0.,0.), 80 initialCenterOfMass(0.,0.,0.),
81 remnant(true), 81 remnant(true),
82 initialEnergy(0.), 82 initialEnergy(0.),
83 tryCN(false), 83 tryCN(false),
84 theUniverseRadius(universeRadius), 84 theUniverseRadius(universeRadius),
85 isNucleusNucleus(false), 85 isNucleusNucleus(false),
86 theProjectileRemnant(NULL), 86 theProjectileRemnant(NULL),
87 theDensity(NULL), 87 theDensity(NULL),
88 thePotential(NULL), 88 thePotential(NULL),
89 theAType(AType) 89 theAType(AType)
90 { 90 {
91 PotentialType potentialType; 91 PotentialType potentialType;
92 G4bool pionPotential; 92 G4bool pionPotential;
93 if(conf) { 93 if(conf) {
94 potentialType = conf->getPotentialType() 94 potentialType = conf->getPotentialType();
95 pionPotential = conf->getPionPotential() 95 pionPotential = conf->getPionPotential();
96 } else { // By default we don't use energy 96 } else { // By default we don't use energy dependent
97 // potential. This is convenient for som 97 // potential. This is convenient for some tests.
98 potentialType = IsospinPotential; 98 potentialType = IsospinPotential;
99 pionPotential = true; 99 pionPotential = true;
100 } 100 }
101 101
102 thePotential = NuclearPotential::createPot 102 thePotential = NuclearPotential::createPotential(potentialType, theA, theZ, pionPotential);
103 103
104 ParticleTable::setProtonSeparationEnergy(t 104 ParticleTable::setProtonSeparationEnergy(thePotential->getSeparationEnergy(Proton));
105 ParticleTable::setNeutronSeparationEnergy( 105 ParticleTable::setNeutronSeparationEnergy(thePotential->getSeparationEnergy(Neutron));
106 106
107 if (theAType==PType) theDensity = NuclearD 107 if (theAType==PType) theDensity = NuclearDensityFactory::createDensity(theA+1, theZ+1, theS);
108 else if (theAType==NType) theDensity = Nuc 108 else if (theAType==NType) theDensity = NuclearDensityFactory::createDensity(theA+1, theZ, theS);
109 else 109 else
110 theDensity = NuclearDensityFactory::create 110 theDensity = NuclearDensityFactory::createDensity(theA, theZ, theS);
111 111
112 theParticleSampler->setPotential(thePotent 112 theParticleSampler->setPotential(thePotential);
113 theParticleSampler->setDensity(theDensity) 113 theParticleSampler->setDensity(theDensity);
114 114
115 if(theUniverseRadius<0) 115 if(theUniverseRadius<0)
116 theUniverseRadius = theDensity->getMaxim 116 theUniverseRadius = theDensity->getMaximumRadius();
117 theStore = new Store(conf); 117 theStore = new Store(conf);
118 } 118 }
119 119
120 Nucleus::~Nucleus() { 120 Nucleus::~Nucleus() {
121 delete theStore; 121 delete theStore;
122 deleteProjectileRemnant(); 122 deleteProjectileRemnant();
123 /* We don't delete the potential and the d 123 /* We don't delete the potential and the density here any more -- Factories
124 * are caching them 124 * are caching them
125 delete thePotential; 125 delete thePotential;
126 delete theDensity;*/ 126 delete theDensity;*/
127 } 127 }
128 128
129 AnnihilationType Nucleus::getAType() const { 129 AnnihilationType Nucleus::getAType() const {
130 return theAType; 130 return theAType;
131 } 131 }
132 132
133 void Nucleus::setAType(AnnihilationType type 133 void Nucleus::setAType(AnnihilationType type) {
134 theAType = type; 134 theAType = type;
135 } 135 }
136 136
137 void Nucleus::initializeParticles() { 137 void Nucleus::initializeParticles() {
138 // Reset the variables connected with the 138 // Reset the variables connected with the projectile remnant
139 delete theProjectileRemnant; 139 delete theProjectileRemnant;
140 theProjectileRemnant = NULL; 140 theProjectileRemnant = NULL;
141 Cluster::initializeParticles(); 141 Cluster::initializeParticles();
142 142
143 for(ParticleIter i=particles.begin(), e=pa 143 for(ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) {
144 updatePotentialEnergy(*i); 144 updatePotentialEnergy(*i);
145 } 145 }
146 theStore->add(particles); 146 theStore->add(particles);
147 particles.clear(); 147 particles.clear();
148 initialInternalEnergy = computeTotalEnergy 148 initialInternalEnergy = computeTotalEnergy();
149 initialCenterOfMass = thePosition; 149 initialCenterOfMass = thePosition;
150 } 150 }
151 151
152 void Nucleus::applyFinalState(FinalState *fi 152 void Nucleus::applyFinalState(FinalState *finalstate) {
153 if(!finalstate) // do nothing if no final 153 if(!finalstate) // do nothing if no final state was returned
154 return; 154 return;
155 155
156 G4double totalEnergy = 0.0; 156 G4double totalEnergy = 0.0;
157 157
158 FinalStateValidity const validity = finals 158 FinalStateValidity const validity = finalstate->getValidity();
159 if(validity == ValidFS) { 159 if(validity == ValidFS) {
160 160
161 ParticleList const &created = finalstate 161 ParticleList const &created = finalstate->getCreatedParticles();
162 for(ParticleIter iter=created.begin(), e 162 for(ParticleIter iter=created.begin(), e=created.end(); iter!=e; ++iter) {
163 theStore->add((*iter)); 163 theStore->add((*iter));
164 if(!(*iter)->isOutOfWell()) { 164 if(!(*iter)->isOutOfWell()) {
165 totalEnergy += (*iter)->getEnergy() 165 totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
166 } 166 }
167 } 167 }
168 168
169 ParticleList const &deleted = finalstate 169 ParticleList const &deleted = finalstate->getDestroyedParticles();
170 for(ParticleIter iter=deleted.begin(), e 170 for(ParticleIter iter=deleted.begin(), e=deleted.end(); iter!=e; ++iter) {
171 theStore->particleHasBeenDestroyed(*it 171 theStore->particleHasBeenDestroyed(*iter);
172 } 172 }
173 173
174 ParticleList const &modified = finalstat 174 ParticleList const &modified = finalstate->getModifiedParticles();
175 for(ParticleIter iter=modified.begin(), 175 for(ParticleIter iter=modified.begin(), e=modified.end(); iter!=e; ++iter) {
176 theStore->particleHasBeenUpdated(*iter 176 theStore->particleHasBeenUpdated(*iter);
177 totalEnergy += (*iter)->getEnergy() - 177 totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
178 } 178 }
179 179
180 ParticleList const &out = finalstate->ge 180 ParticleList const &out = finalstate->getOutgoingParticles();
181 for(ParticleIter iter=out.begin(), e=out 181 for(ParticleIter iter=out.begin(), e=out.end(); iter!=e; ++iter) {
182 if((*iter)->isCluster()) { 182 if((*iter)->isCluster()) {
183 Cluster *clusterOut = dynamic_cast<Clust 183 Cluster *clusterOut = dynamic_cast<Cluster*>((*iter));
184 // assert(clusterOut); 184 // assert(clusterOut);
185 #ifdef INCLXX_IN_GEANT4_MODE 185 #ifdef INCLXX_IN_GEANT4_MODE
186 if(!clusterOut) 186 if(!clusterOut)
187 continue; 187 continue;
188 #endif 188 #endif
189 ParticleList const &components = clu 189 ParticleList const &components = clusterOut->getParticles();
190 for(ParticleIter in=components.begin 190 for(ParticleIter in=components.begin(), end=components.end(); in!=end; ++in)
191 theStore->particleHasBeenEjected(* 191 theStore->particleHasBeenEjected(*in);
192 } else { 192 } else {
193 theStore->particleHasBeenEjected(*it 193 theStore->particleHasBeenEjected(*iter);
194 } 194 }
195 totalEnergy += (*iter)->getEnergy(); 195 totalEnergy += (*iter)->getEnergy(); // No potential here because the particle is gone
196 theA -= (*iter)->getA(); 196 theA -= (*iter)->getA();
197 theZ -= (*iter)->getZ(); 197 theZ -= (*iter)->getZ();
198 theS -= (*iter)->getS(); 198 theS -= (*iter)->getS();
199 theStore->addToOutgoing(*iter); 199 theStore->addToOutgoing(*iter);
200 (*iter)->setEmissionTime(theStore->get 200 (*iter)->setEmissionTime(theStore->getBook().getCurrentTime());
201 } 201 }
202 202
203 ParticleList const &entering = finalstat 203 ParticleList const &entering = finalstate->getEnteringParticles();
204 for(ParticleIter iter=entering.begin(), 204 for(ParticleIter iter=entering.begin(), e=entering.end(); iter!=e; ++iter) {
205 insertParticle(*iter); 205 insertParticle(*iter);
206 totalEnergy += (*iter)->getEnergy() - 206 totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
207 } 207 }
208 208
209 // actually perform the removal of the s 209 // actually perform the removal of the scheduled avatars
210 theStore->removeScheduledAvatars(); 210 theStore->removeScheduledAvatars();
211 } else if(validity == ParticleBelowFermiFS 211 } else if(validity == ParticleBelowFermiFS || validity == ParticleBelowZeroFS) {
212 INCL_DEBUG("A Particle is entering below 212 INCL_DEBUG("A Particle is entering below the Fermi sea:" << '\n' << finalstate->print() << '\n');
213 tryCN = true; 213 tryCN = true;
214 ParticleList const &entering = finalstat 214 ParticleList const &entering = finalstate->getEnteringParticles();
215 for(ParticleIter iter=entering.begin(), 215 for(ParticleIter iter=entering.begin(), e=entering.end(); iter!=e; ++iter) {
216 insertParticle(*iter); 216 insertParticle(*iter);
217 } 217 }
218 } 218 }
219 219
220 if(validity==ValidFS && 220 if(validity==ValidFS &&
221 std::abs(totalEnergy - finalstate->get 221 std::abs(totalEnergy - finalstate->getTotalEnergyBeforeInteraction()) > 0.1) {
222 INCL_ERROR("Energy nonconservation! Ener 222 INCL_ERROR("Energy nonconservation! Energy at the beginning of the event = "
223 << finalstate->getTotalEnergyBeforeI 223 << finalstate->getTotalEnergyBeforeInteraction()
224 <<" and after interaction = " 224 <<" and after interaction = "
225 << totalEnergy << '\n' 225 << totalEnergy << '\n'
226 << finalstate->print()); 226 << finalstate->print());
227 } 227 }
228 } 228 }
229 229
230 void Nucleus::propagateParticles(G4double /* 230 void Nucleus::propagateParticles(G4double /*step*/) {
231 INCL_WARN("Useless Nucleus::propagateParti 231 INCL_WARN("Useless Nucleus::propagateParticles -method called." << '\n');
232 } 232 }
233 233
234 G4double Nucleus::computeTotalEnergy() const 234 G4double Nucleus::computeTotalEnergy() const {
235 G4double totalEnergy = 0.0; 235 G4double totalEnergy = 0.0;
236 ParticleList const &inside = theStore->get 236 ParticleList const &inside = theStore->getParticles();
237 for(ParticleIter p=inside.begin(), e=insid 237 for(ParticleIter p=inside.begin(), e=inside.end(); p!=e; ++p) {
238 if((*p)->isNucleon()) // Ugly: we should 238 if((*p)->isNucleon()) // Ugly: we should calculate everything using total energies!
239 totalEnergy += (*p)->getKineticEnergy( 239 totalEnergy += (*p)->getKineticEnergy() - (*p)->getPotentialEnergy();
240 else if((*p)->isResonance()) 240 else if((*p)->isResonance())
241 totalEnergy += (*p)->getEnergy() - (*p 241 totalEnergy += (*p)->getEnergy() - (*p)->getPotentialEnergy() - ParticleTable::effectiveNucleonMass;
242 else if((*p)->isHyperon()) 242 else if((*p)->isHyperon())
243 totalEnergy += (*p)->getEnergy() - (*p 243 totalEnergy += (*p)->getEnergy() - (*p)->getPotentialEnergy() - ParticleTable::getRealMass((*p)->getType());
244 else if((*p)->isAntiNucleon()) 244 else if((*p)->isAntiNucleon())
245 totalEnergy += (*p)->getEnergy() - (*p 245 totalEnergy += (*p)->getEnergy() - (*p)->getPotentialEnergy() + ParticleTable::getINCLMass(Proton) - ParticleTable::getProtonSeparationEnergy();
246 else if((*p)->isAntiLambda()) 246 else if((*p)->isAntiLambda())
247 totalEnergy += (*p)->getEnergy() - (*p 247 totalEnergy += (*p)->getEnergy() - (*p)->getPotentialEnergy() + ParticleTable::getRealMass((*p)->getType()) - ParticleTable::getSeparationEnergyINCL(Lambda, theA, theZ);
248 //std::cout << ParticleTable::getRealM 248 //std::cout << ParticleTable::getRealMass((*p)->getType()) << std::endl;}
249 else 249 else
250 totalEnergy += (*p)->getEnergy() - (*p 250 totalEnergy += (*p)->getEnergy() - (*p)->getPotentialEnergy();
251 } 251 }
252 252
253 return totalEnergy; 253 return totalEnergy;
254 } 254 }
255 255
256 void Nucleus::computeRecoilKinematics() { 256 void Nucleus::computeRecoilKinematics() {
257 // If the remnant consists of only one nuc 257 // If the remnant consists of only one nucleon, we need to apply a special
258 // procedure to put it on mass shell. 258 // procedure to put it on mass shell.
259 if(theA==1) { 259 if(theA==1) {
260 emitInsidePions(); 260 emitInsidePions();
261 computeOneNucleonRecoilKinematics(); 261 computeOneNucleonRecoilKinematics();
262 remnant=false; 262 remnant=false;
263 return; 263 return;
264 } 264 }
265 265
266 // Compute the recoil momentum and angular 266 // Compute the recoil momentum and angular momentum
267 theMomentum = incomingMomentum; 267 theMomentum = incomingMomentum;
268 theSpin = incomingAngularMomentum; 268 theSpin = incomingAngularMomentum;
269 269
270 ParticleList const &outgoing = theStore->g 270 ParticleList const &outgoing = theStore->getOutgoingParticles();
271 for(ParticleIter p=outgoing.begin(), e=out 271 for(ParticleIter p=outgoing.begin(), e=outgoing.end(); p!=e; ++p) {
272 theMomentum -= (*p)->getMomentum(); 272 theMomentum -= (*p)->getMomentum();
273 theSpin -= (*p)->getAngularMomentum(); 273 theSpin -= (*p)->getAngularMomentum();
274 } 274 }
275 if(theProjectileRemnant) { 275 if(theProjectileRemnant) {
276 theMomentum -= theProjectileRemnant->get 276 theMomentum -= theProjectileRemnant->getMomentum();
277 theSpin -= theProjectileRemnant->getAngu 277 theSpin -= theProjectileRemnant->getAngularMomentum();
278 } 278 }
279 279
280 // Subtract orbital angular momentum 280 // Subtract orbital angular momentum
281 thePosition = computeCenterOfMass(); 281 thePosition = computeCenterOfMass();
282 theSpin -= (thePosition-initialCenterOfMas 282 theSpin -= (thePosition-initialCenterOfMass).vector(theMomentum);
283 283
284 setMass(ParticleTable::getTableMass(theA,t 284 setMass(ParticleTable::getTableMass(theA,theZ,theS) + theExcitationEnergy);
285 adjustEnergyFromMomentum(); 285 adjustEnergyFromMomentum();
286 remnant=true; 286 remnant=true;
287 } 287 }
288 288
289 ThreeVector Nucleus::computeCenterOfMass() c 289 ThreeVector Nucleus::computeCenterOfMass() const {
290 ThreeVector cm(0.,0.,0.); 290 ThreeVector cm(0.,0.,0.);
291 G4double totalMass = 0.0; 291 G4double totalMass = 0.0;
292 ParticleList const &inside = theStore->get 292 ParticleList const &inside = theStore->getParticles();
293 for(ParticleIter p=inside.begin(), e=insid 293 for(ParticleIter p=inside.begin(), e=inside.end(); p!=e; ++p) {
294 const G4double mass = (*p)->getMass(); 294 const G4double mass = (*p)->getMass();
295 cm += (*p)->getPosition() * mass; 295 cm += (*p)->getPosition() * mass;
296 totalMass += mass; 296 totalMass += mass;
297 } 297 }
298 cm /= totalMass; 298 cm /= totalMass;
299 return cm; 299 return cm;
300 } 300 }
301 301
302 G4double Nucleus::computeExcitationEnergy() 302 G4double Nucleus::computeExcitationEnergy() const {
303 const G4double totalEnergy = computeTotalE 303 const G4double totalEnergy = computeTotalEnergy();
304 const G4double separationEnergies = comput 304 const G4double separationEnergies = computeSeparationEnergyBalance();
305 305
306 G4double eSep = 0; 306 G4double eSep = 0;
307 if (getAType() == AnnihilationType::Def) { 307 if (getAType() == AnnihilationType::Def) {
308 } else if (getAType() == AnnihilationType::PT 308 } else if (getAType() == AnnihilationType::PType) {
309 } else if (getAType() == AnnihilationType::NT 309 } else if (getAType() == AnnihilationType::NType) {
310 } else if (getAType() == AnnihilationType::PT 310 } else if (getAType() == AnnihilationType::PTypeInFlight) {
311 eSep = ParticleTable::getProtonSeparationEne 311 eSep = ParticleTable::getProtonSeparationEnergy();
312 } else if (getAType() == AnnihilationType::NT 312 } else if (getAType() == AnnihilationType::NTypeInFlight) {
313 eSep = ParticleTable::getNeutronSeparationEn 313 eSep = ParticleTable::getNeutronSeparationEnergy();
314 } else if (getAType() == AnnihilationType::Nb 314 } else if (getAType() == AnnihilationType::NbarPTypeInFlight) {
315 eSep = ParticleTable::getProtonSeparationEne 315 eSep = ParticleTable::getProtonSeparationEnergy();
316 } else if (getAType() == AnnihilationType::Nb 316 } else if (getAType() == AnnihilationType::NbarNTypeInFlight) {
317 eSep = ParticleTable::getNeutronSeparationEn 317 eSep = ParticleTable::getNeutronSeparationEnergy();
318 } 318 }
319 319
320 if (eSep > 0. && (totalEnergy - initialInter 320 if (eSep > 0. && (totalEnergy - initialInternalEnergy - separationEnergies - eSep) < 0.) {
321 INCL_DEBUG("Negative Excitation Energy du 321 INCL_DEBUG("Negative Excitation Energy due to a Nbar Annihilation process (separation energy of the nucleon annihilated...); E* = " << (totalEnergy - initialInternalEnergy - separationEnergies - eSep) << '\n');
322 } 322 }
323 323
324 return totalEnergy - initialInternalEnergy - 324 return totalEnergy - initialInternalEnergy - separationEnergies - eSep;
325 325
326 } 326 }
327 327
328 //thePotential->getSeparationEnergy(Proton) 328 //thePotential->getSeparationEnergy(Proton)
329 329
330 std::string Nucleus::print() 330 std::string Nucleus::print()
331 { 331 {
332 std::stringstream ss; 332 std::stringstream ss;
333 ss << "Particles in the nucleus:" << '\n' 333 ss << "Particles in the nucleus:" << '\n'
334 << "Inside:" << '\n'; 334 << "Inside:" << '\n';
335 G4int counter = 1; 335 G4int counter = 1;
336 ParticleList const &inside = theStore->get 336 ParticleList const &inside = theStore->getParticles();
337 for(ParticleIter p=inside.begin(), e=insid 337 for(ParticleIter p=inside.begin(), e=inside.end(); p!=e; ++p) {
338 ss << "index = " << counter << '\n' 338 ss << "index = " << counter << '\n'
339 << (*p)->print(); 339 << (*p)->print();
340 counter++; 340 counter++;
341 } 341 }
342 ss <<"Outgoing:" << '\n'; 342 ss <<"Outgoing:" << '\n';
343 ParticleList const &outgoing = theStore->g 343 ParticleList const &outgoing = theStore->getOutgoingParticles();
344 for(ParticleIter p=outgoing.begin(), e=out 344 for(ParticleIter p=outgoing.begin(), e=outgoing.end(); p!=e; ++p)
345 ss << (*p)->print(); 345 ss << (*p)->print();
346 346
347 return ss.str(); 347 return ss.str();
348 } 348 }
349 349
350 G4bool Nucleus::decayOutgoingDeltas() { 350 G4bool Nucleus::decayOutgoingDeltas() {
351 ParticleList const &out = theStore->getOut 351 ParticleList const &out = theStore->getOutgoingParticles();
352 ParticleList deltas; 352 ParticleList deltas;
353 for(ParticleIter i=out.begin(), e=out.end( 353 for(ParticleIter i=out.begin(), e=out.end(); i!=e; ++i) {
354 if((*i)->isDelta()) deltas.push_back((*i 354 if((*i)->isDelta()) deltas.push_back((*i));
355 } 355 }
356 if(deltas.empty()) return false; 356 if(deltas.empty()) return false;
357 357
358 for(ParticleIter i=deltas.begin(), e=delta 358 for(ParticleIter i=deltas.begin(), e=deltas.end(); i!=e; ++i) {
359 INCL_DEBUG("Decay outgoing delta particl 359 INCL_DEBUG("Decay outgoing delta particle:" << '\n'
360 << (*i)->print() << '\n'); 360 << (*i)->print() << '\n');
361 const ThreeVector beta = -(*i)->boostVec 361 const ThreeVector beta = -(*i)->boostVector();
362 const G4double deltaMass = (*i)->getMass 362 const G4double deltaMass = (*i)->getMass();
363 363
364 // Set the delta momentum to zero and sa 364 // Set the delta momentum to zero and sample the decay in the CM frame.
365 // This makes life simpler if we are usi 365 // This makes life simpler if we are using real particle masses.
366 (*i)->setMomentum(ThreeVector()); 366 (*i)->setMomentum(ThreeVector());
367 (*i)->setEnergy((*i)->getMass()); 367 (*i)->setEnergy((*i)->getMass());
368 368
369 // Use a DecayAvatar 369 // Use a DecayAvatar
370 IAvatar *decay = new DecayAvatar((*i), 0 370 IAvatar *decay = new DecayAvatar((*i), 0.0, NULL);
371 FinalState *fs = decay->getFinalState(); 371 FinalState *fs = decay->getFinalState();
372 Particle * const pion = fs->getCreatedPa 372 Particle * const pion = fs->getCreatedParticles().front();
373 Particle * const nucleon = fs->getModifi 373 Particle * const nucleon = fs->getModifiedParticles().front();
374 374
375 // Adjust the decay momentum if we are u 375 // Adjust the decay momentum if we are using the real masses
376 const G4double decayMomentum = Kinematic 376 const G4double decayMomentum = KinematicsUtils::momentumInCM(deltaMass,
377 nucleon->getTableMass(), 377 nucleon->getTableMass(),
378 pion->getTableMass()); 378 pion->getTableMass());
379 ThreeVector newMomentum = pion->getMomen 379 ThreeVector newMomentum = pion->getMomentum();
380 newMomentum *= decayMomentum / newMoment 380 newMomentum *= decayMomentum / newMomentum.mag();
381 381
382 pion->setTableMass(); 382 pion->setTableMass();
383 pion->setMomentum(newMomentum); 383 pion->setMomentum(newMomentum);
384 pion->adjustEnergyFromMomentum(); 384 pion->adjustEnergyFromMomentum();
385 pion->setEmissionTime(nucleon->getEmissi 385 pion->setEmissionTime(nucleon->getEmissionTime());
386 pion->boost(beta); 386 pion->boost(beta);
387 pion->setBiasCollisionVector(nucleon->ge 387 pion->setBiasCollisionVector(nucleon->getBiasCollisionVector());
388 388
389 nucleon->setTableMass(); 389 nucleon->setTableMass();
390 nucleon->setMomentum(-newMomentum); 390 nucleon->setMomentum(-newMomentum);
391 nucleon->adjustEnergyFromMomentum(); 391 nucleon->adjustEnergyFromMomentum();
392 nucleon->boost(beta); 392 nucleon->boost(beta);
393 393
394 theStore->addToOutgoing(pion); 394 theStore->addToOutgoing(pion);
395 395
396 delete fs; 396 delete fs;
397 delete decay; 397 delete decay;
398 } 398 }
399 399
400 return true; 400 return true;
401 } 401 }
402 402
403 G4bool Nucleus::decayInsideDeltas() { 403 G4bool Nucleus::decayInsideDeltas() {
404 /* If there is a pion potential, do nothin 404 /* If there is a pion potential, do nothing (deltas will be counted as
405 * excitation energy). 405 * excitation energy).
406 * If, however, the remnant is unphysical 406 * If, however, the remnant is unphysical (Z<0 or Z>A), force the deltas to
407 * decay and get rid of all the pions. In 407 * decay and get rid of all the pions. In case you're wondering, you can
408 * end up with Z<0 or Z>A if the remnant c 408 * end up with Z<0 or Z>A if the remnant contains more pi- than protons or
409 * more pi+ than neutrons, respectively. 409 * more pi+ than neutrons, respectively.
410 */ 410 */
411 const G4bool unphysicalRemnant = (theZ<0 | 411 const G4bool unphysicalRemnant = (theZ<0 || theZ>theA);
412 if(thePotential->hasPionPotential() && !un 412 if(thePotential->hasPionPotential() && !unphysicalRemnant)
413 return false; 413 return false;
414 414
415 // Build a list of deltas (avoid modifying 415 // Build a list of deltas (avoid modifying the list you are iterating on).
416 ParticleList const &inside = theStore->get 416 ParticleList const &inside = theStore->getParticles();
417 ParticleList deltas; 417 ParticleList deltas;
418 for(ParticleIter i=inside.begin(), e=insid 418 for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i)
419 if((*i)->isDelta()) deltas.push_back((*i 419 if((*i)->isDelta()) deltas.push_back((*i));
420 420
421 // Loop over the deltas, make them decay 421 // Loop over the deltas, make them decay
422 for(ParticleIter i=deltas.begin(), e=delta 422 for(ParticleIter i=deltas.begin(), e=deltas.end(); i!=e; ++i) {
423 INCL_DEBUG("Decay inside delta particle: 423 INCL_DEBUG("Decay inside delta particle:" << '\n'
424 << (*i)->print() << '\n'); 424 << (*i)->print() << '\n');
425 // Create a forced-decay avatar. Note th 425 // Create a forced-decay avatar. Note the last boolean parameter. Note
426 // also that if the remnant is unphysica 426 // also that if the remnant is unphysical we more or less explicitly give
427 // up energy conservation and CDPP by pa 427 // up energy conservation and CDPP by passing a NULL pointer for the
428 // nucleus. 428 // nucleus.
429 IAvatar *decay; 429 IAvatar *decay;
430 if(unphysicalRemnant) { 430 if(unphysicalRemnant) {
431 INCL_WARN("Forcing delta decay inside 431 INCL_WARN("Forcing delta decay inside an unphysical remnant (A=" << theA
432 << ", Z=" << theZ << "). Might le 432 << ", Z=" << theZ << "). Might lead to energy-violation warnings."
433 << '\n'); 433 << '\n');
434 decay = new DecayAvatar((*i), 0.0, NUL 434 decay = new DecayAvatar((*i), 0.0, NULL, true);
435 } else 435 } else
436 decay = new DecayAvatar((*i), 0.0, thi 436 decay = new DecayAvatar((*i), 0.0, this, true);
437 FinalState *fs = decay->getFinalState(); 437 FinalState *fs = decay->getFinalState();
438 438
439 // The pion can be ejected only if we ma 439 // The pion can be ejected only if we managed to satisfy energy
440 // conservation and if pion emission doe 440 // conservation and if pion emission does not lead to negative excitation
441 // energies. 441 // energies.
442 if(fs->getValidity()==ValidFS) { 442 if(fs->getValidity()==ValidFS) {
443 // Apply the final state to the nucleu 443 // Apply the final state to the nucleus
444 applyFinalState(fs); 444 applyFinalState(fs);
445 } 445 }
446 delete fs; 446 delete fs;
447 delete decay; 447 delete decay;
448 } 448 }
449 449
450 // If the remnant is unphysical, emit all 450 // If the remnant is unphysical, emit all the pions
451 if(unphysicalRemnant) { 451 if(unphysicalRemnant) {
452 INCL_DEBUG("Remnant is unphysical: Z=" < 452 INCL_DEBUG("Remnant is unphysical: Z=" << theZ << ", A=" << theA << ", emitting all the pions" << '\n');
453 emitInsidePions(); 453 emitInsidePions();
454 } 454 }
455 455
456 return true; 456 return true;
457 } 457 }
458 458
459 G4bool Nucleus::decayInsideStrangeParticles( 459 G4bool Nucleus::decayInsideStrangeParticles() {
460 460
461 /* Transform each strange particles into a 461 /* Transform each strange particles into a lambda
462 * Every Kaon (KPlus and KZero) are emited 462 * Every Kaon (KPlus and KZero) are emited
463 */ 463 */
464 const G4bool unphysicalRemnant = (theZ<0 | 464 const G4bool unphysicalRemnant = (theZ<0 || theZ>theA);
465 if(unphysicalRemnant){ 465 if(unphysicalRemnant){
466 emitInsideStrangeParticles(); 466 emitInsideStrangeParticles();
467 INCL_WARN("Remnant is unphysical: Z=" 467 INCL_WARN("Remnant is unphysical: Z=" << theZ << ", A=" << theA << ", too much strange particles? -> all emit" << '\n');
468 return false; 468 return false;
469 } 469 }
470 470
471 /* Build a list of particles with a strang 471 /* Build a list of particles with a strangeness == -1 except Lambda,
472 * and two other one for proton and neutro 472 * and two other one for proton and neutron
473 */ 473 */
474 ParticleList const &inside = theStore->get 474 ParticleList const &inside = theStore->getParticles();
475 ParticleList stranges; 475 ParticleList stranges;
476 ParticleList protons; 476 ParticleList protons;
477 ParticleList neutrons; 477 ParticleList neutrons;
478 for(ParticleIter i=inside.begin(), e=insid 478 for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i){
479 if((*i)->isSigma() || (*i)->isAntiKaon 479 if((*i)->isSigma() || (*i)->isAntiKaon()) stranges.push_back((*i));
480 else if((*i)->isNucleon() && (*i)->get 480 else if((*i)->isNucleon() && (*i)->getZ() == 1) protons.push_back((*i));
481 else if((*i)->isNucleon() && (*i)->get 481 else if((*i)->isNucleon() && (*i)->getZ() == 0) neutrons.push_back((*i));
482 } 482 }
483 483
484 if((stranges.size() > protons.size()) || ( 484 if((stranges.size() > protons.size()) || (stranges.size() > neutrons.size())){
485 INCL_WARN("Remnant is unphysical: Npro 485 INCL_WARN("Remnant is unphysical: Nproton=" << protons.size() << ", Nneutron=" << neutrons.size() << ", Strange particles : " << stranges.size() << '\n');
486 emitInsideStrangeParticles(); 486 emitInsideStrangeParticles();
487 return false; 487 return false;
488 } 488 }
489 489
490 // Loop over the strange particles, make t 490 // Loop over the strange particles, make them absorbe
491 ParticleIter protonIter = protons.begin(); 491 ParticleIter protonIter = protons.begin();
492 ParticleIter neutronIter = neutrons.begin( 492 ParticleIter neutronIter = neutrons.begin();
493 for(ParticleIter i=stranges.begin(), e=str 493 for(ParticleIter i=stranges.begin(), e=stranges.end(); i!=e; ++i) {
494 INCL_DEBUG("Absorbe inside strange parti 494 INCL_DEBUG("Absorbe inside strange particles:" << '\n'
495 << (*i)->print() << '\n'); 495 << (*i)->print() << '\n');
496 IAvatar *decay; 496 IAvatar *decay;
497 if((*i)->getType() == SigmaMinus){ 497 if((*i)->getType() == SigmaMinus){
498 decay = new DecayAvatar((*protonIter 498 decay = new DecayAvatar((*protonIter), (*i), 0.0, this, true);
499 ++protonIter; 499 ++protonIter;
500 } 500 }
501 else if((*i)->getType() == SigmaPlus){ 501 else if((*i)->getType() == SigmaPlus){
502 decay = new DecayAvatar((*neutronIte 502 decay = new DecayAvatar((*neutronIter), (*i), 0.0, this, true);
503 ++neutronIter; 503 ++neutronIter;
504 } 504 }
505 else if(Random::shoot()*(protons.size() 505 else if(Random::shoot()*(protons.size() + neutrons.size()) < protons.size()){
506 decay = new DecayAvatar((*protonIter 506 decay = new DecayAvatar((*protonIter), (*i), 0.0, this, true);
507 ++protonIter; 507 ++protonIter;
508 } 508 }
509 else { 509 else {
510 decay = new DecayAvatar((*neutronIte 510 decay = new DecayAvatar((*neutronIter), (*i), 0.0, this, true);
511 ++neutronIter; 511 ++neutronIter;
512 } 512 }
513 FinalState *fs = decay->getFinalState(); 513 FinalState *fs = decay->getFinalState();
514 applyFinalState(fs); 514 applyFinalState(fs);
515 delete fs; 515 delete fs;
516 delete decay; 516 delete decay;
517 } 517 }
518 518
519 return true; 519 return true;
520 } 520 }
521 521
522 G4bool Nucleus::decayOutgoingPionResonances( 522 G4bool Nucleus::decayOutgoingPionResonances(G4double timeThreshold) {
523 ParticleList const &out = theStore->ge 523 ParticleList const &out = theStore->getOutgoingParticles();
524 ParticleList pionResonances; 524 ParticleList pionResonances;
525 for(ParticleIter i=out.begin(), e=out. 525 for(ParticleIter i=out.begin(), e=out.end(); i!=e; ++i) {
526 // if((*i)->isEta() || (*i)->isOmeg 526 // if((*i)->isEta() || (*i)->isOmega()) pionResonances.push_back((*i));
527 if(((*i)->isEta() && timeThreshold 527 if(((*i)->isEta() && timeThreshold > ParticleTable::getWidth(Eta)) || ((*i)->isOmega() && timeThreshold > ParticleTable::getWidth(Omega))) pionResonances.push_back((*i));
528 } 528 }
529 if(pionResonances.empty()) return fals 529 if(pionResonances.empty()) return false;
530 530
531 for(ParticleIter i=pionResonances.begi 531 for(ParticleIter i=pionResonances.begin(), e=pionResonances.end(); i!=e; ++i) {
532 INCL_DEBUG("Decay outgoing pionRes 532 INCL_DEBUG("Decay outgoing pionResonances particle:" << '\n'
533 << (*i)->print() << '\n 533 << (*i)->print() << '\n');
534 const ThreeVector beta = -(*i)->bo 534 const ThreeVector beta = -(*i)->boostVector();
535 const G4double pionResonanceMass = 535 const G4double pionResonanceMass = (*i)->getMass();
536 536
537 // Set the pionResonance momentum 537 // Set the pionResonance momentum to zero and sample the decay in the CM frame.
538 // This makes life simpler if we a 538 // This makes life simpler if we are using real particle masses.
539 (*i)->setMomentum(ThreeVector()); 539 (*i)->setMomentum(ThreeVector());
540 (*i)->setEnergy((*i)->getMass()); 540 (*i)->setEnergy((*i)->getMass());
541 541
542 // Use a DecayAvatar 542 // Use a DecayAvatar
543 IAvatar *decay = new DecayAvatar(( 543 IAvatar *decay = new DecayAvatar((*i), 0.0, NULL);
544 FinalState *fs = decay->getFinalSt 544 FinalState *fs = decay->getFinalState();
545 545
546 Particle * const theModifiedPartic 546 Particle * const theModifiedParticle = fs->getModifiedParticles().front();
547 ParticleList const &created = fs-> 547 ParticleList const &created = fs->getCreatedParticles();
548 Particle * const theCreatedParticl 548 Particle * const theCreatedParticle1 = created.front();
549 549
550 if (created.size() == 1) { 550 if (created.size() == 1) {
551 551
552 // Adjust the decay momentum i 552 // Adjust the decay momentum if we are using the real masses
553 const G4double decayMomentum = 553 const G4double decayMomentum = KinematicsUtils::momentumInCM(pionResonanceMass,theModifiedParticle->getTableMass(),theCreatedParticle1->getTableMass());
554 ThreeVector newMomentum = theC 554 ThreeVector newMomentum = theCreatedParticle1->getMomentum();
555 newMomentum *= decayMomentum / 555 newMomentum *= decayMomentum / newMomentum.mag();
556 556
557 theCreatedParticle1->setTableM 557 theCreatedParticle1->setTableMass();
558 theCreatedParticle1->setMoment 558 theCreatedParticle1->setMomentum(newMomentum);
559 theCreatedParticle1->adjustEne 559 theCreatedParticle1->adjustEnergyFromMomentum();
560 //theCreatedParticle1->setEmis 560 //theCreatedParticle1->setEmissionTime(nucleon->getEmissionTime());
561 theCreatedParticle1->boost(bet 561 theCreatedParticle1->boost(beta);
562 theCreatedParticle1->setBiasCo 562 theCreatedParticle1->setBiasCollisionVector(theModifiedParticle->getBiasCollisionVector());
563 563
564 theModifiedParticle->setTableM 564 theModifiedParticle->setTableMass();
565 theModifiedParticle->setMoment 565 theModifiedParticle->setMomentum(-newMomentum);
566 theModifiedParticle->adjustEne 566 theModifiedParticle->adjustEnergyFromMomentum();
567 theModifiedParticle->boost(bet 567 theModifiedParticle->boost(beta);
568 568
569 theStore->addToOutgoing(theCre 569 theStore->addToOutgoing(theCreatedParticle1);
570 } 570 }
571 else if (created.size() == 2) { 571 else if (created.size() == 2) {
572 Particle * const theCreatedPar 572 Particle * const theCreatedParticle2 = created.back();
573 573
574 theCreatedParticle1->boost(bet 574 theCreatedParticle1->boost(beta);
575 theCreatedParticle1->setBiasCo 575 theCreatedParticle1->setBiasCollisionVector(theModifiedParticle->getBiasCollisionVector());
576 theCreatedParticle2->boost(bet 576 theCreatedParticle2->boost(beta);
577 theCreatedParticle2->setBiasCo 577 theCreatedParticle2->setBiasCollisionVector(theModifiedParticle->getBiasCollisionVector());
578 theModifiedParticle->boost(bet 578 theModifiedParticle->boost(beta);
579 579
580 theStore->addToOutgoing(theCre 580 theStore->addToOutgoing(theCreatedParticle1);
581 theStore->addToOutgoing(theCre 581 theStore->addToOutgoing(theCreatedParticle2);
582 } 582 }
583 else { 583 else {
584 INCL_ERROR("Wrong number (< 2) 584 INCL_ERROR("Wrong number (< 2) of created particles during the decay of a pion resonance");
585 } 585 }
586 delete fs; 586 delete fs;
587 delete decay; 587 delete decay;
588 } 588 }
589 589
590 return true; 590 return true;
591 } 591 }
592 592
593 G4bool Nucleus::decayOutgoingSigmaZero(G4dou 593 G4bool Nucleus::decayOutgoingSigmaZero(G4double timeThreshold) {
594 ParticleList const &out = theStore->ge 594 ParticleList const &out = theStore->getOutgoingParticles();
595 ParticleList neutralsigma; 595 ParticleList neutralsigma;
596 for(ParticleIter i=out.begin(), e=out. 596 for(ParticleIter i=out.begin(), e=out.end(); i!=e; ++i) {
597 if((*i)->getType() == SigmaZero && 597 if((*i)->getType() == SigmaZero && timeThreshold > ParticleTable::getWidth(SigmaZero)) neutralsigma.push_back((*i));
598 } 598 }
599 if(neutralsigma.empty()) return false; 599 if(neutralsigma.empty()) return false;
600 600
601 for(ParticleIter i=neutralsigma.begin( 601 for(ParticleIter i=neutralsigma.begin(), e=neutralsigma.end(); i!=e; ++i) {
602 INCL_DEBUG("Decay outgoing neutral 602 INCL_DEBUG("Decay outgoing neutral sigma:" << '\n'
603 << (*i)->print() << '\n 603 << (*i)->print() << '\n');
604 const ThreeVector beta = -(*i)->bo 604 const ThreeVector beta = -(*i)->boostVector();
605 const G4double neutralsigmaMass = 605 const G4double neutralsigmaMass = (*i)->getMass();
606 606
607 // Set the neutral sigma momentum 607 // Set the neutral sigma momentum to zero and sample the decay in the CM frame.
608 // This makes life simpler if we a 608 // This makes life simpler if we are using real particle masses.
609 (*i)->setMomentum(ThreeVector()); 609 (*i)->setMomentum(ThreeVector());
610 (*i)->setEnergy((*i)->getMass()); 610 (*i)->setEnergy((*i)->getMass());
611 611
612 // Use a DecayAvatar 612 // Use a DecayAvatar
613 IAvatar *decay = new DecayAvatar(( 613 IAvatar *decay = new DecayAvatar((*i), 0.0, NULL);
614 FinalState *fs = decay->getFinalSt 614 FinalState *fs = decay->getFinalState();
615 615
616 Particle * const theModifiedPartic 616 Particle * const theModifiedParticle = fs->getModifiedParticles().front();
617 ParticleList const &created = fs-> 617 ParticleList const &created = fs->getCreatedParticles();
618 Particle * const theCreatedParticl 618 Particle * const theCreatedParticle = created.front();
619 619
620 if (created.size() == 1) { 620 if (created.size() == 1) {
621 621
622 // Adjust the decay momentum i 622 // Adjust the decay momentum if we are using the real masses
623 const G4double decayMomentum = 623 const G4double decayMomentum = KinematicsUtils::momentumInCM(neutralsigmaMass,theModifiedParticle->getTableMass(),theCreatedParticle->getTableMass());
624 ThreeVector newMomentum = theC 624 ThreeVector newMomentum = theCreatedParticle->getMomentum();
625 newMomentum *= decayMomentum / 625 newMomentum *= decayMomentum / newMomentum.mag();
626 626
627 theCreatedParticle->setTableMa 627 theCreatedParticle->setTableMass();
628 theCreatedParticle->setMomentu 628 theCreatedParticle->setMomentum(newMomentum);
629 theCreatedParticle->adjustEner 629 theCreatedParticle->adjustEnergyFromMomentum();
630 theCreatedParticle->boost(beta 630 theCreatedParticle->boost(beta);
631 theCreatedParticle->setBiasCol 631 theCreatedParticle->setBiasCollisionVector(theModifiedParticle->getBiasCollisionVector());
632 632
633 theModifiedParticle->setTableM 633 theModifiedParticle->setTableMass();
634 theModifiedParticle->setMoment 634 theModifiedParticle->setMomentum(-newMomentum);
635 theModifiedParticle->adjustEne 635 theModifiedParticle->adjustEnergyFromMomentum();
636 theModifiedParticle->boost(bet 636 theModifiedParticle->boost(beta);
637 637
638 theStore->addToOutgoing(theCre 638 theStore->addToOutgoing(theCreatedParticle);
639 } 639 }
640 else { 640 else {
641 INCL_ERROR("Wrong number (!= 1 641 INCL_ERROR("Wrong number (!= 1) of created particles during the decay of a sigma zero");
642 } 642 }
643 delete fs; 643 delete fs;
644 delete decay; 644 delete decay;
645 } 645 }
646 646
647 return true; 647 return true;
648 } 648 }
649 649
650 G4bool Nucleus::decayOutgoingNeutralKaon() { 650 G4bool Nucleus::decayOutgoingNeutralKaon() {
651 ParticleList const &out = theStore->ge 651 ParticleList const &out = theStore->getOutgoingParticles();
652 ParticleList neutralkaon; 652 ParticleList neutralkaon;
653 for(ParticleIter i=out.begin(), e=out. 653 for(ParticleIter i=out.begin(), e=out.end(); i!=e; ++i) {
654 if((*i)->getType() == KZero || (* 654 if((*i)->getType() == KZero || (*i)->getType() == KZeroBar) neutralkaon.push_back((*i));
655 } 655 }
656 if(neutralkaon.empty()) return false; 656 if(neutralkaon.empty()) return false;
657 657
658 for(ParticleIter i=neutralkaon.begin() 658 for(ParticleIter i=neutralkaon.begin(), e=neutralkaon.end(); i!=e; ++i) {
659 INCL_DEBUG("Transform outgoing neu 659 INCL_DEBUG("Transform outgoing neutral kaon:" << '\n'
660 << (*i)->print() << '\n 660 << (*i)->print() << '\n');
661 661
662 // Use a DecayAvatar 662 // Use a DecayAvatar
663 IAvatar *decay = new DecayAvatar(( 663 IAvatar *decay = new DecayAvatar((*i), 0.0, NULL);
664 FinalState *fs = decay->getFinalSt 664 FinalState *fs = decay->getFinalState();
665 665
666 delete fs; 666 delete fs;
667 delete decay; 667 delete decay;
668 } 668 }
669 669
670 return true; 670 return true;
671 } 671 }
672 672
673 G4bool Nucleus::decayOutgoingClusters() { 673 G4bool Nucleus::decayOutgoingClusters() {
674 ParticleList const &out = theStore->getOut 674 ParticleList const &out = theStore->getOutgoingParticles();
675 ParticleList clusters; 675 ParticleList clusters;
676 for(ParticleIter i=out.begin(), e=out.end( 676 for(ParticleIter i=out.begin(), e=out.end(); i!=e; ++i) {
677 if((*i)->isCluster()) clusters.push_back 677 if((*i)->isCluster()) clusters.push_back((*i));
678 } 678 }
679 if(clusters.empty()) return false; 679 if(clusters.empty()) return false;
680 680
681 for(ParticleIter i=clusters.begin(), e=clu 681 for(ParticleIter i=clusters.begin(), e=clusters.end(); i!=e; ++i) {
682 Cluster *cluster = dynamic_cast<Cluster* 682 Cluster *cluster = dynamic_cast<Cluster*>(*i); // Can't avoid using a cast here
683 // assert(cluster); 683 // assert(cluster);
684 #ifdef INCLXX_IN_GEANT4_MODE 684 #ifdef INCLXX_IN_GEANT4_MODE
685 if(!cluster) 685 if(!cluster)
686 continue; 686 continue;
687 #endif 687 #endif
688 cluster->deleteParticles(); // Don't nee 688 cluster->deleteParticles(); // Don't need them
689 ParticleList decayProducts = ClusterDeca 689 ParticleList decayProducts = ClusterDecay::decay(cluster);
690 for(ParticleIter j=decayProducts.begin() 690 for(ParticleIter j=decayProducts.begin(), end=decayProducts.end(); j!=end; ++j){
691 (*j)->setBiasCollisionVector(cluster-> 691 (*j)->setBiasCollisionVector(cluster->getBiasCollisionVector());
692 theStore->addToOutgoing(*j); 692 theStore->addToOutgoing(*j);
693 } 693 }
694 } 694 }
695 return true; 695 return true;
696 } 696 }
697 697
698 G4bool Nucleus::decayMe() { 698 G4bool Nucleus::decayMe() {
699 // Do the phase-space decay only if Z=0 or 699 // Do the phase-space decay only if Z=0 or N=0
700 if(theA<=1 || (theZ!=0 && (theA+theS)!=the 700 if(theA<=1 || (theZ!=0 && (theA+theS)!=theZ))
701 return false; 701 return false;
702 702
703 ParticleList decayProducts = ClusterDecay: 703 ParticleList decayProducts = ClusterDecay::decay(this);
704 for(ParticleIter j=decayProducts.begin(), 704 for(ParticleIter j=decayProducts.begin(), e=decayProducts.end(); j!=e; ++j){
705 (*j)->setBiasCollisionVector(this->getBi 705 (*j)->setBiasCollisionVector(this->getBiasCollisionVector());
706 theStore->addToOutgoing(*j); 706 theStore->addToOutgoing(*j);
707 } 707 }
708 708
709 return true; 709 return true;
710 } 710 }
711 711
712 void Nucleus::emitInsidePions() { 712 void Nucleus::emitInsidePions() {
713 /* Forcing emissions of all pions in the n 713 /* Forcing emissions of all pions in the nucleus. This probably violates
714 * energy conservation (although the compu 714 * energy conservation (although the computation of the recoil kinematics
715 * might sweep this under the carpet). 715 * might sweep this under the carpet).
716 */ 716 */
717 INCL_WARN("Forcing emissions of all pions 717 INCL_WARN("Forcing emissions of all pions in the nucleus." << '\n');
718 718
719 // Emit the pions with this kinetic energy 719 // Emit the pions with this kinetic energy
720 const G4double tinyPionEnergy = 0.1; // Me 720 const G4double tinyPionEnergy = 0.1; // MeV
721 721
722 // Push out the emitted pions 722 // Push out the emitted pions
723 ParticleList const &inside = theStore->get 723 ParticleList const &inside = theStore->getParticles();
724 ParticleList toEject; 724 ParticleList toEject;
725 for(ParticleIter i=inside.begin(), e=insid 725 for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i) {
726 if((*i)->isPion()) { 726 if((*i)->isPion()) {
727 Particle * const thePion = *i; 727 Particle * const thePion = *i;
728 INCL_DEBUG("Forcing emission of the fo 728 INCL_DEBUG("Forcing emission of the following particle: "
729 << thePion->print() << '\n' 729 << thePion->print() << '\n');
730 thePion->setEmissionTime(theStore->get 730 thePion->setEmissionTime(theStore->getBook().getCurrentTime());
731 // Correction for real masses 731 // Correction for real masses
732 const G4double theQValueCorrection = t 732 const G4double theQValueCorrection = thePion->getEmissionQValueCorrection(theA,theZ,theS);
733 const G4double kineticEnergyOutside = 733 const G4double kineticEnergyOutside = thePion->getKineticEnergy() - thePion->getPotentialEnergy() + theQValueCorrection;
734 thePion->setTableMass(); 734 thePion->setTableMass();
735 if(kineticEnergyOutside > 0.0) 735 if(kineticEnergyOutside > 0.0)
736 thePion->setEnergy(thePion->getMass( 736 thePion->setEnergy(thePion->getMass()+kineticEnergyOutside);
737 else 737 else
738 thePion->setEnergy(thePion->getMass( 738 thePion->setEnergy(thePion->getMass()+tinyPionEnergy);
739 thePion->adjustMomentumFromEnergy(); 739 thePion->adjustMomentumFromEnergy();
740 thePion->setPotentialEnergy(0.); 740 thePion->setPotentialEnergy(0.);
741 theZ -= thePion->getZ(); 741 theZ -= thePion->getZ();
742 toEject.push_back(thePion); 742 toEject.push_back(thePion);
743 } 743 }
744 } 744 }
745 for(ParticleIter i=toEject.begin(), e=toEj 745 for(ParticleIter i=toEject.begin(), e=toEject.end(); i!=e; ++i) {
746 theStore->particleHasBeenEjected(*i); 746 theStore->particleHasBeenEjected(*i);
747 theStore->addToOutgoing(*i); 747 theStore->addToOutgoing(*i);
748 (*i)->setParticleBias(Particle::getTotal 748 (*i)->setParticleBias(Particle::getTotalBias());
749 } 749 }
750 } 750 }
751 751
752 void Nucleus::emitInsideStrangeParticles() { 752 void Nucleus::emitInsideStrangeParticles() {
753 /* Forcing emissions of Sigmas and antiKao 753 /* Forcing emissions of Sigmas and antiKaons.
754 * This probably violates energy conservat 754 * This probably violates energy conservation
755 * (although the computation of the recoil 755 * (although the computation of the recoil kinematics
756 * might sweep this under the carpet). 756 * might sweep this under the carpet).
757 */ 757 */
758 INCL_DEBUG("Forcing emissions of all stran 758 INCL_DEBUG("Forcing emissions of all strange particles in the nucleus." << '\n');
759 759
760 // Emit the strange particles with this ki 760 // Emit the strange particles with this kinetic energy
761 const G4double tinyEnergy = 0.1; // MeV 761 const G4double tinyEnergy = 0.1; // MeV
762 762
763 // Push out the emitted strange particles 763 // Push out the emitted strange particles
764 ParticleList const &inside = theStore->get 764 ParticleList const &inside = theStore->getParticles();
765 ParticleList toEject; 765 ParticleList toEject;
766 for(ParticleIter i=inside.begin(), e=insid 766 for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i) {
767 if((*i)->isSigma() || (*i)->isAntiKaon() 767 if((*i)->isSigma() || (*i)->isAntiKaon()) {
768 Particle * const theParticle = *i; 768 Particle * const theParticle = *i;
769 INCL_DEBUG("Forcing emission of the fo 769 INCL_DEBUG("Forcing emission of the following particle: "
770 << theParticle->print() << 770 << theParticle->print() << '\n');
771 theParticle->setEmissionTime(theStore- 771 theParticle->setEmissionTime(theStore->getBook().getCurrentTime());
772 // Correction for real masses 772 // Correction for real masses
773 const G4double theQValueCorrection = t 773 const G4double theQValueCorrection = theParticle->getEmissionQValueCorrection(theA,theZ,theS); // Does it work for strange particles? should be check
774 const G4double kineticEnergyOutside = 774 const G4double kineticEnergyOutside = theParticle->getKineticEnergy() - theParticle->getPotentialEnergy() + theQValueCorrection;
775 theParticle->setTableMass(); 775 theParticle->setTableMass();
776 if(kineticEnergyOutside > 0.0) 776 if(kineticEnergyOutside > 0.0)
777 theParticle->setEnergy(theParticle-> 777 theParticle->setEnergy(theParticle->getMass()+kineticEnergyOutside);
778 else 778 else
779 theParticle->setEnergy(theParticle-> 779 theParticle->setEnergy(theParticle->getMass()+tinyEnergy);
780 theParticle->adjustMomentumFromEnergy( 780 theParticle->adjustMomentumFromEnergy();
781 theParticle->setPotentialEnergy(0.); 781 theParticle->setPotentialEnergy(0.);
782 theA -= theParticle->getA(); 782 theA -= theParticle->getA();
783 theZ -= theParticle->getZ(); 783 theZ -= theParticle->getZ();
784 theS -= theParticle->getS(); 784 theS -= theParticle->getS();
785 toEject.push_back(theParticle); 785 toEject.push_back(theParticle);
786 } 786 }
787 } 787 }
788 for(ParticleIter i=toEject.begin(), e=toEj 788 for(ParticleIter i=toEject.begin(), e=toEject.end(); i!=e; ++i) {
789 theStore->particleHasBeenEjected(*i); 789 theStore->particleHasBeenEjected(*i);
790 theStore->addToOutgoing(*i); 790 theStore->addToOutgoing(*i);
791 (*i)->setParticleBias(Particle::getTotal 791 (*i)->setParticleBias(Particle::getTotalBias());
792 } 792 }
793 } 793 }
794 794
795 G4int Nucleus::emitInsideLambda() { 795 G4int Nucleus::emitInsideLambda() {
796 /* Forcing emissions of all Lambda in the 796 /* Forcing emissions of all Lambda in the nucleus.
797 * This probably violates energy conservat 797 * This probably violates energy conservation
798 * (although the computation of the recoil 798 * (although the computation of the recoil kinematics
799 * might sweep this under the carpet). 799 * might sweep this under the carpet).
800 */ 800 */
801 INCL_DEBUG("Forcing emissions of all Lambd 801 INCL_DEBUG("Forcing emissions of all Lambda in the nucleus." << '\n');
802 802
803 // Emit the Lambda with this kinetic energ 803 // Emit the Lambda with this kinetic energy
804 const G4double tinyEnergy = 0.1; // MeV 804 const G4double tinyEnergy = 0.1; // MeV
805 805
806 // Push out the emitted Lambda 806 // Push out the emitted Lambda
807 ParticleList const &inside = theStore->get 807 ParticleList const &inside = theStore->getParticles();
808 ParticleList toEject; 808 ParticleList toEject;
809 for(ParticleIter i=inside.begin(), e=insid 809 for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i) {
810 if((*i)->isLambda()) { 810 if((*i)->isLambda()) {
811 Particle * const theLambda = *i; 811 Particle * const theLambda = *i;
812 INCL_DEBUG("Forcing emission of the fo 812 INCL_DEBUG("Forcing emission of the following particle: "
813 << theLambda->print() << '\ 813 << theLambda->print() << '\n');
814 theLambda->setEmissionTime(theStore->g 814 theLambda->setEmissionTime(theStore->getBook().getCurrentTime());
815 // Correction for real masses 815 // Correction for real masses
816 const G4double theQValueCorrection = t 816 const G4double theQValueCorrection = theLambda->getEmissionQValueCorrection(theA,theZ,theS); // Does it work for strange particles? Should be check
817 const G4double kineticEnergyOutside = 817 const G4double kineticEnergyOutside = theLambda->getKineticEnergy() - theLambda->getPotentialEnergy() + theQValueCorrection;
818 theLambda->setTableMass(); 818 theLambda->setTableMass();
819 if(kineticEnergyOutside > 0.0) 819 if(kineticEnergyOutside > 0.0)
820 theLambda->setEnergy(theLambda->getM 820 theLambda->setEnergy(theLambda->getMass()+kineticEnergyOutside);
821 else 821 else
822 theLambda->setEnergy(theLambda->getM 822 theLambda->setEnergy(theLambda->getMass()+tinyEnergy);
823 theLambda->adjustMomentumFromEnergy(); 823 theLambda->adjustMomentumFromEnergy();
824 theLambda->setPotentialEnergy(0.); 824 theLambda->setPotentialEnergy(0.);
825 theA -= theLambda->getA(); 825 theA -= theLambda->getA();
826 theS -= theLambda->getS(); 826 theS -= theLambda->getS();
827 toEject.push_back(theLambda); 827 toEject.push_back(theLambda);
828 } 828 }
829 } 829 }
830 for(ParticleIter i=toEject.begin(), e=toEj 830 for(ParticleIter i=toEject.begin(), e=toEject.end(); i!=e; ++i) {
831 theStore->particleHasBeenEjected(*i); 831 theStore->particleHasBeenEjected(*i);
832 theStore->addToOutgoing(*i); 832 theStore->addToOutgoing(*i);
833 (*i)->setParticleBias(Particle::getTotal 833 (*i)->setParticleBias(Particle::getTotalBias());
834 } 834 }
835 return (G4int)toEject.size(); 835 return (G4int)toEject.size();
836 } 836 }
837 837
838 G4bool Nucleus::emitInsideKaon() { 838 G4bool Nucleus::emitInsideKaon() {
839 /* Forcing emissions of all Kaon (not anti 839 /* Forcing emissions of all Kaon (not antiKaons) in the nucleus.
840 * This probably violates energy conservat 840 * This probably violates energy conservation
841 * (although the computation of the recoil 841 * (although the computation of the recoil kinematics
842 * might sweep this under the carpet). 842 * might sweep this under the carpet).
843 */ 843 */
844 INCL_DEBUG("Forcing emissions of all Kaon 844 INCL_DEBUG("Forcing emissions of all Kaon in the nucleus." << '\n');
845 845
846 // Emit the Kaon with this kinetic energy 846 // Emit the Kaon with this kinetic energy (not supposed to append
847 const G4double tinyEnergy = 0.1; // MeV 847 const G4double tinyEnergy = 0.1; // MeV
848 848
849 // Push out the emitted kaon 849 // Push out the emitted kaon
850 ParticleList const &inside = theStore->get 850 ParticleList const &inside = theStore->getParticles();
851 ParticleList toEject; 851 ParticleList toEject;
852 for(ParticleIter i=inside.begin(), e=insid 852 for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i) {
853 if((*i)->isKaon()) { 853 if((*i)->isKaon()) {
854 Particle * const theKaon = *i; 854 Particle * const theKaon = *i;
855 INCL_DEBUG("Forcing emission of the fo 855 INCL_DEBUG("Forcing emission of the following particle: "
856 << theKaon->print() << '\n' 856 << theKaon->print() << '\n');
857 theKaon->setEmissionTime(theStore->get 857 theKaon->setEmissionTime(theStore->getBook().getCurrentTime());
858 // Correction for real masses 858 // Correction for real masses
859 const G4double theQValueCorrection = t 859 const G4double theQValueCorrection = theKaon->getEmissionQValueCorrection(theA,theZ,theS);
860 const G4double kineticEnergyOutside = 860 const G4double kineticEnergyOutside = theKaon->getKineticEnergy() - theKaon->getPotentialEnergy() + theQValueCorrection;
861 theKaon->setTableMass(); 861 theKaon->setTableMass();
862 if(kineticEnergyOutside > 0.0) 862 if(kineticEnergyOutside > 0.0)
863 theKaon->setEnergy(theKaon->getMass( 863 theKaon->setEnergy(theKaon->getMass()+kineticEnergyOutside);
864 else 864 else
865 theKaon->setEnergy(theKaon->getMass( 865 theKaon->setEnergy(theKaon->getMass()+tinyEnergy);
866 theKaon->adjustMomentumFromEnergy(); 866 theKaon->adjustMomentumFromEnergy();
867 theKaon->setPotentialEnergy(0.); 867 theKaon->setPotentialEnergy(0.);
868 theZ -= theKaon->getZ(); 868 theZ -= theKaon->getZ();
869 theS -= theKaon->getS(); 869 theS -= theKaon->getS();
870 toEject.push_back(theKaon); 870 toEject.push_back(theKaon);
871 } 871 }
872 } 872 }
873 for(ParticleIter i=toEject.begin(), e=toEj 873 for(ParticleIter i=toEject.begin(), e=toEject.end(); i!=e; ++i) {
874 theStore->particleHasBeenEjected(*i); 874 theStore->particleHasBeenEjected(*i);
875 theStore->addToOutgoing(*i); 875 theStore->addToOutgoing(*i);
876 (*i)->setParticleBias(Particle::getTotal 876 (*i)->setParticleBias(Particle::getTotalBias());
877 } 877 }
878 theNKaon -= 1; 878 theNKaon -= 1;
879 return toEject.size() != 0; 879 return toEject.size() != 0;
880 } 880 }
881 881
882 G4bool Nucleus::isEventTransparent() const { 882 G4bool Nucleus::isEventTransparent() const {
883 883
884 Book const &theBook = theStore->getBook(); 884 Book const &theBook = theStore->getBook();
885 const G4int nEventCollisions = theBook.get 885 const G4int nEventCollisions = theBook.getAcceptedCollisions();
886 const G4int nEventDecays = theBook.getAcce 886 const G4int nEventDecays = theBook.getAcceptedDecays();
887 const G4int nEventClusters = theBook.getEm 887 const G4int nEventClusters = theBook.getEmittedClusters();
888 if(nEventCollisions==0 && nEventDecays==0 888 if(nEventCollisions==0 && nEventDecays==0 && nEventClusters==0)
889 return true; 889 return true;
890 890
891 return false; 891 return false;
892 892
893 } 893 }
894 894
895 void Nucleus::computeOneNucleonRecoilKinemat 895 void Nucleus::computeOneNucleonRecoilKinematics() {
896 // We should be here only if the nucleus c 896 // We should be here only if the nucleus contains only one nucleon
897 // assert(theStore->getParticles().size()==1); 897 // assert(theStore->getParticles().size()==1);
898 898
899 // No excitation energy! 899 // No excitation energy!
900 theExcitationEnergy = 0.0; 900 theExcitationEnergy = 0.0;
901 901
902 // Move the nucleon to the outgoing list 902 // Move the nucleon to the outgoing list
903 Particle *remN = theStore->getParticles(). 903 Particle *remN = theStore->getParticles().front();
904 theA -= remN->getA(); 904 theA -= remN->getA();
905 theZ -= remN->getZ(); 905 theZ -= remN->getZ();
906 theS -= remN->getS(); 906 theS -= remN->getS();
907 theStore->particleHasBeenEjected(remN); 907 theStore->particleHasBeenEjected(remN);
908 theStore->addToOutgoing(remN); 908 theStore->addToOutgoing(remN);
909 remN->setEmissionTime(theStore->getBook(). 909 remN->setEmissionTime(theStore->getBook().getCurrentTime());
910 910
911 // Treat the special case of a remaining d 911 // Treat the special case of a remaining delta
912 if(remN->isDelta()) { 912 if(remN->isDelta()) {
913 IAvatar *decay = new DecayAvatar(remN, 0 913 IAvatar *decay = new DecayAvatar(remN, 0.0, NULL);
914 FinalState *fs = decay->getFinalState(); 914 FinalState *fs = decay->getFinalState();
915 // Eject the pion 915 // Eject the pion
916 ParticleList const &created = fs->getCre 916 ParticleList const &created = fs->getCreatedParticles();
917 for(ParticleIter j=created.begin(), e=cr 917 for(ParticleIter j=created.begin(), e=created.end(); j!=e; ++j)
918 theStore->addToOutgoing(*j); 918 theStore->addToOutgoing(*j);
919 delete fs; 919 delete fs;
920 delete decay; 920 delete decay;
921 } 921 }
922 922
923 // Do different things depending on how ma 923 // Do different things depending on how many outgoing particles we have
924 ParticleList const &outgoing = theStore->g 924 ParticleList const &outgoing = theStore->getOutgoingParticles();
925 if(outgoing.size() == 2) { 925 if(outgoing.size() == 2) {
926 926
927 INCL_DEBUG("Two particles in the outgoin 927 INCL_DEBUG("Two particles in the outgoing channel, applying exact two-body kinematics" << '\n');
928 928
929 // Can apply exact 2-body kinematics her 929 // Can apply exact 2-body kinematics here. Keep the CM emission angle of
930 // the first particle. 930 // the first particle.
931 Particle *p1 = outgoing.front(), *p2 = o 931 Particle *p1 = outgoing.front(), *p2 = outgoing.back();
932 const ThreeVector aBoostVector = incomin 932 const ThreeVector aBoostVector = incomingMomentum / initialEnergy;
933 // Boost to the initial CM 933 // Boost to the initial CM
934 p1->boost(aBoostVector); 934 p1->boost(aBoostVector);
935 const G4double sqrts = std::sqrt(initial 935 const G4double sqrts = std::sqrt(initialEnergy*initialEnergy - incomingMomentum.mag2());
936 const G4double pcm = KinematicsUtils::mo 936 const G4double pcm = KinematicsUtils::momentumInCM(sqrts, p1->getMass(), p2->getMass());
937 const G4double scale = pcm/(p1->getMomen 937 const G4double scale = pcm/(p1->getMomentum().mag());
938 // Reset the momenta 938 // Reset the momenta
939 p1->setMomentum(p1->getMomentum()*scale) 939 p1->setMomentum(p1->getMomentum()*scale);
940 p2->setMomentum(-p1->getMomentum()); 940 p2->setMomentum(-p1->getMomentum());
941 p1->adjustEnergyFromMomentum(); 941 p1->adjustEnergyFromMomentum();
942 p2->adjustEnergyFromMomentum(); 942 p2->adjustEnergyFromMomentum();
943 // Unboost 943 // Unboost
944 p1->boost(-aBoostVector); 944 p1->boost(-aBoostVector);
945 p2->boost(-aBoostVector); 945 p2->boost(-aBoostVector);
946 946
947 } else { 947 } else {
948 948
949 INCL_DEBUG("Trying to adjust final-state 949 INCL_DEBUG("Trying to adjust final-state momenta to achieve energy and momentum conservation" << '\n');
950 950
951 const G4int maxIterations=8; 951 const G4int maxIterations=8;
952 G4double totalEnergy, energyScale; 952 G4double totalEnergy, energyScale;
953 G4double val=1.E+100, oldVal=1.E+100, ol 953 G4double val=1.E+100, oldVal=1.E+100, oldOldVal=1.E+100, oldOldOldVal;
954 ThreeVector totalMomentum, deltaP; 954 ThreeVector totalMomentum, deltaP;
955 std::vector<ThreeVector> minMomenta; // 955 std::vector<ThreeVector> minMomenta; // use it to store the particle momenta that minimize the merit function
956 956
957 // Reserve the vector size 957 // Reserve the vector size
958 minMomenta.reserve(outgoing.size()); 958 minMomenta.reserve(outgoing.size());
959 959
960 // Compute the initial total momentum 960 // Compute the initial total momentum
961 totalMomentum.setX(0.0); 961 totalMomentum.setX(0.0);
962 totalMomentum.setY(0.0); 962 totalMomentum.setY(0.0);
963 totalMomentum.setZ(0.0); 963 totalMomentum.setZ(0.0);
964 for(ParticleIter i=outgoing.begin(), e=o 964 for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
965 totalMomentum += (*i)->getMomentum(); 965 totalMomentum += (*i)->getMomentum();
966 966
967 // Compute the initial total energy 967 // Compute the initial total energy
968 totalEnergy = 0.0; 968 totalEnergy = 0.0;
969 for(ParticleIter i=outgoing.begin(), e=o 969 for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
970 totalEnergy += (*i)->getEnergy(); 970 totalEnergy += (*i)->getEnergy();
971 971
972 // Iterative algorithm starts here: 972 // Iterative algorithm starts here:
973 for(G4int iterations=0; iterations < max 973 for(G4int iterations=0; iterations < maxIterations; ++iterations) {
974 974
975 // Save the old merit-function values 975 // Save the old merit-function values
976 oldOldOldVal = oldOldVal; 976 oldOldOldVal = oldOldVal;
977 oldOldVal = oldVal; 977 oldOldVal = oldVal;
978 oldVal = val; 978 oldVal = val;
979 979
980 if(iterations%2 == 0) { 980 if(iterations%2 == 0) {
981 INCL_DEBUG("Momentum step" << '\n'); 981 INCL_DEBUG("Momentum step" << '\n');
982 // Momentum step: modify all the par 982 // Momentum step: modify all the particle momenta
983 deltaP = incomingMomentum - totalMom 983 deltaP = incomingMomentum - totalMomentum;
984 G4double pOldTot = 0.0; 984 G4double pOldTot = 0.0;
985 for(ParticleIter i=outgoing.begin(), 985 for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
986 pOldTot += (*i)->getMomentum().mag 986 pOldTot += (*i)->getMomentum().mag();
987 for(ParticleIter i=outgoing.begin(), 987 for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
988 const ThreeVector mom = (*i)->getM 988 const ThreeVector mom = (*i)->getMomentum();
989 (*i)->setMomentum(mom + deltaP*mom 989 (*i)->setMomentum(mom + deltaP*mom.mag()/pOldTot);
990 (*i)->adjustEnergyFromMomentum(); 990 (*i)->adjustEnergyFromMomentum();
991 } 991 }
992 } else { 992 } else {
993 INCL_DEBUG("Energy step" << '\n'); 993 INCL_DEBUG("Energy step" << '\n');
994 // Energy step: modify all the parti 994 // Energy step: modify all the particle momenta
995 energyScale = initialEnergy/totalEne 995 energyScale = initialEnergy/totalEnergy;
996 for(ParticleIter i=outgoing.begin(), 996 for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
997 const ThreeVector mom = (*i)->getM 997 const ThreeVector mom = (*i)->getMomentum();
998 G4double pScale = ((*i)->getEnergy 998 G4double pScale = ((*i)->getEnergy()*energyScale - std::pow((*i)->getMass(),2))/mom.mag2();
999 if(pScale>0) { 999 if(pScale>0) {
1000 (*i)->setEnergy((*i)->getEnergy 1000 (*i)->setEnergy((*i)->getEnergy()*energyScale);
1001 (*i)->adjustMomentumFromEnergy( 1001 (*i)->adjustMomentumFromEnergy();
1002 } 1002 }
1003 } 1003 }
1004 } 1004 }
1005 1005
1006 // Compute the current total momentum 1006 // Compute the current total momentum and energy
1007 totalMomentum.setX(0.0); 1007 totalMomentum.setX(0.0);
1008 totalMomentum.setY(0.0); 1008 totalMomentum.setY(0.0);
1009 totalMomentum.setZ(0.0); 1009 totalMomentum.setZ(0.0);
1010 totalEnergy = 0.0; 1010 totalEnergy = 0.0;
1011 for(ParticleIter i=outgoing.begin(), 1011 for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
1012 totalMomentum += (*i)->getMomentum( 1012 totalMomentum += (*i)->getMomentum();
1013 totalEnergy += (*i)->getEnergy(); 1013 totalEnergy += (*i)->getEnergy();
1014 } 1014 }
1015 1015
1016 // Merit factor 1016 // Merit factor
1017 val = std::pow(totalEnergy - initialE 1017 val = std::pow(totalEnergy - initialEnergy,2) +
1018 0.25*(totalMomentum - incomingMomen 1018 0.25*(totalMomentum - incomingMomentum).mag2();
1019 INCL_DEBUG("Merit function: val=" << 1019 INCL_DEBUG("Merit function: val=" << val << ", oldVal=" << oldVal << ", oldOldVal=" << oldOldVal << ", oldOldOldVal=" << oldOldOldVal << '\n');
1020 1020
1021 // Store the minimum 1021 // Store the minimum
1022 if(val < oldVal) { 1022 if(val < oldVal) {
1023 INCL_DEBUG("New minimum found, stor 1023 INCL_DEBUG("New minimum found, storing the particle momenta" << '\n');
1024 minMomenta.clear(); 1024 minMomenta.clear();
1025 for(ParticleIter i=outgoing.begin() 1025 for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
1026 minMomenta.push_back((*i)->getMom 1026 minMomenta.push_back((*i)->getMomentum());
1027 } 1027 }
1028 1028
1029 // Stop the algorithm if the search d 1029 // Stop the algorithm if the search diverges
1030 if(val > oldOldVal && oldVal > oldOld 1030 if(val > oldOldVal && oldVal > oldOldOldVal) {
1031 INCL_DEBUG("Search is diverging, br 1031 INCL_DEBUG("Search is diverging, breaking out of the iteration loop: val=" << val << ", oldVal=" << oldVal << ", oldOldVal=" << oldOldVal << ", oldOldOldVal=" << oldOldOldVal << '\n');
1032 break; 1032 break;
1033 } 1033 }
1034 } 1034 }
1035 1035
1036 // We should have made at least one suc 1036 // We should have made at least one successful iteration here
1037 // assert(minMomenta.size()==outgoing.size()) 1037 // assert(minMomenta.size()==outgoing.size());
1038 1038
1039 // Apply the optimal momenta 1039 // Apply the optimal momenta
1040 INCL_DEBUG("Applying the solution" << ' 1040 INCL_DEBUG("Applying the solution" << '\n');
1041 std::vector<ThreeVector>::const_iterato 1041 std::vector<ThreeVector>::const_iterator v = minMomenta.begin();
1042 for(ParticleIter i=outgoing.begin(), e= 1042 for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i, ++v) {
1043 (*i)->setMomentum(*v); 1043 (*i)->setMomentum(*v);
1044 (*i)->adjustEnergyFromMomentum(); 1044 (*i)->adjustEnergyFromMomentum();
1045 INCL_DATABLOCK((*i)->print()); 1045 INCL_DATABLOCK((*i)->print());
1046 } 1046 }
1047 1047
1048 } 1048 }
1049 1049
1050 } 1050 }
1051 1051
1052 void Nucleus::fillEventInfo(EventInfo *even 1052 void Nucleus::fillEventInfo(EventInfo *eventInfo) {
1053 eventInfo->nParticles = 0; 1053 eventInfo->nParticles = 0;
1054 G4bool isNucleonAbsorption = false; 1054 G4bool isNucleonAbsorption = false;
1055 G4bool isPionAbsorption = false; 1055 G4bool isPionAbsorption = false;
1056 // It is possible to have pion absorption 1056 // It is possible to have pion absorption event only if the
1057 // projectile is pion. 1057 // projectile is pion.
1058 if(eventInfo->projectileType == PiPlus || 1058 if(eventInfo->projectileType == PiPlus ||
1059 eventInfo->projectileType == PiMinus | 1059 eventInfo->projectileType == PiMinus ||
1060 eventInfo->projectileType == PiZero) { 1060 eventInfo->projectileType == PiZero) {
1061 isPionAbsorption = true; 1061 isPionAbsorption = true;
1062 } 1062 }
1063 1063
1064 // Forced CN 1064 // Forced CN
1065 eventInfo->forcedCompoundNucleus = tryCN; 1065 eventInfo->forcedCompoundNucleus = tryCN;
1066 1066
1067 // Outgoing particles 1067 // Outgoing particles
1068 ParticleList const &outgoingParticles = g 1068 ParticleList const &outgoingParticles = getStore()->getOutgoingParticles();
1069 1069
1070 // Check if we have a nucleon absorption 1070 // Check if we have a nucleon absorption event: nucleon projectile
1071 // and no ejected particles. 1071 // and no ejected particles.
1072 if(outgoingParticles.size() == 0 && 1072 if(outgoingParticles.size() == 0 &&
1073 (eventInfo->projectileType == Proton | 1073 (eventInfo->projectileType == Proton ||
1074 eventInfo->projectileType == Neutron)) { 1074 eventInfo->projectileType == Neutron)) {
1075 isNucleonAbsorption = true; 1075 isNucleonAbsorption = true;
1076 } 1076 }
1077 1077
1078 // Reset the remnant counter 1078 // Reset the remnant counter
1079 eventInfo->nRemnants = 0; 1079 eventInfo->nRemnants = 0;
1080 eventInfo->history.clear(); 1080 eventInfo->history.clear();
1081 1081
1082 for(ParticleIter i=outgoingParticles.begi 1082 for(ParticleIter i=outgoingParticles.begin(), e=outgoingParticles.end(); i!=e; ++i ) {
1083 // We have a pion absorption event only 1083 // We have a pion absorption event only if the projectile is
1084 // pion and there are no ejected pions. 1084 // pion and there are no ejected pions.
1085 if(isPionAbsorption) { 1085 if(isPionAbsorption) {
1086 if((*i)->isPion()) { 1086 if((*i)->isPion()) {
1087 isPionAbsorption = false; 1087 isPionAbsorption = false;
1088 } 1088 }
1089 } 1089 }
1090 1090
1091 eventInfo->ParticleBias[eventInfo->nPar 1091 eventInfo->ParticleBias[eventInfo->nParticles] = (*i)->getParticleBias();
1092 1092
1093 #ifdef INCLXX_IN_GEANT4_MODE 1093 #ifdef INCLXX_IN_GEANT4_MODE
1094 eventInfo->A[eventInfo->nParticles] = ( 1094 eventInfo->A[eventInfo->nParticles] = (G4INCL::Short_t)(*i)->getA();
1095 eventInfo->Z[eventInfo->nParticles] = ( 1095 eventInfo->Z[eventInfo->nParticles] = (G4INCL::Short_t)(*i)->getZ();
1096 eventInfo->S[eventInfo->nParticles] = ( 1096 eventInfo->S[eventInfo->nParticles] = (G4INCL::Short_t)(*i)->getS();
1097 #else 1097 #else
1098 eventInfo->A[eventInfo->nParticles] = ( 1098 eventInfo->A[eventInfo->nParticles] = (Short_t)(*i)->getA();
1099 eventInfo->Z[eventInfo->nParticles] = ( 1099 eventInfo->Z[eventInfo->nParticles] = (Short_t)(*i)->getZ();
1100 eventInfo->S[eventInfo->nParticles] = ( 1100 eventInfo->S[eventInfo->nParticles] = (Short_t)(*i)->getS();
1101 #endif 1101 #endif
1102 eventInfo->emissionTime[eventInfo->nPar 1102 eventInfo->emissionTime[eventInfo->nParticles] = (*i)->getEmissionTime();
1103 eventInfo->EKin[eventInfo->nParticles] 1103 eventInfo->EKin[eventInfo->nParticles] = (*i)->getKineticEnergy();
1104 ThreeVector mom = (*i)->getMomentum(); 1104 ThreeVector mom = (*i)->getMomentum();
1105 eventInfo->px[eventInfo->nParticles] = 1105 eventInfo->px[eventInfo->nParticles] = mom.getX();
1106 eventInfo->py[eventInfo->nParticles] = 1106 eventInfo->py[eventInfo->nParticles] = mom.getY();
1107 eventInfo->pz[eventInfo->nParticles] = 1107 eventInfo->pz[eventInfo->nParticles] = mom.getZ();
1108 eventInfo->theta[eventInfo->nParticles] 1108 eventInfo->theta[eventInfo->nParticles] = Math::toDegrees(mom.theta());
1109 eventInfo->phi[eventInfo->nParticles] = 1109 eventInfo->phi[eventInfo->nParticles] = Math::toDegrees(mom.phi());
1110 eventInfo->origin[eventInfo->nParticles 1110 eventInfo->origin[eventInfo->nParticles] = -1;
1111 #ifdef INCLXX_IN_GEANT4_MODE 1111 #ifdef INCLXX_IN_GEANT4_MODE
1112 eventInfo->parentResonancePDGCode[event 1112 eventInfo->parentResonancePDGCode[eventInfo->nParticles] = (*i)->getParentResonancePDGCode();
1113 eventInfo->parentResonanceID[eventInfo- 1113 eventInfo->parentResonanceID[eventInfo->nParticles] = (*i)->getParentResonanceID();
1114 #endif 1114 #endif
1115 eventInfo->history.push_back(""); 1115 eventInfo->history.push_back("");
1116 if ((*i)->getType() != Composite) { 1116 if ((*i)->getType() != Composite) {
1117 ParticleSpecies pt((*i)->getType()); 1117 ParticleSpecies pt((*i)->getType());
1118 eventInfo->PDGCode[eventInfo->nPartic 1118 eventInfo->PDGCode[eventInfo->nParticles] = pt.getPDGCode();
1119 } 1119 }
1120 else { 1120 else {
1121 ParticleSpecies pt((*i)->getA(), (*i) 1121 ParticleSpecies pt((*i)->getA(), (*i)->getZ(), (*i)->getS());
1122 eventInfo->PDGCode[eventInfo->nPartic 1122 eventInfo->PDGCode[eventInfo->nParticles] = pt.getPDGCode();
1123 } 1123 }
1124 eventInfo->nParticles++; 1124 eventInfo->nParticles++;
1125 } 1125 }
1126 eventInfo->nucleonAbsorption = isNucleonA 1126 eventInfo->nucleonAbsorption = isNucleonAbsorption;
1127 eventInfo->pionAbsorption = isPionAbsorpt 1127 eventInfo->pionAbsorption = isPionAbsorption;
1128 eventInfo->nCascadeParticles = eventInfo- 1128 eventInfo->nCascadeParticles = eventInfo->nParticles;
1129 1129
1130 // Projectile-like remnant characteristic 1130 // Projectile-like remnant characteristics
1131 if(theProjectileRemnant && theProjectileR 1131 if(theProjectileRemnant && theProjectileRemnant->getA()>0) {
1132 #ifdef INCLXX_IN_GEANT4_MODE 1132 #ifdef INCLXX_IN_GEANT4_MODE
1133 eventInfo->ARem[eventInfo->nRemnants] = 1133 eventInfo->ARem[eventInfo->nRemnants] = (G4INCL::Short_t)theProjectileRemnant->getA();
1134 eventInfo->ZRem[eventInfo->nRemnants] = 1134 eventInfo->ZRem[eventInfo->nRemnants] = (G4INCL::Short_t)theProjectileRemnant->getZ();
1135 eventInfo->SRem[eventInfo->nRemnants] = 1135 eventInfo->SRem[eventInfo->nRemnants] = (G4INCL::Short_t)theProjectileRemnant->getS();
1136 #else 1136 #else
1137 eventInfo->ARem[eventInfo->nRemnants] = 1137 eventInfo->ARem[eventInfo->nRemnants] = (Short_t)theProjectileRemnant->getA();
1138 eventInfo->ZRem[eventInfo->nRemnants] = 1138 eventInfo->ZRem[eventInfo->nRemnants] = (Short_t)theProjectileRemnant->getZ();
1139 eventInfo->SRem[eventInfo->nRemnants] = 1139 eventInfo->SRem[eventInfo->nRemnants] = (Short_t)theProjectileRemnant->getS();
1140 #endif 1140 #endif
1141 G4double eStar = theProjectileRemnant-> 1141 G4double eStar = theProjectileRemnant->getExcitationEnergy();
1142 if(std::abs(eStar)<1E-10) 1142 if(std::abs(eStar)<1E-10)
1143 eStar = 0.0; // blame rounding and se 1143 eStar = 0.0; // blame rounding and set the excitation energy to zero
1144 eventInfo->EStarRem[eventInfo->nRemnant 1144 eventInfo->EStarRem[eventInfo->nRemnants] = eStar;
1145 if(eventInfo->EStarRem[eventInfo->nRemn 1145 if(eventInfo->EStarRem[eventInfo->nRemnants]<0.) {
1146 INCL_WARN("Negative excitation energy in 1146 INCL_WARN("Negative excitation energy in projectile-like remnant! EStarRem = " << eventInfo->EStarRem[eventInfo->nRemnants] << '\n');
1147 } 1147 }
1148 const ThreeVector &spin = theProjectile 1148 const ThreeVector &spin = theProjectileRemnant->getSpin();
1149 if(eventInfo->ARem[eventInfo->nRemnants 1149 if(eventInfo->ARem[eventInfo->nRemnants]%2==0) { // even-A nucleus
1150 eventInfo->JRem[eventInfo->nRemnants] = ( 1150 eventInfo->JRem[eventInfo->nRemnants] = (G4int) (spin.mag()/PhysicalConstants::hc + 0.5);
1151 } else { // odd-A nucleus 1151 } else { // odd-A nucleus
1152 eventInfo->JRem[eventInfo->nRemnants] = ( 1152 eventInfo->JRem[eventInfo->nRemnants] = ((G4int) (spin.mag()/PhysicalConstants::hc)) + 0.5;
1153 } 1153 }
1154 eventInfo->EKinRem[eventInfo->nRemnants 1154 eventInfo->EKinRem[eventInfo->nRemnants] = theProjectileRemnant->getKineticEnergy();
1155 const ThreeVector &mom = theProjectileR 1155 const ThreeVector &mom = theProjectileRemnant->getMomentum();
1156 eventInfo->pxRem[eventInfo->nRemnants] 1156 eventInfo->pxRem[eventInfo->nRemnants] = mom.getX();
1157 eventInfo->pyRem[eventInfo->nRemnants] 1157 eventInfo->pyRem[eventInfo->nRemnants] = mom.getY();
1158 eventInfo->pzRem[eventInfo->nRemnants] 1158 eventInfo->pzRem[eventInfo->nRemnants] = mom.getZ();
1159 eventInfo->jxRem[eventInfo->nRemnants] 1159 eventInfo->jxRem[eventInfo->nRemnants] = spin.getX() / PhysicalConstants::hc;
1160 eventInfo->jyRem[eventInfo->nRemnants] 1160 eventInfo->jyRem[eventInfo->nRemnants] = spin.getY() / PhysicalConstants::hc;
1161 eventInfo->jzRem[eventInfo->nRemnants] 1161 eventInfo->jzRem[eventInfo->nRemnants] = spin.getZ() / PhysicalConstants::hc;
1162 eventInfo->thetaRem[eventInfo->nRemnant 1162 eventInfo->thetaRem[eventInfo->nRemnants] = Math::toDegrees(mom.theta());
1163 eventInfo->phiRem[eventInfo->nRemnants] 1163 eventInfo->phiRem[eventInfo->nRemnants] = Math::toDegrees(mom.phi());
1164 eventInfo->nRemnants++; 1164 eventInfo->nRemnants++;
1165 } 1165 }
1166 1166
1167 // Target-like remnant characteristics 1167 // Target-like remnant characteristics
1168 if(hasRemnant()) { 1168 if(hasRemnant()) {
1169 #ifdef INCLXX_IN_GEANT4_MODE 1169 #ifdef INCLXX_IN_GEANT4_MODE
1170 eventInfo->ARem[eventInfo->nRemnants] = 1170 eventInfo->ARem[eventInfo->nRemnants] = (G4INCL::Short_t)getA();
1171 eventInfo->ZRem[eventInfo->nRemnants] = 1171 eventInfo->ZRem[eventInfo->nRemnants] = (G4INCL::Short_t)getZ();
1172 eventInfo->SRem[eventInfo->nRemnants] = 1172 eventInfo->SRem[eventInfo->nRemnants] = (G4INCL::Short_t)getS();
1173 #else 1173 #else
1174 eventInfo->ARem[eventInfo->nRemnants] = 1174 eventInfo->ARem[eventInfo->nRemnants] = (Short_t)getA();
1175 eventInfo->ZRem[eventInfo->nRemnants] = 1175 eventInfo->ZRem[eventInfo->nRemnants] = (Short_t)getZ();
1176 eventInfo->SRem[eventInfo->nRemnants] = 1176 eventInfo->SRem[eventInfo->nRemnants] = (Short_t)getS();
1177 #endif 1177 #endif
1178 eventInfo->EStarRem[eventInfo->nRemnant 1178 eventInfo->EStarRem[eventInfo->nRemnants] = getExcitationEnergy();
1179 if(eventInfo->EStarRem[eventInfo->nRemn 1179 if(eventInfo->EStarRem[eventInfo->nRemnants]<0.) {
1180 INCL_WARN("Negative excitation energy in 1180 INCL_WARN("Negative excitation energy in target-like remnant! EStarRem = " << eventInfo->EStarRem[eventInfo->nRemnants] << " eventNumber=" << eventInfo->eventNumber << '\n');
1181 } 1181 }
1182 const ThreeVector &spin = getSpin(); 1182 const ThreeVector &spin = getSpin();
1183 if(eventInfo->ARem[eventInfo->nRemnants 1183 if(eventInfo->ARem[eventInfo->nRemnants]%2==0) { // even-A nucleus
1184 eventInfo->JRem[eventInfo->nRemnants] = ( 1184 eventInfo->JRem[eventInfo->nRemnants] = (G4int) (spin.mag()/PhysicalConstants::hc + 0.5);
1185 } else { // odd-A nucleus 1185 } else { // odd-A nucleus
1186 eventInfo->JRem[eventInfo->nRemnants] = ( 1186 eventInfo->JRem[eventInfo->nRemnants] = ((G4int) (spin.mag()/PhysicalConstants::hc)) + 0.5;
1187 } 1187 }
1188 eventInfo->EKinRem[eventInfo->nRemnants 1188 eventInfo->EKinRem[eventInfo->nRemnants] = getKineticEnergy();
1189 const ThreeVector &mom = getMomentum(); 1189 const ThreeVector &mom = getMomentum();
1190 eventInfo->pxRem[eventInfo->nRemnants] 1190 eventInfo->pxRem[eventInfo->nRemnants] = mom.getX();
1191 eventInfo->pyRem[eventInfo->nRemnants] 1191 eventInfo->pyRem[eventInfo->nRemnants] = mom.getY();
1192 eventInfo->pzRem[eventInfo->nRemnants] 1192 eventInfo->pzRem[eventInfo->nRemnants] = mom.getZ();
1193 eventInfo->jxRem[eventInfo->nRemnants] 1193 eventInfo->jxRem[eventInfo->nRemnants] = spin.getX() / PhysicalConstants::hc;
1194 eventInfo->jyRem[eventInfo->nRemnants] 1194 eventInfo->jyRem[eventInfo->nRemnants] = spin.getY() / PhysicalConstants::hc;
1195 eventInfo->jzRem[eventInfo->nRemnants] 1195 eventInfo->jzRem[eventInfo->nRemnants] = spin.getZ() / PhysicalConstants::hc;
1196 eventInfo->thetaRem[eventInfo->nRemnant 1196 eventInfo->thetaRem[eventInfo->nRemnants] = Math::toDegrees(mom.theta());
1197 eventInfo->phiRem[eventInfo->nRemnants] 1197 eventInfo->phiRem[eventInfo->nRemnants] = Math::toDegrees(mom.phi());
1198 eventInfo->nRemnants++; 1198 eventInfo->nRemnants++;
1199 } 1199 }
1200 1200
1201 // Global counters, flags, etc. 1201 // Global counters, flags, etc.
1202 Book const &theBook = theStore->getBook() 1202 Book const &theBook = theStore->getBook();
1203 eventInfo->nCollisions = theBook.getAccep 1203 eventInfo->nCollisions = theBook.getAcceptedCollisions();
1204 eventInfo->nBlockedCollisions = theBook.g 1204 eventInfo->nBlockedCollisions = theBook.getBlockedCollisions();
1205 eventInfo->nDecays = theBook.getAcceptedD 1205 eventInfo->nDecays = theBook.getAcceptedDecays();
1206 eventInfo->nBlockedDecays = theBook.getBl 1206 eventInfo->nBlockedDecays = theBook.getBlockedDecays();
1207 eventInfo->firstCollisionTime = theBook.g 1207 eventInfo->firstCollisionTime = theBook.getFirstCollisionTime();
1208 eventInfo->firstCollisionXSec = theBook.g 1208 eventInfo->firstCollisionXSec = theBook.getFirstCollisionXSec();
1209 eventInfo->firstCollisionSpectatorPositio 1209 eventInfo->firstCollisionSpectatorPosition = theBook.getFirstCollisionSpectatorPosition();
1210 eventInfo->firstCollisionSpectatorMomentu 1210 eventInfo->firstCollisionSpectatorMomentum = theBook.getFirstCollisionSpectatorMomentum();
1211 eventInfo->firstCollisionIsElastic = theB 1211 eventInfo->firstCollisionIsElastic = theBook.getFirstCollisionIsElastic();
1212 eventInfo->nReflectionAvatars = theBook.g 1212 eventInfo->nReflectionAvatars = theBook.getAvatars(SurfaceAvatarType);
1213 eventInfo->nCollisionAvatars = theBook.ge 1213 eventInfo->nCollisionAvatars = theBook.getAvatars(CollisionAvatarType);
1214 eventInfo->nDecayAvatars = theBook.getAva 1214 eventInfo->nDecayAvatars = theBook.getAvatars(DecayAvatarType);
1215 eventInfo->nEnergyViolationInteraction = 1215 eventInfo->nEnergyViolationInteraction = theBook.getEnergyViolationInteraction();
1216 } 1216 }
1217 1217
1218 1218
1219 1219
1220 Nucleus::ConservationBalance Nucleus::getCo 1220 Nucleus::ConservationBalance Nucleus::getConservationBalance(const EventInfo &theEventInfo, const G4bool afterRecoil) const {
1221 ConservationBalance theBalance; 1221 ConservationBalance theBalance;
1222 // Initialise balance variables with the 1222 // Initialise balance variables with the incoming values
1223 INCL_DEBUG("theEventInfo " << theEventInf 1223 INCL_DEBUG("theEventInfo " << theEventInfo.Zt << " " << theEventInfo.At << '\n');
1224 theBalance.Z = theEventInfo.Zp + theEvent 1224 theBalance.Z = theEventInfo.Zp + theEventInfo.Zt;
1225 theBalance.A = theEventInfo.Ap + theEvent 1225 theBalance.A = theEventInfo.Ap + theEventInfo.At;
1226 theBalance.S = theEventInfo.Sp + theEvent 1226 theBalance.S = theEventInfo.Sp + theEventInfo.St;
1227 INCL_DEBUG("theBalance Z and A " << theBa 1227 INCL_DEBUG("theBalance Z and A " << theBalance.Z << " " << theBalance.A << '\n');
1228 theBalance.energy = getInitialEnergy(); 1228 theBalance.energy = getInitialEnergy();
1229 theBalance.momentum = getIncomingMomentum 1229 theBalance.momentum = getIncomingMomentum();
1230 1230
1231 // Process outgoing particles 1231 // Process outgoing particles
1232 ParticleList const &outgoingParticles = t 1232 ParticleList const &outgoingParticles = theStore->getOutgoingParticles();
1233 for(ParticleIter i=outgoingParticles.begi 1233 for(ParticleIter i=outgoingParticles.begin(), e=outgoingParticles.end(); i!=e; ++i ) {
1234 theBalance.Z -= (*i)->getZ(); 1234 theBalance.Z -= (*i)->getZ();
1235 theBalance.A -= (*i)->getA(); 1235 theBalance.A -= (*i)->getA();
1236 theBalance.S -= (*i)->getS(); 1236 theBalance.S -= (*i)->getS();
1237 // For outgoing clusters, the total ene 1237 // For outgoing clusters, the total energy automatically includes the
1238 // excitation energy 1238 // excitation energy
1239 theBalance.energy -= (*i)->getEnergy(); 1239 theBalance.energy -= (*i)->getEnergy(); // Note that outgoing particles should have the real mass
1240 theBalance.momentum -= (*i)->getMomentu 1240 theBalance.momentum -= (*i)->getMomentum();
1241 } 1241 }
1242 1242
1243 // Projectile-like remnant contribution, 1243 // Projectile-like remnant contribution, if present
1244 if(theProjectileRemnant && theProjectileR 1244 if(theProjectileRemnant && theProjectileRemnant->getA()>0) {
1245 theBalance.Z -= theProjectileRemnant->g 1245 theBalance.Z -= theProjectileRemnant->getZ();
1246 theBalance.A -= theProjectileRemnant->g 1246 theBalance.A -= theProjectileRemnant->getA();
1247 theBalance.S -= theProjectileRemnant->g 1247 theBalance.S -= theProjectileRemnant->getS();
1248 theBalance.energy -= ParticleTable::get 1248 theBalance.energy -= ParticleTable::getTableMass(theProjectileRemnant->getA(),theProjectileRemnant->getZ(),theProjectileRemnant->getS()) +
1249 theProjectileRemnant->getExcitationEn 1249 theProjectileRemnant->getExcitationEnergy();
1250 theBalance.energy -= theProjectileRemna 1250 theBalance.energy -= theProjectileRemnant->getKineticEnergy();
1251 theBalance.momentum -= theProjectileRem 1251 theBalance.momentum -= theProjectileRemnant->getMomentum();
1252 } 1252 }
1253 1253
1254 // Target-like remnant contribution, if p 1254 // Target-like remnant contribution, if present
1255 if(hasRemnant()) { 1255 if(hasRemnant()) {
1256 theBalance.Z -= getZ(); 1256 theBalance.Z -= getZ();
1257 theBalance.A -= getA(); 1257 theBalance.A -= getA();
1258 theBalance.S -= getS(); 1258 theBalance.S -= getS();
1259 theBalance.energy -= ParticleTable::get 1259 theBalance.energy -= ParticleTable::getTableMass(getA(),getZ(),getS()) +
1260 getExcitationEnergy(); 1260 getExcitationEnergy();
1261 if(afterRecoil) 1261 if(afterRecoil)
1262 theBalance.energy -= getKineticEnergy 1262 theBalance.energy -= getKineticEnergy();
1263 theBalance.momentum -= getMomentum(); 1263 theBalance.momentum -= getMomentum();
1264 } 1264 }
1265 1265
1266 return theBalance; 1266 return theBalance;
1267 } 1267 }
1268 1268
1269 void Nucleus::useFusionKinematics() { 1269 void Nucleus::useFusionKinematics() {
1270 setEnergy(initialEnergy); 1270 setEnergy(initialEnergy);
1271 setMomentum(incomingMomentum); 1271 setMomentum(incomingMomentum);
1272 setSpin(incomingAngularMomentum); 1272 setSpin(incomingAngularMomentum);
1273 theExcitationEnergy = std::sqrt(theEnergy 1273 theExcitationEnergy = std::sqrt(theEnergy*theEnergy-theMomentum.mag2()) - getTableMass();
1274 setMass(getTableMass() + theExcitationEne 1274 setMass(getTableMass() + theExcitationEnergy);
1275 } 1275 }
1276 1276
1277 void Nucleus::finalizeProjectileRemnant(con 1277 void Nucleus::finalizeProjectileRemnant(const G4double anEmissionTime) {
1278 // Deal with the projectile remnant 1278 // Deal with the projectile remnant
1279 const G4int prA = theProjectileRemnant->g 1279 const G4int prA = theProjectileRemnant->getA();
1280 if(prA>=1) { 1280 if(prA>=1) {
1281 // Set the mass 1281 // Set the mass
1282 const G4double aMass = theProjectileRem 1282 const G4double aMass = theProjectileRemnant->getInvariantMass();
1283 theProjectileRemnant->setMass(aMass); 1283 theProjectileRemnant->setMass(aMass);
1284 1284
1285 // Compute the excitation energy from t 1285 // Compute the excitation energy from the invariant mass
1286 const G4double anExcitationEnergy = aMa 1286 const G4double anExcitationEnergy = aMass
1287 - ParticleTable::getTableMass(prA, th 1287 - ParticleTable::getTableMass(prA, theProjectileRemnant->getZ(), theProjectileRemnant->getS());
1288 1288
1289 // Set the excitation energy 1289 // Set the excitation energy
1290 theProjectileRemnant->setExcitationEner 1290 theProjectileRemnant->setExcitationEnergy(anExcitationEnergy);
1291 1291
1292 // No spin! 1292 // No spin!
1293 theProjectileRemnant->setSpin(ThreeVect 1293 theProjectileRemnant->setSpin(ThreeVector());
1294 1294
1295 // Set the emission time 1295 // Set the emission time
1296 theProjectileRemnant->setEmissionTime(a 1296 theProjectileRemnant->setEmissionTime(anEmissionTime);
1297 } 1297 }
1298 } 1298 }
1299 1299
1300 } 1300 }
1301 1301
1302 #endif 1302 #endif
1303 1303