Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLProjectileRemnant.cc
Go to the documentation of this file.
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  for(ParticleIter p=pL.begin(), e=pL.end(); p!=e; ++p) {
151 // assert((*p)->isNucleon());
152  // Add the initial (off-shell) momentum and energy to the projectile remnant
153  theNewMomentum += getStoredMomentum(*p);
154  theNewEnergy += (*p)->getEnergy();
155  theNewA += (*p)->getA();
156  theNewZ += (*p)->getZ();
157  }
158 
159  // Check that the excitation energy of the new projectile remnant is non-negative
160  const G4double theNewMass = ParticleTable::getTableMass(theNewA,theNewZ);
161  const G4double theNewExcitationEnergy = computeExcitationEnergyWith(pL);
162  const G4double theNewEffectiveMass = theNewMass + theNewExcitationEnergy;
163 
164  // If this condition is satisfied, there is no solution. Fall back on the
165  // "most" method
166  if(theNewEnergy<theNewEffectiveMass) {
167  INCL_WARN("Could not add all the dynamical spectators back into the projectile remnant."
168  << " Falling back to the \"most\" method." << '\n');
169  return addMostDynamicalSpectators(pL);
170  }
171 
172  // Add all the participants to the projectile remnant
173  for(ParticleIter p=pL.begin(), e=pL.end(); p!=e; ++p) {
174  particles.push_back(*p);
175  }
176 
177  // Rescale the momentum of the projectile remnant so that sqrt(s) has the
178  // correct value
179  const G4double scalingFactorSquared = (theNewEnergy*theNewEnergy-theNewEffectiveMass*theNewEffectiveMass)/theNewMomentum.mag2();
180  const G4double scalingFactor = std::sqrt(scalingFactorSquared);
181  INCL_DEBUG("Scaling factor for the projectile-remnant momentum = " << scalingFactor << '\n');
182 
183  theA = theNewA;
184  theZ = theNewZ;
185  theMomentum = theNewMomentum * scalingFactor;
186  theEnergy = theNewEnergy;
187 
188  return ParticleList();
189  }
190 
192  // Try as hard as possible to add back all the dynamical spectators.
193  // Don't add spectators that lead to negative excitation energies. Start by
194  // adding all of them, and repeatedly remove the most troublesome one until
195  // the excitation energy becomes non-negative.
196 
197  // Put all the spectators in the projectile
198  ThreeVector theNewMomentum = theMomentum;
199  G4double theNewEnergy = theEnergy;
200  G4int theNewA = theA;
201  G4int theNewZ = theZ;
202  for(ParticleIter p=pL.begin(), e=pL.end(); p!=e; ++p) {
203 // assert((*p)->isNucleon());
204  // Add the initial (off-shell) momentum and energy to the projectile remnant
205  theNewMomentum += getStoredMomentum(*p);
206  theNewEnergy += (*p)->getEnergy();
207  theNewA += (*p)->getA();
208  theNewZ += (*p)->getZ();
209  }
210 
211  // Check that the excitation energy of the new projectile remnant is non-negative
212  const G4double theNewMass = ParticleTable::getTableMass(theNewA,theNewZ);
213  const G4double theNewInvariantMassSquared = theNewEnergy*theNewEnergy-theNewMomentum.mag2();
214 
215  G4bool positiveExcitationEnergy = false;
216  if(theNewInvariantMassSquared>=0.) {
217  const G4double theNewInvariantMass = std::sqrt(theNewInvariantMassSquared);
218  positiveExcitationEnergy = (theNewInvariantMass-theNewMass>-1.e-5);
219  }
220 
221  // Keep removing nucleons from the projectile remnant until we achieve a
222  // non-negative excitation energy.
223  ParticleList rejected;
224  while(!positiveExcitationEnergy && !pL.empty()) { /* Loop checking, 10.07.2015, D.Mancusi */
225  G4double maxExcitationEnergy = -1.E30;
226  ParticleMutableIter best = pL.end();
227  ThreeVector bestMomentum;
228  G4double bestEnergy = -1.;
229  G4int bestA = -1, bestZ = -1;
230  for(ParticleList::iterator p=pL.begin(), e=pL.end(); p!=e; ++p) {
231  // Subtract the initial (off-shell) momentum and energy from the new
232  // projectile remnant
233  const ThreeVector theNewerMomentum = theNewMomentum - getStoredMomentum(*p);
234  const G4double theNewerEnergy = theNewEnergy - (*p)->getEnergy();
235  const G4int theNewerA = theNewA - (*p)->getA();
236  const G4int theNewerZ = theNewZ - (*p)->getZ();
237 
238  const G4double theNewerMass = ParticleTable::getTableMass(theNewerA,theNewerZ);
239  const G4double theNewerInvariantMassSquared = theNewerEnergy*theNewerEnergy-theNewerMomentum.mag2();
240 
241  if(theNewerInvariantMassSquared>=-1.e-5) {
242  const G4double theNewerInvariantMass = std::sqrt(std::max(0.,theNewerInvariantMassSquared));
243  const G4double theNewerExcitationEnergy = ((theNewerA>1) ? theNewerInvariantMass-theNewerMass : 0.);
244  // Pick the nucleon that maximises the excitation energy of the
245  // ProjectileRemnant
246  if(theNewerExcitationEnergy>maxExcitationEnergy) {
247  best = p;
248  maxExcitationEnergy = theNewerExcitationEnergy;
249  bestMomentum = theNewerMomentum;
250  bestEnergy = theNewerEnergy;
251  bestA = theNewerA;
252  bestZ = theNewerZ;
253  }
254  }
255  }
256 
257  // If we couldn't even calculate the excitation energy, fail miserably
258  if(best==pL.end())
259  return pL;
260 
261  rejected.push_back(*best);
262  pL.erase(best);
263  theNewMomentum = bestMomentum;
264  theNewEnergy = bestEnergy;
265  theNewA = bestA;
266  theNewZ = bestZ;
267 
268  if(maxExcitationEnergy>0.) {
269  // Stop here
270  positiveExcitationEnergy = true;
271  }
272  }
273 
274  // Add the accepted participants to the projectile remnant
275  for(ParticleIter p=pL.begin(), e=pL.end(); p!=e; ++p) {
276  particles.push_back(*p);
277  }
278  theA = theNewA;
279  theZ = theNewZ;
280  theMomentum = theNewMomentum;
281  theEnergy = theNewEnergy;
282 
283  return rejected;
284  }
285 
286  G4bool ProjectileRemnant::addDynamicalSpectator(Particle * const p) {
287 // assert(p->isNucleon());
288 
289  // Add the initial (off-shell) momentum and energy to the projectile remnant
290  ThreeVector const &oldMomentum = getStoredMomentum(p);
291  const ThreeVector theNewMomentum = theMomentum + oldMomentum;
292  const G4double oldEnergy = p->getEnergy();
293  const G4double theNewEnergy = theEnergy + oldEnergy;
294 
295  // Check that the excitation energy of the new projectile remnant is non-negative
296  const G4double theNewMass = ParticleTable::getTableMass(theA+p->getA(),theZ+p->getZ());
297  const G4double theNewInvariantMassSquared = theNewEnergy*theNewEnergy-theNewMomentum.mag2();
298 
299  if(theNewInvariantMassSquared<0.)
300  return false;
301 
302  const G4double theNewInvariantMass = std::sqrt(theNewInvariantMassSquared);
303 
304  if(theNewInvariantMass-theNewMass<-1.e-5)
305  return false; // negative excitation energy here
306 
307  // Add the spectator to the projectile remnant
308  theA += p->getA();
309  theZ += p->getZ();
310  theMomentum = theNewMomentum;
311  theEnergy = theNewEnergy;
312  particles.push_back(p);
313  return true;
314  }
315 
317  const EnergyLevels theEnergyLevels = getPresentEnergyLevelsExcept(exceptID);
318  return computeExcitationEnergy(theEnergyLevels);
319  }
320 
322  const EnergyLevels theEnergyLevels = getPresentEnergyLevelsWith(pL);
323  return computeExcitationEnergy(theEnergyLevels);
324  }
325 
326  G4double ProjectileRemnant::computeExcitationEnergy(const EnergyLevels &levels) const {
327  // The ground-state energy is the sum of the A smallest initial projectile
328  // energies.
329  // For the last nucleon, return 0 so that the algorithm will just put it on
330  // shell.
331  const unsigned theNewA = levels.size();
332 // assert(theNewA>0);
333  if(theNewA==1)
334  return 0.;
335 
336  const G4double groundState = theGroundStateEnergies.at(theNewA-1);
337 
338  // Compute the sum of the presently occupied energy levels
339  const G4double excitedState = std::accumulate(
340  levels.begin(),
341  levels.end(),
342  0.);
343 
344  return excitedState-groundState;
345  }
346 
347  ProjectileRemnant::EnergyLevels ProjectileRemnant::getPresentEnergyLevelsExcept(const long exceptID) const {
348  EnergyLevels theEnergyLevels;
349  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
350  if((*p)->getID()!=exceptID) {
351  EnergyLevelMap::const_iterator i = theInitialEnergyLevels.find((*p)->getID());
352 // assert(i!=theInitialEnergyLevels.end());
353  theEnergyLevels.push_back(i->second);
354  }
355  }
356 // assert(theEnergyLevels.size()==particles.size()-1);
357  return theEnergyLevels;
358  }
359 
360  ProjectileRemnant::EnergyLevels ProjectileRemnant::getPresentEnergyLevelsWith(const ParticleList &pL) const {
361  EnergyLevels theEnergyLevels;
362  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
363  EnergyLevelMap::const_iterator i = theInitialEnergyLevels.find((*p)->getID());
364 // assert(i!=theInitialEnergyLevels.end());
365  theEnergyLevels.push_back(i->second);
366  }
367  for(ParticleIter p=pL.begin(), e=pL.end(); p!=e; ++p) {
368  EnergyLevelMap::const_iterator i = theInitialEnergyLevels.find((*p)->getID());
369 // assert(i!=theInitialEnergyLevels.end());
370  theEnergyLevels.push_back(i->second);
371  }
372 
373 // assert(theEnergyLevels.size()==particles.size()+pL.size());
374  return theEnergyLevels;
375  }
376 
377 }
378 
G4int getA() const
Returns the baryon number.
const char * p
Definition: xmltok.h:285
const G4INCL::ThreeVector & getMomentum() const
Class for constructing a projectile-like remnant.
std::string print() const
G4double getEnergy() const
#define INCL_WARN(x)
void reset()
Reset the projectile remnant to the state at the beginning of the cascade.
int G4int
Definition: G4Types.hh:78
G4double mag2() const
G4double computeExcitationEnergyWith(const ParticleList &pL) const
Compute the excitation energy if some nucleons are put back.
std::string print() const
bool G4bool
Definition: G4Types.hh:79
ParticleList particles
G4int getZ() const
Returns the charge number.
G4double computeExcitationEnergyExcept(const long exceptID) const
Compute the excitation energy when a nucleon is removed.
G4INCL::ThreeVector thePosition
void removeParticle(Particle *const p)
Remove a particle from the cluster components.
ParticleList addDynamicalSpectators(ParticleList pL)
Add back dynamical spectators to the projectile remnant.
ParticleList addAllDynamicalSpectators(ParticleList const &pL)
Add back all dynamical spectators to the projectile remnant.
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4ThreadLocal NuclearMassFn getTableMass
Static pointer to the mass function for nuclei.
void setTableMass()
Set the mass of the Particle to its table mass.
ParticleList addMostDynamicalSpectators(ParticleList pL)
Add back dynamical spectators to the projectile remnant.
void addParticle(Particle *const p)
std::vector< G4double > EnergyLevels
void removeParticle(Particle *const p, const G4double theProjectileCorrection)
Remove a nucleon from the projectile remnant.
G4INCL::ThreeVector theMomentum
double G4double
Definition: G4Types.hh:76
#define INCL_DEBUG(x)
ParticleList::iterator ParticleMutableIter
G4double thePotentialEnergy
ParticleList::const_iterator ParticleIter
long getID() const