1// -*- C++ -*-
  2#include "Rivet/Analyses/MC_JETS_BASE.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/VetoedFinalState.hh"
  5#include "Rivet/Projections/LeptonFinder.hh"
  6#include "Rivet/Projections/MissingMomentum.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief MC validation analysis for e+mu W^+W^- + jets events
 12  class MC_WWJETS : public MC_JETS_BASE {
 13  public:
 14
 15    /// Default constructor
 16    MC_WWJETS()
 17      : MC_JETS_BASE("MC_WWJETS", 4, "Jets")
 18    {    }
 19
 20
 21    /// @name Analysis methods
 22    /// @{
 23
 24    /// Book histograms
 25    void init() {
 26
 27      declare("MET", MissingMomentum());
 28
 29      // Find electrons with cuts from input options
 30      const double ETAECUT = getOption<double>("ABSETAEMAX", 3.5);
 31      const double PTECUT = getOption<double>("PTEMIN", 25.);
 32      const Cut cut_e = Cuts::abseta < ETAECUT && Cuts::pT > PTECUT*GeV;
 33      LeptonFinder ef(0.2, cut_e && Cuts::abspid == PID::ELECTRON);
 34      declare(ef, "Elecs");
 35
 36      // Find muons with cuts from input options
 37      const double ETAMUCUT = getOption<double>("ABSETAMUMAX", 3.5);
 38      const double PTMUCUT = getOption<double>("PTMUMIN", 25.);
 39      const Cut cut_m = Cuts::abseta < ETAMUCUT && Cuts::pT > PTMUCUT*GeV;
 40      LeptonFinder mf(0.2, cut_m && Cuts::abspid == PID::MUON);
 41      declare(mf, "Muons");
 42
 43      // Set pT cut, jet alg, and clustering radius from input options
 44      const double jetptcut = getOption<double>("PTJMIN", 20.0);
 45      _jetptcut = jetptcut * GeV;
 46      const double R = getOption<double>("R", 0.4);
 47      JetAlg clusterAlgo;
 48      const string algoopt = getOption("ALGO", "ANTIKT");
 49      if ( algoopt == "KT" ) {
 50        clusterAlgo = JetAlg::KT;
 51      } else if ( algoopt == "CA" ) {
 52        clusterAlgo = JetAlg::CA;
 53      } else if ( algoopt == "ANTIKT" ) {
 54        clusterAlgo = JetAlg::ANTIKT;
 55      } else {
 56        MSG_WARNING("Unknown jet clustering algorithm option " + algoopt + ". Defaulting to anti-kT");
 57        clusterAlgo = JetAlg::ANTIKT;
 58      }
 59
 60      // Find jets with clustering radius from input option
 61      VetoedFinalState jetinput;
 62      jetinput
 63      .addVetoOnThisFinalState(ef)
 64      .addVetoOnThisFinalState(mf);
 65      FastJets fj(jetinput, clusterAlgo, R);
 66      declare(fj, "Jets");
 67
 68
 69      // Histogram correlations with jets and global quantities
 70      book(_h_WW_jet1_deta ,"WW_jet1_deta", 70, -7.0, 7.0);
 71      book(_h_WW_jet1_dR ,"WW_jet1_dR", 25, 1.5, 7.0);
 72      book(_h_We_jet1_dR ,"We_jet1_dR", 25, 0.0, 7.0);
 73      //
 74      book(_h_HT ,"HT", logspace(100, 100.0, 0.5*(sqrtS()>0.?sqrtS():14000.)));
 75      book(_h_jets_m_12 ,"jets_m_12", logspace(100, 1.0, 0.25*(sqrtS()>0.?sqrtS():14000.)));
 76
 77      MC_JETS_BASE::init();
 78    }
 79
 80
 81
 82    /// Do the analysis
 83    void analyze(const Event& event) {
 84
 85      // MET cut
 86      const P4& pmiss = apply<MissingMom>(event, "MET").missingMom();
 87      if (pmiss.pT() < 25*GeV) vetoEvent;
 88
 89      // Identify the closest-matching l+MET to m == mW
 90      /// @note Dubious strategy, given there are two neutrinos...
 91      const Particles& es = apply<LeptonFinder>(event, "Elecs").particles();
 92      const int iefound = closestMatchIndex(es, pmiss, Kin::mass, 80.4*GeV, 60*GeV, 100*GeV);
 93      const Particles& mus = apply<LeptonFinder>(event, "Muons").particles();
 94      const int imfound = closestMatchIndex(mus, pmiss, Kin::mass, 80.4*GeV, 60*GeV, 100*GeV);
 95
 96      // Require two valid W candidates
 97      if (iefound < 0 || imfound < 0) vetoEvent;
 98
 99      // Get momenta
100      FourMomentum pe = es[iefound].mom();
101      FourMomentum pm = mus[imfound].mom();
102      FourMomentum pww = pe + pm + pmiss;
103
104      // Get jets and compute obs relative to the WW objects
105      const Jets& jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > _jetptcut);
106      if (jets.size() > 0) {
107        _h_WW_jet1_deta->fill(pww.eta()-jets[0].eta());
108        _h_WW_jet1_dR->fill(deltaR(pww, jets[0].momentum()));
109        _h_We_jet1_dR->fill(deltaR(pe, jets[0].momentum()));
110      }
111      const double HT = sum(jets, Kin::pT, pe.pT() + pm.pT() + pmiss.pT());
112      if (HT > 0.0) _h_HT->fill(HT/GeV);
113
114      if (jets.size() > 1) {
115        const FourMomentum jet1 = jets[0].momentum();
116        const FourMomentum jet2 = jets[1].momentum();
117        _h_jets_m_12->fill((jet1+jet2).mass()/GeV);
118      }
119
120      MC_JETS_BASE::analyze(event);
121    }
122
123
124    /// Finalize
125    void finalize() {
126      const double norm = crossSection()/picobarn/sumOfWeights();
127      scale(_h_WW_jet1_deta, norm);
128      scale(_h_WW_jet1_dR, norm);
129      scale(_h_We_jet1_dR, norm);
130      scale(_h_jets_m_12, norm);
131      scale(_h_HT, norm);
132      MC_JETS_BASE::finalize();
133    }
134
135    /// @}
136
137
138  private:
139
140    /// @name Histograms
141    /// @{
142    Histo1DPtr _h_WW_jet1_deta;
143    Histo1DPtr _h_WW_jet1_dR;
144    Histo1DPtr _h_We_jet1_dR;
145    Histo1DPtr _h_jets_m_12;
146    Histo1DPtr _h_HT;
147    /// @}
148
149  };
150
151
152  RIVET_DECLARE_PLUGIN(MC_WWJETS);
153
154}