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