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