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