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: I2703254
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
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/DressedLeptons.hh"
#include "Rivet/Projections/PromptFinalState.hh"
#include "Rivet/Projections/JetAlg.hh"

namespace Rivet {


  /// @brief ttbb in lepton+jets at 13 TeV
  class CMS_2023_I2703254 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2023_I2703254);


    /// @name Analysis methods
    ///@{

    /// Book histograms and initialise projections before the run
    void init() {

      Cut eta_cut = (Cuts::abseta < 5.0);
      const FinalState fs(eta_cut);

      PromptFinalState photons      (eta_cut && Cuts::abspid == PID::PHOTON, true, true);
      PromptFinalState electrons    (eta_cut && Cuts::abspid == PID::ELECTRON, true, true);
      PromptFinalState muons        (eta_cut && Cuts::abspid == PID::MUON, true, true);

      Cut electron_cut     = (Cuts::abseta < 2.4) && (Cuts::pT > 29*GeV);
      Cut muon_cut         = (Cuts::abseta < 2.4) && (Cuts::pT > 26*GeV);
      Cut vetoelectron_cut = (Cuts::abseta < 2.5) && (Cuts::pT > 15*GeV);
      Cut vetomuon_cut     = (Cuts::abseta < 2.4) && (Cuts::pT > 15*GeV);

      DressedLeptons dressedelectrons     (photons, electrons,     0.1, electron_cut,     true);
      DressedLeptons dressedmuons         (photons, muons,         0.1, muon_cut,         true);
      DressedLeptons dressedvetoelectrons (photons, electrons,     0.1, vetoelectron_cut, true);
      DressedLeptons dressedvetomuons     (photons, muons,         0.1, vetomuon_cut,     true);

      declare(dressedelectrons,     "elecs");
      declare(dressedmuons,         "muons");
      declare(dressedvetoelectrons, "vetoelecs");
      declare(dressedvetomuons,     "vetomuons");

      FastJets jetfs(fs, FastJets::ANTIKT, 0.4, JetAlg::Muons::ALL, JetAlg::Invisibles::NONE);
      declare(jetfs, "jets");

      // Book histograms
      book(_h["xsec"],                         1, 1, 1);
      book(_h["average_DR_6j_4b"],            11, 1, 1);
      book(_h["bjet3_abseta_5j_3b"],          12, 1, 1);
      book(_h["bjet3_abseta_6j_4b"],          13, 1, 1);
      book(_h["bjet3_pt_5j_3b"],              14, 1, 1);
      book(_h["bjet3_pt_6j_4b"],              15, 1, 1);
      book(_h["bjet4_abseta_6j_4b"],          16, 1, 1);
      book(_h["bjet4_pt_6j_4b"],              17, 1, 1);
      book(_h["bJet_Ht_5j_3b"],               18, 1, 1);
      book(_h["bJet_Ht_6j_4b"],               19, 1, 1);
      book(_h["dPhi_lj_bsoft_6j_3b_3lj"],     20, 1, 1);
      book(_h["dPhi_lj_bsoft_7j_4b_3lj"],     21, 1, 1);
      book(_h["extra_jet1_abseta_6j_4b"],     22, 1, 1);
      book(_h["extra_jet1_pt_6j_4b"],         23, 1, 1);
      book(_h["extra_jet2_abseta_6j_4b"],     24, 1, 1);
      book(_h["extra_jet2_pt_6j_4b"],         25, 1, 1);
      book(_h["extra_jet_abseta_6j_4b"],      26, 1, 1);
      book(_h["extra_jet_DR_6j_4b"],          27, 1, 1);
      book(_h["extra_jet_M_6j_4b"],           28, 1, 1);
      book(_h["extra_jet_pt_6j_4b"],          29, 1, 1);
      book(_h["extra_lightJet_pt_6j_3b_3lj"], 30, 1, 1);
      book(_h["extra_lightJet_pt_7j_4b_3lj"], 31, 1, 1);
      book(_h["jet_Ht_5j_3b"],                32, 1, 1);
      book(_h["jet_Ht_6j_4b"],                33, 1, 1);
      book(_h["largest_Mbb_6j_4b"],           34, 1, 1);
      book(_h["lightJet_Ht_6j_3b_3lj"],       35, 1, 1);
      book(_h["lightJet_Ht_7j_4b_3lj"],       36, 1, 1);
      book(_h["Njets_5j_3b"],                 37, 1, 1);
      book(_h["Njets_6j_4b"],                 38, 1, 1);
      book(_h["n_M_tagged_jets_5j_3b"],       39, 1, 1);
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {

      vector<DressedLepton> leptons;
      for (auto &lep : apply<DressedLeptons>(event, "elecs").dressedLeptons()) { leptons.push_back(lep); }
      for (auto &lep : apply<DressedLeptons>(event, "muons").dressedLeptons()) { leptons.push_back(lep); }

      vector<DressedLepton> Electrons     = apply<DressedLeptons>(event, "elecs").dressedLeptons();
      vector<DressedLepton> Muons         = apply<DressedLeptons>(event, "muons").dressedLeptons();
      vector<DressedLepton> vetoElectrons = apply<DressedLeptons>(event, "vetoelecs").dressedLeptons();
      vector<DressedLepton> vetoMuons     = apply<DressedLeptons>(event, "vetomuons").dressedLeptons();


      // Retrieve clustered jets, sorted by pT, with a minimum pT cut
      const Jets nonisojets = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.4);
      const Jets jets = discardIfAnyDeltaRLess(nonisojets, leptons, 0.4);

      Jets bjets;
      Jets lfjets;

      for (const Jet& jet : jets) {
        if (jet.bTagged())  bjets += jet;
        else lfjets += jet;
      }

      size_t njets   = jets.size();
      size_t nbjets  = bjets.size();
      size_t nlfjets = lfjets.size();

      bool pass_ljets   = (Electrons.size() == 1 && vetoElectrons.size() == 1 && vetoMuons.size()     == 0)
                           || (Muons.size() == 1 && vetoMuons.size()     == 1 && vetoElectrons.size() == 0);
      if (!(pass_ljets))            vetoEvent;
      if (nbjets < 3 || njets < 5)  vetoEvent;


      double ht_jets   = sum(  jets, Kin::pT, 0.0);
      double ht_bjets  = sum( bjets, Kin::pT, 0.0);
      double ht_lfjets = sum(lfjets, Kin::pT, 0.0);

      size_t ind1 {0}, ind2 {0}; double mindr = 999.;
      double maxmbb  =   0.;
      double sum_dr = 0.0;
      size_t sum_n_dr = 0;
      for (size_t i = 0; i < bjets.size(); ++i) {
        for (size_t j = i + 1; j < bjets.size(); ++j) {
          double dr = deltaR(bjets[i], bjets[j]);
          double mbb = (bjets[i].momentum() + bjets[j].momentum()).mass();

          if (dr < mindr) {
            ind1 = i;
            ind2 = j;
            mindr = dr;
          }
          sum_dr += dr;
          sum_n_dr += 1;

          if (mbb > maxmbb) maxmbb = mbb;
        }
      }

      FourMomentum bb_closest = bjets[ind1].momentum() + bjets[ind2].momentum();
      double dr_closest = deltaR(bjets[ind1], bjets[ind2]);

      // 6j3b3l and 7j4b3l:
      double gen_w_mass = 79.6;
      size_t lfind1 {0}, lfind2 {0}; double minwmasdiff = 9999.;
      size_t lfind = 0;
      double phi_first_ljet_extra_bjet_soft = 0.;
      if (nlfjets >= 3) {
        for (size_t i = 0; i < nlfjets; ++i) {
          for (size_t j = i + 1; j < nlfjets; ++j) {
            double wmasdiff = abs(gen_w_mass - (lfjets[i].momentum()+lfjets[j].momentum()).mass());
            if (wmasdiff < minwmasdiff) {
              lfind1 = i;
              lfind2 = j;
              minwmasdiff = wmasdiff;
            }
          }
        }
        for (size_t i = 0; i < 3; ++i){
            if (!(( i == lfind1) || ( i == lfind2))){
            lfind = i;
            break;
            }
        }
        phi_first_ljet_extra_bjet_soft = deltaPhi(lfjets[lfind].momentum(),bjets[nbjets - 1].momentum());
      }

      // Fill histograms
      // 5j3b category
      if (njets >= 5 && nbjets >= 3) {
        _h["xsec"]                          ->fill(4);
        _h["Njets_5j_3b"]                   ->fill(njets);
        _h["n_M_tagged_jets_5j_3b"]         ->fill(nbjets);
        _h["bjet3_pt_5j_3b"]                ->fill(bjets[2].pT()/GeV);
        _h["bjet3_abseta_5j_3b"]            ->fill(bjets[2].abseta()/GeV);
        _h["jet_Ht_5j_3b"]                  ->fill(ht_jets/GeV);
        _h["bJet_Ht_5j_3b"]                 ->fill(ht_bjets/GeV);
      }

      // 6j4b category
      if (njets >= 6 && nbjets >= 4) {
        _h["xsec"]                          ->fill(2);
        _h["Njets_6j_4b"]                   ->fill(njets);
        _h["bjet3_pt_6j_4b"]                ->fill(bjets[2].pT()/GeV);
        _h["bjet3_abseta_6j_4b"]            ->fill(bjets[2].abseta()/GeV);
        _h["bjet4_pt_6j_4b"]                ->fill(bjets[3].pT()/GeV);
        _h["bjet4_abseta_6j_4b"]            ->fill(bjets[3].abseta()/GeV);
        _h["jet_Ht_6j_4b"]                  ->fill(ht_jets/GeV);
        _h["bJet_Ht_6j_4b"]                 ->fill(ht_bjets/GeV);
        _h["average_DR_6j_4b"]              ->fill(sum_dr/sum_n_dr);
        _h["largest_Mbb_6j_4b"]             ->fill(maxmbb);
        _h["extra_jet_DR_6j_4b"]            ->fill(dr_closest);
        _h["extra_jet_M_6j_4b"]             ->fill(bb_closest.mass()/GeV);
        _h["extra_jet_pt_6j_4b"]            ->fill(bb_closest.pT()/GeV);
        _h["extra_jet_abseta_6j_4b"]        ->fill(bb_closest.abseta()/GeV);
        _h["extra_jet1_pt_6j_4b"]           ->fill(bjets[ind1].pT()/GeV);
        _h["extra_jet1_abseta_6j_4b"]       ->fill(bjets[ind1].abseta()/GeV);
        _h["extra_jet2_pt_6j_4b"]           ->fill(bjets[ind2].pT()/GeV);
        _h["extra_jet2_abseta_6j_4b"]       ->fill(bjets[ind2].abseta()/GeV);
      }
      // 6j3b3l category
      if (njets >= 6 && nbjets >= 3 && nlfjets >= 3) {
        _h["xsec"]                          ->fill(3);
        _h["extra_lightJet_pt_6j_3b_3lj"]   ->fill(lfjets[lfind].pT()/GeV);
        _h["dPhi_lj_bsoft_6j_3b_3lj"]       ->fill(phi_first_ljet_extra_bjet_soft);
        _h["lightJet_Ht_6j_3b_3lj"]         ->fill(ht_lfjets/GeV);
      }

      // 7j4b3l category
      if (njets >= 7 && nbjets >= 4 && nlfjets >= 3) {
        _h["xsec"]                          ->fill(1);
        _h["extra_lightJet_pt_7j_4b_3lj"]   ->fill(lfjets[lfind].pT()/GeV);
        _h["dPhi_lj_bsoft_7j_4b_3lj"]       ->fill(phi_first_ljet_extra_bjet_soft);
        _h["lightJet_Ht_7j_4b_3lj"]         ->fill(ht_lfjets/GeV);
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      const double sf = crossSection() / femtobarn / sumOfWeights();

      // move the overflow into the last "true" bin of the distribution _just_ for the following
      std::list<string> overflows =
      {
        "Njets_5j_3b",
        "bJet_Ht_5j_3b",
        "bjet3_pt_5j_3b",
        "jet_Ht_5j_3b",
        "n_M_tagged_jets_5j_3b",
        "extra_lightJet_pt_6j_3b_3lj",
        "lightJet_Ht_6j_3b_3lj",
        "Njets_6j_4b",
        "bJet_Ht_6j_4b",
        "bjet3_pt_6j_4b",
        "bjet4_pt_6j_4b",
        "extra_jet1_pt_6j_4b",
        "extra_jet2_pt_6j_4b",
        "extra_jet_M_6j_4b",
        "extra_jet_abseta_6j_4b",
        "extra_jet_pt_6j_4b",
        "jet_Ht_6j_4b",
        "largest_Mbb_6j_4b",
        "extra_lightJet_pt_7j_4b_3lj",
        "lightJet_Ht_7j_4b_3lj"
      };
      for (auto of = overflows.begin(); of != overflows.end(); of++) {
        _h[*of]->fillBin(_h[*of]->numBins()-1, _h[*of]->overflow().sumW());
        _h[*of]->overflow().reset();
      }

      for (auto const& h : _h) {
        scale(h.second, sf);
        if (h.first == "xsec")
          continue;
        normalize(h.second);
      }

    }

    ///@}


    /// @name Histograms
    ///@{
    map<string, Histo1DPtr> _h;
    ///@}


  };


  RIVET_DECLARE_PLUGIN(CMS_2023_I2703254);

}