| 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.hh
39 * \brief Class for constructing a projectile-like remnant.
40 *
41 * \date 20 March 2012
42 * \author Davide Mancusi
43 */
44
45 #ifndef G4INCLPROJECTILEREMNANT_HH_
46 #define G4INCLPROJECTILEREMNANT_HH_
47
48 #include "G4INCLCluster.hh"
49 #include "G4INCLRandom.hh"
50 #include "G4INCLAllocationPool.hh"
51 #include <vector>
52 #include <map>
53 #include <numeric>
54 #include <functional>
55
56 namespace G4INCL {
57
58 class ProjectileRemnant : public Cluster {
59 public:
60
61 // typedefs for the calculation of the projectile excitation energy
62 typedef std::vector<G4double> EnergyLevels;
63 typedef std::map<long, G4double> EnergyLevelMap;
64
65 ProjectileRemnant(ParticleSpecies const &species, const G4double kineticEnergy)
66 : Cluster(species.theZ, species.theA, species.theS) {
67
68 // Use the table mass
69 setTableMass();
70
71 // Set the kinematics
72 const G4double projectileMass = getMass();
73 const G4double energy = kineticEnergy + projectileMass;
74 const G4double momentumZ = std::sqrt(energy*energy - projectileMass*projectileMass);
75
76 // Initialise the particles
77 initializeParticles();
78 internalBoostToCM();
79 putParticlesOffShell();
80
81 // Store the energy levels of the ProjectileRemnant (used to compute its
82 // excitation energy)
83 storeEnergyLevels();
84
85 // Boost the whole thing
86 const ThreeVector aBoostVector = ThreeVector(0.0, 0.0, momentumZ / energy);
87 boost(-aBoostVector);
88
89 // Freeze the internal motion of the particles
90 freezeInternalMotion();
91
92 // Set as projectile spectator
93 makeProjectileSpectator();
94 }
95
96 ~ProjectileRemnant() {
97 deleteStoredComponents();
98 // The ProjectileRemnant owns its particles
99 deleteParticles();
100 clearEnergyLevels();
101 }
102
103 /// \brief Reset the projectile remnant to the state at the beginning of the cascade
104 void reset();
105
106 /** \brief Remove a nucleon from the projectile remnant
107 *
108 * \param p particle to be removed
109 * \param theProjectileCorrection correction to be given to the projectile total energy
110 */
111 void removeParticle(Particle * const p, const G4double theProjectileCorrection);
112
113 /** \brief Add back dynamical spectators to the projectile remnant
114 *
115 * Try to add the dynamical spectators back to the projectile remnant.
116 * Refuse to do so if this leads to a negative projectile excitation
117 * energy.
118 *
119 * Return a list of rejected dynamical spectators.
120 */
121 ParticleList addDynamicalSpectators(ParticleList pL);
122
123 /** \brief Add back dynamical spectators to the projectile remnant
124 *
125 * Try as hard as possible to add back all the dynamical spectators. Don't
126 * add spectators that lead to negative excitation energies. Start by
127 * adding all of them, and repeatedly remove the most troublesome one until
128 * the excitation energy becomes non-negative.
129 *
130 * Return a list of rejected dynamical spectators.
131 */
132 ParticleList addMostDynamicalSpectators(ParticleList pL);
133
134 /** \brief Add back all dynamical spectators to the projectile remnant
135 *
136 * Return a list of rejected dynamical spectators.
137 */
138 ParticleList addAllDynamicalSpectators(ParticleList const &pL);
139
140 /// \brief Clear the stored projectile components and delete the particles
141 void deleteStoredComponents() {
142 for(std::map<long,Particle*>::const_iterator p=storedComponents.begin(), e=storedComponents.end(); p!=e; ++p)
143 delete p->second;
144 clearStoredComponents();
145 }
146
147 /// \brief Clear the stored projectile components
148 void clearStoredComponents() {
149 storedComponents.clear();
150 }
151
152 /// \brief Clear the stored energy levels
153 void clearEnergyLevels() {
154 theInitialEnergyLevels.clear();
155 theGroundStateEnergies.clear();
156 }
157
158 /** \brief Compute the excitation energy when a nucleon is removed
159 *
160 * Compute the excitation energy of the projectile-like remnant as the
161 * difference between the initial and the present configuration. This
162 * follows the algorithm proposed by A. Boudard in INCL4.2-HI, as
163 * implemented in Geant4.
164 *
165 * \return the excitation energy
166 */
167 G4double computeExcitationEnergyExcept(const long exceptID) const;
168
169 /** \brief Compute the excitation energy if some nucleons are put back
170 *
171 * \return the excitation energy
172 */
173 G4double computeExcitationEnergyWith(const ParticleList &pL) const;
174
175 /// \brief Store the projectile components
176 void storeComponents() {
177 for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
178 // Store the particles (needed for forced CN)
179 storedComponents[(*p)->getID()]=new Particle(**p);
180 }
181 }
182
183 /// \brief Get the number of the stored components
184 G4int getNumberStoredComponents() const {
185 return (G4int)storedComponents.size();
186 }
187
188 /// \brief Store the energy levels
189 void storeEnergyLevels() {
190 EnergyLevels energies;
191
192 for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
193 const G4double theCMEnergy = (*p)->getEnergy();
194 // Store the CM energy in the EnergyLevels map
195 theInitialEnergyLevels[(*p)->getID()] = theCMEnergy;
196 energies.push_back(theCMEnergy);
197 }
198
199 std::sort(energies.begin(), energies.end());
200 // assert(energies.size()==(unsigned int)theA);
201 theGroundStateEnergies.resize(energies.size());
202 // Compute the partial sums of the CM energies -- they are our reference
203 // ground-state energies for any number of nucleons
204 std::partial_sum(energies.begin(), energies.end(), theGroundStateEnergies.begin());
205 }
206
207 EnergyLevels const &getGroundStateEnergies() const {
208 return theGroundStateEnergies;
209 }
210
211 private:
212
213 /** \brief Compute the excitation energy for a given configuration
214 *
215 * The function that does the real job of calculating the excitation energy
216 * for a given configuration of energy levels.
217 *
218 * \param levels a configuration of energy levels
219 * \return the excitation energy
220 */
221 G4double computeExcitationEnergy(const EnergyLevels &levels) const;
222
223 EnergyLevels getPresentEnergyLevelsExcept(const long exceptID) const;
224
225 EnergyLevels getPresentEnergyLevelsWith(const ParticleList &pL) const;
226
227 /// \brief Shuffle the list of stored projectile components
228 ParticleList shuffleStoredComponents() {
229 ParticleList pL = getStoredComponents();
230 std::shuffle(pL.begin(), pL.end(), Random::getAdapter());
231 return pL;
232 }
233
234 ParticleList getStoredComponents() const {
235 ParticleList pL;
236 for(std::map<long,Particle*>::const_iterator p=storedComponents.begin(), e=storedComponents.end(); p!=e; ++p)
237 pL.push_back(p->second);
238 return pL;
239 }
240
241 /// \brief Return the stored momentum of a given projectile component
242 ThreeVector const &getStoredMomentum(Particle const * const p) const {
243 std::map<long,Particle*>::const_iterator i = storedComponents.find(p->getID());
244 if(i==storedComponents.end()) {
245 INCL_ERROR("Couldn't find particle " << p->getID() << " in the list of projectile components" << '\n');
246 return p->getMomentum();
247 } else {
248 return i->second->getMomentum();
249 }
250 }
251
252 /** \brief Add back a nucleon to the projectile remnant
253 *
254 * Try to add a dynamical spectator back to the projectile remnant. Refuse
255 * to do so if this leads to a negative projectile excitation energy.
256 * Return true on success, false on failure.
257 */
258 G4bool addDynamicalSpectator(Particle * const p);
259
260 /// \brief Return the stored energy of a given projectile component
261 /* G4double getStoredEnergy(Particle const * const p) {
262 std::map<long,Particle*>::const_iterator i = initialProjectileComponents.find(p->getID());
263 if(i==initialProjectileComponents.end()) {
264 INCL_ERROR("Couldn't find particle " << p->getID() << " in the list of projectile components" << '\n');
265 return 0.;
266 } else {
267 return i->second->getEnergy();
268 }
269 }*/
270
271 /** \brief Stored projectile components
272 *
273 * These particles are owned by the ProjectileRemnant.
274 */
275 std::map<long, Particle*> storedComponents;
276
277 /// \brief Initial energy levels of the projectile
278 EnergyLevelMap theInitialEnergyLevels;
279
280 /// \brief Ground-state energies for any number of nucleons
281 EnergyLevels theGroundStateEnergies;
282
283 INCL_DECLARE_ALLOCATION_POOL(ProjectileRemnant)
284 };
285 }
286
287 #endif // G4INCLPROJECTILEREMNANT_HH_
288
289