Measurement of double-parton scattering in inclusive production of four jets with low transverse momentum in proton-proton collisions at $\sqrt{s}$ = 13 TeV

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 $\lvert\eta\rvert$ $\lt$ 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.


Introduction
Quantum chromodynamics (QCD), the theory of strong interactions, provides a good description of the production of hadron jets with large transverse momentum (p T ) in high-energy proton-proton (pp) collisions. This is achieved by factorizing the cross section into a perturbatively calculable matrix element describing the scattering between partons, and parton distribution functions (PDFs) that provide the probability to find a parton with given properties within the proton. The PDFs cannot be perturbatively calculated and are obtained by fitting available data. This fitting process includes nonperturbative effects such as the underlying event, hadronization, and parton showering. Measurements of the cross section for the production of inclusive high-p T jets have been performed by the CMS collaboration at various center-of-mass energies and show agreement with perturbative QCD predictions at next-toleading-order (NLO) accuracy [1-3]. However, final states with multiple jets are not as well understood [4], suggesting a need for additional theoretical treatments of the strong interaction.
Multijet final states can be produced in a single-parton scattering (SPS). Depending on the order of the matrix element in the strong coupling, two or more jets can be produced in SPS. Radiation before and/or after the interaction between the partons, as described by parton shower models, can contribute additional jets to the final state. Thus, predictions for multijet processes in SPS provide an important test of the matching between fixed-order matrix element calculations and the parton-shower formalism. A different approach introduces an additional hard scattering in the pp collision, which also contributes a number of jets to the final state. Such processes are in general referred to as double-parton scattering (DPS), and they represent the simplest case of multiple-parton interactions (MPI). A schematic depiction of inclusive four-jet production through SPS and DPS is shown in Fig. 1.
The cross section of a DPS process, σ DPS A,B , where A and B denote two processes with their own respective cross sections σ A and σ B , can be expressed as: The factor m is a combinatorial factor, which is equal to 1 for identical processes and 2 for nonidentical processes. The effective cross section (σ eff ) reflects how strongly the occurrence of A and B is correlated [5]. For fully uncorrelated production of A and B, σ eff tends to the total inelastic pp cross section, whereas a small σ eff indicates an enhanced simultaneous occurrence of processes A and B. For multijet production, SPS processes often exhibit strong kinematic correlations between all jets, whereas DPS processes will manifest a distinctly different behavior. Indeed, the jets resulting from DPS are more often produced in two independent pairs, each in a back-to-back configuration in the transverse plane. The relevance of DPS rises with increasing center-of-mass energy; at higher energy and for fixed p T of the jets, smaller values of the momentum fraction of the protons carried by the partons are probed, resulting in a strong increase of the gluon density and a larger probability for DPS. A study of the extent to which DPS processes can supplement various SPS models is therefore beneficial for a complete description of hadronic interactions.
Various This paper presents an analysis of the inclusive production of four-jet events in pp collisions at a center-of-mass energy of 13 TeV. The data correspond to an integrated luminosity of 42 nb −1 and were collected with the CMS detector at the CERN LHC in 2016 during a special datataking period with a low probability for several pp interactions occurring within the same or neighbouring bunch crossings (pileup). This avoids the challenges posed by pileup and enables us to include jets with low p T . As a result of the low-p T jets, a custom calibration of the jet energy scale is required. Data are corrected for detector efficiency and resolution effects by means of an unfolding procedure.
Several aspects of multijet production are studied by comparing the distributions of DPSsensitive observables predicted by various Monte Carlo (MC) event generators with the distributions measured in data. These observables all exploit differences in the kinematic correlations between the jets expected for SPS and DPS. The DPS cross section is extracted with a template method. A pure DPS signal template is reconstructed from data by randomly mixing two inclusive single-jet events into one inclusive DPS four-jet event. This is then fitted together with several SPS-only background MC models to the distributions obtained from inclusive four-jet data. Finally, the effective cross section is computed using Eq. (1), with σ A and σ B measured from data.
Tabulated results are provided in the HEPData record for this analysis [27].
This paper is organized as follows. In Section 2, the observables of interest are defined. Section 3 gives a brief overview of the CMS detector. The MC models used in the comparison with data are detailed in Section 4, whereas the data samples, event selection, and correction procedure are discussed in Section 5. The strategy for the extraction of the DPS cross section and σ eff is detailed in Section 6. Systematic uncertainties for each of the unfolded observables are discussed in Section 7. Section 8 contains a discussion of the results, which are summarized in Section 9.
• The azimuthal angular difference between the two softest jets: The two softest jets are more likely to be in a back-to-back configuration when produced by DPS since there is an increased probability that the two softest jets are produced in an independent scattering fron the two hardest jets and since the momentum should be conserved in the transverse plane. The increased probability leads to an enhanced DPS contribution around ∆φ Soft = π. • The minimal combined azimuthal angular range of three jets: In DPS, at least two out of three jets are more likely to be in a back-to-back configuration, while SPS processes have a more random distribution in their azimuthal angular difference. Therefore, a DPS process is prone to yield larger values of ∆φ min 3j [33]. • The maximum η difference between two jets: As the maximum separation in η between two jets becomes larger, the probability for the two jets to originate from two different parton interactions increases.
• The azimuthal angular difference between the jets with the largest η separation: Since the jets with the largest η separation are more likely to be produced in separate DPS subprocesses, a decorrelation in the distribution of the azimuthal angular difference of these jets is expected, whereas the jets will show stronger correlations in a SPS event.
• The transverse momentum balance of the two softest jets: When the two softest jets originate from a DPS process, they are more likely to be in a back-to-back configuration rendering the value for ∆p T,Soft small. In SPS processes, the two softest jets do not necessarily balance.
• The azimuthal angular difference between the hard and the soft jet pairs: In a SPS process, the four jets must balance so that the ∆S distribution peaks around π, while in DPS the two jet pairs are more likely to be independently produced, yielding a less correlated ∆S distribution. Thus we anticipate that DPS events tend toward lower values of ∆S.

