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