Study of hard double-parton scattering in four-jet events in $pp$ collisions at $\sqrt{s} = 7$ TeV with the ATLAS experiment

Inclusive four-jet events produced in proton-proton collisions at a centre-of-mass energy of $\sqrt{s} = 7$ TeV are analysed for the presence of hard double-parton scattering using data corresponding to an integrated luminosity of 37.3 pb$^{-1}$, collected with the ATLAS detector at the LHC. The contribution of hard double-parton scattering to the production of four-jet events is extracted using an artificial neural network, assuming that hard double-parton scattering can be approximated by an uncorrelated overlaying of dijet events. For events containing at least four jets with transverse momentum $p_{\mathrm{T}} \geq 20$ GeV and pseudorapidity $\eta \leq 4.4$, and at least one having $p_{\mathrm{T}} \geq 42.5$ GeV, the contribution of hard double-parton scattering is estimated to be $f_{\mathrm{DPS}} = 0.092 ^{+0.005}_{-0.011} (\mathrm{stat.}) ^{+0.033}_{-0.037} (\mathrm{syst.})$. After combining this measurement with those of the inclusive dijet and four-jet cross-sections in the appropriate phase space regions, the effective overlap area between the interacting protons, $\sigma_{\mathrm{eff}}$, was determined to be $\sigma_{\mathrm{eff}} = 14.9 ^{+1.2}_{-1.0} (\mathrm{stat.}) ^{+5.1}_{-3.8} (\mathrm{syst.})$ mb. This result is consistent within the quoted uncertainties with previous measurements of $\sigma_{\mathrm{eff}}$, performed at centre-of-mass energies between 63 GeV and 8 TeV using various final states, and it corresponds to $21^{+7}_{-6}$% of the total inelastic cross-section measured at $\sqrt{s} = 7$ TeV. The distributions of the observables sensitive to the contribution of hard double-parton scattering, corrected for detector effects, are also provided.


Introduction
Interactions involving more than one pair of incident partons in the same collision have been discussed on theoretical grounds since the introduction of the parton model to the description of particle production in hadron-hadron collisions [1][2][3]. These first studies were followed by the generalization of the Altarelli-Parisi evolution equations to the case of multi-parton states in Refs. [4,5] and a discussion of possible correlations in the colour and spin degrees of freedom of the incident partons [6]. In the first phenomenological studies of such effects, the most prominent role was played by processes known as double-parton scattering (DPS), which is the simplest case of multi-parton interactions (MPI), leading to final states such as four leptons, four jets, three jets plus a photon, or a leptonically decaying gauge boson accompanied by two jets [7][8][9][10][11][12][13][14][15]. These studies have been supplemented by experimental measurements of DPS effects in hadron collisions at different centre-of-mass energies, which now range over two orders of magnitude, from 63 GeV to , and which have firmly established the existence of this mechanism. The abundance of MPI phenomena at the LHC and their importance for the full picture of hadronic collisions have reignited the phenomenological interest in DPS and have led to a deepening of its theoretical understanding [31][32][33][34][35][36][37][38][39]. Despite this progress, quantitative measurements of the effect of DPS on distributions of observables sensitive to it are affected by large systematic uncertainties. This is a clear indication of the experimental challenges and of the complexity of the analysis related to such measurements. Therefore, the cross-section of DPS continues to be estimated by ignoring the likely existence of complicated correlation effects. For a process in which a final state A + B is produced at a hadronic centre-of-mass energy √ s, the simplified formalism of Refs. [12,13] yields The quantity δ AB is the Kronecker delta used to construct a symmetry factor such that for identical final states with identical phase space, the DPS cross-section is divided by two. The σ eff , usually referred to as the effective cross-section, is a purely phenomenological parameter describing the effective overlap of the spatial distribution of partons in the plane perpendicular to the direction of motion. In hadronic collisions it was typically found to range between 10 and 25 mb [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30]. In Eq. (1), the variousσ are the parton-level cross-sections, either for the DPS events, indicated by the subscript A + B, or for the production of a final state A or B in a single parton scatter (SPS), given by dσ A (s) = 1 2s i j dx 1 dx 2 f i (x 1 , µ F ) f j (x 2 , µ F ) dΦ A |M i j→A (x 1 x 2 s, µ F , µ R )| 2 .
Here the functions f i (x, µ F ) are the single parton distribution functions (PDFs) which at leading order parameterize the probability of finding a parton i at a momentum fraction x at a given factorization scale µ F in the incident hadron; dΦ A is the invariant differential phase-space element for the final state A; M is the perturbative matrix element for the process i j → A; and µ R is the renormalization scale at which the couplings are evaluated. To constrain the phase space to that allowed by the energy of each incoming proton, a simple two-parton PDF is defined as where Θ(x) is the Heaviside step function, Γ(b) the area overlap function, and the x and scale dependence of the PDF are assumed to be independent of the impact parameter b. Equation (3) reflects the omission of correlations between the partons in the proton. At high energy, Eq. (1) can be derived using Eq. (3) by neglecting the contribution of the step function.
Typically, the main challenge in measurements of DPS is to determine if the A+B final state was produced in an SPS via the 2 → 4 process or in DPS through two independent 2 → 2 interactions. In one of the first studies of DPS in four-jet production at hadron colliders [10] the kinematic configuration in which there is a pairwise balance of the transverse momenta (p T ) of the jets was identified as increasing the contribution of the DPS mechanism relative to the perturbative QCD production of four jets in SPS. The idea is that in typical 2 → 2 scattering processes the two outgoing particles -here the partons identified as jets -are oriented back-to-back in transverse plane such that their net transverse momentum is zero. Corrections to this simple picture include initial-and final-state radiation as well as fragmentation and hadronization. In addition, recoil against the underlying event can modify the four-momentum of the overall final-state particle configuration. In attempting to describe all of these features, Monte Carlo (MC) event generators form an integral part, providing a link between the experimentally observed jets and the simple partonic picture of DPS as two almost independent 2 → 2 scatters.
An analysis of inclusive four-jet events produced in proton-proton collisions at a centre-of-mass energy of √ s = 7 TeV at the LHC and collected during 2010 with the ATLAS detector is presented here. The topology of the four jets is exploited to construct observables sensitive to the DPS contribution. The DPS contribution to the four-jet final state is estimated and combined with the measured inclusive dijet and four-jet cross-sections in the appropriate phase space regions to determine σ eff . The normalized differential four-jet cross-sections as a function of DPS-sensitive observables are measured and presented here as well.

