//
// ********************************************************************
// * 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.0_rc3
//
#define INCLXX_IN_GEANT4_MODE 1

#include "globals.hh"

/*
 * Particle.hh
 *
 *  Created on: Jun 5, 2009
 *      Author: Pekka Kaitaniemi
 */

#ifndef PARTICLE_HH_
#define PARTICLE_HH_

#include "G4INCLThreeVector.hh"
#include "G4INCLParticleTable.hh"
#include "G4INCLParticleType.hh"
#include "G4INCLLogger.hh"
#include <list>
#include <sstream>
#include <string>

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();

    /**
     * Get the particle type.
     * @see G4INCL::ParticleType
     */
    G4INCL::ParticleType getType() const {
      return 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);
          break;
        case UnknownParticle:
          ERROR("Trying to set particle type to Unknown!" << std::endl);
          break;
      }

      if( !isResonance() && t!=Composite )
        theMass = ParticleTable::getMass(theType);
    }

    /**
     * Is this a nucleon?
     */
    G4bool isNucleon() const {
      if(theType == G4INCL::Proton || theType == G4INCL::Neutron)
	return true;
      else
	return false;
    };

    G4bool isParticipant() const {
      return participant;
    }

    void makeParticipant() {
      participant = true;
    }

    void makeSpectator() {
      participant = false;
    }

    /** \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 &boostVector) {
      const G4double beta2 = boostVector.mag2();
      const G4double gamma = 1.0 / std::sqrt(1.0 - beta2);
      const G4double bp = theMomentum.dot(boostVector);
      const G4double alpha = (gamma*gamma)/(1.0 + gamma);

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

    /** \brief Get the cached particle mass. */
    inline G4double getMass() const { return theMass; }

    /** \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 */
    G4INCL::ThreeVector getAngularMomentum() const
    {
      return thePosition.vector(theMomentum);
    };

    /**
     * Set the momentum vector.
     */
    void setMomentum(const G4INCL::ThreeVector &momentum)
    {
      this->theMomentum = momentum;
    };

    /**
     * Set the position vector.
     */
    const G4INCL::ThreeVector &getPosition() const
    {
      return thePosition;
    };

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

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

    void propagate(G4double step) {
      thePosition += (theMomentum*(step/theEnergy));
    };

    /** \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() {
      return thePosition - theMomentum *
        (thePosition.dot(theMomentum)/theMomentum.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 {
      for(ParticleIter i=l.begin(); i!=l.end(); ++i)
        if((*i)->getID()==ID) return true;
      return false;
    }

    G4bool isCluster() const {
      if(theType == Composite) return true;
      else return false;
    }

    std::string prG4int() const {
      std::stringstream ss;
      ss << "Particle (ID = " << ID << ") type = ";
      ss << ParticleTable::getName(theType);
      ss << std::endl
        << "   energy = " << theEnergy << std::endl
        << "   momentum = "
        << theMomentum.prG4int()
        << std::endl
        << "   position = "
        << thePosition.prG4int()
        << 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 poG4inter
     */
    ParticleList const *getParticles() const {
      WARN("Particle::getParticles() method was called on a Particle object" << std::endl);
      return 0;
    }

  protected:
    G4int theZ, theA;
    G4bool participant;
    G4INCL::ParticleType theType;
    G4double theEnergy;
    G4INCL::ThreeVector theMomentum;
    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_ */