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