// -*- C++ -*-
#ifndef RIVET_ThrustGamma_HH
#define RIVET_ThrustGamma_HH

#include "Rivet/Projection.hh"
#include "Rivet/Projections/AxesDefinition.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Event.hh"
#include "Rivet/Particle.hh"
#include "Rivet/Projections/DISKinematics.hh"
#include "Rivet/Projections/DISLepton.hh"
#include "Rivet/Projections/Beam.hh"
#include "Rivet/Rivet.hh"
#include "ThrustGamma.hh"
#include "Rivet/Tools/Logging.hh"
#include "Rivet/Cmp.hh"
#include "Rivet/Projections/DISKinematics.hh"
#include "Rivet/Math/Constants.hh"
#include "Rivet/Tools/ParticleIdUtils.hh"
#include "math.h"
#include "FinalStateBreitCurrent.hh"



namespace Rivet {
  
  //Projection to find the Thrust along an axis defined by the exchange photon. 

  class ThrustGamma : public AxesDefinition {
  public:
    
    /// Constructor.
    ThrustGamma(const DISKinematics& kinematicsp, const FinalState& fs);
    
    /// Clone on the heap.
    virtual const Projection* clone() const {
      return new ThrustGamma(*this);
    }
    
  protected:
    
    /// Perform the projection on the Event
    void project(const Event& e);
    
    /// Compare projections
    int compare(const Projection& p)const;
    
    
  public:
    
    ///@{ ThrustGamma scalar accessors
    /// The thrustgamma scalar, \f$ T \f$, (maximum thrustgamma).
    const double thrustgamma() const { return _thrustgammas[0]; }
    /// The thrustgamma major scalar, \f$ M \f$, (thrustgamma along thrustgamma major axis).
    const double thrustgammaMajor() const { return _thrustgammas[1]; }
    /// The thrustgamma minor scalar, \f$ m \f$, (thrustgamma along thrustgamma minor axis).
    const double thrustgammaMinor() const { return _thrustgammas[2]; }
    /// The oblateness, \f$ O = M - m \f$ .
    const double oblateness() const { return _thrustgammas[1] - _thrustgammas[2]; }
    ///@}
 
    ///@{ ThrustGamma axis accessors
    /// The thrustgamma axis.
    const Vector3& thrustgammaAxis() const { return _thrustgammaAxes[0]; }
    /// The thrustgamma major axis (axis of max thrustgamma perpendicular to thrustgamma axis).
    const Vector3& thrustgammaMajorAxis() const { return _thrustgammaAxes[1]; }
    /// The thrustgamma minor axis (axis perpendicular to thrustgamma and thrustgamma major).
    const Vector3& thrustgammaMinorAxis() const { return _thrustgammaAxes[2]; }
    ///@}
    
    ///@{ AxesDefinition axis accessors.
    const Vector3& axis1() const { return thrustgammaAxis(); }
    const Vector3& axis2() const { return thrustgammaMajorAxis(); }
    const Vector3& axis3() const { return thrustgammaMinorAxis(); }
    ///@}
    
    
  public:
    
    /// @name Direct methods
    /// Ways to do the calculation directly, without engaging the caching system
    //@{
 
    
    /// Manually calculate the thrustgamma, without engaging the caching system
    void calcg(const vector<Particle>& fsparticles, FourMomentum pgamma);
    
    /// Manually calculate the thrustgamma, without engaging the caching system
    void calcg(const vector<FourMomentum>& fsmomenta, FourMomentum pgamma);
    
    /// Manually calculate the thrustgamma, without engaging the caching system
    void calcg(const vector<Vector3>& threeMomenta, FourMomentum pgamma);
    
    //@}
    
    
  private:
    
    /// The thrustgamma scalars.
    vector<double> _thrustgammas;
    
    /// The thrustgamma axes.
    vector<Vector3> _thrustgammaAxes;
    
    /// Caching flag to avoid costly recalculations.
    bool _calbgulatedThrustGamma;
    
  private:
    
    /// Explicitly calculate the thrustgamma values.
    void _calcgThrustGamma(const vector<Vector3>& fsmomenta, FourMomentum pgamma);
    
  };
  
 }
 
#endif