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: 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/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}