| 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 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics
28 // Joseph Cugnon, University of Liege, Belgium << 28 // Davide Mancusi, CEA
29 // Jean-Christophe David, CEA-Saclay, France << 29 // Alain Boudard, CEA
30 // Pekka Kaitaniemi, CEA-Saclay, France, and H << 30 // Sylvie Leray, CEA
31 // Sylvie Leray, CEA-Saclay, France << 31 // Joseph Cugnon, University of Liege
32 // Davide Mancusi, CEA-Saclay, France << 32 //
>> 33 // INCL++ revision: v5.1.8
33 // 34 //
34 #define INCLXX_IN_GEANT4_MODE 1 35 #define INCLXX_IN_GEANT4_MODE 1
35 36
36 #include "globals.hh" 37 #include "globals.hh"
37 38
38 /* 39 /*
39 * G4INCLParticle.hh << 40 * Particle.hh
40 * 41 *
41 * \date Jun 5, 2009 42 * \date Jun 5, 2009
42 * \author Pekka Kaitaniemi 43 * \author Pekka Kaitaniemi
43 */ 44 */
44 45
45 #ifndef PARTICLE_HH_ 46 #ifndef PARTICLE_HH_
46 #define PARTICLE_HH_ 47 #define PARTICLE_HH_
47 48
48 #include "G4INCLThreeVector.hh" 49 #include "G4INCLThreeVector.hh"
49 #include "G4INCLParticleTable.hh" 50 #include "G4INCLParticleTable.hh"
50 #include "G4INCLParticleType.hh" 51 #include "G4INCLParticleType.hh"
51 #include "G4INCLParticleSpecies.hh" 52 #include "G4INCLParticleSpecies.hh"
52 #include "G4INCLLogger.hh" 53 #include "G4INCLLogger.hh"
53 #include "G4INCLUnorderedVector.hh" << 54 #include <list>
54 #include "G4INCLAllocationPool.hh" <<
55 #include <sstream> 55 #include <sstream>
56 #include <string> 56 #include <string>
>> 57 #include <algorithm>
57 58
58 namespace G4INCL { 59 namespace G4INCL {
59 60
60 class Particle; 61 class Particle;
61 62
62 class ParticleList : public UnorderedVector< << 63 typedef std::list<G4INCL::Particle*> ParticleList;
63 public: << 64 typedef std::list<G4INCL::Particle*>::const_iterator ParticleIter;
64 void rotatePositionAndMomentum(const G4d <<
65 void rotatePosition(const G4double angle <<
66 void rotateMomentum(const G4double angle <<
67 void boost(const ThreeVector &b) const; <<
68 G4double getParticleListBias() const; <<
69 std::vector<G4int> getParticleListBiasVe <<
70 }; <<
71 <<
72 typedef ParticleList::const_iterator Particl <<
73 typedef ParticleList::iterator Particl <<
74 65
75 class Particle { 66 class Particle {
76 public: 67 public:
77 Particle(); 68 Particle();
78 Particle(ParticleType t, G4double energy, << 69 Particle(ParticleType t, G4double energy, ThreeVector momentum, ThreeVector position);
79 Particle(ParticleType t, ThreeVector const << 70 Particle(ParticleType t, ThreeVector momentum, ThreeVector position);
80 virtual ~Particle() {} 71 virtual ~Particle() {}
81 72
82 /** \brief Copy constructor 73 /** \brief Copy constructor
83 * 74 *
84 * Does not copy the particle ID. 75 * Does not copy the particle ID.
85 */ 76 */
86 Particle(const Particle &rhs) : 77 Particle(const Particle &rhs) :
87 theZ(rhs.theZ), 78 theZ(rhs.theZ),
88 theA(rhs.theA), 79 theA(rhs.theA),
89 theS(rhs.theS), <<
90 theParticipantType(rhs.theParticipantTyp 80 theParticipantType(rhs.theParticipantType),
91 theType(rhs.theType), 81 theType(rhs.theType),
92 theEnergy(rhs.theEnergy), 82 theEnergy(rhs.theEnergy),
93 theFrozenEnergy(rhs.theFrozenEnergy), 83 theFrozenEnergy(rhs.theFrozenEnergy),
94 theMomentum(rhs.theMomentum), 84 theMomentum(rhs.theMomentum),
95 theFrozenMomentum(rhs.theFrozenMomentum) 85 theFrozenMomentum(rhs.theFrozenMomentum),
96 thePosition(rhs.thePosition), 86 thePosition(rhs.thePosition),
97 nCollisions(rhs.nCollisions), 87 nCollisions(rhs.nCollisions),
98 nDecays(rhs.nDecays), 88 nDecays(rhs.nDecays),
99 thePotentialEnergy(rhs.thePotentialEnerg 89 thePotentialEnergy(rhs.thePotentialEnergy),
100 rpCorrelated(rhs.rpCorrelated), <<
101 uncorrelatedMomentum(rhs.uncorrelatedMom <<
102 theParticleBias(rhs.theParticleBias), <<
103 theNKaon(rhs.theNKaon), <<
104 #ifdef INCLXX_IN_GEANT4_MODE <<
105 theParentResonancePDGCode(rhs.theParentR <<
106 theParentResonanceID(rhs.theParentResona <<
107 #endif <<
108 theHelicity(rhs.theHelicity), 90 theHelicity(rhs.theHelicity),
109 emissionTime(rhs.emissionTime), 91 emissionTime(rhs.emissionTime),
110 outOfWell(rhs.outOfWell), 92 outOfWell(rhs.outOfWell),
111 theMass(rhs.theMass) 93 theMass(rhs.theMass)
112 { 94 {
113 if(rhs.thePropagationEnergy == &(rhs.t 95 if(rhs.thePropagationEnergy == &(rhs.theFrozenEnergy))
114 thePropagationEnergy = &theFrozenEne 96 thePropagationEnergy = &theFrozenEnergy;
115 else 97 else
116 thePropagationEnergy = &theEnergy; 98 thePropagationEnergy = &theEnergy;
117 if(rhs.thePropagationMomentum == &(rhs 99 if(rhs.thePropagationMomentum == &(rhs.theFrozenMomentum))
118 thePropagationMomentum = &theFrozenM 100 thePropagationMomentum = &theFrozenMomentum;
119 else 101 else
120 thePropagationMomentum = &theMomentu 102 thePropagationMomentum = &theMomentum;
121 // ID intentionally not copied 103 // ID intentionally not copied
122 ID = nextID++; 104 ID = nextID++;
123 <<
124 theBiasCollisionVector = rhs.theBiasCo <<
125 } 105 }
126 106
127 protected: 107 protected:
128 /// \brief Helper method for the assignmen 108 /// \brief Helper method for the assignment operator
129 void swap(Particle &rhs) { 109 void swap(Particle &rhs) {
130 std::swap(theZ, rhs.theZ); 110 std::swap(theZ, rhs.theZ);
131 std::swap(theA, rhs.theA); 111 std::swap(theA, rhs.theA);
132 std::swap(theS, rhs.theS); <<
133 std::swap(theParticipantType, rhs.thePar 112 std::swap(theParticipantType, rhs.theParticipantType);
134 std::swap(theType, rhs.theType); 113 std::swap(theType, rhs.theType);
135 if(rhs.thePropagationEnergy == &(rhs.the 114 if(rhs.thePropagationEnergy == &(rhs.theFrozenEnergy))
136 thePropagationEnergy = &theFrozenEnerg 115 thePropagationEnergy = &theFrozenEnergy;
137 else 116 else
138 thePropagationEnergy = &theEnergy; 117 thePropagationEnergy = &theEnergy;
139 std::swap(theEnergy, rhs.theEnergy); 118 std::swap(theEnergy, rhs.theEnergy);
140 std::swap(theFrozenEnergy, rhs.theFrozen 119 std::swap(theFrozenEnergy, rhs.theFrozenEnergy);
141 if(rhs.thePropagationMomentum == &(rhs.t 120 if(rhs.thePropagationMomentum == &(rhs.theFrozenMomentum))
142 thePropagationMomentum = &theFrozenMom 121 thePropagationMomentum = &theFrozenMomentum;
143 else 122 else
144 thePropagationMomentum = &theMomentum; 123 thePropagationMomentum = &theMomentum;
145 std::swap(theMomentum, rhs.theMomentum); 124 std::swap(theMomentum, rhs.theMomentum);
146 std::swap(theFrozenMomentum, rhs.theFroz 125 std::swap(theFrozenMomentum, rhs.theFrozenMomentum);
147 std::swap(thePosition, rhs.thePosition); 126 std::swap(thePosition, rhs.thePosition);
148 std::swap(nCollisions, rhs.nCollisions); 127 std::swap(nCollisions, rhs.nCollisions);
149 std::swap(nDecays, rhs.nDecays); 128 std::swap(nDecays, rhs.nDecays);
150 std::swap(thePotentialEnergy, rhs.thePot 129 std::swap(thePotentialEnergy, rhs.thePotentialEnergy);
151 // ID intentionally not swapped 130 // ID intentionally not swapped
152 131
153 #ifdef INCLXX_IN_GEANT4_MODE <<
154 std::swap(theParentResonancePDGCode, rhs <<
155 std::swap(theParentResonanceID, rhs.theP <<
156 #endif <<
157 <<
158 std::swap(theHelicity, rhs.theHelicity); 132 std::swap(theHelicity, rhs.theHelicity);
159 std::swap(emissionTime, rhs.emissionTime 133 std::swap(emissionTime, rhs.emissionTime);
160 std::swap(outOfWell, rhs.outOfWell); 134 std::swap(outOfWell, rhs.outOfWell);
161 135
162 std::swap(theMass, rhs.theMass); 136 std::swap(theMass, rhs.theMass);
163 std::swap(rpCorrelated, rhs.rpCorrelated <<
164 std::swap(uncorrelatedMomentum, rhs.unco <<
165 <<
166 std::swap(theParticleBias, rhs.thePartic <<
167 std::swap(theBiasCollisionVector, rhs.th <<
168 <<
169 } 137 }
170 138
171 public: 139 public:
172 140
173 /** \brief Assignment operator 141 /** \brief Assignment operator
174 * 142 *
175 * Does not copy the particle ID. 143 * Does not copy the particle ID.
176 */ 144 */
177 Particle &operator=(const Particle &rhs) { 145 Particle &operator=(const Particle &rhs) {
178 Particle temporaryParticle(rhs); 146 Particle temporaryParticle(rhs);
179 swap(temporaryParticle); 147 swap(temporaryParticle);
180 return *this; 148 return *this;
181 } 149 }
182 150
183 /** 151 /**
184 * Get the particle type. 152 * Get the particle type.
185 * @see G4INCL::ParticleType 153 * @see G4INCL::ParticleType
186 */ 154 */
187 G4INCL::ParticleType getType() const { 155 G4INCL::ParticleType getType() const {
188 return theType; 156 return theType;
189 }; 157 };
190 158
191 /// \brief Get the particle species 159 /// \brief Get the particle species
192 virtual G4INCL::ParticleSpecies getSpecies 160 virtual G4INCL::ParticleSpecies getSpecies() const {
193 return ParticleSpecies(theType); 161 return ParticleSpecies(theType);
194 }; 162 };
195 163
196 void setType(ParticleType t) { 164 void setType(ParticleType t) {
197 theType = t; 165 theType = t;
198 switch(theType) 166 switch(theType)
199 { 167 {
200 case DeltaPlusPlus: 168 case DeltaPlusPlus:
201 theA = 1; 169 theA = 1;
202 theZ = 2; 170 theZ = 2;
203 theS = 0; <<
204 break; 171 break;
205 case Proton: 172 case Proton:
206 case DeltaPlus: 173 case DeltaPlus:
207 theA = 1; 174 theA = 1;
208 theZ = 1; 175 theZ = 1;
209 theS = 0; <<
210 break; 176 break;
211 case Neutron: 177 case Neutron:
212 case DeltaZero: 178 case DeltaZero:
213 theA = 1; 179 theA = 1;
214 theZ = 0; 180 theZ = 0;
215 theS = 0; <<
216 break; 181 break;
217 case DeltaMinus: 182 case DeltaMinus:
218 theA = 1; 183 theA = 1;
219 theZ = -1; 184 theZ = -1;
220 theS = 0; <<
221 break; 185 break;
222 case PiPlus: 186 case PiPlus:
223 theA = 0; 187 theA = 0;
224 theZ = 1; 188 theZ = 1;
225 theS = 0; <<
226 break; 189 break;
227 case PiZero: 190 case PiZero:
228 case Eta: <<
229 case Omega: <<
230 case EtaPrime: <<
231 case Photon: <<
232 theA = 0; 191 theA = 0;
233 theZ = 0; 192 theZ = 0;
234 theS = 0; <<
235 break; 193 break;
236 case PiMinus: 194 case PiMinus:
237 theA = 0; 195 theA = 0;
238 theZ = -1; 196 theZ = -1;
239 theS = 0; <<
240 break; <<
241 case Lambda: <<
242 theA = 1; <<
243 theZ = 0; <<
244 theS = -1; <<
245 break; <<
246 case SigmaPlus: <<
247 theA = 1; <<
248 theZ = 1; <<
249 theS = -1; <<
250 break; <<
251 case SigmaZero: <<
252 theA = 1; <<
253 theZ = 0; <<
254 theS = -1; <<
255 break; <<
256 case SigmaMinus: <<
257 theA = 1; <<
258 theZ = -1; <<
259 theS = -1; <<
260 break; <<
261 case antiProton: <<
262 theA = -1; <<
263 theZ = -1; <<
264 theS = 0; <<
265 break; <<
266 case XiMinus: <<
267 theA = 1; <<
268 theZ = -1; <<
269 theS = -2; <<
270 break; <<
271 case XiZero: <<
272 theA = 1; <<
273 theZ = 0; <<
274 theS = -2; <<
275 break; <<
276 case antiNeutron: <<
277 theA = -1; <<
278 theZ = 0; <<
279 theS = 0; <<
280 break; <<
281 case antiLambda: <<
282 theA = -1; <<
283 theZ = 0; <<
284 theS = 1; <<
285 break; <<
286 case antiSigmaMinus: <<
287 theA = -1; <<
288 theZ = 1; <<
289 theS = 1; <<
290 break; <<
291 case antiSigmaPlus: <<
292 theA = -1; <<
293 theZ = -1; <<
294 theS = 1; <<
295 break; <<
296 case antiSigmaZero: <<
297 theA = -1; <<
298 theZ = 0; <<
299 theS = 1; <<
300 break; <<
301 case antiXiMinus: <<
302 theA = -1; <<
303 theZ = 1; <<
304 theS = 2; <<
305 break; <<
306 case antiXiZero: <<
307 theA = -1; <<
308 theZ = 0; <<
309 theS = 2; <<
310 break; <<
311 case KPlus: <<
312 theA = 0; <<
313 theZ = 1; <<
314 theS = 1; <<
315 break; <<
316 case KZero: <<
317 theA = 0; <<
318 theZ = 0; <<
319 theS = 1; <<
320 break; <<
321 case KZeroBar: <<
322 theA = 0; <<
323 theZ = 0; <<
324 theS = -1; <<
325 break; <<
326 case KShort: <<
327 theA = 0; <<
328 theZ = 0; <<
329 // theS should not be defined <<
330 break; <<
331 case KLong: <<
332 theA = 0; <<
333 theZ = 0; <<
334 // theS should not be defined <<
335 break; <<
336 case KMinus: <<
337 theA = 0; <<
338 theZ = -1; <<
339 theS = -1; <<
340 break; 197 break;
341 case Composite: 198 case Composite:
342 // INCL_ERROR("Trying to set particle << 199 // ERROR("Trying to set particle type to Composite! Construct a Cluster object instead" << std::endl);
343 theA = 0; 200 theA = 0;
344 theZ = 0; 201 theZ = 0;
345 theS = 0; << 202 break;
346 break; <<
347 case UnknownParticle: 203 case UnknownParticle:
348 theA = 0; 204 theA = 0;
349 theZ = 0; 205 theZ = 0;
350 theS = 0; << 206 ERROR("Trying to set particle type to Unknown!" << std::endl);
351 INCL_ERROR("Trying to set particle t <<
352 break; 207 break;
353 } 208 }
354 209
355 if( !isResonance() && t!=Composite ) 210 if( !isResonance() && t!=Composite )
356 setINCLMass(); 211 setINCLMass();
357 } 212 }
358 213
359 /** 214 /**
360 * Is this a nucleon? 215 * Is this a nucleon?
361 */ 216 */
362 G4bool isNucleon() const { 217 G4bool isNucleon() const {
363 if(theType == G4INCL::Proton || theType 218 if(theType == G4INCL::Proton || theType == G4INCL::Neutron)
364 return true; << 219 return true;
365 else 220 else
366 return false; << 221 return false;
367 }; 222 };
368 223
369 ParticipantType getParticipantType() const 224 ParticipantType getParticipantType() const {
370 return theParticipantType; 225 return theParticipantType;
371 } 226 }
372 227
373 void setParticipantType(ParticipantType co 228 void setParticipantType(ParticipantType const p) {
374 theParticipantType = p; 229 theParticipantType = p;
375 } 230 }
376 231
377 G4bool isParticipant() const { 232 G4bool isParticipant() const {
378 return (theParticipantType==Participant) 233 return (theParticipantType==Participant);
379 } 234 }
380 235
381 G4bool isTargetSpectator() const { 236 G4bool isTargetSpectator() const {
382 return (theParticipantType==TargetSpecta 237 return (theParticipantType==TargetSpectator);
383 } 238 }
384 239
385 G4bool isProjectileSpectator() const { 240 G4bool isProjectileSpectator() const {
386 return (theParticipantType==ProjectileSp 241 return (theParticipantType==ProjectileSpectator);
387 } 242 }
388 243
389 virtual void makeParticipant() { 244 virtual void makeParticipant() {
390 theParticipantType = Participant; 245 theParticipantType = Participant;
391 } 246 }
392 247
393 virtual void makeTargetSpectator() { 248 virtual void makeTargetSpectator() {
394 theParticipantType = TargetSpectator; 249 theParticipantType = TargetSpectator;
395 } 250 }
396 251
397 virtual void makeProjectileSpectator() { 252 virtual void makeProjectileSpectator() {
398 theParticipantType = ProjectileSpectator 253 theParticipantType = ProjectileSpectator;
399 } 254 }
400 255
401 /** \brief Is this a pion? */ 256 /** \brief Is this a pion? */
402 G4bool isPion() const { return (theType == 257 G4bool isPion() const { return (theType == PiPlus || theType == PiZero || theType == PiMinus); }
403 258
404 /** \brief Is this an eta? */ <<
405 G4bool isEta() const { return (theType == <<
406 <<
407 /** \brief Is this an omega? */ <<
408 G4bool isOmega() const { return (theType = <<
409 <<
410 /** \brief Is this an etaprime? */ <<
411 G4bool isEtaPrime() const { return (theTyp <<
412 <<
413 /** \brief Is this a photon? */ <<
414 G4bool isPhoton() const { return (theType <<
415 <<
416 /** \brief Is it a resonance? */ 259 /** \brief Is it a resonance? */
417 inline G4bool isResonance() const { return 260 inline G4bool isResonance() const { return isDelta(); }
418 261
419 /** \brief Is it a Delta? */ 262 /** \brief Is it a Delta? */
420 inline G4bool isDelta() const { 263 inline G4bool isDelta() const {
421 return (theType==DeltaPlusPlus || theTyp 264 return (theType==DeltaPlusPlus || theType==DeltaPlus ||
422 theType==DeltaZero || theType==Delta << 265 theType==DeltaZero || theType==DeltaMinus);
423 << 266 }
424 /** \brief Is this a Sigma? */ <<
425 G4bool isSigma() const { return (theType = <<
426 <<
427 /** \brief Is this a Kaon? */ <<
428 G4bool isKaon() const { return (theType == <<
429 <<
430 /** \brief Is this an antiKaon? */ <<
431 G4bool isAntiKaon() const { return (theTyp <<
432 <<
433 /** \brief Is this a Lambda? */ <<
434 G4bool isLambda() const { return (theType <<
435 <<
436 /** \brief Is this a Nucleon or a Lambda? <<
437 G4bool isNucleonorLambda() const { return <<
438 <<
439 /** \brief Is this an Hyperon? */ <<
440 G4bool isHyperon() const { return (isLambd <<
441 <<
442 /** \brief Is this a Meson? */ <<
443 G4bool isMeson() const { return (isPion() <<
444 <<
445 /** \brief Is this a Baryon? */ <<
446 G4bool isBaryon() const { return (isNucleo <<
447 <<
448 /** \brief Is this a Strange? */ <<
449 G4bool isStrange() const { return (isKaon( <<
450 <<
451 /** \brief Is this a Xi? */ <<
452 G4bool isXi() const { return (theType == X <<
453 <<
454 /** \brief Is this an antinucleon? */ <<
455 G4bool isAntiNucleon() const { return (the <<
456 <<
457 /** \brief Is this an antiSigma? */ <<
458 G4bool isAntiSigma() const { return (theTy <<
459 <<
460 /** \brief Is this an antiXi? */ <<
461 G4bool isAntiXi() const { return (theType <<
462 <<
463 /** \brief Is this an antiLambda? */ <<
464 G4bool isAntiLambda() const { return (theT <<
465 <<
466 /** \brief Is this an antiHyperon? */ <<
467 G4bool isAntiHyperon() const { return (isA <<
468 <<
469 /** \brief Is this an antiBaryon? */ <<
470 G4bool isAntiBaryon() const { return (isAn <<
471 <<
472 /** \brief Is this an antiNucleon or an an <<
473 G4bool isAntiNucleonorAntiLambda() const { <<
474 267
475 /** \brief Returns the baryon number. */ 268 /** \brief Returns the baryon number. */
476 G4int getA() const { return theA; } 269 G4int getA() const { return theA; }
477 270
478 /** \brief Returns the charge number. */ 271 /** \brief Returns the charge number. */
479 G4int getZ() const { return theZ; } 272 G4int getZ() const { return theZ; }
480 <<
481 /** \brief Returns the strangeness number. <<
482 G4int getS() const { return theS; } <<
483 273
484 G4double getBeta() const { 274 G4double getBeta() const {
485 const G4double P = theMomentum.mag(); 275 const G4double P = theMomentum.mag();
486 return P/theEnergy; 276 return P/theEnergy;
487 } 277 }
488 278
489 /** 279 /**
490 * Returns a three vector we can give to t 280 * Returns a three vector we can give to the boost() -method.
491 * 281 *
492 * In order to go to the particle rest fra 282 * In order to go to the particle rest frame you need to multiply
493 * the boost vector by -1.0. 283 * the boost vector by -1.0.
494 */ 284 */
495 ThreeVector boostVector() const { 285 ThreeVector boostVector() const {
496 return theMomentum / theEnergy; 286 return theMomentum / theEnergy;
497 } 287 }
498 288
499 /** 289 /**
500 * Boost the particle using a boost vector 290 * Boost the particle using a boost vector.
501 * 291 *
502 * Example (go to the particle rest frame) 292 * Example (go to the particle rest frame):
503 * particle->boost(particle->boostVector() 293 * particle->boost(particle->boostVector());
504 */ 294 */
505 void boost(const ThreeVector &aBoostVector 295 void boost(const ThreeVector &aBoostVector) {
506 const G4double beta2 = aBoostVector.mag2 296 const G4double beta2 = aBoostVector.mag2();
507 const G4double gamma = 1.0 / std::sqrt(1 297 const G4double gamma = 1.0 / std::sqrt(1.0 - beta2);
508 const G4double bp = theMomentum.dot(aBoo 298 const G4double bp = theMomentum.dot(aBoostVector);
509 const G4double alpha = (gamma*gamma)/(1. 299 const G4double alpha = (gamma*gamma)/(1.0 + gamma);
510 300
511 theMomentum = theMomentum + aBoostVector 301 theMomentum = theMomentum + aBoostVector * (alpha * bp - gamma * theEnergy);
512 theEnergy = gamma * (theEnergy - bp); 302 theEnergy = gamma * (theEnergy - bp);
513 } 303 }
514 304
515 /** \brief Lorentz-contract the particle p 305 /** \brief Lorentz-contract the particle position around some center
516 * 306 *
517 * Apply Lorentz contraction to the positi 307 * Apply Lorentz contraction to the position component along the
518 * direction of the boost vector. 308 * direction of the boost vector.
519 * 309 *
520 * \param aBoostVector the boost vector (v 310 * \param aBoostVector the boost vector (velocity) [c]
521 * \param refPos the reference position 311 * \param refPos the reference position
522 */ 312 */
523 void lorentzContract(const ThreeVector &aB 313 void lorentzContract(const ThreeVector &aBoostVector, const ThreeVector &refPos) {
524 const G4double beta2 = aBoostVector.mag2 314 const G4double beta2 = aBoostVector.mag2();
525 const G4double gamma = 1.0 / std::sqrt(1 315 const G4double gamma = 1.0 / std::sqrt(1.0 - beta2);
526 const ThreeVector theRelativePosition = 316 const ThreeVector theRelativePosition = thePosition - refPos;
527 const ThreeVector transversePosition = t 317 const ThreeVector transversePosition = theRelativePosition - aBoostVector * (theRelativePosition.dot(aBoostVector) / aBoostVector.mag2());
528 const ThreeVector longitudinalPosition = 318 const ThreeVector longitudinalPosition = theRelativePosition - transversePosition;
529 319
530 thePosition = refPos + transversePositio 320 thePosition = refPos + transversePosition + longitudinalPosition / gamma;
531 } 321 }
532 322
533 /** \brief Get the cached particle mass. * 323 /** \brief Get the cached particle mass. */
534 inline G4double getMass() const { return t 324 inline G4double getMass() const { return theMass; }
535 325
536 /** \brief Get the INCL particle mass. */ 326 /** \brief Get the INCL particle mass. */
537 inline G4double getINCLMass() const { 327 inline G4double getINCLMass() const {
538 switch(theType) { 328 switch(theType) {
539 case Proton: 329 case Proton:
540 case Neutron: 330 case Neutron:
541 case PiPlus: 331 case PiPlus:
542 case PiMinus: 332 case PiMinus:
543 case PiZero: 333 case PiZero:
544 case Lambda: <<
545 case SigmaPlus: <<
546 case SigmaZero: <<
547 case SigmaMinus: <<
548 case antiProton: <<
549 case XiZero: <<
550 case XiMinus: <<
551 case antiNeutron: <<
552 case antiLambda: <<
553 case antiSigmaPlus: <<
554 case antiSigmaZero: <<
555 case antiSigmaMinus: <<
556 case antiXiZero: <<
557 case antiXiMinus: <<
558 case KPlus: <<
559 case KZero: <<
560 case KZeroBar: <<
561 case KShort: <<
562 case KLong: <<
563 case KMinus: <<
564 case Eta: <<
565 case Omega: <<
566 case EtaPrime: <<
567 case Photon: <<
568 return ParticleTable::getINCLMass(th 334 return ParticleTable::getINCLMass(theType);
569 break; 335 break;
570 336
571 case DeltaPlusPlus: 337 case DeltaPlusPlus:
572 case DeltaPlus: 338 case DeltaPlus:
573 case DeltaZero: 339 case DeltaZero:
574 case DeltaMinus: 340 case DeltaMinus:
575 return theMass; 341 return theMass;
576 break; 342 break;
577 343
578 case Composite: 344 case Composite:
579 return ParticleTable::getINCLMass(th << 345 return ParticleTable::getINCLMass(theA,theZ);
580 break; 346 break;
581 347
582 default: 348 default:
583 INCL_ERROR("Particle::getINCLMass: U << 349 ERROR("Particle::getINCLMass: Unknown particle type." << std::endl);
584 return 0.0; 350 return 0.0;
585 break; 351 break;
586 } 352 }
587 } 353 }
588 354
589 /** \brief Get the tabulated particle mass 355 /** \brief Get the tabulated particle mass. */
590 inline virtual G4double getTableMass() con 356 inline virtual G4double getTableMass() const {
591 switch(theType) { 357 switch(theType) {
592 case Proton: 358 case Proton:
593 case Neutron: 359 case Neutron:
594 case PiPlus: 360 case PiPlus:
595 case PiMinus: 361 case PiMinus:
596 case PiZero: 362 case PiZero:
597 case Lambda: <<
598 case SigmaPlus: <<
599 case SigmaZero: <<
600 case SigmaMinus: <<
601 case antiProton: <<
602 case XiZero: <<
603 case XiMinus: <<
604 case antiNeutron: <<
605 case antiLambda: <<
606 case antiSigmaPlus: <<
607 case antiSigmaZero: <<
608 case antiSigmaMinus: <<
609 case antiXiZero: <<
610 case antiXiMinus: <<
611 case KPlus: <<
612 case KZero: <<
613 case KZeroBar: <<
614 case KShort: <<
615 case KLong: <<
616 case KMinus: <<
617 case Eta: <<
618 case Omega: <<
619 case EtaPrime: <<
620 case Photon: <<
621 return ParticleTable::getTablePartic 363 return ParticleTable::getTableParticleMass(theType);
622 break; 364 break;
623 365
624 case DeltaPlusPlus: 366 case DeltaPlusPlus:
625 case DeltaPlus: 367 case DeltaPlus:
626 case DeltaZero: 368 case DeltaZero:
627 case DeltaMinus: 369 case DeltaMinus:
628 return theMass; 370 return theMass;
629 break; 371 break;
630 372
631 case Composite: 373 case Composite:
632 return ParticleTable::getTableMass(t << 374 return ParticleTable::getTableMass(theA,theZ);
633 break; 375 break;
634 376
635 default: 377 default:
636 INCL_ERROR("Particle::getTableMass: << 378 ERROR("Particle::getTableMass: Unknown particle type." << std::endl);
637 return 0.0; 379 return 0.0;
638 break; 380 break;
639 } 381 }
640 } 382 }
641 383
642 /** \brief Get the real particle mass. */ 384 /** \brief Get the real particle mass. */
643 inline G4double getRealMass() const { 385 inline G4double getRealMass() const {
644 switch(theType) { 386 switch(theType) {
645 case Proton: 387 case Proton:
646 case Neutron: 388 case Neutron:
647 case PiPlus: 389 case PiPlus:
648 case PiMinus: 390 case PiMinus:
649 case PiZero: 391 case PiZero:
650 case Lambda: <<
651 case SigmaPlus: <<
652 case SigmaZero: <<
653 case SigmaMinus: <<
654 case antiProton: <<
655 case XiZero: <<
656 case XiMinus: <<
657 case antiNeutron: <<
658 case antiLambda: <<
659 case antiSigmaPlus: <<
660 case antiSigmaZero: <<
661 case antiSigmaMinus: <<
662 case antiXiZero: <<
663 case antiXiMinus: <<
664 case KPlus: <<
665 case KZero: <<
666 case KZeroBar: <<
667 case KShort: <<
668 case KLong: <<
669 case KMinus: <<
670 case Eta: <<
671 case Omega: <<
672 case EtaPrime: <<
673 case Photon: <<
674 return ParticleTable::getRealMass(th 392 return ParticleTable::getRealMass(theType);
675 break; 393 break;
676 394
677 case DeltaPlusPlus: 395 case DeltaPlusPlus:
678 case DeltaPlus: 396 case DeltaPlus:
679 case DeltaZero: 397 case DeltaZero:
680 case DeltaMinus: 398 case DeltaMinus:
681 return theMass; 399 return theMass;
682 break; 400 break;
683 401
684 case Composite: 402 case Composite:
685 return ParticleTable::getRealMass(th << 403 return ParticleTable::getRealMass(theA,theZ);
686 break; 404 break;
687 405
688 default: 406 default:
689 INCL_ERROR("Particle::getRealMass: U << 407 ERROR("Particle::getRealMass: Unknown particle type." << std::endl);
690 return 0.0; 408 return 0.0;
691 break; 409 break;
692 } 410 }
693 } 411 }
694 412
695 /// \brief Set the mass of the Particle to 413 /// \brief Set the mass of the Particle to its real mass
696 void setRealMass() { setMass(getRealMass() 414 void setRealMass() { setMass(getRealMass()); }
697 415
698 /// \brief Set the mass of the Particle to 416 /// \brief Set the mass of the Particle to its table mass
699 void setTableMass() { setMass(getTableMass 417 void setTableMass() { setMass(getTableMass()); }
700 418
701 /// \brief Set the mass of the Particle to 419 /// \brief Set the mass of the Particle to its table mass
702 void setINCLMass() { setMass(getINCLMass() 420 void setINCLMass() { setMass(getINCLMass()); }
703 421
704 /**\brief Computes correction on the emiss 422 /**\brief Computes correction on the emission Q-value
705 * 423 *
706 * Computes the correction that must be ap 424 * Computes the correction that must be applied to INCL particles in
707 * order to obtain the correct Q-value for 425 * order to obtain the correct Q-value for particle emission from a given
708 * nucleus. For absorption, the correction 426 * nucleus. For absorption, the correction is obviously equal to minus
709 * the value returned by this function. 427 * the value returned by this function.
710 * 428 *
711 * \param AParent the mass number of the e 429 * \param AParent the mass number of the emitting nucleus
712 * \param ZParent the charge number of the 430 * \param ZParent the charge number of the emitting nucleus
713 * \return the correction 431 * \return the correction
714 */ 432 */
715 G4double getEmissionQValueCorrection(const 433 G4double getEmissionQValueCorrection(const G4int AParent, const G4int ZParent) const {
716 const G4int SParent = 0; <<
717 const G4int ADaughter = AParent - theA; 434 const G4int ADaughter = AParent - theA;
718 const G4int ZDaughter = ZParent - theZ; 435 const G4int ZDaughter = ZParent - theZ;
719 const G4int SDaughter = 0; <<
720 436
721 // Note the minus sign here 437 // Note the minus sign here
722 G4double theQValue; 438 G4double theQValue;
723 if(isCluster()) 439 if(isCluster())
724 theQValue = -ParticleTable::getTableQV << 440 theQValue = -ParticleTable::getTableQValue(theA, theZ, ADaughter, ZDaughter);
725 else { 441 else {
726 const G4double massTableParent = Parti << 442 const G4double massTableParent = ParticleTable::getTableMass(AParent,ZParent);
727 const G4double massTableDaughter = Par << 443 const G4double massTableDaughter = ParticleTable::getTableMass(ADaughter,ZDaughter);
728 const G4double massTableParticle = get 444 const G4double massTableParticle = getTableMass();
729 theQValue = massTableParent - massTabl 445 theQValue = massTableParent - massTableDaughter - massTableParticle;
730 } 446 }
731 447
732 const G4double massINCLParent = Particle << 448 const G4double massINCLParent = ParticleTable::getINCLMass(AParent,ZParent);
733 const G4double massINCLDaughter = Partic << 449 const G4double massINCLDaughter = ParticleTable::getINCLMass(ADaughter,ZDaughter);
734 const G4double massINCLParticle = getINC 450 const G4double massINCLParticle = getINCLMass();
735 451
736 // The rhs corresponds to the INCL Q-val 452 // The rhs corresponds to the INCL Q-value
737 return theQValue - (massINCLParent-massI 453 return theQValue - (massINCLParent-massINCLDaughter-massINCLParticle);
738 } 454 }
739 455
740 G4double getEmissionPbarQvalueCorrection(c <<
741 G4int SParent = 0; <<
742 G4int SDaughter = 0; <<
743 G4int ADaughter = AParent - 1; <<
744 G4int ZDaughter; <<
745 G4bool isProton = Victim; <<
746 if(isProton){ //proton is annihilate <<
747 ZDaughter = ZParent - 1; <<
748 } <<
749 else { //neutron is annihilated <<
750 ZDaughter = ZParent; <<
751 } <<
752 <<
753 G4double theQValue; //same procedure as <<
754 <<
755 const G4double massTableParent = Particl <<
756 const G4double massTableDaughter = Parti <<
757 const G4double massTableParticle = getTa <<
758 theQValue = massTableParent - massTableD <<
759 <<
760 const G4double massINCLParent = Particle <<
761 const G4double massINCLDaughter = Partic <<
762 const G4double massINCLParticle = getINC <<
763 <<
764 return theQValue - (massINCLParent-massI <<
765 } <<
766 <<
767 /**\brief Computes correction on the trans 456 /**\brief Computes correction on the transfer Q-value
768 * 457 *
769 * Computes the correction that must be ap 458 * Computes the correction that must be applied to INCL particles in
770 * order to obtain the correct Q-value for 459 * order to obtain the correct Q-value for particle transfer from a given
771 * nucleus to another. 460 * nucleus to another.
772 * 461 *
773 * Assumes that the receving nucleus is IN 462 * Assumes that the receving nucleus is INCL's target nucleus, with the
774 * INCL separation energy. 463 * INCL separation energy.
775 * 464 *
776 * \param AFrom the mass number of the don 465 * \param AFrom the mass number of the donating nucleus
777 * \param ZFrom the charge number of the d 466 * \param ZFrom the charge number of the donating nucleus
778 * \param ATo the mass number of the recei 467 * \param ATo the mass number of the receiving nucleus
779 * \param ZTo the charge number of the rec 468 * \param ZTo the charge number of the receiving nucleus
780 * \return the correction 469 * \return the correction
781 */ 470 */
782 G4double getTransferQValueCorrection(const 471 G4double getTransferQValueCorrection(const G4int AFrom, const G4int ZFrom, const G4int ATo, const G4int ZTo) const {
783 const G4int SFrom = 0; <<
784 const G4int STo = 0; <<
785 const G4int AFromDaughter = AFrom - theA 472 const G4int AFromDaughter = AFrom - theA;
786 const G4int ZFromDaughter = ZFrom - theZ 473 const G4int ZFromDaughter = ZFrom - theZ;
787 const G4int SFromDaughter = 0; <<
788 const G4int AToDaughter = ATo + theA; 474 const G4int AToDaughter = ATo + theA;
789 const G4int ZToDaughter = ZTo + theZ; 475 const G4int ZToDaughter = ZTo + theZ;
790 const G4int SToDaughter = 0; << 476 const G4double theQValue = ParticleTable::getTableQValue(AToDaughter,ZToDaughter,AFromDaughter,ZFromDaughter,AFrom,ZFrom);
791 const G4double theQValue = ParticleTable <<
792 477
793 const G4double massINCLTo = ParticleTabl << 478 const G4double massINCLTo = ParticleTable::getINCLMass(ATo,ZTo);
794 const G4double massINCLToDaughter = Part << 479 const G4double massINCLToDaughter = ParticleTable::getINCLMass(AToDaughter,ZToDaughter);
795 /* Note that here we have to use the tab 480 /* Note that here we have to use the table mass in the INCL Q-value. We
796 * cannot use theMass, because at this s 481 * cannot use theMass, because at this stage the particle is probably
797 * still off-shell; and we cannot use ge 482 * still off-shell; and we cannot use getINCLMass(), because it leads to
798 * violations of global energy conservat 483 * violations of global energy conservation.
799 */ 484 */
800 const G4double massINCLParticle = getTab 485 const G4double massINCLParticle = getTableMass();
801 486
802 // The rhs corresponds to the INCL Q-val 487 // The rhs corresponds to the INCL Q-value for particle absorption
803 return theQValue - (massINCLToDaughter-m 488 return theQValue - (massINCLToDaughter-massINCLTo-massINCLParticle);
804 } 489 }
805 490
806 /**\brief Computes correction on the emiss <<
807 * <<
808 * Computes the correction that must be ap <<
809 * order to obtain the correct Q-value for <<
810 * nucleus. For absorption, the correction <<
811 * the value returned by this function. <<
812 * <<
813 * \param AParent the mass number of the e <<
814 * \param ZParent the charge number of the <<
815 * \param SParent the strangess number of <<
816 * \return the correction <<
817 */ <<
818 G4double getEmissionQValueCorrection(const <<
819 const G4int ADaughter = AParent - theA; <<
820 const G4int ZDaughter = ZParent - theZ; <<
821 const G4int SDaughter = SParent - theS; <<
822 <<
823 // Note the minus sign here <<
824 G4double theQValue; <<
825 if(isCluster()) <<
826 theQValue = -ParticleTable::getTableQV <<
827 else { <<
828 const G4double massTableParent = Parti <<
829 const G4double massTableDaughter = Par <<
830 const G4double massTableParticle = get <<
831 theQValue = massTableParent - massTabl <<
832 } <<
833 <<
834 const G4double massINCLParent = Particle <<
835 const G4double massINCLDaughter = Partic <<
836 const G4double massINCLParticle = getINC <<
837 <<
838 // The rhs corresponds to the INCL Q-val <<
839 return theQValue - (massINCLParent-massI <<
840 } <<
841 <<
842 /**\brief Computes correction on the trans <<
843 * <<
844 * Computes the correction that must be ap <<
845 * order to obtain the correct Q-value for <<
846 * nucleus to another. <<
847 * <<
848 * Assumes that the receving nucleus is IN <<
849 * INCL separation energy. <<
850 * <<
851 * \param AFrom the mass number of the don <<
852 * \param ZFrom the charge number of the d <<
853 * \param SFrom the strangess number of th <<
854 * \param ATo the mass number of the recei <<
855 * \param ZTo the charge number of the rec <<
856 * \param STo the strangess number of the <<
857 * \return the correction <<
858 */ <<
859 G4double getTransferQValueCorrection(const <<
860 const G4int AFromDaughter = AFrom - theA <<
861 const G4int ZFromDaughter = ZFrom - theZ <<
862 const G4int SFromDaughter = SFrom - theS <<
863 const G4int AToDaughter = ATo + theA; <<
864 const G4int ZToDaughter = ZTo + theZ; <<
865 const G4int SToDaughter = STo + theS; <<
866 const G4double theQValue = ParticleTable <<
867 <<
868 const G4double massINCLTo = ParticleTabl <<
869 const G4double massINCLToDaughter = Part <<
870 /* Note that here we have to use the tab <<
871 * cannot use theMass, because at this s <<
872 * still off-shell; and we cannot use ge <<
873 * violations of global energy conservat <<
874 */ <<
875 const G4double massINCLParticle = getTab <<
876 <<
877 // The rhs corresponds to the INCL Q-val <<
878 return theQValue - (massINCLToDaughter-m <<
879 } <<
880 <<
881 <<
882 <<
883 /** \brief Get the the particle invariant 491 /** \brief Get the the particle invariant mass.
884 * 492 *
885 * Uses the relativistic invariant 493 * Uses the relativistic invariant
886 * \f[ m = \sqrt{E^2 - {\vec p}^2}\f] 494 * \f[ m = \sqrt{E^2 - {\vec p}^2}\f]
887 **/ 495 **/
888 G4double getInvariantMass() const { 496 G4double getInvariantMass() const {
889 const G4double mass = std::pow(theEnergy 497 const G4double mass = std::pow(theEnergy, 2) - theMomentum.dot(theMomentum);
890 if(mass < 0.0) { 498 if(mass < 0.0) {
891 INCL_ERROR("E*E - p*p is negative." << << 499 ERROR("E*E - p*p is negative." << std::endl);
892 return 0.0; 500 return 0.0;
893 } else { 501 } else {
894 return std::sqrt(mass); 502 return std::sqrt(mass);
895 } 503 }
896 }; 504 };
897 505
898 /// \brief Get the particle kinetic energy 506 /// \brief Get the particle kinetic energy.
899 inline G4double getKineticEnergy() const { 507 inline G4double getKineticEnergy() const { return theEnergy - theMass; }
900 508
901 /// \brief Get the particle potential ener 509 /// \brief Get the particle potential energy.
902 inline G4double getPotentialEnergy() const 510 inline G4double getPotentialEnergy() const { return thePotentialEnergy; }
903 511
904 /// \brief Set the particle potential ener 512 /// \brief Set the particle potential energy.
905 inline void setPotentialEnergy(G4double v) 513 inline void setPotentialEnergy(G4double v) { thePotentialEnergy = v; }
906 514
907 /** 515 /**
908 * Get the energy of the particle in MeV. 516 * Get the energy of the particle in MeV.
909 */ 517 */
910 G4double getEnergy() const 518 G4double getEnergy() const
911 { 519 {
912 return theEnergy; 520 return theEnergy;
913 }; 521 };
914 522
915 /** 523 /**
916 * Set the mass of the particle in MeV/c^2 524 * Set the mass of the particle in MeV/c^2.
917 */ 525 */
918 void setMass(G4double mass) 526 void setMass(G4double mass)
919 { 527 {
920 this->theMass = mass; 528 this->theMass = mass;
921 } 529 }
922 530
923 /** 531 /**
924 * Set the energy of the particle in MeV. 532 * Set the energy of the particle in MeV.
925 */ 533 */
926 void setEnergy(G4double energy) 534 void setEnergy(G4double energy)
927 { 535 {
928 this->theEnergy = energy; 536 this->theEnergy = energy;
929 }; 537 };
930 538
931 /** 539 /**
932 * Get the momentum vector. 540 * Get the momentum vector.
933 */ 541 */
934 const G4INCL::ThreeVector &getMomentum() c 542 const G4INCL::ThreeVector &getMomentum() const
935 { 543 {
936 return theMomentum; 544 return theMomentum;
937 }; 545 };
938 546
939 /** Get the angular momentum w.r.t. the or 547 /** Get the angular momentum w.r.t. the origin */
940 virtual G4INCL::ThreeVector getAngularMome 548 virtual G4INCL::ThreeVector getAngularMomentum() const
941 { 549 {
942 return thePosition.vector(theMomentum); 550 return thePosition.vector(theMomentum);
943 }; 551 };
944 552
945 /** 553 /**
946 * Set the momentum vector. 554 * Set the momentum vector.
947 */ 555 */
948 virtual void setMomentum(const G4INCL::Thr 556 virtual void setMomentum(const G4INCL::ThreeVector &momentum)
949 { 557 {
950 this->theMomentum = momentum; 558 this->theMomentum = momentum;
951 }; 559 };
952 560
953 /** 561 /**
954 * Set the position vector. 562 * Set the position vector.
955 */ 563 */
956 const G4INCL::ThreeVector &getPosition() c 564 const G4INCL::ThreeVector &getPosition() const
957 { 565 {
958 return thePosition; 566 return thePosition;
959 }; 567 };
960 568
961 virtual void setPosition(const G4INCL::Thr 569 virtual void setPosition(const G4INCL::ThreeVector &position)
962 { 570 {
963 this->thePosition = position; 571 this->thePosition = position;
964 }; 572 };
965 573
966 G4double getHelicity() { return theHelicit 574 G4double getHelicity() { return theHelicity; };
967 void setHelicity(G4double h) { theHelicity 575 void setHelicity(G4double h) { theHelicity = h; };
968 576
969 void propagate(G4double step) { 577 void propagate(G4double step) {
970 thePosition += ((*thePropagationMomentum 578 thePosition += ((*thePropagationMomentum)*(step/(*thePropagationEnergy)));
971 }; 579 };
972 580
973 /** \brief Return the number of collisions 581 /** \brief Return the number of collisions undergone by the particle. **/
974 G4int getNumberOfCollisions() const { retu 582 G4int getNumberOfCollisions() const { return nCollisions; }
975 583
976 /** \brief Set the number of collisions un 584 /** \brief Set the number of collisions undergone by the particle. **/
977 void setNumberOfCollisions(G4int n) { nCol 585 void setNumberOfCollisions(G4int n) { nCollisions = n; }
978 586
979 /** \brief Increment the number of collisi 587 /** \brief Increment the number of collisions undergone by the particle. **/
980 void incrementNumberOfCollisions() { nColl 588 void incrementNumberOfCollisions() { nCollisions++; }
981 589
982 /** \brief Return the number of decays und 590 /** \brief Return the number of decays undergone by the particle. **/
983 G4int getNumberOfDecays() const { return n 591 G4int getNumberOfDecays() const { return nDecays; }
984 592
985 /** \brief Set the number of decays underg 593 /** \brief Set the number of decays undergone by the particle. **/
986 void setNumberOfDecays(G4int n) { nDecays 594 void setNumberOfDecays(G4int n) { nDecays = n; }
987 595
988 /** \brief Increment the number of decays 596 /** \brief Increment the number of decays undergone by the particle. **/
989 void incrementNumberOfDecays() { nDecays++ 597 void incrementNumberOfDecays() { nDecays++; }
990 598
991 /** \brief Mark the particle as out of its 599 /** \brief Mark the particle as out of its potential well
992 * 600 *
993 * This flag is used to control pions crea 601 * This flag is used to control pions created outside their potential well
994 * in delta decay. The pion potential chec 602 * in delta decay. The pion potential checks it and returns zero if it is
995 * true (necessary in order to correctly e 603 * true (necessary in order to correctly enforce energy conservation). The
996 * Nucleus::applyFinalState() method uses 604 * Nucleus::applyFinalState() method uses it to determine whether new
997 * avatars should be generated for the par 605 * avatars should be generated for the particle.
998 */ 606 */
999 void setOutOfWell() { outOfWell = true; } 607 void setOutOfWell() { outOfWell = true; }
1000 608
1001 /// \brief Check if the particle is out o 609 /// \brief Check if the particle is out of its potential well
1002 G4bool isOutOfWell() const { return outOf 610 G4bool isOutOfWell() const { return outOfWell; }
1003 611
1004 void setEmissionTime(G4double t) { emissi 612 void setEmissionTime(G4double t) { emissionTime = t; }
1005 G4double getEmissionTime() { return emiss 613 G4double getEmissionTime() { return emissionTime; };
1006 614
1007 /** \brief Transverse component of the po 615 /** \brief Transverse component of the position w.r.t. the momentum. */
1008 ThreeVector getTransversePosition() const 616 ThreeVector getTransversePosition() const {
1009 return thePosition - getLongitudinalPos 617 return thePosition - getLongitudinalPosition();
1010 } 618 }
1011 619
1012 /** \brief Longitudinal component of the 620 /** \brief Longitudinal component of the position w.r.t. the momentum. */
1013 ThreeVector getLongitudinalPosition() con 621 ThreeVector getLongitudinalPosition() const {
1014 return *thePropagationMomentum * (thePo 622 return *thePropagationMomentum * (thePosition.dot(*thePropagationMomentum)/thePropagationMomentum->mag2());
1015 } 623 }
1016 624
1017 /** \brief Rescale the momentum to match 625 /** \brief Rescale the momentum to match the total energy. */
1018 const ThreeVector &adjustMomentumFromEner 626 const ThreeVector &adjustMomentumFromEnergy();
1019 627
1020 /** \brief Recompute the energy to match 628 /** \brief Recompute the energy to match the momentum. */
1021 G4double adjustEnergyFromMomentum(); 629 G4double adjustEnergyFromMomentum();
1022 630
>> 631 /** \brief Check if the particle belongs to a given list **/
>> 632 G4bool isInList(ParticleList const &l) const {
>> 633 return (std::find(l.begin(), l.end(), this)!=l.end());
>> 634 }
>> 635
1023 G4bool isCluster() const { 636 G4bool isCluster() const {
1024 return (theType == Composite); 637 return (theType == Composite);
1025 } 638 }
1026 639
1027 /// \brief Set the frozen particle moment 640 /// \brief Set the frozen particle momentum
1028 void setFrozenMomentum(const ThreeVector 641 void setFrozenMomentum(const ThreeVector &momentum) { theFrozenMomentum = momentum; }
1029 642
1030 /// \brief Set the frozen particle moment 643 /// \brief Set the frozen particle momentum
1031 void setFrozenEnergy(const G4double energ 644 void setFrozenEnergy(const G4double energy) { theFrozenEnergy = energy; }
1032 645
1033 /// \brief Get the frozen particle moment 646 /// \brief Get the frozen particle momentum
1034 ThreeVector getFrozenMomentum() const { r 647 ThreeVector getFrozenMomentum() const { return theFrozenMomentum; }
1035 648
1036 /// \brief Get the frozen particle moment 649 /// \brief Get the frozen particle momentum
1037 G4double getFrozenEnergy() const { return 650 G4double getFrozenEnergy() const { return theFrozenEnergy; }
1038 651
1039 /// \brief Get the propagation velocity o 652 /// \brief Get the propagation velocity of the particle
1040 ThreeVector getPropagationVelocity() cons 653 ThreeVector getPropagationVelocity() const { return (*thePropagationMomentum)/(*thePropagationEnergy); }
1041 654
1042 /** \brief Freeze particle propagation 655 /** \brief Freeze particle propagation
1043 * 656 *
1044 * Make the particle use theFrozenMomentu 657 * Make the particle use theFrozenMomentum and theFrozenEnergy for
1045 * propagation. The normal state can be r 658 * propagation. The normal state can be restored by calling the
1046 * thawPropagation() method. 659 * thawPropagation() method.
1047 */ 660 */
1048 void freezePropagation() { 661 void freezePropagation() {
1049 thePropagationMomentum = &theFrozenMome 662 thePropagationMomentum = &theFrozenMomentum;
1050 thePropagationEnergy = &theFrozenEnergy 663 thePropagationEnergy = &theFrozenEnergy;
1051 } 664 }
1052 665
1053 /** \brief Unfreeze particle propagation 666 /** \brief Unfreeze particle propagation
1054 * 667 *
1055 * Make the particle use theMomentum and 668 * Make the particle use theMomentum and theEnergy for propagation. Call
1056 * this method to restore the normal prop 669 * this method to restore the normal propagation if the
1057 * freezePropagation() method has been ca 670 * freezePropagation() method has been called.
1058 */ 671 */
1059 void thawPropagation() { 672 void thawPropagation() {
1060 thePropagationMomentum = &theMomentum; 673 thePropagationMomentum = &theMomentum;
1061 thePropagationEnergy = &theEnergy; 674 thePropagationEnergy = &theEnergy;
1062 } 675 }
1063 676
1064 /** \brief Rotate the particle position a 677 /** \brief Rotate the particle position and momentum
1065 * 678 *
1066 * \param angle the rotation angle 679 * \param angle the rotation angle
1067 * \param axis a unit vector representing 680 * \param axis a unit vector representing the rotation axis
1068 */ 681 */
1069 virtual void rotatePositionAndMomentum(co << 682 virtual void rotate(const G4double angle, const ThreeVector &axis) {
1070 rotatePosition(angle, axis); <<
1071 rotateMomentum(angle, axis); <<
1072 } <<
1073 <<
1074 /** \brief Rotate the particle position <<
1075 * <<
1076 * \param angle the rotation angle <<
1077 * \param axis a unit vector representing <<
1078 */ <<
1079 virtual void rotatePosition(const G4doubl <<
1080 thePosition.rotate(angle, axis); 683 thePosition.rotate(angle, axis);
1081 } <<
1082 <<
1083 /** \brief Rotate the particle momentum <<
1084 * <<
1085 * \param angle the rotation angle <<
1086 * \param axis a unit vector representing <<
1087 */ <<
1088 virtual void rotateMomentum(const G4doubl <<
1089 theMomentum.rotate(angle, axis); 684 theMomentum.rotate(angle, axis);
1090 theFrozenMomentum.rotate(angle, axis); 685 theFrozenMomentum.rotate(angle, axis);
1091 } 686 }
1092 687
1093 std::string print() const { 688 std::string print() const {
1094 std::stringstream ss; 689 std::stringstream ss;
1095 ss << "Particle (ID = " << ID << ") typ 690 ss << "Particle (ID = " << ID << ") type = ";
1096 ss << ParticleTable::getName(theType); 691 ss << ParticleTable::getName(theType);
1097 ss << '\n' << 692 ss << std::endl
1098 << " energy = " << theEnergy << '\n << 693 << " energy = " << theEnergy << std::endl
1099 << " momentum = " 694 << " momentum = "
1100 << theMomentum.print() 695 << theMomentum.print()
1101 << '\n' << 696 << std::endl
1102 << " position = " 697 << " position = "
1103 << thePosition.print() 698 << thePosition.print()
1104 << '\n'; << 699 << std::endl;
1105 return ss.str(); 700 return ss.str();
1106 }; 701 };
1107 702
1108 std::string dump() const { 703 std::string dump() const {
1109 std::stringstream ss; 704 std::stringstream ss;
1110 ss << "(particle " << ID << " "; 705 ss << "(particle " << ID << " ";
1111 ss << ParticleTable::getName(theType); 706 ss << ParticleTable::getName(theType);
1112 ss << '\n' << 707 ss << std::endl
1113 << thePosition.dump() 708 << thePosition.dump()
1114 << '\n' << 709 << std::endl
1115 << theMomentum.dump() 710 << theMomentum.dump()
1116 << '\n' << 711 << std::endl
1117 << theEnergy << ")" << '\n'; << 712 << theEnergy << ")" << std::endl;
1118 return ss.str(); 713 return ss.str();
1119 }; 714 };
1120 715
1121 long getID() const { return ID; }; 716 long getID() const { return ID; };
1122 717
1123 /** 718 /**
1124 * Return a NULL pointer 719 * Return a NULL pointer
1125 */ 720 */
1126 ParticleList const *getParticles() const 721 ParticleList const *getParticles() const {
1127 INCL_WARN("Particle::getParticles() met << 722 WARN("Particle::getParticles() method was called on a Particle object" << std::endl);
1128 return 0; 723 return 0;
1129 } 724 }
1130 725
1131 /** \brief Return the reflection momentum <<
1132 * <<
1133 * The reflection momentum is used by cal <<
1134 * the radius of the sphere where the nuc <<
1135 * introduce fuzzy r-p correlations. <<
1136 */ <<
1137 G4double getReflectionMomentum() const { <<
1138 if(rpCorrelated) <<
1139 return theMomentum.mag(); <<
1140 else <<
1141 return uncorrelatedMomentum; <<
1142 } <<
1143 <<
1144 /// \brief Set the uncorrelated momentum <<
1145 void setUncorrelatedMomentum(const G4doub <<
1146 <<
1147 /// \brief Make the particle follow a str <<
1148 void rpCorrelate() { rpCorrelated = true; <<
1149 <<
1150 /// \brief Make the particle not follow a <<
1151 void rpDecorrelate() { rpCorrelated = fal <<
1152 <<
1153 /// \brief Get the cosine of the angle be <<
1154 G4double getCosRPAngle() const { <<
1155 const G4double norm = thePosition.mag2( <<
1156 if(norm>0.) <<
1157 return thePosition.dot(*thePropagatio <<
1158 else <<
1159 return 1.; <<
1160 } <<
1161 <<
1162 /// \brief General bias vector function <<
1163 static G4double getTotalBias(); <<
1164 static void setINCLBiasVector(std::vector <<
1165 static void FillINCLBiasVector(G4double n <<
1166 static G4double getBiasFromVector(std::ve <<
1167 <<
1168 static std::vector<G4int> MergeVectorBias <<
1169 static std::vector<G4int> MergeVectorBias <<
1170 <<
1171 /// \brief Get the particle bias. <<
1172 G4double getParticleBias() const { return <<
1173 <<
1174 /// \brief Set the particle bias. <<
1175 void setParticleBias(G4double ParticleBia <<
1176 <<
1177 /// \brief Get the vector list of biased <<
1178 std::vector<G4int> getBiasCollisionVector <<
1179 <<
1180 /// \brief Set the vector list of biased <<
1181 void setBiasCollisionVector(std::vector<G <<
1182 this->theBiasCollisionVector = BiasCollis <<
1183 this->setParticleBias(Particle::getBiasFr <<
1184 } <<
1185 <<
1186 /** \brief Number of Kaon inside de nucle <<
1187 * <<
1188 * Put in the Particle class in order to <<
1189 * "correct" mass of composit particle. <<
1190 * <<
1191 */ <<
1192 <<
1193 G4int getNumberOfKaon() const { return th <<
1194 void setNumberOfKaon(const G4int NK) { th <<
1195 <<
1196 #ifdef INCLXX_IN_GEANT4_MODE <<
1197 G4int getParentResonancePDGCode() const { <<
1198 void setParentResonancePDGCode(const G4in <<
1199 G4int getParentResonanceID() const { retu <<
1200 void setParentResonanceID(const G4int par <<
1201 #endif <<
1202 <<
1203 public: <<
1204 /** \brief Time ordered vector of all bia <<
1205 * <<
1206 * /!\ Caution /!\ <<
1207 * methods Assotiated to G4VectorCache<T> <<
1208 * Push_back(…), <<
1209 * operator[], <<
1210 * Begin(), <<
1211 * End(), <<
1212 * Clear(), <<
1213 * Size() and <<
1214 * Pop_back() <<
1215 * <<
1216 */ <<
1217 #ifdef INCLXX_IN_GEANT4_MODE <<
1218 static std::vector<G4double> INCLBiasVe <<
1219 //static G4VectorCache<G4double> INCLBi <<
1220 #else <<
1221 static G4ThreadLocal std::vector<G4doub <<
1222 //static G4VectorCache<G4double> INCLBi <<
1223 #endif <<
1224 static G4ThreadLocal G4int nextBiasedColl <<
1225 <<
1226 protected: 726 protected:
1227 G4int theZ, theA, theS; << 727 G4int theZ, theA;
1228 ParticipantType theParticipantType; 728 ParticipantType theParticipantType;
1229 G4INCL::ParticleType theType; 729 G4INCL::ParticleType theType;
1230 G4double theEnergy; 730 G4double theEnergy;
1231 G4double *thePropagationEnergy; 731 G4double *thePropagationEnergy;
1232 G4double theFrozenEnergy; 732 G4double theFrozenEnergy;
1233 G4INCL::ThreeVector theMomentum; 733 G4INCL::ThreeVector theMomentum;
1234 G4INCL::ThreeVector *thePropagationMoment 734 G4INCL::ThreeVector *thePropagationMomentum;
1235 G4INCL::ThreeVector theFrozenMomentum; 735 G4INCL::ThreeVector theFrozenMomentum;
1236 G4INCL::ThreeVector thePosition; 736 G4INCL::ThreeVector thePosition;
1237 G4int nCollisions; 737 G4int nCollisions;
1238 G4int nDecays; 738 G4int nDecays;
1239 G4double thePotentialEnergy; 739 G4double thePotentialEnergy;
1240 long ID; 740 long ID;
1241 741
1242 G4bool rpCorrelated; <<
1243 G4double uncorrelatedMomentum; <<
1244 <<
1245 G4double theParticleBias; <<
1246 /// \brief The number of Kaons inside the <<
1247 G4int theNKaon; <<
1248 <<
1249 #ifdef INCLXX_IN_GEANT4_MODE <<
1250 G4int theParentResonancePDGCode; <<
1251 G4int theParentResonanceID; <<
1252 #endif <<
1253 <<
1254 private: 742 private:
1255 G4double theHelicity; 743 G4double theHelicity;
1256 G4double emissionTime; 744 G4double emissionTime;
1257 G4bool outOfWell; 745 G4bool outOfWell;
1258 <<
1259 /// \brief Time ordered vector of all bia <<
1260 std::vector<G4int> theBiasCollisionVector <<
1261 746
1262 G4double theMass; 747 G4double theMass;
1263 static G4ThreadLocal long nextID; << 748 static long nextID;
1264 749
1265 INCL_DECLARE_ALLOCATION_POOL(Particle) <<
1266 }; 750 };
1267 } 751 }
1268 752
1269 #endif /* PARTICLE_HH_ */ 753 #endif /* PARTICLE_HH_ */
1270 754