G4INCLParticle.hh

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLParticle.hh

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
// INCL++ intra-nuclear cascade model
// Alain Boudard, CEA-Saclay, France
// Joseph Cugnon, University of Liege, Belgium
// Jean-Christophe David, CEA-Saclay, France
// Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
// Sylvie Leray, CEA-Saclay, France
// Davide Mancusi, CEA-Saclay, France
//
#define INCLXX_IN_GEANT4_MODE 1

#include "globals.hh"

/*
 * G4INCLParticle.hh
 *
 *  \date Jun 5, 2009
 * \author Pekka Kaitaniemi
 */

#ifndef PARTICLE_HH_
#define PARTICLE_HH_

#include "G4INCLThreeVector.hh"
#include "G4INCLParticleTable.hh"
#include "G4INCLParticleType.hh"
#include "G4INCLParticleSpecies.hh"
#include "G4INCLLogger.hh"
#include "G4INCLUnorderedVector.hh"
#include "G4INCLAllocationPool.hh"
#include <sstream>
#include <string>

namespace G4INCL {

  class Particle;

  class ParticleList : public UnorderedVector<Particle*> {
    public:
      void rotatePositionAndMomentum(const G4double angle, const ThreeVector &axis) const;
      void rotatePosition(const G4double angle, const ThreeVector &axis) const;
      void rotateMomentum(const G4double angle, const ThreeVector &axis) const;
      void boost(const ThreeVector &b) const;
      G4double getParticleListBias() const;
      std::vector<G4int> getParticleListBiasVector() const;
  };

  typedef ParticleList::const_iterator ParticleIter;
  typedef ParticleList::iterator       ParticleMutableIter;

  class Particle {
  public:
    Particle();
    Particle(ParticleType t, G4double energy, ThreeVector const &momentum, ThreeVector const &position);
    Particle(ParticleType t, ThreeVector const &momentum, ThreeVector const &position);
    virtual ~Particle() {}

    Particle(const Particle &rhs) :
      theZ(rhs.theZ),
      theA(rhs.theA),
      theS(rhs.theS),
      theParticipantType(rhs.theParticipantType),
      theType(rhs.theType),
      theEnergy(rhs.theEnergy),
      theFrozenEnergy(rhs.theFrozenEnergy),
      theMomentum(rhs.theMomentum),
      theFrozenMomentum(rhs.theFrozenMomentum),
      thePosition(rhs.thePosition),
      nCollisions(rhs.nCollisions),
      nDecays(rhs.nDecays),
      thePotentialEnergy(rhs.thePotentialEnergy),
      rpCorrelated(rhs.rpCorrelated),
      uncorrelatedMomentum(rhs.uncorrelatedMomentum),
      theParticleBias(rhs.theParticleBias),
      theNKaon(rhs.theNKaon),
      theHelicity(rhs.theHelicity),
      emissionTime(rhs.emissionTime),
      outOfWell(rhs.outOfWell),
      theMass(rhs.theMass)
      {
        if(rhs.thePropagationEnergy == &(rhs.theFrozenEnergy))
          thePropagationEnergy = &theFrozenEnergy;
        else
          thePropagationEnergy = &theEnergy;
        if(rhs.thePropagationMomentum == &(rhs.theFrozenMomentum))
          thePropagationMomentum = &theFrozenMomentum;
        else
          thePropagationMomentum = &theMomentum;
        // ID intentionally not copied
        ID = nextID++;
        
        theBiasCollisionVector = rhs.theBiasCollisionVector;
      }

