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