Optimizing the parton shower model in PYTHIA with pp collision data at $\sqrt{s}$ = 13 TeV

TL;DR

This paper describes the optimization of the PYTHIA8 parton shower model using recent 13 TeV proton-proton collision data from the LHC, improving the simulation of quark and gluon evolution into jets.

Contribution

The paper introduces a new tuning of the PYTHIA8 parton shower model based on recent experimental data, enhancing its accuracy in simulating hadron collision events.

Findings

Optimized PYTHIA8 model better matches experimental data

Improved simulation of event shape variables

Enhanced prediction of jet cross-sections

Abstract

Production of quarks and gluons in hadron collisions tests Quantum Chromodynamics (QCD) over a wide range of energy. Models of QCD are implemented in event generators to simulate hadron collisions and evolution of quarks and gluons into jets of hadrons. PYTHIA8 uses the parton shower model for simulating particle collisions and is optimized using experimental observations. Recent measurements of event shape variables and jet cross-sections in pp collisions at $\sqrt{s}$ = 13 TeV at the Large Hadron Collider have been used to optimize the parton shower model as used in PYTHIA8.

Figures

Tables

Table 1. Table 1 : Optimization of three parameters of pythia8 is shown along with their default values in the Monash tune and the sampling range.

pythia8	Monash	Sampling range	Optimized
Parameters set	values		values
SpaceShower:alphaSvalue	0.1365	$0.1092-0.1638$	${0.11409}_{-0.00073}^{+0.00078}$
TimeShower:alphaSvalue	0.1365	$0.1092-0.1638$	${0.15052}_{-0.00076}^{+0.00077}$
SpaceShower:PTmaxFudge	1.0	$0.6-1.4$	${0.9323}_{-0.0064}^{+0.0065}$

Table 2. Table 2 : The goodness-of-fit ( χ 2 / NDF superscript 𝜒 2 NDF \mathrm{\chi^{2}/NDF} ) for each ESV corresponding to the Monash tune (left column) and the present optimization (right column) are shown for each H T,2 subscript 𝐻 T,2 H_{\textrm{T,2}} range.

$H_{\text{T,2}}$ range	${\tau }_{⟂}$		${\rho }_{\text{Tot}}$		$B_{\text{T}}$		${\rho }_{\text{Tot}}^{\text{T}}$
in GeV	Monash	This	Monash	This	Monash	This	Monash	This
	tune	study	tune	study	tune	study	tune	study
73-93	8.49	16.40	38.11	16.86	69.84	22.12	12.32	20.19
93-165	3.90	9.96	19.44	7.31	23.79	6.49	4.16	7.82
165- 225	1.14	5.39	13.63	4.38	21.30	7.17	2.56	4.60
225-298	1.81	4.59	11.44	3.58	23.51	7.42	2.86	3.43
298-365	1.54	3.53	6.29	1.68	12.19	3.78	2.03	3.17
365-452	1.36	3.95	7.42	2.29	8.53	5.15	1.07	3.68
452-557	1.72	2.40	6.29	2.09	10.10	4.56	2.61	1.88
$>557$	1.16	3.88	6.38	2.76	8.46	6.23	3.43	3.91

Table 3. Table 3 : Optimized pythia shows better agreement in terms of χ 2 / NDF superscript 𝜒 2 NDF \mathrm{\chi^{2}/NDF} values with the CMS measurement of inclusive jet cross section for anti- ​ k T anti- subscript 𝑘 T \textrm{anti-}k_{\textrm{T}} jets with R=0.4, 0.7 and also for ALTAS measurement of inclusive jet cross section with R=0.4

	CMS(R$=0.4$)		CMS(R$=0.7$)		ATLAS(R$=0.4$)
