Rivet analyses referenceCMS_2023_I2703254Inclusive and differential cross section measurements of ttbb production in the lepton+jets channel at 13 TeVExperiment: CMS (LHC) Inspire ID: I2703254 Status: VALIDATED Authors:
Beam energies: (6500.0, 6500.0) GeV Run details:
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}
|