Analysis strategy
To extract σ eff in the four-jet final state, Eq. (1) is rearranged as follows. The differential cross-sections in Eq. (1) are rewritten for the four-jet and dijet final states and integrated over the phase space defined by the selection requirements of the dijet phase space regions A and B. This yields the following expression for the DPS cross-section in the four-jet final state: where σ A 2j and σ B 2j are the cross-sections for dijet events in the phase space regions labelled A and B respectively. The assumed dependence of the cross-sections and σ eff on s is omitted for simplicity. The DPS cross-section may be expressed as where σ 4j is the inclusive cross-section for four-jet events in the phase-space region A ⊕ B, including all four-jet final states, namely both the SPS and DPS topologies, and where f DPS represents the fraction of DPS events in these four-jet final states. The expression for σ eff then becomes, To extract σ eff , it is therefore necessary to measure three cross-sections, σ A 2j , σ B 2j and σ 4j , and estimate f DPS .
The four-jet and dijet final states are defined inclusively [40,41] such that at least four jets or two jets respectively are required in the event, while no restrictions are applied to additional jets. When measuring the cross-section of n-jet events, the leading (highest-p T ) n jets in the event are considered. The general expression for the measured four-jet and dijet cross-sections may be written as where the subscript nj denotes either dijet (2j) or four-jet (4j) topologies. For each nj channel, N nj is the number of observed events, C nj is the correction for detector effects, particularly due to the jet energy scale and resolution, and L nj is the corresponding proton-proton integrated luminosity.
The DPS model contributes in two ways to the production of events with at least four jets, leading to two separate event classifications. In one contribution, the secondary scatter produces two of the four leading jets in the event; such events are classified as complete-DPS (cDPS). In the second contribution of DPS to four-jet production, three of the four leading jets are produced in the hardest scatter, and the fourth jet is produced in the secondary scatter; such events are classified as semi-DPS (sDPS). The DPS fraction is therefore rewritten as f DPS = f cDPS + f sDPS , and f cDPS and f sDPS are both determined from data. The dijet cross-sections in Eq. (6) do not require any modification since they are all inclusive cross-sections, i.e., the three-jet cross-section accounting for the production of an sDPS event is already included in the dijet cross-sections.
Denoting the observed cross-section at the detector level by and the ratio of the corrections for detector effects by yields the expression from which σ eff is determined, The main challenge of the measurement is the extraction of f DPS = f cDPS + f sDPS from optimally selected measured observables. An artificial neural network (NN) is used for the classification of events [42], using as input various observables sensitive to the contribution of DPS. The differential distributions of these observables are also presented here.

The ATLAS detector
The ATLAS detector is described in detail in Ref. [43]. In this analysis, the tracking detectors are used to define candidate collision events by constructing vertices from tracks, and the calorimeters are used to reconstruct jets.
GeV. This SPS sample was compared to the SPS sample extracted from the AHJ sample for validation purposes.
In addition, a sample of multi-jet events was generated with Pythia 6.425 [57] using a 2 → 2 matrix element at leading order with additional radiation modelled in the leading-logarithmic approximation by p T -ordered parton showers. The sample was generated utilizing the modified leading-order PDF set MRST LO* [58] with the AMBT1 [59] tune.
To account for the effects of multiple proton-proton interactions in the LHC (pile-up), the multi-jet events were overlaid with inelastic soft QCD events generated with Pythia 6.423 using the MRST LO* PDF set with the AMBT1 tune. All the events were processed through the ATLAS detector simulation framework [60], which is based on Geant4 [61]. They were then reconstructed and analysed by the same program chain used for the data.

