| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35
36 #include "globals.hh"
37
38 /*
39 * 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<Particle*> {
63 public:
64 void rotatePositionAndMomentum(const G4double angle, const ThreeVector &axis) const;
65 void rotatePosition(const G4double angle, const ThreeVector &axis) const;
66 void rotateMomentum(const G4double angle, const ThreeVector &axis) const;
67 void boost(const ThreeVector &b) const;
68 G4double getParticleListBias() const;
69 std::vector<G4int> getParticleListBiasVector() const;
70 };
71
72 typedef ParticleList::const_iterator ParticleIter;
73 typedef ParticleList::iterator ParticleMutableIter;
74
75 class Particle {
76 public:
77 Particle();
78 Particle(ParticleType t, G4double energy, ThreeVector const &momentum, ThreeVector const &position);
79 Particle(ParticleType t, ThreeVector const &momentum, ThreeVector const &position);
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.theParticipantType),
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.thePotentialEnergy),
100 rpCorrelated(rhs.rpCorrelated),
101 uncorrelatedMomentum(rhs.uncorrelatedMomentum),
102 theParticleBias(rhs.theParticleBias),
103 theNKaon(rhs.theNKaon),
104 #ifdef INCLXX_IN_GEANT4_MODE
105 theParentResonancePDGCode(rhs.theParentResonancePDGCode),
106 theParentResonanceID(rhs.theParentResonanceID),
107 #endif
108 theHelicity(rhs.theHelicity),
109 emissionTime(rhs.emissionTime),
110 outOfWell(rhs.outOfWell),
111 theMass(rhs.theMass)
112 {
113 if(rhs.thePropagationEnergy == &(rhs.theFrozenEnergy))
114 thePropagationEnergy = &theFrozenEnergy;
115 else
116 thePropagationEnergy = &theEnergy;
117 if(rhs.thePropagationMomentum == &(rhs.theFrozenMomentum))
118 thePropagationMomentum = &theFrozenMomentum;
119 else
120 thePropagationMomentum = &theMomentum;
121 // ID intentionally not copied
122 ID = nextID++;
123
124 theBiasCollisionVector = rhs.theBiasCollisionVector;
125 }
126
127 protected:
128 /// \brief Helper method for the assignment operator
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.theParticipantType);
134 std::swap(theType, rhs.theType);
135 if(rhs.thePropagationEnergy == &(rhs.theFrozenEnergy))
136 thePropagationEnergy = &theFrozenEnergy;
137 else
138 thePropagationEnergy = &theEnergy;
139 std::swap(theEnergy, rhs.theEnergy);
140 std::swap(theFrozenEnergy, rhs.theFrozenEnergy);
141 if(rhs.thePropagationMomentum == &(rhs.theFrozenMomentum))
142 thePropagationMomentum = &theFrozenMomentum;
143 else
144 thePropagationMomentum = &theMomentum;
145 std::swap(theMomentum, rhs.theMomentum);
146 std::swap(theFrozenMomentum, rhs.theFrozenMomentum);
147 std::swap(thePosition, rhs.thePosition);
148 std::swap(nCollisions, rhs.nCollisions);
149 std::swap(nDecays, rhs.nDecays);
150 std::swap(thePotentialEnergy, rhs.thePotentialEnergy);
151 // ID intentionally not swapped
152
153 #ifdef INCLXX_IN_GEANT4_MODE
154 std::swap(theParentResonancePDGCode, rhs.theParentResonancePDGCode);
155 std::swap(theParentResonanceID, rhs.theParentResonanceID);
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.uncorrelatedMomentum);
165
166 std::swap(theParticleBias, rhs.theParticleBias);
167 std::swap(theBiasCollisionVector, rhs.theBiasCollisionVector);
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() const {
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 type to Composite! Construct a Cluster object instead" << '\n');
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 type to Unknown!" << '\n');
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 == G4INCL::Neutron)
364 return true;
365 else
366 return false;
367 };
368
369 ParticipantType getParticipantType() const {
370 return theParticipantType;
371 }
372
373 void setParticipantType(ParticipantType const p) {
374 theParticipantType = p;
375 }
376
377 G4bool isParticipant() const {
378 return (theParticipantType==Participant);
379 }
380
381 G4bool isTargetSpectator() const {
382 return (theParticipantType==TargetSpectator);
383 }
384
385 G4bool isProjectileSpectator() const {
386 return (theParticipantType==ProjectileSpectator);
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 == PiPlus || theType == PiZero || theType == PiMinus); }
403
404 /** \brief Is this an eta? */
405 G4bool isEta() const { return (theType == Eta); }
406
407 /** \brief Is this an omega? */
408 G4bool isOmega() const { return (theType == Omega); }
409
410 /** \brief Is this an etaprime? */
411 G4bool isEtaPrime() const { return (theType == EtaPrime); }
412
413 /** \brief Is this a photon? */
414 G4bool isPhoton() const { return (theType == Photon); }
415
416 /** \brief Is it a resonance? */
417 inline G4bool isResonance() const { return isDelta(); }
418
419 /** \brief Is it a Delta? */
420 inline G4bool isDelta() const {
421 return (theType==DeltaPlusPlus || theType==DeltaPlus ||
422 theType==DeltaZero || theType==DeltaMinus); }
423
424 /** \brief Is this a Sigma? */
425 G4bool isSigma() const { return (theType == SigmaPlus || theType == SigmaZero || theType == SigmaMinus); }
426
427 /** \brief Is this a Kaon? */
428 G4bool isKaon() const { return (theType == KPlus || theType == KZero); }
429
430 /** \brief Is this an antiKaon? */
431 G4bool isAntiKaon() const { return (theType == KZeroBar || theType == KMinus); }
432
433 /** \brief Is this a Lambda? */
434 G4bool isLambda() const { return (theType == Lambda); }
435
436 /** \brief Is this a Nucleon or a Lambda? */
437 G4bool isNucleonorLambda() const { return (isNucleon() || isLambda()); }
438
439 /** \brief Is this an Hyperon? */
440 G4bool isHyperon() const { return (isLambda() || isSigma() ); } //|| isXi()
441
442 /** \brief Is this a Meson? */
443 G4bool isMeson() const { return (isPion() || isKaon() || isAntiKaon() || isEta() || isEtaPrime() || isOmega()); }
444
445 /** \brief Is this a Baryon? */
446 G4bool isBaryon() const { return (isNucleon() || isResonance() || isHyperon()); }
447
448 /** \brief Is this a Strange? */
449 G4bool isStrange() const { return (isKaon() || isAntiKaon() || isHyperon()); }
450
451 /** \brief Is this a Xi? */
452 G4bool isXi() const { return (theType == XiZero || theType == XiMinus); }
453
454 /** \brief Is this an antinucleon? */
455 G4bool isAntiNucleon() const { return (theType == antiProton || theType == antiNeutron); }
456
457 /** \brief Is this an antiSigma? */
458 G4bool isAntiSigma() const { return (theType == antiSigmaPlus || theType == antiSigmaZero || theType == antiSigmaMinus); }
459
460 /** \brief Is this an antiXi? */
461 G4bool isAntiXi() const { return (theType == antiXiZero || theType == antiXiMinus); }
462
463 /** \brief Is this an antiLambda? */
464 G4bool isAntiLambda() const { return (theType == antiLambda); }
465
466 /** \brief Is this an antiHyperon? */
467 G4bool isAntiHyperon() const { return (isAntiLambda() || isAntiSigma() || isAntiXi()); }
468
469 /** \brief Is this an antiBaryon? */
470 G4bool isAntiBaryon() const { return (isAntiNucleon() || isAntiHyperon()); }
471
472 /** \brief Is this an antiNucleon or an antiLambda? */
473 G4bool isAntiNucleonorAntiLambda() const { return (isAntiNucleon() || isAntiLambda()); }
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 the boost() -method.
491 *
492 * In order to go to the particle rest frame you need to multiply
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.0 - beta2);
508 const G4double bp = theMomentum.dot(aBoostVector);
509 const G4double alpha = (gamma*gamma)/(1.0 + gamma);
510
511 theMomentum = theMomentum + aBoostVector * (alpha * bp - gamma * theEnergy);
512 theEnergy = gamma * (theEnergy - bp);
513 }
514
515 /** \brief Lorentz-contract the particle position around some center
516 *
517 * Apply Lorentz contraction to the position component along the
518 * direction of the boost vector.
519 *
520 * \param aBoostVector the boost vector (velocity) [c]
521 * \param refPos the reference position
522 */
523 void lorentzContract(const ThreeVector &aBoostVector, const ThreeVector &refPos) {
524 const G4double beta2 = aBoostVector.mag2();
525 const G4double gamma = 1.0 / std::sqrt(1.0 - beta2);
526 const ThreeVector theRelativePosition = thePosition - refPos;
527 const ThreeVector transversePosition = theRelativePosition - aBoostVector * (theRelativePosition.dot(aBoostVector) / aBoostVector.mag2());
528 const ThreeVector longitudinalPosition = theRelativePosition - transversePosition;
529
530 thePosition = refPos + transversePosition + longitudinalPosition / gamma;
531 }
532
533 /** \brief Get the cached particle mass. */
534 inline G4double getMass() const { return theMass; }
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(theType);
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(theA,theZ,theS);
580 break;
581
582 default:
583 INCL_ERROR("Particle::getINCLMass: Unknown particle type." << '\n');
584 return 0.0;
585 break;
586 }
587 }
588
589 /** \brief Get the tabulated particle mass. */
590 inline virtual G4double getTableMass() const {
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::getTableParticleMass(theType);
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(theA,theZ,theS);
633 break;
634
635 default:
636 INCL_ERROR("Particle::getTableMass: Unknown particle type." << '\n');
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(theType);
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(theA,theZ,theS);
686 break;
687
688 default:
689 INCL_ERROR("Particle::getRealMass: Unknown particle type." << '\n');
690 return 0.0;
691 break;
692 }
693 }
694
695 /// \brief Set the mass of the Particle to its real mass
696 void setRealMass() { setMass(getRealMass()); }
697
698 /// \brief Set the mass of the Particle to its table mass
699 void setTableMass() { setMass(getTableMass()); }
700
701 /// \brief Set the mass of the Particle to its table mass
702 void setINCLMass() { setMass(getINCLMass()); }
703
704 /**\brief Computes correction on the emission Q-value
705 *
706 * Computes the correction that must be applied to INCL particles in
707 * order to obtain the correct Q-value for particle emission from a given
708 * nucleus. For absorption, the correction is obviously equal to minus
709 * the value returned by this function.
710 *
711 * \param AParent the mass number of the emitting nucleus
712 * \param ZParent the charge number of the emitting nucleus
713 * \return the correction
714 */
715 G4double getEmissionQValueCorrection(const G4int AParent, const G4int ZParent) 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::getTableQValue(theA, theZ, theS, ADaughter, ZDaughter, SDaughter);
725 else {
726 const G4double massTableParent = ParticleTable::getTableMass(AParent,ZParent,SParent);
727 const G4double massTableDaughter = ParticleTable::getTableMass(ADaughter,ZDaughter,SDaughter);
728 const G4double massTableParticle = getTableMass();
729 theQValue = massTableParent - massTableDaughter - massTableParticle;
730 }
731
732 const G4double massINCLParent = ParticleTable::getINCLMass(AParent,ZParent,SParent);
733 const G4double massINCLDaughter = ParticleTable::getINCLMass(ADaughter,ZDaughter,SDaughter);
734 const G4double massINCLParticle = getINCLMass();
735
736 // The rhs corresponds to the INCL Q-value
737 return theQValue - (massINCLParent-massINCLDaughter-massINCLParticle);
738 }
739
740 G4double getEmissionPbarQvalueCorrection(const G4int AParent, const G4int ZParent, const G4bool Victim) const {
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 annihilated
747 ZDaughter = ZParent - 1;
748 }
749 else { //neutron is annihilated
750 ZDaughter = ZParent;
751 }
752
753 G4double theQValue; //same procedure as for normal case
754
755 const G4double massTableParent = ParticleTable::getTableMass(AParent,ZParent,SParent);
756 const G4double massTableDaughter = ParticleTable::getTableMass(ADaughter,ZDaughter,SDaughter);
757 const G4double massTableParticle = getTableMass();
758 theQValue = massTableParent - massTableDaughter - massTableParticle;
759
760 const G4double massINCLParent = ParticleTable::getINCLMass(AParent,ZParent,SParent);
761 const G4double massINCLDaughter = ParticleTable::getINCLMass(ADaughter,ZDaughter,SDaughter);
762 const G4double massINCLParticle = getINCLMass();
763
764 return theQValue - (massINCLParent-massINCLDaughter-massINCLParticle);
765 }
766
767 /**\brief Computes correction on the transfer Q-value
768 *
769 * Computes the correction that must be applied to INCL particles in
770 * order to obtain the correct Q-value for particle transfer from a given
771 * nucleus to another.
772 *
773 * Assumes that the receving nucleus is INCL's target nucleus, with the
774 * INCL separation energy.
775 *
776 * \param AFrom the mass number of the donating nucleus
777 * \param ZFrom the charge number of the donating nucleus
778 * \param ATo the mass number of the receiving nucleus
779 * \param ZTo the charge number of the receiving nucleus
780 * \return the correction
781 */
782 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;
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::getTableQValue(AToDaughter,ZToDaughter,SToDaughter,AFromDaughter,ZFromDaughter,SFromDaughter,AFrom,ZFrom,SFrom);
792
793 const G4double massINCLTo = ParticleTable::getINCLMass(ATo,ZTo,STo);
794 const G4double massINCLToDaughter = ParticleTable::getINCLMass(AToDaughter,ZToDaughter,SToDaughter);
795 /* Note that here we have to use the table mass in the INCL Q-value. We
796 * cannot use theMass, because at this stage the particle is probably
797 * still off-shell; and we cannot use getINCLMass(), because it leads to
798 * violations of global energy conservation.
799 */
800 const G4double massINCLParticle = getTableMass();
801
802 // The rhs corresponds to the INCL Q-value for particle absorption
803 return theQValue - (massINCLToDaughter-massINCLTo-massINCLParticle);
804 }
805
806 /**\brief Computes correction on the emission Q-value for hypernuclei
807 *
808 * Computes the correction that must be applied to INCL particles in
809 * order to obtain the correct Q-value for particle emission from a given
810 * nucleus. For absorption, the correction is obviously equal to minus
811 * the value returned by this function.
812 *
813 * \param AParent the mass number of the emitting nucleus
814 * \param ZParent the charge number of the emitting nucleus
815 * \param SParent the strangess number of the emitting nucleus
816 * \return the correction
817 */
818 G4double getEmissionQValueCorrection(const G4int AParent, const G4int ZParent, const G4int SParent) 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::getTableQValue(theA, theZ, theS, ADaughter, ZDaughter, SDaughter);
827 else {
828 const G4double massTableParent = ParticleTable::getTableMass(AParent,ZParent,SParent);
829 const G4double massTableDaughter = ParticleTable::getTableMass(ADaughter,ZDaughter,SDaughter);
830 const G4double massTableParticle = getTableMass();
831 theQValue = massTableParent - massTableDaughter - massTableParticle;
832 }
833
834 const G4double massINCLParent = ParticleTable::getINCLMass(AParent,ZParent,SParent);
835 const G4double massINCLDaughter = ParticleTable::getINCLMass(ADaughter,ZDaughter,SDaughter);
836 const G4double massINCLParticle = getINCLMass();
837
838 // The rhs corresponds to the INCL Q-value
839 return theQValue - (massINCLParent-massINCLDaughter-massINCLParticle);
840 }
841
842 /**\brief Computes correction on the transfer Q-value for hypernuclei
843 *
844 * Computes the correction that must be applied to INCL particles in
845 * order to obtain the correct Q-value for particle transfer from a given
846 * nucleus to another.
847 *
848 * Assumes that the receving nucleus is INCL's target nucleus, with the
849 * INCL separation energy.
850 *
851 * \param AFrom the mass number of the donating nucleus
852 * \param ZFrom the charge number of the donating nucleus
853 * \param SFrom the strangess number of the donating nucleus
854 * \param ATo the mass number of the receiving nucleus
855 * \param ZTo the charge number of the receiving nucleus
856 * \param STo the strangess number of the receiving nucleus
857 * \return the correction
858 */
859 G4double getTransferQValueCorrection(const G4int AFrom, const G4int ZFrom, const G4int SFrom, const G4int ATo, const G4int ZTo , const G4int STo) 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::getTableQValue(AToDaughter,ZToDaughter,SFromDaughter,AFromDaughter,ZFromDaughter,SToDaughter,AFrom,ZFrom,SFrom);
867
868 const G4double massINCLTo = ParticleTable::getINCLMass(ATo,ZTo,STo);
869 const G4double massINCLToDaughter = ParticleTable::getINCLMass(AToDaughter,ZToDaughter,SToDaughter);
870 /* Note that here we have to use the table mass in the INCL Q-value. We
871 * cannot use theMass, because at this stage the particle is probably
872 * still off-shell; and we cannot use getINCLMass(), because it leads to
873 * violations of global energy conservation.
874 */
875 const G4double massINCLParticle = getTableMass();
876
877 // The rhs corresponds to the INCL Q-value for particle absorption
878 return theQValue - (massINCLToDaughter-massINCLTo-massINCLParticle);
879 }
880
881
882
883 /** \brief Get the the particle invariant mass.
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, 2) - theMomentum.dot(theMomentum);
890 if(mass < 0.0) {
891 INCL_ERROR("E*E - p*p is negative." << '\n');
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 { return theEnergy - theMass; }
900
901 /// \brief Get the particle potential energy.
902 inline G4double getPotentialEnergy() const { return thePotentialEnergy; }
903
904 /// \brief Set the particle potential energy.
905 inline void setPotentialEnergy(G4double v) { thePotentialEnergy = 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() const
935 {
936 return theMomentum;
937 };
938
939 /** Get the angular momentum w.r.t. the origin */
940 virtual G4INCL::ThreeVector getAngularMomentum() const
941 {
942 return thePosition.vector(theMomentum);
943 };
944
945 /**
946 * Set the momentum vector.
947 */
948 virtual void setMomentum(const G4INCL::ThreeVector &momentum)
949 {
950 this->theMomentum = momentum;
951 };
952
953 /**
954 * Set the position vector.
955 */
956 const G4INCL::ThreeVector &getPosition() const
957 {
958 return thePosition;
959 };
960
961 virtual void setPosition(const G4INCL::ThreeVector &position)
962 {
963 this->thePosition = position;
964 };
965
966 G4double getHelicity() { return theHelicity; };
967 void setHelicity(G4double h) { theHelicity = h; };
968
969 void propagate(G4double step) {
970 thePosition += ((*thePropagationMomentum)*(step/(*thePropagationEnergy)));
971 };
972
973 /** \brief Return the number of collisions undergone by the particle. **/
974 G4int getNumberOfCollisions() const { return nCollisions; }
975
976 /** \brief Set the number of collisions undergone by the particle. **/
977 void setNumberOfCollisions(G4int n) { nCollisions = n; }
978
979 /** \brief Increment the number of collisions undergone by the particle. **/
980 void incrementNumberOfCollisions() { nCollisions++; }
981
982 /** \brief Return the number of decays undergone by the particle. **/
983 G4int getNumberOfDecays() const { return nDecays; }
984
985 /** \brief Set the number of decays undergone by the particle. **/
986 void setNumberOfDecays(G4int n) { nDecays = n; }
987
988 /** \brief Increment the number of decays undergone by the particle. **/
989 void incrementNumberOfDecays() { nDecays++; }
990
991 /** \brief Mark the particle as out of its potential well
992 *
993 * This flag is used to control pions created outside their potential well
994 * in delta decay. The pion potential checks it and returns zero if it is
995 * true (necessary in order to correctly enforce energy conservation). The
996 * Nucleus::applyFinalState() method uses it to determine whether new
997 * avatars should be generated for the particle.
998 */
999 void setOutOfWell() { outOfWell = true; }
1000
1001 /// \brief Check if the particle is out of its potential well
1002 G4bool isOutOfWell() const { return outOfWell; }
1003
1004 void setEmissionTime(G4double t) { emissionTime = t; }
1005 G4double getEmissionTime() { return emissionTime; };
1006
1007 /** \brief Transverse component of the position w.r.t. the momentum. */
1008 ThreeVector getTransversePosition() const {
1009 return thePosition - getLongitudinalPosition();
1010 }
1011
1012 /** \brief Longitudinal component of the position w.r.t. the momentum. */
1013 ThreeVector getLongitudinalPosition() const {
1014 return *thePropagationMomentum * (thePosition.dot(*thePropagationMomentum)/thePropagationMomentum->mag2());
1015 }
1016
1017 /** \brief Rescale the momentum to match the total energy. */
1018 const ThreeVector &adjustMomentumFromEnergy();
1019
1020 /** \brief Recompute the energy to match the momentum. */
1021 G4double adjustEnergyFromMomentum();
1022
1023 G4bool isCluster() const {
1024 return (theType == Composite);
1025 }
1026
1027 /// \brief Set the frozen particle momentum
1028 void setFrozenMomentum(const ThreeVector &momentum) { theFrozenMomentum = momentum; }
1029
1030 /// \brief Set the frozen particle momentum
1031 void setFrozenEnergy(const G4double energy) { theFrozenEnergy = energy; }
1032
1033 /// \brief Get the frozen particle momentum
1034 ThreeVector getFrozenMomentum() const { return theFrozenMomentum; }
1035
1036 /// \brief Get the frozen particle momentum
1037 G4double getFrozenEnergy() const { return theFrozenEnergy; }
1038
1039 /// \brief Get the propagation velocity of the particle
1040 ThreeVector getPropagationVelocity() const { return (*thePropagationMomentum)/(*thePropagationEnergy); }
1041
1042 /** \brief Freeze particle propagation
1043 *
1044 * Make the particle use theFrozenMomentum and theFrozenEnergy for
1045 * propagation. The normal state can be restored by calling the
1046 * thawPropagation() method.
1047 */
1048 void freezePropagation() {
1049 thePropagationMomentum = &theFrozenMomentum;
1050 thePropagationEnergy = &theFrozenEnergy;
1051 }
1052
1053 /** \brief Unfreeze particle propagation
1054 *
1055 * Make the particle use theMomentum and theEnergy for propagation. Call
1056 * this method to restore the normal propagation if the
1057 * freezePropagation() method has been called.
1058 */
1059 void thawPropagation() {
1060 thePropagationMomentum = &theMomentum;
1061 thePropagationEnergy = &theEnergy;
1062 }
1063
1064 /** \brief Rotate the particle position and momentum
1065 *
1066 * \param angle the rotation angle
1067 * \param axis a unit vector representing the rotation axis
1068 */
1069 virtual void rotatePositionAndMomentum(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 the rotation axis
1078 */
1079 virtual void rotatePosition(const G4double angle, const ThreeVector &axis) {
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 the rotation axis
1087 */
1088 virtual void rotateMomentum(const G4double angle, const ThreeVector &axis) {
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 << ") type = ";
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() method was called on a Particle object" << '\n');
1128 return 0;
1129 }
1130
1131 /** \brief Return the reflection momentum
1132 *
1133 * The reflection momentum is used by calls to getSurfaceRadius to compute
1134 * the radius of the sphere where the nucleon moves. It is necessary to
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 G4double p) { uncorrelatedMomentum = p; }
1146
1147 /// \brief Make the particle follow a strict r-p correlation
1148 void rpCorrelate() { rpCorrelated = true; }
1149
1150 /// \brief Make the particle not follow a strict r-p correlation
1151 void rpDecorrelate() { rpCorrelated = false; }
1152
1153 /// \brief Get the cosine of the angle between position and momentum
1154 G4double getCosRPAngle() const {
1155 const G4double norm = thePosition.mag2()*thePropagationMomentum->mag2();
1156 if(norm>0.)
1157 return thePosition.dot(*thePropagationMomentum) / std::sqrt(norm);
1158 else
1159 return 1.;
1160 }
1161
1162 /// \brief General bias vector function
1163 static G4double getTotalBias();
1164 static void setINCLBiasVector(std::vector<G4double> NewVector);
1165 static void FillINCLBiasVector(G4double newBias);
1166 static G4double getBiasFromVector(std::vector<G4int> VectorBias);
1167
1168 static std::vector<G4int> MergeVectorBias(Particle const * const p1, Particle const * const p2);
1169 static std::vector<G4int> MergeVectorBias(std::vector<G4int> p1, Particle const * const p2);
1170
1171 /// \brief Get the particle bias.
1172 G4double getParticleBias() const { return theParticleBias; };
1173
1174 /// \brief Set the particle bias.
1175 void setParticleBias(G4double ParticleBias) { this->theParticleBias = ParticleBias; }
1176
1177 /// \brief Get the vector list of biased vertices on the particle path.
1178 std::vector<G4int> getBiasCollisionVector() const { return theBiasCollisionVector; }
1179
1180 /// \brief Set the vector list of biased vertices on the particle path.
1181 void setBiasCollisionVector(std::vector<G4int> BiasCollisionVector) {
1182 this->theBiasCollisionVector = BiasCollisionVector;
1183 this->setParticleBias(Particle::getBiasFromVector(std::move(BiasCollisionVector)));
1184 }
1185
1186 /** \brief Number of Kaon inside de nucleus
1187 *
1188 * Put in the Particle class in order to calculate the
1189 * "correct" mass of composit particle.
1190 *
1191 */
1192
1193 G4int getNumberOfKaon() const { return theNKaon; };
1194 void setNumberOfKaon(const G4int NK) { theNKaon = NK; }
1195
1196 #ifdef INCLXX_IN_GEANT4_MODE
1197 G4int getParentResonancePDGCode() const { return theParentResonancePDGCode; };
1198 void setParentResonancePDGCode(const G4int parentPDGCode) { theParentResonancePDGCode = parentPDGCode; };
1199 G4int getParentResonanceID() const { return theParentResonanceID; };
1200 void setParentResonanceID(const G4int parentID) { theParentResonanceID = parentID; };
1201 #endif
1202
1203 public:
1204 /** \brief Time ordered vector of all bias applied
1205 *
1206 * /!\ Caution /!\
1207 * methods Assotiated to G4VectorCache<T> are:
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> INCLBiasVector;
1219 //static G4VectorCache<G4double> INCLBiasVector;
1220 #else
1221 static G4ThreadLocal std::vector<G4double> INCLBiasVector;
1222 //static G4VectorCache<G4double> INCLBiasVector;
1223 #endif
1224 static G4ThreadLocal G4int nextBiasedCollisionID;
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 *thePropagationMomentum;
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 nucleus (update during the cascade)
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 biased vertices on the particle path
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