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/FinalState.hh"
#include "Rivet/Projections/JetFinder.hh"
#include "Rivet/Tools/RivetFastJet.hh"
 
#include "fastjet/ATLASConePlugin.hh"
#include "fastjet/CDFJetCluPlugin.hh"
#include "fastjet/CDFMidPointPlugin.hh"
#include "fastjet/CMSIterativeConePlugin.hh"
#include "fastjet/D0RunIIConePlugin.hh"
#include "fastjet/JadePlugin.hh"
#include "fastjet/SISConePlugin.hh"
#include "fastjet/TrackJetPlugin.hh"
#include "fastjet/contrib/VariableRPlugin.hh"
 
#include "Rivet/Projections/PxConePlugin.hh"
#include "Rivet/Tools/TypeTraits.hh"
 
namespace Rivet {
 
  typedef std::shared_ptr<fastjet::JetDefinition::Plugin> FJPluginPtr;

  class FastJets : public JetFinder {
  public:
 
    using JetFinder::operator=;
 


    FastJets(const FinalState& fsp,
             const fastjet::JetDefinition& jdef,
             JetMuons usemuons = JetMuons::ALL,
             JetInvisibles useinvis = JetInvisibles::NONE,
             fastjet::AreaDefinition* adef = nullptr)
        : JetFinder(fsp, usemuons, useinvis), _jdef(jdef), _adef(adef), _cuts(Cuts::OPEN) {
      _initBase();
    }

    FastJets(const FinalState& fsp,
             const fastjet::JetDefinition& jdef,
             const Cut& c,
             JetMuons usemuons = JetMuons::ALL,
             JetInvisibles useinvis = JetInvisibles::NONE,
             fastjet::AreaDefinition* adef = nullptr)
        : JetFinder(fsp, usemuons, useinvis), _jdef(jdef), _adef(adef), _cuts(c) {
      _initBase();
    }

    FastJets(const FinalState& fsp,
             const fastjet::JetDefinition& jdef,
             fastjet::AreaDefinition* adef,
             JetMuons usemuons = JetMuons::ALL,
             JetInvisibles useinvis = JetInvisibles::NONE)
        : FastJets(fsp, jdef, usemuons, useinvis, adef) { }

    FastJets(const FinalState& fsp,
             const fastjet::JetDefinition& jdef,
             fastjet::AreaDefinition* adef,
             const Cut& c,
             JetMuons usemuons = JetMuons::ALL,
             JetInvisibles useinvis = JetInvisibles::NONE)
        : FastJets(fsp, jdef, c, usemuons, useinvis, adef) { }
 

    FastJets(const FinalState& fsp,
             fastjet::JetAlgorithm type,
             fastjet::RecombinationScheme recom,
             double rparameter,
             JetMuons usemuons = JetMuons::ALL,
             JetInvisibles useinvis = JetInvisibles::NONE,
             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,
             const Cut& c,
             JetMuons usemuons = JetMuons::ALL,
             JetInvisibles useinvis = JetInvisibles::NONE,
             fastjet::AreaDefinition* adef = nullptr)
        : FastJets(fsp, fastjet::JetDefinition(type, rparameter, recom), c, usemuons, useinvis, adef) { }

    FastJets(const FinalState& fsp,
             fastjet::JetAlgorithm type,
             fastjet::RecombinationScheme recom,
             double rparameter,
             fastjet::AreaDefinition* adef,
             JetMuons usemuons = JetMuons::ALL,
             JetInvisibles useinvis = JetInvisibles::NONE)
        : FastJets(fsp, type, recom, rparameter, usemuons, useinvis, adef) { }

    FastJets(const FinalState& fsp,
             fastjet::JetAlgorithm type,
             fastjet::RecombinationScheme recom,
             double rparameter,
             fastjet::AreaDefinition* adef,
             const Cut& c,
             JetMuons usemuons = JetMuons::ALL,
             JetInvisibles useinvis = JetInvisibles::NONE)
        : FastJets(fsp, type, recom, rparameter, c, usemuons, useinvis, adef) { }

    FastJets(const FinalState& fsp,
             FJPluginPtr plugin,
             JetMuons usemuons = JetMuons::ALL,
             JetInvisibles useinvis = JetInvisibles::NONE,
             fastjet::AreaDefinition* adef = nullptr)
        : FastJets(fsp, fastjet::JetDefinition(plugin.get()), usemuons, useinvis, adef) {
      _plugin = plugin;
    }

    FastJets(const FinalState& fsp,
             FJPluginPtr plugin,
             const Cut& c,
             JetMuons usemuons = JetMuons::ALL,
             JetInvisibles useinvis = JetInvisibles::NONE,
             fastjet::AreaDefinition* adef = nullptr)
        : FastJets(fsp, fastjet::JetDefinition(plugin.get()), c, usemuons, useinvis, adef) {
      _plugin = plugin;
    }