Data set and event selection
The measurement presented here is based on the full ATLAS 2010 data sample from proton-proton collisions at √ s = 7 TeV. The trigger conditions evolved during the year with changing thresholds and prescales. A full description of the trigger strategy, developed and used for the measurement of the dijet cross-section using 2010 data, is given in Ref. [62]. For the events in the samples used in this study, the trigger was fully efficient. In total, the data used correspond to a luminosity of 37.3 pb −1 , with a systematic uncertainty of 3.5% [63]. This data set was chosen because it has a low number of protonproton interactions per bunch crossing, averaging to approximately 0.4. It was therefore possible to collect multi-jet events with low p T thresholds and to efficiently select events with exactly one reconstructed vertex (single-vertex events), thereby removing any contribution from pile-up collisions to the four-jet final-state topologies.
To reject events initiated by cosmic-ray muons and other non-collision backgrounds, events were required to have at least one reconstructed primary vertex, defined as a vertex that is consistent with the beam spot and is associated with at least five tracks with transverse momentum p track T > 150 MeV. The efficiency for collision events to pass these requirements was over 99%, while the contribution from fake vertices was negligible [62,64].
Jets were identified using the anti-k t jet algorithm [65], implemented in the FastJet [66] package, with radius parameter R = 0.6. The inputs to jet reconstruction are the energies in three-dimensional topological clusters [67,68] built from calorimeter cells, calibrated at the electromagnetic (EM) scale. 2 A jet energy calibration was subsequently applied at the jet level, relating the jet energy measured with the ATLAS calorimeter to the true energy of the stable particles entering the detector. A full description of the jet energy calibration is given in Ref. [64]. For the MC samples, particle jets were built from particles with a lifetime longer than 30 ps in the Monte Carlo event record, excluding muons and neutrinos.
For the purpose of measuring σ eff in the four-jet final state, three samples of events were selected, two dijet samples and one four-jet sample. The former two samples have at least two, and the latter at least four, jets in the final state, where each jet was required to have p T ≥ 20 GeV and |η| ≤ 4.4. In each event, jets were sorted in decreasing order of their transverse momenta. The transverse momentum of the i th jet is denoted by p i T and the jet with the highest p T (p 1 T ) is referred to as the leading jet. To ensure 100% trigger efficiency, the leading jet in four-jet events was required to have p 1 T ≥ 42.5 GeV. The selection requirements for the dijet samples were dictated by those used to select four-jet events. In one class of dijet events, the requirement on the transverse momentum of the leading jet must be equivalent to the requirement on the leading jet in four-jet events, p 1 T ≥ 42.5 GeV. The other type of dijet event corresponds to the sub-leading pair of jets in the four-jet event, with a requirement of p T ≥ 20 GeV. In the following, the cross-section for dijets selected with p 1 T ≥ 20 GeV is denoted by σ A 2j and the crosssection for dijets with p 1 T ≥ 42.5 GeV is denoted by σ B 2j . To summarize, the measurement was performed using the dijet A sample and its two subsamples (dijet B and four-jet), selected using the following requirements: where N jet denotes the number of reconstructed jets. All of the selected events were corrected for jet reconstruction and trigger inefficiencies, the corrections ranging from 2%-4% for low-p T jets to less than 1% for jets with p T ≥ 60 GeV. The observed distributions of the p T and y of the four leading jets in the events are shown in Figures 1(a) and 1(b) respectively.

Correction for detector effects
The correction for detector effects was estimated separately for each class of events using the Pythia6 MC sample. The same restrictions on the phase space of reconstructed jets, defined in Eq. (11), were applied to particle jets. The correction is given by where N A,B reco nj (N A,B particle nj ) is the number of n-jet events passing the A-or-B selection requirements using reconstructed (particle) jets.
This correction is sensitive to the migration of events into and out of the phase space of the measurement. Due to the very steep jet-p T spectrum in dijet and four-jet events, it is crucial to have good agreement between the jet p T spectra in data and in MC simulation close to the selection threshold before calculating the correction. Therefore, the jet p T threshold was lowered to 10 GeV and the fiducial |η| range was increased to 4.5 for both the reconstructed and particle jets, and the MC events were reweighted such that the jet p T -y distributions reproduced those measured in data. The value of α 4j 2j (see Eq. (9)), as determined from the reweighted MC events, is where the uncertainty is statistical. The systematic uncertainties are discussed in Section 7.

