| 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 /* 38 /*
39 * Particle.cc 39 * Particle.cc
40 * 40 *
41 * \date Jun 5, 2009 41 * \date Jun 5, 2009
42 * \author Pekka Kaitaniemi 42 * \author Pekka Kaitaniemi
43 */ 43 */
44 44
45 #include "G4INCLParticle.hh" 45 #include "G4INCLParticle.hh"
46 #include "G4INCLParticleTable.hh" 46 #include "G4INCLParticleTable.hh"
47 47
48 namespace G4INCL { 48 namespace G4INCL {
49 49
50 #ifdef INCLXX_IN_GEANT4_MODE << 50 G4ThreadLocal long Particle::nextID = 1;
51 std::vector<G4double> Particle::INCLBiasVe <<
52 #else <<
53 G4ThreadLocal std::vector<G4double> Partic <<
54 //G4VectorCache<G4double> Particle::INCLBias <<
55 #endif <<
56 G4ThreadLocal long Particle::nextID = 1; <<
57 G4ThreadLocal G4int Particle::nextBiasedColl <<
58 51
59 Particle::Particle() 52 Particle::Particle()
60 : theZ(0), theA(0), theS(0), << 53 : theZ(0), theA(0),
61 theParticipantType(TargetSpectator), 54 theParticipantType(TargetSpectator),
62 theType(UnknownParticle), 55 theType(UnknownParticle),
63 theEnergy(0.0), 56 theEnergy(0.0),
64 thePropagationEnergy(&theEnergy), 57 thePropagationEnergy(&theEnergy),
65 theFrozenEnergy(theEnergy), 58 theFrozenEnergy(theEnergy),
66 theMomentum(ThreeVector(0.,0.,0.)), 59 theMomentum(ThreeVector(0.,0.,0.)),
67 thePropagationMomentum(&theMomentum), 60 thePropagationMomentum(&theMomentum),
68 theFrozenMomentum(theMomentum), 61 theFrozenMomentum(theMomentum),
69 thePosition(ThreeVector(0.,0.,0.)), 62 thePosition(ThreeVector(0.,0.,0.)),
70 nCollisions(0), 63 nCollisions(0),
71 nDecays(0), 64 nDecays(0),
72 thePotentialEnergy(0.0), 65 thePotentialEnergy(0.0),
73 rpCorrelated(false), 66 rpCorrelated(false),
74 uncorrelatedMomentum(0.), 67 uncorrelatedMomentum(0.),
75 theParticleBias(1.), <<
76 theNKaon(0), <<
77 #ifdef INCLXX_IN_GEANT4_MODE <<
78 theParentResonancePDGCode(0), <<
79 theParentResonanceID(0), <<
80 #endif <<
81 theHelicity(0.0), 68 theHelicity(0.0),
82 emissionTime(0.0), 69 emissionTime(0.0),
83 outOfWell(false), 70 outOfWell(false),
84 theMass(0.) 71 theMass(0.)
85 { 72 {
86 ID = nextID; 73 ID = nextID;
87 nextID++; 74 nextID++;
88 } 75 }
89 76
90 Particle::Particle(ParticleType t, G4double 77 Particle::Particle(ParticleType t, G4double energy,
91 ThreeVector const &momentum, ThreeVector 78 ThreeVector const &momentum, ThreeVector const &position)
92 : theEnergy(energy), 79 : theEnergy(energy),
93 thePropagationEnergy(&theEnergy), 80 thePropagationEnergy(&theEnergy),
94 theFrozenEnergy(theEnergy), 81 theFrozenEnergy(theEnergy),
95 theMomentum(momentum), 82 theMomentum(momentum),
96 thePropagationMomentum(&theMomentum), 83 thePropagationMomentum(&theMomentum),
97 theFrozenMomentum(theMomentum), 84 theFrozenMomentum(theMomentum),
98 thePosition(position), 85 thePosition(position),
99 nCollisions(0), nDecays(0), 86 nCollisions(0), nDecays(0),
100 thePotentialEnergy(0.), << 87 thePotentialEnergy(0.),
101 rpCorrelated(false), << 88 rpCorrelated(false),
102 uncorrelatedMomentum(theMomentum.mag()), << 89 uncorrelatedMomentum(theMomentum.mag()),
103 theParticleBias(1.), << 90 theHelicity(0.0),
104 theNKaon(0), << 91 emissionTime(0.0), outOfWell(false)
105 #ifdef INCLXX_IN_GEANT4_MODE <<
106 theParentResonancePDGCode(0), <<
107 theParentResonanceID(0), <<
108 #endif <<
109 theHelicity(0.0), <<
110 emissionTime(0.0), outOfWell(false) <<
111 { 92 {
112 theParticipantType = TargetSpectator; 93 theParticipantType = TargetSpectator;
113 ID = nextID; 94 ID = nextID;
114 nextID++; 95 nextID++;
115 if(theEnergy <= 0.0) { 96 if(theEnergy <= 0.0) {
116 INCL_WARN("Particle with energy " << the 97 INCL_WARN("Particle with energy " << theEnergy << " created." << '\n');
117 } 98 }
118 setType(t); 99 setType(t);
119 setMass(getInvariantMass()); 100 setMass(getInvariantMass());
120 } 101 }
121 102
122 Particle::Particle(ParticleType t, 103 Particle::Particle(ParticleType t,
123 ThreeVector const &momentum, ThreeVector 104 ThreeVector const &momentum, ThreeVector const &position)
124 : thePropagationEnergy(&theEnergy), 105 : thePropagationEnergy(&theEnergy),
125 theMomentum(momentum), 106 theMomentum(momentum),
126 thePropagationMomentum(&theMomentum), 107 thePropagationMomentum(&theMomentum),
127 theFrozenMomentum(theMomentum), 108 theFrozenMomentum(theMomentum),
128 thePosition(position), 109 thePosition(position),
129 nCollisions(0), nDecays(0), 110 nCollisions(0), nDecays(0),
130 thePotentialEnergy(0.), << 111 thePotentialEnergy(0.),
131 rpCorrelated(false), << 112 rpCorrelated(false),
132 uncorrelatedMomentum(theMomentum.mag()), << 113 uncorrelatedMomentum(theMomentum.mag()),
133 theParticleBias(1.), << 114 theHelicity(0.0),
134 theNKaon(0), << 115 emissionTime(0.0), outOfWell(false)
135 #ifdef INCLXX_IN_GEANT4_MODE <<
136 theParentResonancePDGCode(0), <<
137 theParentResonanceID(0), <<
138 #endif <<
139 theHelicity(0.0), <<
140 emissionTime(0.0), outOfWell(false) <<
141 { 116 {
142 theParticipantType = TargetSpectator; 117 theParticipantType = TargetSpectator;
143 ID = nextID; 118 ID = nextID;
144 nextID++; 119 nextID++;
145 setType(t); 120 setType(t);
146 if( isResonance() ) { 121 if( isResonance() ) {
147 INCL_ERROR("Cannot create resonance with 122 INCL_ERROR("Cannot create resonance without specifying its momentum four-vector." << '\n');
148 } 123 }
149 G4double energy = std::sqrt(theMomentum.ma 124 G4double energy = std::sqrt(theMomentum.mag2() + theMass*theMass);
150 theEnergy = energy; 125 theEnergy = energy;
151 theFrozenEnergy = theEnergy; 126 theFrozenEnergy = theEnergy;
152 } 127 }
153 128
154 const ThreeVector &Particle::adjustMomentumF 129 const ThreeVector &Particle::adjustMomentumFromEnergy() {
155 const G4double p2 = theMomentum.mag2(); 130 const G4double p2 = theMomentum.mag2();
156 G4double newp2 = theEnergy*theEnergy - the 131 G4double newp2 = theEnergy*theEnergy - theMass*theMass;
157 if( newp2<0.0 ) { 132 if( newp2<0.0 ) {
158 INCL_ERROR("Particle has E^2 < m^2." << 133 INCL_ERROR("Particle has E^2 < m^2." << '\n' << print());
159 newp2 = 0.0; 134 newp2 = 0.0;
160 theEnergy = theMass; 135 theEnergy = theMass;
161 } 136 }
162 137
163 theMomentum *= std::sqrt(newp2/p2); 138 theMomentum *= std::sqrt(newp2/p2);
164 return theMomentum; 139 return theMomentum;
165 } 140 }
166 141
167 G4double Particle::adjustEnergyFromMomentum( 142 G4double Particle::adjustEnergyFromMomentum() {
168 theEnergy = std::sqrt(theMomentum.mag2() + 143 theEnergy = std::sqrt(theMomentum.mag2() + theMass*theMass);
169 return theEnergy; 144 return theEnergy;
170 } 145 }
171 146
172 void ParticleList::rotatePositionAndMomentum 147 void ParticleList::rotatePositionAndMomentum(const G4double angle, const ThreeVector &axis) const {
173 for(const_iterator i=begin(), e=end(); i!= 148 for(const_iterator i=begin(), e=end(); i!=e; ++i) {
174 (*i)->rotatePositionAndMomentum(angle, a 149 (*i)->rotatePositionAndMomentum(angle, axis);
175 } 150 }
176 } 151 }
177 152
178 void ParticleList::rotatePosition(const G4do 153 void ParticleList::rotatePosition(const G4double angle, const ThreeVector &axis) const {
179 for(const_iterator i=begin(), e=end(); i!= 154 for(const_iterator i=begin(), e=end(); i!=e; ++i) {
180 (*i)->rotatePosition(angle, axis); 155 (*i)->rotatePosition(angle, axis);
181 } 156 }
182 } 157 }
183 158
184 void ParticleList::rotateMomentum(const G4do 159 void ParticleList::rotateMomentum(const G4double angle, const ThreeVector &axis) const {
185 for(const_iterator i=begin(), e=end(); i!= 160 for(const_iterator i=begin(), e=end(); i!=e; ++i) {
186 (*i)->rotateMomentum(angle, axis); 161 (*i)->rotateMomentum(angle, axis);
187 } 162 }
188 } 163 }
189 164
190 void ParticleList::boost(const ThreeVector & 165 void ParticleList::boost(const ThreeVector &b) const {
191 for(const_iterator i=begin(), e=end(); i!= 166 for(const_iterator i=begin(), e=end(); i!=e; ++i) {
192 (*i)->boost(b); 167 (*i)->boost(b);
193 } 168 }
194 } <<
195 <<
196 G4double ParticleList::getParticleListBias() <<
197 if(G4int((*this).size())==0) return 1.; <<
198 std::vector<G4int> MergedVector; <<
199 for(ParticleIter i = (*this).begin(), e = <<
200 MergedVector = Particle::MergeVectorBi <<
201 } <<
202 return Particle::getBiasFromVector(std::mo <<
203 } <<
204 <<
205 std::vector<G4int> ParticleList::getParticle <<
206 std::vector<G4int> MergedVector; <<
207 if(G4int((*this).size())==0) return Merged <<
208 for(ParticleIter i = (*this).begin(), e = <<
209 MergedVector = Particle::MergeVectorBi <<
210 } <<
211 return MergedVector; <<
212 } <<
213 <<
214 void Particle::FillINCLBiasVector(G4double n <<
215 // assert(G4int(Particle::INCLBiasVector.size( <<
216 //assert(G4int(Particle::INCLBiasVector.Si <<
217 // assert(std::fabs(newBias - 1.) > 1E-6); <<
218 Particle::INCLBiasVector.push_back(newBias); <<
219 //Particle::INCLBiasVector.Push_back(newBias <<
220 Particle::nextBiasedCollisionID++; <<
221 } <<
222 <<
223 G4double Particle::getBiasFromVector(std::ve <<
224 if(VectorBias.empty()) return 1.; <<
225 <<
226 G4double ParticleBias = 1.; <<
227 <<
228 for(G4int i=0; i<G4int(VectorBias.size()); <<
229 ParticleBias *= Particle::INCLBiasVect <<
230 } <<
231 <<
232 return ParticleBias; <<
233 } <<
234 <<
235 std::vector<G4int> Particle::MergeVectorBias <<
236 std::vector<G4int> MergedVectorBias; <<
237 std::vector<G4int> VectorBias1 = p1->getBi <<
238 std::vector<G4int> VectorBias2 = p2->getBi <<
239 G4int i = 0; <<
240 G4int j = 0; <<
241 if(VectorBias1.size()==0 && VectorBias2.si <<
242 else if(VectorBias1.size()==0) return Vect <<
243 else if(VectorBias2.size()==0) return Vect <<
244 <<
245 while(i < G4int(VectorBias1.size()) || j < <<
246 if(VectorBias1[i]==VectorBias2[j]){ <<
247 MergedVectorBias.push_back(VectorB <<
248 i++; <<
249 j++; <<
250 if(i == G4int(VectorBias1.size())) <<
251 for(;j<G4int(VectorBias2.size( <<
252 } <<
253 else if(j == G4int(VectorBias2.siz <<
254 for(;i<G4int(VectorBias1.size( <<
255 } <<
256 } else if(VectorBias1[i]<VectorBias2[j <<
257 MergedVectorBias.push_back(VectorB <<
258 i++; <<
259 if(i == G4int(VectorBias1.size())) <<
260 for(;j<G4int(VectorBias2.size( <<
261 } <<
262 } <<
263 else { <<
264 MergedVectorBias.push_back(VectorB <<
265 j++; <<
266 if(j == G4int(VectorBias2.size())) <<
267 for(;i<G4int(VectorBias1.size( <<
268 } <<
269 } <<
270 } <<
271 return MergedVectorBias; <<
272 } <<
273 <<
274 std::vector<G4int> Particle::MergeVectorBias <<
275 std::vector<G4int> MergedVectorBias; <<
276 std::vector<G4int> VectorBias = p2->getBia <<
277 G4int i = 0; <<
278 G4int j = 0; <<
279 if(p1.size()==0 && VectorBias.size()==0) r <<
280 else if(p1.size()==0) return VectorBias; <<
281 else if(VectorBias.size()==0) return p1; <<
282 <<
283 while(i < G4int(p1.size()) || j < G4int(Ve <<
284 if(p1[i]==VectorBias[j]){ <<
285 MergedVectorBias.push_back(p1[i]); <<
286 i++; <<
287 j++; <<
288 if(i == G4int(p1.size())){ <<
289 for(;j<G4int(VectorBias.size() <<
290 } <<
291 else if(j == G4int(VectorBias.size <<
292 for(;i<G4int(p1.size());i++) M <<
293 } <<
294 } else if(p1[i]<VectorBias[j]){ <<
295 MergedVectorBias.push_back(p1[i]); <<
296 i++; <<
297 if(i == G4int(p1.size())){ <<
298 for(;j<G4int(VectorBias.size() <<
299 } <<
300 } <<
301 else { <<
302 MergedVectorBias.push_back(VectorB <<
303 j++; <<
304 if(j == G4int(VectorBias.size())){ <<
305 for(;i<G4int(p1.size());i++) M <<
306 } <<
307 } <<
308 } <<
309 return MergedVectorBias; <<
310 } <<
311 <<
312 G4double Particle::getTotalBias() { <<
313 G4double TotalBias = 1.; <<
314 for(G4int i=0; i<G4int(INCLBiasVector.si <<
315 return TotalBias; <<
316 } <<
317 <<
318 void Particle::setINCLBiasVector(std::vector <<
319 Particle::INCLBiasVector = std::move(New <<
320 } 169 }
321 } 170 }
322 171