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