Study of the underlying event in top quark pair production in p p collisions at 13 Te .

Measurements of normalized differential cross sections as functions of the multiplicity and kinematic variables of charged-particle tracks from the underlying event in top quark and antiquark pair production are presented. The measurements are performed in proton-proton collisions at a center-of-mass energy of 13 Te , and are based on data collected by the CMS experiment at the LHC in 2016 corresponding to an integrated luminosity of 35.9 fb - 1 . Events containing one electron, one muon, and two jets from the hadronization and fragmentation of b quarks are used. These measurements characterize, for the first time, properties of the underlying event in top quark pair production and show no deviation from the universality hypothesis at energy scales typically above twice the top quark mass.

Introduction

At the LHC, top quark and antiquark pairs () are dominantly produced in the scattering of the proton constituents via quantum chromodynamics (QCD) at an energy scale (Q) of about two times the quark mass (). The properties of the quark can be studied directly from its decay products, as it decays before hadronizing. Mediated by the electroweak interaction, the quark decay yields a boson and a quark, the latter carrying the QCD color charge of the mother particle. Given the large branching fraction for the decay into a bottom quark,  [1], in this analysis we assume that each or quark yields a corresponding bottom () or antibottom () quark in its decay. Other quarks may also be produced, in a color-singlet state, if a decay occurs. Being colored, these quarks will fragment and hadronize giving rise to an experimental signature with jets. Thus, when performing precision measurements of quark properties at hadron colliders, an accurate description of the fragmentation and hadronization of the quarks from the hard scatter process as well as of the “underlying event” (UE), defined below, is essential. First studies of the fragmentation and hadronization of the quarks in events have been reported in Refs. [2, 3]. In this paper, we present the first measurement of the properties of the UE in events at a scale .

The UE is defined as any hadronic activity that cannot be attributed to the particles stemming from the hard scatter, and in this case from decays. Because of energy-momentum conservation, the UE constitutes the recoil against the system. In this study, the hadronization products of initial- and final-state radiation (ISR and FSR) that cannot be associated to the particles from the decays are probed as part of the UE, even if they can be partially modeled by perturbative QCD. The main contribution to the UE comes from the color exchanges between the beam particles and is modeled in terms of multiparton interactions (MPI), color reconnection (CR), and beam-beam remnants (BBR), whose model parameters can be tuned to minimum bias and Drell–Yan (DY) data.

The study of the UE in events provides a direct test of its universality at higher energy scales than those probed in minimum bias or DY events. This is relevant as a direct probe of CR, which is needed to confine the initial QCD color charge of the quark into color-neutral states. The CR mainly occurs between one of the products of the fragmentation of the quark from the quark decay and the proton remnants. This is expected to induce an ambiguity in the origin of some of the final states present in a bottom quark jet [4-6]. The impact of these ambiguities in the measurement of quark properties is evaluated through phenomenological models that need to be tuned to the data. Recent examples of the impact that different model parameters have on can be found in Refs. [7, 8].

The analysis is performed using final states where both of the bosons decay to leptons, yielding one electron and one muon with opposite charge sign, and the corresponding neutrinos. In addition, two jets are required in the selection, as expected from the decay. This final state is chosen because of its expected high purity and because the products of the hard process can be distinguished with high efficiency and small contamination from objects not associated with quark decays, e.g., jets from ISR.

After discussing the experimental setup in Sect. 2, and the signal and background modeling in Sect. 3, we present the strategy employed to select the events in Sect. 4 and to measure the UE contribution in each selected event in Sect. 5. The measurements are corrected to a particle-level definition using the method described in Sect. 6 and the associated systematic uncertainties are discussed in Sect. 7. Finally, in Sect. 8, the results are discussed and compared to predictions from different Monte Carlo (MC) simulations. The measurements are summarized in Sect. 9.

The CMS detector

The central feature of the CMS apparatus is a superconducting solenoid of 6 internal diameter, providing a magnetic field of 3.8 parallel to the beam direction.

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. A preshower detector, consisting of two planes of silicon sensors interleaved with about three radiation lengths of lead, is located in front of the endcap regions of the ECAL. Hadron forward calorimeters, using steel as an absorber and quartz fibers as the sensitive material, extend the pseudorapidity coverage provided by the barrel and endcap detectors from to 5.2. Muons are detected in the window in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid.

Charged-particle trajectories with are measured by the tracker system. The particle-flow algorithm [9] is used to reconstruct and identify individual particles in an event, with an optimized combination of information from the various elements of the CMS detector. The energy of the photons is directly obtained from the ECAL measurement, corrected for zero-suppression effects. The energy of the electrons is determined from a combination of the electron momentum at the primary interaction vertex as determined by the tracker, the energy of the corresponding ECAL cluster, and the energy sum of all bremsstrahlung photons spatially compatible with originating from the electron track. The energy of the muons is obtained from the curvature of the corresponding track. The energy of charged hadrons is determined from a combination of their momentum measured in the tracker and the matching ECAL and HCAL energy deposits, corrected for zero-suppression effects and for the response function of the calorimeters to hadronic showers. Finally, the energy of neutral hadrons is obtained from the corresponding corrected ECAL and HCAL energies.

Events of interest are selected using a two-tiered trigger system [10]. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 within a time interval of less than 4. The second level, known as the high-level trigger, 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 before data storage.

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

Signal and background modeling

This analysis is based on proton-proton ( ) collision data at a center-of-mass energy , collected by the CMS detector in 2016 and corresponds to an integrated luminosity of  [12].