Rapidity	Monash	This	Monash	This	Monash	This
range	tune	study	tune	study	tune	study
	2.10	0.63	2.17	0.68	9.38	1.13
	1.73	0.41	2.11	0.50	8.20	0.96
	2.03	0.46	3.02	0.97	9.33	2.05
	0.89	0.25	1.30	0.40	6.92	1.56
	0.72	0.58	1.21	0.76	5.64	1.39
	0.36	0.31	0.27	0.15	3.75	1.23
	0.28	0.14	0.28	0.23	-	-

Equations

$$d{\mathcal{P}}_a=\frac{d{p}_{\perp }^2}{p_{\perp }^2}\sum _{b,c}\frac{{\alpha }_s(p_{\perp }^2)}{2\pi }{P}_{a\to bc}(z)dz$$

$$p_{{\perp }_{\text{evol}}}^2$$

$$\text{d}{\mathcal{P}}_b$$

$$\text{d}{\mathcal{P}}_a$$

Full text

1 Introduction

High energy collisions of hadrons produce quarks and gluons, collectively known as partons, which evolve into a large number of stable hadrons. Successive stages of this process involve lower energy and hence increasing value of the strong coupling $\mathrm{\alpha_{s}}$. This makes analytical calculations based on perturbation theory prohibitively difficult as well as unreliable. Thus, a combination of analytical calculations in the early stages, including the hard scattering, and approximate numerical calculations later, is used to describe such events.

A number of event shape variables (ESVs) [1, 2] have been defined to reflect the different aspects of QCD over a wide range of energy as well as theoretical complexity. Their measurements[3, 4, 5, 6, 7, 8, 9, 10] at the colliders have been used to improve our understanding of QCD - parton distribution function (PDF), multi-parton interaction (MPI), initial state radiation(ISR), final state radiation (FSR), fragmentation, hadronization, etc. Although data from LEP, Tevatron and RunI of the Large Hadron Collider(LHC) have been used to this end, RunII data of LHC has so far not been used exhaustively. This study attempts to improve the understanding of QCD using four ESVs measured in RunII data [3].

Monte Carlo (MC) simulation programs such as pythia8 [11], herwig++ [12], MadGraph5_amc@nlo [13], etc. combine analytical calculations and approximate models to describe the collisions. pythia8 is a popular MC event generator and is extensively used by the particle physics community. Event generators which use matrix element (ME) calculations for the hard scattering process often use pythia8 for emulating the subsequent fragmentation and hadronization process. Data from colliders have been used to optimize pythia [14, 15, 16]and it is imperative that data from the RunII of LHC is used to optimize it.

The plan of the paper is as follows. Section 2 briefly reviews the theory relevant to parton shower. Section 3 examines the agreement of pythia8 with measurement of ESVs by the CMS experiment [3] in pp collisions at $\mathrm{\surd{s}=13}$ TeV. Section 4 describes the studies of the different parameters of pythia8 and optimization of the parameter set. It also describes the application of the improved parameter set to compare with results of other analyses of CMS [17] and ATLAS [18]. Section 5 contains the summary and outlook.

2 Parton Shower description of hadron collisions

In high energy collisions, it is intensely complicated to achieve an exclusive picture of the final state partons using matrix elements calculation and only fixed order treatment is not sufficient. However, for comparison with experimental analyses, a fully exclusive description of the final states using calculations based on the shower evolution and hadronization is more suitable. Such methods are described through phenomenological models embedded in the shower MC codes. MC event generators such as pythia8, herwig++ and sherpa use their own parton shower treatment which have their own merits. herwig++ uses the coherent branching algorithm [19] for angular order branching whereas Catani-Seymour dipole factorisation [20] is used in sherpa.

pythia8 uses leading order(LO) calculations for generating the $2\rightarrow 2$ hard scattering processes. It uses ‘transverse momentum’ ordered parton shower[21] with $p^{2}_{\perp}$ as evolution variable for the generation of $2\rightarrow n$ ($n\geq 2$) final states by taking account ISR and FSR. pythia8 also emulates MPI and evolution of the partons into hadrons. In this showering scheme, resummation is done to all orders with a certain logarithmic accuracy.

