1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/VetoedFinalState.hh"
  4#include "Rivet/Projections/IdentifiedFinalState.hh"
  5#include "Rivet/Projections/PromptFinalState.hh"
  6#include "Rivet/Projections/LeptonFinder.hh"
  7#include "Rivet/Projections/FastJets.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief Lepton differential ttbar analysis at 8 TeV
 13  class ATLAS_2017_I1626105 : public Analysis {
 14  public:
 15
 16    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2017_I1626105);
 17
 18    void init() {
 19
 20      Cut eta_full = Cuts::abseta < 5.0 && Cuts::pT > 1.0*MeV;
 21
 22      // All final state particles
 23      const FinalState fs;
 24
 25      // Get photons to dress leptons
 26      IdentifiedFinalState photons(fs, PID::PHOTON);
 27      PromptFinalState prompt_photons(photons);
 28
 29      // Projection to find the electrons
 30      PromptFinalState prompt_el(Cuts::abspid == PID::ELECTRON, TauDecaysAs::PROMPT);
 31      LeptonFinder elecs(prompt_el, photons, 0.1, Cuts::abseta < 2.5 && Cuts::pT > 25*GeV);
 32      LeptonFinder veto_elecs(prompt_el, prompt_photons, 0.1, eta_full);
 33      declare(elecs, "elecs");
 34
 35      // Projection to find the muons
 36      PromptFinalState prompt_mu(Cuts::abspid == PID::MUON, TauDecaysAs::PROMPT);
 37      LeptonFinder muons(prompt_mu, photons, 0.1, Cuts::abseta < 2.5 && Cuts::pT > 25*GeV);
 38      LeptonFinder veto_muons(prompt_mu, prompt_photons, 0.1, eta_full);
 39      declare(muons, "muons");
 40
 41      // Jet clustering.
 42      VetoedFinalState vfs;
 43      vfs.addVetoOnThisFinalState(veto_elecs);
 44      vfs.addVetoOnThisFinalState(veto_muons);
 45      FastJets jets(vfs, JetAlg::ANTIKT, 0.4);
 46      jets.useInvisibles();
 47      declare(jets, "jets");
 48
 49      // Book histograms
 50      bookHistos("lep_pt",       1);
 51      bookHistos("lep_eta",      3);
 52      bookHistos("dilep_pt",     5);
 53      bookHistos("dilep_mass",   7);
 54      bookHistos("dilep_rap",    9);
 55      bookHistos("dilep_dphi",  11);
 56      bookHistos("dilep_sumpt", 13);
 57      bookHistos("dilep_sumE",  15);
 58    }
 59
 60
 61    void analyze(const Event& event) {
 62      DressedLeptons elecs = sortByPt(apply<LeptonFinder>(event, "elecs").dressedLeptons());
 63      DressedLeptons muons = sortByPt(apply<LeptonFinder>(event, "muons").dressedLeptons());
 64      Jets jets = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5);
 65
 66      // Check overlap of jets/leptons.
 67      for (const Jet& jet : jets) {
 68        idiscard(elecs, deltaRLess(jet, 0.4));
 69        idiscard(muons, deltaRLess(jet, 0.4));
 70      }
 71      if (elecs.empty() || muons.empty()) vetoEvent;
 72      if (elecs[0].charge() == muons[0].charge()) vetoEvent;
 73
 74      FourMomentum el = elecs[0].momentum();
 75      FourMomentum mu = muons[0].momentum();
 76      FourMomentum ll = elecs[0].momentum() + muons[0].momentum();
 77
 78      // Fill histograms
 79      fillHistos("lep_pt",      el.pT()/GeV);
 80      fillHistos("lep_pt",      mu.pT()/GeV);
 81      fillHistos("lep_eta",     el.abseta());
 82      fillHistos("lep_eta",     mu.abseta());
 83      fillHistos("dilep_pt",    ll.pT()/GeV);
 84      fillHistos("dilep_mass",  ll.mass()/GeV);
 85      fillHistos("dilep_rap",   ll.absrap());
 86      fillHistos("dilep_dphi",  deltaPhi(el, mu));
 87      fillHistos("dilep_sumpt", (el.pT() + mu.pT())/GeV);
 88      fillHistos("dilep_sumE",  (el.E() + mu.E())/GeV);
 89    }
 90
 91
 92    void finalize() {
 93      // Normalize to cross-section
 94      const double sf = crossSection()/femtobarn/sumOfWeights();
 95      for (auto& hist : _h) {
 96        const double norm = 1.0 / hist.second->integral();
 97        // add overflow to last bin
 98        const size_t nBins = hist.second->numBins();
 99        const double overflow = hist.second->bin(nBins+1).effNumEntries();
100        hist.second->fill(hist.second->bin(nBins).xMid(), overflow);
101        // histogram normalisation
102        if (hist.first.find("norm") != string::npos)  scale(hist.second, norm);
103        else  scale(hist.second, sf);
104      }
105
106    }
107
108
109  private:
110
111    /// @name Histogram helper functions
112    /// @{
113    void bookHistos(const std::string name, unsigned int index) {
114      book(_h[name], index, 1, 1);
115      book(_h["norm_" + name], index + 1, 1, 1);
116    }
117
118    void fillHistos(const std::string name, double value) {
119      _h[name]->fill(value);
120      _h["norm_" + name]->fill(value);
121    }
122
123    map<string, Histo1DPtr> _h;
124    /// @}
125
126  };
127
128
129  RIVET_DECLARE_PLUGIN(ATLAS_2017_I1626105);
130
131}