Measurements of jet multiplicity and jet transverse momentum in multijet events in proton–proton collisions at ${\sqrt{s}=13\, \text {TeV}}$

Abstract

Multijet events at large transverse momentum ($p_{\textrm{T}}$) are measured at $\sqrt{s}=13\,\text {TeV}$ using data recorded with the CMS detector at the LHC, corresponding to an integrated luminosity of $36.3{\,\text {fb}^{-1}}$. The multiplicity of jets with $p_{\textrm{T}} >50\,\text {GeV}$ that are produced in association with a high-$p_{\textrm{T}}$ dijet system is measured in various ranges of the $p_{\textrm{T}}$ of the jet with the highest transverse momentum and as a function of the azimuthal angle difference $\varDelta \phi _{1,2}$ between the two highest $p_{\textrm{T}}$ jets in the dijet system. The differential production cross sections are measured as a function of the transverse momenta of the four highest $p_{\textrm{T}}$ jets. The measurements are compared with leading and next-to-leading order matrix element calculations supplemented with simulations of parton shower, hadronization, and multiparton interactions. In addition, the measurements are compared with next-to-leading order matrix element calculations combined with transverse-momentum dependent parton densities and transverse-momentum dependent parton shower.

1 Introduction

The production of jets, which are reconstructed from a stream of hadrons coming from the fragmentation of energetic partons, is described by the theory of strong interactions, quantum chromodynamics (QCD). In proton–proton ($\text {p}\text {p}$) collisions, at leading order (LO) in the strong coupling $\alpha _\textrm{S}$, two colliding partons from the incident protons scatter and produce two high transverse-momentum ($p_{\textrm{T}}$) partons in the final state. The jets that originate from such a process are strongly correlated in the transverse plane, and the azimuthal angle difference between them, $\varDelta \phi _{1,2}$, should be close to $\pi$. However, higher-order corrections to the lowest order process will result in a decorrelation in the azimuthal plane, and $\varDelta \phi _{1,2}$ will significantly deviate from $\pi$. These corrections can be due to either hard parton radiation, calculated at the matrix element (ME) level at next-to-leading order (NLO), or softer multiple parton radiation described by parton showers. In a recent approach [1], transverse-momentum dependent (TMD) parton densities are obtained with the parton-branching method [2, 3] (PB-TMDs). These PB-TMDs were combined with NLO ME calculations [4] supplemented with PB initial-state parton showers [5], leading to predictions where the initial-state parton shower is determined by the PB-TMD densities without tunable parameters. In Drell-Yan production at the LHC this approach leads to a good description of the measurements [6], whereas other approaches need specific tunes. Therefore, it is interesting to study the PB prediction in a jet environment, especially since the PB-TMD initial-state shower also becomes important. Although $\varDelta \phi _{1,2}$ is an inclusive observable, it is interesting for the theoretical understanding of the complete event to measure the multiplicity of jets in different regions of $\varDelta \phi _{1,2}$ and the transverse momenta of the first four jets. The $\varDelta \phi _{1,2}$ measurement is mainly sensitive to initial-state parton showers at an inclusive level, whereas the measurement of the jet multiplicity in different regions of $\varDelta \phi _{1,2}$ illustrates how many high-$p_{\textrm{T}}$ jets contribute to the $\varDelta \phi _{1,2}$ decorrelation.

The azimuthal correlation in high-$p_{\textrm{T}}$ dijet events was measured previously at: the Fermilab Tevatron in proton-antiproton collisions by the D0 Collaboration at $\sqrt{s}=1.96\,\text {TeV}$ [7, 8]; and at the CERN LHC in $\text {p}\text {p}$ collisions by both the ATLAS Collaboration at $\sqrt{s}=7\,\text {TeV}$ [9] and the CMS Collaboration at $\sqrt{s}=7$, 8, and $13\,\text {TeV}$ [10,11,12,13].

In this paper, we describe new measurements of dijet events with rapidity $|y | < 2.5$ and with transverse momenta of the leading jet $p_{\textrm{T1}} > 200\,\text {GeV}$ and the subleading jet $p_{\textrm{T2}} > 100\,\text {GeV}$. The multiplicity of jets with $p_{\textrm{T}} > 50\,\text {GeV}$ is measured in bins of $p_{\textrm{T1}}$ and $\varDelta \phi _{1,2}$. The jet multiplicity in bins of $\varDelta \phi _{1,2}$ provides information on the $\varDelta \phi _{1,2}$ decorrelation. The cross sections for the four leading jets are measured as a function of $p_{\textrm{T}}$ of each jet, which can give additional information on the structure of the higher-order corrections.

This paper is organized as follows; in Sect. 2, a brief summary of the CMS detector and the relevant components is given. In Sect. 3, the theoretical models for comparison at detector level, as well as with the final results are described. Section 4 gives an overview of the analysis, with the event selection, data correction, and a discussion of the uncertainties. The final results and comparison with theoretical predictions are discussed in Sect. 5.

2 The CMS detector and event reconstruction

The central feature of the CMS apparatus is a superconducting solenoid of 6$\,\text {m}$ internal diameter, providing a magnetic field of 3.8$\,\text {T}$. Within the solenoid volume are silicon pixel and strip tracker detectors, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel part and two endcap sections.