The process is simulated with the powheg (v2) generator in the heavy quark production (hvq) mode [13-15]. The NNPDF3.0 next-to-leading-order (NLO) parton distribution functions (PDFs) with the strong coupling parameter at the boson mass scale () [16] are utilized in the matrix-element (ME) calculation. The renormalization and factorization scales, and , are set to , where and is the transverse momentum in the rest frame. Parton showering is simulated using pythia8 (v8.219) [17] and the CUETP8M2T4 UE tune [18]. The CUETP8M2T4 tune is based on the CUETP8M1 tune [19] but uses a lower value of in the parton shower (PS); this value leads to a better description of jet multiplicities in events at  [20]. The leading-order (LO) version of the same NNPDF3.0 is used in the PS and MPI simulation in the CUETP8M2T4 tune. The cross section used for the simulation is , computed at the next-to-next-to-leading-order (NNLO) plus next-to-next-to-leading-logarithmic accuracy  [21].

Throughout this paper, data are compared to the predictions of different generator settings for the process. Table 1 summarizes the main characteristics of the setups and abbreviations used in the paper. Among other UE properties, CR and MPI are modeled differently in the alternative setups considered, hence the interest in comparing them to the data. Three different signal ME generators are used: powheg, MadGraph 5_amc@nlo (v2.2.2) with the FxFx merging scheme [22, 23] for jets from the ME calculations and PS, and sherpa (v2.2.4) [24]. The latter is used in combination with OpenLoops (v1.3.1) [25], and with the CS parton shower based on the Catani–Seymour dipole subtraction scheme [26]. In addition, two different herwig PS versions are used and interfaced with powheg: herwig++  [27] with the EE5C UE tune [28] and the CTEQ6 (L1) [29] PDF set, and herwig 7 [27, 30] with its default tune and the MMHT2014 (LO) [31] PDF set.

Table 1

MC simulation settings used for the comparisons with the differential cross section measurements of the UE. The table lists the main characteristics and values used for the most relevant parameters of the generators. The row labeled “Setup designation” shows the definitions of the abbreviations used throughout this paper

Event generator	powheg (v2)	MadGraph 5_amc@nlo  (v2.2.2)	sherpa (v2.2.4)
Matrix element characteristics
Mode	hvq	FxFx Merging	OpenLoops
Scales (μ_R,μ_F)	m_T	∑_t,t¯m_T/2	METS+QCD
α_S(M_Z)	0.118	0.118	0.118
PDF	NNPDF3.0 NLO	NNPDF3.0 NLO	NNPDF3.0 NNLO
Accuracy	tt¯ [NLO]	tt¯+0,1,2jets [NLO]	tt¯ [NLO]
	1 jet [LO]	3 jets [LO]
Parton shower
Setup designation	Pw+Py8	MG5_aMC	sherpa
PS	pythia (v8.219)		CS
Tune	CUETP8M2T4		default
PDF	NNPDF2.3 LO		NNPDF3.0 NNLO
(αSISR(M_Z),αSFSR(M_Z))	(0.1108, 0.1365)		(0.118, 0.118)
ME corrections	on		—
Setup designation	Pw+Hw++	Pw+Hw7
PS	herwig++	herwig 7
Tune	EE5C	Default
PDF	CTEQ6 (L1)	MMHT2014 LO
(αSISR(M_Z),αSFSR(M_Z))	(0.1262, 0.1262)	(0.1262, 0.1262)
ME corrections	Off	On

Additional variations of the Pw+Py8 sample are used to illustrate the sensitivity of the measurements to different parameters of the UE model. A supplementary table, presented in the appendix, details the parameters that have been changed with respect to the CUETP8M2T4 tune in these additional variations. The variations include extreme models that highlight separately the contributions of MPI and CR to the UE, fine-grained variations of different CR models [5, 32], an alternative MPI model based on the “Rope hadronization” framework describing Lund color strings overlapping in the same area [33, 34], variations of the choice of in the parton shower, and a variation of the values that constitute the CUETP8M2T4 tune, according to their uncertainties.