    FastJets(const FinalState& fsp,
             FJPluginPtr plugin,
             fastjet::AreaDefinition* adef,
             JetMuons usemuons = JetMuons::ALL,
             JetInvisibles useinvis = JetInvisibles::NONE)
        : FastJets(fsp, plugin, usemuons, useinvis, adef) { }

    FastJets(const FinalState& fsp,
             FJPluginPtr plugin,
             fastjet::AreaDefinition* adef,
             const Cut& c,
             JetMuons usemuons = JetMuons::ALL,
             JetInvisibles useinvis = JetInvisibles::NONE)
        : FastJets(fsp, plugin, c, usemuons, useinvis, adef) { }

    FastJets(const FinalState& fsp,
             JetAlg alg,
             double rparameter = -1,
             JetMuons usemuons = JetMuons::ALL,
             JetInvisibles useinvis = JetInvisibles::NONE,
             fastjet::AreaDefinition* adef = nullptr)
        : JetFinder(fsp, usemuons, useinvis), _adef(adef), _cuts(Cuts::OPEN) {
      _initBase();
      _initJdef(alg, rparameter);
    }

    FastJets(const FinalState& fsp,
             JetAlg alg,
             double rparameter,
             const Cut& c,
             JetMuons usemuons = JetMuons::ALL,
             JetInvisibles useinvis = JetInvisibles::NONE,
             fastjet::AreaDefinition* adef = nullptr)
        : JetFinder(fsp, usemuons, useinvis), _adef(adef), _cuts(c) {
      _initBase();
      _initJdef(alg, rparameter);
    }
 

    RIVET_DEFAULT_PROJ_CLONE(FastJets);

 

    using Projection::operator=;
 


    static fastjet::JetDefinition mkJetDef(JetAlg alg, double rparameter);

    template <JetAlg JETALG, typename... Args>
    static FJPluginPtr mkPlugin(Args&&... args) {
      return std::make_shared<mapJetAlg2Plugin_t<JETALG>>(std::forward<Args>(args)...);
    }

    static FJPluginPtr mkPlugin(JetAlg alg, double rparameter = -1);

    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 = Particles(),
                       const Particles& tagparticles = Particles());

 

 
    typedef std::pair<PseudoJets, Particles> PJetsParts;

    template <typename CONTAINER, typename = std::enable_if_t<is_citerable_v<CONTAINER>, Jet>>
    static CONTAINER reclusterJets(const CONTAINER& jetsIn, const fastjet::JetDefinition& jDef) {
      PJetsParts reclusteredConsts = reclusterJetsParts(jetsIn, jDef);
      return mkTaggedJets(jetsIn, reclusteredConsts);
    }

    template <typename CONTAINER, typename = std::enable_if_t<is_citerable_v<CONTAINER>, Jet>>
    static CONTAINER reclusterJets(const CONTAINER& jetsIn,
                                   const fastjet::JetDefinition& jDef,
                                   const fastjet::Filter& filter) {
      PJetsParts reclusteredConsts = reclusterJetsParts(jetsIn, jDef);
      ifilterPseudoJets(reclusteredConsts.first, filter);
      return mkTaggedJets(jetsIn, reclusteredConsts);
    }

    template <typename CONTAINER, typename = std::enable_if_t<is_citerable_v<CONTAINER>, Jet>>
    static CONTAINER reclusterJets(const CONTAINER& jetsIn, const FJPluginPtr& jetAlg) {
      const fastjet::JetDefinition jDef(jetAlg.get());
      return reclusterJets(jetsIn, jDef);
    }

    template <typename CONTAINER, typename = std::enable_if_t<is_citerable_v<CONTAINER>, Jet>>
    static CONTAINER reclusterJets(const CONTAINER& jetsIn,
                                   const FJPluginPtr& jetAlg,
                                   const fastjet::Filter& filter) {
      const fastjet::JetDefinition jDef(jetAlg.get());
      return reclusterJets(jetsIn, jDef, filter);
    }

    template <JetAlg JETALG,
              typename... Args,
              typename CONTAINER,
              typename = std::enable_if_t<is_citerable_v<CONTAINER>, Jet>>
    static CONTAINER reclusterJets(const CONTAINER& jetsIn, Args&&... args) {
      if constexpr (JETALG < JetAlg::SISCONE) { // JetDefinition
        const fastjet::JetDefinition jDef = mkJetDef(JETALG, std::forward<Args>(args)...);
        return reclusterJets(jetsIn, jDef);
      }
      else { // Plugin
        const FJPluginPtr plugin = mkPlugin<JETALG>(std::forward<Args>(args)...);
        return reclusterJets(jetsIn, plugin);
      }
      throw std::invalid_argument("Unknown jet algorithm: " + to_string(int(JETALG)));
    }