Events of interest are selected using a two-tiered trigger system. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100$\,\text {kHz}$ within a fixed latency of about 4$\,\upmu \text {s}$  [14]. The second level, known as the high-level trigger (HLT), consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, while reduces the event rate to around 1$\,\text {kHz}$ for data storage [15].

During the 2016 data-taking period, a gradual shift in the timing of the inputs of the ECAL first level trigger in the pseudorapidity region $|\eta | > 2.0$, also known as “prefiring”, caused some trigger inefficiencies [14]. For events containing a jet with $p_{\textrm{T}} > 100\,\text {GeV}$, in the region $2.5< |\eta |< 3.0$ the efficiency loss is 10–20%, depending on $p_{\textrm{T}}$, $\eta$, and the data-taking period. Correction factors were computed from data and applied to the acceptance evaluated by simulation.

The particle-flow (PF) algorithm [16] reconstructs and identifies each individual particle in an event, with an optimized combination of information from the various elements of the CMS detector. The energy of charged hadrons is determined from a combination of their momentum measured in the tracker and the matching ECAL and HCAL energy deposits, corrected for the response function of the calorimeters to hadronic showers. The energy of neutral hadrons is obtained from the corresponding corrected ECAL and HCAL energies.

The candidate vertex with the largest value of summed physics-object $p_{\textrm{T}} ^2$ is taken to be the primary vertex (PV) of $\text {p}\text {p}$ interactions as described in Section 9.4.1 of Ref. [17]. The physics objects are the jets, clustered using the jet finding algorithm [18, 19] with the tracks assigned to candidate vertices as inputs, and the associated missing transverse momentum.

Jets are reconstructed from PF objects, clustered using the anti-$k_{\textrm{T}}$ algorithm [18, 19] with a distance parameter of $R=0.4$. The jet momentum is determined as the vectorial sum of all particle momenta in the jet, and is found from simulation to be, typically, within 5–10% of the true momentum over the entire $p_{\textrm{T}}$ spectrum and detector acceptance. Additional $\text {p}\text {p}$ interactions within the same or nearby bunch crossings can contribute additional tracks and calorimetric energy deposits, increasing the apparent jet momentum. To mitigate this effect, tracks identified as originating from pileup vertices are discarded and an offset correction based on the average amount of transverse energy in the event per unit area is applied to correct for the remaining contributions [20]. Jet energy corrections are derived from simulation studies so that the average measured energy of jets becomes identical to that of particle-level jets. However, the selective ECAL readout leads to a bias in the jet energy scale. In situ measurements of the momentum balance in dijet, photon + jet, Z + jet, and multijet events are used to determine any residual differences between the jet energy scale (JES) in data and in simulation, and appropriate corrections are made [21]. Additional selection criteria are applied to each jet to remove jets potentially dominated by instrumental effects or reconstruction failures. The jet energy resolution (JER) amounts typically to 15–20% at 30$\,\text {GeV}$, 10% at 100$\,\text {GeV}$, and 5% at 1$\,\text {TeV}$  [21].

The missing transverse momentum vector ${\vec p}_{\textrm{T}}^{\text {miss}}$ is computed as the negative vector sum of the transverse momenta of all the PF candidates in an event, and its magnitude is denoted as $p_{\textrm{T}} ^\text {miss}$  [22]. The ${\vec p}_{\textrm{T}}^{\text {miss}}$ is modified to include corrections to the energy scale of the reconstructed jets in the event. Anomalous high-$p_{\textrm{T}} ^\text {miss}$ events can be due to a variety of reconstruction failures, detector malfunctions, or noncollision backgrounds. Such events are rejected by event filters that are designed to identify more than 85–90% of the spurious high-$p_{\textrm{T}} ^\text {miss}$ events with a mistagging rate less than 0.1% [22].

A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [23].

3 Theoretical predictions

Theoretical predictions from Monte Carlo (MC) event generators at LO and NLO are used for comparison with measurements of the jet multiplicities as well as with the $p_{\textrm{T}}$ spectra in multijet final states.

We use the following predictions at LO:

• pythia 8 [24] (version 8.219) simulates LO $2\rightarrow 2$ hard processes. The parton shower is generated in a phase space ordered in transverse momentum and longitudinal momentum of the emitted partons, and the colored strings are hadronized using the Lund string fragmentation model. The CUETP8M1 [25] tune (with the parton distribution function (PDF) set NNPDF2.3LO [26]) gives the parameters for multiparton interactions (MPI).

• herwig++ [27] (version 2.7.1) simulates LO $2\rightarrow 2$ hard processes. The emitted partons in the parton shower follow angular ordering conditions, and the cluster fragmentation model is used to transform colored partons into observable hadrons. The CUETHppS1 [25] tune (with the PDF set CTEQ6L1 [28]) is applied.

• MadGraph 5_amc@nlo [4] (version 2.3.3) event generator, labeled MadGraph+Py8, is used in the LO mode, with up to four noncollinear high-$p_{\textrm{T}}$ partons included in the matrix element (ME), supplemented with parton showering and MPI using pythia 8 with the CUETP8M1 tune and merged according to the $k_{\textrm{T}}$-MLM matching procedure [29] with a matching scale of $10\,\text {GeV}$.

