# A Comparison of Deep Inelastic Scattering Monte Carlo Event Generators to HERA Data

## 1 Introduction

Monte Carlo generators are an essential tool in modern day experimental High Energy Physics. They play a crucial rôle in the analysis of the data, often in assessing the systematic errors of a measurement. For that reason it is of great importance that the Monte Carlo programs give results that agree closely with the experimental data. This paper aims to describe the agreement, deficiencies and tuning of the Monte Carlo models with the neutral current deep inelastic scattering (DIS) data at HERA. Extensive use is made of the utility package HzTool , which is a FORTRAN library containing a collection of experimental results from the H1 and ZEUS collaborations.

The work described here is part of an ongoing program. During the workshop a forum was established between the H1 and ZEUS collaborations for a joint coordinated investigation of the generators working closely with the programs’ authors.

## 2 Monte Carlo Models

The ARIADNE , HERWIG  and LEPTO  Monte Carlo generators for DIS data have been investigated during the course of the workshop. Other programs such as RAPGAP  and those developed over the duration of the workshop will be examined as part of the ongoing program of work. In the following sections a brief introduction to each of the three generators studied is given.

### 2.1 ARIADNE

In ARIADNE the QCD cascade is modelled by emitting gluons from a chain of independently radiating dipoles spanning colour connected partons , correcting the first emission to reproduce the first order matrix elements . The hadronisation of the partons into final state particles is performed by the Lund string model  as incorporated in JETSET . Since the proton remnant at one endpoint of the parton chain is treated as an extended object, the coherence condition allows only a fraction of this source to be involved in gluon radiation. Since the photon probing the proton only resolves the struck quark to a distance $`\lambda 1/Q,`$ the struck quark is also treated as an extended object. As a consequence gluon emissions in the proton and photon directions are suppressed. This phase space restriction is governed by $`a=(\mu /k_T)^\alpha `$ where $`k_T`$ is the transverse momentum of the emission, $`a`$ is the fraction of the colour antenna involved in the radiation, $`\mu `$ is a parameter related to a typical inverse size of a hadron and $`\alpha `$ governs the distribution of the energy along the dipole.

In the default version of ARIADNE, the mechanism for soft suppression of radiation due to the extended source of the proton remnant results in a suppression of radiation in the current region of the Breit frame at high $`Q^2.`$ In the course of the workshop a high $`Q^2`$ modification was developed  where this suppression in the current region was removed.

### 2.2 HERWIG

HERWIG relies on a coherent parton branching algorithm with additional first order matrix element corrections  to populate the extremities of phase space which the partons from the conventional QCD cascade fail to occupy. The partons are transformed into hadrons using the cluster fragmentation model , whereby the primary hadrons are produced from an isotropic two body decay of colour-singlet clusters formed from partonic constituents.

Since the Monte Carlo tuning to HERA data at the ‘Future Physics at HERA workshop’ , a new version of HERWIG (version 5.9) has become available. This version includes the modified remnant treatment of version 5.8d whereby the fragmentation of the cluster containing the hadronic remnant is treated differently to that containing the perturbative parton from the incident hadron. In addition, the particle decay tables have been updated and now contain a large amount of information on additional resonance decays.

The default version of HERWIG implements the next-to-leading order (NLO) running of the QCD coupling constant $`\alpha _s`$. The HERWIG philosophy is to incorporate as much perturbative QCD behaviour as possible, so even though the generator only uses a leading order (LO) parton shower cascade, a NLO $`\alpha _s`$ behaviour is implemented. This can be justified, to some degree, because of the HERWIG implementation of angular ordering in the QCD cascades. The H1 collaboration have modified HERWIG to allow a LO behaviour of $`\alpha _s`$ .

### 2.3 LEPTO

In LEPTO the hard parton processes are described by a leading order matrix element (ME). The soft and collinear divergences are regulated with a lower and upper cut in $`z_p`$ where $`z_p=pj_1/pq`$ where $`p`$ (q) is the proton (photon) four-vector and $`j_1`$ the four-vector of one of the partons produced in the hard subprocess. In addition, the invariant mass squared of the two hard partons is required to exceed a minimal value, $`\widehat{s}_{min}`$. Below the ME cut-offs, parton emissions are treated by parton showers based on the DGLAP evolution equations . The amount of parton radiation depends on the virtuality chosen between a lower cut-off ($`Q_0^2`$) and a maximum given by the scale of the hard process or the ME cut-off. LEPTO uses JETSET for the hadronization of the partons. In addition to this non-perturbative phase, LEPTO introduces another non-perturbative mechanism. This is a soft (i.e. at a scale below $`Q_0^2`$) colour interaction which assumes that the colour configuration of the partonic system can be changed whilst traversing the colour field of the proton remnant. This was introduced in order to reproduce the rapidity gap events observed at HERA. 

