1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/PromptFinalState.hh"
  4#include "Rivet/Projections/VetoedFinalState.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6#include "Rivet/Projections/LeptonFinder.hh"
  7#include "Rivet/Projections/MissingMomentum.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// W + jets jet multiplicities and pT
 13  class ATLAS_2010_I882534 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2010_I882534);
 18
 19
 20    /// @name Analysis methods
 21    /// @{
 22
 23    /// Book histograms and initialise projections before the run
 24    void init() {
 25
 26      declare("MET", MissingMomentum());
 27
 28      // Find electrons and muons
 29      const Cut cut_e = (Cuts::abseta < 1.37 || Cuts::absetaIn(1.52, 2.47)) && Cuts::pT > 20*GeV;
 30      LeptonFinder ef(0.1, cut_e && Cuts::abspid == PID::ELECTRON);
 31      declare(ef, "Elecs");
 32      const Cut cut_m = Cuts::abseta < 2.4 && Cuts::pT > 20*GeV;
 33      LeptonFinder mf(0.1, cut_m && Cuts::abspid == PID::MUON);
 34      declare(mf, "Muons");
 35
 36      // Input for the jets: no neutrinos, no muons, and no electron which passed the electron cuts
 37      VetoedFinalState veto;
 38      veto.addVetoOnThisFinalState(ef);
 39      veto.addVetoOnThisFinalState(mf);
 40      FastJets jets(veto, JetAlg::ANTIKT, 0.4, JetMuons::NONE);
 41      declare(jets, "Jets");
 42
 43      /// Book histograms
 44      book(_h_el_njet_inclusive,1,1,1);
 45      book(_h_mu_njet_inclusive,2,1,1);
 46      book(_h_el_pT_jet1,5,1,1);
 47      book(_h_mu_pT_jet1,6,1,1);
 48      book(_h_el_pT_jet2,7,1,1);
 49      book(_h_mu_pT_jet2,8,1,1);
 50    }
 51
 52
 53    /// Perform the per-event analysis
 54    void analyze(const Event& event) {
 55
 56      if (_edges.empty())  _edges = _h_mu_njet_inclusive->xEdges();
 57
 58      const Jets& jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 20*GeV);
 59
 60      // MET cut
 61      const P4& pmiss = apply<MissingMom>(event, "MET").missingMom();
 62      if (pmiss.Et() < 25*GeV) vetoEvent;
 63
 64      // Identify the closest-matching l+MET to m == mW
 65      const Particles& es = apply<LeptonFinder>(event, "Elecs").particles();
 66      const int iefound = closestMatchIndex(es, pmiss, Kin::mass, 80.4*GeV, 0*GeV, 1000*GeV);
 67      const Particles& mus = apply<LeptonFinder>(event, "Muons").particles();
 68      const int imfound = closestMatchIndex(mus, pmiss, Kin::mass, 80.4*GeV, 0*GeV, 1000*GeV);
 69
 70      // Require two valid W candidates
 71      if (iefound < 0 && imfound < 0) vetoEvent; //< no W's
 72      if (iefound >= 0 && imfound >= 0) vetoEvent; //< multi-W
 73
 74      // Histogramming
 75      if (iefound >= 0) {
 76        const Particle& e = es[iefound].constituents()[0];
 77        if (mT(pmiss, e) > 40*GeV) {
 78          const Jets js = select(jets, [&](const Jet& j) { return j.abseta() < 2.8 && deltaR(e, j) > 0.5; });
 79          _h_el_njet_inclusive->fill(_edges[0]);
 80          if (js.size() >= 1) {
 81            _h_el_njet_inclusive->fill(_edges[1]);
 82            _h_el_pT_jet1->fill(js[0].pT());
 83          }
 84          if (js.size() >= 2) {
 85            _h_el_njet_inclusive->fill(_edges[2]);
 86            _h_el_pT_jet2->fill(js[1].pT());
 87          }
 88          if (js.size() >= 3) {
 89            _h_el_njet_inclusive->fill(_edges[3]);
 90          }
 91        }
 92      }
 93
 94      if (imfound >= 0) {
 95        const Particle& mu = mus[imfound];
 96        if (mT(pmiss, mu) > 40*GeV) {
 97          const Jets js = select(jets, [&](const Jet& j) { return j.abseta() < 2.8 && deltaR(mu, j) > 0.5; });
 98          _h_mu_njet_inclusive->fill(_edges[0]);
 99          if (js.size() >= 1) {
100            _h_mu_njet_inclusive->fill(_edges[1]);
101            _h_mu_pT_jet1->fill(js[0].pT());
102          }
103          if (js.size() >= 2) {
104            _h_mu_njet_inclusive->fill(_edges[2]);
105            _h_mu_pT_jet2->fill(js[1].pT());
106          }
107          if (js.size() >= 3) {
108            _h_mu_njet_inclusive->fill(_edges[3]);
109          }
110          if (js.size() >= 4) {
111            _h_mu_njet_inclusive->fill(_edges[4]);
112          }
113        }
114      }
115
116    }
117
118
119    /// Normalise histograms etc., after the run
120    void finalize() {
121      double normfac = crossSection()/picobarn/sumOfWeights();
122      scale(_h_el_njet_inclusive, normfac/nanobarn);
123      scale(_h_mu_njet_inclusive, normfac/nanobarn);
124      scale(_h_el_pT_jet1, normfac);
125      scale(_h_mu_pT_jet1, normfac);
126      scale(_h_el_pT_jet2, normfac);
127      scale(_h_mu_pT_jet2, normfac);
128    }
129
130    ///@}
131
132
133  private:
134
135    /// @name Histograms
136    /// @{
137    BinnedHistoPtr<string> _h_el_njet_inclusive;
138    BinnedHistoPtr<string> _h_mu_njet_inclusive;
139    Histo1DPtr _h_el_pT_jet1;
140    Histo1DPtr _h_mu_pT_jet1;
141    Histo1DPtr _h_el_pT_jet2;
142    Histo1DPtr _h_mu_pT_jet2;
143    vector<string> _edges;
144    /// @}
145
146  };
147
148
149  RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2010_I882534, ATLAS_2010_S8919674);
150
151}