FastJets.hh

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

#include "Rivet/Jet.hh"
#include "Rivet/Particle.hh"
#include "Rivet/Projection.hh"
#include "Rivet/Projections/JetAlg.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Tools/RivetFastJet.hh"

#include "fastjet/SISConePlugin.hh"
#include "fastjet/ATLASConePlugin.hh"
#include "fastjet/CMSIterativeConePlugin.hh"
#include "fastjet/CDFJetCluPlugin.hh"
#include "fastjet/CDFMidPointPlugin.hh"
#include "fastjet/D0RunIIConePlugin.hh"
#include "fastjet/TrackJetPlugin.hh"
#include "fastjet/JadePlugin.hh"
#include "Rivet/Projections/PxConePlugin.hh"
//#include "fastjet/PxConePlugin.hh"

namespace Rivet {


  class FastJets : public JetAlg {
  public:

    enum JetAlgName { KT, CAM, SISCONE, ANTIKT, PXCONE,
                      ATLASCONE, CMSCONE,
                      CDFJETCLU, CDFMIDPOINT, D0ILCONE,
                      JADE, DURHAM, TRACKJET, GENKTEE };




    FastJets(const FinalState& fsp,
             const fastjet::JetDefinition& jdef,
             JetAlg::MuonsStrategy usemuons=JetAlg::ALL_MUONS,
             JetAlg::InvisiblesStrategy useinvis=JetAlg::NO_INVISIBLES,
             fastjet::AreaDefinition* adef=nullptr)
      : JetAlg(fsp, usemuons, useinvis), _jdef(jdef), _adef(adef)
    {
      _initBase();
    }

    FastJets(const FinalState& fsp,
             const fastjet::JetDefinition& jdef,
             fastjet::AreaDefinition* adef,
             JetAlg::MuonsStrategy usemuons=JetAlg::ALL_MUONS,
             JetAlg::InvisiblesStrategy useinvis=JetAlg::NO_INVISIBLES)
      : FastJets(fsp, jdef, usemuons, useinvis, adef)
    {    }

    FastJets(const FinalState& fsp,
             fastjet::JetAlgorithm type,
             fastjet::RecombinationScheme recom, double rparameter,
             JetAlg::MuonsStrategy usemuons=JetAlg::ALL_MUONS,
             JetAlg::InvisiblesStrategy useinvis=JetAlg::NO_INVISIBLES,
             fastjet::AreaDefinition* adef=nullptr)
      : FastJets(fsp, fastjet::JetDefinition(type, rparameter, recom), usemuons, useinvis, adef)
    {    }

    FastJets(const FinalState& fsp,
             fastjet::JetAlgorithm type,
             fastjet::RecombinationScheme recom, double rparameter,
             fastjet::AreaDefinition* adef,
             JetAlg::MuonsStrategy usemuons=JetAlg::ALL_MUONS,
             JetAlg::InvisiblesStrategy useinvis=JetAlg::NO_INVISIBLES)
      : FastJets(fsp, type, recom, rparameter, usemuons, useinvis, adef)
    {    }

    FastJets(const FinalState& fsp,
             fastjet::JetDefinition::Plugin* plugin,
             JetAlg::MuonsStrategy usemuons=JetAlg::ALL_MUONS,
             JetAlg::InvisiblesStrategy useinvis=JetAlg::NO_INVISIBLES,
             fastjet::AreaDefinition* adef=nullptr)
      : FastJets(fsp, fastjet::JetDefinition(plugin), usemuons, useinvis, adef)
    {
      _plugin.reset(plugin);
    }

    FastJets(const FinalState& fsp,
             fastjet::JetDefinition::Plugin* plugin,
             fastjet::AreaDefinition* adef,
             JetAlg::MuonsStrategy usemuons=JetAlg::ALL_MUONS,
             JetAlg::InvisiblesStrategy useinvis=JetAlg::NO_INVISIBLES)
      : FastJets(fsp, plugin, usemuons, useinvis, adef)
    {    }

