Geant4  10.02.p01
G4INCLParticleEntryChannel.cc
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 
39 #include "G4INCLRootFinder.hh"
40 #include "G4INCLIntersection.hh"
41 #include <algorithm>
42 
43 namespace G4INCL {
44 
46  :theNucleus(n), theParticle(p)
47  {}
48 
50  {}
51 
53  // Behaves slightly differency if a third body (the projectile) is present
55 
56  /* Corrections to the energy of the entering nucleon
57  *
58  * In particle-nucleus reactions, the goal of this correction is to satisfy
59  * energy conservation in particle-nucleus reactions using real particle
60  * and nuclear masses.
61  *
62  * In nucleus-nucleus reactions, in addition to the above, the correction
63  * is determined by a model for the excitation energy of the
64  * quasi-projectile (QP). The energy of the entering nucleon is such that
65  * the QP excitation energy, as determined by conservation, is what given
66  * by our model.
67  *
68  * Possible choices for the correction (or, equivalently, for the QP
69  * excitation energy):
70  *
71  * 1. the correction is 0. (same as in particle-nucleus);
72  * 2. the correction is the separation energy of the entering nucleon in
73  * the current QP;
74  * 3. the QP excitation energy is given by A. Boudard's algorithm, as
75  * implemented in INCL4.2-HI/Geant4.
76  * 4. the QP excitation energy vanishes.
77  *
78  * Ideally, the QP excitation energy should always be >=0. Algorithms 1.
79  * and 2. do not guarantee this, although violations to the rule seem to be
80  * more severe for 1. than for 2.. Algorithms 3. and 4., by construction,
81  * yields non-negative QP excitation energies.
82  */
83  G4double theCorrection;
84  if(isNN) {
85 // assert(theParticle->isNucleon());
86  ProjectileRemnant * const projectileRemnant = theNucleus->getProjectileRemnant();
87 // assert(projectileRemnant);
88 
89  // No correction (model 1. above)
90  /*
91  theCorrection = theParticle->getEmissionQValueCorrection(
92  theNucleus->getA() + theParticle->getA(),
93  theNucleus->getZ() + theParticle->getZ())
94  + theParticle->getTableMass() - theParticle->getINCLMass();
95  const G4double theProjectileCorrection = 0.;
96  */
97 
98  // Correct the energy of the entering particle for the Q-value of the
99  // emission from the projectile (model 2. above)
100  /*
101  theCorrection = theParticle->getTransferQValueCorrection(
102  projectileRemnant->getA(), projectileRemnant->getZ(),
103  theNucleus->getA(), theNucleus->getZ());
104  G4double theProjectileCorrection;
105  if(projectileRemnant->getA()>theParticle->getA()) { // if there are any particles left
106  // Compute the projectile Q-value (to be used as a correction to the
107  // other components of the projectile remnant)
108  theProjectileCorrection = ParticleTable::getTableQValue(
109  projectileRemnant->getA() - theParticle->getA(),
110  projectileRemnant->getZ() - theParticle->getZ(),
111  theParticle->getA(),
112  theParticle->getZ());
113  } else
114  theProjectileCorrection = 0.;
115  */
116 
117  // Fix the correction in such a way that the quasi-projectile excitation
118  // energy is given by A. Boudard's INCL4.2-HI model (model 3. above).
119  const G4double theProjectileExcitationEnergy =
120  (projectileRemnant->getA()-theParticle->getA()>1) ?
121  (projectileRemnant->computeExcitationEnergyExcept(theParticle->getID())) :
122  0.;
123  // Set the projectile excitation energy to zero (cold quasi-projectile,
124  // model 4. above).
125  // const G4double theProjectileExcitationEnergy = 0.;
126  // The part that follows is common to model 3. and 4.
127  const G4double theProjectileEffectiveMass =
128  ParticleTable::getTableMass(projectileRemnant->getA() - theParticle->getA(), projectileRemnant->getZ() - theParticle->getZ())
129  + theProjectileExcitationEnergy;
130  const ThreeVector &theProjectileMomentum = projectileRemnant->getMomentum() - theParticle->getMomentum();
131  const G4double theProjectileEnergy = std::sqrt(theProjectileMomentum.mag2() + theProjectileEffectiveMass*theProjectileEffectiveMass);
132  const G4double theProjectileCorrection = theProjectileEnergy - (projectileRemnant->getEnergy() - theParticle->getEnergy());
133  theCorrection = theParticle->getEmissionQValueCorrection(
137  + theProjectileCorrection;
138  // end of part common to model 3. and 4.
139 
140 
141  projectileRemnant->removeParticle(theParticle, theProjectileCorrection);
142  } else {
143  const G4int ACN = theNucleus->getA() + theParticle->getA();
144  const G4int ZCN = theNucleus->getZ() + theParticle->getZ();
145  // Correction to the Q-value of the entering particle
146  theCorrection = theParticle->getEmissionQValueCorrection(ACN,ZCN);
147  INCL_DEBUG("The following Particle enters with correction " << theCorrection << '\n'
148  << theParticle->print() << '\n');
149  }
150 
151  const G4double energyBefore = theParticle->getEnergy() - theCorrection;
152  G4bool success = particleEnters(theCorrection);
154 
155  if(!success) {
156  fs->makeParticleBelowZero();
157  } else if(theParticle->isNucleon() &&
159  // If the participant is a nucleon entering below its Fermi energy, force a
160  // compound nucleus
162  }
163 
164  fs->setTotalEnergyBeforeInteraction(energyBefore);
165  }
166 
168 
169  // \todo{this is the place to add refraction}
170 
171  theParticle->setINCLMass(); // Will automatically put the particle on shell
172 
173  // Add the nuclear potential to the kinetic energy when entering the
174  // nucleus
175 
176  class IncomingEFunctor : public RootFunctor {
177  public:
178  IncomingEFunctor(Particle * const p, Nucleus const * const n, const G4double correction) :
179  RootFunctor(0., 1E6),
180  theParticle(p),
181  thePotential(n->getPotential()),
182  theEnergy(theParticle->getEnergy()),
183  theMass(theParticle->getMass()),
184  theQValueCorrection(correction),
185  refraction(n->getStore()->getConfig()->getRefraction()),
186  theMomentumDirection(theParticle->getMomentum())
187  {
188  if(refraction) {
190  const G4double r2 = position.mag2();
191  if(r2>0.)
192  normal = - position / std::sqrt(r2);
193  G4double cosIncidenceAngle = theParticle->getCosRPAngle();
194  if(cosIncidenceAngle < -1.)
195  sinIncidenceAnglePOut = 0.;
196  else
197  sinIncidenceAnglePOut = theMomentumDirection.mag()*std::sqrt(1.-cosIncidenceAngle*cosIncidenceAngle);
198  } else {
199  sinIncidenceAnglePOut = 0.;
200  }
201  }
202  ~IncomingEFunctor() {}
203  G4double operator()(const G4double v) const {
204  G4double energyInside = std::max(theMass, theEnergy + v - theQValueCorrection);
205  theParticle->setEnergy(energyInside);
207  if(refraction) {
208  // Compute the new direction of the particle momentum
209  const G4double pIn = std::sqrt(energyInside*energyInside-theMass*theMass);
210  const G4double sinRefractionAngle = sinIncidenceAnglePOut/pIn;
211  const G4double cosRefractionAngle = (sinRefractionAngle>1.) ? 0. : std::sqrt(1.-sinRefractionAngle*sinRefractionAngle);
212  const ThreeVector momentumInside = theMomentumDirection - normal * normal.dot(theMomentumDirection) + normal * (pIn * cosRefractionAngle);
213  theParticle->setMomentum(momentumInside);
214  } else {
215  theParticle->setMomentum(theMomentumDirection); // keep the same direction
216  }
217  // Scale the particle momentum
219  return v - thePotential->computePotentialEnergy(theParticle);
220  }
221  void cleanUp(const G4bool /*success*/) const {}
222  private:
224  NuclearPotential::INuclearPotential const *thePotential;
225  const G4double theEnergy;
226  const G4double theMass;
227  const G4double theQValueCorrection;
228  const G4bool refraction;
229  const ThreeVector theMomentumDirection;
231  G4double sinIncidenceAnglePOut;
232  } theIncomingEFunctor(theParticle,theNucleus,theQValueCorrection);
233 
235  if(theParticle->getKineticEnergy()+v-theQValueCorrection<0.) { // Particle entering below 0. Die gracefully
236  INCL_DEBUG("Particle " << theParticle->getID() << " is trying to enter below 0" << '\n');
237  return false;
238  }
239  const RootFinder::Solution theSolution = RootFinder::solve(&theIncomingEFunctor, v);
240  if(theSolution.success) { // Apply the solution
241  theIncomingEFunctor(theSolution.x);
242  INCL_DEBUG("Particle successfully entered:\n" << theParticle->print() << '\n');
243  } else {
244  INCL_WARN("Couldn't compute the potential for incoming particle, root-finding algorithm failed." << '\n');
245  }
246  return theSolution.success;
247  }
248 
249 }
250 
G4int getA() const
Returns the baryon number.
G4double getCosRPAngle() const
Get the cosine of the angle between position and momentum.
G4double getMass() const
Get the cached particle mass.
G4bool particleEnters(const G4double theQValueCorrection)
Modify particle that enters the nucleus.
G4bool isNucleusNucleusCollision() const
Is it a nucleus-nucleus collision?
const G4INCL::ThreeVector & getMomentum() const
Get the momentum vector.
Config const * getConfig()
Get the config object.
Definition: G4INCLStore.hh:273
Store * getStore() const
std::string print() const
G4double getEnergy() const
Get the energy of the particle in MeV.
void setINCLMass()
Set the mass of the Particle to its table mass.
#define INCL_WARN(x)
G4double getINCLMass() const
Get the INCL particle mass.
ParticleEntryChannel(Nucleus *n, Particle *p)
int G4int
Definition: G4Types.hh:78
G4double mag2() const
Get the square of the length.
virtual ~ParticleEntryChannel()
G4double getFermiEnergy(const Particle *const p) const
Return the Fermi energy for a particle.
void setEnergy(G4double energy)
Set the energy of the particle in MeV.
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
Final state of an interaction.
#define position
Definition: xmlparse.cc:622
bool G4bool
Definition: G4Types.hh:79
G4int getZ() const
Returns the charge number.
void setPotentialEnergy(G4double v)
Set the particle potential energy.
G4double computeExcitationEnergyExcept(const long exceptID) const
Compute the excitation energy when a nucleon is removed.
void fillFinalState(FinalState *fs)
const G4int n
void setTotalEnergyBeforeInteraction(G4double E)
NuclearPotential::INuclearPotential const * getPotential() const
Getter for thePotential.
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4ThreadLocal NuclearMassFn getTableMass
Static pointer to the mass function for nuclei.
const G4INCL::ThreeVector & getPosition() const
Set the position vector.
Solution solve(RootFunctor const *const f, const G4double x0)
Numerically solve a one-dimensional equation.
G4bool isNucleon() const
Is this a nucleon?
Nucleus * theNucleus
void removeParticle(Particle *const p, const G4double theProjectileCorrection)
Remove a nucleon from the projectile remnant.
Particle * theParticle
virtual G4double getTableMass() const
Get the tabulated particle mass.
ProjectileRemnant * getProjectileRemnant() const
Get the projectile remnant.
G4double getKineticEnergy() const
Get the particle kinetic energy.
virtual G4double computePotentialEnergy(const Particle *const p) const =0
double G4double
Definition: G4Types.hh:76
void addEnteringParticle(Particle *p)
G4double getEmissionQValueCorrection(const G4int AParent, const G4int ZParent) const
Computes correction on the emission Q-value.
#define INCL_DEBUG(x)
const ThreeVector & adjustMomentumFromEnergy()
Rescale the momentum to match the total energy.
Simple class for computing intersections between a straight line and a sphere.
Static root-finder algorithm.
long getID() const
G4bool getRefraction() const
True if we should use refraction.
virtual void setMomentum(const G4INCL::ThreeVector &momentum)
Set the momentum vector.