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