Rivet analyses referenceCMS_2017_I1497519Measurements of differential production cross sections for a Z boson in association with jets in pp collisions at 8 TeVExperiment: CMS (LHC) Inspire ID: 1497519 Status: VALIDATED Authors:
Beam energies: (4000.0, 4000.0) GeV Run details:
Cross sections for the production of a Z boson in association with jets in proton-proton collisions at a centre-of-mass energy of $\sqrt{s} = 8\,$TeV are measured using a data sample collected by the CMS experiment at the LHC corresponding to 19.6 fb$^{-1}$. Differential cross sections are presented as functions of up to three observables that describe the jet kinematics and the jet activity. Correlations between the azimuthal directions and the rapidities of the jets and the Z boson are studied in detail. The predictions of a number of multileg generators with leading or next-to-leading order accuracy are compared with the measurements. The comparison shows the importance of including multi-parton contributions in the matrix elements and the improvement in the predictions when next-to-leading order terms are included. Source code: CMS_2017_I1497519.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FastJets.hh"
4#include "Rivet/Projections/FinalState.hh"
5#include "Rivet/Projections/PromptFinalState.hh"
6#include "Rivet/Projections/DressedLeptons.hh"
7
8namespace Rivet {
9
10/// @brief Measurements of differential production cross sections for a Z boson in association with jets in pp collisions at 8 TeV
11
12 class CMS_2017_I1497519 : public Analysis {
13 private:
14 enum histIds {
15 //Normalized differential cross sections
16 kYZ, kYZmidPt, kYZhighPt, //Figs. 9a, 9b, 9c
17 kYdiff, kYdiffMidPt, kYdiffHighPt, //Figs. 10a, 10b, 10c
18 kYdiff2jZj1, kYdiff2jZj2, kYdiff2jZdijet, kYdiff2jJJ, //Figs. 11a, 11b, 11c, 12a
19 kYsum, kYsumMidPt, kYsumHighPt, //Figs. 13a, 13b, 13c
20 kYsum2jZj1, kYsum2jZj2, kYsum2jZdijet, kYsum2jJJ, //Figs. 14a, 14b, 14c, 12b
21
22 //Absolute differential cross sections. Figure numbers refer to doi:10.1007/JHEP04(2017)022
23 kNjets_exc, kNjets_inc, //Figs. 2a, 2b
24 kPtj1, kPtj2, kPtj3, kPtj4, kPtj5, //Figs. 3a, 3b, 4a, 4b, 5
25 kYj1, kYj2, kYj3, kYj4, kYj5, //Figs. 6a, 6b, 7a, 7b, 8
26 kHt1, kHt2, kHt3, kHt4, kHt5, //Figs. 15a, 15b, 16a, 16b, 17
27 kDphiZj1, kDphiZj1MidPt, kDphiZj1HighPt, // Figs. 18a, 20a, 22a
28 kDphi2jZj1, kDphi2jZj1MidPt, kDphi2jZj1HighPt, //Figs. 18b, 20b, 22b
29 kDphi3jZj1, kDphi3jZj1MidPt, kDphi3jZj1HighPt,// Figs. 18c, 20c, 22c
30 kDphi3jZj2, kDphi3jZj2MidPt, kDphi3jZj2HighPt, // Figs. 19a, 21a, 23a
31 kDphi3jZj3, kDphi3jZj3MidPt, kDphi3jZj3HighPt, // Figs. 19b, 21b, 23b
32 kDphi3jZj1HighHT, kDphi3jZj2HighHT, kDphi3jZj3HighHT, // Figs. 24a, 24b, 24c
33 kDphi3jJ1j2, kDphi3jJ1j2MidPt, kDphi3jJ1j2HighPt, //Figs. 25a, 26a, 27a
34 kDphi3jJ1j3, kDphi3jJ1j3MidPt, kDphi3jJ1j3HighPt, //Figs. 25b, 26b, 27b
35 kDphi3jJ2j3, kDphi3jJ2j3MidPt, kDphi3jJ2j3HighPt, //Figs. 25c, 26c, 27c
36 kDijetMass, //Fig. 28
37 kPtjYjBin0, kPtjYjBin1, kPtjYjBin2, kPtjYjBin3, kPtjYjBin4, kPtjYjBin5, kPtjYjBin6, //Fig. 29
38 kYjYZBin0, kYjYZBin1, kYjYZBin2, kYjYZBin3, //Fig. 33
39 kPtjYjYZBin0, kPtjYjYZBin1, kPtjYjYZBin2, //Fig. 37 top left frame (Z and j on same side)
40 kPtjYjYZBin3, kPtjYjYZBin4, kPtjYjYZBin5, //Fig. 37 bottom left frame (same side)
41 kPtjYjYZBin6, kPtjYjYZBin7, kPtjYjYZBin8, //Fig. 37 top right frame (opposite sides)
42 kPtjYjYZBin9, kPtjYjYZBin10, kPtjYjYZBin11, //Fig. 37 bottom right frame (opposite sides)
43 nHistos
44 };
45
46 unsigned nNormalized = kYsum2jJJ + 1;
47
48 template<typename T>
49 inline double ydiff(const ParticleBase& p1, const T& p2){ return 0.5*fabs(p1.rapidity()-p2.rapidity()); }
50
51 template<typename T>
52 inline double ysum(const ParticleBase& p1, const T& p2){ return 0.5*fabs(p1.rapidity()+p2.rapidity()); }
53
54 /// Fill histograms _h[ih+0..2] with x with different lower thresholds on tocut:
55 /// no threshold, threshold cut2, threshold cut2
56 void fill3cuts(int ih, double tocut, double cut2, double cut3, double x){
57 _h[ih]->fill(x);
58 if (tocut > cut2) _h[ih+1]->fill(x);
59 if (tocut > cut3) _h[ih+2]->fill(x);
60 }
61
62 /// Fills a set of N 1-D histograms _h[ih + 0...N-1] that represents a 2-D distribution
63 /// @param bins: boundaries of the y bins covered by each 1-D histogram
64 void fill2D(int ih, std::vector<double>& ybins, double x, double y, double w){
65 int iybin = -1;
66 double lowEdge;
67 for(auto highEdge: ybins){
68 if(y < highEdge){
69 if (iybin >= 0) _h[ih + iybin]->fill(x, w / (highEdge - lowEdge));
70 break;
71 }
72 lowEdge = highEdge;
73 ++iybin;
74 }
75 }
76
77 void fill3D(int ih, std::vector<double>& ybins, std::vector<double> zbins, double x, double y, double z, double w){
78 int izbin = -1;
79 double lowEdge = 0;
80 int nybins = ybins.size() - 1;
81 for(auto highEdge: zbins){
82 if(z < highEdge){
83 if (izbin >= 0) fill2D(ih + izbin*nybins, ybins, x, y, w / (highEdge - lowEdge));
84 break;
85 }
86 lowEdge = highEdge;
87 ++izbin;
88 }
89 }
90
91 public:
92
93 /// Constructor
94 CMS_2017_I1497519()
95 : Analysis("CMS_2017_I1497519")
96 {
97 _histListInPaperOrder = { /*number NN in comment is the id from the histogram name dNN-x01-y01*/
98 /*1*/kNjets_exc, /*2*/ kNjets_inc, //Figs. 2a, 2b
99 /*3*/kPtj1, /*4*/kPtj2, /*5*/kPtj3, /*6*/kPtj4, /*7*/kPtj5, //Figs. 3a, 3b, 4a, 4b, 5
100 /*8*/kYj1, /*9*/kYj2, /*10*/kYj3, /*11*/kYj4, /*12*/kYj5, //Figs. 6a, 6b, 7a, 7b, 8
101 /*13*/kYZ, /*14*/kYZmidPt, /*15*/kYZhighPt, //Figs. 9a, 9b, 9c
102 /*16*/kYdiff, /*17*/kYdiffMidPt, /*18*/kYdiffHighPt, //Figs. 10a, 10b, 10c
103 /*19*/kYdiff2jZj1, /*20*/kYdiff2jZj2, /*21*/kYdiff2jZdijet, /*22*/kYdiff2jJJ, //Figs. 11a, 11b, 11c, 12a
104 /*23*/kYsum2jJJ, //12b
105 /*24*/kYsum, /*25*/kYsumMidPt, /*26*/kYsumHighPt, //Figs. 13a, 13b, 13c
106 /*27*/kYsum2jZj1, /*28*/kYsum2jZj2, /*29*/kYsum2jZdijet, //Figs. 14a, 14b, 14c
107 /*30*/kHt1, /*31*/kHt2, /*32*/kHt3, /*33*/kHt4, /*34*/kHt5, //Figs. 15a, 15b, 16a, 16b, 17
108 /*35*/kDphiZj1, /*36*/kDphi2jZj1, /*37*/kDphi3jZj1, // Figs. 18a, 18b, 18c
109 /*38*/kDphi3jZj2, /*39*/kDphi3jZj3, // Figs. 19a, 19b
110 /*40*/kDphiZj1MidPt, /*41*/kDphi2jZj1MidPt, /*42*/kDphi3jZj1MidPt, // Figs. 20a, 20b, 20c
111 /*43*/kDphi3jZj2MidPt, /*44*/kDphi3jZj3MidPt, // Figs. 21a, 21b
112 /*45*/kDphiZj1HighPt, /*46*/kDphi2jZj1HighPt, /*47*/kDphi3jZj1HighPt,// Figs. 22a, 22b, 22c
113 /*48*/kDphi3jZj2HighPt, /*49*/kDphi3jZj3HighPt, // Figs. 23a, 23b
114 /*50*/kDphi3jZj1HighHT, /*51*/kDphi3jZj2HighHT, /*52*/kDphi3jZj3HighHT, // Figs. 24a, 24b, 24c
115 /*53*/kDphi3jJ1j2, /*54*/kDphi3jJ1j3, /*55*/kDphi3jJ2j3, //Figs. 25a, 25b, 25c
116 /*56*/kDphi3jJ1j2MidPt, /*57*/kDphi3jJ1j3MidPt, /*58*/kDphi3jJ2j3MidPt, //Figs. 26a, 26b, 26c
117 /*59*/kDphi3jJ1j2HighPt, /*60*/kDphi3jJ1j3HighPt, /*61*/kDphi3jJ2j3HighPt, //Figs. 27a, 27b, 27c
118 /*62*/kDijetMass, //Fig. 28
119 /*63*/kPtjYjBin0, /*64*/kPtjYjBin1, /*65*/kPtjYjBin2, /*66*/kPtjYjBin3, /*67*/kPtjYjBin4, /*68*/kPtjYjBin5, /*69*/kPtjYjBin6, //Fig. 29
120 /*70*/kYjYZBin0, /*71*/kYjYZBin1, /*72*/kYjYZBin2, /*73*/kYjYZBin3, //Fig. 33
121 /*74*/kPtjYjYZBin0, /*75*/kPtjYjYZBin1, /*76*/kPtjYjYZBin2, //Fig. 37 top left frame (Z and j on same side)
122 /*77*/kPtjYjYZBin6, /*78*/kPtjYjYZBin7, /*79*/kPtjYjYZBin8, //Fig. 37 top right frame (opposite sides)
123 /*80*/kPtjYjYZBin3, /*81*/kPtjYjYZBin4, /*82*/kPtjYjYZBin5, //Fig. 37 bottom left frame (same side)
124 /*83*/kPtjYjYZBin9, /*84*/kPtjYjYZBin10, /*85*/kPtjYjYZBin11, //Fig. 37 bottom right frame (opposite sides)
125 };
126 }
127
128
129 /// Book histograms and initialise projections before the run
130 void init() {
131
132 // Get options from the new option system
133 // default to combined.
134 _mode = 2;
135 if ( getOption("LMODE") == "EL" ) _mode = 0;
136 if ( getOption("LMODE") == "MU" ) _mode = 1;
137 if ( getOption("LMODE") == "EMU" ) _mode = 2;
138
139 FinalState fs;
140
141 PromptFinalState bareMuons(Cuts::abspid == PID::MUON);
142 declare(DressedLeptons(fs, bareMuons, /*dRmax = */0.1,
143 Cuts::pT > 20*GeV && Cuts::abseta < 2.4,
144 /*useDecayPhotons = */ true),
145 "muons");
146
147 PromptFinalState bareElectrons(Cuts::abspid == PID::ELECTRON);
148 declare(DressedLeptons(fs, bareElectrons, /*dRmax =*/ 0.1,
149 Cuts::pT > 20*GeV && Cuts::abseta < 2.4,
150 /*useDecayPhotons = */ true),
151 "electrons");
152
153 FastJets jets(fs, FastJets::ANTIKT, 0.5);
154 declare(jets, "jets");
155
156 _h = std::vector<Histo1DPtr>(nHistos);
157 for(int ih = 0; ih < nHistos; ++ih){
158 book(_h[_histListInPaperOrder[ih]], ih + 1, 1, 1);
159 }
160
161 _ptjYjBins = {0, 0.5, 1., 1.5, 2., 2.5, 3.2, 4.7};
162 _yJyZbins = {0, 0.5, 1., 1.5, 2.5};
163 _ptJyJyZbinsYj = {0., 1.5, 2.5, 4.7};
164 _ptJyJyZbinsYZ = {0., 1., 2.5};
165 }
166
167 /// Z boson finder.
168 /// Note: we don't use the standard ZFinder class in order to stick to
169 /// the definition of the publication that is simpler than the ZFinder
170 /// algorithm
171 /// @param leptons pt-ordered of electron or muon collection to use to build
172 /// the Z boson
173 std::unique_ptr<Particle> zfinder(const std::vector<DressedLepton>& leptons){
174 if(leptons.size() < 2) return 0;
175 if(leptons[0].charge()*leptons[1].charge() > 0) return 0;
176 std::unique_ptr<Particle> cand(new Particle(PID::ZBOSON, leptons[0].mom()
177 + leptons[1].mom()));
178 if (cand->mass() < 71.*GeV || cand->mass() > 111.*GeV) return 0;
179 return cand;
180 }
181
182 /// Perform the per-event analysis
183 void analyze(const Event& event) {
184
185 std::vector<DressedLepton> muons = apply<DressedLeptons>(event, "muons").dressedLeptons();
186 std::vector<DressedLepton> electrons = apply<DressedLeptons>(event, "electrons").dressedLeptons();
187
188 //Look for Z->ee
189 std::unique_ptr<Particle> z = zfinder(electrons);
190
191 const std::vector<DressedLepton>* dressedLeptons = 0;
192
193 //Look for Z->ee
194 if(z.get() != nullptr && _mode != 1) {
195 dressedLeptons = &electrons;
196 } else{ //look for Z->mumu
197 z = zfinder(muons);
198 if(z.get() != nullptr && _mode != 0){
199 dressedLeptons = &muons;
200 } else{ //no Z boson found
201 vetoEvent;
202 }
203 }
204
205 // Cluster jets
206 const FastJets& fj = apply<FastJets>(event, "jets");
207 const Jets& jets = fj.jetsByPt(Cuts::absrap < 4.7 && Cuts::pT > 30*GeV);
208
209 // Remove jets overlapping with any of the two selected leptons
210 Jets goodjets47 = filter_discard(jets, [dressedLeptons](const ParticleBase& j){
211 return deltaR(j, (*dressedLeptons)[0]) < 0.5
212 || deltaR(j, (*dressedLeptons)[1]) < 0.5;
213 });
214
215 // Jets in the CMS tracker acceptance
216 Jets goodjets24 = filter_select(goodjets47, [](const ParticleBase& j){
217 return j.absrapidity() < 2.4;
218 });
219
220 goodjets24 = sortByPt(goodjets24);
221
222 // Compute jet pt scalar sum, H_T:
223 double ht = sum(goodjets24, Kin::pT, 0.);
224
225 _h[kNjets_exc]->fill(goodjets24.size());
226
227 // Fill jet number integral histograms
228 /// @todo Could be better computed by toIntegral transform on exclusive histo
229 for (size_t iJet = 0; iJet <= goodjets24.size(); iJet++ )
230 _h[kNjets_inc]->fill(iJet);
231
232 int offset = 0;
233 for(const auto& j: goodjets24){
234 int nj = 1 + offset;
235 if(nj > 5) break;
236
237 _h[kPtj1 + offset]->fill(j.pT() / GeV);
238 _h[kYj1 + offset]->fill(j.absrapidity());
239 _h[kHt1 + offset]->fill(ht);
240 ++offset;
241 }
242
243 ///////////////////////////////
244 /// Nj >=1 in |y| < 4.7
245 if (goodjets47.size() < 1) return;
246
247 const Jet& j1_47 = goodjets47[0];
248
249 //note: 0.5 factors in the four following fill statements is required
250 //because the cross section is differiated in y(j1) while the binning
251 //is in abs(y(j1)): bin filled twice for the same y(j1) value.
252 //An extra factor 0.5 is included for differential cross including jet and
253 //Z boson rapidities to matches the definition used in the publication.
254 fill2D(kPtjYjBin0, _ptjYjBins, j1_47.pt(), j1_47.absrapidity(), 0.5);
255 fill2D(kYjYZBin0, _yJyZbins, j1_47.absrapidity()*sign(z->rapidity()*j1_47.rapidity()), z->absrapidity(), 0.25);
256
257 if(j1_47.rapidity()*z->rapidity() >= 0){
258 fill3D(kPtjYjYZBin0, _ptJyJyZbinsYj, _ptJyJyZbinsYZ, j1_47.pt(), j1_47.absrapidity(), z->absrapidity(), 0.25);
259 } else{
260 fill3D(kPtjYjYZBin6, _ptJyJyZbinsYj, _ptJyJyZbinsYZ, j1_47.pt(), j1_47.absrapidity(), z->absrapidity(), 0.25);
261 }
262 ///////////////////////////////
263
264 ////////////////////////////////
265 /// Nj >= 1 in |y| < 2.4
266 if (goodjets24.size() < 1) return;
267
268 const Jet& j1 = goodjets24[0];
269
270 fill3cuts(kYZ, z->pt(), 150*GeV, 300*GeV, z->absrapidity());
271
272 fill3cuts(kYdiff, z->pt(), 150*GeV, 300*GeV, ydiff(*z, j1));
273 fill3cuts(kYsum, z->pt(), 150*GeV, 300*GeV, ysum(*z, j1));
274
275 fill3cuts(kDphiZj1, z->pt(), 150*GeV, 300*GeV, deltaPhi(*z, j1));
276 //////////////////////////////
277
278 ////////////////////////////////
279 /// Nj >= 2 in |y| < 2.4
280 if (goodjets24.size() < 2) return;
281
282 const Jet& j2 = goodjets24[1];
283 _h[kYdiff2jZj1]->fill(ydiff(*z, j1));
284 _h[kYdiff2jZj2]->fill(ydiff(*z, j2));
285 _h[kYdiff2jZdijet]->fill(ydiff(*z, j1.mom() + j2.mom()));
286 _h[kYdiff2jJJ]->fill(ydiff(j1, j2));
287
288 _h[kYsum2jZj1]->fill(ysum(*z, j1));
289 _h[kYsum2jZj2]->fill(ysum(*z, j2));
290 _h[kYsum2jZdijet]->fill(ysum(*z, j1.mom() + j2.mom()));
291 _h[kYsum2jJJ]->fill(ysum(j1, j2));
292
293 fill3cuts(kDphi2jZj1, z->pt(), 150*GeV, 300*GeV, deltaPhi(*z, j1));
294 _h[kDijetMass]->fill((j1.mom() + j2.mom()).mass());
295 //////////////////////////////
296
297 ////////////////////////////////
298 /// Nj >= 3 in |y| < 2.4
299 if (goodjets24.size() < 3) return;
300
301 const Jet& j3 = goodjets24[2];
302
303 fill3cuts(kDphi3jZj1, z->pt(), 150*GeV, 300*GeV, deltaPhi(*z, j1));
304 fill3cuts(kDphi3jZj2, z->pt(), 150*GeV, 300*GeV, deltaPhi(*z, j2));
305 fill3cuts(kDphi3jZj3, z->pt(), 150*GeV, 300*GeV, deltaPhi(*z, j3));
306 fill3cuts(kDphi3jJ1j2, z->pt(), 150*GeV, 300*GeV, deltaPhi(j1, j2));
307 fill3cuts(kDphi3jJ1j3, z->pt(), 150*GeV, 300*GeV, deltaPhi(j1, j3));
308 fill3cuts(kDphi3jJ2j3, z->pt(), 150*GeV, 300*GeV, deltaPhi(j2, j3));
309
310 //Although not specified in the paper, the H^{jet}_{T} varible used in
311 //Fig. 24 differs from H_T and is defined as the scalar sum
312 //of the pt of the three leading jets
313 double ht_3jets = goodjets24[0].pt() + goodjets24[1].pt() + goodjets24[2].pt();
314 if(z->pt() > 150*GeV && ht_3jets > 300*GeV){
315 _h[kDphi3jZj1HighHT]->fill(deltaPhi(*z, j1));
316 _h[kDphi3jZj2HighHT]->fill(deltaPhi(*z, j2));
317 _h[kDphi3jZj3HighHT]->fill(deltaPhi(*z, j3));
318 }
319 //////////////////////////////
320 }
321
322
323 /// Normalise histograms etc., after the run
324 void finalize() {
325
326 double norm = (sumOfWeights() != 0) ? crossSection()/sumOfWeights() : 1.0;
327
328 // when running in combined mode, need to average to get lepton xsec
329 if (_mode == 2) norm /= 2.;
330
331 MSG_DEBUG("Cross section = " << std::setfill(' ') << std::setw(14) << std::fixed << std::setprecision(3) << crossSection() << " pb");
332 MSG_DEBUG("# Events = " << std::setfill(' ') << std::setw(14) << std::fixed << std::setprecision(3) << numEvents() );
333 MSG_DEBUG("SumW = " << std::setfill(' ') << std::setw(14) << std::fixed << std::setprecision(3) << sumOfWeights());
334 MSG_DEBUG("Norm factor = " << std::setfill(' ') << std::setw(14) << std::fixed << std::setprecision(6) << norm);
335
336 unsigned ih = 0;
337 for(auto& h: _h){
338 if(ih < nNormalized) normalize(h);
339 else scale(h, norm);
340 ++ih;
341 }
342
343 }
344
345 protected:
346
347 size_t _mode;
348
349 private:
350
351 /// Histograms
352 std::vector<Histo1DPtr> _h;
353
354 /// List of histogram in the order of appearance
355 /// in the paper. Histograms are numbered according
356 /// to this order.
357 std::vector<enum histIds> _histListInPaperOrder;
358
359 //@name Binning of pseudo-2D/3D histograms
360 std::vector<double> _ptjYjBins;
361 std::vector<double> _yJyZbins;
362 std::vector<double> _ptJyJyZbinsYj;
363 std::vector<double> _ptJyJyZbinsYZ;
364
365 };
366
367 RIVET_DECLARE_PLUGIN(CMS_2017_I1497519);
368}
|