RivetHTT.hh

// -*- C++ -*-
#ifndef RIVET_RIVETHTT_HH
#define RIVET_RIVETHTT_HH
 
#include "Rivet/Jet.hh"
#include "HEPTopTagger/HEPTopTagger.hh"
 
namespace Rivet {
 
 
  enum class HTTMode {
    EARLY_MASSRATIO_SORT_MASS,     // applies 2D mass plane requirements then select the candidate which minimizes |m_cand-mt|
    LATE_MASSRATIO_SORT_MASS,      // selects the candidate which minimizes |m_cand-mt|
    EARLY_MASSRATIO_SORT_MODDJADE, // applies the 2D mass plane requirements then select the candidate with highest jade distance
    LATE_MASSRATIO_SORT_MODDJADE,  // selects the candidate with highest modified jade distance
    TWO_STEP_FILTER                // only analyzes the candidate built with the highest pT(t) after unclustering
  };
  
  
  class HTT {
  public:
 
    struct InputParameters {
 
      Mode mode = HTTMode::EARLY_MASSRATIO_SORT_MASS;
 
      bool do_optimalR = true;
      double optimalR_min = 0.5; // min jet size
      double optimalR_step = 0.1; // step size
      double optimalR_threshold = 0.2; // step size
 
      double mass_drop = 0.8;
      double max_subjet_mass = 30*GeV;
 
      unsigned int filt_N = 5; // set_nfilt
      double filtering_R = 0.3; // max subjet distance for filtering
      double filtering_minpt = 0.; // min subjet pt for filtering
 
      JetAlg filtering_algorithm = JetAlg::CA;
 
      JetAlg reclustering_algorithm = JetAlg::CA;
 
      double top_mass = 172.3*GeV;
      double W_mass = 80.4*GeV;
      double Mtop_min = 150*GeV;
      double Mtop_max = 200*GeV; //set_top_range(min,max)
 
      double fw = 0.15;
      double mass_ratio_range_min = (1.-fw)*W_mass/top_mass;
      double mass_ratio_range_max = (1.+fw)*W_mass/top_mass;
 
      double m23cut = 0.35;
      double m13cutmin = 0.2;
      double m13cutmax = 1.3;
 
      double prune_zcut = 0.1;
      double prune_rcut = 0.5;
 
    };
 
 
 
    HTT() {}
 
    HTT(HTT::InputParameters& params) {
      setParams(params);
    }
 
    // /// Destructor
    // ~HTT() {}
 
    void setParams(HTT::InputParameters& params);
 
    void calc(Jet& jet);
 
    const Jet topJet() const;
    const Jet bJet() const;
    const Jet wJet() const;
    const Jet w1Jet() const;
    const Jet w2Jet() const;
 
    const PseudoJet& topJet() const;
    const PseudoJet& bJet() const;
    const PseudoJet& wJet() const;
    const PseudoJet& w1Jet() const;
    const PseudoJet& w2Jet() const;
 
    const Jets& subjets() const;
 
    // /// Print tagger information
    // void info() const;
    // /// Print tagger settings
    // void settings() const;
 
    double prunedMass() const;
 
    // The unfiltered mass
    double unfilteredMass() const;
 
    double deltaTopMass() const;
 
    bool isTopTagged() const;
 
    bool passedMassCutTop() const;
 
    bool passedMassCut2D() const;
 
 
  private:
 
    fastjet::HEPTopTagger::HEPTopTagger _tagger;
 
  };
 
 
 