  protected:
    void swap(Particle &rhs) {
      std::swap(theZ, rhs.theZ);
      std::swap(theA, rhs.theA);
      std::swap(theS, rhs.theS);
      std::swap(theParticipantType, rhs.theParticipantType);
      std::swap(theType, rhs.theType);
      if(rhs.thePropagationEnergy == &(rhs.theFrozenEnergy))
        thePropagationEnergy = &theFrozenEnergy;
      else
        thePropagationEnergy = &theEnergy;
      std::swap(theEnergy, rhs.theEnergy);
      std::swap(theFrozenEnergy, rhs.theFrozenEnergy);
      if(rhs.thePropagationMomentum == &(rhs.theFrozenMomentum))
        thePropagationMomentum = &theFrozenMomentum;
      else
        thePropagationMomentum = &theMomentum;
      std::swap(theMomentum, rhs.theMomentum);
      std::swap(theFrozenMomentum, rhs.theFrozenMomentum);
      std::swap(thePosition, rhs.thePosition);
      std::swap(nCollisions, rhs.nCollisions);
      std::swap(nDecays, rhs.nDecays);
      std::swap(thePotentialEnergy, rhs.thePotentialEnergy);
      // ID intentionally not swapped

      std::swap(theHelicity, rhs.theHelicity);
      std::swap(emissionTime, rhs.emissionTime);
      std::swap(outOfWell, rhs.outOfWell);

      std::swap(theMass, rhs.theMass);
      std::swap(rpCorrelated, rhs.rpCorrelated);
      std::swap(uncorrelatedMomentum, rhs.uncorrelatedMomentum);
      
      std::swap(theParticleBias, rhs.theParticleBias);
      std::swap(theBiasCollisionVector, rhs.theBiasCollisionVector);

    }

  public:

    Particle &operator=(const Particle &rhs) {
      Particle temporaryParticle(rhs);
      swap(temporaryParticle);
      return *this;
    }

    G4INCL::ParticleType getType() const {
      return theType;
    };

    virtual G4INCL::ParticleSpecies getSpecies() const {
      return ParticleSpecies(theType);
    };

    void setType(ParticleType t) {
      theType = t;
      switch(theType)
      {
        case DeltaPlusPlus:
          theA = 1;
          theZ = 2;
          theS = 0;
          break;
        case Proton:
        case DeltaPlus:
          theA = 1;
          theZ = 1;
          theS = 0;
          break;
        case Neutron:
        case DeltaZero:
          theA = 1;
          theZ = 0;
          theS = 0;
          break;
        case DeltaMinus:
          theA = 1;
          theZ = -1;
          theS = 0;
          break;
        case PiPlus:
          theA = 0;
          theZ = 1;
          theS = 0;
          break;
        case PiZero:
        case Eta:
        case Omega:
        case EtaPrime:
        case Photon:
          theA = 0;
          theZ = 0;
          theS = 0;
          break;
        case PiMinus:
          theA = 0;
          theZ = -1;
          theS = 0;
          break;
        case Lambda:
          theA = 1;
          theZ = 0;
          theS = -1;
          break;
        case SigmaPlus:
          theA = 1;
          theZ = 1;
          theS = -1;
          break;
        case SigmaZero:
          theA = 1;
          theZ = 0;
          theS = -1;
          break;
        case SigmaMinus:
          theA = 1;
          theZ = -1;
          theS = -1;
          break;
        case KPlus:
          theA = 0;
          theZ = 1;
          theS = 1;
          break;
        case KZero:
          theA = 0;
          theZ = 0;
          theS = 1;
          break;
        case KZeroBar:
          theA = 0;
          theZ = 0;
          theS = -1;
          break;
        case KShort:
          theA = 0;
          theZ = 0;
//        theS should not be defined
          break;
        case KLong:
          theA = 0;
          theZ = 0;
//        theS should not be defined
          break;
        case KMinus:
          theA = 0;
          theZ = -1;
          theS = -1;
          break;
        case Composite:
         // INCL_ERROR("Trying to set particle type to Composite! Construct a Cluster object instead" << '\n');
          theA = 0;
          theZ = 0;
          theS = 0;
          break;
        case UnknownParticle:
          theA = 0;
          theZ = 0;
          theS = 0;
          INCL_ERROR("Trying to set particle type to Unknown!" << '\n');
          break;
      }

      if( !isResonance() && t!=Composite )
        setINCLMass();
    }