Background processes are simulated with several generators. The , +jets, and (where denotes any of the charged leptons e//) processes are simulated at NLO, using MadGraph 5_amc@nlo with the FxFx merging. Drell–Yan production, with dilepton invariant mass, , greater than 50, is simulated at LO with MadGraph 5_amc@nlo using the so-called MLM matching scheme [35] for jet merging. The powheg (v2) program is furthermore used to simulate the , and processes [36, 37], while powheg (v1) is used to simulate the process [38]. The single top quark t-channel background is simulated at NLO using powheg (v2) and MadSpin contained in MadGraph 5_amc@nlo (v2.2.2) [39, 40]. The residual +V backgrounds, where or , are generated at NLO using MadGraph 5_amc@nlo. The cross sections of the DY and +jets processes are normalized to the NNLO prediction, computed using fewz (v3.1.b2) [41], and single top quark processes are normalized to the approximate NNLO prediction [42]. Processes containing two vector bosons (hereafter referred to as dibosons) are normalized to the NLO predictions computed with MadGraph 5_amc@nlo, with the exception of the process, for which the NNLO prediction [43] is used.

All generated events are processed through the Geant4-based [44-46] CMS detector simulation and the standard CMS event reconstruction. Additional collisions per bunch crossing (pileup) are included in the simulations. These simulate the effect of pileup in the events, with the same multiplicity distribution as that observed in data, i.e., about 23 simultaneous interactions, on average, per bunch crossing.

Event reconstruction and selection

The selection targets events in which each boson decays to a charged lepton and a neutrino. Data are selected online with single-lepton and dilepton triggers. The particle flow (PF) algorithm [9] is used for the reconstruction of final-state objects. The offline event selection is similar to the one described in Ref. [47]. At least one PF charged lepton candidate with and another one with , both having , are required. The two leptons must have opposite charges and an invariant mass . When extra leptons are present in the event, the dilepton candidate is built from the highest leptons in the event. Events with in the final state are used for the main analysis, while and events are used to derive the normalization of the DY background. The simulated events are corrected for the differences between data and simulation in the efficiencies of the trigger, lepton identification, and lepton isolation criteria. The corrections are derived with and events using the “tag-and-probe” method [48] and are parameterized as functions of the and of the leptons.

Jets are clustered using the anti- jet finding algorithm [49, 50] with a distance parameter of 0.4 and all the reconstructed PF candidates in the event. The charged hadron subtraction algorithm is used to mitigate the contribution from pileup to the jets [51]. At least two jets with , and identified by a -tagging algorithm are required. The -tagging is based on a “combined secondary vertex” algorithm  [52] characterized by an efficiency of about 66%, corresponding to misidentification probabilities for light quark and quark jets of 1.5 and 18%, respectively. A -dependent scale factor is applied to the simulations in order to reproduce the efficiency of this algorithm, as measured in data.

The reconstructed vertex with the largest value of summed physics-object is taken to be the primary interaction vertex. The physics objects are the jets, clustered using the jet finding algorithm [49, 50] with the tracks assigned to the vertex as inputs, and the associated missing transverse momentum, , taken as the negative vector sum of the of those jets. The latter is defined as the magnitude of the negative vector sum of the momenta of all reconstructed PF candidates in an event, projected onto the plane perpendicular to the direction of the proton beams.

All backgrounds are estimated from simulation, with the exception of the DY background normalization. The latter is estimated making use of the so-called method [53], in which events with same-flavor leptons are used to normalize the yield of pairs from DY production of lepton pairs. The normalization of the simulation is estimated from the number of events in the data within a 15 window around the boson mass [53]. For events, we use the geometric mean of the scale factors determined for and events. With respect to the simulated predictions, a scale factor is obtained from this method, with statistical and systematic uncertainties added in quadrature. The systematic uncertainty is estimated from the differences found in the scale factor for events with 0 or 1 -tagged jets, in the same-flavor channels.

We select a total of 52 645 events with an expected purity of 96%. The data agree with the expected yields within 2.2%, a value smaller than the uncertainty in the integrated luminosity alone, 2.5% [12]. The events are expected to constitute 90% of the total background.

In the simulation, the selection is mimicked at the particle level with the techniques described in Ref. [54]. Jets and leptons are defined at the particle level with the same conventions as adopted by the rivet framework [55]. The procedure ensures that the selections and definitions of the objects at particle level are consistent with those used in the rivet routines. A brief description of the particle-level definitions follows:Using these definitions, the fiducial region of this analysis is specified by the same requirements that are applied offline (reconstruction level) for leptons and jets. Simulated events are categorized as fulfilling only the reconstruction-based, only the particle-based, or both selection requirements. If a simulated event passes only the reconstruction-level selection, it is considered in the “misidentified signal” category, i.e., it does not contribute to the fiducial region defined in the analysis and thus is considered as a background process. In the majority of the bins of each of the distributions analyzed, the fraction of signal events passing both the reconstruction- and particle-level selections is estimated to be about 80%, while the fraction of misidentified signal events is estimated to be less than 10%.

prompt charged leptons (i.e., not produced as a result of hadron decays) are reconstructed as “dressed” leptons with nearby photon recombination using the anti- algorithm with a distance parameter of 0.1;

jets are clustered with the anti- algorithm with a distance parameter of 0.4 using all particles remaining after removing both the leptons from the hard process and the neutrinos;

the flavor of a jet is identified by including hadrons in the clustering.

Characterization of the underlying event

In order to isolate the UE activity in data, the contribution from both pileup and the hard process itself must be identified and excluded from the analysis. The contamination from pileup events is expected to yield soft particles in time with the hard process, as well as tails in the energy deposits from out-of-time interactions. The contamination from the hard process is expected to be associated with the two charged leptons and two jets originating from the decay chain.

In order to minimize the contribution from these sources, we use the properties of the reconstructed PF candidates in each event. The track associated to the charged PF candidate is required to be compatible with originating from the primary vertex. This condition reduces to a negligible amount the contamination from pileup in the charged particle collection. A simple association by proximity in z with respect to the primary vertex of the event is expected to yield a pileup-robust, high-purity selection. For the purpose of this analysis all charged PF candidates are required to satisfy the following requirements:After performing the selection of the charged PF candidates we check which ones have been used in the clustering of the two -tagged jets and which ones match the two charged lepton candidates within a cone, where is the azimuthal angle in radians. All PF candidates failing the kinematic requirements, being matched to another primary vertex in the event, or being matched to the charged leptons and -tagged jets, are removed from the analysis. The UE analysis proceeds by using the remaining charged PF candidates. Figure 1 shows, in a simulated event, the contribution from charged and neutral PF candidates, the charged component of the pileup, and the hard process. The charged PF candidates that are used in the study of the UE are represented after applying the selection described above.

Fig. 1

Distribution of all PF candidates reconstructed in a Pw+Py8 simulated event in the – plane. Only particles with are shown, with a marker whose area is proportional to the particle . The fiducial region in is represented by the dashed lines

and ;

the associated track needs to be either used in the fit of the primary vertex or to be closer to it in z than with respect to other reconstructed vertices in the event.

Various characteristics, such as the multiplicity of the selected charged particles, the flux of momentum, and the topology or shape of the event have different sensitivity to the modeling of the recoil, the contribution from MPI and CR, and other parameters.

The first set of observables chosen in this analysis is related to the multiplicity and momentum flux in the event:The second set of observables characterizes the UE shape and it is computed from the so-called linearized sphericity tensor [56, 57]:where the i index runs over the particles associated with the UE, as for the previous variables, and the and indices refer to one of the (x, y, z) components of the momentum of the particles. The eigenvalues () of are in decreasing order, i.e., with the largest one, and are used to compute the following observables [58]:Further insight can be gained by studying the evolution of the two sets of observables in different categories of the system kinematic quantities. The categories chosen below are sensitive to the recoil or the scale of the energy of the hard process, and are expected to be reconstructed with very good resolution. Additionally, these variables minimize the effect of migration of events between categories due to resolution effects.

charged-particle multiplicity: ;

magnitude of the of the charged particle recoil system: ;

scalar sum of the (or ) of charged particles: , where or z;

average (or ) per charged particle: computed from the ratio between the scalar sum and the charged multiplicity: (or ).

Aplanarity: measures the component out of the event plane, defined by the two leading eigenvectors. Isotropic (planar) events are expected to have .

Sphericity: measures the with respect to the axis of the event. An isotropic (dijet) event is expected to have .

identifies 3 jet events (tends to be 0 for dijet events).

identifies 4 jet events (tends to be 0 otherwise).

The dependence of the UE on the recoil system is studied in categories that are defined according to the multiplicity of additional jets with and , excluding the two selected -tagged jets. The categories with 0, 1, or more than 1 additional jet are used for this purpose. The additional jet multiplicity is partially correlated with the charged-particle multiplicity and helps to factorize the contribution from ISR. The distribution of the number of additional jets is shown in Fig. 2 (upper left).

Fig. 2

Distributions of the variables used to categorize the study of the UE. Upper left: multiplicity of additional jets (). Upper right: . Lower: . The distributions in data are compared to the sum of the expectations for the signal and backgrounds. The shaded band represents the uncertainty associated to the integrated luminosity and the theoretical value of the cross section

In addition to these categories, the transverse momentum of the dilepton system, , is used as it preserves some correlation with the transverse momentum of the system and, consequently, with the recoil of the system. The direction is used to define three regions in the transverse plane of each event. The regions are denoted as “transverse” (), “away” (), and “toward” (). Each reconstructed particle in an event is assigned to one of these regions, depending on the difference of their azimuthal angle with respect to the vector. Figure 3 illustrates how this classification is performed on a typical event. This classification is expected to enhance the sensitivity of the measurements to the contributions from ISR, MPI and CR in different regions. In addition, the magnitude, , is used to categorize the events and its distribution is shown in Fig. 2 (upper right). The variable is estimated with a resolution better than 3%.

Fig. 3

Display of the transverse momentum of the selected charged particles, the two leptons, and the dilepton pair in the transverse plane corresponding to the same event as in Fig. 1. The of the particles is proportional to the length of the arrows and the dashed lines represent the regions that are defined relative to the direction. For clarity, the of the leptons has been rescaled by a factor of 0.5

Lastly, the dependence of the UE on the energy scale of the hard process is characterized by measuring it in different categories of the variable. This variable is correlated with the invariant mass of the system, but not with its . The  distribution is shown in Fig. 2 (lower). A resolution better than 2% is expected in the measurement of .

Although both and are only partially correlated with the kinematic quantities, they are expected to be reconstructed with very good resolution. Because of the two escaping neutrinos, the kinematics of the pair can only be reconstructed by using the measurement, which has poorer experimental resolution when compared to the leptons. In addition, given that is correlated with the UE activity, as it stems from the balance of all PF candidates in the transverse plane, it could introduce a bias in the definition of the categories and the observables studied to characterize the UE. Hence the choice to use only dilepton-related variables.

Corrections to the particle level

Inefficiencies of the track reconstruction due to the residual contamination from pileup, nuclear interactions in the tracker material, and accidental splittings of the primary vertex [59] are expected to cause a slight bias in the observables described above. The correction for these biases is estimated from simulation and applied to the data by means of an unfolding procedure.

At particle (generator) level, the distributions of the observables of interest are binned according to the resolutions expected from simulation. Furthermore, we require that each bin contains at least 2% of the total number of events. The migration matrix (K), used to map the reconstruction- to particle-level distributions, is constructed using twice the number of bins at the reconstruction level than the ones used at particle level. This procedure ensures almost diagonal matrices, which have a numerically stable inverse. The matrix is extended with an additional row that is used to count the events failing the reconstruction-level requirements, but found in the fiducial region of the analysis, i.e., passing the particle-level requirements. The inversion of the migration matrix is made using a Tikhonov regularization procedure [60], as implemented in the TUnfoldDensity package [61]. The unfolded distribution is found by minimizing a functionwhere y are the observations, is an estimate of the covariance of y (calculated using the simulated signal sample), is the particle-level expectation, is a penalty function (with being estimated from the simulated samples), and is the so-called regularization parameter. The latter regulates how strongly the penalty term should contribute to the minimization of . In our setup we choose the function L to be the curvature, i.e., the second derivative, of the output distribution. The chosen value of the parameter is optimized for each distribution by minimizing its average global correlation coefficient [61]. Small values, i.e., , are found for all the distributions; the global correlation coefficients are around 50%. After unfolding, the distributions are normalized to unity.

The statistical coverage of the unfolding procedure is checked by means of pseudo-experiments based on independent Pw+Py8 samples. The pull of each bin in each distribution is found to be consistent with that of a standard normal distribution. The effect of the regularization term in the unfolding is checked in the data by folding the measured distributions and comparing the outcome to the originally-reconstructed data. In general the folded and the original distributions agree within 1–5% in each bin, with the largest differences observed in bins with low yield.

Systematic uncertainties

The impact of different sources of uncertainty is evaluated by unfolding the data with alternative migration matrices, which are obtained after changing the settings in the simulations as explained below. The effect of a source of uncertainty in non-fiducial events is included in this estimate, by updating the background prediction. The observed bin-by-bin differences are used as estimates of the uncertainty. The impact of the uncertainty in the background normalization is the only exception to this procedure, as detailed below. The covariance matrices associated to each source of uncertainty are built using the procedure described in detail in [62]. In case several sub-contributions are used to estimate a source of uncertainty, the corresponding differences in each bin are treated independently, symmetrized, and used to compute individual covariance matrices, which preserve the normalization. Variations on the event yields are fully absorbed by normalizing the measured cross sections. Thus, only the sources of uncertainty that yield variations in the shapes have a non-negligible impact.

Experimental uncertainties

The following experimental sources of uncertainty are considered:

Although pileup is included in the simulation, there is an intrinsic uncertainty in modeling its multiplicity. An uncertainty of in the inelastic cross section is used and propagated to the event weights [63].

The scale factors used to correct the simulation for different trigger and lepton selection efficiencies in data and simulation are varied up or down, according to their uncertainty. The uncertainties in the muon track and electron reconstruction efficiencies are included in this category and added in quadrature.

The corrections applied to the electron energy and muon momentum scales are varied separately, according to their uncertainties. The corrections and uncertainties are obtained using methods similar to those described in Refs. [64, 65]. These variations lead to a small migration of events between the different or categories used in the analysis.

A - and -dependent parameterization of the jet energy scale is used to vary the calibration of the jets in the simulation. The corrections and uncertainties are obtained using methods similar to those described in Ref. [51]. The effect of these variations is similar to that described for the lepton energy scale uncertainty; in this case the migration of events occurs between different jet multiplicity categories.

Each jet is further smeared up or down depending on its and , with respect to the central value measured in data. The difference with respect to data is measured using methods similar to those described in Ref. [51]. The main effect induced in the analysis from altering the jet energy resolution is similar to that described for the jet energy scale uncertainty.

The scale factors used to correct for the difference in performance between data and simulation are varied according to their uncertainties and depending on the flavor of the jet [52]. The main effect of this variation is to move jets into the candidate jets sample or remove them from it.

The impact of the uncertainty in the normalization of the backgrounds is estimated by computing the difference obtained with respect to the nominal result when these contributions are not subtracted from data. This difference is expected to cover the uncertainty in the normalization of the main backgrounds, i.e., DY and the process, and the uncertainty in the normalization of the events that are expected to pass the reconstruction-level requirements but fail the generator-level ones. The total expected background contribution is at the level of 8–10%, depending on the bin. The impact from this uncertainty is estimated to be .

The efficiency of track reconstruction is found to be more than 90%. It is monitored using muon candidates from decays, and the ratio of the four-body final decay to the two-body decay. The latter is used to determine a data-to-simulation scale factor () as a function of the pseudorapidity of the tracks, and for different periods of the data taking used in this analysis. The envelope of the values, with uncertainties included, ranges from 0.89 to 1.17 [66], and it provides an adequate coverage for the residual variations observed in the charged-particle multiplicity between different data taking periods. The impact of the variation of by its uncertainty is estimated by using the value of for the probability to remove a reconstructed track from the event or to promote an unmatched generator-level charged particle to a reconstructed track, depending on whether or , respectively. Different migration matrices, reflecting the different tracking efficiencies obtained from varying the uncertainty in , are obtained by this method and used to unfold the data. Although the impact is nonnegligible on variables such as or , it has very small impact () on variables such as and .

Theoretical uncertainties

The following theoretical uncertainties are considered:

and are varied individually in the ME by factors between 0.5 and 2, excluding the extreme cases and , according to the prescription described in Refs. [67, 68].

In powheg, the real emission cross section is scaled by a damping function, parameterized by the so-called variable [13-15]. This parameter controls the ME-PS matching and regulates the high- radiation by reducing real emissions generated by powheg with a factor of . In the simulation used to derive the migration matrices, and the uncertainty in this value is evaluated by changing it by +42 or -37%, a range that is determined from the jet multiplicity measurements in at  [20]. Likewise, the uncertainty associated with the choice of for space-like and for time-like showers in the CUETP8M2T4 tune is evaluated by varying the scale at which it is computed, , by a factor of 2 or 1/2.

The dependence of the migration matrix on the UE model assumed in the simulation is tested by varying the parameters that model the MPI and CR in the range of values corresponding to the uncertainty envelope associated to the CUETP8M2T4 tune. The uncertainty envelope has been determined using the same methods as described in Ref. [19]. In the following, these will be referred to as UE up/down variations. The dependence on the CR model is furthermore tested using other models besides the nominal one, which is the MPI-based CR model where the decay products are excluded from reconnections to the UE. A dedicated sample where the reconnections to resonant decay products are enabled (hereafter designated as ERDon) is used to evaluate possible differences in the unfolded results. In addition, alternative models for the CR are tested. One sample utilizing the “gluon move” model [5], in which gluons can be moved to another string, and another utilizing the “QCD-based” model with string formation beyond LO [32] are used for this purpose. In both samples, the reconnections to the decay resonant processes are enabled. The envelope of the differences is considered as a systematic uncertainty.

The effect of reweighting of the simulated quark () distribution to match the one reconstructed from data [69, 70] is added as an additional uncertainty. This has the most noticeable effect on the fraction of events that do not pass the reconstruction-level requirements and migrate out of the fiducial phase space.

An additional uncertainty is considered, related to the value of used in the simulations, by varying this parameter by  [71].

Any possible uncertainty from the choice of the hadronization model is expected to be significantly smaller than the theory uncertainties described above. This has been explicitly tested by comparing the results at reconstruction level and after unfolding the data with the Pw+Py8 and Pw+Hw++ migration matrices. The latter relies on a different hadronization model, but it folds other modelling differences such as the underlying event tune or the parton shower as well. Thus it can only be used as a test setup to validate the measurement.

Summary of systematic uncertainties

The uncertainties on the measurement of the normalized differential cross sections are dominated by the systematic uncertainties, although in some bins of the distributions the statistical uncertainties are a large component. The experimental uncertainties have, in general, small impact; the most relevant are the tracking reconstruction efficiency for the , , , and observables. Other observables are affected at a sub-percent level by this uncertainty. Theory uncertainties affect the measurements more significantly, a fact that underlines the need of better tuning of the model parameters.

Event shape observables are found to be the most robust against this uncertainty, while , , and are the ones that suffer more from it. Other sources of theoretical uncertainty typically have a smaller effect.

To further illustrate the impact of different sources on the observables considered, we list in Table 2 the uncertainties on the average of each observable. In the table, only systematic uncertainties that impact the average of one of the observables by at least 0.5% are included. The total uncertainty on the average of a given quantity ranges from 1 to 8%, and hence the comparison with model predictions can be carried out in a discrete manner.

Table 2

Uncertainties affecting the measurement of the average of the UE observables. The values are expressed in % and the last row reports the quadratic sum of the individual contributions

Source	% Uncertainty
	N_ch	∑p_T	∑p_z	p_T¯	p_z¯	|p→_T|	S	A	C	D
Statistical	0.1	0.2	0.3	0.2	0.2	0.3	0.1	0.1	0.1	0.1
Experimental
Background	1.2	1.6	1.8	0.4	0.7	1.6	0.4	0.7	0.3	0.7
Tracking eff.	4.4	4.2	4.9	0.8	0.4	4.0	0.4	0.6	0.2	0.6
Theory
μ_R/μ_F	0.5	0.8	1.0	0.3	0.3	1.0	0.1	0.1	0.1	0.2
Resummation scale	0.2	0.8	0.5	1.1	0.2	1.6	0.8	0.4	0.2	0.7
αSFSR(M_Z)	0.5	0.7	0.7	0.8	1.7	0.7	0.2	1.0	0.2	1.2
αSISR(M_Z)	0.1	0.3	1.1	1.2	0.7	0.4	0.2	0.5	0.1	1.3
UE model	0.1	0.1	0.2	1.0	0.4	0.5	0.2	0.2	0.1	0.9
m_t	0.4	0.7	1.5	0.6	0.9	0.5	0.1	0.1	0.1	0.7
p_T(t)	1.4	4.4	4.5	2.8	2.1	6.7	0.2	0.5	0.2	0.3
Total	4.9	6.5	7.3	3.7	3.1	8.2	1.1	1.6	0.6	2.4

Results

Inclusive distributions

The normalized differential cross sections measured as functions of , , , , , , sphericity, aplanarity, C, and D are shown in Figs. 4, 5, 6, 7, 8, 9, 10, 11, 12 and 13, respectively. The distributions are obtained after unfolding the background-subtracted data and normalizing the result to unity. The result is compared to the simulations, whose settings are summarized in Table 1 and in the appendix. For the predictions, the statistical uncertainty is represented as an error bar. In the specific case of the Pw+Py8 setup, the error bar represents the envelope obtained by varying the main parameters of the CUETP8M2T4 tune, according to their uncertainties. The envelope includes the variation of the CR model, , , the parameter, and the scales at the ME level. Thus, the uncertainty band represented for the Pw+Py8 setup should be interpreted as the theory uncertainty in that prediction. For each distribution we give, in addition, the ratio between different predictions and the data.

Fig. 4

The normalized differential cross section as a function of is shown on the upper panel. The data (colored boxes) are compared to the nominal Pw+Py8 predictions and to the expectations obtained from varied or Pw+Py8 setups (markers). The different panels on the lower display show the ratio between each model tested (see text) and the data. In both cases the shaded (hatched) band represents the total (statistical) uncertainty of the data, while the error bars represent either the total uncertainty of the Pw+Py8 setup, computed as described in the text, or the statistical uncertainty of the other MC simulation setups

Fig. 5

Normalized differential cross section as function of , compared to the predictions of different models. The conventions of Fig. 4 are used

Fig. 6

Normalized differential cross section as function of , compared to the predictions of different models. The conventions of Fig. 4 are used

Fig. 7

Normalized differential cross section as function of , compared to the predictions of different models. The conventions of Fig. 4 are used

Fig. 8

Normalized differential cross section as function of , compared to the predictions of different models. The conventions of Fig. 4 are used

Fig. 9

Normalized differential cross section as function of , compared to the predictions of different models. The conventions of Fig. 4 are used

Fig. 10

Normalized differential cross section as function of the sphericity variable, compared to the predictions of different models. The conventions of Fig. 4 are used

Fig. 11

Normalized differential cross section as function of the aplanarity variable, compared to the predictions of different models. The conventions of Fig. 4 are used

Fig. 12

Normalized differential cross section as function of the C variable, compared to the predictions of different models. The conventions of Fig. 4 are used

Fig. 13

Normalized differential cross section as function of the D variable, compared to the predictions of different models. The conventions of Fig. 4 are used

In events the UE contribution is determined to have typically charged particles with , vectorially summing to a recoil of about 10. The distribution of the UE activity is anisotropic (as the sphericity is ), close to planar (as the aplanarity peaks at low values of 0.1), and peaks at around 0.75 in the C variable, which identifies three-jet topologies. The D variable, which identifies the four-jet topology, is instead found to have values closer to 0. The three-prong configuration in the energy flux of the UE described by the C variable can be identified with two of the eigenvectors of the linearized sphericity tensor being correlated with the direction of the -tagged jets, and the third one being determined by energy conservation. When an extra jet with is selected, we measure a change in the profile of the event shape variables, with average values lower by 20–40% with respect to the distributions in which no extra jet is found. Thus when an extra jet is present, the event has a dijet-like topology instead of an isotropic shape.

The results obtained with pythia8 for the parton shower simulation show negligible dependence on the ME generator with which it is interfaced, i.e., Pw+Py8 and MG5_aMC yield similar results. In all distributions the contribution from MPI is strong: switching off this component in the simulation has a drastic effect on the predictions of all the variables analyzed. Color reconnection effects are more subtle to identify in the data. In the inclusive distributions, CR effects are needed to improve the theory accuracy for or . The differences between the CR models tested (as discussed in detail in Sect. 3) are nevertheless small and almost indistinguishable in the inclusive distributions. In general the Pw+Py8 setup is found to be in agreement with the data, when the total theory uncertainty is taken into account. In most of the distributions it is the variation of that dominates the theory uncertainty, as this variation leads to the most visible changes in the UE. The other parton shower setups tested do not describe the data as accurately, but they were not tuned to the same level of detail as Pw+Py8. The Pw+Hw++ and Pw+Hw7-based setups show distinct trends with respect to the data from those observed in any of the pythia8-based setups. While describing fairly well the UE event shape variables, herwig++ and herwig 7 disagree with the , , and measurements. The sherpa predictions disagree with data in most of the observables.

For each distribution the level of agreement between theory predictions and data is quantified by means of a variable defined as:where () are the differences between the data and the model in the i-th (j-th) bin; here n represents the total number of bins, and is the inverse of the covariance matrix. Given that the distributions are normalized, the covariance matrix becomes invertible after removing its first row and column. In the calculation of Eq. 3 we assume that the theory uncertainties are uncorrelated with those assigned to the measurements. Table 3 summarizes the values obtained for the main models with respect to each of the distributions shown in Figs. 4, 5, 6, 7, 8, 9, 10, 11, 12 and 13. The values presented in the table quantify the level of agreement of each model with the measurements. Low per number of degrees of freedom (dof) values are obtained for the Pw+Py8 setup when all the theory uncertainties of the model are taken into account, in particular for the event shape variables. This indicates that the theory uncertainty envelope is conservative.

Table 3

Comparison between the measured distributions at particle level and the predictions of different generator setups. We list the results of the tests together with dof. For the comparison no uncertainties in the predictions are taken into account, except for the Pw+Py8 setup for which the comparison including the theoretical uncertainties is quoted separately in parenthesis

Observable	χ^2/dof
	Pw+Py8	Pw+Hw++	Pw+Hw7	MG5_aMC	sherpa
N_ch	30/11 (15/11)	33/11	17/11	34/11	95/11
∑p_T	24/13 (13/13)	129/13	56/13	30/13	37/13
∑p_z	8/11 (4/11)	34/11	20/11	9/11	18/11
p_T¯	12/9 (1/9)	40/9	56/9	6/9	56/9
p_z¯	2/9 (1/9)	9/9	32/9	1/9	36/9
|p→_T|	17/11 (7/11)	102/11	49/11	20/11	34/11
S	29/7 (3/7)	7/7	17/7	36/7	194/7
A	18/7 (1/7)	8/7	13/7	26/7	167/7
C	34/7 (4/7)	7/7	27/7	38/7	187/7
D	7/7 (1/7)	5/7	8/7	11/7	83/7

Profile of the UE in different categories

The differential cross sections as functions of different observables are measured in different event categories introduced in Sect. 5. We report the profile, i.e., the average of the measured differential cross sections in different event categories, and compare it to the expectations from the different simulation setups. Figures 14, 15, 16, 17, 18, 19, 20, 21, 22 and 23 summarize the results obtained. Additional results for , profiled in different categories of and/or jet multiplicity, are shown in Figs. 24 and 25, respectively. In all figures, the pull of the simulation distributions with respect to data, defined as the difference between the model and the data divided by the total uncertainty, is used to quantify the level of agreement.

Fig. 14

Average in different event categories. The mean observed in data (boxes) is compared to the predictions from different models (markers), which are superimposed in the upper figure. The total (statistical) uncertainty of the data is represented by a shaded (hatched) area and the statistical uncertainty of the models is represented with error bars. In the specific case of the Pw+Py8 model the error bars represent the total uncertainty (see text). The lower figure displays the pull between different models and the data, with the different panels corresponding to different sets of models. The bands represent the interval where . The error bar for the Pw+Py8 model represents the range of variation of the pull for the different configurations described in the text

Fig. 15

Average in different event categories. The conventions of Fig. 14 are used

Fig. 16

Average in different categories. The conventions of Fig. 14 are used

Fig. 17

Average in different categories. The conventions of Fig. 14 are used

Fig. 18

Average in different categories. The conventions of Fig. 14 are used

Fig. 19

Average in different categories. The conventions of Fig. 14 are used

Fig. 20

Average sphericity in different categories. The conventions of Fig. 14 are used

Fig. 21

Average aplanarity in different categories. The conventions of Fig. 14 are used

Fig. 22

Average C in different categories. The conventions of Fig. 14 are used

Fig. 23

Average D in different categories. The conventions of Fig. 14 are used

Fig. 24

Average in different categories. The conventions of Fig. 14 are used

Fig. 25

Average in different jet multiplicity categories. The conventions of Fig. 14 are used

The average charged-particle multiplicity and the average of the momentum flux observables vary significantly when extra jets are found in the event or for higher values. The same set of variables varies very slowly as a function of . Event shape variables are mostly affected by the presence of extra jets in the event, while varying slowly as a function of or . The average sphericity increases significantly when no extra jets are present in the event showing that the UE is slightly more isotropic in these events. A noticeable change is also observed for the other event shape variables in the same categories.

For all observables, the MPI contribution is crucial: most of the pulls are observed to be larger than 5 when MPI is switched off in the simulation. Color reconnection effects are on the other hand more subtle and are more relevant for , specifically when no additional jet is present in the event. This is illustrated by the fact that the pulls of the setup without CR are larger for events belonging to these categories. Event shape variables also show sensitivity to CR. All other variations of the UE and CR models tested yield smaller variations of the pulls, compared to the ones discussed.

Although a high pull value of the Pw+Py8 simulation is obtained for several categories, when different theory variations are taken into account, the envelope encompasses the data. The variations of and account for the largest contribution to this envelope. As already noted in the previous section, the Pw+Hw++, Pw+Hw7, and sherpa models tend to be in worse agreement with data than Pw+Py8, indicating that further tuning of the first two is needed.

Sensitivity to the choice of in the parton shower

The sensitivity of these results to the choice of in the parton shower is tested by performing a scan of the value defined by Eq. (3), as a function of or . The is scanned fixing all the other parameters of the generator. A more complete treatment could only be achieved with a fully tuned UE, which lies beyond the scope of this paper. While no sensitivity is found to , most observables are influenced by the choice of . The most sensitive variable is found to be and the corresponding variation of the function is reported in Fig. 26. A polynomial interpolation is used to determine the minimum of the scan (best fit), and the points at which the function increases by one unit are used to derive the 68% confidence interval (CI). The degree of the polynomial is selected by a stepwise regression based on an F-test statistics [72]. A value of is obtained, which is lower than the one assumed in the Monash tune [73] and used in the CUETP8M2T4 tune. The value obtained is compatible with the one obtained from the differential cross sections measured as a function of in different regions or in events with different additional jet multiplicities. Table 4 summarizes the results obtained. From the inclusive results, we conclude that the range of the energy scale that corresponds to the 5% uncertainty attained in the determination of can be approximated by a variation, improving considerably over the canonical [2, 0.5] scale variations.

Fig. 26

Scan of the as a function of the value of employed in the Pw+Py8 simulation, when the inclusive or the distribution measured in different regions is used. The curves result from a fourth-order polynomial interpolation between the simulated points

Table 4

The first rows give the best fit values for for the Pw+Py8 setup, obtained from the inclusive distribution of different observables and the corresponding 68 and 95% confidence intervals. The last two rows give the preferred value of the renormalization scale in units of , and the associated interval that can be used as an estimate of its variation to encompass the differences between data and the Pw+Py8 setup

p_T(ℓℓ) region	Inclusive	Away	Toward	Transverse
Best fit αSFSR(M_Z)	0.120	0.119	0.116	0.119
68% CI	[− 0.006, + 0.006]	[− 0.011, + 0.010]	[− 0.013, + 0.011]	[− 0.006, + 0.006]
95% CI	[− 0.013, + 0.011]	[− 0.022, + 0.019]	[− 0.030, + 0.021]	[− 0.013, + 0.012]
μ_R/M_Z	2.3	2.4	2.9	2.4
68% CI	[1.7, 3.3]	[1.4, 4.9]	[1.6, 7.4]	[1.7, 3.5]

Summary

The first measurement of the underlying event (UE) activity in dilepton events produced in hadron colliders has been reported. The measurement makes use of proton-proton collision data collected by the CMS experiment in 2016, and corresponding to 35.9. Using particle-flow reconstruction, the contribution from the UE has been isolated by removing charged particles associated with the decay products of the event candidates as well as with pileup interactions from the set of reconstructed charged particles per event. The measurements performed are expected to be valid for other final states, and can be used as a reference for complementary studies, e.g., of how different color reconnection (CR) models compare to data in the description of the jets from decays. The chosen observables and categories enhance the sensitivity to the modeling of multiparton interactions (MPI), CR and the choice of strong coupling parameter at the mass of boson () in the pythia8 parton shower Monte Carlo simulation. These parameters have significant impact on the modeling of production at the LHC. In particular, the compatibility of the data with different choices of the parameter in pythia8 has been quantified, resulting in a lower value than the one considered in Ref. [73].

The majority of the distributions analyzed indicate a fair agreement between the data and the powheg +pythia8 setup with the CUETP8M2T4 tune [18], but disfavor the setups in which MPI and CR are switched off, or in which is increased. The data also disfavor the default configurations in powheg +herwig++, powheg +herwig7, and sherpa. It has been furthermore verified that, as expected, the choice of the next-to-leading-order matrix-element generator does not impact significantly the expected characteristics of the UE by comparing predictions from powheg and MadGraph 5_amc@nlo, both interfaced with pythia8.

The present results test the hypothesis of universality in UE at an energy scale typically higher than the ones at which models have been studied. The UE model is tested up to a scale of two times the top quark mass, and the measurements in categories of dilepton invariant mass indicate that it should be valid at even higher scales. In addition, they can be used to improve the assessment of systematic uncertainties in future top quark analyses. The results obtained in this study show that a value of is consistent with the data. The corresponding uncertainties translate to a variation of the renormalization scale by a factor of .

Table 5

Variations of the Pw+Py8 setup used for the comparison with the measurements. The values changed with respect to the CUETP8M2T4 tune are given in the columns corresponding to each model. Further details on parameters or specificities of the models can be found in Refs. [4, 5, 17, 32, 32–34, 73]. For the Rope hadronization model two variations are considered: one with no CR and the other with the default CR model. The settings for the former are denoted in parenthesis in the last column

Parameter	Pw+Py8 simulation setups
	CUETP8M2T4	Extreme variations		Fine grain variations
				MPI/CR	Parton shower scale		CR including tt¯
		no	no	UE	ISR	FSR	ERD	QCD	Gluon	Rope (no CR)
		MPI	CR	up/down	up/down	up/down	on	based [32]	move [5]	[33, 34]
PartonLevel
MPI	on	off
SpaceShower
renormMultFac	1.0				4/0.25
alphaSvalue	0.1108
TimeShower
renormMultFac	1.0					4/0.25
alphaSvalue	0.1365
MultipartonInteractions
pT0Ref	2.2			2.20/2.128				2.174	2.3
ecmPow	0.2521							0.2521
expPow	1.6			1.711/1.562				1.312	1.35
ColorReconnection
reconnect	on		off							(off)
range	6.59			6.5/8.7
mode	0							1	2
junctionCorrection								0.1222
timeDilationPar								15.86
m0								1.204
flipMode									0
m2Lambda									1.89
fracGluon									1
dLambdaCut									0
PartonVertex
setVertex										on
Ropewalk
RopeHadronization										on
doShoving										on
doFlavour										on
PartonLevel
earlyResDec	off						on	on	on	on