Extraction and validation of a new set of CMS PYTHIA8 tunes from underlying-event measurements

New sets of CMS underlying-event parameters ("tunes") are presented for the PYTHIA8 event generator. These tunes use the NNPDF3.1 parton distribution functions (PDFs) at leading (LO), next-to-leading (NLO), or next-to-next-to-leading (NNLO) orders in perturbative quantum chromodynamics, and the strong coupling evolution at LO or NLO. Measurements of charged-particle multiplicity and transverse momentum densities at various hadron collision energies are fit simultaneously to determine the parameters of the tunes. Comparisons of the predictions of the new tunes are provided for observables sensitive to the event shapes at LEP, global underlying event, soft multiparton interactions, and double-parton scattering contributions. In addition, comparisons are made for observables measured in various specific processes, such as multijet, Drell-Yan, and top quark-antiquark pair production including jet substructure observables. The simulation of the underlying event provided by the new tunes is interfaced to a higher-order matrix-element calculation. For the first time, predictions from PYTHIA8 obtained with tunes based on NLO or NNLO PDFs are shown to reliably describe minimum-bias and underlying-event data with a similar level of agreement to predictions from tunes using LO PDF sets.


Introduction
Monte Carlo (MC) simulation codes describe hadron-hadron collisions with models based on several components. The hard scattering component of the event consists of particles from the hadronization of partons whose kinematics are predicted using perturbative matrix elements (MEs), along with partons from initial-state radiation (ISR) and final-state radiation (FSR) that are simulated using a showering algorithm. The underlying event (UE) consists of the beambeam remnants (BBR) and the particles that arise from multiple-parton interactions (MPI). The BBR are what remains after a parton is scattered out of each of the two initial beam hadrons. The MPI are additional soft or semi-hard parton-parton scatterings that occur within the same hadron-hadron collision. Generally, observables sensitive to the UE also receive contributions from the hard-scattering components. Accurately describing observables that are sensitive to the UE not only requires a good description of BBR and MPI, but also a good modeling of hadronization, ISR, and FSR. Standard MC event generators, such as PYTHIA8 [1], HERWIG [2,3], and SHERPA [4] have adjustable parameters to control the behavior of their event modeling. A set of these parameters, which has been adjusted to better fit some aspects of the data, is referred to as a tune.
In a previous study [5], we presented several PYTHIA8 and HERWIG++ UE tunes constructed for a center-of-mass energy √ s lower than 13 TeV. The CMS PYTHIA8 tune CUETP8M1 is based on the Monash tune [6], using both the NNPDF2.3LO parton distribution function (PDF) set [7]. The CMS PYTHIA8 tune CUETP8S1-CTEQ6L1 is based instead on the tune 4C [8]. Both tunes CUETP8M1 and CUETP8S1-CTEQ6L1 were constructed by fitting the CDF UE data at √ s = 900 GeV and 1.96 TeV [9] together with CMS UE data at √ s = 7 TeV [10]. A similar procedure was used for the determination of the HERWIG++ tune (CUETHppS1) with the CTEQ6L1 PDF set [11]. A collection of previously published tunes is documented in [6,8,[12][13][14][15].
In this paper, a new set of tunes for the UE simulation in the PYTHIA8 (version 8.226) event generator is obtained by fitting various measurements sensitive to soft and semi-hard MPI at different hadron collision energies [9,10], including data from √ s = 13 TeV [16]. These tunes are constructed with the leading order (LO), next-to-leading order (NLO), and next-tonext-to-leading order (NNLO) versions of the NNPDF3.1 PDF set [17] for the simulation of all UE components. Typically, the values of strong coupling used for the simulation of the hard scattering are chosen consistent with the order of the PDF set used.
The new tunes are obtained by fitting CDF UE data at √ s = 1.96 TeV [9], together with CMS UE data at √ s = 7 TeV [10] and at 13 TeV [16,18]. For the first time, we show that predictions obtained with tunes based on higher-order PDF sets are able to give a reliable description of minimum-bias (MB) and UE measurements with a similar level of agreement to predictions from tunes using LO PDF sets. We also compare the predictions for multijet, Drell-Yan, and top-antiquark (tt) processes from PYTHIA8 with new tunes in ME-parton shower (PS) merged configurations.
In Section 2 we describe observables that are sensitive to MB and UE: diffractive processes [19], where one or both protons remain intact after the collision; and double-parton scattering (DPS), where two hard scatterings occur within the same collision. In Section 3, we compare the tunes that were constructed before the data at √ s = 13 TeV were available ("Pre-13 TeV" tunes) with UE data measured at 13 TeV. Section 4 is dedicated to a general discussion of the choice of PDF sets and strong coupling values for the UE simulation. In Section 5 we describe the new tunes. Section 6 shows the validation of the new CMS PYTHIA8 tunes for multijet, Drell-Yan, tt, and DPS processes. Section 7 is the summary and conclusions.

Observables for characterizing minimum bias, underlying event, and double-parton scattering
Minimum bias is a generic term that refers to inelastic events that are collected with a loose event selection that has the smallest bias possible. The MB observables are constructed from data with little or no additional selection requirements. The majority of MB collisions are soft, with a typical transverse momentum scale p T 2 GeV. The UE is defined as the activity that is not associated with the particles originating from the hard scattering of two partons and is generally studied in events that contain a hard scattering with p T 2 GeV. The main contribution to the UE comes from color exchanges between the beam partons and is modeled in terms of MPI, BBR, and color reconnection (CR). The MB and UE observables have quite different kinematic properties because they are affected by different mixtures of hard and soft scattering processes.
As illustrated in Fig. 1, one can use the topological structure of a typical hard hadron-hadron collision to study the UE experimentally. On an event-by-event basis, a leading object is used to define regions of η-φ space that are sensitive to the modeling of the UE, where η is the pseudorapidity and φ is the azimuthal scattering angle defined in the xy plane. The azimuthal separation between charged particles and the leading object, ∆φ = φ − φ max , is used to define the UE-sensitive regions. Here φ max is the azimuth of the leading object and φ is the azimuth angle of an outgoing charged particle. The regions are labelled as 'toward' (|∆φ| ≤ 60 • ), 'away' (|∆φ| < 120 • ), and 'transverse' (60 • < |∆φ| ≤ 120 • ). The transverse region can further be Published UE studies used the charged-particle jet with the largest p T [16], the dilepton system in DY [20,21], or tt [22] events as the leading (i.e., the highest p T ) objects. The tunes from CDF and CMS data [9,10] made use of the charged particle with the largest p T (p max T ) as the "leading object", and use only charged particles with p T > 0.5 GeV and |η| < 0.8 to characterize the UE. The toward region contains the leading object, and the away region is expected to include the object recoiling against the leading one. Most of the UE contributions, i.e., PS and MPI, are contained in the two transverse regions. For events with multiple ISR or FSR emissions, transMAX often contains a third hard jet, while both transMAX and transMIN receive contributions from the MPI and BBR components. Typically, the transMIN observables are more sensitive to the MPI and BBR components of the UE.
Observables sensitive to UE contributions are the charged-particle multiplicity and the charged-particle scalar-p T sum densities in the η-φ space, measured in transMIN and trans-MAX. The tunes that are constructed by fitting such UE-sensitive observables are referred to as "UE tunes".
The PYTHIA8 MC event generator also simulates single-diffractive (SD) dissociation, doublediffractive (DD) dissociation, central-diffractive (CD), and nondiffractive (ND) processes [23], which contribute to the inelastic cross section in hadron-hadron collisions. In SD, CD, and DD events, one or both of the beam particles are excited into color singlet states, which then decay. The SD and DD processes correspond to color singlet exchanges between the beam hadrons, while CD corresponds to double color singlet exchange with a diffractive system produced centrally. For ND processes, color exchanges occur, the outgoing remnants are no longer color singlets, and a multitude of particles is produced. All processes except SD are defined as nonSD (NSD) processes. An NSD-enhanced sample is required to have an energy deposit in both the backward (−5 < η < −3) and the forward (3 < η < 5) regions of the detector. The details of the selection for different types of diffractive events can be found in Ref. [24].
Generally, MC models such as PYTHIA8 regularize the contributions of the primary hardscattering processes and MPI to the differential cross section by using a threshold parameter p 0 T . The primary hard-scattering processes and the MPI are regularized in the same way with this parameter. This threshold is expected to have a dependence on the center-of-mass energy of the hadron-hadron collision, √ s. The threshold at a reference center-of-mass energy √ s = 7 TeV is called pT0Ref. In PYTHIA8 the energy dependence is parameterized using a power law function with a reference energy parameter s 0 and an exponent . At a given center-of-mass energy, the amount of MPI depends on the threshold p 0 T , the PDF, and the overlap of the matter distributions of the two colliding hadrons. Smaller values of p 0 T result in larger MPI contributions because of a higher MPI cross section. Each MPI adds colored partons to the final state, creating a dense net of color lines that spatially overlap with the fields produced by the partons of the hard scattering and with each other. All the generated color lines may connect to each other according to the CR model.
Since PYTHIA8 regularizes both the cross section for MPI and the cross section of collisions with low-p T exchange using the p 0 T parameter, one can model the overall ND cross section by letting the p T of the primary hard scattering become small. In this simple approach, the UE in a hardscattering process is related to MB collisions. At the same center-of-mass energy, the activity in the UE of a hard-scattering process is greater than that of an average MB collision. In PYTHIA8, this is caused by the higher MPI activity in hard-scattering processes compared to a typical MB collision. By demanding a hard scattering, one forces the collision to be more central, i.e., with a small impact parameter between the protons, and this increases the probability of MPI. For MB collisions, peripheral collisions, where the impact parameter between the two colliding protons is large, are most common.
Typically MPI interactions contain particles with substantially lower p T ("softer"). However, occasionally two hard 2-to-2 parton scatterings can occur within the same hadron-hadron collision. This is referred to as DPS. Tunes that are constructed by fitting DPS-sensitive observables are referred to as "DPS tunes". Ultimately, one universal tune that simultaneously accurately describes observables in hard scattering events, as well as MB collisions, is desirable.
The goals of this paper are to produce improved 13 TeV PYTHIA8 tunes with well-motivated parameters, and to provide an investigation of the possible choices that can be made in PYTHIA8 which simultaneously describe a wide range of UE and MB measurements and are suitable for merged configurations, where a ME calculation is interfaced to the simulation of UE contributions.

Comparisons of predictions for UE observables from previous tunes to measurements at 1TeV
In this section, comparisons are presented between data collected at √ s = 13 TeV and predictions from tunes obtained using fits to measurements performed at lower center-of-mass energies. Figure 2 displays comparisons of CMS data at 13 TeV [16] for the transMIN and transMAX charged-particle p sum T densities, as functions of the leading charged-particle p max T . The data are compared with predictions from the PYTHIA8 tunes CUETP8S1-CTEQ6L1 [5], CUETP8M1 [5], and Monash [6].
The CMS Monash-based tune CUETP8M1 does not describe the central values of the data at √ s = 13 TeV well, nor does the original Monash tune. For example, CUETP8M1 and Monash tunes do not predict enough UE activity in the region with p max T > 5 GeV (the "plateau" region) of transMIN at 13 TeV, with a disagreement of ≈10% and ≈5%, respectively. The transMIN observables are very sensitive to MPI, which suggests that tune CUETP8M1 does not produce enough charged particles at 13 TeV. In addition, CUETP8M1 does not provide a good fit to the jet multiplicity in tt production either at 8 TeV or at 13 TeV [25,26]. High jet multiplicity tt events are sensitive to the modeling of the ISR. Hence, CUETP8M1 may not have the proper mixture of MPI and ISR. The ATLAS collaboration has also observed some discrepancies between the predictions of the A14 tune [12], used as standard tune for analyses of 7 and 8 TeV data, and the data at 13 TeV [27]. The CMS UE tunes were constructed by fitting CDF UE data at √ s = 900 GeV and 1.96 TeV, together with CMS UE data at √ s = 7 TeV.
Since no currently available tune is able to optimally reproduce the UE data at √ s = 13 TeV, we aim to produce improved PYTHIA8 UE tunes.

PDF and strong coupling values for the tunes
Two of the basic input parameters to the predictions are the choice of the order of the PDF sets and values of the strong coupling α S . These appear in the hard partonic MEs, the PS model, and the MPI model. The α S values used in simulations at LO or NLO are typically different. Traditionally, the perturbative order of the PDF is matched to the order of the ME calculation. Merged calculations capture some higher-order corrections with respect to the formal order of the ME calculation. Merging schemes, such as the k T -MLM [31] or CKKW [32,33], allow the combination of predictions of jet production using ME calculations with those from PS emissions for soft and collinear parton radiation at leading-log accuracy without double counting or dead regions. Merging can be applied also for processes generated at NLO. Using the same PDF set and α S value in the ME calculations and in the simulation of the various components of the PS is advocated in Ref. [34], and by the HERWIG7 and SHERPA Collaborations, especially when the PS simulation is merged with calculations of higher-order MEs. The PDF used for the hard process is constrained by the accuracy of the ME calculation. If we require the PDF to match between the ME and PS, simulations with a (N)NLO ME will also require a (N)NLO PDF in the PS. Depending on the process, this may not have a significant effect. For PS MC event generators, different strategies are adopted; CMS [5] and ATLAS [12] tunes are traditionally based on LO PDFs, PYTHIA8 [6] tunes are mostly based on LO PDFs, new SHERPA [4] tunes are based on NNLO PDFs, and HERWIG7 [30] provide tunes based on NLO PDFs. The usage of a LO PDF set in the UE simulation is motivated by the fact that MPI processes occur at very low energy scales, where a physical (positive) gluon distribution is required by the parton shower. However, there is no consensus on the choice of the order of the PDF. For example, in the NNPDF3.1 set at NNLO, the gluon distribution remains physical even at very low scales.
In the PYTHIA8 tunes produced prior to this paper, the values used for α S were often not the same as those used in the PDFs. For example, in the Monash tune, the FSR α S (m Z ), set to 0.1365, is obtained by fitting PYTHIA8 predictions to LEP event-shape measurements [6], the ISR α S (m Z ) is assumed to be equal to FSR α S (m Z ), and the hard scattering and MPI α S (m Z ) is set to 0.13 according to the value used in the LO PDF set. Even though the α S values are free parameters in event generators and various possibilities are viable, the usual course is to choose them consistent with the value used by the PDF set.
In this paper, a collection of new tunes is presented for PDF sets that are evaluated at different accuracies and tested against observables of MB, UE, and hard processes. The NNPDF3.1 PDF sets at the LO, NLO, and NNLO accuracy are used [17]. The LO PDF set uses an α S (m Z ) value of 0.13, while 0.118 [35] is the α S (m Z ) value used for the NLO and NNLO PDF sets. None of the central values of the PDF sets have negative values for any parton flavor in the phase space relevant for comparisons. Special care is required when applying these tunes at highx regions, where the parton distributions in NNPDF3.1 NLO and NNLO PDF may become negative, which implies an unphysical (negative) value of the calculated cross sections.
The UE simulation is performed by PYTHIA8, together with PS merged with a calculation of a higher-order or a multileg ME provided by external programs, such as POWHEG [36] or MAD-GRAPH5 aMC@NLO (MG5 aMC) [37]. The issue of combining external ME calculations with PS contributions is addressed by the merging procedure. The procedures considered in this paper are the "FxFx" [38] or the "POWHEG" [39] methods for merging higher-order (NLO) MEs to PS and the "MLM" method [31].  (upper row), and multiplicity (lower row) densities for particles with p T > 0.5 GeV in |η| < 2.0, as a function of the transverse momentum of the leading charged particle (p max T ), from the CMS √ s = 13 TeV analysis [16]. The data are compared with the PYTHIA8 tune Monash, the CMS PYTHIA8 tunes CUETP8S1-CTEQ6L1 and CUETP8M1, and the HERWIG7 (labelled as "H7") tune UE-MMHT. The ratios of the simulations to the data (MC/Data) are also shown, where the shaded band indicates the total experimental uncertainty in the data. Vertical lines drawn on the data points refer to the total uncertainty in the data. Vertical lines drawn on the MC points refer to the statistical uncertainty in the predictions. Horizontal bars indicate the associated bin width.
During this study, we also investigated the effect of imposing an additional rapidity (y) ordering to ISR in these merging calculations. The PYTHIA8 Monash tune includes a rapidity ordering for both ISR and MPI. The rapidity ordering acts as an extra constraint on the p T -ordered emissions, thus reducing the phase space for parton emission.

New CMS PYTHIA8 tunes at 13 TeV
In the following, a set of new 13 TeV PYTHIA8.226 tunes is presented with different choices of values of the strong coupling used in the modeling of the ISR, FSR, hard scattering, and MPI, as well as the order of its evolution as a function of the four-momentum squared Q 2 . We distinguish the new tunes according to the order of the PDF set used: LO-PDF, NLO-PDF, or NNLO-PDF. The tunes are labeled as CPi, where CP stands for "CMS PYTHIA8" and i is a progressive number from 1 to 5. Only five parameters related to the simulation of MPI, to the overlap matter distribution function [40], and to the amount of CR are constrained for the new CMS tunes. In all tunes, we use the MPI-based CR model [41]. The CP tunes are multipurpose tunes, aiming for a consistent description of UE and MB observables at several collision energies and a reliable prediction of the UE simulation in various processes when merged with higher-order ME calculations.
The settings, used in the determination of the new CMS PYTHIA8 UE tunes, are as follows: • Tune CP1 uses the NNPDF3.1 PDF set at LO, with α S values used for the simulation of MPI, hard scattering, FSR, and ISR equal to, respectively, 0.13, 0.13, 0.1365, and 0.1365, and running according to an LO evolution.
• Tune CP2 is a slight variation with respect to CP1, uses the NNPDF3.1 PDF set at LO, with α S values used for the simulation of MPI, hard scattering, FSR, and ISR contributions equal to 0.13, and running according to an LO evolution.
• Tune CP3 uses the NNPDF3.1 PDF set at NLO, with α S values used for the simulation of MPI, hard scattering, FSR, and ISR contributions equal to 0.118, and running according to an NLO evolution.
• Tune CP4 uses the NNPDF3.1 PDF set at NNLO, with α S values used for the simulation of MPI, hard scattering, FSR, and ISR contributions equal to 0.118, and running according to an NLO evolution.
• Tune CP5 has the same settings as CP4, but with the ISR emissions ordered according to rapidity.
The parameters related to the simulation of the hadronization and beam remnants are not varied in the fits and are kept fixed to the values of the Monash tune. The overlap distribution between the two colliding protons is modeled according to a double-Gaussian functional form with the parameters coreRadius and coreFraction. This parametrization of the transverse partonic overlap of two protons identifies an inner, denser part, the so-called core, and an outer less dense part. The coreRadius parameter represents the width of the core and the coreFraction, the fraction of quark and gluon content enclosed in the core. A double-Gaussian function is preferred for modeling the proton overlap over the negative exponential used in some previous tunes. Tunes using a double-Gaussian function tend to better reproduce the cross sections measured by the CMS experiment at √ s = 7 TeV [10], simultaneously as a function of charged-particle multiplicity and transverse momenta.
The parameter that determines the amount of simulated CR in the MPI-based model is varied in the fits. A small (large) value of the final-state CR parameter tends to increase (reduce) the final particle multiplicities.
The new CMS PYTHIA8 tunes are extracted by varying the parameters listed in Table 1 and by fitting UE observables at various collision energies. In the fitting procedure, we use the charged-particle and p sum T densities, measured in transMIN and transMAX regions as a function of p max T , as well as the charged-particle multiplicity as a function of pseudorapidity η, measured by CMS at √ s = 13 TeV [16,18]. In addition, we also use the charged-particle and p sum T densities as a function of the leading charged-particle p T , measured in transMIN and transMAX by CMS at √ s = 7 TeV [10] and by CDF at √ s = 1.96 TeV [9]. The predictions form a grid in the fivedimensional parameter space which is fitted using a third-order polynomial function. The uncertainty introduced in the fitted parameters due to the interpolation procedure is negligible compared with the quoted tune uncertainty. Results are found to be stable if one decreases this number to 100 or increases to 200, or uses a fourth-order polynomial function for the grid interpolation. The generated inelastic events include ND and diffractive (DD+SD+CD) contributions. The UE observables used to determine the tunes are sensitive to diffractive contributions only at very small p max T values (<3 GeV). The ND component is dominant for p max T values greater than ≈3.0 GeV, since the cross section of the diffractive components rapidly decreases as a function of the exchanged p T . Minimum-bias observables, such as the inclusive charged-particle multiplicity as a function of η, are sensitive to all contributions over the whole spectrum.
The fit is performed by minimizing the χ 2 function where the sum runs over each bin i of every observable O j . The f i (p) functions represent a parametrization of the dependence of the predictions in bin i on the tuning parameters, R i is the value of the measured observable in bin i, and ∆ i is the total experimental uncertainty of R i . The best fit values of the tuned parameters are shown in Table 2 for CP1 and CP2, i.e., the tunes using LO PDF sets, and in Table 3 for CP3, CP4, CP5, i.e., the tunes using NLO or NNLO PDF sets. Uncertainties in the parameters of these tunes are discussed in Appendix A.
No correlation across bins is included in the minimized χ 2 function.
The value of pT0Ref and its energy dependence is very different between tunes based on LO PDF sets and tunes based on NLO or NNLO PDFs. While pT0Ref is around 2.3-2.4 GeV for CP1 and CP2 tunes with ≈ 0.14-0.15, CP3, CP4, and CP5 tunes prefer much lower values for both pT0Ref (≈1.4-1.5) and (≈0.03-0.04). A value of of ≈0.03-0.04 corresponds to a very weak energy dependence of the threshold of the MPI cross section. These results can be understood by considering the shapes of the gluon densities at small x for the different PDF sets. In order to describe the UE observables, the rapidly increasing gluon densities at small x in LO PDF sets favor large values of pT0Ref. Meanwhile NLO and NNLO PDF sets, whose gluon densities are more flat at low x, need higher contributions of MPI, i.e., a small value of pT0Ref. Figure 3 shows the number of MPI observed for the various tunes and the gluon distribution at a reference scale of µ = 3 GeV for various NNPDF versions. The larger number of simulated MPI for NLO and NNLO tunes with respect to LO tunes is apparent.        We have found that the values of pT0Ref and also depend on the order of the running used for α S . In particular, fits based on NLO or NNLO PDF sets, i.e., CP3, CP4, or CP5, with an LO α S running prefer even smaller values for both pT0Ref and than the ones in the tunes obtained with an NLO α S running. This is because α S runs faster at NLO than at LO. When α S is run from the same value at the same scale (m Z ), the effective coupling at low scales is larger for NLO running than for LO running. Therefore, a lower pT0Ref is needed for NLO α S running than for LO α S running to obtain a similar number of MPI.
For tunes based on NLO and NNLO PDF sets, the value of pT0Ref is as low as the initial scale of the PDF Q 2 min . For interactions occurring at Q 2 which are lower than Q 2 min , the value of the PDF is left frozen to the value assumed at the initial scale.
The contribution from CR also changes among the different tunes and depends on the choice of PDF and its order. In particular, the amount of CR is also affected by the shape of the PDFs at small fractional momenta x.
Parameters related to the overlap matter distribution function differ between the different tunes. They are strongly correlated with the other UE parameters governing the MPI and CR contributions. In general, for a given value of the matter fraction (coreFraction), MPI contributions increase for decreasing values of the core radius (coreRadius). The inclusion of the Table 2: CMS PYTHIA8 LO-PDF tunes CP1 and CP2. Both the values at Q = m Z and the order of running with Q 2 of the strong coupling α S (m Z ) are listed. In these tunes, we use the Schuler-Sjöstrand diffraction model [44] and also include the simulation of CD processes. The number of degrees of freedom for tunes CP1 and CP2 is 63.
SpaceShower:rapidityOrder off off  Table 3: CMS PYTHIA8 NLO-PDF tune CP3 and NNLO-PDF tunes CP4 and CP5. Both the values at Q = m Z and the order of running with Q 2 of the strong coupling α S are listed. In these tunes, we use the Schuler-Sjöstrand diffraction model [44] and also include the simulation of CD processes. The number of degrees of freedom for tunes CP3, CP4, and CP5 is 63.  densities, as a function of the transverse momentum of the leading charged particle, p max T , from the CMS √ s = 13 TeV analysis [16]. Charged hadrons are measured with p T > 0.5 GeV in |η| < 2.0. The transMIN densities are more sensitive to the MPI, whereas the transMAX densities are more sensitive to ISR and FSR. The data are compared with the CMS PYTHIA8 LO-PDF tunes CP1 and CP2. The ratios of the simulations to the data (MC/Data) are also shown, where the shaded band indicates the total experimental uncertainty in the data. Vertical lines drawn on the data points refer to the total uncertainty in the data. Vertical lines drawn on the MC points refer to the statistical uncertainty in the predictions. Horizontal bars indicate the associated bin width. rapidity ordering for ISR in tune CP5 impacts the UE observables by reducing the number of charged particles, and needs to be compensated by a larger amount of MPI contributions.
The χ 2 per degree of freedom (dof) listed in Tables 2 and 3 refers to the quantity χ 2 (p) in Eq. (1), divided by dof in the fit. The eigentunes (Appendix A) correspond to the tunes in which the changes in the χ 2 (∆χ 2 ) of the fit relative to the best fit value equals the χ 2 value obtained in the tune, i.e., ∆χ 2 min = χ 2 . Such a variation of the χ 2 produces a tune whose uncertainty bands are roughly the same as the uncertainties in the fitted data points. This is the main motivation why this choice of variation was considered. For all tunes in Tables 2 and 3, the fit quality is good,    MC/Data Figure 6: The transMIN (upper left) charged-particle and (upper right) charged p sum T densities and the transMAX (lower left) charged-particle and (lower right) charged p sum T densities, as a function of the transverse momentum of the leading charged particle, p max T , from the CMS √ s = 7 TeV analysis [10]. Charged hadrons are measured with p T > 0.5 GeV in |η| < 0.8.
The data are compared with the CMS PYTHIA8 LO-PDF tunes CP1 and CP2. The ratios of simulations to the data (MC/Data) are also shown, where the shaded band indicates the total experimental uncertainty in the data. Vertical lines drawn on the data points refer to the total uncertainty in the data. Vertical lines drawn on the MC points refer to the statistical uncertainty in the predictions. Horizontal bars indicate the associated bin width. transMIN and transMAX regions. All predictions reproduce well the UE observables at √ s = 1.96, 7, and 13 TeV. Predictions from LO tunes are slightly better than the higher-order tunes in describing the energy dependence of the considered UE measurements.
In the region of small p max T values (p max T < 3 GeV), where contributions from diffractive processes are relevant, the predictions do not always reproduce the measurements and exhibit discrepancies up to 20%. Predictions from all of the new tunes cannot reproduce the UE data measured at √ s = 300 and 900 GeV [9]. Figure 10 shows the charged-particle multiplicity as a function of pseudorapidity for charged  densities, as a function of the transverse momentum of the leading charged particle, p max T , from the CMS √ s = 7 TeV analysis [10]. Charged hadrons are measured with p T > 0.5 GeV in |η| < 0.8.
The data are compared with the CMS PYTHIA8 (N)NLO-PDF tunes CP3, CP4, and CP5. The ratios of simulations to the data (MC/Data) are also shown, where the shaded band indicates the total experimental uncertainty in the data. Vertical lines drawn on the data points refer to the total uncertainty in the data. Vertical lines drawn on the MC points refer to the statistical uncertainty in the predictions. Horizontal bars indicate the associated bin width. particles in |η| < 2 measured by the CMS experiment at √ s = 13 TeV [18] in MB events.
These events were recorded with no magnetic field, so all particles irrespective of their p T are measured. Data are compared with the predictions of the new PYTHIA8 tunes. All of them are able to reproduce the measurement at the same level of agreement, independently of the PDF used for the UE simulation. We could not find any MB or UE observable where the level of agreement between data and predictions from the different tunes is significantly different.  densities, as a function of the transverse momentum of the leading charged particle, p max T , from the CDF √ s = 1.96 TeV analysis [9]. Charged hadrons are measured with p T > 0.5 GeV in |η| < 0.8.
The data are compared with the CMS PYTHIA8 LO-PDF tunes CP1 and CP2. The ratios of simulations to the data (MC/Data) are also shown, where the shaded band indicates the total experimental uncertainty in the data. Vertical lines drawn on the data points refer to the total uncertainty in the data. Vertical lines drawn on the MC points refer to the statistical uncertainty in the predictions. Horizontal bars indicate the associated bin width.

Validation of the new PYTHIA8 tunes
In this section, comparisons of the predictions obtained with the new tunes to various experimental measurements performed by the CMS experiment are provided. Unless otherwise stated, the comparisons are made at √ s = 13 TeV. We compare the CMS UE tunes with MB and UE data measured at central and forward pseudorapidities that are not used in the fits. We examine how well multijet, Drell-Yan, and top quark observables are predicted by MC simulations using higher-order ME generators merged with PYTHIA8 with the various new tunes. MC/Data Figure 9: The transMIN (upper left) charged-particle and (upper right) charged p sum T densities and the transMAX (lower left) charged-particle and (lower right) charged p sum T densities, as a function of the transverse momentum of the leading charged particle, p max T , from the CDF √ s = 1.96 TeV analysis [9]. Charged hadrons are measured with p T > 0.5 GeV in |η| < 0.8.
The data are compared with the CMS PYTHIA8 (N)NLO-PDF tunes CP3, CP4, and CP5. The ratios of simulations to the data (MC/Data) are also shown, where the shaded band indicates the total experimental uncertainty in the data. Vertical lines drawn on the data points refer to the total uncertainty in the data. Vertical lines drawn on the MC points refer to the statistical uncertainty in the predictions. Horizontal bars indicate the associated bin width.

Comparisons using event-shape observables
In this subsection, predictions of the new tunes are compared to event-shape observables measured at LEP, in electron-positron collisions. These observables are particularly sensitive to the value of α FSR S (m Z ). Given the leptonic initial state, there is no effect coming from the values of the MPI, color reconnection, and ISR parameters.
When predictions with PYTHIA 8 are used, an optimal value of α FSR S (m Z ) ∼ 0.13 is found, which best describes these observables, independent of the PDF used for the modeling of the PS evolution. Predictions obtained with MG5 aMC with up to 4 partons in the final state, and interfaced with the UE from the tune CUETP8M1 and the new PYTHIA 8 tunes CP2, CP3, and CP5 are considered ( Figure 11). Predictions using the tune CP2 do not describe the event-shape observables very well, with discrepancies with the data up to 30% in the T and T major . In particular, tune CP2 predicts too many isotropical events. A similar description is obtained for predictions of MG5 aMC+PYTHIA 8 with the tune CUETP8M1. A better agreement in the event-shape variables is observed for predictions using tune CP3 and CP5. A correct description of event-shape observables strongly depends on the value of the FSR strong coupling. The observations above indicate that when merged configurations are considered, i.e., MG5 aMC + PYTHIA8, where partons at higher multiplicities in the final state are simulated at the ME level, the description of event-shape observables degrades. A value of α FSR S (m Z ) ∼ 0.13 generally overestimates the number of final-state partons, while a lower α FSR S (m Z ) ∼ 0.12 performs better. At large values of T, where the hadronization effects become relevant, we observe a large difference between predictions from tunes using a small α FSR S (CP3 and CP5) and tunes using a large α FSR S (CP2 and CUETP8M1). These differences may be due to the interplay between the value of the strong coupling and the hadronization. Analyses particularly sensitive to hadronization should carefully evaluate the corresponding systematic uncertainties. In some cases retuning hadronization parameters may be desired.

Comparisons using MB and other UE observables
In this subsection, predictions of the new tunes are compared to the observables measured in MB collisions that are sensitive to contributions from soft emissions and MPI. Figure 13 shows the charged-particle multiplicity as a function of pseudorapidity [24] in NSD-enhanced and SD-enhanced event samples. The details of the selections can be found in Ref. [24]. These observables are sensitive to SD, CD, and DD dissociation. It is observed that predictions from all of the tunes are similar to each other and describe well the measurements for both considered + ++ + + + + + ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ ++ + ++ ++ ++ + + + + + + + + + + +  selections. This shows that the number of charged particles produced in diffractive processes and inelastic collisions is simultaneously described by the new CMS tunes. Figure 13 also demonstrates that tunes based on NNPDF3.1 PDF sets at orders higher than LO adequately describe the MB data. Figure 14 shows the UE observables, i.e., charged-particle multiplicity and p sum T densities [16], as a function of the p T of the leading jet reconstructed using just charged particles. The observables shown in Fig. 14 are from events selected without requiring any NSD-or SD-enhanced selections. The CMS UE tunes describe well UE-sensitive data measured using the leading charged-particle jet for p T jet > 10 GeV. Tunes based on NLO or NNLO PDF sets, i.e., CP3, CP4, and CP5, describe the region at lower p T jet better than CP1 and CP2, which are based on LO PDF sets. Predictions obtained with CP1 and CP2 underestimate the UE observables by about ≈15-20%. Predictions obtained with CP3, CP4, and CP5 describe the UE in events characterized using the leading charged particle, as well as those characterized by the leading charged-particle jet, quite well.
Predictions for observables measured in the forward region are compared with data and shown in Figs. 15 and 16. The energy flow, defined as the average energy per event [47], as a function of η with the Hadron Forward (HF) calorimeter [48] covering 3.15 < |η| < 5.20 and the CAS-TOR calorimeter [48] covering −6.6 < η < −5.2, is well reproduced by all tunes. A different level of agreement is achieved for predictions from the new CMS tunes for the spectrum of the total energy E measured in the CASTOR calorimeter at √ s = 13 TeV [49], displayed in Fig. 16.
In particular, the tunes based on LO PDF sets reproduce the energy spectrum well at large values (E > 2000 GeV), but have differences of up to 40% at low values (E < 800 GeV). The tunes using higher-order PDF sets are closer to the data at low energy values, with differences up to 20%, but tend to overestimate the energy at large values. This dissimilar behaviour is driven by the different pT0Ref values of the tunes. The fiducial inelastic cross sections [50], when two different selections are applied in the forward region, are not well reproduced by any of the new tunes or by CUETP8M1, with differences up to 10%. This might be because of the Schüler-Sjöstrand [44] diffraction model used in the simulation, which might have a suboptimal description of the low-mass diffractive components. A better description might be provided by tunes using the Donnachie-Landshoff [51] or minimum-bias Rockefeller [52] diffractive models.   Figure 14: The transMIN charged-particle multiplicity (left column) and p T sum densities (right column) for particles with p T > 0.5 GeV in |η| < 2.0 as a function of the transverse momentum of the leading charged-particle jet, p jet T , from the CMS √ s = 13 TeV analysis [16]. The upperrow plots show the LO tunes, while the lower-row plots show the higher-order tunes. The ratio of the simulations to the data (MC/Data) is also shown, where the shaded band indicates the total experimental uncertainty in the data. Vertical lines drawn on the data points refer to the total uncertainty in the data. Vertical lines drawn on the MC points refer to the statistical uncertainty in the predictions. Horizontal bars indicate the associated bin width.

Comparisons using observables in multijet final states
In this subsection, we present comparisons of observables measured in multijet final states. For these studies, the NLO dijet MEs implemented in the POWHEG event generator merged with the PYTHIA8 simulation of the PS and UE are used. The merging between the POWHEG ME calculations and the PYTHIA8 UE simulation is performed using the shower-veto procedure, which rejects showers if their transverse momentum is greater than the minimal p T of all final-state partons simulated in the ME (parameter p T hard = 2 GeV [53]). Variables in multijet events, such as jet transverse momenta or azimuthal dijet correlations, are expected to be less affected by MPI contributions, since jets at high p T (>100 GeV) mainly originate from the MC/Data Figure 16: The total energy spectrum measured in the pseudorapidity interval −6.6 < η < −5.2, from the CMS √ s = 13 TeV analysis [49]. The data are compared with (left) the CMS PYTHIA8 LO-PDF tunes CP1 and CP2, and with (right) the CMS PYTHIA8 NLO-PDF tune CP3 and the CMS PYTHIA8 NNLO-PDF tunes CP4 and CP5. The ratio of the simulations to the data (MC/Data) is also shown, where the shaded band indicates the total experimental uncertainty in the data. Vertical lines drawn on the data points refer to the total uncertainty in the data. Vertical lines drawn on the MC points refer to the statistical uncertainty in the predictions.
Horizontal bars indicate the associated bin width.
hard scattering or additional hard emissions, which are simulated in the POWHEG calculation by the ME formalism. However, the MPI contribution still has some impact because it adds an average energy offset to the event, which is then included in the jet reconstruction [54,55]. The predictions reproduce well inclusive jet cross sections as a function of jet p T at both central and forward jet rapidities, irrespective of the cone size (0.4 or 0.7) used for the jet clustering algorithm [56].  < 300 GeV and (right) 300 < p lead T < 400 GeV, from the CMS √ s = 13 TeV analysis [57]. The jets are reconstructed using the anti-k T jet finding algorithm [58,59] with a distance parameter of 0.4. The data are compared with predictions of the NLO dijet ME calculation from POWHEG, interfaced to the PYTHIA8 tunes CUETP8M1, CP2, CP4, and CP5. Tunes CP1 and CP3 are not shown in the plot but present a similar behavior as tunes CP2 and CP4. The ratios of simulations to the data (MC/Data) are also shown, where the shaded band indicates the total experimental uncertainty in the data. Vertical lines drawn on the data points refer to the total uncertainty in the data. Vertical lines drawn on the MC points refer to the statistical uncertainty in the predictions. Horizontal bars indicate the associated bin width. Figure 17 shows the normalized cross section [57] as a function of the azimuthal difference ∆φ 1,2 between the two leading jets for two different selections on the leading jet p T : 200 < p T < 300 GeV and 300 < p T < 400 GeV. The results indicate that UE tunes based on an NLO evaluation of α S (m Z ) describe the data better than UE tunes based on LO evolution. In particular, the better agreement is driven by the lower value of α ISR S (m Z ). In fact, predictions obtained with POWHEG merged with PYTHIA8 with the CUETP8M1 or CP2 tune exhibit a strong jet decorrelation, due to a large contribution from emissions simulated from the PS, and they overestimate the cross sections at small and medium ∆φ 1,2 values (∆φ 1,2 < 2.4). The PS component is reduced by the lower value of α S (m Z ), which increases the degree of correlation between the selected jets, resulting in a better description of the data by predictions of the CP4 and CP5 tunes. A similar outcome was also observed for an analogous measurement performed at the D0 experiment at √ s = 1.96 TeV [60]. In general, predictions obtained from POWHEG + PYTHIA8 tend to differ from the data at low and intermediate ∆φ 1,2 values (∆φ 1,2 < 2.7) by 10-40%.

Comparisons using observables sensitive to double-parton scattering
In this subsection, we present comparisons of predictions of the new tunes to DPS-sensitive observables measured by the CMS experiment at √ s = 7 TeV in final states with four jets (4j) [61], and with two jets originating from bottom quarks (b jets) [62] and two other jets (2b2j) [63].
The topology in the transverse plane of the physics objects measured in the final state is sensitive to contributions from DPS. In particular, the 4j analysis performed by the CMS experiment requires two jets at high p T (hard jets) and two jets at low p T (soft jets); the 2b2j measurement selects two jets originating from b quarks and two other jets (light-flavor jets). Both of them measured the ∆S observable, defined as: where p T,1 refers to the momentum of the hard-jet or bottom jet pair system and p T,2 to that of the soft-jet or light-flavor jet pair system. This variable relates the production planes of the hard (bottom) jet and soft-(light-flavor) jet pairs. Details of the event selection and of the specific analyses can be found in Refs. [61] and [63].
Assuming that the two hard scatterings occurring within the same collision are completely independent of each other, the DPS cross section for a given process can be expressed through the inclusive partonic cross sections of the two single scatterings and an effective cross section, σ eff . In a geometrical approach, this cross section is related to the transverse size of the proton and to the total inelastic proton-proton (pp) cross section [64,65]. When no correlations among the partons inside the proton are present, σ eff is similar to the inelastic pp cross section. In this simple factorized approach, one expects σ eff to be independent of the partonic final states of the two hard processes occurring within the same collision. In PYTHIA8, the value of σ eff is calculated by dividing the ND cross section by the so-called "enhancement factor", which depends on the parameters of the overlap matter distribution function and on pT0Ref [40]. For central pp collisions, the enhancement factor tends to be large, translating to a lower value of σ eff and a larger DPS contribution. For peripheral interactions, enhancement factors are small, giving large values of σ eff and a small DPS contribution. Table 4 shows the values of σ eff published by the CMS Collaboration for the 4j and the 2b2j measurements. A previous study [5] concluded that observables sensitive to semi-hard MPI and those sensitive to DPS cannot be described by a single set of parameters. Table 5 displays the σ eff values obtained from the new CMS UE tunes. The central values of σ eff are consistent among the new tunes and are slightly larger than the values of the DPS-based tunes [5]. Figure 18 shows the comparisons of predictions obtained from PYTHIA8 with tunes CUETP8M1, CP2, CP4, and CP5 to the DPS observables measured in the 4j and 2b2j final states. Predictions from the CP2 tune based on a LO PDF set describe the central values better than the CP4 and CP5 tunes based on an NNLO PDF set or the old tune CUETP8M1. This is due to the different pT0Ref value used by CP2, CP4, and CP5, which determines the amount of simulated MPI. The value of the pT0Ref parameter is driven by the distribution of the gluon distribution function at low x, which is very different in LO and NNLO PDF sets. Additionally, predictions obtained with CP4 describe the DPS-sensitive observables better than CP5. This is due to the different rapidity ordering used for the PS emissions in the two tunes. By removing the rapidity ordering for the PS emissions (CP4), the simulation produces more radiation and decreases the correlation between the selected jet pairs compared to CP5. This reduced jet correlation tends to mimic a DPS event by producing low values of ∆S. We have checked that the observables sensitive to color coherence, which were measured by the CMS experiment at √ s = 7 TeV [66], are well described by predictions from both CP4 and CP5 tunes, despite the difference in the rapidity ordering of the PS simulation between the two tunes.  Figure 18: The correlation observable ∆S measured in 4j (left) and 2b2j (right) production, compared to predictions of PYTHIA8 tunes CUETP8M1, CP2, CP4, and CP5, from the CMS √ s = 7 TeV analyses [61,63]. Tunes CP1 and CP3 are not shown in the plot but show a similar behaviour as, respectively, tunes CP2 and CP4. The ratios of simulations to the data (MC/Data) are also shown, where the shaded band indicates the total experimental uncertainty in the data. Vertical lines drawn on the data points refer to the total uncertainty in the data. Vertical lines drawn on the MC points refer to the statistical uncertainty in the predictions. Horizontal bars indicate the associated bin width.

Comparisons using observables in top quark production
In the following, we investigate how the new PYTHIA8 tunes describe the CMS tt data when different ME generators, namely POWHEG and MG5 aMC, are employed. Both ME configurations use the NNPDF3.1 NNLO PDF with α S (m Z ) = 0.118 and assume a top quark mass (m t ) value of 172.5 GeV.
In the POWHEG configuration, the ME heavy quark production mode [36,39,68] is used. In this configuration, POWHEG simulates inclusive tt production at NLO, where the first additional jet is computed at LO, while MG5 aMC performs the calculation with up to 2 additional jets at NLO, with a third jet simulated at LO. The POWHEG generator scales the real emission cross section by a damping function that controls the ME-PS merging and that regulates the high-p T radiation. The damping variable used in the POWHEG simulation is set to 1.379 times m t , a value derived from data at √ s = 8 TeV in the dilepton channel using a similar ME calculation and assuming the CP5 tune. The factorization and renormalization scales are assumed equal to the transverse mass of the top quark, m t T = √ m 2 t + p 2 T . The minimum p T for the emission of light quarks in POWHEG is 0.8 GeV. The pThard parameter is set to 0 and the POWHEG hardness criterion, defined by the pTdef option, is set to 1. The merging scale in MG5 aMC is set to 40 GeV, and the threshold applied to regulate multijet MEs in the MG5 aMC FxFx merging procedure, is 20 GeV.
Distributions [69] in the lepton+jets channel are compared to predictions from different tunes using various settings, namely, POWHEG + PYTHIA8, and MG5 aMC+PYTHIA8 with FxFx merging [38], referred to as MG5 aMC [FxFx] hereafter, with the CUETP8M1, CP2, CP4, and CP5 tunes. Figure 19 (upper panel) displays the normalized tt cross section in bins of p T of the top quark decaying leptonically (t ), in data and simulation. For all tunes, POWHEG + PYTHIA8 predictions have deviations below 10% with respect to the central values of the data. The central values of predictions from MG5 aMC [FxFx] and data agree within ≈10% for p T (t ) < 400 GeV and within ≈20% for higher p T . . The central values predicted by POWHEG + PYTHIA8 are in good agreement with data when CP5 tune is used. The value of α ISR S (m Z ) in combination with the rapidity ordering for ISR in the PYTHIA8 simulation affects the additional jet distribution in tt events. Predictions obtained from POWHEG + PYTHIA8 overestimate the data when a high value of α ISR S (m Z ) ≈ 0.13 is used (CUETP8M1 and CP2 tunes) irrespective of rapidity ordering for ISR. It is observed that even when α ISR S (m Z ) = 0.118 is used, predictions from the CP4 tune overshoot the data at high jet multiplicities. A much better agreement of central values is obtained only when rapidity ordering for ISR is switched on in the PYTHIA8 simulation and α ISR S (m Z ) = 0.118 is used as in the CP5 tune. Predictions from MG5 aMC [FxFx]+PYTHIA8 with CUETP8M1, CP2, CP4, and CP5 tunes describe the central values of the data reasonably well.
Comparisons are also made using jet substructure observables in tt events in the lepton+jets channel using measurements by CMS at √ s = 13 TeV [70]. Figure 20 displays the comparisons using the distribution of the angle between two groomed subjets, ∆R g , which is found to be the most sensitive to α FSR S (m Z ) [70]. The data are compared to simulations with the tunes CUETP8M1, CP2, CP4, and CP5, as well as CP5 FSR up (α FSR S (m Z ) = 0.122), CP5 FSR down (α FSR S (m Z ) = 0.115), and CP5 with CMW rescaling. It is observed that tunes with higher α FSR S (m Z ) (CUETP8M1, CP2, and CP5 FSR up) describe the data better. Tune CP5 with CMW rescaling resolves the discrepancy of CP5 at high ∆R g , but worsens the description at ∆R g ∼ 0.27 compared to CP5. It should be noted that a fit to the ∆R g distribution using a b-enriched sample yields α FSR S (m Z ) = 0.130 +0.016 −0.020 [70] without CMW rescaling, while a fit to the distirubtion of the UE observable p T measured in tt events yields α FSR S (m Z ) = 0.120 ± 0.006 [22]. Therefore, in tt events, UE and jet substructure observables prefer different central α FSR S (m Z ) values, but they are compatible within uncertainties.   [70] of the angle between two groomed subjets, ∆R g in tt events predicted by POWHEG + PYTHIA8 for the different tunes. The data are compared to tunes CUETP8M1, CP2, CP4, and CP5 (left). Tunes CP1 and CP3 are not displayed but they present a similar behavior as tunes CP2 and CP4, respectively. The data are also compared to CP5, CP5 FSR up, CP5 FSR down, and CP5 with CMW rescaling (right). The ratios of simulations to the data (MC/Data) are also shown, where the shaded band indicates the total experimental uncertainty in data. Vertical lines drawn on the data points refer to the total uncertainty in the data. Horizontal bars indicate the associated bin width.

Comparisons using observables in W and Z boson production
In this subsection, we present a validation of the new CMS UE tunes for observables measured in events with a W or Z boson in the final state at √ s = 13 TeV. For the comparisons, we use predictions obtained with MG5 aMC + PYTHIA8 at LO using the k T -MLM merging scheme, and at NLO using the FxFx merging scheme. The k T -MLM merging scale is set to 19 GeV, while for FxFx the corresponding scale is set to 30 GeV. In both cases the MEs include the final states with 0, 1, 2, and 3 partons, and up to 2 partons are calculated at NLO precision in the FxFx case. To ease the comparison of the different tunes, the same PDF, NNPDF3.1 NNLO, and α S (m Z ) = 0.118 are used for the ME calculation independently of the tune.
First, UE observables [21] in Drell-Yan events in an invariant mass window of 81-101 GeV around the Z boson peak for muonic decays are studied. The charged-particle density and transverse momentum sum are measured as a function of the Z boson p T in the three regions introduced in Section 6.2: toward, away, and transverse. The regions are defined with respect to the Z boson direction. The measurements are compared with FxFx predictions obtained with the CUETP8M1, CP2, CP4, and CP5 tunes in Fig. 21. The measurements are, in general, well-described by all tunes.
The description of the cross section as a function of the jet multiplicity is also investigated in Z +jets [71] and W +jets [72] final states. The Z +jets measurement is restricted to the phase space where the two leptons have p T > 20 GeV and |y| < 2.4 and the dilepton mass lies in a ±20 GeV window around 91 GeV. The momenta of the photons inside a cone of ∆R < 0.1 are added to the lepton momentum in order to partly recover the energy lost by FSR. Jets are clustered using the anti-k T algorithm with R = 0.4 and must satisfy the criteria p T > 30 GeV and |y| < 2.4. The distance between the selected leptons and the leading jet ∆R( , j) must be greater than 0.4. For the W +jets measurement, the phase space is restricted by a transverse mass requirement,  √ s = 13 TeV [21], and compared with the predictions obtained by an inclusive NLO ME calculated by MG5 aMC, interfaced to the UE simulation of PYTHIA8 with the CUETP8M1, CP2, CP4, and CP5 tunes. Tunes CP1 and CP3 are not shown in the plot but present a similar behaviour as, respectively, tunes CP2 and CP4. The ratios of simulations to the data (MC/Data) are also shown, where the shaded band indicates the total experimental uncertainty of the data. Vertical lines drawn on the data points refer to the total uncertainty in the data. Vertical lines drawn on the MC points refer to the statistical uncertainty in the predictions. Horizontal bars indicate the associated bin width.
m T > 50 GeV, and by requirements on the muon, p T > 25 GeV and |y| < 2.4. In the Z +jets measurements the same clustering algorithm, the FSR recovery prescription described above, and the lepton jet separation requirement are applied.
The comparisons of the jet multiplicities to various predictions are shown in Fig. 22. The measurement of the cross section inclusive in the number of jets, N, is not available for the W +jets analysis and the lower plots start at N = 1. The k T -MLM predictions of the jet multiplicity have little sensitivity to the UE and PS tunes, so all the tunes provide a good description this observable, with a slightly better agreement observed for the CP2 tune. In the case of the FxFx sample, the CP5 tune predicts fewer events with a jet multiplicity of more than four with respect to the measurement. The deficit increases for increasing jet multiplicities. The CUETP8M1 tune shows a similar behaviour, though.
Predictions using the new CMS UE tunes are also compared with the p T balance between the Z boson and the jets with p T > 30 GeV and |y| < 2.4 using the variable p T bal = | p T (Z) + ∑ jets p T (j i )| [71]. This variable is sensitive to PS and UE. The comparison is shown in Fig. 23 for events with at least one jet. Differences between the tunes are significant only in the region below ≈20 GeV. The discrepancy in this region for the FxFx samples indicates that the distribution peaks at lower values for CP4 and CP5 than in data.
Results of Ref. [71] are also used to validate the description of the transverse momentum of the weak vector boson in Z + ≥ 1 jet events. The comparison is shown in Fig. 24. The new tunes provide similar descriptions for this distribution. Predictions using k T -MLM achieve a poor agreement with the data, independently of the UE tune, with respect to FxFx, which is able to describe the transverse momentum of the Z boson at p T > 10 GeV. The region below 10 GeV is poorly described for both FxFx and k T -MLM and the new tunes, but is well-described by predictions using the CUETP8M1 tune.
To summarize the study of weak vector boson production, the CP2, CP4, and CP5 tunes provide similar descriptions of the UE observables with a reasonable agreement with the data. In general, the CP2 tune performs better in describing variables such as p T bal and p T (Z). For the jet multiplicity, the CP2 and CP4 tunes are equally good in describing the measurement, whereas CP5 tends to undershoot the PS dominated region with at least five jets with a significance of 3.5 standard deviations.

Summary and conclusions
A new set of tunes for the underlying-event (UE) simulation in the PYTHIA8 event generator is obtained by fitting various measurements sensitive to soft and semihard multipartonic interactions at different collision energies. To derive these tunes, the leading order (LO), nextto-leading order (NLO), or next-to-next-to-leading order (NNLO) versions of the NNPDF3.1 parton distribution function (PDF) set for the simulation of the underlying-event components are used. In these tunes, the values of the strong coupling, α S (m Z ), used for the simulation of hard scattering, initial-and final-state radiation, and multiple-parton interactions are chosen consistent with the order of the PDF used. In the LO NNPDF3.1 set, α S (m Z ) = 0.130, whereas for the NLO and NNLO NNPDF3.1 sets, α S (m Z ) = 0.118. In general, the combination of contributions from multiple-parton interactions and parton-shower emissions is crucial to give a good description of variables measured in soft-collision events. The infrared threshold is relatively independent of center-of-mass energy when using NLO or NNLO PDF sets. Irrespective of the specific PDF used, predictions from the new tunes reproduce well the UE measurements at center-of-mass energies √ s = 1.96 and 7 TeV. A significant improvement in the description of UE measurements at 13 TeV is observed with respect to predictions from old tunes that  Figure 22: Comparison with the measurement [71,72] of the inclusive jet multiplicity in Z +jets (upper) and W +jets (lower) events predicted by MG5 aMC + PYTHIA8 with k T -MLM merging (left) and FxFx merging (right) for the different tunes. Tunes CP1 and CP3 are not shown in the plot but present a similar behaviour as, respectively, tunes CP2 and CP4. The ratios of simulations to the data (MC/Data) are also shown, where the shaded band indicates the total experimental uncertainty in the data. Vertical lines drawn on the data points refer to the total uncertainty in the data. Vertical lines drawn on the MC points refer to the statistical uncertainty in the predictions. Horizontal bars indicate the associated bin width.    Comparison with the measurement [71] of the p T (Z) predicted by MG5 aMC + PYTHIA8 with k T -MLM merging (left) and FxFx merging (right) for the different tunes. Tunes CP1 and CP3 are not shown in the plot but they present a similar behaviour as tunes CP2 and CP4, respectively. The ratios of simulations to the data (MC/Data) are also shown, where the shaded band indicates the total experimental uncertainty in data. Vertical lines drawn on the data points refer to the total uncertainty in the data. Vertical lines drawn on the MC points refer to the statistical uncertainty in the predictions. Horizontal bars indicate the associated bin width.
were extracted using data at lower collision energies. For the first time, predictions based on higher-order PDF sets are shown to give a reliable description of minimum-bias (MB) and UE measurements, with a similar level of agreement as predictions from tunes using LO PDF sets.
Predictions of the new tunes agree well with the data for MB observables measured at pseudorapidities in the central (|η| < 2.4) and forward (3.2 < |η| < 4.7) regions. The new CMS tunes simultaneously describe the number of charged particles produced in diffractive processes and MB collisions. Neither the new CMS tunes nor the CUETP8M1 tune describe the very forward region (−6.6 < η < −5.2) well.
Measurements sensitive to double-parton scattering contributions are reproduced better by predictions using the LO PDF set in the UE simulation, without rapidity ordering of the initialstate shower.
The UE simulation provided by the new tunes can be interfaced to higher-order and multileg matrix element generators, such as POWHEG and MG5 aMC, without degrading the good description of UE observables. Such predictions also reproduce well observables measured in multijet final states, Drell-Yan, and top quark production processes. The central values of the normalized tt cross section in bins of the number of additional jets predicted by POWHEG +PYTHIA8 overestimate the data when a high value of α ISR S (m Z ) 0.130 is used (CMS PYTHIA8 CP1 and CP2 tunes). Even when α ISR S (m Z ) = 0.118 is used, the CP4 tune overestimates the data at high jet multiplicities. This is cured by the rapidity ordering of the initialstate shower (CP5 tune). Measurements of azimuthal dijet correlations are also better described when a value of α ISR S (m Z ) = 0.118 is used in predictions obtained with POWHEG merged with PYTHIA8.
Comparisons with LEP event-shape observables and the distribution of the angle between two groomed subjets (∆R g ) in tt events at the LHC show that in ME-PS merged configurations CMW rescaling is disfavored. It is also found that ∆R g is better described by tunes with α FSR S (m Z ) higher than ∼0.120 while LEP event-shape observables and UE event observables in tt events prefer a central value ∼0.120 [22].
All of the new CMS tunes are supplied with their eigentunes, which can also be used to determine the uncertainties associated with the theoretical predictions. We show that predictions using the new tunes based on PDFs determined at LO, NLO, and NNLO agree reasonably well with the measurements, and that the new tunes can also be applied to LO and NLO calculations merged with parton showers, multiple-parton interactions, and hadronization.

A Tables of tune uncertainties
This section provides the values of the parameters corresponding to the uncertainties when the new CMS PYTHIA8 tunes are used. The tune uncertainty is obtained by extracting the eigentunes, which are defined by a change in the χ 2 of the fit that equals the absolute χ 2 value obtained in the tune. The eigentunes refer to the variations of the tunes along each of the maximally independent directions in the parameter space, obtained by using the covariance matrix in the region of the best tune. The number of directions defined in the parameter space equals the number of free parameters n used in the fit and results into 2n parameter variations, i.e., eigentunes. These variations represent a good set of systematic uncertainties in the given tune.
The estimations of the tune uncertainties, which have 2n parameter variations, i.e., 10 for the new CMS PYTHIA8 tunes, are very time consuming in analyses, since for each variation separate samples must be produced. Therefore, a lower number of variations is preferred. Hence, two variations, one "up" and one "down", are defined. For the definition of the two variations, predictions using the parameters of the eigentunes are implemented for the UE observables at √ s = 13 TeV and their differences with respect to the central predictions are added in quadrature. This procedure is applied in each bin and tune uncertainties are estimated without any correlation across the different bins. Positive differences between central predictions and tune variations define the upper edge of the bin-by-bin uncertainty, while negative differences define the lower edge of the bin-by-bin uncertainty. By following the same approach used for the extraction of the central values of the new CMS PYTHIA8 tunes, the upper edge is fitted to obtain the "up variation", while the "down variation" is obtained by fitting the lower edge. The parameters of the up-and down-variations are listed in Table 6 for the tunes using LO PDF sets, and in Table 7 for the tunes using (N)NLO PDF sets. We checked that for a wide range of MB and UE observables at √ s = 13 TeV predictions from up-and down-variations, obtained by including the full set of eigenvalues, reproduce well the upper and the lower edge of the predictions. Hence, tune uncertainties estimated by evaluating predictions of up-and downvariations represent a reliable way of estimating the systematic uncertainties in the tunes. The correlation matrix for the fit of the CP5 tune is displayed in Table 8. It is retrieved by evaluating the correlation of the parameter variations obtained in the eigentunes. Variations of the values of the ISR and FSR are also studied, in order to check the consistency of the selected α ISR S (m Z ) and α FSR S (m Z ) values selected for the tunes and to estimate the allowed range of α ISR S (m Z ) and α FSR S (m Z ) values in the PS using the CP5 tune. Starting from tune CP5, the value of α ISR S (m Z ) is fitted to UE observables measured by CMS at √ s = 13 TeV. The same procedure is repeated when α FSR S (m Z ) is fitted. The parameters obtained from the fits are shown in Table 9, along with the up and down variation. Figure 25 shows the predictions of the CP5 tune, with the corresponding variation bands rela-   Figure 25: The variations allowed by the CP5 tune when α ISR S (m Z ) (blue band) and α FSR S (m Z ) (red band) are left free in the fit for charged-particle (left) and charged p sum T (right) density in the transMIN region at √ s = 13 TeV. Vertical lines drawn on the data points refer to the total uncertainty in the data. The grey band represents the total UE uncertainty for the tune CP5.
Horizontal bars indicate the associated bin width.