    G4bool isNucleon() const {
      if(theType == G4INCL::Proton || theType == G4INCL::Neutron)
    return true;
      else
    return false;
    };

    ParticipantType getParticipantType() const {
      return theParticipantType;
    }

    void setParticipantType(ParticipantType const p) {
      theParticipantType = p;
    }

    G4bool isParticipant() const {
      return (theParticipantType==Participant);
    }

    G4bool isTargetSpectator() const {
      return (theParticipantType==TargetSpectator);
    }

    G4bool isProjectileSpectator() const {
      return (theParticipantType==ProjectileSpectator);
    }

    virtual void makeParticipant() {
      theParticipantType = Participant;
    }

    virtual void makeTargetSpectator() {
      theParticipantType = TargetSpectator;
    }

    virtual void makeProjectileSpectator() {
      theParticipantType = ProjectileSpectator;
    }

    G4bool isPion() const { return (theType == PiPlus || theType == PiZero || theType == PiMinus); }

    G4bool isEta() const { return (theType == Eta); }

    G4bool isOmega() const { return (theType == Omega); }

    G4bool isEtaPrime() const { return (theType == EtaPrime); }

    G4bool isPhoton() const { return (theType == Photon); }

    inline G4bool isResonance() const { return isDelta(); }

    inline G4bool isDelta() const {
      return (theType==DeltaPlusPlus || theType==DeltaPlus ||
          theType==DeltaZero || theType==DeltaMinus); }
    
    G4bool isSigma() const { return (theType == SigmaPlus || theType == SigmaZero || theType == SigmaMinus); }
    
    G4bool isKaon() const { return (theType == KPlus || theType == KZero); }
    
    G4bool isAntiKaon() const { return (theType == KZeroBar || theType == KMinus); }
    
    G4bool isLambda() const { return (theType == Lambda); }

    G4bool isNucleonorLambda() const { return (isNucleon() || isLambda()); }
    
    G4bool isHyperon() const { return (isLambda() || isSigma()); }
    
    G4bool isMeson() const { return (isPion() || isKaon() || isAntiKaon() || isEta() || isEtaPrime() || isOmega()); }
    
    G4bool isBaryon() const { return (isNucleon() || isResonance() || isHyperon()); }
    
    G4bool isStrange() const { return (isKaon() || isAntiKaon() || isHyperon()); }

    G4int getA() const { return theA; }

    G4int getZ() const { return theZ; }
    
    G4int getS() const { return theS; }

    G4double getBeta() const {
      const G4double P = theMomentum.mag();
      return P/theEnergy;
    }

    ThreeVector boostVector() const {
      return theMomentum / theEnergy;
    }

    void boost(const ThreeVector &aBoostVector) {
      const G4double beta2 = aBoostVector.mag2();
      const G4double gamma = 1.0 / std::sqrt(1.0 - beta2);
      const G4double bp = theMomentum.dot(aBoostVector);
      const G4double alpha = (gamma*gamma)/(1.0 + gamma);

      theMomentum = theMomentum + aBoostVector * (alpha * bp - gamma * theEnergy);
      theEnergy = gamma * (theEnergy - bp);
    }

    void lorentzContract(const ThreeVector &aBoostVector, const ThreeVector &refPos) {
      const G4double beta2 = aBoostVector.mag2();
      const G4double gamma = 1.0 / std::sqrt(1.0 - beta2);
      const ThreeVector theRelativePosition = thePosition - refPos;
      const ThreeVector transversePosition = theRelativePosition - aBoostVector * (theRelativePosition.dot(aBoostVector) / aBoostVector.mag2());
      const ThreeVector longitudinalPosition = theRelativePosition - transversePosition;

      thePosition = refPos + transversePosition + longitudinalPosition / gamma;
    }

    inline G4double getMass() const { return theMass; }

