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