rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

CMS_2018_I1690148

Measurement of jet substructure observables in $t\bar{t}$ events from $pp$ collisions at 13 TeV
Experiment: CMS (LHC)
Inspire ID: 1690148
Status: VALIDATED
Authors:
  • Markus Seidel
References:
  • arXiv: 1808.07340
  • CMS-TOP-17-013
  • Phys. Rev. D 98, 092014 (2018)
Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • $t\bar{t}$ events at $\sqrt{s} = 13 \text{TeV}$, lepton+jets selection at particle level.

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/DressedLeptons.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      DressedLeptons dressed_leptons(prompt_photons, prompt_leptons, 0.1, lepton_cut, true);
 69      declare(dressed_leptons, "DressedLeptons");
 70
 71      // Projection for jets
 72      VetoedFinalState fsForJets(fs);
 73      fsForJets.addVetoOnThisFinalState(dressed_leptons);
 74      declare(FastJets(fsForJets, FastJets::ANTIKT, 0.4,
 75                       JetAlg::Muons::ALL, JetAlg::Invisibles::NONE), "Jets");
 76
 77      // Booking of histograms
 78      int d = 0;
 79      for (int r = 0; r < 2; ++r) { // reconstruction (charged, all)
 80        for (int o = 0; o < 33; ++o) { // observable
 81          d += 1;
 82          for (int f = 0; f < 4; ++f) { // flavor
 83            char buffer [20];
 84            sprintf(buffer, "d%02d-x01-y%02d", d, f+1);
 85            book(_h[r][o][f], buffer);
 86          }
 87        }
 88      }
 89    }
 90
 91
 92    void analyze(const Event& event) {
 93
 94      // select ttbar -> lepton+jets
 95      const vector<DressedLepton>& leptons = applyProjection<DressedLeptons>(event, "DressedLeptons").dressedLeptons();
 96      int nsel_leptons = 0;
 97      for (const DressedLepton& lepton : leptons) {
 98        if (lepton.pt() > 26.) nsel_leptons += 1; else vetoEvent; // found veto lepton
 99      }
100      if (nsel_leptons != 1) vetoEvent;
101
102      const Jets all_jets = applyProjection<FastJets>(event, "Jets").jetsByPt(jet_cut);
103      if (all_jets.size() < 4) vetoEvent;
104
105      // categorize jets
106      int nsel_bjets = 0;
107      int nsel_wjets = 0;
108      Jets jets[4];
109      for (const Jet& jet : all_jets) {
110        // check for jet-lepton overlap -> do not consider for selection
111        if (deltaR(jet, leptons[0]) < 0.4) continue;
112
113        bool overlap = false;
114        bool w_jet   = false;
115        for (const Jet& jet2 : all_jets) {
116          if (jet.momentum() == jet2.momentum()) continue;
117          // check for jet-jet overlap -> do not consider for analysis
118          if (deltaR(jet, jet2) < 0.8)
119            overlap = true;
120          // check for W candidate
121          if (jet.bTagged() or jet2.bTagged()) continue;
122          FourMomentum w_cand = jet.momentum() + jet2.momentum();
123          if (abs(w_cand.mass() - 80.4) < 15.) w_jet = true;
124        }
125
126        // count jets for event selection
127        if (jet.bTagged()) nsel_bjets += 1;
128        if (w_jet) nsel_wjets += 1;
129
130        // jets for analysis
131        if (jet.abseta() > 2. or overlap) continue;
132
133        if (jet.bTagged()) {
134          jets[BOTTOM].push_back(jet);
135        } else if (w_jet) {
136          jets[QUARK].push_back(jet);
137        } else {
138          jets[GLUON].push_back(jet);
139        }
140      }
141
142      if (nsel_bjets != 2) vetoEvent;
143      if (nsel_wjets < 2) vetoEvent;
144
145      // substructure analysis
146      // no loop over incl jets -> more loc but faster
147      for (int f = 1; f < 4; ++f) {
148        for (const Jet& jet : jets[f]) {
149          // apply cuts on constituents
150          vector<PseudoJet> particles[2];
151          for (const Particle& p : jet.particles(Cuts::pT > 1.*GeV)) {
152            particles[ALL].push_back( PseudoJet(p.px(), p.py(), p.pz(), p.energy()) );
153            if (p.charge3() != 0)
154              particles[CHARGED].push_back( PseudoJet(p.px(), p.py(), p.pz(), p.energy()) );
155          }
156
157          if (particles[CHARGED].size() == 0) continue;
158
159          // recluster with C/A and anti-kt+WTA
160          PseudoJet ca_jet[2];
161          JetDefinition ca_def(fastjet::cambridge_algorithm, fastjet::JetDefinition::max_allowable_R);
162          ClusterSequence ca_charged(particles[CHARGED], ca_def);
163          ClusterSequence ca_all(particles[ALL], ca_def);
164          ca_jet[CHARGED] = ca_charged.exclusive_jets(1)[0];
165          ca_jet[ALL] = ca_all.exclusive_jets(1)[0];
166
167          PseudoJet akwta_jet[2];
168          JetDefinition akwta_def(fastjet::antikt_algorithm, fastjet::JetDefinition::max_allowable_R, fastjet::RecombinationScheme::WTA_pt_scheme);
169          ClusterSequence akwta_charged(particles[CHARGED], akwta_def);
170          ClusterSequence akwta_all(particles[ALL], akwta_def);
171          akwta_jet[CHARGED] = akwta_charged.exclusive_jets(1)[0];
172          akwta_jet[ALL]     = akwta_all.exclusive_jets(1)[0];
173
174          // calculate observables
175          for (int r = 0; r < 2; ++r) {
176            int mult = akwta_jet[r].constituents().size();
177            // generalized angularities
178            _h[r][MULT][INCL]->fill(mult);
179            _h[r][MULT][f]->fill(mult);
180            if (mult > 1) {
181              double ptds = getPtDs(akwta_jet[r]);
182              double ga_lha = calcGA(0.5, 1., akwta_jet[r]);
183              double ga_width = calcGA(1., 1., akwta_jet[r]);
184              double ga_thrust = calcGA(2., 1., akwta_jet[r]);
185              _h[r][PTDS][INCL]->fill(ptds);
186              _h[r][PTDS][f]->fill(ptds);
187              _h[r][GA_LHA][INCL]->fill(ga_lha);
188              _h[r][GA_LHA][f]->fill(ga_lha);
189              _h[r][GA_WIDTH][INCL]->fill(ga_width);
190              _h[r][GA_WIDTH][f]->fill(ga_width);
191              _h[r][GA_THRUST][INCL]->fill(ga_thrust);
192              _h[r][GA_THRUST][f]->fill(ga_thrust);
193            }
194            // eccentricity
195            if (mult > 3) {
196              double ecc = getEcc(akwta_jet[r]);
197              _h[r][ECC][INCL]->fill(ecc);
198              _h[r][ECC][f]->fill(ecc);
199            }
200            // N-subjettiness
201            if (mult > 2) {
202              double tau21 = getTau(2, 1, ca_jet[r]);
203              _h[r][TAU21][INCL]->fill(tau21);
204              _h[r][TAU21][f]->fill(tau21);
205            }
206            if (mult > 3) {
207              double tau32 = getTau(3, 2, ca_jet[r]);
208              _h[r][TAU32][INCL]->fill(tau32);
209              _h[r][TAU32][f]->fill(tau32);
210            }
211            if (mult > 4) {
212              double tau43 = getTau(4, 3, ca_jet[r]);
213              _h[r][TAU43][INCL]->fill(tau43);
214              _h[r][TAU43][f]->fill(tau43);
215            }
216            // soft drop
217            if (mult > 1) {
218              vector<double> sd_results = getZg(ca_jet[r]);
219              if (sd_results[0] > 0.) {
220                _h[r][ZG][INCL]->fill(sd_results[0]);
221                _h[r][ZG][f]->fill(sd_results[0]);
222                _h[r][ZGDR][INCL]->fill(sd_results[1]);
223                _h[r][ZGDR][f]->fill(sd_results[1]);
224              }
225            }
226            int nsd = getNSD(0.007, -1., ca_jet[r]);
227            _h[r][NSD][INCL]->fill(nsd);
228            _h[r][NSD][f]->fill(nsd);
229            // C-series energy correlation ratios
230            if (mult > 1) {
231              double cn_00 = getC(1, 0.0, ca_jet[r]);
232              double cn_02 = getC(1, 0.2, ca_jet[r]);
233              double cn_05 = getC(1, 0.5, ca_jet[r]);
234              double cn_10 = getC(1, 1.0, ca_jet[r]);
235              double cn_20 = getC(1, 2.0, ca_jet[r]);
236              _h[r][C1_00][INCL]->fill(cn_00);
237              _h[r][C1_02][INCL]->fill(cn_02);
238              _h[r][C1_05][INCL]->fill(cn_05);
239              _h[r][C1_10][INCL]->fill(cn_10);
240              _h[r][C1_20][INCL]->fill(cn_20);
241              _h[r][C1_00][f]->fill(cn_00);
242              _h[r][C1_02][f]->fill(cn_02);
243              _h[r][C1_05][f]->fill(cn_05);
244              _h[r][C1_10][f]->fill(cn_10);
245              _h[r][C1_20][f]->fill(cn_20);
246            }
247            if (mult > 2) {
248              double cn_00 = getC(2, 0.0, ca_jet[r]);
249              double cn_02 = getC(2, 0.2, ca_jet[r]);
250              double cn_05 = getC(2, 0.5, ca_jet[r]);
251              double cn_10 = getC(2, 1.0, ca_jet[r]);
252              double cn_20 = getC(2, 2.0, ca_jet[r]);
253              _h[r][C2_00][INCL]->fill(cn_00);
254              _h[r][C2_02][INCL]->fill(cn_02);
255              _h[r][C2_05][INCL]->fill(cn_05);
256              _h[r][C2_10][INCL]->fill(cn_10);
257              _h[r][C2_20][INCL]->fill(cn_20);
258              _h[r][C2_00][f]->fill(cn_00);
259              _h[r][C2_02][f]->fill(cn_02);
260              _h[r][C2_05][f]->fill(cn_05);
261              _h[r][C2_10][f]->fill(cn_10);
262              _h[r][C2_20][f]->fill(cn_20);
263            }
264            if (mult > 3) {
265              double cn_00 = getC(3, 0.0, ca_jet[r]);
266              double cn_02 = getC(3, 0.2, ca_jet[r]);
267              double cn_05 = getC(3, 0.5, ca_jet[r]);
268              double cn_10 = getC(3, 1.0, ca_jet[r]);
269              double cn_20 = getC(3, 2.0, ca_jet[r]);
270              _h[r][C3_00][INCL]->fill(cn_00);
271              _h[r][C3_02][INCL]->fill(cn_02);
272              _h[r][C3_05][INCL]->fill(cn_05);
273              _h[r][C3_10][INCL]->fill(cn_10);
274              _h[r][C3_20][INCL]->fill(cn_20);
275              _h[r][C3_00][f]->fill(cn_00);
276              _h[r][C3_02][f]->fill(cn_02);
277              _h[r][C3_05][f]->fill(cn_05);
278              _h[r][C3_10][f]->fill(cn_10);
279              _h[r][C3_20][f]->fill(cn_20);
280            }
281            // M/N-series energy correlation ratios
282            if (mult > 2) {
283              double m2_b1 = getM(2, 1., ca_jet[r]);
284              double m2_b2 = getM(2, 2., ca_jet[r]);
285              double n2_b1 = getN(2, 1., ca_jet[r]);
286              double n2_b2 = getN(2, 2., ca_jet[r]);
287              _h[r][M2_B1][INCL]->fill(m2_b1);
288              _h[r][M2_B2][INCL]->fill(m2_b2);
289              _h[r][N2_B1][INCL]->fill(n2_b1);
290              _h[r][N2_B2][INCL]->fill(n2_b2);
291              _h[r][M2_B1][f]->fill(m2_b1);
292              _h[r][M2_B2][f]->fill(m2_b2);
293              _h[r][N2_B1][f]->fill(n2_b1);
294              _h[r][N2_B2][f]->fill(n2_b2);
295            }
296            if (mult > 3) {
297              double n3_b1 = getN(3, 1., ca_jet[r]);
298              double n3_b2 = getN(3, 2., ca_jet[r]);
299              _h[r][N3_B1][INCL]->fill(n3_b1);
300              _h[r][N3_B2][INCL]->fill(n3_b2);
301              _h[r][N3_B1][f]->fill(n3_b1);
302              _h[r][N3_B2][f]->fill(n3_b2);
303            }
304          }
305        }
306      }
307    }
308
309
310    void finalize() {
311      for (int r = 0; r < 2; ++r) { // reconstruction (charged, all)
312        for (int o = 0; o < 33; ++o) { // observable
313          for (int f = 0; f < 4; ++f) { // flavor
314            normalize(_h[r][o][f], 1.0, false);
315          }
316        }
317      }
318    }
319
320    //@}
321
322
323  private:
324
325    double deltaR(PseudoJet j1, PseudoJet j2) {
326      double deta = j1.eta() - j2.eta();
327      double dphi = j1.delta_phi_to(j2);
328      return sqrt(deta*deta + dphi*dphi);
329    }
330
331
332    double getPtDs(PseudoJet jet) {
333      double mult   = jet.constituents().size();
334      double sumpt  = 0.; // would be jet.pt() in WTA scheme but better keep it generic
335      double sumpt2 = 0.;
336      for (auto p : jet.constituents()) {
337        sumpt  += p.pt();
338        sumpt2 += pow(p.pt(), 2);
339      }
340      double ptd = sumpt2/pow(sumpt,2);
341      return max(0., sqrt((ptd-1./mult) * mult/(mult-1.)));
342    }
343
344
345    double calcGA(double beta, double kappa, PseudoJet jet) {
346      double sumpt = 0.;
347      for (const auto& p : jet.constituents()) {
348        sumpt += p.pt();
349      }
350      double ga = 0.;
351      for (auto p : jet.constituents()) {
352        ga += pow(p.pt()/sumpt, kappa) * pow(deltaR(jet, p)/0.4, beta);
353      }
354      return ga;
355    }
356
357
358    double getEcc(PseudoJet jet) {
359      // Covariance matrix
360      Matrix<2> M;
361      for (const auto& p : jet.constituents()) {
362        Matrix<2> MPart;
363        MPart.set(0, 0, (p.eta() - jet.eta()) * (p.eta() - jet.eta()));
364        MPart.set(0, 1, (p.eta() - jet.eta()) * mapAngleMPiToPi(p.phi() - jet.phi()));
365        MPart.set(1, 0, mapAngleMPiToPi(p.phi() - jet.phi()) * (p.eta() - jet.eta()));
366        MPart.set(1, 1, mapAngleMPiToPi(p.phi() - jet.phi()) * mapAngleMPiToPi(p.phi() - jet.phi()));
367        M += MPart * p.e();
368      }
369      // Calculate eccentricity from eigenvalues
370      // Check that the matrix is symmetric.
371      const bool isSymm = M.isSymm();
372      if (!isSymm) {
373        MSG_ERROR("Error: energy tensor not symmetric:");
374        MSG_ERROR("[0,1] vs. [1,0]: " << M.get(0,1) << ", " << M.get(1,0));
375      }
376      // If not symmetric, something's wrong (we made sure the error msg appeared first).
377
378      assert(isSymm);
379      const double a = M.get(0,0);
380      const double b = M.get(1,1);
381      const double c = M.get(1,0);
382
383      const double l1 = 0.5*(a+b+sqrt( (a-b)*(a-b) + 4 *c*c));
384      const double l2 = 0.5*(a+b-sqrt( (a-b)*(a-b) + 4 *c*c));
385
386      return 1. - l2/l1;
387    }
388
389
390    double getTau(int N, int M, PseudoJet jet) {
391      fjcontrib::NsubjettinessRatio tau_ratio(N, M, fjcontrib::OnePass_WTA_KT_Axes(),
392                                              fjcontrib::NormalizedMeasure(1.0, 0.4));
393      return tau_ratio(jet);
394    }
395
396
397    vector<double> getZg(PseudoJet jet) {
398      PseudoJet jet0 = jet;
399      PseudoJet jet1, jet2;
400      double zg = 0.;
401      while (zg < 0.1 && jet0.has_parents(jet1, jet2)) {
402        zg   = jet2.pt()/jet0.pt();
403        jet0 = jet1;
404      }
405      if (zg < 0.1) return {-1., -1.};
406      return {zg, jet1.delta_R(jet2)};
407    }
408
409
410    int getNSD(double zcut, double beta, PseudoJet jet) {
411      PseudoJet jet0 = jet;
412      PseudoJet jet1, jet2;
413      int nsd = 0.;
414      double zg = 0.;
415      while (jet0.has_parents(jet1, jet2)) {
416        zg = jet2.pt()/jet0.pt();
417        if (zg > zcut * pow(jet1.delta_R(jet2)/0.4, beta))
418          nsd += 1;
419        jet0 = jet1;
420      }
421      return nsd;
422    }
423
424
425    double getC(int N, double beta, PseudoJet jet) {
426      fjcontrib::EnergyCorrelatorDoubleRatio C(N, beta);
427      return C(jet);
428    }
429
430
431    double getM(int N, double beta, PseudoJet jet) {
432      fjcontrib::EnergyCorrelatorMseries CM(N, beta);
433      return CM(jet);
434    }
435
436
437    double getN(int N, double beta, PseudoJet jet) {
438      fjcontrib::EnergyCorrelatorNseries CN(N, beta);
439      return CN(jet);
440    }
441
442
443  private:
444
445    // @name Histogram data members
446    //@{
447
448    Cut particle_cut, lepton_cut, jet_cut;
449    Histo1DPtr _h[2][33][4];
450
451    //@}
452
453  };
454
455
456
457  // The hook for the plugin system
458  RIVET_DECLARE_PLUGIN(CMS_2018_I1690148);
459
460}