    inline G4double getINCLMass() const {
      switch(theType) {
        case Proton:
        case Neutron:
        case PiPlus:
        case PiMinus:
        case PiZero:
        case Lambda:
        case SigmaPlus:
        case SigmaZero:
        case SigmaMinus:
        case KPlus:
        case KZero:
        case KZeroBar:
        case KShort:
        case KLong:
        case KMinus:
        case Eta:
        case Omega:
        case EtaPrime:
        case Photon:
          return ParticleTable::getINCLMass(theType);
          break;

        case DeltaPlusPlus:
        case DeltaPlus:
        case DeltaZero:
        case DeltaMinus:
          return theMass;
          break;

        case Composite:
          return ParticleTable::getINCLMass(theA,theZ,theS);
          break;

        default:
          INCL_ERROR("Particle::getINCLMass: Unknown particle type." << '\n');
          return 0.0;
          break;
      }
    }

    inline virtual G4double getTableMass() const {
      switch(theType) {
        case Proton:
        case Neutron:
        case PiPlus:
        case PiMinus:
        case PiZero:
        case Lambda:
        case SigmaPlus:
        case SigmaZero:
        case SigmaMinus:
        case KPlus:
        case KZero:
        case KZeroBar:
        case KShort:
        case KLong:
        case KMinus:
        case Eta:
        case Omega:
        case EtaPrime:
        case Photon:
          return ParticleTable::getTableParticleMass(theType);
          break;

        case DeltaPlusPlus:
        case DeltaPlus:
        case DeltaZero:
        case DeltaMinus:
          return theMass;
          break;

        case Composite:
          return ParticleTable::getTableMass(theA,theZ,theS);
          break;

        default:
          INCL_ERROR("Particle::getTableMass: Unknown particle type." << '\n');
          return 0.0;
          break;
      }
    }

    inline G4double getRealMass() const {
      switch(theType) {
        case Proton:
        case Neutron:
        case PiPlus:
        case PiMinus:
        case PiZero:
        case Lambda:
        case SigmaPlus:
        case SigmaZero:
        case SigmaMinus:
        case KPlus:
        case KZero:
        case KZeroBar:
        case KShort:
        case KLong:
        case KMinus:
        case Eta:
        case Omega:
        case EtaPrime:
        case Photon:
          return ParticleTable::getRealMass(theType);
          break;

        case DeltaPlusPlus:
        case DeltaPlus:
        case DeltaZero:
        case DeltaMinus:
          return theMass;
          break;

        case Composite:
          return ParticleTable::getRealMass(theA,theZ,theS);
          break;

        default:
          INCL_ERROR("Particle::getRealMass: Unknown particle type." << '\n');
          return 0.0;
          break;
      }
    }

    void setRealMass() { setMass(getRealMass()); }

    void setTableMass() { setMass(getTableMass()); }

    void setINCLMass() { setMass(getINCLMass()); }

    G4double getEmissionQValueCorrection(const G4int AParent, const G4int ZParent) const {
      const G4int SParent = 0;
      const G4int ADaughter = AParent - theA;
      const G4int ZDaughter = ZParent - theZ;
      const G4int SDaughter = 0;

      // Note the minus sign here
      G4double theQValue;
      if(isCluster())
        theQValue = -ParticleTable::getTableQValue(theA, theZ, theS, ADaughter, ZDaughter, SDaughter);
      else {
        const G4double massTableParent = ParticleTable::getTableMass(AParent,ZParent,SParent);
        const G4double massTableDaughter = ParticleTable::getTableMass(ADaughter,ZDaughter,SDaughter);
        const G4double massTableParticle = getTableMass();
        theQValue = massTableParent - massTableDaughter - massTableParticle;
      }

      const G4double massINCLParent = ParticleTable::getINCLMass(AParent,ZParent,SParent);
      const G4double massINCLDaughter = ParticleTable::getINCLMass(ADaughter,ZDaughter,SDaughter);
      const G4double massINCLParticle = getINCLMass();

      // The rhs corresponds to the INCL Q-value
      return theQValue - (massINCLParent-massINCLDaughter-massINCLParticle);
    }

