Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4INCLParticleEntryChannel.cc
이 파일의 문서화 페이지로 가기
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 
NuclearPotential::INuclearPotential const * getPotential() const
Getter for thePotential.
T max(const T t1, const T t2)
brief Return the largest of the two arguments
virtual G4double getTableMass() const
Get the tabulated particle mass.
const G4INCL::ThreeVector & getPosition() const
G4double getCosRPAngle() const
Get the cosine of the angle between position and momentum.
Static root-finder algorithm.
G4double getEnergy() const
const ThreeVector & adjustMomentumFromEnergy()
Rescale the momentum to match the total energy.
ParticleEntryChannel(Nucleus *n, Particle *p)
G4double getKineticEnergy() const
Get the particle kinetic energy.
const G4INCL::ThreeVector & getMomentum() const
const char * p
Definition: xmltok.h:285
long getID() const
ProjectileRemnant * getProjectileRemnant() const
Get the projectile remnant.
void setINCLMass()
Set the mass of the Particle to its table mass.
G4double computeExcitationEnergyExcept(const long exceptID) const
Compute the excitation energy when a nucleon is removed.
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
G4double getMass() const
Get the cached particle mass.
#define INCL_DEBUG(x)
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
Solution solve(RootFunctor const *const f, const G4double x0)
Numerically solve a one-dimensional equation.
G4ThreadLocal NuclearMassFn getTableMass
Static pointer to the mass function for nuclei.
G4bool getRefraction() const
True if we should use refraction.
void setTotalEnergyBeforeInteraction(G4double E)
void setPotentialEnergy(G4double v)
Set the particle potential energy.
std::string print() const
G4int getZ() const
Returns the charge number.
G4int getA() const
Returns the baryon number.
G4double getINCLMass() const
Get the INCL particle mass.
virtual void setMomentum(const G4INCL::ThreeVector &momentum)
G4bool particleEnters(const G4double theQValueCorrection)
Modify particle that enters the nucleus.
#define INCL_WARN(x)
int G4int
Definition: G4Types.hh:78
Config const * getConfig()
Definition: G4INCLStore.hh:273
Simple class for computing intersections between a straight line and a sphere.
G4double mag2() const
G4double getEmissionQValueCorrection(const G4int AParent, const G4int ZParent) const
Computes correction on the emission Q-value.
void setEnergy(G4double energy)
void addEnteringParticle(Particle *p)
void removeParticle(Particle *const p, const G4double theProjectileCorrection)
Remove a nucleon from the projectile remnant.
G4double getFermiEnergy(const Particle *const p) const
Return the Fermi energy for a particle.
Char_t n[5]
G4bool isNucleon() const
virtual G4double computePotentialEnergy(const Particle *const p) const =0
Store * getStore() const
G4bool isNucleusNucleusCollision() const
Is it a nucleus-nucleus collision?