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

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

namespace Rivet {


  class Spherocity : public AxesDefinition {
  public:

    // Default Constructor
    Spherocity() {}

    Spherocity(const FinalState& fsp) {
      setName("Spherocity");
      declare(fsp, "FS");
    }

    DEFAULT_RIVET_PROJ_CLONE(Spherocity);


  protected:

    void project(const Event& e) {
      const vector<Particle> ps
        = applyProjection<FinalState>(e, "FS").particles();
      calc(ps);
    }


    CmpState compare(const Projection& p) const {
      return mkNamedPCmp(p, "FS");
    }


  public:



    double spherocity() const { return _spherocities[0]; }




    const Vector3& spherocityAxis() const { return _spherocityAxes[0]; }
    const Vector3& spherocityMajorAxis() const { return _spherocityAxes[1]; }
    const Vector3& spherocityMinorAxis() const { return _spherocityAxes[2]; }



    const Vector3& axis1() const { return spherocityAxis(); }
    const Vector3& axis2() const { return spherocityMajorAxis(); }
    const Vector3& axis3() const { return spherocityMinorAxis(); }


  public:



    void calc(const FinalState& fs);

    void calc(const vector<Particle>& fsparticles);

    void calc(const vector<FourMomentum>& fsmomenta);

    void calc(const vector<Vector3>& threeMomenta);



  private:

    vector<double> _spherocities;

    vector<Vector3> _spherocityAxes;


  private:

    void _calcSpherocity(const vector<Vector3>& fsmomenta);

  };

}

#endif