    G4double getTransferQValueCorrection(const G4int AFrom, const G4int ZFrom, const G4int ATo, const G4int ZTo) const {
      const G4int SFrom = 0;
      const G4int STo = 0;
      const G4int AFromDaughter = AFrom - theA;
      const G4int ZFromDaughter = ZFrom - theZ;
      const G4int SFromDaughter = 0;
      const G4int AToDaughter = ATo + theA;
      const G4int ZToDaughter = ZTo + theZ;
      const G4int SToDaughter = 0;
      const G4double theQValue = ParticleTable::getTableQValue(AToDaughter,ZToDaughter,SToDaughter,AFromDaughter,ZFromDaughter,SFromDaughter,AFrom,ZFrom,SFrom);

      const G4double massINCLTo = ParticleTable::getINCLMass(ATo,ZTo,STo);
      const G4double massINCLToDaughter = ParticleTable::getINCLMass(AToDaughter,ZToDaughter,SToDaughter);
      /* Note that here we have to use the table mass in the INCL Q-value. We
       * cannot use theMass, because at this stage the particle is probably
       * still off-shell; and we cannot use getINCLMass(), because it leads to
       * violations of global energy conservation.
       */
      const G4double massINCLParticle = getTableMass();

      // The rhs corresponds to the INCL Q-value for particle absorption
      return theQValue - (massINCLToDaughter-massINCLTo-massINCLParticle);
    }

    G4double getEmissionQValueCorrection(const G4int AParent, const G4int ZParent, const G4int SParent) const {
      const G4int ADaughter = AParent - theA;
      const G4int ZDaughter = ZParent - theZ;
      const G4int SDaughter = SParent - theS;

      // Note the minus sign here
      G4double theQValue;
      if(isCluster())
        theQValue = -ParticleTable::getTableQValue(theA, theZ, theS, ADaughter, ZDaughter, SDaughter);
      else {
        const G4double massTableParent = ParticleTable::getTableMass(AParent,ZParent,SParent);
        const G4double massTableDaughter = ParticleTable::getTableMass(ADaughter,ZDaughter,SDaughter);
        const G4double massTableParticle = getTableMass();
        theQValue = massTableParent - massTableDaughter - massTableParticle;
      }

      const G4double massINCLParent = ParticleTable::getINCLMass(AParent,ZParent,SParent);
      const G4double massINCLDaughter = ParticleTable::getINCLMass(ADaughter,ZDaughter,SDaughter);
      const G4double massINCLParticle = getINCLMass();

      // The rhs corresponds to the INCL Q-value
      return theQValue - (massINCLParent-massINCLDaughter-massINCLParticle);
    }

    G4double getTransferQValueCorrection(const G4int AFrom, const G4int ZFrom, const G4int SFrom, const G4int ATo, const G4int ZTo , const G4int STo) const {
      const G4int AFromDaughter = AFrom - theA;
      const G4int ZFromDaughter = ZFrom - theZ;
      const G4int SFromDaughter = SFrom - theS;
      const G4int AToDaughter = ATo + theA;
      const G4int ZToDaughter = ZTo + theZ;
      const G4int SToDaughter = STo + theS;
      const G4double theQValue = ParticleTable::getTableQValue(AToDaughter,ZToDaughter,SFromDaughter,AFromDaughter,ZFromDaughter,SToDaughter,AFrom,ZFrom,SFrom);

      const G4double massINCLTo = ParticleTable::getINCLMass(ATo,ZTo,STo);
      const G4double massINCLToDaughter = ParticleTable::getINCLMass(AToDaughter,ZToDaughter,SToDaughter);
      /* Note that here we have to use the table mass in the INCL Q-value. We
       * cannot use theMass, because at this stage the particle is probably
       * still off-shell; and we cannot use getINCLMass(), because it leads to
       * violations of global energy conservation.
       */
      const G4double massINCLParticle = getTableMass();

      // The rhs corresponds to the INCL Q-value for particle absorption
      return theQValue - (massINCLToDaughter-massINCLTo-massINCLParticle);
    }



