| 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 /*
39 * Particle.cc
40 *
41 * \date Jun 5, 2009
42 * \author Pekka Kaitaniemi
43 */
44
45 #include "G4INCLParticle.hh"
46 #include "G4INCLParticleTable.hh"
47
48 namespace G4INCL {
49
50 #ifdef INCLXX_IN_GEANT4_MODE
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
59 Particle::Particle()
60 : theZ(0), theA(0), theS(0),
61 theParticipantType(TargetSpectator),
62 theType(UnknownParticle),
63 theEnergy(0.0),
64 thePropagationEnergy(&theEnergy),
65 theFrozenEnergy(theEnergy),
66 theMomentum(ThreeVector(0.,0.,0.)),
67 thePropagationMomentum(&theMomentum),
68 theFrozenMomentum(theMomentum),
69 thePosition(ThreeVector(0.,0.,0.)),
70 nCollisions(0),
71 nDecays(0),
72 thePotentialEnergy(0.0),
73 rpCorrelated(false),
74 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),
82 emissionTime(0.0),
83 outOfWell(false),
84 theMass(0.)
85 {
86 ID = nextID;
87 nextID++;
88 }
89
90 Particle::Particle(ParticleType t, G4double
91 ThreeVector const &momentum, ThreeVector
92 : theEnergy(energy),
93 thePropagationEnergy(&theEnergy),
94 theFrozenEnergy(theEnergy),
95 theMomentum(momentum),
96 thePropagationMomentum(&theMomentum),
97 theFrozenMomentum(theMomentum),
98 thePosition(position),
99 nCollisions(0), nDecays(0),
100 thePotentialEnergy(0.),
101 rpCorrelated(false),
102 uncorrelatedMomentum(theMomentum.mag()),
103 theParticleBias(1.),
104 theNKaon(0),
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 {
112 theParticipantType = TargetSpectator;
113 ID = nextID;
114 nextID++;
115 if(theEnergy <= 0.0) {
116 INCL_WARN("Particle with energy " << the
117 }
118 setType(t);
119 setMass(getInvariantMass());
120 }
121
122 Particle::Particle(ParticleType t,
123 ThreeVector const &momentum, ThreeVector
124 : thePropagationEnergy(&theEnergy),
125 theMomentum(momentum),
126 thePropagationMomentum(&theMomentum),
127 theFrozenMomentum(theMomentum),
128 thePosition(position),
129 nCollisions(0), nDecays(0),
130 thePotentialEnergy(0.),
131 rpCorrelated(false),
132 uncorrelatedMomentum(theMomentum.mag()),
133 theParticleBias(1.),
134 theNKaon(0),
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 {
142 theParticipantType = TargetSpectator;
143 ID = nextID;
144 nextID++;
145 setType(t);
146 if( isResonance() ) {
147 INCL_ERROR("Cannot create resonance with
148 }
149 G4double energy = std::sqrt(theMomentum.ma
150 theEnergy = energy;
151 theFrozenEnergy = theEnergy;
152 }
153
154 const ThreeVector &Particle::adjustMomentumF
155 const G4double p2 = theMomentum.mag2();
156 G4double newp2 = theEnergy*theEnergy - the
157 if( newp2<0.0 ) {
158 INCL_ERROR("Particle has E^2 < m^2." <<
159 newp2 = 0.0;
160 theEnergy = theMass;
161 }
162
163 theMomentum *= std::sqrt(newp2/p2);
164 return theMomentum;
165 }
166
167 G4double Particle::adjustEnergyFromMomentum(
168 theEnergy = std::sqrt(theMomentum.mag2() +
169 return theEnergy;
170 }
171
172 void ParticleList::rotatePositionAndMomentum
173 for(const_iterator i=begin(), e=end(); i!=
174 (*i)->rotatePositionAndMomentum(angle, a
175 }
176 }
177
178 void ParticleList::rotatePosition(const G4do
179 for(const_iterator i=begin(), e=end(); i!=
180 (*i)->rotatePosition(angle, axis);
181 }
182 }
183
184 void ParticleList::rotateMomentum(const G4do
185 for(const_iterator i=begin(), e=end(); i!=
186 (*i)->rotateMomentum(angle, axis);
187 }
188 }
189
190 void ParticleList::boost(const ThreeVector &
191 for(const_iterator i=begin(), e=end(); i!=
192 (*i)->boost(b);
193 }
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 }
321 }
322