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
//
// INCL++ revision: v5.1.8
//
#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 <list>
#include <sstream>
#include <string>
#include <algorithm>
 
namespace G4INCL {
 
  class Particle;
 
  typedef std::list<G4INCL::Particle*> ParticleList;
  typedef std::list<G4INCL::Particle*>::const_iterator ParticleIter;
 
  class Particle {
  public:
    Particle();
    Particle(ParticleType t, G4double energy, ThreeVector momentum, ThreeVector position);
    Particle(ParticleType t, ThreeVector momentum, ThreeVector position);
    virtual ~Particle() {}
 
    /** \brief Copy constructor
     *
     * Does not copy the particle ID.
     */
    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),
      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:
    /// \brief Helper method for the assignment operator
    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);
    }
 
  public:
 
    /** \brief Assignment operator
     *
     * Does not copy the particle ID.
     */
    Particle &operator=(const Particle &rhs) {
      Particle temporaryParticle(rhs);
      swap(temporaryParticle);
      return *this;
    }
 
    /**
     * Get the particle type.
     * @see G4INCL::ParticleType
     */
    G4INCL::ParticleType getType() const {
      return theType;
    };
 
    /// \brief Get the particle species
    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:
         // 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;
          ERROR("Trying to set particle type to Unknown!" << std::endl);
          break;
      }
 
      if( !isResonance() && t!=Composite )
        setINCLMass();
    }
 
    /**
     * Is this a nucleon?
     */
    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;
    }
 
    /** \brief Is this a pion? */
    G4bool isPion() const { return (theType == PiPlus || theType == PiZero || theType == PiMinus); }
 
    /** \brief Is it a resonance? */
    inline G4bool isResonance() const { return isDelta(); }
 
    /** \brief Is it a Delta? */
    inline G4bool isDelta() const {
      return (theType==DeltaPlusPlus || theType==DeltaPlus ||
          theType==DeltaZero || theType==DeltaMinus);
    }
 
    /** \brief Returns the baryon number. */
    G4int getA() const { return theA; }
 
    /** \brief Returns the charge number. */
    G4int getZ() const { return theZ; }
 
    G4double getBeta() const {
      const G4double P = theMomentum.mag();
      return P/theEnergy;
    }
 
    /**
     * Returns a three vector we can give to the boost() -method.
     *
     * In order to go to the particle rest frame you need to multiply
     * the boost vector by -1.0.
     */
    ThreeVector boostVector() const {
      return theMomentum / theEnergy;
    }
 
    /**
     * Boost the particle using a boost vector.
     *
     * Example (go to the particle rest frame):
     * particle->boost(particle->boostVector());
     */
    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);
    }
 
    /** \brief Lorentz-contract the particle position around some center
     *
     * Apply Lorentz contraction to the position component along the
     * direction of the boost vector.
     *
     * \param aBoostVector the boost vector (velocity) [c]
     * \param refPos the reference position
     */
    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;
    }
 
    /** \brief Get the cached particle mass. */
    inline G4double getMass() const { return theMass; }
 
    /** \brief Get the INCL particle mass. */
    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:
          ERROR("Particle::getINCLMass: Unknown particle type." << std::endl);
          return 0.0;
          break;
      }
    }
 
    /** \brief Get the tabulated particle mass. */
    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:
          ERROR("Particle::getTableMass: Unknown particle type." << std::endl);
          return 0.0;
          break;
      }
    }
 
    /** \brief Get the real particle mass. */
    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:
          ERROR("Particle::getRealMass: Unknown particle type." << std::endl);
          return 0.0;
          break;
      }
    }
 
    /// \brief Set the mass of the Particle to its real mass
    void setRealMass() { setMass(getRealMass()); }
 
    /// \brief Set the mass of the Particle to its table mass
    void setTableMass() { setMass(getTableMass()); }
 
    /// \brief Set the mass of the Particle to its table mass
    void setINCLMass() { setMass(getINCLMass()); }
 
    /**\brief Computes correction on the emission Q-value
     *
     * Computes the correction that must be applied to INCL particles in
     * order to obtain the correct Q-value for particle emission from a given
     * nucleus. For absorption, the correction is obviously equal to minus
     * the value returned by this function.
     *
     * \param AParent the mass number of the emitting nucleus
     * \param ZParent the charge number of the emitting nucleus
     * \return the correction
     */
    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);
    }
 
    /**\brief Computes correction on the transfer Q-value
     *
     * Computes the correction that must be applied to INCL particles in
     * order to obtain the correct Q-value for particle transfer from a given
     * nucleus to another.
     *
     * Assumes that the receving nucleus is INCL's target nucleus, with the
     * INCL separation energy.
     *
     * \param AFrom the mass number of the donating nucleus
     * \param ZFrom the charge number of the donating nucleus
     * \param ATo the mass number of the receiving nucleus
     * \param ZTo the charge number of the receiving nucleus
     * \return the correction
     */
    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);
    }
 
    /** \brief Get the the particle invariant mass.
     *
     * Uses the relativistic invariant
     * \f[ m = \sqrt{E^2 - {\vec p}^2}\f]
     **/
    G4double getInvariantMass() const {
      const G4double mass = std::pow(theEnergy, 2) - theMomentum.dot(theMomentum);
      if(mass < 0.0) {
        ERROR("E*E - p*p is negative." << std::endl);
        return 0.0;
      } else {
        return std::sqrt(mass);
      }
    };
 
    /// \brief Get the particle kinetic energy.
    inline G4double getKineticEnergy() const { return theEnergy - theMass; }
 
    /// \brief Get the particle potential energy.
    inline G4double getPotentialEnergy() const { return thePotentialEnergy; }
 
    /// \brief Set the particle potential energy.
    inline void setPotentialEnergy(G4double v) { thePotentialEnergy = v; }
 
    /**
     * Get the energy of the particle in MeV.
     */
    G4double getEnergy() const
    {
      return theEnergy;
    };
 
    /**
     * Set the mass of the particle in MeV/c^2.
     */
    void setMass(G4double mass)
    {
      this->theMass = mass;
    }
 
    /**
     * Set the energy of the particle in MeV.
     */
    void setEnergy(G4double energy)
    {
      this->theEnergy = energy;
    };
 
    /**
     * Get the momentum vector.
     */
    const G4INCL::ThreeVector &getMomentum() const
    {
      return theMomentum;
    };
 
    /** Get the angular momentum w.r.t. the origin */
    virtual G4INCL::ThreeVector getAngularMomentum() const
    {
      return thePosition.vector(theMomentum);
    };
 
    /**
     * Set the momentum vector.
     */
    virtual void setMomentum(const G4INCL::ThreeVector &momentum)
    {
      this->theMomentum = momentum;
    };
 
    /**
     * Set the position vector.
     */
    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)));
    };
 
    /** \brief Return the number of collisions undergone by the particle. **/
    G4int getNumberOfCollisions() const { return nCollisions; }
 
    /** \brief Set the number of collisions undergone by the particle. **/
    void setNumberOfCollisions(G4int n) { nCollisions = n; }
 
    /** \brief Increment the number of collisions undergone by the particle. **/
    void incrementNumberOfCollisions() { nCollisions++; }
 
    /** \brief Return the number of decays undergone by the particle. **/
    G4int getNumberOfDecays() const { return nDecays; }
 
    /** \brief Set the number of decays undergone by the particle. **/
    void setNumberOfDecays(G4int n) { nDecays = n; }
 
    /** \brief Increment the number of decays undergone by the particle. **/
    void incrementNumberOfDecays() { nDecays++; }
 
    /** \brief Mark the particle as out of its potential well
     *
     * This flag is used to control pions created outside their potential well
     * in delta decay. The pion potential checks it and returns zero if it is
     * true (necessary in order to correctly enforce energy conservation). The
     * Nucleus::applyFinalState() method uses it to determine whether new
     * avatars should be generated for the particle.
     */
    void setOutOfWell() { outOfWell = true; }
 
    /// \brief Check if the particle is out of its potential well
    G4bool isOutOfWell() const { return outOfWell; }
 
    void setEmissionTime(G4double t) { emissionTime = t; }
    G4double getEmissionTime() { return emissionTime; };
 
    /** \brief Transverse component of the position w.r.t. the momentum. */
    ThreeVector getTransversePosition() const {
      return thePosition - getLongitudinalPosition();
    }
 
    /** \brief Longitudinal component of the position w.r.t. the momentum. */
    ThreeVector getLongitudinalPosition() const {
      return *thePropagationMomentum * (thePosition.dot(*thePropagationMomentum)/thePropagationMomentum->mag2());
    }
 
    /** \brief Rescale the momentum to match the total energy. */
    const ThreeVector &adjustMomentumFromEnergy();
 
    /** \brief Recompute the energy to match the momentum. */
    G4double adjustEnergyFromMomentum();
 
    /** \brief Check if the particle belongs to a given list **/
    G4bool isInList(ParticleList const &l) const {
      return (std::find(l.begin(), l.end(), this)!=l.end());
    }
 
    G4bool isCluster() const {
      return (theType == Composite);
    }
 
    /// \brief Set the frozen particle momentum
    void setFrozenMomentum(const ThreeVector &momentum) { theFrozenMomentum = momentum; }
 
    /// \brief Set the frozen particle momentum
    void setFrozenEnergy(const G4double energy) { theFrozenEnergy = energy; }
 
    /// \brief Get the frozen particle momentum
    ThreeVector getFrozenMomentum() const { return theFrozenMomentum; }
 
    /// \brief Get the frozen particle momentum
    G4double getFrozenEnergy() const { return theFrozenEnergy; }
 
    /// \brief Get the propagation velocity of the particle
    ThreeVector getPropagationVelocity() const { return (*thePropagationMomentum)/(*thePropagationEnergy); }
 
    /** \brief Freeze particle propagation
     *
     * Make the particle use theFrozenMomentum and theFrozenEnergy for
     * propagation. The normal state can be restored by calling the
     * thawPropagation() method.
     */
    void freezePropagation() {
      thePropagationMomentum = &theFrozenMomentum;
      thePropagationEnergy = &theFrozenEnergy;
    }
 
    /** \brief Unfreeze particle propagation
     *
     * Make the particle use theMomentum and theEnergy for propagation. Call
     * this method to restore the normal propagation if the
     * freezePropagation() method has been called.
     */
    void thawPropagation() {
      thePropagationMomentum = &theMomentum;
      thePropagationEnergy = &theEnergy;
    }
 
    /** \brief Rotate the particle position and momentum
     *
     * \param angle the rotation angle
     * \param axis a unit vector representing the rotation axis
     */
    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; };
 
    /**
     * Return a NULL pointer
     */
    ParticleList const *getParticles() const {
      WARN("Particle::getParticles() method was called on a Particle object" << std::endl);
      return 0;
    }
 
  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;
 
  private:
    G4double theHelicity;
    G4double emissionTime;
    G4bool outOfWell;
 
    G4double theMass;
    static long nextID;
 
  };
}
 
#endif /* PARTICLE_HH_ */