Determination of the fraction of DPS events
The main challenge in the measurement of σ eff is to estimate the DPS contribution to the four-jet data sample. It is impossible to extract cDPS and sDPS candidate events on an event-by-event basis. Therefore, the usual approach adopted is to fit the distributions of variables sensitive to cDPS and sDPS in the data to a combination of templates for the expected SPS, cDPS and sDPS contributions. The template for the SPS contribution is extracted from the AHJ MC sample, while the cDPS and sDPS templates are obtained by overlaying two events from the data. In addition to assuming that the two interactions producing the four-jet final state in a DPS event are kinematically decoupled, the analysis relies on the assumption that the SPS template from AHJ properly describes the expected topology of four-jet production in a single interaction. The latter assumption is supported by the observation of good agreement between various distributions in the SPS samples in AHJ and in Sherpa. To exploit the full spectrum of variables sensitive to the various contributions and their correlations, the classification was performed with an artificial neural network.

Template samples
Differences were observed when comparing the p T and y distributions in data with those in AHJ. Therefore, before extracting template samples, the events in the four-jet AHJ sample selected with the requirements detailed in Eq. (11) are reweighted such that they reproduce the distributions in data.
In events generated in AHJ, the outgoing partons can be assigned to the primary interaction from the Alpgen generator or to a secondary interaction, generated by Jimmy, based on the MC generator's event record. The former are referred to as primary-scatter partons and the latter as secondary-scatter partons.
The p T of secondary-scatter partons was required to be p T ≥ 15 GeV in order to match the minimum p T of primary-scatter partons set by the MLM matching scale in AHJ. Once the outgoing partons were classified, the jets in the event were matched to outgoing partons and the event was classified as an SPS, cDPS or sDPS event.
The matching of jets to partons is done in the φ-y plane by calculating the angular distance, ∆R parton−jet , between the jet and the outgoing parton as For 99% of the primary-scatter partons, the parton can be matched to a jet within ∆R parton−jet ≤ 1.0, which was therefore used as a requirement for the matching of jets and partons. Jets were first matched to primary-scatter partons and the remaining jets were matched to secondary-scatter partons.
Events in which none of the leading four jets match a secondary-scatter parton were assigned to the SPS sample. This method of obtaining an SPS sample is preferred over turning off the MPI module in the generator since it retains all of the soft MPI and underlying activity in the selected SPS events. Events were classified as cDPS events if two of the four leading jets match primary-scatter partons and the other two match secondary-scatter partons. Events in which three of the leading jets match primary-scatter partons and the fourth jet matches a secondary-scatter parton were classified as sDPS events.
Four-jet DPS events were modelled by overlaying two different events. To reduce any dependence of the measurement on the modelling of jet production, this construction used events from data rather than MC simulation. Complete-DPS events were built using dijet events from the A and B samples selected from data (see Eq. (11)). To build sDPS events, two other samples were selected with the following requirements: The overlay was performed at the reconstructed jet level. When constructing cDPS and sDPS events the following conditions were imposed for a given pair of events to be overlaid: • none of the four jets contains the axis of one of the other jets, i.e., ∆R jet−jet > 0.6; • the vertices of the two overlaid events are no more than 10 mm apart in the z direction; • when building cDPS events, each of the overlaid events contributes two jets to the four leading jets in the constructed event; • when building sDPS events, one of the overlaid events contributes three jets to the four leading jets in the constructed event and the other contributes one jet.
The first condition ensures that none of the jets would be merged if the four-jet event had been reconstructed as a real event; the second condition avoids possible kinematic bias due to events where two jet pairs originate from far-away vertices; the last two conditions enforce the appropriate composition of the four leading jets in the constructed event.
As is discussed in Section 6.4, the topology of cDPS and sDPS events constructed by overlaying two events is compared to the topology of cDPS and sDPS events extracted from the AHJ sample respectively.