During the course of the workshop a new version of LEPTO was released (version 6.5.2$`\beta `$). This version introduced a new scheme for dealing with SCI events, in which the probability of soft colour interactions is suppressed depending on the difference in area spanned by the possible string configurations (after or before a soft colour interaction) . This means at high $`Q^2`$ there are effectively no soft colour interactions.

## 3 Model comparison with the data

### 3.1 ARIADNE

These studies closely followed those of the previous tuning exercise  performed at the ‘Future Physics at HERA’ workshop. However, they have been extended: jet data were now available for inclusion in the comparisons; and, the behaviour of the parameter PARA(25) was considered for the first time. ARIADNE version 4.10 has been investigated including the modified treatment of high $`Q^2`$ DIS events (MHAR(151) = 2.)

The four model parameters considered are listed in Table 1, which includes a short description of their influence. The two parameters PARA(10) and PARA(15) govern the slope of the suppression line (in the phase space available for gluon emission) for the proton and the struck quark, respectively. PARA(25) governs the probability of emissions outside the soft suppression cut-off, while PARA(27) corresponds to the square root of the primordial $`k_T^2`$ in the proton.

#### 3.1.1 Approach 1

This approach concentrated on HERA data at $`Q^2>80\mathrm{GeV}^2.`$ The motivation for this was to minimize the theoretical uncertainties in the generator associated with parton evolution in the low $`(x,Q^2)`$ region. The distributions that were most sensitive to the parameters under investigation were first ascertained. Next a combined $`\chi ^2`$ was calculated for each parameter setting according to

$$\chi _{Comb}^2=\frac{1}{nsets}\underset{i=1}{\overset{nsets}{}}\chi _i^2$$
(1)   
where $`\chi _i^2`$ represents the total (average) $`\chi ^2`$ per degree of freedom (d.o.f) of data set $`i`$. The parameter combination that yields the minimum of the overall $`\chi _{Comb}^2`$ corresponds to the tuned result.

The following distributions for $`Q^2>80\mathrm{GeV}^2`$ were investigated:

* scaled momentum $`x_p`$ distributions in the current region of the Breit frame ;
* flows of transverse energy in the hadronic centre of mass system ;
* differential distributions and evolution of mean of event shape variables thrust $`T_c`$ and $`T_z`$, jet broadening $`B_c`$ and jet mass $`\rho _c`$ in the current region of the Breit frame ;
* fragmentation functions and charged particle multiplicities in the current region of the Breit frame ;
* differential and integrated jet shapes as a function of pseudorapidity $`\eta `$ and transverse energy $`E_T`$ ;
* (2+1) jet event rate as a function of the transferred momentum squared,$`Q^2`$ .

Table 2 summarises the total $`\chi ^2`$ for each of the six sets of data, together with the combined $`\chi _{Comb}^2`$ given by equation 1. The results of tuning the parameters, set 1 Table 1, agree very well with those previously obtained . Both transverse energy flows and the jet data strongly favour high values for PARA(10), in contrast to the lower value favoured by the charged particle momentum distributions . Again the value PARA(15)=0.5 yields a better overall description of data, as compared to its default value of 1.0. The jet variables are particularly sensitive to PARA(15). The behaviour of PARA(25) was studied for the first time. Although its influence is not large in general, it clearly has a significant effect on the predictions related to the (2+1) jet event rate and on the transverse energy flows. The results suggest a lower value compared to the default one. Parameter PARA(27) is a relatively insensitive parameter, but the data disfavour high values, such as 0.8-1.0 GeV, when describing transverse energy flows and jet shapes. The tuning has been performed using the GRV94 parton density function . The use of the parton density function CTEQ4M  results in a slightly worse value for $`\chi _{Comb}^2`$.

The improvement achieved with the tuned parameters is mostly due to a better description of jet shapes and (2+1) jet event rates (see Figs. 1 and 2). This improved agreement with the jet data leads to a slightly worse description of other distributions such as fragmentation functions and transverse energy flows (see Figs. 3 and 4). Conversely, the new treatment at high $`Q^2`$ of ARIADNE describes much better than before transverse energy flows and event shape variables, but lessens the agreement with the data on jet shapes. The current study seems to suggest that a simultaneous description of jet and charged particle distributions is difficult.