Measuring jets with the CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections. Forward calorimeters extend the η coverage provided by the barrel and endcap detectors. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid.
The silicon tracker measures charged particles within the range |η| < 2.5. During the LHC running period when the data used in this paper were recorded, the silicon tracker consisted of 1440 silicon pixel and 15 148 silicon strip detector modules.
The ECAL consists of 75 848 lead tungstate crystals, which provide coverage in |η| < 1.48 in a barrel region and 1.48 < |η| < 3.0 in the two endcap regions. Preshower detectors consisting of two planes of silicon sensors interleaved with a total of 3X 0 of lead are located in front of each endcap detector.
In the η-φ plane, and for |η| < 1.48, the HCAL cells map on to 5 × 5 arrays of ECAL crystals to form calorimeter towers projecting radially outwards from close to the nominal interaction point. For |η| > 1. 74, the coverage of the towers increases progressively to a maximum of 0.174 in ∆η and ∆φ. Within each tower, the energy deposits in ECAL and HCAL cells are summed to define the calorimeter tower energies.
The forward hadron (HF) calorimeter uses steel as an absorber and quartz fibers as the sensitive material. The two halves of the HF are located 11.2 m from the interaction region, one on each end, and together they provide coverage in the range 3.0 < |η| < 5.2. They also serve as luminosity monitors.
Events of interest are selected using a two-tiered trigger system. The first level (L1), composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz within a fixed latency of about 4 µs [34]. 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, and reduces the event rate to around 1 kHz before data storage [35].
A global event reconstruction (particle flow) algorithm [36] reconstructs and identifies each individual particle in an event, with an optimized combination of all subdetector information. Jets are clustered from these reconstructed particles using the infrared and collinear safe antik T algorithm [37,38] with a distance parameter of 0.4. The jet momentum is determined as the vectorial sum of the momenta of all particles in the jet, and is typically within 5 to 10% of the true momentum over the whole p T spectrum and detector acceptance, based on simulation.
The CH3 tune has been obtained by the CMS Collaboration from a study of underlying event measurements [58]. It is used in combination with the NNPDF2.3 NNLO PDF.

Multijet models
A second group of models, referred to as the multijet models, uses higher-order matrix elements to produce more than two jets in the hard parton scattering. These MC event generators are interfaced with PYTHIA, HERWIG, or CASCADE [59] to include a description of the underlying event. Details of the generated event samples are given below.
• MADGRAPH5 aMC@NLO (version 2.6.5) [60] is a generator with the ability to compute tree-level and NLO matrix elements for arbitrary processes. Two LO samples and one NLO sample are generated, as listed below.
The LO samples combine a 2 → 2, a 2 → 3, and a 2 → 4 matrix element, referred to as 2 → 2, 3, 4. An H T > 50 GeV generation condition is used, where H T is defined as the scalar sum of the transverse momenta of the produced partons, and all partons must have |η| < 5. For one sample the description of the underlying event, parton shower, and hadronization is handled by PYTHIA8.240, using the CP5 tune, whereas the other sample is interfaced to PYTHIA8.301 with VINCIA. The former uses the NNPDF2.3 NNLO PDFs, whereas the latter uses the NNPDF2.3 LO PDFs. The MLM scheme [61] is used to match jets produced via matrix-element calculations with those from parton showers, using the matching p T scale of 18 GeV, which was optimized by analyzing the differential jet-rate distributions. A first sample is generated with a 2 → 2 NLO matrix element [65]. The factorization and renormalization scales are set to the p T of the underlying Born configuration.
A second sample was generated with a 2 → 3 NLO matrix element [66]. The generator-level minimal p T of the underlying Born configuration is set to 10 GeV, and the factorization and renormalization scales are set to H T /2.
with the CP5 and CH3 tunes, respectively. Both samples use the NNPDF2.3 NNLO PDFs.

