| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 H
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(N
46 :theNucleus(n), theParticle(p)
47 {}
48
49 ParticleEntryChannel::~ParticleEntryChannel(
50 {}
51
52 void ParticleEntryChannel::fillFinalState(Fi
53 // Behaves slightly differency if a third
54 G4bool isNN = theNucleus->isNucleusNucleus
55
56 /* Corrections to the energy of the enteri
57 *
58 * In particle-nucleus reactions, the goal
59 * energy conservation in particle-nucleus
60 * and nuclear masses.
61 *
62 * In nucleus-nucleus reactions, in additi
63 * is determined by a model for the excita
64 * quasi-projectile (QP). The energy of th
65 * the QP excitation energy, as determined
66 * by our model.
67 *
68 * Possible choices for the correction (or
69 * excitation energy):
70 *
71 * 1. the correction is 0. (same as in par
72 * 2. the correction is the separation ene
73 * the current QP;
74 * 3. the QP excitation energy is given by
75 * implemented in INCL4.2-HI/Geant4.
76 * 4. the QP excitation energy vanishes.
77 *
78 * Ideally, the QP excitation energy shoul
79 * and 2. do not guarantee this, although
80 * more severe for 1. than for 2.. Algorit
81 * yields non-negative QP excitation energ
82 */
83 G4double theCorrection;
84 if(isNN) {
85 // assert(theParticle->isNucleonorLambda()); /
86 ProjectileRemnant * const projectileRemn
87 // assert(projectileRemnant);
88
89 // No correction (model 1. above)
90 /*
91 theCorrection = theParticle->getEmission
92 theNucleus->getA() + theParticle->ge
93 theNucleus->getZ() + theParticle->ge
94 + theParticle->getTableMass() - thePar
95 const G4double theProjectileCorrection =
96 */
97
98 // Correct the energy of the entering pa
99 // emission from the projectile (model 2
100 /*
101 theCorrection = theParticle->getTransfer
102 projectileRemnant->getA(), projectil
103 theNucleus->getA(), theNucleus->getZ
104 G4double theProjectileCorrection;
105 if(projectileRemnant->getA()>theParticle
106 // Compute the projectile Q-value (to
107 // other components of the projectile
108 theProjectileCorrection = ParticleTabl
109 projectileRemnant->getA() - thePar
110 projectileRemnant->getZ() - thePar
111 theParticle->getA(),
112 theParticle->getZ());
113 } else
114 theProjectileCorrection = 0.;
115 */
116
117 // Fix the correction in such a way that
118 // energy is given by A. Boudard's INCL4
119 const G4double theProjectileExcitationEn
120 (projectileRemnant->getA()-theParticle
121 (projectileRemnant->computeExcitationE
122 0.;
123 // Set the projectile excitation energy
124 // model 4. above).
125 // const G4double theProjectileExcitatio
126 // The part that follows is common to mo
127 const G4double theProjectileEffectiveMas
128 ParticleTable::getTableMass(projectile
129 + theProjectileExcitationEnergy;
130 const ThreeVector &theProjectileMomentum
131 const G4double theProjectileEnergy = std
132 const G4double theProjectileCorrection =
133 theCorrection = theParticle->getEmission
134 theNucleus->getA() + theParticle->ge
135 theNucleus->getZ() + theParticle->ge
136 theNucleus->getS() + theParticle->ge
137 + theParticle->getTableMass() - thePar
138 + theProjectileCorrection;
139 // end of part common to model 3. and 4.
140
141
142 projectileRemnant->removeParticle(thePar
143 } else {
144 const G4int ACN = theNucleus->getA() + t
145 const G4int ZCN = theNucleus->getZ() + t
146 const G4int SCN = theNucleus->getS() + t
147 // Correction to the Q-value of the ente
148 if(theParticle->isKaon()) theCorrection
149 else theCorrection = theParticle->getEmi
150 INCL_DEBUG("The following Particle enter
151 << theParticle->print() << '\n');
152 }
153
154 const G4double energyBefore = theParticle-
155 G4bool success = particleEnters(theCorrect
156 fs->addEnteringParticle(theParticle);
157
158 if(!success) {
159 fs->makeParticleBelowZero();
160 } else if(theParticle->isNucleonorLambda()
161 theParticle->getKineticEnergy()<theNuc
162 // If the participant is a nucleon enter
163 // compound nucleus
164 fs->makeParticleBelowFermi();
165 } else if(theParticle->isKaon()) theNucleu
166
167 fs->setTotalEnergyBeforeInteraction(energy
168 }
169
170 G4bool ParticleEntryChannel::particleEnters(
171
172 // \todo{this is the place to add refracti
173
174 theParticle->setINCLMass(); // Will automa
175
176 // Add the nuclear potential to the kineti
177 // nucleus
178
179 class IncomingEFunctor : public RootFuncto
180 public:
181 IncomingEFunctor(Particle * const p, N
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(
189 theMomentumDirection(theParticle->ge
190 {
191 if(refraction) {
192 const ThreeVector &position = th
193 const G4double r2 = position.mag
194 if(r2>0.)
195 normal = - position / std::sqr
196 G4double cosIncidenceAngle = the
197 if(cosIncidenceAngle < -1.)
198 sinIncidenceAnglePOut = 0.;
199 else
200 sinIncidenceAnglePOut = theMom
201 } else {
202 sinIncidenceAnglePOut = 0.;
203 }
204 }
205 ~IncomingEFunctor() {}
206 G4double operator()(const G4double v)
207 G4double energyInside = std::max(the
208 theParticle->setEnergy(energyInside)
209 theParticle->setPotentialEnergy(v);
210 if(refraction) {
211 // Compute the new direction of th
212 const G4double pIn = std::sqrt(ene
213 const G4double sinRefractionAngle
214 const G4double cosRefractionAngle
215 const ThreeVector momentumInside =
216 theParticle->setMomentum(momentumI
217 } else {
218 theParticle->setMomentum(theMoment
219 }
220 // Scale the particle momentum
221 theParticle->adjustMomentumFromEnerg
222 return v - thePotential->computePote
223 }
224 void cleanUp(const G4bool /*success*/)
225 private:
226 Particle *theParticle;
227 NuclearPotential::INuclearPotential co
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,theNucle
236
237 G4double v = theNucleus->getPotential()->c
238 if(theParticle->getKineticEnergy()+v-theQV
239 INCL_DEBUG("Particle " << theParticle->g
240 return false;
241 }
242 const RootFinder::Solution theSolution = R
243 if(theSolution.success) { // Apply the sol
244 theIncomingEFunctor(theSolution.x);
245 INCL_DEBUG("Particle successfully entere
246 } else {
247 INCL_WARN("Couldn't compute the potentia
248 }
249 return theSolution.success;
250 }
251
252 }
253
254