| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // INCL++ intra-nuclear cascade model 26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France 27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium 28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France 29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and H 30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France 31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France 32 // Davide Mancusi, CEA-Saclay, France
33 // 33 //
34 #define INCLXX_IN_GEANT4_MODE 1 34 #define INCLXX_IN_GEANT4_MODE 1
35 35
36 #include "globals.hh" 36 #include "globals.hh"
37 37
38 #include "G4INCLParticleEntryChannel.hh" 38 #include "G4INCLParticleEntryChannel.hh"
39 #include "G4INCLRootFinder.hh" 39 #include "G4INCLRootFinder.hh"
40 #include "G4INCLIntersection.hh" 40 #include "G4INCLIntersection.hh"
41 #include <algorithm> 41 #include <algorithm>
42 42
43 namespace G4INCL { 43 namespace G4INCL {
44 44
45 ParticleEntryChannel::ParticleEntryChannel(N 45 ParticleEntryChannel::ParticleEntryChannel(Nucleus *n, Particle *p)
46 :theNucleus(n), theParticle(p) 46 :theNucleus(n), theParticle(p)
47 {} 47 {}
48 48
49 ParticleEntryChannel::~ParticleEntryChannel( 49 ParticleEntryChannel::~ParticleEntryChannel()
50 {} 50 {}
51 51
52 void ParticleEntryChannel::fillFinalState(Fi 52 void ParticleEntryChannel::fillFinalState(FinalState *fs) {
53 // Behaves slightly differency if a third 53 // Behaves slightly differency if a third body (the projectile) is present
54 G4bool isNN = theNucleus->isNucleusNucleus 54 G4bool isNN = theNucleus->isNucleusNucleusCollision();
55 55
56 /* Corrections to the energy of the enteri << 56 /* Corrections to the energy of the entering nucleon
57 * 57 *
58 * In particle-nucleus reactions, the goal 58 * In particle-nucleus reactions, the goal of this correction is to satisfy
59 * energy conservation in particle-nucleus 59 * energy conservation in particle-nucleus reactions using real particle
60 * and nuclear masses. 60 * and nuclear masses.
61 * 61 *
62 * In nucleus-nucleus reactions, in additi 62 * In nucleus-nucleus reactions, in addition to the above, the correction
63 * is determined by a model for the excita 63 * is determined by a model for the excitation energy of the
64 * quasi-projectile (QP). The energy of th 64 * quasi-projectile (QP). The energy of the entering nucleon is such that
65 * the QP excitation energy, as determined 65 * the QP excitation energy, as determined by conservation, is what given
66 * by our model. 66 * by our model.
67 * 67 *
68 * Possible choices for the correction (or 68 * Possible choices for the correction (or, equivalently, for the QP
69 * excitation energy): 69 * excitation energy):
70 * 70 *
71 * 1. the correction is 0. (same as in par 71 * 1. the correction is 0. (same as in particle-nucleus);
72 * 2. the correction is the separation ene 72 * 2. the correction is the separation energy of the entering nucleon in
73 * the current QP; 73 * the current QP;
74 * 3. the QP excitation energy is given by 74 * 3. the QP excitation energy is given by A. Boudard's algorithm, as
75 * implemented in INCL4.2-HI/Geant4. 75 * implemented in INCL4.2-HI/Geant4.
76 * 4. the QP excitation energy vanishes. 76 * 4. the QP excitation energy vanishes.
77 * 77 *
78 * Ideally, the QP excitation energy shoul 78 * Ideally, the QP excitation energy should always be >=0. Algorithms 1.
79 * and 2. do not guarantee this, although 79 * and 2. do not guarantee this, although violations to the rule seem to be
80 * more severe for 1. than for 2.. Algorit 80 * more severe for 1. than for 2.. Algorithms 3. and 4., by construction,
81 * yields non-negative QP excitation energ 81 * yields non-negative QP excitation energies.
82 */ 82 */
83 G4double theCorrection; 83 G4double theCorrection;
84 if(isNN) { 84 if(isNN) {
85 // assert(theParticle->isNucleonorLambda()); / << 85 // assert(theParticle->isNucleon());
86 ProjectileRemnant * const projectileRemn 86 ProjectileRemnant * const projectileRemnant = theNucleus->getProjectileRemnant();
87 // assert(projectileRemnant); 87 // assert(projectileRemnant);
88 88
89 // No correction (model 1. above) 89 // No correction (model 1. above)
90 /* 90 /*
91 theCorrection = theParticle->getEmission 91 theCorrection = theParticle->getEmissionQValueCorrection(
92 theNucleus->getA() + theParticle->ge 92 theNucleus->getA() + theParticle->getA(),
93 theNucleus->getZ() + theParticle->ge 93 theNucleus->getZ() + theParticle->getZ())
94 + theParticle->getTableMass() - thePar 94 + theParticle->getTableMass() - theParticle->getINCLMass();
95 const G4double theProjectileCorrection = 95 const G4double theProjectileCorrection = 0.;
96 */ 96 */
97 97
98 // Correct the energy of the entering pa 98 // Correct the energy of the entering particle for the Q-value of the
99 // emission from the projectile (model 2 99 // emission from the projectile (model 2. above)
100 /* 100 /*
101 theCorrection = theParticle->getTransfer 101 theCorrection = theParticle->getTransferQValueCorrection(
102 projectileRemnant->getA(), projectil 102 projectileRemnant->getA(), projectileRemnant->getZ(),
103 theNucleus->getA(), theNucleus->getZ 103 theNucleus->getA(), theNucleus->getZ());
104 G4double theProjectileCorrection; 104 G4double theProjectileCorrection;
105 if(projectileRemnant->getA()>theParticle 105 if(projectileRemnant->getA()>theParticle->getA()) { // if there are any particles left
106 // Compute the projectile Q-value (to 106 // Compute the projectile Q-value (to be used as a correction to the
107 // other components of the projectile 107 // other components of the projectile remnant)
108 theProjectileCorrection = ParticleTabl 108 theProjectileCorrection = ParticleTable::getTableQValue(
109 projectileRemnant->getA() - thePar 109 projectileRemnant->getA() - theParticle->getA(),
110 projectileRemnant->getZ() - thePar 110 projectileRemnant->getZ() - theParticle->getZ(),
111 theParticle->getA(), 111 theParticle->getA(),
112 theParticle->getZ()); 112 theParticle->getZ());
113 } else 113 } else
114 theProjectileCorrection = 0.; 114 theProjectileCorrection = 0.;
115 */ 115 */
116 116
117 // Fix the correction in such a way that 117 // Fix the correction in such a way that the quasi-projectile excitation
118 // energy is given by A. Boudard's INCL4 118 // energy is given by A. Boudard's INCL4.2-HI model (model 3. above).
119 const G4double theProjectileExcitationEn 119 const G4double theProjectileExcitationEnergy =
120 (projectileRemnant->getA()-theParticle 120 (projectileRemnant->getA()-theParticle->getA()>1) ?
121 (projectileRemnant->computeExcitationE 121 (projectileRemnant->computeExcitationEnergyExcept(theParticle->getID())) :
122 0.; 122 0.;
123 // Set the projectile excitation energy 123 // Set the projectile excitation energy to zero (cold quasi-projectile,
124 // model 4. above). 124 // model 4. above).
125 // const G4double theProjectileExcitatio 125 // const G4double theProjectileExcitationEnergy = 0.;
126 // The part that follows is common to mo 126 // The part that follows is common to model 3. and 4.
127 const G4double theProjectileEffectiveMas 127 const G4double theProjectileEffectiveMass =
128 ParticleTable::getTableMass(projectile << 128 ParticleTable::getTableMass(projectileRemnant->getA() - theParticle->getA(), projectileRemnant->getZ() - theParticle->getZ())
129 + theProjectileExcitationEnergy; 129 + theProjectileExcitationEnergy;
130 const ThreeVector &theProjectileMomentum 130 const ThreeVector &theProjectileMomentum = projectileRemnant->getMomentum() - theParticle->getMomentum();
131 const G4double theProjectileEnergy = std 131 const G4double theProjectileEnergy = std::sqrt(theProjectileMomentum.mag2() + theProjectileEffectiveMass*theProjectileEffectiveMass);
132 const G4double theProjectileCorrection = 132 const G4double theProjectileCorrection = theProjectileEnergy - (projectileRemnant->getEnergy() - theParticle->getEnergy());
133 theCorrection = theParticle->getEmission 133 theCorrection = theParticle->getEmissionQValueCorrection(
134 theNucleus->getA() + theParticle->ge 134 theNucleus->getA() + theParticle->getA(),
135 theNucleus->getZ() + theParticle->ge << 135 theNucleus->getZ() + theParticle->getZ())
136 theNucleus->getS() + theParticle->ge <<
137 + theParticle->getTableMass() - thePar 136 + theParticle->getTableMass() - theParticle->getINCLMass()
138 + theProjectileCorrection; 137 + theProjectileCorrection;
139 // end of part common to model 3. and 4. 138 // end of part common to model 3. and 4.
140 139
141 140
142 projectileRemnant->removeParticle(thePar 141 projectileRemnant->removeParticle(theParticle, theProjectileCorrection);
143 } else { 142 } else {
144 const G4int ACN = theNucleus->getA() + t 143 const G4int ACN = theNucleus->getA() + theParticle->getA();
145 const G4int ZCN = theNucleus->getZ() + t 144 const G4int ZCN = theNucleus->getZ() + theParticle->getZ();
146 const G4int SCN = theNucleus->getS() + t <<
147 // Correction to the Q-value of the ente 145 // Correction to the Q-value of the entering particle
148 if(theParticle->isKaon()) theCorrection << 146 theCorrection = theParticle->getEmissionQValueCorrection(ACN,ZCN);
149 else theCorrection = theParticle->getEmi <<
150 INCL_DEBUG("The following Particle enter 147 INCL_DEBUG("The following Particle enters with correction " << theCorrection << '\n'
151 << theParticle->print() << '\n'); 148 << theParticle->print() << '\n');
152 } 149 }
153 150
154 const G4double energyBefore = theParticle- 151 const G4double energyBefore = theParticle->getEnergy() - theCorrection;
155 G4bool success = particleEnters(theCorrect 152 G4bool success = particleEnters(theCorrection);
156 fs->addEnteringParticle(theParticle); 153 fs->addEnteringParticle(theParticle);
157 154
158 if(!success) { 155 if(!success) {
159 fs->makeParticleBelowZero(); 156 fs->makeParticleBelowZero();
160 } else if(theParticle->isNucleonorLambda() << 157 } else if(theParticle->isNucleon() &&
161 theParticle->getKineticEnergy()<theNuc 158 theParticle->getKineticEnergy()<theNucleus->getPotential()->getFermiEnergy(theParticle)) {
162 // If the participant is a nucleon enter 159 // If the participant is a nucleon entering below its Fermi energy, force a
163 // compound nucleus 160 // compound nucleus
164 fs->makeParticleBelowFermi(); 161 fs->makeParticleBelowFermi();
165 } else if(theParticle->isKaon()) theNucleu << 162 }
166 163
167 fs->setTotalEnergyBeforeInteraction(energy 164 fs->setTotalEnergyBeforeInteraction(energyBefore);
168 } 165 }
169 166
170 G4bool ParticleEntryChannel::particleEnters( 167 G4bool ParticleEntryChannel::particleEnters(const G4double theQValueCorrection) {
171 168
172 // \todo{this is the place to add refracti 169 // \todo{this is the place to add refraction}
173 170
174 theParticle->setINCLMass(); // Will automa 171 theParticle->setINCLMass(); // Will automatically put the particle on shell
175 172
176 // Add the nuclear potential to the kineti 173 // Add the nuclear potential to the kinetic energy when entering the
177 // nucleus 174 // nucleus
178 175
179 class IncomingEFunctor : public RootFuncto 176 class IncomingEFunctor : public RootFunctor {
180 public: 177 public:
181 IncomingEFunctor(Particle * const p, N 178 IncomingEFunctor(Particle * const p, Nucleus const * const n, const G4double correction) :
182 RootFunctor(0., 1E6), 179 RootFunctor(0., 1E6),
183 theParticle(p), 180 theParticle(p),
184 thePotential(n->getPotential()), 181 thePotential(n->getPotential()),
185 theEnergy(theParticle->getEnergy()), 182 theEnergy(theParticle->getEnergy()),
186 theMass(theParticle->getMass()), 183 theMass(theParticle->getMass()),
187 theQValueCorrection(correction), 184 theQValueCorrection(correction),
188 refraction(n->getStore()->getConfig( 185 refraction(n->getStore()->getConfig()->getRefraction()),
189 theMomentumDirection(theParticle->ge 186 theMomentumDirection(theParticle->getMomentum())
190 { 187 {
191 if(refraction) { 188 if(refraction) {
192 const ThreeVector &position = th 189 const ThreeVector &position = theParticle->getPosition();
193 const G4double r2 = position.mag 190 const G4double r2 = position.mag2();
194 if(r2>0.) 191 if(r2>0.)
195 normal = - position / std::sqr 192 normal = - position / std::sqrt(r2);
196 G4double cosIncidenceAngle = the 193 G4double cosIncidenceAngle = theParticle->getCosRPAngle();
197 if(cosIncidenceAngle < -1.) 194 if(cosIncidenceAngle < -1.)
198 sinIncidenceAnglePOut = 0.; 195 sinIncidenceAnglePOut = 0.;
199 else 196 else
200 sinIncidenceAnglePOut = theMom 197 sinIncidenceAnglePOut = theMomentumDirection.mag()*std::sqrt(1.-cosIncidenceAngle*cosIncidenceAngle);
201 } else { 198 } else {
202 sinIncidenceAnglePOut = 0.; 199 sinIncidenceAnglePOut = 0.;
203 } 200 }
204 } 201 }
205 ~IncomingEFunctor() {} 202 ~IncomingEFunctor() {}
206 G4double operator()(const G4double v) 203 G4double operator()(const G4double v) const {
207 G4double energyInside = std::max(the 204 G4double energyInside = std::max(theMass, theEnergy + v - theQValueCorrection);
208 theParticle->setEnergy(energyInside) 205 theParticle->setEnergy(energyInside);
209 theParticle->setPotentialEnergy(v); 206 theParticle->setPotentialEnergy(v);
210 if(refraction) { 207 if(refraction) {
211 // Compute the new direction of th 208 // Compute the new direction of the particle momentum
212 const G4double pIn = std::sqrt(ene 209 const G4double pIn = std::sqrt(energyInside*energyInside-theMass*theMass);
213 const G4double sinRefractionAngle 210 const G4double sinRefractionAngle = sinIncidenceAnglePOut/pIn;
214 const G4double cosRefractionAngle 211 const G4double cosRefractionAngle = (sinRefractionAngle>1.) ? 0. : std::sqrt(1.-sinRefractionAngle*sinRefractionAngle);
215 const ThreeVector momentumInside = 212 const ThreeVector momentumInside = theMomentumDirection - normal * normal.dot(theMomentumDirection) + normal * (pIn * cosRefractionAngle);
216 theParticle->setMomentum(momentumI 213 theParticle->setMomentum(momentumInside);
217 } else { 214 } else {
218 theParticle->setMomentum(theMoment 215 theParticle->setMomentum(theMomentumDirection); // keep the same direction
219 } 216 }
220 // Scale the particle momentum 217 // Scale the particle momentum
221 theParticle->adjustMomentumFromEnerg 218 theParticle->adjustMomentumFromEnergy();
222 return v - thePotential->computePote 219 return v - thePotential->computePotentialEnergy(theParticle);
223 } 220 }
224 void cleanUp(const G4bool /*success*/) 221 void cleanUp(const G4bool /*success*/) const {}
225 private: 222 private:
226 Particle *theParticle; 223 Particle *theParticle;
227 NuclearPotential::INuclearPotential co 224 NuclearPotential::INuclearPotential const *thePotential;
228 const G4double theEnergy; 225 const G4double theEnergy;
229 const G4double theMass; 226 const G4double theMass;
230 const G4double theQValueCorrection; 227 const G4double theQValueCorrection;
231 const G4bool refraction; 228 const G4bool refraction;
232 const ThreeVector theMomentumDirection 229 const ThreeVector theMomentumDirection;
233 ThreeVector normal; 230 ThreeVector normal;
234 G4double sinIncidenceAnglePOut; 231 G4double sinIncidenceAnglePOut;
235 } theIncomingEFunctor(theParticle,theNucle 232 } theIncomingEFunctor(theParticle,theNucleus,theQValueCorrection);
236 233
237 G4double v = theNucleus->getPotential()->c 234 G4double v = theNucleus->getPotential()->computePotentialEnergy(theParticle);
238 if(theParticle->getKineticEnergy()+v-theQV 235 if(theParticle->getKineticEnergy()+v-theQValueCorrection<0.) { // Particle entering below 0. Die gracefully
239 INCL_DEBUG("Particle " << theParticle->g 236 INCL_DEBUG("Particle " << theParticle->getID() << " is trying to enter below 0" << '\n');
240 return false; 237 return false;
241 } 238 }
242 const RootFinder::Solution theSolution = R 239 const RootFinder::Solution theSolution = RootFinder::solve(&theIncomingEFunctor, v);
243 if(theSolution.success) { // Apply the sol 240 if(theSolution.success) { // Apply the solution
244 theIncomingEFunctor(theSolution.x); 241 theIncomingEFunctor(theSolution.x);
245 INCL_DEBUG("Particle successfully entere 242 INCL_DEBUG("Particle successfully entered:\n" << theParticle->print() << '\n');
246 } else { 243 } else {
247 INCL_WARN("Couldn't compute the potentia 244 INCL_WARN("Couldn't compute the potential for incoming particle, root-finding algorithm failed." << '\n');
248 } 245 }
249 return theSolution.success; 246 return theSolution.success;
250 } 247 }
251 248
252 } 249 }
253 250
254 251