SPS+DPS samples
The PYTHIA8.240 and KATIE MC event generators both produce two 2 → 2 matrix elements per event, resulting in a pure DPS sample. MPIs are also present as part of the underlying event. In KATIE, σ eff is a parameter that directly determines the size of the DPS contribution relative to the SPS cross section. A value of 21.3 mb for σ eff is adopted, as in [51]. For PYTHIA8, it is not possible to set σ eff , because it is determined by the underlying event parameters, and therefore the second 2 → 2 process is simply added with the same kinematic requirements, without any additional scaling of the cross section.
Four samples with an explicit DPS contribution are used in this paper.
• A PYTHIA8.240 sample is generated with the CP5 tune and the NNPDF2.3 NNLO PDFs. It is the aforementioned PYTHIA8 sample to which a pure DPS sample, obtained by overlaying two 2 → 2 matrix elements, is added. • The second PYTHIA8 sample, which includes an explicit DPS contribution, is the one already mentioned in Section 4.1, since the CDPSTP8S1-4j tune has been fitted to DPS-sensitive observables.
• An on-shell KATIE sample is generated with an explicit DPS contribution that is obtained by overlaying two 2 → 2 matrix elements, with the exact same generation parameters as the on-shell KATIE LO sample from the multijet models. Showering and hadronization are handled by PYTHIA8.240 with the CP5 tune and the NNPDF2.3 NNLO PDFs.
• Two off-shell KATIE samples with an explicit DPS contribution are generated using the same TMD PDFs as for the multijet samples. Since CASCADE cannot handle two 2 → 2 matrix elements, nonperturbative corrections have been derived from the onshell SPS and DPS KATIE samples, and are applied to the off-shell DPS KATIE parton level sample. The nonperturbative corrections range from 1-4% for all observables, except for the ∆S observable for which corrections up to 11% were found.

Event selection and unfolding
This analysis uses data from pp collisions at a center-of-mass energy of 13 TeV, collected during a data taking period at low luminosity, with an average pileup of 1.3 and an integrated luminosity of 42 nb −1 . The online selection of multi-jet events was based on four single-jet triggers each requiring the presence of at least one jet with a p T above 30, 50, 80, or 100 GeV, and within |η| < 4.7. Because the triggers have been prescaled, they are used in disjoint p T ranges. Offline requirements are imposed to ensure that the triggers are fully efficient, except for the trigger with the lowest p T threshold. In this last case, a correction as a function of the jet p T is applied to the selected event, effectively altering its weight. The trigger efficiency was determined by comparing the performance of the jet trigger with a minimum-bias trigger serving as an unbiased reference.
Events are selected offline by requiring exactly one primary vertex, so effects of pileup can be neglected. The correction of the event yield is based on the measurement of the average pileup and has negligible uncertainty. A systematic uncertainty due to a possible contamination of events containing two or more pp collisions is nevertheless included in the results. Two phase space regions defined by selections on jet p T are used. In Region I, the four leading jets within |η| < 4.7 are required to exceed p T thresholds of 35,30,25, and 20 GeV. Asymmetric thresholds have been chosen over symmetric ones because the latter tend to dampen the DPS contribution with respect to the SPS fraction, according to higher-order calculations or calculations performed in the k T -factorization framework [31,33]. The ∆S distribution is obtained in Region II, with p T thresholds of 50, 30, 30, and 30 GeV. On the one hand, the resolution of the ∆S observable is improved by imposing higher p T cuts. On the other hand, the ∆S observable can now be used to perform the extraction of σ eff , using the lowest jet-p T trigger threshold of 30 GeV. The second set of selections is needed to obtain the cross sections σ A and σ B , as detailed in Section 6.
The measured distributions are corrected for detector effects with the TUnfold program [73,74], which is based on a least squares fit and Tikhonov regularization [75]. The regularization is necessary to avoid possible instabilities in the inversion of the matrix describing the migrations within the phase space. Bin-to-bin migrations are kept to a minimum by choosing a bin width that is two times larger than the resolution of the considered variable, as obtained from a simulation study with PYTHIA8. The migration matrix, as well as the probabilities for migration into and out of the phase space, are obtained from the PYTHIA8 and HERWIG++ MC models.