    template <JetAlg JETALG,
              typename... Args,
              typename CONTAINER,
              typename = std::enable_if_t<is_citerable_v<CONTAINER>, Jet>>
    static CONTAINER reclusterJets(const CONTAINER& jetsIn, const fastjet::Filter& filter, Args&&... args) {
      if constexpr (JETALG < JetAlg::SISCONE) { // JetDefinition
        const fastjet::JetDefinition jDef = mkJetDef(JETALG, std::forward<Args>(args)...);
        return reclusterJets(jetsIn, jDef, filter);
      }
      else { // Plugin
        const FJPluginPtr plugin = mkPlugin<JETALG>(std::forward<Args>(args)...);
        return reclusterJets(jetsIn, plugin, filter);
      }
      throw std::invalid_argument("Unknown jet algorithm: " + to_string(int(JETALG)));
    }

    template <typename T, typename U, typename... Args>
    static std::map<T, U> reclusterJets(const std::map<T, U>& jetsMap, Args&&... args) {
      std::map<T, U> rtn;
      for (const auto& [key, jetsIn] : jetsMap) rtn[key] = reclusterJets(jetsIn, std::forward<Args>(args)...);
      return rtn;
    }

    template <JetAlg JETALG, typename T, typename U, typename... Args>
    static std::map<T, U> reclusterJets(const std::map<T, U>& jetsMap, Args&&... args) {
      std::map<T, U> rtn;
      for (const auto& [key, jetsIn] : jetsMap)
        rtn[key] = reclusterJets<JETALG>(jetsIn, std::forward<Args>(args)...);
      return rtn;
    }

 

    void reset();
 


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

    void clearJetArea() {
      _adef.reset();
    }

 


    void addTrf(fastjet::Transformer* trf) {
      _trfs.push_back(shared_ptr<fastjet::Transformer>(trf));
    }

    template <typename TRFS, typename TRF = typename TRFS::value_type>
    typename std::enable_if<Derefable<TRF>::value, void>::type addTrfs(const TRFS& trfs) {
      for (auto& trf : trfs) addTrf(trf);
    }

    void addFilter(fastjet::Filter* filter) {
      addTrf(filter);
    }

    template <typename FILTERS, typename FILTER = typename FILTERS::value_type>
    typename std::enable_if<Derefable<FILTER>::value, void>::type addFilters(const FILTERS& filters) {
      addTrfs(filters);
    }

    void clearTrfs() {
      _trfs.clear();
    }

 


    Jets _jets() const;

    PseudoJets pseudojets(double ptmin = 0.0) const;

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

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

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

 


    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;
    }

 
 
  protected:

    void _initBase();
 
 
    void _initJdef(JetAlg alg, double rparameter) {
      MSG_DEBUG("JetAlg = " << static_cast<int>(alg));
      MSG_DEBUG("R parameter = " << rparameter);
      if (alg < JetAlg::SISCONE) { // fastjet::JetDefinitions
        _jdef = mkJetDef(alg, rparameter);
      }
      else { // fastjet::JetDefinition::Plugins
        _plugin = mkPlugin(alg, rparameter);
        _jdef = fastjet::JetDefinition(_plugin.get());
      }
    }

    static Jets mkTaggedJets(const Jets& jetsIn, const PJetsParts& pJetsParts);

    static PJetsParts reclusterJetsParts(const Jets& jetsIn, const fastjet::JetDefinition& jDef);
 
  protected:

    void project(const Event& e);

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

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

    fastjet::JetDefinition _jdef;

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

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

    Cut _cuts;

    FJPluginPtr _plugin;

    std::vector<std::shared_ptr<fastjet::Transformer>> _trfs;

    static const std::map<JetAlg, std::pair<fastjet::JetAlgorithm, fastjet::RecombinationScheme>> jetAlgMap;

    template <JetAlg, typename X>
    struct mapJetAlg2Plugin;

    template <typename X>
    struct mapJetAlg2Plugin<JetAlg::SISCONE, X> {
      using type = fastjet::SISConePlugin;
    };
    template <typename X>
    struct mapJetAlg2Plugin<JetAlg::PXCONE, X> {
      using type = Rivet::PxConePlugin;
    };
    template <typename X>
    struct mapJetAlg2Plugin<JetAlg::CDFJETCLU, X> {
      using type = fastjet::CDFJetCluPlugin;
    };
    template <typename X>
    struct mapJetAlg2Plugin<JetAlg::CDFMIDPOINT, X> {
      using type = fastjet::CDFMidPointPlugin;
    };
    template <typename X>
    struct mapJetAlg2Plugin<JetAlg::D0ILCONE, X> {
      using type = fastjet::D0RunIIConePlugin;
    };
    template <typename X>
    struct mapJetAlg2Plugin<JetAlg::JADE, X> {
      using type = fastjet::JadePlugin;
    };
    template <typename X>
    struct mapJetAlg2Plugin<JetAlg::TRACKJET, X> {
      using type = fastjet::TrackJetPlugin;
    };
    template <typename X>
    struct mapJetAlg2Plugin<JetAlg::VARIABLER, X> {
      using type = fastjet::contrib::VariableRPlugin;
    };

    template <JetAlg JETALG>
    using mapJetAlg2Plugin_t = typename mapJetAlg2Plugin<JETALG, void>::type;

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

    Particles _fsparticles, _tagparticles;
  };
 
}
 
#endif