rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

CMS_2023_I2703254

Inclusive and differential cross section measurements of ttbb production in the lepton+jets channel at 13 TeV
Experiment: CMS (LHC)
Inspire ID: 2703254
Status: VALIDATED
Authors:
  • Seo Hyeon An
  • Florencia Canelli
  • Ji Eun Choi
  • Kyle Cormier
  • Tae Jeong Kim
  • Umberto Molinatti
  • Emanuel Pfeffer
  • Jan van der Linden
  • Sebastien Wertz
References: Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • ttbar events at sqrt(s) = 13 TeV

Measurement of inclusive and differential cross sections of the associated production of top quark and b quark pairs, $t\bar{t}b\bar{b}$. The cross sections are measured in the semileptonic decay channel of the top quark pair, using events containing exactly one isolated electron or muon. Four fiducial phase spaces are defined. "6j4b" requires the presence of at least six jets and at least four b jets. "5j3b" requires at least five jets and at least three b jets. "6j3b3l" requires at least six jets, including at least three b jets, and at least three light jets "7j4b3l" phase space, requiring at least seven jets, including at least four b jets and at least three light jets. All jets are required to have $p_\perp > 25 \text{GeV}$ and $|\eta| < 2.4$. B jets are defined using ghost-matching of B hadrons. No requirement is applied on the origin of of the b jets, i.e. the phase space definition is independent from simulated parton content.

