# Parton shower Monte Carlo event generators

Post-publication activity

Curator: Bryan Webber Figure 1: Schematic depicting the stages of simulation of a hadron-hadron collisions by a parton shower Monte Carlo event generator.

## Monte Carlo method

The Monte Carlo method is fundamentally a technique for numerical integration using random numbers, or, in practice, numbers from a pseudorandom number generator. To predict the scattering cross section for a collider process, one must integrate the probability density, given by a quantum mechanical matrix element-squared, over the phase space of the process. For many processes of interest, the number of final-state particles is very large and so the phase space has many dimensions. However, the production of many particles can often be factored into a series of lower-dimensional integrations.

Schematically, we have a probability density $$f(x)$$ which is a function of an $$n$$-component vector $$x\ ,$$ related to the momenta of the particles involved, and we want to integrate it over some region in $$x$$-space of volume $$V$$, $I[f] = \int_V d^n x\,f(x)\;.$ Standard methods of integration (Simpson's, Gaussian, ...) are too laborious and/or inaccurate for $$n$$ large (say $$n>3$$). However, if $$N$$ points $$\{x_i,i=1,\ldots,N\}$$ are distributed randomly in $$V\ ,$$ then the central limit theorem of statistics tells us that the mean value of $$f$$ on those points is an unbiased estimator of the integral, $I[f]\simeq\langle f\rangle = \frac VN \sum_{i=1}^N f(x_i)\;,$ and that the estimated error $$E[f]$$ on this evaluation is given by the variance of $$f\ ,$$ $\mbox{Var}(f)=\langle (f-\langle f\rangle )^2\rangle= \langle f^2\rangle-\langle f\rangle ^2\;,$ as $E[f]= V\sqrt{\frac{\mbox{Var}(f)}{N-1}}\;.$ Thus the error decreases as the inverse square root of the number of points, independent of the dimensionality of the integral. Furthermore, for a given number of points, the error will be less if the variance of the integrand is small.

The variance can be reduced by a change of variables that "flattens" the integrand. Consider the mapping $$x\to y(x)$$ with Jacobian $\left|\frac{\partial y}{\partial x}\right| = g(x)\;.$ Then $I[f] = \int_{V'} d^n y\,\frac{f(x)}{g(x)}\;,$ where $$V'$$ is the region in $$y$$-space corresponding to $$V$$ in $$x$$-space. If $$h=f/g$$ is a function with less variance than $$f$$ itself then the error will be reduced by distributing points uniformly in $$y$$-space. For example, consider the following transformation: $I[f] = \int \frac{dx}{x} = \int dy,$ with $$y=\ln(x)\ .$$ This is known as importance sampling.

A further advantage of the Monte Carlo method is that we can use the points generated to represent simulated events, since each point corresponds to a set of momenta for the particles involved. These simulated events can be used to histogram the distributions of observables of interest or passed to a detector simulation to study the response of the detector.

To obtain a set of points distributed according to the predicted probability distribution $$f(x)\ ,$$ as desired for an event generator, we can apply the hit-or-miss method, accepting points with probability $$h/h_{lim}\ ,$$ where $$h_{\rm lim}$$ is an upper bound on the value of the flattened distribution $$h$$ in $$V'\ .$$ The Monte Carlo efficiency, as measured by the fraction of points accepted, $$\langle h\rangle/h_{\rm lim}\ ,$$ will usually also be increased as a result of the variance reduction.

The above procedure corresponds to the generation of unweighted events, which can be treated like real events when computing observables. If it is difficult to arrive at an acceptable efficiency by reducing the variance of the integrand and/or finding a good upper bound on it, one may wish to resort to generating weighted events. In that case the phase-space points $$\{x_i\}$$ (or $$\{y_i\}$$ if some variance reduction has been achieved) are used to represent events, but each event has a different weight $$f_i$$ (or $$h_i$$) when contributions to observables are computed. In that case one has to take account of the variance of the weights when computing error on observables. That means, for example, that one must keep track of the sum of the squared weights as well as the weights contributing to each bin of a histogram. In contrast the error for the unweighted events obtained from hit-or-miss is just given by the square root of the number of events in the bin.

## Hard subprocess

Event simulation normally begins with a relatively simple subprocess resulting from a highly energetic collision of constituents of colliding particle beams. For example, at the CERN Large Hadron Collider (LHC) a top quark-antiquark pair can be created in the collision of a pair of gluons or a light quark-antiquark pair from the incoming protons. The momenta of the colliding constituents are selected by sampling the parton distribution functions of the proton at the energy scale of the subprocess. These distributions have been measured at lower energies and in other processes and are evolved to higher scales using the QCD evolution equations for parton densities. Convolution with the differential cross section of the subprocess and integration over phase space gives the relevant production cross section (for top quark pair production in this case). Using the Monte Carlo method for the convolution and phase-space integrations means that we obtain at the same time the momentum distributions of the primary produced objects, in the form of a set of representative phase-space points.

Any of the incoming or outgoing objects of the hard subprocess that carry QCD colour will be eligible for parton showering in the next stage of the simulation. For this purpose, and for the subsequent hadronization, they need to carry information about the colour flow in the subprocess. This is normally assigned in the so-called large-$$N$$ approximation, where $$N$$ is the number of colours. In this approximation, any subprocess can be represented as a sum of distinct colour flows. For example, in the subprocess $$gg\to t\bar t$$ in this approximation, the colour of the top quark comes from one of the incoming gluons and the anticolour of the antiquark comes from the other. In general, if $$M_N$$ represents the matrix element-squared for $$N$$ colours, we have $M_\infty = \sum_i M_\infty^i$ where $$i$$ labels distinct colour flows. Then colour flow $$i$$ is assigned with probability corresponding to the matrix element-squared $$M^i$$ where $M^i= \frac{M_\infty^i}{M_\infty} M_3\;.$ This ensures that the sum over colour flows gives the exact $$N=3$$ probability distribution.

Depending on the level of precision required, the hard subprocess may be treated at the lowest relevant order of perturbation theory (leading order, LO), or to higher order. At present, event generators do not go beyond next-to-leading order, NLO.

## Parton showers

The hard subprocess, by definition, involves large momentum transfers and therefore the partons involved in it are violently accelerated. Just as accelerated electric charges emit QED radiation (photons), the accelerated coloured partons will emit QCD radiation in the form of gluons. Unlike the uncharged photons, the gluons themselves carry colour charges and can therefore emit further radiation, leading to parton showers. In principle, the showers represent higher-order corrections to the hard subprocess. However, it is not feasible to calculate these corrections exactly. Instead, an approximation scheme is used, in which the dominant contributions are included in each order. These dominant contributions are associated with collinear parton splitting or soft (low-energy) gluon emission.

We consider first the almost-collinear splitting of a parton of type $$i$$ into $$j+k\ ,$$ for example $$q\to q+g\ .$$ If the $$n$$-parton differential cross section before splitting is $$d\sigma_n\ ,$$ then after the splitting it becomes $\tag{1} d\sigma_{n+1} \approx d\sigma_n \frac{\alpha_S}{2\pi} \, \frac{d\theta^2}{\theta^2} \, dz\,d\phi\,P_{ji}(z,\phi)$

at the leading order in perturbation theory, where $$\alpha_S$$ is the strong coupling constant, $$\theta$$ and $$\phi$$ are the opening angle and azimuthal angle of the splitting, and $$P_{ji}$$ is the $$i\to j$$ splitting function, which describes the distribution of the fraction $$z$$ of the energy of $$i$$ carried by $$j\ .$$ This formula takes into account only the collinear-dominant contribution, i.e. the part that diverges as $$\theta\to 0\ ,$$ which factorizes in this way. After averaging over the azimuth $$\phi\ ,$$ the leading-order splitting functions are the same as those appearing in the QCD evolution equations for parton densities.

By sequential application of Eq. (1), using the Monte Carlo method to generate values of $$z,\,\theta$$ and $$\phi$$ for each splitting, a parton shower is developed from each coloured parton of the hard subprocess. An important feature of the showering algorithm, not apparent from Eq. (1), is the evolution variable. The simplest evolution variable to understand is the virtual mass-squared (the virtuality, $$q^2$$) of the partons in the shower. The splitting process $$i\to j+k$$ cannot take place with all the partons on their mass-shells. The dominant contributions will come from configurations in which the virtualities are strongly ordered, with the parton nearest to the hard subprocess farthest from its mass shell and the virtualities falling sharply as the shower evolves away from it. The upper limit on the initial virtuality is set by some momentum transfer scale $$Q$$of the hard subprocess, $$q^2<Q^2\ ,$$ and the shower is terminated when the virtualities have fallen to the hadronization scale, $$q^2= Q_0^2\sim 1\,\mathrm{GeV}^2\ .$$

It might appear that the parton shower approximation takes into account only the collinear-enhanced real parton emissions at each order in perturbation theory and neglects virtual (quantum loop) effects of the same order. However, this is not the case: such effects are included in the probability of not splitting during evolution from scale $$q_1^2$$ to $$q_2^2\ ,$$ which is given by the Sudakov form factor $\tag{2} \Delta_i(q_1^2,q_2^2) = \exp\left\{-\int_{q_2^2}^{q_1^2}\frac{dq^2}{q^2}\,\frac{\alpha_S}{2\pi} \, \int_{Q_0^2/q^2}^{1-Q_0^2/q^2} dz\int_0^{2\pi} d\phi \, P_{ji}(z,\phi)\right\}\,.$

Here we have used the fact that $$d\theta^2/\theta^2=dq^2/q^2$$ for any valid definition of the evolution variable. However, the limits of the $$z$$-integration depend on the precise definition of this quantity and of the evolution variable. They specify the range of $$z$$ in which the splitting is resolvable. An emission that would lie outside this range is too soft (i.e. with too little energy) or at too small an angle to be detected: it is declared unresolvable and is not included in the shower. The virtual and unresolvable contributions are combined in the Sudakov form factor (2). The situation is analogous to that in radioactive decay, where the probability of decay not occurring in a time interval $$t$$ is $$\exp(-\lambda t)$$ where $$\lambda$$ is the decay probability per unit time.

The generation of a parton shower thus proceeds as follows. Given the initial scale $$Q^2\ ,$$ one solves the equation $$\Delta_i(Q^2,q_1^2) = R_1\ ,$$ where $$R_1$$ is a random number uniform on the interval [0,1], for the scale $$q_1^2$$ of the first splitting. If $$q_1^2<Q_0^2\ ,$$ then the splitting was unresolvable and the showering of that parton is terminated. Otherwise, if the splitting was $$i\to j+k\ ,$$ then the procedure is repeated on parton $$j$$ using $$\Delta_j(q_1^2,q_2^2) = R_2$$ to determine the scale $$q_2^2$$ for the splitting of that parton, and similarly for parton $$k\ ,$$ and so on, until all attempted splittings have fallen below the resolvable scale $$Q_0^2\ .$$ At each splitting, the variables $$z$$ and $$\phi$$ are chosen according to the distribution $$P_{ji}(z,\phi)$$ using the Monte Carlo method, with $$z$$ in the resolvable region specified by the limits of integration in the Sudakov form factor.

### Final-state showers

A final-state shower is one that develops from an outgoing parton of the hard subprocess. In this case the evolution of the shower proceeds as described above: the primary parton starts at a high energy and a large time-like virtuality scale $$Q^2$$ set by the subprocess, and it loses energy and virtuality until it and all its descendant partons have fallen to the scale $$Q_0^2$$ at which splitting is terminated. At this point the final configuration of parton momenta can be passed to one of the hadronization models described below.

A complication that has to be addressed is the fact that the hard subprocess was generated with on-mass-shell partons, whereas any partons that initiated showers are now off mass-shell. This necessitates some momentum reshuffling to restore momentum conservation. For example, the outgoing parton three-momenta in the subprocess rest frame can be rescaled by a common factor and the showers boosted to those three-momenta.

### Initial-state showers

An initial-state shower is one that develops on an incoming parton of the hard subprocess, of which there are two at a hadron collider. Here a constituent parton from each of the incoming hadrons starts at a high energy and low virtuality and evolves to a higher space-like virtuality by emitting partons and losing energy. The showering of these partons terminates when they collide to initiate the hard subprocess, which sets the scale that limits the endpoint virtualities of the showers. For example, in the process $$q \bar q\to Z^0$$ the limit on the quark and antiquark virtualities is of the order of the $$Z^0$$ boson mass squared. As a result of showering the incoming partons, which started out in the beam directions, also acquire transverse momenta, and the vector sum of these is communicated to the hard subprocess. The partons emitted in the initial-state showers each initiate secondary showers that evolve in the same way as final-state showers.

The sequence of initial-state showering just described is not suitable for Monte Carlo event generation. If we want to simulate $$Z^0$$ production, for example, then the quark and antiquark at the end of their showering must have precisely the right 4-momenta to combine to form a system with the $$Z^0$$ mass. If their momenta at the start of showering are chosen according to the parton distribution functions at the initial low scale, there is very little probability of this and so the Monte Carlo efficiency will be very low. A better procedure, used by all the major event generators, is backward evolution. First the longitudinal momentum fractions $$x_1$$ and $$x_2$$ of the incoming partons of the hard subprocess are chosen using the parton distribution functions at the high hard-subprocess scale, subject to the kinematic constraint that $$x_1 x_2 S = s$$ where $$S$$ and $$s$$ are the collider and subprocess centre-of-mass energies-squared respectively. The incoming partons are then evolved backwards, gaining energy in each emission, to the low scale appropriate for constituents of the incoming hadrons. The virtualities and transverse momenta of the incoming partons of the hard subprocess follow from momentum conservation at the successive splittings in the showers. The only complication is that the no-splitting probability is no longer given by the Sudakov form factor (2) alone, but rather by that factor modified by a ratio of parton distribution functions: $\tag{3} \Delta'_i(q_1^2,q_2^2) = \Delta_i(q_1^2,q_2^2)\frac{f_i(x,q_2^2)}{f_i(x,q_1^2)}\,.$

An extra factor of $$f_j(x/z,q_2^2)$$ in the distribution of the momentum fraction $$z$$ for the splitting $$j\to i+k$$ cancels the corresponding term in the denominator of $$\Delta'_j(q_2^2,q_3^2)$$ for the next step of evolution and ensures that the parton distribution evolves correctly to the lower scale.

### Coherent showering

Many choices of the evolution variable are equivalent as far as collinear-enhanced contributions are concerned but they differ in the treatment of soft (low-energy) gluon emission, which is also enhanced. The preferred variable for this purpose is related to the opening angle of splitting, $$\theta\ ,$$ as written in Eq. (1), leading to an angular-ordered parton shower. This treats the coherence of soft radiation correctly and is therefore often called coherent showering. The reason may be seen most simply in QED. Consider the radiation of soft photons from an electron-positron pair moving apart at a small opening angle $$\theta\ .$$ At angles larger than $$\theta\ ,$$ the pair appears as an unresolvable neutral object and therefore there is no significant radiation. Another way to see this is to boost the dipole radiation pattern of the pair from its rest frame into a frame where the opening angle is $$\theta\ :$$ the radiation is essentially confined to two cones of opening angle $$\theta$$ around the directions of the two charges, corresponding to angular-ordered showering.

In QCD the same argument applies except that two colour charges may combine in different ways. For example a quark-antiquark pair may come from the splitting of a virtual photon, in which case they form a colour singlet and only radiate soft gluons inside the above cones, or else they come from a virtual gluon and form a colour octet. In that case when unresolved they appear as a single octet source and therefore the radiation outside the cones looks like that from a gluon. In the angular-ordered parton shower, this radiation is literally generated from the parent gluon before the splitting, although in reality it is the coherent emission from the pair. Similarly for the gluon splitting $$g\to gg$$ the coherent emission outside the cones of the produced gluons is generated as if it came from the parent gluon. Inside the cones the two colour charges are resolved and their subsequent splitting can be treated separately.

The main effect of angular ordering is to reduce the amount of soft gluon emission, relative to a virtuality-ordered shower. After hadronization, this is reflected in the number and energy distribution of hadrons produced in hard processes. For example, the energy distribution of hadrons produced in electron-positron annihilation assumes a Gaussian form in the variable $$\xi = \ln(\sqrt{S}/2E)$$ at high centre-of-mass energy $$\sqrt{S}\ ,$$ with a peak that moves as $$\xi_{\rm peak}\sim (\ln S)/4\ ,$$ corresponding to a suppression of low energies relative to the prediction neglecting coherence, $$\xi_{\rm peak}\sim (\ln S)/2\ .$$ A fuller discussion may be found in Ellis, Stirling and Webber (1996).

### Dipole showering

An interesting alternative treatment of parton showers, which also takes soft gluon coherence into account, is dipole showering. In this approach gluon emission is generated not by splitting a parton but according to the dipole radiation pattern of a pair of partons. As explained earlier, in the approximation of a large number of colours each quark or antiquark is uniquely connected to a colour partner and each gluon to two colour partners. Each pair of colour partners forms a dipole which, in this approximation, splits into two dipoles when it emits a gluon. The masses and momenta of the produced dipoles are determined by the kinematics of the gluon emission. These dipoles can split again, and so on until some minimum scale of dipole mass or relative transverse momentum is reached. The dipole radiation pattern can be adapted to take account of the different colour factors and collinear radiation distributions for quarks and gluons, thus correctly reproducing the leading soft and collinear enhanced terms to all orders in perturbation theory.

From the partonic viewpoint, dipole splitting is a $$2\to 3$$ process rather than $$1\to 2\ .$$ This has the advantage that momentum conservation can be satisfied with the partons on mass-shell at all stages of the shower. However, the question arises of how to share the recoil against the emitted gluon between the other two produced partons. Furthermore if one of those partons is a gluon then it also belongs to another dipole. This question can be resolved by requiring the parton distributions to satisfy the appropriate evolution equations, including enhanced higher-order terms. The recoils in initial-state showering must also be handled in a way that correctly generates the transverse momentum distribution of the hard subprocess.

The coupling 'constant' $$\alpha_S$$ appearing in Eqs. (1) and (2) is in reality a scale-dependent quantity, the QCD running coupling, which increases at low values of the shower evolution scale. Thus at some point in the evolution perturbation theory becomes invalid and the dynamics enter a non-perturbative phase, which leads to the formation of the observed final-state hadrons. This hadronization process is not amenable to the currently available non-perturbative techniques for calculation, and therefore event generators have to rely on models based on general features of QCD.

### String model

The string model for hadronization, depicted schematically in Figure 2, is based on the observation, from lattice simulations of QCD, that at large distances the potential energy of colour sources, such as a heavy quark-antiquark pair, increases linearly with their separation, corresponding to a distance-independent force of attraction. This is thought to be due to the self-attraction of the gluonic field, causing it to collapse into a string or tube configuration with thickness of the order of 1 fm when the separation of the sources becomes much larger than this. Figure 3: Space-time picture of string hadronization. $$A$$ is the world-sheet of the string; $$h_1,\ldots,h_n$$ represent produced hadrons.

The simplest case to consider first is the production of a quark-antiquark pair in electron-positron annihilation (see Figure 3). When produced at the annihilation point, the quark and antiquark are moving rapidly apart, each with half the total energy in the centre-of-mass frame. As the gluonic string is stretched between them, its potential energy grows at the expense of their kinetic energy. When the potential energy becomes of the order of hadron masses, it becomes energetically favourable for the string to break at some point along its length, through the creation of a new quark-antiquark pair, the new antiquark being at the end of the string segment connected to the original quark and similarly the new quark being connected to the antiquark. The two string segments then begin to stretch and break again, and so on until all the energy has been converted into quark-antiquark pairs connected by short string segments, which can be identified with hadrons. Because the quarks and antiquarks are moving almost at the speed of light, the string breaking events have space-like separation and so the sequence of hadron formation is frame-dependent. One can equally well view it as starting at one end of the string or the other, or in the middle. The types and momentum distributions of the hadrons formed depend on the quark-antiquark pairs produced in string breaking. These can be regarded as arising from vacuum fluctuations, so there is a preference for production of light flavours with low transverse momentum relative to the string axis.

The hadronization of a system more complicated than a single quark-antiquark pair depends on its colour structure, which is specified in the large-$$N$$ limit discussed earlier. Then each parton in the system has a unique colour partner, connected to it by a string segment which stretches and breaks as described above in its rest frame. Thus, for example, a quark-antiquark-gluon system formed in electron-positron annihilation will have a string stretching from the quark to the gluon and from there to the antiquark: the gluon produces a kink on the string, which is sharper for a high-momentum gluon. The result in the overall centre-of-mass frame is that there is more hadron production in the direction of the gluon and in the quark-gluon and antiquark-gluon angular regions, and less in the region between the quark and the antiquark, as observed experimentally. This contrasts with the situation in quark-antiquark-photon production, where the dragging of the string by the quark and antiquark results in more hadron production between them and less in the direction of the photon. This phenomenon is known as the string effect.

The development of parton showers means that in general one has a system of many quarks and gluons at the hadronization scale $$Q_0\ ,$$ which can however be assigned a large-$$N$$ colour structure specifying how the strings are to be stretched, as illustrated in Figure 2.

### Cluster model

The cluster hadronization model is based on the so-called preconfinement property of QCD discovered by Amati and Veneziano. They showed that at evolution scales much less than the hard subprocess scale, $$q\ll Q\ ,$$ the partons in a shower are clustered in colourless groups with an invariant mass distribution that is independent of the nature and scale of the hard subprocess, depending only on $$q$$ and the fundamental QCD scale $$\Lambda\ .$$ It is then natural to identify these clusters at the hadronization scale $$Q_0$$ as proto-hadrons that decay into the observed final-state hadrons. Figure 5: Cluster mass distribution. The distribution of colour singlet masses tends to a universal curve, whereas that of random quark-antiquark combinations does not.

In practical implementations of the model, the colour structure is again specified in the large-$$N$$ limit. At a scale at or near the cutoff $$Q_0$$ the gluons in the shower are forced to split into quark-antiquark pairs, as illustrated in Figure 4, and these form clusters with the corresponding colour partners. For $$Q_0\sim 1$$ GeV, most clusters have masses below 3 GeV (see Figure 5) and can be decayed into hadrons using a simple isotropic quasi-two-body phase space model. For these, the limited cluster mass limits the amount of transverse momentum generated in hadronization. For clusters in the higher-mass tail of the distribution, a string-like model of sequential cluster decay is adopted. As may be seen by comparing Figure 2 and Figure 4, a difference from the string model is that here the string is always broken at a gluon, rather than just having a kink there.

## Underlying event

In hadron collider events that contain a hard subprocess, there is extra hadron production that cannot be ascribed to showering from the coloured partons participating in the subprocess. Furthermore this extra activity, known as the underlying event is greater than that in so-called minimum-bias events, i.e. collisions that do not yield an identifiable hard subprocess. The underlying event is believed to arise from collisions between those partons in the incoming hadrons that do not directly participate in the hard subprocess.

### Multiple parton interactions

The most common hard subprocess at a high-energy hadron collider, such as the LHC, is elastic gluon-gluon scattering, $$gg\to gg\ .$$ The leading-order differential cross section for this subprocess diverges at zero momentum transfer, due to the exchange of a massless virtual gluon. This divergence is presumably regularized below some momentum transfer $$t_{\rm min}$$ by higher-order and non-perturbative effects. Nevertheless, for reasonable values of $$t_{\rm min}$$ the integrated gluon-gluon scattering cross section is very large, larger even than the total proton-proton scattering cross section. This result is not absurd: it simply indicates that the average number of gluon-gluon collisions per proton-proton collisions is greater than one. Including the cross sections for elastic scattering of quarks, antiquarks and gluons in all possible combinations, all of which diverge in leading order, we see that multiple parton interactions, each with relatively small momentum transfer, are highly probable in hadron-hadron collisions. This is the basis on which modern event generators model both minimum-bias collisions and the underlying event.

To account for the extra hadron production when a hard subprocess is present, an event generator must model the impact parameter structure of hadron-hadron collisions. The partons in each incoming hadron are distributed over a transverse area of the order of 1 fm$$^2$$. The impact parameter of a collision is the transverse distance between the centroids of these areas before the collision. When the impact parameter is large, the areas overlap little and the collision is peripheral, with a low probability of a hard parton-parton interaction and few multiple interactions. On the other hand at small impact parameter the collision is central, with a large overlap, several multiple interactions and a higher probability of a hard interaction. Thus the presence of a hard subprocess is correlated with more multiple interactions and a higher level of underlying event activity.

The final stage of event generation is the sequential decay of any unstable hadrons produced in the hadronization process. The experimental data indicate that a large fraction of observed final-state particles come from the decays of excited hadronic states, so most of the states listed in the Review of Particle Physics need to be included, together with their decay modes. To deal with the hadronization of heavy-flavour (charm and bottom) quarks, and to avoid flavour biases in the final states generated, complete flavour multiplets must be included for each spin and parity. In some cases, this involves modelling the properties of states that are not yet well established experimentally. In addition, for event generation the sum of branching fractions for all decays of a given state must be unity, which may require the modelling of unlisted decay modes, particularly for heavy-flavour hadrons. The momentum distributions and correlations in multiparticle decays need to be included or modelled, and spin correlations should also be taken into account. All this makes the implementation of decays a laborious but crucial part of event generator development.

## General-purpose event generators

There are three main programs in use at present for the generation of simulated collider events. They incorporate different combinations of the approaches and models described above, as outlined very briefly below. For fuller details and references, see Buckley et al. (2011). For regularly updated comparisons with data, see MCPLOTS.

### HERWIG

The HERWIG event generator was originally developed in Fortran but that version, while still in use and maintained at a low level, has been superseded by the C++ version, Herwig++. A range of hard subprocesses are hard-coded in both, and Herwig++ includes a generator for new two-to-two subprocesses given the relevant Feynman rules. Both versions are based on angular-ordered parton showers, although with different definitions of the shower variables, and a cluster model of hadronization. A multiple-interaction model for the underlying event is built into Herwig++ but a separate module (JIMMY) has to be interfaced to the Fortran version.

### PYTHIA

PYTHIA is the most senior and established general-purpose event generator. The Fortran version PYTHIA 6 is still very widely used, but development concentrates on the C++ version, PYTHIA 8. They both have a wide range of hard-coded subprocesses. PYTHIA 6 has the options of parton showers with virtuality or transverse momentum as the evolution variable, whereas PYTHIA 8 is based on dipole showering. Both versions use the Lund string model for hadronization and a highly developed multiple-interaction model for the underlying event.

### SHERPA

SHERPA is the most recent contender amongst general-purpose event generators and has been coded in C++ from the start. It has built-in subprocess generators for the Standard Model (both two-to-two and higher multiplicity) and for new models given the Feynman rules. A dipole formulation is used for parton showering, and a cluster model for hadronization. For the underlying event SHERPA uses a multiple-interaction model based on that of PYTHIA but differing in some respects.

All three generators have their own hadron decay modules with extensive tables of particle properties, branching fractions and decay distributions. While these have much in common, they are not identical since, as explained above, significant physics choices and modelling are involved. For the same reason they are not interchangeable, since these features influence the values of parameters relevant to other parts of the program that are needed to give the best agreement with experimental data.

## Acknowledgements

The work of B.W. was partly performed at the Kavli Institute for Theoretical Physics, University of California, Santa Barbara, supported in part by the National Science Foundation under Grant No. NSF PHY05-51164 and in part by a Leverhulme Trust Emeritus Fellowship.