#### 3.1.2 Approach 2

The 2nd approach has additionally investigated the behaviour of ARIADNE for $`Q^2<80\mathrm{GeV}^2`$ and has also concentrated on different data sets than those used in approach 1. In particular new preliminary data from H1 on dijet production  has been used along with charged particle distributions in $`\gamma ^{}P`$ centre of mass frame . In addition, the $`E_T`$ flows and the charged particle distributions in the Breit frame have been considered but, again, at lower $`Q^2`$ values than approach 1. ARIADNE 4.10, with the high $`Q^2`$ modifications, has been studied using CTEQ4L for the parton density parametrisation.

Investigations showed that parameter PARA(10) was very sensitive to the dijet cross section, especially at low $`Q^2,`$ and was also sensitive, to a lesser degree, to the $`E_T`$ flows. Parameter PARA(25) also influenced the agreement with the dijet measurement and to the rapidity distribution of hard $`p_T`$ charged particles but otherwise displayed little sensitivity to the data. The other two parameters, PARA(15) and PARA(27), displayed a lesser sensitivity to the data, though the hard $`p_T`$ particles and the dijet distributions proved the most affected to changes in these parameters.

Figure 5 shows the sensitivity of the dijet cross section to PARA(10). The default ARIADNE produces $`E_T`$ spectra for the dijets that are too hard, with the discrepancy predominantly occurring in the forward region, $`\eta _{fwd,lab}>1.0.`$ At low $`Q^2`$ there is a large variation in the $`E_T`$ spectrum but PARA(10) affects the distribution at both low and high $`E_T`$, which results in this parameter alone not being able to describe the complete $`E_T`$ spectra. This problem can be circumvented by varying PARA(25) in conjunction with PARA(10). Variation in PARA(25) alone gives larger fractional changes in the cross section at large $`E_T`$ than at smaller values of $`E_T`$, see Fig. 6.

The influence of PARA(10) on the $`E_T`$ flows can be seen in Fig. 7. The increase of this parameter suppresses $`E_T`$ production across the whole $`\eta `$ range. A similar effect is seen in the charged particle rapidity distribution, particularly for particles with $`p_T>1\mathrm{GeV}.`$ As can be seen from the $`x_p`$ spectra, Fig. 7, the current region of the Breit frame seems relatively insensitive to this parameter. The $`E_T`$ flows are less sensitive to PARA(25) than PARA(10), see Fig. 8. However the data seem to prefer values of PARA(25) smaller than the default. This preference is also true for the charged particle rapidity distributions regardless of any $`p_T`$ selection, for the default value of PARA(10).

The average $`\chi ^2`$ for the low and high $`Q^2`$ region, as well as the combined $`Q^2`$ regions, are shown in Table 3 at various settings of PARA(10) and PARA(25). An improved fit to the data was found for all distributions for the parameters listed as set 2 in table 1. It should be noted though that a comparison with the ZEUS jet shapes was not included in the data sets investigated in this approach.

### 3.2 HERWIG

HERWIG overall has fewer tunable parameters than the Lund family of generators . In particular the cluster fragmentation model has far fewer tunable parameters than the Lund string model. Many of the parameters are well constrained by e<sup>+</sup>e<sup>-</sup> annihilation data. Consequently, those involved with the hard subprocess and the perturbative QCD evolution of the final state parton shower were not varied for this study. It was found previously , that of the remaining parameters only a small number were seen to have any sensitivity to the distributions under study in DIS. Therefore it was decided to limit this study primarily to the effects of the CLMAX and PSPLT parameters, where CLMAX relates to the maximum allowed cluster mass and PSPLT is the exponent in generating the mass distribution of split clusters.

#### 3.2.1 Approach 1

The data studied in this approach corresponds to the same data sets considered in Approach 1 for ARIADNE but also extended to lower $`Q^2.`$ In addition to the parameters CLMAX and PSPLT, this study investigated the parameter DECWT, which provides the relative weight between decuplet and octet baryon production relevant to the new decay tables. The dependence on the parton density parameterisation (pdf) of the proton has also been investigated by studying CTEQ4L  and MRSD-  pdfs. Even though the MRSD- has in principle been retracted by the authors and is known to be too high at low $`x`$ it was used here to provide a more significant variation of the underlying distribution.

The three parameters were studied over the following ranges:

