| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
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 Helsinki Institute of Physics, Finland
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::INCLBiasVector;
52 #else
53 G4ThreadLocal std::vector<G4double> Particle::INCLBiasVector;
54 //G4VectorCache<G4double> Particle::INCLBiasVector;
55 #endif
56 G4ThreadLocal long Particle::nextID = 1;
57 G4ThreadLocal G4int Particle::nextBiasedCollisionID = 0;
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 energy,
91 ThreeVector const &momentum, ThreeVector const &position)
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 " << theEnergy << " created." << '\n');
117 }
118 setType(t);
119 setMass(getInvariantMass());
120 }
121
122 Particle::Particle(ParticleType t,
123 ThreeVector const &momentum, ThreeVector const &position)
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 without specifying its momentum four-vector." << '\n');
148 }
149 G4double energy = std::sqrt(theMomentum.mag2() + theMass*theMass);
150 theEnergy = energy;
151 theFrozenEnergy = theEnergy;
152 }
153
154 const ThreeVector &Particle::adjustMomentumFromEnergy() {
155 const G4double p2 = theMomentum.mag2();
156 G4double newp2 = theEnergy*theEnergy - theMass*theMass;
157 if( newp2<0.0 ) {
158 INCL_ERROR("Particle has E^2 < m^2." << '\n' << print());
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() + theMass*theMass);
169 return theEnergy;
170 }
171
172 void ParticleList::rotatePositionAndMomentum(const G4double angle, const ThreeVector &axis) const {
173 for(const_iterator i=begin(), e=end(); i!=e; ++i) {
174 (*i)->rotatePositionAndMomentum(angle, axis);
175 }
176 }
177
178 void ParticleList::rotatePosition(const G4double angle, const ThreeVector &axis) const {
179 for(const_iterator i=begin(), e=end(); i!=e; ++i) {
180 (*i)->rotatePosition(angle, axis);
181 }
182 }
183
184 void ParticleList::rotateMomentum(const G4double angle, const ThreeVector &axis) const {
185 for(const_iterator i=begin(), e=end(); i!=e; ++i) {
186 (*i)->rotateMomentum(angle, axis);
187 }
188 }
189
190 void ParticleList::boost(const ThreeVector &b) const {
191 for(const_iterator i=begin(), e=end(); i!=e; ++i) {
192 (*i)->boost(b);
193 }
194 }
195
196 G4double ParticleList::getParticleListBias() const {
197 if(G4int((*this).size())==0) return 1.;
198 std::vector<G4int> MergedVector;
199 for(ParticleIter i = (*this).begin(), e = (*this).end(); i!=e; ++i){
200 MergedVector = Particle::MergeVectorBias(MergedVector,(*i));
201 }
202 return Particle::getBiasFromVector(std::move(MergedVector));
203 }
204
205 std::vector<G4int> ParticleList::getParticleListBiasVector() const {
206 std::vector<G4int> MergedVector;
207 if(G4int((*this).size())==0) return MergedVector;
208 for(ParticleIter i = (*this).begin(), e = (*this).end(); i!=e; ++i){
209 MergedVector = Particle::MergeVectorBias(MergedVector,(*i));
210 }
211 return MergedVector;
212 }
213
214 void Particle::FillINCLBiasVector(G4double newBias){
215 // assert(G4int(Particle::INCLBiasVector.size())==nextBiasedCollisionID);
216 //assert(G4int(Particle::INCLBiasVector.Size())==nextBiasedCollisionID);
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::vector<G4int> VectorBias) {
224 if(VectorBias.empty()) return 1.;
225
226 G4double ParticleBias = 1.;
227
228 for(G4int i=0; i<G4int(VectorBias.size()); i++){
229 ParticleBias *= Particle::INCLBiasVector[G4int(VectorBias[i])];
230 }
231
232 return ParticleBias;
233 }
234
235 std::vector<G4int> Particle::MergeVectorBias(Particle const * const p1, Particle const * const p2){
236 std::vector<G4int> MergedVectorBias;
237 std::vector<G4int> VectorBias1 = p1->getBiasCollisionVector();
238 std::vector<G4int> VectorBias2 = p2->getBiasCollisionVector();
239 G4int i = 0;
240 G4int j = 0;
241 if(VectorBias1.size()==0 && VectorBias2.size()==0) return MergedVectorBias;
242 else if(VectorBias1.size()==0) return VectorBias2;
243 else if(VectorBias2.size()==0) return VectorBias1;
244
245 while(i < G4int(VectorBias1.size()) || j < G4int(VectorBias2.size())){
246 if(VectorBias1[i]==VectorBias2[j]){
247 MergedVectorBias.push_back(VectorBias1[i]);
248 i++;
249 j++;
250 if(i == G4int(VectorBias1.size())){
251 for(;j<G4int(VectorBias2.size());j++) MergedVectorBias.push_back(VectorBias2[j]);
252 }
253 else if(j == G4int(VectorBias2.size())){
254 for(;i<G4int(VectorBias1.size());i++) MergedVectorBias.push_back(VectorBias1[i]);
255 }
256 } else if(VectorBias1[i]<VectorBias2[j]){
257 MergedVectorBias.push_back(VectorBias1[i]);
258 i++;
259 if(i == G4int(VectorBias1.size())){
260 for(;j<G4int(VectorBias2.size());j++) MergedVectorBias.push_back(VectorBias2[j]);
261 }
262 }
263 else {
264 MergedVectorBias.push_back(VectorBias2[j]);
265 j++;
266 if(j == G4int(VectorBias2.size())){
267 for(;i<G4int(VectorBias1.size());i++) MergedVectorBias.push_back(VectorBias1[i]);
268 }
269 }
270 }
271 return MergedVectorBias;
272 }
273
274 std::vector<G4int> Particle::MergeVectorBias(std::vector<G4int> p1, Particle const * const p2){
275 std::vector<G4int> MergedVectorBias;
276 std::vector<G4int> VectorBias = p2->getBiasCollisionVector();
277 G4int i = 0;
278 G4int j = 0;
279 if(p1.size()==0 && VectorBias.size()==0) return MergedVectorBias;
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(VectorBias.size())){
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());j++) MergedVectorBias.push_back(VectorBias[j]);
290 }
291 else if(j == G4int(VectorBias.size())){
292 for(;i<G4int(p1.size());i++) MergedVectorBias.push_back(p1[i]);
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());j++) MergedVectorBias.push_back(VectorBias[j]);
299 }
300 }
301 else {
302 MergedVectorBias.push_back(VectorBias[j]);
303 j++;
304 if(j == G4int(VectorBias.size())){
305 for(;i<G4int(p1.size());i++) MergedVectorBias.push_back(p1[i]);
306 }
307 }
308 }
309 return MergedVectorBias;
310 }
311
312 G4double Particle::getTotalBias() {
313 G4double TotalBias = 1.;
314 for(G4int i=0; i<G4int(INCLBiasVector.size());i++) TotalBias *= Particle::INCLBiasVector[i];
315 return TotalBias;
316 }
317
318 void Particle::setINCLBiasVector(std::vector<G4double> NewVector) {
319 Particle::INCLBiasVector = std::move(NewVector);
320 }
321 }
322