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