| 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() || p->isLambda());
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, S, 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(); 85 theS -= p->getS();
86 86
87 ThreeVector const &oldMomentum = p->getMom 87 ThreeVector const &oldMomentum = p->getMomentum();
88 const G4double oldEnergy = p->getEnergy(); 88 const G4double oldEnergy = p->getEnergy();
89 Cluster::removeParticle(p); 89 Cluster::removeParticle(p);
90 90
91 #if !defined(NDEBUG) && !defined(INCLXX_IN_GEA 91 #if !defined(NDEBUG) && !defined(INCLXX_IN_GEANT4_MODE)
92 ThreeVector theTotalMomentum; 92 ThreeVector theTotalMomentum;
93 G4double theTotalEnergy = 0.; 93 G4double theTotalEnergy = 0.;
94 const G4double theThreshold = 0.1; 94 const G4double theThreshold = 0.1;
95 #endif 95 #endif
96 96
97 if(getA()>0) { // if there are any particl 97 if(getA()>0) { // if there are any particles left
98 // assert((unsigned int)getA()==particles.size 98 // assert((unsigned int)getA()==particles.size());
99 99
100 const G4double theProjectileCorrectionPe 100 const G4double theProjectileCorrectionPerNucleon = theProjectileCorrection / particles.size();
101 101
102 // Update the kinematics of the componen 102 // Update the kinematics of the components
103 for(ParticleIter i=particles.begin(), e= 103 for(ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) {
104 (*i)->setEnergy((*i)->getEnergy() + th 104 (*i)->setEnergy((*i)->getEnergy() + theProjectileCorrectionPerNucleon);
105 (*i)->setMass((*i)->getInvariantMass() 105 (*i)->setMass((*i)->getInvariantMass());
106 #if !defined(NDEBUG) && !defined(INCLXX_IN_GEA 106 #if !defined(NDEBUG) && !defined(INCLXX_IN_GEANT4_MODE)
107 theTotalMomentum += (*i)->getMomentum( 107 theTotalMomentum += (*i)->getMomentum();
108 theTotalEnergy += (*i)->getEnergy(); 108 theTotalEnergy += (*i)->getEnergy();
109 #endif 109 #endif
110 } 110 }
111 } 111 }
112 112
113 theMomentum -= oldMomentum; 113 theMomentum -= oldMomentum;
114 theEnergy -= oldEnergy - theProjectileCorr 114 theEnergy -= oldEnergy - theProjectileCorrection;
115 115
116 // assert(std::abs((theTotalMomentum-theMoment 116 // assert(std::abs((theTotalMomentum-theMomentum).mag())<theThreshold);
117 // assert(std::abs(theTotalEnergy-theEnergy)<t 117 // assert(std::abs(theTotalEnergy-theEnergy)<theThreshold);
118 INCL_DEBUG("After Particle removal, the Pr 118 INCL_DEBUG("After Particle removal, the ProjectileRemnant looks like this:"
119 << '\n' << print()); 119 << '\n' << print());
120 } 120 }
121 121
122 ParticleList ProjectileRemnant::addDynamical 122 ParticleList ProjectileRemnant::addDynamicalSpectators(ParticleList pL) {
123 // Try as hard as possible to add back all 123 // Try as hard as possible to add back all the dynamical spectators.
124 // Don't add spectators that lead to negat 124 // Don't add spectators that lead to negative excitation energies, but
125 // iterate over the spectators as many tim 125 // iterate over the spectators as many times as possible, until
126 // absolutely sure that all of them were r 126 // absolutely sure that all of them were rejected.
127 unsigned int accepted; 127 unsigned int accepted;
128 unsigned long loopCounter = 0; 128 unsigned long loopCounter = 0;
129 const unsigned long maxLoopCounter = 10000 129 const unsigned long maxLoopCounter = 10000000;
130 do { 130 do {
131 accepted = 0; 131 accepted = 0;
132 ParticleList toBeAdded = pL; 132 ParticleList toBeAdded = pL;
133 for(ParticleIter p=toBeAdded.begin(), e= 133 for(ParticleIter p=toBeAdded.begin(), e=toBeAdded.end(); p!=e; ++p) {
134 G4bool isAccepted = addDynamicalSpecta 134 G4bool isAccepted = addDynamicalSpectator(*p);
135 if(isAccepted) { 135 if(isAccepted) {
136 pL.remove(*p); 136 pL.remove(*p);
137 accepted++; 137 accepted++;
138 } 138 }
139 } 139 }
140 ++loopCounter; 140 ++loopCounter;
141 } while(loopCounter<maxLoopCounter && acce 141 } while(loopCounter<maxLoopCounter && accepted > 0); /* Loop checking, 10.07.2015, D.Mancusi */
142 return pL; 142 return pL;
143 } 143 }
144 144
145 ParticleList ProjectileRemnant::addAllDynami 145 ParticleList ProjectileRemnant::addAllDynamicalSpectators(ParticleList const &pL) {
146 // Put all the spectators in the projectil 146 // Put all the spectators in the projectile
147 ThreeVector theNewMomentum = theMomentum; 147 ThreeVector theNewMomentum = theMomentum;
148 G4double theNewEnergy = theEnergy; 148 G4double theNewEnergy = theEnergy;
149 G4int theNewA = theA; 149 G4int theNewA = theA;
150 G4int theNewZ = theZ; 150 G4int theNewZ = theZ;
151 G4int theNewS = theS; 151 G4int theNewS = theS;
152 for(ParticleIter p=pL.begin(), e=pL.end(); 152 for(ParticleIter p=pL.begin(), e=pL.end(); p!=e; ++p) {
153 // assert((*p)->isNucleonorLambda()); 153 // assert((*p)->isNucleonorLambda());
154 // Add the initial (off-shell) momentum 154 // Add the initial (off-shell) momentum and energy to the projectile remnant
155 theNewMomentum += getStoredMomentum(*p); 155 theNewMomentum += getStoredMomentum(*p);
156 theNewEnergy += (*p)->getEnergy(); 156 theNewEnergy += (*p)->getEnergy();
157 theNewA += (*p)->getA(); 157 theNewA += (*p)->getA();
158 theNewZ += (*p)->getZ(); 158 theNewZ += (*p)->getZ();
159 theNewS += (*p)->getS(); 159 theNewS += (*p)->getS();
160 } 160 }
161 161
162 // Check that the excitation energy of the 162 // Check that the excitation energy of the new projectile remnant is non-negative
163 const G4double theNewMass = ParticleTable: 163 const G4double theNewMass = ParticleTable::getTableMass(theNewA,theNewZ,theNewS);
164 const G4double theNewExcitationEnergy = co 164 const G4double theNewExcitationEnergy = computeExcitationEnergyWith(pL);
165 const G4double theNewEffectiveMass = theNe 165 const G4double theNewEffectiveMass = theNewMass + theNewExcitationEnergy;
166 166
167 // If this condition is satisfied, there i 167 // If this condition is satisfied, there is no solution. Fall back on the
168 // "most" method 168 // "most" method
169 if(theNewEnergy<theNewEffectiveMass) { 169 if(theNewEnergy<theNewEffectiveMass) {
170 INCL_WARN("Could not add all the dynamic 170 INCL_WARN("Could not add all the dynamical spectators back into the projectile remnant."
171 << " Falling back to the \"most\" m 171 << " Falling back to the \"most\" method." << '\n');
172 return addMostDynamicalSpectators(pL); 172 return addMostDynamicalSpectators(pL);
173 } 173 }
174 174
175 // Add all the participants to the project 175 // Add all the participants to the projectile remnant
176 for(ParticleIter p=pL.begin(), e=pL.end(); 176 for(ParticleIter p=pL.begin(), e=pL.end(); p!=e; ++p) {
177 particles.push_back(*p); 177 particles.push_back(*p);
178 } 178 }
179 179
180 // Rescale the momentum of the projectile 180 // Rescale the momentum of the projectile remnant so that sqrt(s) has the
181 // correct value 181 // correct value
182 const G4double scalingFactorSquared = (the 182 const G4double scalingFactorSquared = (theNewEnergy*theNewEnergy-theNewEffectiveMass*theNewEffectiveMass)/theNewMomentum.mag2();
183 const G4double scalingFactor = std::sqrt(s 183 const G4double scalingFactor = std::sqrt(scalingFactorSquared);
184 INCL_DEBUG("Scaling factor for the project 184 INCL_DEBUG("Scaling factor for the projectile-remnant momentum = " << scalingFactor << '\n');
185 185
186 theA = theNewA; 186 theA = theNewA;
187 theZ = theNewZ; 187 theZ = theNewZ;
188 theS = theNewS; 188 theS = theNewS;
189 theMomentum = theNewMomentum * scalingFact 189 theMomentum = theNewMomentum * scalingFactor;
190 theEnergy = theNewEnergy; 190 theEnergy = theNewEnergy;
191 191
192 return ParticleList(); 192 return ParticleList();
193 } 193 }
194 194
195 ParticleList ProjectileRemnant::addMostDynam 195 ParticleList ProjectileRemnant::addMostDynamicalSpectators(ParticleList pL) {
196 // Try as hard as possible to add back all 196 // Try as hard as possible to add back all the dynamical spectators.
197 // Don't add spectators that lead to negat 197 // Don't add spectators that lead to negative excitation energies. Start by
198 // adding all of them, and repeatedly remo 198 // adding all of them, and repeatedly remove the most troublesome one until
199 // the excitation energy becomes non-negat 199 // the excitation energy becomes non-negative.
200 200
201 // Put all the spectators in the projectil 201 // Put all the spectators in the projectile
202 ThreeVector theNewMomentum = theMomentum; 202 ThreeVector theNewMomentum = theMomentum;
203 G4double theNewEnergy = theEnergy; 203 G4double theNewEnergy = theEnergy;
204 G4int theNewA = theA; 204 G4int theNewA = theA;
205 G4int theNewZ = theZ; 205 G4int theNewZ = theZ;
206 G4int theNewS = theS; 206 G4int theNewS = theS;
207 for(ParticleIter p=pL.begin(), e=pL.end(); 207 for(ParticleIter p=pL.begin(), e=pL.end(); p!=e; ++p) {
208 // assert((*p)->isNucleonorLambda()); 208 // assert((*p)->isNucleonorLambda());
209 // Add the initial (off-shell) momentum 209 // Add the initial (off-shell) momentum and energy to the projectile remnant
210 theNewMomentum += getStoredMomentum(*p); 210 theNewMomentum += getStoredMomentum(*p);
211 theNewEnergy += (*p)->getEnergy(); 211 theNewEnergy += (*p)->getEnergy();
212 theNewA += (*p)->getA(); 212 theNewA += (*p)->getA();
213 theNewZ += (*p)->getZ(); 213 theNewZ += (*p)->getZ();
214 theNewS += (*p)->getS(); 214 theNewS += (*p)->getS();
215 } 215 }
216 216
217 // Check that the excitation energy of the 217 // Check that the excitation energy of the new projectile remnant is non-negative
218 const G4double theNewMass = ParticleTable: 218 const G4double theNewMass = ParticleTable::getTableMass(theNewA,theNewZ,theNewS);
219 const G4double theNewInvariantMassSquared 219 const G4double theNewInvariantMassSquared = theNewEnergy*theNewEnergy-theNewMomentum.mag2();
220 220
221 G4bool positiveExcitationEnergy = false; 221 G4bool positiveExcitationEnergy = false;
222 if(theNewInvariantMassSquared>=0.) { 222 if(theNewInvariantMassSquared>=0.) {
223 const G4double theNewInvariantMass = std 223 const G4double theNewInvariantMass = std::sqrt(theNewInvariantMassSquared);
224 positiveExcitationEnergy = (theNewInvari 224 positiveExcitationEnergy = (theNewInvariantMass-theNewMass>-1.e-5);
225 } 225 }
226 226
227 // Keep removing nucleons from the project 227 // Keep removing nucleons from the projectile remnant until we achieve a
228 // non-negative excitation energy. 228 // non-negative excitation energy.
229 ParticleList rejected; 229 ParticleList rejected;
230 while(!positiveExcitationEnergy && !pL.emp 230 while(!positiveExcitationEnergy && !pL.empty()) { /* Loop checking, 10.07.2015, D.Mancusi */
231 G4double maxExcitationEnergy = -1.E30; 231 G4double maxExcitationEnergy = -1.E30;
232 ParticleMutableIter best = pL.end(); 232 ParticleMutableIter best = pL.end();
233 ThreeVector bestMomentum; 233 ThreeVector bestMomentum;
234 G4double bestEnergy = -1.; 234 G4double bestEnergy = -1.;
235 G4int bestA = -1, bestZ = -1, bestS = 0; 235 G4int bestA = -1, bestZ = -1, bestS = 0;
236 for(ParticleList::iterator p=pL.begin(), 236 for(ParticleList::iterator p=pL.begin(), e=pL.end(); p!=e; ++p) {
237 // Subtract the initial (off-shell) mo 237 // Subtract the initial (off-shell) momentum and energy from the new
238 // projectile remnant 238 // projectile remnant
239 const ThreeVector theNewerMomentum = t 239 const ThreeVector theNewerMomentum = theNewMomentum - getStoredMomentum(*p);
240 const G4double theNewerEnergy = theNew 240 const G4double theNewerEnergy = theNewEnergy - (*p)->getEnergy();
241 const G4int theNewerA = theNewA - (*p) 241 const G4int theNewerA = theNewA - (*p)->getA();
242 const G4int theNewerZ = theNewZ - (*p) 242 const G4int theNewerZ = theNewZ - (*p)->getZ();
243 const G4int theNewerS = theNewS - (*p) 243 const G4int theNewerS = theNewS - (*p)->getS();
244 244
245 const G4double theNewerMass = Particle 245 const G4double theNewerMass = ParticleTable::getTableMass(theNewerA,theNewerZ,theNewerS);
246 const G4double theNewerInvariantMassSq 246 const G4double theNewerInvariantMassSquared = theNewerEnergy*theNewerEnergy-theNewerMomentum.mag2();
247 247
248 if(theNewerInvariantMassSquared>=-1.e- 248 if(theNewerInvariantMassSquared>=-1.e-5) {
249 const G4double theNewerInvariantMass 249 const G4double theNewerInvariantMass = std::sqrt(std::max(0.,theNewerInvariantMassSquared));
250 const G4double theNewerExcitationEne 250 const G4double theNewerExcitationEnergy = ((theNewerA>1) ? theNewerInvariantMass-theNewerMass : 0.);
251 // Pick the nucleon that maximises t 251 // Pick the nucleon that maximises the excitation energy of the
252 // ProjectileRemnant 252 // ProjectileRemnant
253 if(theNewerExcitationEnergy>maxExcit 253 if(theNewerExcitationEnergy>maxExcitationEnergy) {
254 best = p; 254 best = p;
255 maxExcitationEnergy = theNewerExci 255 maxExcitationEnergy = theNewerExcitationEnergy;
256 bestMomentum = theNewerMomentum; 256 bestMomentum = theNewerMomentum;
257 bestEnergy = theNewerEnergy; 257 bestEnergy = theNewerEnergy;
258 bestA = theNewerA; 258 bestA = theNewerA;
259 bestZ = theNewerZ; 259 bestZ = theNewerZ;
260 bestS = theNewerS; 260 bestS = theNewerS;
261 } 261 }
262 } 262 }
263 } 263 }
264 264
265 // If we couldn't even calculate the exc 265 // If we couldn't even calculate the excitation energy, fail miserably
266 if(best==pL.end()) 266 if(best==pL.end())
267 return pL; 267 return pL;
268 268
269 rejected.push_back(*best); 269 rejected.push_back(*best);
270 pL.erase(best); 270 pL.erase(best);
271 theNewMomentum = bestMomentum; 271 theNewMomentum = bestMomentum;
272 theNewEnergy = bestEnergy; 272 theNewEnergy = bestEnergy;
273 theNewA = bestA; 273 theNewA = bestA;
274 theNewZ = bestZ; 274 theNewZ = bestZ;
275 theNewS = bestS; 275 theNewS = bestS;
276 276
277 if(maxExcitationEnergy>0.) { 277 if(maxExcitationEnergy>0.) {
278 // Stop here 278 // Stop here
279 positiveExcitationEnergy = true; 279 positiveExcitationEnergy = true;
280 } 280 }
281 } 281 }
282 282
283 // Add the accepted participants to the pr 283 // Add the accepted participants to the projectile remnant
284 for(ParticleIter p=pL.begin(), e=pL.end(); 284 for(ParticleIter p=pL.begin(), e=pL.end(); p!=e; ++p) {
285 particles.push_back(*p); 285 particles.push_back(*p);
286 } 286 }
287 theA = theNewA; 287 theA = theNewA;
288 theZ = theNewZ; 288 theZ = theNewZ;
289 theS = theNewS; 289 theS = theNewS;
290 theMomentum = theNewMomentum; 290 theMomentum = theNewMomentum;
291 theEnergy = theNewEnergy; 291 theEnergy = theNewEnergy;
292 292
293 return rejected; 293 return rejected;
294 } 294 }
295 295
296 G4bool ProjectileRemnant::addDynamicalSpecta 296 G4bool ProjectileRemnant::addDynamicalSpectator(Particle * const p) {
297 // assert(p->isNucleon()); 297 // assert(p->isNucleon());
298 298
299 // Add the initial (off-shell) momentum an 299 // Add the initial (off-shell) momentum and energy to the projectile remnant
300 ThreeVector const &oldMomentum = getStored 300 ThreeVector const &oldMomentum = getStoredMomentum(p);
301 const ThreeVector theNewMomentum = theMome 301 const ThreeVector theNewMomentum = theMomentum + oldMomentum;
302 const G4double oldEnergy = p->getEnergy(); 302 const G4double oldEnergy = p->getEnergy();
303 const G4double theNewEnergy = theEnergy + 303 const G4double theNewEnergy = theEnergy + oldEnergy;
304 304
305 // Check that the excitation energy of the 305 // Check that the excitation energy of the new projectile remnant is non-negative
306 const G4double theNewMass = ParticleTable: 306 const G4double theNewMass = ParticleTable::getTableMass(theA+p->getA(),theZ+p->getZ(),theS+p->getS());
307 const G4double theNewInvariantMassSquared 307 const G4double theNewInvariantMassSquared = theNewEnergy*theNewEnergy-theNewMomentum.mag2();
308 308
309 if(theNewInvariantMassSquared<0.) 309 if(theNewInvariantMassSquared<0.)
310 return false; 310 return false;
311 311
312 const G4double theNewInvariantMass = std:: 312 const G4double theNewInvariantMass = std::sqrt(theNewInvariantMassSquared);
313 313
314 if(theNewInvariantMass-theNewMass<-1.e-5) 314 if(theNewInvariantMass-theNewMass<-1.e-5)
315 return false; // negative excitation ene 315 return false; // negative excitation energy here
316 316
317 // Add the spectator to the projectile rem 317 // Add the spectator to the projectile remnant
318 theA += p->getA(); 318 theA += p->getA();
319 theZ += p->getZ(); 319 theZ += p->getZ();
320 theMomentum = theNewMomentum; 320 theMomentum = theNewMomentum;
321 theEnergy = theNewEnergy; 321 theEnergy = theNewEnergy;
322 particles.push_back(p); 322 particles.push_back(p);
323 return true; 323 return true;
324 } 324 }
325 325
326 G4double ProjectileRemnant::computeExcitatio 326 G4double ProjectileRemnant::computeExcitationEnergyExcept(const long exceptID) const {
327 const EnergyLevels theEnergyLevels = getPr 327 const EnergyLevels theEnergyLevels = getPresentEnergyLevelsExcept(exceptID);
328 return computeExcitationEnergy(theEnergyLe 328 return computeExcitationEnergy(theEnergyLevels);
329 } 329 }
330 330
331 G4double ProjectileRemnant::computeExcitatio 331 G4double ProjectileRemnant::computeExcitationEnergyWith(const ParticleList &pL) const {
332 const EnergyLevels theEnergyLevels = getPr 332 const EnergyLevels theEnergyLevels = getPresentEnergyLevelsWith(pL);
333 return computeExcitationEnergy(theEnergyLe 333 return computeExcitationEnergy(theEnergyLevels);
334 } 334 }
335 335
336 G4double ProjectileRemnant::computeExcitatio 336 G4double ProjectileRemnant::computeExcitationEnergy(const EnergyLevels &levels) const {
337 // The ground-state energy is the sum of t 337 // The ground-state energy is the sum of the A smallest initial projectile
338 // energies. 338 // energies.
339 // For the last nucleon, return 0 so that 339 // For the last nucleon, return 0 so that the algorithm will just put it on
340 // shell. 340 // shell.
341 const std::size_t theNewA = levels.size(); 341 const std::size_t theNewA = levels.size();
342 // assert(theNewA>0); 342 // assert(theNewA>0);
343 if(theNewA==1) 343 if(theNewA==1)
344 return 0.; 344 return 0.;
345 345
346 const G4double groundState = theGroundStat 346 const G4double groundState = theGroundStateEnergies.at(theNewA-1);
347 347
348 // Compute the sum of the presently occupi 348 // Compute the sum of the presently occupied energy levels
349 const G4double excitedState = std::accumul 349 const G4double excitedState = std::accumulate(
350 levels.cbegin(), 350 levels.cbegin(),
351 levels.cend(), 351 levels.cend(),
352 0.); 352 0.);
353 353
354 return excitedState-groundState; 354 return excitedState-groundState;
355 } 355 }
356 356
357 ProjectileRemnant::EnergyLevels ProjectileRe 357 ProjectileRemnant::EnergyLevels ProjectileRemnant::getPresentEnergyLevelsExcept(const long exceptID) const {
358 EnergyLevels theEnergyLevels; 358 EnergyLevels theEnergyLevels;
359 for(ParticleIter p=particles.begin(), e=pa 359 for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
360 if((*p)->getID()!=exceptID) { 360 if((*p)->getID()!=exceptID) {
361 EnergyLevelMap::const_iterator i = the 361 EnergyLevelMap::const_iterator i = theInitialEnergyLevels.find((*p)->getID());
362 // assert(i!=theInitialEnergyLevels.end()); 362 // assert(i!=theInitialEnergyLevels.end());
363 theEnergyLevels.push_back(i->second); 363 theEnergyLevels.push_back(i->second);
364 } 364 }
365 } 365 }
366 // assert(theEnergyLevels.size()==particles.si 366 // assert(theEnergyLevels.size()==particles.size()-1);
367 return theEnergyLevels; 367 return theEnergyLevels;
368 } 368 }
369 369
370 ProjectileRemnant::EnergyLevels ProjectileRe 370 ProjectileRemnant::EnergyLevels ProjectileRemnant::getPresentEnergyLevelsWith(const ParticleList &pL) const {
371 EnergyLevels theEnergyLevels; 371 EnergyLevels theEnergyLevels;
372 for(ParticleIter p=particles.begin(), e=pa 372 for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
373 EnergyLevelMap::const_iterator i = theIn 373 EnergyLevelMap::const_iterator i = theInitialEnergyLevels.find((*p)->getID());
374 // assert(i!=theInitialEnergyLevels.end()); 374 // assert(i!=theInitialEnergyLevels.end());
375 theEnergyLevels.push_back(i->second); 375 theEnergyLevels.push_back(i->second);
376 } 376 }
377 for(ParticleIter p=pL.begin(), e=pL.end(); 377 for(ParticleIter p=pL.begin(), e=pL.end(); p!=e; ++p) {
378 EnergyLevelMap::const_iterator i = theIn 378 EnergyLevelMap::const_iterator i = theInitialEnergyLevels.find((*p)->getID());
379 // assert(i!=theInitialEnergyLevels.end()); 379 // assert(i!=theInitialEnergyLevels.end());
380 theEnergyLevels.push_back(i->second); 380 theEnergyLevels.push_back(i->second);
381 } 381 }
382 382
383 // assert(theEnergyLevels.size()==particles.si 383 // assert(theEnergyLevels.size()==particles.size()+pL.size());
384 return theEnergyLevels; 384 return theEnergyLevels;
385 } 385 }
386 386
387 } 387 }
388 388
389 389