rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

CMS_2021_I1932460

Measurement of double-parton scattering in inclusive production of four jets with low transverse momentum in proton-proton collisions at $\sqrt{s}$ = 13 TeV.
Experiment: CMS (LHC)
Inspire ID: 1932460
Status: VALIDATED
Authors:
  • cms-pag-conveners-smp@cern.ch
  • Maxim Pieters
  • Hans Van Haevermaet
  • Pierre Van Mechelen
References: Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • Event types Hard QCD; Cuts jet $p_\perp > 20$ GeV, $|\eta|$ < 4.7.

A measurement of inclusive four-jet production in proton-proton collisions at a center-of-mass energy of 13 TeV is presented. The transverse momenta of jets within $| \eta | < 4.7$ reach down to 35, 30, 25, and 20 GeV for the first-, second-, third-, and fourth- leading jet, respectively. Differential cross sections are measured as functions of the jet transverse momentum, jet pseudorapidity, and several other observables that describe the angular correlations between the jets. The measured distributions show sensitivity to different aspects of the underlying event, parton shower, and matrix element calculations. In particular, the interplay between angular correlations caused by parton shower and double-parton scattering contributions is shown to be important. The double-parton scattering contribution is extracted by means of a template fit to the data, using distributions for single-parton scattering obtained from Monte Carlo event generators and a double-parton scattering distribution constructed from inclusive single-jet events in data. The effective double-parton scattering cross section is calculated and discussed in view of previous measurements and of its dependence on the models used to describe the single-parton scattering background.