Kinematic characteristics of event classes
In cDPS, double dijet production should result in pairwise p T -balanced jets with a distance |φ 1 − φ 2 | ≈ π between the jets in each pair. In addition, the azimuthal angle between the two planes of interactions is expected to have a uniform random distribution. In SPS, the pairwise p T balancing of jets is not as likely; therefore the topology of the four jets is expected to be different for cDPS and SPS.
The topology of three of the jets in sDPS events would resemble the topology of the jets in SPS interactions. The fourth jet initiated by the primary interaction in an SPS is expected to be closer, in the φ-y plane, to the other three jets originating from that interaction. In an sDPS event, the jet produced in the secondary interaction would be emitted in a random direction relative to the other three jets.
In constructing possible differentiating variables, three guiding principles were followed: 1. use pairwise relations that have the potential to differentiate between SPS and cDPS topologies; 2. include angular relations between all jets in light of the expected topology of sDPS events; 3. attempt to construct variables least sensitive to systematic uncertainties.
The first two guidelines encapsulate the different characteristics of SPS and DPS events. The third guideline led to the usage of ratios of p T in order to avoid large dependencies on the jet energy scale (JES) uncertainty. Various studies, including the use of a principal component analysis [69], led to the following list of candidate variables for distinguishing event topologies: where p i T , p i T , y i and φ i stand for the scalar and vectorial transverse momentum, the rapidity and the azimuthal angle of jet i respectively, with i = 1, 2, 3, 4. The variables with the subscript i j are calculated for all possible jet combinations. The term φ i+ j denotes the azimuthal angle of the four-vector obtained by the sum of jets i and j.
In the following, the pairing notation { i, j k, l } is used to describe a cDPS event in which jets i and j originate from one interaction and jets k and l originate from the other. In around 85% of cDPS events, the two leading jets originate from one interaction and jets 3 and 4 originate from the other. indicating that both the leading and the sub-leading jet pairs are balanced in p T . The small peak around unity is due to events in which the appropriate pairing of the jets is { 1, 3 2, 4 } or { 1, 4 2, 3 }. In the SPS and sDPS samples, the leading jet-pair exhibits a wider peak at higher values of ∆ p T 12 compared to that in the cDPS sample. This indicates that the two leading jets are not well balanced in p T since a significant fraction of the hard-scatter momentum is carried by additional jets.
The ∆φ 34 distributions in the three samples are shown in Figure 2(c). The p T balance between the jets seen in the ∆  and ∆φ 34 variables can be readily understood through the following approximation: The peak around unity observed in the ∆ p T 34 distributions in the SPS and sDPS samples is thus a direct consequence of the Jacobian of the relation between ∆ p T 34 and ∆φ 34 . The set of variables quantifying the distance between jets in rapidity, ∆y i j , is particularly important for the sDPS topology. The colour flow is different in SPS leading to the four-jet final state and results in smaller angles between the sub-leading jets. Hence, on average, smaller distances between non-leading jets are expected in the SPS sample compared to the sDPS sample. This is observed in the comparison of the ∆y 34 distributions shown in Figure 2(d), where the distribution in the sDPS sample is slightly wider than in the other two samples.
The study of the various distributions in the three samples is summed up as follows: • Strong correlations between all variables are observed. The ∆ p T i j and ∆φ i j variables are correlated in a non-linear way, while geometrical constraints correlate the ∆y i j and ∆φ i j variables. Transverse momentum conservation correlates the φ i+ j − φ k+l variables with the ∆ p T i j and ∆φ i j variables.
• None of the variables displays a clear separation between all three samples. The variables in which a large difference is observed between the SPS and cDPS distributions, e.g., ∆ p T 34 , do not provide any differentiating power between SPS and sDPS.
• All variables are important -in cDPS events, where the pairing of the jets is different from { 1, 2 3, 4 }, variables relating the other possible pairs, e.g., ∆φ 13 , may indicate which is the correct pairing.
These conclusions led to the decision to use a multivariate technique in the form of an NN to perform event classification.

Extraction of the fraction of DPS events using an artificial neural network
For the purpose of training the NN, events from each sample were divided into two statistically independent subsamples, the training sample and the test sample. The former was used to train the NN and the latter to test the robustness of the result. To avoid bias during training, the events in the SPS, cDPS and sDPS training samples were reweighted such that each sample contributed a third of the total sum of weights. In all subsequent figures, only the test samples are shown.
The NN used is a feed-forward multilayer perceptron with two hidden layers, implemented in the ROOT analysis framework [70]. The input layer has 21 neurons, corresponding to the variables defined in Eq. (16), and the first and second hidden layers have 42 and 12 neurons respectively. These choices represent the product of a study conducted to optimize the performance of the NN and balance the complexity of the network with the computation time of the training. The output of the NN consists of three variables, which are interpreted as the probability for an event to be more like SPS (ξ SPS ), cDPS (ξ cDPS ) or sDPS (ξ sDPS ). During training, each event is marked as belonging to one of the samples; e.g., an event from the cDPS sample is marked as For each event, the three outputs are plotted as a single point inside an equilateral triangle (ternary plot) using the constraint ξ SPS + ξ cDPS + ξ sDPS = 1. A point in the triangle expresses the three probabilities as three distances from each of the sides of the triangle. The vertices would therefore be populated by events with high probability to belong to a single sample. Figure 3 shows an illustration of the ternary plot, where the horizontal axis corresponds to 1   Figure 4(a). However, a ridge of SPS events extending towards the sDPS corner is observed as well. A contribution from SPS events is also visible in the bottom right corner. The clearest peak is seen for events from the cDPS sample in the bottom right corner in Figure 4(b). A visible cluster of sDPS events is seen in Figure 4(c) concentrated around ξ sDPS ∼ 0.75 and there is a tail of events along the side connecting the SPS and sDPS corners. The NN output distribution in the data, shown in Figure 4(d), is visually consistent with a superposition of the three components, SPS, cDPS and sDPS.
Based on these observations, it is clear that event classification on an event-by-event basis is not possible. However, the differences between the SPS, cDPS and sDPS distributions suggest that an estimation of the different contributions can be performed. To estimate the cDPS and sDPS fractions in four-jet events, the ternary distribution in data (D) is fitted to a weighted sum of the ternary distributions in the SPS (M SPS ), cDPS (M cDPS ) and sDPS (M sDPS ) samples, each normalized to the measured four-jet cross-section in data, with the fractions as free parameters. The optimal fractions were obtained using a fit of the form, where a χ 2 minimization was performed, as implemented in the Minuit [71] package in ROOT, taking into account the statistical uncertainties of all the samples in each bin. The results of the fit are presented in Section 8, after the methodology validation and discussion of systematic uncertainties.

