| Geant4 Cross Reference |
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 #include "G4INCLParticleEntryChannel.hh"
39 #include "G4INCLRootFinder.hh"
40 #include "G4INCLIntersection.hh"
41 #include <algorithm>
42
43 namespace G4INCL {
44
45 ParticleEntryChannel::ParticleEntryChannel(Nucleus *n, Particle *p)
46 :theNucleus(n), theParticle(p)
47 {}
48
49 ParticleEntryChannel::~ParticleEntryChannel()
50 {}
51
52 void ParticleEntryChannel::fillFinalState(FinalState *fs) {
53 // Behaves slightly differency if a third body (the projectile) is present
54 G4bool isNN = theNucleus->isNucleusNucleusCollision();
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(
134 theNucleus->getA() + theParticle->getA(),
135 theNucleus->getZ() + theParticle->getZ(),
136 theNucleus->getS() + theParticle->getS())
137 + theParticle->getTableMass() - theParticle->getINCLMass()
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);
156 fs->addEnteringParticle(theParticle);
157
158 if(!success) {
159 fs->makeParticleBelowZero();
160 } else if(theParticle->isNucleonorLambda() &&
161 theParticle->getKineticEnergy()<theNucleus->getPotential()->getFermiEnergy(theParticle)) {
162 // If the participant is a nucleon entering below its Fermi energy, force a
163 // compound nucleus
164 fs->makeParticleBelowFermi();
165 } else if(theParticle->isKaon()) theNucleus->setNumberOfKaon(theNucleus->getNumberOfKaon()+1);
166
167 fs->setTotalEnergyBeforeInteraction(energyBefore);
168 }
169
170 G4bool ParticleEntryChannel::particleEnters(const G4double theQValueCorrection) {
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) {
192 const ThreeVector &position = theParticle->getPosition();
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);
209 theParticle->setPotentialEnergy(v);
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
221 theParticle->adjustMomentumFromEnergy();
222 return v - thePotential->computePotentialEnergy(theParticle);
223 }
224 void cleanUp(const G4bool /*success*/) const {}
225 private:
226 Particle *theParticle;
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;
233 ThreeVector normal;
234 G4double sinIncidenceAnglePOut;
235 } theIncomingEFunctor(theParticle,theNucleus,theQValueCorrection);
236
237 G4double v = theNucleus->getPotential()->computePotentialEnergy(theParticle);
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
254