Rivet analyses referenceCMS_2018_I1690148Measurement of jet substructure observables in $t\bar{t}$ events from $pp$ collisions at 13 TeVExperiment: CMS (LHC) Inspire ID: 1690148 Status: VALIDATED Authors:
Beam energies: (6500.0, 6500.0) GeV Run details:
A measurement of jet substructure observables is presented using $t\bar{t}$ events in the lepton+jets channel from proton-proton collisions at $\sqrt{s}=13$ TeV recorded by the CMS experiment at the LHC, corresponding to an integrated luminosity of 35.9/fb. Multiple jet substructure observables are measured for jets identified as bottom, light-quark, and gluon jets, as well as for inclusive jets (no flavor information). The results are unfolded to the particle level and compared to next-to-leading-order predictions from POWHEG interfaced with the parton shower generators PYTHIA 8 and HERWIG 7, as well as from SHERPA 2 and DIRE 2. A value of the strong coupling at the Z boson mass, $\alpha_{S}(m_{Z}) = 0.115^{+0.015}_{-0.013}$, is extracted from the substructure data at leading-order plus leading-log accuracy. Source code: CMS_2018_I1690148.cc 1#include "Rivet/Analysis.hh"
2#include "Rivet/Projections/FinalState.hh"
3#include "Rivet/Projections/FastJets.hh"
4#include "Rivet/Projections/ChargedLeptons.hh"
5#include "Rivet/Projections/LeptonFinder.hh"
6#include "Rivet/Projections/IdentifiedFinalState.hh"
7#include "Rivet/Projections/PromptFinalState.hh"
8#include "Rivet/Projections/VetoedFinalState.hh"
9#include "Rivet/Projections/InvMassFinalState.hh"
10#include "Rivet/Projections/MissingMomentum.hh"
11#include "Rivet/Math/MatrixN.hh"
12
13#include "fastjet/contrib/Nsubjettiness.hh"
14#include "fastjet/contrib/EnergyCorrelator.hh"
15
16namespace Rivet {
17
18
19 /// Measurement of jet substructure observables in $t\bar{t}$ events from $pp$ collisions at 13~TeV
20 class CMS_2018_I1690148 : public Analysis {
21 public:
22
23 enum Reconstruction { CHARGED=0, ALL=1 };
24 enum Observable { MULT=0, PTDS=1, GA_LHA=2, GA_WIDTH=3,
25 GA_THRUST=4, ECC=5, ZG=6, ZGDR=7, NSD=8,
26 TAU21=9, TAU32=10, TAU43=11, C1_00=12,
27 C1_02=13, C1_05=14, C1_10=15, C1_20=16,
28 C2_00=17, C2_02=18, C2_05=19, C2_10=20,
29 C2_20=21, C3_00=22, C3_02=23, C3_05=24,
30 C3_10=25, C3_20=26, M2_B1=27, N2_B1=28,
31 N3_B1=29, M2_B2=30, N2_B2=31, N3_B2=32 };
32 enum Flavor { INCL=0, BOTTOM=1, QUARK=2, GLUON=3 };
33
34
35 /// Minimal constructor
36 RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2018_I1690148);
37
38
39 /// @name Analysis methods
40 /// @{
41
42
43 /// Set up projections and book histograms
44 void init() {
45
46 // Cuts
47 particle_cut = (Cuts::abseta < 5.0 && Cuts::pT > 0.*GeV);
48 lepton_cut = (Cuts::abseta < 2.4 && Cuts::pT > 15.*GeV);
49 jet_cut = (Cuts::abseta < 2.4 && Cuts::pT > 30.*GeV);
50
51 // Generic final state
52 FinalState fs(particle_cut);
53
54 // Dressed leptons
55 ChargedLeptons charged_leptons(fs);
56 IdentifiedFinalState photons(fs);
57 photons.acceptIdPair(PID::PHOTON);
58
59 PromptFinalState prompt_leptons(charged_leptons);
60 prompt_leptons.acceptMuonDecays(true);
61 prompt_leptons.acceptTauDecays(true);
62
63 PromptFinalState prompt_photons(photons);
64 prompt_photons.acceptMuonDecays(true);
65 prompt_photons.acceptTauDecays(true);
66
67 // NB. useDecayPhotons=true allows for photons with tau ancestor; photons from hadrons are vetoed by the PromptFinalState;
68 LeptonFinder dressed_leptons(prompt_leptons, prompt_photons, 0.1, lepton_cut);
69 declare(dressed_leptons, "LeptonFinder");
70
71 // Projection for jets
72 VetoedFinalState fsForJets(fs);
73 fsForJets.addVetoOnThisFinalState(dressed_leptons);
74 declare(FastJets(fsForJets, JetAlg::ANTIKT, 0.4, JetMuons::ALL, JetInvisibles::NONE), "Jets");
75
76 // Booking of histograms
77 int d = 0;
78 for (int r = 0; r < 2; ++r) { // reconstruction (charged, all)
79 for (int o = 0; o < 33; ++o) { // observable
80 d += 1;
81 for (int f = 0; f < 4; ++f) { // flavor
82 book(_h[r][o][f], d, 1, f+1);
83 }
84 }
85 }
86 }
87
88
89 void analyze(const Event& event) {
90
91 // select ttbar -> lepton+jets
92 const DressedLeptons& leptons = apply<LeptonFinder>(event, "LeptonFinder").dressedLeptons();
93 int nsel_leptons = 0;
94 for (const DressedLepton& lepton : leptons) {
95 if (lepton.pt() > 26.) nsel_leptons += 1; else vetoEvent; // found veto lepton
96 }
97 if (nsel_leptons != 1) vetoEvent;
98
99 const Jets all_jets = apply<FastJets>(event, "Jets").jetsByPt(jet_cut);
100 if (all_jets.size() < 4) vetoEvent;
101
102 // categorize jets
103 int nsel_bjets = 0;
104 int nsel_wjets = 0;
105 Jets jets[4];
106 for (const Jet& jet : all_jets) {
107 // check for jet-lepton overlap -> do not consider for selection
108 if (deltaR(jet, leptons[0]) < 0.4) continue;
109
110 bool overlap = false;
111 bool w_jet = false;
112 for (const Jet& jet2 : all_jets) {
113 if (jet.momentum() == jet2.momentum()) continue;
114 // check for jet-jet overlap -> do not consider for analysis
115 if (deltaR(jet, jet2) < 0.8)
116 overlap = true;
117 // check for W candidate
118 if (jet.bTagged() or jet2.bTagged()) continue;
119 FourMomentum w_cand = jet.momentum() + jet2.momentum();
120 if (abs(w_cand.mass() - 80.4) < 15.) w_jet = true;
121 }
122
123 // count jets for event selection
124 if (jet.bTagged()) nsel_bjets += 1;
125 if (w_jet) nsel_wjets += 1;
126
127 // jets for analysis
128 if (jet.abseta() > 2. or overlap) continue;
129
130 if (jet.bTagged()) {
131 jets[BOTTOM].push_back(jet);
132 } else if (w_jet) {
133 jets[QUARK].push_back(jet);
134 } else {
135 jets[GLUON].push_back(jet);
136 }
137 }
138
139 if (nsel_bjets != 2) vetoEvent;
140 if (nsel_wjets < 2) vetoEvent;
141
142 // substructure analysis
143 // no loop over incl jets -> more loc but faster
144 for (int f = 1; f < 4; ++f) {
145 for (const Jet& jet : jets[f]) {
146 // apply cuts on constituents
147 vector<PseudoJet> particles[2];
148 for (const Particle& p : jet.particles(Cuts::pT > 1.*GeV)) {
149 particles[ALL].push_back( PseudoJet(p.px(), p.py(), p.pz(), p.energy()) );
150 if (p.charge3() != 0)
151 particles[CHARGED].push_back( PseudoJet(p.px(), p.py(), p.pz(), p.energy()) );
152 }
153
154 if (particles[CHARGED].size() == 0) continue;
155
156 // recluster with C/A and anti-kt+WTA
157 PseudoJet ca_jet[2];
158 JetDefinition ca_def(fastjet::cambridge_algorithm, fastjet::JetDefinition::max_allowable_R);
159 ClusterSequence ca_charged(particles[CHARGED], ca_def);
160 ClusterSequence ca_all(particles[ALL], ca_def);
161 ca_jet[CHARGED] = ca_charged.exclusive_jets(1)[0];
162 ca_jet[ALL] = ca_all.exclusive_jets(1)[0];
163
164 PseudoJet akwta_jet[2];
165 JetDefinition akwta_def(fastjet::antikt_algorithm, fastjet::JetDefinition::max_allowable_R, fastjet::RecombinationScheme::WTA_pt_scheme);
166 ClusterSequence akwta_charged(particles[CHARGED], akwta_def);
167 ClusterSequence akwta_all(particles[ALL], akwta_def);
168 akwta_jet[CHARGED] = akwta_charged.exclusive_jets(1)[0];
169 akwta_jet[ALL] = akwta_all.exclusive_jets(1)[0];
170
171 // calculate observables
172 for (int r = 0; r < 2; ++r) {
173 int mult = akwta_jet[r].constituents().size();
174 // generalized angularities
175 _h[r][MULT][INCL]->fill(mult);
176 _h[r][MULT][f]->fill(mult);
177 if (mult > 1) {
178 double ptds = getPtDs(akwta_jet[r]);
179 double ga_lha = calcGA(0.5, 1., akwta_jet[r]);
180 double ga_width = calcGA(1., 1., akwta_jet[r]);
181 double ga_thrust = calcGA(2., 1., akwta_jet[r]);
182 _h[r][PTDS][INCL]->fill(ptds);
183 _h[r][PTDS][f]->fill(ptds);
184 _h[r][GA_LHA][INCL]->fill(ga_lha);
185 _h[r][GA_LHA][f]->fill(ga_lha);
186 _h[r][GA_WIDTH][INCL]->fill(ga_width);
187 _h[r][GA_WIDTH][f]->fill(ga_width);
188 _h[r][GA_THRUST][INCL]->fill(ga_thrust);
189 _h[r][GA_THRUST][f]->fill(ga_thrust);
190 }
191 // eccentricity
192 if (mult > 3) {
193 double ecc = getEcc(akwta_jet[r]);
194 _h[r][ECC][INCL]->fill(ecc);
195 _h[r][ECC][f]->fill(ecc);
196 }
197 // N-subjettiness
198 if (mult > 2) {
199 double tau21 = getTau(2, 1, ca_jet[r]);
200 _h[r][TAU21][INCL]->fill(tau21);
201 _h[r][TAU21][f]->fill(tau21);
202 }
203 if (mult > 3) {
204 double tau32 = getTau(3, 2, ca_jet[r]);
205 _h[r][TAU32][INCL]->fill(tau32);
206 _h[r][TAU32][f]->fill(tau32);
207 }
208 if (mult > 4) {
209 double tau43 = getTau(4, 3, ca_jet[r]);
210 _h[r][TAU43][INCL]->fill(tau43);
211 _h[r][TAU43][f]->fill(tau43);
212 }
213 // soft drop
214 if (mult > 1) {
215 vector<double> sd_results = getZg(ca_jet[r]);
216 if (sd_results[0] > 0.) {
217 _h[r][ZG][INCL]->fill(sd_results[0]);
218 _h[r][ZG][f]->fill(sd_results[0]);
219 _h[r][ZGDR][INCL]->fill(sd_results[1]);
220 _h[r][ZGDR][f]->fill(sd_results[1]);
221 }
222 }
223 int nsd = getNSD(0.007, -1., ca_jet[r]);
224 _h[r][NSD][INCL]->fill(nsd);
225 _h[r][NSD][f]->fill(nsd);
226 // C-series energy correlation ratios
227 if (mult > 1) {
228 double cn_00 = getC(1, 0.0, ca_jet[r]);
229 double cn_02 = getC(1, 0.2, ca_jet[r]);
230 double cn_05 = getC(1, 0.5, ca_jet[r]);
231 double cn_10 = getC(1, 1.0, ca_jet[r]);
232 double cn_20 = getC(1, 2.0, ca_jet[r]);
233 _h[r][C1_00][INCL]->fill(cn_00);
234 _h[r][C1_02][INCL]->fill(cn_02);
235 _h[r][C1_05][INCL]->fill(cn_05);
236 _h[r][C1_10][INCL]->fill(cn_10);
237 _h[r][C1_20][INCL]->fill(cn_20);
238 _h[r][C1_00][f]->fill(cn_00);
239 _h[r][C1_02][f]->fill(cn_02);
240 _h[r][C1_05][f]->fill(cn_05);
241 _h[r][C1_10][f]->fill(cn_10);
242 _h[r][C1_20][f]->fill(cn_20);
243 }
244 if (mult > 2) {
245 double cn_00 = getC(2, 0.0, ca_jet[r]);
246 double cn_02 = getC(2, 0.2, ca_jet[r]);
247 double cn_05 = getC(2, 0.5, ca_jet[r]);
248 double cn_10 = getC(2, 1.0, ca_jet[r]);
249 double cn_20 = getC(2, 2.0, ca_jet[r]);
250 _h[r][C2_00][INCL]->fill(cn_00);
251 _h[r][C2_02][INCL]->fill(cn_02);
252 _h[r][C2_05][INCL]->fill(cn_05);
253 _h[r][C2_10][INCL]->fill(cn_10);
254 _h[r][C2_20][INCL]->fill(cn_20);
255 _h[r][C2_00][f]->fill(cn_00);
256 _h[r][C2_02][f]->fill(cn_02);
257 _h[r][C2_05][f]->fill(cn_05);
258 _h[r][C2_10][f]->fill(cn_10);
259 _h[r][C2_20][f]->fill(cn_20);
260 }
261 if (mult > 3) {
262 double cn_00 = getC(3, 0.0, ca_jet[r]);
263 double cn_02 = getC(3, 0.2, ca_jet[r]);
264 double cn_05 = getC(3, 0.5, ca_jet[r]);
265 double cn_10 = getC(3, 1.0, ca_jet[r]);
266 double cn_20 = getC(3, 2.0, ca_jet[r]);
267 _h[r][C3_00][INCL]->fill(cn_00);
268 _h[r][C3_02][INCL]->fill(cn_02);
269 _h[r][C3_05][INCL]->fill(cn_05);
270 _h[r][C3_10][INCL]->fill(cn_10);
271 _h[r][C3_20][INCL]->fill(cn_20);
272 _h[r][C3_00][f]->fill(cn_00);
273 _h[r][C3_02][f]->fill(cn_02);
274 _h[r][C3_05][f]->fill(cn_05);
275 _h[r][C3_10][f]->fill(cn_10);
276 _h[r][C3_20][f]->fill(cn_20);
277 }
278 // M/N-series energy correlation ratios
279 if (mult > 2) {
280 double m2_b1 = getM(2, 1., ca_jet[r]);
281 double m2_b2 = getM(2, 2., ca_jet[r]);
282 double n2_b1 = getN(2, 1., ca_jet[r]);
283 double n2_b2 = getN(2, 2., ca_jet[r]);
284 _h[r][M2_B1][INCL]->fill(m2_b1);
285 _h[r][M2_B2][INCL]->fill(m2_b2);
286 _h[r][N2_B1][INCL]->fill(n2_b1);
287 _h[r][N2_B2][INCL]->fill(n2_b2);
288 _h[r][M2_B1][f]->fill(m2_b1);
289 _h[r][M2_B2][f]->fill(m2_b2);
290 _h[r][N2_B1][f]->fill(n2_b1);
291 _h[r][N2_B2][f]->fill(n2_b2);
292 }
293 if (mult > 3) {
294 double n3_b1 = getN(3, 1., ca_jet[r]);
295 double n3_b2 = getN(3, 2., ca_jet[r]);
296 _h[r][N3_B1][INCL]->fill(n3_b1);
297 _h[r][N3_B2][INCL]->fill(n3_b2);
298 _h[r][N3_B1][f]->fill(n3_b1);
299 _h[r][N3_B2][f]->fill(n3_b2);
300 }
301 }
302 }
303 }
304 }
305
306
307 void finalize() {
308 for (int r = 0; r < 2; ++r) { // reconstruction (charged, all)
309 for (int o = 0; o < 33; ++o) { // observable
310 for (int f = 0; f < 4; ++f) { // flavor
311 normalize(_h[r][o][f], 1.0, false);
312 }
313 }
314 }
315 }
316
317 /// @}
318
319
320 private:
321
322 double deltaR(PseudoJet j1, PseudoJet j2) {
323 double deta = j1.eta() - j2.eta();
324 double dphi = j1.delta_phi_to(j2);
325 return sqrt(deta*deta + dphi*dphi);
326 }
327
328
329 double getPtDs(PseudoJet jet) {
330 double mult = jet.constituents().size();
331 double sumpt = 0.; // would be jet.pt() in WTA scheme but better keep it generic
332 double sumpt2 = 0.;
333 for (auto p : jet.constituents()) {
334 sumpt += p.pt();
335 sumpt2 += pow(p.pt(), 2);
336 }
337 double ptd = sumpt2/pow(sumpt,2);
338 return max(0., sqrt((ptd-1./mult) * mult/(mult-1.)));
339 }
340
341
342 double calcGA(double beta, double kappa, PseudoJet jet) {
343 double sumpt = 0.;
344 for (const auto& p : jet.constituents()) {
345 sumpt += p.pt();
346 }
347 double ga = 0.;
348 for (auto p : jet.constituents()) {
349 ga += pow(p.pt()/sumpt, kappa) * pow(deltaR(jet, p)/0.4, beta);
350 }
351 return ga;
352 }
353
354
355 double getEcc(PseudoJet jet) {
356 // Covariance matrix
357 Matrix<2> M;
358 for (const auto& p : jet.constituents()) {
359 Matrix<2> MPart;
360 MPart.set(0, 0, (p.eta() - jet.eta()) * (p.eta() - jet.eta()));
361 MPart.set(0, 1, (p.eta() - jet.eta()) * mapAngleMPiToPi(p.phi() - jet.phi()));
362 MPart.set(1, 0, mapAngleMPiToPi(p.phi() - jet.phi()) * (p.eta() - jet.eta()));
363 MPart.set(1, 1, mapAngleMPiToPi(p.phi() - jet.phi()) * mapAngleMPiToPi(p.phi() - jet.phi()));
364 M += MPart * p.e();
365 }
366 // Calculate eccentricity from eigenvalues
367 // Check that the matrix is symmetric.
368 const bool isSymm = M.isSymm();
369 if (!isSymm) {
370 MSG_ERROR("Error: energy tensor not symmetric:");
371 MSG_ERROR("[0,1] vs. [1,0]: " << M.get(0,1) << ", " << M.get(1,0));
372 }
373 // If not symmetric, something's wrong (we made sure the error msg appeared first).
374
375 assert(isSymm);
376 const double a = M.get(0,0);
377 const double b = M.get(1,1);
378 const double c = M.get(1,0);
379
380 const double l1 = 0.5*(a+b+sqrt( (a-b)*(a-b) + 4 *c*c));
381 const double l2 = 0.5*(a+b-sqrt( (a-b)*(a-b) + 4 *c*c));
382
383 return 1. - l2/l1;
384 }
385
386
387 double getTau(int N, int M, PseudoJet jet) {
388 fjcontrib::NsubjettinessRatio tau_ratio(N, M, fjcontrib::OnePass_WTA_KT_Axes(),
389 fjcontrib::NormalizedMeasure(1.0, 0.4));
390 return tau_ratio(jet);
391 }
392
393
394 vector<double> getZg(PseudoJet jet) {
395 PseudoJet jet0 = jet;
396 PseudoJet jet1, jet2;
397 double zg = 0.;
398 while (zg < 0.1 && jet0.has_parents(jet1, jet2)) {
399 zg = jet2.pt()/jet0.pt();
400 jet0 = jet1;
401 }
402 if (zg < 0.1) return {-1., -1.};
403 return {zg, jet1.delta_R(jet2)};
404 }
405
406
407 int getNSD(double zcut, double beta, PseudoJet jet) {
408 PseudoJet jet0 = jet;
409 PseudoJet jet1, jet2;
410 int nsd = 0.;
411 double zg = 0.;
412 while (jet0.has_parents(jet1, jet2)) {
413 zg = jet2.pt()/jet0.pt();
414 if (zg > zcut * pow(jet1.delta_R(jet2)/0.4, beta))
415 nsd += 1;
416 jet0 = jet1;
417 }
418 return nsd;
419 }
420
421
422 double getC(int N, double beta, PseudoJet jet) {
423 fjcontrib::EnergyCorrelatorDoubleRatio C(N, beta);
424 return C(jet);
425 }
426
427
428 double getM(int N, double beta, PseudoJet jet) {
429 fjcontrib::EnergyCorrelatorMseries CM(N, beta);
430 return CM(jet);
431 }
432
433
434 double getN(int N, double beta, PseudoJet jet) {
435 fjcontrib::EnergyCorrelatorNseries CN(N, beta);
436 return CN(jet);
437 }
438
439
440 private:
441
442 // @name Histogram data members
443 /// @{
444 Cut particle_cut, lepton_cut, jet_cut;
445 Histo1DPtr _h[2][33][4];
446 /// @}
447
448 };
449
450
451 RIVET_DECLARE_PLUGIN(CMS_2018_I1690148);
452
453}
|