Methodology validation
A sizeable discrepancy was found in the ∆ p T 34 and ∆φ 34 distributions between the data and AHJ (See Section 9 for details), suggesting that there are more sub-leading jets (jets 3 and 4) that are back-to-back in AHJ than in the data. In order to test that the discrepancies are not from mis-modelling of SPS in AHJ, the ∆ p T 34 and ∆φ 34 distributions in the SPS sample extracted from AHJ were compared to the distributions in the SPS sample generated in Sherpa. Good agreement between the shapes of the distributions was observed for both variables. This and further studies indicate that the excess of events with jets 3 and 4 in the back-to-back topology is due to an excess of DPS events in the AHJ sample compared to the data.
In order to verify that the topologies of cDPS and sDPS events can be reproduced by overlaying two events, the overlay samples are compared to the cDPS and sDPS samples extracted from AHJ. An extensive comparison between the distributions of the variables used as input to the NN in the overlay samples and in AHJ was performed and good agreement was observed. This can be summarized by comparing the NN output distributions. The NN is applied to the cDPS and sDPS samples extracted from AHJ and the output distributions are compared to the output distributions in the corresponding samples constructed by overlaying events selected from data. Normalized distributions of the projection of the full ternary plot on the horizontal axis are shown in Figures 5(a) and 5(b) for the cDPS and sDPS samples respectively. Good agreement is observed between the distributions. Based on these results, it is concluded that the topology of the four jets in the overlaid events is comparable to that of the four leading jets in DPS events extracted from AHJ. The added advantage of using overlaid events from data to construct the DPS samples is that the jets are at the same JES as the jets in four-jet events in data, leading to a smaller systematic uncertainty in the final result. As an additional validation step, the NN is applied to the inclusive AHJ sample and the resulting distribution is fitted with the NN output distributions of the SPS, cDPS and sDPS samples. The fraction obtained from the fit, f (MC) DPS , is compared to the fraction at parton level, f (P) DPS , extracted from the event record, Fair agreement is observed between the value obtained from the fit and that at parton level. The larger statistical uncertainty in f (MC) DPS compared to f (P) DPS reflects the loss of statistical power due to the use of a template fit to estimate the former.

Systematic uncertainties
For jets with 20 ≤ p T < 30 GeV, the fractional JES uncertainty is about 4.5% in the central region of the detector, rising to about 10% in the forward region [64]. The overall impact of the JES on the distributions, f DPS and α 4j 2j was estimated by shifting the jet energy upwards and downwards in the MC samples by the JES uncertainty and repeating the analysis. Similarly, the overall impact of the jet energy and angular resolution was determined by varying the jet energy and angular resolution in the MC samples by the corresponding resolution uncertainty [72].
The systematic uncertainties in the measured cross-sections due to the integrated luminosity measurement uncertainty (±3.5%), the jet reconstruction efficiency uncertainty (±2%) and the uncertainty as a result of selecting single-vertex events (±0.5%) were propagated to the uncertainty in σ eff .
The statistical uncertainty in the AHJ sample was translated to a systematic uncertainty in f DPS by varying the reweighting function used to reweight AHJ and repeating the analysis.
The stability of the value of σ eff relative to the various parameter values used in the measurement was studied. Parameters such as p parton T and ∆R jet−jet were varied and the requirement ∆R parton−jet ≤ 0.6 was applied, leading to a relative change in σ eff of the order of a few percent. Since the observed relative changes are small compared to the statistical uncertainty in σ eff , no systematic uncertainty was assigned due to these parameters.
The relative systematic uncertainties in f DPS , α 4j 2j and σ eff are summarized in Table 1. The dominant systematic uncertainty on f DPS originates from the JES variation. A variation in the JES results in a modification of the NN output distribution for the SPS template used in the fit, which directly impacts the value of f DPS .

