rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

CMS_2017_I1497519

Measurements of differential production cross sections for a Z boson in association with jets in pp collisions at 8 TeV
Experiment: CMS (LHC)
Inspire ID: 1497519
Status: VALIDATED
Authors:
  • cms-pag-conveners-smp@cern.ch
  • Philippe Gras
References: Beams: p+ p+
Beam energies: (4000.0, 4000.0) GeV
Run details:
  • Run MC generators with Z decaying leptonically into both electrons and muons at 8 TeV CoM energy. If only one of the two decay channels is included, multiply the provided cross-section value by two for a proper normalisation.

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}