| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 H
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-
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_itera
62 Particle *p = new Particle(*(i->second))
63 EnergyLevelMap::iterator energyIter = th
64 // assert(energyIter!=theInitialEnergyLevels.e
65 const G4double energyLevel = energyIter-
66 theInitialEnergyLevels.erase(energyIter)
67 theInitialEnergyLevels[p->getID()] = ene
68 addParticle(p);
69 }
70 if(theA>0)
71 thePosition /= theA;
72 setTableMass();
73 INCL_DEBUG("ProjectileRemnant object was r
74 }
75
76 void ProjectileRemnant::removeParticle(Parti
77 // assert(p->isNucleon() || p->isLambda());
78
79 INCL_DEBUG("The following Particle is abou
80 << '\n' << p->print()
81 << "theProjectileCorrection=" << thePr
82 // Update A, Z, S, momentum, and energy of
83 theA -= p->getA();
84 theZ -= p->getZ();
85 theS -= p->getS();
86
87 ThreeVector const &oldMomentum = p->getMom
88 const G4double oldEnergy = p->getEnergy();
89 Cluster::removeParticle(p);
90
91 #if !defined(NDEBUG) && !defined(INCLXX_IN_GEA
92 ThreeVector theTotalMomentum;
93 G4double theTotalEnergy = 0.;
94 const G4double theThreshold = 0.1;
95 #endif
96
97 if(getA()>0) { // if there are any particl
98 // assert((unsigned int)getA()==particles.size
99
100 const G4double theProjectileCorrectionPe
101
102 // Update the kinematics of the componen
103 for(ParticleIter i=particles.begin(), e=
104 (*i)->setEnergy((*i)->getEnergy() + th
105 (*i)->setMass((*i)->getInvariantMass()
106 #if !defined(NDEBUG) && !defined(INCLXX_IN_GEA
107 theTotalMomentum += (*i)->getMomentum(
108 theTotalEnergy += (*i)->getEnergy();
109 #endif
110 }
111 }
112
113 theMomentum -= oldMomentum;
114 theEnergy -= oldEnergy - theProjectileCorr
115
116 // assert(std::abs((theTotalMomentum-theMoment
117 // assert(std::abs(theTotalEnergy-theEnergy)<t
118 INCL_DEBUG("After Particle removal, the Pr
119 << '\n' << print());
120 }
121
122 ParticleList ProjectileRemnant::addDynamical
123 // Try as hard as possible to add back all
124 // Don't add spectators that lead to negat
125 // iterate over the spectators as many tim
126 // absolutely sure that all of them were r
127 unsigned int accepted;
128 unsigned long loopCounter = 0;
129 const unsigned long maxLoopCounter = 10000
130 do {
131 accepted = 0;
132 ParticleList toBeAdded = pL;
133 for(ParticleIter p=toBeAdded.begin(), e=
134 G4bool isAccepted = addDynamicalSpecta
135 if(isAccepted) {
136 pL.remove(*p);
137 accepted++;
138 }
139 }
140 ++loopCounter;
141 } while(loopCounter<maxLoopCounter && acce
142 return pL;
143 }
144
145 ParticleList ProjectileRemnant::addAllDynami
146 // Put all the spectators in the projectil
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();
153 // assert((*p)->isNucleonorLambda());
154 // Add the initial (off-shell) momentum
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
163 const G4double theNewMass = ParticleTable:
164 const G4double theNewExcitationEnergy = co
165 const G4double theNewEffectiveMass = theNe
166
167 // If this condition is satisfied, there i
168 // "most" method
169 if(theNewEnergy<theNewEffectiveMass) {
170 INCL_WARN("Could not add all the dynamic
171 << " Falling back to the \"most\" m
172 return addMostDynamicalSpectators(pL);
173 }
174
175 // Add all the participants to the project
176 for(ParticleIter p=pL.begin(), e=pL.end();
177 particles.push_back(*p);
178 }
179
180 // Rescale the momentum of the projectile
181 // correct value
182 const G4double scalingFactorSquared = (the
183 const G4double scalingFactor = std::sqrt(s
184 INCL_DEBUG("Scaling factor for the project
185
186 theA = theNewA;
187 theZ = theNewZ;
188 theS = theNewS;
189 theMomentum = theNewMomentum * scalingFact
190 theEnergy = theNewEnergy;
191
192 return ParticleList();
193 }
194
195 ParticleList ProjectileRemnant::addMostDynam
196 // Try as hard as possible to add back all
197 // Don't add spectators that lead to negat
198 // adding all of them, and repeatedly remo
199 // the excitation energy becomes non-negat
200
201 // Put all the spectators in the projectil
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();
208 // assert((*p)->isNucleonorLambda());
209 // Add the initial (off-shell) momentum
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
218 const G4double theNewMass = ParticleTable:
219 const G4double theNewInvariantMassSquared
220
221 G4bool positiveExcitationEnergy = false;
222 if(theNewInvariantMassSquared>=0.) {
223 const G4double theNewInvariantMass = std
224 positiveExcitationEnergy = (theNewInvari
225 }
226
227 // Keep removing nucleons from the project
228 // non-negative excitation energy.
229 ParticleList rejected;
230 while(!positiveExcitationEnergy && !pL.emp
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(),
237 // Subtract the initial (off-shell) mo
238 // projectile remnant
239 const ThreeVector theNewerMomentum = t
240 const G4double theNewerEnergy = theNew
241 const G4int theNewerA = theNewA - (*p)
242 const G4int theNewerZ = theNewZ - (*p)
243 const G4int theNewerS = theNewS - (*p)
244
245 const G4double theNewerMass = Particle
246 const G4double theNewerInvariantMassSq
247
248 if(theNewerInvariantMassSquared>=-1.e-
249 const G4double theNewerInvariantMass
250 const G4double theNewerExcitationEne
251 // Pick the nucleon that maximises t
252 // ProjectileRemnant
253 if(theNewerExcitationEnergy>maxExcit
254 best = p;
255 maxExcitationEnergy = theNewerExci
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 exc
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 pr
284 for(ParticleIter p=pL.begin(), e=pL.end();
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::addDynamicalSpecta
297 // assert(p->isNucleon());
298
299 // Add the initial (off-shell) momentum an
300 ThreeVector const &oldMomentum = getStored
301 const ThreeVector theNewMomentum = theMome
302 const G4double oldEnergy = p->getEnergy();
303 const G4double theNewEnergy = theEnergy +
304
305 // Check that the excitation energy of the
306 const G4double theNewMass = ParticleTable:
307 const G4double theNewInvariantMassSquared
308
309 if(theNewInvariantMassSquared<0.)
310 return false;
311
312 const G4double theNewInvariantMass = std::
313
314 if(theNewInvariantMass-theNewMass<-1.e-5)
315 return false; // negative excitation ene
316
317 // Add the spectator to the projectile rem
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::computeExcitatio
327 const EnergyLevels theEnergyLevels = getPr
328 return computeExcitationEnergy(theEnergyLe
329 }
330
331 G4double ProjectileRemnant::computeExcitatio
332 const EnergyLevels theEnergyLevels = getPr
333 return computeExcitationEnergy(theEnergyLe
334 }
335
336 G4double ProjectileRemnant::computeExcitatio
337 // The ground-state energy is the sum of t
338 // energies.
339 // For the last nucleon, return 0 so that
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 = theGroundStat
347
348 // Compute the sum of the presently occupi
349 const G4double excitedState = std::accumul
350 levels.cbegin(),
351 levels.cend(),
352 0.);
353
354 return excitedState-groundState;
355 }
356
357 ProjectileRemnant::EnergyLevels ProjectileRe
358 EnergyLevels theEnergyLevels;
359 for(ParticleIter p=particles.begin(), e=pa
360 if((*p)->getID()!=exceptID) {
361 EnergyLevelMap::const_iterator i = the
362 // assert(i!=theInitialEnergyLevels.end());
363 theEnergyLevels.push_back(i->second);
364 }
365 }
366 // assert(theEnergyLevels.size()==particles.si
367 return theEnergyLevels;
368 }
369
370 ProjectileRemnant::EnergyLevels ProjectileRe
371 EnergyLevels theEnergyLevels;
372 for(ParticleIter p=particles.begin(), e=pa
373 EnergyLevelMap::const_iterator i = theIn
374 // assert(i!=theInitialEnergyLevels.end());
375 theEnergyLevels.push_back(i->second);
376 }
377 for(ParticleIter p=pL.begin(), e=pL.end();
378 EnergyLevelMap::const_iterator i = theIn
379 // assert(i!=theInitialEnergyLevels.end());
380 theEnergyLevels.push_back(i->second);
381 }
382
383 // assert(theEnergyLevels.size()==particles.si
384 return theEnergyLevels;
385 }
386
387 }
388
389