Determination of σ eff
To determine f DPS and σ eff and their statistical uncertainties taking into account all of the correlations, many replica fits were performed by random sampling from the NN output distributions. uncertainties were obtained by propagating the expected variations into this analysis, and the resulting shifts were added in quadrature. The result for f DPS is where the contribution of f sDPS to f DPS was found to be about 40%. The fraction of DPS estimated in data is 65 +23 −27 % of the fraction in AHJ as extracted from the event record (see Eq. (20)). Taking into account the systematic uncertainties in the calculation of the goodness-of-fit χ 2 , a value for χ 2 /N DF of 112/84 = 1.3 is obtained, where N DF is the number of degrees of freedom in the fit.
In order to visualize the results of the fit, the ternary distribution is divided into five slices, A comparison of the fit distributions with the distributions in data in the five slices of the ternary plot is shown in Figure 6. Considering the systematic uncertainties, the most significant difference between the data and the fit is seen for the two left-most bins in the range 0.0 ≤ ξ sDPS < 0.1 (Figure 6(a)) of the ternary plot. These bins are dominated by the SPS contribution. Thus, a discrepancy between the data and the fit result in these bins is expected to have a negligible effect on the measurement of the DPS rate. A discrepancy between the data and the fit result is also observed in the three rightmost bins in Figure 6(a). These bins have about a 30% contribution from cDPS. To test the effect of this discrepancy on the description of observables in data, the distributions of the various variables in data were compared to a combination of the distributions in the SPS, cDPS and sDPS samples, normalizing the latter three distributions to their respective fractions in the data as obtained in the fit. This comparison for the ∆ p T 34 and ∆φ 34 variables is shown in Figure 7, where a good description of the data is observed. The same level of agreement is seen for all the variables.
Before calculating σ eff , the symmetry factor in Eq. (6) has to be adjusted because there is an overlap in the cross-sections σ A 2j and σ B 2j when the leading jet in sample A has p T ≥ 42.5 GeV (see Eq. (11)). The  adjusted symmetry factor is as determined from the measured dijet cross-sections. This factor was also determined using Pythia6 and good agreement was observed between the two values. The relative difference in the value of σ eff obtained by using the symmetry factors extracted from the data and from Pythia6 was of the order of 0.2%, a negligible difference compared to the statistical uncertainty of σ eff .
An additional correction of +4% is applied to the measured DPS cross-section due to the probability of jets from the secondary interaction overlapping with jets from the primary interaction. In this configuration, the anti-k t algorithm merges the two overlapping jets into one, and hence the event cannot pass the four-jet requirement. The value of this correction was determined from the fraction of phase space occupied by a jet. It was also determined directly in AHJ and good agreement between the two values was observed.
Finally, the measurements of the dijet and four-jet cross-sections can be used to calculate the effective cross-section, yielding This value is consistent within the quoted uncertainties with previous measurements, performed by the ATLAS collaboration and by other experiments [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30], all of which are summarized in Figure 8. Figure 9 shows σ eff as a function of √ s, where the AFS result and some of the LHCb results are omitted for clarity. Within the large uncertainties, the measurements are consistent with no √ s dependence of σ eff . The σ eff value obtained is 21 +7 −6 % of the inelastic cross-section, σ inel , measured by ATLAS at √ s = 7 TeV [73].

Normalized differential cross-sections
To allow the results of this study to be used in future comparisons with MPI models, the distributions of the variables used as input to the NN were corrected for detector effects. The corrections were derived using an iterative unfolding, producing an unfolding matrix for each observable, relating the particle-level and reconstructed-level quantities. These matrices were derived using samples of four-jet events selected from the AHJ and Pythia6 samples by imposing the cuts detailed in Eq. (11) on particle jets. The AHJ sample generated with the AUET1 tune was used to derive the unfolding matrix. The distributions were unfolded with the Bayesian unfolding algorithm, implemented in the RooUnfold package [74], using two iterations.
The unfolding matrices derived from AHJ were taken as the nominal matrices and the differences observed when using the matrices derived from Pythia6 were used as an additional systematic uncertainty, typically of the order of a few percent in each bin. The total systematic uncertainty of the differential distributions in data was obtained by summing in quadrature the uncertainty due to MC modelling in a given bin with the systematic uncertainties in this bin due to the JES and jet energy and angular resolution uncertainties, while preserving correlations between bins. Figure 10 shows the normalized differential cross-section distribution in data for the ∆ p T 34 and ∆φ 34 variables compared to the particle-level distributions in the AHJ samples generated with the AUET1 and AUET2 tunes. The particle-level distributions in the AUET2 AHJ sample overestimate the normalized differential cross-section distributions in data in the regions ∆ p T 34 ≤ 0.15 and ∆φ 34 ≥ 2.8, demonstrating the excess of the DPS contribution in this sample compared to the data. On the other hand, the DPS contribution in the data is underestimated by the prediction obtained with the AUET1 tune. These comparisons demonstrate the power of these distributions to constrain MPI models and tunes. In Appendix A, the normalized differential cross-sections in data for the remaining variables are compared to the particle-level distributions in the AHJ samples generated using the AUET1 and AUET2 tunes.

