1// -*- C++ -*-
  2#include "Rivet/Analyses/MC_JETS_BASE.hh"
  3#include "Rivet/Projections/DileptonFinder.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief MC validation analysis for Higgs [-> tau tau] + jets events
 10  class MC_HJETS : public MC_JETS_BASE {
 11  public:
 12
 13    /// Default constructor
 14    MC_HJETS()
 15      : MC_JETS_BASE("MC_HJETS", 4, "Jets")
 16    {    }
 17
 18
 19    /// @name Analysis methods
 20    /// @{
 21
 22    /// Book histograms
 23    void init() {
 24      // set FS cuts from input options
 25      const double etacut = getOption<double>("ABSETATAUMAX", 3.5);
 26      const double ptcut = getOption<double>("PTTAUMIN", 25.);
 27
 28      Cut cut = Cuts::abseta < etacut && Cuts::pT > ptcut*GeV;
 29      /// @todo Hmm, FS taus??
 30      DileptonFinder hfinder(125*GeV, 0.0, cut && Cuts::abspid == PID::TAU, Cuts::massIn(115*GeV, 135*GeV));
 31
 32      declare(hfinder, "Hfinder");
 33
 34      // set ptcut from input option
 35      _jetptcut = getOption<double>("PTJMIN", 20.0) * GeV;
 36
 37      // set clustering radius from input option
 38      const double R = getOption<double>("R", 0.4);
 39
 40      // set clustering algorithm from input option
 41      JetAlg clusterAlgo;
 42      const string algoopt = getOption("ALGO", "ANTIKT");
 43      if ( algoopt == "KT" ) {
 44        clusterAlgo = JetAlg::KT;
 45      } else if ( algoopt == "CA" ) {
 46        clusterAlgo = JetAlg::CA;
 47      } else if ( algoopt == "ANTIKT" ) {
 48        clusterAlgo = JetAlg::ANTIKT;
 49      } else {
 50        MSG_WARNING("Unknown jet clustering algorithm option " + algoopt + ". Defaulting to anti-kT");
 51        clusterAlgo = JetAlg::ANTIKT;
 52      }
 53
 54      FastJets jetpro(hfinder.remainingFinalState(), clusterAlgo, R);
 55      declare(jetpro, "Jets");
 56
 57      book(_h_H_jet1_deta ,"H_jet1_deta", 50, -5.0, 5.0);
 58      book(_h_H_jet1_dR ,"H_jet1_dR", 25, 0.5, 7.0);
 59
 60      MC_JETS_BASE::init();
 61    }
 62
 63
 64
 65    /// Do the analysis
 66    void analyze(const Event & e) {
 67      const DileptonFinder& hfinder = apply<DileptonFinder>(e, "Hfinder");
 68      if (hfinder.bosons().size() != 1) vetoEvent;
 69
 70      FourMomentum hmom(hfinder.bosons()[0].momentum());
 71      const Jets& jets = apply<FastJets>(e, "Jets").jetsByPt(Cuts::pT > _jetptcut);
 72      if (jets.size() > 0) {
 73        _h_H_jet1_deta->fill(hmom.eta()-jets[0].eta());
 74        _h_H_jet1_dR->fill(deltaR(hmom, jets[0].momentum()));
 75      }
 76
 77      MC_JETS_BASE::analyze(e);
 78    }
 79
 80
 81    /// Finalize
 82    void finalize() {
 83      normalize(_h_H_jet1_deta, crossSection()/picobarn);
 84      normalize(_h_H_jet1_dR, crossSection()/picobarn);
 85      MC_JETS_BASE::finalize();
 86    }
 87
 88    /// @}
 89
 90
 91  private:
 92
 93    /// @name Histograms
 94    /// @{
 95    Histo1DPtr _h_H_jet1_deta;
 96    Histo1DPtr _h_H_jet1_dR;
 97    /// @}
 98
 99  };
100
101
102
103  RIVET_DECLARE_PLUGIN(MC_HJETS);
104
105}