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