ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLParticleEntryChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 particle
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->isNucleonorLambda()); // Possible hypernucleus projectile of inverse kinematic
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(), projectileRemnant->getS() - theParticle->getS())
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(
138  + theProjectileCorrection;
139  // end of part common to model 3. and 4.
140 
141 
142  projectileRemnant->removeParticle(theParticle, theProjectileCorrection);
143  } else {
144  const G4int ACN = theNucleus->getA() + theParticle->getA();
145  const G4int ZCN = theNucleus->getZ() + theParticle->getZ();
146  const G4int SCN = theNucleus->getS() + theParticle->getS();
147  // Correction to the Q-value of the entering particle
148  if(theParticle->isKaon()) theCorrection = theParticle->getEmissionQValueCorrection(ACN,ZCN,theNucleus->getS());
149  else theCorrection = theParticle->getEmissionQValueCorrection(ACN,ZCN,SCN);
150  INCL_DEBUG("The following Particle enters with correction " << theCorrection << '\n'
151  << theParticle->print() << '\n');
152  }
153 
154  const G4double energyBefore = theParticle->getEnergy() - theCorrection;
155  G4bool success = particleEnters(theCorrection);
157 
158  if(!success) {
159  fs->makeParticleBelowZero();
160  } else if(theParticle->isNucleonorLambda() &&
162  // If the participant is a nucleon entering below its Fermi energy, force a
163  // compound nucleus
166 
167  fs->setTotalEnergyBeforeInteraction(energyBefore);
168  }
169 
171 
172  // \todo{this is the place to add refraction}
173 
174  theParticle->setINCLMass(); // Will automatically put the particle on shell
175 
176  // Add the nuclear potential to the kinetic energy when entering the
177  // nucleus
178 
179  class IncomingEFunctor : public RootFunctor {
180  public:
181  IncomingEFunctor(Particle * const p, Nucleus const * const n, const G4double correction) :
182  RootFunctor(0., 1E6),
183  theParticle(p),
184  thePotential(n->getPotential()),
185  theEnergy(theParticle->getEnergy()),
186  theMass(theParticle->getMass()),
187  theQValueCorrection(correction),
188  refraction(n->getStore()->getConfig()->getRefraction()),
189  theMomentumDirection(theParticle->getMomentum())
190  {
191  if(refraction) {
193  const G4double r2 = position.mag2();
194  if(r2>0.)
195  normal = - position / std::sqrt(r2);
196  G4double cosIncidenceAngle = theParticle->getCosRPAngle();
197  if(cosIncidenceAngle < -1.)
198  sinIncidenceAnglePOut = 0.;
199  else
200  sinIncidenceAnglePOut = theMomentumDirection.mag()*std::sqrt(1.-cosIncidenceAngle*cosIncidenceAngle);
201  } else {
202  sinIncidenceAnglePOut = 0.;
203  }
204  }
205  ~IncomingEFunctor() {}
206  G4double operator()(const G4double v) const {
207  G4double energyInside = std::max(theMass, theEnergy + v - theQValueCorrection);
208  theParticle->setEnergy(energyInside);
210  if(refraction) {
211  // Compute the new direction of the particle momentum
212  const G4double pIn = std::sqrt(energyInside*energyInside-theMass*theMass);
213  const G4double sinRefractionAngle = sinIncidenceAnglePOut/pIn;
214  const G4double cosRefractionAngle = (sinRefractionAngle>1.) ? 0. : std::sqrt(1.-sinRefractionAngle*sinRefractionAngle);
215  const ThreeVector momentumInside = theMomentumDirection - normal * normal.dot(theMomentumDirection) + normal * (pIn * cosRefractionAngle);
216  theParticle->setMomentum(momentumInside);
217  } else {
218  theParticle->setMomentum(theMomentumDirection); // keep the same direction
219  }
220  // Scale the particle momentum
222  return v - thePotential->computePotentialEnergy(theParticle);
223  }
224  void cleanUp(const G4bool /*success*/) const {}
225  private:
227  NuclearPotential::INuclearPotential const *thePotential;
228  const G4double theEnergy;
229  const G4double theMass;
230  const G4double theQValueCorrection;
231  const G4bool refraction;
232  const ThreeVector theMomentumDirection;
234  G4double sinIncidenceAnglePOut;
235  } theIncomingEFunctor(theParticle,theNucleus,theQValueCorrection);
236 
238  if(theParticle->getKineticEnergy()+v-theQValueCorrection<0.) { // Particle entering below 0. Die gracefully
239  INCL_DEBUG("Particle " << theParticle->getID() << " is trying to enter below 0" << '\n');
240  return false;
241  }
242  const RootFinder::Solution theSolution = RootFinder::solve(&theIncomingEFunctor, v);
243  if(theSolution.success) { // Apply the solution
244  theIncomingEFunctor(theSolution.x);
245  INCL_DEBUG("Particle successfully entered:\n" << theParticle->print() << '\n');
246  } else {
247  INCL_WARN("Couldn't compute the potential for incoming particle, root-finding algorithm failed." << '\n');
248  }
249  return theSolution.success;
250  }
251 
252 }
253