Extraction of the effective cross section
The DPS formula (1) allows the calculation of σ eff if the DPS cross section, as well as the cross sections of the two Processes A and B, are known. In the simplest case, the Processes A and B would both be dijet production, resulting in a four-jet final state with uncorrelated jet pairs, as depicted in Fig. 1. However, initial-and final-state radiation, and higher-order interactions can produce final states with more than two jets, yielding additional possibilities to form a four-jet topology.
To avoid additional model dependencies, a DPS signal template is constructed from data, following an approach similar to the one laid out in Ref.
[25]. The Processes A and B are both defined as inclusive single-jet production. Combining two events of Type A and B will result in a multijet final state. Whenever at least four jets are found in the combined final state originating from the same process, the combined event is labeled as SPS process, otherwise the event is labeled as a DPS process.
Region II is used for the extraction of σ eff . The choice is motivated by Ref. [31] where it is suggested that such asymmetric cuts could boost the DPS signature. The cross sections of the Processes A and B are defined as inclusive single jet production with Combining two inclusive single-jet processes results in final states with at least four jets in only a fraction of the cases. A "four-jet efficiency" ( 4j ) has been obtained from the combined sample as detailed below.
From the event sample with at least one jet with p T > 30 GeV, two events are drawn at random with the second event containing at least one jet with p T > 50 GeV. The two selected events are combined to form one single event. A combined event is discarded whenever two or more jet axes spatially coincide. This veto condition is formulated as where the indices i and j indicate jets belonging to an event from the first and second data sample, respectively. The newly constructed combined event sample is then subjected to the four-jet selection criteria of Region II. The four-jet efficiency was estimated to be where the statistical uncertainty is negligible and the systematic error is detailed in the next section.
The ∆S observable is chosen for the extraction of the DPS cross section and σ eff because it is the least affected by parton shower effects. The signal template is taken from the combined data sample, which is used to extract the ∆S distribution in exactly the same manner as before, including the correction for detector effects by means of unfolding.
The SPS MC models are taken as background templates. To avoid contamination by MPI, additional samples are provided where an event is omitted if it contains a generator level parton with a p T > 20 GeV that originates from a MPI. This selection ensures that no hard jets originating from MPI enter the four-jet analysis, and will be referred to as "hard MPI removed".
The fraction of DPS events, f DPS , is extracted by performing a template fit to the unfolded ∆S distribution, obtained from the original inclusive four-jet sample. The DPS signal and SPS background ∆S distributions are both normalized to the integral of the ∆S distribution obtained from the four-jet events in data. The optimal value of the DPS fraction f DPS is determined with a maximum likelihood technique using Poisson statistics [76]: The cross section σ DPS A,B , needed for the extraction of σ eff , is then given by the integral of the ∆S distribution measured in data, scaled with the DPS fraction f DPS : Because of the overlapping p T ranges, the Processes A and B are not always distinguishable and the cross section for Process B has therefore to be rewritten as the sum of the cross section for Process A (σ A ) and the difference between the cross sections for Processes B and A (σ B − σ A ). Taking into account the correct combinatorial factor along with the four-jet efficiency, the DPS formula (1) can be reformulated as: Jet energy resolution (JER) uncertainty: The JER obtained from MC simulation differs from the one estimated for data, which would lead to a wrong estimation of the bin-to-bin migrations. An additional smearing of the jet p T at detector level is therefore applied to both MC samples used in the unfolding. To estimate the JER uncertainty, the data-to-simulation smearing factor is varied up and down with its own uncertainty, resulting in migration matrices that differ from the nominal ones. The newly obtained migration matrices are used to unfold the distributions, which are then compared with the nominal distributions for all observables. The JER uncertainty is less than 9%, except for the p T spectrum of the leading jet where it reaches a maximum of 26%.
Trigger uncertainty: An event weight as a function of jet p T is applied to data to correct for the efficiency of the trigger with the lowest threshold. These weights are obtained by fitting the trigger efficiency curve determined in data using a least-squares minimization. Varying the fit parameters by their uncertainty leads to a trigger uncertainty that never exceeds 1%.

Model uncertainty:
The data distributions are unfolded using migration matrices from the PYTHIA8 and HERWIG++ models. The averages of the two unfolded distributions are taken as the nominal unfolded distributions and the systematic error is estimated as half of the difference. The model uncertainty varies between 1% and 16%, depending on the observable.
Pileup: Events with two pp collisions in the same bunch crossing may be reconstructed with only one vertex if the collision points are separated by less than 0.12 cm along the beam axis. Taking into account the spread of vertices and the relative yields of events with 1 and 2 collisions, the pileup contamination in the data sample is estimated to be 0.28%. No further correction is applied, and a systematic uncertainty is included by reproducing all distributions with a sample of events containing two vertices, normalizing these to 0.28% of the nominal distributions, and estimating the effect of such pileup correction on the data. The systematic uncertainty is smaller than 1% in all bins for all observables.

Integrated luminosity uncertainty:
The uncertainty in the integrated luminosity for data collected in 2016 is 1.2% [77].
Because the four-jet efficiency is determined using uncorrected data, it has neither JER or model systematic uncertainty. However, an additional systematic uncertainty is included to cover a possible difference with respect to the true efficiency to be applied on the corrected cross sections. This uncertainty is determined by examining the four-jet efficiency obtained with PYTHIA8 and HERWIG++ at both the detector and generator levels, after a 2-dimensional reweighting as a function of leading-jet p T and jet multiplicity to obtain a better description of the data. A four-jet efficiency of 0.404 and 0.412 is found with PYTHIA8 on detector and generator level, respectively, whereas values of 0.403 and 0.392 are obtained with HERWIG++, with negligible statistical uncertainty. The systematic uncertainty is therefore conservatively estimated to be smaller than 2%.

Results
The total cross sections in the two phase space regions defined by thresholds on the p T of the four leading jets, Region I and Region II, are obtained by integrating the differential cross section as a function of the leading jet η and the ∆S observable, respectively, yielding: Table 2: Systematic uncertainties, along with the statistical and the total uncertainties for the cross sections of the two phase space regions, along with the observables needed for the extraction of σ DPS A,B , in percent. The JES uncertainty leads to asymmetric uncertainties (an upper and a lower error): all other systematic uncertainties, as well as the statistical uncertainty, are symmetric. An additional uncertainty in 4j because of possible differences between generatorand detector-level events, is estimated to be 2%.
Tables 3-5 compares the values measured in data with the ones obtained from MC event generators. A discussion of the total and differential cross sections for each of the sets of models introduced in Section 4 is presented below.
The cross sections of the inclusive single-jet Processes A and B are determined by integrating the leading jet η distribution for both processes, and are: A large increase in cross section is observed when lowering the p T threshold from 50 GeV to 30 GeV, as expected. These cross sections will be used as input to Eq. (1) for the determination of the DPS cross section and σ eff along with the four-jet efficiency from Eq. (2).

