rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

CMS_2013_I1224539

CMS jet mass measurement in $W, Z$ and dijet events
Experiment: CMS (LHC)
Inspire ID: 1224539
Status: VALIDATED
Authors:
  • Salvatore Rappoccio
  • Nhan Tran
  • Kalanand Mishra
  • dlopes@mail.cern.ch
References: Beams: p+ p+
Beam energies: (3500.0, 3500.0) GeV
Run details:
  • 7 TeV $pp$ collisions. Events are required to have a electron channel $Z$ with $p_\perp > 120$ GeV, and at least 1 jet opposed to the $Z$ with $p_\perp > 50$ GeV and $|y| < 2.4$.

Measurements of the mass spectra of the jets in dijet and $W/Z$+jet events from proton--proton collisions at a centre-of-mass energy of 7 TeV using a data sample of integrated luminosity of 5 fb$^-1$. The jets are reconstructed using the both the anti-$k_T$ algorithm with $R=0.7$ (AK7) and the Cambridge-Aachen algorithm with $R=0.8$ (CA8) and $R=1.2$ (CA12) with several grooming techniques applied (ungroomed, filtered, pruned and trimmed). See the text of the paper for more details. For the dijet events the distributions are presented as a function of the mean Mass of the two leading jets in bins of the mean p_\perp of the two jets.

