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