HERWIG 6.5
An event generator for Hadron Emission Reactions With Interfering Gluons (including supersymmetric processes)1
Abstract: HERWIG is a general-purpose Monte Carlo event generator, which includes the simulation of hard lepton-lepton, lepton-hadron and hadron-hadron scattering and soft hadron-hadron collisions in one package. It uses the parton-shower approach for initial- and final-state QCD radiation, including colour coherence effects and azimuthal correlations both within and between jets. This article updates the description of HERWIG published in 1992, emphasising the new features incorporated since then. These include, in particular, the matching of first-order matrix elements with parton showers, a more correct treatment of spin correlations and heavy quark decays, and a wide range of new processes, including many predicted by the Minimal Supersymmetric Standard Model, with the option of R-parity violation. At the same time we offer a brief review of the physics underlying HERWIG, together with details of the input and control parameters and the output data, to provide a self-contained guide for prospective users of the program. This version of the manual (version 3) is updated to HERWIG version 6.5, which is expected to be the last major release of Fortran HERWIG. Future developments will be implemented in a new C++ event generator, HERWIG++.

Table of Contents

1  Introduction

HERWIG is a general-purpose event generator for high-energy processes, with particular emphasis on the detailed simulation of QCD parton showers. The program provides a full simulation of hard lepton-lepton, lepton-hadron and hadron-hadron scattering and soft hadron-hadron collisions in a single package, and has the following special features: Several of these features were already present in HERWIG version 5.1 and were described accordingly in some detail in ref. [1]. This was in turn based on the work of refs. [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. In the present article we concentrate therefore on the new features incorporated into HERWIG since 1992. The most important of these are the matching of first-order matrix elements with parton showers [15, 16, 17, 18, 19, 20], a more correct treatment of spin correlations [21] and heavy quark decays, again including matrix element matching [19], and a wide range of new hard processes. In particular HERWIG now includes the production and decay of supersymmetric particles, with or without the assumption of R-parity conservation [22].

The precise HERWIG version described here is 6.500, which we shall generally refer to as ``version 6'' in the following.

Let us recall briefly the main features of a generic hard, i.e. high momentum transfer, process of the type simulated by HERWIG. It can be divided notionally into four components, corresponding roughly to increasing scales of distance and time:
  1. Elementary hard subprocess. A pair of incoming beam particles or their constituents interact to produce one or more primary outgoing fundamental objects. This can be computed exactly to lowest order in perturbation theory. The hard momentum transfer scale Q together with the colour flow of the subprocess set the boundary conditions for the initial- and final-state parton showers, if there are any.
  2. Initial- and final-state parton showers. A parton constituent of an incident beam hadron with low spacelike virtuality can radiate timelike partons. In the process it decreases its energy to a fraction x of that of the beam and increases its spacelike virtuality, which is bounded in absolute value by the scale Q of the hard subprocess. This initial-state emission process leads to the evolution of the structure function F(x,Q) of the incident hadron. On a similar time-scale, an outgoing parton with large timelike virtuality can generate a shower of partons with lower virtuality. The amount of emission depends on the upper limit on the virtuality of the initiating parton, which is again controlled by the momentum transfer scale Q of the hard subprocess. Timelike partons from the initial-state emission may in turn initiate parton showering. The coherence of soft gluon emission from different parton showers is controlled by the colour flow of the subprocess. Within the showers, it is simulated by angular ordering of successive emissions.
  3. Heavy object decays. Massive produced objects such as top quarks, electroweak gauge and Higgs bosons, and possibly non-Standard Model (e.g. SUSY) particles, can decay on time-scales that may be shorter than or comparable to that of the QCD parton showers. Depending on their nature and the decay mode, they may also initiate parton showers before and/or after decaying.
  4. Hadronization process. In order to construct a realistic simulation one needs to combine the partons into hadrons. This applies to the constituent partons of incoming hadronic beams as well as to the outgoing products of parton showering, which give rise to hadronic jets. This hadronization process takes place at a low momentum transfer scale, for which the strong coupling aS is large and perturbation theory is not applicable. In the absence of a firm theoretical understanding of non-perturbative processes, it must be described by a phenomenological model. The model adopted in HERWIG is intended to disrupt as little as possible the event structure established in the parton showering phase. Showering is terminated at a low scale, Q0<1 GeV, and the preconfinement property of perturbative QCD [23, 24] is exploited to form colour-neutral clusters [4] which decay into the observed hadrons. Initial-state partons are incorporated into the incoming hadron beams through a soft, non-perturbative ``forced branching'' phase of spacelike showering. The remnants of incoming hadron, i.e. those constituent partons which do not participate in the hard subprocess, undergo a soft ``underlying event'' interaction modelled on soft minimum bias hadron-hadron collisions.
After a brief, practical section on the use of HERWIG, in the following sections we discuss each of these components in turn, concentrating on the new features since version 5.1.

2  Using HERWIG

The latest version of the program, together with all relevant information, is always available via the official HERWIG web page:
http://hepwww.rl.ac.uk/theory/seymour/herwig/
which is also mirrored at CERN:
http://home.cern.ch/seymour/herwig/

The program is written in Fortran and the user has to modify the main program HWIGPR to generate the type and number of events required. See section 8.1 for a sample main program. The program operates by setting up parameters in common blocks and then calling a sequence of subroutines to generate an event. Parameters not set explicitly in the block data HWUDAT or in HWIGPR are set to default values in the main initialisation routine HWIGIN. Output data are delivered in the LEP standard common block HEPEVT  [25, 26]. Note that all real variables accessible to the user, including those in HEPEVT, are of type DOUBLE PRECISION.

Since version 6.3, to take account of the increased energy and complexity of interactions at LHC and future colliders, the default value of the parameter NMXHEP, which sets the array sizes in the standard /HEPEVT/ common block, has been increased to 4000.

To generate events the user must first set up the beam particle names PART1, PART2 (type CHARACTER*8) and momenta PBEAM1, PBEAM2 (in GeV/c), a process code IPROC and the number of events required MAXEV.

See section 4 for beams and processes available.

All analysis of generated events (histogramming, etc.) should be performed by the user-provided routines HWABEG (to initialise the run), HWANAL (to analyse an event) and HWAEND (to terminate the run).

A detailed event summary is printed out for the first MAXPR events (default MAXPR=1). Setting IPRINT=2 lists the particle identity codes, properties and decay schemes used in the program.

The programming language is standard Fortran 77 as far as possible. However, the following may require modification for running on some computers In the INCLUDE file, common blocks have been regrouped so that all commons that were present in version 6.1 are unchanged. All new variables since version 6.1 are in separate common blocks, which are frozen after each version.

It is anticipated that version 6.5, with minor modifications in the series 6.5xx, will be the last Fortran version of HERWIG. Substantial physics improvements will be included in the new C++ event generator HERWIG++ [27], to be released in 2003. Watch for progress on the HERWIG++ web page:
http://www.hep.phy.cam.ac.uk/~gieseke/Herwig++/

3  Physics simulated by HERWIG

3.1  Elementary hard subprocess

In HERWIG version 6 there is a fairly large library of QCD, electroweak and supersymmetric elementary subprocesses. These are listed and discussed in section 4.

The hard subprocess plays an important rôle in defining the phase space of any associated initial- and final-state parton showers. As discussed in ref. [1] and references therein, the parton showers are branching processes in which the branchings are ordered in angle from a maximum to a minimum value determined by the cutoff Q0. The maximum value is determined by the elementary subprocess and is due to interference among soft gluons. The general result [7, 6, 28] is that the initial and final branchings are approximately confined within cones around the incoming and outgoing partons from the elementary subprocess. For the branching of parton i, the aperture of the cone is defined by the direction of the other parton j that is colour-connected to i. The relation between soft gluon interference and the colour connection structure of the elementary subprocess leads to many detectable effects in hadronic final states. For recent examples, see e.g. refs. [29, 30].

For a general process there are various contributions with different colour connections. The HERWIG library of elementary subprocesses includes the separate colour connection contributions. In general there is some ambiguity in the separation of contributions which are suppressed by inverse powers of the number of colours Nc. In earlier versions of HERWIG, these sub-leading terms were divided amongst the various colour connection contributions according to the recipe in ref. [7]. In version 6 the following improved prescription [31] has been followed, for both the QCD and SUSY QCD subprocesses. The matrix element squared | Mi|2 for each colour connection is computed in the limit Nc®¥ and the corresponding contribution is defined as
| Mi|2
 
å
j
 | Mj|2
 | M|2 ,
where | M|2 represents the sum over all colour connections including terms sub-leading in Nc. This ensures that each colour connection contribution is positive-definite and has the correct pole structure, while the sum of contributions yields the exact (Born-level) cross section.

Parton emission into phase space regions which are outside the above-mentioned angular cones, called in the following ``dead zones'', do not contribute to leading order and often not even to next-to-leading order. However, for a more complete description of the event it is also necessary to take into account emission into these ``dead zones'', which may be done using exact fixed-order matrix elements (see the following sections).

Another important function of the elementary subprocess is to set up the polarizations of any electroweak bosons, heavy (e.g. SUSY) particles or gluon jets that may be involved. These polarizations give rise to angular asymmetries and correlations in boson decays and jet fragmentation. They are included in HERWIG for many of the subprocesses provided, using the approach of refs. [9, 10, 21] to generate all correlations in decays, and in jet fragmentation to leading-logarithmic accuracy.

3.2  Parton showers

3.2.1  Final-state showers

Final-state parton showering in HERWIG is generated by a coherent branching algorithm with the following properties:
  1. The energy fractions are distributed according to the leading-order Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) splitting functions.
  2. The full available phase space is restricted to an angular-ordered region. Such a restriction is the result of interference and correctly takes important infrared singularities into account. At each branching, the angle between the two emitted partons is smaller than that of the previous branching.
  3. The emission angles are distributed according to the Sudakov form factors, which sum the virtual corrections and unresolved real emissions. The Sudakov form factor normalizes the branching distributions to give the probabilistic interpretation needed for a Monte Carlo simulation. This fact is a consequence of unitarity and of the infrared finiteness of inclusive quantities.
  4. The azimuthal angular distribution in each branching is determined by two effects:
    1. for a soft emitted gluon the azimuth is distributed according to the eikonal dipole distribution [7];
    2. for non-soft emission one finds azimuthal correlations due to spin effects. See [9, 10] for the method used to implement these correlations in full, to leading collinear logarithmic accuracy, in HERWIG.
  5. In each branching the scale of aS is the relative transverse momentum of the two daughter partons.

  6. In the case of heavy flavour production the mass of the quark modifies the angular-ordered phase space. The most important effect is that the soft radiation in the direction of the heavy quark is depleted. One finds that emission within an angle of order M/E is suppressed, M and E being the mass and energy of the heavy quark: this is known as the ``dead cone'' [32].
Specifically, the HERWIG parton shower evolution is done in terms of the parton energy fraction z and an angular variable x. In the parton splitting i® jk, zj = Ej/Ei and xjk = (pj· pk)/(Ej Ek). Thus xjk~1/2qjk2 for massless partons at small angles. The values of z are chosen according to the DGLAP splitting functions and the distribution of x values is determined by the Sudakov form factors. Angular ordering implies that each x value must be smaller than the x value for the previous branching of the parent parton.

The parton showers are terminated as follows. For partons there is a cutoff of the form Qi = mi + Q0, where mi (i=1,...,6 for d, u, s, c, b, t) is set by the relevant mass parameter RMASS(i) and Q0 is set by the quark and gluon virtuality cutoff parameters VQCUT and VGCUT (see section 5). Showering from any parton stops when a value of x below Qi2/Ei2 is selected for the next branching. For heavy quarks, the condition x > Qi2/Ei2 corresponds to the ``dead cone" mentioned above. At this point the parton is put on mass-shell or given a small non-zero effective mass in the case of gluons.2 Working backwards from these on-shell partons, one can now construct the virtual masses of all the internal lines of the shower and the overall jet mass, from the energies and opening angles of the branchings. Finally one can assign the azimuthal angles of the branchings, including EPR-type correlations (from Einstein-Podolski-Rosen [33]), and deduce completely all the 4-momenta in the shower.

Next the parton showers are used to replace the (on mass-shell) partons that were generated in the original hard process. This is done in such a way that the jet 3-momenta have the same directions as the original partons in the c.m. frame of the hard process, but they are boosted to conserve 4-momentum taking into account their extra masses.

The main improvements in the final-state emission algorithm of HERWIG version 6, relative to version 5.1, are as follows.

The Sudakov form factors can be calculated using the one-loop or two-loop aS, according to the variable SUDORD (default = 1). The parton showering still incorporates the two-loop aS in either case but if SUDORD = 1 this is done using a veto algorithm, whereas if SUDORD = 2 no vetoes are used in the final-state evolution. The usefulness of this option is discussed briefly in section 8.2.

Matrix element corrections have been introduced into final-state parton showers in e+e- and deep inelastic processes [15, 16], in heavy flavour decays [19] and in Drell-Yan processes [20] (see section 3.2.3).

