include/Geant4/G4INCLParticle.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // INCL++ intra-nuclear cascade model
0027 // Alain Boudard, CEA-Saclay, France
0028 // Joseph Cugnon, University of Liege, Belgium
0029 // Jean-Christophe David, CEA-Saclay, France
0030 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
0031 // Sylvie Leray, CEA-Saclay, France
0032 // Davide Mancusi, CEA-Saclay, France
0033 //
0034 #define INCLXX_IN_GEANT4_MODE 1
0035
0036 #include "globals.hh"
0037
0038 /*
0039  * G4INCLParticle.hh
0040  *
0041  *  \date Jun 5, 2009
0042  * \author Pekka Kaitaniemi
0043  */
0044
0045 #ifndef PARTICLE_HH_
0046 #define PARTICLE_HH_
0047
0048 #include "G4INCLThreeVector.hh"
0049 #include "G4INCLParticleTable.hh"
0050 #include "G4INCLParticleType.hh"
0051 #include "G4INCLParticleSpecies.hh"
0052 #include "G4INCLLogger.hh"
0053 #include "G4INCLUnorderedVector.hh"
0054 #include "G4INCLAllocationPool.hh"
0055 #include <sstream>
0056 #include <string>
0057
0058 namespace G4INCL {
0059
0060   class Particle;
0061
0062   class ParticleList : public UnorderedVector<Particle*> {
0063     public:
0064       void rotatePositionAndMomentum(const G4double angle, const ThreeVector &axis) const;
0065       void rotatePosition(const G4double angle, const ThreeVector &axis) const;
0066       void rotateMomentum(const G4double angle, const ThreeVector &axis) const;
0067       void boost(const ThreeVector &b) const;
0068       G4double getParticleListBias() const;
0069       std::vector<G4int> getParticleListBiasVector() const;
0070   };
0071
0072   typedef ParticleList::const_iterator ParticleIter;
0073   typedef ParticleList::iterator       ParticleMutableIter;
0074
0075   class Particle {
0076   public:
0077     Particle();
0078     Particle(ParticleType t, G4double energy, ThreeVector const &momentum, ThreeVector const &position);
0079     Particle(ParticleType t, ThreeVector const &momentum, ThreeVector const &position);
0080     virtual ~Particle() {}
0081
0082     /** \brief Copy constructor
0083      *
0084      * Does not copy the particle ID.
0085      */
0086     Particle(const Particle &rhs) :
0087       theZ(rhs.theZ),
0088       theA(rhs.theA),
0089       theS(rhs.theS),
0090       theParticipantType(rhs.theParticipantType),
0091       theType(rhs.theType),
0092       theEnergy(rhs.theEnergy),
0093       theFrozenEnergy(rhs.theFrozenEnergy),
0094       theMomentum(rhs.theMomentum),
0095       theFrozenMomentum(rhs.theFrozenMomentum),
0096       thePosition(rhs.thePosition),
0097       nCollisions(rhs.nCollisions),
0098       nDecays(rhs.nDecays),
0099       nSrcPair(rhs.nSrcPair),
0100       thePotentialEnergy(rhs.thePotentialEnergy),
0101       rpCorrelated(rhs.rpCorrelated),
0102       uncorrelatedMomentum(rhs.uncorrelatedMomentum),
0103       theParticleBias(rhs.theParticleBias),
0104       theNKaon(rhs.theNKaon),
0105 #ifdef INCLXX_IN_GEANT4_MODE
0106       theParentResonancePDGCode(rhs.theParentResonancePDGCode),
0107       theParentResonanceID(rhs.theParentResonanceID),
0108 #endif
0109       theHelicity(rhs.theHelicity),
0110       emissionTime(rhs.emissionTime),
0111       outOfWell(rhs.outOfWell),
0112       theSrcPartner(rhs.theSrcPartner),
0113       theMass(rhs.theMass)
0114       {
0115         if(rhs.thePropagationEnergy == &(rhs.theFrozenEnergy))
0116           thePropagationEnergy = &theFrozenEnergy;
0117         else
0118           thePropagationEnergy = &theEnergy;
0119         if(rhs.thePropagationMomentum == &(rhs.theFrozenMomentum))
0120           thePropagationMomentum = &theFrozenMomentum;
0121         else
0122           thePropagationMomentum = &theMomentum;
0123         // ID intentionally not copied
0124         ID = nextID++;
0125
0126         theBiasCollisionVector = rhs.theBiasCollisionVector;
0127       }
0128
0129   protected:
0130     /// \brief Helper method for the assignment operator
0131     void swap(Particle &rhs) {
0132       std::swap(theZ, rhs.theZ);
0133       std::swap(theA, rhs.theA);
0134       std::swap(theS, rhs.theS);
0135       std::swap(theParticipantType, rhs.theParticipantType);
0136       std::swap(theType, rhs.theType);
0137       if(rhs.thePropagationEnergy == &(rhs.theFrozenEnergy))
0138         thePropagationEnergy = &theFrozenEnergy;
0139       else
0140         thePropagationEnergy = &theEnergy;
0141       std::swap(theEnergy, rhs.theEnergy);
0142       std::swap(theFrozenEnergy, rhs.theFrozenEnergy);
0143       if(rhs.thePropagationMomentum == &(rhs.theFrozenMomentum))
0144         thePropagationMomentum = &theFrozenMomentum;
0145       else
0146         thePropagationMomentum = &theMomentum;
0147       std::swap(theMomentum, rhs.theMomentum);
0148       std::swap(theFrozenMomentum, rhs.theFrozenMomentum);
0149       std::swap(thePosition, rhs.thePosition);
0150       std::swap(nCollisions, rhs.nCollisions);
0151       std::swap(nDecays, rhs.nDecays);
0152       std::swap(nSrcPair, rhs.nSrcPair),
0153       std::swap(thePotentialEnergy, rhs.thePotentialEnergy);
0154       // ID intentionally not swapped
0155
0156 #ifdef INCLXX_IN_GEANT4_MODE
0157       std::swap(theParentResonancePDGCode, rhs.theParentResonancePDGCode);
0158       std::swap(theParentResonanceID, rhs.theParentResonanceID);
0159 #endif
0160
0161       std::swap(theHelicity, rhs.theHelicity);
0162       std::swap(emissionTime, rhs.emissionTime);
0163       std::swap(outOfWell, rhs.outOfWell);
0164       std::swap(theSrcPartner, rhs.theSrcPartner);
0165
0166       std::swap(theMass, rhs.theMass);
0167       std::swap(rpCorrelated, rhs.rpCorrelated);
0168       std::swap(uncorrelatedMomentum, rhs.uncorrelatedMomentum);
0169
0170       std::swap(theParticleBias, rhs.theParticleBias);
0171       std::swap(theBiasCollisionVector, rhs.theBiasCollisionVector);
0172
0173     }
0174
0175   public:
0176
0177     /** \brief Assignment operator
0178      *
0179      * Does not copy the particle ID.
0180      */
0181     Particle &operator=(const Particle &rhs) {
0182       Particle temporaryParticle(rhs);
0183       swap(temporaryParticle);
0184       return *this;
0185     }
0186
0187     /**
0188      * Get the particle type.
0189      * @see G4INCL::ParticleType
0190      */
0191     G4INCL::ParticleType getType() const {
0192       return theType;
0193     };
0194
0195     /// \brief Get the particle species
0196     virtual G4INCL::ParticleSpecies getSpecies() const {
0197       return ParticleSpecies(theType);
0198     };
0199
0200     void setType(ParticleType t) {
0201       theType = t;
0202       switch(theType)
0203       {
0204         case DeltaPlusPlus:
0205           theA = 1;
0206           theZ = 2;
0207           theS = 0;
0208           break;
0209         case Proton:
0210         case DeltaPlus:
0211           theA = 1;
0212           theZ = 1;
0213           theS = 0;
0214           break;
0215         case Neutron:
0216         case DeltaZero:
0217           theA = 1;
0218           theZ = 0;
0219           theS = 0;
0220           break;
0221         case DeltaMinus:
0222           theA = 1;
0223           theZ = -1;
0224           theS = 0;
0225           break;
0226         case PiPlus:
0227           theA = 0;
0228           theZ = 1;
0229           theS = 0;
0230           break;
0231         case PiZero:
0232         case Eta:
0233         case Omega:
0234         case EtaPrime:
0235         case Photon:
0236           theA = 0;
0237           theZ = 0;
0238           theS = 0;
0239           break;
0240         case PiMinus:
0241           theA = 0;
0242           theZ = -1;
0243           theS = 0;
0244           break;
0245         case Lambda:
0246           theA = 1;
0247           theZ = 0;
0248           theS = -1;
0249           break;
0250         case SigmaPlus:
0251           theA = 1;
0252           theZ = 1;
0253           theS = -1;
0254           break;
0255         case SigmaZero:
0256           theA = 1;
0257           theZ = 0;
0258           theS = -1;
0259           break;
0260         case SigmaMinus:
0261           theA = 1;
0262           theZ = -1;
0263           theS = -1;
0264           break;
0265         case antiProton:
0266           theA = -1;
0267           theZ = -1;
0268           theS = 0;
0269           break;
0270         case XiMinus:
0271           theA = 1;
0272           theZ = -1;
0273           theS = -2;
0274           break;
0275         case XiZero:
0276           theA = 1;
0277           theZ = 0;
0278           theS = -2;
0279           break;
0280         case antiNeutron:
0281           theA = -1;
0282           theZ = 0;
0283           theS = 0;
0284           break;
0285         case antiLambda:
0286           theA = -1;
0287           theZ = 0;
0288           theS = 1;
0289           break;
0290         case antiSigmaMinus:
0291           theA = -1;
0292           theZ = 1;
0293           theS = 1;
0294           break;
0295         case antiSigmaPlus:
0296           theA = -1;
0297           theZ = -1;
0298           theS = 1;
0299           break;
0300         case antiSigmaZero:
0301           theA = -1;
0302           theZ = 0;
0303           theS = 1;
0304           break;
0305         case antiXiMinus:
0306           theA = -1;
0307           theZ = 1;
0308           theS = 2;
0309           break;
0310         case antiXiZero:
0311           theA = -1;
0312           theZ = 0;
0313           theS = 2;
0314           break;
0315         case KPlus:
0316           theA = 0;
0317           theZ = 1;
0318           theS = 1;
0319           break;
0320         case KZero:
0321           theA = 0;
0322           theZ = 0;
0323           theS = 1;
0324           break;
0325         case KZeroBar:
0326           theA = 0;
0327           theZ = 0;
0328           theS = -1;
0329           break;
0330         case KShort:
0331           theA = 0;
0332           theZ = 0;
0333 //        theS should not be defined
0334           break;
0335         case KLong:
0336           theA = 0;
0337           theZ = 0;
0338 //        theS should not be defined
0339           break;
0340         case KMinus:
0341           theA = 0;
0342           theZ = -1;
0343           theS = -1;
0344           break;
0345         case Composite:
0346          // INCL_ERROR("Trying to set particle type to Composite! Construct a Cluster object instead" << '\n');
0347           theA = 0;
0348           theZ = 0;
0349           theS = 0;
0350           break;
0351         case antiComposite:
0352           theA = 0;
0353           theZ = 0;
0354           theS = 0;
0355           break;
0356         case UnknownParticle:
0357           theA = 0;
0358           theZ = 0;
0359           theS = 0;
0360           INCL_ERROR("Trying to set particle type to Unknown!" << '\n');
0361           break;
0362       }
0363
0364       if( !isResonance() && t!=Composite && t!=antiComposite )
0365         setINCLMass();
0366     }
0367
0368     /**
0369      * Is this a nucleon?
0370      */
0371     G4bool isNucleon() const {
0372       if(theType == G4INCL::Proton || theType == G4INCL::Neutron)
0373     return true;
0374       else
0375     return false;
0376     };
0377
0378     ParticipantType getParticipantType() const {
0379       return theParticipantType;
0380     }
0381
0382     void setParticipantType(ParticipantType const p) {
0383       theParticipantType = p;
0384     }
0385
0386     G4bool isParticipant() const {
0387       return (theParticipantType==Participant);
0388     }
0389
0390     G4bool isTargetSpectator() const {
0391       return (theParticipantType==TargetSpectator);
0392     }
0393
0394     G4bool isProjectileSpectator() const {
0395       return (theParticipantType==ProjectileSpectator);
0396     }
0397
0398     virtual void makeParticipant() {
0399       theParticipantType = Participant;
0400     }
0401
0402     virtual void makeTargetSpectator() {
0403       theParticipantType = TargetSpectator;
0404     }
0405
0406     virtual void makeProjectileSpectator() {
0407       theParticipantType = ProjectileSpectator;
0408     }
0409
0410     /** \brief Is this a pion? */
0411     G4bool isPion() const { return (theType == PiPlus || theType == PiZero || theType == PiMinus); }
0412
0413     /** \brief Is this an eta? */
0414     G4bool isEta() const { return (theType == Eta); }
0415
0416     /** \brief Is this an omega? */
0417     G4bool isOmega() const { return (theType == Omega); }
0418
0419     /** \brief Is this an etaprime? */
0420     G4bool isEtaPrime() const { return (theType == EtaPrime); }
0421
0422     /** \brief Is this a photon? */
0423     G4bool isPhoton() const { return (theType == Photon); }
0424
0425     /** \brief Is it a resonance? */
0426     inline G4bool isResonance() const { return isDelta(); }
0427
0428     /** \brief Is it a Delta? */
0429     inline G4bool isDelta() const {
0430       return (theType==DeltaPlusPlus || theType==DeltaPlus ||
0431           theType==DeltaZero || theType==DeltaMinus); }
0432
0433     /** \brief Is this a Sigma? */
0434     G4bool isSigma() const { return (theType == SigmaPlus || theType == SigmaZero || theType == SigmaMinus); }
0435
0436     /** \brief Is this a Kaon? */
0437     G4bool isKaon() const { return (theType == KPlus || theType == KZero); }
0438
0439     /** \brief Is this an antiKaon? */
0440     G4bool isAntiKaon() const { return (theType == KZeroBar || theType == KMinus); }
0441
0442     /** \brief Is this a Lambda? */
0443     G4bool isLambda() const { return (theType == Lambda); }
0444
0445     /** \brief Is this a Nucleon or a Lambda? */
0446     G4bool isNucleonorLambda() const { return (isNucleon() || isLambda()); }
0447
0448     /** \brief Is this an Hyperon? */
0449     G4bool isHyperon() const { return (isLambda() || isSigma() ); } //|| isXi()
0450
0451     /** \brief Is this a Meson? */
0452     G4bool isMeson() const { return (isPion() || isKaon() || isAntiKaon() || isEta() || isEtaPrime() || isOmega()); }
0453
0454     /** \brief Is this a Baryon? */
0455     G4bool isBaryon() const { return (isNucleon() || isResonance() || isHyperon()); }
0456
0457     /** \brief Is this a Strange? */
0458     G4bool isStrange() const { return (isKaon() || isAntiKaon() || isHyperon()); }
0459
0460     /** \brief Is this a Xi? */
0461     G4bool isXi() const { return (theType == XiZero || theType == XiMinus); }
0462
0463     /** \brief Is this an antinucleon? */
0464     G4bool isAntiNucleon() const { return (theType == antiProton || theType == antiNeutron); }
0465
0466     /** \brief Is this an antiSigma? */
0467     G4bool isAntiSigma() const { return (theType == antiSigmaPlus || theType == antiSigmaZero || theType == antiSigmaMinus); }
0468
0469     /** \brief Is this an antiXi? */
0470     G4bool isAntiXi() const { return (theType == antiXiZero || theType == antiXiMinus); }
0471
0472     /** \brief Is this an antiLambda? */
0473     G4bool isAntiLambda() const { return (theType == antiLambda); }
0474
0475     /** \brief Is this an antiHyperon? */
0476     G4bool isAntiHyperon() const { return (isAntiLambda() || isAntiSigma() || isAntiXi()); }
0477
0478     /** \brief Is this an antiBaryon? */
0479     G4bool isAntiBaryon() const { return (isAntiNucleon() || isAntiHyperon()); }
0480
0481     /** \brief Is this an antiNucleon or an antiLambda? */
0482     G4bool isAntiNucleonorAntiLambda() const { return (isAntiNucleon() || isAntiLambda()); }
0483
0484     /** \brief Returns the baryon number. */
0485     G4int getA() const { return theA; }
0486
0487     /** \brief Returns the charge number. */
0488     G4int getZ() const { return theZ; }
0489
0490     /** \brief Returns the strangeness number. */
0491     G4int getS() const { return theS; }
0492
0493     /** \brief Returns the strangeness number. */
0494     G4int getSrcPair() const { return nSrcPair; }
0495
0496     G4double getBeta() const {
0497       const G4double P = theMomentum.mag();
0498       return P/theEnergy;
0499     }
0500
0501     /**
0502      * Returns a three vector we can give to the boost() -method.
0503      *
0504      * In order to go to the particle rest frame you need to multiply
0505      * the boost vector by -1.0.
0506      */
0507     ThreeVector boostVector() const {
0508       return theMomentum / theEnergy;
0509     }
0510
0511     /**
0512      * Boost the particle using a boost vector.
0513      *
0514      * Example (go to the particle rest frame):
0515      * particle->boost(particle->boostVector());
0516      */
0517     void boost(const ThreeVector &aBoostVector) {
0518       const G4double beta2 = aBoostVector.mag2();
0519       const G4double gamma = 1.0 / std::sqrt(1.0 - beta2);
0520       const G4double bp = theMomentum.dot(aBoostVector);
0521       const G4double alpha = (gamma*gamma)/(1.0 + gamma);
0522
0523       theMomentum = theMomentum + aBoostVector * (alpha * bp - gamma * theEnergy);
0524       theEnergy = gamma * (theEnergy - bp);
0525     }
0526
0527     /** \brief Lorentz-contract the particle position around some center
0528      *
0529      * Apply Lorentz contraction to the position component along the
0530      * direction of the boost vector.
0531      *
0532      * \param aBoostVector the boost vector (velocity) [c]
0533      * \param refPos the reference position
0534      */
0535     void lorentzContract(const ThreeVector &aBoostVector, const ThreeVector &refPos) {
0536       const G4double beta2 = aBoostVector.mag2();
0537       const G4double gamma = 1.0 / std::sqrt(1.0 - beta2);
0538       const ThreeVector theRelativePosition = thePosition - refPos;
0539       const ThreeVector transversePosition = theRelativePosition - aBoostVector * (theRelativePosition.dot(aBoostVector) / aBoostVector.mag2());
0540       const ThreeVector longitudinalPosition = theRelativePosition - transversePosition;
0541
0542       thePosition = refPos + transversePosition + longitudinalPosition / gamma;
0543     }
0544
0545     /** \brief Get the cached particle mass. */
0546     inline G4double getMass() const { return theMass; }
0547
0548     /** \brief Get the INCL particle mass. */
0549     inline G4double getINCLMass() const {
0550       switch(theType) {
0551         case Proton:
0552         case Neutron:
0553         case PiPlus:
0554         case PiMinus:
0555         case PiZero:
0556         case Lambda:
0557         case SigmaPlus:
0558         case SigmaZero:
0559         case SigmaMinus:
0560         case antiProton:
0561         case XiZero:
0562         case XiMinus:
0563         case antiNeutron:
0564         case antiLambda:
0565         case antiSigmaPlus:
0566         case antiSigmaZero:
0567         case antiSigmaMinus:
0568         case antiXiZero:
0569         case antiXiMinus:
0570         case KPlus:
0571         case KZero:
0572         case KZeroBar:
0573         case KShort:
0574         case KLong:
0575         case KMinus:
0576         case Eta:
0577         case Omega:
0578         case EtaPrime:
0579         case Photon:
0580           return ParticleTable::getINCLMass(theType);
0581           break;
0582
0583         case DeltaPlusPlus:
0584         case DeltaPlus:
0585         case DeltaZero:
0586         case DeltaMinus:
0587           return theMass;
0588           break;
0589
0590         case Composite:
0591           return ParticleTable::getINCLMass(theA,theZ,theS);
0592           break;
0593         case antiComposite:
0594           return ParticleTable::getINCLMass(-theA,-theZ,theS);
0595           break;
0596
0597         default:
0598           INCL_ERROR("Particle::getINCLMass: Unknown particle type." << '\n');
0599           return 0.0;
0600           break;
0601       }
0602     }
0603
0604     /** \brief Get the tabulated particle mass. */
0605     inline virtual G4double getTableMass() const {
0606       switch(theType) {
0607         case Proton:
0608         case Neutron:
0609         case PiPlus:
0610         case PiMinus:
0611         case PiZero:
0612         case Lambda:
0613         case SigmaPlus:
0614         case SigmaZero:
0615         case SigmaMinus:
0616         case antiProton:
0617         case XiZero:
0618         case XiMinus:
0619         case antiNeutron:
0620         case antiLambda:
0621         case antiSigmaPlus:
0622         case antiSigmaZero:
0623         case antiSigmaMinus:
0624         case antiXiZero:
0625         case antiXiMinus:
0626         case KPlus:
0627         case KZero:
0628         case KZeroBar:
0629         case KShort:
0630         case KLong:
0631         case KMinus:
0632         case Eta:
0633         case Omega:
0634         case EtaPrime:
0635         case Photon:
0636           return ParticleTable::getTableParticleMass(theType);
0637           break;
0638
0639         case DeltaPlusPlus:
0640         case DeltaPlus:
0641         case DeltaZero:
0642         case DeltaMinus:
0643           return theMass;
0644           break;
0645
0646         case Composite:
0647           return ParticleTable::getTableMass(theA,theZ,theS);
0648           break;
0649         case antiComposite:
0650           return ParticleTable::getTableMass(-theA,-theZ,theS);
0651           break;
0652
0653         default:
0654           INCL_ERROR("Particle::getTableMass: Unknown particle type." << '\n');
0655           return 0.0;
0656           break;
0657       }
0658     }
0659
0660     /** \brief Get the real particle mass. */
0661     inline G4double getRealMass() const {
0662       switch(theType) {
0663         case Proton:
0664         case Neutron:
0665         case PiPlus:
0666         case PiMinus:
0667         case PiZero:
0668         case Lambda:
0669         case SigmaPlus:
0670         case SigmaZero:
0671         case SigmaMinus:
0672         case antiProton:
0673         case XiZero:
0674         case XiMinus:
0675         case antiNeutron:
0676         case antiLambda:
0677         case antiSigmaPlus:
0678         case antiSigmaZero:
0679         case antiSigmaMinus:
0680         case antiXiZero:
0681         case antiXiMinus:
0682         case KPlus:
0683         case KZero:
0684         case KZeroBar:
0685         case KShort:
0686         case KLong:
0687         case KMinus:
0688         case Eta:
0689         case Omega:
0690         case EtaPrime:
0691         case Photon:
0692           return ParticleTable::getRealMass(theType);
0693           break;
0694
0695         case DeltaPlusPlus:
0696         case DeltaPlus:
0697         case DeltaZero:
0698         case DeltaMinus:
0699           return theMass;
0700           break;
0701
0702         case Composite:
0703           return ParticleTable::getRealMass(theA,theZ,theS);
0704           break;
0705         case antiComposite:
0706           return ParticleTable::getRealMass(-theA,-theZ,theS);
0707           break;
0708
0709         default:
0710           INCL_ERROR("Particle::getRealMass: Unknown particle type." << '\n');
0711           return 0.0;
0712           break;
0713       }
0714     }
0715
0716     /// \brief Set the mass of the Particle to its real mass
0717     void setRealMass() { setMass(getRealMass()); }
0718
0719     /// \brief Set the mass of the Particle to its table mass
0720     void setTableMass() { setMass(getTableMass()); }
0721
0722     /// \brief Set the mass of the Particle to its table mass
0723     void setINCLMass() { setMass(getINCLMass()); }
0724
0725     /**\brief Computes correction on the emission Q-value
0726      *
0727      * Computes the correction that must be applied to INCL particles in
0728      * order to obtain the correct Q-value for particle emission from a given
0729      * nucleus. For absorption, the correction is obviously equal to minus
0730      * the value returned by this function.
0731      *
0732      * \param AParent the mass number of the emitting nucleus
0733      * \param ZParent the charge number of the emitting nucleus
0734      * \return the correction
0735      */
0736     G4double getEmissionQValueCorrection(const G4int AParent, const G4int ZParent) const {
0737       const G4int SParent = 0;
0738       const G4int ADaughter = AParent - theA;
0739       const G4int ZDaughter = ZParent - theZ;
0740       const G4int SDaughter = 0;
0741
0742       // Note the minus sign here
0743       G4double theQValue;
0744       if(isCluster())
0745         theQValue = -ParticleTable::getTableQValue(theA, theZ, theS, ADaughter, ZDaughter, SDaughter);
0746       else {
0747         const G4double massTableParent = ParticleTable::getTableMass(AParent,ZParent,SParent);
0748         const G4double massTableDaughter = ParticleTable::getTableMass(ADaughter,ZDaughter,SDaughter);
0749         const G4double massTableParticle = getTableMass();
0750         theQValue = massTableParent - massTableDaughter - massTableParticle;
0751       }
0752
0753       const G4double massINCLParent = ParticleTable::getINCLMass(AParent,ZParent,SParent);
0754       const G4double massINCLDaughter = ParticleTable::getINCLMass(ADaughter,ZDaughter,SDaughter);
0755       const G4double massINCLParticle = getINCLMass();
0756
0757       // The rhs corresponds to the INCL Q-value
0758       return theQValue - (massINCLParent-massINCLDaughter-massINCLParticle);
0759     }
0760
0761     /**\brief Computes correction on the transfer Q-value
0762      *
0763      * Computes the correction that must be applied to INCL particles in
0764      * order to obtain the correct Q-value for particle transfer from a given
0765      * nucleus to another.
0766      *
0767      * Assumes that the receving nucleus is INCL's target nucleus, with the
0768      * INCL separation energy.
0769      *
0770      * \param AFrom the mass number of the donating nucleus
0771      * \param ZFrom the charge number of the donating nucleus
0772      * \param ATo the mass number of the receiving nucleus
0773      * \param ZTo the charge number of the receiving nucleus
0774      * \return the correction
0775      */
0776     G4double getTransferQValueCorrection(const G4int AFrom, const G4int ZFrom, const G4int ATo, const G4int ZTo) const {
0777       const G4int SFrom = 0;
0778       const G4int STo = 0;
0779       const G4int AFromDaughter = AFrom - theA;
0780       const G4int ZFromDaughter = ZFrom - theZ;
0781       const G4int SFromDaughter = 0;
0782       const G4int AToDaughter = ATo + theA;
0783       const G4int ZToDaughter = ZTo + theZ;
0784       const G4int SToDaughter = 0;
0785       const G4double theQValue = ParticleTable::getTableQValue(AToDaughter,ZToDaughter,SToDaughter,AFromDaughter,ZFromDaughter,SFromDaughter,AFrom,ZFrom,SFrom);
0786
0787       const G4double massINCLTo = ParticleTable::getINCLMass(ATo,ZTo,STo);
0788       const G4double massINCLToDaughter = ParticleTable::getINCLMass(AToDaughter,ZToDaughter,SToDaughter);
0789       /* Note that here we have to use the table mass in the INCL Q-value. We
0790        * cannot use theMass, because at this stage the particle is probably
0791        * still off-shell; and we cannot use getINCLMass(), because it leads to
0792        * violations of global energy conservation.
0793        */
0794       const G4double massINCLParticle = getTableMass();
0795
0796       // The rhs corresponds to the INCL Q-value for particle absorption
0797       return theQValue - (massINCLToDaughter-massINCLTo-massINCLParticle);
0798     }
0799
0800     /**\brief Computes correction on the emission Q-value for hypernuclei
0801      *
0802      * Computes the correction that must be applied to INCL particles in
0803      * order to obtain the correct Q-value for particle emission from a given
0804      * nucleus. For absorption, the correction is obviously equal to minus
0805      * the value returned by this function.
0806      *
0807      * \param AParent the mass number of the emitting nucleus
0808      * \param ZParent the charge number of the emitting nucleus
0809      * \param SParent the strangess number of the emitting nucleus
0810      * \return the correction
0811      */
0812     G4double getEmissionQValueCorrection(const G4int AParent, const G4int ZParent, const G4int SParent) const {
0813       const G4int ADaughter = AParent - theA;
0814       const G4int ZDaughter = ZParent - theZ;
0815       const G4int SDaughter = SParent - theS;
0816
0817       // Note the minus sign here
0818       G4double theQValue;
0819       if(isCluster())
0820         theQValue = -ParticleTable::getTableQValue(theA, theZ, theS, ADaughter, ZDaughter, SDaughter);
0821       else {
0822         const G4double massTableParent = ParticleTable::getTableMass(AParent,ZParent,SParent);
0823         const G4double massTableDaughter = ParticleTable::getTableMass(ADaughter,ZDaughter,SDaughter);
0824         const G4double massTableParticle = getTableMass();
0825         theQValue = massTableParent - massTableDaughter - massTableParticle;
0826       }
0827
0828       const G4double massINCLParent = ParticleTable::getINCLMass(AParent,ZParent,SParent);
0829       const G4double massINCLDaughter = ParticleTable::getINCLMass(ADaughter,ZDaughter,SDaughter);
0830       const G4double massINCLParticle = getINCLMass();
0831
0832       // The rhs corresponds to the INCL Q-value
0833       return theQValue - (massINCLParent-massINCLDaughter-massINCLParticle);
0834     }
0835
0836     /**\brief Computes correction on the transfer Q-value for hypernuclei
0837      *
0838      * Computes the correction that must be applied to INCL particles in
0839      * order to obtain the correct Q-value for particle transfer from a given
0840      * nucleus to another.
0841      *
0842      * Assumes that the receving nucleus is INCL's target nucleus, with the
0843      * INCL separation energy.
0844      *
0845      * \param AFrom the mass number of the donating nucleus
0846      * \param ZFrom the charge number of the donating nucleus
0847      * \param SFrom the strangess number of the donating nucleus
0848      * \param ATo the mass number of the receiving nucleus
0849      * \param ZTo the charge number of the receiving nucleus
0850      * \param STo the strangess number of the receiving nucleus
0851      * \return the correction
0852      */
0853     G4double getTransferQValueCorrection(const G4int AFrom, const G4int ZFrom, const G4int SFrom, const G4int ATo, const G4int ZTo , const G4int STo) const {
0854       const G4int AFromDaughter = AFrom - theA;
0855       const G4int ZFromDaughter = ZFrom - theZ;
0856       const G4int SFromDaughter = SFrom - theS;
0857       const G4int AToDaughter = ATo + theA;
0858       const G4int ZToDaughter = ZTo + theZ;
0859       const G4int SToDaughter = STo + theS;
0860       const G4double theQValue = ParticleTable::getTableQValue(AToDaughter,ZToDaughter,SFromDaughter,AFromDaughter,ZFromDaughter,SToDaughter,AFrom,ZFrom,SFrom);
0861
0862       const G4double massINCLTo = ParticleTable::getINCLMass(ATo,ZTo,STo);
0863       const G4double massINCLToDaughter = ParticleTable::getINCLMass(AToDaughter,ZToDaughter,SToDaughter);
0864       /* Note that here we have to use the table mass in the INCL Q-value. We
0865        * cannot use theMass, because at this stage the particle is probably
0866        * still off-shell; and we cannot use getINCLMass(), because it leads to
0867        * violations of global energy conservation.
0868        */
0869       const G4double massINCLParticle = getTableMass();
0870
0871       // The rhs corresponds to the INCL Q-value for particle absorption
0872       return theQValue - (massINCLToDaughter-massINCLTo-massINCLParticle);
0873     }
0874
0875
0876
0877     /** \brief Get the the particle invariant mass.
0878      *
0879      * Uses the relativistic invariant
0880      * \f[ m = \sqrt{E^2 - {\vec p}^2}\f]
0881      **/
0882     G4double getInvariantMass() const {
0883       const G4double mass = std::pow(theEnergy, 2) - theMomentum.dot(theMomentum);
0884       if(mass < 0.0) {
0885         INCL_ERROR("E*E - p*p is negative." << '\n');
0886         return 0.0;
0887       } else {
0888         return std::sqrt(mass);
0889       }
0890     };
0891
0892     /// \brief Get the particle kinetic energy.
0893     inline G4double getKineticEnergy() const { return theEnergy - theMass; }
0894
0895     /// \brief Get the particle potential energy.
0896     inline G4double getPotentialEnergy() const { return thePotentialEnergy; }
0897
0898     /// \brief Set the particle potential energy.
0899     inline void setPotentialEnergy(G4double v) { thePotentialEnergy = v; }
0900
0901     /**
0902      * Get the energy of the particle in MeV.
0903      */
0904     G4double getEnergy() const
0905     {
0906       return theEnergy;
0907     };
0908
0909     /**
0910      * Set the mass of the particle in MeV/c^2.
0911      */
0912     void setMass(G4double mass)
0913     {
0914       this->theMass = mass;
0915     }
0916
0917     /**
0918      * Set the energy of the particle in MeV.
0919      */
0920     void setEnergy(G4double energy)
0921     {
0922       this->theEnergy = energy;
0923     };
0924
0925     /**
0926      * Get the momentum vector.
0927      */
0928     const G4INCL::ThreeVector &getMomentum() const
0929     {
0930       return theMomentum;
0931     };
0932
0933     /** Get the angular momentum w.r.t. the origin */
0934     virtual G4INCL::ThreeVector getAngularMomentum() const
0935     {
0936       return thePosition.vector(theMomentum);
0937     };
0938
0939     /**
0940      * Set the momentum vector.
0941      */
0942     virtual void setMomentum(const G4INCL::ThreeVector &momentum)
0943     {
0944       this->theMomentum = momentum;
0945     };
0946
0947     /**
0948      * Set the position vector.
0949      */
0950     const G4INCL::ThreeVector &getPosition() const
0951     {
0952       return thePosition;
0953     };
0954
0955     virtual void setPosition(const G4INCL::ThreeVector &position)
0956     {
0957       this->thePosition = position;
0958     };
0959
0960     G4double getHelicity() { return theHelicity; };
0961     void setHelicity(G4double h) { theHelicity = h; };
0962
0963     void propagate(G4double step) {
0964       thePosition += ((*thePropagationMomentum)*(step/(*thePropagationEnergy)));
0965     };
0966
0967     /** \brief Return the number of collisions undergone by the particle. **/
0968     G4int getNumberOfCollisions() const { return nCollisions; }
0969
0970     /** \brief Set the number of collisions undergone by the particle. **/
0971     void setNumberOfCollisions(G4int n) { nCollisions = n; }
0972
0973     /** \brief Increment the number of collisions undergone by the particle. **/
0974     void incrementNumberOfCollisions() { nCollisions++; }
0975
0976     /** \brief Return the number of decays undergone by the particle. **/
0977     G4int getNumberOfDecays() const { return nDecays; }
0978
0979     /** \brief Set the number of decays undergone by the particle. **/
0980     void setNumberOfDecays(G4int n) { nDecays = n; }
0981
0982     /** \brief Increment the number of decays undergone by the particle. **/
0983     void incrementNumberOfDecays() { nDecays++; }
0984
0985     /** \brief Set the number of srcpairs. **/
0986     void setNumberOfSrcPair(int n) { nSrcPair = n; }
0987
0988     /** \brief Mark the particle as out of its potential well
0989      *
0990      * This flag is used to control pions created outside their potential well
0991      * in delta decay. The pion potential checks it and returns zero if it is
0992      * true (necessary in order to correctly enforce energy conservation). The
0993      * Nucleus::applyFinalState() method uses it to determine whether new
0994      * avatars should be generated for the particle.
0995      */
0996     void setOutOfWell() { outOfWell = true; }
0997
0998     /// \brief Check if the particle is out of its potential well
0999     G4bool isOutOfWell() const { return outOfWell; }
1000
1001     /// \brief Set and reset src partner
1002     void setSrcPartner() { theSrcPartner = true; }
1003     void resetSrcPartner() { theSrcPartner = false; nSrcPair=0; }
1004
1005     /// \brief Check if the particle is a src partner
1006     G4bool isSrcPartner() const { return theSrcPartner; }
1007
1008     void setEmissionTime(G4double t) { emissionTime = t; }
1009     G4double getEmissionTime() { return emissionTime; };
1010
1011     /** \brief Transverse component of the position w.r.t. the momentum. */
1012     ThreeVector getTransversePosition() const {
1013       return thePosition - getLongitudinalPosition();
1014     }
1015
1016     /** \brief Longitudinal component of the position w.r.t. the momentum. */
1017     ThreeVector getLongitudinalPosition() const {
1018       return *thePropagationMomentum * (thePosition.dot(*thePropagationMomentum)/thePropagationMomentum->mag2());
1019     }
1020
1021     /** \brief Rescale the momentum to match the total energy. */
1022     const ThreeVector &adjustMomentumFromEnergy();
1023
1024     /** \brief Recompute the energy to match the momentum. */
1025     G4double adjustEnergyFromMomentum();
1026
1027     G4bool isCluster() const {
1028       return ((theType == Composite || theType == antiComposite));
1029     }
1030
1031     /// \brief Set the frozen particle momentum
1032     void setFrozenMomentum(const ThreeVector &momentum) { theFrozenMomentum = momentum; }
1033
1034     /// \brief Set the frozen particle momentum
1035     void setFrozenEnergy(const G4double energy) { theFrozenEnergy = energy; }
1036
1037     /// \brief Get the frozen particle momentum
1038     ThreeVector getFrozenMomentum() const { return theFrozenMomentum; }
1039
1040     /// \brief Get the frozen particle momentum
1041     G4double getFrozenEnergy() const { return theFrozenEnergy; }
1042
1043     /// \brief Get the propagation velocity of the particle
1044     ThreeVector getPropagationVelocity() const { return (*thePropagationMomentum)/(*thePropagationEnergy); }
1045
1046     /** \brief Freeze particle propagation
1047      *
1048      * Make the particle use theFrozenMomentum and theFrozenEnergy for
1049      * propagation. The normal state can be restored by calling the
1050      * thawPropagation() method.
1051      */
1052     void freezePropagation() {
1053       thePropagationMomentum = &theFrozenMomentum;
1054       thePropagationEnergy = &theFrozenEnergy;
1055     }
1056
1057     /** \brief Unfreeze particle propagation
1058      *
1059      * Make the particle use theMomentum and theEnergy for propagation. Call
1060      * this method to restore the normal propagation if the
1061      * freezePropagation() method has been called.
1062      */
1063     void thawPropagation() {
1064       thePropagationMomentum = &theMomentum;
1065       thePropagationEnergy = &theEnergy;
1066     }
1067
1068     /** \brief Rotate the particle position and momentum
1069      *
1070      * \param angle the rotation angle
1071      * \param axis a unit vector representing the rotation axis
1072      */
1073     virtual void rotatePositionAndMomentum(const G4double angle, const ThreeVector &axis) {
1074       rotatePosition(angle, axis);
1075       rotateMomentum(angle, axis);
1076     }
1077
1078     /** \brief Rotate the particle position
1079      *
1080      * \param angle the rotation angle
1081      * \param axis a unit vector representing the rotation axis
1082      */
1083     virtual void rotatePosition(const G4double angle, const ThreeVector &axis) {
1084       thePosition.rotate(angle, axis);
1085     }
1086
1087     /** \brief Rotate the particle momentum
1088      *
1089      * \param angle the rotation angle
1090      * \param axis a unit vector representing the rotation axis
1091      */
1092     virtual void rotateMomentum(const G4double angle, const ThreeVector &axis) {
1093       theMomentum.rotate(angle, axis);
1094       theFrozenMomentum.rotate(angle, axis);
1095     }
1096
1097     std::string print() const {
1098       std::stringstream ss;
1099       ss << "Particle (ID = " << ID << ") type = ";
1100       ss << ParticleTable::getName(theType);
1101       ss << ", SRC pair = " << nSrcPair;
1102       ss << ", Potential energy = " << thePotentialEnergy;
1103       ss << '\n'
1104         << "   energy = " << theEnergy << '\n'
1105         << "   momentum = "
1106         << theMomentum.print()
1107         << '\n'
1108         << "   position = "
1109         << thePosition.print()
1110         << '\n';
1111       return ss.str();
1112     };
1113
1114     std::string dump() const {
1115       std::stringstream ss;
1116       ss << "(particle " << ID << " ";
1117       ss << ParticleTable::getName(theType);
1118       ss << nSrcPair << " ";
1119       ss << '\n'
1120         << thePosition.dump()
1121         << '\n'
1122         << theMomentum.dump()
1123         << '\n'
1124         << theEnergy << ")" << '\n';
1125       return ss.str();
1126     };
1127
1128     long getID() const { return ID; };
1129
1130     /**
1131      * Return a NULL pointer
1132      */
1133     ParticleList const *getParticles() const {
1134       INCL_WARN("Particle::getParticles() method was called on a Particle object" << '\n');
1135       return 0;
1136     }
1137
1138     /** \brief Return the reflection momentum
1139      *
1140      * The reflection momentum is used by calls to getSurfaceRadius to compute
1141      * the radius of the sphere where the nucleon moves. It is necessary to
1142      * introduce fuzzy r-p correlations.
1143      */
1144     G4double getReflectionMomentum() const {
1145       if(rpCorrelated)
1146         return theMomentum.mag();
1147       else
1148         return uncorrelatedMomentum;
1149     }
1150
1151     /// \brief Set the uncorrelated momentum
1152     void setUncorrelatedMomentum(const G4double p) { uncorrelatedMomentum = p; }
1153
1154     /// \brief Make the particle follow a strict r-p correlation
1155     void rpCorrelate() { rpCorrelated = true; }
1156
1157     /// \brief Make the particle not follow a strict r-p correlation
1158     void rpDecorrelate() { rpCorrelated = false; }
1159
1160     /// \brief Get the cosine of the angle between position and momentum
1161     G4double getCosRPAngle() const {
1162       const G4double norm = thePosition.mag2()*thePropagationMomentum->mag2();
1163       if(norm>0.)
1164         return thePosition.dot(*thePropagationMomentum) / std::sqrt(norm);
1165       else
1166         return 1.;
1167     }
1168
1169     /// \brief General bias vector function
1170     static G4double getTotalBias();
1171     static void setINCLBiasVector(std::vector<G4double> NewVector);
1172     static void FillINCLBiasVector(G4double newBias);
1173     static G4double getBiasFromVector(std::vector<G4int> VectorBias);
1174
1175     static std::vector<G4int> MergeVectorBias(Particle const * const p1, Particle const * const p2);
1176     static std::vector<G4int> MergeVectorBias(std::vector<G4int> p1, Particle const * const p2);
1177
1178     /// \brief Get the particle bias.
1179     G4double getParticleBias() const { return theParticleBias; };
1180
1181     /// \brief Set the particle bias.
1182     void setParticleBias(G4double ParticleBias) { this->theParticleBias = ParticleBias; }
1183
1184     /// \brief Get the vector list of biased vertices on the particle path.
1185     std::vector<G4int> getBiasCollisionVector() const { return theBiasCollisionVector; }
1186
1187     /// \brief Set the vector list of biased vertices on the particle path.
1188     void setBiasCollisionVector(std::vector<G4int> BiasCollisionVector) {
1189       this->theBiasCollisionVector = BiasCollisionVector;
1190       this->setParticleBias(Particle::getBiasFromVector(std::move(BiasCollisionVector)));
1191       }
1192
1193     /** \brief Number of Kaon inside de nucleus
1194      *
1195      * Put in the Particle class in order to calculate the
1196      * "correct" mass of composit particle.
1197      *
1198      */
1199
1200     G4int getNumberOfKaon() const { return theNKaon; };
1201     void setNumberOfKaon(const G4int NK) { theNKaon = NK; }
1202
1203 #ifdef INCLXX_IN_GEANT4_MODE
1204     G4int getParentResonancePDGCode() const { return theParentResonancePDGCode; };
1205     void setParentResonancePDGCode(const G4int parentPDGCode) { theParentResonancePDGCode = parentPDGCode; };
1206     G4int getParentResonanceID() const { return theParentResonanceID; };
1207     void setParentResonanceID(const G4int parentID) { theParentResonanceID = parentID; };
1208 #endif
1209
1210   public:
1211     /** \brief Time ordered vector of all bias applied
1212      *
1213      * /!\ Caution /!\
1214      * methods Assotiated to G4VectorCache<T> are:
1215      * Push_back(…),
1216      * operator[],
1217      * Begin(),
1218      * End(),
1219      * Clear(),
1220      * Size() and
1221      * Pop_back()
1222      *
1223      */
1224 #ifdef INCLXX_IN_GEANT4_MODE
1225       static std::vector<G4double> INCLBiasVector;
1226       //static G4VectorCache<G4double> INCLBiasVector;
1227 #else
1228       static G4ThreadLocal std::vector<G4double> INCLBiasVector;
1229       //static G4VectorCache<G4double> INCLBiasVector;
1230 #endif
1231     static G4ThreadLocal G4int nextBiasedCollisionID;
1232
1233   protected:
1234     G4int theZ, theA, theS;
1235     ParticipantType theParticipantType;
1236     G4INCL::ParticleType theType;
1237     G4double theEnergy;
1238     G4double *thePropagationEnergy;
1239     G4double theFrozenEnergy;
1240     G4INCL::ThreeVector theMomentum;
1241     G4INCL::ThreeVector *thePropagationMomentum;
1242     G4INCL::ThreeVector theFrozenMomentum;
1243     G4INCL::ThreeVector thePosition;
1244     G4int nCollisions;
1245     G4int nDecays;
1246     G4int nSrcPair;
1247     G4double thePotentialEnergy;
1248     long ID;
1249
1250     G4bool rpCorrelated;
1251     G4double uncorrelatedMomentum;
1252
1253     G4double theParticleBias;
1254     /// \brief The number of Kaons inside the nucleus (update during the cascade)
1255     G4int theNKaon;
1256
1257 #ifdef INCLXX_IN_GEANT4_MODE
1258     G4int theParentResonancePDGCode;
1259     G4int theParentResonanceID;
1260 #endif
1261
1262   private:
1263     G4double theHelicity;
1264     G4double emissionTime;
1265     G4bool outOfWell;
1266     G4bool theSrcPartner;
1267
1268     /// \brief Time ordered vector of all biased vertices on the particle path
1269     std::vector<G4int> theBiasCollisionVector;
1270
1271     G4double theMass;
1272     static G4ThreadLocal long nextID;
1273
1274     INCL_DECLARE_ALLOCATION_POOL(Particle)
1275   };
1276 }
1277
1278 #endif /* PARTICLE_HH_ */