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