In HERWIG, the angular-ordering constraint, which is derived for soft gluon emission, is applied to all parton shower vertices, including g® qq_. In versions before 6.1, this resulted in a severe suppression (an absence in fact) of configurations in which the gluon energy is very unevenly shared between the quarks. For light quarks this is irrelevant, because in this region one is dominated by gluon emissions, which are correctly treated. However, for heavy quarks this energy sharing (or equivalently the quarks' angular distribution in their rest frame) is a directly measurable quantity and was badly described. Related to this was an inconsistency in the calculation of the Sudakov form factor for g® qq_. This was calculated using the entire allowed kinematic range (with massless kinematics) for the energy fraction, 0£ z£ 1, while the z distribution generated was actually confined to the angular-ordered region, z,1-z ³ m/E x1/2.

In version 6, these defects are corrected as follows. We generate the E,x and z values for the shower as before. We then apply an a posteriori adjustment to the kinematics of the g® qq_ vertex during the kinematic reconstruction. At this stage, the masses of the q and q_ showers are known. We can therefore guarantee to stay within the kinematically allowed region. In fact, the adjustment we perform is purely of the angular distribution of the q and q_ showers in the g rest frame, preserving all the masses and the gluon four-momentum. Therefore we do not disturb the kinematics of the rest of the shower at all.

Although this cures the inconsistency above, it actually introduces a new one: the upper limit for subsequent emission is calculated from the generated E,x and x values, rather than from the finally-used kinematics. This correlation is of NNLL importance, so we can formally neglect it. It would be manifested as an incorrect correlation between the masses and directions of the produced q and q_ jets. This is, in principle, physically measurable, but it seems less important than getting the angular distribution itself right. In fact the solution we propose maps the old angular distribution smoothly onto the new, so the sign of the correlation will still be preserved, even if the magnitude is wrong. Even with this modification, the HERWIG kinematic reconstruction can only cope with particles that are emitted into the forward hemisphere in the showering frame. Thus one cannot populate the whole of kinematically-allowed phase space. Nevertheless, we find that this is usually a rather weak condition and that most of phase space is actually populated.

Using this procedure, we find that the predicted angular distribution for secondary b quarks at LEP energies is well-behaved, i.e. it looks reasonably similar to the leading-order result (1+cos2q*), and has relatively small hadronization corrections.

Real photon emission is included in timelike parton showering. The infrared photon cutoff is VPCUT, which defaults to 0.4 GeV. Agreement with LEP data is satisfactory if showering is used together with the matrix element correction to produce photons in the back-to-back region. The results are insensitive to VPCUT variations in the range 0.1-1.0 GeV. Setting VPCUT greater than the total c.m. energy switches off such emission. As an expedient way of improving the efficiency of photon final-state radiation studies, the electromagnetic coupling ALPHEM can be multiplied by a factor ALPFAC (default = 1) for all quark-photon vertices in jets, and in the `dead zone' in e+e-. Results at small photon-jet separation become sensitive to ALPFAC above about 5.

3.2.2  Initial-state showers

The theoretical analysis of initial-state showering is more complex than the final-state case. The most relevant parameters of the hard subprocess are the hard scale Q and the energy fraction x of the incoming parton after the emission of initial state radiation. For lepton-hadron processes x corresponds to the Bjorken variable, while for hadron-hadron processes x is related to Q2/s where s is the c.m. energy squared.

The main result is that for any value of x, even for x small [34], the initial-state emission process factorizes and can be described as a coherent branching process suitable for Monte Carlo simulations. The properties which characterize this process include all those discussed above for the final-state emission. However, in the initial-state case the angular-ordering restriction on the phase space applies to the angles qi between the directions of the incoming hadron and the emitted partons i.

For large x, the coherent branching algorithm sums correctly [13] not only the leading but also the next-to-leading contributions. This accuracy allows us to identify the relation between the QCD scale used in the Monte Carlo program and the fundamental parameter LMS. This is achieved by using the one-loop Altarelli-Parisi splitting functions and the two-loop expression for aS with the following universal relation between the scale parameter Lphys  [13] used in the simulation and LMS (here, Nf is the number of flavours)
Lphys = exp æ
ç
ç
è
67 - 3 p2 - 10 Nf/3
2 (33 - 2 Nf)
 ö
÷
÷
ø		 L
 
MS
    ~ 1.569 L
 
MS
      for Nf=5 .     (1)
Therefore a Monte Carlo simulation with next-to-leading accuracy can be used to determine LMS from semi-inclusive data at large momentum fractions.3

In the case of small values of x, the initial state branching process has additional properties, which are not yet included fully in HERWIG. This was discussed in ref. [1] and the situation remains unchanged since version 5.1.

The initial-state branching algorithm in HERWIG is of the backward evolution type. It proceeds from the elementary subprocess, at a hard scale set by colour coherence (see section 3.1), back to the hadron scale, set here by the spacelike cutoff parameter QSPAC. At this point there is a forced non-perturbative stage of branching which ensures that the emitting parton fits smoothly with the valence parton distributions of the incoming hadron.

Matrix-element corrections have been introduced into initial-state parton showers in deep inelastic [16] and Drell-Yan processes [20], as discussed in the following subsection.

To avoid double-counting of hard parton emission, all radiation at transverse momenta greater than the hard process scale EMSCA is vetoed. In the case of initial-state radiation, this affects all events, while for final-state radiation it only affects those events in which the two jets have a rapidity difference of more than about 3.4.

In the backward evolution of initial-state radiation for photons the ``anomalous'' branching qq_¬g is included. Variables ANOMSC(1,IBEAM) and ANOMSC(2,IBEAM) record the evolution scale and transverse momentum, respectively, at which an anomalous splitting was generated in the backward evolution of beam IBEAM. If zero, then no such splitting was generated.

The treatment of forced branching of gluons and sea (anti-)quarks in backward evolution has been improved, by allowing it to occur at a random scale between the spacelike cutoff QSPAC and the infrared cutoff, instead of exactly at QSPAC as before. A new option ISPAC = 2 allows the freezing of structure functions at the scale QSPAC, while evolution continues to the infrared cutoff. The default, ISPAC = 0 is equivalent to previous versions, in which perturbative evolution stops at QSPAC.

The width of the (gaussian) intrinsic transverse momentum distribution of valence partons in incoming hadrons at scale QSPAC is set by the parameter PTRMS (default value zero). The intrinsic transverse momentum is chosen before the initial state cascade is performed and is held fixed even if the generated cascade is rejected. This is done to avoid correlation between the amount of perturbative and non-perturbative transverse momentum generated.

It is possible to completely switch off initial-state emission, by setting NOSPAC = .TRUE., in which case only the forced splitting of non-valence partons is generated.

3.2.3  Matrix-element corrections

One of the new features of HERWIG 6 is the matching of first-order matrix elements with parton showers.

The HERWIG parton showers are performed in the soft or collinear approximation and emission is allowed only in regions of the phase space satisfying the condition x<1 or x<z2, for the final- (timelike branching) and initial-state (spacelike branching) radiation respectively, where x and z are the showering variables defined above.

The emission is entirely suppressed inside the so-called dead zones (x >1 or x>z2), corresponding to hard and/or large-angle parton radiation. According to the exact matrix elements, the radiation in the dead zones is suppressed, since it is not soft or collinear logarithmically enhanced, but it is not completely absent as happens in the HERWIG standard shower algorithm. The HERWIG parton cascades need to be supplemented by matrix-element corrections for a full description of the physical phase space.

The method of matrix-element corrections to the HERWIG parton showers is discussed in [16, 17]. The radiation in the dead zones is generated according to the exact first-order matrix element (`hard correction'); the shower in the already-populated region of the phase space is corrected by the use of the exact O(aS) amplitude any time an emission is capable of being the `hardest so far' (`soft correction'4).

By `hardest-so-far', we mean the radiation of a parton whose transverse momentum relative to the splitting one is larger than all those previously emitted. This is not always the first emission, as angular ordering does not necessarily imply ordering in transverse momentum. As shown in [16], if we corrected only the first emission, we would have problems in the implementation of the Sudakov form factor whenever a subsequent harder emission occurs, as we would find that the probability of hard radiation would depend on the infrared cutoff, which is clearly unphysical. Using the O (aS) result for the hardest-so-far emission in the already-filled phase space as well as in the dead zone allows one to have matching over the boundary of the dead zone itself.

Since the fraction of events which receive a hard correction is typically small, we neglect multiple hard emissions in the dead zones and rely on the first-order result plus showering in those regions.

Our method is quite different from the one used to implement matrix-element corrections in JETSET [35], where the parton shower probability is applied over the whole phase space and the first-order amplitude is used only to correct the first emission.

Following these general prescriptions, matrix-element corrections have been implemented in some e+e- processes [15] (including e+e-® WW/ZZ), deep inelastic lepton scattering [18], top quark decay [19], and Drell-Yan processes [20]. Those for the process gg® Higgs are now in progress [36].

The variables HARDME (default = .TRUE.) and SOFTME (default = .TRUE.) allow respectively the application of hard and soft matrix-element corrections to the HERWIG parton cascades.

3.3  Heavy flavour production and decay

Heavy quark decays are treated as secondary hard subprocesses. Top quarks and any hypothetical heavier quarks always decay before hadronization. Heavy-flavoured hadrons are split into collinear heavy quark and spectator and the former decays independently. After decay, parton showers may be generated from coloured decay products, in the usual way. See ref. [11] for details of the treatment of colour coherence in these showers.

In HERWIG version 6 matrix-element corrections to the simulation of top quark decays are available. The routine HWBTOP implements the hard corrections; HWBRAN has been modified to implement the soft corrections. Since the dead zone includes part of the soft singularity, a cutoff is required: only gluons with energy above GCUTME (default value 2 GeV) in the top rest frame are corrected. Physical quantities are not strongly dependent on GCUTME in the range 1 to 5 GeV, after the typical experimental cuts are applied. For more details see ref. [19].

The structure of the program has been altered so that the secondary hard subprocess and subsequent fragmentation associated with each partonic heavy hadron decay appear separately in the event record. Thus top quark decays are treated individually as are any subsequent bottom hadron partonic decays. Note that the statement CALL HWDHOB, which deals with the decays of all heavy objects (including SUSY particles), must appear in the main program between the calls to HWBGEN and HWCFOR, in order to carry out any decays before hadronization.

The partonic decay fractions of heavy quarks are specified in the decay tables like the decay modes of other particles. This permits different decays to be given to individual heavy hadrons. Changes to the decay table entries can be made on an event by event basis if desired. Partonic decays of charm hadrons and quarkonium states are also now supported. The order of the products in a partonic decay mode is significant. For example, if the decay is Q ® W+q ® (f+f_')+q occurring inside a Qs_ hadron, the required orderings are:
Q+
_
s
 
 ®
(f+
_
f
 
 ')+(q+
_
s
 
 )
  or
(q+
_
f
 
 ')+(f+
_
s
 
 )      (`colour rearranged') .
In both cases the (V-A)2 matrix element-squared is proportional to (pQ· p2)(p1· p3), where p1 etc. correspond to the ordering given. Decays of heavy-flavoured hadrons to exclusive non-partonic final states are also supported. No check is made against double counting from partonic modes. However this is not expected to be a major problem for the semi-leptonic and two-body hadronic modes supplied.

The default masses of the c and b quarks have been lowered to 1.55 and 4.95 respectively: this corresponds to the mass of the lightest meson minus the u or d quark mass. This increases the number of heavy mesons, and hence total multiplicities, and slightly softens their momentum spectrum. The rate of photoproduced charm states increases and B-p momentum correlations become smoother. The default top quark mass is 174.3 GeV/c2. The same value is used in the production and decay matrix elements and for all kinematics. Note that higher-order corrections are not fully included, and so the HERWIG top mass does not necessarily correspond to that defined in any particular rigorous scheme (e.g. the pole mass or the MS running mass). However, since it is probably the decay kinematics that are most sensitive to this parameter, it should be close to the pole mass. See subsection 4.2.1 for notes on the treatment of quark masses in various processes.

3.4  Gauge and Higgs boson decays

MODBOS(i) W± Decay Z0 Decay
0 all all
1 qq_ qq_
2 en e+e-
3 µn µ+µ-
4 tn t+t-
5 en + µn e+e- + µ+µ-
6 all nn_
7 all bb_
>7 all all

Table 1: W,Z decays.

The total decay widths of the electroweak gauge bosons V=W,Z are specified by the input parameters GAMW and GAMZ. Their branching fractions to various final states are computed automatically from the other SM input parameters. Which decays actually occur is controlled as follows. The variable MODBOS(i) controls the decay of the ith gauge boson per event (table 1).

All entries of MODBOS default to 0. Bosons which are produced in pairs (i.e. from VV pair production, or Higgs decay) are symmetrized in MODBOS(i) and MODBOS(i+1). For processes which directly produce gauge bosons, the event weight includes the branching fraction to the requested decay, but this is only true for Higgs production if decay to W+W-/Z0Z0 is forced (IPROC = 310, 311 but not 399, etc.). Users can thus force Z ® bb_ decays, with MODBOS(i)=7. For example, IPROC = 250, MODBOS(1) = 7, MODBOS(2) = 0 gives Z0Z0 production with one Z0 decaying to bb_.

The spin correlations in the decays are handled in one of two ways:
  1. The diagonal members of the spin density matrix are stored in RHOHEP(i,IHEP), where i=1,2,3 for helicity=i-2 in the centre-of-mass frame of their production, for processes where this matrix is diagonal (i.e. there is no interference between spin states).
  2. The correlations in the decay are handled directly by the production routine where (1) is not possible.
The processing of the parton showers in hadronic W and Z decays is handled in the rest frame of the vector boson if WZRFR is .TRUE. (the default), otherwise in the lab frame. In the latter case, which was the default in earlier versions, the initial cone angles of the showers depend on the velocity of the boson, which leads to a slight Lorentz non-invariance of decay distributions.

The total decay width of the SM Higgs boson is computed from its input mass RMASS(201) and stored as GAMH. Its decay branching fractions are also computed and stored in BRHIG(I): I = 1-6 for dd_,..., tt_; I = 7-9 for e+e-,..., t+t-; I = 10,11,12 for W+W-, Z0Z0, gg. Non-SM Higgs bosons, on the other hand, such as those in supersymmetric models, have to have their widths and decay tables provided as input data (see section 3.5.1). To avoid any ambiguity, the SM Higgs boson has a distinct identity code in HERWIG and is represented by the special symbol HSM0.

There are two choices for the treatment of the SM Higgs width, both controlled by the variable IOPHIG:
IOPHIG = 2I + J ,
where I and J are both zero or one. Whenever a Higgs boson is generated, its mass is chosen from a distribution that, for heavy SM Higgs bosons, can be rather broad. The choice of I makes a significant difference to the physical meaning of the distribution generated: for I=0, the cross section corresponds to the tree level process containing an s-channel Breit-Wigner resonance for the Higgs boson with a running Higgs width. As discussed in [37], this neglects important contributions from interference with non-resonant diagrams and can violate unitarity at high energy, so I=1 (the default) uses the improved prescription of [37]. This replaces the s-channel propagator by an effective propagator that sums the interference terms to all orders. This increases the cross section below resonance and decreases it above, causing an overall increase in cross section. More details can be found in [37].

The variable J is a more technical parameter that does not affect the physical results, only the method by which they are generated: J=1 (the default) generates the mass according to a fixed-width Breit-Wigner resonance, while J=0 biases the distribution more towards higher masses. In either case, the appropriate jacobian factor is included in the event weight, so that the physical cross section is independent of J.

In all the above cases, the SM Higgs mass distribution is restricted to the range [mH-GAMMAX×GH , mH+GAMMAX×GH]. GAMMAX defaults to 10, but in the non-perturbative region mH³ 1 TeV should be reduced to protect against poor weight distributions. These considerations do not affect the distribution noticeably for mH£500 GeV, and GAMMAX can safely be increased if necessary.

For a SUSY Higgs, the width is never large enough for unitarity to be violated and these issues are unimportant. In this case, the mass distribution is chosen according to a fixed-width Breit-Wigner resonance, like that of any other SUSY particle.

The SM Higgs decays that can occur are normally controlled by the process code IPROC, as in IPROC=300+ID for example: ID= 1-6 for quarks, 7-9 for leptons, 10/11 for W+W-/Z0Z0 pairs, and 12 for photons. In addition ID= 0 gives quarks of all flavours, and ID= 99 gives all decays. For each process, the average event weight is the cross section in nb times the branching fraction to the requested decay. The branching ratios to quarks use the next-to-leading logarithm corrections, those to W+W-/Z0Z0 pairs allow for one or both bosons being off mass-shell.

All Higgs vertices include an optional enhancement factor to account for non-SM and non-MSSM couplings. The amplitudes for all Higgs vertices are multiplied by the factor ENHANC(ID) where ID is the same as in IPROC=300+ID except the gg H `vertex' which is calculated from ENHANC(6) and ENHANC(10) for the top and W± loops. This allows the simulation of the production of any chargeless scalar Higgs-like particle. Note however that pseudoscalar and charged Higgs bosons, and processes involving more than one Higgs particle (e.g. the decay H0® h0Z) are not included this way (see section 4.7).

The array ENHANC(ID) is initialised as usual in HWIGIN. Note, however, that it will be overwritten if MSSM Higgs production is required by IPROC. In that case, as mentioned earlier, the Higgs widths and decay modes are simply read from an input particle data file (see section 3.5.1).

3.5  Supersymmetry

Particle   Spin Particle   Spin
quark q 1/2 squarks qL,R~ 0
charged lepton l 1/2 charged sleptons lL,R~ 0
neutrino n 1/2 sneutrino n~ 0
gluon g 1 gluino g~ 1/2
photon g 1 photino g~ 1/2
neutral gauge boson Z0 1 zino Z~ 1/2
neutral Higgs bosons h0,H0,A0 0 neutral Higgsinos H~1,20 1/2
charged gauge boson W± 1 wino W~± 1/2
charged Higgs boson H± 0 charged Higgsino H~± 1/2
graviton G 2 gravitino G~ 3/2
W~±, H~± mix to form 2 chargino mass eigenstates c~1±, c~2±
g~, Z~, H~1,20 mix to form 4 neutralino mass eigenstates c~10,c~20,c~30,c~40
t~L,t~R (and similarly b~, t~) mix to form the mass eigenstates t~1, t~2

Table 2: SUSY particle content.

HERWIG now includes the production and decay of superparticles, as given by the Minimal Supersymmetric Standard Model (MSSM) [22]. The mass spectrum and decay modes, being read from input files (see below), are completely general. The particle content, listed in the following table (table 2), includes the gravitino/goldstino. For sparticles that mix, the subscripts label the mass eigenstates in the ascending order of mass. The two Higgs Doublet Model (2HDM) Higgs sector, intrinsic to the MSSM, is also included. The three neutral Higgs bosons are denoted by h0, H0 and A0.

The SUSY particle names in HERWIG are as shown in table 3. IDHW is the HERWIG identity code and IDPDG is the corresponding Particle Data Group code [38].
IDHW   NAME IDPDG IDHW   NAME IDPDG
401		 dL~ 'SSDL ' 1000001 413 dR~ 'SSDR ' 2000001
402 uL~ 'SSUL ' 1000002 414 uR~ 'SSUR ' 2000002
403 sL~ 'SSSL ' 1000003 415 sR~ 'SSSR ' 2000003
404 cL~ 'SSCL ' 1000004 416 cR~ 'SSCR ' 2000004
405 b1~ 'SSB1 ' 1000005 417 b2~ 'SSB2 ' 2000005
406 t1~ 'SST1 ' 1000006 418 t2~ 'SST2 ' 2000006
               
425 eL~ 'SSEL- ' 1000011 437 eR~ 'SSER- ' 2000011
426 ne~ 'SSNUEL ' 1000012        
427 µL~ 'SSMUL- ' 1000013 439 µR~ 'SSMUR- ' 2000013
428 nµ~ 'SSNUMUL ' 1000014        
429 t1~ 'SSTAU1- ' 1000015 441 t2~ 'SSTAU2- ' 2000015
430 nt~ 'SSNUTL ' 1000016        
               
449 g~ 'GLUINO ' 1000021 458 G~ 'GRAVTINO' 1000039
450 c~10 'NTLINO1 ' 1000022 451 c~20 'NTLINO2 ' 1000023
452 c~30 'NTLINO3 ' 1000025 453 c~40 'NTLINO4 ' 1000035
454 c~1+ 'CHGINO1+' 1000024 455 c~2+ 'CHGINO2+' 1000037
203 h0 'HIGGSL0 ' 26 204 H0 'HIGGSH0 ' 35
205 A0 'HIGGSA0 ' 36 206 H+ 'HIGGS+ ' 37

Table 3: SUSY particle names.

Antiparticles generally appear in sequence after the corresponding particles, e.g. antisquarks dL~*-t1~* at IDHW = 407-412, dR~*-t2~* at 419-424. They have 'BR' added to the name, e.g. 'SSDLBR ', or opposite charge, and negative PDG codes. A full list can be obtained using the print option IPRINT = 2 (see section 6).

Note that the HERWIG particle labelling of the lightest MSSM Higgs boson departs from the PDG recommendation: it is given PDG code 26, to avoid confusion with the SM Higgs boson (PDG code 25) in our implementation (specifically, in our use of the array ENHANC for the MSSM processes: see the relevant Higgs sections for more details).

HERWIG does not contain any built-in models for SUSY scenarios beyond the MSSM, such as, Supergravity (SUGRA) or Gauge Mediated Symmetry Breaking (GMSB). In all cases the SUSY particle spectrum and decay tables must be provided just like those for any other particles. The subroutine HWISSP, if called, reads these from an input file. The production subprocesses are then generated by HWHESP, in lepton-antilepton collisions, HWHSSP, in hadron-hadron collisions, or by one of the ¬  Rp production routines. The decays of the sparticles produced, as well as any top quarks or Higgs bosons, are then performed by HWDHOB.

3.5.1  Data input

A package ISAWIG, see section 9.4, has been created to work with ISAJET [39] to produce a file containing the SUSY particle masses, lifetimes and decay modes. This package takes the outputs of the ISAJET SUGRA or general MSSM programs and produces a data file in a format that can be read into HERWIG by the subroutine HWISSP. In principle the user can produce a similar file provided that the correct format is used, as explained below.

For the mixing terms of the MSSM lagrangian we follow the Haber-Kane [40, 41] conventions, so that we differ from ISAJET on the sign for gaugino masses, the ordering and signs of the gaugino current eigenstates, the interchange of the rows and columns of the gaugino mixing matrices, and the sign of the neutral Higgs mixing angle a.

In addition to the decay modes included in the ISAJET package ISAWIG allows for the possibility of violating R-parity and includes the calculation of all 2-body squark and slepton, and 3-body gaugino and gluino R-parity violating (¬  Rp) decay modes.

It can happen that some of the SUSY particle decay modes generated by ISAJET are found to be kinematically forbidden in HERWIG, owing to the slightly different values assumed for the light quark masses. In this case a warning message is printed by HERWIG and these modes are deleted, the other branching ratios being rescaled accordingly. Such modes normally have negligible ISAJET branching ratios anyway, because of their tiny phase space.

The input file organisation expected by HWISSP is as follows. First the SUSY particle and top quark masses and lifetimes (in seconds) are given according to their HERWIG identity codes IDHW, for example: 
         65
         401 927.3980    0.74510E-25
         402 925.3307    0.74009E-25
         ....etc.
That is,
NSUSY = Number of SUSY + top particles
IDHW, RMASS(IDHW), RLTIM(IDHW)
repeated NSUSY times.
Next each particle's decay modes together with their branching ratios and matrix element codes are given as, for example: 
         6
         401  0.18842796E-01     0   450     1     0     0     0
          :         :            :    :      :     :     :     :
         401  0.32755006E-02     0   457     2     0     0     0
         6
         402  0.94147678E-02     0   450     2     0     0     0
         ....etc.
That is,
Number of decay modes for a given particle IDK
IDK(IM), BRFRAC(IM), NME(IM), IDKPRD(1-5,IM)
repeated for each mode IM
all repeated for each particle (NSUSY times).
The order in which the decay products appear is important in order to obtain appropriate showering and hadronization. The correct orderings are indicated in the table below (table 4).
Decaying No. of Type of mode Order of decay products
Particle products   1st 2nd 3rd
top 2 two-body to Higgs Higgs bottom  
  3 three-body via Higgs/W quarks or leptons bottom
      from W/Higgs  
gluinos 2 without gluon any order  
    with gluon gluon colour  
        neutral  
  3 R-parity conserved colour quark or antiquark
      neutral  
squark 2 gaugino/gluino gaugino quark  
or slepton   quark/lepton gluino lepton  
  3 weak sparticle particles from
        W decay
squarks 2 lepton number violated quark lepton  
    baryon number violated quark quark  
sleptons 2 lepton number violated quark or antiquark  
Higgs 2 (s)quark-anti(s)quark (s)quark or anti(s)quark  
    (s)lepton-anti(s)lepton (s)lepton or anti(s)lepton  
  3 all three-body colour quark or antiquark
      neutral lepton or antilepton
gaugino 2 squark-quark quark or squark  
    slepton-lepton lepton or slepton  
  3 R-parity conserved colour quark or antiquark
      neutral lepton or antilepton
gaugino 3 R-parity violating particles in the order i,j,k as
or gluino     in the superpotential

Table 4: SUSY decay product ordering.

New matrix element codes have been added for SUSY and Higgs decays: The indices i,j,k in ¬  Rp gaugino/gluino decays refer to the ordering of the indices in the ¬  Rp couplings in the superpotential. The convention is as in ref. [42].

Next a number of parameters derived from the SUSY lagrangian must be given. These are: the ratio of Higgs VEVs, tanb, and the scalar Higgs mixing angle, a; the mixing parameters for the Higgses, gauginos and the sleptons; the trilinear couplings; and the Higgsino mass parameter µ.

Finally the logical variable RPARTY must be set .FALSE. if R-parity is violated, and the ¬  Rp couplings must also be supplied; otherwise not.

Details of the FORMAT statements employed can be found by examining the subroutine HWISSP.

HWISSP reads the data from UNIT = LRSUSY (default LRSUSY = 66). If the data are stored in a fort.LRSUSY file on a UNIX system5 no further action is required, but if the data are to be read from a file named fname.dat then appropriate OPEN and CLOSE statements must be added by hand: 
OPEN(UNIT=LRSUSY,FORM='FORMATTED',STATUS='UNKNOWN',FILE='fname.dat')
CALL HWISSP CLOSE(UNIT=LRSUSY)
A number of sets of SUSY parameter files, produced using ISAJET, for the standard LHC SUGRA and GMSB points are available from the ISAWIG home page: http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/

3.5.2  Processes

The implementation of supersymmetric particle production processes in lepton-antilepton and hadron-hadron collisions is described in sections 4.4 and 4.7, respectively. As in SM processes, we do not include any higher-order QCD corrections to the relevant matrix elements, even though these are now known in many cases. This is in order to avoid double counting of corrections generated approximately by parton showering.

The procedure by which the MSSM matrix elements describing the elementary hard subprocesses are interfaced to the initial- and final-state parton showers is similar to that described in sections 3.1 and 3.2 for the SM case. However, at present showers are generated from partons but not from spartons, whose short lifetimes should make this a reasonable approximation.

We include matrix elements and spin correlations in all SUSY decays as described in [21]. Breit-Wigner mass distributions have been implemented for all unstable SUSY particles, according to their widths as given in the input file. For a list of the most relevant decay modes of MSSM particles, see ref. [22].

One difference between the SM and MSSM implementations of the Higgs decay channels should be mentioned. Whereas for the SM Higgs boson all decay rates are calculated by HERWIG itself, in terms of the other (known) parameters of the SM, in the case of the MSSM scalars these are passed to the generator through the data files.

A more detailed description of the new SUSY reactions in HERWIG, along with the relevant formulae for the hard scattering processes involved, can be found in [22].

We include the possibility of violating R-parity both in the decays of the sparticles and in the initial hard subprocess. The ¬  Rp model we consider is specified by a spectrum which can be given in either the SUGRA, GMSB or a more general MSSM scenario, and a set of ¬  Rp couplings at the weak scale. We include only the tri-linear couplings and not the bi-linear terms which mix the leptons and gauginos; a recent review of these models can be found in [43]. However we do include the possibility that more than one of the ¬  Rp couplings can be non-zero. All the two-body squark and slepton and three-body gaugino and gluino decay modes, and resonant production processes in hadron-hadron collisions are included, as well as a range of production processes in e+e- collisions. A decay matrix element is implemented for the three-body decays. The colour structure of these events is very different from that of the MSSM [42], due to the presence of the baryon-number violating (¬  B) vertex. This means that a new subroutine HWBRCN is required to handle the colour connections between jets in this case. A full discussion of the colour connections in these processes and the matrix elements for the cross sections and decays can be found in [42].

3.5.3  Related changes

A large number of changes have been made in HERWIG to enable SUSY processes to be included in hadron-hadron collisions. The main changes are:

3.6  Spin correlations

In addition to the spin correlations in parton showers and vector boson decays discussed above, in versions 6.4 and higher spin correlation effects are included in processes where top quarks, t leptons and SUSY particles are produced, as described in [21]. At the moment the effects are only calculated for the production of these particles in the following processes: IPROC=100-199, 700-799, 800-899, 1300-1399, 1400-1499, 1500-1599, 1700-1799, 2000-2099, 2800-2825, 3000-3030, 4000-4199. However, if these particles are produced in other processes, the spin correlation algorithm will still be used to perform their decays. The correlations are also included for the decay of the MSSM Higgs bosons, regardless of how they are produced.

The spin correlations are controlled by the logical variable6 SYSPIN [.TRUE.] which switches the correlations on. If required the correlations are initialised by the new routine HWISPN. This routine initialises the two, three and four body matrix elements.

The three and four body matrix elements can be used separately to generate the decay distributions without spin correlation effects. These are switched on by the switches THREEB [.TRUE.] for three body decay and FOURB [.FALSE.] for four body decays. The four body decays are only important in SUSY Higgs studies, and have small branching ratios. However, they take some time to initialise and are therefore switched off by default.

The initialisation of the spin correlations and/or decay matrix elements can be time consuming and we have therefore included an option to read/write the information. The information is written to unit LWDEC [88] and read from LRDEC [0]. If either are zero the data is not written/read. If IPRINT=2 then information on the branching ratios for the decay modes and the maximum weights for the matrix elements is outputted.

If the spin correlation (SYSPIN) or matrix element switches (THREEB, FOURB) are .TRUE., then the matrix element codes (NME entries) for the decays concerned are not used; the calculated matrix elements are used instead.

When we included spin correlations in HERWIG6.4 [21] we did not include either R-parity violating decays or decays producing gravitinos in the algorithm. This led to HERWIG stopping when such decays were included. This of course could be stopped by switching the spin correlations off, i.e. SYSPIN=.FALSE.. In version 6.5, we have included the relevant matrix elements for R-parity violating decays and hard processes and decays producing gravitinos. At the same time we have made changes so that at both the initialisation and event generation stages many of the terminal warnings which were caused by the code not having the correct matrix elements are now information-only warnings. If you still get terminal error messages from any of the spin correlation routines please let us know.

The effect of polarization for incoming leptonic beams in MSSM and ¬  Rp SUSY processes has also been included. These effects are included both in the production of SUSY particles and via the spin correlation algorithm in their decays.

3.7  Hadronization

For a general hard process in hadron-hadron collisions, we have to consider: (a) the representation of the incoming partons as constituents of the incident hadrons; (b) the conversion of the emitted partons into outgoing hadrons; (c) the `underlying soft event' associated with the presence of spectator partons.

The first of these is dealt with through the use of non-perturbative parton distribution functions, which are discussed below in section 4.1.1, and by the remnant hadronization model. The cluster model for hadron formation, remnant hadronization and the underlying event is as follows.

3.7.1  Cluster model

The preconfinement property mentioned in section 1 is used by HERWIG as the basis for a simple hadronization model which is local in colour and independent of the hard process and the energy [4, 7].

After the perturbative parton showering, all outgoing gluons are split non-perturbatively, into light quark-antiquark or diquark-antidiquark pairs (the default option is to disallow diquark splitting). At this point, each jet consists of a set of outgoing quarks and antiquarks (also possibly some diquarks and antidiquarks) and, in the case of spacelike jets, a single incoming valence quark or antiquark. The latter is replaced by an outgoing spectator carrying the opposite colour and the residual flavour and momentum of the corresponding beam hadron.

In the limit of a large number of colours, each final-state colour line can now be followed from a quark/anti-diquark to an antiquark/diquark with which it can form a colour-singlet cluster.7 By virtue of pre-confinement, these clusters have a distribution of mass and spatial size that peaks at low values, falls rapidly for large cluster masses and sizes, and is asymptotically independent of the hard subprocess type and scale.

The clusters thus formed are fragmented into hadrons. If a cluster is too light to decay into two hadrons, it is taken to represent the lightest single hadron of its flavour. Its mass is shifted to the appropriate value by an exchange of 4-momentum with a neighbouring cluster in the jet. Similarly, any diquark-antidiquark clusters with masses below threshold for decay into a baryon-antibaryon pair are shifted to the threshold via a transfer of 4-momentum to a neighbouring cluster.

Those clusters massive enough to decay into two hadrons, but below a fission threshold to be specified below, decay isotropically 8 into pairs of hadrons selected in the following way. A flavour f is chosen at random from among u, d, s, the six corresponding diquark flavour combinations, and c. For a cluster of flavour f1 f2_, this specifies the flavours f1 f_ and f f2_ of the decay products, which are then selected at random from tables of hadrons of those flavours. See section 7 for details of the hadrons included. The selected choice of decay products is accepted in proportion to the density of states (phase space times spin degeneracy) for that channel. Otherwise, f is rejected and the procedure is repeated.

The above method of selection for cluster decays is simple and fast but does not automatically satisfy constraints such as strong isospin symmetry. The decay rate into hadrons of a certain flavour depends on the average phase space for channels involving that flavour. Thus, for example, the existence of the h or h', with the same quark content as the p0, leads to a slight reduction of direct p0 production relative to p+ and p-. Quantitatively, the effect is too small to be observed even with the high statistics of the LEP1 data. However, the method can give rise to strange effects if the particle data tables are extended, and modifications to avoid this have been proposed [45].

In the decays of clusters to h or h', the parameter ETAMIX gives the h8/h0 mixing angle in degrees (default = -20). Rates are not very sensitive to its exact value, as the h'/h suppression is dominated by mass effects in the cluster model. See section 7 for more details.

A fraction of clusters have masses too high for isotropic two-body decay to be a reasonable ansatz, even though the cluster mass spectrum falls rapidly (faster than any power) at high masses. These are fragmented using an iterative fission model until the masses of the fission products fall below the fission threshold. In the fission model the produced flavour f is limited to u, d or s and the product clusters f1 f_ and f f2_ move in the directions of the original constituents f1 and f2_ in their c.m. frame. Thus the fission mechanism is not unlike string fragmentation [46].

In HERWIG there are three main fission parameters, CLMAX, CLPOW and PSPLT. The maximum cluster mass parameter CLMAX and CLPOW specify the fission threshold Mf according to the formula
MfCLPOW = CLMAXCLPOW+ (m1 + m2)CLPOW ,
where m1 and m2 are the quark mass parameters RMASS(i) for flavours f1 and f2 (see section 3.2). The parameter PSPLT specifies the mass spectrum of the produced clusters, which is taken to be MPSPLT within the allowed phase space. Provided the parameter CLMAX is not chosen too small (the default value is 3.35 GeV), the gross features of events are insensitive to the details of the fission model, since only a small fraction of clusters undergo fission. However, the production rates of high-pt or heavy particles (especially baryons) are affected, because they are sensitive to the tail of the cluster mass distribution. The default value of the power CLPOW is 2. Smaller values will increase the yield of heavier clusters (and hence of baryons) for heavy quarks, without affecting light quarks much. For example, the default value gives no b-baryons (for the default value of CLMAX) whereas CLPOW = 1.0 makes the ratio of b-baryons to b-hadrons about 1/4.

There is also a switch CLDIR for cluster decays. If CLDIR = 1 (the default) then a cluster that contains a `perturbative' quark, i.e. one coming from the perturbative stage of the event (the hard process or perturbative gluon splitting) `remembers' its direction. Thus when the cluster decays, the hadron carrying its flavour continues in the same direction (in the cluster c.m. frame) as the quark. This considerably hardens the spectrum of heavy hadrons, particularly of c- and b-flavoured hadrons. It also introduces a tendency for baryon-antibaryon pairs preferentially to align themselves with the event axis (the `TPC/2g string effect' [47]). CLDIR = 0 turns off this option, treating clusters containing quarks of perturbative and non-perturbative origin equivalently. In the CLDIR = 1 option, the parameter CLSMR (default = 0.0) allows for a gaussian smearing of the direction of the perturbative quark's momentum. The smearing is actually exponential in 1-cosq with mean CLSMR. Thus increasing CLSMR decorrelates the cluster decay from the initial quark direction.

The process of b-quark hadronization requires special treatment and the results obtained using HERWIG are still not fully satisfactory. Generally speaking, it is difficult to obtain a sufficiently hard B-hadron spectrum and the observed b-meson/b-baryon ratio. These depend not only on the perturbative subprocess and parton shower but also on non-perturbative issues such as the fraction of b-flavoured clusters that become a single B meson, the fractions that decay into a B meson and another meson, or into a b-baryon and an antibaryon, and the fraction that are split into more clusters. Thus the properties of b-jets depend on the parameters RMASS(5), CLMAX, CLPOW and PSPLT in a rather complicated way. In practice these parameters are tuned to global final-state properties and one needs extra parameters to describe b-jets.

A parameter B1LIM has therefore been introduced to allow clusters somewhat above the Bp threshold mass Mth to form a single B meson if
M < Mlim = (1+B1LIM) Mth .
The probability of such single-meson clustering is assumed to decrease linearly for Mth < M < Mlim. This has the effect of hardening the B spectrum if B1LIM is increased from the default value of zero. In addition, in version 6, the parameters PSPLT, CLDIR and CLSMR have been converted into two-dimensional arrays, with the first element controlling clusters that do not contain a b-quark and the second those that do. Thus tuning of b-fragmentation can now be performed separately from other flavours, by setting CLDIR(2) = 1 and varying PSPLT(2) and CLSMR(2). By reducing the value of PSPLT(2), further hardening of the B-hadron spectrum can be achieved.

3.7.2  Underlying soft event

In hadron-hadron and lepton-hadron collisions there are `beam clusters' containing the spectators from the incoming hadrons. In the formation of beam clusters, the colour connection between the spectators and the initial-state parton showers is cut by the forced emission of a soft quark-antiquark pair. The underlying soft event in a hard hadron-hadron collision is then assumed to be a soft collision between these two beam clusters. In a lepton-hadron collision the corresponding `soft hadronic remnant' is represented by a soft collision between the beam cluster and the adjacent cluster, i.e. the one produced by the forced emission mentioned above.

The model used for the underlying event is based on the minimum-bias pp_ event generator of the UA5 Collaboration [48], modified to make use of our cluster fragmentation algorithm. This model is explained in the following subsection.

Adding 10000 to the HERWIG process code IPROC suppresses the underlying event, in which case the beam clusters are simply fragmented like other clusters, without any soft collision. The parameter PRSOF enables one to produce an underlying event in only a fraction PRSOF of events (default = 1.0). Adding 10000 to IPROC is thus equivalent to setting PRSOF = 0.

A parameter BTCLM is available to users to adjust the mass parameter equivalent to PMBM1 (see below) in remnant cluster formation. Its default value, 1.0, is identical to previous versions. There is also an option for the special treatment of the splitting of clusters containing hadron (or photon) remnants. IOPREM = 0 gives the fragments a gaussian mass spectrum typical of soft processes. When IOPREM = 1 (default), the child containing the remnant is treated as before but the other cluster, containing a perturbative parton, is treated as a normal cluster, with mass spectrum MPSPLT.

Two special remnant `particles' have been defined: 'REMG ' with IDHW = 71, IDHEP = 98 and 'REMN ' with IDHW = 72, IDHEP = 99. These are remnant photons and nucleons respectively. They are identical to photons and nucleons, except that gluons are labelled as valence partons and, for the nucleon, valence quark distributions are set to zero. They are used by an external package for simulating multi-parton interactions, called JIMMY [49]. See section 9.7 for further details.

3.7.3  Minimum bias processes

The minimum-bias event generator of the UA5 Collaboration [48] starts from a parametrization of the pp_ inelastic charged multiplicity distribution as a negative binomial distribution. In HERWIG version 6, the relevant parameters are made available to the user for modification, the default values being the UA5 ones as used in previous versions. These parameters are given in table 5.
Name Description Default
     
PMBN1 a in n_ =asb+c 9.110
PMBN2 b in n_ =asb+c 0.115
PMBN3 c in n_ =asb+c -9.500
     
PMBK1 a in 1/k =alns+b 0.029
PMBK2 b in 1/k =alns+b -0.104
     
PMBM1 a in (M-m1-m2-a)e-bM 0.4
PMBM2 b in (M-m1-m2-a)e-bM 2.0
     
PMBP1 pt slope for d,u 5.2
PMBP2 pt slope for s,c 3.0
PMBP3 pt slope for qq 5.2
     

Table 5: Soft/min.bias parameters.

The first three parameters control the mean charged multiplicity n_ at c.m. energy s1/2 as indicated. The next two specify the parameter k in the negative binomial charged multiplicity distribution,
P(n) =
G(n+k)
n! G(k)
(
_
n
 
 /k)n
(1+
_
n
 
 /k)n+k
  .
The parameters PMBM1 and PMBM2 describe the distribution of cluster masses M in the soft collision. These soft clusters are generated using a flat rapidity distribution with gaussian shoulders. The transverse momentum distribution of soft clusters has the form
P(pt)µ ptexp æ
ç
ç
è		 -b ( pt2+M2 )
1/2
 
 
 ö
÷
÷
ø		  ,
where the slope parameter b depends as indicated on the flavour of the quark or diquark pair created when the cluster was produced. As an option, for underlying events the value of s1/2 used to choose the multiplicity n may be increased by a factor ENSOF to allow for an enhanced underlying activity in hard events. The actual charged multiplicity is taken to be n plus the sum of the moduli of the charges of the colliding hadrons or clusters.

There is now also an interface to the multiple-interaction model JIMMY [49]. For this purpose, several routines have been added or modified. New are HWHREM for identifying and cleaning up the beam remnants and HWHSCT to administer the extra scatters. Minor modifications to HWBGEN and HWSBRN suppress energy conservation errors when ISLENT = -1; HWSSPC has an improved approximation for remnant mass at high energies; and HWUPCM improves safety against negative square roots.

3.8  Spacetime structure

The space-time structure of events is available for all types of subprocess. The production vertex of each parton, cluster, unstable resonance and final-state particle is supplied in the VHEP array of /HEPEVT/. Set PRVTX = .TRUE. to include this information when printing the event record (120 column format). The units are: x,y,z in mm and t in mm/c. In the case of partons and clusters the production points are always given in a local coordinate system with its origin at the relevant hard subprocess. This helps to separate the fermi-scale partonic showers from millimetre-scale distances possible in particle decays, for example the partonic decays of heavy (c,b) hadrons. The vertices of hadrons produced in cluster decays are always corrected back into the laboratory coordinate system.

It is possible to vary the principal interaction point, assigned to the c.m. frame entry in /HEPEVT/ (with ISTHEP=103), by setting PIPSMR=.TRUE. The smearing is generated by the routine HWRPIP according to a triple gaussian given by parameters VIPWID(I) (I = 1,2,3 for x,y,z widths): the default values correspond to LEP1.

It is also possible to veto particle decays that would occur outside a specified volume by setting MAXDKL=.TRUE. Each putative decay is tested in HWDXLM and if the particle would have decayed outside the chosen volume it is frozen and labelled as final state. Using IOPDKL = 1,2 selects a cylindrical or spherical allowed region (centred about the origin): then parameters DXRCYL, DXZMAX or DXRSPH specify the dimensions of the region.

3.8.1  Particle decays

Lepton and hadron lifetimes (in seconds) are supplied in the array RLTIM. In the case of MSSM (s)particles, including Higgs states, RLTIM values are entered through the input files (see discussion in section 3.5.2). The lifetimes of heavy quarks (top and any hypothetical extra generations) and weak bosons (including the SM Higgs) are derived from their calculated or specified widths in HWUDKS, whilst light quarks and gluons are given an effective minimum width that acts as a lifetime cutoff - see below. All particles whose lifetimes are larger than PLTCUT are set stable.

The proper (i.e. rest-frame) time t* at which an unstable lepton or hadron decays is generated according to the exponential decay law with mean lifetime át*ñ=tºRLTIM: Prob(proper time> t*) = exp(-t*/t) . The laboratory-frame decay time t and distance travelled d are obtained by applying a boost: t=g t*, d=bg t* where b=v/c and g = 1/(1-b2)1/2. The production vertices of the daughter particles are then calculated by adding the distance travelled by the mother particle as given above to its production vertex. The mean lifetime t of a particle is set, taking into account its width and virtuality, by:
t(q2)=
_
h
 
 ( q2 )
1/2
 
 
( (q2-M2)2 + (G q2/M)2 )
1/2
 
 
  .     (2)
This formula is used for all particles: light partons; heavy quarks and weak bosons, which have appreciable widths; resonances; and unstable leptons. It interpolates between t=h_/G for a particle that is on mass-shell and t=h_(q2)1/2/(q2-M2) for one that is far off mass-shell.

3.8.2  Parton showers

The above prescription, based on an exponential proper lifetime distribution, is also used to describe parton showers. For light quarks and gluons, whose natural widths are small, this could lead to unreasonably large distances being generated in the final, low virtuality steps of showering. To avoid this they are given a width G=VMIN2/M; the parameter VMIN2 (default value 0.1 GeV2) acts effectively as a lower limit on a parton's virtuality. This is particularly important for the forced splitting of gluons (see section 3.7.1), which uses t=h_ RMASS(13)/VMIN2.

3.8.3  Hadronization

In the case of a cluster its initial production vertex is taken as the midpoint of a line perpendicular to the cluster's direction of travel and with its two ends on the trajectories of the constituent quark-antiquark pair. If such a cluster undergoes a forced splitting to two clusters the string picture is adopted. The vertex of the light quark pair is positioned so that the masses of the two daughter clusters would be the same as those for two equivalent string fragments. The production vertices of the daughter clusters are given by the first crossing of their constituent qq_ pairs. The production positions of primary hadrons from cluster decays are smeared, relative to the cluster position, according to a gaussian distribution of width 1/(cluster mass).

3.8.4  Colour rearrangement

HERWIG version 6 contains a colour rearrangement model based on the space-time structure of an event at the end of the parton shower. This is illustrated in the simple example shown below where showering results in a colour-neutral qggq_ final state. In the conventional HERWIG hadronization model (corresponding to the default value of the reconnection parameter, CLRECO=.FALSE.), after a non-perturbative splitting of the final-state gluons, colour singlet clusters are formed from neighbouring qq_ pairs: (ij)(pq)(kl). However when CLRECO=.TRUE. the program first creates colour singlet clusters as normal but then checks all (non-neighbouring) pairs of clusters to test if a colour rearrangement lowers the sum of the clusters' spatial sizes added in quadrature. A cluster's size dij is defined to be the Lorentz-invariant space-time distance between the production points of its constituent quark qi and antiquark qj_. If an allowed alternative is found, that is, (ij)(kl)®(il)(jk) such that |dij|2+|dkl|2 > |dil|2+|dkl|2, then it is accepted with a probability given by the parameter PRECO (default value 1/9).



Note that not all colour rearrangements are allowed, for instance in the example shown (ij)(pq)® (iq)(jp) is forbidden since the cluster (jp) is a colour octet - it contains both products from a non-perturbative gluon splitting.

Multiple colour rearrangements are considered by the program, as are those between clusters in jets arising from a single, colour neutral source, for example Z0 decay (as shown above), or due to more than one source, for example e+e-® W+W- ® 4 jets. In the latter case a new parameter, EXAG, is available to exaggerate the lifetime of the W± or any other weak boson, so that any dependence of rearrangement effects on source separation can be investigated.

The CLRECO option can be used for all the processes available in HERWIG. Note, however, that before using the program with CLRECO=.TRUE. for detailed physics analyses the default parameters should be retuned to LEP data with this option switched on.

3.8.5   B-B mixing

When MIXING=.TRUE., particle-antiparticle mixing for Bd,s0 mesons is implemented. The probability that a meson is mixed when it decays is given in terms of its lab-frame decay time t by:
Pmix(t) =
1
2
 -
cos(Xtm/ct E)
2cosh(Ytm/ct E)
  ,
where X=D M/G, Y=DG/2G and m, t, E are the B0 mass, lifetime and energy. The ratios X and Y are stored in XMIX(I) and YMIX(I), I=1,2 for q=s,d. Whenever a neutral B meson occurs in an event, a copy of the original entry is always added to the event record, with ISTHEP=200, which gives the particle's flavour at the production (or cluster decay) time. This is in addition to the usual decaying particle entry with ISTHEP=199.

4  Processes

4.1  Beams

Name Description Default
PART1 Type of particle in beam 1 'PBAR '
PART2 Type of particle in beam 2 'P '
PBEAM1 Momentum of beam 1 900.0
PBEAM2 Momentum of beam 2 900.0
IPROC Type of process to generate 1500
MAXEV Number of events to generate    100

Table 6: Main program variables.

As indicated in table 6, a number of variables must be set in the main program HWIGPR to specify what is to be simulated. The beam particle types PART1, PART2 can take any of the values NAME listed in table 7.
  NAME   NAME
e+ 'E+ ' e- 'E- '
µ+ 'MU+ ' µ- 'MU- '
ne 'NU_E ' ne_ 'NU_EBAR '
nµ 'NU_MU ' nµ_ 'NU_MUBAR'
nt 'NU_TAU ' nt_ 'NU_TAUBR'
p 'P ' p_ 'PBAR '
n 'N ' n_ 'NBAR '
p+ 'PI+ ' p- 'PI- '
g 'GAMMA '    

Table 7: Beam particles.

In the case of point-like photon/QCD processes, IPROC = 5000-5999, the first particle must be the photon or a lepton. In addition, beams 'K+ ' and 'K- ' are supported for minimum bias non-diffractive soft hadronic events (IPROC=8000) only. In the case that the beam momenta PBEAM1 and PBEAM2 are not equal, the default procedure (USECMF=.TRUE.) is to generate events in the beam-beam centre-of-mass frame and boost them back to the laboratory frame afterwards.

In hadronic processes with lepton beams (e.g. photoproduction in ep), the lepton ® lepton + photon vertex uses the full transverse-momentum dependent splitting function, with exact light-cone kinematics, i.e the Equivalent Photon Approximation (EPA). This means that the photon-hadron collision has a transverse momentum in the lepton-hadron frame and must be boosted to a frame where it has no transverse momentum. Thus the c.m.f. boost described above is always used in these processes, regardless of the value of USECMF. The correct lower energy cutoff appropriate to the hadronic process is applied to the photon. The Q2 of the photon is generated within the kinematically allowed limits, or the user-defined limits Q2WWMN and Q2WWMX (defaults 0 and 4) whichever is more restrictive.9 Similarly for the photon's light-cone momentum fraction, with user-defined limits YWWMIN and YWWMAX (default 0 and 1). Together with the Bjorken y-variable limits YBMIN and YBMAX, this allows different ranges for the tagged and untagged photons in two-photon DIS.

4.1.1  Parton distributions

The parton momentum fraction distributions of the beam particles are used in the generation of initial-state parton showers and also in the non-perturbative process of linking the shower with the beam hadron and its remnant. Since the parton showering is done in leading-logarithmic order, there is no strong motivation to use next-to-leading order parton distributions, although this has become customary since the most up-to-date distributions are deduced from next-to-leading order fits to (inclusive) data. Thus the most common option is to use the interface to the PDFLIB parton distribution library [50].

The HERWIG interface is compatible with PDFLIB version 4. AUTPDF should be set to the author group as listed in the PDFLIB manual, e.g. 'MRS', and MODPDF to the set number in the new convention. It is permissible to choose the PDFLIB set independently for each of the two beams. For example, to use MRS D- for the proton and Gordon-Storrow set 1 for the photon in g-hadron or lepton-hadron collisions, one sets: 
          AUTPDF(2)='MRS'
          MODPDF(2)=28
          AUTPDF(1)='GS'
          MODPDF(1)=2
If the PDFLIB interface is not used, the parton distributions are chosen from the HERWIG internal sets according to the value of the parameter NSTRU. The default parton distributions in HERWIG versions prior to 6.3 were very old and did not include fits to any of the HERA data. Therefore several new PDFs have been included in versions 6.3 and higher. These are shown in table 8.

It should be noted that we have only added leading-order fits because the evolution algorithms in HERWIG, in particular the backward-evolution algorithm for initial-state parton showering, are only leading-order and therefore inconsistencies could occur with next-to-leading-order distributions.

The new default structure function set NSTRU=8 is the average of two of the published fits [51], because this has been found [52] to be closer to the central value of more recent next-to-leading-order fits. The other fits can then be used to assess the effects of varying the high-x gluon.

These new HERWIG parton distributions are only available for nucleons. For pion beams, either the old NSTRU=1,2 pion sets or PDFLIB should be used.

For photons, the default is to use the Drees-Grassie parton distributions [53]. The heavy quark content of the photon uses the corrections to the Drees-Grassie distribution functions for light quarks, calculated by Drees and Kim [54]. There is also an interface to the Schuler-Sjöstrand [55] parton distribution functions for the photon, version 2. These appear as PDFLIB sets with author group `SaSph', but are actually implemented via a call to their SASGAM code. The value in MODPDF specifies the set (1-4 for 1D [recommended set],1M, 2D,2M), whether the Bethe-Heitler process is used for heavy flavours (add 10), whether the P2-dependence is included (add 20), and which of their P2 models is used (add 100 times their IP2 parameter).
NSTRU Description
6 Central aS and gluon leading-order fit of [51]
7 Higher gluon leading-order fit of [51]
8 Average of central and higher gluon leading-order fits of [51]

Table 8: New internal MRST parton distributions.

An option to damp the parton distributions of off mass-shell photons relative to on-shell photons, according to the scheme of Drees and Godbole [56] has been introduced. The adjustable parameter PHOMAS defines the crossover from the non-suppressed to suppressed regimes. Recommended values lie in the range from QCDLAM to 1 GeV. The default value PHOMAS = 0 corresponds to no suppression, as in previous versions.

4.2  Summary of subprocesses

We give in table 9 a list of the currently available hard subprocesses IPROC. More detailed descriptions are given in sections 4.3-4.12, and then in section 4.13 there are instructions to users on how to add a new process.
IPROC Process
100 l+ l- ® q q_(g) (all q flavours)
100+IQ l+ l- ® q q_(g) (IQ=1,2,3,4,5,6 for q=d,u,s,c,b,t)
107 l+ l- ® g g (g) (fictitious process)
110 l+ l- ® q q_ g (all flavours)
110+IQ l+ l- ® q q_ g (IQ as above)
120 l+ l- ® q q_ (all flavours, no hard gluon correction)
120+IQ l+ l- ® q q_ (IQ as above, no hard gluon correction)
127 l+ l- ® g g (fictitious process, no hard gluon correction)
150+IL l+ l- ® l' l_' (IL=1,2,3 for l'=e,µ,t, N.B. l¹l')
200 l+ l- ® W+ W- (see sect. 4.3.2 on control of W/Z decays)
250 l+ l- ® Z0 Z0 (see sect. 4.3.2 on control of W/Z decays)
300 l+ l- ® Z0HSM0 ® Z0 q q_ (all flavours)
300+IQ l+ l- ® Z0HSM0 ® Z0 q q_ (IQ as above)
306+IL l+ l- ® Z0HSM0 ® Z0 l l_ (IL as above)
310, 311 l+ l- ® Z0HSM0 ® Z0 W+ W-, Z0 Z0 Z0
312 l+ l- ® Z0HSM0 ® Z0 g g
399 l+ l- ® Z0HSM0 ® Z0 anything
400+ID l+ l- ® n n_HSM0 + l+ l-HSM0 (ID as in IPROC=300+ID)
500+ID
l+ l- ® l+l-gg® l+l -q
_
q
 
 /l
_
l
 
 /W+W-
(ID=0-10 as in IPROC=300+ID)
550+ID l+ l- ®lnl g W® lnl qq_'/ll_' (ID=0-9 as in IPROC=300+ID)
600 l+ l-® qq_ gg, qq_ q'q_' (all q flavours)
600+IQ l+ l-® qq_ gg, qq_ q'q_' (IQ as above)
  After generation, IHPRO is subprocess (see sect. 4.3.5)
700-99 Minimal Supersymmetric Standard Model (MSSM) processes
700 l+ l- ® 2-sparticle processes (sum of 710, 730, 740 and 760)
710 l+ l- ® neutralino pairs (all neutralinos)
706+4IN1+IN2 l+ l- ® c~IN10 c~IN20 (IN1,2=neutralino mass eigenstate)
730 l+ l- ® chargino pairs (all charginos)
728+2IC1+IC2 l+ l- ® c~IC1+ c~IC2- (IC1,2=chargino mass eigenstate)
740 l+ l- ® slepton pairs (all flavours)
736+5IL l+ l- ® lL,R~ lL,R~* (IL=1,2,3 for l~=e~,µ~,t~)
737+5IL l+ l- ® lL~ lL~* (IL as above)
738+5IL l+ l- ® lL~ lR~* (IL as above)
739+5IL l+ l- ® lR~ lR~* (IL as above)
740+5IL l+ l- ® nL~ nL~* (IL=1,2,3 for ne~, nµ~, nt~)
760 l+ l- ® squark pairs (all flavours)
757+4IQ l+ l- ® qL,R~ q~L,R* (IQ=1,2,3,4,5,6 for q~=d~,u~,s~, c~,b~,t~)
758+4IQ l+ l- ® qL~ q~L* (IQ as above)
759+4IQ l+ l- ® qL~ q~R* (IQ as above)
760+4IQ l+ l- ® qR~ q~R* (IQ as above)
800-99 R-parity violating supersymmetric processes
800 Single sparticle production, sum of 810-840
810 l+ l- ® c~0 ni, (all neutralinos)
810+IN l+ l- ® c~IN0 ni, (IN=neutralino mass state)
820 l+ l- ® c~- ei+ (all charginos)
820+IC l+ l- ® c~IC- ei+, (IC=chargino mass state)
830 l+ l- ® ni~ Z0 and l+ l- ® l~i+ W-
840 l+ l- ® ni~ h0/H0/A0 and l+ l- ® l~i+ H-
850 l+ l- ® ni~ g
860 Sum of 870 and 880
870 l+ l- ® l+ l-, via LLE only
867+3IL1+IL2 l+ l- ® lIL1+ lIL2- (IL1,2=1,2,3 for e,µ,t)
880 l+ l- ® d_ d, via LLE and LQD
877+3IQ1+IQ2 l+ l- ® dIL1 dIL2_ (IQ1,2=1,2,3 for d,s,b)
910 l+ l- ® ne ne_ h0 + e+ e- h0
920 l+ l- ® ne ne_ H0 + e+ e- H0
960 l+ l- ® Z0 h0
970 l+ l- ® Z0 H0
955 l+ l- ® H+ H-
965 l+ l- ® A0 h0
965 l+ l- ® A0 H0
1000+ID l+l-® t t_ HSM0 (ID as in IPROC=300+ID)
1110+IQ l+l-® q q_  h0 (IQ as in IPROC=100+IQ)
1116+IL l+l-® l+l- h0 (IL=1,2,3 for e,µ,t)
1120+IQ l+l-® q q_  H0 (IQ as in IPROC=100+IQ)
1126+IL l+l-® l+l- H0 (IL=1,2,3 for e,µ,t)
1130+IQ l+l-® q q_  A0 (IQ as in IPROC=100+IQ)
1136+IL l+l-® l+l- A0 (IL=1,2,3 for e,µ,t)
1140 l+l-® d u_  H+ + ch. conj.
1141 l+l-® s c_  H+ + ch. conj.
1142 l+l-® b t_  H+ + ch. conj.
1143 l+l-® e ne_ H+ + ch. conj.
1144 l+l-® µ nµ_H+ + ch. conj.
1145 l+l-® t nt_H+ + ch. conj.
1200-99 Reserved for other l+l- processes
1300 q q_ ® Z0/g ® q' q_' (all flavours)
1300+IQ q q_ ® Z0/g ® q' q_' (IQ=1,2,3,4,5,6 for q=d,u,s,c,b,t)
1350 q q_ ® Z0/g ® l l_ (all lepton species)
1350+IL q q_ ® Z0/g ® l l_ (IL=1-6 for l=e,ne,µ,nµ, etc.)
1399 q q_ ® Z0/g ® anything
1400 q q_ ® W±® q' q_'' (all flavours)
1400+IQ q q_ ® W±® q' q_'' (q' or q'' as above)
1450 q q_ ® W±® l nl (all lepton species)
1450+IL q q_ ® W±® l nl (IL=1,2,3 for l=e,µ,t)
1499 q q_ ® W±® anything
1500 QCD 2 ® 2 hard parton scattering
  After generation, IHPRO is subprocess (see sect. 4.6.2)
1600+ID g g/qq_ ®HSM0 (ID as in IPROC=300+ID)
1700+IQ QCD heavy quark production (IQ as above)
  After generation, IHPRO is subprocess (see sect. 4.6.2)
1800 QCD direct photon + jet production
  After generation, IHPRO is subprocess (see sect. 4.6.5)
1900+ID qq_® q'q_'W+W-/Z0Z0® q'q_'HSM0 (ID as in IPROC=300+ID)
2000 t production via W± exchange (sum of 2001-2008)
2001-4 u_ b_ ® d_ t_ ,    d b_ ® u t_ ,   d_ b_ ® u_ t_ ,    u b® d t
2005-8 c_ b_ ® s_ t_ ,    sb_ ® ct_  ,   s_ b ® c_ t ,    c b® s t
2100 W± + jet production
2110 W± + jet production (Compton only: g q ® W q)
2120 W± + jet production (annihilation only: q q_ ® W g)
2150 Z0 + jet production
2160 Z0 + jet production (Compton only: g q ® Z q)
2170 Z0 + jet production (annihilation only: q q_ ® Z g)
2200 QCD direct photon pair production
  After generation, IHPRO is subprocess (see sect. 4.6.5)
2300+ID QCD SM Higgs + jet production (ID as in IPROC=300+ID)
  After generation, IHPRO is subprocess (see sect. 4.6.10)
2400 Mueller-Tang colour singlet exchange
2450 Quark scattering via photon exchange
2500+ID gg/qq_® tt_HSM0 (ID as in IPROC=300+ID)
2600+ID qq_' ® W±HSM0 (ID as in IPROC=300+ID)
2700+ID qq_ ® Z0HSM0 (ID as in IPROC=300+ID)
2800 W+W- production in hadron-hadron collisions
2810 Z0Z0 production in hadron-hadron collisions (including photon terms)
2815 Z0Z0 production in hadron-hadron collisions (Z0 only)
2820 W±Z0 production in hadron-hadron collisions (including photon terms)
2825 W±Z0 production in hadron-hadron collisions (Z0 only)
2850 hadron-hadron ® W+ W- X using MC@NLO
2860 hadron-hadron ® Z0 Z0 X using MC@NLO
2870 hadron-hadron ® W+ Z0 X using MC@NLO
2880 hadron-hadron ® W- Z0 X using MC@NLO
2900+IQ gg+qq_® QQ_ Z0 for massless Q and Q_ (IQ=1...6 for Q=d... t)
2910+IQ gg+qq_® QQ_ Z0, for massive Q and Q_ (IQ=1...6 for Q=d... t)
3000-3999 Minimal Supersymmetric Standard Model (MSSM) processes
3000 2-parton ® 2-sparticle processes (sum of those below)
3010 2-parton ® 2-sparton processes
3020 2-parton ® 2-gaugino processes
3030 2-parton ® 2-slepton processes
3100+ISQ gg/qq_® q~q~'* H± (ISQ=IPROC-3100 as from table 15)
3200+ISQ gg/qq_® q~q~'* h,H,A (ISQ=IPROC-3200 as from table 16)
3310,3315
q
_
q
 
 ' ® W±h0,H±h0 (all
q,q'flavours - gauge bosons mediated only)
3320,3325 qq_' ® W±H0,H±H0 ('')
       3335 qq_' ® H±A0 ('')
3350        qq_ ® W±H± (Higgstrahlung and Higgs mediated)
       3355 qq_ ® H±H± (all q flavours - gauge boson mediated only)
3360,3365 qq_ ® Z0 h0,A0 h0 ('')
3370,3375 qq_ ® Z0 H0,A0 H0 ('')
3410 bg ® b h0 + ch. conj.
3420 bg ® b H0 + ch. conj.
3430 bg ® b A0 + ch. conj.
3450 bg ® t H- + ch. conj.
3500 b q ® b q' H± + ch. conj.
3610 qq_/gg ® h0 (light scalar Higgs)
3620 qq_/gg ® H0  (heavy scalar Higgs)
3630 qq_/gg ® A0  (pseudoscalar Higgs)
3710 qq_® q'q_'W+W-/Z0Z0® q'q_'h0
3720 qq_® q'q_'W+W-/Z0Z0® q'q_'H0
3810+IQ gg+qq_® QQ_ h0 (all q flavours in s-channel, IQ as usual for Q flavour)
3820+IQ gg+qq_® QQ_ H0 ('')
3830+IQ gg+qq_® QQ_ A0 ('')
3839        gg+qq_® bt_ H+ + ch. conjg. (all q flavours in s-channel)
3840+IQ gg ® QQ_ h0 (IQ as above)
3850+IQ gg ® QQ_ H0 ('')
3860+IQ gg ® QQ_ A0 ('')
3869        gg ® bt_ H+ + ch. conjg.
3870+IQ qq_ ® QQ_ h0 (all q flavours in s-channel, IQ as above)
3880+IQ qq_ ® QQ_ H0 ('')
3890+IQ qq_ ® QQ_ A0 ('')
3899        qq_ ® bt_ H+ + ch. conjg. (all q flavours in s-channel)
3900-99 Reserved for other hadron-hadron MSSM processes
4000-99 R-parity violating supersymmetric processes via LQD
4000 single sparticle production, sum of 4010-4050
4010 uj_ dk ® c~0 li-, dj_ dk ® c~0 ni (all neutralinos)
4010+IN uj_ dk ® c~IN0 li-, dj_ dk ® c~IN0 ni (IN=neutralino mass state)
4020 uj_ dk ® c~- ni, dj_ dk ® c~- ei+ (all charginos)
4020+IC uj_ dk ® c~IC- ni, dj_ dk ® c~IC- ei+ (IC=chargino mass state)
4040 uj dk_ ® t~i+ Z0, uj dk_ ® ni~ W+ and dj dk_ ® l~i+ W-
4050 uj dk_ ® l~i+ h0/H0/A0, uj dk_ ® ni~ H+ and dj dk_ ® l~i+ H-
4060 Sum of 4070 and 4080
4070 uj_ dk ® ul_ dm and dj_ dk ® dl_ dm , via LQD only
4080 uj_ dk ® nj lk- and dj_ dk ® lj+ lk- , via LQD and LLE
4100-99 R-parity violating supersymmetric processes via UDD
4100 single sparticle production, sum of 4110-4150
4110 ui dj ® c~0 dk_, dj dk ® c~0 ui_ (all neutralinos)
4110 +IN ui dj ® c~IN0 dk_, dj dk ® c~IN0 ui_(IN as above)
4120 ui dj ® c~+ uk_, dj dk ® c~- di_ (all charginos)
4120 +IC ui dj ® c~IC+ uk_, dj dk ® c~IC- di_ (IC as above)
4130 ui dj ® g~ dk_, dj dk ® g~ ui_
4140 ui dj ® b~1* Z0, dj dk ® t~1* Z0, ui dj ® t~i* W+ and dj dk ® b~i* W-
4150
ui dj ®
~
d
 
 k1* h0/H0/A0, dj dk ®
~
u
 
 i1* h0/H0/A0, ui dj ®
~
u
 
 ka* H+,
dj dk ®
~
d
 
 ia* H-
4160 ui dj ® ul dm, dj dk ® dl dm via UDD.
4200-99 Graviton resonance production
4200 Sum of 4210, 4250 and 4270
4210 gg/qq_® G ® gg/qq_ (all partons)
4210+IQ gg/qq_® G ® qq_ (IQ as above)
4220 gg/qq_® G ® gg
4250 gg/qq_® G ® ll_ (all leptons)
4250+IL gg/qq_® G ® ll_ (IL=1-6 for l=e,ne,µ,nµ, etc.)
4260 gg/qq_® G ® g g
4270 gg/qq_® G ® W+W-/Z0Z0/HSM0HSM0
4271 gg/qq_® G ® W+W-
4272 gg/qq_® G ® Z0Z0
4273 gg/qq_® G ® HSM0HSM0
5000 Pointlike photon-hadron jet production (all flavours)
5100+IQ Pointlike photon heavy flavour pair production (IQ as above)
5200+IQ Pointlike photon heavy flavour single excitation (IQ as above)
  After generation, IHPRO is subprocess (see sect. 4.6.5)
5300 Quark-photon Compton scattering
5500 Pointlike photon production of light (u,d,s) L=0 mesons
5510,20 S=0 mesons only, S=1 mesons only
  After generation, IHPRO is subprocess (see sect. 4.6.5)
6000 gg® qq_ (all flavours)
6000+IQ gg® qq_ (IQ as above)
6006+IL gg® ll_ (IL=1,2,3 for l=e,µ,t)
6010 gg® W+W-
7000 - Baryon-number violating and other multi-W± processes
7999 generated by HERBVI package
8000 Minimum bias soft hadron-hadron event
9000 Deep inelastic lepton scattering (all neutral current)
9000+IQ Deep inelastic lepton scattering (NC on flavour IQ)
9010 Deep inelastic lepton scattering (all charged current)
9010+IQ Deep inelastic lepton scattering (CC on flavour IQ)
9100 Boson-gluon fusion in neutral current DIS (all flavours)
9100+IQ Boson-gluon fusion in neutral current DIS (IQ as above)
9107 J/y + gluon production by boson-gluon fusion
9110 QCD Compton process in neutral current DIS (all flavours)
9110+IP QCD Compton process in NC DIS (IP=1-12 for d-t,d_-t_)
9130 All O(aS) NC processes (i.e. 9100+9110)
9140+IP Heavy quark production by charged-current boson-gluon fusion
  IP: 1 = s c_, 2 = b c_, 3 = s t_, 4 = b t_ (+ ch. conj.)
9500+ID W+W-/Z0Z0®HSM0 in DIS (ID as in IPROC=300+ID)
10000+IP as IPROC=IP but with soft underlying event
  (soft remnant fragmentation in lepton-hadron) suppressed

Table 9: Process codes.

4.2.1  Treatment of quark masses

The extent to which quark mass effects are included in the hard process cross section is different in different processes. In many processes, they are always treated as massless: IPROC = 1300, 1800, 1900, 2100, 2300, 2400, 5300, 9000. In two processes they are all treated as massless except the top quark, for which the mass is correctly incorporated: 1400, 2000. In the case of massless pair production, only quark flavours that are kinematically allowed are produced. In all cases the event kinematics incorporate the quark mass, even when it is not used to calculate the cross section. In two processes, quarks are always treated as massive: 500, 9100. Finally, in several processes, the behaviour is different depending on whether a specific quark flavour is requested, in which case its mass is included, or not, in which case all quarks are treated as massless. These are: IPROC = 100, 110, 120, QCD 2® 2 scattering (1500 vs. 1700+IQ), jets in direct photoproduction (5000 vs. 5100+IQ and 5200+IQ). In the case of IPROC = 2900,2910 one has the option of using the massless or massive matrix element.

These differences can cause inconsistencies between different ways of generating the same process. The most noticeable example is in direct photoproduction, where one can use process 9130, which uses the exact 2® 3 matrix element e+g ® e+q+q_, or process 5000, which uses the Equivalent Photon Approximation (EPA) for e ® e+g and the 2® 2 matrix element for g+g ® q+q_. For typical HERA kinematics, the EPA is valid to a few per cent, but the difference between the two processes is much larger, about 20% for PTMIN=2 GeV. This is entirely due to the difference in quark mass treatments, as can be checked by comparing process 9130 with processes 5100+IQ and 5200+IQ summed over IQ.

4.2.2  Couplings

The two-loop QCD coupling at scale Q is given by subroutine HWUALF with arguments IOPT=1 and SCALE=Q. Threshold matching is performed at the quark mass scales Q=RMASS(i). Setting IOPT=0 initialises the coupling using the 5-flavour value LMS=QCDLAM. Other values of IOPT are for internal use only.

The electromagnetic coupling is given by HWUAEM(Q2)=e2/(4p); it runs according to the prescription in ref. [57] with the hadronic term as given in ref. [58]. The parameter ALPHEMºHWUAEM(0), default value 0.0072993, provides the normalization at the Thomson limit; it is used for all processes involving real photons. Photon emission in parton showers and in the `dead-zone' in e+e- processes can be enhanced by a factor of ALPFAC (default = 1). The normalised electric charges of the fundamental fermions are stored in the array QFCH(I), where I = 1-6 for the quarks d,u,s,c,b,t (e.g. QFCH(4)=2./3.) and 11-16 for the leptons e,ne,µ,nµ, t,nt.

The weak neutral current is taken to be of the form e(vf+afg5)gµ, where the electric charge is evaluated at a scale appropriate to the process. The arrays VFCH(I,J) and AFCH(I,J) store the couplings: I as before, J = 1 for the minimal Standard Model and 2 for possible Z' couplings (only used if ZPRIME=.TRUE.). Note that universality is not assumed - couplings can be arbitrarily set separately for each fermion species. The default couplings are given in terms of of SWEIN=sin2qW, default value 0.2319, as:
vf=(T3/2-Qsin2qW)/(cosqWsinqW) ,      af=T3/(2cosqWsinqW) .

The weak charged current is given in terms of g=e/sinqW and the Cabbibo-Kobayashi-Maskawa mixing matrix, the elements squared of which are stored in VCKM(K,L), K=1,2,3 for u,c,t, L=1,2,3 for d,s,b. The variable SCABI=sin2qCabibbo is however also retained for the present. Note the Fermi constant GFermi is eliminated from all cross sections.

The overall scale for all cross sections, given in nanobarns, is set by GEV2NB=(h_ c/e)2, default value 389379.

We now give more detailed descriptions of the various subprocesses, concentrating again on the new features since ref. [1].

4.3  Lepton-antilepton Standard Model processes

Lepton beam polarization effects are included in e+e-® 2/3 jet production and the Bjorken process (ZH production). Incoming lepton and antilepton beam polarizations are specified by setting the two 3-vectors EPOLN and PPOLN: component 3 is longitudinal and 1,2 transverse.

Photon initial-state radiation (ISR) in e+e- annihilation events is allowed. The parameter TMNISR sets the minimum s^/s value (default = 10-4), ZMXISR sets the (arbitrary) separation between unresolved and resolved emission (default = 1-10-6). Setting ZMXISR = 0 switches off photon ISR.

Where processes are listed for l+ l- they are available for e+e- and µ+µ- annihilation. Many of the processes listed for e+e- will also work for µ+µ-, but we have not been systematic in ensuring this.

4.3.1  IPROC=100-127: hadron production

A correction to hard gluon emission in e+e- events has been added and is now the default process for IPROC=100+IQ. The O(aS) matrix element is used to add events in the `dead zone' of phase-space corresponding to a quark-antiquark pair recoiling from a hard gluon [16]. Although this is asymptotically negligible, and cannot be produced within the shower itself, it has a significant effect at LEP1 energies. As a result, the default parameters have been retuned, and show a marked improvement in agreement with e+e- data for event shapes sensitive to three-jet configurations.

The routine HWBDED implements this hard correction while HWBRAN has been modified to include the soft matrix-element corrections described in section 3.2.3.

When IPROC=100+IQ, hard gluons emitted into the dead zone are assigned to the quark or antiquark shower and do not appear explicitly in the event record.

The qq_ g process alone, generated according to the O(aS) qq_ g matrix element with a maximum thrust cutoff THMAX (default 0.9), is given by IPROC=110+IQ.

The uncorrected qq_ process has been retained for comparative purposes and is available as IPROC=120+IQ.

The fictional e+e- processes e+e- ® g+g(+g), IPROC=107 and 127, is treated just like e+e- ® qq_, summed over light quark flavours, for direct comparisons between quark and gluon jets.

4.3.2  IPROC=150-250: lepton and electroweak boson production

In IPROC=150, only the s-channel process, mediated by a virtual photon or Z0, is included, so the final-state leptons must be different from the initial ones.

The processes of W+W- and Z0Z0 pair production, IPROC=200 and 250, are based on a program kindly supplied by Zoltan Kunszt, which fully includes decay correlations. The QCD O(aS) matrix element correction for hard gluon emission in hadronic W and Z decays has also been implemented in these processes, according to the method described in section 3.2.3. In contrast to IPROC=100+IQ, any hard gluons emitted into the dead zone are shown explicitly in the event record.

4.3.3  IPROC=300-499: Higgs boson production

HERWIG generates SM Higgs bosons in lepton-antilepton collisions through the Bjorken process Z(*) ® Z(*)HSM0 with one or both Z0's off-shell (IPROC=300+ID) and W+W-/Z0Z0 fusion (IPROC=400+ID). See section 3.4 for explanation of how the Higgs decay is controlled by the value of ID.

4.3.4  IPROC=500-559: two-photon/photon-boson processes

In the e+e- two-photon processes, IPROC = 500+ID, ID=0-10 is the same as in Higgs processes for qq_, it ll_ and W+W-. The Equivalent Photon Approximation (EPA) is used for the e® eg vertices. The phase space is controlled by EMMIN and EMMAX for the two-photon centre-of-mass frame (CMF) mass, PTMIN and PTMAX for the transverse momentum of the CMF in the lab, and CTMAX for the c.m. angle of the outgoing particles. The additional phase-space variable WHMIN sets the minimum allowed hadronic mass and affects photoproduction reactions (g-hadron and g-g) and DIS.

In photon-W± fusion, IPROC = 550+ID, ID=0-9 is also the same as in Higgs processes, except that ID=1 or 2 both give the sum of du_ and ud_ etc. The EPA is used for the e® eg vertex. The phase space is controlled by EMMIN and EMMAX only. The full 2® 3 matrix elements for g e® ff_'n are used, so the cross section for real W± production is correctly included. In the case of gg® W W the decay correlations are not yet correctly included: the W's currently decay isotropically.

4.3.5  IPROC=600-656: four jet production

Electron-positron annihilation to four jets is provided by IPROC=600+IQ, where a non-zero value for IQ guarantees production of quark flavour IQ whilst IQ=0 corresponds to the natural flavour mix. IPROC=650+IQ is as above but without those terms in the matrix element which orient the event w.r.t. the lepton beam direction. The matrix elements are based on those of Ellis, Ross and Terrano [59] with orientation terms from Catani and Seymour [60]. The soft and collinear divergences are avoided by imposing a minimum ycut, Y4JT (default 0.01), on the initial four partons. The interjet distance ycut is calculated using either the Durham or JADE metrics, as selected by the logical variable DURHAM (default .TRUE.). In order to improve efficiency parameterizations of the volume of four-body phase space are used: these are accurate up to a few percent for ycut values less than 0.14. Note also that the phase space is for massless partons, as are the matrix elements, though a mass threshold cut is applied.

The argument of the strong coupling is set equal to EMSCA, the scale for the parton showers. This is taken to be the smallest of the ycut values times the c.m. energy squared if FIX4JT=.FALSE. (default); otherwise the argument is fixed at Y4JT times the c.m. energy squared.
IHPRO g*® 1+2+3+4 c/f conn.
91 q+q_+g+g 3 1 4 2
92 q+q_+g+g 4 1 2 3
93 q+q_+q+q_ 4 1 2 3
94 q+q_+q+q_ 2 1 4 3
95 q+q_+q'+q_' 4 1 2 3

Table 10: Four jet subprocesses.

The matrix elements for the qq_ gg and qq_ qq_ (same flavour quark) final states receive contributions from two colour flows each. The actual process and colour flow generated is indicated by IHPRO as shown in table 10. The meaning of `c/f conn.' is discussed in section 4.6.2 below. The treatment of the interference terms between the two colour flows is controlled by the array IOP4JT(1) for qq_ gg and IOP4JT(2) for qq_ qq_ (identical quark flavour):
IOP4JT(1)= ì
í
î
0 neglect
1 extreme 3142
2 extreme 4123
 IOP4JT(2)= ì
í
î
0 neglect
1 extreme 4123
2 extreme 2143
In both instances the default value is 0.

See ref. [61] for some applications and discussions of the new four-jet implementation in e+e- annihilations.

4.4  Lepton-antilepton supersymmetric processes (MSSM)

The R-parity conserving lepton-antilepton SUSY processes have IPROC=700-799 and ¬  Rp processes have IPROC=800-899. Lepton beam polarization effects are included for all of the SUSY production processes. The processes have all been implemented in such a way as to allow either e+e- or µ+µ- as the initial state. As with the SM lepton-antilepton processes, ISR is allowed for the SUSY production processes.

As, by probing the individual thresholds, it may be possible to study the production of a given sparticle pair in lepton-antilepton collisions, we have provided more control over the sparticles produced than for the hadron-hadron SUSY production processes described in section 4.7.

We remind the reader here that all SUSY particle data have to be read from an input file before event generation (see section 3.5.1).

4.4.1  IPROC=700-799: gaugino, slepton and/or squark production

With IPROC=700 one obtains the four processes IPROC=710,730,740,760 in the correct proportions. The matrix elements have been derived independently and the cross sections are in good agreement with those from SUSYGEN [62]. For all these processes the hard process scale EMSCA has been set to the centre-of-mass energy.

The gaugino and sfermion mixing conventions of Haber and Kane [40] are adopted in all cases. In addition, the neutralino mixing matrix ZMIXSS is defined internally in terms of the photino, zino and current eigenstate neutral higgsino components, instead of the bino, W3-ino and higgsino components adopted for ZMXNSS, equivalent to Nij in [40].

IPROC=710 gives lepton-antilepton®neutralino pair production. A number of additional IPROC codes have been provided to enable the user to produce given neutralino mass eigenstates. It should be noted that in order to provide a compact code for these processes there is more than one possible IPROC number for some processes. For example the final state c~10 c~30 can be produced using either the IPROC codes 713 or 721.

IPROC=730 gives lepton-antilepton®chargino pair production. As with the neutralino production there are codes to allow a given pair of charginos to be produced.

IPROC=740 gives lepton-antilepton®slepton pair production. In these processes for the first two generations the left/right eigenstates are produced, while for the third generation the mass eigenstates are produced. In the processes producing given slepton pairs for t~ production processes the left eigenstate is replaced by the lighter mass eigenstate and the right eigenstate by the heavier mass eigenstate.

IPROC=760 gives lepton-antilepton®squark pair production. As with the other processes additional codes are provided to allow the production of given q~q~* pairs. For stop and sbottom production the mass eigenstates are produced, while for the first two generations the left/right eigenstates are generated. As with the slepton production for the third generation when a given squark pair is requested the left eigenstate is replaced by the lighter mass eigenstate and the right eigenstate by the heavier mass eigenstate.

4.5  Other lepton-antilepton non-Standard-Model processes

4.5.1  IPROC=800-899: R-parity violating SUSY processes

A range of possible ¬  Rp production processes in lepton-antilepton collisions is included. Unlike the case of ¬  Rp production in hadron-hadron collisions, section 4.8.1, we have included processes for which there is either no s-channel resonance or the resonance is not kinematically accessible.

All the possible single sparticle production mechanisms which occur via the first term in the superpotential given in [42] are included. This includes the process l+l-®gn~ for which there is no s-channel resonance. As the cross section for this process diverges in the limit that the photon is collinear with the incoming lepton beams a cut on the pT of the outgoing particles pT>PTMIN has been imposed for this process (IPROC=850). The ISR is switched off for this process as including it would lead to a double counting of the photon radiation. For this reason this process is not included in the code IPROC=800 which generates all the other single sparticle production mechanisms.

We have also included some processes for the production of SM particles via s-channel sneutrino exchange. In addition to the s-channel sneutrino diagrams the t-channel sparticle exchanges, SM diagrams and all the interference terms are included. This uses a generalization of the formulae of [63]. A cut pT>PTMIN is used in the process l+l-®l+l- (IPROC=870) to avoid the divergence as t®0 in the Bhabha scattering cross section.

Except where stated explicitly above, no PTMIN cut is applied.

4.5.2  IPROC=900-999: MSSM Higgs production processes

For MSSM Higgs pair production (IPROC=955-965) a new subroutine has been introduced, HWHIHH, whereas for the other processes we make use of the implementation of their SM counterparts, which are based on the subroutines HWHIGW and HWHIGZ.

4.5.3  IPROC=1000-1199: SM and MSSM Higgs associated production

SM and MSSM Higgs production in association with fermion pairs in lepton-lepton collisions is available in version 6.5. These processes were introduced in [64] and their phenomenological relevance was discussed in [65].

Fermion masses are retained in the final state according to the HERWIG defaults. The same values appear in the Yukawa couplings. Notice that in the case of charged Higgs boson production the Cabibbo-Kobayashi-Maskawa mixing matrix has been assumed to be diagonal. Furthermore, due to the rather different phase space distribution of the final state products, all processes can only be produced separately, not collectively. Initial- and final-state radiation (both QED and QCD) and beamsstrahlung are included via the usual HERWIG algorithms. Finally notice that the use of the IPROC series 1000 and 1100 for l+l- processes required some internal modification to HERWIG, which was originally designed to generate leptonic processes only for IPROC< 1000. Now we assume an l+l- process whenever IPROC< 1300. These modifications have no implications for the traditional user, but may affect more knowledgeable ones who have edited previous versions of the main HERWIG code.

4.6  Hadron-hadron Standard Model processes

4.6.1  IPROC=1300-1499: Drell-Yan processes

The Drell-Yan code is extended to the production of all fermion pairs; 1300 gives all quark flavours; 1300+IQ a specific quark flavour, 1350 all leptons (including neutrinos) 1350+IL a specific lepton flavour. The s-channel component of the interference with like-flavour qq_ scattering is included here.

The initial-state parton showers in Drell-Yan processes are matched to the exact O (aS) matrix-element result as discussed in section 3.2.3. The routine HWBDYP implements the hard corrections whilst HWSBRN includes the soft corrections to the initial-state radiation. For further details see ref. [20].

4.6.2  IPROC=1500: QCD 2® 2 processes

IHPRO 1 + 2 ® 3 + 4 c/f conn.
1 q + q ® q + q 3 4 2 1
2 q + q ® q + q 4 3 1 2
3 q + q' ® q + q' 3 4 2 1
4 q + q_ ® q'+ q_' 2 4 1 3
5 q + q_ ® q + q_ 3 1 4 2
6 q + q_ ® q + q_ 2 4 1 3
7 q + q_ ® g + g 2 4 1 3
8 q + q_ ® g + g 2 3 4 1
9 q + q_' ® q + q_' 3 1 4 2
10 q + g ® q + g 3 1 4 2
11 q + g ® q + g 3 4 2 1
12 q_ + q ® q_' +q' 3 1 4 2
13 q_ + q ® q_ + q 2 4 1 3
14 q_ + q ® q_ + q 3 1 4 2
15 q_ + q ® g + g 3 1 4 2
16 q_ + q ® g + g 4 1 2 3
17 q_ + q' ® q_ + q' 2 4 1 3
18 q_ + q_ ® q_ + q_ 4 3 1 2
19 q_ + q_ ® q_ + q_ 3 4 2 1
20 q_ + q_' ® q_ + q_' 4 3 1 2
21 q_ + g ® q_ + g 2 4 1 3
22 q_ + g ® q_ + g 4 3 1 2
23 g + q ® g + q 2 4 1 3
24 g + q ® g + q 3 4 2 1
25 g + q_ ® g + q_ 3 1 4 2
26 g + q_ ® g + q_ 4 3 1 2
27 g + g ® q + q_ 2 4 1 3
28 g + g ® q + q_ 4 1 2 3
29 g + g ® g + g 4 1 2 3
30 g + g ® g + g 4 3 1 2
31 g + g ® g + g 2 4 1 3

Table 11: QCD subprocesses.

At present only 2®2 subprocesses are implemented. They are classified in table 11. Here and in other subprocess tables, `c/f conn.' refers to the colour/flavour connections between the partons: `i j k l' means that the colour of parton 1 comes from parton i, that of 2 from j, etc. For antiquarks, which have no colour (only anticolour), the label shows instead to which parton the flavour is connected. For this colour/flavour labelling all partons are defined as outgoing. Thus, for example, process 10 has colour connections 3 1 4 2, corresponding to the colour flow diagram:



When different colour flows are possible, they are listed as separate subprocesses. This separation is not exact but is normally a good approximation [6, 7]. The separation is now performed using the improved method proposed in [31], as outlined in section 3.1. The sum of the colour flows is the exact lowest-order cross section.

4.6.3  IPROC=1600-1699: Higgs boson production by parton fusion

IPROC = 1600+ID gives the sum of gg and qq_ fusion. The lowest-order formulae that we have used can be found in [41]. The hard processes are implemented in the subroutine HWHIGS.

4.6.4  IPROC=1700-1706: heavy quark production

The separation of colour flows is now performed using the improved method proposed in [31], as outlined in section 3.1. The classification of subprocesses according to IHPRO is as for IPROC=1500.

4.6.5  IPROC=1800: QCD direct photon plus jet production

The relevant IHPRO codes are 41-47 in table 12. For future reference we also collect here the codes for other processes that involve outgoing direct photons or incoming pointlike photons (IPROC=2200, 5000-5520). Note that the photon is colour/flavour-connected to itself. In the cases IHPRO=71-76, M represents an L=0 meson (see IPROC=5500).
  IHPRO   1 + 2    ®    3 + 4   c/f conn.  
  41   q + q_    ®    g + g   2 3 1 4  
  42   q + g    ®    q + g   3 1 2 4  
  43   q_ + q    ®    g + g   3 1 2 4  
  44   q_ + g    ®    q_ + g   2 3 1 4  
  45   g + q    ®    q + g   2 3 1 4  
  46   g + q_    ®    q_ + g   3 1 2 4  
  47   g + g    ®    g + g   2 3 1 4  
  51   g+ q    ®    g+ q   1 4 2 3  
  52   g+ q_    ®    g+ q_   1 3 4 2  
  53   g+ g    ®    q + q_   1 4 2 3  
  61   q + q_    ®    g+g   2 1 3 4  
  62   q_ + q    ®    g+g   2 1 3 4  
  63   g + g    ®    g+g   2 1 3 4  
  71   g+ q    ®    M(S=0) +q'   1 4 3 2  
  72   g+ q    ®    M(S=1)L+q'   1 4 3 2  
  73   g+ q    ®    M(S=1)T+q'   1 4 3 2  
  74   g+ q_    ®    M(S=0)+q_'   1 4 3 2  
  75   g+ q_    ®    M(S=1)L+q_'   1 4 3 2  
  76   g+ q_    ®    M(S=1)T+q_'   1 4 3 2  

Table 12: Direct photon subprocesses.

4.6.6  IPROC=1900-1999: Higgs boson production by weak boson fusion

The qq_® q(')q_(') VV® q(')q_ (')HSM0 subprocesses, for VV=W+W-,Z0Z0, summed over initial- and final-state quarks can be invoked by setting IPROC=1900+ID, with ID used to identify the Higgs decay.

The formulae used are well known in the literature and can be found e.g. in [22]. This process is administered by the subroutine HWHIGW, which also handles the similar cases initiated by e+e- and e±p collisions.

4.6.7  IPROC=2000-2008: single top production

The process of single top quark production by W-boson exchange includes so far only those processes initiated by a b-quark (or antiquark) and a first- or second-generation quark or antiquark. Note that this requires b-quarks to be available in the parton distribution functions.

4.6.8  IPROC=2100-2170: electro-
weak boson plus jet production

The electroweak boson decay correlations and width are now correctly included in these processes.

4.6.9  IPROC=2200: direct photon pair production

See section 4.6.5 for IHPRO codes (61-63).

4.6.10  IPROC=2300-2399: Higgs boson plus jet production

High transverse momentum, scalar Higgs production via one-loop diagrams, in association with a jet, is available as IPROC=2300, within the SM. Only the top quark is included in the loops with IAPHIG controlling the approximation used: IAPHIG=0 gives the zero top mass limit, 1 (default) the exact result, 2 the infinite top mass limit. The various subprocesses are illustrated in table 13.
IHPRO 1 + 2 ® 3 + 4 c/f conn.
81 q + q_ ® g +HSM0 2 3 1 4
82 q + g ® q +HSM0 3 1 2 4
83 q_ + q ® g +HSM0 3 1 2 4
84 q_ + g ® q_ +HSM0 2 3 1 4
85 g + q ® q +HSM0 2 3 1 4
86 g + q_ ® q_ +HSM0 3 1 2 4
87 g + g ® g +HSM0 2 3 1 4

Table 13: Higgs plus jet subprocesses.

Note that the Higgs boson is colour/flavour connected to itself.

The relevant routines HWHGJ1, HWHGJA/B/C/D, HWUCI2 and HWULI2 use (non-standard Fortran 77) DOUBLE COMPLEX variables which may not be accepted by some compilers and are called COMPLEX*16 by others. Users can change to COMPLEX variables, but this involves a risk of rounding errors spoiling numerical cancellations.

4.6.11  IPROC=2400-2450: colour singlet exchange

IPROC=2400: Two-to-two parton scattering via exchange of a colour singlet, Mueller-Tang pomeron [66]. The fixed aS and w0 are given by ASFIXD (default 0.25) and OMEGA0 (0.3) respectively.
IPROC=2450: Photon exchange, for like-flavour qq_ pairs including the t-channel component of the interference with qq_® qq_ via an s-channel photon or Z0.

4.6.12  IPROC=2500-2599: Higgs boson plus top quark pair production

The SM 2®3 Higgs production subprocesses of the type gg® tt_HSM0 and qq_® tt_HSM0, for any flavour of initial state quarks q, are new to HERWIG version 6 and are handled by the subroutines HWHIGQ and HWH2QH. They are invoked by setting IPROC=2500+ID (both gg and qq_), with ID administering the Higgs decay channels as in IPROC=300+ID. The initial state quark flavours q are always summed over. Notice that, given the size of the Yukawa couplings of the Higgs boson to quarks, in practice only the associated production with top quarks is of phenomenological relevance in the SM. The matrix elements used in the implementation can be found in ref. [22]. The treatment of the Higgs width here is as described in section 3.4.

4.6.13  IPROC=2600-2799: Higgs plus weak boson production

The associated production of SM Higgs scalars with W± (IPROC=2600-2699) and Z0 (IPROC=2700-2799) gauge bosons initiated by quark-antiquark fusion via the 2®2 processes qq_® W±HSM0 and qq_® Z0HSM0 is now available. A summation is as usual intended on the incoming quarks. The formulae given in [22] are used. The production processes are generated by the two new subroutines HWHIGV and HWH2VH whereas the Higgs decays are administered through the ID increment, as in IPROC=300+ID. Again, the treatment of the Higgs boson width here is as discussed in section 3.4.

4.6.14  IPROC=2800-2825: Gauge boson pair production

In version 6.3, the code already included in HERWIG for e+e-® WW/ZZ [67] was adapted for hadron-hadron collisions, including photons and the photon/Z interference for the resonant diagrams.

All of these processes use a cut EMMIN (default value 20 GeV) on the mass of the gauge bosons produced. The cut PTMIN (default 10 GeV) on the transverse momentum of the bosons is also used. Both these cuts should not be taken to zero simultaneously if photon terms are included. The phase space for these processes contains a number of peaks and it was therefore necessary to use an adaptive multi-channel phase-space integration method which is described below.

In versions 6.4 and above, the hard process scale EMSCA for this process has been changed to the parton-parton centre-of-mass energy from the average of the produced boson masses, which was used previously.

4.6.15  IPROC=2850-2880: Gauge boson pairs using MC@NLO

These codes activate the interface to the program MC@NLO for diboson production at next-to-leading order (see sect. 9.6).

4.6.16  IPROC=2900-2916: QQ_ Z production

The matrix elements of [68] were used for the massless case and an independent calculation, using the approach of [69], which was checked both for gauge invariance and against the massless case for the massive result.

In both cases the decay of the Z is fully included and is selected using MODBOS. PTMIN controls the minimum transverse momentum of the outgoing quarks. As with gauge boson pair production it was necessary to use an optimized multi-channel phase-space integrator which is described below.

The phase space for both gauge boson pair production and QQ_ Z is complicated, as these processes are both treated as 2®4 processes. In order to obtain a reasonable efficiency it was necessary to adopt a multi-channel approach based on that described in [70, 71].

In each case a number of different channels are included which attempt to map the phase space for the different processes. The default weights for these different channels have been chosen to optimize the efficiency for the Tevatron and LHC; a choice of which to use is made based on the beam energy. However, the choice is affected by the phase space cuts applied. Therefore if these are significantly altered the weights for the different channels need re-optimizing.

This is controlled by the new variable OPTM (default .FALSE.). If OPTM=.TRUE., before performing the initial search for the maximum weight HERWIG will attempt to optimize the efficiency using the procedure suggested in [71].

This is done by generating IOPSTP (default 10) iterations of IOPSH (default 1000) events. The choice of IOPSTP and IOPSH is a compromise between run time and accuracy. The value of IOPSH should not be significantly reduced because the procedure attempts to minimize the error on the Monte Carlo evaluation of the cross section and if IOPSH is small the error on the error can be significant. If you need to re-optimize the weights we would recommend a long run just to optimize the weights which can then be used in all the runs to generate events. The new subroutine HWIPHS was added to initialise the phase space.

4.7  Hadron-hadron supersymmetric processes (MSSM)

The R-parity conserving SUSY processes occupy the IPROC entries of the series 3000 (sparticle processes) and 3300-3600 (Higgs boson production), while ¬  Rp processes have IPROC=4000-4199.

As with the lepton-antilepton SUSY processes the SUSY particle data must be read in from an input file before event generation as described in section 3.5.1. In particular, unlike those of the SM Higgs boson HSM0, the widths and decay modes of the SUSY Higgs bosons are not computed by HERWIG.

A large number of final states involving the production of both neutral and charged Higgs bosons of the Minimal Supersymmetric Standard Model (MSSM) have been made available in version 6.3. They all proceed via 2®3 body hard scattering subprocesses. They are listed below, with corresponding process numbers (IQ and ISQ are as detailed in the following subsections). Further details of their implementation can be found in [22]

The new subroutines introduced to administer the following processes are HWHIBQ and HWH2BH for IPROC=3500, plus HWHISQ and HWH2SH for IPROC=3100, 3200. In addition, HWHIGQ has been modified to accommodate the IPROC=3800 series.

4.7.1  IPROC=3000-3030: sparton, gaugino and/or slepton production

With IPROC=3000 one obtains the three following processes, IPROC=3010, 3020, 3030, in the correct proportions. The variable IHPRO gives the subprocess actually generated (table 14).

The matrix elements have been derived independently and the cross sections are in good agreement with those from ISAJET [39].

The hard process scale EMSCA has to be chosen globally for all sparton processes, e.g. for the argument of the QCD coupling. This is done using the kinematics appropriate to production of the lightest supersymmetric particle (LSP) and
EMSCA= æ
ç
ç
ç
ç
ç
è
2
^
s
 
^
t
 
 '
^
u
 
 '
^
s
 
 2+
^
t
 
 '2+
^
u
 
 '2
 ö
÷
÷
÷
÷
÷
ø
1/2






 
 .     (3)
with t^'=t^-m2, u^'=u^-m2 where m is the LSP mass.

IPROC=3010 gives 2-parton ® 2-sparton processes. All QCD sparton, i.e. squark and gluino, pair production processes are implemented. The matrix elements and the scheme for separating different colour flow parts are as given in [31].

IPROC=3020 gives 2-parton ® 2-gaugino or gaugino+sparton processes. All gaugino, i.e. chargino and neutralino, pair production processes and gaugino-sparton associated production processes are implemented.

The gaugino and sfermion mixing conventions of Haber and Kane  [40] are used as described in section 4.4.

The various subprocesses, which include the 1«2 and charge conjugate reactions omitted below for brevity, are shown in table 14.
IHPRO 1 + 2 ® 3 + 4 c/f conn.
3021 q+q_ ® c~a±+c~b± 2 1 3 4
3022 q+q_ ® c~i0 +c~j0 2 1 3 4
3023 q+q_' ® c~a±+c~i0 2 1 3 4
3024 q+q_ ® c~i0 +g~ 2 4 3 1
3025 q+q_' ® c~a±+g~ 2 4 3 1
3026 g+ q ® c~i0 +q~ 2 4 3 1
3027 g+ q ® c~a±+q~' 2 4 3 1

Table 14: SUSY subprocesses.

Note that the gauginos connect to themselves. The indices a,b,i,j label gauginos in the order of increasing mass and take values a,b=1-2 and i,j=1-4.

Gaugino mixing matrices are implemented in all subprocesses. The associated production subprocesses IHPRO=3026, 3027 include stop and sbottom left-right mixings. CKM mixing is implemented in subprocesses IHPRO=3023, 3025, 3027 but neglected in subprocess IHPRO=3021.

IPROC=3030 gives 2-parton ® 2-slepton processes. All Drell-Yan slepton production processes are implemented. The formulae agree with those of refs. [72, 73] in the limit of no stau left-right mixing.

4.7.2  IPROC=3110-3178: charged Higgs plus squark pair production

The production of charged Higgs bosons of the MSSM in association with squark pairs, of bottom and top flavours only, is implemented via the 3100 series of IPROC numbers (table 15). Their phenomenological relevance has been discussed in [74].
IPROC partons ® spartons Higgs
3110 gg+qq_ ® q~iq~j'* H±
3111 gg+qq_ ® b~1t~1* H+
3112 gg+qq_ ® b~1t~2* H+
3113 gg+qq_ ® b~2t~1* H+
3114 gg+qq_ ® b~2t~2* H+
3115 gg+qq_ ® t~1b~1* H-
3116 gg+qq_ ® t~1b~2* H-
3117 gg+qq_ ® t~2b~1* H-
3118 gg+qq_ ® t~2b~2* H-

Table 15: Processes for IPROC=3110-3178.

One should add 30(60) to IPROC for gg(qq_)-only initiated processes.

4.7.3  IPROC=3210-3298: neutral
Higgs plus squark pair production

The production of neutral Higgs bosons of the MSSM in association with squark pairs, of bottom and top flavours only, is implemented via the 3200 series of IPROC numbers (table 16). Their phenomenological relevance has been discussed in [74, 75].

One should add 30(60) to IPROC for gg(qq_)-only initiated processes.

4.7.4  IPROC=3310-3375: Higgs-Higgs and Higgs-gauge boson pair
production

The production of Higgs-Higgs and Higgs-gauge boson pairs of the MSSM is implemented at tree level. We include gauge boson mediated contributions but not Higgstrahlung from the initial state, the only exception being W±H± production (IPROC=3350), which does include diagrams where the Higgs boson couples to the initial partons as well as those mediated by neutral Higgs states [76]. Some of these processes are the MSSM equivalent of IPROC=2600 and 2700 described earlier for the case of the SM. Here, the cases IPROC=3310(3320) correspond to W±h0(W±H0) and IPROC=3360(3370) to Z0 h0(Z0 H0) final states. (No similar A0 production can occur at leading order.) The array ENHANC is used to implement the MSSM couplings of Higgs scalars to gauge bosons.
IPROC partons ® spartons Higgs
3210(3220)[3230] gg+qq_ ® q~iq~j* h(H)[A]
3211(3221)[3231] gg+qq_ ® b~1b~1* h(H)[A]
3212(3222)[3232] gg+qq_ ® b~1b~2* h(H)[A]
3213(3223)[3233] gg+qq_ ® b~2b~1* h(H)[A]
3214(3224)[3234] gg+qq_ ® b~2b~2* h(H)[A]
3215(3225)[3235] gg+qq_ ® t~1t~1* h(H)[A]
3216(3226)[3236] gg+qq_ ® t~1t~2* h(H)[A]
3217(3227)[3237] gg+qq_ ® t~2t~1* h(H)[A]
3218(3228)[3238] gg+qq_ ® t~2t~2* h(H)[A]

Table 16: Processes for IPROC=3210-3298.

4.7.5  IPROC=3410-3450: Higgs boson plus heavy quark production

We have included so far only those processes initiated by a b-quark (or antiquark) and a gluon. Note that this requires b-quarks to be available in the parton distribution functions.

4.7.6  IPROC=3500: charged Higgs boson from bq-initiated processes

This process is relevant for charged Higgs scalar production at large tanb values, see [77].

4.7.7  IPROC=3610-3630: neutral Higgs production by parton fusion

These processes are the MSSM analogues of the SM processes 1600 etc. Recall however that the MSSM Higgs decay modes are controlled by the SUSY input data (see section 3.5.1) and not by the value of IPROC. The subroutines HWHIGS and HWHIGT have been modified to take account of squark loop contributions and parity-violating Higgs-fermion couplings in the MSSM case.

4.7.8  IPROC=3710-3720: neutral Higgs production via weak boson fusion

These processes are the MSSM counterparts of the SM process of weak vector-vector fusion in hadronic collisions (IPROC=1900). They are computed using the same subroutines and setting the F0 VV couplings to sin(b-a) for F0=h0 (IPROC=3710) and to cos(b-a) for F0=H0 (IPROC=3720), respectively, where V=W±,Z0. (There is no A0VV coupling at tree level.)

These reactions are two of the major direct production channels of neutral CP-even Higgs bosons of the MSSM at hadron colliders, such as the LHC (see e.g. [78]).

4.7.9  IPROC=3811-3899: Higgs boson plus heavy quark pair production

For the case of neutral Higgs states, IPROC=3810+IQ corresponds to h0 production, IPROC=3820+IQ to H0 and IPROC=3830+IQ to A0. (For the last case, the variable PARITY is set to -1.) Note also the production of charged Higgs states, via IPROC=3839, 3869 and 3899, in association with pairs of top-bottom quarks.

In the usage of the IPROC numbers corresponding to neutral Higgs states, when b-quarks are involved in gg-fusion modes (IPROC(+30)=3845, 3855 or 3865), the user should take care to avoid double-counting the chosen process with the corresponding 2®1 and 2® 2 cases of IPROC=3610-3630 and IPROC=3410-3430 initiated by quark-antiquark annihilation, i.e. b b_® Higgs, and (anti)quark-gluon scattering, i.e. bg® b Higgs, respectively: see [79]. Similar arguments hold for the charged Higgs states, as the gg-induced process (IPROC(+30)=3869) is an alternative implementation of IPROC=3450 [80].

The associate production of neutral Higgs bosons (both CP-even and CP-odd) of the MSSM with heavy QQ_ pairs (Q=b and t) is of extreme phenomenological relevance as a discovery mode of Higgs scalars, both at the Tevatron (Run 2) and the LHC (see e.g. [78, 81]), as is the case of the charged Higgs channel [80, 82].

4.8  Other hadron-hadron non-Standard-Model processes

4.8.1  IPROC=4000-4199: R-parity violating SUSY processes

We include all the possible production processes of resonant sleptons and squarks in hadron-hadron collisions, for arbitrary numbers of non-zero ¬  Rp couplings. These processes are implemented as two-to-two processes, i.e. with the decay of the resonant particle included. This allows us to include the t-channel diagrams where these occur. However we have not implemented those processes which can only occur via a t-channel diagram, or where the resonance will never be accessible. So for example while we include the process ui dj ® b~1* Z0, which can occur via a resonant b~2*, we do not include the process ui dj ® b~2* Z0, which cannot occur via a resonant diagram. In all cases both the processes listed and their charge conjugates are included.

The scale choice is s^1/2 rather than the conventional transverse mass, due to the large number of different processes which must be calculated. The colour connection structure of these processes and their matrix elements can be found in ref. [42].

4.8.2  IPROC=4200-4299 : graviton resonance production

In some models with extra dimensions, Kaluza-Klein excitations of the graviton can be produced with significant cross sections at TeV-scale colliders. When the scale of the extra dimensions is not large, the excitations are manifest as discrete resonances.

The production of a resonant excitation of the graviton is implemented as a 2®2 process including the decay of the resonance. The process is treated in a model-independent way, assuming only that there is a universal coupling of the graviton resonance to the SM fields. The effective lagrangian is given by
LI = -
1
Lp
 hµnTµn ,
where hµn is the spin-2 field and Tµn is the energy-momentum tensor of the SM fields. Although in these models there are usually many resonances, we have implemented only one. Others can be studied by making appropriate changes in the parameters. Graviton resonance production is described in more detail in [83].

The process is controlled by the coupling GRVLAM=Lp, with dimensions of mass and default value 10 TeV, and by the mass EMGRV (default 1 TeV) and width GAMGRV of the resonance. If the width is set to zero (the default), the subroutine HWHGRV which calculates the cross section also calculates the width.

The parton-level cross section for this process is non-unitary and is proportional to s^/EMGRV4 at high energies. The fall-off of the parton distribution functions is not sufficient to suppress this bad high energy behaviour. Hence the parameters EMMIN and EMMAX controlling the minimum and maximum values of s^1/2 must be set. The default is to set these to 90% and 110% of the graviton resonance mass respectively. If the width of the resonance is more than a few percent of its mass then these limits should be reset.

After event generation, IHPRO is set to 50 for qq_ initiated subprocesses and to 51 for gg initiated subprocesses.

4.9  Photon-hadron and photon-photon processes

4.9.1  IPROC=5000-5520: pointlike photon-hadron processes

Pointlike photon-hadron scattering to produce QCD jets is available as IPROC=5000-5206. This is suitable for fixed-target photoproduction, provided events are generated in a frame in which the target has high momentum, and then boosted back to the lab. This is done if USECMF=.TRUE., the default, in which case the frame for event generation is the beam-target c.m. frame. IPROC=5100+IQ gives flavour IQ pair production, g+g® QQ_, and IPROC=5200+IQ gives flavour IQ single excitation, g+Q® g+Q, both including quark masses. IPROC=5000 gives a sum over all processes and flavours, 5100 and 5200, with massless quark kinematics. In all cases, after event generation the code IHPRO is set to 51, 52 or 53 according to the hard subprocess, as specified in section 4.6.5. IPROC=5300 gives Compton scattering, g + q®g + q.

The direct, higher twist, production of light (u,d,s) L=0 mesons by point-like photons is also available: IPROC = 5500 for all spin =0 and 1 mesons; 5510 for only S=0 mesons; and 5520 for only S=1 mesons. The vector mesons are produced with transverse or longitudinal polarization and decayed accordingly. The corresponding IHPRO codes (71-76) are also listed in section 4.6.5.

All these processes are available with lepton as well as hadron beams, using the Equivalent Photon Approximation. The phase-space variable WHMIN sets the minimum allowed hadronic mass and affects photoproduction reactions (g-hadron and g-g) and DIS. In lepton-hadron DIS it is largely irrelevant since there is already a cut on Bjorken y which at fixed s is almost the same, but for lepton-gamma DIS it makes a big difference. Direct g+g*® q+q_ is included in the hard correction for lepton-gamma DIS.

4.9.2  IPROC=6000-6010: pointlike photon-photon processes

Direct gg® charged particle pairs has been implemented with IPROC=6000+IQ: if IQ=1-6 then only quark flavour IQ is produced, if IQ=7,8 or 9 then only lepton flavour e, µ or t is produced and if IQ=10 then only W pairs are produced: in these cases particle mass effects are included. If IQ=0, the natural mix of quark pairs is produced using massless matrix elements but including a mass threshold cut. The range of allowed transverse momenta is controlled by PTMIN and PTMAX as usual.

4.10  Baryon number violating processes

4.10.1  IPROC=7000-7999: generated by the HERBVI package

See section 9.8 for details.

4.11  Minimum bias soft hadron-hadron collisions

4.11.1  IPROC=8000: minimum bias soft hadron-hadron event

Non-diffractive, soft hadronic, minimum bias events (IPROC=8000) can be generated for the following combinations of beam and target: p,p_,p±,K±,e±,µ±,g on target p (or vice versa); p,p_ on target n (or vice versa); or g on g. The event weight is the estimated cross section based on the parameterizations of Donnachie and Landshoff [84]. The non-diffractive cross section is assumed to be 70% of the total. For lepton beams a photon is first generated using the Effective Photon Approximation (see section 4.1) and then the on-shell photon cross section is used. See section 3.7.3 for discussion of the model used and the relevant parameters.

4.12  Deep inelastic lepton scattering

Deep inelastic (DIS) processes are broadly divided into those that start at O(aS0) (IPROC=90**) and those like heavy quark and dijet production which start at O(aS1) (IPROC=91**). Note that the DIS O(aS) jet production processes, IPROC=92**, have been withdrawn since they are subsumed within IPROC=91**.

The default limits on Q2 in DIS processes (Q2MIN, Q2MAX) have been set very small/large (0, 1010 GeV) and are reset to the kinematic limits unless changed by the user. This means the default Q2MIN is not suitable for simple neutral current DIS (IPROC=9000 etc), but is appropriate for jet and heavy quark photoproduction. The range of the Bjorken-y variable (y=Q2/xs) can be limited by YBMIN and YBMAX.

The kinematic reconstruction of DIS processes can now take place in the Breit frame, if BREIT=.TRUE. (the default value). Previous versions used the lab frame. Although the reconstruction is fully invariant under Lorentz boosts along the incoming hadron's direction, it is not under transverse boosts, so there should be some difference between the two frames. The boost is not performed for very small Q2 (<10-4) to avoid numerical instabilities, but the two frames are in any case equivalent for such small values.

The phase space for boson-gluon fusion is controlled by the parameters EMMIN, EMMAX (see section 4.3.4). The default values (0, s1/2) correspond to the behaviour of version 5.1.

Lepton beam polarization effects are included in all DIS processes apart from J/y production. The polarization is specified as in lepton-antilepton processes, i.e. by setting the 3-vectors EPOLN and PPOLN: component 3 is longitudinal and 1,2 transverse. Transverse only occurs in e+e- routines; recall that two transverse `measurements' are needed to see an effect so it should not arise elsewhere. Note that in DIS processes one has to set either EPOLN if it is a lepton or (exclusive) PPOLN if an antilepton.

All the DIS processes IPROC=9000-9599 are available in e+e- as well as lepton-hadron collisions. The program generates a photon from the second beam (only) in the Equivalent Photon Approximation and the default is to use Drees-Grassie [53, 54] structure functions for DIS on the photon.

The parameter WHMIN sets the minimum allowed hadronic mass in DIS. In lepton-hadron DIS it is largely irrelevant since there is already a cut on Bjorken y which at fixed s is almost the same but for lepton-gamma DIS it makes a big difference.

In addition to the processes listed here, a full simulation of QCD instanton-induced processes in DIS [85] is available through an interface to the program QCDINS [86]. For details, see the web page http://www.desy.de/ t00fri/qcdins/qcdins.html

4.12.1  IPROC=9000-9006: neutral current

Matrix-element corrections to DIS are available, following the general method described in section 3.2.3. The hard correction is implemented in HWBDIS and the soft correction is included in routines HWBRAN and HWSBRN for the final- and initial-state radiation respectively.

4.12.2  IPROC=9010-9016: charged current

These are the charged current processes corresponding to those above, with the same treatment of hard and soft matrix element corrections.

4.12.3  IPROC=9100-9130: O(aS) neutral current processes

The photoproduction processes have been extended from the original heavy quark production program, to include all quark pair production (IPROC=9100-9106) and QCD Compton (IPROC=9110-9122), as well as the sum of the two (IPROC=9130). The possible flavours for processes 9100, 9110 and 9130 are limited by the input parameters IFLMIN and IFLMAX (defaults are 1 and 3, i.e. only u,d,s flavours).

A sign error has been corrected that led to the incorrect sign for the lepton-jet azimuthal correlation in QCD Compton processes in versions before 5.7.

J/y production (IPROC=9107) now uses the Equivalent Photon Approximation instead of Weizsacker-Williams, with the same phase-space cuts as hadronic processes with lepton beams (see section 4.1).

The argument of the running coupling is controlled by the parameter BGSHAT - see section 6.

4.12.4  IPROC=9140-9144: charged-current heavy quark production

At present only W±g fusion processes are implemented.

4.12.5  IPROC=9500-9599: Higgs production by weak boson fusion

The process of W+ W-/Z0Z0 fusion into the SM Higgs boson is also available in e±p collisions, as IPROC = 9500+ID. For details of the implementation, see section 3.4 and the corresponding processes initiated by lepton-lepton and hadron-hadron collisions, IPROC=400+ID and IPROC=1900+ID, respectively.

4.13  Including new subprocesses

The procedure for including further subprocesses remains substantially as described in ref. [1] but is repeated here for completeness. However, we now recommend that users make use of the Les Houches interface, see Section 9.1.

The parton and hard subprocess 4-momenta, masses and identity codes need to be entered in COMMON/HEPEVT/ with the appropriate status codes ISTHEP(I)=110-114 to tell the program which is which (see the table in section 8.3.1). The colour/flavour structure should be specified by the second mother and daughter pointers as explained in section 4.6.2.

The HERWIG identity codes IDHW(I) in COMMON/HWEVNT/ also need to be set correctly. The IDHW codes can be listed in a run with IPRINT=2: the most important are the quarks 1-6 (as IDHEP), antiquarks 7-12, gluon 13, overall centre-of-mass 14, hard centre-of-mass 15, soft centre-of-mass 16, photon 59, leptons 121-126, antileptons 127-132.

The utility subroutine HWUIDT(IOPT,IPDG,IHWG,NAME) is provided to translate between Particle Data Group code IPDG, HERWIG code IHWG, and HERWIG CHARACTER*8 NAME, with IOPT=1,2,3 depending on which of IPDG, IHWG and NAME is the input argument.

Consider for example the process of virtual photon-gluon fusion to make b+b_ in proton-electron collisions (in fact this process is included as IPROC=9105). We assume the user provides a subroutine to generate the momenta PHEP for the hard subprocess e+g® e b b_. The colour structure is



Thus the momenta generated, together with those of the initial beams and the overall centre of mass, could be entered in the sequence shown in table 17.
IHEP Entry ISTHEP IDHEP JMOHEP JDAHEP IDHW
1 e beam 101 11 0 0 0 0 121
2 p beam 102 2212 0 0 0 0 73
3 ep c.m. 103 0 0 0 0 0 14
4 e in 111 11 6 7 0 7 121
5 gluon 112 21 6 9 0 8 13
6 hard cm 110 0 4 5 7 9 15
7 e out 113 11 6 4 0 4 121
8 b 114 5 6 5 0 9 5
9 b_ 114 -5 6 8 0 5 11

Table 17: Event record entries for eg® ebb_.

Note that if there are more than two outgoing partons, the first has status 113 and all the others 114. Each parton has JMOHEP(1,I)=6 to indicate the location of the hard c.m. for this subprocess, while JMOHEP(2,I) gives the location of the colour mother (treating the incoming gluon as outgoing) or the connected electron. JDAHEP(1,I) will be set by the jet generator HWBGEN, while JDAHEP(2,I) points to the anticolour mother (or connected electron). Finally the HERWIG identifiers IDHW(I) could be set to the indicated values by means of the translation subroutine HWUIDT as follows: 
                CHARACTER*8 NAME
                .....
                NHEP=9
                IDHEP(1)=11
                IDHEP(2)=2212
                .....
                IDHEP(9)=-5
                DO 10 I=1,NHEP
             10 CALL HWUIDT(1,IDHEP(I),IDHW(I),NAME)
                IDHW(6)=15
The last statement is needed because IDPDG(I)=0 returns IDHW(I)=14. If subroutine HWBGEN is now called, it will find the coloured partons and generate QCD jets from them. Subsequent calls to HWCFOR etc. can then be used to form clusters and hadronize them.

If the hard subprocess routine is called from HWEPRO, like those already provided, it must have two options controlled by the logical variable GENEV in COMMON/HWHARD/. For GENEV=.FALSE., an event weight (normally the cross section in nanobarns) is generated and stored as EVWGT in COMMON/HWEVNT/. If this weight is accepted by HWEPRO, the subroutine is called a second time with GENEV=.TRUE. and the corresponding event data should then be generated and stored as explained above. On certain computers it will be necessary to SAVE those variables that determine event characteristics between the two subroutine calls.

The parameter NMXJET sets the maximum number of outgoing partons in a hard subprocess (default 200).

5  Parameters

The quantities that may be regarded as adjustable parameters are indicated in table 18. Notes on parameters are given below.
Name Description Default
QCDLAM LQCD (see below) 0.18
RMASS(1) Down quark mass 0.32
RMASS(2) Up quark mass 0.32
RMASS(3) Strange quark mass 0.50
RMASS(4) Charmed quark mass 1.55
RMASS(5) Bottom quark mass 4.95
RMASS(6) Top quark mass 174.3
RMASS(13) Gluon effective mass 0.75
VQCUT Quark virtuality cutoff (added to 0.48
  quark masses in parton showers)  
VGCUT Gluon virtuality cutoff (added to 0.10
  effective mass in parton showers)  
VPCUT Photon virtuality cutoff 0.40
CLMAX Maximum cluster mass parameter 3.35
CLPOW Power in maximum cluster mass 2.00
PSPLT(1) Split cluster spectrum parameter 1.00
PSPLT(2) 1: light cluster, 2 heavy b-cluster PSPLT(1)
QDIQK Maximum scale for gluon®diquarks 0.00
PDIQK Gluon®diquarks rate parameter 5.00
QSPAC Cutoff for spacelike evolution 2.50
PTRMS Intrinsic pT in incoming hadrons 0.00

Table 18: Adjustable parameters.

In practice, the parameters that have been found most effective in fitting data are QCDLAM, the gluon effective mass RMASS(13), and the cluster mass parameter CLMAX. Note that QSPAC, PTRMS and ENSOF do not affect lepton-lepton collisions.

The default parameter values are based on those that were found to give good agreement when comparing earlier versions with event shape distributions at LEP. However, the substantial changes in this version mean that a re-tuning of parameters would be very worthwhile.

Up-to-date details of HERWIG parameter tunes can be found via the official web page cited in section 2.

6  Control switches, constants and options

A number of quantities can be reset to control the program and various options:
Name Description Default
NEVHEP Current number of events 0
NHEP Current number of entries in /HEPEVT/ 0
IPRINT Information to include in print out 1
MAXPR Number of events to print out 1
PRVTX Include vertex information in print out .TRUE.
NPRFMT Controls number of sig. figs. in print out 1
PRNDEC Use decimal/hexadecimal in print out .TRUE.
PRNDEF Produce ASCII (stout) version of print out .TRUE.
PRNTEX Produce LATEX version of print out .FALSE.
PRNWEB Produce html version of print out .FALSE.
MAXER Maximum number of errors to tolerate 10
LWEVT Unit for writing output events 0
LRSUD Unit for reading Sudakov table 0
LWSUD Unit for writing Sudakov table 77
SUDORD aS order in Sudakov table 1
INTER Order of interpolation in Sudakov tables 3
NRN(1) Random number seed 1 17673
NRN(2) Random number seed 2 63565
WGTMAX Max. weight (0 to search for it) 0.0
NOWGT Generate unweighted events with EVWGT=AVWGT .TRUE.
AVWGT Mean event weight 1.0
EFFMIN Min. acceptable Monte Carlo efficiency 0.001
NEGWTS Whether or not to allow negative weight events .FALSE.
AZSOFT Include soft gluon azimuthal correlations .TRUE.
AZSPIN Include gluon spin azimuthal correlations .TRUE.
HARDME Use hard matrix-element corrections .TRUE.
SOFTME Use soft matrix-element corrections .TRUE.
GCUTME Gluon energy cut in top M.E. correction 2.0
NCOLO Number of colours 3
NFLAV Number of (producible) flavours 6
MODPDF(I) PDFLIB parton set and author group for beam -1
AUTPDF(I) I (=1,2) (if MODPDF<0 do not use PDFLIB) 'MRS'
NSTRU Input parton set (1,2 = Duke-Owens sets 1,2;  
  3,4 = EHLQ sets 1,2; 5 = Owens set 1.1, 8
  6,7,8 = MRST, see table 8)  
PRSOF Probability of soft underlying event 1.0
ENSOF Multiplicity enhancement for SUE:  
n=á npp_ñ(ENSOFs1/2) 1.0
PMBN1 Mean multiplicity in SUE/Min. bias event +9.110
PMBN2 á npp_ñ(s1/2)= PMBN1sPMBN2+PMBN3 +0.115
PMBN3   -9.500
PMBK1 Negative binomial param.  
  k-1=PMBK1 loge(s)+PMBK2 +0.029
PMBK2   -0.014
PMBM1 Soft cluster mass spectrum:  
(M-M1-M2 -PMBM1)e-PMBM2M 0.2
PMBM2   2.0
PMBP1 Soft cluster PT spectrum: pTe-PMBPi
( pT2+M2 )
1/2
 
 
,		  
  d,u quarks 5.2
PMBP2 s,c quarks 3.0
PMBP3 diquarks 5.2
IOPREM Options for treatment of remnant clusters 1
BTCLM Mass parameter in remnant fragmentation 1.0
VMIN2 Min. parton virtuality2 in distance calcs. 0.1
CLRECO Include colour rearrangement .FALSE.
PRECO Probability for rearrangement 1/9
EXAG Lifetime scaling for weak bosons 1.0
ETAMIX h/h' mixing angle in degrees -23
PHIMIX f/w mixing angle in degrees +36
H1MIX h1(1380)/h1(1170) mixing angle in degrees tan-1(1/21/2)
F0MIX -/f0(1370) mixing angle in degrees tan-1(1/21/2)
F1MIX f1(1420)/f1(1285) mixing angle in degrees tan-1(1/21/2)
F2MIX f'2/f2 mixing angle in degrees +26
ET2MIX h2(1645)/h2(1870) mixing angle in degrees tan-1(1/21/2)
OMHMIX -/w(1600) mixing angle in degrees tan-1(1/21/2)
PH3MIX f3/w3 mixing angle in degrees +28
B1LIM B cluster ® 1 hadron parameter 0.0
CLDIR(I) Decay orientation of perturbative clusters, 1,1
  0: isotropic, 1: along quark direction  
CLSMR(I) Width of gaussian angle smearing, 0.0,0.0
  (I=1: light cluster, I=2: heavy b-cluster)  
PWT(I) A priori weights for ff_-pairs in cluster decay, 1.0
  I=1-6: f=d,u,s,c,b,t I=7: f=qq'  
REPWT(L,J,N) A priori weight for n(2S+1)LJ mesons 1.0
SNGWT A priori weight for singlet baryons 1.0
DECWT A priori weight for decuplet baryons 1.0
PLTCUT Minimum lifetime for particle to be set stable 1.0×10-8
VTOCDK(I) Veto decay of clusters to hadron I .FALSE.
VTOCDK(I) Veto decay of resonances to hadron I .FALSE.
  I=290-293, f0(980),a0(980) .TRUE.
PIPSMR Smear the primary vertex .FALSE.
VIPWID(1) x width (mm) 0.25
VIPWID(2) y width (mm) 0.015
VIPWID(3) z width (mm) 1.8
MAXDKL Veto decays outside given volume .FALSE.
IOPDKL Option for volume: 1=cylinder, 2=sphere 1
DXRCYL Radius for cylindrical option (mm) 20
DXZMAX Length for cylindrical option (mm) 500
DXRSPH Radius for spherical option (mm) 100
BDECAY Controls which B Decay package is used. 'HERW'
  Allowed values are: 'HERW', 'EURO' or 'CLEO'  
MIXING Include neutral B meson mixing .TRUE.
XMIX(1) D M/G for Bs0 10.0
XMIX(2) D M/G for Bd0 0.7
YMIX(1) DG/2G for Bs0 0.2
YMIX(2) DG/2G for Bd0 0.0
RMASS(198) W+ mass 80.42
RMASS(199) W- mass RMASS(198)
GAMW W± width 2.12
RMASS(200) Z0 mass 91.188
GAMZ Z0 width 2.495
WZRFR Use W/Z rest frame for decay parton showers .TRUE.
MODBOS(I) Force decay modes for weak bosons,  
  see sect. 3.4 0
RMASS(201) SM Higgs mass 115
IOPHIG Options for large Higgs mass distribution 3
GAMMAX Limit on range of Higgs mass distribution 10.
ENHANC(I) Enhancement factor for SM Higgs decay mode I 1.0
RMASS(209) Hypothetical 4th generation `bottom' quark mass 200.
RMASS(215) corresponding antiquark mass RMASS(209)
ALPHEM Thompson limit value of aem(0) 0.0072993
SWEIN Value of sin2qW 0.2319
QFCH(I) Fermion electric charge I=1-6: d,..,t  
AFCH(I,J) Fermion weak axial charge I=10-16: e,..,nt see sect. 4.2.2
VFCH(I,J) Fermion weak vector charge J=1: Z, J=2: Z'  
ZPRIME Include a Z' in g*/Z0 processes .FALSE
RMASS(202) Mass of the Z' 500.
GAMZP Width of the Z' 5.0
VCKM(I,J) Cabibbo-Kobayashi-Maskawa matrix elements:
æ
ç
ç
è
0.9512 0.0488 0
0.0488 0.9492 0.002
0 0.002 0.998
 ö
÷
÷
ø
  VKL2,K=1-3: u,c,t L=1-3: d,s,b  
SCABI Value of sin2qCabibbo 0.0488
EPOLN(1) Electron and positron beam 0.0
EPOLN(2) polarizations in DIS and e+e- 0.0
EPOLN(3) annihilation. First two cmpts are 0.0
PPOLN(1) transverse and only used in e+e-, 0.0
PPOLN(2) 3rd cmpt is longitudinal, and is 0.0
PPOLN(3) +/-1 for fully rh/lh polarized 0.0
QLIM Upper limit on hard process scale 108
THMAX Max. value of thrust in IPROC=110-116 0.9
Y4JT Min. jet separation in IPROC=600-656 0.01
DURHAM Use Durham/JADE algorithm in IPROC=600-656 .TRUE.
  Colour interferences in IPROC=600-656:  
IOP4JT(1) qq_ gg 0: neglect, 1: extreme 3142, 2: extreme 4123 0
IOP4JT(2) qq_ qq_ 0: neglect, 1: extreme 4123, 2: extreme 2143 0
BGSHAT Boson-gluon fusion scale (see below) .TRUE.
BREIT Use Breit frame for DIS kinematics .TRUE.
USECMF Use hadron-hadron c.m. frame .TRUE.
NOSPAC Switch off spacelike showers .FALSE.
ISPAC Changes meaning of QSPAC 0
  (see the earlier notes on QSPAC)  
TMNISR Min. value of s^/S for photon ISR 10-4
ZMXISR Max. momentum fraction for photon ISR 1-10-6
ASFIXD Values of fixed as and w=12loge(2)as/p 0.25
OMEGA for Mueller-Tang cross section 0.3
IAPHIG Approx. used in Higgs+jet production 1
  IPROC=2300-2312  
PHOMAS Damp structure functions for off mass-shell 0.0
  photons (0 for no damping)  
PRESPL Preserve longitudinal momentum of hard c.m. .TRUE.
PTMIN Min. pT in hadronic jet production 10.0
PTMAX Max. pT in hadronic jet production 108
PTPOW 1/pTPTPOW for jet sampling 4.0
YJMIN Min. jet rapidity -8.0
YJMAX Max. jet rapidity +8.0
EMMIN Min. dilepton mass in Drell-Yan 10.0
EMMAX Max. dilepton mass in Drell-Yan 108
EMPOW 1/mEMPOW for Drell-Yan sampling 4.0
Q2MIN Min. Q2 in deep inelastic scattering 0
Q2MAX Max. Q2 in deep inelastic scattering 1010
Q2POW 1/Q2Q2POW for DIS sampling 2.5
YBMIN Min. and max. Bjorken-y in DIS 0.0
YBMAX   1.0
WHMIN Min. hadronic mass in 0.0
  g-induced processes (inc. DIS)  
ZJMAX Max. z in J/y production 0.9
Q2WWMN Min. and max. Q2 in 0.0
Q2WWMX Equivalent Photon Approximation 4.0
YWWMIN Min. and max. photon light-cone fraction 0.0
YWWMAX in Equiv. Photon Approx. 1.0
CSPEED Speed of light in vacuum (mm/s) 2.99792×1011
GEV2NB Value of (h_ c/e)2 389 379
IBSH Number of shots for initial max. weight search 10 000
IBRN(1) 1st random number seed for max. weight search 1246579
IBRN(2) 2nd random number seed for max. weight search 8447766
NQEV Number of entries in Sudakov FF  
  look-up table 1024
ZBINM Max. bin size for z in spacelike branching 0.05
NZBIN Max. number of z bins in spacelike branching 100
NBTRY Max. number of attempts to branch a parton 200
NCTRY Max. number of attempts to decay a cluster 200
NETRY Max. number of attempts to generate a mass 200
NSTRY Max. number of attempts at soft subprocess 200
ACCUR Precision for soft gaussian integration 10-6
RPARTY R-parity conservation in SUSY .TRUE.
SUSYIN Check to see if SUSY data are already loaded .FALSE.
LRSUSY Unit for reading SUSY data (if needed) 66
SYSPIN Spin correlations in decays .TRUE.
THREEB SUSY three body decays .TRUE.
FOURB SUSY four body decays .FALSE.
TAUDEC Tau decay package (HERWIG or TAUOLA) HERWIG
LHSOFT Generation of soft event for Les Houches interface .TRUE.
LHGLSF Self-connected gluons for Les Houches interface .FALSE.
OPTM Optimisation of phase space .FALSE.
IOPSTP Number of steps for phase space optimisation 10
IOPSH Number of weights for phase space optimisation 1000

Table 19: Control switches, constants and options.

Printout options are listed in table 20.
IPRINT = 0 Print program title only
1 Print selected input parameters
2 1 + table of particle codes and properties
3 2 + tables of Sudakov form factors

Table 20: Printout options.

The contents of /HEPEVT/ can by printed by calling HWUEPR, those of /HWPART/ (the last parton shower) by calling HWUBPR. The logical variable PRNDEC (default .TRUE. unless NMXHEP>9999) causes track numbers in event listings to be printed in decimal, or hexadecimal if false. The latter is necessary for very large events such as those generated by the HERBVI package (see above).

The maximum number of errors MAXER refers to errors from which the program cannot recover without killing an event and starting a new one. Such errors are not necessarily a cause for grave concern because the phase space for backward evolution of initial-state showers is complicated and the program may occasionally step outside it (in which case the event weight should be zero anyway). When generating large numbers of events, it is advisable to increase MAXER in proportion, e.g. to MAXEV/100.

See section 8.2 on form factors for details of LRSUD, LWSUD and SUDORD.

The parameter EFFMIN sets the minimum allowed efficiency for the generation of unweighted events. A warning is printed once in every 10/EFFMIN weights if the efficiency is below 10×EFFMIN, and running is stopped if the efficiency is below EFFMIN. See sect. 6.1 for details of NEGWTS.

Variables HARDME and SOFTME invoke hard and soft matrix-element corrections respectively, as described in subsection 3.2.3.

If BGSHAT is .TRUE., the scale used for heavy quark production via boson-gluon fusion in lepton-hadron collisions will be the hard subprocess c.m. energy s^. If it is .FALSE., the scale used will be
2
^
s
 
^
t
 
^
u
 
^
s
 
 2+
^
t
 
 2+
^
u
 
 2
  ,
except in the case of J/y+g production, where u^ is used.

If BREIT is true, the kinematic reconstruction of deep inelastic events takes place in the Breit frame (i.e. the frame where the exchanged boson is purely spacelike, and collinear with the incoming hadron). In fact the reconstruction procedure is invariant under longitudinal boosts, so any frame in which the boson and hadron are collinear would be equivalent, and it is only the transverse part of the boost that has an effect. The BREIT frame option becomes very inaccurate for very small Q2. It is therefore only used if Q2 > 10-4 (the lab and Breit frames are anyway equivalent for such small Q2). If BREIT is false, reconstruction takes place in the lab frame.

If USECMF is true, the entire event record is boosted to the hadron-hadron c.m. frame before event processing, and boosted back afterwards. This means that fixed-target simulation can be done in the lab frame, i.e. with PBEAM2=0. For hadronic processes with lepton beams, this boosting is always done, regardless of the value of USECMF.

In version 6.5, a new logical input variable, PRESPL [.TRUE.], has been introduced to control whether the longitudinal momentum (PRESPL = .TRUE.), or rapidity (PRESPL = .FALSE.), of the hard process centre-of-mass is preserved in hadron collisions after initial-state parton showering. At present the only function of this variable is to allow users to study the effects of momentum reshuffling, which is necessary after showering to compensate for jet masses. In future, it is anticipated that setting PRESPL=.FALSE. will simplify the treatment of other processes in MC@NLO (see sect. 9.6).

The quantities from PTMIN to ZJMAX control the region of phase space in which events are generated and importance sampling inside those regions. See section 8.3.2 on event weights for further details on these quantities and the use of WGTMAX and NOWGT.

If hadronic processes with lepton beams are requested, the photon emission vertex includes the full transverse-momentum-dependent kinematics (the Equivalent Photon Approximation). The variables Q2WWMN and Q2WWMX set the minimum and maximum virtualities generated respectively. For normal simulation, Q2WWMN should be zero, and Q2WWMX should be the largest Q2 through which the lepton can be scattered without being detected. The variables YWWMIN and YWWMAX control the range of lightcone momentum fraction generated.

In addition there are options to give different weights to the various flavours of quarks and diquarks, and to resonances of different spins. So far, these options have not been used. See the comments in the initialisation routine HWIGIN for details.

6.1  Negative weights option

In a number of new applications such as MC@NLO (see sect. 9.6), parton configurations with negative weight are used to produce more accurate predictions, and therefore the possibility of negative event weights has to be considered. In general, a Monte Carlo program generates Nw weights {wi} such that the estimated cross section is
s =
1
Nw
Nw
å
i=1
 wi º
w
  .     (4)
The corresponding error depends on the width of the weight distribution:
ds
s
 =
1
( Nw )
1/2
 
 
d wrms
w
  .     (5)

If only positive weights are generated, and there exists a maximum weight wmax, then unweighted events can be generated by `hit and miss': w'i = 0 or 1 with probability P(w'i = 1) = wi/wmax. Then
ds
s
 = æ
ç
ç
ç
è
wmax-
w
Nw
w
 ö
÷
÷
÷
ø
1/2




 
 = æ
ç
ç
è
wmax-
w
Ne wmax
 ö
÷
÷
ø
1/2



 
 ~
1
( Ne )
1/2
 
 
     (6)
where Ne=Nw w/wmax is the number of events generated. The time needed (especially for detector simulation) depends mainly on the number of events. Hence the inefficiency of `hit and miss' is not necessarily a disaster. This is the usual approach adopted in Monte Carlo event generators.

Negative weights can be generated by subtraction procedures for matrix element corrections. These are not a problem of principle but prevent naive `hit and miss'. To generalize `hit and miss', one can generate unweighted events (w'i = 1) and `antievents' (w'i = -1) with
  sign(w'i) = sign(wi) ,     (7)
  P(w'i = ± 1) = |wi|/|w|max .     (8)
Then
ds
s
 = æ
ç
ç
ç
è
|w|
 |w|max -
w
 2
Nw
w
 2
 ö
÷
÷
÷
ø
1/2




 
 ~
1
( Ne )
1/2
 
 
|w|
w
     (9)
where Ne=Nw |w|/|w|max is now the total number of events+antievents generated. Again, the time needed is almost proportional to Ne, so this is tolerable as long as |w|~w. The cross section after any cuts that may be applied is
  sc =
|w|
Ne
 (N+ - N-)     (10)
where N+ events and N- antievents pass the cuts.

To allow for the possibility of negative weights, a new logical parameter NEGWTS has been introduced. The default (NEGWTS=.FALSE.) is as before: negative weights are forbidden. If one is detected, a non-fatal warning is issued and the event weight is set to zero.

If NEGWTS=.TRUE., negative weights are allowed. Statistics are computed and printed accordingly. If unweighted events are requested (NOWGT=.TRUE.), the initial search stores the maximum and mean absolute weights, |w|max and |w|. Events and antievents are selected according to Eqs. (7,8) and EVWGT is reset to |w|sign(wi), so that the numerator in Eq. (10) is the sum of EVWGTs for contributing (anti)events.

7  Particle data

From HERWIG version 5.9 onwards, new 8-character particle names have been introduced and the revised 7 digit PDG numbering scheme, as advocated in the LEP2 report [26], has been adopted. All hadron and lepton masses are given to five significant figures whenever possible.

Unstable hadrons from clusters produced in both the hard and soft components of the event decay according to simplified decay schemes, which can be tabulated by specifying the print option IPRINT=2. Decays modes are `invented' where necessary to make the branching ratios add up to 100%. Phase space distributions are assumed except where stated otherwise. See section 3.3 for the treatment of heavy quark decays. After a t, partonic b or quarkonium decay, secondary parton showers are produced by outgoing partons as discussed in ref. [11]; these are hadronized in the same way as primary jets.

There have been a number of additions/changes to the default hadrons included via HWUDAT. Here the identification of hadrons follows the PDG  [38, table 13.2], numbered according to their section 30.

All S and P wave mesons are present including the 1P0 and 3P1 states and many new, excited B**, Bc and quarkonium states. Also all D wave kaons and some `light' I=3 states [p2, r(1700) and r3]. All the baryons (singlet/octet/decuplet) containing up to one heavy (c,b) quark are included.

New isoscalars states have been added to try to complete the 13D3, 11D2 and 13D1 multiplets (table 21).
IDHW RNAME IDPDG IDHW RNAME IDPDG
395 OMEGA_3 227 396 PHI_3 337
397 ETA_2(L) 10225 398 ETA_2(H) 10335
399 OMEGA(H) 30223      

Table 21: New isoscalar states.

Also the states in table 22 have been redefined.
IDHW RNAME IDPDG IDHW RNAME IDPDG
57 FH_1 20333      
293 F0P0 9010221 294 FH_00 10221
62 A_0(H)0 10111 290 A_00 9000111
63 A_0(H)+ 10211 291 A_0+ 9000211
64 A_0(H)- -10211 292 A_0- -9000211

Table 22: Re-identified/replaced states.

The f1(1420) state completely replaces the f1(1520) in the 13P0 multiplet, taking over 57. The f0(1370) (294) replaces the f0(980) (293) in the 13P0 multiplet; the latter is retained as it appears in the decays of several other states. The new a0(1450) states (62-64) replace the three old a0(980) states (290-292) in the 13P0 multiplet; the latter are kept as they appear in f1(1285) decays.
ETAMIX h/h',
PHIMIX w/f,
H1MIX h1(1170)/h1(1380),
F0MIX f0(1300)/f0(980),
F1MIX f1(1285)/f1(1510),
F2MIX f2/f2'.

Table 23: Mixing angles.

By default production of the f0(980) and a0(980) states in cluster decays is vetoed.

The mixing angles (in degrees) of all the light, I=0 mesons can now be set using table 23.

There were previously some inconsistencies and ambiguities in our conventions for the mixing of flavour `octet' and `singlet' mesons. They are now as in table 24.
Multiplet Octet Singlet Mixing Angle
11S0 h h' ETAMIX=-23.
13S1 f w PHIMIX=+36.
11P1 h1(1380) h1(1170) H1MIX=ANGLE
13P0 missing f0(1370) F0MIX=ANGLE
13P1 f1(1420) f1(1285) F1MIX=ANGLE
13P2 f'2 f2 F2MIX=+26.
11D2 h2(1645) h2(1870) ET2MIX=ANGLE
13D1 missing w(1600) OMHMIX=ANGLE
13D3 f3 w3 PH3MIX=+28.

Table 24: Mixed meson states.

After mixing, the quark content of the physical states is given in terms of the mixing angle, q, by table 25 where tanq0=
( 2 )
1/2
 
 
.
State (dd_+uu_)/
( 2 )
1/2
 
 
 ss_
Octet cos(q+q0) -sin(q+q0)
Singlet sin(q+q0) cos(q+q0)

Table 25: Quark content of mixed states.

Hence, using the default value of ANGLE=arctan(1/
( 2 )
1/2
 
 
)=+35.3 for q gives ideal mixing, that is, the `octet' state=ss_ and the `singlet'=(dd_+uu_)/
( 2 )
1/2
 
 
. This choice is important to avoid large isospin violations in the 13P0 and 13D1 multiplets in which the octet member is unknown.

Since version 6 contains a large number of supersymmetry processes many new hypothetical particles have been added - see section 3.5. These new states do not interfere with the user's ability to add new particles as described below.

The resonance decay tables supplied in the program have also been largely revised. Measured/expected modes with branching fraction at or above 1 per mille are given, including 4 and 5 body decays. To print the new tables call HWUDPR.

The layout of HWUDAT has been altered to make it easier to identify and modify particle properties. Three new arrays have been introduced RLTIM, RSPIN and IFLAV. These are: the particle's lifetime (s), spin, and a code which specifies the flavour content of each hadron - used (in HWURES) to create sets of iso-flavour hadrons for cluster decay. Using the standard numbering of quark flavours the convention is: Some parts of the program have been automated so that it is possible for the user to add new particles by specifying their properties via the arrays in /HWPROP/ and /HWUNAM/ and increasing NRES appropriately: this should be done before a call to HWUINC.

As an example, the following lines add an isoscalar, spin-2 state 'STAN' and a (very light) stable toponium state 'BEER' with the decay mode: STAN®BEER+BEER+BEER. 
      NRES=NRES+1
      RNAME(NRES)='STAN    '
      IDPDG(NRES)=666
      IFLAV(NRES)=11
      ICHRG(NRES)=0.
      RMASS(NRES)=0.5
      RLTIM(NRES)=1.000D-10
      RSPIN(NRES)=2.0
      NRES=NRES+1
      RNAME(NRES)='BEER    '
      IDPDG(NRES)=66
      IFLAV(NRES)=66
      ICHRG(NRES)=0.
      RMASS(NRES)=0.1
      RLTIM(NRES)=1.000D+30
      RSPIN(NRES)=0.0
      CALL HWMODK(666,1.D0,0,66,66,66,0,0)
Using the logical arrays VTOCDK and VTORDK the production of specified particles can be stopped in both cluster decays and via the decay of other unstable resonances.

A priori weights for the relative production rates in cluster decays of mesons and baryons differing only via their S and L quantum numbers can be supplied using SNGWT and DECWT for singlet (i.e. L-like) and decuplet baryons and REPWT for mesons. The old VECWT now corresponds to REPWT(0,1,0) and TENWT to REPWT(0,2,0).

The arrays FBTM, FTOP and FHVY which stored the branching fractions of the bottom, top and heavier quarks' `partonic' decays are now no longer used. Such decays are specified in the same way as all other decay modes: this permits different decays to be given to individual heavy hadrons. Partonic decays of charm hadrons and quarkonium states are also now supported. As already mentioned, the products' order in a partonic decay mode is significant: see discussion in section 3.3.

The structure of the program has been altered so that the secondary hard subprocess and subsequent fragmentation associated with each partonic heavy hadron decay appears separately. Thus pre-hadronization top quark decays are treated individually, as are any subsequent bottom hadron partonic decays.

Additionally decays of heavy hadrons to exclusive non-partonic final states are supported. No check against double counting from partonic modes is included. However this is not expected to be a major problem for the semi-leptonic and 2-body hadronic modes supplied.

B decays can also be performed by the EURODEC or CLEO Monte Carlo packages. The new variable BDECAY controls which package is used: 'HERW' for HERWIG; 'EURO' for EURODEC; 'CLEO' for CLEO. The EURODEC package can be obtained from the CERN library. The CLEO package is available by kind permission of the CLEO collaboration.

An array NME has been introduced to enable a possible matrix element to be specified for each decay mode. The list of matrix elements currently supported is modest; users are urged to contact an author to have others implemented. It should be noted however that a number of additional matrix elements are now available via the spin correlation algorithm.

A Z' has been introduced with PDG code 32, HERWIG identifier 202, default mass 500 GeV, width GAMZP (default 5 GeV) and name 'Z0PR '. It is invoked by setting ZPRIME=.TRUE. (default .FALSE.).

The decay tables can be written to/read from a file by using HWIODK, adopting the format advocated in the LEP2 report [26]. In addition to the PDG numbering of particles the HERWIG numbers or character names can be used. This permits easy alteration of the decay tables. In HWUINC a call is made to HWUDKS which sets up HERWIG internal pointers and performs some basic checks of the decay tables. Each decay mode must conserve charge and be kinematically allowed and not contain vetoed decay products. The sum of all branching ratios is set to 1 for all particles. Also a warning is printed if an antiparticle does not have all the charge conjugate decays modes of the particle.

HWMODK enables changes to the decay tables to be made by altering or adding single decay modes including on an event-by-event basis. This can be done before calling HWUINC, in which case when altering the branching ratio and/or matrix element code of an existing mode a warning is given of a duplicate second mode which superseeds the first. Branching ratios set below 10-6 are eliminated, whilst if one mode is within 10-6 of unity all other modes are removed. Note that some forethought is required if the branching ratios of two modes of the same particle are changed since the operation of rescaling the branching ratio sum to unity causes a non-commutativity in the order of the calls.

It is possible to create particle property and event listings in any combination of 3 formats - standard ASCII, LATEX or html. These options are controlled by the logical variables PRNDEF (default .TRUE.) PRNTEX (default .FALSE.) and PRNWEB (default .FALSE.). The ASCII output is directed to stout (screen/log file) as in previous versions. When a listing of particle properties is requested (IPRINT.GE.2 or HWUDPR is called explicitly) then the following files are produced:
If (PRNTEX): HW_decays.tex
If (PRNWEB): HW_decays/index.html
/PART0000001.html etc.
The HW_decays.tex file is written to the working directory whilst the many **.html files appear in the sub-directory HW_decays/ which must have been created previously. Paper sizes and offsets for the LATEX output are stored at the top of the block data file HWUDAT: they may need modifying to suit a particular printer. When event listings are requested (NEVHEP.LE.MAXPR or HWUEPR is called explicitly) the following files are created in the current working directory:
If (PRNTEX): HWEV_*******.tex
If (PRNWEB): HWEV_*******.html
where *******=0000001 etc. is the event number.

Note that the html file automatically makes links to the index.html file of particle properties, assumed to be in the HW_decays sub-directory.

A new integer variable NPRFMT (default 1) has been introduced to control how many significant figures are shown in each of the 3 event outputs. Basically NPRFMT=1 gives short compact outputs whilst NPRFMT=2 gives long formats.

Note that all the LATEX files use the package longtable.sty to format the tables. Also if NPRFMT=2 or PRVTX=.TRUE. then the LATEX files are designed to be printed in landscape mode.

8  Structure and output

8.1  Main program

The main program HWIGPR has the following form: 
      PROGRAM HWIGPR
C---COMMON BLOCKS ARE INCLUDED AS FILE HERWIG65.INC
      INCLUDE 'HERWIG65.INC'
      INTEGER N
      EXTERNAL HWUDAT
C---MAX NUMBER OF EVENTS THIS RUN
      MAXEV=100
C---BEAM PARTICLES
      PART1='P'
      PART2='P'
C---BEAM MOMENTA
      PBEAM1=7000.
      PBEAM2=PBEAM1
C---PROCESS
      IPROC=3000
C---INITIALISE OTHER COMMON BLOCKS
      CALL HWIGIN
C---USER CAN RESET PARAMETERS AT
C   THIS POINT, OTHERWISE DEFAULT
C   VALUES IN HWIGIN WILL BE USED.
      PRVTX=.FALSE.
      MAXER=MAXEV/100
      MAXPR=0
      PTMIN=100.
C   N.B. TO READ SUDAKOV FORM FACTOR FILE ON UNIT 77
C   INSERT THE FOLLOWING TWO LINES IN SUBSEQUENT RUNS
C      LRSUD=77
C      LWSUD=0
C---READ IN SUSY INPUT FILE, IN THIS CASE LHC SUGRA POINT 2
      OPEN(UNIT=LRSUSY,FORM='FORMATTED',STATUS='OLD',ERR=999,
     &     FILE='sugra_pt2.1200.in')
      CALL HWISSP
      CLOSE(UNIT=LRSUSY)
C---COMPUTE PARAMETER-DEPENDENT CONSTANTS
      CALL HWUINC
C---CALL HWUSTA TO MAKE ANY PARTICLE STABLE
      CALL HWUSTA('PI0     ')
C---USER'S INITIAL CALCULATIONS
      CALL HWABEG
C---INITIALISE ELEMENTARY PROCESS
      CALL HWEINI
C---LOOP OVER EVENTS
      DO 100 N=1,MAXEV
C---INITIALISE EVENT
      CALL HWUINE
C---GENERATE HARD SUBPROCESS
      CALL HWEPRO
C---GENERATE PARTON CASCADES
      CALL HWBGEN
C---DO HEAVY OBJECT DECAYS
      CALL HWDHOB
C---DO CLUSTER FORMATION
      CALL HWCFOR
C---DO CLUSTER DECAYS
      CALL HWCDEC
C---DO UNSTABLE PARTICLE DECAYS
      CALL HWDHAD
C---DO HEAVY FLAVOUR HADRON DECAYS
      CALL HWDHVY
C---ADD SOFT UNDERLYING EVENT IF NEEDED
      CALL HWMEVT
C---FINISH EVENT
      CALL HWUFNE
C---USER'S EVENT ANALYSIS
      CALL HWANAL
  100 CONTINUE
C---TERMINATE ELEMENTARY PROCESS
      CALL HWEFIN
C---USER'S TERMINAL CALCULATIONS
      CALL HWAEND
      STOP
 999  WRITE (6,*)
      WRITE (6,*) 'SUSY input file did not open correctly.'
      WRITE (6,*) 'Please check that it is in the right place.'
      WRITE (6,*) 'Examples can be obtained from the ISAWIG web page.'
      WRITE (6,*)
      END
The declaration EXTERNAL HWUDAT is recommended to help the linker with finding the block data on some systems.

Various phases of the simulation can be suppressed by deleting the corresponding subroutine calls, or different subroutines may be substituted. For example, in non-SUSY studies the call to HWISSP should be omitted, and in studies at the parton level everything from CALL HWCFOR to CALL HWMEVT can be omitted.

Note that the functionality of the routine HWUINE in earlier versions has now been split between it and a new routine, HWUFNE. A call to the latter must be made between the calls to HWMEVT and HWANAL. A check is built in to prevent execution if this is not done.

The analysis routine HWANAL should always begin with the line 
      IF (IERROR.NE.0) RETURN
since if an event is cancelled, each of the routines is still called in turn until reaching the end of the main loop.

8.2  Form factor file

HERWIG uses look-up tables of Sudakov form factors for the evolution of initial- and final-state parton showers. These can be read from an input file rather than being recomputed each time. The reading, writing and computing of form factor tables is controlled by integer parameters LRSUD and LWSUD (table 26).
LRSUD=N>0 Read form factors for this run from unit N
LRSUD = 0 Compute new form factor tables for this run
LRSUD < 0 Form factor tables are already loaded
LWSUD=N>0 Write form factors on unit N for future use
LWSUD = 0 Do not write new form factor tables

Table 26: Form factor read/write options.

The option LRSUD<0 allows the program to be initialised several times in the same run (e.g. to generate various event types) without recomputing or rereading form factors.

Note that the Sudakov form factors depend on the parameters QCDLAM, VQCUT, VGCUT, NCOLO, NFLAV, RMASS(13) and RMASS(i) for i=1,...,NFLAV. Consequently form factor tables must be recomputed every time any of these parameters is changed. These parameters are written/read with the form factor tables and checks are performed to ensure consistency.

The parton showering algorithm uses the two-loop running coupling, with matching at each flavour threshold. However, the Sudakov table can be computed with either the one-loop or two-loop form, according to the variable SUDORD (= 1 or 2 respectively, default=1). If SUDORD=1 the two-loop value is recovered using the veto algorithm in the shower, whereas if SUDORD=2 no vetoes are used in the final-state evolution. This means that the relative weight of any shower configuration can be calculated in a closed form, hence that showers can be `forced'.

To next-to-leading order the two possibilities should be identical, but they differ at beyond-NLO, so some results may change a little. The most noticeable difference is that the form factor table takes a factor of about five times longer to compute with SUDORD=2 than with 1.

When SUDORD=2, no veto is needed for gluon splitting to quarks. This means that no vetoes are needed for final state showering, except for the previously-mentioned transverse momentum cut. The removal of vetoes allows preselection of the flavours that a jet will contain, giving a huge increase in the efficiency of rare process simulation  [87].

8.3  Event data

/HEPEVT/ is the LEP standard common block containing current event data (table 27).
NEVHEP event number
NHEP number of entries for this event
ISTHEP(I) status of entry I (see below)
IDHEP(I) identity of entry I (Particle Data Group code)
JMOHEP(1,I) pointer to first mother of entry I (see below)
JMOHEP(2,I) pointer to second mother of entry I (see below)
JDAHEP(1,I) pointer to first daughter of entry I (see below)
JDAHEP(2,I) pointer to last daughter of entry I (see below)
PHEP(*,I) (px,py,pz,E,M) of entry I: M=sign
æ
ç
ç
è		 ( abs(m2) )
1/2
 
 
 ,m2 ö
÷
÷
ø
VHEP(*,I) (x,y,z,t) of production vertex of entry I (see sect. 3.8)

Table 27: Current event data.

All momenta are given in GeV/c in the laboratory frame, in which the input beam momenta are PBEAM1 and PBEAM2 as specified by the user and point along the +z and -z directions respectively. Final state particles have ISTHEP(I) = 1. See section 8.3.1 for a complete list of the special status codes used by HERWIG.

The identity codes IDHEP are almost as recommended by the LEP Working Group [26], i.e. the revised PDG [38] numbers where defined, IDHEP = 91 for clusters, 94 for jets, and 0 for objects with no PDG code. The only exception is our use of IDHEP = 26 for the lightest MSSM Higgs boson, to distinguish it from the SM Higgs boson (PDG code 25). In addition, the `generator-specific' (IDHEP=81-100) codes 98 and 99 are used for remnant photons and nucleons, respectively (see section 3.7.2).

HERWIG also has its own internal identity codes IDHW(I), stored in /HWEVNT/. The utility subroutine HWUIDT translates between HERWIG and PDG identity codes. See section 4.13 for further details.

The mother and daughter pointers are standard, except that JMOHEP(2,I) and JDAHEP(2,I) for a parton are its colour mother and colour daughter, i.e. the partons to which its colour and anticolour are connected, respectively. For this purpose the primary partons from a hard subprocess are all regarded as outgoing (see examples in section 4.6.2, 4.13 and 8.5.1). Since a quark has no anticolour, JDAHEP(2,I) is used to point to its flavour partner. Similarly for JMOHEP(2,I) in the case of an antiquark.

In addition to entries representing partons, particles, clusters etc, /HEPEVT/ contains purely informational entries representing the total centre-of-mass momentum, hard and soft subprocess momenta, etc. See section 8.3.1 for the corresponding status codes.

Information from all stages of event processing is retained in /HEPEVT/ so the same particle may appear several times with different status codes. For example, an outgoing parton from a hard scattering (entered initially with status 113 or 114) will appear after processing as an on-mass-shell parton before QCD branching (status 123,124), an off-mass-shell entry representing the flavour and momentum of the outgoing jet (status 143,144), and a jet constituent (157). It might also appear again in other contexts, e.g. as a spectator in a heavy flavour decay (status 154,160).

Incoming partons (entered with status 111 or 112, changed to 121 or 122 after branching) give rise to spacelike jets (status 141 or 142) with m2<0, indicated by PHEP(5,IHEP)<0, due to the loss of momentum via initial-state bremsstrahlung. The same applies in principle to incoming leptons, but in that case emission of at most a single photon is permitted from each initial-state lepton and the off-shell lepton is given status code 3 (see section 8.3.1).

Each parton jet begins with a status 141-144 jet entry giving the total flavour and momentum of the jet. The first mother pointer of this entry gives the location of the parent hard parton, while the second gives that of the subprocess centre-of-mass momentum. If QCD branching has occurred, this is followed by a lightlike CONE entry, which fixes the angular extent of the jet and its azimuthal orientation relative to the parton with which it interferes. The interfering parton is listed as the second mother of the cone. Next come the actual constituents of the jet. If no branching has occurred, there is no cone and the single jet constituent is the same as the jet.

Since version 5.1, the event record has been modified to retain entries for all partons before hadronization (with status ISTHEP=2). During hadronization, the gluons are split into quark-antiquark, while other partons are copied to a location (indicated by JDAHEP(1,*)) where their momenta may be shifted slightly, to conserve momentum, during heavy cluster splitting. Previously the original momenta were shifted, so momentum appeared not to be conserved at the parton level.

Since version 6.3, to take account of the increased energy and complexity of interactions at LHC and future colliders, the default value of the parameter NMXHEP, which sets the array sizes in the standard /HEPEVT/ common block, has been increased to 4000.

8.3.1  Status codes

A complete list of currently-used HERWIG status codes is given below (table 28). Many are used only in intermediate stages of event processing. The most important for users are probably 1 (final-state particle), 101-3 (initial state), 141-4 (jets), and 199 (decayed b-flavoured hadrons).

For technical reasons, some HERWIG status codes ISTHEP between 153 and 165 have changed their meanings since version 5.1.

The event status ISTAT in common /HWEVNT/ is roughly ISTHEP-100 where ISTHEP is the status of entries being processed. For completed events, ISTAT=100.
ISTHEP Description
1 final state particle
2 parton before hadronization
3 documentation line
100 cone limiting jet evolution
101 `beam' (beam 1)
102 `target' (beam 2)
103 overall centre of mass
110 unprocessed hard process c.m.
111 unprocessed beam parton
112 unprocessed target parton
113 unproc. first outgoing parton
114 unproc. other outgoing parton
115 unprocessed spectator parton
120-25 as 110-15, after processing
130 lepton in jet (unboosted)
131-34 as 141-44, unboosted to c.m.
135 spacelike parton (beam, unboosted)
136 spacelike parton (target,unboosted)
137 spectator (beam, unboosted)
138 spectator (target, unboosted)
139 parton from branching (unboosted)
140 parton from gluon splitting (unboosted)
141-44 jet from parton type 111-14
145-50 as 135-40 boosted, unclustered
151 as 159, not yet clustered
152 as 160, not yet clustered
153 spectator from beam
154 spectator from target
155 heavy quark before decay
156 spectator before heavy decay
157 parton from QCD branching
158 parton from gluon splitting
159 parton from cluster splitting
160 spectator after heavy decay
161 beam spectator after gluon splitting
162 target spectator after gluon splitting
163 other cluster before soft process
164 beam cluster before soft process
165 target cluster before soft process
167 unhadronized beam cluster
168 unhadronized target cluster
170 soft process centre of mass
171 soft cluster (beam, unhadronized)
172 soft cluster (target, unhadronized)
173 soft cluster (other, unhadronized)
181 beam cluster (no soft process)
182 target cluster (no soft process)
183 hard process cluster (hadronized)
184 soft cluster (beam, hadronized)
185 soft cluster (target, hadronized)
186 soft cluster (other, hadronized)
190-93 as 195-98, before decays
195 direct unstable non-hadron
196 direct unstable hadron (1-body clus.)
197 direct unstable hadron (2-body clus.)
198 indirect unstable hadron or lepton
199 decayed heavy flavour hadron
200 neutral B meson, flavour at prod'n

Table 28: Status codes.

8.3.2  Event weights

The default is to generate unweighted events (EVWGT=AVWGT). Then event distributions are generated by computing a weight proportional to the cross section and comparing it with a random number times the maximum weight. Set WGTMAX to the maximum weight, or to zero for the program to compute it. If a weight greater than WGTMAX is generated during execution, a warning is printed and WGTMAX is reset. If this occur too often, output event distributions could be distorted.

To generate weighted events, set NOWGT=.FALSE. in common /HWEVNT/.

In QCD hard scattering and heavy flavour and direct photon production (IPROC= 1500-1800) the transverse energy distribution of weighted events (or the efficiency for unweighted events) can be varied using the parameters PTMIN, PTMAX and PTPOW.

Similarly in Drell-Yan processes (IPROC = 1350 etc.) the lepton pair mass distribution is controlled by the parameters EMMIN, EMMAX and EMPOW, and in deep inelastic scattering the Q2 distribution depends on Q2MIN, Q2MAX and Q2POW.

Data on weights generated are output at the end of the run. The mean weight is an estimate of the cross section (in nanobarns) integrated over the region used for event generation. Note that the mean weight is the sum of weights divided by the total number of weights generated, not the total number of events.

The maximum weight is now always printed in full precision. This is needed to be sure of generating the same events in repeated runs.

In versions 6.3 and higher, the option of generating negative weights has been available. This is controlled by the input parameter NEGWTS [.FALSE.], see sect. 6.1 for details.

8.4  Error conditions

Certain combinations of input parameters may lead to problems in execution. HERWIG tries to detect these and print a warning. Errors during execution are dealt with by HWWARN, which prints the calling subprogram and a code and takes appropriate action. In general, the larger the code the more serious the problem. Refer to the source code to find out why HWWARN was called. It is important to note the subprogram from which the call was issued: many different subprograms use the same error code, but each code is unique within a given subprogram.

Events can be rerun by setting the random number seeds NRN(1) and NRN(2) to the values given in the error message or event dump, and MAXWGT to the maximum weight encountered in the run. The contents of /HEPEVT/ can by printed by calling HWUEPR, those of /HWPART/ (the last parton shower) by calling HWUBPR.

Note that if WGTMAX is increased during event generation, so that this type of message is printed: 
   HWWARN CALLED FROM SUBPROGRAM HWEPRO: CODE =   1
   EVENT      21:   SEEDS =  836291635 & 1823648329  WEIGHT = 0.3893E-08
   EVENT SURVIVES. EXECUTION CONTINUES
            NEW MAXIMUM WEIGHT = 0.428217360829367E-08
then to regenerate any later events, WGTMAX must be set to the printed value, as well as setting NRN to the appropriate seeds.

Examples of error messages are: 
   HWWARN CALLED FROM SUBPROGRAM HWSBRN: CODE = 101
   EVENT      31:   SEEDS =  422399901 &  771980111  WEIGHT = 0.3893E-08
   EVENT KILLED.   EXECUTION CONTINUES
Spacelike (initial-state) parton branching had no phase space. This can happen due to cutoffs which are slightly different in the hard subprocess and the parton shower.
Action taken: program throws away this event and starts a new one. 
   HWWARN CALLED FROM SUBPROGRAM HWCHAD: CODE = 102
   EVENT      51:   SEEDS = 1033784787 & 1428957533  WEIGHT = 0.3893E-08
   EVENT KILLED.   EXECUTION CONTINUES
A cluster has been formed with too low a mass to represent any hadron of the correct flavour, and there is no colour-connected cluster from which the necessary additional mass could be transferred.
Action taken: program throws away this event and starts a new one. 
   HWWARN CALLED FROM SUBPROGRAM HWUINE: CODE= 200
   EVENT SURVIVES.  RUN ENDS GRACEFULLY
CPU time limit liable to be reached before generating MAXEV events.
Action taken: skips to terminal calculations using existing events. 
   HWWARN CALLED FROM SUBPROGRAM HWBSUD: CODE= 500
   RUN CANNOT CONTINUE
The table of Sudakov form factors read on unit LRSUD does not extend to the maximum momentum scale QLIM specified for this run.
Action taken: run aborted. The user must either reduce QLIM or set LRSUD to zero to make a bigger table (set LWSUD non-zero to write it). 
   HWWARN CALLED FROM SUBPROGRAM HWBSUD: CODE= 515
   RUN CANNOT CONTINUE
The table of Sudakov form factors read on unit LRSUD is for a different value of a relevant parameter (in this case the b quark mass).
Action taken: run aborted. The user must make a new table (set LWSUD non-zero to write it).

8.5  Sample output

This is the output from the main program listed in section 8.1, with no event printout or user analysis. 
          HERWIG 6.500    16 Oct 2002 

          Please reference:  G. Marchesini, B.R. Webber,
          G.Abbiendi, I.G.Knowles, M.H.Seymour & L.Stanco
          Computer Physics Communications 67 (1992) 465
                             and
          G.Corcella, I.G.Knowles, G.Marchesini, S.Moretti,
          K.Odagiri, P.Richardson, M.H.Seymour & B.R.Webber,
          JHEP 0101 (2001) 010


                                      

          Since SUSY processes are called,
          please also reference: S.Moretti, K.Odagiri,
          P.Richardson, M.H.Seymour & B.R.Webber,
          JHEP 0204 (2002) 028

          Reading in SUSY data from unit 66

          INPUT CONDITIONS FOR THIS RUN

          BEAM 1 (P       ) MOM. =   7000.00
          BEAM 2 (P       ) MOM. =   7000.00
          PROCESS CODE (IPROC)   =    3000
          NUMBER OF FLAVOURS     =    6
          STRUCTURE FUNCTION SET =    8
          AZIM SPIN CORRELATIONS =    T
          AZIM SOFT CORRELATIONS =    T
          QCD LAMBDA (GEV)       =    0.1800
          DOWN     QUARK  MASS   =    0.3200
          UP       QUARK  MASS   =    0.3200
          STRANGE  QUARK  MASS   =    0.5000
          CHARMED  QUARK  MASS   =    1.5500
          BOTTOM   QUARK  MASS   =    4.9500
          TOP      QUARK  MASS   =  175.0000
          GLUON EFFECTIVE MASS   =    0.7500
          EXTRA SHOWER CUTOFF (Q)=    0.4800
          EXTRA SHOWER CUTOFF (G)=    0.1000
          PHOTON SHOWER CUTOFF   =    0.4000
          CLUSTER MASS PARAMETER =    3.3500
          SPACELIKE EVOLN CUTOFF =    2.5000
          INTRINSIC P-TRAN (RMS) =    0.0000
          DECAY SPIN CORRELATIONS=    T
          SUSY THREE BODY ME     =    T
          SUSY FOUR  BODY ME     =    F

          NO EVENTS WILL BE WRITTEN TO DISK

          B_d: Delt-M/Gam =0.7000 Delt-Gam/2*Gam =0.0000
          B_s: Delt-M/Gam = 10.00 Delt-Gam/2*Gam =0.2000

          PDFLIB NOT USED FOR BEAM 1
          PDFLIB NOT USED FOR BEAM 2


          Checking consistency of particle properties


          Checking consistency of decay tables

 Line    1 is the same as line 2634
 Take BR 0.333 and ME code 100 from second entry
 Line    2 is the same as line 2635
 Take BR 0.333 and ME code 100 from second entry
 Line    3 is the same as line 2636
 Take BR 0.111 and ME code 100 from second entry
 Line    4 is the same as line 2637
 Take BR 0.111 and ME code 100 from second entry
 Line    5 is the same as line 2638
 Take BR 0.111 and ME code 100 from second entry
 Line    6 is the same as line 2639
 Take BR 0.333 and ME code 100 from second entry
 Line    7 is the same as line 2640
 Take BR 0.333 and ME code 100 from second entry
 Line    8 is the same as line 2641
 Take BR 0.111 and ME code 100 from second entry
 Line    9 is the same as line 2642
 Take BR 0.111 and ME code 100 from second entry
 Line   10 is the same as line 2643
 Take BR 0.111 and ME code 100 from second entry

          WRITING SUDAKOV TABLE ON UNIT  77

          WRITING MATRIX ELEMENT TABLE ON UNIT  88

          CHECKING SUSY DECAY MATRIX ELEMENTS

          PARTICLE TYPE  21=PI0      SET STABLE

          INITIAL SEARCH FOR MAX WEIGHT

          PROCESS CODE IPROC =        3000
          RANDOM NO. SEED 1  =     1246579
                     SEED 2  =     8447766
          NUMBER OF SHOTS    =       10000
          NEW MAXIMUM WEIGHT =  2.4359288434851706E-03
          NEW MAXIMUM WEIGHT =  3.7017893759328114E-03
          NEW MAXIMUM WEIGHT =  1.9578185739344733E-02

 HWWARN CALLED FROM SUBPROGRAM HWSMRS: CODE =   5
 EVENT       0:   SEEDS =      17673 &      63565  WEIGHT = 0.0000E+00
 EVENT SURVIVES. EXECUTION CONTINUES
 WARNING: MRST98 CALLED WITH Q OUTSIDE ALLOWED RANGE!
 Q VALUE=3.869E+03, MINIMUM=1.118E+00, MAXIMUM=3.162E+03
 NO FURTHER WARNINGS WILL BE ISSUED
          NEW MAXIMUM WEIGHT =  2.2550328399448680E-02
          NEW MAXIMUM WEIGHT =  5.5633613959185035E-02
          NEW MAXIMUM WEIGHT =  6.9166986913034176E-02
          NEW MAXIMUM WEIGHT =  0.1197405440791255    
          NEW MAXIMUM WEIGHT =  0.1488818465259009    

          INITIAL SEARCH FINISHED

          OUTPUT ON ELEMENTARY PROCESS

          N.B. NEGATIVE WEIGHTS NOT ALLOWED

          NUMBER OF EVENTS   =           0
          NUMBER OF WEIGHTS  =       10000
          MEAN VALUE OF WGT  =  4.1213E-03
          RMS SPREAD IN WGT  =  1.2862E-02
          ACTUAL MAX WEIGHT  =  1.3755E-01
          ASSUMED MAX WEIGHT =  1.4888E-01

          PROCESS CODE IPROC =        3000
          CROSS SECTION (PB) =   4.121    
          ERROR IN C-S  (PB) =  0.1286    
          EFFICIENCY PERCENT =   2.768    

          SUBROUTINE TIMEL CALLED BUT NOT LINKED.
          DUMMY TIMEL WILL BE USED. DELETE DUMMY
          AND LINK CERNLIB FOR CPU TIME REMAINING.

          OUTPUT ON ELEMENTARY PROCESS

          N.B. NEGATIVE WEIGHTS NOT ALLOWED

          NUMBER OF EVENTS   =         100
          NUMBER OF WEIGHTS  =        3699
          MEAN VALUE OF WGT  =  3.7473E-03
          RMS SPREAD IN WGT  =  1.2205E-02
          ACTUAL MAX WEIGHT  =  1.3900E-01
          ASSUMED MAX WEIGHT =  1.4888E-01

          PROCESS CODE IPROC =        3000
          CROSS SECTION (PB) =   3.747    
          ERROR IN C-S  (PB) =  0.2007    
          EFFICIENCY PERCENT =   2.517

8.5.1  Guide to sample output

See ref. [1] for a full discussion of the basic features of HERWIG output, including a listing of a sample event. Here we point out only some new features in comparison to version 5.1.

The beam particles, their energies and the process code IPROC=3000 indicate the SUSY process of squark and/or gluino production at LHC. The call to HWISSP triggers a request for a SUSY particle data file, in the format specified in section 3.5, which is read in from the default unit. In this case the file corresponds to the second LHC SUGRA point discussed in section 3.5.1.

First the program lists the main relevant input parameters, including Bd and Bs mixing parameters. Parton distributions were not requested from the PDFLIB library and therefore the default MRST set [51] is used (see sect. 4.1.1.)

Next the program performs some basic checks on the particle data provided. Here it finds that the input file read by HWISSP contains top quark decay modes which duplicate the default modes stored in HWUDAT. The branching ratios and matrix element codes are accordingly updated to those in the input file.

After an initial search for the maximum weight, the program prints its estimate of the relevant cross section and the expected Monte Carlo efficiency for event generation. During the search for the maximum weight the default parton distributions are called with a parameter, in this case the scale, outside the allowed range and a warning is issued. The maximum or minimum value of the parameter is used instead, depending on whether the value requested is too large or too small.

In the course of event generation a new cross section estimate is generated, which is printed together with the actual Monte Carlo efficiency at the end of the run.

8.6  Subroutine descriptions

We give in table 29 a list of all subroutines with their functions. Note that the third letter of the name usually follows a rough classification scheme.
Name Description
  Main program and initialisation
HWIGPR Main program
HWIGIN Default initialisations
  Reading/writing/altering decay modes
HWIGUP Initialisation for Les Houches interface
HWIMDE Adds modes for SUSY four body decays
HWIODK Inputs/outputs formatted decay tables
HWISPC Initialisation of couplings for spin correlations and SUSY decays
HWISPN Initialisation for spin correlations and SUSY decays
HWISP2 Initialisation of two body decays for spin correlations and SUSY decays
HWISP3 Initialisation of three body decays for spin correlations and SUSY decays
HWISP4 Initialisation of four body decays for spin correlations and SUSY decays
HWIPHS Initialises optimized phase space
HWISSP Inputs supersymmetric particle data
HWMODK Modifies or adds an individual decay mode
  User-provided analysis routines
HWABEG Initialises user's analysis
HWAEND Terminates user's analysis
HWANAL Performs user's analysis on event
  Parton branching with interfering gluons
HWBAZF Computes azimuthal correlation functions
HWBCON Makes colour connections between jets
HWBDED Correction to the `dead zone' in e+e-
HWBDIS Correction to the `dead zone' in DIS
HWBDYP Correction to the `dead zone' in Drell-Yan
HWBFIN Transfers external lines of jet to /HEPEVT/
HWBGEN Finds unevolved partons and generates jets
HWBGUP Makes colour and flavour connections for Les Houches events
HWBJCO Combines jets with correct kinematics
HWBMAS Computes masses and trans. momenta in jet
HWBRAN Generates a timelike parton branching
HWBRCN Replaces HWBCON if R-parity is violated
HWBRC1 Finds colour partner in gluino decay
HWBRC2 Finds colour partner in jet
HWBSPA Computes momenta in spacelike jet
HWBSPN Computes spin density/decay matrices
HWBSU1 First term in quark Sudakov form factor
HWBSU2 Second term in quark Sudakov form factor
HWBSUD Computes (or reads) Sudakov form factors
HWBSUG Integrand in gluon Sudakov form factor
HWBSUL Logarithmic part of Sudakov form factor
HWBTIM Computes momenta in timelike jet
HWBTOP Correction to the `dead zone' in top decay
HWBVMC Virtual mass cutoff for parton type ID
  Cluster hadronization model
HWCBCT Cuts a massive baryon cluster in two
HWCBVI Clusters quarks from a ¬  B interaction
HWCBVT Finds which ¬  B interaction partons came from
HWCCCC Finds colour connections after gluon splitting if ¬  B
HWCCUT Cuts a massive cluster in two
HWCDEC Decays clusters into primary hadrons
HWCFLA Sets up flavours for HWCHAD
HWCFOR Forms clusters
HWCGSP Splits gluons
HWCHAD Decays a cluster into one or two hadrons
  Particle and heavy quark decays
HWD2ME Finds maximum weight for a two body decay for spin correlations
HWD3ME SUSY three body decays and spin correlations master routine
HWD3M0 Calculation of SUSY three body matrix element
HWD3M1 Helicity amplitude for f® fff_ via vector boson exchange
HWD3M2 Helicity amplitude for f® fff_ via Higgs exchange
HWD3M3 Helicity amplitude for f® fff_ via antisfermion exchange
HWD3M4 Helicity amplitude for f® fff_ via sfermion exchange
HWD3M5 Helicity amplitude for f_®f_ff_ via vector boson exchange
HWD3M6 Helicity amplitude for f®f ff_ via vector boson exchange
HWD3M7 Helicity amplitude for f®G~ff_ via vector boson exchange
HWD3M8 Helicity amplitude for f®f_f_f via scalar exchange (1st diagram)
HWD3M9 Helicity amplitude for f®f_f_f via scalar exchange (2nd diagram)
HWD3MA Helicity amplitude for f®f_f_f via scalar exchange (3rd diagram)
HWD3MB Helicity amplitude for f® fff via scalar exchange (1st diagram)
HWD3MC Helicity amplitude for f® fff via scalar exchange (2nd diagram)
HWD3MD Helicity amplitude for f® fff via scalar exchange (3rd diagram)
HWD3MF Helicity amplitude for f®f_f_f_ via scalar exchange (1st diagram)
HWD3MG Helicity amplitude for f®f_f_f_ via scalar exchange (2nd diagram)
HWD3MH Helicity amplitude for f®f_f_f_ via scalar exchange (3rd diagram)
HWD3MI Helicity amplitude for f_®f_ff_ via scalar exchange
HWD4ME Master routine for SUSY four body decays
HWD4M0 Matrix Element for f® V*V*® ff_ff_
HWDBOS Finds and decays W and Z bosons
HWDBOZ Chooses decay mode of W and Z bosons
HWDBZ2 Copy of HWDBOZ used by hadronic WW, WZ and ZZ production
HWDCHK Checks given decay mode is self-consistent
HWDCLE Interface to CLEO package for B decays
HWDEUR Interface to EURODEC package for B decays
HWDFIV Generates a five-body decay
HWDFOR Generates a four-body decay
HWDHAD Generates decays of unstable hadrons
HWDHGC Higgs ®gg decay
HWDHGF Higgs ® W+ W- decay
HWDHIG Finds and decays Higgs bosons
HWDHOB Finds and decays heavy objects
HWDHO1 Selects decay mode for heavy object
HWDHO2 Calculates momenta of heavy object's decay products
HWDHO3 Makes colour connects for heavy object decay
HWDHO4 Performs parton shower for heavy object decay
HWDHO5 Checks for colour disconnections in heavy object decay
HWDHO6 Perform ¬  Rp colour connections in heavy object decay
HWDHVY Finds and decays heavy flavours
HWDHWT Subroutine for top decay via a virtual H±
HWDRCL Colour connections for a ¬  B SUSY decay
HWDRME Main ¬  Rp 3-body decay matrix element subroutine
HWDRM1 ¬  Rp 3-body decay matrix element subroutine
HWDRM2 ¬  Rp 3-body decay matrix element subroutine
HWDRM3 ¬  Rp 3-body decay matrix element subroutine
HWDRM4 ¬  Rp 3-body decay matrix element subroutine
HWDRM5 ¬  Rp 3-body decay matrix element subroutine
HWDPWT Phase space three-body decay weight
HWDSIN Performs decays with spin correlations
HWDSI1 Picks next particle to decay for spin correlations
HWDSI2 Pick moment of decay products with spin correlations
HWDSI3 Selects t polarization and passes it to the decay routine
HWDSM2 Two body decay with spin correlations
HWDSM3 Three body decay with spin correlations
HWDSM4 Four body decay with spin correlations
HWDTAU Interface to TAUOLA for t decays
HWDTHR Generates a three-body decay
HWDTOP Decides whether to decay top quark
HWDTWO Generates a two-body decay
HWDWWT Weak (V-A) three-body decay weight
HWDXLM Tests if decay vertex lies in given volume
  Elementary subprocess generation
HWECIR Interface to CIRCE
HWEFIN Final calculations on elementary subprocess
HWEGAM Generates incoming photon
HWEINI Initialises elementary subprocess
HWEISR Generates photon emission from e or µ
HWEONE Sets up a 2®1 hard subprocess
HWEPRO Generates elementary subprocess
HWETWO Sets up a 2®2 hard subprocess
  Individual hard subprocesses
HWH2BK Matrix element for bb_® W±H±
HWH2BH Matrix element for H± production via bq-fusion
HWH2DD Function to return the D function of [68]
HWH2F1 Subroutine to return the F function of [69] for a fixed first momentum
HWH2F2 Subroutine to return the F function of [69] for a fixed second momentum
HWH2F3 Subroutine to return the F function of [69] for all first and second momenta
HWH2HE Matrix element for Higgs associated production
HWH2M0 Subroutine to compute the massless matrix element for QQ_ Z
HWH2MQ Subroutine to compute the massive matrix element for QQ_ Z
HWH2PS Subroutine to perform the phase-space for Z+two jets
HWH2P1 Subroutine to select quark masses for HWH2PS
HWH2P2 Subroutine to select quark masses for HWH2PS
HWH2QH Matrix element for qq_,gg® QQ_(') Higgs
HWH2SH Matrix element for squark pair plus Higgs production
HWH2SS Subroutine to return the S function of [68]
HWH2T1 Function to return the T1 function of [68]
HWH2T2 Function to return the T2 function of [68]
HWH2T3 Function to return the T3 function of [68]
HWH2T4 Function to return the T4 function of [68]
HWH2T5 Function to return the T5 function of [68]
HWH2T6 Function to return the T6 function of [68]
HWH2T7 Function to return the T7 function of [68]
HWH2T8 Function to return the T8 function of [68]
HWH2T9 Function to return the T9 function of [68]
HWH2T0 Function to return the T10 function of [68]
HWH2VH Matrix element for qq_(') ® V Higgs, V=W±,Z0
HWH4JT Hard subprocess: e+e-® 4 jets
HWH4J1 Matrix element for e+e-® 4 jets
HWHBGF Hard subprocess: boson-gluon fusion (BGF)
HWHBKI Computes kinematics for BGF
HWHBRN Returns a phase-space point for BGF
HWHBSG Computes cross section for BGF
HWHDIS Hard subprocess: deep inelastic e/µ quark
HWHDYP Hard subprocess: Drell-Yan Z0/g production
HWHDYQ Subroutine for QQ_ Z
HWHEGG Hard subprocess: two-photon processes in e+e-
HWHEGW Hard subprocess: g W processes in e+e-
HWHEGX Calculates cross section for HWHEGW
HWHEPA Hard subprocess: e+e- ® q q_
HWHEPG Hard subprocess: e+e- ® q q_ g
HWHESG Gaugino pair production in e+e- collisions
HWHESL Slepton pair production in e+e- collisions
HWHESP Sparticle pair production in e+e- collisions
HWHESQ Squark pair production in e+e- collisions
HWHEW0 e+e- ® W+ W- subroutine
HWHEW1 e+e- ® W+ W- subroutine
HWHEW2 e+e- ® W+ W- subroutine
HWHEW3 e+e- ® W+ W- subroutine
HWHEW4 e+e- ® W+ W- subroutine
HWHEW5 e+e- ® W+ W- subroutine
HWHEWW Hard subprocess: e+e- ® W+ W-
HWHGBP Main routine for gauge boson pair production in hadron-hadron
HWHGBS Phase space for gauge boson pair production in hadron-hadron
HWHGB1 Selects gauge boson mass for HWHGBS
HWHGB2 Matrix element for WW in hadron-hadron
HWHGB3 Matrix element for ZZ in hadron-hadron
HWHGB4 Matrix element for WZ in hadron-hadron
HWHGB5 Selects t and u for HWHGBS
HWHGRV Graviton resonance production
HWHGUP External hard process using Les Houches interface
HWHHVY Hard subprocess: heavy quark production
HWHIBG Hard subprocess: for bg ® Q Higgs, with Q=t,b
HWHIBK Hard subprocess: bb_® W±H±
HWHIBQ Subroutine for H± production via bq-fusion
HWHIG1 Matrix elements for Higgs + jet production
HWHIGA Amplitudes squared for Higgs + jet
HWHIGB Loop integrals for Higgs + jet
HWHIGE Hard subprocess: Higgs associated production
HWHIGH Hard subprocess: qq_® Higgs1 + Higgs2
HWHIGJ QCD Higgs + jet production
HWHIGM Choose any Higgs mass for production routines
HWHIGQ Hard subprocess: qq_,gg® QQ_(') Higgs
HWHIGS Hard subprocess: g g/qq_ ® Higgs
HWHIGT Computes gg ® Higgs cross section
HWHIGV Hard subprocess: qq_(') ® V Higgs, V=W±,Z0
HWHIGW Hard subprocess: W+W-/Z0Z0® Higgs
HWHIGY Computes e+e-® Z0 ® Z0 HSM0 cross section
HWHIGZ Hard subprocess: l+ l- ® Z0 ® Z0 HSM0
HWHIHH Hard subprocess: l+l- SUSY Higgs pair production
HWHISQ Subroutine for squark pair plus Higgs production
HWHPH2 Hard subprocess: direct photon pairs
HWHPHO Hard subprocess: direct photon production
HWHPPB Box contribution to gg®gg
HWHPPE Pointlike photon-parton (fixed flavour)
HWHPPH Pointlike photon-parton (fixed pair flavour)
HWHPPM Pointlike photon-parton direct light meson
HWHPPT Pointlike photon-parton (all flavours)
HWHPQS Pointlike photon-quark (Compton) scattering
HWHQCD QCD 2® 2 hard subprocesses
HWHQCP Identifies QCD 2® 2 hard subprocess
HWHQCM Hard subprocess: gg® qq_, l+l-, W+W-
HWHRBB ¬  Rp resonant squark to SM particles
HWHRBS ¬  Rp resonant squark to SUSY particles
HWHREE ¬  Rp SM particle production in e+e- collisions
HWHREM Treats hard scattering remnants
HWHREP Decides which ¬  Rp subroutine to use in e+e-
HWHRES ¬  Rp single sparticle production in e+e- collisions
HWHRLL ¬  Rp resonant slepton to SM particles
HWHRLS ¬  Rp resonant slepton to SUSY particles
HWHRSP Decides which ¬  Rp subroutine to use in hadron-hadron
HWHRSS Identifies ¬  Rp process
HWHSCT Process extra hard scatterings
HWHSNG Colour singlet parton scattering
HWHSNM Colour singlet parton scattering matrix element
HWHSPN Main routine for spin correlations in the hard process
HWHS01 Helicity amplitude for ff_® ff via s-channel vector boson exchange
HWHS02 Helicity amplitude for ff_® ff via t-channel scalar exchange
HWHS03 Helicity amplitude for ff_® ff via u-channel scalar exchange
HWHS04 Helicity amplitude for ff_® ff_ via s-channel vector boson exchange
HWHS05 Helicity amplitude for gg® ff_ via t-channel fermion exchange
HWHS06 Helicity amplitude for gg® ff_ via u-channel fermion exchange
HWHS07 Helicity amplitude for gg® ff_ via s-channel gluon exchange
HWHS08 Helicity amplitude for gf® ff via t-channel scalar exchange
HWHS09 Helicity amplitude for gf_® f f via t-channel scalar exchange
HWHS10 Helicity amplitude for gf® ff via s-channel fermion exchange
HWHS11 Helicity amplitude for gf_® f f via s-channel antifermion exchange
HWHS12 Helicity amplitude for gf® ff via u-channel fermion exchange
HWHS13 Helicity amplitude for gf_® f f via u-channel fermion exchange
HWHS14 Helicity amplitude for gg® ff via t-channel fermion exchange
HWHS15 Helicity amplitude for gg® ff via u-channel fermion exchange
HWHS16 Helicity amplitude for gg® ff via s-channel gluon exchange
HWHS17 Helicity amplitude for ff® ff via t-channel gauge boson exchange
HWHS18 Helicity amplitude for ff_® ff_ via t-channel gauge boson exchange
HWHS19 Helicity amplitude for f_f®f_f via t-channel gauge boson exchange
HWHS20 Helicity amplitude for f_f_®f_f_ via t-channel gauge boson exchange
HWHS21 Helicity amplitude for ff_® ff_ via s-channel scalar exchange
HWHS22 Helicity amplitude for ff_® ff_ via t-channel scalar exchange
HWHS23 Helicity amplitude for ff_® ff_ via u-channel scalar exchange
HWHS24 Helicity amplitude for f_ f ® ff via s-channel scalar exchange
HWHS25 Helicity amplitude for f_ f ® ff via t-channel scalar exchange
HWHS26 Helicity amplitude for f_ f ® ff via u-channel scalar exchange
HWHS27 Helicity amplitude for ff® ff_ via s-channel scalar exchange
HWHS28 Helicity amplitude for ff® ff_ via t-channel scalar exchange
HWHS29 Helicity amplitude for ff® ff_ via u-channel scalar exchange
HWHS30 Helicity amplitude for f_f_® ff via s-channel scalar exchange
HWHS31 Helicity amplitude for f_f_® ff via t-channel scalar exchange
HWHS32 Helicity amplitude for f_f_® ff via u-channel scalar exchange
HWHS33 Helicity amplitude for ff® ff via s-channel scalar exchange
HWHS34 Helicity amplitude for f_f_®f_f_ via s-channel scalar exchange
HWHSS1 Gaugino-gaugino production matrix element
HWHSS2 Gaugino-gaugino production matrix element with polarizations
HWHSSG Gaugino-gaugino/gaugino-sparton production
HWHSSL Slepton pair production
HWHSSP Combines MSSM subprocesses
HWHSSQ SQCD 2® 2 hard subprocesses
HWHSSS Identifies MSSM hard subprocess
HWHV1J Hard subprocess W/Z + jet production
HWHV2J Master subroutine for all gauge boson + two jet processes
HWHVVJ Dummy WW,WZ,ZZ production subroutine (see sect. 9.6)
HWHWEX Top production by W exchange
HWHWPR Hard subprocess: W production
  Soft minimum-bias or underlying event
HWMEVT Generates min bias or soft underlying event
HWMLPS Generates longitudinal phase space
HWMNBI Computes negative binomial probability
HWMULT Chooses min bias charged multiplicity
HWMWGT Calculates weight for minimum bias events
  Random number generators
HWRAZM Randomly rotated azimuth
HWREXP Random number: exponential distribution
HWREXQ Random number: exp. dist. with cutoff
HWREXT Random number: exponential transverse mass
HWRGAU Random number: gaussian
HWRGEN Random number generator (l'Ecuyer's method)
HWRINT Random integer
HWRLOG Random logical
HWRPIP Random primary interaction point
HWRPOW Random number: power distribution
HWRUNG Random number: uniform + gaussian tails
HWRUNI Random number: uniform
  Spacelike branching of incoming partons
HWSBRN Generates spacelike parton branching
HWSDGG Drees-Grassie gluon distribution in photon
HWSDGQ Drees-Grassie quark distribution in photon
HWSFBR Chooses a spacelike branching
HWSFUN Hadron structure functions
HWSGAM Gamma function (for structure functions)
HWSGEN Generates x values for spacelike partons
HWSGQQ Inserts g® qq_ part of gluon form factor
HWSMRS Subroutine for MRST PDFs
HWSSPC Replaces spacelike partons by spectators
HWSSUD Sudakov form factor/structure function
HWSTAB Interpolates in function table (for HWSSUD)
HWSVAL Checks for valence parton
  Miscellaneous utilities
HWUAEM Running electromagnetic coupling constant
HWUAER Real part of photon self-energy
HWUALF Two-loop QCD running coupling constant
HWUANT Finds a particle's antiparticle
HWUATS Replaces & with
HWUBPR Prints branching data for last parton shower
HWUBST Boost event record to/from hadron-hadron c.m.f.
HWUCFF Coefficients for e+e- and DIS cross sections
HWUCI2 Logarithmic integral Ci2
HWUDAT Particle properties (N.B. BLOCK DATA)
HWUDKL Generates decay vertex of unstable particle
HWUDKS Converts decay modes into internal format
HWUDPR Prints particle properties and decay modes
HWUECM Centre-of-mass energy
HWUEDT Insert or delete entries in the event record
HWUEEC Computes coefficients for e+e- cross section
HWUEMV Moves entries within the event record
HWUEPR Prints event data
HWUFNE Finishes an event
HWUGAU Adaptive gaussian integration
HWUGUP Run termination for Les Houches Accord
HWUIDT Translates particle identity codes
HWUINC Initial parameter-dependent calculations
HWUINE Initialises an event
HWULB4 Boost: rest frame to lab, no masses assumed
HWULDO Lorentz 4-vector dot product
HWULF4 Boost: lab frame to rest, no masses assumed
HWULI2 Logarithmic integral Li2 (Spence function)
HWULOB Lorentz transformation: rest frame ® lab
HWULOF Lorentz transformation: lab ® rest frame
HWULOR Multiplies by Lorentz matrix
HWUMAS Puts mass in 5th component of vector
HWUMBW Generates mass (Breit-Wigner distribution)
HWUMPO Spinor routine
HWUMPP Spinor routine
HWUNST Converts integer to character
HWUPCM Centre-of-mass momentum
HWUPUP Prints contents of Les Houches common block
HWURAP Rapidity
HWURES Computes/prints resonance data
HWUROB Rotation by inverse of matrix R
HWUROF Rotation by matrix R
HWUROT Computes rotation R from vector to z-axis
HWUSOR Sorts an array in ascending order
HWUSPR Print contents of spin correlations common block
HWUSQR Square root with sign retention
HWUSTA Makes a particle type stable
HWUTAB Interpolates in a table
HWUTIM Checks time remaining
  Vector manipulation
HWVDIF Vector difference
HWVDOT Vector dot product
HWVEQU Vector equality
HWVSCA Vector times scalar
HWVSUM Vector sum
HWVZRI Vector zero (integer)
HWVZRO Vector zero
HWWARN Issues warnings and deals with errors

Table 29: Subroutines.

In addition there are the routines for generating the Schuler-Sjöstrand parton distributions of the photon: 
   SASANO  SASBEH  SASDIR  SASGAM  SASVMD
and for putting the output of TAUOLA [88] into the /HEPEVT/ common block: 
   FILHEP
Finally, there are dummy versions of external routines, which should be deleted if the relevant packages are used:

9  Interfaces

9.1  Les Houches

Version 6.5 includes support for the interface between parton level generators and HERWIG using the Les Houches Accord as described in [89]. In general we have tried to code the interface in such a way that from the user point of view it behaves in the same way as that already included in PYTHIA [90].

The interface operates in the following way. If the IPROC code is set negative then HERWIG assumes that the user wants an external hard process using the Les Houches accord.

If this option is used the initialisation will call the routine UPINIT to initialise the external hard process; this name is the same as that used by PYTHIA. This routine should set the values of the run parameters in the Les Houches run common block.

After the initialisation during the event generation phase the routine UPEVNT (again this name is the same as that used by PYTHIA) is called. This routine should fill the event common block as described in [89].

Dummy copies of both these routines are supplied with HERWIG and should be deleted and replaced if you are using this option. Due to the internal structure of HERWIG two new parameters are needed to control the interface, in addition to those in the Les Houches common block. The logical input variable LHSOFT [.TRUE] controls the generation of the soft underlying event; the default is to generate a soft underlying event and so LHSOFT must be set .FALSE. for lepton-lepton processes. The second variable LHGLSF [.FALSE.] controls the treatment of colour self-connected gluons, which may occur with some methods of colour decomposition. The default is to kill events with self-connected gluons whereas if LHGLSF = .TRUE. an information-only warning is issued instead.

We have not included full support for the interface. In particular These restrictions are the same as those imposed by PYTHIA. If you have any problems with the interface or need the options we have not yet supported please let us know.

It should also be noted that while the interface has been tested with Standard Model processes, for example all the processes in ALPGEN [91], it is less well tested for SUSY processes.

9.2  CIRCE

Simulation of beamstrahlung is now available, via an interface to the CIRCE program [92]. It is implemented as a modification of the standard bremsstrahlung implementation, so is only available for processes for which that is available. The produced radiation is treated in the collinear limit, i.e. COLISR and the option WWA in routine HWEGAM are forcibly set .TRUE..

For both types of distribution function, fe/e and fg/e, we use a simplifying approximation: where the correct result should use the convolution of beamstrahlung and bremsstrahlung, we actually use the sum of the two. For the photon distribution this is quite reasonable, but for the electron it is perhaps more questionable. It boils down to the replacement:
fe/e(x) = ó
õ
1


x
dz
z
 fe/e,brem ( z ) fe/e,beam æ
ç
ç
è
x
z
 ö
÷
÷
ø		 ® fe/e,brem(x) + fe/e,beam(x) - d(1-x).     (11)
This is good to the extent that both distributions are dominated by the x®1 region. A small utility program can be obtained from the HERWIG web page to test this approximation. For example, for TESLA running at 500 GeV, the integral over fe/e(x) for large x, relevant for threshold shapes, is accurate to 1% for x>0.95 and to 10% for x>0.99, and the value of fe/e((90/500)2), relevant for radiative return to the Z, is accurate to 2%.

In order for the distribution to be probabilistic, the coefficient of the delta function has to remain positive. If it is not, HERWIG will terminate. This can be prevented by adjusting the ZMXISR parameter: the default value, ZMXISR=0.999999 is too large and HERWIG will fail, but reducing it to ZMXISR=0.99999 works fine for all the parameter sets currently available. It is straightforward to check that ZMXISR values in this range do not have an influence on any physical observables, by rerunning with it further reduced.

The interface is controlled by five new variables, which are given in the following table together with their default values:
CIRCOP 0
CIRCAC 2
CIRCVR 7
CIRCRV 9999 12 31
CIRCCH 0
CIRCOP is the main control option: CIRCOP=0 switches off beamstrahlung and uses standard HERWIG; CIRCOP=1 switches to collinear kinematics, but still leaves beamstrahlung switched off; CIRCOP=2 uses only beamstrahlung; and CIRCOP=3 uses both beamstrahlung and bremsstrahlung. CIRCOP=0 and CIRCOP=3 should therefore be regarded as `off' and `on', respectively, with the other two options mainly for cross-checking purposes. The variables CIRCAC, CIRCVR, CIRCRV and CIRCCH are simply passed to CIRCE as its input variables acc, ver, rev and chat, as described in its documentation. The default values correspond to the most up-to-date revision of version 7 of the TESLA parametrization, with minimal output.

CIRCE can be obtained from
http://heplix.ikp.physik.tu-darmstadt.de/nlc/beam.html

9.3  TAUOLA

An interface to the TAUOLA decay package [88] has been added. To use this interface, the dummy subroutines DEXAY, INIETC, INIMAS, INIPHX, INITDK, PHOINI, PHOTOS must be deleted. You should then link to both TAUOLA and PHOTOS. The easiest way to do this is to obtain the TAUOLA/PHOTOS versioning system from
http://wasm.home.cern.ch/wasm/test.html
and produce a version with the correct size of the HEPEVT common block. This system will provide a JETSET demo. If you link HERWIG with the code produced by this apart from the JETSET interface and main program everything should work.

This interface uses the information from the spin density matrices to select the helicity of the decaying t leptons if the spin correlation algorithm is being used; otherwise the helicity of the decaying t is averaged over.

9.4  ISAWIG

ISAWIG is an interface which allows SUSY spectra and decay tables generated by ISAJET [39] to be read into HERWIG. The web address is
http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/

To coincide with the release of HERWIG 6.5 we have produced a new version of ISAWIG. Since the original version of ISAWIG there have been a number of new versions of ISAJET, often with changes in the sizes of the common blocks from which we extract the information we need. This has necessitated periodical updating of the code to run with the most recent version of ISAJET, with the result that the code could no longer be run with older versions. However, the Snowmass points and slopes [93] are defined with ISAJET version 7.58 and so we can no longer continue in this way and still be able to generate these points. Therefore, starting with the new ISAWIG version 1.2, we are using C preprocessing so that users can define the version of ISAJET they are using at compile time. When compiling the main ISAWIG code and the modified SUGRUN and SSRUN programs the following compiler options should be specified: The default at the moment is to compile code to run with ISAJET7.64. In the future this will change so that the most recent version of ISAJET becomes the default but we will continue to support the older versions.

9.5  HDECAY

There has also been a new release of the HDECAY package [94], version 3. In order to support this version, and the previous version 2.0, we are also using C preprocessing to control the version of HDECAY used. In order to achieve this the HDECAY interface code has been merged with the main ISAWIG program. The following options should be used when compiling ISAWIG if you are using HDECAY: As before the default is to compile a dummy routine. If you are using HDECAY you must link with a version of the HDECAY code which has the main HDECAY program removed.

9.6  MC@NLO

The program MC@NLO [95] generates events with At present only the processes of gauge boson pair production in hadron (h=p or p_) collisions are implemented in MC@NLO, making use of the HERWIG process codes IPROC=2850-2880 (see table 9).

To use these codes, one must delete the dummy HERWIG subroutine HWHVVJ and replace it by the one in the MC@NLO package. This will then read and process an input file of parton configurations prepared by MC@NLO. For further details see the MC@NLO web page:
http://www.hep.phy.cam.ac.uk/theory/webber/MCatNLO/

Implementation of other processes, such as heavy quark production, in MC@NLO is in progress [96]. These future development will make use of the Les Houches interface discussed in sect. 9.1, rather than special process codes.

9.7  JIMMY

Contrary to our original plans, the JIMMY generator [49] for multiple parton-parton interactions in hadron-hadron, hadron-photon and photon-photon collisions is still not integrated into HERWIG, and is still a separate add-on package. It can be obtained from the JIMMY web page:
http://jetweb.hep.ucl.ac.uk/JIMMY/index.html
where references to the relevant papers can also be found.

9.8  HERBVI

An interface is provided to the HERBVI package for electroweak baryon number violation (¬  B), and other multi-W± production processes [44]. For full details, see

http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/HERBVI/

Acknowledgments

We thank G. Abbiendi, J. Chýla and L. Stanco for their collaboration on earlier versions of HERWIG and for their continuing interest. We also thank U. Baur, M.J. Gibbs, E.W.N. Glover, Z. Kunszt and D.R. Ward for contributing code to be incorporated into the program.

We are grateful to the many experimental colleagues who have made valuable criticisms and suggestions over the years, especially K. Hamacher, R. Hemingway and G. Rudolph.

Since this is expected to be the last (Fortran) HERWIG version, we would like to take this opportunity to thank again all those colleagues and users who have contributed to the development of the program, whether by providing code, suggesting improvements or reporting problems. Thanks also to those who worked so hard to establish the Les Houches accord, which should make it possible to expand the application of Fortran HERWIG to new processes without (much) further intervention by the authors. Meanwhile, many of us will be transferring our main efforts to HERWIG++, in preparation for the new era of object-oriented event generation.

References

[1]
G. Marchesini, B.R. Webber, G. Abbiendi, I.G. Knowles, M. H. Seymour and L. Stanco, HERWIG: a Monte Carlo event generator for simulating hadron emission reactions with interfering gluons. Version 5.1 - april 1991, Comput. Phys. Commun. 67 (1992) 465.

[2]
A. Bassetto, M. Ciafaloni and G. Marchesini, Jet structure and infrared sensitive quantities in perturbative QCD, Phys. Rep. 100 (1983) 201.

[3]
G. Marchesini and B. R. Webber, Simulation of QCD jets including soft gluon interference, Nucl. Phys. B 238 (1984) 1.

[4]
B.R. Webber, A QCD model for jet fragmentation including soft gluon interference, Nucl. Phys. B 238 (1984) 492.

[5]
B.R. Webber, Monte Carlo simulation of hard hadronic processes, Ann. Rev. Nucl. Part. Sci. 36 (1986) 253.

[6]
R.K. Ellis, G. Marchesini and B.R. Webber, Soft radiation in parton parton scattering, Nucl. Phys. B 286 (1987) 643.

[7]
G. Marchesini and B.R. Webber, Monte Carlo simulation of general hard processes with coherent QCD radiation, Nucl. Phys. B 310 (1988) 461.

[8]
G. Marchesini and B.R. Webber, Associated transverse energy in hadronic jet production, Phys. Rev. D 38 (1988) 3419.

[9]
I.G. Knowles, Spin correlations in parton-parton scattering, Nucl. Phys. B 310 (1988) 571.

[10]
I.G. Knowles, A linear algorithm for calculating spin correlations in hadronic collisions, Comput. Phys. Commun. 58 (1990) 271.

[11]
G. Marchesini and B.R. Webber, Simulation of QCD coherence in heavy quark production and decay, Nucl. Phys. B 330 (1990) 261.

[12]
G. Marchesini and B.R. Webber, Simulation of QCD initial state radiation at small x, Nucl. Phys. B 349 (1991) 617.

[13]
S. Catani, B.R. Webber and G. Marchesini, QCD coherent branching and semiinclusive processes at large x, Nucl. Phys. B 349 (1991) 635.

[14]
G. Abbiendi and L. Stanco, Heavy flavor production at HERA. Simulation with a new Monte Carlo event generator, Z. Physik C 51 (1991) 81; A new heavy flavor generator in e-p collisions, Comput. Phys. Commun. 66 (1991) 16.

[15]
M.H. Seymour, Photon radiation in final state parton showering, Z. Physik C 56 (1992) 161.

[16]
M.H. Seymour, Matrix element corrections to parton shower algorithms, Comput. Phys. Commun. 90 (1995) 95.

[17]
M.H. Seymour, A simple prescription for first order corrections to quark scattering and annihilation processes, Nucl. Phys. B 436 (1995) 443 [ hep-ph/9410244 ].

[18]
M.H. Seymour, Matrix element corrections to parton shower simulation of deep inelastic scattering, contributed to 27th International Conference on High Energy Physics (ICHEP), Glasgow 1994, Lund preprint LU-TP-94-12, unpublished.

[19]
G. Corcella and M.H. Seymour, Matrix element corrections to parton shower simulations of heavy quark decay, Phys. Lett. B 442 (1998) 417 [ hep-ph/9809451 ].

[20]
G. Corcella, M.H. Seymour, Initial state radiation in simulations of vector boson production at hadron colliders, Nucl. Phys. B 565 (2000) 227 [ hep-ph/9908388 ].

[21]
P. Richardson, Spin correlations in Monte Carlo simulations, J. High Energy Phys. 0111 (2001) 029 [ hep-ph/0110108 ].

[22]
S. Moretti, K. Odagiri, P. Richardson, M. H. Seymour and B. R. Webber, Implementation of supersymmetric processes in the HERWIG event generator, J. High Energy Phys. 0204 (2002) 028 [ hep-ph/0204123 ].

[23]
D. Amati and G. Veneziano, Preconfinement as a property of perturbative QCD, Phys. Lett. B 83 (1979) 87.

[24]
Yu.L. Dokshitzer and S.I. Troyan, Nonleading perturbative corrections to the dynamics of quark-gluon cascades and soft hadron spectra in e+ e- annihilation, Leningrad Nuclear Physics Institute preprint N 922 (1984);
Y.I. Azimov, Yu.L. Dokshitzer, V.A. Khoze and S.I. Troian, The string effect and QCD coherence, Phys. Lett. B 165 (1985) 147; Similarity of parton and hadron spectra in QCD jets, Z. Physik C 27 (1985) 65.

[25]
T. Sjöstrand et al., QCD event generators for LEP, in Z Physics at LEP, G. Altarelli, R. Kleiss and C. Verzegnassi eds., CERN 89-08.

[26]
I.G. Knowles, T. Sjöstrand et al. Report of QCD Event Generators Working Group, in Physics at LEP2, ed. G. Altarelli, T. Sjöstrand and F. Zwirner, CERN 96-01.

[27]
S. Gieseke, Event generators - new developments, to appear in Proc. 14th Topical Conference on Hadron Collider Physics (Karlsruhe, October 2002).

[28]
Yu.L. Dokshitzer, V.A. Khoze, S.I. Troyan and A.H. Müller, QCD coherence in high-energy reactions, Rev. Mod. Phys. 60 (1988) 373.

[29]
CDF Collaboration, Evidence for color coherence in pp_ collisions at s1/2 = 1.8 TeV, Phys. Rev. D 50 (1994) 556.

[30]
DØ collaboration,Color coherent radiation in multijet events from ppbar collisions at s1/2=1.8 TeV Phys. Lett. B 414 (1997) 419 [ hep-ex/9706012 ];
Evidence of color coherence effects in W+ jets events from pp_ collisions at s1/2 = 1.8 TeV, Phys. Lett. B 464 (1999) 145 [ hep-ex/9908017 ].

[31]
K. Odagiri, Color connection structure of (supersymmetric) QCD (2®2) processes, J. High Energy Phys. 10 (1998) 006 [ hep-ph/9806531 ].

[32]
Yu.L. Dokshitzer, V.A. Khoze and S.I. Troian, New perturbative results in hadron jet physics, in Proc. Physics in Collision 6, Chicago 1986, p.417.

[33]
A. Einstein, B. Podolski and N. Rosen, Can quantum mechanical description of physical reality be considered complete?, Phys. Rev. 47 (1935) 777.

[34]
G. Marchesini, QCD coherence in the structure function and associated distributions at small x, Nucl. Phys. B 445 (1995) 49 [ hep-ph/9412327 ].

[35]
T. Sjöstrand, High-energy physics event generation with PYTHIA 5.7 and JETSET 7.4, Comput. Phys. Commun. 82 (1994) 74.

[36]
G. Corcella and S. Moretti, in preparation.

[37]
M.H. Seymour, The Higgs boson line shape and perturbative unitarity, Phys. Lett. B 354 (1995) 409 [ hep-ph/9505211 ].

[38]
D. E. Groom et al., Review of particle physics, Eur. Phys. J. C 15 (2000) 1.

[39]
F.E. Paige, S.D. Protopopescu, H. Baer and X. Tata, ISAJET 7.37: a Monte Carlo event generator for pp, p_ p and e+e- reactions, [ hep-ph/9804321 ]; ISAJET 7.40: a Monte Carlo event generator for pp, p_p, and e+e- reactions preprint BNL-HET-98/39, FSU-HEP-981016, UH-511-917-98, October 1998, [ hep-ph/9810440 ]. H. Baer, F.E. Paige, S.D. Protopopescu, X. Tata, preprint BNL-HET-99-43, FSU-HEP-991218, UH-511-952-00, [ hep-ph/0001086 ]. Latest version available from ftp://penguin.phy.bnl.gov:/pub/isajet/.

[40]
H.E. Haber and G.L. Kane, The search for supersymmetry: probing physics beyond the Standard Model, Phys. Rep. 117 (1985) 75.

[41]
J.F. Gunion, H. E. Haber, G.L. Kane and S. Dawson, The Higgs hunter guide, Addison-Wesley, Reading MA 1990.

[42]
H. Dreiner, P. Richardson and M.H. Seymour, Parton-shower simulations of R-parity violating supersymmetric models, J. High Energy Phys. 04 (2000) 008 [ hep-ph/9912407 ].

[43]
H. Dreiner, An introduction to explicit R-parity violation in Perspectives on Supersymmetry, G.L. Kane ed., World Scientific, pp. 462-479 [ hep-ph/9707435 ].

[44]
M.J. Gibbs, A. Ringwald, B.R. Webber and J.T. Zadrozny, Monte Carlo simulation of baryon and lepton number violating processes at high-energies, Z. Physik C 66 (1995) 285 [ hep-ph/9406266 ];
M.J. Gibbs, B.R. Webber, HERBVI: a program for simulation of baryon and lepton number violating processes, Comput. Phys. Commun. 90 (1995) 369 [ hep-ph/9504232 ].

[45]
A. Kupco, Cluster hadronization in HERWIG 5.9, in Proc. Workshop on Monte Carlo Generators for HERA Physics, Hamburg 1998, p.292 [ hep-ph/9906412 ].

[46]
B. Andersson, G. Gustafsson, G. Ingelman and T. Sjöstrand, Parton fragmentation and string dynamics, Phys. Rep. 97 (1983) 33.

[47]
TPC/Two Gamma collaboration, H. Aihara et al., Study of baryon correlations in e+e- annihilation at 29 GeV, Phys. Rev. Lett. 57 (1986) 3140.

[48]
UA5 collaboration, G.J. Alner et al., The UA5 high-energy p_ p simulation program, Nucl. Phys. B 291 (1987) 445.

[49]
J.M. Butterworth and J.R. Forshaw, Photoproduction of multi-jet events at HERA: a Monte Carlo simulation, J. Phys. G 19 (1993) 1657;
J.M. Butterworth, J.R. Forshaw and M.H. Seymour, Multiparton interactions in photoproduction at HERA, Z. Physik C 72 (1996) 637 [ hep-ph/9601371 ].

[50]
H. Plothow-Besch, PDFLIB: a library of all available parton density functions of the nucleon, the pion and the photon and the corresponding alpha-s calculations, Comput. Phys. Commun. 75 (1993) 396.

[51]
A.D. Martin, R.G. Roberts, W.J. Stirling and R.S. Thorne, Scheme dependence, leading order and higher twist studies of MRST partons, Phys. Lett. B 443 (1998) 301 [ hep-ph/9808371 ].

[52]
R.S. Thorne, private communication.

[53]
M. Drees and K. Grassie, Parametrizations of the photon structure and applications to supersymmetric particle production at HERA, Z. Physik C 28 (1985) 451.

[54]
M. Drees and C.S. Kim, Associate production of a heavy quark and a gauge boson at e+ p colliders, Z. Physik C 52 (1991) 503.

[55]
G. A. Schuler and T. Sjostrand, Parton distributions of the virtual photon, Phys. Lett. B 376 (1996) 193 [ hep-ph/9601282 ].

[56]
M. Drees and R.M. Godbole, Virtual photon structure functions and the parton content of the electron, Phys. Rev. D 50 (1994) 3124 [ hep-ph/9403229 ].

[57]
Z Physics at LEP1, CERN 89-09, vol. 3, p. 129.

[58]
H. Burkhardt, F. Jegerlehner, G. Penso and C. Verzegnassi, Uncertainties in the hadronic contribution to the QED vacuum polarization, Z. Physik C 43 (1989) 497.

[59]
R. K. Ellis, D. A. Ross and A. E. Terrano, The perturbative calculation of jet structure in e+e- annihilation, Nucl. Phys. B 178 (1981) 421.

[60]
S. Catani and M.H. Seymour, A general algorithm for calculating jet cross sections in NLO QCD, Nucl. Phys. B 485 (1997) 291 [ hep-ph/9605323 ];
erratum ibid. 510 (1997) 503.

[61]
S. Moretti and W.J. Stirling, Spin correlations in e+e-® 4 jets, Eur. Phys. J. C 9 (1999) 81 [ hep-ph/9808429 ];
S. Moretti, The implementation of the four-jet matrix elements in HERWIG and elsewhere, preprint RAL-TR-2000-039, August 2000 [ hep-ph/0008197 ].

[62]
S. Katsanevas and P. Morawitz, SUSYGEN 2.2: a Monte Carlo event generator for MSSM sparticle production at e+e- colliders, Comput. Phys. Commun. 112 (1998) 227 [ hep-ph/9711417 ].

[63]
J. Kalinowski, R. Ruckl, H. Spiesberger and P.M. Zerwas, Supersymmetry with R-parity breaking: contact interactions and resonance formation in leptonic processes at LEP2, Phys. Lett. B 406 (1997) 314 [ hep-ph/9703436 ].

[64]
A. Djouadi, J. Kalinowski and P. M. Zerwas, Higgs radiation off top quarks in high-energy e+e- colliders, Z. Physik C 54 (1992) 255; S. Komamiya, Searching For Charged Higgs Bosons At O(1/2-TeV To 1-TeV) e+e- Colliders, Phys. Rev. D 38 (1988) 2158; S. Kanemura, S. Moretti and K. Odagiri, Single charged Higgs boson production at next generation linear colliders, J. High Energy Phys. 0102 (2001) 011 [ hep-ph/0012030 ]; Single production of charged Higgs bosons at linear colliders, hep-ph/0101354.
[65]
A. Djouadi, J. Kalinowski and P. M. Zerwas, Measuring the H t t_ coupling in e+e- collisions, Mod. Phys. Lett. A 7 (1992) 1765; A. Juste and G. Merino, Top-Higgs Yukawa coupling measurement at a linear e+e- collider, [ hep-ph/9910301 ] (and references therein); M. Battaglia, A. Ferrari, A. Kiiskinen and T. Maki, Pair production of charged Higgs bosons at future linear e+e- colliders, in Proc. of the APS/DPF/DPB Summer Study on the Future of Particle Physics (Snowmass 2001) ed. R. Davidson and C. Quigg, [ hep-ex/0112015 ]; S. Moretti, Detection of heavy charged Higgs bosons at future linear colliders, [ hep-ph/0206208 ]; Higgs radiation off top-antitop pairs at future Linear Colliders: a background study, [ hep-ph/9911501 ]; S. Moretti, The process e+e-® H t t_ and its backgrounds at future electron positron colliders, Phys. Lett. B 452 (1999) 338 [ hep-ph/9902214 ].

[66]
A.H. Mueller and W.K. Tang, High-energy parton-parton elastic scattering in QCD, Phys. Lett. B 284 (1992) 123.

[67]
J. F. Gunion and Z. Kunszt, Lepton correlations in gauge boson pair production and decay, Phys. Rev. D 33 (1986) 665.
[68]
R. Kleiss and W. J. Stirling, Spinor techniques for calculating pp_® W±/Z0 + jets, Nucl. Phys. B 262 (1985) 235.
[69]
B. van Eijk and R. Kleiss, On the calculation of the exact gg® Z bb_ cross-section including Z decay and b quark mass effects, In Proceeding of the LHC Workshop, Aachen 1990, vol. 2 183-194.
[70]
F. A. Berends, R. Pittau and R. Kleiss, Excalibur: a Monte Carlo program to evaluate all four fermion processes at LEP-200 and beyond, Comput. Phys. Commun. 85 (1995) 437 [ hep-ph/9409326 ].
[71]
R. Kleiss and R. Pittau, Weight optimization in multichannel Monte Carlo, Comput. Phys. Commun. 83 (1994) 141 [ hep-ph/9405257 ].

[72]
E. Eichten, I. Hinchliffe, K. Lane and C. Quigg, Super collider physics, Rev. Mod. Phys. 56 (1984) 579; erratum ibid. 58 (1986) 1065.

[73]
S. Dawson, E. Eichten and C. Quigg, Search for supersymmetric particles in hadron-hadron collisions, Phys. Rev. D 31 (1985) 1581.

[74]
A. Dedes and S. Moretti, Eur. Phys. J. C 10 (1999) 515; preprint RAL-TR-1999-067, June 1999, [ hep-ph/9909526 ];
A. Djouadi et al., preprint February 2000, [ hep-ph/0002258 ].

[75]
A. Dedes and S. Moretti, Phys. Rev. D 60 (1999) 015007;
A. Djouadi, J.-L. Kneur and G. Moultaka, Phys. Rev. Lett. 80 (1998) 1830; Nucl. Phys. B 569 (2000) 53;
G. Bélanger, F. Boudjema and K. K. Sridhar, Nucl. Phys. B 568 (2000) 3.

[76]
A.A. Barrientos-Bendezù and B.A. Kniehl, W± H± associated production at the Large Hadron Collider, Phys. Rev. D 59 (1999) 015009 [ hep-ph/9807480 ];
S. Moretti and K. Odagiri, The phenomenology of W± H± production at the Large Hadron Collider, Phys. Rev. D 59 (1999) 055008 [ hep-ph/9809244 ].

[77]
S. Moretti and K. Odagiri, Phys. Rev. D 55 (1997) 5627.

[78]
M. Spira, Fortsch. Phys. 46 (1998) 203.

[79]
D. Dicus, T. Stelzer, Z. Sullivan and S. Willenbrock, Phys. Rev. D 59 (1999) 094016.

[80]
F. Borzumati, J.-L. Kneur and N. Polonsky, Phys. Rev. D 60 (1999) 115011;
S. Moretti and D.P. Roy, Phys. Lett. B 470 (1999) 209.

[81]
M. Carena el al., Report of the Tevatron Higgs working group, hep-ph/0010338 .

[82]
D.J. Miller, S. Moretti, D.P. Roy and W.J. Stirling, Detecting heavy charged Higgs bosons at the LHC with four b quark tags, Phys. Rev. D 61 (2000) 055011.

[83]
B.C. Allanach, K. Odagiri, M.A. Parker and B.R. Webber, Searching for narrow graviton resonances with the ATLAS detector at the large hadron collider, J. High Energy Phys. 09 (2000) 019.

[84]
A. Donnachie and P.V. Landshoff, Total cross sections, Phys. Lett. B 296 (1992) 227 [ hep-ph/9209205 ].

[85]
S. Moch, A. Ringwald and F. Schrempp, Instantons in deep-inelastic scattering. The simplest process, Nucl. Phys. B 507 (1997) 134 [ hep-ph/9609445 ]; Instanton induced cross-sections in deep inelastic scattering, Phys. Lett. B 438 (1998) 217 [ hep-ph/9806528 ].

[86]
M. Gibbs, A. Ringwald and F. Schrempp, QCD-instanton induced final states in deep inelastic scattering, in Proc. DIS 1995 (Paris), pp. 341-344 [ hep-ph/9506392 ;
A. Ringwald and F. Schrempp, QCDINS 2.0:a Monte Carlo generator for instanton-induced processes in deep-inelastic scattering, Comput. Phys. Commun. 132 (2000) 267 [ hep-ph/9911516 ]

[87]
M.H. Seymour, Heavy quark production in jets, Z. Physik C 63 (1994) 99.

[88]
S. Jadach, Z. Was, R. Decker and J. H. Kuhn, The tau decay library TAUOLA: Version 2.4, Comput. Phys. Commun. 76 (1993) 361; Z. Was,TAUOLA the library for tau lepton decay, and KKMC/KORALB/KORALZ/ status report, Nucl. Phys. 98 (Proc. Suppl.) (2001) 96. [ hep-ph/0011305 ].

[89]
E. Boos et al., Generic user process interface for event generators, [ hep-ph/0109068 ].
[90]
T. Sjostrand, P. Eden, C. Friberg, L. Lonnblad, G. Miu, S. Mrenna and E. Norrbin, High-energy-physics event generation with PYTHIA 6.1, Comput. Phys. Commun. 135 (2001) 238 [ hep-ph/0010017 ]; T. Sjostrand, L. Lonnblad and S. Mrenna, PYTHIA 6.2: Physics and manual, [ hep-ph/0108264 ].
[91]
M. L. Mangano, M. Moretti, F. Piccinini, R. Pittau and A. D. Polosa, ALPGEN, a generator for hard multiparton processes in hadronic collisions, hep-ph/0206293.

[92]
T. Ohl, CIRCE version 1.0: Beam spectra for simulating linear collider physics, Comput. Phys. Commun. 101 (1997) 269 [ hep-ph/9607454 ].



[93]
B. C. Allanach et al., The Snowmass points and slopes: Benchmarks for SUSY searches, in Proc. of the APS/DPF/DPB Summer Study on the Future of Particle Physics (Snowmass 2001) ed. R. Davidson and C. Quigg, Eur. Phys. J. C 25 (2002) 113 [ hep-ph/0202233 ].

[94]
A. Djouadi, J. Kalinowski and M. Spira, HDECAY: A program for Higgs boson decays in the standard model and its supersymmetric extension, Comput. Phys. Commun. 108 (1998) 56 [ hep-ph/9704448 ].

[95]
S. Frixione and B. R. Webber, Matching NLO QCD computations and parton shower simulations, J. High Energy Phys. 0206 (2002) 029 [ hep-ph/0204244 ]; The MC@NLO event generator, [ hep-ph/0207182 ].

[96]
S. Frixione, P. Nason and B. R. Webber, Matching NLO QCD computations and parton showers in heavy quark production, in preparation.
1
Work supported in part by the U.K. Particle Physics and Astronomy Research Council, by the Italian Ministero della Università e Ricerca Scientifica, and by the E.U. Fourth Framework Programme `Training and Mobility of Researchers', Network `Quantum Chromodynamics and the Deep Structure of Elementary Particles', contract FMRX-CT98-0194 (DG 12 - MIHT). Research of G.C. supported by grant number DE-FG02-91ER40685 from the U.S. Department of Energy.
2
The quark mass parameters should also be thought of as effective or constituent masses rather than current quark masses.
3
This applies also to final-state emission, i.e. to jet fragmentation at large values of the jet momentum fraction.
4
We point out that in the expression `soft correction', `soft' refers to the phase space where such corrections are applied and not to the amplitude, since we still use the `hard' exact matrix element for the soft correction as well.
5
Or the equivalent, e.g. fortLRSUSY.dat on a VAX system.
6
Default values for input variables are shown in square brackets.
7
The situation when baryon number is violated is more complicated and is discussed in [44] for the Standard Model and in [42] for R-parity violating SUSY models.
8
Except for those containing a `perturbative' quark when CLDIR = 1 - see below.
9
The WW in parameter names is a relic from earlier versions that used the less accurate Weizsacker-Williams approximation.
This document was translated from LATEX by HEVEA.