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 $t\bar{t}$ + jets at 13 TeV
 13  class ATLAS_2017_I1495243 : public Analysis {
 14  public:
 15
 16    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2017_I1495243);
 17
 18
 19    void init() {
 20
 21      Cut eta_full = Cuts::abseta < 5.0 && Cuts::pT > 1.0*MeV;
 22      Cut eta_lep = Cuts::abseta < 2.5;
 23
 24      // Collect final state particles
 25      FinalState FS(eta_full);
 26
 27      // Get photons to dress leptons
 28      IdentifiedFinalState photons(FS);
 29      photons.acceptIdPair(PID::PHOTON);
 30
 31      // Projection to find the electrons
 32      IdentifiedFinalState el_id(FS);
 33      el_id.acceptIdPair(PID::ELECTRON);
 34      PromptFinalState electrons(el_id);
 35      electrons.acceptTauDecays(false);
 36      LeptonFinder dressedelectrons(electrons, photons, 0.1, Cuts::abseta < 2.5 && Cuts::pT > 25*GeV);
 37      declare(dressedelectrons, "electrons");
 38      LeptonFinder fulldressedelectrons(electrons, photons, 0.1, eta_full);
 39
 40      // Projection to find the muons
 41      IdentifiedFinalState mu_id(FS);
 42      mu_id.acceptIdPair(PID::MUON);
 43      PromptFinalState muons(mu_id);
 44      muons.acceptTauDecays(false);
 45      LeptonFinder dressedmuons(muons, photons, 0.1, Cuts::abseta < 2.5 && Cuts::pT > 25*GeV);
 46      declare(dressedmuons, "muons");
 47      LeptonFinder fulldressedmuons(muons, photons, 0.1, eta_full);
 48
 49      // Projection to find neutrinos to exclude from jets
 50      IdentifiedFinalState nu_id;
 51      nu_id.acceptNeutrinos();
 52      PromptFinalState neutrinos(nu_id);
 53      neutrinos.acceptTauDecays(false);
 54
 55      // Jet clustering
 56      VetoedFinalState vfs;
 57      vfs.addVetoOnThisFinalState(fulldressedelectrons);
 58      vfs.addVetoOnThisFinalState(fulldressedmuons);
 59      vfs.addVetoOnThisFinalState(neutrinos);
 60      FastJets jets(vfs, JetAlg::ANTIKT, 0.4, JetMuons::ALL, JetInvisibles::DECAY);
 61      declare(jets, "jets");
 62
 63      // Book Histograms
 64      book(_h["bjet_pt"] , 5,1,1);
 65      book(_h["2bjet_pt"], 6,1,1);
 66      book(_h["ljet_pt"] , 7,1,1);
 67
 68      for (size_t i = 0; i < 4; ++i) {
 69        book(_d["njet"  + to_str(i)], i+1, 1, 1);
 70        book(_h["Q0"    + to_str(i)], "_Q0"    + to_str(i+ 7), refData((i>1?"d":"d0") + to_str(i+ 8) + "-x01-y01"));
 71        book(_h["MQ0"   + to_str(i)], "_MQ0"   + to_str(i+12), refData("d" + to_str(i+12) + "-x01-y01"));
 72        book(_h["Qsum"  + to_str(i)], "_Qsum"  + to_str(i+16), refData("d" + to_str(i+16) + "-x01-y01"));
 73        book(_h["MQsum" + to_str(i)], "_MQsum" + to_str(i+20), refData("d" + to_str(i+20) + "-x01-y01"));
 74        book(_s["gapFracQ0"    + to_str(i)],  8+i, 1 ,1);
 75        book(_s["gapFracMQ0"   + to_str(i)], 12+i, 1, 1);
 76        book(_s["gapFracQsum"  + to_str(i)], 16+i, 1, 1);
 77        book(_s["gapFracMQsum" + to_str(i)], 20+i, 1, 1);
 78      }
 79    }
 80
 81
 82    void analyze(const Event& event) {
 83      // Get the selected objects, using the projections.
 84      Jets all_jets = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5);
 85
 86      const DressedLeptons electrons = discard(apply<LeptonFinder>(event, "electrons").dressedLeptons(),
 87        [&](const DressedLepton &e) {
 88          return any(all_jets, deltaRLess(e, 0.4));
 89        });
 90
 91      const DressedLeptons muons = discard(apply<LeptonFinder>(event, "muons").dressedLeptons(),
 92        [&](const DressedLepton &m) {
 93          return any(all_jets, deltaRLess(m, 0.4));
 94        });
 95
 96      if (electrons.size() != 1 || muons.size() != 1)  vetoEvent;
 97      if (electrons[0].charge() == muons[0].charge())  vetoEvent;
 98
 99      Jets bjets, extrajets;