• MadGraph 5_amc@nlo [4] (version 2.6.3) event generator, labeled MadGraph+CA3, is used in the LO mode to generate up to four noncollinear high-$p_{\textrm{T}}$ partons included in the ME. The PB-method describes the DGLAP evolution as a process of successive branching processes. It has been shown [2, 3] that on an inclusive level, the DGLAP evolution is exactly reproduced, whereas the simulation of each branching determines the transverse momentum distribution obtained during evolution. The PB-TMDs are obtained from a fit to inclusive HERA deep-inelastic electron proton collision data [1]. Since the PB-TMDs come from a (forward) branching evolution, the initial-state parton shower (in a backward evolution) follows directly from the PB-TMD distributions, and therefore no free parameters are left for the initial-state shower. We use the TMD merging [30] procedure for combining the TMD parton shower with the ME calculations. The NLO PB-TMD set 2 [1] with $\alpha _\textrm{S} (m_{\text {Z}}) = 0.118$ is used. The inclusion of the transverse momentum $k_{\textrm{T}}$ and initial-state PB-TMD parton shower is performed with Cascade3  [5] (version 3.2.3). The initial-state parton shower follows the PB-TMD distribution, and has no free parameters left. The final-state radiation (since not constrained by TMDs) as well as hadronization is performed with pythia 6 (version 6.428) [31] with an angular ordering veto imposed. A merging scale value of 30$\,\text {GeV}$ is used, since it provides a smooth transition between ME and PS computations. MPI effects are not simulated.

At NLO, different theoretical predictions are used. The factorization and renormalization scales are set to half the sum over the scalar transverse momenta of all produced partons, $1/2 \sum H_{\textrm{T}}$. However, we have explicitly checked that the distributions do not change when choosing $1/4 \sum H_{\textrm{T}}$ or even $1/6 \sum H_{\textrm{T}}$. The uncertainty bands of the NLO predictions are determined from the variation of the factorization and renormalization scales by a factor of two up and down, avoiding the largest scale differences.

• MadGraph 5_amc@nlo [4] (version 2.6.3) interfaced with pythia 8 (version 8.306), labeled MG5_aMC+Py8 (jj), is used with MEs computed at NLO for the process $\text {p}\text {p}\rightarrow \textrm{jj}$. The NNPDF 3.0 NLO PDF set [32] is used and $\alpha _\textrm{S} (m_{\text {Z}})$ is set to 0.118. This calculation uses the CUETP8M1 tune to study the effect of MPI.

• MadGraph 5_amc@nlo [4] (version 2.6.3) interfaced with Cascade3 (version 3.2.2) [5], labeled MG5_aMC+CA3 (jj), is used with MEs computed at NLO for the process $\text {p}\text {p}\rightarrow \textrm{jj}$. The Herwig6 [33] subtraction terms are used in mc@nlo  [34, 35], since they are closest to the needs for applying PB-TMD parton densities, as described in Ref. [5]. The NLO PB-TMD set 2 [1] with $\alpha _\textrm{S} (m_{\text {Z}}) = 0.118$ is used. MPI are not simulated in this approach.

• MadGraph 5_amc@nlo [4] (version 2.6.3) interfaced with Cascade3, labeled MG5_aMC+CA3 (jjj) NLO, is used with MEs computed at NLO for the process $\text {p}\text {p}\rightarrow \textrm{jjj}$. The same PB-TMD distribution and parton shower calculation procedure as for MG5_aMC+CA3 (jj) are used.

In Table 1, the Monte Carlo event generators utilized in this analysis are summarized. Events generated by pythia 8, MadGraph+Py8, and herwig++ are passed through a full detector simulation based on Geant4  [36]. The simulated events are reconstructed the same way as the observed data events.

Table 1 Description of the simulated samples used in the analysis

Table 2 The integrated luminosity for each trigger sample considered in this analysis with the $p_{\textrm{T}}$ thresholds for HLT (PF) reconstruction

4 Data analysis

The measured data samples recorded in 2016, corresponding to an integrated luminosity of $36.3{\,\text {fb}^{-1}}$, were collected with single-jet HLTs. For each single-jet HLT at least one jet with $p_{\textrm{T}}$ higher than the $p_{\textrm{T}}^{\textrm{HLT}}$ trigger threshold is required. All triggers, except the one with the highest $p_{\textrm{T}}^{\textrm{HLT}}$ threshold, were prescaled. We consider events only if the leading jet, reconstructed with the PF algorithm, can be matched with an HLT jet. In Table 2, the integrated luminosity $\mathcal {L}$ for each trigger is shown. The trigger efficiency ($>99.5\%$) is estimated using triggers with lower $p_{\textrm{T}}^{\textrm{HLT}}$ thresholds; for the trigger with lowest $p_{\textrm{T}}^{\textrm{HLT}}$ threshold a tag-and-probe [14] method is used to determine the $p_{\textrm{T}}$ threshold. To address the issue of trigger inefficiency resulting from prefiring (discussed in Sect. 2), the simulated event samples are corrected on an event-by-event basis prior to unfolding. Maps of the prefiring probability in the region of $2.0< |\eta | < 3.0$ are utilized, taking into account the $p_{\textrm{T}}$ and $\eta$ of the jets. The total event weight is then calculated as the product of the nonprefiring probability of all jets. For each $p_{\textrm{T}}$ and $\eta$ bin in the prefiring maps, the uncertainty per jet is determined by taking the maximum value between 20% of the prefiring probability and the statistical uncertainty.