For the splitting of a parton $a\rightarrow bc$, pythia8 uses the branching probability expressed by the DGLAP evolution equations:

[TABLE]

where $\mathrm{P_{a\rightarrow bc}}$ is the DGLAP splitting function and $p^{2}_{\perp}$ represents the scale of the branching; $z$ represents the sharing of $p_{\perp}$ of $a$ between the two daughters, with $b$ taking a fraction $z$ and $c$ the rest, $1-z$. Here the summation goes over all allowed branching, e.g. $q\rightarrow qg$ and $\mathrm{q\rightarrow q\gamma}$ and so on. Now, these probabilities become larger than unity due to divergence when $p^{2}_{\perp}\rightarrow 0$ which is taken care by introducing a term $\mathcal{P}^{no}_{a}(p^{2}_{{\perp}_{\textrm{max}}},p^{2}_{{\perp}_{\textrm{evol}}})$ known as Sudakov form factor [22]. This Sudakov factor ensures that there will be no emission between scale $p^{2}_{{\perp}_{\textrm{max}}}$ to a given $p^{2}_{{\perp}_{\textrm{evol}}}$.

Considering lightcone kinematics, evolution variables for $a\rightarrow bc$ at virtuality scale $Q^{2}$ for space-like branching (ISR) and time-like branching (FSR) are given by equations 2.2 and 2.3, respectively.

[TABLE]

Finally, equations 2.4 and 2.5 describe the evolutions for ISR and FSR respectively [21].

[TABLE]

Currently both the running re-normalisation and factorisation shower scales, i.e. the scales at which $\alpha_{s}$ and the PDFs are evaluated, are chosen to be $p^{2}_{\perp_{\textrm{evol}}}$ [23]. The general methodology of pythia8 for ISR, FSR, MPI is to start from some maximum scale $p^{2}_{{\perp}_{\textrm{max}}}$ and evolve downward in energies towards next branching untill the daughter partons reach some cut-off.

3 Event Shape Variables at RunII of LHC

Event shape variables[1] are defined in terms of the four-momenta of the objects with multi-jet final states. They are sensitive to the topology of the primary hard scattering and also, the evolution of the primary partons into stable hadrons. With a judicious choice of the ESVs it is possible to use experimental data to confront QCD predictions - from perturbative analytical calculation in the early phase to the nonperturbative models in the later phase of the event. These variables also have the advantage that, they are safe from collinear and infrared(IRC) divergences [2]. The ESVs may also be used to look for signature of physics beyond the Standard Model [24, 25]. LHC experiments ATLAS, CMS, and ALICE have studied various ESVs, evaluating them with jets and charged particles in proton-proton (pp) collision [3, 4, 5, 6, 7, 8, 9, 10].

Most recently the CMS has evaluated four ESVs in multi-jet events (with atleast three jets) in pp collisions at $\surd{s}=13$ TeV [3]. These are - the complement of transverse thrust ($\mathrm{\tau_{\perp}}$), total jet mass ($\mathrm{\rho_{\text{Tot}}}$), total transverse jet mass ($\mathrm{\rho^{\textrm{T}}_{\text{Tot}}}$) and total jet broadening ($\mathrm{B_{\textrm{T}}}$) and defined as ratios of momenta of the jets in an event and hence many uncertainties cancel out. These ESVs have higher values for multijet, spherical events and lower values for two-jet, pencil like events. The complement of transverse thrust ($\mathrm{\tau_{\perp}}$) is expected to be sensitive to the initial hard scattering, multi-parton interaction, and emissions of high $p_{\textrm{T}}$ gluons in the form of ISR and FSR. The other three - $\mathrm{\rho_{\text{Tot}}}$, $\mathrm{\rho^{\textrm{T}}_{\text{Tot}}}$ and $\mathrm{B_{\textrm{T}}}$ - have additional dependence on the details of fragmentation and hadronization.

