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

#include "Rivet/Projection.hh"
#include "Rivet/Projections/AxesDefinition.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Event.hh"

namespace Rivet {
  
  //Projection to find the broadening for a given event along the Thrust axis.
  
  class Broadening : public AxesDefinition {
  public:
    
    /// Constructor.
    Broadening(const FinalState& fsp)
      : _calculatedBroadening(false)
    {
      setName("Broadening");
      addProjection(fsp, "FS");
    }
    
    /// Clone on the heap.
    virtual const Projection* clone() const {
      return new Broadening(*this);
    }
    
  protected:
    
    /// Perform the projection on the Event
    void project(const Event& e) {
      const vector<Particle> ps
	= applyProjection<FinalState>(e, "FS").particles();
      calc(ps);
    }
    
    /// Compare projections
    int compare(const Projection& p) const {
      return mkNamedPCmp(p, "FS");
     }
    
    
  public:
    
    ///@{ Broadening scalar accessors
    /// The broadening scalar, \f$ T \f$, (maximum broadening).
    const double broadening() const { return _broadenings[0]; }
    /// The broadening major scalar, \f$ M \f$, (broadening along broadening major axis).
    const double broadeningMajor() const { return _broadenings[1]; }
    /// The broadening minor scalar, \f$ m \f$, (broadening along broadening minor axis).
    const double broadeningMinor() const { return _broadenings[2]; }
    /// The oblateness, \f$ O = M - m \f$ .
    const double oblateness() const { return _broadenings[1] - _broadenings[2]; }
    ///@}
 
    ///@{ Broadening axis accessors
    /// The broadening axis.
    const Vector3& broadeningAxis() const { return _broadeningAxes[0]; }
    /// The broadening major axis (axis of max broadening perpendicular to broadening axis).
    const Vector3& broadeningMajorAxis() const { return _broadeningAxes[1]; }
    /// The broadening minor axis (axis perpendicular to broadening and broadening major).
    const Vector3& broadeningMinorAxis() const { return _broadeningAxes[2]; }
    ///@}
    
    ///@{ AxesDefinition axis accessors.
    const Vector3& axis1() const { return broadeningAxis(); }
    const Vector3& axis2() const { return broadeningMajorAxis(); }
    const Vector3& axis3() const { return broadeningMinorAxis(); }
    ///@}
    
    
  public:
    
    /// @name Direct methods
    /// Ways to do the calculation directly, without engaging the caching system
    //@{
 
    /// Manually calculate the broadening, without engaging the caching system
    void calc(const FinalState& fs);
    
    /// Manually calculate the broadening, without engaging the caching system
    void calc(const vector<Particle>& fsparticles);
    
    /// Manually calculate the broadening, without engaging the caching system
    void calc(const vector<FourMomentum>& fsmomenta);
    
    /// Manually calculate the broadening, without engaging the caching system
    void calc(const vector<Vector3>& threeMomenta);
    
    //@}
    
    
  private:
    
    /// The broadening scalars.
    vector<double> _broadenings;
    
    /// The broadening axes.
    vector<Vector3> _broadeningAxes;
    
    /// Caching flag to avoid costly recalculations.
    bool _calculatedBroadening;
    
  private:
    
    /// Explicitly calculate the broadening values.
    void _calcBroadening(const vector<Vector3>& fsmomenta);
    
  };
  
 }
 
#endif