ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLCascade.hh
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLCascade.hh
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"
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 
51 namespace G4INCL {
52  class INCL {
53  public:
54  INCL(Config const * const config);
55 
56  ~INCL();
57 
59  INCL(const INCL &rhs);
60 
62  INCL &operator=(const INCL &rhs);
63 
64  G4bool prepareReaction(const ParticleSpecies &projectileSpecies, const G4double kineticEnergy, const G4int A, const G4int Z, const G4int S);
65  G4bool initializeTarget(const G4int A, const G4int Z, const G4int S);
66  inline const EventInfo &processEvent() {
67  return processEvent(
73  );
74  }
75  const EventInfo &processEvent(
76  ParticleSpecies const &projectileSpecies,
77  const G4double kineticEnergy,
78  const G4int targetA,
79  const G4int targetZ,
80  const G4int targetS
81  );
82 
83  void finalizeGlobalInfo(Random::SeedVector const &initialSeeds);
84  const GlobalInfo &getGlobalInfo() const { return theGlobalInfo; }
85 
86 
87  private:
96  Config const * const theConfig;
99 
102 
105 
107  class RecoilFunctor : public RootFunctor {
108  public:
113  RecoilFunctor(Nucleus * const n, const EventInfo &ei) :
114  RootFunctor(0., 1E6),
115  nucleus(n),
116  outgoingParticles(n->getStore()->getOutgoingParticles()),
117  theEventInfo(ei) {
118  for(ParticleIter p=outgoingParticles.begin(), e=outgoingParticles.end(); p!=e; ++p) {
119  particleMomenta.push_back((*p)->getMomentum());
120  particleKineticEnergies.push_back((*p)->getKineticEnergy());
121  }
122  ProjectileRemnant * const aPR = n->getProjectileRemnant();
123  if(aPR && aPR->getA()>0) {
124  particleMomenta.push_back(aPR->getMomentum());
125  particleKineticEnergies.push_back(aPR->getKineticEnergy());
126  outgoingParticles.push_back(aPR);
127  }
128  }
129  virtual ~RecoilFunctor() {}
130 
136  G4double operator()(const G4double x) const {
139  }
140 
142  void cleanUp(const G4bool success) const {
143  if(!success)
145  }
146 
147  private:
152  // \brief Reference to the EventInfo object
155  std::list<ThreeVector> particleMomenta;
157  std::list<G4double> particleKineticEnergies;
158 
163  void scaleParticleEnergies(const G4double rescale) const {
164  // Rescale the energies (and the momenta) of the outgoing particles.
166  std::list<ThreeVector>::const_iterator iP = particleMomenta.begin();
167  std::list<G4double>::const_iterator iE = particleKineticEnergies.begin();
168  for( ParticleIter i = outgoingParticles.begin(), e = outgoingParticles.end(); i!=e; ++i, ++iP, ++iE)
169  {
170  const G4double mass = (*i)->getMass();
171  const G4double newKineticEnergy = (*iE) * rescale;
172 
173  (*i)->setMomentum(*iP);
174  (*i)->setEnergy(mass + newKineticEnergy);
175  (*i)->adjustMomentumFromEnergy();
176 
177  pBalance -= (*i)->getMomentum();
178  }
179 
180  nucleus->setMomentum(pBalance);
182  const G4double pRem2 = pBalance.mag2();
183  const G4double recoilEnergy = pRem2/
184  (std::sqrt(pRem2+remnantMass*remnantMass) + remnantMass);
185  nucleus->setEnergy(remnantMass + recoilEnergy);
186  }
187  };
188 
190  class RecoilCMFunctor : public RootFunctor {
191  public:
196  RecoilCMFunctor(Nucleus * const n, const EventInfo &ei) :
197  RootFunctor(0., 1E6),
198  nucleus(n),
199  theIncomingMomentum(nucleus->getIncomingMomentum()),
200  outgoingParticles(n->getStore()->getOutgoingParticles()),
201  theEventInfo(ei) {
203  for(ParticleIter p=outgoingParticles.begin(), e=outgoingParticles.end(); p!=e; ++p) {
204  (*p)->boost(thePTBoostVector);
205  particleCMMomenta.push_back((*p)->getMomentum());
206  }
207  ProjectileRemnant * const aPR = n->getProjectileRemnant();
208  if(aPR && aPR->getA()>0) {
209  aPR->boost(thePTBoostVector);
210  particleCMMomenta.push_back(aPR->getMomentum());
211  outgoingParticles.push_back(aPR);
212  }
213  }
214  virtual ~RecoilCMFunctor() {}
215 
221  G4double operator()(const G4double x) const {
224  }
225 
227  void cleanUp(const G4bool success) const {
228  if(!success)
230  }
231 
232  private:
244  std::list<ThreeVector> particleCMMomenta;
245 
250  void scaleParticleCMMomenta(const G4double rescale) const {
251  // Rescale the CM momenta of the outgoing particles.
252  ThreeVector remnantMomentum = theIncomingMomentum;
253  std::list<ThreeVector>::const_iterator iP = particleCMMomenta.begin();
254  for( ParticleIter i = outgoingParticles.begin(), e = outgoingParticles.end(); i!=e; ++i, ++iP)
255  {
256  (*i)->setMomentum(*iP * rescale);
257  (*i)->adjustEnergyFromMomentum();
258  (*i)->boost(-thePTBoostVector);
259 
260  remnantMomentum -= (*i)->getMomentum();
261  }
262 
263  nucleus->setMomentum(remnantMomentum);
265  const G4double pRem2 = remnantMomentum.mag2();
266  const G4double recoilEnergy = pRem2/
267  (std::sqrt(pRem2+remnantMass*remnantMass) + remnantMass);
268  nucleus->setEnergy(remnantMass + recoilEnergy);
269  }
270  };
271 
278 
279 #ifndef INCLXX_IN_GEANT4_MODE
280 
289  void globalConservationChecks(G4bool afterRecoil);
290 #endif
291 
298 
317 
325  void makeCompoundNucleus();
326 
328  G4bool preCascade(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy);
329 
331  void cascade();
332 
334  void postCascade();
335 
340  void initMaxInteractionDistance(ParticleSpecies const &p, const G4double kineticEnergy);
341 
347  void initUniverseRadius(ParticleSpecies const &p, const G4double kineticEnergy, const G4int A, const G4int Z);
348 
350  void updateGlobalInfo();
351  };
352 }
353 
354 #endif