The PYTHIA8 and HERWIG models
The models based on LO 2 → 2 matrix elements, PYTHIA8, HERWIG++, and HERWIG7, respectively labeled as P8, H++, and H7 in the figures, are compared with data. Table 3 gives an overview of the event generators, tunes, and PDF sets, along with their respective cross sections. All LO 2 → 2 models predict cross sections that are much larger than the measured ones; especially PYTHIA8 with the CDPSTP8S1-4j tune predicts a cross section that is roughly 2.5 times larger than the one observed in data. Figures 2-5 show a comparison of the data to various MC models as a function of p T , η, and the DPS-sensitive observables. Three PYTHIA8 models and one HERWIG7 model are shown in direct comparison with the data. These models employ the most recent CP5 and CH3 tunes, the dedicated DPS tune or combine PYTHIA8 with a dipole-antenna shower provided by VINCIA, while all of the models are represented in the ratio plots.
The p T spectra in Fig. 2, obtained for Region I, show that the much larger integrated cross section of the MC models, compared with the data, comes from an abundance of low-p T jets, whereas for p T 100 GeV the models show agreement within 50% of the data; for HERWIG7 it is even within the total uncertainty. Fig. 3, the η spectra, shows that a large part of the excess of low-p T jets is located in the forward η regions. and ∆Y are normalized to the average of their first four bins. Normalizing to the average of multiple bins reduces the effect of statistical fluctuations. The distributions in φ ij , ∆p T,Soft , and ∆S are all normalized to their last bin, since these bins already have a small relative statistical uncertainty.
The ∆φ Soft and ∆p T,Soft distributions are relatively well described by all LO 2 → 2 models. Deviations from data never exceed 20%, albeit being larger than the total uncertainty in the data points in certain bins for some of the models. Similar results are observed for the predictions of the ∆p T,Soft observable from the DPS tune in [15]. Deviations of 10-20% occur between a model employing a similar DPS tune (CDPSTP8S1-WJ) and the data in the Z+jets final state.
The shape of the ∆Y distribution predicted by all the LO 2 → 2 models differs significantly from data, and the MC-to-data ratio increases towards higher values of ∆Y. The overshoot at large ∆Y is consistent with the excess of low-p T forward jets.
A distinction between two classes of models becomes apparent in the ∆φ min 3j and φ ij distributions. The models implementing a p T -ordered parton shower describe the distribution in ∆φ min 3j well, although yielding a distribution in φ ij that is more uncorrelated than observed in data. The slope of the ∆φ min 3j distribution obtained from the models with a p T -ordered shower algorithm going to zero, overshoots the slope of the distribution obtained from data. For models that use an angular-ordered or dipole-antenna parton shower, the ∆φ min 3j correlation is too strong. The slope of the distributions ∆φ min 3j distributions obtained from the models with an angular-ordered shower overshoot the distribution obtained from data when going to π, whereas the shape of the data is described more accurately by φ ij . The PYTHIA8+VINCIA model confirms that the parton shower algorithm is responsible for the different tendencies observed for the two classes of models. This observation makes these observables less suitable to untangle SPS and DPS contributions to the cross section.       The ∆S distribution is less affected by different parton shower implementations. The DPS tune CDPSTP8S1-4j agrees very well with the shape of the data, but lies slightly above the data at low values of ∆S, pointing to a potential overestimation of the DPS contribution. All other models underestimate the data at low ∆S, indicating a possible need for more DPS to obtain a proper description of the shape of the ∆S distribution.