Source code: CMS_2023_I2703254.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/LeptonFinder.hh"
  6#include "Rivet/Projections/PromptFinalState.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief ttbb in lepton+jets at 13 TeV
 12  class CMS_2023_I2703254 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2023_I2703254);
 17
 18
 19    /// @name Analysis methods
 20    ///@{
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24
 25      Cut eta_cut = (Cuts::abseta < 5.0);
 26      const FinalState fs(eta_cut);
 27
 28      PromptFinalState photons      (eta_cut && Cuts::abspid == PID::PHOTON, TauDecaysAs::PROMPT, MuDecaysAs::PROMPT);
 29      PromptFinalState electrons    (eta_cut && Cuts::abspid == PID::ELECTRON, TauDecaysAs::PROMPT, MuDecaysAs::PROMPT);
 30      PromptFinalState muons        (eta_cut && Cuts::abspid == PID::MUON, TauDecaysAs::PROMPT, MuDecaysAs::PROMPT);
 31
 32      Cut electron_cut     = (Cuts::abseta < 2.4) && (Cuts::pT > 29*GeV);
 33      Cut muon_cut         = (Cuts::abseta < 2.4) && (Cuts::pT > 26*GeV);
 34      Cut vetoelectron_cut = (Cuts::abseta < 2.5) && (Cuts::pT > 15*GeV);
 35      Cut vetomuon_cut     = (Cuts::abseta < 2.4) && (Cuts::pT > 15*GeV);
 36
 37      LeptonFinder dressedelectrons    (electrons, photons, 0.1, electron_cut);
 38      LeptonFinder dressedmuons        (muons,     photons, 0.1, muon_cut);
 39      LeptonFinder dressedvetoelectrons(electrons, photons, 0.1, vetoelectron_cut);
 40      LeptonFinder dressedvetomuons    (muons,     photons, 0.1, vetomuon_cut);
 41
 42      declare(dressedelectrons,     "elecs");
 43      declare(dressedmuons,         "muons");
 44      declare(dressedvetoelectrons, "vetoelecs");
 45      declare(dressedvetomuons,     "vetomuons");
 46
 47      FastJets jetfs(fs, JetAlg::ANTIKT, 0.4, JetMuons::ALL, JetInvisibles::NONE);
 48      declare(jetfs, "jets");
 49
 50      // Book histograms
 51      book(_xsec, 1, 1, 1);
 52      book(_h["average_DR_6j_4b"],            11, 1, 1);
 53      book(_h["bjet3_abseta_5j_3b"],          12, 1, 1);
 54      book(_h["bjet3_abseta_6j_4b"],          13, 1, 1);
 55      book(_h["bjet3_pt_5j_3b"],              14, 1, 1);
 56      book(_h["bjet3_pt_6j_4b"],              15, 1, 1);
 57      book(_h["bjet4_abseta_6j_4b"],          16, 1, 1);
 58      book(_h["bjet4_pt_6j_4b"],              17, 1, 1);
 59      book(_h["bJet_Ht_5j_3b"],               18, 1, 1);
 60      book(_h["bJet_Ht_6j_4b"],               19, 1, 1);
 61      book(_h["dPhi_lj_bsoft_6j_3b_3lj"],     20, 1, 1);
 62      book(_h["dPhi_lj_bsoft_7j_4b_3lj"],     21, 1, 1);
 63      book(_h["extra_jet1_abseta_6j_4b"],     22, 1, 1);
 64      book(_h["extra_jet1_pt_6j_4b"],         23, 1, 1);
 65      book(_h["extra_jet2_abseta_6j_4b"],     24, 1, 1);
 66      book(_h["extra_jet2_pt_6j_4b"],         25, 1, 1);
 67      book(_h["extra_jet_abseta_6j_4b"],      26, 1, 1);
 68      book(_h["extra_jet_DR_6j_4b"],          27, 1, 1);
 69      book(_h["extra_jet_M_6j_4b"],           28, 1, 1);
 70      book(_h["extra_jet_pt_6j_4b"],          29, 1, 1);
 71      book(_h["extra_lightJet_pt_6j_3b_3lj"], 30, 1, 1);
 72      book(_h["extra_lightJet_pt_7j_4b_3lj"], 31, 1, 1);
 73      book(_h["jet_Ht_5j_3b"],                32, 1, 1);
 74      book(_h["jet_Ht_6j_4b"],                33, 1, 1);
 75      book(_h["largest_Mbb_6j_4b"],           34, 1, 1);
 76      book(_h["lightJet_Ht_6j_3b_3lj"],       35, 1, 1);
 77      book(_h["lightJet_Ht_7j_4b_3lj"],       36, 1, 1);
 78      book(_h["Njets_5j_3b"],                 37, 1, 1);
 79      book(_h["Njets_6j_4b"],                 38, 1, 1);
 80      book(_h["n_M_tagged_jets_5j_3b"],       39, 1, 1);
 81    }
 82
 83
 84    /// Perform the per-event analysis
 85    void analyze(const Event& event) {
 86
 87      if (_edges.empty())  _edges = _xsec->xEdges();
 88
 89      Particles leptons;
 90      for (auto &lep : apply<LeptonFinder>(event, "elecs").dressedLeptons()) { leptons.push_back(lep); }
 91      for (auto &lep : apply<LeptonFinder>(event, "muons").dressedLeptons()) { leptons.push_back(lep); }
 92
 93      DressedLeptons Electrons     = apply<LeptonFinder>(event, "elecs").dressedLeptons();
 94      DressedLeptons Muons         = apply<LeptonFinder>(event, "muons").dressedLeptons();
 95      DressedLeptons vetoElectrons = apply<LeptonFinder>(event, "vetoelecs").dressedLeptons();
 96      DressedLeptons vetoMuons     = apply<LeptonFinder>(event, "vetomuons").dressedLeptons();
 97
 98
 99      // Retrieve clustered jets, sorted by pT, with a minimum pT cut
100      const Jets nonisojets = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.4);
101      const Jets jets = discardIfAnyDeltaRLess(nonisojets, leptons, 0.4);
102
103      Jets bjets;
104      Jets lfjets;
105
106      for (const Jet& jet : jets) {
107        if (jet.bTagged())  bjets += jet;
108        else lfjets += jet;
109      }
110
111      size_t njets   = jets.size();
112      size_t nbjets  = bjets.size();
113      size_t nlfjets = lfjets.size();
114
115      bool pass_ljets   = (Electrons.size() == 1 && vetoElectrons.size() == 1 && vetoMuons.size()     == 0)
116                           || (Muons.size() == 1 && vetoMuons.size()     == 1 && vetoElectrons.size() == 0);
117      if (!(pass_ljets))            vetoEvent;
118      if (nbjets < 3 || njets < 5)  vetoEvent;
119
120
121      double ht_jets   = sum(  jets, Kin::pT, 0.0);
122      double ht_bjets  = sum( bjets, Kin::pT, 0.0);
123      double ht_lfjets = sum(lfjets, Kin::pT, 0.0);
124
125      size_t ind1{0}, ind2{0}; double mindr = 999.;
126      double maxmbb  =   0.;
127      double sum_dr = 0.0;
128      size_t sum_n_dr = 0;
129      for (size_t i = 0; i < bjets.size(); ++i) {
130        for (size_t j = i + 1; j < bjets.size(); ++j) {
131          double dr = deltaR(bjets[i], bjets[j]);
132          double mbb = (bjets[i].momentum() + bjets[j].momentum()).mass();
133
134          if (dr < mindr) {
135            ind1 = i;
136            ind2 = j;
137            mindr = dr;
138          }
139          sum_dr += dr;
140          sum_n_dr += 1;
141
142          if (mbb > maxmbb) maxmbb = mbb;
143        }
144      }
145
146      FourMomentum bb_closest = bjets[ind1].momentum() + bjets[ind2].momentum();
147      double dr_closest = deltaR(bjets[ind1], bjets[ind2]);
148
149      // 6j3b3l and 7j4b3l:
150      double gen_w_mass = 79.6;
151      size_t lfind1{0}, lfind2{0}; double minwmasdiff = 9999.;
152      size_t lfind = 0;
153      double phi_first_ljet_extra_bjet_soft = 0.;
154      if (nlfjets >= 3) {
155        for (size_t i = 0; i < nlfjets; ++i) {
156          for (size_t j = i + 1; j < nlfjets; ++j) {
157            double wmasdiff = abs(gen_w_mass - (lfjets[i].momentum()+lfjets[j].momentum()).mass());
158            if (wmasdiff < minwmasdiff) {
159              lfind1 = i;
160              lfind2 = j;
161              minwmasdiff = wmasdiff;
162            }
163          }
164        }
165        for (size_t i = 0; i < 3; ++i){
166            if (!(( i == lfind1) || ( i == lfind2))){
167            lfind = i;
168            break;
169            }
170        }
171        phi_first_ljet_extra_bjet_soft = deltaPhi(lfjets[lfind].momentum(),bjets[nbjets - 1].momentum());
172      }
173
174      // Fill histograms
175      // 5j3b category
176      if (njets >= 5 && nbjets >= 3) {
177        _xsec->fill(_edges[3]);
178        _h["Njets_5j_3b"]                   ->fill(njets);
179        _h["n_M_tagged_jets_5j_3b"]         ->fill(nbjets);
180        _h["bjet3_pt_5j_3b"]                ->fill(bjets[2].pT()/GeV);
181        _h["bjet3_abseta_5j_3b"]            ->fill(bjets[2].abseta()/GeV);
182        _h["jet_Ht_5j_3b"]                  ->fill(ht_jets/GeV);
183        _h["bJet_Ht_5j_3b"]                 ->fill(ht_bjets/GeV);
184      }
185
186      // 6j4b category
187      if (njets >= 6 && nbjets >= 4) {
188        _xsec->fill(_edges[1]);
189        _h["Njets_6j_4b"]                   ->fill(njets);
190        _h["bjet3_pt_6j_4b"]                ->fill(bjets[2].pT()/GeV);
191        _h["bjet3_abseta_6j_4b"]            ->fill(bjets[2].abseta()/GeV);
192        _h["bjet4_pt_6j_4b"]                ->fill(bjets[3].pT()/GeV);
193        _h["bjet4_abseta_6j_4b"]            ->fill(bjets[3].abseta()/GeV);
194        _h["jet_Ht_6j_4b"]                  ->fill(ht_jets/GeV);
195        _h["bJet_Ht_6j_4b"]                 ->fill(ht_bjets/GeV);
196        _h["average_DR_6j_4b"]              ->fill(sum_dr/sum_n_dr);
197        _h["largest_Mbb_6j_4b"]             ->fill(maxmbb);
198        _h["extra_jet_DR_6j_4b"]            ->fill(dr_closest);
199        _h["extra_jet_M_6j_4b"]             ->fill(bb_closest.mass()/GeV);
200        _h["extra_jet_pt_6j_4b"]            ->fill(bb_closest.pT()/GeV);
201        _h["extra_jet_abseta_6j_4b"]        ->fill(bb_closest.abseta()/GeV);
202        _h["extra_jet1_pt_6j_4b"]           ->fill(bjets[ind1].pT()/GeV);
203        _h["extra_jet1_abseta_6j_4b"]       ->fill(bjets[ind1].abseta()/GeV);
204        _h["extra_jet2_pt_6j_4b"]           ->fill(bjets[ind2].pT()/GeV);
205        _h["extra_jet2_abseta_6j_4b"]       ->fill(bjets[ind2].abseta()/GeV);
206      }
207      // 6j3b3l category
208      if (njets >= 6 && nbjets >= 3 && nlfjets >= 3) {
209        _xsec->fill(_edges[2]);
210        _h["extra_lightJet_pt_6j_3b_3lj"]   ->fill(lfjets[lfind].pT()/GeV);
211        _h["dPhi_lj_bsoft_6j_3b_3lj"]       ->fill(phi_first_ljet_extra_bjet_soft);
212        _h["lightJet_Ht_6j_3b_3lj"]         ->fill(ht_lfjets/GeV);
213      }
214
215      // 7j4b3l category
216      if (njets >= 7 && nbjets >= 4 && nlfjets >= 3) {
217        _xsec->fill(_edges[0]);
218        _h["extra_lightJet_pt_7j_4b_3lj"]   ->fill(lfjets[lfind].pT()/GeV);
219        _h["dPhi_lj_bsoft_7j_4b_3lj"]       ->fill(phi_first_ljet_extra_bjet_soft);
220        _h["lightJet_Ht_7j_4b_3lj"]         ->fill(ht_lfjets/GeV);
221      }
222    }
223
224
225    /// Normalise histograms etc., after the run
226    void finalize() {
227      const double sf = crossSection() / femtobarn / sumOfWeights();
228
229      // move the overflow into the last "true" bin of the distribution _just_ for the following
230      const std::list<string> overflows =
231      {
232        "Njets_5j_3b",
233        "bJet_Ht_5j_3b",
234        "bjet3_pt_5j_3b",
235        "jet_Ht_5j_3b",
236        "n_M_tagged_jets_5j_3b",
237        "extra_lightJet_pt_6j_3b_3lj",
238        "lightJet_Ht_6j_3b_3lj",
239        "Njets_6j_4b",
240        "bJet_Ht_6j_4b",
241        "bjet3_pt_6j_4b",
242        "bjet4_pt_6j_4b",
243        "extra_jet1_pt_6j_4b",
244        "extra_jet2_pt_6j_4b",
245        "extra_jet_M_6j_4b",
246        "extra_jet_abseta_6j_4b",
247        "extra_jet_pt_6j_4b",
248        "jet_Ht_6j_4b",
249        "largest_Mbb_6j_4b",
250        "extra_lightJet_pt_7j_4b_3lj",
251        "lightJet_Ht_7j_4b_3lj"
252      };
253      for (const string& of : overflows) {
254        auto& of_bin = _h[of]->bin(_h[of]->numBins()+1);
255        _h[of]->bin(_h[of]->numBins()) += of_bin;
256        of_bin.reset();
257      }
258
259      scale(_xsec, sf);
260      normalize(_h);
261
262    }
263
264    ///@}
265
266
267    /// @name Histograms
268    ///@{
269    map<string, Histo1DPtr> _h;
270    BinnedHistoPtr<string> _xsec;
271    vector<string> _edges;
272    ///@}
273
274
275  };
276
277
278  RIVET_DECLARE_PLUGIN(CMS_2023_I2703254);
279
280}