Rivet analyses referenceCMS_2013_I1224539CMS jet mass measurement in $W, Z$ and dijet eventsExperiment: CMS (LHC) Inspire ID: 1224539 Status: VALIDATED Authors:
Beam energies: (3500.0, 3500.0) GeV Run details:
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}
|