Multijet models
Data distributions are also compared with the multijet samples that are obtained from models based on LO 2 → n(n ≥ 2) and NLO 2 → 2 and 2 → 3 matrix elements. The group of multijet samples includes 2 → 4 on-shell and off-shell predictions made by KATIE, two MADGRAPH5 aMC@NLO LO samples for which 2 → 2, 3, 4 matrix element are all included, a MADGRAPH5 aMC@NLO NLO 2 → 2 sample, and two POWHEG NLO samples that use a 2 → 2 and a 2 → 3 matrix element. Table 4 gives a complete overview of all models, tunes, and PDF sets, along with their respective cross sections. In the figures, the labels KT, PW, and MG5 are used for KATIE, POWHEG, and MADGRAPH5 aMC@NLO, respectively. Figures 6-9 show a comparison of the data to various MC models as a function of p T , η, and the DPS-sensitive observables. Four models are shown in direct comparison with the data. These models include one of the two on-and off-shell KATIE models, the MADGRAPH5 aMC@NLO LO sample interfaced with the CP5 tune and the POWHEG NLO 2 → 2 sample, while all of the models are represented in the ratio plots. The predicted cross sections obtained from the on-shell KATIE samples interfaced with PYTHIA8 and HERWIG7 are larger than the cross sections obtained from data. They sharply decrease when an off-shell matrix element is used, showing agreement within the data uncertainty for Region I. The MADGRAPH5 aMC@NLO LO samples and all the NLO samples predict cross sections that are roughly in agreement with the cross sections obtained from data in Region I, but are larger than those from data in Region II, as are all the KATIE cross sections in the same region. Fig. 6 compares the p T spectra of the various models with the data. The on-shell KATIE predictions agree with the data in the first bin of each of the p T spectra, but are above the data at higher p T . This may be explained because most jets originate from the 2 → 4 matrix element and not from the parton shower. The same effect, but less pronounced, is observed for the off-shell KATIE curves. The different PDF sets used with off-shell KATIE result in small variations. A better description of the p T spectra is given by the MADGRAPH5 aMC@NLO LO sample, with a p T -ordered parton shower. In this case, some jets must originate from the parton shower, yielding a softer spectrum. The combination of the MADGRAPH5 aMC@NLO LO 2 → 2, 3, 4 sample with the dipole-antenna showering from PYTHIA8+VINCIA, results in a lowering of the total cross section. All NLO models give a similar description as the MADGRAPH5 aMC@NLO sample; the higher-order matrix element, including virtual corrections, contributes to a lower cross section. A comparison of the multi jet data samples to the standard PYTHIA8 and HER-WIG curves from the previous section demonstrates that NLO corrections and the inclusion of multi-leg matrix elements improve the description of the p T spectra.
The η spectra are shown in Fig. 7. The central region is consistently described by all models, even the on-shell KATIE models; the overall cross section is too large but the ratio remains flat for |η| ≤ 3.0. An excess of jets is observed in the forward region, although this is less pronounced than in the case of the PYTHIA8 and HERWIG models. The excess is also strongest for the leading jet and diminishes for the second, third and fourth leading jet, yielding a good description of the shape for the latter.
Differential cross sections for all other observables are shown in Figs. 8 and 9. As before, these distributions have been normalized to a region with a much reduced DPS contribution.
The distributions in ∆φ Soft and ∆p T,Soft demonstrate that most multijet models leave room for an additional DPS contribution in the regions where this can be expected. The exceptions are the MADGRAPH5 aMC@NLO LO 2 → 2, 3, 4 distributions that describe the shape of the distributions from data reasonably well. This last model contains a LO 2 → 2 contribution, and, as observed in Sec 8.1, the LO 2 → 2 models describe the ∆φ Soft and ∆p T,Soft distributions well.
The ∆Y distributions show that some multijet models still have an excess of jets at large rapidity separation.
The distributions in φ ij confirm the strong dependence on the parton shower implementation observed with the LO 2 → 2 models (as shown, e.g., the MADGRAPH5 aMC@NLO LO 2 → 2, 3, 4 predictions for φ ij ). Angular-ordered parton showers, reproduced by VINCIA and KATIE, describe the shape of the data distributions better. The off-shell KATIE predictions show a depletion in the region where additional DPS is expected. The effect of the parton shower on the ∆φ min 3j distributions is less pronounced in the multijet models compared to the LO 2 → 2 models.
The ∆S distributions again show a more robust behavior with respect to the parton shower implementation. The LO MADGRAPH5 aMC@NLO samples leave less room for an additional DPS contribution compared with the NLO models. Due to the sole use of a 2 → 4 matrix element, the p T spectra are far too hard, resulting in a large overestimation of the slope for all the KATIE models.

SPS+DPS Models
The last group of model predictions that are compared with data are those including an explicit DPS contribution, as described in Section 4.3. The PYTHIA8 CDPSTP8S1-4j predictions presented in this section are the same as in previous sections; there is no need to add a DPS contribution explicitly, because the underlying event parameters are specifically tuned to in-

KT Onshell+CP5
KT Onshell+CH3  A comparison of the integrated cross sections obtained from the SPS+DPS samples in Table 5 with those from pure SPS samples, shows the expected increase. Only the values of the off-shell KATIE models agree with the measured cross section in Region I within the uncertainty. The cross sections predicted for Region II are all too high.  Figures 10 and 11 show the p T and η spectra of the four leading jets, respectively. A comparison with the spectra of the pure SPS samples from the previous sections shows that the DPS samples contribute in the low-p T region, and mostly in the forward regions of the η spectra.

KT Offshell+MRW KT Offshell+PBTMD
Normalized distributions in the DPS-sensitive observables are shown in Fig. 12. Small differences in shape, typically ∼ 5%, occur in the DPS-sensitive regions when comparing the SPS+DPS PYTHIA8 sample, interfaced with the CP5 tune, with its pure SPS counterpart in Figs. 4 and 5. The PYTHIA8 samples give a good description of shape of the distributions in the observables ∆φ Soft , ∆φ min 3j , ∆p T,Soft and ∆S, with the CDPSTP8S1-4j tune performing slightly better than the PYTHIA8 SPS+DPS sample with the CP5 tune. This might be expected since the CDPSTP8S1-4j tune was obtained by fitting the PYTHIA8 underlying event parameters to distributions in ∆φ Soft , ∆p T,Soft , and ∆S measured with data at a center-of-mass energy of 7 TeV, as detailed in [51].
The KATIE curves show a more noticeable increase where a DPS contribution is expected, especially in the case of the off-shell samples. This is in contrast with the performance of the SPS-only predictions (see Figs. 8 and 9) where the off-shell KATIE models underestimate the data. This shows that a DPS contribution is needed in the KATIE model to improve the description of the data, but that the value of σ eff = 21.3 mb (taken from Ref. [51]), is too small for this model based on a 2 → 4 matrix element.
All models fail to describe the shape of the ∆Y observable, and the φ ij observable shows the same tendencies found in the previous sections. The LO 2 → 2PYTHIA8 models using a p Tordered parton shower show too large a decorrelation, whereas predictions using higher-order matrix element calculations perform better, demonstrated by the KATIE distributions.

