| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // INCL++ intra-nuclear cascade model 26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France 27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium 28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France 29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and H 30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France 31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France 32 // Davide Mancusi, CEA-Saclay, France
33 // 33 //
34 #define INCLXX_IN_GEANT4_MODE 1 34 #define INCLXX_IN_GEANT4_MODE 1
35 35
36 #include "globals.hh" 36 #include "globals.hh"
37 37
38 /** \file G4INCLProjectileRemnant.cc 38 /** \file G4INCLProjectileRemnant.cc
39 * \brief Class for constructing a projectile- 39 * \brief Class for constructing a projectile-like remnant.
40 * 40 *
41 * \date 20 March 2012 41 * \date 20 March 2012
42 * \author Davide Mancusi 42 * \author Davide Mancusi
43 */ 43 */
44 44
45 #include "G4INCLProjectileRemnant.hh" 45 #include "G4INCLProjectileRemnant.hh"
46 #include <algorithm> 46 #include <algorithm>
47 #include <numeric> 47 #include <numeric>
48 48
49 namespace G4INCL { 49 namespace G4INCL {
50 50
51 void ProjectileRemnant::reset() { 51 void ProjectileRemnant::reset() {
52 deleteParticles(); 52 deleteParticles();
53 thePosition = ThreeVector(); 53 thePosition = ThreeVector();
54 theMomentum = ThreeVector(); 54 theMomentum = ThreeVector();
55 theEnergy = 0.0; 55 theEnergy = 0.0;
56 thePotentialEnergy = 0.0; 56 thePotentialEnergy = 0.0;
57 theA = 0; 57 theA = 0;
58 theZ = 0; 58 theZ = 0;
59 nCollisions = 0; 59 nCollisions = 0;
60 60
61 for(std::map<long, Particle*>::const_itera 61 for(std::map<long, Particle*>::const_iterator i=storedComponents.begin(); i!=storedComponents.end(); ++i) {
62 Particle *p = new Particle(*(i->second)) 62 Particle *p = new Particle(*(i->second));
63 EnergyLevelMap::iterator energyIter = th 63 EnergyLevelMap::iterator energyIter = theInitialEnergyLevels.find(i->first);
64 // assert(energyIter!=theInitialEnergyLevels.e 64 // assert(energyIter!=theInitialEnergyLevels.end());
65 const G4double energyLevel = energyIter- 65 const G4double energyLevel = energyIter->second;
66 theInitialEnergyLevels.erase(energyIter) 66 theInitialEnergyLevels.erase(energyIter);
67 theInitialEnergyLevels[p->getID()] = ene 67 theInitialEnergyLevels[p->getID()] = energyLevel;
68 addParticle(p); 68 addParticle(p);
69 } 69 }
70 if(theA>0) 70 if(theA>0)
71 thePosition /= theA; 71 thePosition /= theA;
72 setTableMass(); 72 setTableMass();
73 INCL_DEBUG("ProjectileRemnant object was r 73 INCL_DEBUG("ProjectileRemnant object was reset:" << '\n' << print());
74 } 74 }
75 75
76 void ProjectileRemnant::removeParticle(Parti 76 void ProjectileRemnant::removeParticle(Particle * const p, const G4double theProjectileCorrection) {
77 // assert(p->isNucleon() || p->isLambda()); << 77 // assert(p->isNucleon());
78 78
79 INCL_DEBUG("The following Particle is abou 79 INCL_DEBUG("The following Particle is about to be removed from the ProjectileRemnant:"
80 << '\n' << p->print() 80 << '\n' << p->print()
81 << "theProjectileCorrection=" << thePr 81 << "theProjectileCorrection=" << theProjectileCorrection << '\n');
82 // Update A, Z, S, momentum, and energy of << 82 // Update A, Z, momentum and energy of the projectile remnant
83 theA -= p->getA(); 83 theA -= p->getA();
84 theZ -= p->getZ(); 84 theZ -= p->getZ();
85 theS -= p->getS(); <<
86 85
87 ThreeVector const &oldMomentum = p->getMom 86 ThreeVector const &oldMomentum = p->getMomentum();
88 const G4double oldEnergy = p->getEnergy(); 87 const G4double oldEnergy = p->getEnergy();
89 Cluster::removeParticle(p); 88 Cluster::removeParticle(p);
90 89
91 #if !defined(NDEBUG) && !defined(INCLXX_IN_GEA 90 #if !defined(NDEBUG) && !defined(INCLXX_IN_GEANT4_MODE)
92 ThreeVector theTotalMomentum; 91 ThreeVector theTotalMomentum;
93 G4double theTotalEnergy = 0.; 92 G4double theTotalEnergy = 0.;
94 const G4double theThreshold = 0.1; 93 const G4double theThreshold = 0.1;
95 #endif 94 #endif
96 95
97 if(getA()>0) { // if there are any particl 96 if(getA()>0) { // if there are any particles left
98 // assert((unsigned int)getA()==particles.size 97 // assert((unsigned int)getA()==particles.size());
99 98
100 const G4double theProjectileCorrectionPe 99 const G4double theProjectileCorrectionPerNucleon = theProjectileCorrection / particles.size();
101 100
102 // Update the kinematics of the componen 101 // Update the kinematics of the components
103 for(ParticleIter i=particles.begin(), e= 102 for(ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) {
104 (*i)->setEnergy((*i)->getEnergy() + th 103 (*i)->setEnergy((*i)->getEnergy() + theProjectileCorrectionPerNucleon);
105 (*i)->setMass((*i)->getInvariantMass() 104 (*i)->setMass((*i)->getInvariantMass());
106 #if !defined(NDEBUG) && !defined(INCLXX_IN_GEA 105 #if !defined(NDEBUG) && !defined(INCLXX_IN_GEANT4_MODE)
107 theTotalMomentum += (*i)->getMomentum( 106 theTotalMomentum += (*i)->getMomentum();
108 theTotalEnergy += (*i)->getEnergy(); 107 theTotalEnergy += (*i)->getEnergy();
109 #endif 108 #endif
110 } 109 }
111 } 110 }
112 111
113 theMomentum -= oldMomentum; 112 theMomentum -= oldMomentum;
114 theEnergy -= oldEnergy - theProjectileCorr 113 theEnergy -= oldEnergy - theProjectileCorrection;
115 114
116 // assert(std::abs((theTotalMomentum-theMoment 115 // assert(std::abs((theTotalMomentum-theMomentum).mag())<theThreshold);
117 // assert(std::abs(theTotalEnergy-theEnergy)<t 116 // assert(std::abs(theTotalEnergy-theEnergy)<theThreshold);
118 INCL_DEBUG("After Particle removal, the Pr 117 INCL_DEBUG("After Particle removal, the ProjectileRemnant looks like this:"
119 << '\n' << print()); 118 << '\n' << print());
120 } 119 }
121 120
122 ParticleList ProjectileRemnant::addDynamical 121 ParticleList ProjectileRemnant::addDynamicalSpectators(ParticleList pL) {
123 // Try as hard as possible to add back all 122 // Try as hard as possible to add back all the dynamical spectators.
124 // Don't add spectators that lead to negat 123 // Don't add spectators that lead to negative excitation energies, but
125 // iterate over the spectators as many tim 124 // iterate over the spectators as many times as possible, until
126 // absolutely sure that all of them were r 125 // absolutely sure that all of them were rejected.
127 unsigned int accepted; 126 unsigned int accepted;
128 unsigned long loopCounter = 0; 127 unsigned long loopCounter = 0;
129 const unsigned long maxLoopCounter = 10000 128 const unsigned long maxLoopCounter = 10000000;
130 do { 129 do {
131 accepted = 0; 130 accepted = 0;
132 ParticleList toBeAdded = pL; 131 ParticleList toBeAdded = pL;
133 for(ParticleIter p=toBeAdded.begin(), e= 132 for(ParticleIter p=toBeAdded.begin(), e=toBeAdded.end(); p!=e; ++p) {
134 G4bool isAccepted = addDynamicalSpecta 133 G4bool isAccepted = addDynamicalSpectator(*p);
135 if(isAccepted) { 134 if(isAccepted) {
136 pL.remove(*p); 135 pL.remove(*p);
137 accepted++; 136 accepted++;
138 } 137 }
139 } 138 }
140 ++loopCounter; 139 ++loopCounter;
141 } while(loopCounter<maxLoopCounter && acce 140 } while(loopCounter<maxLoopCounter && accepted > 0); /* Loop checking, 10.07.2015, D.Mancusi */
142 return pL; 141 return pL;
143 } 142 }
144 143
145 ParticleList ProjectileRemnant::addAllDynami 144 ParticleList ProjectileRemnant::addAllDynamicalSpectators(ParticleList const &pL) {
146 // Put all the spectators in the projectil 145 // Put all the spectators in the projectile
147 ThreeVector theNewMomentum = theMomentum; 146 ThreeVector theNewMomentum = theMomentum;
148 G4double theNewEnergy = theEnergy; 147 G4double theNewEnergy = theEnergy;
149 G4int theNewA = theA; 148 G4int theNewA = theA;
150 G4int theNewZ = theZ; 149 G4int theNewZ = theZ;
151 G4int theNewS = theS; 150 G4int theNewS = theS;
152 for(ParticleIter p=pL.begin(), e=pL.end(); 151 for(ParticleIter p=pL.begin(), e=pL.end(); p!=e; ++p) {
153 // assert((*p)->isNucleonorLambda()); 152 // assert((*p)->isNucleonorLambda());
154 // Add the initial (off-shell) momentum 153 // Add the initial (off-shell) momentum and energy to the projectile remnant
155 theNewMomentum += getStoredMomentum(*p); 154 theNewMomentum += getStoredMomentum(*p);
156 theNewEnergy += (*p)->getEnergy(); 155 theNewEnergy += (*p)->getEnergy();
157 theNewA += (*p)->getA(); 156 theNewA += (*p)->getA();
158 theNewZ += (*p)->getZ(); 157 theNewZ += (*p)->getZ();
159 theNewS += (*p)->getS(); 158 theNewS += (*p)->getS();
160 } 159 }
161 160
162 // Check that the excitation energy of the 161 // Check that the excitation energy of the new projectile remnant is non-negative
163 const G4double theNewMass = ParticleTable: 162 const G4double theNewMass = ParticleTable::getTableMass(theNewA,theNewZ,theNewS);
164 const G4double theNewExcitationEnergy = co 163 const G4double theNewExcitationEnergy = computeExcitationEnergyWith(pL);
165 const G4double theNewEffectiveMass = theNe 164 const G4double theNewEffectiveMass = theNewMass + theNewExcitationEnergy;
166 165
167 // If this condition is satisfied, there i 166 // If this condition is satisfied, there is no solution. Fall back on the
168 // "most" method 167 // "most" method
169 if(theNewEnergy<theNewEffectiveMass) { 168 if(theNewEnergy<theNewEffectiveMass) {
170 INCL_WARN("Could not add all the dynamic 169 INCL_WARN("Could not add all the dynamical spectators back into the projectile remnant."
171 << " Falling back to the \"most\" m 170 << " Falling back to the \"most\" method." << '\n');
172 return addMostDynamicalSpectators(pL); 171 return addMostDynamicalSpectators(pL);
173 } 172 }
174 173
175 // Add all the participants to the project 174 // Add all the participants to the projectile remnant
176 for(ParticleIter p=pL.begin(), e=pL.end(); 175 for(ParticleIter p=pL.begin(), e=pL.end(); p!=e; ++p) {
177 particles.push_back(*p); 176 particles.push_back(*p);
178 } 177 }
179 178
180 // Rescale the momentum of the projectile 179 // Rescale the momentum of the projectile remnant so that sqrt(s) has the
181 // correct value 180 // correct value
182 const G4double scalingFactorSquared = (the 181 const G4double scalingFactorSquared = (theNewEnergy*theNewEnergy-theNewEffectiveMass*theNewEffectiveMass)/theNewMomentum.mag2();
183 const G4double scalingFactor = std::sqrt(s 182 const G4double scalingFactor = std::sqrt(scalingFactorSquared);
184 INCL_DEBUG("Scaling factor for the project 183 INCL_DEBUG("Scaling factor for the projectile-remnant momentum = " << scalingFactor << '\n');
185 184
186 theA = theNewA; 185 theA = theNewA;
187 theZ = theNewZ; 186 theZ = theNewZ;
188 theS = theNewS; 187 theS = theNewS;
189 theMomentum = theNewMomentum * scalingFact 188 theMomentum = theNewMomentum * scalingFactor;
190 theEnergy = theNewEnergy; 189 theEnergy = theNewEnergy;
191 190
192 return ParticleList(); 191 return ParticleList();
193 } 192 }
194 193
195 ParticleList ProjectileRemnant::addMostDynam 194 ParticleList ProjectileRemnant::addMostDynamicalSpectators(ParticleList pL) {
196 // Try as hard as possible to add back all 195 // Try as hard as possible to add back all the dynamical spectators.
197 // Don't add spectators that lead to negat 196 // Don't add spectators that lead to negative excitation energies. Start by
198 // adding all of them, and repeatedly remo 197 // adding all of them, and repeatedly remove the most troublesome one until
199 // the excitation energy becomes non-negat 198 // the excitation energy becomes non-negative.
200 199
201 // Put all the spectators in the projectil 200 // Put all the spectators in the projectile
202 ThreeVector theNewMomentum = theMomentum; 201 ThreeVector theNewMomentum = theMomentum;
203 G4double theNewEnergy = theEnergy; 202 G4double theNewEnergy = theEnergy;
204 G4int theNewA = theA; 203 G4int theNewA = theA;
205 G4int theNewZ = theZ; 204 G4int theNewZ = theZ;
206 G4int theNewS = theS; 205 G4int theNewS = theS;
207 for(ParticleIter p=pL.begin(), e=pL.end(); 206 for(ParticleIter p=pL.begin(), e=pL.end(); p!=e; ++p) {
208 // assert((*p)->isNucleonorLambda()); 207 // assert((*p)->isNucleonorLambda());
209 // Add the initial (off-shell) momentum 208 // Add the initial (off-shell) momentum and energy to the projectile remnant
210 theNewMomentum += getStoredMomentum(*p); 209 theNewMomentum += getStoredMomentum(*p);
211 theNewEnergy += (*p)->getEnergy(); 210 theNewEnergy += (*p)->getEnergy();
212 theNewA += (*p)->getA(); 211 theNewA += (*p)->getA();
213 theNewZ += (*p)->getZ(); 212 theNewZ += (*p)->getZ();
214 theNewS += (*p)->getS(); 213 theNewS += (*p)->getS();
215 } 214 }
216 215
217 // Check that the excitation energy of the 216 // Check that the excitation energy of the new projectile remnant is non-negative
218 const G4double theNewMass = ParticleTable: 217 const G4double theNewMass = ParticleTable::getTableMass(theNewA,theNewZ,theNewS);
219 const G4double theNewInvariantMassSquared 218 const G4double theNewInvariantMassSquared = theNewEnergy*theNewEnergy-theNewMomentum.mag2();
220 219
221 G4bool positiveExcitationEnergy = false; 220 G4bool positiveExcitationEnergy = false;
222 if(theNewInvariantMassSquared>=0.) { 221 if(theNewInvariantMassSquared>=0.) {
223 const G4double theNewInvariantMass = std 222 const G4double theNewInvariantMass = std::sqrt(theNewInvariantMassSquared);
224 positiveExcitationEnergy = (theNewInvari 223 positiveExcitationEnergy = (theNewInvariantMass-theNewMass>-1.e-5);
225 } 224 }
226 225
227 // Keep removing nucleons from the project 226 // Keep removing nucleons from the projectile remnant until we achieve a
228 // non-negative excitation energy. 227 // non-negative excitation energy.
229 ParticleList rejected; 228 ParticleList rejected;
230 while(!positiveExcitationEnergy && !pL.emp 229 while(!positiveExcitationEnergy && !pL.empty()) { /* Loop checking, 10.07.2015, D.Mancusi */
231 G4double maxExcitationEnergy = -1.E30; 230 G4double maxExcitationEnergy = -1.E30;
232 ParticleMutableIter best = pL.end(); 231 ParticleMutableIter best = pL.end();
233 ThreeVector bestMomentum; 232 ThreeVector bestMomentum;
234 G4double bestEnergy = -1.; 233 G4double bestEnergy = -1.;
235 G4int bestA = -1, bestZ = -1, bestS = 0; 234 G4int bestA = -1, bestZ = -1, bestS = 0;
236 for(ParticleList::iterator p=pL.begin(), 235 for(ParticleList::iterator p=pL.begin(), e=pL.end(); p!=e; ++p) {
237 // Subtract the initial (off-shell) mo 236 // Subtract the initial (off-shell) momentum and energy from the new
238 // projectile remnant 237 // projectile remnant
239 const ThreeVector theNewerMomentum = t 238 const ThreeVector theNewerMomentum = theNewMomentum - getStoredMomentum(*p);
240 const G4double theNewerEnergy = theNew 239 const G4double theNewerEnergy = theNewEnergy - (*p)->getEnergy();
241 const G4int theNewerA = theNewA - (*p) 240 const G4int theNewerA = theNewA - (*p)->getA();
242 const G4int theNewerZ = theNewZ - (*p) 241 const G4int theNewerZ = theNewZ - (*p)->getZ();
243 const G4int theNewerS = theNewS - (*p) 242 const G4int theNewerS = theNewS - (*p)->getS();
244 243
245 const G4double theNewerMass = Particle 244 const G4double theNewerMass = ParticleTable::getTableMass(theNewerA,theNewerZ,theNewerS);
246 const G4double theNewerInvariantMassSq 245 const G4double theNewerInvariantMassSquared = theNewerEnergy*theNewerEnergy-theNewerMomentum.mag2();
247 246
248 if(theNewerInvariantMassSquared>=-1.e- 247 if(theNewerInvariantMassSquared>=-1.e-5) {
249 const G4double theNewerInvariantMass 248 const G4double theNewerInvariantMass = std::sqrt(std::max(0.,theNewerInvariantMassSquared));
250 const G4double theNewerExcitationEne 249 const G4double theNewerExcitationEnergy = ((theNewerA>1) ? theNewerInvariantMass-theNewerMass : 0.);
251 // Pick the nucleon that maximises t 250 // Pick the nucleon that maximises the excitation energy of the
252 // ProjectileRemnant 251 // ProjectileRemnant
253 if(theNewerExcitationEnergy>maxExcit 252 if(theNewerExcitationEnergy>maxExcitationEnergy) {
254 best = p; 253 best = p;
255 maxExcitationEnergy = theNewerExci 254 maxExcitationEnergy = theNewerExcitationEnergy;
256 bestMomentum = theNewerMomentum; 255 bestMomentum = theNewerMomentum;
257 bestEnergy = theNewerEnergy; 256 bestEnergy = theNewerEnergy;
258 bestA = theNewerA; 257 bestA = theNewerA;
259 bestZ = theNewerZ; 258 bestZ = theNewerZ;
260 bestS = theNewerS; 259 bestS = theNewerS;
261 } 260 }
262 } 261 }
263 } 262 }
264 263
265 // If we couldn't even calculate the exc 264 // If we couldn't even calculate the excitation energy, fail miserably
266 if(best==pL.end()) 265 if(best==pL.end())
267 return pL; 266 return pL;
268 267
269 rejected.push_back(*best); 268 rejected.push_back(*best);
270 pL.erase(best); 269 pL.erase(best);
271 theNewMomentum = bestMomentum; 270 theNewMomentum = bestMomentum;
272 theNewEnergy = bestEnergy; 271 theNewEnergy = bestEnergy;
273 theNewA = bestA; 272 theNewA = bestA;
274 theNewZ = bestZ; 273 theNewZ = bestZ;
275 theNewS = bestS; 274 theNewS = bestS;
276 275
277 if(maxExcitationEnergy>0.) { 276 if(maxExcitationEnergy>0.) {
278 // Stop here 277 // Stop here
279 positiveExcitationEnergy = true; 278 positiveExcitationEnergy = true;
280 } 279 }
281 } 280 }
282 281
283 // Add the accepted participants to the pr 282 // Add the accepted participants to the projectile remnant
284 for(ParticleIter p=pL.begin(), e=pL.end(); 283 for(ParticleIter p=pL.begin(), e=pL.end(); p!=e; ++p) {
285 particles.push_back(*p); 284 particles.push_back(*p);
286 } 285 }
287 theA = theNewA; 286 theA = theNewA;
288 theZ = theNewZ; 287 theZ = theNewZ;
289 theS = theNewS; 288 theS = theNewS;
290 theMomentum = theNewMomentum; 289 theMomentum = theNewMomentum;
291 theEnergy = theNewEnergy; 290 theEnergy = theNewEnergy;
292 291
293 return rejected; 292 return rejected;
294 } 293 }
295 294
296 G4bool ProjectileRemnant::addDynamicalSpecta 295 G4bool ProjectileRemnant::addDynamicalSpectator(Particle * const p) {
297 // assert(p->isNucleon()); 296 // assert(p->isNucleon());
298 297
299 // Add the initial (off-shell) momentum an 298 // Add the initial (off-shell) momentum and energy to the projectile remnant
300 ThreeVector const &oldMomentum = getStored 299 ThreeVector const &oldMomentum = getStoredMomentum(p);
301 const ThreeVector theNewMomentum = theMome 300 const ThreeVector theNewMomentum = theMomentum + oldMomentum;
302 const G4double oldEnergy = p->getEnergy(); 301 const G4double oldEnergy = p->getEnergy();
303 const G4double theNewEnergy = theEnergy + 302 const G4double theNewEnergy = theEnergy + oldEnergy;
304 303
305 // Check that the excitation energy of the 304 // Check that the excitation energy of the new projectile remnant is non-negative
306 const G4double theNewMass = ParticleTable: 305 const G4double theNewMass = ParticleTable::getTableMass(theA+p->getA(),theZ+p->getZ(),theS+p->getS());
307 const G4double theNewInvariantMassSquared 306 const G4double theNewInvariantMassSquared = theNewEnergy*theNewEnergy-theNewMomentum.mag2();
308 307
309 if(theNewInvariantMassSquared<0.) 308 if(theNewInvariantMassSquared<0.)
310 return false; 309 return false;
311 310
312 const G4double theNewInvariantMass = std:: 311 const G4double theNewInvariantMass = std::sqrt(theNewInvariantMassSquared);
313 312
314 if(theNewInvariantMass-theNewMass<-1.e-5) 313 if(theNewInvariantMass-theNewMass<-1.e-5)
315 return false; // negative excitation ene 314 return false; // negative excitation energy here
316 315
317 // Add the spectator to the projectile rem 316 // Add the spectator to the projectile remnant
318 theA += p->getA(); 317 theA += p->getA();
319 theZ += p->getZ(); 318 theZ += p->getZ();
320 theMomentum = theNewMomentum; 319 theMomentum = theNewMomentum;
321 theEnergy = theNewEnergy; 320 theEnergy = theNewEnergy;
322 particles.push_back(p); 321 particles.push_back(p);
323 return true; 322 return true;
324 } 323 }
325 324
326 G4double ProjectileRemnant::computeExcitatio 325 G4double ProjectileRemnant::computeExcitationEnergyExcept(const long exceptID) const {
327 const EnergyLevels theEnergyLevels = getPr 326 const EnergyLevels theEnergyLevels = getPresentEnergyLevelsExcept(exceptID);
328 return computeExcitationEnergy(theEnergyLe 327 return computeExcitationEnergy(theEnergyLevels);
329 } 328 }
330 329
331 G4double ProjectileRemnant::computeExcitatio 330 G4double ProjectileRemnant::computeExcitationEnergyWith(const ParticleList &pL) const {
332 const EnergyLevels theEnergyLevels = getPr 331 const EnergyLevels theEnergyLevels = getPresentEnergyLevelsWith(pL);
333 return computeExcitationEnergy(theEnergyLe 332 return computeExcitationEnergy(theEnergyLevels);
334 } 333 }
335 334
336 G4double ProjectileRemnant::computeExcitatio 335 G4double ProjectileRemnant::computeExcitationEnergy(const EnergyLevels &levels) const {
337 // The ground-state energy is the sum of t 336 // The ground-state energy is the sum of the A smallest initial projectile
338 // energies. 337 // energies.
339 // For the last nucleon, return 0 so that 338 // For the last nucleon, return 0 so that the algorithm will just put it on
340 // shell. 339 // shell.
341 const std::size_t theNewA = levels.size(); << 340 const unsigned theNewA = levels.size();
342 // assert(theNewA>0); 341 // assert(theNewA>0);
343 if(theNewA==1) 342 if(theNewA==1)
344 return 0.; 343 return 0.;
345 344
346 const G4double groundState = theGroundStat 345 const G4double groundState = theGroundStateEnergies.at(theNewA-1);
347 346
348 // Compute the sum of the presently occupi 347 // Compute the sum of the presently occupied energy levels
349 const G4double excitedState = std::accumul 348 const G4double excitedState = std::accumulate(
350 levels.cbegin(), << 349 levels.begin(),
351 levels.cend(), << 350 levels.end(),
352 0.); 351 0.);
353 352
354 return excitedState-groundState; 353 return excitedState-groundState;
355 } 354 }
356 355
357 ProjectileRemnant::EnergyLevels ProjectileRe 356 ProjectileRemnant::EnergyLevels ProjectileRemnant::getPresentEnergyLevelsExcept(const long exceptID) const {
358 EnergyLevels theEnergyLevels; 357 EnergyLevels theEnergyLevels;
359 for(ParticleIter p=particles.begin(), e=pa 358 for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
360 if((*p)->getID()!=exceptID) { 359 if((*p)->getID()!=exceptID) {
361 EnergyLevelMap::const_iterator i = the 360 EnergyLevelMap::const_iterator i = theInitialEnergyLevels.find((*p)->getID());
362 // assert(i!=theInitialEnergyLevels.end()); 361 // assert(i!=theInitialEnergyLevels.end());
363 theEnergyLevels.push_back(i->second); 362 theEnergyLevels.push_back(i->second);
364 } 363 }
365 } 364 }
366 // assert(theEnergyLevels.size()==particles.si 365 // assert(theEnergyLevels.size()==particles.size()-1);
367 return theEnergyLevels; 366 return theEnergyLevels;
368 } 367 }
369 368
370 ProjectileRemnant::EnergyLevels ProjectileRe 369 ProjectileRemnant::EnergyLevels ProjectileRemnant::getPresentEnergyLevelsWith(const ParticleList &pL) const {
371 EnergyLevels theEnergyLevels; 370 EnergyLevels theEnergyLevels;
372 for(ParticleIter p=particles.begin(), e=pa 371 for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
373 EnergyLevelMap::const_iterator i = theIn 372 EnergyLevelMap::const_iterator i = theInitialEnergyLevels.find((*p)->getID());
374 // assert(i!=theInitialEnergyLevels.end()); 373 // assert(i!=theInitialEnergyLevels.end());
375 theEnergyLevels.push_back(i->second); 374 theEnergyLevels.push_back(i->second);
376 } 375 }
377 for(ParticleIter p=pL.begin(), e=pL.end(); 376 for(ParticleIter p=pL.begin(), e=pL.end(); p!=e; ++p) {
378 EnergyLevelMap::const_iterator i = theIn 377 EnergyLevelMap::const_iterator i = theInitialEnergyLevels.find((*p)->getID());
379 // assert(i!=theInitialEnergyLevels.end()); 378 // assert(i!=theInitialEnergyLevels.end());
380 theEnergyLevels.push_back(i->second); 379 theEnergyLevels.push_back(i->second);
381 } 380 }
382 381
383 // assert(theEnergyLevels.size()==particles.si 382 // assert(theEnergyLevels.size()==particles.size()+pL.size());
384 return theEnergyLevels; 383 return theEnergyLevels;
385 } 384 }
386 385
387 } 386 }
388 387
389 388