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