The jets are corrected using the JES correction procedure in CMS [21], and an additional smoothing procedure (described in Ref. [37]) is applied to the JES correction. The simulated samples are corrected to take into account the JER by spreading the $p_{\textrm{T}}$ of the jets according to the resolution extracted from data. Jets reconstructed in regions of the detector corresponding to noisy towers in the calorimeters are excluded from the measurement [38].

The pileup profile used in simulation is corrected to reproduce the one in data.

4.1 Event selection

Each event is required to have a reconstructed PV. The PV must satisfy $|z_\textrm{PV} |<24\,\text {cm}$ and $\rho _\textrm{PV}<2\,\text {cm}$, where $z_\textrm{PV}$ ($\rho _\textrm{PV}$) is the longitudinal (radial) distance from the nominal interaction point.

All events that contain jets clustered using the anti-$k_{\textrm{T}}$ algorithm [18, 19] with a distance parameter of $R=0.4$ and reconstructed with $|\eta |<3.2$ and transverse momentum $p_{\textrm{T}} >20\,\text {GeV}$ are preselected. From these events, the ones with at least a pair of jets with $p_{\textrm{T1}} >200\,\text {GeV}$, $p_{\textrm{T2}} >100\,\text {GeV}$ and $|y^{{1,2}} |<2.5$ are selected (events with one of the leading two jets with $|y |>2.5$ are vetoed). Additional jets must have $p_{\textrm{T}} > 50\,\text {GeV}$ and $|y | < 2.5$. In addition, jets must satisfy quality criteria based on the jet constituents, in order to reject misidentified jets [39]. The selected events must have a missing transverse energy fraction smaller than 0.1 (more details in Sect. 4.3). The selection at particle-level includes only the $p_{\textrm{T}}$ and $|y |$ constraints on the jets.

4.2 Observables and phase space

In the following, the observables will be described at particle-level. The particle-level is defined after all the partons have hadronized to form stable particles with a proper lifetime above 10$\,\text {mm}$/c.

The exclusive jet multiplicity ($N_{\text {jets}}$) of up to 6 jets (inclusive for $N_{\text {jets}} \ge 7$), in three bins of $p_{\textrm{T1}}$ ($200< p_{\textrm{T1}} < 400$, $400< p_{\textrm{T1}} < 800$, and $p_{\textrm{T1}} > 800\,\text {GeV}$) and for three different regions in $\varDelta \phi _{1,2}$ ($0< \varDelta \phi _{1,2} < 150$, $150< \varDelta \phi _{1,2} < 170$, $170< \varDelta \phi _{1,2} < 180^\circ$) is measured:

$$\begin{aligned} \frac{\text {d}\sigma _\text {dijet}}{\text {d}N^\textrm{i}_\text {jets}\text {d}p^\textrm{j}_\textrm{T1}\text {d}(\varDelta \phi _{1,2} ^\textrm{k})} \end{aligned}$$

(1)

where $\textrm{i},\textrm{j},\textrm{k}$ corresponds to the binning in $N_{\text {jets}}$, $p_{\textrm{T1}}$, and $\varDelta \phi _{1,2}$. We measure exclusive jet multiplicities since we are interested in how the hard and soft multijet radiation is modeled by including higher-order corrections and parton shower effects.

In addition, the differential cross section as a function of the $p_{\textrm{T}}$ of each of the four leading jets are measured:

$$\begin{aligned} \frac{\text {d} \sigma _{\text {p}\text {p}\rightarrow \textrm{jj}}}{\text {d} p_{\textrm{T1}}}, \quad \frac{\text {d} \sigma _{\text {p}\text {p}\rightarrow \textrm{jj}}}{\text {d} p_{\textrm{T2}}}, \quad \frac{\text {d} \sigma _{\text {p}\text {p}\rightarrow \textrm{jjj}}}{\text {d} p_{\textrm{T3}}}, \quad \frac{\text {d} \sigma _{\text {p}\text {p}\rightarrow \textrm{jjjj}}}{\text {d} p_{\textrm{T4}}}, \end{aligned}$$

(2)

where $\sigma _{\text {p}\text {p}\rightarrow \textrm{jj}}$, $\sigma _{\text {p}\text {p}\rightarrow \textrm{jjj}}$, $\sigma _{\text {p}\text {p}\rightarrow \textrm{jjjj}}$ correspond to the dijet, three-jet, and four-jet cross sections in proton–proton collision.

Fig. 1

Distribution of $E_{\textrm{T}}^{\text {miss}}/\sum {E_{\textrm{T}}}$ for data and simulated jet production for three regions of $\varDelta \phi _{1,2}$. Shown are the contributions from QCD, $\text {W}/\text {Z}$ and ${\text {t}{}\overline{\text {t}}}$ events. The main contributions of events with large $E_{\textrm{T}}^{\text {miss}}$ in the final state come from processes like $\text {Z}\rightarrow \upnu \overline{\upnu }$ and $\text {W}\rightarrow \text {l}\upnu$. The data (MC prediction) statistical uncertainty is shown as a vertical line (grey shaded bar in the ratio)

4.3 Background treatment

