1#include "Rivet/Analysis.hh"
  2#include "Rivet/Projections/VetoedFinalState.hh"
  3#include "Rivet/Projections/PromptFinalState.hh"
  4#include "Rivet/Projections/InvisibleFinalState.hh"
  5#include "Rivet/Projections/LeptonFinder.hh"
  6#include "Rivet/Projections/FastJets.hh"
  7#include "Rivet/Tools/HistoGroup.hh"
  8
  9namespace Rivet {
 10
 11  /// @brief lepton differential ttbar analysis at 13 TeV
 12  class ATLAS_2023_I2648096 : public Analysis {
 13  public:
 14
 15    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2023_I2648096);
 16
 17    void init() {
 18
 19      Cut eta_full = Cuts::abseta < 5.0 && Cuts::pT > 1.0*MeV;
 20
 21      // Get photons to dress leptons
 22      FinalState photons(Cuts::abspid == PID::PHOTON);
 23
 24      // Projection to find the electrons
 25      PromptFinalState prompt_el(Cuts::abspid == PID::ELECTRON, TauDecaysAs::PROMPT);
 26      LeptonFinder elecs(prompt_el, photons, 0.1, (Cuts::abseta < 2.5 && Cuts::pT > 25*GeV));
 27      LeptonFinder veto_elecs(prompt_el, photons, 0.1, eta_full);
 28      declare(elecs, "elecs");
 29
 30      // Projection to find the muons
 31      PromptFinalState prompt_mu(Cuts::abspid == PID::MUON, TauDecaysAs::PROMPT);
 32      LeptonFinder muons(prompt_mu, photons, 0.1, (Cuts::abseta < 2.5 && Cuts::pT > 25*GeV));
 33      LeptonFinder veto_muons(prompt_mu, photons, 0.1, eta_full);
 34      declare(muons, "muons");
 35
 36      const InvisibleFinalState invis(OnlyPrompt::YES, TauDecaysAs::PROMPT);
 37      VetoedFinalState vfs;
 38      vfs.addVetoOnThisFinalState(veto_elecs);
 39      vfs.addVetoOnThisFinalState(veto_muons);
 40      vfs.addVetoOnThisFinalState(invis);
 41      FastJets jets(vfs, JetAlg::ANTIKT, 0.4, JetMuons::ALL, JetInvisibles::ALL);
 42      declare(jets, "jets");
 43
 44      // Book histograms
 45      bookHistos("lep_pt",       6);
 46      bookHistos("lep_eta",      9);
 47      bookHistos("dilep_sumE",  12);
 48      bookHistos("dilep_mass",  15);
 49      bookHistos("dilep_sumpt", 18);
 50      bookHistos("dilep_pt",    21);
 51      bookHistos("dilep_dphi",  24);
 52      bookHistos("dilep_rap",   27);
 53
 54      // unrolled 2D distributions - 2nd-dim bin edges must be specified
 55      const vector<double> b_mass_2D = { 0., 70., 100., 130., 200., 800. };
 56      const vector<double> b_pt_2D = { 0., 40., 65., 100. };
 57      const vector<double> b_sumE_2D = { 0., 110., 140., 200., 250., 900. };
 58
 59      bookHisto2D("dilep_rap_mass",  78, b_mass_2D);
 60      bookHisto2D("dilep_dphi_mass", 79, b_mass_2D);
 61      bookHisto2D("dilep_dphi_pt",   80, b_pt_2D);
 62      bookHisto2D("dilep_dphi_sumE", 81, b_sumE_2D);
 63    }
 64
 65    void analyze(const Event& event) {
 66      DressedLeptons elecs = apply<LeptonFinder>(event, "elecs").dressedLeptons();
 67      DressedLeptons muons = apply<LeptonFinder>(event, "muons").dressedLeptons();
 68      Jets jets = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5);
 69
 70      // Check overlap of jets/leptons.
 71      for (const Jet& jet : jets) {
 72        idiscard(elecs, deltaRLess(jet, 0.4));
 73        idiscard(muons, deltaRLess(jet, 0.4));
 74      }
 75      if (elecs.empty() || muons.empty())  vetoEvent;
 76      if (elecs[0].charge() == muons[0].charge())  vetoEvent;
 77
 78      FourMomentum el = elecs[0].momentum();
 79      FourMomentum mu = muons[0].momentum();
 80      FourMomentum ll = elecs[0].momentum() + muons[0].momentum();
 81
 82      if (max(el.pT(), mu.pT()) < 27*GeV)  vetoEvent;
 83
 84      const double dphi = deltaPhi(el,mu);
 85      const double drap = ll.absrap();
 86      const double pt   = ll.pT()/GeV;
 87      const double sumE = min((el.E()+mu.E())/GeV, 899);
 88      const double mass = min(ll.mass()/GeV, 799.);
 89
 90      // Fill histograms
 91      // include explicit overflow protection as last bins are inclusive
 92      fillHistos("lep_pt",      min(el.pT()/GeV, 349.));
 93      fillHistos("lep_pt",      min(mu.pT()/GeV, 349.));
 94      fillHistos("lep_eta",     el.abseta());
 95      fillHistos("lep_eta",     mu.abseta());
 96      fillHistos("dilep_pt",    min(pt, 299.));
 97      fillHistos("dilep_mass",  mass);
 98      fillHistos("dilep_rap",   drap);
 99      fillHistos("dilep_dphi",  dphi);
100      fillHistos("dilep_sumpt", min((el.pT()+mu.pT())/GeV, 599));
101      fillHistos("dilep_sumE",  sumE);
102
103      // Fill unrolled 2D histograms vs mass
104      fillHisto2D("dilep_rap_mass",  mass, drap);
105      fillHisto2D("dilep_dphi_mass", mass, dphi);
106      fillHisto2D("dilep_dphi_pt",   min(pt, 99.), dphi);
107      fillHisto2D("dilep_dphi_sumE", sumE, dphi);
108    }
109
110    void finalize() {
111      // Normalize to cross-section
112      const double sf = crossSection()/femtobarn/sumOfWeights();
113
114      // finalisation of 1D histograms
115      for (auto& hist : _h) {
116        const double norm = 1.0 / hist.second->integral();
117        // histogram normalisation
118        if (hist.first.find("norm") != string::npos)  scale(hist.second, norm);
119        else  scale(hist.second, sf);
120      }
121
122      // finalisation of 2D histograms
123      for (auto& hist : _h_multi) {
124        if (hist.first.find("_norm") != std::string::npos) {
125          const double sf = safediv(1.0, hist.second->integral(false));
126          scale(hist.second, sf);
127        }
128        else {
129          scale(hist.second, sf);
130        }
131        divByGroupWidth(hist.second);
132      }
133    }
134
135
136  private:
137
138    /// @name Histogram helper functions
139    //@{
140    void bookHistos(const string& name, const unsigned int id) {
141      book(_h[name], id, 1, 1);
142      book(_h["norm_" + name], id+36, 1, 1);
143    }
144
145    void fillHistos(const string& name, const double value) {
146      _h[name]->fill(value);
147      _h["norm_" + name]->fill(value);
148    }
149
150    void bookHisto2D(const string& name, const unsigned int id, const vector<double>& axis1) {
151      book(_h_multi[name], axis1);
152      book(_h_multi[name + "_norm"], axis1);
153      for (size_t i = 1; i < axis1.size(); ++i) {
154        book(_h_multi[name]->bin(i), id, 1, i);
155        book(_h_multi[name + "_norm"]->bin(i), id+4, 1, i);
156      }
157    }
158
159
160
161    void fillHisto2D(const string& name, const double val1, const double val2) {
162      _h_multi[name]->fill(val1, val2);
163      _h_multi[name+"_norm"]->fill(val1, val2);
164    }
165
166
167    // pointers to 1D and 2D histograms
168    map<string, Histo1DPtr> _h;
169    map<string, Histo1DGroupPtr> _h_multi;
170    //@}
171    // acceptance counter
172
173  };
174
175  // Declare the class as a hook for the plugin system
176  RIVET_DECLARE_PLUGIN(ATLAS_2023_I2648096);
177}