1#include "Rivet/Analysis.hh"
  2#include "Rivet/Projections/VetoedFinalState.hh"
  3#include "Rivet/Projections/IdentifiedFinalState.hh"
  4#include "Rivet/Projections/PromptFinalState.hh"
  5#include "Rivet/Projections/LeptonFinder.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief Lepton differential ttbar analysis at 13 TeV
 11  class ATLAS_2019_I1759875 : public Analysis {
 12  public:
 13
 14    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2019_I1759875);
 15
 16    void init() {
 17
 18      Cut eta_full = Cuts::abseta < 5.0 && Cuts::pT > 1.0*MeV;
 19
 20      // Get photons to dress leptons
 21      PromptFinalState photons(Cuts::pid == PID::PHOTON);
 22
 23      // Projection to find the electrons
 24      PromptFinalState prompt_el(Cuts::abspid == PID::ELECTRON, TauDecaysAs::PROMPT);
 25      LeptonFinder elecs(prompt_el, photons, 0.1, Cuts::abseta < 2.5 && Cuts::pT > 20*GeV);
 26      LeptonFinder veto_elecs(prompt_el, photons, 0.1, eta_full);
 27      declare(elecs, "elecs");
 28
 29      // Projection to find the muons
 30      PromptFinalState prompt_mu(Cuts::abspid == PID::MUON, TauDecaysAs::PROMPT);
 31      LeptonFinder muons(prompt_mu, photons, 0.1, Cuts::abseta < 2.5 && Cuts::pT > 20*GeV);
 32      LeptonFinder veto_muons(prompt_mu, photons, 0.1, eta_full);
 33      declare(muons, "muons");
 34
 35      VetoedFinalState vfs;
 36      vfs.addVetoOnThisFinalState(veto_elecs);
 37      vfs.addVetoOnThisFinalState(veto_muons);
 38
 39      // Book histograms
 40      bookHistos("lep_pt",       1);
 41      bookHistos("lep_eta",      3);
 42      bookHistos("dilep_pt",     5);
 43      bookHistos("dilep_mass",   7);
 44      bookHistos("dilep_rap",    9);
 45      bookHistos("dilep_dphi",  11);
 46      bookHistos("dilep_sumpt", 13);
 47      bookHistos("dilep_sumE",  15);
 48
 49      // unrolled 2D distributions - 2nd-dim bin edges must be specified
 50      std::vector<double> massbins={0.,80.,120.,200.,500.};
 51
 52      bookHisto2D("lep_eta_mass",17,massbins);
 53      bookHisto2D("dilep_rap_mass",19,massbins);
 54      bookHisto2D("dilep_dphi_mass",21,massbins);
 55    }
 56
 57    void analyze(const Event& event) {
 58      DressedLeptons elecs = apply<LeptonFinder>(event, "elecs").dressedLeptons();
 59      DressedLeptons muons = apply<LeptonFinder>(event, "muons").dressedLeptons();
 60
 61      if (elecs.empty() || muons.empty())  vetoEvent;
 62      if (elecs[0].charge() == muons[0].charge())  vetoEvent;
 63
 64      FourMomentum el = elecs[0].momentum();
 65      FourMomentum mu = muons[0].momentum();
 66      FourMomentum ll = elecs[0].momentum() + muons[0].momentum();
 67
 68      // Fill histograms
 69      // include explicit overflow protection as last bins are inclusive
 70      fillHistos("lep_pt",      min(el.pT()/GeV,299.));
 71      fillHistos("lep_pt",      min(mu.pT()/GeV,299.));
 72      fillHistos("lep_eta",     el.abseta());
 73      fillHistos("lep_eta",     mu.abseta());
 74      fillHistos("dilep_pt",    min(ll.pT()/GeV,299.));
 75      fillHistos("dilep_mass",  min(ll.mass()/GeV,499.));
 76      fillHistos("dilep_rap",   ll.absrap());
 77      fillHistos("dilep_dphi",  deltaPhi(el, mu));
 78      fillHistos("dilep_sumpt", min((el.pT()+mu.pT())/GeV,399.));
 79      fillHistos("dilep_sumE",  min((el.E()+mu.E())/GeV,699.));
 80
 81      // find mass bin variable, with overflow protection
 82      float massv=ll.mass()/GeV;
 83      if (massv>499.) massv=499.;
 84      // Fill unrolled 2D histograms vs mass
 85      fillHisto2D("lep_eta_mass",el.abseta(),massv);
 86      fillHisto2D("lep_eta_mass",mu.abseta(),massv);
 87      fillHisto2D("dilep_rap_mass",ll.absrap(),massv);
 88      fillHisto2D("dilep_dphi_mass",deltaPhi(el,mu),massv);
 89    }
 90
 91    void finalize() {
 92      // Normalize to cross-section
 93      const double sf = crossSection()/femtobarn/sumOfWeights();
 94
 95      // finalisation of 1D histograms
 96      for (auto& hist : _h) {
 97        const double norm = 1.0 / hist.second->integral();
 98        // histogram normalisation
 99        if (hist.first.find("norm") != string::npos)  scale(hist.second, norm);
100        else  scale(hist.second, sf);
101      }
102
103      // finalisation of 2D histograms
104      for (auto& hist : _h_multi) {
105        if (hist.first.find("_norm") != std::string::npos) {
106          // scaling for normalised distribution according integral of whole set
107          hist.second->normalizeGroup(1.0, false);
108        }
109        else {
110          // scaling for non-normalised distribution
111          scale(hist.second, sf);
112        }
113      }
114      divByGroupWidth(_h_multi);
115    }
116
117
118  private:
119
120    /// @name Histogram helper functions
121    /// @{
122    void bookHistos(const std::string name, unsigned int index) {
123      book(_h[name], index, 1, 1);
124      book(_h["norm_" + name],index + 1, 1, 1);
125    }
126
127    void fillHistos(const std::string name, double value) {
128      _h[name]->fill(value);
129      _h["norm_" + name]->fill(value);
130    }
131
132    void bookHisto2D(const std::string& name, unsigned int index, const std::vector<double>& massbins) {
133      book(_h_multi[name], massbins);
134      book(_h_multi[name+"_norm"], massbins);
135      for (size_t i=1; i < _h_multi[name]->numBins()+1; ++i) {
136        book(_h_multi[name]->bin(i), index, 1, i);
137        book(_h_multi[name+"_norm"]->bin(i), index+1, 1, i);
138      }
139    }
140
141
142
143    void fillHisto2D(const std::string& name, double val, double massval) {
144      _h_multi[name]->fill(massval, val);
145      _h_multi[name+"_norm"]->fill(massval, val);
146    }
147
148
149    // pointers to 1D and 2D histograms
150    map<string, Histo1DPtr> _h;
151    map<string, Histo1DGroupPtr> _h_multi;
152    /// @}
153    // acceptance counter
154
155  };
156
157  RIVET_DECLARE_PLUGIN(ATLAS_2019_I1759875);
158}