To remove background contributions from ${\text {t}{}\overline{\text {t}}} +\text {jets}$, $\text {W}/\text {Z}+\text {jets}$, $\text {Z}\rightarrow \upnu \overline{\upnu }$ and $\text {W}\rightarrow \text {l}\upnu$, events with $E_{\textrm{T}}^{\text {miss}}/\sum {E_{\textrm{T}}} > 0.1$ are rejected. The missing transverse energy is calculated from $p_{\textrm{T}} ^\text {miss}$ and the sum $\sum {E_{\textrm{T}}}$ runs over all PF objects in the event. Less than 1% of the events are rejected in the whole sample. For $N_{\text {jets}} = 2$ this constraint becomes important and the backgournd contribution is reduced from 20% to the percent level for $p_{\textrm{T1}} > 800\,\text {GeV}$ in the low $\varDelta \phi _{1,2}$. In Fig. 1, the measured fraction $E_{\textrm{T}}^{\text {miss}}/\sum {E_{\textrm{T}}}$ is compared with the simulation of signal and background processes in bins of $\varDelta \phi _{1,2}$.

Fig. 2

Probability matrix (condition number: 3.0) for the jet multiplicity distribution constructed with the MadGraph+Py8 sample. The global 3$\times$3 sectors (separated by the thick black lines) correspond to the $p_{\textrm{T1}}$ bins, indicated by the labels in the x (lower) and y (left) axes; the smaller 3$\times$3 structures correspond to the $\varDelta \phi _{1,2}$ bins, indicated in the leftmost row and lowest column, the x(y) axis of these $\varDelta \phi _{1,2}$ cells corresponds to the jet multiplicity at particle (detector) level. The z axis covers a range from $10^{-6}$ to 1 indicating the probability of migrations from the particle-level bin to the corresponding detector-level bin. The lower right corner of the matrix, which describes (extreme) migrations between hight $p_{\textrm{T1}}$ at particle level and low $p_{\textrm{T1}}$ at dectector level, is characterized by statistical fluctuations

4.4 Correction for detector effects

The measured cross sections are corrected for misidentified events, detector resolution, and inefficiencies for comparison with particle-level predictions through the procedure of unfolding. A response matrix (RM) mapping the “true” distribution onto the measured one is constructed using simulated MC samples from pythia 8, MadGraph+Py8, and herwig++. The MadGraph+Py8 sample, which has the smallest statistical uncertainty, is used as default for constructing the RM, whereas the herwig++ and pythia 8 samples are used to investigate the model dependence. The RM is constructed by matching the detector and particle-level distributions. If events (or jets) cannot be matched, they contribute to the background or inefficiency distributions. For the jet multiplicity, the dijet system (leading and subleading jets) is matched if the jets coincide within $\varDelta R < 0.2$ (half of the jet radius of $R = 0.4$). For the $p_{\textrm{T}}$ distributions, jets are matched with $\varDelta R < 0.2$, and from the matched candidates the one highest in $p_{\textrm{T}}$ is selected (only events with at least two jets are considered in the matching). The TUnfold (version 17.9) package [40] is used to perform the unfolding.

The jet multiplicity is obtained by matrix inversion:

$$\begin{aligned} \textbf{A}\alpha + \beta = \gamma \end{aligned}$$

(3)

where $\gamma (\alpha )$ is the distribution at detector (particle) level, $\textbf{A}$ is the probability matrix (PM), and $\beta$ is the contribution from background (or noise). The PM is obtained by normalizing the RM to the particle-level axis (row-by-row normalization), and describes the probability that a particle-level jet (or event) generated in a bin is reconstructed in (migrates to) another bin at detector-level. The condition number of the PM is defined as the ratio between the largest and smallest singular value of the matrix. The condition number of the PM for the jet multiplicity is 3.0, and matrix inversion is possible.

For the $p_{\textrm{T}}$ distributions, the condition number of the PM is 4.9, and least-square minimisation is applied:

$$\begin{aligned} \chi ^2&= \min _{\alpha } \left[ (\gamma - \beta - \textbf{A} \alpha )^\intercal \, {\textbf{V}_{\gamma \gamma }}^{-1}\, (\gamma - \beta - \textbf{A}\alpha ) \right] , \end{aligned}$$

(4)

where $\textbf{V}_{\gamma \gamma }$ represents the statistical covariance matrix of $\gamma$ and $\beta$, with twice the number of bins at detector level compared to particle-level. No additional regularization is needed.

Fig. 3

Probability matrix (condition number: 4.9) for the $p_{\textrm{T}}$ of the four leading jets constructed with the MadGraph+Py8 sample. The global 4$\times$4 sectors correspond to the $p_{\textrm{T}}$ distributions each of the first four jets, the x axis corresponds to the particle (gen) level and y axis corresponds to the detector (rec) level. The z axis covers a range from $10^{-6}$ to 1 indicating the probability of migrations from the particle-level bin to the corresponding detector-level bin

In Figs. 2 and 3, the PMs are shown for the multiplicity distributions and the $p_{\textrm{T}}$ distributions of the first four jets.

4.5 Uncertainties

The statistical uncertainties in the measured spectra and response matrices are propagated to the final results. In Figs. 4 and 5, the statistical correlations are shown for the jet multiplicity and for the $p_{\textrm{T}}$ spectra of the four leading jets.

Fig. 4

