JetShape.hh

// -*- C++ -*-
#ifndef RIVET_JetShape_HH
#define RIVET_JetShape_HH
 
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/Projection.hh"
#include "Rivet/Projections/JetFinder.hh"
#include "Rivet/Particle.hh"
#include "Rivet/Event.hh"
#include "Rivet/Tools/Utils.hh"
 
namespace Rivet {
 

  class JetShape : public Projection {
  public:


    JetShape(const JetFinder& jetalg,
             double rmin, double rmax, size_t nbins,
             double ptmin=0, double ptmax=DBL_MAX,
             double absrapmin=-DBL_MAX, double absrapmax=-DBL_MAX,
             RapScheme rapscheme=RAPIDITY);

    JetShape(const JetFinder& jetalg, vector<double> binedges,
             double ptmin=0, double ptmax=DBL_MAX,
             double absrapmin=-DBL_MAX, double absrapmax=-DBL_MAX,
             RapScheme rapscheme=RAPIDITY);

    RIVET_DEFAULT_PROJ_CLONE(JetShape);


    using Projection::operator =;
 

    void clear();
 

    void calc(const Jets& jets);
 
 
  public:
 

    size_t numBins() const {
      return _binedges.size() - 1;
    }

    size_t numJets() const {
      return _diffjetshapes.size();
    }

    double rMin() const {
      return _binedges.front();
    }

    double rMax() const {
      return _binedges.back();
    }

    double ptMin() const {
      return _ptcuts.first;
    }

    double ptMax() const {
      return _ptcuts.second;
    }

    double rBinMin(size_t rbin) const {
      assert(inRange(rbin, 0u, numBins()));
      return _binedges[rbin];
    }

    double rBinMax(size_t rbin) const {
      assert(inRange(rbin, 0u, numBins()));
      return _binedges[rbin+1];
    }

    double rBinMid(size_t rbin) const {
      assert(inRange(rbin, 0u, numBins()));
      //cout << _binedges << '\n';
      return (_binedges[rbin] + _binedges[rbin+1])/2.0;
    }

    double diffJetShape(size_t ijet, size_t rbin) const {
      assert(inRange(ijet, 0u, numJets()));
      assert(inRange(rbin, 0u, numBins()));
      return _diffjetshapes[ijet][rbin];
    }

    double intJetShape(size_t ijet, size_t rbin) const {
      assert(inRange(ijet, 0u, numJets()));
      assert(inRange(rbin, 0u, numBins()));
      double rtn  = 0;
      for (size_t i = 0; i <= rbin; ++i) {
        rtn += _diffjetshapes[ijet][i];
      }
      return rtn;
    }

 
    // /// Return value of \f$ \Psi \f$ (integrated jet shape) at given radius for a \f$ p_T \f$ bin.
    // /// @todo Remove this external indexing thing
    // double psi(size_t pTbin) const {
    //   return _PsiSlot[pTbin];
    // }
 
 
  protected:

    void project(const Event& e);

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


    vector<double> _binedges;

    pair<double, double> _ptcuts;

    pair<double, double> _rapcuts;

    RapScheme _rapscheme;

 


    vector< vector<double> > _diffjetshapes;

 
  };
 
 
}
 
#endif