The distributions of the ESVs observed in data are unfolded to remove the effects of the efficiencies and acceptances of the CMS detector and a broad agreement between data and the predictions of pythia8, herwig++ and MadGraph5_amc@nlo is observed. In case of pythia8, it is noted that the agreement is better for the $\mathrm{\tau_{\perp}}$ and $\mathrm{\rho^{\textrm{T}}_{\text{Tot}}}$ which are computed using $p_{\textrm{T}}$ of the jets. However, in case of $\mathrm{B_{\textrm{T}}}$ and $\mathrm{\rho_{\text{Tot}}}$, which are computed using $\vec{p}$ of the jets, pythia8 overestimates the multijet nature of the events. This indicates that the flow of energy in the transverse plane is better modelled by pythia8 while the overall three-dimensional modelling is not adequate.

CMS also observes [3] that ISR has a large effect on the distributions of the ESVs and increases the spherical nature of the multijet events. The effect of FSR is much smaller and that of MPI is negligible.

For all the four ESVs, the overall agreement improves with the energy scale of the event. The strong coupling $\mathrm{\alpha_{s}}$ is smaller for higher energy and reduces the emission of hard gluons in the early stage of the evolution of the partons. This makes the higher order calculations less important for describing such event. CMS has used average $p_{\textrm{T}}$ of the two leading jets of an event, $H_{\textrm{T,2}}=(p_{\textrm{T,\, jet1}}+p_{\textrm{T,\, jet2}})/2$, to represent its energy scale.

4 Optimizing the Parton Shower Model of PYTHIA8

Recently, both CMS and ATLAS have done several optimizations of pythia8 around its Monash tune [14]. CMS has used underlying event (UE) and charged particle data from the CDF ($\mathrm{p\bar{p}}$ collisions, Fermilab) and CMS (pp collisions, LHC) for tuning the strong coupling and MPI related parameters for different PDF sets [26, 27]. ATLAS has optimized ISR, FSR and MPI related parameters using a number of observables [28]. The present study also uses the Monash tune, pythia v8.235 with NNPDF2.3 PDF (LO) set, as the default. For each parameter, variation has been done around the default value of the Monash tune.

Since Monash tune overestimates the multijet nature of the events [3], the role of ISR and FSR need to be examined along with the choice of a suitable value of $\mathrm{\alpha_{s}}$. pythia8 has the provision to use separate values of $\mathrm{\alpha_{s}(M_{Z})}$ for ISR and FSR and the corresponding parameters are SpaceShower:alphaSvalue and TimeShower:alphaSvalue respectively. These will be referred to as $\mathrm{\alpha_{s}^{\textrm{ISR}}}$ and $\mathrm{\alpha_{s}^{\textrm{FSR}}}$. The Monash tune uses $\mathrm{\alpha_{s}(M_{Z})=0.1365}$ for both. The maximum evolution scale involved in the showering is set to match the scale of the hard process itself. In pythia8 it set equal to the factorization scale, but allows its modification by multiplicative factors SpaceShower:PTmaxFudge for ISR and TimeShower:PTmaxFudge for FSR. Both parameters in the Monash tune have 1 as their default values.

For each point in the parameter space, $10^{6}$ events have been generated to ensure that the statistical uncertainty is smaller than the experimental uncertainty. The resulting distributions have been compared with data and $\mathrm{\chi^{2}/NDF}$ has been calculated to check the goodness of fit. The optimization has been done using professor v2.3.0 [29] along with RIVET v2.6 [30] with the complete set of ESVs distributions from pythia8 as available from [3]. Post optimization, the new parameter set is checked using other relevant results from the CMS [17] and ATLAS [18].

4.1 Optimization of the Strong Coupling for ISR and FSR

