| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and H
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35
36 #include "globals.hh"
37
38 /*
39 * G4INCLParticle.hh
40 *
41 * \date Jun 5, 2009
42 * \author Pekka Kaitaniemi
43 */
44
45 #ifndef PARTICLE_HH_
46 #define PARTICLE_HH_
47
48 #include "G4INCLThreeVector.hh"
49 #include "G4INCLParticleTable.hh"
50 #include "G4INCLParticleType.hh"
51 #include "G4INCLParticleSpecies.hh"
52 #include "G4INCLLogger.hh"
53 #include "G4INCLUnorderedVector.hh"
54 #include "G4INCLAllocationPool.hh"
55 #include <sstream>
56 #include <string>
57
58 namespace G4INCL {
59
60 class Particle;
61
62 class ParticleList : public UnorderedVector<
63 public:
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
75 class Particle {
76 public:
77 Particle();
78 Particle(ParticleType t, G4double energy,
79 Particle(ParticleType t, ThreeVector const
80 virtual ~Particle() {}
81
82 /** \brief Copy constructor
83 *
84 * Does not copy the particle ID.
85 */
86 Particle(const Particle &rhs) :
87 theZ(rhs.theZ),
88 theA(rhs.theA),
89 theS(rhs.theS),
90 theParticipantType(rhs.theParticipantTyp
91 theType(rhs.theType),
92 theEnergy(rhs.theEnergy),
93 theFrozenEnergy(rhs.theFrozenEnergy),
94 theMomentum(rhs.theMomentum),
95 theFrozenMomentum(rhs.theFrozenMomentum)
96 thePosition(rhs.thePosition),
97 nCollisions(rhs.nCollisions),
98 nDecays(rhs.nDecays),
99 thePotentialEnergy(rhs.thePotentialEnerg
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),
109 emissionTime(rhs.emissionTime),
110 outOfWell(rhs.outOfWell),
111 theMass(rhs.theMass)
112 {
113 if(rhs.thePropagationEnergy == &(rhs.t
114 thePropagationEnergy = &theFrozenEne
115 else
116 thePropagationEnergy = &theEnergy;
117 if(rhs.thePropagationMomentum == &(rhs
118 thePropagationMomentum = &theFrozenM
119 else
120 thePropagationMomentum = &theMomentu
121 // ID intentionally not copied
122 ID = nextID++;
123
124 theBiasCollisionVector = rhs.theBiasCo
125 }
126
127 protected:
128 /// \brief Helper method for the assignmen
129 void swap(Particle &rhs) {
130 std::swap(theZ, rhs.theZ);
131 std::swap(theA, rhs.theA);
132 std::swap(theS, rhs.theS);
133 std::swap(theParticipantType, rhs.thePar
134 std::swap(theType, rhs.theType);
135 if(rhs.thePropagationEnergy == &(rhs.the
136 thePropagationEnergy = &theFrozenEnerg
137 else
138 thePropagationEnergy = &theEnergy;
139 std::swap(theEnergy, rhs.theEnergy);
140 std::swap(theFrozenEnergy, rhs.theFrozen
141 if(rhs.thePropagationMomentum == &(rhs.t
142 thePropagationMomentum = &theFrozenMom
143 else
144 thePropagationMomentum = &theMomentum;
145 std::swap(theMomentum, rhs.theMomentum);
146 std::swap(theFrozenMomentum, rhs.theFroz
147 std::swap(thePosition, rhs.thePosition);
148 std::swap(nCollisions, rhs.nCollisions);
149 std::swap(nDecays, rhs.nDecays);
150 std::swap(thePotentialEnergy, rhs.thePot
151 // ID intentionally not swapped
152
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);
159 std::swap(emissionTime, rhs.emissionTime
160 std::swap(outOfWell, rhs.outOfWell);
161
162 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 }
170
171 public:
172
173 /** \brief Assignment operator
174 *
175 * Does not copy the particle ID.
176 */
177 Particle &operator=(const Particle &rhs) {
178 Particle temporaryParticle(rhs);
179 swap(temporaryParticle);
180 return *this;
181 }
182
183 /**
184 * Get the particle type.
185 * @see G4INCL::ParticleType
186 */
187 G4INCL::ParticleType getType() const {
188 return theType;
189 };
190
191 /// \brief Get the particle species
192 virtual G4INCL::ParticleSpecies getSpecies
193 return ParticleSpecies(theType);
194 };
195
196 void setType(ParticleType t) {
197 theType = t;
198 switch(theType)
199 {
200 case DeltaPlusPlus:
201 theA = 1;
202 theZ = 2;
203 theS = 0;
204 break;
205 case Proton:
206 case DeltaPlus:
207 theA = 1;
208 theZ = 1;
209 theS = 0;
210 break;
211 case Neutron:
212 case DeltaZero:
213 theA = 1;
214 theZ = 0;
215 theS = 0;
216 break;
217 case DeltaMinus:
218 theA = 1;
219 theZ = -1;
220 theS = 0;
221 break;
222 case PiPlus:
223 theA = 0;
224 theZ = 1;
225 theS = 0;
226 break;
227 case PiZero:
228 case Eta:
229 case Omega:
230 case EtaPrime:
231 case Photon:
232 theA = 0;
233 theZ = 0;
234 theS = 0;
235 break;
236 case PiMinus:
237 theA = 0;
238 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;
341 case Composite:
342 // INCL_ERROR("Trying to set particle
343 theA = 0;
344 theZ = 0;
345 theS = 0;
346 break;
347 case UnknownParticle:
348 theA = 0;
349 theZ = 0;
350 theS = 0;
351 INCL_ERROR("Trying to set particle t
352 break;
353 }
354
355 if( !isResonance() && t!=Composite )
356 setINCLMass();
357 }
358
359 /**
360 * Is this a nucleon?
361 */
362 G4bool isNucleon() const {
363 if(theType == G4INCL::Proton || theType
364 return true;
365 else
366 return false;
367 };
368
369 ParticipantType getParticipantType() const
370 return theParticipantType;
371 }
372
373 void setParticipantType(ParticipantType co
374 theParticipantType = p;
375 }
376
377 G4bool isParticipant() const {
378 return (theParticipantType==Participant)
379 }
380
381 G4bool isTargetSpectator() const {
382 return (theParticipantType==TargetSpecta
383 }
384
385 G4bool isProjectileSpectator() const {
386 return (theParticipantType==ProjectileSp
387 }
388
389 virtual void makeParticipant() {
390 theParticipantType = Participant;
391 }
392
393 virtual void makeTargetSpectator() {
394 theParticipantType = TargetSpectator;
395 }
396
397 virtual void makeProjectileSpectator() {
398 theParticipantType = ProjectileSpectator
399 }
400
401 /** \brief Is this a pion? */
402 G4bool isPion() const { return (theType ==
403
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? */
417 inline G4bool isResonance() const { return
418
419 /** \brief Is it a Delta? */
420 inline G4bool isDelta() const {
421 return (theType==DeltaPlusPlus || theTyp
422 theType==DeltaZero || theType==Delta
423
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
475 /** \brief Returns the baryon number. */
476 G4int getA() const { return theA; }
477
478 /** \brief Returns the charge number. */
479 G4int getZ() const { return theZ; }
480
481 /** \brief Returns the strangeness number.
482 G4int getS() const { return theS; }
483
484 G4double getBeta() const {
485 const G4double P = theMomentum.mag();
486 return P/theEnergy;
487 }
488
489 /**
490 * Returns a three vector we can give to t
491 *
492 * In order to go to the particle rest fra
493 * the boost vector by -1.0.
494 */
495 ThreeVector boostVector() const {
496 return theMomentum / theEnergy;
497 }
498
499 /**
500 * Boost the particle using a boost vector
501 *
502 * Example (go to the particle rest frame)
503 * particle->boost(particle->boostVector()
504 */
505 void boost(const ThreeVector &aBoostVector
506 const G4double beta2 = aBoostVector.mag2
507 const G4double gamma = 1.0 / std::sqrt(1
508 const G4double bp = theMomentum.dot(aBoo
509 const G4double alpha = (gamma*gamma)/(1.
510
511 theMomentum = theMomentum + aBoostVector
512 theEnergy = gamma * (theEnergy - bp);
513 }
514
515 /** \brief Lorentz-contract the particle p
516 *
517 * Apply Lorentz contraction to the positi
518 * direction of the boost vector.
519 *
520 * \param aBoostVector the boost vector (v
521 * \param refPos the reference position
522 */
523 void lorentzContract(const ThreeVector &aB
524 const G4double beta2 = aBoostVector.mag2
525 const G4double gamma = 1.0 / std::sqrt(1
526 const ThreeVector theRelativePosition =
527 const ThreeVector transversePosition = t
528 const ThreeVector longitudinalPosition =
529
530 thePosition = refPos + transversePositio
531 }
532
533 /** \brief Get the cached particle mass. *
534 inline G4double getMass() const { return t
535
536 /** \brief Get the INCL particle mass. */
537 inline G4double getINCLMass() const {
538 switch(theType) {
539 case Proton:
540 case Neutron:
541 case PiPlus:
542 case PiMinus:
543 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
569 break;
570
571 case DeltaPlusPlus:
572 case DeltaPlus:
573 case DeltaZero:
574 case DeltaMinus:
575 return theMass;
576 break;
577
578 case Composite:
579 return ParticleTable::getINCLMass(th
580 break;
581
582 default:
583 INCL_ERROR("Particle::getINCLMass: U
584 return 0.0;
585 break;
586 }
587 }
588
589 /** \brief Get the tabulated particle mass
590 inline virtual G4double getTableMass() con
591 switch(theType) {
592 case Proton:
593 case Neutron:
594 case PiPlus:
595 case PiMinus:
596 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
622 break;
623
624 case DeltaPlusPlus:
625 case DeltaPlus:
626 case DeltaZero:
627 case DeltaMinus:
628 return theMass;
629 break;
630
631 case Composite:
632 return ParticleTable::getTableMass(t
633 break;
634
635 default:
636 INCL_ERROR("Particle::getTableMass:
637 return 0.0;
638 break;
639 }
640 }
641
642 /** \brief Get the real particle mass. */
643 inline G4double getRealMass() const {
644 switch(theType) {
645 case Proton:
646 case Neutron:
647 case PiPlus:
648 case PiMinus:
649 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
675 break;
676
677 case DeltaPlusPlus:
678 case DeltaPlus:
679 case DeltaZero:
680 case DeltaMinus:
681 return theMass;
682 break;
683
684 case Composite:
685 return ParticleTable::getRealMass(th
686 break;
687
688 default:
689 INCL_ERROR("Particle::getRealMass: U
690 return 0.0;
691 break;
692 }
693 }
694
695 /// \brief Set the mass of the Particle to
696 void setRealMass() { setMass(getRealMass()
697
698 /// \brief Set the mass of the Particle to
699 void setTableMass() { setMass(getTableMass
700
701 /// \brief Set the mass of the Particle to
702 void setINCLMass() { setMass(getINCLMass()
703
704 /**\brief Computes correction on the emiss
705 *
706 * Computes the correction that must be ap
707 * order to obtain the correct Q-value for
708 * nucleus. For absorption, the correction
709 * the value returned by this function.
710 *
711 * \param AParent the mass number of the e
712 * \param ZParent the charge number of the
713 * \return the correction
714 */
715 G4double getEmissionQValueCorrection(const
716 const G4int SParent = 0;
717 const G4int ADaughter = AParent - theA;
718 const G4int ZDaughter = ZParent - theZ;
719 const G4int SDaughter = 0;
720
721 // Note the minus sign here
722 G4double theQValue;
723 if(isCluster())
724 theQValue = -ParticleTable::getTableQV
725 else {
726 const G4double massTableParent = Parti
727 const G4double massTableDaughter = Par
728 const G4double massTableParticle = get
729 theQValue = massTableParent - massTabl
730 }
731
732 const G4double massINCLParent = Particle
733 const G4double massINCLDaughter = Partic
734 const G4double massINCLParticle = getINC
735
736 // The rhs corresponds to the INCL Q-val
737 return theQValue - (massINCLParent-massI
738 }
739
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
768 *
769 * Computes the correction that must be ap
770 * order to obtain the correct Q-value for
771 * nucleus to another.
772 *
773 * Assumes that the receving nucleus is IN
774 * INCL separation energy.
775 *
776 * \param AFrom the mass number of the don
777 * \param ZFrom the charge number of the d
778 * \param ATo the mass number of the recei
779 * \param ZTo the charge number of the rec
780 * \return the correction
781 */
782 G4double getTransferQValueCorrection(const
783 const G4int SFrom = 0;
784 const G4int STo = 0;
785 const G4int AFromDaughter = AFrom - theA
786 const G4int ZFromDaughter = ZFrom - theZ
787 const G4int SFromDaughter = 0;
788 const G4int AToDaughter = ATo + theA;
789 const G4int ZToDaughter = ZTo + theZ;
790 const G4int SToDaughter = 0;
791 const G4double theQValue = ParticleTable
792
793 const G4double massINCLTo = ParticleTabl
794 const G4double massINCLToDaughter = Part
795 /* Note that here we have to use the tab
796 * cannot use theMass, because at this s
797 * still off-shell; and we cannot use ge
798 * violations of global energy conservat
799 */
800 const G4double massINCLParticle = getTab
801
802 // The rhs corresponds to the INCL Q-val
803 return theQValue - (massINCLToDaughter-m
804 }
805
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
884 *
885 * Uses the relativistic invariant
886 * \f[ m = \sqrt{E^2 - {\vec p}^2}\f]
887 **/
888 G4double getInvariantMass() const {
889 const G4double mass = std::pow(theEnergy
890 if(mass < 0.0) {
891 INCL_ERROR("E*E - p*p is negative." <<
892 return 0.0;
893 } else {
894 return std::sqrt(mass);
895 }
896 };
897
898 /// \brief Get the particle kinetic energy
899 inline G4double getKineticEnergy() const {
900
901 /// \brief Get the particle potential ener
902 inline G4double getPotentialEnergy() const
903
904 /// \brief Set the particle potential ener
905 inline void setPotentialEnergy(G4double v)
906
907 /**
908 * Get the energy of the particle in MeV.
909 */
910 G4double getEnergy() const
911 {
912 return theEnergy;
913 };
914
915 /**
916 * Set the mass of the particle in MeV/c^2
917 */
918 void setMass(G4double mass)
919 {
920 this->theMass = mass;
921 }
922
923 /**
924 * Set the energy of the particle in MeV.
925 */
926 void setEnergy(G4double energy)
927 {
928 this->theEnergy = energy;
929 };
930
931 /**
932 * Get the momentum vector.
933 */
934 const G4INCL::ThreeVector &getMomentum() c
935 {
936 return theMomentum;
937 };
938
939 /** Get the angular momentum w.r.t. the or
940 virtual G4INCL::ThreeVector getAngularMome
941 {
942 return thePosition.vector(theMomentum);
943 };
944
945 /**
946 * Set the momentum vector.
947 */
948 virtual void setMomentum(const G4INCL::Thr
949 {
950 this->theMomentum = momentum;
951 };
952
953 /**
954 * Set the position vector.
955 */
956 const G4INCL::ThreeVector &getPosition() c
957 {
958 return thePosition;
959 };
960
961 virtual void setPosition(const G4INCL::Thr
962 {
963 this->thePosition = position;
964 };
965
966 G4double getHelicity() { return theHelicit
967 void setHelicity(G4double h) { theHelicity
968
969 void propagate(G4double step) {
970 thePosition += ((*thePropagationMomentum
971 };
972
973 /** \brief Return the number of collisions
974 G4int getNumberOfCollisions() const { retu
975
976 /** \brief Set the number of collisions un
977 void setNumberOfCollisions(G4int n) { nCol
978
979 /** \brief Increment the number of collisi
980 void incrementNumberOfCollisions() { nColl
981
982 /** \brief Return the number of decays und
983 G4int getNumberOfDecays() const { return n
984
985 /** \brief Set the number of decays underg
986 void setNumberOfDecays(G4int n) { nDecays
987
988 /** \brief Increment the number of decays
989 void incrementNumberOfDecays() { nDecays++
990
991 /** \brief Mark the particle as out of its
992 *
993 * This flag is used to control pions crea
994 * in delta decay. The pion potential chec
995 * true (necessary in order to correctly e
996 * Nucleus::applyFinalState() method uses
997 * avatars should be generated for the par
998 */
999 void setOutOfWell() { outOfWell = true; }
1000
1001 /// \brief Check if the particle is out o
1002 G4bool isOutOfWell() const { return outOf
1003
1004 void setEmissionTime(G4double t) { emissi
1005 G4double getEmissionTime() { return emiss
1006
1007 /** \brief Transverse component of the po
1008 ThreeVector getTransversePosition() const
1009 return thePosition - getLongitudinalPos
1010 }
1011
1012 /** \brief Longitudinal component of the
1013 ThreeVector getLongitudinalPosition() con
1014 return *thePropagationMomentum * (thePo
1015 }
1016
1017 /** \brief Rescale the momentum to match
1018 const ThreeVector &adjustMomentumFromEner
1019
1020 /** \brief Recompute the energy to match
1021 G4double adjustEnergyFromMomentum();
1022
1023 G4bool isCluster() const {
1024 return (theType == Composite);
1025 }
1026
1027 /// \brief Set the frozen particle moment
1028 void setFrozenMomentum(const ThreeVector
1029
1030 /// \brief Set the frozen particle moment
1031 void setFrozenEnergy(const G4double energ
1032
1033 /// \brief Get the frozen particle moment
1034 ThreeVector getFrozenMomentum() const { r
1035
1036 /// \brief Get the frozen particle moment
1037 G4double getFrozenEnergy() const { return
1038
1039 /// \brief Get the propagation velocity o
1040 ThreeVector getPropagationVelocity() cons
1041
1042 /** \brief Freeze particle propagation
1043 *
1044 * Make the particle use theFrozenMomentu
1045 * propagation. The normal state can be r
1046 * thawPropagation() method.
1047 */
1048 void freezePropagation() {
1049 thePropagationMomentum = &theFrozenMome
1050 thePropagationEnergy = &theFrozenEnergy
1051 }
1052
1053 /** \brief Unfreeze particle propagation
1054 *
1055 * Make the particle use theMomentum and
1056 * this method to restore the normal prop
1057 * freezePropagation() method has been ca
1058 */
1059 void thawPropagation() {
1060 thePropagationMomentum = &theMomentum;
1061 thePropagationEnergy = &theEnergy;
1062 }
1063
1064 /** \brief Rotate the particle position a
1065 *
1066 * \param angle the rotation angle
1067 * \param axis a unit vector representing
1068 */
1069 virtual void rotatePositionAndMomentum(co
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);
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);
1090 theFrozenMomentum.rotate(angle, axis);
1091 }
1092
1093 std::string print() const {
1094 std::stringstream ss;
1095 ss << "Particle (ID = " << ID << ") typ
1096 ss << ParticleTable::getName(theType);
1097 ss << '\n'
1098 << " energy = " << theEnergy << '\n
1099 << " momentum = "
1100 << theMomentum.print()
1101 << '\n'
1102 << " position = "
1103 << thePosition.print()
1104 << '\n';
1105 return ss.str();
1106 };
1107
1108 std::string dump() const {
1109 std::stringstream ss;
1110 ss << "(particle " << ID << " ";
1111 ss << ParticleTable::getName(theType);
1112 ss << '\n'
1113 << thePosition.dump()
1114 << '\n'
1115 << theMomentum.dump()
1116 << '\n'
1117 << theEnergy << ")" << '\n';
1118 return ss.str();
1119 };
1120
1121 long getID() const { return ID; };
1122
1123 /**
1124 * Return a NULL pointer
1125 */
1126 ParticleList const *getParticles() const
1127 INCL_WARN("Particle::getParticles() met
1128 return 0;
1129 }
1130
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:
1227 G4int theZ, theA, theS;
1228 ParticipantType theParticipantType;
1229 G4INCL::ParticleType theType;
1230 G4double theEnergy;
1231 G4double *thePropagationEnergy;
1232 G4double theFrozenEnergy;
1233 G4INCL::ThreeVector theMomentum;
1234 G4INCL::ThreeVector *thePropagationMoment
1235 G4INCL::ThreeVector theFrozenMomentum;
1236 G4INCL::ThreeVector thePosition;
1237 G4int nCollisions;
1238 G4int nDecays;
1239 G4double thePotentialEnergy;
1240 long ID;
1241
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:
1255 G4double theHelicity;
1256 G4double emissionTime;
1257 G4bool outOfWell;
1258
1259 /// \brief Time ordered vector of all bia
1260 std::vector<G4int> theBiasCollisionVector
1261
1262 G4double theMass;
1263 static G4ThreadLocal long nextID;
1264
1265 INCL_DECLARE_ALLOCATION_POOL(Particle)
1266 };
1267 }
1268
1269 #endif /* PARTICLE_HH_ */
1270