// -*- C++ -*-
#ifndef RIVET_BroadeningGamma_HH
#define RIVET_BroadeningGamma_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/Math/Constants.hh"
#include "Rivet/Tools/ParticleIdUtils.hh"
#include "math.h"
#include "FinalStateBreit.hh"



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

  class BroadeningGamma : public AxesDefinition {
  public:
    
    /// Constructor.
    BroadeningGamma(const DISKinematics& kinematicsp, const FinalState& fs){

    setName("BroadeningGamma");
    addProjection(kinematicsp, "DK");
    
    const FinalStateBreit fsb = FinalStateBreit(kinematicsp, 1, 0.25);
    addProjection (fsb, "FSB");
      
  }
    /// Clone on the heap.
    virtual const Projection* clone() const {
      return new BroadeningGamma(*this);
    }
    
  protected:
    
    /// Perform the projection on the Event
    void project(const Event& e){
      
      const vector<Particle> ps = applyProjection<FinalStateBreit>(e, "FSB").fsbparticles();
      const FourMomentum pgam = applyProjection<FinalStateBreit>(e, "FSB").pgamma();

      calcbg(ps, pgam);

    }
    
    /// Compare projections
    int compare(const Projection& p)const;
    
    
  public:
    
    ///@{ BroadeningGamma scalar accessors
    /// The broadeinggamma scalar, \f$ T \f$, (maximum broadeninggamma).
    const double broadeninggamma() const { return _broadeninggammas[0]; }

    /// The Broadeninggamma major scalar, \f$ M \f$, (broadeninggamma along briadeninggamma major axis).
    const double broadeninggammaMajor() const { return _broadeninggammas[1]; }

    /// The broadeninggamma minor scalar, \f$ m \f$, (broadeninggamma along broadeninggamma minor axis).
    const double broadeninggammaMinor() const { return _broadeninggammas[2]; }

    /// The oblateness, \f$ O = M - m \f$ .
    const double oblateness() const { return _broadeninggammas[1] - _broadeninggammas[2]; }
    ///@}
 
    ///@{ BroadeningGamma axis accessors
    /// The Broadeninggamma axis.
    const Vector3& broadeninggammaAxis() const { return _broadeninggammaAxes[0]; }

    /// The Broadeninggamma major axis (axis of max broadeninggamma perpendicular to broadeninggamma axis).
    const Vector3& broadeninggammaMajorAxis() const { return _broadeninggammaAxes[1]; }

    /// The broadeninggamma minor axis (axis perpendicular to broadeninggamma and broadeninggamma major).
    const Vector3& broadeninggammaMinorAxis() const { return _broadeninggammaAxes[2]; }
    ///@}
    
    ///@{ AxesDefinition axis accessors.
    const Vector3& axis1() const { return broadeninggammaAxis(); }
    const Vector3& axis2() const { return broadeninggammaMajorAxis(); }
    const Vector3& axis3() const { return broadeninggammaMinorAxis(); }
    ///@}
    
    
  public:
    
    /// @name Direct methods
    /// Ways to do the calculation directly, without engaging the caching system
    //@{
 
    
    /// Manually calculate the broadeninggamma, without engaging the caching system
    void calcbg(const vector<Particle>& fsparticles, FourMomentum pgamma);
    
    /// Manually calculate the broadeninggamma, without engaging the caching system
    void calcbg(const vector<FourMomentum>& fsmomenta, FourMomentum pgamma);
    
    /// Manually calculate the broadeninggamma, without engaging the caching system
    void calcbg(const vector<Vector3>& threeMomenta, FourMomentum pgamma);
    
    //@}
    
    
  private:
     
    /// The broadeninggamma scalars.
    vector<double> _broadeninggammas;
    
    /// The broadeninggamma axes.
    vector<Vector3> _broadeninggammaAxes;
    
    /// Caching flag to avoid costly recalculations.
    bool _calbgulatedbroadeningGamma;
    
  private:
    
    /// Explicitly calculate the broadeninggamma values.
    void _calcBroadeningGamma(const vector<Vector3>& fsmomenta, FourMomentum pgamma);
    
  };
  
 }
 
#endif