Extraction of the effective cross section
As demonstrated in the above sections, the ∆S observable exhibits the most robust sensitivity to DPS. Other observables are either less sensitive to DPS or suffer from large variations induced by different parton shower models. Therefore, the ∆S distribution is used in the extraction of σ eff . The DPS cross section is determined and the effective cross section, σ eff , is extracted using different SPS MC event samples with and without the hard MPI removed following the template method laid out in Section 6. The PYTHIA8 sample with the CUETP8M1 tune yields a DPS  Figure 12: Comparison of the distributions in DPS-sensitive observables obtained from data to different SPS+DPS KATIE (KT) and PYTHIA8 (P8) models. All distributions have been determined in Region I, except for the ∆S distribution which has been measured in Region II. All distributions have been normalized to the region where a reduced DPS contribution is expected. The error bars represent the statistical uncertainty, and the yellow band indicates the total (sta-tistical+systematic) uncertainty on the measurement.
fraction lower compared to all other tunes, while the HERWIG++ sample with the CUETHS1 tune gives similar results as when HERWIG7 with the CH3 tune is used, and therefore both are omitted below. The PYTHIA8 sample interfaced with the CDPSTP8S1-4j tune without the hard MPI removed already contains a DPS contribution, scaled with an effective cross section equal to 21.3 mb [51]. Therefore, the extracted DPS fraction for this model is expected to be close to zero. The extraction is still performed as a consistency check and to test the performance of the tune at a center-of-mass energy of √ s = 13 TeV. The multijet KATIE models are not considered; they significantly overshoot the DPS-sensitive slope of the ∆S observable and using them would result in a negative DPS contribution.
The mixed event sample obtained from data, as described in Section 6, is used to get the ∆S distributions for pure DPS events and is corrected by means of unfolding in the same way as the other distributions. This corrected ∆S distribution is shown in Fig. 13, along with the ∆S distributions obtained from the DPS component of the PYTHIA8 and on-shell KATIE generators, both interfaced with the CP5 tune. All distributions have been normalized to unity. The distributions show a maximum at ∆S = π because of jet pairs with overlapping p T values such that the two softest and the two hardest jets do not coincide with the jet pairs from separate events or parton collisions. The DPS data sample exhibits a larger decorrelation as compared to the DPS MC samples, which can be attributed to disparities in the p T spectra observed in data and for MC events. The template fitter program [76] takes the SPS MC distributions along with the DPS data sample as input, and uses the template method to determine the fraction of DPS events, f DPS . The results for the extracted values of f DPS are shown in Table 6. Equations (5) and (4) are used to determine the DPS cross sections and values for σ eff shown in Tables 7 and 8, respectively. The results for both sets of samples (with and without the hard MPI removed) are shown, along with the net difference between both cross sections since this can be interpreted as the amount of DPS inherent to the tune. An example of the template fit is presented in Fig. 14, where the fitted distributions using the POWHEG NLO 2 → 2 model without the hard MPI removed along with its statistical and systematic uncertainties are shown. The DPS cross section obtained for all models with and without the hard MPI removed range from 14.6 to 70 mb, yielding values for σ eff between 7.7 and 34.8 mb. The exception is PYTHIA8 with the CDPSTP8S1-4j tune without the hard MPI removed where an excess of DPS events is found, resulting in a negative DPS cross section. Therefore, the effective cross section for this model has not been calculated because it would yield a nonphysical result.

