1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/PromptFinalState.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6#include "Rivet/Projections/Sphericity.hh"
  7#include "Rivet/Projections/SmearedParticles.hh"
  8#include "Rivet/Projections/SmearedJets.hh"
  9#include "Rivet/Projections/SmearedMET.hh"
 10
 11namespace Rivet {
 12
 13
 14  /// @brief ATLAS 2016 1-lepton + many jets SUSY search, from 14.8/fb CONF note
 15  class ATLAS_2016_CONF_2016_094 : public Analysis {
 16  public:
 17
 18    /// Constructor
 19    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2016_CONF_2016_094);
 20
 21
 22    /// @name Analysis methods
 23    /// @{
 24
 25    /// Book histograms and initialise projections before the run
 26    void init() {
 27
 28      // Initialise and register projections
 29      FinalState calofs(Cuts::abseta < 4.9);
 30      FastJets fj(calofs, JetAlg::ANTIKT, 0.4);
 31      declare(fj, "TruthJets");
 32      declare(SmearedJets(fj, JET_SMEAR_ATLAS_RUN2, [](const Jet& j) {
 33            if (j.abseta() > 2.5) return 0.;
 34            return j.bTagged(Cuts::pT > 5*GeV) ? 0.80 :
 35              j.cTagged(Cuts::pT > 5*GeV) ? 1/6. : 1/106.; }), "Jets");
 36
 37      // MissingMomentum mm(calofs);
 38      // declare(mm, "TruthMET");
 39      // declare(SmearedMET(mm, MET_SMEAR_ATLAS_RUN2), "MET");
 40
 41      FinalState es(Cuts::abspid == PID::ELECTRON && Cuts::abseta < 2.47 && !Cuts::absetaIn(1.37, 1.52) && Cuts::pT > 10*GeV);
 42      declare(es, "TruthElectrons");
 43      declare(SmearedParticles(es, ELECTRON_RECOEFF_ATLAS_RUN2, ELECTRON_SMEAR_ATLAS_RUN2), "Electrons");
 44
 45      FinalState mus(Cuts::abspid == PID::MUON && Cuts::abseta < 2.4 && Cuts::pT > 10*GeV);
 46      declare(mus, "TruthMuons");
 47      declare(SmearedParticles(mus, MUON_EFF_ATLAS_RUN2, MUON_SMEAR_ATLAS_RUN2), "Muons");
 48
 49
 50      // Book histograms/counters
 51      book(_h_08j40_0b,"08j40_0b");
 52      book(_h_09j40_0b,"09j40_0b");
 53      book(_h_10j40_0b,"10j40_0b");
 54      book(_h_08j40_3b,"08j40_3b");
 55      book(_h_09j40_3b,"09j40_3b");
 56      book(_h_10j40_3b,"10j40_3b");
 57      book(_h_08j60_0b,"08j60_0b");
 58      book(_h_09j60_0b,"09j60_0b");
 59      book(_h_10j60_0b,"10j60_0b");
 60      book(_h_08j60_3b,"08j60_3b");
 61      book(_h_09j60_3b,"09j60_3b");
 62      book(_h_10j60_3b,"10j60_3b");
 63
 64    }
 65
 66
 67    /// Perform the per-event analysis
 68    void analyze(const Event& event) {
 69
 70      // Get baseline electrons, muons, and jets
 71      // NB. for electrons, we don't apply the loose ID here, since we don't want to double-count effs with later use of tight ID
 72      Particles elecs = apply<ParticleFinder>(event, "Electrons").particles();
 73      Particles muons = apply<ParticleFinder>(event, "Muons").particles();
 74      Jets jets = apply<JetFinder>(event, "Jets").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 2.4);
 75      iselect(jets, JetEffFilter([](const Jet& j) { return j.pT() > 60*GeV ? 1.0 : 0.94; }));
 76
 77
 78      // Jet/electron/muon overlap removal and selection
 79      // Remove any untagged jet within dR = 0.2 of an electron
 80      for (const Particle& e : elecs)
 81        idiscard(jets, [&](const Jet& j) { return !j.bTagged(Cuts::pT > 5*GeV) && deltaR(e, j, RAPIDITY) < 0.2; });
 82      // Remove any untagged low-multiplicity/muon-dominated jet within dR = 0.4 of a muon
 83      for (const Particle& m : muons)
 84        idiscard(jets, [&](const Jet& j) {
 85            if (j.bTagged(Cuts::pT > 5*GeV)) return false; /// @note A different b-tag working point, 85%, was actually used here *sigh*
 86            if (deltaR(m, j, RAPIDITY) > 0.4) return false;
 87            if (j.particles(Cuts::abscharge != 0).size() < 3) return true;
 88            return m.pT()/j.pT() > 0.5;
 89          });
 90      // Removing leptons within dR = 0.4 of remaining jets
 91      for (const Jet& j : jets) {
 92        idiscard(elecs, deltaRLess(j, 0.4, RAPIDITY));
 93        idiscard(muons, deltaRLess(j, 0.4, RAPIDITY));
 94      }
 95
 96      // Signal jet and lepton selection
 97      const Jets sigjets40 = select(jets, Cuts::pT > 40*GeV);
 98      const Jets sigjets60 = select(sigjets40, Cuts::pT > 60*GeV);
 99      const Jets sigbjets40 = select(sigjets40, [](const Jet& j) { return j.bTagged(Cuts::pT > 5*GeV); });
100      const Jets sigbjets60 = select(sigjets60, [](const Jet& j) { return j.bTagged(Cuts::pT > 5*GeV); });
101      const Particles sigmuons = select(muons, Cuts::pT > 35*GeV);
102      Particles sigelecs = select(elecs, Cuts::pT > 35*GeV);
103      iselect(sigelecs, ParticleEffFilter(ELECTRON_EFF_ATLAS_RUN2_TIGHT));
104
105
106      //////////////////
107
108
109      // Event selection cuts
110      if (sigelecs.size() + sigmuons.size() != 1) vetoEvent;
111      const Particle siglepton = sigelecs.empty() ? sigmuons.front() : sigelecs.front();
112
113      /// @note The note describes Nj = 5, 6, 7, 8, 9, >= 10 and Nb = 0, 1, 2, 3, >= 4 = 30 2D bins
114      ///  for each jet cut... but only provides data for six Nj = >= 8, 9, 10, Nb = 0, >= 3 bins.
115      /// We just implement the latter for now.
116
117      // Fill counters
118      if (sigjets40.size() >= 8  && sigbjets40.empty()) _h_08j40_0b->fill();
119      if (sigjets40.size() >= 9  && sigbjets40.empty()) _h_09j40_0b->fill();
120      if (sigjets40.size() >= 10 && sigbjets40.empty()) _h_10j40_0b->fill();
121      if (sigjets40.size() >= 8  && sigbjets40.size() >= 3) _h_08j40_3b->fill();
122      if (sigjets40.size() >= 9  && sigbjets40.size() >= 3) _h_09j40_3b->fill();
123      if (sigjets40.size() >= 10 && sigbjets40.size() >= 3) _h_10j40_3b->fill();
124
125      if (sigjets60.size() >= 8  && sigbjets60.empty()) _h_08j60_0b->fill();
126      if (sigjets60.size() >= 9  && sigbjets60.empty()) _h_09j60_0b->fill();
127      if (sigjets60.size() >= 10 && sigbjets60.empty()) _h_10j60_0b->fill();
128      if (sigjets60.size() >= 8  && sigbjets60.size() >= 3) _h_08j60_3b->fill();
129      if (sigjets60.size() >= 9  && sigbjets60.size() >= 3) _h_09j60_3b->fill();
130      if (sigjets60.size() >= 10 && sigbjets60.size() >= 3) _h_10j60_3b->fill();
131
132    }
133
134
135    /// Normalise counters after the run
136    void finalize() {
137
138      const double sf = 14.8*crossSection()/femtobarn/sumOfWeights();
139      scale(_h_08j40_0b, sf); scale(_h_09j40_0b, sf); scale(_h_10j40_0b, sf);
140      scale(_h_08j40_3b, sf); scale(_h_09j40_3b, sf); scale(_h_10j40_3b, sf);
141      scale(_h_08j60_0b, sf); scale(_h_09j60_0b, sf); scale(_h_10j60_0b, sf);
142      scale(_h_08j60_3b, sf); scale(_h_09j60_3b, sf); scale(_h_10j60_3b, sf);
143
144    }
145
146    /// @}
147
148
149  private:
150
151    /// @name Histograms
152    /// @{
153    CounterPtr _h_08j40_0b, _h_09j40_0b, _h_10j40_0b, _h_08j40_3b, _h_09j40_3b, _h_10j40_3b;
154    CounterPtr _h_08j60_0b, _h_09j60_0b, _h_10j60_0b, _h_08j60_3b, _h_09j60_3b, _h_10j60_3b;
155    /// @}
156
157
158  };
159
160
161
162  RIVET_DECLARE_PLUGIN(ATLAS_2016_CONF_2016_094);
163
164
165}