1#include "Rivet/Analysis.hh"
  2#include "Rivet/Projections/FinalState.hh"
  3#include "Rivet/Projections/PromptFinalState.hh"
  4#include "Rivet/Projections/LeptonFinder.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// WWbb at 13 TeV
 11  class ATLAS_2018_I1677498 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2018_I1677498);
 16
 17
 18    /// Book cuts and projections
 19    void init() {
 20
 21      // All final state particles
 22      FinalState fs(Cuts::abseta < 5.0);
 23
 24      PromptFinalState photons(Cuts::abspid == PID::PHOTON, TauDecaysAs::PROMPT);
 25      PromptFinalState bare_el(Cuts::abspid == PID::ELECTRON, TauDecaysAs::PROMPT);
 26      LeptonFinder elecs(bare_el, photons, 0.1, Cuts::pT > 7*GeV && Cuts::abseta < 2.47);
 27      declare(elecs, "elecs");
 28
 29      PromptFinalState bare_mu(Cuts::abspid == PID::MUON, TauDecaysAs::PROMPT);
 30      LeptonFinder muons(bare_mu, photons, 0.1, Cuts::pT > 6*GeV && Cuts::abseta < 2.5);
 31      declare(muons, "muons");
 32
 33      FastJets jets(fs, JetAlg::ANTIKT, 0.4, JetMuons::NONE, JetInvisibles::NONE);
 34      declare(jets, "jets");
 35
 36      book(_h, 3, 1, 1);
 37    }
 38
 39
 40    void analyze(const Event& event) {
 41
 42      // Identify bjets (signal), light jets (OR) and soft b-jets (veto)
 43      size_t soft_bjets = 0;
 44      Jets bjets, lightjets;
 45      for (const Jet& jet : apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 5*GeV)) {
 46        bool isBjet = jet.bTagged(Cuts::pT > 5*GeV);
 47        if (isBjet) soft_bjets += 1;
 48        if (jet.abseta() < 2.5) {
 49          if ( isBjet && jet.pT() > 25*GeV) bjets += jet;
 50          if (!isBjet && jet.pT() > 20*GeV) lightjets += jet;
 51        }
 52      }
 53      if (soft_bjets != 2) vetoEvent;
 54      if (bjets.size() != 2) vetoEvent;
 55
 56      // Get dressed leptons
 57      DressedLeptons leptons;
 58      for (auto& lep : apply<LeptonFinder>(event, "muons").dressedLeptons()) { leptons.push_back(lep); }
 59      for (auto& lep : apply<LeptonFinder>(event, "elecs").dressedLeptons()) { leptons.push_back(lep); }
 60
 61      // 1. Find which light jets survive OR
 62      for (const auto& lep : leptons) {
 63        idiscard(lightjets, [&](const Jet& jet) {
 64          return deltaR(jet, lep) < 0.2 && (lep.abspid() == PID::ELECTRON || lep.pT()/jet.pT() > 0.7);
 65        });
 66      }
 67
 68      // 2. Find which leptons survive the OR and apply signal selection
 69      for (const auto& jet : (lightjets + bjets)) {
 70        idiscard(leptons, [&](const DressedLepton& lep) {
 71          return lep.pT() < 28*GeV || deltaR(jet, lep) < min(0.4, 0.04+10*GeV/lep.pT());
 72        });
 73      }
 74
 75      if (leptons.size() != 2) vetoEvent;
 76      std::sort(leptons.begin(), leptons.end(), cmpMomByPt);
 77
 78      // Z veto
 79      const size_t nEl = count(leptons, [](const DressedLepton& lep) { return  lep.abspid() == PID::ELECTRON; });
 80      const double mll = (leptons[0].mom() + leptons[1].mom()).mass();
 81      if (nEl != 1 && !(fabs(mll - 91*GeV) > 15*GeV && mll > 10*GeV)) vetoEvent;
 82
 83      const double m00 = (leptons[0].mom() + bjets[0].mom()).mass();
 84      const double m10 = (leptons[1].mom() + bjets[0].mom()).mass();
 85      const double m01 = (leptons[0].mom() + bjets[1].mom()).mass();
 86      const double m11 = (leptons[1].mom() + bjets[1].mom()).mass();
 87      double minimax = min( max(m00,m11), max(m10,m01) );
 88      minimax = min(minimax, 419.); // put overflow in the last bin
 89      _h->fill(minimax/GeV);
 90    }
 91
 92
 93    /// Finalise
 94    void finalize() {
 95      normalize(_h);
 96    }
 97
 98
 99  private:
100
101    /// Histogram
102    Histo1DPtr _h;
103
104  };
105
106
107  RIVET_DECLARE_PLUGIN(ATLAS_2018_I1677498);
108
109}