    FastJets(const FinalState& fsp,
             JetAlgName alg, double rparameter,
             JetAlg::MuonsStrategy usemuons=JetAlg::ALL_MUONS,
             JetAlg::InvisiblesStrategy useinvis=JetAlg::NO_INVISIBLES,
             fastjet::AreaDefinition* adef=nullptr,
             double seed_threshold=1.0)
      : JetAlg(fsp, usemuons, useinvis)
    {
      _initBase();
      _initJdef(alg, rparameter, seed_threshold);
    }


    // /// Same thing as above, but without an FS (for when we want to pass the particles directly to the calc method)
    // /// @todo Does this work properly, without internal HeavyQuarks etc.?
    // FastJets(JetAlgName alg, double rparameter, double seed_threshold=1.0) { _initJdef(alg, rparameter, seed_threshold); }
    // /// Same thing as above, but without an FS (for when we want to pass the particles directly to the calc method)
    // /// @todo Does this work properly, without internal HeavyQuarks etc.?
    // FastJets(fastjet::JetAlgorithm type, fastjet::RecombinationScheme recom, double rparameter) { _initJdef(type, recom, rparameter); }
    // /// Same thing as above, but without an FS (for when we want to pass the particles directly to the calc method)
    // /// @todo Does this work properly, without internal HeavyQuarks etc.?
    // FastJets(fastjet::JetDefinition::Plugin* plugin) : _jdef(plugin), _plugin(plugin) {
    //   // _plugin.reset(plugin);
    //   // _jdef = fastjet::JetDefinition(plugin);
    // }


    DEFAULT_RIVET_PROJ_CLONE(FastJets);





    static PseudoJets mkClusterInputs(const Particles& fsparticles, const Particles& tagparticles=Particles());
    static Jet mkJet(const PseudoJet& pj, const Particles& fsparticles, const Particles& tagparticles=Particles());
    static Jets mkJets(const PseudoJets& pjs, const Particles& fsparticles, const Particles& tagparticles=Particles());



    void reset();


    void useJetArea(fastjet::AreaDefinition* adef) {
      _adef.reset(adef);
    }




    Jets _jets() const;

    PseudoJets pseudoJets(double ptmin=0.0) const;
    PseudoJets pseudojets(double ptmin=0.0) const { return pseudoJets(ptmin); }

    PseudoJets pseudoJetsByPt(double ptmin=0.0) const {
      return sorted_by_pt(pseudoJets(ptmin));
    }
    PseudoJets pseudojetsByPt(double ptmin=0.0) const { return pseudoJetsByPt(ptmin); }

    PseudoJets pseudoJetsByE(double ptmin=0.0) const {
      return sorted_by_E(pseudoJets(ptmin));
    }
    PseudoJets pseudojetsByE(double ptmin=0.0) const { return pseudoJetsByE(ptmin); }

    PseudoJets pseudoJetsByRapidity(double ptmin=0.0) const {
      return sorted_by_rapidity(pseudoJets(ptmin));
    }
    PseudoJets pseudojetsByRapidity(double ptmin=0.0) const { return pseudoJetsByRapidity(ptmin); }

    Jet trimJet(const Jet& input, const fastjet::Filter& trimmer) const;





    const shared_ptr<fastjet::ClusterSequence> clusterSeq() const {
      return _cseq;
    }

    const shared_ptr<fastjet::ClusterSequenceArea> clusterSeqArea() const {
      return areaDef() ? dynamic_pointer_cast<fastjet::ClusterSequenceArea>(_cseq) : nullptr;
    }

    const fastjet::JetDefinition& jetDef() const {
      return _jdef;
    }

    const shared_ptr<fastjet::AreaDefinition> areaDef() const {
      return _adef;
    }



  private:

    void _initBase();
    void _initJdef(JetAlgName alg, double rparameter, double seed_threshold);

  protected:

    void project(const Event& e);

    int compare(const Projection& p) const;

  public:

    void calc(const Particles& fsparticles, const Particles& tagparticles=Particles());


  private:

    fastjet::JetDefinition _jdef;

    std::shared_ptr<fastjet::AreaDefinition> _adef;

    std::shared_ptr<fastjet::ClusterSequence> _cseq;

    std::shared_ptr<fastjet::JetDefinition::Plugin> _plugin;

    mutable std::map<int, vector<double> > _yscales;

    Particles _fsparticles, _tagparticles;

  };

}

#endif