1#include "Rivet/Analysis.hh"
  2#include "Rivet/Projections/FinalState.hh"
  3#include "Rivet/Projections/IdentifiedFinalState.hh"
  4#include "Rivet/Projections/VisibleFinalState.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6#include "Rivet/Projections/MissingMomentum.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  /// ATLAS 13 TeV monojet search with 3.2/fb of pp data
 15  class ATLAS_2016_I1452559 : public Analysis {
 16  public:
 17
 18    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2016_I1452559);
 19
 20    void init() {
 21
 22      FastJets jets(FinalState(Cuts::abseta < 4.9), JetAlg::ANTIKT, 0.4);
 23      SmearedJets recojets(jets, JET_SMEAR_ATLAS_RUN1);
 24      declare(recojets, "Jets");
 25
 26      FinalState electrons(Cuts::abspid == PID::ELECTRON && Cuts::abseta < 2.47 && Cuts::pT > 20*GeV);
 27      SmearedParticles recoelectrons(electrons, ELECTRON_EFF_ATLAS_RUN1_MEDIUM);
 28      declare(recoelectrons, "Electrons");
 29
 30      FinalState muons(Cuts::abspid == PID::MUON && Cuts::abseta < 2.50 && Cuts::pT > 10*GeV);
 31      SmearedParticles recomuons(muons, MUON_EFF_ATLAS_RUN1);
 32      declare(recomuons, "Muons");
 33
 34      VisibleFinalState calofs(Cuts::abseta < 4.9 && Cuts::abspid != PID::MUON);
 35      MissingMomentum met(calofs);
 36      SmearedMET recomet(met, MET_SMEAR_ATLAS_RUN1);
 37      declare(recomet, "MET");
 38
 39
 40      /// Book histograms
 41      for (size_t i = 0; i < 7; ++i)
 42        book(_count_IM[i], "count_IM" + toString(i+1));
 43      for (size_t i = 0; i < 6; ++i)
 44        book(_count_EM[i], "count_EM" + toString(i+1));
 45
 46    }
 47
 48
 49    void analyze(const Event& event) {
 50
 51      const Jets jets = apply<JetFinder>(event, "Jets").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 2.8);
 52      const Particles elecs = apply<ParticleFinder>(event, "Electrons").particlesByPt();
 53      const Particles mus = apply<ParticleFinder>(event, "Muons").particlesByPt();
 54      MSG_DEBUG("Number of raw jets, electrons, muons = "
 55                << jets.size() << ", " << elecs.size() << ", " << mus.size());
 56
 57      // Discard jets very close to electrons, or with low track multiplicity and close to muons
 58      const Jets isojets = discard(jets, [&](const Jet& j) {
 59          /// @todo Add track efficiency random filtering
 60          if (any(elecs, deltaRLess(j, 0.2))) return true;
 61          if (j.particles(Cuts::abscharge > 0 && Cuts::pT > 0.4*GeV).size() < 3 &&
 62              any(mus, deltaRLess(j, 0.4))) return true;
 63          return false;
 64        });
 65
 66      // Discard electrons close to remaining jets
 67      const Particles isoelecs = discard(elecs, [&](const Particle& e) {
 68          return any(isojets, deltaRLess(e, 0.4));
 69        });
 70
 71      // Discard muons close to remaining jets
 72      const Particles isomus = discard(mus, [&](const Particle& m) {
 73          for (const Jet& j : isojets) {
 74            if (deltaR(j,m) > 0.4) continue;
 75            if (j.particles(Cuts::abscharge > 0 && Cuts::pT > 0.4*GeV).size() > 3) return true;
 76          }
 77          return false;
 78        });
 79
 80      // Calculate ETmiss
 81      //const Vector3& vet = apply<MissingMomentum>(event, "MET").vectorEt();
 82      const Vector3& vet = apply<SmearedMET>(event, "MET").vectorEt();
 83      const double etmiss = vet.perp();
 84
 85
 86      // Event selection cuts
 87      if (etmiss < 250*GeV) vetoEvent;
 88      // Require at least one jet with pT > 250 GeV and |eta| < 2.4
 89      if (select(isojets, Cuts::pT > 250*GeV && Cuts::abseta < 2.4).empty()) vetoEvent;
 90      // Require at most 4 jets with pT > 30 GeV and |eta| < 2.8
 91      if (select(isojets, Cuts::pT > 30*GeV).size() > 4) vetoEvent;
 92      // Require no isolated jets within |dphi| < 0.4 of the MET vector
 93      if (any(isojets, deltaPhiLess(-vet, 0.4))) vetoEvent;
 94      // Require no isolated electrons or muons
 95      if (!isoelecs.empty() || !isomus.empty()) vetoEvent;
 96
 97
 98      ////////////////////
 99
100
101      // Get ETmiss bin number and fill counters
102      const int i_etmiss = binIndex(etmiss/GeV, ETMISS_CUTS);
103      // Inclusive ETmiss bins
104      for (int ibin = 0; ibin < 7; ++ibin)
105        if (i_etmiss >= ibin) _count_IM[ibin]->fill();
106      // Exclusive ETmiss bins
107      if (inRange(i_etmiss, 0, 6)) _count_EM[i_etmiss]->fill();
108
109    }
110
111
112    void finalize() {
113      const double norm = 3.2*crossSection()/femtobarn;
114      scale(_count_IM, norm/sumOfWeights());
115      scale(_count_EM, norm/sumOfWeights());
116    }
117
118
119  private:
120
121    const vector<double> ETMISS_CUTS = { 250, 300, 350, 400, 500, 600, 700, 13000 };
122    CounterPtr _count_IM[7], _count_EM[6];
123
124  };
125
126
127  RIVET_DECLARE_PLUGIN(ATLAS_2016_I1452559);
128
129}