Correlation matrix at the particle-level for the jet multiplicity distribution. It contains contributions from the data and from the limited-size MadGraph+Py8 sample. The global 3$\times$3 sectors (separated by the thick black lines) correspond to the $p_{\textrm{T1}}$ bins, indicated by the labels next to the x (lower) and y (left) axes; the smaller 3$\times$3 structures correspond to the $\varDelta \phi _{1,2}$ bins, indicated in the leftmost row and lowest column, the x and y axes correspond to the jet multiplicity. The z axis covers a range from $-\,1$ to 1 indicating the correlations in blue shades and anticorrelations in red shades, the values between $-\,0.1$ and 0.1 are represented in white

Fig. 5

Correlation matrix for the particle-level $p_{\textrm{T}}$ of the four leading jets. It contains contributions from the data and from the limited-size MadGraph+Py8 sample. Here each one of the 4$\times$4 sectors corresponds to one of the $p_{\textrm{T}}$ spectra measured, indicated by the x and y axis labels. The z axis covers a range from $-\,1$ to 1 indicating the correlations in blue shades and anticorrelations in red shades, the values between $-\,0.1$ and 0.1 are represented in white

The systematic uncertainties originate from the following sources:

• JES: The JES uncertainty is estimated from variations of one standard deviation in the JES corrections applied to data (at detector-level) and repeating the whole unfolding procedure for each variation.

• JER: The JER uncertainty is estimated by varying the resolution by one standard deviation in the simulated sample, and repeating the unfolding for each variation.

• Integrated luminosity: The uncertainty in the integrated luminosity is 1.2% [41] and is applied as a global scaling factor to the cross section.

• Pileup: The pileup distribution in the simulated samples is reweighted to match that of the data. The corresponding uncertainty is estimated by varying the total inelastic cross section by $\pm \,5\%$ [42], affecting the measurement by less than 1%.

• Trigger prefiring uncertainty: The trigger prefiring uncertainty is estimated by varying the correction applied to the simulated samples and repeating the unfolding for each variation, resulting in an uncertainty of 1–3%.

• Model uncertainty: The model uncertainty is estimated by varying the distributions of the factorization and renormalization scales in the MC sample such as to still describe the detector level distributions. Additionally background and inefficiencies are varied separately by 15%. The final uncertainty is the quadratic sum of each of the uncertainties. It was validated with a set of closure tests performed using pythia 8 and herwig++ samples as pseudodata unfolded with MadGraph+Py8 RM. This uncertainty ranges from 1 to 7%.

The total systematic uncertainty is obtained by adding all the systematic uncertainties in quadrature, assuming independent sources.

In Fig. 6, the relative uncertainties for the jet multiplicity in bins of $p_{\textrm{T1}}$ and $\varDelta \phi _{1,2}$ are shown. The dominant uncertainty is JES. The total statistical uncertainty (stat. unc.) is mainly driven by the limited event counts in data. The total experimental uncertainty (Total) is typically about 10–15%.

In Fig. 7, the relative uncertainties as a function of the jet $p_{\textrm{T}}$ for the four leading jets are shown. The dominant uncertainty is JES. The measurement is limited by the systematic uncertainty and the total experimental uncertainty ranges from 5 to 10%.

Fig. 6

Relative uncertainties for JES, JER, “Other” and total statistical uncertainty for the jet multiplicity distribution in bins of $p_{\textrm{T1}}$ and $\varDelta \phi _{1,2}$. Here “Other” indicates luminosity, pileup, prefiring, and model uncertainty

Fig. 7

Relative uncertainties for JES, JER, “Other” and total statistical uncertainty for the $p_{\textrm{T}}$ distributions of the four leading jets. Here “Other” indicates luminosity, pileup, prefiring, and model uncertainty

5 Results

The phase space of the measurements at particle-level (particles with a proper lifetime above 10$\,\text {mm}$/c) is defined by jets clustered using the anti-$k_{\textrm{T}}$ algorithm [18, 19] with a distance parameter of $R=0.4$ within $|y | < 3.2$. Events are selected if the two highest $p_{\textrm{T}}$ jets with $p_{\textrm{T1}} > 200\,\text {GeV}$, $p_{\textrm{T2}} > 100\,\text {GeV}$ have $|y |< 2.5$ (i.e., events are not counted where one of the leading jets has $|y |> 2.5$). For the additional jets $p_{\textrm{T}} >50\,\text {GeV}$ and $|y |< 2.5$ is required. All predictions (LO and NLO) are normalized to the measured dijet cross section in the measurement phase space, with the normalization factors explicitly given in the figures.

5.1 Jet multiplicity distribution

The multiplicity of jets with $p_{\textrm{T}} > 50\,\text {GeV}$ in $|y |< 2.5$ is measured for various regions of the transverse momentum of the leading jet, $p_{\textrm{T1}}$, and the azimuthal angle $\varDelta \phi _{1,2}$ between the two leading jets as shown in Fig. 8.

As a characterization of the jet multiplicity we compare the production rate for 3 jets with the one for 7 jets. In the region of low $p_{\textrm{T1}}$ ($200< p_{\textrm{T1}} < 400\,\text {GeV}$), a large number of additional jets is observed at low $\varDelta \phi _{1,2}$ ($0< \varDelta \phi _{1,2} < 150^\circ$), the production rate between 3 and 7 jets changes by two orders of magnitude. In the large-$\varDelta \phi _{1,2}$ region ($170< \varDelta \phi _{1,2} < 180^\circ$), where the leading jets are nearly back-to-back, the production rate for 3 to 7 jets changes by three orders of magnitude. We note that even in this back-to-back region a large number of additional jets are observed that do not contribute significantly to the momentum imbalance of the two leading jets.

