Rivet analyses referenceCMS_2018_I1663958ttbar lepton+jets 13 TeVExperiment: CMS (LHC) Inspire ID: 1663958 Status: VALIDATED Authors:
Beam energies: (6500.0, 6500.0) GeV Run details:
Abstract: Differential and double-differential cross sections for the production of top quark pairs in proton-proton collisions at $\sqrt{s} = 13$\,TeV are measured as a function of kinematic variables of the top quarks and the top quark-antiquark ($\text{t}\bar{\text{t}}$) system. In addition, kinematic variables and multiplicities of jets associated with the $\text{t}\bar{\text{t}}$ production are measured. This analysis is based on data collected by the CMS experiment at the LHC in 2016 corresponding to an integrated luminosity of 35.8\,fb$^{-1}$. The measurements are performed in the lepton+jets decay channels with a single muon or electron and jets in the final state. The differential cross sections are presented at the particle level, within a phase space close to the experimental acceptance, and at the parton level in the full phase space. The results are compared to several standard model predictions that use different methods and approximations. The kinematic variables of the top quarks and the $\text{t}\bar{\text{t}}$ system are reasonably described in general, though none predict all the measured distributions. In particular, the transverse momentum distribution of the top quarks is more steeply falling than predicted. The kinematic distributions and multiplicities of jets are adequately modeled by certain combinations of next-to-leading-order calculations and parton shower models. Rivet: This analysis is to be run on $\text{t}\bar{\text{t}}$ Monte Carlo. The particle-level phase space is defined using the following definitions: \begin{description} \item[lepton]: an electron or muon with $p_\text{T}>30\,\text{GeV}$ and $|\eta|<2.4$, dressed within a cone of radius 0.1, \item[jet]: a jet is reconstructed with the anti-$k_\text{T}$ algorithm with a radius of 0.4, after removing the neutrinos and dressed leptons, with $p_\text{T]>25\,\text{GeV}$ and $|\eta|<2.4$, \item[b-jet]: a jet that contains a B-hadron. \end{description} A W boson is reconstructed from a lepton and the sum of the neutrino energies, while another W boson is reconstructed from a light jet pair. The two top quarks are reconstructed by combining b jets to these W boson. A check based on the W boson and top quark masses is performed to choose the proper combinations. Source code: CMS_2018_I1663958.cc 1#include "Rivet/Analysis.hh"
2#include "Rivet/Projections/FinalState.hh"
3#include "Rivet/Projections/VisibleFinalState.hh"
4#include "Rivet/Projections/VetoedFinalState.hh"
5#include "Rivet/Projections/FastJets.hh"
6#include "Rivet/Projections/LeptonFinder.hh"
7
8namespace Rivet {
9
10
11 /// @brief ttbar lepton+jets 13 TeV
12 class CMS_2018_I1663958 : public Analysis {
13 public:
14
15 RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2018_I1663958);
16
17 void init() {
18
19 const FinalState fs(Cuts::abseta < 6.);
20 const VisibleFinalState vfs(Cuts::abseta < 6.);
21
22 VetoedFinalState invisibles(fs);
23 invisibles.addVetoOnThisFinalState(vfs);
24 declare(invisibles, "Invisibles");
25
26 FinalState all_photons(vfs, Cuts::abspid == PID::PHOTON);
27 FinalState leptons(vfs, Cuts::abspid == PID::ELECTRON || Cuts::abspid == PID::MUON);
28 LeptonFinder dressed_leptons(leptons, all_photons, 0.1, Cuts::abseta < 2.4 && Cuts::pT > 15*GeV);
29 declare(dressed_leptons, "MyLeptons");
30
31 VetoedFinalState photons(all_photons);
32 photons.addVetoOnThisFinalState(dressed_leptons);
33 declare(photons, "MyPhotons");
34
35 VetoedFinalState isolationparticles(vfs);
36 isolationparticles.addVetoOnThisFinalState(dressed_leptons);
37 declare(isolationparticles, "IsoParticles");
38
39 declare(FastJets(vfs, JetAlg::ANTIKT, 0.4), "Jets");
40
41 book(_h["thadpt"], 1, 1, 1);
42 book(_h["thady"], 3, 1, 1);
43 book(_h["tleppt"], 5, 1, 1);
44 book(_h["tlepy"], 7, 1, 1);
45 book(_h["ttm"], 9, 1, 1);
46 book(_h["ttpt"], 11, 1, 1);
47 book(_h["tty"], 13, 1, 1);
48 book(_h["njet"], 15, 1, 1);
49
50 const vector<double> njetbins{-0.5, 0.5, 1.5, 2.5, 3.5};
51 book(_b["njet_ttm"], njetbins);
52 book(_b["njet_ttm_norm"], njetbins);
53 book(_b["njet_thadpt"], njetbins);
54 book(_b["njet_thadpt_norm"], njetbins);
55 book(_b["njet_ttpt"], njetbins);
56 book(_b["njet_ttpt_norm"], njetbins);
57 for (size_t i = 1; i < _b["njet_ttm"]->numBins()+1; ++i) {
58 book(_b["njet_ttm"]->bin(i), 16 + i, 1, 1);
59 book(_b["njet_ttm_norm"]->bin(i), 98 + i, 1, 1);
60 book(_b["njet_thadpt"]->bin(i), 21 + i, 1, 1);
61 book(_b["njet_thadpt_norm"]->bin(i), 103 + i, 1, 1);
62 book(_b["njet_ttpt"]->bin(i), 26 + i, 1, 1);
63 book(_b["njet_ttpt_norm"]->bin(i), 108 + i, 1, 1);
64 }
65
66 book(_b["thady_thadpt"], {0.0, 0.5, 1.0, 1.5, 2.5}, {"d32-x01-y01", "d33-x01-y01", "d34-x01-y01", "d35-x01-y01"});
67 book(_b["thady_thadpt_norm"], {0.0, 0.5, 1.0, 1.5, 2.5}, {"d114-x01-y01", "d115-x01-y01", "d116-x01-y01", "d117-x01-y01"});
68
69 book(_b["ttm_tty"], {300., 450., 625., 850., 2000.}, {"d37-x01-y01", "d38-x01-y01", "d39-x01-y01", "d40-x01-y01"});
70 book(_b["ttm_tty_norm"], {300., 450., 625., 850., 2000.}, {"d119-x01-y01", "d120-x01-y01", "d121-x01-y01", "d122-x01-y01"});
71
72 book(_b["thadpt_ttm"], {0., 90., 180., 270., 800.}, {"d42-x01-y01", "d43-x01-y01", "d44-x01-y01", "d45-x01-y01"});
73 book(_b["thadpt_ttm_norm"], {0., 90., 180., 270., 800.}, {"d124-x01-y01", "d125-x01-y01", "d126-x01-y01", "d127-x01-y01"});
74
75 const vector<double> jetbins{-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5};
76 book(_b["jetspt"], jetbins);
77 book(_b["jetspt_norm"], jetbins);
78 book(_b["jetseta"], jetbins);
79 book(_b["jetseta_norm"], jetbins);
80 book(_b["jetsdr"], jetbins);
81 book(_b["jetsdr_norm"], jetbins);
82 book(_b["jetsdrtops"], jetbins);
83 book(_b["jetsdrtops_norm"], jetbins);
84
85 for (size_t i = 1; i < _b["jetspt"]->numBins()+1; ++i) {
86 book(_b["jetspt"]->bin(i), 46 + i, 1, 1);
87 book(_b["jetspt_norm"]->bin(i), 128 + i, 1, 1);
88 book(_b["jetseta"]->bin(i), 55 + i, 1, 1);
89 book(_b["jetseta_norm"]->bin(i), 137 + i, 1, 1);
90 book(_b["jetsdr"]->bin(i), 64 + i, 1, 1);
91 book(_b["jetsdr_norm"]->bin(i), 146 + i, 1, 1);
92 book(_b["jetsdrtops"]->bin(i), 73 + i, 1, 1);
93 book(_b["jetsdrtops_norm"]->bin(i), 155 + i, 1, 1);
94 }
95
96 book(_b["njetspt"], {0., 40., 60., 80., 120.}, {"d169-x01-y01", "d170-x01-y01", "d171-x01-y01", "d172-x01-y01"});
97
98 book(_h["thadpt_norm"], 83, 1, 1);
99 book(_h["thady_norm"], 85, 1, 1);
100 book(_h["tleppt_norm"], 87, 1, 1);
101 book(_h["tlepy_norm"], 89, 1, 1);
102 book(_h["ttm_norm"], 91, 1, 1);
103 book(_h["ttpt_norm"], 93, 1, 1);
104 book(_h["tty_norm"], 95, 1, 1);
105 book(_h["njet_norm"], 97, 1, 1);
106 /// @todo Memory leak?
107 book(m_hist_gap1, "d165-x01-y01");
108 book(m_hist_gap2, "d167-x01-y01");
109 }
110
111
112 void analyze(const Event& event) {
113
114 Jets bjets, ljets;
115 Particles leptons, vetoleptons, additionalobjects, additionaljets;
116
117 const Particles& isopars = apply<VetoedFinalState>(event, "IsoParticles").particles();
118 const Particles& dressedleptons = apply<LeptonFinder>(event, "MyLeptons").particles();
119 for (const Particle& lep : dressedleptons) {
120 double isolation = sum(select(isopars, deltaRLess(lep, 0.4)), Kin::pT, 0.);
121 isolation = isolation/lep.pt();
122 if (isolation > 0.35) continue;
123 if (lep.pt() > 30*GeV && lep.abseta() < 2.4) leptons += lep;
124 else if(lep.pt() > 15*GeV && lep.abseta() < 2.4) vetoleptons += lep;
125 }
126
127 const Particles& photons = apply<VetoedFinalState>(event, "MyPhotons").particles(Cuts::abseta < 2.4 && Cuts::pT > 15*GeV);
128 for (const Particle& ph : photons) {
129 double isolation = sum(select(isopars, deltaRLess(ph, 0.4)), Kin::pT, 0.);
130 isolation = isolation/ph.pt() - 1.;
131 if (isolation > 0.25) continue;
132 additionalobjects += ph;
133 }
134
135 FourMomentum invsum(0.,0.,0.,0.);
136 const Particles& invfspars = apply<FinalState>(event, "Invisibles").particles();
137 for (const Particle& par : invfspars) invsum += par.mom();
138
139 Jets allJets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::abseta < 2.4 && Cuts::pT > 25*GeV);
140 idiscardIfAnyDeltaRLess(allJets, leptons, 0.4);
141 idiscardIfAnyDeltaRLess(allJets, vetoleptons, 0.4);
142 idiscardIfAnyDeltaRLess(allJets, additionalobjects, 0.4);
143 for (const Jet& jet : allJets) {
144 if (jet.bTagged()) bjets += jet;
145 else ljets += jet;
146 }
147
148 //Semi-leptonic reconstruction
149 if (leptons.size() != 1 || vetoleptons.size() || bjets.size() < 2 || ljets.size() < 2) vetoEvent;
150
151 FourMomentum wl = invsum + leptons[0].mom();
152
153 Particle thad, tlep;
154 Particles tlep_decay(3), thad_decay(3);
155 double Kmin = numeric_limits<double>::max();
156 for (size_t a = 0 ; a < ljets.size() ; ++a){
157 const Jet& lja = ljets[a];
158 for (size_t b = 0 ; b < a ; ++b) {
159 const Jet& ljb = ljets[b];
160 FourMomentum wh(lja.momentum() + ljb.momentum());
161 for (const Jet& bjh : bjets) {
162 FourMomentum th(wh + bjh.momentum());
163 for (const Jet& bjl : bjets) {
164 if (&bjh == &bjl) continue;
165 FourMomentum tl(wl + bjl.momentum());
166
167 double K = pow(wh.mass() - 80.4, 2) + pow(th.mass() - 172.5, 2) + pow(tl.mass() - 172.5, 2);
168 if (K < Kmin)
169 {
170 Kmin = K;
171 thad = Particle(6, th);
172 thad_decay[0] = Particle(5, bjh);
173 thad_decay[1] = lja.pt() > ljb.pt() ? Particle(1, lja) : Particle(1, ljb);
174 thad_decay[2] = lja.pt() <= ljb.pt() ? Particle(1, lja) : Particle(1, ljb);
175 tlep = Particle(-6, tl);
176 tlep_decay[0] = Particle(5, bjl);
177 tlep_decay[1] = leptons[0];
178 tlep_decay[2] = Particle(-1*(leptons[0].pid()+1), invsum);
179 }
180 }
181 }
182 }
183 }
184
185 Particles tt_jets({ tlep_decay[0], thad_decay[0], thad_decay[1], thad_decay[2] });
186
187 const double eps = 1E-5;
188 for (const Jet& jet : bjets) {
189 if(jet.pt() < 30*GeV || jet.abseta() > 2.4) continue;
190 if(find_if(tt_jets.begin(), tt_jets.end(), [&](const Particle& par){return deltaR(jet, par) < eps;}) != tt_jets.end()) continue;
191 additionaljets += Particle(5, jet.mom());
192 }
193 for (const Jet& jet : ljets) {
194 if(jet.pt() < 30*GeV || jet.abseta() > 2.4) continue;
195 if(find_if(tt_jets.begin(), tt_jets.end(), [&](const Particle& par){return deltaR(jet, par) < eps;}) != tt_jets.end()) continue;
196 if(jet.cTagged()) additionaljets += Particle(4, jet.mom());
197 else additionaljets += Particle(1, jet.mom());
198 }
199
200 sort(additionaljets.begin(), additionaljets.end(), cmpMomByPt);
201
202 FourMomentum tt(thad.mom() + tlep.mom());
203
204 dualfill("thadpt", thad.pt()/GeV);
205 dualfill("thady", thad.absrap());
206 dualfill("tleppt", tlep.pt()/GeV);
207 dualfill("tlepy", tlep.absrap());
208 dualfill("ttm", tt.mass()/GeV);
209 dualfill("ttpt", tt.pt()/GeV);
210 dualfill("tty", tt.absrap());
211 dualfill("njet", min(additionaljets.size(), (size_t)5));
212 double njet = double(min((size_t)3, additionaljets.size()));
213 dualfill("njet_ttm", njet, tt.mass()/GeV);
214 dualfill("njet_thadpt", njet, thad.pt()/GeV);
215 dualfill("njet_ttpt", njet, tt.pt()/GeV);
216 dualfill("thady_thadpt", thad.absrap(), thad.pt()/GeV);
217 dualfill("ttm_tty", tt.mass()/GeV, tt.absrap());
218 dualfill("thadpt_ttm", thad.pt()/GeV, tt.mass()/GeV);
219 int jpos = -1;
220 for (const Particles& jets : {tt_jets, additionaljets}) {
221 for (const Particle& jet : jets) {
222 jpos++;
223 dualfill("jetspt", jpos, jet.pt()/GeV);
224 dualfill("jetseta", jpos, jet.abseta());
225 double drmin = 1E10;
226 for (const Particle& tjet : tt_jets) {
227 double dr = deltaR(jet, tjet);
228 if(dr > eps && dr < drmin) drmin = dr;
229 }
230 dualfill("jetsdr", jpos, drmin);
231 dualfill("jetsdrtops", jpos, min(deltaR(jet, thad), deltaR(jet, tlep)));
232 }
233 }
234 for (double ptcut : {30, 50, 75, 100}) {
235 _b["njetspt"]->fill(ptcut , count_if(additionaljets.begin(), additionaljets.end(),
236 [&ptcut](const Particle& j) {return j.pt() > ptcut;}));
237 }
238 }
239
240 void dualfill(const string& tag, const double value) {
241 _h[tag]->fill(value); _h[tag + "_norm"]->fill(value);
242 }
243
244 void dualfill(const string& tag, const double val1, const double val2) {
245 _b[tag]->fill(val1, val2); _b[tag + "_norm"]->fill(val1, val2);
246 }
247
248 void gapfractionfromjetpt(const string& tag, Estimate1DPtr hgap, size_t njet) {
249 size_t hn = njet+4;
250 const double total = _b[tag]->bin(1)->integral();
251 const double totalj = _b[tag]->bin(hn)->integral();
252 double acc = total-totalj;
253 double gf = acc/total;
254 if (!std::isnan(gf)) {
255 for (size_t nb = 1 ; nb < _b[tag]->bin(hn)->numBins()+1; ++nb) {
256 hgap->bin(nb).set(gf, 0.);
257 acc += _b[tag]->bin(hn)->bin(nb).sumW();
258 gf = acc/total;
259 }
260 } else MSG_WARNING("Gap fraction is NaN. Histogram not produced.");
261 }
262
263
264 void finalize() {
265
266 const double sf = crossSection()/picobarn/sumOfWeights();
267
268 gapfractionfromjetpt("jetspt", m_hist_gap1, 1);
269 gapfractionfromjetpt("jetspt", m_hist_gap2, 2);
270
271 for (auto& item : _h) {
272 if (item.first.find("_norm") != string::npos) {
273 normalize(item.second, 1.0, false);
274 }
275 else scale(item.second, sf);
276 }
277
278 for (auto& item : _b) {
279 if (item.first.find("_norm") != string::npos) {
280 const double area = item.second->sumW(false);
281 //normalize(item.second, 1.0, false);
282 if (area) scale(item.second, 1/area);
283 divByGroupWidth(item.second);
284 }
285 else if (item.first.find("njetspt") != string::npos) {
286 // skip division by bin width for secondary axis
287 scale(item.second, sf);
288 }
289 else {
290 scale(item.second, sf);
291 divByGroupWidth(item.second);
292 }
293 }
294
295 }
296
297 map<string,Histo1DPtr> _h;
298 map<string,Histo1DGroupPtr> _b;
299
300 Estimate1DPtr m_hist_gap1;
301 Estimate1DPtr m_hist_gap2;
302
303 };
304
305
306 RIVET_DECLARE_PLUGIN(CMS_2018_I1663958);
307}
|