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