In the region of large $p_{\textrm{T1}}$ ($p_{\textrm{T1}} >800\,\text {GeV}$), we observe that the rate of additional jets at low $\varDelta \phi _{1,2}$ is essentially constant and the rate between 3 and 7 jets only shows small changes, indicating that many jets participate in the compensation of the $\varDelta \phi _{1,2}$ decorrelation. In the large-$\varDelta \phi _{1,2}$ region ($170< \varDelta \phi _{1,2} < 180^\circ$), the rate between 3 and 7 jets changes by less than 2 orders of magnitude, in contrast to the low-$p_{\textrm{T1}}$ region. A multiplicity of more than two or three additional jets at large $p_{\textrm{T1}}$ is observed in all regions of $\varDelta \phi _{1,2}$.

In Fig. 8, predictions from the LO $2\rightarrow 2$ generators pythia 8 and herwig++ including parton showering and MPI are shown. The shape of the prediction coming from pythia 8 is different from what is observed in the measurement. The shape of the prediction from herwig++ agrees rather well with the measurement, especially in the large $\varDelta \phi _{1,2}$ region. The difference between pythia 8 and herwig++ in jet multiplicity comes from the different treatment of the parton shower. In addition the predictions from MadGraph+Py8 and MadGraph+CA3 with additional noncollinear high-$p_{\textrm{T}}$ partons, supplemented with different parton showering approaches and MPI are shown. We find that the merged prediction from MadGraph+Py8 does not agree in shape with the measurement, whereas the MadGraph+CA3 merged prediction (for $N_{\text {jets}} > 2$ and $p_{\textrm{T1}} < 800 \,\text {GeV}$) does, similarly to the case of herwig++.

Fig. 8

Differential cross section as a function of the exclusive jet multiplicity (inclusive for 7 jets) in bins of $p_{\textrm{T1}}$ and $\varDelta \phi _{1,2}$. The data are compared with LO predictions of pythia 8, herwig++, MadGraph+Py8 and MadGraph+CA3. The predictions are normalized to the measured dijet cross section using the scaling factors shown in the legend. The vertical error bars correspond to the statistical uncertainty, the yellow band shows the total experimental uncertainty

The calculations with NLO MEs matched with parton shower compared with the measurements are shown in Fig. 9. The uncertainty bands of the predictions come from the variation of the factorization and renormalization scales by a factor of two up and down, avoiding the largest scale differences. The normalization of the MG5_aMC+Py8 (jj) calculation is in reasonable agreement with the measured cross section even for three jets. For higher jet multiplicities the prediction falls below the measurement. MG5_aMC+CA3 (jj) predicts a smaller cross section for more than three jets compared with the measurement. The MG5_aMC+CA3 (jjj) NLO calculation (using the same normalization factor as for MG5_aMC+CA3 (jj)) predicts a larger three- and four-jet cross section, whereas the higher jet multiplicities are still underestimated.

Fig. 9

Differential cross section as a function of the exclusive jet multiplicity (inclusive for 7 jets) in bins of $p_{\textrm{T1}}$ and $\varDelta \phi _{1,2}$. The data are compared with NLO dijet predictions MG5_aMC+Py8 (jj) and MG5_aMC+CA3 (jj) as well as the NLO three-jet prediction of MG5_aMC+CA3 (jjj). The vertical error bars correspond to the statistical uncertainty, the yellow band shows the total experimental uncertainty. The shaded bands show the uncertainty from a variation of the renormalization and factorization scales. The predictions are normalized to the measured inclusive dijet cross section using the scaling factors shown in the legend

5.2 Transverse momenta of the four leading jets

The measured differential jet cross section as a function of the jet transverse momentum, $p_{\textrm{T}}$, for the four leading $p_{\textrm{T}}$ jets is shown in Fig. 10 and compared with the predictions of pythia 8 and MG5_aMC+Py8 (jj). The $p_{\textrm{T}}$ values of the two leading-$p_{\textrm{T}}$ jets reaches up to 2$\,\text {TeV}$ and the third and fourth jets $p_{\textrm{T}}$ reaches about 1$\,\text {TeV}$. We observe that the shape of the $p_{\textrm{T}}$ spectrum for the third and fourth leading jets is qualitatively similar to the one of the two leading jets, whereas the cross section is different. The rise of the cross section for the second jet between 100$\,\text {GeV}$ and 200$\,\text {GeV}$ is a consequence of the higher $p_{\textrm{T}}$ requirement ($p_{\textrm{T1}} > 200\,\text {GeV}$) applied to the leading jet in the event selection.

Fig. 10

Transverse momenta of the four leading jets, with the yellow band representing the total experimental uncertainty. The data are compared with LO (pythia 8) and NLO (MG5_aMC+Py8) predictions. The red band in the NLO prediction represents the renormalization and factorization scale uncertainty