Summary and conclusions
A measurement of the rate of hard double-parton scattering in four-jet events was performed using a sample of data collected with the ATLAS experiment at the LHC in 2010, with an average of approximately 0.4 proton-proton interactions per bunch crossing, corresponding to an integrated luminosity of 37.3 ± 1.3 pb −1 . Three different samples were selected, all consisting of single-vertex events from proton-proton collisions at a centre-of-mass energy of √ s = 7 TeV. Four-jet events were defined as The contribution of hard double-parton scattering to the production of four-jet events was extracted using an artificial neural network. The four-jet topology originating from hard double-parton scattering was represented by a random combination of events selected in data. The fraction of events corresponding to the contribution made by hard double-parton scattering in four-jet events was determined to be, After combining this result with measurements of the dijet and four-jet cross-sections in the appropriate phase space regions, the effective cross-section was determined to be This value is 21 +7 −6 % of the measured value of σ inel at √ s = 7 TeV and is consistent with previous measurements performed at various centre-of-mass energies and in various final states. It is compatible with a model in which σ eff is a universal parameter that does not depend on the process or phase space. To facilitate future studies of the dynamics of multi-parton interactions, distributions of observables sensitive to the presence of hard double-parton scattering are also presented.  16), in data after unfolding to particle level, compared to the MC prediction from AHJ at the particle level, generated using the AUET1 and AUET2 tunes, as indicated in the legend. The hatched area represents the sum in quadrature of the statistical and systematic uncertainties in the normalized differential cross-sections and all histograms are normalized to unity. The ratio of the particle-level distribution to the normalized differential cross-section is shown in the bottom panels, where the shaded areas represent statistical uncertainties. Appendix A. Normalized differential cross-sections  show the normalized differential cross-sections in data for all the observables used as input to the NN, compared to the particle-level distributions in the AHJ samples generated using the AUET1 and AUET2 tunes. The hatched areas in the distributions represent the total uncertainty of the normalized differential cross-section, obtained by adding in quadrature the statistical and systematic uncertainties. , in data after unfolding to particle level, compared to the MC prediction from AHJ at the particle level, generated using the AUET1 and AUET2 tunes, as indicated in the legend. The hatched areas represent the sum in quadrature of the statistical and systematic uncertainties in the normalized differential cross-sections and all histograms are normalized to unity. The ratio of the particle-level distribution to the normalized differential cross-section is shown in the bottom panels, where the shaded areas represent statistical uncertainties.  (d) Figure 12: Distributions of the variables (a) ∆ p T 24 , (b) ∆φ 12 , (c) ∆φ 13 and (d) ∆φ 23 , defined in Eq. (16), in data after unfolding to particle level, compared to the MC prediction from AHJ at the particle level, generated using the AUET1 and AUET2 tunes, as indicated in the legend. The hatched areas represent the sum in quadrature of the statistical and systematic uncertainties in the normalized differential cross-sections and all histograms are normalized to unity. The ratio of the particle-level distribution to the normalized differential cross-section is shown in the bottom panels, where the shaded areas represent statistical uncertainties.  (d) Figure 13: Distributions of the variables (a) ∆φ 14 , (b) ∆φ 24 , (c) ∆y 12 and (d) ∆y 34 , defined in Eq. (16), in data after unfolding to particle level, compared to the MC prediction from AHJ at the particle level, generated using the AUET1 and AUET2 tunes, as indicated in the legend. The hatched areas represent the sum in quadrature of the statistical and systematic uncertainties in the normalized differential cross-sections and all histograms are normalized to unity. The ratio of the particle-level distribution to the normalized differential cross-section is shown in the bottom panels, where the shaded areas represent statistical uncertainties.  (d) Figure 14: Distributions of the variables (a) ∆y 13 , (b) ∆y 23 , (c) ∆y 14 and (d) ∆y 24 defined in Eq. (16), in data after unfolding to particle level, compared to the MC prediction from AHJ at the particle level, generated using the AUET1 and AUET2 tunes, as indicated in the legend. The hatched areas represent the sum in quadrature of the statistical and systematic uncertainties in the normalized differential cross-sections and all histograms are normalized to unity. The ratio of the particle-level distribution to the normalized differential cross-section is shown in the bottom panels, where the shaded areas represent statistical uncertainties.  16), in data after unfolding to particle level, compared to the MC prediction from AHJ at the particle level, generated using the AUET1 and AUET2 tunes, as indicated in the legend. The hatched areas represent the sum in quadrature of the statistical and systematic uncertainties in the normalized differential cross-sections and all histograms are normalized to unity. The ratio of the differential distribution to the particle-level distributions is shown in the bottom panels, where the shaded areas represent statistical uncertainties.