Source code: CMS_2013_I1224539.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/DileptonFinder.hh"
  7#include "Rivet/Projections/MissingMomentum.hh"
  8#include "fastjet/tools/Filter.hh"
  9#include "fastjet/tools/Pruner.hh"
 10
 11namespace Rivet {
 12
 13
 14  /// CMS jet mass measurement in W, Z and dijet events
 15  class CMS_2013_I1224539 : public Analysis {
 16  public:
 17
 18    /// Constructor
 19    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2013_I1224539);
 20
 21
 22    /// @name Analysis methods
 23    /// @{
 24
 25    /// Book histograms and initialise projections before the run
 26    void init() {
 27
 28      // Get options
 29      WJET = true;
 30      ZJET = true;
 31      DIJET = true;
 32      if ( getOption("JMODE") == "W" ) {
 33        ZJET = false;
 34        DIJET = false;
 35      }
 36      if ( getOption("JMODE") == "Z" ) {
 37        WJET = false;
 38        DIJET = false;
 39      }
 40      if ( getOption("JMODE") == "DIJET" ) {
 41        WJET = false;
 42        ZJET = false;
 43      }
 44
 45      if (WJET) {
 46        // filling W+jet histos
 47
 48        // Find W's with pT > 120, MET > 50
 49        LeptonFinder lf(Cuts::abseta < 2.4 && Cuts::pT > 80*GeV && Cuts::abspid == PID::ELECTRON, 0.2);
 50        declare(lf, "Leptons");
 51        declare(MissingMomentum(), "MET");
 52
 53        // W+jet jet collections
 54        VetoedFinalState rfs;
 55        rfs.vetoFinalState(lf);
 56        declare(FastJets(rfs, JetAlg::ANTIKT, 0.7), "JetsAK7_wj");
 57        declare(FastJets(rfs, JetAlg::CAM, 0.8), "JetsCA8_wj");
 58        declare(FastJets(rfs, JetAlg::CAM, 1.2), "JetsCA12_wj");
 59
 60        // Histograms
 61        /// @note These are 2D histos rendered into slices
 62        const int wjetsOffset = 51;
 63        for (size_t i = 0; i < N_PT_BINS_vj; ++i) {
 64          book(_h_ungroomedJetMass_AK7_wj[i] ,wjetsOffset+i+1+0*N_PT_BINS_vj, 1, 1);
 65          book(_h_filteredJetMass_AK7_wj[i] ,wjetsOffset+i+1+1*N_PT_BINS_vj, 1, 1);
 66          book(_h_trimmedJetMass_AK7_wj[i] ,wjetsOffset+i+1+2*N_PT_BINS_vj, 1, 1);
 67          book(_h_prunedJetMass_AK7_wj[i] ,wjetsOffset+i+1+3*N_PT_BINS_vj, 1, 1);
 68          book(_h_prunedJetMass_CA8_wj[i] ,wjetsOffset+i+1+4*N_PT_BINS_vj, 1, 1);
 69          if (i > 0) book(_h_filteredJetMass_CA12_wj[i] ,wjetsOffset+i+5*N_PT_BINS_vj, 1, 1);
 70        }
 71      }
 72
 73      if (ZJET) {
 74        // filling Z+jet histos
 75
 76        // Find Zs with pT > 120 GeV
 77        DileptonFinder zfinder(91.2*GeV, 0.2, Cuts::abseta < 2.4 && Cuts::pT > 30*GeV &&
 78                               Cuts::abspid == PID::ELECTRON, Cuts::massIn(80*GeV, 100*GeV));
 79        declare(zfinder, "DileptonFinder");
 80
 81        // Z+jet jet collections
 82        VetoedFinalState rfs;
 83        rfs.vetoFinalState(zfinder);
 84        declare(FastJets(rfs, JetAlg::ANTIKT, 0.7), "JetsAK7_zj");
 85        declare(FastJets(rfs, JetAlg::CAM, 0.8), "JetsCA8_zj");
 86        declare(FastJets(rfs, JetAlg::CAM, 1.2), "JetsCA12_zj");
 87
 88        // Histograms
 89        /// @note These are 2D histos rendered into slices
 90        const int zjetsOffset = 28;
 91        for (size_t i = 0; i < N_PT_BINS_vj; ++i ) {
 92          book(_h_ungroomedJetMass_AK7_zj[i] ,zjetsOffset+i+1+0*N_PT_BINS_vj, 1, 1);
 93          book(_h_filteredJetMass_AK7_zj[i] ,zjetsOffset+i+1+1*N_PT_BINS_vj,1,1);
 94          book(_h_trimmedJetMass_AK7_zj[i] ,zjetsOffset+i+1+2*N_PT_BINS_vj,1,1);
 95          book(_h_prunedJetMass_AK7_zj[i] ,zjetsOffset+i+1+3*N_PT_BINS_vj,1,1);
 96          book(_h_prunedJetMass_CA8_zj[i] ,zjetsOffset+i+1+4*N_PT_BINS_vj,1,1);
 97          if (i > 0) book(_h_filteredJetMass_CA12_zj[i] ,zjetsOffset+i+5*N_PT_BINS_vj,1,1);
 98        }
 99
100      }
101
102      if (DIJET){
103
104        // Jet collections
105        FinalState fs;
106        declare(FastJets(fs, JetAlg::ANTIKT, 0.7), "JetsAK7");
107        declare(FastJets(fs, JetAlg::CAM, 0.8), "JetsCA8");
108        declare(FastJets(fs, JetAlg::CAM, 1.2), "JetsCA12");
109
110        // Histograms
111        for (size_t i = 0; i < N_PT_BINS_dj; ++i ) {
112          book(_h_ungroomedAvgJetMass_dj[i] ,i+1+0*N_PT_BINS_dj, 1, 1);
113          book(_h_filteredAvgJetMass_dj[i] ,i+1+1*N_PT_BINS_dj, 1, 1);
114          book(_h_trimmedAvgJetMass_dj[i] ,i+1+2*N_PT_BINS_dj, 1, 1);
115          book(_h_prunedAvgJetMass_dj[i] ,i+1+3*N_PT_BINS_dj, 1, 1);
116        }
117
118      }
119    }
120
121    bool isBackToBack_wj(const Particle& l, const P4& pmiss, const fastjet::PseudoJet& psjet) {
122      const FourMomentum w = l + pmiss;
123      /// @todo We should make FourMomentum know how to construct itself from a PseudoJet
124      const FourMomentum jmom(psjet.e(), psjet.px(), psjet.py(), psjet.pz());
125      return (deltaPhi(w, jmom) > 2.0 && deltaR(l, jmom) > 1.0 && deltaPhi(pmiss, jmom) > 0.4);
126    }
127
128    bool isBackToBack_zj(const DileptonFinder& zf, const fastjet::PseudoJet& psjet) {
129      const FourMomentum& z = zf.bosons()[0].momentum();
130      const FourMomentum& l1 = zf.constituents()[0].momentum();
131      const FourMomentum& l2 = zf.constituents()[1].momentum();
132      /// @todo We should make FourMomentum know how to construct itself from a PseudoJet
133      const FourMomentum jmom(psjet.e(), psjet.px(), psjet.py(), psjet.pz());
134      return (deltaPhi(z, jmom) > 2.0 && deltaR(l1, jmom) > 1.0 && deltaR(l2, jmom) > 1.0);
135    }
136
137
138    // Find the pT histogram bin index for value pt (in GeV), to hack a 2D histogram equivalent
139    /// @todo Use a YODA axis/finder alg when available
140    size_t findPtBin_vj(double ptJ) {
141      const double ptBins_vj[N_PT_BINS_vj+1] = { 125.0, 150.0, 220.0, 300.0, 450.0 };
142      for (size_t ibin = 0; ibin < N_PT_BINS_vj; ++ibin) {
143        if (inRange(ptJ, ptBins_vj[ibin], ptBins_vj[ibin+1])) return ibin;
144      }
145      return N_PT_BINS_vj;
146    }
147
148    // Find the pT histogram bin index for value pt (in GeV), to hack a 2D histogram equivalent
149    /// @todo Use a YODA axis/finder alg when available
150    size_t findPtBin_jj(double ptJ) {
151      const double ptBins_dj[N_PT_BINS_dj+1] = { 220.0, 300.0, 450.0, 500.0, 600.0, 800.0, 1000.0, 1500.0};
152      for (size_t ibin = 0; ibin < N_PT_BINS_dj; ++ibin) {
153        if (inRange(ptJ, ptBins_dj[ibin], ptBins_dj[ibin+1])) return ibin;
154      }
155      return N_PT_BINS_dj;
156    }
157
158
159    /// Perform the per-event analysis
160    void analyze(const Event& event) {
161
162      if (WJET) {
163        do { //< use do..while to allow breaks
164
165          // Get the W
166          /// @note Closest-matching mu+MET to mT == mW: not a great strategy!
167          const P4& pmiss = apply<MissingMom>(event, "MET").missingMom();
168          if (pmiss.pT() < 50*GeV) break;
169          const Particles& es = apply<LeptonFinder>(event, "Leptons").particles();
170          const int iefound = closestMatchIndex(es, pmiss, Kin::mT, 80.4*GeV, 50*GeV, 1000*GeV);
171
172          // Require a fairly high-pT W and charged lepton (lepton cut already applied)
173          if (iefound < 0) break;
174          const Particle& e = es[iefound];
175          const P4 pW = e.mom() + pmiss;
176          if (pW.pT() < 120*GeV) break;
177
178          // Get the pseudojets.
179          const PseudoJets psjetsCA8_wj = apply<FastJets>(event, "JetsCA8_wj").pseudojetsByPt( 50.0*GeV );
180          const PseudoJets psjetsCA12_wj = apply<FastJets>(event, "JetsCA12_wj").pseudojetsByPt( 50.0*GeV );
181
182          // AK7 jets
183          const PseudoJets psjetsAK7_wj = apply<FastJets>(event, "JetsAK7_wj").pseudojetsByPt( 50.0*GeV );
184          if (!psjetsAK7_wj.empty()) {
185            // Get the leading jet and make sure it's back-to-back with the W
186            const fastjet::PseudoJet& j0 = psjetsAK7_wj[0];
187            if (isBackToBack_wj(e, pmiss, j0)) {
188              const size_t njetBin = findPtBin_vj(j0.pt()/GeV);
189              if (njetBin < N_PT_BINS_vj) {
190                fastjet::PseudoJet filtered0 = _filter(j0);
191                fastjet::PseudoJet trimmed0 = _trimmer(j0);
192                fastjet::PseudoJet pruned0 = _pruner(j0);
193                _h_ungroomedJetMass_AK7_wj[njetBin]->fill(j0.m()/GeV);
194                _h_filteredJetMass_AK7_wj[njetBin]->fill(filtered0.m()/GeV);
195                _h_trimmedJetMass_AK7_wj[njetBin]->fill(trimmed0.m()/GeV);
196                _h_prunedJetMass_AK7_wj[njetBin]->fill(pruned0.m()/GeV);
197              }
198            }
199          }
200
201          // CA8 jets
202          if (!psjetsCA8_wj.empty()) {
203            // Get the leading jet and make sure it's back-to-back with the W
204            const fastjet::PseudoJet& j0 = psjetsCA8_wj[0];
205            if (isBackToBack_wj(e, pmiss, j0)) {
206              const size_t njetBin = findPtBin_vj(j0.pt()/GeV);
207              if (njetBin < N_PT_BINS_vj) {
208                fastjet::PseudoJet pruned0 = _pruner(j0);
209                _h_prunedJetMass_CA8_wj[njetBin]->fill(pruned0.m()/GeV);
210              }
211            }
212          }
213
214          // CA12 jets
215          if (!psjetsCA12_wj.empty()) {
216            // Get the leading jet and make sure it's back-to-back with the W
217            const fastjet::PseudoJet& j0 = psjetsCA12_wj[0];
218            if (isBackToBack_wj(e, pmiss, j0)) {
219              const size_t njetBin = findPtBin_vj(j0.pt()/GeV);
220              if (njetBin < N_PT_BINS_vj&&njetBin>0) {
221                fastjet::PseudoJet filtered0 = _filter(j0);
222                _h_filteredJetMass_CA12_wj[njetBin]->fill( filtered0.m() / GeV);
223              }
224            }
225          }
226
227          break;
228        } while (true);
229      }
230
231      if (ZJET) {
232        // Get the Z
233        const DileptonFinder& zfinder = apply<DileptonFinder>(event, "DileptonFinder");
234        if (zfinder.bosons().size() == 1) {
235          const Particle& z = zfinder.bosons()[0];
236          if (z.constituents().size() < 2) {
237            MSG_WARNING("Found a Z with less than 2 constituents.");
238            vetoEvent;
239          }
240          const Particle& l1 = z.constituents()[0];
241          const Particle& l2 = z.constituents()[1];
242          MSG_DEBUG(l1.pT() << " " << l2.pT());
243
244          // Require a high-pT Z (and constituents)
245          if (l1.pT() >= 30*GeV && l2.pT() >= 30*GeV && z.pT() >= 120*GeV) {
246
247            // AK7 jets
248            const PseudoJets& psjetsAK7_zj = apply<FastJets>(event, "JetsAK7_zj").pseudojetsByPt(50.0*GeV);
249            if (!psjetsAK7_zj.empty()) {
250              // Get the leading jet and make sure it's back-to-back with the Z
251              const fastjet::PseudoJet& j0 = psjetsAK7_zj[0];
252              if (isBackToBack_zj(zfinder, j0)) {
253                const size_t njetBin = findPtBin_vj(j0.pt()/GeV);
254                if (njetBin < N_PT_BINS_vj) {
255                  fastjet::PseudoJet filtered0 = _filter(j0);
256                  fastjet::PseudoJet trimmed0 = _trimmer(j0);
257                  fastjet::PseudoJet pruned0 = _pruner(j0);
258                  _h_ungroomedJetMass_AK7_zj[njetBin]->fill(j0.m()/GeV);
259                  _h_filteredJetMass_AK7_zj[njetBin]->fill(filtered0.m()/GeV);
260                  _h_trimmedJetMass_AK7_zj[njetBin]->fill(trimmed0.m()/GeV);
261                  _h_prunedJetMass_AK7_zj[njetBin]->fill(pruned0.m()/GeV);
262                }
263              }
264            }
265
266            // CA8 jets
267            const PseudoJets& psjetsCA8_zj = apply<FastJets>(event, "JetsCA8_zj").pseudojetsByPt(50.0*GeV);
268            if (!psjetsCA8_zj.empty()) {
269              // Get the leading jet and make sure it's back-to-back with the Z
270              const fastjet::PseudoJet& j0 = psjetsCA8_zj[0];
271              if (isBackToBack_zj(zfinder, j0)) {
272                const size_t njetBin = findPtBin_vj(j0.pt()/GeV);
273                if (njetBin < N_PT_BINS_vj) {
274                  fastjet::PseudoJet pruned0 = _pruner(j0);
275                  _h_prunedJetMass_CA8_zj[njetBin]->fill(pruned0.m()/GeV);
276                }
277              }
278            }
279
280            // CA12 jets
281            const PseudoJets& psjetsCA12_zj = apply<FastJets>(event, "JetsCA12_zj").pseudojetsByPt(50.0*GeV);
282            if (!psjetsCA12_zj.empty()) {
283              // Get the leading jet and make sure it's back-to-back with the Z
284              const fastjet::PseudoJet& j0 = psjetsCA12_zj[0];
285              if (isBackToBack_zj(zfinder, j0)) {
286                const size_t njetBin = findPtBin_vj(j0.pt()/GeV);
287                if (njetBin>0 && njetBin < N_PT_BINS_vj) {
288                  fastjet::PseudoJet filtered0 = _filter(j0);
289                  _h_filteredJetMass_CA12_zj[njetBin]->fill( filtered0.m() / GeV);
290                }
291              }
292            }
293          }
294        }
295      }
296
297      if (DIJET) {
298        // Look at events with >= 2 jets
299        const PseudoJets& psjetsAK7 = apply<FastJets>(event, "JetsAK7").pseudojetsByPt( 50.0*GeV );
300        if (psjetsAK7.size() >= 2) {
301
302          // Get the leading two jets and find their average pT
303          const fastjet::PseudoJet& j0 = psjetsAK7[0];
304          const fastjet::PseudoJet& j1 = psjetsAK7[1];
305          double ptAvg = 0.5 * (j0.pt() + j1.pt());
306
307          // Find the appropriate mean pT bin and escape if needed
308          const size_t njetBin = findPtBin_jj(ptAvg/GeV);
309          if (njetBin < N_PT_BINS_dj) {
310
311            // Now run the substructure algs...
312            fastjet::PseudoJet filtered0 = _filter(j0);
313            fastjet::PseudoJet filtered1 = _filter(j1);
314            fastjet::PseudoJet trimmed0 = _trimmer(j0);
315            fastjet::PseudoJet trimmed1 = _trimmer(j1);
316            fastjet::PseudoJet pruned0 = _pruner(j0);
317            fastjet::PseudoJet pruned1 = _pruner(j1);
318
319            // ... and fill the histograms
320            _h_ungroomedAvgJetMass_dj[njetBin]->fill(0.5*(j0.m() + j1.m())/GeV);
321            _h_filteredAvgJetMass_dj[njetBin]->fill(0.5*(filtered0.m() + filtered1.m())/GeV);
322            _h_trimmedAvgJetMass_dj[njetBin]->fill(0.5*(trimmed0.m() + trimmed1.m())/GeV);
323            _h_prunedAvgJetMass_dj[njetBin]->fill(0.5*(pruned0.m() + pruned1.m())/GeV);
324          }
325        }
326      }
327    }
328
329
330    /// Normalise histograms etc., after the run
331    void finalize() {
332
333      const double normalizationVal = 1000;
334
335      for (size_t i = 0; i < N_PT_BINS_vj; ++i) {
336        if (WJET) {
337          normalize(_h_ungroomedJetMass_AK7_wj[i], normalizationVal);
338          normalize(_h_filteredJetMass_AK7_wj[i], normalizationVal);
339          normalize(_h_trimmedJetMass_AK7_wj[i], normalizationVal);
340          normalize(_h_prunedJetMass_AK7_wj[i], normalizationVal);
341          normalize(_h_prunedJetMass_CA8_wj[i], normalizationVal);
342          if (i > 0) normalize( _h_filteredJetMass_CA12_wj[i], normalizationVal);
343        }
344        if (ZJET) {
345          normalize( _h_ungroomedJetMass_AK7_zj[i], normalizationVal);
346          normalize( _h_filteredJetMass_AK7_zj[i], normalizationVal);
347          normalize( _h_trimmedJetMass_AK7_zj[i], normalizationVal);
348          normalize( _h_prunedJetMass_AK7_zj[i], normalizationVal);
349          normalize( _h_prunedJetMass_CA8_zj[i], normalizationVal);
350          if (i > 0) normalize( _h_filteredJetMass_CA12_zj[i], normalizationVal);
351        }
352      }
353      if (DIJET) {
354        for (size_t i = 0; i < N_PT_BINS_dj; ++i) {
355          normalize(_h_ungroomedAvgJetMass_dj[i], normalizationVal);
356          normalize(_h_filteredAvgJetMass_dj[i], normalizationVal);
357          normalize(_h_trimmedAvgJetMass_dj[i], normalizationVal);
358          normalize(_h_prunedAvgJetMass_dj[i], normalizationVal);
359        }
360      }
361    }
362
363    /// @}
364
365
366  private:
367
368    bool WJET, ZJET, DIJET;
369
370    /// @name FastJet grooming tools (configured in constructor init list)
371    /// @{
372    const fastjet::Filter _filter{fastjet::Filter(fastjet::JetDefinition(fastjet::cambridge_algorithm, 0.3), fastjet::SelectorNHardest(3))};
373    const fastjet::Filter _trimmer{fastjet::Filter(fastjet::JetDefinition(fastjet::kt_algorithm, 0.2), fastjet::SelectorPtFractionMin(0.03))};
374    const fastjet::Pruner _pruner{fastjet::Pruner(fastjet::cambridge_algorithm, 0.1, 0.5)};
375    /// @}
376
377
378    /// @name Histograms
379    /// @{
380    enum BINS_vj { PT_125_150_vj=0, PT_150_220_vj, PT_220_300_vj, PT_300_450_vj, N_PT_BINS_vj };
381    // W+jet
382    Histo1DPtr _h_ungroomedJetMass_AK7_wj[N_PT_BINS_vj];
383    Histo1DPtr _h_filteredJetMass_AK7_wj[N_PT_BINS_vj];
384    Histo1DPtr _h_trimmedJetMass_AK7_wj[N_PT_BINS_vj];
385    Histo1DPtr _h_prunedJetMass_AK7_wj[N_PT_BINS_vj];
386    Histo1DPtr _h_prunedJetMass_CA8_wj[N_PT_BINS_vj];
387    Histo1DPtr _h_filteredJetMass_CA12_wj[N_PT_BINS_vj];
388    //Z+jet
389    Histo1DPtr _h_ungroomedJetMass_AK7_zj[N_PT_BINS_vj];
390    Histo1DPtr _h_filteredJetMass_AK7_zj[N_PT_BINS_vj];
391    Histo1DPtr _h_trimmedJetMass_AK7_zj[N_PT_BINS_vj];
392    Histo1DPtr _h_prunedJetMass_AK7_zj[N_PT_BINS_vj];
393    Histo1DPtr _h_prunedJetMass_CA8_zj[N_PT_BINS_vj];
394    Histo1DPtr _h_filteredJetMass_CA12_zj[N_PT_BINS_vj];
395    // DIJET
396    enum BINS_dj { PT_220_300_dj=0, PT_300_450_dj, PT_450_500_dj, PT_500_600_dj,
397                   PT_600_800_dj, PT_800_1000_dj, PT_1000_1500_dj, N_PT_BINS_dj };
398    Histo1DPtr _h_ungroomedJet0pt, _h_ungroomedJet1pt;
399    Histo1DPtr _h_ungroomedAvgJetMass_dj[N_PT_BINS_dj];
400    Histo1DPtr _h_filteredAvgJetMass_dj[N_PT_BINS_dj];
401    Histo1DPtr _h_trimmedAvgJetMass_dj[N_PT_BINS_dj];
402    Histo1DPtr _h_prunedAvgJetMass_dj[N_PT_BINS_dj];
403    /// @}
404
405  };
406
407
408  RIVET_DECLARE_PLUGIN(CMS_2013_I1224539);
409
410}