$`2.0<`$ $`\mathrm{CLMAX}`$ $`<5.0`$   
$`0.6<`$ $`\mathrm{PSPLT}`$ $`<0.9`$   
$`0.6<`$ $`\mathrm{DECWT}`$ $`<0.8.`$   
The effect of increasing CLMAX is to increase the $`E_T`$ flow as does increasing PSPLT. Increasing the $`E_T`$ flow with CLMAX has the effect of broadening the jet shapes and producing harder momenta spectra for the charged particles. This is thought to be due to the fact that the clusters are allowed to have more energy before they are forced to split. Reducing DECWT increases the $`E_T`$ flow predictions at low values of PSPLT with a smaller reduction or slight increase for larger values of PSPLT.

An attempt was made to tune the standard HERWIG (using MRSD- parton density functions) and compare with tuned values from LEP data from L3. The best values of the parameters achieved for the HERA data are listed in Table 4. Neither the ‘tuned’ set 1 parameters or the L3 parameters can describe the transverse energy flows at low $`x`$ and $`Q^2`$, see Fig. 9, whilst at higher $`Q^2`$ and $`x`$ the ‘tuned’ values give a better description of the HERA data than the L3 values. The jet shape distributions also prefer the ‘tuned’ values. With the parameters chosen for investigation it was not possible to achieve a consistent description of the data at both low and high $`(x,Q^2).`$

In an attempt to overcome the difficulty in obtaining sufficient $`E_T`$ at higher $`Q^2`$ without using very high values of CLMAX, an investigation of HERWIG with LO running $`\alpha _s`$ was made. Two sets of parameter settings are shown in Table 4 for this modified HERWIG, in conjunction with MRSD-. Set 2 gives the best description of the $`E_T`$ flows, whilst conversely set 3 gives a better description of the jet shape data. Figure 10 compares the HERWIG model predictions for the $`E_T`$ flows with the data. Set 2 describes these distributions well over the whole $`x`$ and $`Q^2`$ range. Set 3 also improves the description of this data in the highest $`Q^2`$ bins, though it underestimates the data in the lowest bins. Figure 11 compares the HERWIG predictions, with the parameter sets, to the jet shape data. Set 3 gives a better description of this data than using the set 2 parameters. Set 2 predicts jets broader than that observed in the data and is in poor agreement with the data.

Investigation of the sensitivity of the data (and the subsequent parameter settings) to the choice of parton density functions was made in the modified HERWIG for the MRSD- and CTEQ4L parametrisations. In particular the $`E_T`$ flows and the jet shapes were sensitive to the choice of parton densities. A 4th parameter set was found using the CTEQ4L parton densities. Again a consistent description of both the $`E_T`$ flows and the jet shapes was not possible. The parameter set listed in table 5 gave a better description of the $`E_T`$ flows than the jet shapes. The $`\chi ^2`$ achieved for the various parameter settings of HERWIG using both MRSD- and CTEQ4L are given in Table 5 for the $`x_p`$ distribution in the current region and the $`E_T`$ flows.

#### 3.2.2 Approach 2

This study considered the same data samples as approach 2 for ARIADNE. Only the H1 modified HERWIG, with the running of $`\alpha _s`$ at leading order, has been considered. The parton densities used in this approach correspond to CTEQ4L.

At high $`Q^2`$ a reasonable description of the dijet data by HERWIG could be obtained only if a larger (than default) value of $`\alpha _s`$ ($`\mathrm{\Lambda }=250\mathrm{M}\mathrm{e}\mathrm{V}`$) was used, Figure 12. At low $`Q^2`$ HERWIG was unable to achieve a good description of the dijet data. The dijet cross sections were relatively insensitive to changes in the hadronisation parameters.

In Figure 12 the DISENT<sup>1</sup><sup>1</sup>1The DISENT program incorporates a NLO calculation for DIS at the parton level. It can also be used to obtain partonic LO predictions.  $`𝒪(\alpha _s)`$ predictions (using $`Q^2`$ as the renormalisation scale) are compared to HERWIG. The HERWIG predictions are in agreement with the DISENT LO calculation. In the same figure, it is also shown that the NLO corrections (K-factors), in particular at low $`Q^2`$ and forward pseudorapidities $`\eta _{\mathrm{fwd},\mathrm{lab}}2,`$ are large. The parton showers used to emulate higher orders in HERWIG are insufficient to account for these large NLO corrections.

