Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLCascade.hh
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 
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);
65  G4bool initializeTarget(const G4int A, const G4int Z);
66  inline const EventInfo &processEvent() {
67  return processEvent(
68  theConfig->getProjectileSpecies(),
69  theConfig->getProjectileKineticEnergy(),
70  theConfig->getTargetA(),
71  theConfig->getTargetZ()
72  );
73  }
74  const EventInfo &processEvent(
75  ParticleSpecies const &projectileSpecies,
76  const G4double kineticEnergy,
77  const G4int targetA,
78  const G4int targetZ
79  );
80 
81  void finalizeGlobalInfo(Random::SeedVector const &initialSeeds);
82  const GlobalInfo &getGlobalInfo() const { return theGlobalInfo; }
83 
84  private:
85  IPropagationModel *propagationModel;
86  G4int theA, theZ;
87  G4bool targetInitSuccess;
88  G4double maxImpactParameter;
89  G4double maxUniverseRadius;
90  G4double maxInteractionDistance;
91  G4double fixedImpactParameter;
92  CascadeAction *cascadeAction;
93  Config const * const theConfig;
94  Nucleus *nucleus;
95  G4bool forceTransparent;
96 
97  EventInfo theEventInfo;
98  GlobalInfo theGlobalInfo;
99 
101  G4int minRemnantSize;
102 
104  class RecoilFunctor : public RootFunctor {
105  public:
110  RecoilFunctor(Nucleus * const n, const EventInfo &ei) :
111  RootFunctor(0., 1E6),
112  nucleus(n),
113  outgoingParticles(n->getStore()->getOutgoingParticles()),
114  theEventInfo(ei) {
115  for(ParticleIter p=outgoingParticles.begin(), e=outgoingParticles.end(); p!=e; ++p) {
116  particleMomenta.push_back((*p)->getMomentum());
117  particleKineticEnergies.push_back((*p)->getKineticEnergy());
118  }
119  ProjectileRemnant * const aPR = n->getProjectileRemnant();
120  if(aPR && aPR->getA()>0) {
121  particleMomenta.push_back(aPR->getMomentum());
122  particleKineticEnergies.push_back(aPR->getKineticEnergy());
123  outgoingParticles.push_back(aPR);
124  }
125  }
126  virtual ~RecoilFunctor() {}
127 
133  G4double operator()(const G4double x) const {
134  scaleParticleEnergies(x);
135  return nucleus->getConservationBalance(theEventInfo,true).energy;
136  }
137 
139  void cleanUp(const G4bool success) const {
140  if(!success)
141  scaleParticleEnergies(1.);
142  }
143 
144  private:
146  Nucleus *nucleus;
148  ParticleList outgoingParticles;
149  // \brief Reference to the EventInfo object
150  EventInfo const &theEventInfo;
152  std::list<ThreeVector> particleMomenta;
154  std::list<G4double> particleKineticEnergies;
155 
160  void scaleParticleEnergies(const G4double rescale) const {
161  // Rescale the energies (and the momenta) of the outgoing particles.
162  ThreeVector pBalance = nucleus->getIncomingMomentum();
163  std::list<ThreeVector>::const_iterator iP = particleMomenta.begin();
164  std::list<G4double>::const_iterator iE = particleKineticEnergies.begin();
165  for( ParticleIter i = outgoingParticles.begin(), e = outgoingParticles.end(); i!=e; ++i, ++iP, ++iE)
166  {
167  const G4double mass = (*i)->getMass();
168  const G4double newKineticEnergy = (*iE) * rescale;
169 
170  (*i)->setMomentum(*iP);
171  (*i)->setEnergy(mass + newKineticEnergy);
172  (*i)->adjustMomentumFromEnergy();
173 
174  pBalance -= (*i)->getMomentum();
175  }
176 
177  nucleus->setMomentum(pBalance);
178  const G4double remnantMass = ParticleTable::getTableMass(nucleus->getA(),nucleus->getZ()) + nucleus->getExcitationEnergy();
179  const G4double pRem2 = pBalance.mag2();
180  const G4double recoilEnergy = pRem2/
181  (std::sqrt(pRem2+remnantMass*remnantMass) + remnantMass);
182  nucleus->setEnergy(remnantMass + recoilEnergy);
183  }
184  };
185 
187  class RecoilCMFunctor : public RootFunctor {
188  public:
193  RecoilCMFunctor(Nucleus * const n, const EventInfo &ei) :
194  RootFunctor(0., 1E6),
195  nucleus(n),
196  theIncomingMomentum(nucleus->getIncomingMomentum()),
197  outgoingParticles(n->getStore()->getOutgoingParticles()),
198  theEventInfo(ei) {
199  thePTBoostVector = nucleus->getIncomingMomentum()/nucleus->getInitialEnergy();
200  for(ParticleIter p=outgoingParticles.begin(), e=outgoingParticles.end(); p!=e; ++p) {
201  (*p)->boost(thePTBoostVector);
202  particleCMMomenta.push_back((*p)->getMomentum());
203  }
204  ProjectileRemnant * const aPR = n->getProjectileRemnant();
205  if(aPR && aPR->getA()>0) {
206  aPR->boost(thePTBoostVector);
207  particleCMMomenta.push_back(aPR->getMomentum());
208  outgoingParticles.push_back(aPR);
209  }
210  }
211  virtual ~RecoilCMFunctor() {}
212 
218  G4double operator()(const G4double x) const {
219  scaleParticleCMMomenta(x);
220  return nucleus->getConservationBalance(theEventInfo,true).energy;
221  }
222 
224  void cleanUp(const G4bool success) const {
225  if(!success)
226  scaleParticleCMMomenta(1.);
227  }
228 
229  private:
231  Nucleus *nucleus;
233  ThreeVector thePTBoostVector;
235  ThreeVector theIncomingMomentum;
237  ParticleList outgoingParticles;
238  // \brief Reference to the EventInfo object
239  EventInfo const &theEventInfo;
241  std::list<ThreeVector> particleCMMomenta;
242 
247  void scaleParticleCMMomenta(const G4double rescale) const {
248  // Rescale the CM momenta of the outgoing particles.
249  ThreeVector remnantMomentum = theIncomingMomentum;
250  std::list<ThreeVector>::const_iterator iP = particleCMMomenta.begin();
251  for( ParticleIter i = outgoingParticles.begin(), e = outgoingParticles.end(); i!=e; ++i, ++iP)
252  {
253  (*i)->setMomentum(*iP * rescale);
254  (*i)->adjustEnergyFromMomentum();
255  (*i)->boost(-thePTBoostVector);
256 
257  remnantMomentum -= (*i)->getMomentum();
258  }
259 
260  nucleus->setMomentum(remnantMomentum);
261  const G4double remnantMass = ParticleTable::getTableMass(nucleus->getA(),nucleus->getZ()) + nucleus->getExcitationEnergy();
262  const G4double pRem2 = remnantMomentum.mag2();
263  const G4double recoilEnergy = pRem2/
264  (std::sqrt(pRem2+remnantMass*remnantMass) + remnantMass);
265  nucleus->setEnergy(remnantMass + recoilEnergy);
266  }
267  };
268 
274  void rescaleOutgoingForRecoil();
275 
276 #ifndef INCLXX_IN_GEANT4_MODE
277 
286  void globalConservationChecks(G4bool afterRecoil);
287 #endif
288 
294  G4bool continueCascade();
295 
313  G4int makeProjectileRemnant();
314 
322  void makeCompoundNucleus();
323 
325  G4bool preCascade(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy);
326 
328  void cascade();
329 
331  void postCascade();
332 
337  void initMaxInteractionDistance(ParticleSpecies const &p, const G4double kineticEnergy);
338 
344  void initUniverseRadius(ParticleSpecies const &p, const G4double kineticEnergy, const G4int A, const G4int Z);
345 
347  void updateGlobalInfo();
348  };
349 }
350 
351 #endif
INCL & operator=(const INCL &rhs)
Dummy assignment operator to silence Coverity warning.
const GlobalInfo & getGlobalInfo() const
G4bool prepareReaction(const ParticleSpecies &projectileSpecies, const G4double kineticEnergy, const G4int A, const G4int Z)
G4int getTargetZ() const
Get the target charge number.
Definition: G4INCLConfig.hh:97
const char * p
Definition: xmltok.h:285
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
Class containing default actions to be performed at intermediate cascade steps.
void finalizeGlobalInfo(Random::SeedVector const &initialSeeds)
INCL(Config const *const config)
Simple container for output of calculation-wide results.
double A(double temperature)
Simple container for output of event results.
bool G4bool
Definition: G4Types.hh:79
G4int getTargetA() const
Get the target mass number.
Definition: G4INCLConfig.hh:94
const EventInfo & processEvent()
const G4int n
G4bool initializeTarget(const G4int A, const G4int Z)
ParticleSpecies getProjectileSpecies() const
Get the projectile species.
G4ThreadLocal NuclearMassFn getTableMass
Static pointer to the mass function for nuclei.
RootFunctor(const G4double x0, const G4double x1)
ProjectileRemnant * getProjectileRemnant() const
Get the projectile remnant.
double G4double
Definition: G4Types.hh:76
G4double getProjectileKineticEnergy() const
Get the projectile kinetic energy.
ParticleList::const_iterator ParticleIter
Static root-finder algorithm.