Source code: CMS_2021_I1932460.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/PromptFinalState.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// Double-parton scattering in inclusive production of four jets with low pT in 13 TeV pp
 11  class CMS_2021_I1932460 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2021_I1932460);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23
 24      // Initialise and register projections
 25
 26      // the basic final-state projection:
 27      // all final-state particles within
 28      // the given eta acceptance: 4.7 (jet range) + 0.4 (cone size)
 29      const FinalState fs(Cuts::abseta < 5.1);
 30
 31      // the final-state particles declared above are clustered using FastJet with
 32      // the anti-kT algorithm and a jet-radius parameter 0.4
 33      // muons and neutrinos are excluded from the clustering
 34      FastJets jetfs(fs, JetAlg::ANTIKT, 0.4, JetMuons::NONE, JetInvisibles::NONE);
 35      declare(jetfs, "jets");
 36
 37      // Book histograms
 38      book(_h["JetPt1"], 1, 1, 1);
 39      book(_h["JetPt2"], 2, 1, 1);
 40      book(_h["JetPt3"], 3, 1, 1);
 41      book(_h["JetPt4"], 4, 1, 1);
 42      book(_h["JetEta1"], 5, 1, 1);
 43      book(_h["JetEta2"], 6, 1, 1);
 44      book(_h["JetEta3"], 7, 1, 1);
 45      book(_h["JetEta4"], 8, 1, 1);
 46
 47      book(_h["DeltaPhiSoft_binNorm"], 9, 1, 1);
 48      book(_h["DeltaPhi3_binNorm"], 10, 1, 1);
 49      book(_h["DeltaY_binNorm"], 11, 1, 1);
 50      book(_h["DeltaPhiY_binNorm"], 12, 1, 1);
 51      book(_h["DeltaPtSoft_binNorm"], 13, 1, 1);
 52      book(_h["DeltaS_binNorm"], 14, 1, 1);
 53
 54      book(_h["DeltaPhiSoft"], 45, 1, 1);
 55      book(_h["DeltaPhi3"], 46, 1, 1);
 56      book(_h["DeltaY"], 47, 1, 1);
 57      book(_h["DeltaPhiY"], 48, 1, 1);
 58      book(_h["DeltaPtSoft"], 49, 1, 1);
 59      book(_h["DeltaS"], 50, 1, 1);
 60
 61    }
 62
 63
 64    /// Perform the per-event analysis
 65    void analyze(const Event& event) {
 66
 67      // retrieve clustered jets, sorted by pT, with a minimum pT cut 10 GeV and eta range 4.7 (similar to PFJet collection)
 68      Jets jets = apply<FastJets>(event, "jets").jetsByPt(Cuts::abseta < 4.7 && Cuts::pT > 10*GeV);
 69
 70      // fill only if there are at least 4 jets
 71      if (jets.size() < 4) vetoEvent;
 72
 73      double pt0 = jets[0].pt();
 74      double pt1 = jets[1].pt();
 75      double pt2 = jets[2].pt();
 76      double pt3 = jets[3].pt();
 77
 78      // pt selection
 79      if (pt0 < 35.0 || pt1 < 30.0 || pt2 < 25.0 || pt3 < 20.0) vetoEvent;
 80
 81      double phi0 = jets[0].phi();
 82      double phi1 = jets[1].phi();
 83      double phi2 = jets[2].phi();
 84      double phi3 = jets[3].phi();
 85
 86      _h["JetPt1"]->fill(pt0);
 87      _h["JetPt2"]->fill(pt1);
 88      _h["JetPt3"]->fill(pt2);
 89      _h["JetPt4"]->fill(pt3);
 90
 91      _h["JetEta1"]->fill(jets[0].eta());
 92      _h["JetEta2"]->fill(jets[1].eta());
 93      _h["JetEta3"]->fill(jets[2].eta());
 94      _h["JetEta4"]->fill(jets[3].eta());
 95
 96      // delta phi and eta of the 2 soft jets
 97      _h["DeltaPhiSoft"]->fill(abs(deltaPhi(phi2, phi3)));
 98      _h["DeltaPhiSoft_binNorm"]->fill(abs(deltaPhi(phi2, phi3)));
 99
100      // delta pt between soft jets
101      double DptSoft = sqrt(pow(pt2*cos(phi2) + pt3*cos(phi3), 2) + pow(pt2*sin(phi2) + pt3*sin(phi3), 2))/(pt2 + pt3);
102      _h["DeltaPtSoft"]->fill(DptSoft);
103      _h["DeltaPtSoft_binNorm"]->fill(DptSoft);
104
105      // delta S
106      if (pt0 > 50.0 && pt1 > 30.0 && pt2 > 30.0 && pt3 > 30.0) {
107        double phiH = atan2(pt0*sin(phi0) + pt1*sin(phi1) , pt0*cos(phi0) + pt1*cos(phi1));
108        double phiS = atan2(pt2*sin(phi2) + pt3*sin(phi3) , pt2*cos(phi2) + pt3*cos(phi3));
109        double DS = abs(deltaPhi(phiH, phiS));
110        _h["DeltaS"]->fill(DS);
111        _h["DeltaS_binNorm"]->fill(DS);
112      }
113
114      // delta Y: most remote jets in rapidity, find min & max eta
115      double mineta = 99999;
116      double maxeta = -99999;
117      int minetapos = -1;
118      int maxetapos = -1;
119
120      for (int i = 0; i < 4; ++i) {
121        if (jets[i].eta() < mineta) {
122          mineta = jets[i].eta();
123          minetapos = i;
124        }
125        if (jets[i].eta() > maxeta) {
126          maxeta = jets[i].eta();
127          maxetapos = i;
128        }
129      }
130
131      _h["DeltaY"]->fill(abs(jets[minetapos].eta() - jets[maxetapos].eta()));
132      _h["DeltaY_binNorm"]->fill(abs(jets[minetapos].eta() - jets[maxetapos].eta()));
133
134      // Delta phi Y: azimuthal angle between most remote jets in eta
135      _h["DeltaPhiY"]->fill(abs(deltaPhi(jets[minetapos].phi(), jets[maxetapos].phi())));
136      _h["DeltaPhiY_binNorm"]->fill(abs(deltaPhi(jets[minetapos].phi(), jets[maxetapos].phi())));
137
138      // delta phi3
139      double minphi3 = 999;
140      for (int iphi1 = 0; iphi1 < 4; ++iphi1) {
141        for (int iphi2 = 0; iphi2 < 4; ++iphi2) {
142          for (int iphi3 = 0; iphi3 < 4; ++iphi3) {
143            if ( !(iphi1 == iphi2 || iphi2 == iphi3 || iphi1 == iphi3) ) {
144              double temp_phi1= jets[iphi1].phi();
145              double temp_phi2= jets[iphi2].phi();
146              double temp_phi3= jets[iphi3].phi();
147              double temp_minphi3 = abs(deltaPhi(temp_phi1, temp_phi2)) + abs(deltaPhi(temp_phi2, temp_phi3));
148              if (temp_minphi3 < minphi3) minphi3 = temp_minphi3;
149            }
150          }
151        }
152      }
153
154      _h["DeltaPhi3"]->fill(minphi3);
155      _h["DeltaPhi3_binNorm"]->fill(minphi3);
156
157    }
158
159
160    /// Normalise histograms etc., after the run
161    void finalize() {
162
163      scale(_h["JetPt1"], crossSection()/picobarn/sumOfWeights()); // norm to cross section
164      scale(_h["JetPt2"], crossSection()/picobarn/sumOfWeights());
165      scale(_h["JetPt3"], crossSection()/picobarn/sumOfWeights());
166      scale(_h["JetPt4"], crossSection()/picobarn/sumOfWeights());
167      scale(_h["JetEta1"], crossSection()/picobarn/sumOfWeights());
168      scale(_h["JetEta2"], crossSection()/picobarn/sumOfWeights());
169      scale(_h["JetEta3"], crossSection()/picobarn/sumOfWeights());
170      scale(_h["JetEta4"], crossSection()/picobarn/sumOfWeights());
171      scale(_h["DeltaPhiSoft"], crossSection()/picobarn/sumOfWeights());
172      scale(_h["DeltaPhi3"], crossSection()/picobarn/sumOfWeights());
173      scale(_h["DeltaY"], crossSection()/picobarn/sumOfWeights());
174      scale(_h["DeltaPhiY"], crossSection()/picobarn/sumOfWeights());
175      scale(_h["DeltaPtSoft"], crossSection()/picobarn/sumOfWeights());
176      scale(_h["DeltaS"], crossSection()/picobarn/sumOfWeights());
177
178      // create bin normalised histograms
179
180      // Correct for binwidths: Rivet automatically normalises histograms to binwidth when plotting AFTER the normalisation executed here.
181      // So we must calculate an extra correction here so that finally our bin-normalised histograms end up around 1 as in the paper.
182      // in YODA bin index starts from 0
183
184      // For DeltaY, which has variable binwidths we need to do following steps
185      // divide histograms by binwidth
186      for (unsigned int i = 1; i < _h["DeltaY_binNorm"]->numBins()+1; ++i) {
187        _h["DeltaY_binNorm"]->bin(i).scaleW(1.0/_h["DeltaY_binNorm"]->bin(i).xWidth());
188      }
189
190      // normalise to average of first 4 bins
191      scale(_h["DeltaY_binNorm"], 1.0/(_h["DeltaY_binNorm"]->integralRange(1,4)/4.0));
192
193      // multiply again with binwidth
194      for (unsigned int i = 1; i < _h["DeltaY_binNorm"]->numBins()+1; ++i) {
195        _h["DeltaY_binNorm"]->bin(i).scaleW(_h["DeltaY_binNorm"]->bin(i).xWidth());
196      }
197
198      // DeltaPhiSoft and DeltaPhi3 histograms have uniform binwidths, so multiply with first binwidth is sufficient
199      scale(_h["DeltaPhiSoft_binNorm"], _h["DeltaPhiSoft_binNorm"]->bin(1).xWidth()/(_h["DeltaPhiSoft_binNorm"]->integralRange(1,5)/5.0));
200      scale(_h["DeltaPhi3_binNorm"], _h["DeltaPhi3_binNorm"]->bin(1).xWidth()/(_h["DeltaPhi3_binNorm"]->integralRange(1,4)/4.0));
201
202      // DeltaPhiY, DeltaPtSoft and DeltaS are normalised to last bin
203      scale(_h["DeltaPhiY_binNorm"], _h["DeltaPhiY_binNorm"]->bin(12).xWidth()/_h["DeltaPhiY_binNorm"]->bin(12).sumW() );
204      scale(_h["DeltaPtSoft_binNorm"], _h["DeltaPtSoft_binNorm"]->bin(8).xWidth()/_h["DeltaPtSoft_binNorm"]->bin(8).sumW() );
205      scale(_h["DeltaS_binNorm"], _h["DeltaS_binNorm"]->bin(7).xWidth()/_h["DeltaS_binNorm"]->bin(7).sumW() );
206    }
207
208    /// @}
209
210    map<string, Histo1DPtr> _h;
211
212  };
213
214
215  RIVET_DECLARE_PLUGIN(CMS_2021_I1932460);
216
217}