In contrast to the dijet data, the $`E_T`$ flows and the rapidity distributions of charged particles exhibit a strong dependence on the fragmentation parameters. As in approach 1, the current region of the Breit frame and the high $`Q^2`$ data prefer different settings of PSPLT and CLMAX parameters than does the low $`Q^2`$ data. The results are summarised in Table 6 and the HERWIG predictions are compared to the data in figure 13. The high $`Q^2`$ and the Breit frame current region data prefers settings of $`\mathrm{CLMAX}=3.0`$ (the default) and $`\mathrm{PSPLT}=1.2`$ whilst the low $`Q^2`$ data favour a higher value of $`\mathrm{CLMAX}=5.0`$ with a slightly lower value of $`\mathrm{PSPLT}=1.0.`$

Although variation of the fragmentation parameters leads to large changes in the prediction of the HERWIG model, the underlying parton dynamics in HERWIG are not sufficient to describe the HERA DIS data.

### 3.3 LEPTO

The new version (6.5.2$`\beta `$) of LEPTO was confronted with preliminary high statistics $`(2+1)`$ jet data from the H1 collaboration  (statistical error only on the data.) This data set consists of DIS events which are all forced to be of a $`(2+1)`$ jet configuration using the modified Durham algorithm. The distributions studied were $`y_2,`$ the cut-off in the algorithm where an event is first defined as $`(2+1)`$ and the angles in the laboratory frame of the forward and backward going jet ($`\theta _{\mathrm{fwd}}`$ and $`\theta _{\mathrm{bwd}}`$). In addition the jet variables $`x_{\mathrm{jet}}`$, defined as $`Q^2/(Q^2+\widehat{s})`$ where $`\widehat{s}`$ is the invariant mass of the jet(parton) pair, and $`z_p`$, defined as $`1/2(1\mathrm{cos}\theta ^{})`$ where $`\theta ^{}`$ is the polar angle of jet in the photon-parton centre of mass system, were investigated.

The following parameters, that control the cut off in the $`𝒪(\alpha _s)`$ matrix element in the generator, were found to have significant impact on the description of the data and have been studied:

* PARL(8) $`z_p^{min}`$ cut off, and
* PARL(9) $`\widehat{s}^{min}`$ cut off.

The new SCI scheme implemented in LEPTO version 6.5.2$`\beta `$ leads to a dramatic improvement in the description of the data compared with version 6.5. The $`\chi ^2`$ is typically reduced by a factor of 5–6, see Table 7. The predictions of LEPTO compared to the data are shown in Figures 14 and 15. Further significant improvement in the description of the data with LEPTO 6.5.2$`\beta `$ was achieved by optimizing the parameters PARL(8) and PARL(9). The results of this optimization are shown in Table 8 (set 1) and the improvement can clearly be seen in the comparison with the data in Figure 15; the corresponding $`\chi ^2`$ values are given in table 7. It should also be noted that LEPTO describes the $`Q^2`$ dependence of the jet distribution well.

A complementary way to optimize LEPTO for jet distributions, instead of applying hard cuts on divergences of the matrix element (set 1), is to loosen these cuts so that LEPTO is forced to find appropriate divergency cuts on an event–by–event basis. The preferred values of PARL(8) and PARL(9) using this approach are listed in Table 8 (set 2) and the corresponding $`\chi ^2`$ values in Table 7.

The variation of the intrinsic $`k_T`$, PARL(3), and the cut–off value of the initial–state parton shower, PYPAR(22), had no effect on the quality of the description of the jet data. Also the jet data were insensitive to the choice of the parton density functions.

Although both approaches to describing the data with LEPTO, via PARL(8) and PARL(9), result in significant improvements, no satisfactory description of the measured 2–jet distributions could be achieved. The parameter sets 1 and 2 were then cross checked against the data samples used in approach 1 for ARIADNE but over the whole $`Q^2`$ range, see Table 9. Besides the $`(2+1)`$ jet rate and the charged particle $`x_p`$ distribution in the current fragmentation region, the default version of LEPTO 6.5.2$`\beta `$ gave a better description of the data.

## 4 Summary

During the course of the workshop new versions of the LEPTO and ARIADNE Monte Carlo generators were made available. These modified versions of the generators were in far better agreement with data.

An attempt was made to find sets of parameters for the ARIADNE, LEPTO and HERWIG generators that would describe the DIS HERA data. It proved difficult to find such a parameter set that would describe the whole range of distributions at both low and high $`Q^2.`$ A number of parameter sets are given for each generator that are optimised for a particular region of phase space.

This paper attempts to summarise a ‘snapshot’ of an ongoing program of work between experimentalists of both the H1 and the ZEUS collaborations and the authors of the event generators. The ultimate aim is to have event generators that are able to describe the complex structure of DIS events at HERA as impressively as they do the LEP data.