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 "G4INCLKinematicsUtils.hh" 38 #include "G4INCLKinematicsUtils.hh" 39 #include "G4INCLParticleTable.hh" 39 #include "G4INCLParticleTable.hh" 40 40 41 namespace G4INCL { 41 namespace G4INCL { 42 42 43 namespace KinematicsUtils { 43 namespace KinematicsUtils { 44 44 45 G4double fiveParFit (const G4double a, const 45 G4double fiveParFit (const G4double a, const G4double b, const G4double c, const G4double d, const G4double e, const G4double x){ 46 return a+b*std::pow(x, c)+d*std::log(x)+e* 46 return a+b*std::pow(x, c)+d*std::log(x)+e*std::log(x)*std::log(x); 47 } 47 } 48 48 49 G4double compute_xs(const std::vector<G4doub 49 G4double compute_xs(const std::vector<G4double> coefficients, const G4double pLab){ 50 G4double sigma = 0.; 50 G4double sigma = 0.; 51 G4double Ethreshold = 0.0; 51 G4double Ethreshold = 0.0; 52 if(coefficients.size() == 6){ 52 if(coefficients.size() == 6){ 53 Ethreshold = coefficients[5]; 53 Ethreshold = coefficients[5]; 54 if(Ethreshold >= 5){ //there are no 54 if(Ethreshold >= 5){ //there are no Ethreshold even close to 5 GeV. 55 if(pLab > Ethreshold){ // E is E 55 if(pLab > Ethreshold){ // E is E cutoff, not threshold, we use it when sigma should be zero. 56 return 0.; 56 return 0.; 57 } 57 } 58 } 58 } 59 else{ 59 else{ 60 if(pLab < Ethreshold){ 60 if(pLab < Ethreshold){ 61 return 0.; 61 return 0.; 62 } 62 } 63 } 63 } 64 } 64 } 65 65 66 sigma = fiveParFit(coefficients[0],coeff 66 sigma = fiveParFit(coefficients[0],coefficients[1],coefficients[2],coefficients[3],coefficients[4], pLab); 67 if(sigma < 0.){ 67 if(sigma < 0.){ 68 return 0.; 68 return 0.; 69 }; 69 }; 70 return sigma; 70 return sigma; 71 } 71 } 72 72 73 void transformToLocalEnergyFrame(Nucleus con 73 void transformToLocalEnergyFrame(Nucleus const * const n, Particle * const p) { 74 // assert(!p->isMeson() && !p->isPhoton() && ! 74 // assert(!p->isMeson() && !p->isPhoton() && !p->isAntiNucleon()); // No local energy for mesons //D nor for photons! 75 const G4double localEnergy = getLocalEnerg 75 const G4double localEnergy = getLocalEnergy(n, p); 76 const G4double localTotalEnergy = p->getEn 76 const G4double localTotalEnergy = p->getEnergy() - localEnergy; 77 p->setEnergy(localTotalEnergy); 77 p->setEnergy(localTotalEnergy); 78 p->adjustMomentumFromEnergy(); 78 p->adjustMomentumFromEnergy(); 79 } 79 } 80 80 81 G4double getLocalEnergy(Nucleus const * cons 81 G4double getLocalEnergy(Nucleus const * const n, Particle * const p) { 82 // assert(!p->isMeson() && !p->isPhoton() && ! 82 // assert(!p->isMeson() && !p->isPhoton() && !p->isAntiNucleon()); // No local energy for mesons //D photons are bad too! 83 G4double vloc = 0.0; 83 G4double vloc = 0.0; 84 const G4double r = p->getPosition().mag(); 84 const G4double r = p->getPosition().mag(); 85 const G4double mass = p->getMass(); 85 const G4double mass = p->getMass(); 86 86 87 // Local energy is constant outside the su 87 // Local energy is constant outside the surface 88 if(r > n->getUniverseRadius()) { 88 if(r > n->getUniverseRadius()) { 89 INCL_WARN("Tried to evaluate local energ 89 INCL_WARN("Tried to evaluate local energy for a particle outside the maximum radius." 90 << '\n' << p->print() << '\n' 90 << '\n' << p->print() << '\n' 91 << "Maximum radius = " << n->getDe 91 << "Maximum radius = " << n->getDensity()->getMaximumRadius() << '\n' 92 << "Universe radius = " << n->getU 92 << "Universe radius = " << n->getUniverseRadius() << '\n'); 93 return 0.0; 93 return 0.0; 94 } 94 } 95 95 96 G4double pfl0 = 0.0; 96 G4double pfl0 = 0.0; 97 const ParticleType t = p->getType(); 97 const ParticleType t = p->getType(); 98 const G4double kinE = p->getKineticEnergy( 98 const G4double kinE = p->getKineticEnergy(); 99 if(kinE <= n->getPotential()->getFermiEner 99 if(kinE <= n->getPotential()->getFermiEnergy(t)) { 100 pfl0 = n->getPotential()->getFermiMoment 100 pfl0 = n->getPotential()->getFermiMomentum(p); 101 } else { 101 } else { 102 const G4double tf0 = p->getPotentialEner 102 const G4double tf0 = p->getPotentialEnergy() - n->getPotential()->getSeparationEnergy(p); 103 if(tf0<0.0) return 0.0; 103 if(tf0<0.0) return 0.0; 104 pfl0 = std::sqrt(tf0*(tf0 + 2.0*mass)); 104 pfl0 = std::sqrt(tf0*(tf0 + 2.0*mass)); 105 } 105 } 106 const G4double pReflection = p->getReflect 106 const G4double pReflection = p->getReflectionMomentum()/pfl0; 107 const G4double reflectionRadius = n->getDe 107 const G4double reflectionRadius = n->getDensity()->getMaxRFromP(p->getType(), pReflection); 108 const G4double pNominal = p->getMomentum() 108 const G4double pNominal = p->getMomentum().mag()/pfl0; 109 const G4double nominalReflectionRadius = n 109 const G4double nominalReflectionRadius = n->getDensity()->getMaxRFromP(p->getType(), pNominal); 110 const G4double pl = pfl0*n->getDensity()-> 110 const G4double pl = pfl0*n->getDensity()->getMinPFromR(t, r*nominalReflectionRadius/reflectionRadius); 111 vloc = std::sqrt(pl*pl + mass*mass) - mass 111 vloc = std::sqrt(pl*pl + mass*mass) - mass; 112 112 113 return vloc; 113 return vloc; 114 } 114 } 115 115 116 ThreeVector makeBoostVector(Particle const * 116 ThreeVector makeBoostVector(Particle const * const p1, Particle const * const p2){ 117 const G4double totalEnergy = p1->getEnergy 117 const G4double totalEnergy = p1->getEnergy() + p2->getEnergy(); 118 return ((p1->getMomentum() + p2->getMoment 118 return ((p1->getMomentum() + p2->getMomentum())/totalEnergy); 119 } 119 } 120 120 121 G4double totalEnergyInCM(Particle const * co 121 G4double totalEnergyInCM(Particle const * const p1, Particle const * const p2){ 122 return std::sqrt(squareTotalEnergyInCM(p1, 122 return std::sqrt(squareTotalEnergyInCM(p1,p2)); 123 } 123 } 124 124 125 G4double squareTotalEnergyInCM(Particle cons 125 G4double squareTotalEnergyInCM(Particle const * const p1, Particle const * const p2) { 126 G4double beta2 = makeBoostVector(p1, p2).m 126 G4double beta2 = makeBoostVector(p1, p2).mag2(); 127 if(beta2 > 1.0) { 127 if(beta2 > 1.0) { 128 INCL_ERROR("squareTotalEnergyInCM: beta2 128 INCL_ERROR("squareTotalEnergyInCM: beta2 == " << beta2 << " > 1.0" << '\n'); 129 beta2 = 0.0; 129 beta2 = 0.0; 130 } 130 } 131 return (1.0 - beta2)*std::pow(p1->getEnerg 131 return (1.0 - beta2)*std::pow(p1->getEnergy() + p2->getEnergy(), 2); 132 } 132 } 133 133 134 G4double momentumInCM(Particle const * const 134 G4double momentumInCM(Particle const * const p1, Particle const * const p2) { 135 const G4double m1sq = std::pow(p1->getMass 135 const G4double m1sq = std::pow(p1->getMass(),2); 136 const G4double m2sq = std::pow(p2->getMass 136 const G4double m2sq = std::pow(p2->getMass(),2); 137 const G4double z = p1->getEnergy()*p2->get 137 const G4double z = p1->getEnergy()*p2->getEnergy() - p1->getMomentum().dot(p2->getMomentum()); 138 G4double pcm2 = (z*z-m1sq*m2sq)/(2*z+m1sq+ 138 G4double pcm2 = (z*z-m1sq*m2sq)/(2*z+m1sq+m2sq); 139 if(pcm2 < 0.0) { 139 if(pcm2 < 0.0) { 140 INCL_ERROR("momentumInCM: pcm2 == " << p 140 INCL_ERROR("momentumInCM: pcm2 == " << pcm2 << " < 0.0" << '\n'); 141 pcm2 = 0.0; 141 pcm2 = 0.0; 142 } 142 } 143 return std::sqrt(pcm2); 143 return std::sqrt(pcm2); 144 } 144 } 145 145 146 G4double momentumInCM(const G4double E, cons 146 G4double momentumInCM(const G4double E, const G4double M1, const G4double M2) { 147 return 0.5*std::sqrt((E*E - std::pow(M1 + 147 return 0.5*std::sqrt((E*E - std::pow(M1 + M2, 2)) 148 *(E*E - std::pow(M1 - M2, 2)))/E; 148 *(E*E - std::pow(M1 - M2, 2)))/E; 149 } 149 } 150 150 151 G4double momentumInLab(const G4double s, con 151 G4double momentumInLab(const G4double s, const G4double m1, const G4double m2) { 152 const G4double m1sq = m1*m1; 152 const G4double m1sq = m1*m1; 153 const G4double m2sq = m2*m2; 153 const G4double m2sq = m2*m2; 154 G4double plab2 = (s*s-2*s*(m1sq+m2sq)+(m1s 154 G4double plab2 = (s*s-2*s*(m1sq+m2sq)+(m1sq-m2sq)*(m1sq-m2sq))/(4*m2sq); 155 if(plab2 < 0.0) { 155 if(plab2 < 0.0) { 156 INCL_ERROR("momentumInLab: plab2 == " << 156 INCL_ERROR("momentumInLab: plab2 == " << plab2 << " < 0.0; m1sq == " << m1sq << "; m2sq == " << m2sq << "; s == " << s << '\n'); 157 plab2 = 0.0; 157 plab2 = 0.0; 158 } 158 } 159 return std::sqrt(plab2); 159 return std::sqrt(plab2); 160 } 160 } 161 161 162 G4double momentumInLab(Particle const * cons 162 G4double momentumInLab(Particle const * const p1, Particle const * const p2) { 163 const G4double m1 = p1->getMass(); 163 const G4double m1 = p1->getMass(); 164 const G4double m2 = p2->getMass(); 164 const G4double m2 = p2->getMass(); 165 const G4double s = squareTotalEnergyInCM(p 165 const G4double s = squareTotalEnergyInCM(p1, p2); 166 return momentumInLab(s, m1, m2); 166 return momentumInLab(s, m1, m2); 167 } 167 } 168 168 169 G4double sumTotalEnergies(const ParticleList 169 G4double sumTotalEnergies(const ParticleList &pl) { 170 G4double E = 0.0; 170 G4double E = 0.0; 171 for(ParticleIter i=pl.begin(), e=pl.end(); 171 for(ParticleIter i=pl.begin(), e=pl.end(); i!=e; ++i) { 172 E += (*i)->getEnergy(); 172 E += (*i)->getEnergy(); 173 } 173 } 174 return E; 174 return E; 175 } 175 } 176 176 177 ThreeVector sumMomenta(const ParticleList &p 177 ThreeVector sumMomenta(const ParticleList &pl) { 178 ThreeVector p(0.0, 0.0, 0.0); 178 ThreeVector p(0.0, 0.0, 0.0); 179 for(ParticleIter i=pl.begin(), e=pl.end(); 179 for(ParticleIter i=pl.begin(), e=pl.end(); i!=e; ++i) { 180 p += (*i)->getMomentum(); 180 p += (*i)->getMomentum(); 181 } 181 } 182 return p; 182 return p; 183 } 183 } 184 184 185 G4double energy(const ThreeVector &p, const 185 G4double energy(const ThreeVector &p, const G4double m) { 186 return std::sqrt(p.mag2() + m*m); 186 return std::sqrt(p.mag2() + m*m); 187 } 187 } 188 188 189 G4double invariantMass(const G4double E, con 189 G4double invariantMass(const G4double E, const ThreeVector & p) { 190 return std::sqrt(squareInvariantMass(E, p) 190 return std::sqrt(squareInvariantMass(E, p)); 191 } 191 } 192 192 193 G4double squareInvariantMass(const G4double 193 G4double squareInvariantMass(const G4double E, const ThreeVector & p) { 194 return E*E - p.mag2(); 194 return E*E - p.mag2(); 195 } 195 } 196 196 197 G4double gammaFromKineticEnergy(const Partic 197 G4double gammaFromKineticEnergy(const ParticleSpecies &p, const G4double EKin) { 198 G4double mass; 198 G4double mass; 199 if(p.theType==Composite) 199 if(p.theType==Composite) 200 mass = ParticleTable::getTableMass(p.the 200 mass = ParticleTable::getTableMass(p.theA, p.theZ, p.theS); 201 else 201 else 202 mass = ParticleTable::getTableParticleMa 202 mass = ParticleTable::getTableParticleMass(p.theType); 203 return (1.+EKin/mass); 203 return (1.+EKin/mass); 204 } 204 } 205 205 206 } 206 } 207 207 208 } 208 } 209 209