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       thePotentialEnergy(rhs.thePotentialEnergy),
0100       rpCorrelated(rhs.rpCorrelated),
0101       uncorrelatedMomentum(rhs.uncorrelatedMomentum),
0102       theParticleBias(rhs.theParticleBias),
0103       theNKaon(rhs.theNKaon),
0104 #ifdef INCLXX_IN_GEANT4_MODE
0105       theParentResonancePDGCode(rhs.theParentResonancePDGCode),
0106       theParentResonanceID(rhs.theParentResonanceID),
0107 #endif
0108       theHelicity(rhs.theHelicity),
0109       emissionTime(rhs.emissionTime),
0110       outOfWell(rhs.outOfWell),
0111       theMass(rhs.theMass)
0112       {
0113         if(rhs.thePropagationEnergy == &(rhs.theFrozenEnergy))
0114           thePropagationEnergy = &theFrozenEnergy;
0115         else
0116           thePropagationEnergy = &theEnergy;
0117         if(rhs.thePropagationMomentum == &(rhs.theFrozenMomentum))
0118           thePropagationMomentum = &theFrozenMomentum;
0119         else
0120           thePropagationMomentum = &theMomentum;
0121         // ID intentionally not copied
0122         ID = nextID++;
0123
0124         theBiasCollisionVector = rhs.theBiasCollisionVector;
0125       }
0126
0127   protected:
0128     /// \brief Helper method for the assignment operator
0129     void swap(Particle &rhs) {
0130       std::swap(theZ, rhs.theZ);
0131       std::swap(theA, rhs.theA);
0132       std::swap(theS, rhs.theS);
0133       std::swap(theParticipantType, rhs.theParticipantType);
0134       std::swap(theType, rhs.theType);
0135       if(rhs.thePropagationEnergy == &(rhs.theFrozenEnergy))
0136         thePropagationEnergy = &theFrozenEnergy;
0137       else
0138         thePropagationEnergy = &theEnergy;
0139       std::swap(theEnergy, rhs.theEnergy);
0140       std::swap(theFrozenEnergy, rhs.theFrozenEnergy);
0141       if(rhs.thePropagationMomentum == &(rhs.theFrozenMomentum))
0142         thePropagationMomentum = &theFrozenMomentum;
0143       else
0144         thePropagationMomentum = &theMomentum;
0145       std::swap(theMomentum, rhs.theMomentum);
0146       std::swap(theFrozenMomentum, rhs.theFrozenMomentum);
0147       std::swap(thePosition, rhs.thePosition);
0148       std::swap(nCollisions, rhs.nCollisions);
0149       std::swap(nDecays, rhs.nDecays);
0150       std::swap(thePotentialEnergy, rhs.thePotentialEnergy);
0151       // ID intentionally not swapped
0152
0153 #ifdef INCLXX_IN_GEANT4_MODE
0154       std::swap(theParentResonancePDGCode, rhs.theParentResonancePDGCode);
0155       std::swap(theParentResonanceID, rhs.theParentResonanceID);
0156 #endif
0157
0158       std::swap(theHelicity, rhs.theHelicity);
0159       std::swap(emissionTime, rhs.emissionTime);
0160       std::swap(outOfWell, rhs.outOfWell);
0161
0162       std::swap(theMass, rhs.theMass);
0163       std::swap(rpCorrelated, rhs.rpCorrelated);
0164       std::swap(uncorrelatedMomentum, rhs.uncorrelatedMomentum);
0165
0166       std::swap(theParticleBias, rhs.theParticleBias);
0167       std::swap(theBiasCollisionVector, rhs.theBiasCollisionVector);
0168
0169     }
0170
0171   public:
0172
0173     /** \brief Assignment operator
0174      *
0175      * Does not copy the particle ID.
0176      */
0177     Particle &operator=(const Particle &rhs) {
0178       Particle temporaryParticle(rhs);
0179       swap(temporaryParticle);
0180       return *this;
0181     }
0182
0183     /**
0184      * Get the particle type.
0185      * @see G4INCL::ParticleType
0186      */
0187     G4INCL::ParticleType getType() const {
0188       return theType;
0189     };
0190
0191     /// \brief Get the particle species
0192     virtual G4INCL::ParticleSpecies getSpecies() const {
0193       return ParticleSpecies(theType);
0194     };
0195
0196     void setType(ParticleType t) {
0197       theType = t;
0198       switch(theType)
0199       {
0200         case DeltaPlusPlus:
0201           theA = 1;
0202           theZ = 2;
0203           theS = 0;
0204           break;
0205         case Proton:
0206         case DeltaPlus:
0207           theA = 1;
0208           theZ = 1;
0209           theS = 0;
0210           break;
0211         case Neutron:
0212         case DeltaZero:
0213           theA = 1;
0214           theZ = 0;
0215           theS = 0;
0216           break;
0217         case DeltaMinus:
0218           theA = 1;
0219           theZ = -1;
0220           theS = 0;
0221           break;
0222         case PiPlus:
0223           theA = 0;
0224           theZ = 1;
0225           theS = 0;
0226           break;
0227         case PiZero:
0228         case Eta:
0229         case Omega:
0230         case EtaPrime:
0231         case Photon:
0232           theA = 0;
0233           theZ = 0;
0234           theS = 0;
0235           break;
0236         case PiMinus:
0237           theA = 0;
0238           theZ = -1;
0239           theS = 0;
0240           break;
0241         case Lambda:
0242           theA = 1;
0243           theZ = 0;
0244           theS = -1;
0245           break;
0246         case SigmaPlus:
0247           theA = 1;
0248           theZ = 1;
0249           theS = -1;
0250           break;
0251         case SigmaZero:
0252           theA = 1;
0253           theZ = 0;
0254           theS = -1;
0255           break;
0256         case SigmaMinus:
0257           theA = 1;
0258           theZ = -1;
0259           theS = -1;
0260           break;
0261         case antiProton:
0262           theA = -1;
0263           theZ = -1;
0264           theS = 0;
0265           break;
0266         case XiMinus:
0267           theA = 1;
0268           theZ = -1;
0269           theS = -2;
0270           break;
0271         case XiZero:
0272           theA = 1;
0273           theZ = 0;
0274           theS = -2;
0275           break;
0276         case antiNeutron:
0277           theA = -1;
0278           theZ = 0;
0279           theS = 0;
0280           break;
0281         case antiLambda:
0282           theA = -1;
0283           theZ = 0;
0284           theS = 1;
0285           break;
0286         case antiSigmaMinus:
0287           theA = -1;
0288           theZ = 1;
0289           theS = 1;
0290           break;
0291         case antiSigmaPlus:
0292           theA = -1;
0293           theZ = -1;
0294           theS = 1;
0295           break;
0296         case antiSigmaZero:
0297           theA = -1;
0298           theZ = 0;
0299           theS = 1;
0300           break;
0301         case antiXiMinus:
0302           theA = -1;
0303           theZ = 1;
0304           theS = 2;
0305           break;
0306         case antiXiZero:
0307           theA = -1;
0308           theZ = 0;
0309           theS = 2;
0310           break;
0311         case KPlus:
0312           theA = 0;
0313           theZ = 1;
0314           theS = 1;
0315           break;
0316         case KZero:
0317           theA = 0;
0318           theZ = 0;
0319           theS = 1;
0320           break;
0321         case KZeroBar:
0322           theA = 0;
0323           theZ = 0;
0324           theS = -1;
0325           break;
0326         case KShort:
0327           theA = 0;
0328           theZ = 0;
0329 //        theS should not be defined
0330           break;
0331         case KLong:
0332           theA = 0;
0333           theZ = 0;
0334 //        theS should not be defined
0335           break;
0336         case KMinus:
0337           theA = 0;
0338           theZ = -1;
0339           theS = -1;
0340           break;
0341         case Composite:
0342          // INCL_ERROR("Trying to set particle type to Composite! Construct a Cluster object instead" << '\n');
0343           theA = 0;
0344           theZ = 0;
0345           theS = 0;
0346           break;
0347         case UnknownParticle:
0348           theA = 0;
0349           theZ = 0;
0350           theS = 0;
0351           INCL_ERROR("Trying to set particle type to Unknown!" << '\n');
0352           break;
0353       }
0354
0355       if( !isResonance() && t!=Composite )
0356         setINCLMass();
0357     }
0358
0359     /**
0360      * Is this a nucleon?
0361      */
0362     G4bool isNucleon() const {
0363       if(theType == G4INCL::Proton || theType == G4INCL::Neutron)
0364     return true;
0365       else
0366     return false;
0367     };
0368
0369     ParticipantType getParticipantType() const {
0370       return theParticipantType;
0371     }
0372
0373     void setParticipantType(ParticipantType const p) {
0374       theParticipantType = p;
0375     }
0376
0377     G4bool isParticipant() const {
0378       return (theParticipantType==Participant);
0379     }
0380
0381     G4bool isTargetSpectator() const {
0382       return (theParticipantType==TargetSpectator);
0383     }
0384
0385     G4bool isProjectileSpectator() const {
0386       return (theParticipantType==ProjectileSpectator);
0387     }
0388
0389     virtual void makeParticipant() {
0390       theParticipantType = Participant;
0391     }
0392
0393     virtual void makeTargetSpectator() {
0394       theParticipantType = TargetSpectator;
0395     }
0396
0397     virtual void makeProjectileSpectator() {
0398       theParticipantType = ProjectileSpectator;
0399     }
0400
0401     /** \brief Is this a pion? */
0402     G4bool isPion() const { return (theType == PiPlus || theType == PiZero || theType == PiMinus); }
0403
0404     /** \brief Is this an eta? */
0405     G4bool isEta() const { return (theType == Eta); }
0406
0407     /** \brief Is this an omega? */
0408     G4bool isOmega() const { return (theType == Omega); }
0409
0410     /** \brief Is this an etaprime? */
0411     G4bool isEtaPrime() const { return (theType == EtaPrime); }
0412
0413     /** \brief Is this a photon? */
0414     G4bool isPhoton() const { return (theType == Photon); }
0415
0416     /** \brief Is it a resonance? */
0417     inline G4bool isResonance() const { return isDelta(); }
0418
0419     /** \brief Is it a Delta? */
0420     inline G4bool isDelta() const {
0421       return (theType==DeltaPlusPlus || theType==DeltaPlus ||
0422           theType==DeltaZero || theType==DeltaMinus); }
0423
0424     /** \brief Is this a Sigma? */
0425     G4bool isSigma() const { return (theType == SigmaPlus || theType == SigmaZero || theType == SigmaMinus); }
0426
0427     /** \brief Is this a Kaon? */
0428     G4bool isKaon() const { return (theType == KPlus || theType == KZero); }
0429
0430     /** \brief Is this an antiKaon? */
0431     G4bool isAntiKaon() const { return (theType == KZeroBar || theType == KMinus); }
0432
0433     /** \brief Is this a Lambda? */
0434     G4bool isLambda() const { return (theType == Lambda); }
0435
0436     /** \brief Is this a Nucleon or a Lambda? */
0437     G4bool isNucleonorLambda() const { return (isNucleon() || isLambda()); }
0438
0439     /** \brief Is this an Hyperon? */
0440     G4bool isHyperon() const { return (isLambda() || isSigma() ); } //|| isXi()
0441
0442     /** \brief Is this a Meson? */
0443     G4bool isMeson() const { return (isPion() || isKaon() || isAntiKaon() || isEta() || isEtaPrime() || isOmega()); }
0444
0445     /** \brief Is this a Baryon? */
0446     G4bool isBaryon() const { return (isNucleon() || isResonance() || isHyperon()); }
0447
0448     /** \brief Is this a Strange? */
0449     G4bool isStrange() const { return (isKaon() || isAntiKaon() || isHyperon()); }
0450
0451     /** \brief Is this a Xi? */
0452     G4bool isXi() const { return (theType == XiZero || theType == XiMinus); }
0453
0454     /** \brief Is this an antinucleon? */
0455     G4bool isAntiNucleon() const { return (theType == antiProton || theType == antiNeutron); }
0456
0457     /** \brief Is this an antiSigma? */
0458     G4bool isAntiSigma() const { return (theType == antiSigmaPlus || theType == antiSigmaZero || theType == antiSigmaMinus); }
0459
0460     /** \brief Is this an antiXi? */
0461     G4bool isAntiXi() const { return (theType == antiXiZero || theType == antiXiMinus); }
0462
0463     /** \brief Is this an antiLambda? */
0464     G4bool isAntiLambda() const { return (theType == antiLambda); }
0465
0466     /** \brief Is this an antiHyperon? */
0467     G4bool isAntiHyperon() const { return (isAntiLambda() || isAntiSigma() || isAntiXi()); }
0468
0469     /** \brief Is this an antiBaryon? */
0470     G4bool isAntiBaryon() const { return (isAntiNucleon() || isAntiHyperon()); }
0471
0472     /** \brief Is this an antiNucleon or an antiLambda? */
0473     G4bool isAntiNucleonorAntiLambda() const { return (isAntiNucleon() || isAntiLambda()); }
0474
0475     /** \brief Returns the baryon number. */
0476     G4int getA() const { return theA; }
0477
0478     /** \brief Returns the charge number. */
0479     G4int getZ() const { return theZ; }
0480
0481     /** \brief Returns the strangeness number. */
0482     G4int getS() const { return theS; }
0483
0484     G4double getBeta() const {
0485       const G4double P = theMomentum.mag();
0486       return P/theEnergy;
0487     }
0488
0489     /**
0490      * Returns a three vector we can give to the boost() -method.
0491      *
0492      * In order to go to the particle rest frame you need to multiply
0493      * the boost vector by -1.0.
0494      */
0495     ThreeVector boostVector() const {
0496       return theMomentum / theEnergy;
0497     }
0498
0499     /**
0500      * Boost the particle using a boost vector.
0501      *
0502      * Example (go to the particle rest frame):
0503      * particle->boost(particle->boostVector());
0504      */
0505     void boost(const ThreeVector &aBoostVector) {
0506       const G4double beta2 = aBoostVector.mag2();
0507       const G4double gamma = 1.0 / std::sqrt(1.0 - beta2);
0508       const G4double bp = theMomentum.dot(aBoostVector);
0509       const G4double alpha = (gamma*gamma)/(1.0 + gamma);
0510
0511       theMomentum = theMomentum + aBoostVector * (alpha * bp - gamma * theEnergy);
0512       theEnergy = gamma * (theEnergy - bp);
0513     }
0514
0515     /** \brief Lorentz-contract the particle position around some center
0516      *
0517      * Apply Lorentz contraction to the position component along the
0518      * direction of the boost vector.
0519      *
0520      * \param aBoostVector the boost vector (velocity) [c]
0521      * \param refPos the reference position
0522      */
0523     void lorentzContract(const ThreeVector &aBoostVector, const ThreeVector &refPos) {
0524       const G4double beta2 = aBoostVector.mag2();
0525       const G4double gamma = 1.0 / std::sqrt(1.0 - beta2);
0526       const ThreeVector theRelativePosition = thePosition - refPos;
0527       const ThreeVector transversePosition = theRelativePosition - aBoostVector * (theRelativePosition.dot(aBoostVector) / aBoostVector.mag2());
0528       const ThreeVector longitudinalPosition = theRelativePosition - transversePosition;
0529
0530       thePosition = refPos + transversePosition + longitudinalPosition / gamma;
0531     }
0532
0533     /** \brief Get the cached particle mass. */
0534     inline G4double getMass() const { return theMass; }
0535
0536     /** \brief Get the INCL particle mass. */
0537     inline G4double getINCLMass() const {
0538       switch(theType) {
0539         case Proton:
0540         case Neutron:
0541         case PiPlus:
0542         case PiMinus:
0543         case PiZero:
0544         case Lambda:
0545         case SigmaPlus:
0546         case SigmaZero:
0547         case SigmaMinus:
0548         case antiProton:
0549         case XiZero:
0550         case XiMinus:
0551         case antiNeutron:
0552         case antiLambda:
0553         case antiSigmaPlus:
0554         case antiSigmaZero:
0555         case antiSigmaMinus:
0556         case antiXiZero:
0557         case antiXiMinus:
0558         case KPlus:
0559         case KZero:
0560         case KZeroBar:
0561         case KShort:
0562         case KLong:
0563         case KMinus:
0564         case Eta:
0565         case Omega:
0566         case EtaPrime:
0567         case Photon:
0568           return ParticleTable::getINCLMass(theType);
0569           break;
0570
0571         case DeltaPlusPlus:
0572         case DeltaPlus:
0573         case DeltaZero:
0574         case DeltaMinus:
0575           return theMass;
0576           break;
0577
0578         case Composite:
0579           return ParticleTable::getINCLMass(theA,theZ,theS);
0580           break;
0581
0582         default:
0583           INCL_ERROR("Particle::getINCLMass: Unknown particle type." << '\n');
0584           return 0.0;
0585           break;
0586       }
0587     }
0588
0589     /** \brief Get the tabulated particle mass. */
0590     inline virtual G4double getTableMass() const {
0591       switch(theType) {
0592         case Proton:
0593         case Neutron:
0594         case PiPlus:
0595         case PiMinus:
0596         case PiZero:
0597         case Lambda:
0598         case SigmaPlus:
0599         case SigmaZero:
0600         case SigmaMinus:
0601         case antiProton:
0602         case XiZero:
0603         case XiMinus:
0604         case antiNeutron:
0605         case antiLambda:
0606         case antiSigmaPlus:
0607         case antiSigmaZero:
0608         case antiSigmaMinus:
0609         case antiXiZero:
0610         case antiXiMinus:
0611         case KPlus:
0612         case KZero:
0613         case KZeroBar:
0614         case KShort:
0615         case KLong:
0616         case KMinus:
0617         case Eta:
0618         case Omega:
0619         case EtaPrime:
0620         case Photon:
0621           return ParticleTable::getTableParticleMass(theType);
0622           break;
0623
0624         case DeltaPlusPlus:
0625         case DeltaPlus:
0626         case DeltaZero:
0627         case DeltaMinus:
0628           return theMass;
0629           break;
0630
0631         case Composite:
0632           return ParticleTable::getTableMass(theA,theZ,theS);
0633           break;
0634
0635         default:
0636           INCL_ERROR("Particle::getTableMass: Unknown particle type." << '\n');
0637           return 0.0;
0638           break;
0639       }
0640     }
0641
0642     /** \brief Get the real particle mass. */
0643     inline G4double getRealMass() const {
0644       switch(theType) {
0645         case Proton:
0646         case Neutron:
0647         case PiPlus:
0648         case PiMinus:
0649         case PiZero:
0650         case Lambda:
0651         case SigmaPlus:
0652         case SigmaZero:
0653         case SigmaMinus:
0654         case antiProton:
0655         case XiZero:
0656         case XiMinus:
0657         case antiNeutron:
0658         case antiLambda:
0659         case antiSigmaPlus:
0660         case antiSigmaZero:
0661         case antiSigmaMinus:
0662         case antiXiZero:
0663         case antiXiMinus:
0664         case KPlus:
0665         case KZero:
0666         case KZeroBar:
0667         case KShort:
0668         case KLong:
0669         case KMinus:
0670         case Eta:
0671         case Omega:
0672         case EtaPrime:
0673         case Photon:
0674           return ParticleTable::getRealMass(theType);
0675           break;
0676
0677         case DeltaPlusPlus:
0678         case DeltaPlus:
0679         case DeltaZero:
0680         case DeltaMinus:
0681           return theMass;
0682           break;
0683
0684         case Composite:
0685           return ParticleTable::getRealMass(theA,theZ,theS);
0686           break;
0687
0688         default:
0689           INCL_ERROR("Particle::getRealMass: Unknown particle type." << '\n');
0690           return 0.0;
0691           break;
0692       }
0693     }
0694
0695     /// \brief Set the mass of the Particle to its real mass
0696     void setRealMass() { setMass(getRealMass()); }
0697
0698     /// \brief Set the mass of the Particle to its table mass
0699     void setTableMass() { setMass(getTableMass()); }
0700
0701     /// \brief Set the mass of the Particle to its table mass
0702     void setINCLMass() { setMass(getINCLMass()); }
0703
0704     /**\brief Computes correction on the emission Q-value
0705      *
0706      * Computes the correction that must be applied to INCL particles in
0707      * order to obtain the correct Q-value for particle emission from a given
0708      * nucleus. For absorption, the correction is obviously equal to minus
0709      * the value returned by this function.
0710      *
0711      * \param AParent the mass number of the emitting nucleus
0712      * \param ZParent the charge number of the emitting nucleus
0713      * \return the correction
0714      */
0715     G4double getEmissionQValueCorrection(const G4int AParent, const G4int ZParent) const {
0716       const G4int SParent = 0;
0717       const G4int ADaughter = AParent - theA;
0718       const G4int ZDaughter = ZParent - theZ;
0719       const G4int SDaughter = 0;
0720
0721       // Note the minus sign here
0722       G4double theQValue;
0723       if(isCluster())
0724         theQValue = -ParticleTable::getTableQValue(theA, theZ, theS, ADaughter, ZDaughter, SDaughter);
0725       else {
0726         const G4double massTableParent = ParticleTable::getTableMass(AParent,ZParent,SParent);
0727         const G4double massTableDaughter = ParticleTable::getTableMass(ADaughter,ZDaughter,SDaughter);
0728         const G4double massTableParticle = getTableMass();
0729         theQValue = massTableParent - massTableDaughter - massTableParticle;
0730       }
0731
0732       const G4double massINCLParent = ParticleTable::getINCLMass(AParent,ZParent,SParent);
0733       const G4double massINCLDaughter = ParticleTable::getINCLMass(ADaughter,ZDaughter,SDaughter);
0734       const G4double massINCLParticle = getINCLMass();
0735
0736       // The rhs corresponds to the INCL Q-value
0737       return theQValue - (massINCLParent-massINCLDaughter-massINCLParticle);
0738     }
0739
0740     G4double getEmissionPbarQvalueCorrection(const G4int AParent, const G4int ZParent, const G4bool Victim) const {
0741       G4int SParent = 0;
0742       G4int SDaughter = 0;
0743       G4int ADaughter = AParent - 1;
0744       G4int ZDaughter;
0745       G4bool isProton = Victim;
0746       if(isProton){     //proton is annihilated
0747         ZDaughter = ZParent - 1;
0748       }
0749       else {       //neutron is annihilated
0750         ZDaughter = ZParent;
0751       }
0752
0753       G4double theQValue; //same procedure as for normal case
0754
0755       const G4double massTableParent = ParticleTable::getTableMass(AParent,ZParent,SParent);
0756       const G4double massTableDaughter = ParticleTable::getTableMass(ADaughter,ZDaughter,SDaughter);
0757       const G4double massTableParticle = getTableMass();
0758       theQValue = massTableParent - massTableDaughter - massTableParticle;
0759
0760       const G4double massINCLParent = ParticleTable::getINCLMass(AParent,ZParent,SParent);
0761       const G4double massINCLDaughter = ParticleTable::getINCLMass(ADaughter,ZDaughter,SDaughter);
0762       const G4double massINCLParticle = getINCLMass();
0763
0764       return theQValue - (massINCLParent-massINCLDaughter-massINCLParticle);
0765     }
0766
0767     /**\brief Computes correction on the transfer Q-value
0768      *
0769      * Computes the correction that must be applied to INCL particles in
0770      * order to obtain the correct Q-value for particle transfer from a given
0771      * nucleus to another.
0772      *
0773      * Assumes that the receving nucleus is INCL's target nucleus, with the
0774      * INCL separation energy.
0775      *
0776      * \param AFrom the mass number of the donating nucleus
0777      * \param ZFrom the charge number of the donating nucleus
0778      * \param ATo the mass number of the receiving nucleus
0779      * \param ZTo the charge number of the receiving nucleus
0780      * \return the correction
0781      */
0782     G4double getTransferQValueCorrection(const G4int AFrom, const G4int ZFrom, const G4int ATo, const G4int ZTo) const {
0783       const G4int SFrom = 0;
0784       const G4int STo = 0;
0785       const G4int AFromDaughter = AFrom - theA;
0786       const G4int ZFromDaughter = ZFrom - theZ;
0787       const G4int SFromDaughter = 0;
0788       const G4int AToDaughter = ATo + theA;
0789       const G4int ZToDaughter = ZTo + theZ;
0790       const G4int SToDaughter = 0;
0791       const G4double theQValue = ParticleTable::getTableQValue(AToDaughter,ZToDaughter,SToDaughter,AFromDaughter,ZFromDaughter,SFromDaughter,AFrom,ZFrom,SFrom);
0792
0793       const G4double massINCLTo = ParticleTable::getINCLMass(ATo,ZTo,STo);
0794       const G4double massINCLToDaughter = ParticleTable::getINCLMass(AToDaughter,ZToDaughter,SToDaughter);
0795       /* Note that here we have to use the table mass in the INCL Q-value. We
0796        * cannot use theMass, because at this stage the particle is probably
0797        * still off-shell; and we cannot use getINCLMass(), because it leads to
0798        * violations of global energy conservation.
0799        */
0800       const G4double massINCLParticle = getTableMass();
0801
0802       // The rhs corresponds to the INCL Q-value for particle absorption
0803       return theQValue - (massINCLToDaughter-massINCLTo-massINCLParticle);
0804     }
0805
0806     /**\brief Computes correction on the emission Q-value for hypernuclei
0807      *
0808      * Computes the correction that must be applied to INCL particles in
0809      * order to obtain the correct Q-value for particle emission from a given
0810      * nucleus. For absorption, the correction is obviously equal to minus
0811      * the value returned by this function.
0812      *
0813      * \param AParent the mass number of the emitting nucleus
0814      * \param ZParent the charge number of the emitting nucleus
0815      * \param SParent the strangess number of the emitting nucleus
0816      * \return the correction
0817      */
0818     G4double getEmissionQValueCorrection(const G4int AParent, const G4int ZParent, const G4int SParent) const {
0819       const G4int ADaughter = AParent - theA;
0820       const G4int ZDaughter = ZParent - theZ;
0821       const G4int SDaughter = SParent - theS;
0822
0823       // Note the minus sign here
0824       G4double theQValue;
0825       if(isCluster())
0826         theQValue = -ParticleTable::getTableQValue(theA, theZ, theS, ADaughter, ZDaughter, SDaughter);
0827       else {
0828         const G4double massTableParent = ParticleTable::getTableMass(AParent,ZParent,SParent);
0829         const G4double massTableDaughter = ParticleTable::getTableMass(ADaughter,ZDaughter,SDaughter);
0830         const G4double massTableParticle = getTableMass();
0831         theQValue = massTableParent - massTableDaughter - massTableParticle;
0832       }
0833
0834       const G4double massINCLParent = ParticleTable::getINCLMass(AParent,ZParent,SParent);
0835       const G4double massINCLDaughter = ParticleTable::getINCLMass(ADaughter,ZDaughter,SDaughter);
0836       const G4double massINCLParticle = getINCLMass();
0837
0838       // The rhs corresponds to the INCL Q-value
0839       return theQValue - (massINCLParent-massINCLDaughter-massINCLParticle);
0840     }
0841
0842     /**\brief Computes correction on the transfer Q-value for hypernuclei
0843      *
0844      * Computes the correction that must be applied to INCL particles in
0845      * order to obtain the correct Q-value for particle transfer from a given
0846      * nucleus to another.
0847      *
0848      * Assumes that the receving nucleus is INCL's target nucleus, with the
0849      * INCL separation energy.
0850      *
0851      * \param AFrom the mass number of the donating nucleus
0852      * \param ZFrom the charge number of the donating nucleus
0853      * \param SFrom the strangess number of the donating nucleus
0854      * \param ATo the mass number of the receiving nucleus
0855      * \param ZTo the charge number of the receiving nucleus
0856      * \param STo the strangess number of the receiving nucleus
0857      * \return the correction
0858      */
0859     G4double getTransferQValueCorrection(const G4int AFrom, const G4int ZFrom, const G4int SFrom, const G4int ATo, const G4int ZTo , const G4int STo) const {
0860       const G4int AFromDaughter = AFrom - theA;
0861       const G4int ZFromDaughter = ZFrom - theZ;
0862       const G4int SFromDaughter = SFrom - theS;
0863       const G4int AToDaughter = ATo + theA;
0864       const G4int ZToDaughter = ZTo + theZ;
0865       const G4int SToDaughter = STo + theS;
0866       const G4double theQValue = ParticleTable::getTableQValue(AToDaughter,ZToDaughter,SFromDaughter,AFromDaughter,ZFromDaughter,SToDaughter,AFrom,ZFrom,SFrom);
0867
0868       const G4double massINCLTo = ParticleTable::getINCLMass(ATo,ZTo,STo);
0869       const G4double massINCLToDaughter = ParticleTable::getINCLMass(AToDaughter,ZToDaughter,SToDaughter);
0870       /* Note that here we have to use the table mass in the INCL Q-value. We
0871        * cannot use theMass, because at this stage the particle is probably
0872        * still off-shell; and we cannot use getINCLMass(), because it leads to
0873        * violations of global energy conservation.
0874        */
0875       const G4double massINCLParticle = getTableMass();
0876
0877       // The rhs corresponds to the INCL Q-value for particle absorption
0878       return theQValue - (massINCLToDaughter-massINCLTo-massINCLParticle);
0879     }
0880
0881
0882
0883     /** \brief Get the the particle invariant mass.
0884      *
0885      * Uses the relativistic invariant
0886      * \f[ m = \sqrt{E^2 - {\vec p}^2}\f]
0887      **/
0888     G4double getInvariantMass() const {
0889       const G4double mass = std::pow(theEnergy, 2) - theMomentum.dot(theMomentum);
0890       if(mass < 0.0) {
0891         INCL_ERROR("E*E - p*p is negative." << '\n');
0892         return 0.0;
0893       } else {
0894         return std::sqrt(mass);
0895       }
0896     };
0897
0898     /// \brief Get the particle kinetic energy.
0899     inline G4double getKineticEnergy() const { return theEnergy - theMass; }
0900
0901     /// \brief Get the particle potential energy.
0902     inline G4double getPotentialEnergy() const { return thePotentialEnergy; }
0903
0904     /// \brief Set the particle potential energy.
0905     inline void setPotentialEnergy(G4double v) { thePotentialEnergy = v; }
0906
0907     /**
0908      * Get the energy of the particle in MeV.
0909      */
0910     G4double getEnergy() const
0911     {
0912       return theEnergy;
0913     };
0914
0915     /**
0916      * Set the mass of the particle in MeV/c^2.
0917      */
0918     void setMass(G4double mass)
0919     {
0920       this->theMass = mass;
0921     }
0922
0923     /**
0924      * Set the energy of the particle in MeV.
0925      */
0926     void setEnergy(G4double energy)
0927     {
0928       this->theEnergy = energy;
0929     };
0930
0931     /**
0932      * Get the momentum vector.
0933      */
0934     const G4INCL::ThreeVector &getMomentum() const
0935     {
0936       return theMomentum;
0937     };
0938
0939     /** Get the angular momentum w.r.t. the origin */
0940     virtual G4INCL::ThreeVector getAngularMomentum() const
0941     {
0942       return thePosition.vector(theMomentum);
0943     };
0944
0945     /**
0946      * Set the momentum vector.
0947      */
0948     virtual void setMomentum(const G4INCL::ThreeVector &momentum)
0949     {
0950       this->theMomentum = momentum;
0951     };
0952
0953     /**
0954      * Set the position vector.
0955      */
0956     const G4INCL::ThreeVector &getPosition() const
0957     {
0958       return thePosition;
0959     };
0960
0961     virtual void setPosition(const G4INCL::ThreeVector &position)
0962     {
0963       this->thePosition = position;
0964     };
0965
0966     G4double getHelicity() { return theHelicity; };
0967     void setHelicity(G4double h) { theHelicity = h; };
0968
0969     void propagate(G4double step) {
0970       thePosition += ((*thePropagationMomentum)*(step/(*thePropagationEnergy)));
0971     };
0972
0973     /** \brief Return the number of collisions undergone by the particle. **/
0974     G4int getNumberOfCollisions() const { return nCollisions; }
0975
0976     /** \brief Set the number of collisions undergone by the particle. **/
0977     void setNumberOfCollisions(G4int n) { nCollisions = n; }
0978
0979     /** \brief Increment the number of collisions undergone by the particle. **/
0980     void incrementNumberOfCollisions() { nCollisions++; }
0981
0982     /** \brief Return the number of decays undergone by the particle. **/
0983     G4int getNumberOfDecays() const { return nDecays; }
0984
0985     /** \brief Set the number of decays undergone by the particle. **/
0986     void setNumberOfDecays(G4int n) { nDecays = n; }
0987
0988     /** \brief Increment the number of decays undergone by the particle. **/
0989     void incrementNumberOfDecays() { nDecays++; }
0990
0991     /** \brief Mark the particle as out of its potential well
0992      *
0993      * This flag is used to control pions created outside their potential well
0994      * in delta decay. The pion potential checks it and returns zero if it is
0995      * true (necessary in order to correctly enforce energy conservation). The
0996      * Nucleus::applyFinalState() method uses it to determine whether new
0997      * avatars should be generated for the particle.
0998      */
0999     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(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_ */