Sphericity.hh

// -*- C++ -*-
#ifndef RIVET_Sphericity_HH
#define RIVET_Sphericity_HH
 
#include "Rivet/Event.hh"
#include "Rivet/Jet.fhh"
#include "Rivet/Projection.hh"
#include "Rivet/Projections/AxesDefinition.hh"
#include "Rivet/Projections/FinalState.hh"
 
namespace Rivet {
 

  class Sphericity : public AxesDefinition {
  public:
 
    using AxesDefinition::operator=;


    Sphericity(double rparam = 2.0)
        : _regparam(rparam) { }
 
    Sphericity(const FinalState& fsp, double rparam = 2.0);

    RIVET_DEFAULT_PROJ_CLONE(Sphericity);


    using Projection::operator=;
 
 
  protected:

    void project(const Event& e);

    CmpState compare(const Projection& p) const;
 
 
  public:

    void clear();

    double sphericity() const {
      return 3.0 / 2.0 * (lambda2() + lambda3());
    }
    double transSphericity() const {
      return 2.0 * lambda2() / (lambda1() + lambda2());
    }

    double planarity() const {
      return 2 * (sphericity() - 2 * aplanarity()) / 3.0;
    }

    double aplanarity() const {
      return 3 / 2.0 * lambda3();
    }

 

    const Vector3& sphericityAxis() const {
      return _sphAxes[0];
    }
    const Vector3& sphericityMajorAxis() const {
      return _sphAxes[1];
    }

    const Vector3& sphericityMinorAxis() const {
      return _sphAxes[2];
    }

 

    const Vector3& axis1() const {
      return sphericityAxis();
    }
    const Vector3& axis2() const {
      return sphericityMajorAxis();
    }
    const Vector3& axis3() const {
      return sphericityMinorAxis();
    }

 

    double lambda1() const {
      return _lambdas[0];
    }
    double lambda2() const {
      return _lambdas[1];
    }
    double lambda3() const {
      return _lambdas[2];
    }
 
    Vector3 mkEigenVector(Matrix3 A, const double& lambda);


    void calc(const FinalState& fs);

    void calc(const Particles& particles);

    void calc(const Jets& jets);

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

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

 
 
  protected:

    vector<double> _lambdas;

    vector<Vector3> _sphAxes;

    const double _regparam;
  };
 
 
}
 
#endif