    G4double getInvariantMass() const {
      const G4double mass = std::pow(theEnergy, 2) - theMomentum.dot(theMomentum);
      if(mass < 0.0) {
        INCL_ERROR("E*E - p*p is negative." << '\n');
        return 0.0;
      } else {
        return std::sqrt(mass);
      }
    };

    inline G4double getKineticEnergy() const { return theEnergy - theMass; }

    inline G4double getPotentialEnergy() const { return thePotentialEnergy; }

    inline void setPotentialEnergy(G4double v) { thePotentialEnergy = v; }

    G4double getEnergy() const
    {
      return theEnergy;
    };

    void setMass(G4double mass)
    {
      this->theMass = mass;
    }

    void setEnergy(G4double energy)
    {
      this->theEnergy = energy;
    };

    const G4INCL::ThreeVector &getMomentum() const
    {
      return theMomentum;
    };

    virtual G4INCL::ThreeVector getAngularMomentum() const
    {
      return thePosition.vector(theMomentum);
    };

    virtual void setMomentum(const G4INCL::ThreeVector &momentum)
    {
      this->theMomentum = momentum;
    };

    const G4INCL::ThreeVector &getPosition() const
    {
      return thePosition;
    };

    virtual void setPosition(const G4INCL::ThreeVector &position)
    {
      this->thePosition = position;
    };

    G4double getHelicity() { return theHelicity; };
    void setHelicity(G4double h) { theHelicity = h; };

    void propagate(G4double step) {
      thePosition += ((*thePropagationMomentum)*(step/(*thePropagationEnergy)));
    };

    G4int getNumberOfCollisions() const { return nCollisions; }

    void setNumberOfCollisions(G4int n) { nCollisions = n; }

    void incrementNumberOfCollisions() { nCollisions++; }

    G4int getNumberOfDecays() const { return nDecays; }

    void setNumberOfDecays(G4int n) { nDecays = n; }

    void incrementNumberOfDecays() { nDecays++; }

    void setOutOfWell() { outOfWell = true; }

    G4bool isOutOfWell() const { return outOfWell; }

    void setEmissionTime(G4double t) { emissionTime = t; }
    G4double getEmissionTime() { return emissionTime; };

    ThreeVector getTransversePosition() const {
      return thePosition - getLongitudinalPosition();
    }

    ThreeVector getLongitudinalPosition() const {
      return *thePropagationMomentum * (thePosition.dot(*thePropagationMomentum)/thePropagationMomentum->mag2());
    }

    const ThreeVector &adjustMomentumFromEnergy();

    G4double adjustEnergyFromMomentum();

    G4bool isCluster() const {
      return (theType == Composite);
    }

    void setFrozenMomentum(const ThreeVector &momentum) { theFrozenMomentum = momentum; }

    void setFrozenEnergy(const G4double energy) { theFrozenEnergy = energy; }

    ThreeVector getFrozenMomentum() const { return theFrozenMomentum; }

    G4double getFrozenEnergy() const { return theFrozenEnergy; }

    ThreeVector getPropagationVelocity() const { return (*thePropagationMomentum)/(*thePropagationEnergy); }

    void freezePropagation() {
      thePropagationMomentum = &theFrozenMomentum;
      thePropagationEnergy = &theFrozenEnergy;
    }

    void thawPropagation() {
      thePropagationMomentum = &theMomentum;
      thePropagationEnergy = &theEnergy;
    }

    virtual void rotatePositionAndMomentum(const G4double angle, const ThreeVector &axis) {
      rotatePosition(angle, axis);
      rotateMomentum(angle, axis);
    }

    virtual void rotatePosition(const G4double angle, const ThreeVector &axis) {
      thePosition.rotate(angle, axis);
    }

    virtual void rotateMomentum(const G4double angle, const ThreeVector &axis) {
      theMomentum.rotate(angle, axis);
      theFrozenMomentum.rotate(angle, axis);
    }