In Fig. 11, the measured differential cross section as a function of the $p_{\textrm{T}}$ for the four leading jets is shown and compared with LO predictions (using the same normalization factors as in Fig. 8). Only the prediction of pythia 8 is able to describe reasonably well the measurement in shape, except for the region $p_{\textrm{T2}} < 200\,\text {GeV}$. Among the LO calculations, pythia 8 and herwig++ predict a rather flat distribution for the third and fourth jet, whereas the other predictions are different in shape. The predictions from herwig++ are not in agreement in shape and rate with the measurements, the differences are up to 50% for the leading and subleading jets at large $p_{\textrm{T}}$. The prediction from MadGraph+Py8 gives a significantly different shape for the $p_{\textrm{T}}$ spectrum for the first 3 jets. Comparing the merged predicitons, MadGraph+CA3 gives a significant improvement in the description of the shape of the $p_{\textrm{T}}$ distribution for the three leading jets, whereas the description of the distribution of the fourth jet is similar to the one with MadGraph+Py8.

The predictions obtained with NLO MEs are shown in Fig. 12 using the same normalization factors as in Fig. 9. The uncertainty bands of the predictions come from the variation of the factorization and renormalization scales by a factor of two up and down, avoiding the largest scale differences. The predictions of MG5_aMC+Py8 (jj) and MG5_aMC+CA3 (jj) describe the normalization and the shape of the first three jets rather well, whereas for the fourth jet (which comes from the parton shower) falls below the measurement by 20–30%. The prediction of MG5_aMC+CA3 (jjj) describes the third and fourth jets rather well within uncertainties (predictions for the first and second jet are meaningless for MG5_aMC+CA3 (jjj) and therefore not shown). The calculations using PB-TMDs together with NLO MEs in the MC@NLO frame provide a good description of jet measurements over a wide range in transverse momentum and jet multiplicities without tunable parameters in the initial-state parton shower.

Fig. 11

Transverse momentum distributions of the four leading jets. The transverse momentum of the leading and subleading (third and fourth leading) $p_{\textrm{T}}$ jets from left to right is shown in the upper (lower) figure. The data are compared with LO predictions of pythia 8, herwig++, MadGraph+Py8 and MadGraph+CA3. The predictions are normalized to the measured dijet cross section using the scaling factors shown in the legend. The vertical error bars correspond to the statistical uncertainty, the yellow band shows the total experimental uncertainty

Fig. 12

Transverse momentum distributions of the four leading jets. The transverse momentum of the leading and subleading (third and fourth leading) $p_{\textrm{T}}$ jets from left to right is shown in the upper (lower) figure. The data are compared with NLO predictions MG5_aMC+Py8 (jj) and MG5_aMC+CA3 (jj) as well as the NLO three-jet prediction of MG5_aMC+CA3 (jjj). The vertical error bars correspond to the statistical uncertainty, and the yellow band to total uncertainty of the measurement. The bands show the uncertainty from a variation of the renormalization and factorization scales. The predictions are normalized to the measured dijet cross section using the scaling factors shown in the legend

6 Summary

A study of multijet events has been performed in proton–proton collisions at a center-of-mass energy of 13$\,\text {TeV}$ using data collected with the CMS detector corresponding to an integrated luminosity of $36.3{\,\text {fb}^{-1}}$. The measurements are performed by selecting a dijet system containing a jet with $p_{\textrm{T}} > 200\,\text {GeV}$ and a subleading jet with $p_{\textrm{T}} >100\,\text {GeV}$ within $|y | < 2.5$.

For the first time, the jet multiplicity in bins of the leading jet $p_{\textrm{T}}$ and the azimuthal angle difference between the two leading jets, $\varDelta \phi _{1,2}$, is measured. The jet multiplicity distributions show that even in the back-to-back region of the dijet system, up to seven jets are measurable. The differential production cross sections are measured for the highest $p_{\textrm{T}}$ jets up to the $\,\text {TeV}$ scale.

The measurement of the differential cross section as a function of the jet $p_{\textrm{T}}$ for the four highest $p_{\textrm{T}}$ jets is an important reference for QCD multijet cross section calculations, and especially for the simulations including parton showers for higher jet multiplicity.

The measured multiplicity distribution of jets with $p_{\textrm{T}} > 50\,\text {GeV}$ and $|y | < 2.5$ is not well described by the leading order MadGraph +pythia 8 simulation. However, in the back-to-back region herwig++ and MadGraph + Cascade3 provide a better description of the shape of the jet multiplicity. The measured differential cross section as a function of the transverse momentum of the four leading $p_{\textrm{T}}$ jets is not described by any of the LO predictions either in normalization or in shape. However, MadGraph + Cascade3 describes the shape of the distribution better than MadGraph +pythia 8.

The predictions using dijet NLO matrix elements, MG5_aMC+pythia 8 (jj) and MG5_aMC+ Cascade3 (jj) describe the lower multiplicity regions, as well as the transverse momenta of the leading jets, reasonably well. The three-jet NLO calculation MG5_aMC+ Cascade3 (jjj) describes very well the cross section of the third and fourth jets.

The measurements presented here provide stringent tests of theoretical predictions in the perturbative high-$p_{\textrm{T}}$ and high-jet multiplicity regions. Although the higher jet multiplicities are not described with either parton shower approach, it is interesting that the lower jet multiplicity cross section is described satisfactorily with NLO dijet calculations supplemented with PB-TMDs and TMD parton shower with fewer tunable parameters than in the case with conventional parton showers.

The measured observables and its statistical correlations are provided in HEPData [43] as tabulated results.