  void HTT::setParams(HTT::InputParameters& params) {
    _tagger = fastjet::HEPTopTagger::HEPTopTagger();
 
    // Optimal R
    _tagger.do_optimalR(params.do_optimalR);
    _tagger.set_optimalR_min(params.optimalR_min);
    _tagger.set_optimalR_step(params.optimalR_step);
    _tagger.set_optimalR_threshold(params.optimalR_threshold);
 
    // Candidate selection
    fastjet::HEPTopTagger::Mode mode;
    if (params.mode == HTTMode::EARLY_MASSRATIO_SORT_MASS) {
      mode = fastjet::HEPTopTagger::EARLY_MASSRATIO_SORT_MASS;
    } else if (params.mode == HTTMode::LATE_MASSRATIO_SORT_MASS) {
      mode = fastjet::HEPTopTagger::LATE_MASSRATIO_SORT_MASS;
    } else if (params.mode == HTTMode::EARLY_MASSRATIO_SORT_MODDJADE) {
      mode = fastjet::HEPTopTagger::EARLY_MASSRATIO_SORT_MODDJADE;
    } else if (params.mode == HTTMode::LATE_MASSRATIO_SORT_MODDJADE) {
      mode = fastjet::HEPTopTagger::LATE_MASSRATIO_SORT_MODDJADE;
    } else {
      mode = fastjet::HEPTopTagger::TWO_STEP_FILTER;
    }
    _tagger.set_mode(mode);
    _tagger.set_mt(params.top_mass);
    _tagger.set_mw(params.W_mass);
    _tagger.set_top_mass_range(params.Mtop_min, params.Mtop_max);
    _tagger.set_fw(params.fw);
    _tagger.set_mass_ratio_range(params.mass_ratio_range_min, params.mass_ratio_range_max);
    _tagger.set_mass_ratio_cut(params.m23cut, params.m13cutmin, params.m13cutmax);
 
    // Filtering
    _tagger.set_filtering_n(params.filt_N);
    _tagger.set_filtering_R(params.filtering_R);
    _tagger.set_filtering_minpt_subjet(params.filtering_minpt);
 
    fastjet::JetAlgorithm algo = fastjet::antikt_algorithm;
    if (params.filtering_algorithm == Algo::CA) algo = fastjet::cambridge_algorithm;
    else if (params.filtering_algorithm == Algo::KT) algo = fastjet::kt_algorithm;
    _tagger.set_filtering_jetalgorithm(algo);
 
    // Reclustering
    algo = fastjet::antikt_algorithm;
    if (params.reclustering_algorithm == Algo::CA) algo = fastjet::cambridge_algorithm;
    else if (params.reclustering_algorithm == Algo::KT) algo = fastjet::kt_algorithm;
    _tagger.set_reclustering_jetalgorithm(algo);
 
    // Mass-drop
    _tagger.set_mass_drop_threshold(params.mass_drop);
    _tagger.set_mass_drop_threshold(params.mass_drop);
 
    // Pruning
    _tagger.set_pruning_rcut_factor(params.prune_rcut);
    _tagger.set_pruning_zcut(params.prune_zcut);
  }
 
 
  void HTT::calc(Jet& jet) {
    _tagger.run(jet.pseudojet());
  }
 
 
  const Jet HTT::topJet() const { return Jet(topPjet()); }
  const Jet HTT::bJet() const { return Jet(bPjet()); }
  const Jet HTT::wJet() const { return Jet(wPjet()); }
  const Jet HTT::w1Jet() const { return Jet(w1Pjet()); }
  const Jet HTT::w2Jet() const { return Jet(w2Pjet()); }
 
 
  const PseudoJet& HTT::topPjet() const { return _tagger.t(); }
  const PseudoJet& HTT::bPjet() const { return _tagger.b(); }
  const PseudoJet& HTT::wPjet() const { return _tagger.W(); }
  const PseudoJet& HTT::w1Pjet() const { return _tagger.W1(); }
  const PseudoJet& HTT::w2Pjet() const { return _tagger.W2(); }
 
 
  Jets HTT::subjets() const {
    Jets rtn;
    rtn.reserve(3);
    rtn.emplace_back(_tagger.j1());
    rtn.emplace_back(_tagger.j2());
    rtn.emplace_back(_tagger.j3());
    return rtn;
  }
 
 
  // void HTT::info() const { _tagger.get_info(); }
  // void HTT::settings() const { _tagger.get_setting(); }
 
  double HTT::prunedMass() const { return _tagger.pruned_mass(); }
 
  double HTT::unfilteredMass() const { return _tagger.unfiltered_mass(); }
 
  double HTT::deltaTopMass() const { return _tagger.delta_top(); }
 
  bool HTT::isTopTagged() const { return _tagger.is_tagged(); }
 
  bool HTT::passedMassCutTop() const { return _tagger.is_maybe_top(); }
 
  bool HTT::passedMassCut2D() const { return _tagger.is_masscut_passed(); }
 
 
 
}
 
#endif