G4INCLParticle.hh

Go to the documentation of this file.

//
// ********************************************************************
// * 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
// Pekka Kaitaniemi, CEA and Helsinki Institute of Physics
// Davide Mancusi, CEA
// Alain Boudard, CEA
// Sylvie Leray, CEA
// Joseph Cugnon, University of Liege
//
#define INCLXX_IN_GEANT4_MODE 1

#include "globals.hh"

/*
 * Particle.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 <vector>
#include <sstream>
#include <string>
#include <algorithm>

namespace G4INCL {

  class Particle;

  template<class T>
    class UnorderedVector : private std::vector<T> {
      public:
        UnorderedVector() {}
        using std::vector<T>::push_back;
        using std::vector<T>::pop_back;
        using std::vector<T>::size;
        using std::vector<T>::begin;
        using std::vector<T>::end;
        using std::vector<T>::rbegin;
        using std::vector<T>::rend;
        using std::vector<T>::front;
        using std::vector<T>::back;
        using std::vector<T>::clear;
        using std::vector<T>::empty;
        using std::vector<T>::insert;
        using std::vector<T>::erase;
        using typename std::vector<T>::iterator;
        using typename std::vector<T>::reverse_iterator;
        using typename std::vector<T>::const_iterator;
        using typename std::vector<T>::const_reverse_iterator;
        void remove(const T &t) {
          const typename std::vector<T>::iterator removeMe = std::find(begin(), end(), t);
// assert(removeMe!=end());
          *removeMe = back();
          pop_back();
        }
    };

  typedef UnorderedVector<Particle*>   ParticleList;
  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),
      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),
      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++;
      }

  protected:
    void swap(Particle &rhs) {
      std::swap(theZ, rhs.theZ);
      std::swap(theA, rhs.theA);
      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);
    }

  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;
          break;
        case Proton:
        case DeltaPlus:
          theA = 1;
          theZ = 1;
          break;
        case Neutron:
        case DeltaZero:
          theA = 1;
          theZ = 0;
          break;
        case DeltaMinus:
          theA = 1;
          theZ = -1;
          break;
        case PiPlus:
          theA = 0;
          theZ = 1;
          break;
        case PiZero:
          theA = 0;
          theZ = 0;
          break;
        case PiMinus:
          theA = 0;
          theZ = -1;
          break;
        case Composite:
         // INCL_ERROR("Trying to set particle type to Composite! Construct a Cluster object instead" << std::endl);
          theA = 0;
          theZ = 0;
          break;
        case UnknownParticle:
          theA = 0;
          theZ = 0;
          INCL_ERROR("Trying to set particle type to Unknown!" << std::endl);
          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); }

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

    inline G4bool isDelta() const {
      return (theType==DeltaPlusPlus || theType==DeltaPlus ||
          theType==DeltaZero || theType==DeltaMinus);
    }

    G4int getA() const { return theA; }

    G4int getZ() const { return theZ; }

    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:
          return ParticleTable::getINCLMass(theType);
          break;

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

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

        default:
          INCL_ERROR("Particle::getINCLMass: Unknown particle type." << std::endl);
          return 0.0;
          break;
      }
    }

    inline virtual G4double getTableMass() const {
      switch(theType) {
        case Proton:
        case Neutron:
        case PiPlus:
        case PiMinus:
        case PiZero:
          return ParticleTable::getTableParticleMass(theType);
          break;

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

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

        default:
          INCL_ERROR("Particle::getTableMass: Unknown particle type." << std::endl);
          return 0.0;
          break;
      }
    }

    inline G4double getRealMass() const {
      switch(theType) {
        case Proton:
        case Neutron:
        case PiPlus:
        case PiMinus:
        case PiZero:
          return ParticleTable::getRealMass(theType);
          break;

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

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

        default:
          INCL_ERROR("Particle::getRealMass: Unknown particle type." << std::endl);
          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 ADaughter = AParent - theA;
      const G4int ZDaughter = ZParent - theZ;

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

      const G4double massINCLParent = ParticleTable::getINCLMass(AParent,ZParent);
      const G4double massINCLDaughter = ParticleTable::getINCLMass(ADaughter,ZDaughter);
      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 AFromDaughter = AFrom - theA;
      const G4int ZFromDaughter = ZFrom - theZ;
      const G4int AToDaughter = ATo + theA;
      const G4int ZToDaughter = ZTo + theZ;
      const G4double theQValue = ParticleTable::getTableQValue(AToDaughter,ZToDaughter,AFromDaughter,ZFromDaughter,AFrom,ZFrom);

      const G4double massINCLTo = ParticleTable::getINCLMass(ATo,ZTo);
      const G4double massINCLToDaughter = ParticleTable::getINCLMass(AToDaughter,ZToDaughter);
      /* 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." << std::endl);
        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 isInList(ParticleList const &l) const {
      return (std::find(l.begin(), l.end(), this)!=l.end());
    }

    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 rotate(const G4double angle, const ThreeVector &axis) {
      thePosition.rotate(angle, 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 << std::endl
        << "   energy = " << theEnergy << std::endl
        << "   momentum = "
        << theMomentum.print()
        << std::endl
        << "   position = "
        << thePosition.print()
        << std::endl;
      return ss.str();
    };

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

    long getID() const { return ID; };

    ParticleList const *getParticles() const {
      INCL_WARN("Particle::getParticles() method was called on a Particle object" << std::endl);
      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.;
    }

  protected:
    G4int theZ, theA;
    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;

  private:
    G4double theHelicity;
    G4double emissionTime;
    G4bool outOfWell;

    G4double theMass;
    static G4ThreadLocal long nextID;

  };
}

#endif /* PARTICLE_HH_ */