| 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 "G4INCLKinematicsUtils.hh"
39 #include "G4INCLParticleTable.hh"
40
41 namespace G4INCL {
42
43 namespace KinematicsUtils {
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
74 // assert(!p->isMeson() && !p->isPhoton() && !
75 const G4double localEnergy = getLocalEnerg
76 const G4double localTotalEnergy = p->getEn
77 p->setEnergy(localTotalEnergy);
78 p->adjustMomentumFromEnergy();
79 }
80
81 G4double getLocalEnergy(Nucleus const * cons
82 // assert(!p->isMeson() && !p->isPhoton() && !
83 G4double vloc = 0.0;
84 const G4double r = p->getPosition().mag();
85 const G4double mass = p->getMass();
86
87 // Local energy is constant outside the su
88 if(r > n->getUniverseRadius()) {
89 INCL_WARN("Tried to evaluate local energ
90 << '\n' << p->print() << '\n'
91 << "Maximum radius = " << n->getDe
92 << "Universe radius = " << n->getU
93 return 0.0;
94 }
95
96 G4double pfl0 = 0.0;
97 const ParticleType t = p->getType();
98 const G4double kinE = p->getKineticEnergy(
99 if(kinE <= n->getPotential()->getFermiEner
100 pfl0 = n->getPotential()->getFermiMoment
101 } else {
102 const G4double tf0 = p->getPotentialEner
103 if(tf0<0.0) return 0.0;
104 pfl0 = std::sqrt(tf0*(tf0 + 2.0*mass));
105 }
106 const G4double pReflection = p->getReflect
107 const G4double reflectionRadius = n->getDe
108 const G4double pNominal = p->getMomentum()
109 const G4double nominalReflectionRadius = n
110 const G4double pl = pfl0*n->getDensity()->
111 vloc = std::sqrt(pl*pl + mass*mass) - mass
112
113 return vloc;
114 }
115
116 ThreeVector makeBoostVector(Particle const *
117 const G4double totalEnergy = p1->getEnergy
118 return ((p1->getMomentum() + p2->getMoment
119 }
120
121 G4double totalEnergyInCM(Particle const * co
122 return std::sqrt(squareTotalEnergyInCM(p1,
123 }
124
125 G4double squareTotalEnergyInCM(Particle cons
126 G4double beta2 = makeBoostVector(p1, p2).m
127 if(beta2 > 1.0) {
128 INCL_ERROR("squareTotalEnergyInCM: beta2
129 beta2 = 0.0;
130 }
131 return (1.0 - beta2)*std::pow(p1->getEnerg
132 }
133
134 G4double momentumInCM(Particle const * const
135 const G4double m1sq = std::pow(p1->getMass
136 const G4double m2sq = std::pow(p2->getMass
137 const G4double z = p1->getEnergy()*p2->get
138 G4double pcm2 = (z*z-m1sq*m2sq)/(2*z+m1sq+
139 if(pcm2 < 0.0) {
140 INCL_ERROR("momentumInCM: pcm2 == " << p
141 pcm2 = 0.0;
142 }
143 return std::sqrt(pcm2);
144 }
145
146 G4double momentumInCM(const G4double E, cons
147 return 0.5*std::sqrt((E*E - std::pow(M1 +
148 *(E*E - std::pow(M1 - M2, 2)))/E;
149 }
150
151 G4double momentumInLab(const G4double s, con
152 const G4double m1sq = m1*m1;
153 const G4double m2sq = m2*m2;
154 G4double plab2 = (s*s-2*s*(m1sq+m2sq)+(m1s
155 if(plab2 < 0.0) {
156 INCL_ERROR("momentumInLab: plab2 == " <<
157 plab2 = 0.0;
158 }
159 return std::sqrt(plab2);
160 }
161
162 G4double momentumInLab(Particle const * cons
163 const G4double m1 = p1->getMass();
164 const G4double m2 = p2->getMass();
165 const G4double s = squareTotalEnergyInCM(p
166 return momentumInLab(s, m1, m2);
167 }
168
169 G4double sumTotalEnergies(const ParticleList
170 G4double E = 0.0;
171 for(ParticleIter i=pl.begin(), e=pl.end();
172 E += (*i)->getEnergy();
173 }
174 return E;
175 }
176
177 ThreeVector sumMomenta(const ParticleList &p
178 ThreeVector p(0.0, 0.0, 0.0);
179 for(ParticleIter i=pl.begin(), e=pl.end();
180 p += (*i)->getMomentum();
181 }
182 return p;
183 }
184
185 G4double energy(const ThreeVector &p, const
186 return std::sqrt(p.mag2() + m*m);
187 }
188
189 G4double invariantMass(const G4double E, con
190 return std::sqrt(squareInvariantMass(E, p)
191 }
192
193 G4double squareInvariantMass(const G4double
194 return E*E - p.mag2();
195 }
196
197 G4double gammaFromKineticEnergy(const Partic
198 G4double mass;
199 if(p.theType==Composite)
200 mass = ParticleTable::getTableMass(p.the
201 else
202 mass = ParticleTable::getTableParticleMa
203 return (1.+EKin/mass);
204 }
205
206 }
207
208 }
209