Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLNucleus.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 /*
39  * G4INCLNucleus.hh
40  *
41  * \date Jun 5, 2009
42  * \author Pekka Kaitaniemi
43  */
44 
45 #ifndef G4INCLNUCLEUS_HH_
46 #define G4INCLNUCLEUS_HH_
47 
48 #include <list>
49 #include <string>
50 
51 #include "G4INCLParticle.hh"
52 #include "G4INCLEventInfo.hh"
53 #include "G4INCLCluster.hh"
54 #include "G4INCLFinalState.hh"
55 #include "G4INCLStore.hh"
56 #include "G4INCLGlobals.hh"
57 #include "G4INCLParticleTable.hh"
58 #include "G4INCLConfig.hh"
59 #include "G4INCLConfigEnums.hh"
60 #include "G4INCLCluster.hh"
62 
63 namespace G4INCL {
64 
65  class Nucleus : public Cluster {
66  public:
67  Nucleus(G4int mass, G4int charge, Config const * const conf, const G4double universeRadius=-1.);
68  virtual ~Nucleus();
69 
71  Nucleus(const Nucleus &rhs);
72 
74  Nucleus &operator=(const Nucleus &rhs);
75 
80  void initializeParticles();
81 
84  theZ += p->getZ();
85  theA += p->getA();
86  theStore->particleHasEntered(p);
87  if(p->isNucleon()) {
88  theNpInitial += Math::heaviside(ParticleTable::getIsospin(p->getType()));
89  theNnInitial += Math::heaviside(-ParticleTable::getIsospin(p->getType()));
90  }
91  if(!p->isTargetSpectator()) theStore->getBook().incrementCascading();
92  };
93 
98 
99  G4int getInitialA() const { return theInitialA; };
100  G4int getInitialZ() const { return theInitialZ; };
101 
107  void propagateParticles(G4double step);
108 
109  G4int getNumberOfEnteringProtons() const { return theNpInitial; };
110  G4int getNumberOfEnteringNeutrons() const { return theNnInitial; };
111 
117  G4double S = 0.0;
118  ParticleList const &outgoing = theStore->getOutgoingParticles();
119  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
120  const ParticleType t = (*i)->getType();
121  switch(t) {
122  case Proton:
123  case Neutron:
124  case DeltaPlusPlus:
125  case DeltaPlus:
126  case DeltaZero:
127  case DeltaMinus:
128  S += thePotential->getSeparationEnergy(*i);
129  break;
130  case Composite:
131  S += (*i)->getZ() * thePotential->getSeparationEnergy(Proton)
132  + ((*i)->getA() - (*i)->getZ()) * thePotential->getSeparationEnergy(Neutron);
133  break;
134  case PiPlus:
135  S += thePotential->getSeparationEnergy(Proton) - thePotential->getSeparationEnergy(Neutron);
136  break;
137  case PiMinus:
138  S += thePotential->getSeparationEnergy(Neutron) - thePotential->getSeparationEnergy(Proton);
139  break;
140  default:
141  break;
142  }
143  }
144 
145  S -= theNpInitial * thePotential->getSeparationEnergy(Proton);
146  S -= theNnInitial * thePotential->getSeparationEnergy(Neutron);
147  return S;
148  }
149 
155 
161 
167 
173 
180  G4bool decayMe();
181 
183  void emitInsidePions();
184 
187 
193 
199 
205 
208  incomingAngularMomentum = j;
209  }
210 
212  const ThreeVector &getIncomingAngularMomentum() const { return incomingAngularMomentum; }
213 
216  incomingMomentum = p;
217  }
218 
221  return incomingMomentum;
222  }
223 
225  void setInitialEnergy(const G4double e) { initialEnergy = e; }
226 
228  G4double getInitialEnergy() const { return initialEnergy; }
229 
235 
238  ParticleList const &inside = theStore->getParticles();
239  for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i)
240  if((*i)->isDelta()) return true;
241  return false;
242  }
243 
245  inline G4bool containsEtas() {
246  ParticleList const &inside = theStore->getParticles();
247  for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i)
248  if((*i)->isEta()) return true;
249  return false;
250  }
251 
254  ParticleList const &inside = theStore->getParticles();
255  for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i)
256  if((*i)->isOmega()) return true;
257  return false;
258  }
259 
263  std::string print();
264 
265  Store* getStore() const {return theStore; };
266  void setStore(Store *s) {
267  delete theStore;
268  theStore = s;
269  };
270 
271  G4double getInitialInternalEnergy() const { return initialInternalEnergy; };
272 
277  G4bool isEventTransparent() const;
278 
283  G4bool hasRemnant() const { return remnant; }
284 
288  void fillEventInfo(EventInfo *eventInfo);
289 
290  G4bool getTryCompoundNucleus() { return tryCN; }
291 
294  const G4double theTransmissionRadius = theDensity->getTransmissionRadius(p);
295  const G4double theParticleZ = p->getZ();
296  return PhysicalConstants::eSquared*(theZ-theParticleZ)*theParticleZ/theTransmissionRadius;
297  }
298 
304  };
305 
307  ConservationBalance getConservationBalance(EventInfo const &theEventInfo, const G4bool afterRecoil) const;
308 
310  void useFusionKinematics();
311 
320  G4double getSurfaceRadius(Particle const * const particle) const {
321 // if(particle->isPion())
322  if(particle->isPion() || particle->isEta() || particle->isOmega() || particle->isEtaPrime())
323  // Temporarily set RPION = RMAX
324  return getUniverseRadius();
325  //return 0.5*(theDensity->getTransmissionRadius(particle)+getUniverseRadius());
326  else {
327  const G4double pr = particle->getReflectionMomentum()/thePotential->getFermiMomentum(particle);
328  if(pr>=1.)
329  return getUniverseRadius();
330  else
331  return theDensity->getMaxRFromP(particle->getType(), pr);
332  }
333  }
334 
336  G4double getUniverseRadius() const { return theUniverseRadius; }
337 
339  void setUniverseRadius(const G4double universeRadius) { theUniverseRadius=universeRadius; }
340 
342  G4bool isNucleusNucleusCollision() const { return isNucleusNucleus; }
343 
345  void setNucleusNucleusCollision() { isNucleusNucleus=true; }
346 
348  void setParticleNucleusCollision() { isNucleusNucleus=false; }
349 
352  delete theProjectileRemnant;
353  theProjectileRemnant = c;
354  }
355 
357  ProjectileRemnant *getProjectileRemnant() const { return theProjectileRemnant; }
358 
361  delete theProjectileRemnant;
362  theProjectileRemnant = NULL;
363  }
364 
373  void finalizeProjectileRemnant(const G4double emissionTime);
374 
376  inline void updatePotentialEnergy(Particle *p) const {
377  p->setPotentialEnergy(thePotential->computePotentialEnergy(p));
378  }
379 
381  void setDensity(NuclearDensity const * const d) {
382  theDensity=d;
384  theParticleSampler->setDensity(theDensity);
385  };
386 
388  NuclearDensity const *getDensity() const { return theDensity; };
389 
391  NuclearPotential::INuclearPotential const *getPotential() const { return thePotential; };
392 
393  private:
399  void computeOneNucleonRecoilKinematics();
400 
401  private:
402  G4int theInitialZ, theInitialA;
404  G4int theNpInitial;
406  G4int theNnInitial;
407  G4double initialInternalEnergy;
408  ThreeVector incomingAngularMomentum, incomingMomentum;
409  ThreeVector initialCenterOfMass;
410  G4bool remnant;
411 
412  G4double initialEnergy;
413  Store *theStore;
414  G4bool tryCN;
415 
417  G4int projectileZ;
419  G4int projectileA;
420 
422  G4double theUniverseRadius;
423 
428  G4bool isNucleusNucleus;
429 
434  ProjectileRemnant *theProjectileRemnant;
435 
437  NuclearDensity const *theDensity;
438 
440  NuclearPotential::INuclearPotential const *thePotential;
441 
443  };
444 
445 }
446 
447 #endif /* G4INCLNUCLEUS_HH_ */
G4bool isEta() const
Is this a eta?
G4int getNumberOfEnteringNeutrons() const
G4int getA() const
Returns the baryon number.
#define INCL_DECLARE_ALLOCATION_POOL(T)
ConservationBalance getConservationBalance(EventInfo const &theEventInfo, const G4bool afterRecoil) const
Compute charge, mass, energy and momentum balance.
void initializeParticles()
G4bool containsDeltas()
Returns true if the nucleus contains any deltas.
void updatePotentialEnergy(Particle *p) const
Update the particle potential energy.
G4double getReflectionMomentum() const
Return the reflection momentum.
void setIncomingAngularMomentum(const ThreeVector &j)
Set the incoming angular-momentum vector.
void setIncomingMomentum(const ThreeVector &p)
Set the incoming momentum vector.
const G4double eSquared
Coulomb conversion factor [MeV*fm].
G4int getInitialZ() const
ParticleList const & getParticles() const
Definition: G4INCLStore.hh:253
std::string print()
double S(double temp)
const ThreeVector & getIncomingMomentum() const
Get the incoming momentum vector.
G4double getFermiMomentum(const Particle *const p) const
Return the Fermi momentum for a particle.
virtual ~Nucleus()
G4bool decayOutgoingClusters()
Force the decay of unstable outgoing clusters.
void setInitialEnergy(const G4double e)
Set the initial energy.
G4int heaviside(G4int n)
G4bool isNucleusNucleusCollision() const
Is it a nucleus-nucleus collision?
const char * p
Definition: xmltok.h:285
G4int getInitialA() const
void applyFinalState(FinalState *)
void setDensity(NuclearDensity const *const d)
Setter for theDensity.
G4bool isTargetSpectator() const
const ThreeVector & getIncomingAngularMomentum() const
Get the incoming angular-momentum vector.
Store * getStore() const
Class for constructing a projectile-like remnant.
void setParticleNucleusCollision()
Set a particle-nucleus collision.
G4bool isEtaPrime() const
Is this a etaprime?
Struct for conservation laws.
G4double getSurfaceRadius(Particle const *const particle) const
Get the maximum allowed radius for a given particle.
G4bool hasRemnant() const
Does the nucleus give a cascade remnant?
Nucleus & operator=(const Nucleus &rhs)
Dummy assignment operator to silence Coverity warning.
int G4int
Definition: G4Types.hh:78
G4double getInitialInternalEnergy() const
const XML_Char * s
Definition: expat.h:262
void useFusionKinematics()
Adjust the kinematics for complete-fusion events.
G4double getInitialEnergy() const
Get the initial energy.
G4bool isOmega() const
Is this a omega?
G4double getSeparationEnergy(const Particle *const p) const
Return the separation energy for a particle.
G4bool decayOutgoingDeltas()
Force the decay of outgoing deltas.
ParticleList const & getOutgoingParticles() const
Definition: G4INCLStore.hh:223
Simple container for output of event results.
void setDensity(NuclearDensity const *const d)
Setter for theDensity.
Book & getBook()
Definition: G4INCLStore.hh:259
bool G4bool
Definition: G4Types.hh:79
G4int getZ() const
Returns the charge number.
void setPotentialEnergy(G4double v)
Set the particle potential energy.
G4double computeSeparationEnergyBalance() const
Outgoing - incoming separation energies.
void setUniverseRadius(const G4double universeRadius)
Setter for theUniverseRadius.
G4double getExcitationEnergy() const
Get the excitation energy of the nucleus.
G4double computeExcitationEnergy() const
Compute the current excitation energy.
G4int getNumberOfEnteringProtons() const
void deleteProjectileRemnant()
Delete the projectile remnant.
G4bool decayMe()
Force the phase-space decay of the Nucleus.
ParticleSampler * theParticleSampler
ThreeVector computeCenterOfMass() const
Compute the current center-of-mass position.
void incrementCascading()
Definition: G4INCLBook.hh:77
G4bool getTryCompoundNucleus()
NuclearPotential::INuclearPotential const * getPotential() const
Getter for thePotential.
G4double theExcitationEnergy
void propagateParticles(G4double step)
G4bool decayInsideDeltas()
Force the decay of deltas inside the nucleus.
void finalizeProjectileRemnant(const G4double emissionTime)
Finalise the projectile remnant.
G4double getTransmissionBarrier(Particle const *const p)
Get the transmission barrier.
G4int getIsospin(const ParticleType t)
Get the isospin of a particle.
void emitInsidePions()
Force emission of all pions inside the nucleus.
G4INCL::ParticleType getType() const
NuclearDensity const * getDensity() const
Getter for theDensity.
void setProjectileRemnant(ProjectileRemnant *const c)
Set the projectile remnant.
void setStore(Store *s)
G4bool isEventTransparent() const
Is the event transparent?
void fillEventInfo(EventInfo *eventInfo)
G4bool isNucleon() const
G4double computeTotalEnergy() const
Compute the current total energy.
ProjectileRemnant * getProjectileRemnant() const
Get the projectile remnant.
G4double getUniverseRadius() const
Getter for theUniverseRadius.
void insertParticle(Particle *p)
Insert a new particle (e.g. a projectile) in the nucleus.
virtual G4double computePotentialEnergy(const Particle *const p) const =0
double G4double
Definition: G4Types.hh:76
Nucleus(G4int mass, G4int charge, Config const *const conf, const G4double universeRadius=-1.)
tuple c
Definition: test.py:13
G4double getTransmissionRadius(Particle const *const p) const
The radius used for calculating the transmission coefficient.
void computeRecoilKinematics()
Compute the recoil momentum and spin of the nucleus.
G4bool containsEtas()
Returns true if the nucleus contains any etas.
G4bool isPion() const
Is this a pion?
void setNucleusNucleusCollision()
Set a nucleus-nucleus collision.
G4bool decayOutgoingPionResonances(G4double timeThreshold)
Force the decay of outgoing PionResonances (eta/omega).
G4double getMaxRFromP(const ParticleType t, const G4double p) const
Get the maximum allowed radius for a given momentum.
ParticleList::const_iterator ParticleIter
void particleHasEntered(Particle *const particle)
Move a particle from incoming to inside.
Definition: G4INCLStore.cc:188
G4bool containsOmegas()
Returns true if the nucleus contains any omegas.