| 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 #ifndef G4INCLCascade_hh
39 #define G4INCLCascade_hh 1
40
41 #include "G4INCLParticle.hh"
42 #include "G4INCLNucleus.hh"
43 #include "G4INCLIPropagationModel.hh"
44 #include "G4INCLCascadeAction.hh"
45 #include "G4INCLEventInfo.hh"
46 #include "G4INCLGlobalInfo.hh"
47 #include "G4INCLLogger.hh"
48 #include "G4INCLConfig.hh"
49 #include "G4INCLRootFinder.hh"
50 #ifndef INCLXX_IN_GEANT4_MODE
51 #include "G4INCLGeant4Compat.hh"
52 #endif
53
54
55 namespace G4INCL {
56 class INCL {
57 public:
58 INCL(Config const * const config);
59
60 ~INCL();
61
62 /// \brief Dummy copy constructor to silence Coverity warning
63 INCL(const INCL &rhs);
64
65 /// \brief Dummy assignment operator to silence Coverity warning
66 INCL &operator=(const INCL &rhs);
67
68 G4bool prepareReaction(const ParticleSpecies &projectileSpecies, const G4double kineticEnergy, const G4int A, const G4int Z, const G4int S);
69
70 G4bool initializeTarget(const G4int A, const G4int Z, const G4int S, AnnihilationType theAType);
71 G4double initUniverseRadiusForAntiprotonAtRest(const G4int A, const G4int Z);
72 inline const EventInfo &processEvent() {
73 return processEvent(
74 theConfig->getProjectileSpecies(),
75 theConfig->getProjectileKineticEnergy(),
76 theConfig->getTargetA(),
77 theConfig->getTargetZ(),
78 theConfig->getTargetS()
79 );
80 }
81 const EventInfo &processEvent(
82 ParticleSpecies const &projectileSpecies,
83 const G4double kineticEnergy,
84 const G4int targetA,
85 const G4int targetZ,
86 const G4int targetS
87 );
88
89 void finalizeGlobalInfo(Random::SeedVector const &initialSeeds);
90 const GlobalInfo &getGlobalInfo() const { return theGlobalInfo; }
91
92
93 private:
94 IPropagationModel *propagationModel;
95 G4int theA, theZ, theS;
96 G4bool targetInitSuccess;
97 G4double maxImpactParameter;
98 G4double maxUniverseRadius;
99 G4double maxInteractionDistance;
100 G4double fixedImpactParameter;
101 CascadeAction *cascadeAction;
102 Config const * const theConfig;
103 Nucleus *nucleus;
104 G4bool forceTransparent;
105
106 EventInfo theEventInfo;
107 GlobalInfo theGlobalInfo;
108
109 /// \brief Remnant size below which cascade stops
110 G4int minRemnantSize;
111
112 /// \brief Class to adjust remnant recoil
113 class RecoilFunctor : public RootFunctor {
114 public:
115 /** \brief Prepare for calling the () operator and scaleParticleEnergies
116 *
117 * The constructor sets the private class members.
118 */
119 RecoilFunctor(Nucleus * const n, const EventInfo &ei) :
120 RootFunctor(0., 1E6),
121 nucleus(n),
122 outgoingParticles(n->getStore()->getOutgoingParticles()),
123 theEventInfo(ei) {
124 for(ParticleIter p=outgoingParticles.begin(), e=outgoingParticles.end(); p!=e; ++p) {
125 particleMomenta.push_back((*p)->getMomentum());
126 particleKineticEnergies.push_back((*p)->getKineticEnergy());
127 }
128 ProjectileRemnant * const aPR = n->getProjectileRemnant();
129 if(aPR && aPR->getA()>0) {
130 particleMomenta.push_back(aPR->getMomentum());
131 particleKineticEnergies.push_back(aPR->getKineticEnergy());
132 outgoingParticles.push_back(aPR);
133 }
134 }
135 virtual ~RecoilFunctor() {}
136
137 /** \brief Compute the energy-conservation violation.
138 *
139 * \param x scale factor for the particle energies
140 * \return the energy-conservation violation
141 */
142 G4double operator()(const G4double x) const {
143 scaleParticleEnergies(x);
144 return nucleus->getConservationBalance(theEventInfo,true).energy;
145 }
146
147 /// \brief Clean up after root finding
148 void cleanUp(const G4bool success) const {
149 if(!success)
150 scaleParticleEnergies(1.);
151 }
152
153 private:
154 /// \brief Pointer to the nucleus
155 Nucleus *nucleus;
156 /// \brief List of final-state particles.
157 ParticleList outgoingParticles;
158 // \brief Reference to the EventInfo object
159 EventInfo const &theEventInfo;
160 /// \brief Initial momenta of the outgoing particles
161 std::list<ThreeVector> particleMomenta;
162 /// \brief Initial kinetic energies of the outgoing particles
163 std::list<G4double> particleKineticEnergies;
164
165 /** \brief Scale the kinetic energies of the outgoing particles.
166 *
167 * \param rescale scale factor
168 */
169 void scaleParticleEnergies(const G4double rescale) const {
170 // Rescale the energies (and the momenta) of the outgoing particles.
171 ThreeVector pBalance = nucleus->getIncomingMomentum();
172 std::list<ThreeVector>::const_iterator iP = particleMomenta.begin();
173 std::list<G4double>::const_iterator iE = particleKineticEnergies.begin();
174 for( ParticleIter i = outgoingParticles.begin(), e = outgoingParticles.end(); i!=e; ++i, ++iP, ++iE)
175 {
176 const G4double mass = (*i)->getMass();
177 const G4double newKineticEnergy = (*iE) * rescale;
178
179 (*i)->setMomentum(*iP);
180 (*i)->setEnergy(mass + newKineticEnergy);
181 (*i)->adjustMomentumFromEnergy();
182
183 pBalance -= (*i)->getMomentum();
184 }
185
186 nucleus->setMomentum(pBalance);
187 const G4double remnantMass = ParticleTable::getTableMass(nucleus->getA(),nucleus->getZ(),nucleus->getS()) + nucleus->getExcitationEnergy();
188 const G4double pRem2 = pBalance.mag2();
189 const G4double recoilEnergy = pRem2/
190 (std::sqrt(pRem2+remnantMass*remnantMass) + remnantMass);
191 nucleus->setEnergy(remnantMass + recoilEnergy);
192 }
193 };
194
195 /// \brief Class to adjust remnant recoil in the reaction CM system
196 class RecoilCMFunctor : public RootFunctor {
197 public:
198 /** \brief Prepare for calling the () operator and scaleParticleEnergies
199 *
200 * The constructor sets the private class members.
201 */
202 RecoilCMFunctor(Nucleus * const n, const EventInfo &ei) :
203 RootFunctor(0., 1E6),
204 nucleus(n),
205 theIncomingMomentum(nucleus->getIncomingMomentum()),
206 outgoingParticles(n->getStore()->getOutgoingParticles()),
207 theEventInfo(ei) {
208 if(theIncomingMomentum.mag() == 0.){ //change the condition
209 thePTBoostVector = {0., 0., 0.};
210 //INCL_WARN("PTBoostVector at rest is zero" << '\n');
211 }
212 else{
213 thePTBoostVector = nucleus->getIncomingMomentum()/(nucleus->getInitialEnergy()); //D
214 //INCL_WARN("PTBoostVector" << '\n');
215 }
216 for(ParticleIter p=outgoingParticles.begin(), e=outgoingParticles.end(); p!=e; ++p) {
217 (*p)->boost(thePTBoostVector);
218 particleCMMomenta.push_back((*p)->getMomentum());
219 }
220 ProjectileRemnant * const aPR = n->getProjectileRemnant();
221 if(aPR && aPR->getA()>0) {
222 aPR->boost(thePTBoostVector);
223 particleCMMomenta.push_back(aPR->getMomentum());
224 outgoingParticles.push_back(aPR);
225 }
226 }
227 virtual ~RecoilCMFunctor() {}
228
229 /** \brief Compute the energy-conservation violation.
230 *
231 * \param x scale factor for the particle energies
232 * \return the energy-conservation violation
233 */
234 G4double operator()(const G4double x) const {
235 scaleParticleCMMomenta(x);
236 return nucleus->getConservationBalance(theEventInfo,true).energy;
237 }
238
239 /// \brief Clean up after root finding
240 void cleanUp(const G4bool success) const {
241 if(!success)
242 scaleParticleCMMomenta(1.);
243 }
244
245 private:
246 /// \brief Pointer to the nucleus
247 Nucleus *nucleus;
248 /// \brief Projectile-target CM boost vector
249 ThreeVector thePTBoostVector;
250 /// \brief Incoming momentum
251 ThreeVector theIncomingMomentum;
252 /// \brief List of final-state particles.
253 ParticleList outgoingParticles;
254 /// \brief Reference to the EventInfo object
255 EventInfo const &theEventInfo;
256 /// \brief Initial CM momenta of the outgoing particles
257 std::list<ThreeVector> particleCMMomenta;
258
259 /** \brief Scale the kinetic energies of the outgoing particles.
260 *
261 * \param rescale scale factor
262 */
263 void scaleParticleCMMomenta(const G4double rescale) const {
264 // Rescale the CM momenta of the outgoing particles.
265 ThreeVector remnantMomentum = theIncomingMomentum;
266 std::list<ThreeVector>::const_iterator iP = particleCMMomenta.begin();
267 for( ParticleIter i = outgoingParticles.begin(), e = outgoingParticles.end(); i!=e; ++i, ++iP)
268 {
269 (*i)->setMomentum(*iP * rescale);
270 (*i)->adjustEnergyFromMomentum();
271 (*i)->boost(-thePTBoostVector);
272
273 remnantMomentum -= (*i)->getMomentum();
274 }
275
276 nucleus->setMomentum(remnantMomentum);
277 const G4double remnantMass = ParticleTable::getTableMass(nucleus->getA(),nucleus->getZ(),nucleus->getS()) + nucleus->getExcitationEnergy();
278 const G4double pRem2 = remnantMomentum.mag2();
279 const G4double recoilEnergy = pRem2/
280 (std::sqrt(pRem2+remnantMass*remnantMass) + remnantMass);
281 nucleus->setEnergy(remnantMass + recoilEnergy);
282 }
283 };
284
285 /** \brief Rescale the energies of the outgoing particles.
286 *
287 * Allow for the remnant recoil energy by rescaling the energy (and
288 * momenta) of the outgoing particles.
289 */
290 void rescaleOutgoingForRecoil();
291
292 #ifndef INCLXX_IN_GEANT4_MODE
293 /** \brief Run global conservation checks
294 *
295 * Check that energy and momentum are correctly conserved. If not, issue
296 * a warning.
297 *
298 * Also feeds the balance variables in theEventInfo.
299 *
300 * \param afterRecoil whether to take into account nuclear recoil
301 */
302 void globalConservationChecks(G4bool afterRecoil);
303 #endif
304
305 /** \brief Stopping criterion for the cascade
306 *
307 * Returns true if the cascade should continue, and false if any of the
308 * stopping criteria is satisfied.
309 */
310 G4bool continueCascade();
311
312 /** \brief Make a projectile pre-fragment out of geometrical spectators
313 *
314 * The projectile pre-fragment is assigned an excitation energy given
315 * by \f$E_\mathrm{sp}-E_\mathrm{i,A}\f$, where \f$E_\mathrm{sp}\f$ is the
316 * sum of the energies of the spectator particles, and \f$E_\mathrm{i,A}\f$
317 * is the sum of the smallest \f$A\f$ particle energies initially present
318 * in the projectile, \f$A\f$ being the mass of the projectile
319 * pre-fragment. This is equivalent to assuming that the excitation
320 * energy is given by the sum of the transitions of all excited
321 * projectile components to the "holes" left by the participants.
322 *
323 * This method can modify the outgoing list and adds a projectile
324 * pre-fragment.
325 *
326 * \return the number of dynamical spectators that were merged back in
327 * the projectile
328 */
329 G4int makeProjectileRemnant();
330
331 /** \brief Make a compound nucleus
332 *
333 * Selects the projectile components that can actually enter their
334 * potential and puts them into the target nucleus. If the CN excitation
335 * energy turns out to be negative, the event is considered a
336 * transparent. This method modifies theEventInfo and theGlobalInfo.
337 */
338 void makeCompoundNucleus();
339
340 /// \brief Initialise the cascade
341 G4bool preCascade(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy);
342
343 /// \brief The actual cascade loop
344 void cascade();
345
346 /// \brief Finalise the cascade and clean up
347 void postCascade(const ParticleSpecies &projectileSpecies, const G4double kineticEnergy);
348
349 /** \brief Initialise the maximum interaction distance.
350 *
351 * Used in forced CN events.
352 */
353 void initMaxInteractionDistance(ParticleSpecies const &p, const G4double kineticEnergy);
354
355 /** \brief Initialize the universe radius.
356 *
357 * Used for determining the energy-dependent size of the volume particles
358 * live in.
359 */
360 void initUniverseRadius(ParticleSpecies const &p, const G4double kineticEnergy, const G4int A, const G4int Z);
361
362 /// \brief Update global counters and other members of theGlobalInfo object
363 void updateGlobalInfo();
364
365 /// \brief Fill probabilities and particle types from datafile and return probability sum for normalization
366 G4double read_file(std::string filename, std::vector<G4double>& probabilities,
367 std::vector<std::vector<G4String>>& particle_types);
368
369 /// \brief This function will tell you the FS line number from the datafile based on your random value
370 G4int findStringNumber(G4double rdm, std::vector<G4double> yields);
371
372 /// \brief Initialise the "cascade" for pbar on H1
373 void preCascade_pbarH1(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy);
374
375 /// \brief Initialise the "cascade" for pbar on H2
376 void preCascade_pbarH2(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy);
377
378 /// \brief Finalise the "cascade" and clean up for pbar on H1
379 void postCascade_pbarH1(ParticleList const &outgoingParticles);
380
381 /// \brief Finalise the "cascade" and clean up for pbar on H2
382 void postCascade_pbarH2(ParticleList const &outgoingParticles, ParticleList const &H2Particles);
383 };
384 }
385
386 #endif
387