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