100      for (Jet j : all_jets) {
101        size_t b_tagged = j.bTags(Cuts::pT > 5*GeV).size();
102        if (bjets.size() < 2 && b_tagged)  bjets += j;
103        else  extrajets += j;
104      }
105      if (bjets.size() < 2)  vetoEvent;
106
107      double bjetpt = bjets[0].pt();
108      if (bjetpt > 250*GeV)  bjetpt = 275*GeV;
109      _h["bjet_pt"]->fill(bjetpt);
110
111      double b2jetpt = bjets[1].pt();
112      if (b2jetpt > 150*GeV)  b2jetpt = 175*GeV;
113      _h["2bjet_pt"]->fill(b2jetpt);
114
115      if (extrajets.size()) {
116        double ljetpt = extrajets[0].pt();
117        if (ljetpt > 250*GeV)  ljetpt = 275*GeV;
118        _h["ljet_pt"]->fill(ljetpt);
119      }
120
121      double Memubb = (electrons[0].momentum() + muons[0].momentum() + bjets[0].momentum() + bjets[1].momentum()).mass();
122      vector<double> leadpt = { 0., 0., 0., 0. }, ptsum = { 0., 0., 0., 0. };
123      vector<size_t> njetcount = { 0, 0, 0, 0 };
124      for (size_t i = 0; i < extrajets.size(); ++i) {
125        double absrap = extrajets[i].absrap(), pt = extrajets[i].pT();
126        if (pt > 25*GeV)  ++njetcount[0];
127        if (pt > 40*GeV)  ++njetcount[1];
128        if (pt > 60*GeV)  ++njetcount[2];
129        if (pt > 80*GeV)  ++njetcount[3];
130
131        if (absrap < 0.8 && pt > leadpt[0])  leadpt[0] = pt;
132        else if (absrap > 0.8 && absrap < 1.5 && pt > leadpt[1])  leadpt[1] = pt;
133        else if (absrap > 1.5 && absrap < 2.1 && pt > leadpt[2])  leadpt[2] = pt;
134        if (absrap < 2.1 && pt > leadpt[3])  leadpt[3] = pt;
135
136        if (absrap < 0.8)  ptsum[0] += pt;
137        else if (absrap > 0.8 && absrap < 1.5)  ptsum[1] += pt;
138        else if (absrap > 1.5 && absrap < 2.1)  ptsum[2] += pt;
139        if (absrap < 2.1)  ptsum[3] += pt;
140      }
141
142
143      for (size_t i = 0; i < 4; ++i) {
144        size_t cutoff = i? 3 : 4;
145        if (njetcount[i] > cutoff)  njetcount[i] = cutoff;
146        _d["njet" + to_str(i)]->fill(discretise(njetcount[i], i));
147
148        if (leadpt[i] > 305*GeV)  leadpt[i] = 305*GeV;
149        _h["Q0" + to_str(i)]->fill(leadpt[i]);
150
151        if (ptsum[i] > 505*GeV)  ptsum[i] = 505*GeV;
152        _h["Qsum" + to_str(i)]->fill(ptsum[i]);
153      }
154
155
156      for (size_t i = 0; i < 4; ++i) {
157        if (i == 0 && !(Memubb < 300*GeV))  continue;
158        if (i == 1 && !(Memubb > 300*GeV && Memubb < 425*GeV))  continue;
159        if (i == 2 && !(Memubb > 425*GeV && Memubb < 600*GeV))  continue;
160        if (i == 3 && !(Memubb > 600*GeV))  continue;
161        _h["MQ0"   + to_str(i)]->fill(leadpt[3]);
162        _h["MQsum" + to_str(i)]->fill(ptsum[3]);
163      }
164    }
165
166
167    void constructGapFraction(Estimate1DPtr out, Histo1DPtr in) {
168      bool hasWeights = in->effNumEntries() != in->numEntries();
169      double denW  = in->sumW();
170      double denW2 = in->sumW2();
171      size_t nEnd  = out->numBins();
172
173      for (auto& b : out->bins()) {
174          double numW = in->sumW(), numW2 = in->sumW2();
175          for (size_t j = b.index(); j <= nEnd; ++j) {
176            numW  -= in->bin(j).sumW();
177            numW2 -= in->bin(j).sumW2();
178          }
179          double yval = safediv(numW, denW);
180          double yerr = sqrt(safediv(yval * (1 - yval), denW));
181          if (hasWeights) { // use F. James's approximation for weighted events
182            yerr = sqrt( safediv((1 - 2 * yval) * numW2 + yval * yval * denW2, denW * denW) );
183          }
184          b.set(yval, yerr);
185      }
186    }
187
188
189    void finalize() {
190
191      // Build gap fraction plots
192      for (size_t i = 0; i < 4; ++i) {
193        constructGapFraction(_s["gapFracQ0"    + to_str(i)], _h["Q0"    + to_str(i)]);
194        constructGapFraction(_s["gapFracMQ0"   + to_str(i)], _h["MQ0"   + to_str(i)]);
195        constructGapFraction(_s["gapFracQsum"  + to_str(i)], _h["Qsum"  + to_str(i)]);
196        constructGapFraction(_s["gapFracMQsum" + to_str(i)], _h["MQsum" + to_str(i)]);
197      }
198
199      // Normalize to cross-section
200      for (map<string, Histo1DPtr>::iterator hit = _h.begin(); hit != _h.end(); ++hit) {
201        if (hit->first.find("jet") != string::npos)  normalize(hit->second);
202      }
203      normalize(_d);
204    }
205
206    string discretise(const size_t n, const size_t axis) const {
207      if (n == 0) return "0"s;
208      if (n == 1) return "1"s;
209      if (n == 2) return "2"s;
210      if (axis) {
211        return ">= 3"s;
212      }
213      else if (n == 3) return "3"s;
214      else if (n <= 8) return "4.0 - 8.0"s;
215      return "OTHER"s;
216    }
217
218
219  private:
220
221    /// @name Histogram helper functions
222    map<string, Histo1DPtr> _h;
223    map<string, Estimate1DPtr> _s;
224    map<string, BinnedHistoPtr<string>> _d;
225  };
226
227
228  RIVET_DECLARE_PLUGIN(ATLAS_2017_I1495243);
229
230
231}