The Monash tune of pythia8 uses the same value of $\mathrm{\alpha_{s}}$ for both ISR and FSR as the default. It is a choice of convenience but not mandatory. It has been argued that the effective scales for ISR and FSR are different which could be absorbed by effective values of the corresponding strong coupling and this is acceptable. This freedom to set nonidentical $\mathrm{\alpha_{s}}$ for ISR and FSR within the validity of the PS approach gives MC tuning studies a space to best describe the observables. The observation of CMS [3] that ISR has a rather large effect on the ESVs compared to FSR adds to the importance of this issue. So it is instructive and important to probe if pythia8 prefers to have different values of $\mathrm{\alpha_{s}}$ for both ISR and FSR.

First,  $\mathrm{\alpha_{s}^{\textrm{ISR}}}$ is varied, keeping  $\mathrm{\alpha_{s}^{\textrm{FSR}}}$ and other parameters fixed at their default values. The value of $\alpha_{s}^{\textrm{ISR}}$ is varied by 20% about the default value in steps of 2%. The results, in terms of $\mathrm{\chi^{2}/NDF}$ are shown in figure 1. It is seen that $\mathrm{\tau_{\perp}}$ and $\mathrm{\rho^{\textrm{T}}_{\text{Tot}}}$ are not sensitive to $\alpha_{s}^{\textrm{ISR}}$ while $\mathrm{\rho_{\text{Tot}}}$ and $\mathrm{B_{\textrm{T}}}$ is mildly sensitive preferring a lower value.

Next, $\mathrm{\alpha_{s}^{\textrm{FSR}}}$ is varied, keeping other parameters including $\mathrm{\alpha_{s}^{\textrm{ISR}}}$ fixed at their default values. The value of $\mathrm{\alpha_{s}^{\textrm{FSR}}}$ is also varied by 20% in steps of 2%. The results, in terms of $\mathrm{\chi^{2}/NDF}$ are shown in figure 2. It is seen that $\mathrm{\tau_{\perp}}$ is not sensitive to $\mathrm{\alpha_{s}^{\textrm{FSR}}}$, $\mathrm{\rho^{\textrm{T}}_{\text{Tot}}}$ prefers the default value, $\mathrm{\rho_{\text{Tot}}}$, and $\mathrm{B_{\textrm{T}}}$ prefer a higher value.

The effect of variation of either $\alpha_{s}^{\textrm{ISR}}$, or, $\alpha_{s}^{\textrm{FSR}}$ on the four ESVs conform to the observations of CMS regarding the effects of ISR and FSR. The disagreements for the first two $H_{\textrm{T,2}}$ ranges for all of them are very large and appear to be beyond the scope of parameter tuning and left out of the next study, where both $\alpha_{s}^{\textrm{ISR}}$ and $\alpha_{s}^{\textrm{FSR}}$ are varied simultaneously and over the same ranges as above. Each plot in figure 3 shows the summed up $\mathrm{\chi^{2}/NDF}$ for the six $H_{\textrm{T,2}}$ ranges for the same variable. This gives us an overall idea of the preferable values of $\alpha_{s}^{\textrm{ISR}}$ and $\alpha_{s}^{\textrm{FSR}}$.

4.2 Optimization of Maximum scale of the Shower Evolution

The sensitivity of the ESVs to maximum scale of the ISR shower evolution SpaceShower:PTmaxFudge ($\mathrm{\tt PTmaxFudge^{ISR}}$) has been checked. Figure 4 shows the variation of $\mathrm{\chi^{2}/NDF}$ for the four ESVs with this parameter. It is observed that $\mathrm{\tau_{\perp}}$ and $\mathrm{\rho^{\textrm{T}}_{\text{Tot}}}$ have similar dependence on $\mathrm{\tt PTmaxFudge^{ISR}}$, preferring a value of $\mathrm{\sim 1}$. On the other hand, $\mathrm{\rho_{\text{Tot}}}$ and $\mathrm{B_{\textrm{T}}}$ prefer a value $\mathrm{\sim 0.6}$. Similar study for TimeShower:PTmaxFudge, related to FSR, shows no significant variation for the ESVs.

4.3 Simultaneous Optimization of the Strong Coupling and Maximum scale of the Shower Evolution