    std::string print() const {
      std::stringstream ss;
      ss << "Particle (ID = " << ID << ") type = ";
      ss << ParticleTable::getName(theType);
      ss << '\n'
        << "   energy = " << theEnergy << '\n'
        << "   momentum = "
        << theMomentum.print()
        << '\n'
        << "   position = "
        << thePosition.print()
        << '\n';
      return ss.str();
    };

    std::string dump() const {
      std::stringstream ss;
      ss << "(particle " << ID << " ";
      ss << ParticleTable::getName(theType);
      ss << '\n'
        << thePosition.dump()
        << '\n'
        << theMomentum.dump()
        << '\n'
        << theEnergy << ")" << '\n';
      return ss.str();
    };

    long getID() const { return ID; };

    ParticleList const *getParticles() const {
      INCL_WARN("Particle::getParticles() method was called on a Particle object" << '\n');
      return 0;
    }

    G4double getReflectionMomentum() const {
      if(rpCorrelated)
        return theMomentum.mag();
      else
        return uncorrelatedMomentum;
    }

    void setUncorrelatedMomentum(const G4double p) { uncorrelatedMomentum = p; }

    void rpCorrelate() { rpCorrelated = true; }

    void rpDecorrelate() { rpCorrelated = false; }

    G4double getCosRPAngle() const {
      const G4double norm = thePosition.mag2()*thePropagationMomentum->mag2();
      if(norm>0.)
        return thePosition.dot(*thePropagationMomentum) / std::sqrt(norm);
      else
        return 1.;
    }

    static G4double getTotalBias();
    static void setINCLBiasVector(std::vector<G4double> NewVector);
    static void FillINCLBiasVector(G4double newBias);
    static G4double getBiasFromVector(std::vector<G4int> VectorBias);

    static std::vector<G4int> MergeVectorBias(Particle const * const p1, Particle const * const p2);
    static std::vector<G4int> MergeVectorBias(std::vector<G4int> p1, Particle const * const p2);

    G4double getParticleBias() const { return theParticleBias; };

    void setParticleBias(G4double ParticleBias) { this->theParticleBias = ParticleBias; }

    std::vector<G4int> getBiasCollisionVector() const { return theBiasCollisionVector; }

    void setBiasCollisionVector(std::vector<G4int> BiasCollisionVector) {
    this->theBiasCollisionVector = BiasCollisionVector;
    this->setParticleBias(Particle::getBiasFromVector(BiasCollisionVector));
    }
    
    G4int getNumberOfKaon() const { return theNKaon; };
    void setNumberOfKaon(const G4int NK) { theNKaon = NK; }
  
  public:
#ifdef INCLXX_IN_GEANT4_MODE
      static std::vector<G4double> INCLBiasVector;
      //static G4VectorCache<G4double> INCLBiasVector;
#else
      static G4ThreadLocal std::vector<G4double> INCLBiasVector;
      //static G4VectorCache<G4double> INCLBiasVector;
#endif
    static G4ThreadLocal G4int nextBiasedCollisionID;
    
  protected:
    G4int theZ, theA, theS;
    ParticipantType theParticipantType;
    G4INCL::ParticleType theType;
    G4double theEnergy;
    G4double *thePropagationEnergy;
    G4double theFrozenEnergy;
    G4INCL::ThreeVector theMomentum;
    G4INCL::ThreeVector *thePropagationMomentum;
    G4INCL::ThreeVector theFrozenMomentum;
    G4INCL::ThreeVector thePosition;
    G4int nCollisions;
    G4int nDecays;
    G4double thePotentialEnergy;
    long ID;

    G4bool rpCorrelated;
    G4double uncorrelatedMomentum;
    
    G4double theParticleBias;
    G4int theNKaon;

  private:
    G4double theHelicity;
    G4double emissionTime;
    G4bool outOfWell;
    
    std::vector<G4int> theBiasCollisionVector;

    G4double theMass;
    static G4ThreadLocal long nextID;

    INCL_DECLARE_ALLOCATION_POOL(Particle)
  };
}

#endif /* PARTICLE_HH_ */