Rivet analyses referenceCMS_2021_I1932460Measurement 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:
Beam energies: (6500.0, 6500.0) GeV Run details:
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}
|