Finally, professor v2.3.0 [29] is used to simultaneously optimize $\mathrm{\alpha_{s}^{\textrm{ISR}}}$, $\mathrm{\alpha_{s}^{\textrm{FSR}}}$, and $\mathrm{\tt PTmaxFudge^{ISR}}$ keeping other parameters fixed to their default values in the Monash tune. CMS data for the four ESVs are made available through RIVET v2.6 [30] and HEP Data [31]. Overall, 120 different combinations of the three parameters are randomly sampled by professor in the ranges mentioned in table 1. professor was instructed to perform a third order polynomial fit to optimize parameter response for the observables, the four ESVs in this case. This goodness-of-fit minimization process gives in return best favourable values of the three parameters as listed in table 1. The uncertainty for the optimized parameters are obtained from Minuit [32] package interfaced with professor.

The values obtained for $\mathrm{\alpha_{s}^{\textrm{ISR}}}$, $\mathrm{\alpha_{s}^{\textrm{FSR}}}$, and $\mathrm{\tt PTmaxFudge^{ISR}}$ are 0.11409, 0.15052 and 0.9323 respectively, see table 1. It is noted that the model of QCD implemented in pythia prefers a lower value of $\mathrm{\alpha_{s}}$ for ISR compared to the default value used in the Monash tune but prefers a higher value of the same for FSR. It is also observed that the optimized values obtained here for $\mathrm{\alpha_{s}^{\textrm{ISR}}}$ and $\mathrm{\alpha_{s}^{\textrm{FSR}}}$ are very close to their values obtained while $\mathrm{\tt PTmaxFudge^{ISR}}$ is fixed to 1.

4.4 Validation of the Optimized Set of Parameters of PYTHIA8

The optimized values of the three parameters (see, table 1) are used to calculate the ESVs. It is seen that the agreement with data deteriorates slightly for $\mathrm{\tau_{\perp}}$ (figure 5) and $\mathrm{\rho^{\textrm{T}}_{\text{Tot}}}$ (figure 8) compared to the good agreement with the Monash tune. But, there is significant improvement in agreement with data for $\mathrm{\rho_{\text{Tot}}}$ (figure 6) and $\mathrm{B_{\textrm{T}}}$ (figure 7) compared to the Monash tune. Since $\mathrm{\rho_{\text{Tot}}}$ and $\mathrm{B_{\textrm{T}}}$ had a rather poor agreement between data and the Monash tune, overall this new set of parameters is better.

It is imperative that the optimized set of parameters is tested against some other QCD results at the LHC. Inclusive jet cross-section is a very important QCD measurement at the LHC and is sensitive to PDF of protons and $\mathrm{\alpha_{s}}$ and both CMS [17] and ATLAS [18] have studied this with the 13 TeV data. The CMS has measurements of inclusive cross-sections for $\textrm{anti-}k_{\textrm{T}}$ jets with R = 0.4, 0.7. Figures 9 and 10 show that the new parameter set improves the agreement with data compared to the Monash tune. Similar improvement is seen for the ATLAS measurement of $\textrm{anti-}k_{\textrm{T}}$ jets with R = 0.4 (figure 11).

5 Summary

The strong coupling and maximum shower evolution scale used in the parton shower model of pythia8 have been optimized using the CMS measurement of four event shape variables over a wide range of energy scale of the events. pythia8 with the optimized parameters shows better agreement with inclusive jet cross-section measurements by CMS and ATLAS experiments. This study suggests that the models of initial and final state radiations in pythia8 can be improved to better represent various Quantum Chromodynamics related studies in the context of the Large Hadron Collider.

Acknowledgments

The authors sincerely thank Prof. Peter Z. Skands (Monash University) for useful discussions and suggestions. Financial support by the Department of Science & Technology, Government of India, through research grants, SR/MF/PS-02/2014-VB(G) and SR/MF/PS-03/2014-VB(G) is thankfully acknowledged.