CMS
In the case of the LO 2 → 2 models, 2 final state jets originate from the ME, and 2 additional jets stem from the parton shower. A distinction between two groups of models using a LO 2 → 2 matrix element can be made. On the one hand, the PYTHIA8 and HERWIG7 models using the CP5 and CH3 tunes, respectively, use the latest underlying event tune and up-to-date PDFs.
On the other hand, both PYTHIA8+VINCIA and HERWIG7 with the SoftTune rely on an older underlying event tune and older PDFs. The extracted values for the DPS cross section for the former two models are 38.5 and 38.2 mb, whereas for the latter two models values of 23.3 and 29.5 mb were obtained where events containing one or more hard partons originating from MPI have been removed. The results indicate that it might be the different underlying event tunes and the usage of older PDFs that are responsible for the more DPS-like topology in both PYTHIA8+VINCIA and HERWIG7 with the SoftTune.
Introducing higher multiplicity MEs reduces the effect of the underlying event tune, parton showers, and the PDFs since more jets will originate from the ME. The dominating contribution to the event sample comes from 2 → 4 ME, where the 4 final state jets all stem from the ME. The extracted DPS cross section with hard MPI removed is 31 mb for both PYTHIA8 with the CP5 tune and PYTHIA8+VINCIA, confirming that the effect of the different underlying event tunes, parton showers, and PDFs is suppressed. In the case of PYTHIA8 with the CP5 tune, the higher multiplicity MEs reduce the need for additional DPS, whereas in the case of PYTHIA8+VINCIA more room for additional DPS exists. For both models, about half of the DPS cross section that is needed to describe the data can be covered by the MPI that are intrinsic to the tune.
The NLO 2 → 2 models contain a combination of exclusive NLO 2 → 2 and LO 2 → 3 ME. One or more parton-shower jets are required to produce a 4-jet final state, with a ME/partonshower matching algorithm applied to avoid double counting. Large DPS cross sections of up to 70 mb are needed to describe the data if events containing hard MPI partons originating from the underlying event description, provided by the CP5 tune, are removed. Compared with the other models, the NLO 2 → 2 models are outliers. If events with hard MPI partons from the underlying event description are included, they can again account for about half of the DPS cross section needed to describe the data.
For the POWHEG NLO 2 → 3 sample, NLO 2 → 3 and LO 2 → 4 MEs are combined, leaving room for zero or one additional jets from the parton shower. The extracted DPS cross section approaches the one obtained for the models with a LO 2 → 2, 3, 4 ME. Contrary to the LO 2 → 2, 3, 4 ME models, however, the NLO 2 → 3 models do not allow events containing hard MPI partons to contribute as much to the DPS cross section. , where a σ eff equal to 14.9 +1.2 −1.0 (stat) +5.1 −3.8 (syst) mb was reported, whereas none agree with the value of 21.3 +1.2 −1.6 mb from the CMS measurement at a center-of-mass energy of 7 TeV [51], which is more in line with the results obtained with some of the models based on older UE tunes.
Two other DPS measurements have been performed at a center-of-mass energy √ s = 13 TeV. A value of σ eff equal to 12.7 +5.0 −2.9 mb has been extracted from the same-sign WW measurement in Ref. [14]. A σ eff of 7.3 ± 0.5 (stat) ± 1.0 (syst) mb has been obtained from the J/ψ pair production measurement in Ref. [20]. It has been shown that σ eff is expected to be process independent for inclusive final states [78], therefore, it is noteworthy that the result the J/ψ meson pair production measurement shows only agreement with the NLO 2 → 2 models. The extracted σ eff of the same-sign WW measurement shows agreement with all models due to the size of the errors in the measurement, indicating that further measurements in this channel are desirable before any further conclusions can be made.

Summary
A study of the inclusive production of four-jet events at low transverse momentum has been presented based on data from proton-proton collisions collected with the CMS detector at a center-of-mass energy of 13 TeV. Various observables sensitive to double-parton scattering (DPS) are studied and values for its effective cross section have been extracted. Table 6: The values of the DPS fraction f DPS extracted from data using different SPS models, along with their statistical and systematic uncertainties. The results are shown for the model where the full tune is used, and for the same models where the hard MPI have been removed. The last column shows the net difference between the two first columns, and is interpreted as the fraction of DPS inherent to the tune.

MC Model
Tune   Models based on leading order (LO) 2 → 2 matrix elements significantly overestimate the absolute four-jet cross section in the phase space domains studied in this paper. This excess is related to an abundance of low-p T and forward jets. The predictions of the absolute cross section generally improve when next-to-leading order (NLO) and/or higher-multiplicity matrix elements are used.
The azimuthal angle between the jets with the largest separation in η, φ ij , has a strong discriminating power for different parton-shower approaches and the data favor the angular-ordered and dipole-antenna parton-shower models over those with a p T -ordered parton shower. The yield of jet pairs with large rapidity separation ∆Y is, however, overestimated by all models, although models based on NLO and/or higher-multiplicity matrix elements are closer to the data.
The distribution of the minimal combined azimuthal angular range of three jets, ∆φ min 3j , also exhibits sensitivity to the parton-shower implementation, with data favoring p T -ordered parton showers with the LO 2 → 2 models for this observable. In the case of models based on NLO and/or higher-multiplicity matrix elements the comparisons are less conclusive.
Other observables, such as the azimuthal angle between the two softest jets, ∆φ Soft , and their transverse momentum balance, ∆p T,Soft , indicate the need for a DPS contribution in the models to various degrees, as confirmed by the extracted values of σ eff . The distribution of the azimuthal angle between the hard and soft jet pairs, ∆S, is the least sensitive to the details of the parton-shower modeling, and it is used for the extraction of the effective cross section, σ eff .
A dependence is observed in the extracted values of σ eff in the model used to describe the SPS contribution. Models based on NLO 2 → 2 matrix elements yield the smallest (∼7 mb) values of σ eff and need the largest DPS contribution. However, models using a 2 → 2 matrix element along with older underlying event descriptions and older PDFs, tend to need the smallest DPS contribution. The sensitivity to the underlying event description, parton showers, and the PDFs is observed to be small when including higher-order matrix elements, since both models using the 2 → 2, 3, 4 matrix elements show agreement with each other.
These results demonstrate the need for further development of models to accurately describe final states with multiple jets in phase space regions with large potential DPS contributions.
[24] I. Sadeh, "Double parton scattering in four-jet events in pp collisions at 7 TeV with the ATLAS experiment at the LHC". PhD thesis, Tel Aviv U., 2013. arXiv:1308.0587.