**********************************************************************
**********************************************************************
**                                                                  **
**                                                September 1991    **
**                                                                  **
**                           A Manual to                            **
**                                                                  **
**           The Lund Monte Carlo for Hadronic Processes            **
**                                                                  **
**                        PYTHIA version 5.6                        **
**                                                                  **
**                        Hans-Uno Bengtsson                        **
**     Department of Theoretical Physics, University of Uppsala     **
**                                                                  **
**                        Torbjorn Sjostrand                        **
**                    CERN/TH, CH-1211 Geneva 23                    **
**                BITNET/EARN address TORSJO@CERNVM                 **
**                       Tel. +22 - 767 28 20                       **
**                                                                  **
**       Copyright Hans-Uno Bengtsson and Torbjorn Sjostrand        **
**                                                                  **
**********************************************************************
**********************************************************************
*                                                                    *
*        Table of Contents                                           *
*                                                                    *
*  1.    Introductory Material                                       *
*  1.1.  Program Objective                                           *
*  1.2.  Update History                                              *
*  1.3.  Major Changes from PYTHIA 4.8 to 5.3                        *
*  1.4.  Major Changes from PYTHIA 5.3 to 5.4                        *
*  1.5.  Major Changes from PYTHIA 5.4 to 5.5                        *
*  1.6.  Major Changes from PYTHIA 5.5 to 5.6                        *
*  1.7.  Installation of Program                                     *
*  1.8.  Programming Philosophy                                      *
*                                                                    *
*  2.    The Program Components                                      *
*  2.1.  The Main Subroutines                                        *
*  2.2.  The Physics Processes                                       *
*  2.3.  Comments on Physics Processes                               *
*  2.4.  Switches for Event Type and Kinematics Selection            *
*  2.5.  The General Switches and Parameters                         *
*  2.6.  General Event Information                                   *
*  2.7.  The Event Record                                            *
*  2.8.  Other Routines and Commonblocks                             *
*  2.9.  The JETSET Routines                                         *
*  2.10. On Cross-Sections                                           *
*  2.11. Examples                                                    *
*  2.12. How to include external processes in PYTHIA                 *
*                                                                    *
*        Acknowledgements                                            *
*        References                                                  *
*                                                                    *
**********************************************************************
* Legend:                                                            *
* >=      larger than or equal to                                    *
* <=      smaller than or equal to                                   *
* /=      not equal to                                               *
* ->      goes to (rightarrow)                                       *
* ^       what follows next is to be read as an upper index          *
* _       what follows next is to be read as a lower index           *
* (D=..)  default value for commonblock parameter                    *
* (R)     commonblock variable which user may read but not change    *
* (I)     commonblock variable for purely internal use               *
**********************************************************************

1.    Introductory Material

PYTHIA is a program intended for the study of hadronic processes.
Historically, this has meant incoming hadron, but today also lepton
beams are included, so that e.g. e+e- and ep events can be generated,
in addition to pp ones. (Further combinations, involving antiprotons,
neutrons, pions, photons, muons, etc., are also possible.) For
completeness, it is also possible to generate a few non-hadronic
final states. In hadronic interactions, the emphasis is on high-pT
physics, but PYTHIA also covers the domain of low-pT interactions
as an integral part of the total cross-section. In its present form it
includes hard scattering matrix elements, structure functions and
initial and final state parton showers. Fragmentation is performed using
the ordinary Lund fragmentation model, JETSET version 7.3, but an
important task for PYTHIA is to set up the correct string configuration,
particularly nontrivial for the low-pT target remnants.
______________________________________________________________________

1.1.  Program Objective

Leaving aside for the moment a discussion of to what extent we can
trust the various pieces that make up a full-fledged Monte Carlo
program like PYTHIA, i.e., assuming that PYTHIA indeed generates
events as they would occur if the underlying theory and assumptions
were true, three important tasks can be singled out as the main
objectives for PYTHIA:

1) to test the underlying theory by comparison with experiments at
   present colliders;
2) to help designing techniques and strategies for probing the
   physics at present and future colliders;
3) to help disentangle possible signals for new physics at colliders
   by giving accurate estimates of standard model backgrounds and
   some scenarios for physics beyond the standard model.

Over the years, PYTHIA has indeed been put to all of these uses; in
particular it has been extensively used to look at the possibilities
of detection of predicted standard model particles like the top quark
and the Higgs boson at LHC and SSC.

It should also be stated what PYTHIA is not. It is not a program for
high-precision tests of the electroweak model, like e.g. Z0 line
shape studies. The program includes an 'improved Born term' approach
(running alpha_em, s-dependent width, etc.) to resonances, and a
modified leading log description of initial and final state radiation,
which should ensure a reasonably accurate representation, but certainly
not without its limitations.
______________________________________________________________________

1.2.  Update History

The PYTHIA program has undergone a contined rapid expansion, in
terms of the number of subprocesses included, and also in terms of
the physics aspects covered. This is a continuous process, with the
official numbered versions little more than snapshots of this process.
Some versions have never been given a larger distribution, but have
instead been handed to a few people for specific purposes. For the
record, we below give a list of the versions, with some brief notes.

no   date   publ.     main new or improved features
1    Dec82  [Ben84]   synthesis of predecessors COMPTON, HIGHPT and
                      KASSANDRA
2    -----
3.1  -----
3.2  -----
3.3  Feb84  [Ben84a]  scale-breaking structure functions
3.4  Sep84  [Ben85]   more efficient kinematics selection
4.1  Dec84            initial and final state parton showers, W and Z
4.2  Jun85            multiple interactions
4.3  Aug85            W W, W Z, Z Z, R
4.4  Nov85            gamma W, gamma Z, gamma gamma
4.5  Jan86            H0 production, diffractive and elastic events
4.6  May86            angular correlation in resonance pair decays
4.7  May86            Z'0, H+
4.8  Jan87  [Ben87]   variable impact parameter in multiple interactions
4.9  May86            g H+
5.1  May86            massive matrix elements for heavy quarks
5.2  Jun87            intermediate boson scattering
5.3  Oct89            new particle and subprocess codes, new commonblock
                      structure, new kinematics selection, some
                      lepton-lepton and lepton-hadron interactions,
                      new subprocesses
5.4  Jun90            s-dependent widths, resonances not on mass-shell,
                      new processes, new structure functions
5.5  Jan91            improved e+e- and ep, several new processes
5.6  Sep91  [this]    reorganized structure functions

Versions preceding 4.8 by now should be considered obsolete. Versions
4.9, 5.1 and 5.2 are all minor expansions on 4.8, and have been given
only limited distribution. Version 5.3 is therefore the first major
revision since 4.8. Since the current version is almost backwards
compatible with 5.3, and since several errors have been found in the
original distribution of 5.3, users are recommended to switch from
5.3 to 5.6.
______________________________________________________________________

1.3.  Major Changes from PYTHIA 4.8 to 5.3

Version 5.3 of PYTHIA, the Lund Monte Carlo for Hadronic Processes,
represented a major break in continuity from a programming point
of view, with no backwards compatibility to PYTHIA 4.8, the former
standard version. Many major rewritings have been prompted by the
adoption of the new particle numbering scheme developed under the aegis
of the Particle Data Group [PDG88]. This scheme is also the standard
adopted for JETSET 7. The internal scheme for subprocess
classification has also been expanded to cover the future inclusion
of further new processes. Also a number of other changes have been
made, to systematize current features and leave space for new ones.

Here is a more extensive list of major changes in PYTHIA 5.3.
- New KF particle codes have been introduced and are consistently used
  in the program.
- The program has, also in other respects, been updated and extended to
  run with JETSET 7.2. In particular, this applies to the commonblock
  LUJETS, which has been expanded to contain more event information.
- The ISUB subprocess classification scheme has been reorganized and
  expanded.
- A number of new subprocesses have been included: heavy quark
  production with massive matrix elements, gauge boson pair scattering,
  etc.
- Many of the subprocesses can also be simulated in lepton-lepton
  interactions (and are prepared for lepton-hadron ones).
- The possibility of event kinematics preselection has been added.
- Switches for event flavour preselection have been reorganized and
  coordinated with JETSET.
- Generation of kinematical variables has been changed and systematized.
- It is possible to generate pileup events, where one single
  bunch crossing gives rise to several hadron-hadron interactions.
______________________________________________________________________

1.4.  Major Changes from PYTHIA 5.3 to 5.4

The updates from version 5.3 to 5.4 are all minor, and just about any
program that ran with version 5.3 will also work with PYTHIA 5.4.
The following changes might give compatibility problems from programs
based on PYTHIA 5.3:
- The routine PYTHIA has been renamed PYEVNT (to avoid confusion
  between the program-as-a-whole and one specific subroutine);
  a dummy routine PYTHIA which calls PYEVNT is provided, however.
- Subprocesses 142 and 143 have been moved to 143 and 144, respectively.
  The related MSEL codes 22 and 23 have also been moved to 23 and 24,
  respectively.
- The non-standard decay channel Z -> H+ + H- has been removed,
  instead the decay Z' -> H+ + H- has been introduced to fill the
  same function.
- The program should be run with JETSET version 7.3 rather than 7.2
  (to take into account the new decay channels and the running of
  alpha_em). Note that many decay channels appear with new numbers.
- Several errors have been corrected, where previously the program gave
  wrong results.

In addition, the following changes have been made, where no
compatibility problems are involved:
- A subprocess 142 has been added; in general the Z'/W' sector has
  been much expanded, and couplings to ordinary fermions and gauge
  bosons appear as free parameters (defined in JETSET).
- Four new structure function parametrization by Morfin and Tung
  [Mor90] come with the program, and two new by Gluck, Reya and Vogt
  [GRV89].
- alpha_em is now taken to be running; if the user so decides it can
  be frozen, however.
- s-channel resonances are generated with shat-dependent widths;
  this should give an improved description of line shapes.
- In the decays H -> Z Z or W W, the contribution to the total
  Higgs width is included also for Higgs so light that either or
  both products are off mass-shell.
- The couplings of H to quarks is taken with running quark masses,
  both for production q + qb -> H0 and decay H0 -> q + qb; see
  MSTP(37) for possibility to use fix masses.
- In a number of 2 -> 2 processes it is now possible to select
  specific mass ranges for resonances in the final state, and in
  particular to generate events below the nominal threshold mass.
- Also for other decays of a resonance into two new resonances it
  is possible to set mass ranges on the decay products. This feature
  is not fully developed yet, however, so some limitations exist.
- Process 22, f + fb -> Z0 + Z0, has been expanded so that the two Z0
  by default represent the correct mixture of gamma*/Z0, with the
  possibility to switch to Z0 only or gamma* only, according to
  MSTP(43) value.
- The process f + q -> f' + Q for single heavy flavour production by
  W exchange in the t channel has been implemented.
- The matrix element for process 73, W_L Z_L -> W_L Z_L, has been
  corrected according to formulae of [Ter90].
- Processes 71, 72, 73, 76 and 77 have been expanded so that they
  alternatively can be used to simulate the strong interaction
  Z_L/W_L models of Dobado, Herrero and Terron [DHT90].
- The process g + g -> Z + q + qbar has been introduced with massive
  matrix elements, as calculated and programmed by Kleiss [Kle90];
  the code used in PYTHIA is a only very slightly modified copy.
- Some processes may now be used also in leptoproduction. This in
  particular includes ZZ, ZW and WW interactions, like in H production.
  PYTEST has been expanded to cover this possibility.
- The table produced with PYSTAT(2) also gives the flavour (KF) codes
  and current masses of resonances, and the decay channel numbers in
  JETSET (IDC) for the channels listed. By default, fourth generation
  channels are not listed.
- A number of internal administrative changes, which should not be
  visible to the outside user.
______________________________________________________________________

1.5.  Major Changes from PYTHIA 5.4 to 5.5

The updates from version 5.4 to 5.5 are all minor, and just about any
program that ran with version 5.4 will also work with PYTHIA 5.5.
The following changes might give compatibility problems from programs
based on PYTHIA 5.4:
- The routine PYTHIA has now been removed; all calls to this routine
  should henceforward be replaced by calls to PYEVNT.
- The expansion of process 11 to also include electroweak interactions
  has been split off and is now found as process 10.
- The size of the XPQ array in PYSTFU and PYSTFE has been expanded
  from -6:6 to -25:25.
- Some COEF coefficients have been moved so as to make space for the
  e+e- motivated new coefficients.
- Some MINT and VINT variables have been moved to leave space for new
  variables.

In addition, the following changes have been made, where no
compatibility problems are involved:
- The use of symmetrization factors for identical final state
  particles has been systematized; see section 2.10.
- It is possible to have a photon as incoming particle, with the
  quark and gluon structure functions of Drees and Grassie [Dre85].
- Initial state gamma radiation in e+e- annihilation has been
  included in the effective structure function approach.
  In an approximate fashion, one can also this way include effects
  of beamstrahlung. See MSTP(11).
- Parton (quark/gluon) structure functions for the electon have
  been implemented by numerical convolution of the photon structure
  function inside an electron with the parton structure functions
  inside a photon. Therefore a number of parton-parton scattering
  processes, otherwise available only for hadron-hadron collisions,
  can also be simulated in e+e- and ep events. See MSTP(12).
- Decays Z' -> Z + H and W' -> W + H have been included in the framework
  of left-right-symmetric models (see [Coc90]). The strength of these
  couplings can be set with PARU(145) and PARU(146), respectively.
- Several processes for leptoquark production have been implemented,
  ISUB numbers 145, 162, 163, and 164. The LQ -> q + l coupling
  strength can be modifed in PARU(151).
- Process 103, gamma + gamma -> H0, has been added.
- A number of photon-induced subprocesses have been included:
  f gamma -> f g (33), f gamma -> f gamma (34), g gamma -> q qbar (54),
  gamma gamma -> f fbar (58), g gamma -> Q Qbar (massive) (84), and
  gamma gamma -> F Fbar (massive) (85).
- Process 10 has been introduced to include neutral and/or charged
  current interactions in the t-channel between partons, e.g. as in
  leptoproduction or Bhabha scattering. For options, see MSTP(21).
- Processes 123 and 124, f f' -> f f' H and f f' -> f" f"' H, have
  beend added. These provide the exact matrix elements for what is
  simulated in the effective W approximation in processes 5 and 8.
- Processes 35, f gamma -> f Z, 36, f gamma -> f' W, 69,
  gamma gamma -> W+ W-, and 70, gamma W -> Z0 W, have been implemented.
  However, only in the first are Z decay angles correctly included,
  while the other decays currently are isotropic.
- Processes 121, g g -> Q Q~ H, and 122, q q~ -> Q Q~ H, have been
  implemented, using the formulae of [Kun84].
- Simulation of the two further neutral Higgs particles, H'0 and A0,
  of the two Higgs doublet scenario. This includes a total of 18 new
  processes, 151 - 153, 156 - 158, 171 - 174, 176 - 179, 181 - 182,
  and 186 - 187, obtained by analogy with the corresponding processes
  for the standard model Higgs. Couplings can either be set freely,
  see PARU(161) - PARU(185), or in accordance with the minimal
  supersymmetric scenario [GHK90], see MSTP(4). This work is still
  in progress: Born level couplings should work fine, while loop
  contributions still are incomplete.
- Possibility of anomalous W magnetic moment in process 20.
- Possibility to simulate contact interaction (quark and lepton
  substructure) in processes 11, 12, 165 and 166.
- Production of excited d* and u* quarks, by q + g -> q*, in
  processes 147 and 148.
- Possibility to have incoming hyperon beams: Lambda, Sigma (-,0,+),
  Xi (-,0), and Omega-. Since no structure functions are available,
  an average valence quark distribution is constructed from the
  proton structure function set used.
______________________________________________________________________

1.6.  Major Changes from PYTHIA 5.5 to 5.6

The updates from version 5.5 to 5.6 are all minor, and just about any
program that ran with version 5.5 will also work with PYTHIA 5.6.
The following changes might give compatibility problems from programs
based on PYTHIA 5.5:
- Structure functions have been reorganized, such that the user
  interface PYSTFE has been replaced by in-code interfaces to the
  PDFLIB [Plo91], PAKPDF [Cha91] and PHOPDF [Cha91a] structure
  function packages. To enable the appropriate interfaces, look for
  lines with C! in the first two columns of the FORTRAN source file.
  Also, the MSTP(51) - MSTP(58) switches for structure functions have
  been changed.
- processes 86 - 89 have been included for g + g -> X + g, where
  X is J/Psi, chi_0c, chi_1c or chi_2c, respectively. PARP(38) and
  PARP(39) are introduced for wave function values at the origin.

In addition, the following changes have been made, where no
compatibility problems are involved:
- Addition of a routine PYEVWT to allow users to modify cross-sections
  and weight events; see section 2.1.
- QCD shower evolution for photon inside electron. Intrinsic k_T
  in photon; see MSTP(93).
- Additional scale choices are possible for Deep Inelastic Scattering
  processes; see MSTP(22) and MSTP(32).
- It is possible to have x and Q2 of the scattered electron agree with
  the originally generated x and Q2 for Deep Inelastic Scattering
  processes, also when including parton showers etc.; see MSTP(23).

______________________________________________________________________

1.7.  Installation of program

The program is written entirely in Fortran 77, and should run on any
machine with such a compiler. The only external program required is
JETSET 7.3. Comments on compiler optimization level, random number
generators and machine precision problems are the same as given in the
JETSET 7.3 manual. (Those wishing to use the IBM AUTODBL complier option
should note that a dummy variable MSEDUM should be inserted after MSEL
in commonblock PYSUBS.)

SAVE statements have been included in accordance with the Fortran
standard. Since most ordinary machines take SAVE for granted, this
part is not particularly well tried out, however. Users on machines
without automatic SAVE are therefore warned to be on the lookout for
any variables which may have been missed.

A test program, PYTEST, is included in the PYTHIA package. It is
disguised as a subroutine, so that the user has to run a main program
      CALL PYTEST(0)
      END
This program will generate some events of different types. If PYTHIA
has not been properly installed, this program is likely to crash, or
at least generate erroneous events. If everything works properly, as
far as the program can judge, a final message to that effect is
produced, else various error messages may appear. If PYTEST(1) is
called instead of PYTEST(0), the same program is run through, but with
more complete initialization and cross-section information, and with
a listing of one event of each type generated. Finally, PYTEST(2) will
in addition give an extensive listing of the initial search for
cross-section coefficients and maxima; normally there is no reason to
use this latter option.
______________________________________________________________________

1.8.  Programming Philosophy

The Monte Carlo program is built as a slave system, i.e. the user
supplies the main program, and from this the various subroutines are
called on to execute specific tasks, after which control is returned to
the main program. Some of these tasks may be very trivial, whereas the
"high-level" routines by themselves may make a large number of
subroutine calls.

It should be noted that, while the physics content is obviously at
the center of attention, the Lund Monte Carlo also contains a more
extensive setup of auxiliary service routines than any other physics
event generator. The hope is that this will provide a comfortable
working environment, where not only events are generated, but users
also linger on to perform a lot of the subsequent studies. (As for
the relatively small attention given to physics in this manual,
the reason is that the physics is documented separately in a series
of papers, but the program pieces only here.)

The general rule is that all routines have names six characters long,
beginning with PY. Also commonblocks have names starting with PY.

Before events can be generated, an initialization call is necessary,
unlike the case of JETSET. Default values and other data are stored
in BLOCK DATA PYDATA. Thus this subprogram, as well as BLOCK DATA LUDATA
in JETSET, must be linked, which does not occur automatically with all
loaders.

Apart from writing initialization information, printing error messages
if need be, and responding to explicit requests for listings, all tasks
of the program are performed silently. All output is directed to unit
MSTU(11), by default 6, and it is up to the user to see to it that this
unit is open for write.

The Lund Monte Carlo is extremely versatile, but the price to be paid
for this is a large number of adjustable parameters and switches for
alternative modes of operation. No single user is ever likely to have
need for more than a fraction of the options available.
Since all these parameters and switches are assigned sensible default
values, there is no reason to worry about them until the need arises.

A number of error checks are performed during execution. Serious errors,
in particular those that may be found already at initialization time,
lead to the program being aborted. Less serious errors may lead to the
treatment of a particular event being cut short; see the JETSET manual.
Consequences are unpredictable when using integer valued input variables
with values not defined, or real-valued variables outside the physically
sensible range. Users beware!

**********************************************************************

2.    The Program Components

It is useful to distinguish three phases in a normal run with PYTHIA.
In the first phase, the initialization, the general character of the
run is determined. At a minimum, this requires the specification of the
incoming hadrons and the energies involved. At the discretion of the
user, it is also possible to select specific final states, and to make
a number of decisions about details in the subsequent generation.
This step is finished by a PYINIT call, at which time several variables
are initialized in accordance with the values set. The second phase
consists of the main loop over the number of events, with each new
event being generated by a PYEVNT call. This event may then be analyzed,
using information stored in some commonblocks, and the statistics
accumulated. In the final phase, results are presented. This may often
be done without the invocation of any PYTHIA routines. From PYSTAT,
however, it is possible to obtain a useful list of cross-sections for
the different subprocesses.

In the presentation in this section, the ordering above is not
strictly adhered to. In particular, information less important for
an efficient use of PYTHIA has been delegated closer to the end.
- In subsection 2.1 the subroutines to be called by the user are
introduced, particularly PYINIT, PYEVNT and PYSTAT.
- The ISUB classification code for subprocesses included in PYTHIA
is tabulated in 2.2. Some comments on the physics of these
subprocesses, and what special possibilities are open, are given
in 2.3.
- The following two subsections, 2.4 and 2.5, deal with variables
that can (with a few exceptions) be set only in the initialization
phase, i.e. before the PYINIT call, if the default values should not
be desirable. In 2.4 is collected the switches for the type of
process to generate and kinematical constraints. Subsection 2.5
covers all the switches and parameters which regulate the details
of the generation, such as choice of structure functions or Q^2
scale, on/off switches for parton showers or fragmentation, etc.
- The information available on each new event is covered in 2.6 and
2.7, in the former the general type of process and some kinematical
variables, in the latter the detailed list of all particles
generated.
- In 2.8 the further subroutines, functions and commonblocks in
PYTHIA are listed, with brief information about their purpose.
- The fragmentation routines of JETSET 7.3 are amply documented
elsewhere [Sjo90], but in 2.9 a brief reminder is given, with
special emphasis on aspects of relevance when running PYTHIA.
- A few comments on cross-section calculations in the program are
collected in 2.10.
- Finally, subsection 2.11 contains some examples on how to run the
program.
_______________________________________________________________________

2.1.  The Main Subroutines

There are two routines that users must know: PYINIT for initialization
and PYEVNT for the subsequent generation of each new event. In addition,
the cross-section and other kinds of information available with PYSTAT
is frequently useful. The three other routines described here, PYFRAM,
PYKCUT, and PYEVWT, are of more specialized interest.

SUBROUTINE PYINIT(FRAME,BEAM,TARGET,WIN)
Purpose: to initialize the generation procedure.
FRAME : a character variable used to specify the frame of the
    experiment. Uppercase and lowercase letters may be freely mixed.
    = 'CMS' : colliding beam experiment in CM frame, with beam momentum
        in +z direction and target momentum in -z direction.
    = 'FIXT' : fixed target experiment, with beam particle momentum
        pointing in +z direction.
    = 'USER' : full freedom to specify frame by giving beam momentum in
        P(1,1), P(1,2) and P(1,3) and target momentum in P(2,1), P(2,2)
        and P(2,3) in common block LUJETS.
    = 'NONE' : there will be no initialization of any processes, but
        only of resonance widths and a few other process-independent
        variables. Subsequent to such a call, PYEVNT can not be used
        to generate events, so this option is mainly intended for those
        who will want to construct their own events afterwards, but
        still want to have access to some of the PYTHIA facilities.
        In this option, the BEAM, TARGET and WIN arguments are dummy.
BEAM, TARGET : character variables to specify beam and target particles.
    Uppercase and lowercase letters may be freely mixed. An
    antiparticle may be denoted either by '~' or 'bar' at the end of
    the name. It is also possible to leave out the underscore ('_')
    directly after 'nu' in neutrino names, and the charge for
    proton and neutron.
    = 'p+' : proton.
    = 'p~-' : antiproton.
    = 'n0' : neutron.
    = 'n~0' : antineutron.
    = 'pi+' : positive pion.
    = 'pi-' : negative pion.
    = 'e-' : electron.
    = 'e+' : positron.
    = 'nu_e' : electron neutrino.
    = 'nu_e~' : electron antineutrino.
    = 'mu-' : muon.
    = 'mu+' : antimuon.
    = 'nu_mu' : muon neutrino.
    = 'nu_mu~' : muon antineutrino.
    = 'gamma' : photon (real, i.e. on mass-shell).
    = 'Lambda0' : Lambda baryon.
    = 'Sigma-' : Sigma- baryon.
    = 'Sigma0' : Sigma0 baryon.
    = 'Sigma+' : Sigma+ baryon.
    = 'Xi-' : Xi- baryon.
    = 'Xi0' : Xi0 baryon.
    = 'Omega-' : Omega- baryon.
WIN : related to energy of system, exact meaning depends on FRAME.
    FRAME='CMS' : total energy of system (in GeV).
    FRAME='FIXT' : momentum of beam particle (in GeV/c).
    FRAME='USER' : dummy (information is taken from the P vectors,
        see above).

SUBROUTINE PYEVNT
Purpose: to generate one event of the type specified by the PYINIT call.
    (This is the main routine, which calls a number of other routines
    for specific tasks.)
Note: Previously this routine was called PYTHIA; it has now been renamed
    PYEVNT to avoid confusion between a subroutine and the program-as-a-
    whole.

SUBROUTINE PYSTAT(MSTAT)
Purpose: to print out cross-sections statistics, decay widths,
    branching ratios, status codes and parameter values. PYSTAT may be
    called at any time, after the PYINIT call, e.g. at the end of the
    run, or not at all.
MSTAT : specification of information desired.
    = 1 : prints a table of how many events of the different kinds that
        have been generated and the corresponding cross-sections. All
        numbers already include the effects of cuts required by the user
        in PYKCUT.
    = 2 : prints a table of the resonances defined in the program, with
        their particle codes (KF), and all allowed decay channels.
        (If the number of generations in MSTP(1) is 3, however, channels
        involving fourth generation particles are not displayed.)
        For each decay channel is shown the sequential channel number
        (IDC) of the JETSET decay tables, the partial decay width,
        branching ratio and effective branching ratio (in the event
        some channels have been excluded by the user).
    = 3 : prints a table with the allowed hard interaction flavours
        KFIN(I,J) for beam and target particles.
    = 4 : prints a table of the kinematical cuts CKIN(I) set by the user
        in the current run.
    = 5 : prints a table with all the values of the status codes
        MSTP(I) and the parameters PARP(I) used in the current run.

SUBROUTINE PYFRAM(IFRAME)
Purpose: to transform event between different frames, if so desired.
IFRAME : specification of frame the event is to be boosted to.
    = 1 : frame specified by user in the PYINIT call.
    = 2 : CM frame of incoming particles.

SUBROUTINE PYKCUT(MCUT)
Purpose: to enable a user to reject a given set of kinematic variables
    at an early stage of the generation procedure (before evaluation
    cross-sections), so as not to spend unnecessary time on the
    generation of events that are not wanted.
    The routine will not be called unless the user requires is by
    setting MSTP(141) = 1, and never if "minimum bias" type events
    (class d of section 2.2) are to be generated as well.
    A dummy routine PYKCUT is included in the program file, so as to
    avoid unresolved external references when the routine is not used.
MCUT : flag to signal effect of user-defined cuts.
    = 0 : event is to be retained and generated in full.
    = 1 : event is to be rejected and a new one generated.
Remark : at the time of selection, several variables in the MINT and
   VINT arrays in the PYINT1 commonblock (see section 2.8) contain
   information that can be used to make the decision. The routine
   provided in the program file explicitly reads the variables that
   have been defined at the time PYKCUT is called, and also calculates
   some derived quantities. The list of information given includes
   subprocess type ISUB, E_CM, s-hat, t-hat, u-hat, p_T-hat,
   x_1, x_2, x_F, tau, y*, tau', cos(theta-hat), and a few more.
   Some of these may not be relevant for the process under study,
   and are then set to zero.

SUBROUTINE PYEVWT(WTXS)
Purpose: to allow a user to reweight event cross-sections, by process
    type and kinematics of the hard scattering. There exists two
    separate modes of usage, described in the following.
    For MSTP(142)=1, it is assumed that the cross-section of the
    process is correctly given by default in PYTHIA, but that one
    wishes to generate events biased to a specific region of phase
    space. While the WTXS factor therefore multiplies the naive
    cross-section in the choice of subprocess type and kinematics,
    the produced event comes with a compensating weight
    PARI(10)=1./WTXS, which should be used when filling histograms
    etc. In the PYSTAT(1) table, the cross-sections are unchanged
    (up to statistical errors) compared to the standard cross-sections,
    but the relative composition of events may be changed and need
    no longer be in proportion to relative cross-sections. A typical
    example of this usage is if one wishes to enhance the production
    of high-p_T events; then a weight like WTXS=(p_T/p_T_0)**2
    (with p_T_0 some fixed number) might be appropriate.
    For MSTP(142)=2, on the other hand, it is assumed that the true
    cross-section is really to be modifed by the multiplicative factor
    WTXS. The events generated therefore come with unit weight, just
    as usual. This option is really equivalent to replacing the basic
    cross-sections coded in PYTHIA, but allows more flexibility: no
    need to recompile the whole of PYTHIA.
    The routine will not be called unless MSTP(142) >= 1, and never if
    "minimum bias" type events (class d of section 2.2) are to be
    generated as well. Further, cross-sections for additional multiple
    interactions or pileup events are never affected.
WTXS: multiplication factor to ordinary event cross-section; to be set
    (by user) in PYEVWT call.
Remark : at the time of selection, several variables in the MINT and
   VINT arrays in the PYINT1 commonblock (see section 2.8) contain
   information that can be used to make the decision. The routine
   provided in the program file explicitly reads the variables that
   have been defined at the time PYEVWT is called, and also calculates
   some derived quantities. The list of information given includes
   subprocess type ISUB, E_CM, s-hat, t-hat, u-hat, p_T-hat,
   x_1, x_2, x_F, tau, y*, tau', cos(theta-hat), and a few more.
   Some of these may not be relevant for the process under study,
   and are then set to zero.
WARNING: the weights only apply to the hard scattering subprocesses.
   There is no way to reweight the shape of initial and final state
   showers, fragmentation, or other aspects of the event.
______________________________________________________________________

2.2.  The Physics Processes

A wide selection of fundamental 2 -> 1 and 2 -> 2 tree processes of the
standard model (electroweak and strong) has been included in PYTHIA,
and slots are provided for those not yet implemented. In addition,
a few "minimum bias" type processes (like elastic scattering), loop
graphs, box graphs, 2 -> 3 tree graphs and some non-standard model
processes are included. The classification is not always unique. A
process that proceeds only via an s-channel state is classified as
a 2 -> 1 process (e.g. q q~ -> gamma* -> e+ e-), but a 2 -> 2
cross-section may well have contributions from s-channel diagrams
(g g -> g g obtains contributions from g g -> g* -> g g). Also, in
the program, 2 -> 1 and 2 -> 2 graphs may sometimes be folded with
two 1 -> 2 splittings to form effective 2 -> 3 or 2 -> 4 processes
(W+ W- -> H0 is folded with q -> q' W to give q q' -> q" q'" H0).

It is possible to select a combination of subprocesses to simulate,
and also afterwards to know which subprocess was actually selected in
each event. For this purpose, all subprocesses are numbered according
to an ISUB code. The list of possible codes is given in this section,
while its detailed use will be made clear in sections 2.4 and 2.6.
Only processes marked with a + sign in the first column have been
implemented in the program to date. The others are still given here,
so that the user will easier understand how the classification works.

In the following f_i represents a fundamental fermion of flavour i,
i.e. either of d, u, s, c, b, t, l, h, e-, nu_e, mu-, nu_mu, tau-,
nu_tau, chi- or nu_chi. A corresponding antifermion is denoted f_i~.
In several cases, some classes of fermions are explicitly excluded,
since they do not couple to the g or gamma (no e+ e- -> g g, e.g.).
When processes have only been included for quarks, while leptons
might also have been possible, the notation q_i is used. In processes
where fermion masses are explicitly included in the matrix elements,
an F is used to denote an arbitrary fermion and a Q a quark.
Flavours appearing already in the initial state are denoted i and j,
whereas new flavours in the final state are denoted k and l.

Charge conjugate channels are always assumed included as well (where
separate), and processes involving a W+ also imply those involving
a W-. Wherever Z0 is written, it is understood that gamma* and
gamma*/Z0 interference should be included as well (with possibilities
to switch off either, if so desired). In practice, the full gamma*/Z0
structure is only included in subprocesses 1 and 22; for the other
processes currently a Z0 does not contain the gamma* piece.
Correspondingly, Z'0 denotes the complete set gamma*/Z0/Z'0 (or some
subset of it). Thus the notation gamma is only used for an
on-mass-shell photon.

In the last column below, references are given to works from which
formulae have been taken. Sometimes these references are to the
original work on the subject, sometimes only to the place where the
formulae are given in the most convenient or accessible form. In
several instances, errata have been obtained from the authors. Often
the formulae given in the literature have been generalized to include
trivial radiative corrections etc.

Subprocess                                    References

a) 2 -> 1, tree
+      1     f_i f_i~  ->  Z0                 EHL84
+      2     f_i f_j~  ->  W+                 EHL84
+      3     f_i f_i~  ->  H0                 EHL84
       4     gamma W+  ->  W+
+      5     Z0 Z0  ->  H0                    EHL84,Cha85
       6     Z0 W+  ->  W+
       7     W+ W-  ->  Z0
+      8     W+ W-  ->  H0                    EHL84,Cha85

b) 2 -> 2, tree
+     10     f_i f_j(~)  ->  f_i f_j(~) (QFD) Ing87
+     11     f_i f_j(~)  ->  f_i f_j(~) (QCD) Cut78,Ben84,Lan91
+     12     f_i f_i~  ->  f_k f_k~           Cut78,Ben84,Lan91
+     13     f_i f_i~  ->  g g                Cut78,Ben84
+     14     f_i f_i~  ->  g gamma            Hal78,Ben84
+     15     f_i f_i~  ->  g Z0               EHL84
+     16     f_i f_j~  ->  g W+               EHL84
      17     f_i f_i~  ->  g H0
+     18     f_i f_i~  ->  gamma gamma        Ber84
+     19     f_i f_i~  ->  gamma Z0           EHL84
+     20     f_i f_j~  ->  gamma W+           EHL84,Sam90
      21     f_i f_i~  ->  gamma H0
+     22     f_i f_i~  ->  Z0 Z0              EHL84,Gun86
+     23     f_i f_j~  ->  Z0 W+              EHL84,Gun86
+     24     f_i f_i~  ->  Z0 H0              EHL84
+     25     f_i f_i~  ->  W+ W-              EHL84,Gun86
+     26     f_i f_j~  ->  W+ H0              EHL84
      27     f_i f_i~  ->  H0 H0
+     28     f_i g  ->  f_i g                 Cut78,Ben84
+     29     f_i g  ->  f_i gamma             Hal78,Ben84
+     30     f_i g  ->  f_i Z0                EHL84
+     31     f_i g  ->  f_k W+                EHL84
      32     f_i g  ->  f_i H0
+     33     f_i gamma  ->  f_i g             Duk82
+     34     f_i gamma  ->  f_i gamma         Duk82
+     35     f_i gamma  ->  f_i Z0            Gab86
+     36     f_i gamma  ->  f_k W+            Gab86
      37     f_i gamma  ->  f_i H0
      38     f_i Z0  ->  f_i g
      39     f_i Z0  ->  f_i gamma
      40     f_i Z0  ->  f_i Z0
      41     f_i Z0  ->  f_k W+
      42     f_i Z0  ->  f_i H0
      43     f_i W+  ->  f_k g
      44     f_i W+  ->  f_k gamma
      45     f_i W+  ->  f_k Z0
      46     f_i W+  ->  f_k W+
      47     f_i W+  ->  f_k H0
      48     f_i H0  ->  f_i g
      49     f_i H0  ->  f_i gamma
      50     f_i H0  ->  f_i Z0
      51     f_i H0  ->  f_k W+
      52     f_i H0  ->  f_i H0
+     53     g g  ->  f_k f_k~                 Cut78,Ben84
+     54     g gamma  ->  f_k f_k~             Duk82
      55     g Z0  ->  f_k f_k~
      56     g W+  ->  f_k f_l~
      57     g H0  ->  f_k f_k~
+     58     gamma gamma  ->  f_k f_k~         Bar90
      59     gamma Z0  ->  f_k f_k~
      60     gamma W+  ->  f_k f_l~
      61     gamma H0  ->  f_k f_k~
      62     Z0 Z0  ->  f_k f_k~
      63     Z0 W+  ->  f_k f_l~
      64     Z0 H0  ->  f_k f_k~
      65     W+ W-  ->  f_k f_k~
      66     W+ H0  ->  f_k f_l~
      67     H0 H0  ->  f_k f_k~
+     68     g g  ->  g g                     Cut78,Ben84
+     69     gamma gamma  ->  W+ W-           Kat83
+     70     gamma W+  ->  Z0 W+              Kun87
+     71     Z0 Z0  ->  Z0 Z0 (longitudinal)
+     72     Z0 Z0  ->  W+ W-       "
+     73     Z0 W+  ->  Z0 W+       "         Ter90
      74     Z0 H0  ->  Z0 H0
      75     W+ W-  ->  gamma gamma
+     76     W+ W-  ->  Z0 Z0       "         BLS87
+     77     W+ W+(-)  ->  W+ W+(-) "         Dun86
      78     W+ H0  ->  W+ H0
      79     H0 H0  ->  H0 H0

c) 2 -> 2, tree, massive final quarks
+     81     f_i f_i~  ->  Q_i Q_i~           Com77
+     82     g g  ->  Q_i Q_i~                Com77
+     83     q_i f_j  ->  Q_k f_l             Zer90
+     84     g gamma  ->  Q_i Q_i~            Fon81
+     85     gamma gamma  ->  F_i F_i~        Bar90
+     86     g g -> J/Psi g                   Bai83
+     87     g g -> chi_0c g                  Gas87
+     88     g g -> chi_1c g                  Gas87
+     89     g g -> chi_2c g                  Gas87

d) "minimum bias"
+     91     elastic scattering               Blo85
+     92     single diffraction               Gou83,Blo85
+     93     double diffraction               Gou83,Blo85
      94     central diffraction
+     95     low-pT production                Sjo87

e) 2 -> 1, loop
     101     g g  ->  Z0
+    102     g g  ->  H0                      EHL84
+    103     gamma gamma  ->  H0              Dre89

f) 2 -> 2, box
+    111     f_i f_i~  ->  g H0               Ell88
+    112     f_i g  ->  f_i H0                Ell88
+    113     g g  ->  g H0                    Ell88
+    114     g g  ->  gamma gamma             Con71,Ber84,Dic88
+    115     g g  ->  g gamma                 Con71,Ber84,Dic88
     116     g g  ->  gamma Z0
     117     g g  ->  Z0 Z0
     118     g g  ->  W+ W-
     119     gamma gamma  ->  g g

g) 2 -> 3, tree
+    121     g g  ->  Q_k Q_k~ H0             Kun84
+    122     q + q~  ->  Q_k Q_k~ H0          Kun84
+    123     f_i f_j  ->  f_i f_j H0          Cah84
+    124     f_i f_j  ->  f_k f_l H0          Cah84
+    131     g g  ->  Z0 Q_k Q_k~             Kle90

h) non-standard model, 2 -> 1
+    141     f_i f_i~  ->  Z'0                Alt89
+    142     f_i f_j~  ->  W'+                Alt89
+    143     f_i f_j~  ->  H+
+    144     f_i f_j~  ->  R                  Ben85a
+    145     q_i l_j  ->  LQ                  Wud86
+    147     d g -> d*                        Bau90
+    148     u g -> u*                        Bau90
+    151     f_i f_i~  ->  H'0                EHL84
+    152     g g  ->  H'0                     EHL84
+    153     gamma gamma  ->  H'0             Dre89
+    156     f_i f_i~  ->  A0                 EHL84
+    157     g g  ->  A0                      EHL84
+    158     gamma gamma  ->  A0              Dre89

i) non-standard model, 2 -> 2 and 2 -> 3
+    161     f_i g  ->  f_k H+
+    162     q g  ->  l LQ                    Hew88
+    163     g g  ->  LQ LQ~                  Hew88,EHL84
+    164     q q~  ->  LQ LQ~                 Hew88
+    165     f_i f_i~ -> f_k f_k~ (gamma/Z)   Lan91
+    166     f_i f_j~ -> f_k f_l~ (W)         Lan91
+    171     f_i f_i~  ->  Z0 H'0             EHL84
+    172     f_i f_j~  ->  W+ H'0             EHL84
+    173     f_i f_j  ->  f_i f_j H'0         Cah84
+    174     f_i f_j  ->  f_k f_l H'0         Cah84
+    176     f_i f_i~  ->  Z0 A0              EHL84
+    177     f_i f_j~  ->  W+ A0              EHL84
+    178     f_i f_j  ->  f_i f_j A0          Cah84
+    179     f_i f_j  ->  f_k f_l A0          Cah84
+    181     g g  ->  Q_k Q_k~ H'0            Kun84
+    182     q q~  ->  Q_k Q_k~ H'0           Kun84
+    186     g g  ->  Q_k Q_k~ A0             Kun84
+    187     q q~  ->  Q_k Q_k~ A0            Kun84

For many of the subprocesses, further notes and qualifications may
be of interest. These are given in the following section.
Also note that some groups of subprocesses are available with the
MSEL switch, see section 2.4.

----------------------------------------------------------------------

2.3.  Comments on Physics Processes

This section contains some useful comments on the processes included
in the program, grouped by physics interest rather than sequentially
by ISUB or MSEL code. The different ISUB and MSEL codes that can be
used to simulate the different groups are given. ISUB codes within
brackets indicate the kind of processes that indirectly involve the
given physics topic, although only as part of a larger whole. Some
obvious examples, like the possibility to produce jets in just about
any process, are not spelled out in detail.

The text at times contains information on which special switches or
parameters are of particular interest to a given process. All these
switches are described in detail in the following sections, but are
alluded to here so as to provide a more complete picture of the
possibilities available for the different subprocesses. However, the
list of possibilities is certainly not exhausted by the text below.
____________________

2.3.1 QCD Jets

MSEL = 1,2
ISUB = 11,12,13,28,53,68

The basic cross-sections are taken from [Cut78]. However, a string-based
fragmentation scheme such as the Lund model needs cross-sections for the
different colour flows; these have been calculated in [Ben84] and differ
from the usual calculations by interference terms of the order 1/N^2. By
default, these interference terms are excluded; however, they can be
introduced by changing MSTP(34). In this case, the interference terms
are distributed on the various colour flows according to the pole
structure of the terms.

As an example, consider subprocess 28, q + g -> q + g. The total
cross-section for this process, obtained by summing and squaring the
Feynman s-, t-, and u-channel graphs, is [Cut78]:
   2*(1 - u*s/t^2) - 4/9*(s/u + u/s) - 1.
(The hats on s, t and u have been suppressed, and an overall factor
alpha_S^2/s^2 ignored.) Using the identity of the Mandelstam variables
for the massless case:
   s + t + u = 0,
this can be rewritten as
   (s^2 + u^2)/t^2 - 4/9*(s/u - u/s).
On the other hand, the cross-sections for the two possible colour flows
of this subprocess are [Ben84]:
   A: 4/9*(2*u^2/t^2 - u/s)
   B: 4/9*(2*s^2/t^2 - s/u),
the sum of which is:
   8/9*(s^2 + u^2)/t^2 - 4/9*(s/u - u/s).
The difference between this expression and that of [Cut78],
corresponding to the interference between the two colour flow
configurations, is then
   1/9*(s^2 + u^2)/t^2,
and can be naturally divided between colour flow A and B:
   A: 1/9*u^2/t^2
   B: 1/9*s^2/t^2.
This procedure is followed also for the other QCD subprocesses.

All the matrix elements in this group are for massless quarks (although
final state quarks are of course put on mass-shell). As a consequence,
some kind of regularization for p_T -> 0 is required. Normally the user
is expected to set the desired pTmin value in CKIN(3).

The new flavour produced in the annihilation processes (ISUB = 12,53)
is determined by the flavours allowed for gluon splitting into
quark-antiquark; see switch MDME.

For production of massive quarks, see below.
____________________

2.3.2 Heavy Flavours

MSEL=4,5,6,7,8,35,36,37,38
ISUB = 81,82,83,84,85

The cross-sections are taken from [Com77] for ISUB = 81,82 and from
[Zer90] for ISUB = 83. The former two processes are pure QCD ones,
and normally dominate. The last process proceeds via t channel W
exchange, and is mainly of interest for the production of very heavy
flavours, where the possibility of producing just one heavy quark
is kinematically favoured over pair production.

The matrix elements in this group differ from the corresponding ones
in the group above in that they correctly take into account the quark
masses. As a consequence, the cross-sections are finite for p_T -> 0.
It is therefore not necessary to introduce any special cuts.

The flavour produced is determined by the heaviest flavour allowed for
gluon splitting into quark-antiquark; see switch MDME. When one of the
MSEL options is used, MDME is automatically set by the program. Note
that only one heavy flavour at a time is allowed; if more than one is
turned on, only the heaviest will be produced (as opposed to the case
for ISUB = 12,53 above, where more than one flavour is allowed
simultaneously). Process 85 is set by the 'decay' channels of the
gamma rather than the gluon, and can therefore be set separately.
Here also lepton production is possible and, since lepton KF codes
come after quark ones, they are counted as being 'heavier', and thus
take precedence if they have been allowed, since again only the
heaviest flavour is simulated.

The lowest order processes above just represent one source of heavy
flavour production. Heavy quarks can also be present in the structure
functions at the Q^2 scale of the hard interaction, leading to
processes like Q + g -> Q + g, so-called flavour excitation, or be
created by gluon splittings g -> Q + Q~ in initial or final state
shower evolution. In fact, as the CM energy is increased, these other
processes gain in importance relative to the lowest order production
graphs above. As as example, only 10% of the b production at LHC
energies come from the lowest order graphs. The figure is even smaller
for charm, while it is at or above 50% for top. At LHC/SSC energies,
the specialized treatment described in this subsection is therefore
only of interest for top (and potential fourth generation quarks) -
the higher order corrections can here be approximated by an effective
K factor, except possibly in some rare corners of phase space. For
charm and bottom, on the other hand, it is necessary to simulate the
full event sample (within the desired kinematics cuts), and then only
keep those events with b/c either from lowest order production, or
flavour excitation, or gluon splitting. Obviously this may be a
time-consuming enterprise - although the probability for a high-p_T
event to contain (at least) one charm or bottom pair is fairly large,
most of these heavy flavours are carrying a small fraction of the total
p_T flow of the jets, and therefore do not survive normal experimental
cuts.

As an aside, it is not only for the lowest order graphs that events
may be generated with a guaranteed heavy flavour content. One may also
generate the flavour excitation process by itself, in the massless
approximation, using ISUB = 28 and setting the KFIN array appropriately.
No trick exists to force the gluon splittings without introducing
undesirable biases, however.

The cross-section for a heavy flavour pair close to threshold can be
modified according to the formulae of [Fad89], see MSTP(35). This
affects the total rate and also kinematical distributions.
____________________

2.3.3 Minimum Bias

MSEL = 1,2
ISUB = 91,92,93,95

The total and elastic cross-sections are given by the parametrizations
in [Blo85]. Several different parametrizations are available; these
can be selected with MSTP(31). For the single and double diffraction
cross-sections, the ansatz in [Gou83] is used. The remaining part of
the cross-section is automatically assigned to low-p_T events.

The simulation of the that variable in elastic or diffractive scattering
is still fairly primitive. Combined with the imprecise kinematics
treatment of small scattering angles, this means that the program should
not be used for studies of the scattered (undiffracted) hadron. In
diffractive scattering, the structure of the hadronic system selected
may be regulated with MSTP(101). No high-p_T jet production in
diffractive events is included so far.

The subprocess 95, low-p_T events, is somewhat unique, in
that no meaningful physical borderline to high-p_T events can be
defined. Even if the QCD 2 -> 2 high-p_T processes are formally
switched off, some of the events generated will be classified as
belonging to this group, with a p_T spectrum of interactions to
match the "minimum bias" event sample. Only with the option
MSTP(82) = 0 will subprocess 95 yield strictly low-p_T  events,
events which will then probably not be compatible with any
experimental event sample. A number of options exists for the detailed
structure of low-p_T events, see in particular MSTP(81) and MSTP(82).
Details of the model(s) used for these events may be found in
[Sjo87].
____________________

2.3.4 Prompt Photons

MSEL = 10
ISUB = 14,18,29,114,115

Processes ISUB = 14,29 give the main source of single gamma production,
with ISUB = 115 giving an additional contribution which in some
kinematics regions may become important. For gamma pair production,
the process ISUB = 18 is often overshadowed in importance by ISUB = 114.
Cross-sections for the Born term graphs 14, 18 and 29 are found e.g. in
[EHL84], while the box graphs 114 and 115 are given in [Con71,Ber84,
Dic88]. For the issue of normalization of single photon vs. diphoton,
see section 2.10.

Another source of photons is bremsstrahlung off incoming or outgoing
quarks. This has to be treated on an equal footing with QCD parton
showering. For timelike parton shower evolution, i.e. in the final
state showering and in the side branches of the initial state showering,
photon emission may be switched on with MSTJ(41) (see JETSET manual).
Photon radiation off the spacelike incoming quark legs is not yet
included, but should be of lesser importance for production at
reasonably large p_T values. Radiation off an incoming electron is
included in a leading log approximation.

WARNING: The cross-sections for the box graphs 114 and 115 become
    very complicated, numerically unstable and slow when the
    full quark mass dependence is included. For quark masses much
    below the shat scale, the simplified massless expressions are
    therefore used - a fairly accurate approximation. However, there
    is another set of subtle numerical cancellations between different
    terms in the massive matrix elements in the region of small-angle
    scattering. The associated problems have not been sorted out yet.
    There are therefore two possible solutions. One is to use the
    massless formulae throughout. The program then becomes faster and
    numerically stable, but does e.g. not give the characteristic dip
    (due to destructive interference) at top threshold. This is the
    current default procedure, with five flavours assumed, but this
    number can be changed in MSTP(38). The other possibility is to
    impose cuts on the scattering angle of the hard process, see
    CKIN(27) and CKIN(28), since the numerically unstable regions
    are when abs(cos(theta-hat)) is close to unity. It is then also
    necessary to change MSTP(38) to 0.
____________________

2.3.5 Single Z/W Production

MSEL = 11,12,13,14,(21)
ISUB = 1,2,15,16,30,31,35,36,131,(141)

This group consists of 2 -> 1 processes, single resonance production,
and 2 -> 2 processes, resonance + jets. With initial state parton
showers turned on, the 2 -> 1 processes also generate resonance + jets;
in order to avoid double counting, the corresponding 2 -> 2 processes
should not be turned on simultaneously. The basic rule is to use the
2 -> 1 processes for inclusive generation of Z/W, whereas the 2 -> 2
processes should be used for study of the subsample with high
transverse momentum.

The Z0 of subprocess 1 includes the full interference structure
Z0/gamma*; via MSTP(43) the user can select to produce only gamma*,
only Z0, or the full Z0/gamma*. The same holds true for the Z'0 of
subprocess 141; via MSTP(44) any combination of gamma*, Z0 and Z'0
can be selected. Thus, subprocess 141 with MSTP(44) = 4 is essentially
equivalent to subprocess 1 with MSTP(43) = 3; however, the former
also includes the possibility of a decay into Higgses. The Z0 of
other processes than 1 and 141 are currently pure Z0 only.

For the 2 -> 1 processes, the Breit-Wigner includes an shat-dependent
width, which should provide an improved description of line shapes.

The process e+e- -> e+e- + Z0 can be simulated in two different ways.
One is to make use of e 'sea' quark structure function inside e, i.e.
have splittings e -> gamma -> e (currently in the collinear limit; and
with an e+e- system in the forward direction represented by a gamma).
This can be obtained, together with ordinary Z production, by using
subprocess 1, with MSTP(11)=1 and MSTP(12)=1. Then the contribution of
the type above is 5.0 pb for a 500 GeV e+e- collider, compared to the
correct 6.2 pb [Hag90]. Alternatively one may use process 35, with
MSTP(11)=1 and MSTP(12)=0. To catch the singularity in the forward
direction, regularized by the electron mass, it is necessary to
set CKIN(3)=CKIN(5)=0.01 - using lower values will only slow down
execution, not significantly increase cross-section. One then obtains
5.1 pb, i.e. 20% below the correct, and also generates a pT distribution
for the Z0; this is therefore to be preferred.

Process 36, f + gamma -> f' + W may have corresponding problems; except
that in e+e- the forward radiation amplitude for e + gamma -> nu W
is killed (radiation zero). It is therefore feasible to use the
default CKIN(3) and CKIN(5) values in e+e-, and one also comes closer
to the correct cross-section.

One single true 2 -> 3 process is included in this class as well;
namely g g -> Z q q~, with full massive matrix elements. The more
complicated phase space and the lengthy matrix element evaluations
makes this process extremely slow. With q = b, it may form an
important background to intermediate mass Higgs searches in the
multilepton channel. The quark flavour is stored in KFPR(131,2);
the default is 5 = b. The kinematics is set up in terms of a Z
recoiling against the q q~ system, and all ordinary kinematics cut
for a 2 -> 2 process can be used on this level, including CKIN(43)
and CKIN(44) to restrict the range of the q q~ invariant mass.
In addition, for this process alone, CKIN(51) - CKIN(54) can be used
to set the p_T range of the two quarks; as is to be expected that of
the Z is set by CKIN(3) - CKIN(4). Since the optimization procedure
is not set up to probe the full multidimensional phase space allowed
in this process, maximum violations may be quite large. It may then
be useful to make a preliminary run to find how big the violations
are in total, and then use the MSTP(121) = 1 option in the full run.
____________________

2.3.6 Neutral Higgs Production

MSEL = 16
ISUB = 3,5,8,24,26,(71,72,73,76,77),102,103,111,112,113,121,122,
       123,124

This section deals with the production of a standard model Higgs
exclusively, as obtained by default for the processes above.
The issue of additional Higgses is addressed in section 2.3.9.
Many different processes can be involved, as seen from the list
above. The proper choice depends on the actual Higgs mass, and
(occasionally) on the desired search strategy.

For a Higgs which can still be handled not unreasonably well in a
narrow width approximation, i.e. with a mass below 700 GeV, say,
the production processes that are involved are ISUB = 3,102,103,123,
124, as obtained for MSEL = 16.

The subprocess t t~ -> H0 (a subset of the more general subprocess
3, but the only subset of importance for heavy Higgs production at
hadronic colliders) is by now known to overestimate the cross-section
for heavy Higgs production as compared to a more careful calculation
based on the subprocess g g -> t t~ H0 (121). This process, as well
as the less important q + q~ -> t t~ H0 (122), are now included,
but the complicated multibody phase space makes it run slower than
most other processes. One should therefore think twice before using
it. As usual, it would be doublecounting to include both 3 and 121.

The complications related to g g -> t t~ H are actually not
major, since heavy Higgs production is anyway dominated by
subprocesses 102,123 and 124. The two subprocesses 5 and 8 are
effective W approximations to the matrix elements now correctly
included in 123 and 124. For light Higgs masses, the effective W
approximation is not good (it overestimates cross-sections), and
subprocesses 5 and 8 should not be used there. For a heavier Higgs,
agreement is better, and either approach may be chosen.

The subprocesses 5 and 8, V V -> H0, which contribute to the
processes V V -> V' V' (V and V' intermediate vector bosons) show
bad high energy behaviour, which can be cured only by the inclusion
of all V V -> V' V' graphs, as is done in subprocesses 71, 72, 73,
76 and 77. In particular, subprocesses 5 and 8 give rise to a
fictitious high-mass tail of the Higgs. If this tail is thrown away,
however, the agreement between the s-channel graphs (subprocesses
5 and 8) and the full set of graphs (subprocesses 71 etc.) is very
good: for a Higgs of nominal mass 300 (800) GeV, a cut at 600 (1200)
GeV retains 95% (84%) of the total cross-section, and differs from
the exact calculation, cut at the same values, by only 2% (11%)
(numbers for SSC energies). With this prescription there is
therefore no need to use subprocesses 71 etc. rather than
subprocesses 5 and 8.

For subprocesses 71, 72, 76 and 77, an option is included (see
MSTP(46)) whereby the user can select only the s-channel Higgs
graph; this will then be essentially equivalent to running
subprocess 5 or 8 with the proper decay channels (i.e. Z0Z0 or
W+W-) set via MDME. The difference is that the Breit-Wigner in
subprocesses 5 and 8 contain an shat-dependent width, whereas the
width in subprocesses 71 - 77 is calculated at nominal Higgs mass;
also, higher order corrections to the widths are treated more
accurately in subprocesses 5 and 8. Further, processes 71 - 77 assume
on-mass-shell incoming W/Z, with associated kinematics factors,
while processes 5 and 8 have W/Z correctly spacelike. All this
leads to differences in the cross-sections by up to a factor 1.5.
For subprocesses 71 - 77, also read comments in next subsection.

A subprocess like 113, with a Higgs recoiling against a gluon jet, is
also effectively generated by initial state corrections to subprocess
102; thus, in order to avoid double counting, just as for the case
of Z0/W+ above, these subprocesses should not be switched on
simultaneously, but 3, 5, 8, and 102 be used for inclusive production
of Higgs, and 111 - 113 for the study of the Higgs subsample with high
transverse momentum.

The Higgs can also be produced in association with a Z0/W,
ISUB = 24,26. In pp collisions, these processes have a lower
cross-section than single Higgs production, but are still of
interest at the lower Higgs mass range because of a potentially
better signal to background ratio. Of course, the process
Z0 -> Z0 + H0 is particularly interesting for e+e- colliders. The
other processes possible here are 103 and, in particular, 123 and 124.
As above, we note that 123 and 124 give much more accurate
cross-sections than the equivalent processes 5 and 8 based on the
effective W approximation.

The branching ratios of the Higgs are very strongly dependent on the
actual mass.In principle, the program is set up to calculate these
correctly at initialization. However, higher order corrections may at
times be important and not fully unambiguous; see e.g. MSTP(37).
____________________

2.3.7 gamma/Z/W Pairs

MSEL = 15
ISUB = 19,20,22,23,25,69,70,71,72,73,76,77

This heading contains two different types of processes:
fermion-antifermion annihilation into gamma/Z/W pairs, and scattering
of intermediate vector bosons. Obviously other sources of gamma/Z/W
pairs also may exist, like the Higgs or new gauge bosons.

The first set contains the standard gamma/Z/W produced by f + fbar
(f + fbar') annihilation. We remind that the notation gamma means a
massless gamma, while Z in principle stands for the full gamma*/Z0
structure. In fact, currently the Z0 appearing in almost all
subprocesses are pure Z0; only the Z0 of process 22 contains the full
interference structure, modifiable with MSTP(43). Currently subprocess
22 does not contain all possible contributions, but only those where
both gamma*/Z0 are produced off the quark line, i.e. excluding e.g.
emission of a gamma* off a Z0 decay product. This approximation is
fully valid for two on-mass-shell Z0:s, but should give the dominant
contribution also in more general situations. We also remind the user
that the mass ranges of the two daughters may be set with the
CKIN(41) - CKIN(44) parameters; this is particularly convenient e.g.
to pick one Z almost on mass-shell and the other not.

The second set is of interest in that only the inclusion of the full
set of VV -> V'V' graphs will cure the bad high energy behaviour of
VV -> H0 (see 2.3.6). There is an option (see MSTP(46)) that selects
only the s-channel Higgs exchange from subprocesses 71,72,76 and 77;
with this option, these subprocesses will be essentially equivalent
to subprocesses 5 and 8 if the Higgs is allowed to decay into Z or W
(see, however, further comments in the Higgs subsection).

For subprocess 77, there is an option (see MSTP(45)) to select the
charge combination of the scattering W's: like-sign, opposite-sign
(relevant for Higgs), or both.

In the absence of a Higgs, the sector of longitudinal Z and W
scattering will become strongly interacting at energies above 1 TeV.
The models proposed by Dobado, Herrero and Terron [DHT90] to describe
this kind of physics have been included as alternative matrix
elements for subprocesses 71,72,73,76,77, selectable by MSTP(46).
(From the point of view of the general classification scheme for
subprocesses, this kind of models should appropriately be included
as separate subprocesses with numbers above 100, but the current
solution allows a more efficient reuse of existing code.)
By a proper choice of parameters, it is also here possible to
simulate the production of a techni-rho.

WARNING: For subprocess 77, the option for like-sign W scattering
    presently gives a non-sensical cross-section and should not be
    used; the default is set for W+W- scattering. The error does not
    affect the Dobado, Herrero and Terron options of this process.
____________________

2.3.8 New Gauge Bosons

MSEL = 21,22,24
ISUB = 141,142,144

The Z' of subprocess 141 contains the full gamma*/Z0/Z'0 interference
structure; the choice between different combinations is made via
MSTP(44). The couplings of Z' to quarks and leptons can be set via
PARU(121) - PARU(128), the coupling to W via PARU(129) - PARU(130), and
to H+ via PARU(143). The coupling of the Z to H+ is set via PARU(142).
By a suitable setting of these parameters, it is possible to simulate
several different physics scenarios. The default values agree with the
'extended gauge model' of [Alt89]. Further description of the coupling
parameters are given in the JETSET 7.3 manual.

The W' of subprocess 142 so far does not contain interference with the
standard model W - in practice this should not be a major limitation.
The couplings of the W' to quarks and leptons are set via PARU(131) -
PARU(134), the coupling to Z + W via PARU(135) - PARU(136). Further
comments as for Z'; in particular, default couplings again agree with
the 'extended gauge model' of [Alt89].

The strenght of couplings Z' -> Z H and W' -> W H can be set via
PARU(145) and PARU(146), respectively. Default values are in agreement
with expectations in a left-right symmmetric model [Coc90].

The R boson (particle code 40) of subprocess 144 represents one
possible scenario for a horizontal gauge boson, i.e. a gauge boson
that couples between the generations, inducing processes like
s + dbar -> mu- + e+. Experimental limits on flavour changing neutral
currents forces such a boson to be fairly heavy. The model implemented
is the one described in [Ben85a].
____________________

2.3.9 Extended Higgs scenario

MSEL = 23
ISUB = 3,24,26,102,103,121,122,123,124,(141),143,151,152,153,156,
       157,158,161,171,172,173,174,176,177,178,179

The particle content of a two Higgs doublet scenario is included:
two neutral scalar particles, 25 and 35, one pseudoscalar one, 36,
and a charged doublet, 37. (Of course, these particles may also be
associated with corresponding Higgs states in larger multiplets).
By convention we choose to call the lighter scalar Higgs H0 and
the heavier H'0 - this differs from the convention in the minimal
supersymmetric extension of the standard model (MSSM), where the
ligher is called h0 and the heavier H0, but allows us to call the
Higgs of the one-Higgs scenario H0. The pseudoscalar is called A0
and the charged H+ and H-.

A number of H0 processes have been duplicated for H'0 and A0.
The correspondence between ISUB numbers is as follows.
                             H0     H'0    A0
f_i f_i~  ->                  3     151   156
g g  ->                     102     152   157
gamma gamma  ->             103     153   158
f_i f_i~  ->  Z +            24     171   176
f_i f_j~  ->  W+ +           26     172   177
g g -> Q_k Q_k~ +           121     181   186
q_i + q_i~ -> Q_k Q_k~ +    122     182   187
f_i f_j  ->  f_i f_j +      123     173   178
f_i f_j  ->  f_k f_l +      124     174   179
Note that several of the processes above are not expected to take
place at all, due to vanishing Born term couplings. We have still
included them for flexibility in simulating arbitrary couplings at
Born or loop level.

A few standard model Higgs processes have no correspondence in the
scheme above. These include: 5 and 8, which anyway have been
superseded by 123 and 124; 71, 72, 73, 76, 77, which deal with
what happens if there is no light Higgs, and so is a scenario
complementary to the one above, where several light Higgses are
assumed; and 111, 112 and 113, which describe the high-p_T tail
of the Higgs production, and are less interesting for most Higgs
studies.

In processes 121,122,181,182,186,187 the recoiling heavy flavour
is assumed to be top, which is the only one of interest in the
standard model, and the one where the structure function approach
invoked in processes 3,151,156 is least reliable. However, it is
possible to change the quark flavour in 121 etc.; for each process
ISUB this flavour is given by KFPR(ISUB,2).

By default, the H0 has the couplings of the standard model Higgs,
while the H'0 and A0 have couplings set in PARU(171) - PARU(177)
and PARU(181) - PARU(187). The default values for the latter have
no deep physics motivation, but are set just so the program will
not crash due to the absence of any couplings whatsoever. The user
should therefore set the couplings above to his desired values if
he wants to simulate either H'0 or A0. Also the couplings of the
H0 particle can be modified, in PARU(161) - PARU(165), provided
that MSTP(4) is set to 1.

For MSTP(4) = 2, the mass of the H0 (in PMAS(25,1)) and the
tan(beta) value (in PARU(141)) are used to derive the masses of
the other Higgses, as well as all Higgs couplings. PMAS(35,1) -
PMAS(37,1) and PARU(161) - PARU(195) are overwritten accordingly.
The relations used are the ones of the Born level MSSM [GHK90].
Note that not all combinations of m_H and tan(beta) are allowed;
the requirement of a finite A mass imposes the constraint
    m_H < m_Z * (tan^2(beta) - 1)/(tan^2(beta) + 1),
or, equivalently,
    tan^2(beta) > (m_Z + m_H)/(m_Z - m_H).
If this condition is not fulfilled, the program will crash.

The basic subprocess for charged Higgs production is ISUB = 143.
However, this process is dominated by t + b~ -> H+, and so
depends on the choice of t structure function. A better
representation is provided by subprocess 161, f + g -> f' + H+-;
i.e. actually b + g -> t + H-. It is therefore recommended to use
161 and not 143. The tan(beta) parameter, which is relevant also
for charged Higgs couplings, is set via PARU(141).

Note in particular, that subprocess 141 is one possibility for
H+ production, in addition to subprocesses 143 and 161. Note also,
that it is only via subprocess 141 that a Z can be made to decay
into charged Higgses. The coupling of the Z to H+ + H- is regulated
by PARU(142), and that of the Z' by PARU(143).

Process 141 can also be used to simulate Z0 -> H0 A0 and
Z0 -> H'0 A0 for associated neutral Higgs production. Finally,
heavier Higgses may decay into lighter ones, if kinematically
allowed, in processes like A -> Z0 H0 or H+ -> W+ H0.
____________________

2.3.10 Leptoquarks

MSEL = 25
ISUB = 145,162,163,164

Several processes that can generate a leptoquark have been included.
Currently only one leptoquark has been implemented, as particle 39.
The leptoquark is assumed to carry specific quark and lepton quantum
numbers, by default u quark + electron. These flavour numbers are
conserved, i.e. a process like u + e -> LQ -> d + nu is not allowed.
This may be a bit restrictive, but represents one of many leptoquark
possibilities.

Although only one leptoquark is implemented, its flavours may be
changed arbitrarily to study the different possibilities. The flavours
of the leptoquark are defined by the quark and lepton flavours in the
decay mode list. Since only one decay channel is allowed, this means
that the quark flavour is stored in KFDP(MDCY(39,2),1) and the lepton
one in KFDP(MDCY(39,2),2). The former must always be a quark, while
the latter could be a lepton or an antilepton; a charge conjugate
partner is automatically defined by the program. At initialization,
the charge is recalculated as a function of the flavours defined; also
the leptoquark name is redefined to be of the type 'LQ_(q)(l)',
where actual quark (q) and lepton (l) flavours are displayed.

The LQ -> q + l vertex contains an undetermined Yukawa coupling
strength, which affects both the width of the leptoquark and the
cross-section for many of the production graphs. This strenght may
be changed in PARU(151). The definition of PARU(151) corresponds to
the k factor of [Hew88], i.e. to lambda**2/(4*pi*alpha_em), where
lambda is the Yukawa coupling strenght of [Wud86]. Note that PARU(151)
thus is quadratic in the coupling.

The leptoquark is likely to be fairly long-lived, in which case it
has time to fragment into a mesonic or baryonic type states, which
would decay later on. This is a bit tedious to handle, so therefore
the leptoquark is always assumed to decay before fragmentation has
to be considered. This may give some imperfections in the event
generation, but should not be off by much in the final analysis.

Inside the program, the leptoquark is treated as a resonance, actually
the only coloured such. This requires some extra care, and in particular
it is not allowed to put the leptoquark stable, by modifying either
MDCY(39,1) or MSTP(41): then the leptoquark would be handed undecayed to
JETSET, which would try to fragment it (as it does with any other
coloured object), and most likely crash.
____________________

2.3.11 Compositeness and Anomalous Couplings

ISUB = 11,12,20,165,166

Some processes have been set up to allow anomalous coupling to be
introduced, in addition to the standard model ones. These can be
switched on by MSTP(5)>=1; the default MSTP(5)=0 corresponds to the
standard model behaviour.

In process 11 and 12, quark substructure is included in the left-left
isoscalar model [Lan91] for MSTP(5)=1; with compositeness scale Lambda
given in PARU(155) (default 1000 GeV) and sign eta of interference
term in PARU(156) (default +1; only other alternative -1). The model
above assumes that only u and d quarks are composite (at least at the
scale studied); with MSTP(5)=2 compositeness terms are included in the
interactions between all quarks.

The processes 165 and 166 are basically equivalent to 1 and 2, i.e.
gamma/Z0 and W+- exchange, respectively, but a bit less fancy (no
s-dependent width etc.). Only one final state flavour at the time is
generated; this flavour should be set in KFPR(165,1) and KFPR(166,1),
respectively. For process 166 one gives the down-type flavour, and the
program will associate the up-type flavour of the same generation.
Defaults are 11 in both cases, i.e. e+e- and e+nu_e/e-nu_e~ final
states. While MSTP(5)=0 gives the standard model results, MSTP(5)=1
contains the left-left model isoscalar model (which does not affect
process 166), and MSTP(5)=3 the helicity-nonconserving model (which
affects both) [Lan91]. Both the models above assume that only u and d
quarks are composite; with MSTP(5)=2 or 4, respectively, contact terms
are included for all quarks in the initial state. Parameters are
PARU(155) and PARU(156), as above.

Note that processes 165 and 166 are bookkept as 2 -> 2 processes,
while 1 and 2 are 2 -> 1 ones. This means that the default Q2 scale in
structure functions is p_T^2 for the former and shat for the latter.
To make contact between the two, it is recommended to set MSTP(32)=4,
so as to use shat as scale also for processes 165 and 166.

In process 20, for W + gamma pair production, it is possible to set
an anomalous magnetic moment for the W in PARU(153) (= eta = kappa-1;
where kappa = 1 is the standard model value). The production process
is affected according to the formulae of [Sam90]; W decay currently
remains unaffected. It is necessary to set MSTP(5)=1 to enable this
extension.
____________________

2.3.12 Excited fermions

ISUB = 147,148

Compositeness scenarios may also give rise to sharp resonances of
excited quarks and leptons. If MSTP(6)=1, then at initialization
the standard fourth generation of fermions will be overwritten,
and made to correspond to an excited copy of the first generation,
consisting of spin 1/2 particles d* (code 7), u* (8), e* (17) and
nu*_e (18). Since the original fourth generation information is lost,
it is then not in the same run possible to generate fourth
generation particles.

The current implementation only contains gauge interactions, i.e.
no contact interactions. The couplings f, f' and f_s to the SU(2),
U(1) and SU(3) groups are stored in PARU(157) - PARU(159), the
scale parameter Lambda in PARU(155); users are also expected to
change f* masses in accordance with what is desired - see [Bau90]
for details on conventions. This means that decay processes are of
the types q* -> q + g, q + gamma, q + Z0 or q' + W+-. Production
is currently only by quark-gluon fusion.
______________________________________________________________________

2.4.  Switches for Event Type and Kinematics Selection

By default, only QCD 2 -> 2 processes are generated by PYTHIA,
composed of hard interactions above p_T-hat_min = PARP(81),
with low-p_T processes added on so as to give the full (parametrized)
inelastic, non-diffractive cross-section. With the help of the
commonblock PYSUBS, it is possible to select for the generation of
another process, or a combination of processes. It is also allowed
to restrict the generation to specific incoming partons/particles
at the hard interaction. This often automatically also restricts
final state flavours but, in processes like resonance production
(Z0, W+, H0, Z'0, H+ or R0) or QCD/QED production of new flavours
(ISUB = 12, 53, 54, 81, 82, 83, 84, 85), switches in the JETSET
program may be used to this end; see section 2.9.

The CKIN array may be used to impose specific kinematics cuts.
The user should here be warned that, if kinematical variables are
too strongly restricted, the generation time per event may become
very long. In extreme cases, where the cuts effectively close the
full phase space, the event generation may run into an infinite
loop. The generation of 2 -> 1 resonance production is performed in
terms of the m-hat and y* variables, and so the ranges CKIN(1) -
CKIN(2) and CKIN(7) - CKIN(8) may be arbitrarily restricted without
a significant loss of speed. For 2 -> 2 processes, cos(theta-hat) is
added as a third generation variable, and so additionally the range
CKIN(27) - CKIN(28) may be restricted without any danger.

In addition to the variables found in PYSUBS, also those in the
PYPARS commonblock may be used to select exactly what one wants
to have simulated. These possibilities will be described in the
following subsection.

The notation used above and in the following is that 'hat' denotes
internal variables in the hard scattering subsystem, while '*' is
for variables in the CM frame of the event-as-a-whole. Effects from
initial and final state radiation are not included, since they are not
known at the time the kinematics at the hard interaction is selected.
The sharp kinematical cutoffs that can be imposed on the generation
process are therefore smeared, both by QCD radiation and by
fragmentation. In a study of events within a given window of
experimentally defined variables, it is up to the user to leave
such liberal margins that no events are missed. In other words, cuts
have to be chosen such that a negligible fraction of events migrate
from outside the simulated region to inside the interesting region.
Often this may lead to low efficiency in terms of what fraction of
the generated events are actually of interest to the user.

COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
Purpose: to allow the user to run the program with any desired subset
    of processes, or restrict flavours or kinematics.

MSEL (D=1) a switch to select between full user control and some
    preprogrammed alternatives.
    = 0 : desired subprocesses have to be switched on in MSUB, i.e.
        full user control.
    = 1 : depending on incoming particles, three different alternatives
        are used.
        Lepton-lepton: Z or W production (ISUB = 1 or 2).
        Lepton-hadron: deep inelastic scattering (ISUB = 10).
        Hadron-hadron: QCD high-p_T processes (ISUB = 11, 12, 13, 28,
            53, 68); additionally low-p_T production if CKIN(3) <
            PARP(81) or PARP(82), depending on MSTP(82) (ISUB = 95).
            If low-p_T is switched on, the other CKIN cuts are not
            used.
    = 2 : as MSEL = 1 for lepton-lepton and lepton-hadron. For
        hadron-hadron all QCD processes, including low-pT, single and
        double diffractive and elastic scattering, are included (ISUB =
        11, 12, 13, 28, 53, 68, 91, 92, 93, 95). The CKIN cuts are
        here not used.
    = 4 : charm (cc~) production with massive matrix elements
        (ISUB = 81, 82, 84, 85).
    = 5 : bottom (bb~) production with massive matrix elements
        (ISUB = 81, 82, 84, 85).
    = 6 : top (tt~) production with massive matrix elements
        (ISUB = 81, 82, 84, 85).
    = 7 : low (ll~) production with massive matrix elements
        (ISUB = 81, 82, 84, 85).
    = 8 : high (hh~) production with massive matrix elements
        (ISUB = 81, 82, 84, 85).
    = 10 : prompt photons (ISUB = 14, 18, 29).
    = 11 : Z0 production (ISUB = 1).
    = 12 : W+/- production (ISUB = 2).
    = 13 : Z0 + jet production (ISUB = 15, 30).
    = 14 : W+/- + jet production (ISUB = 16, 31).
    = 15 : pair production of different combinations of gamma, Z0 and
        W+/- (except gamma + gamma; see MSEL = 10) (ISUB = 19, 20,
        22, 23, 25).
    = 16 : H0 production (ISUB = 3, 102, 103, 123, 124).
    = 17 : H0 + Z0 or H0 + W+/- (ISUB = 24, 26).
    = 18 : H0 production, combination relevant for e+e- annihilation
        (ISUB = 24, 103, 123, 124).
    = 19 : H0, H'0 and A0 production, excepting pair production
        (ISUB = 24, 103, 123, 124, 153, 158, 171, 173, 174, 176,
        178, 179).
    = 21 : Z'0 production (ISUB = 141).
    = 22 : W'+/- production (ISUB = 142).
    = 23 : H+/- production (ISUB = 143).
    = 24 : R0 production (ISUB = 144).
    = 25 : LQ (leptoquark) production (ISUB = 145, 162, 163, 164).
    = 35: single bottom production by W exchange (ISUB = 83).
    = 36: single top production by W exchange (ISUB = 83).
    = 37: single low production by W exchange (ISUB = 83).
    = 38: single high production by W exchange (ISUB = 83).

MSUB : (D=200*0) array to be set when MSEL = 0 (for MSEL >= 1 relevant
    entries are set in PYINIT) to choose which subset of subprocesses to
    include in the generation. The ordering follows the ISUB code given
    in subsection 2.2 (with comments as given there).
    MSUB(ISUB) = 0 : the subprocess is excluded.
    MSUB(ISUB) = 1 : the subprocess is included.
    Note: when MSEL >= 1 the MSUB values set by the user are never
        changed by PYTHIA. If the user wants to combine several
        different 'subruns', each with its own PYINIT call, into one
        single run, it is up to him to remember not only to switch on
        the new processes before each new PYINIT call, but also to
        switch off the old ones that are no longer desired.

KFIN(I,J) : provides an option selectively to switch on and off
    contributions to the cross-sections from the different incoming
    partons/particles at the hard interaction. In combination with the
    JETSET resonance decay switches, this also allows the user to set
    restrictions on flavours appearing in the final state.
    I : is 1 for beam side of event and 2 for target side.
    J : enumerates flavours according to the KF code; see section 2.7,
        or the JETSET manual.
    KFIN(I,J) = 0 : the parton/particle is forbidden.
    KFIN(I,J) = 1 : the parton/particle is allowed.
    By default, everything is on, except for J = 0, which does not
    have a physical meaning.

CKIN : kinematics cuts that can be set by the user before the PYINIT
    call, and that affect the region of phase space within which
    events are generated. Some cuts are 'hardwired' while most are
    'softwired'. The hardwired ones are directly related to the
    kinematical variables used in the event selection procedure,
    and therefore have negligible effects on program efficiency.
    The most important of these are CKIN(1) - CKIN(8), CKIN(27) -
    CKIN(28), and CKIN(31) - CKIN(32). The softwired ones are most
    of the remaining ones, that can not fully be taken into account
    in the kinematical variable selection, so that generation in
    constrained regions of phase space may be slow. In extreme
    cases the phase space may be so small that the maximization
    procedure fails to find any allowed points at all (though some
    small region might still exist somewhere), and therefore switches
    off some subprocesses, or aborts altogether.
CKIN(1), CKIN(2) : (D=2.,-1.) range of allowed m-hat = sqrt(s-hat)
    values. IF CKIN(2) < 0., the upper limit is inactive.
CKIN(3), CKIN(4) : (D=0.,-1.) range of allowed p_T-hat values for
    hard 2 -> 2 processes, with transverse momentum p_T-hat defined
    in the rest frame of the hard interaction. If CKIN(4) < 0.,
    the upper limit is inactive. For processes which are singular
    in the limit p_T-hat -> 0 (see CKIN(6)), CKIN(5) provides
    an additional constraint. The CKIN(3) and CKIN(4) limits can
    also be used in 2 -> 1 -> 2 processes. Here, however, the product
    masses are not known and hence assumed vanishing in the event
    selection. The actual p_T range for massive products is thus
    shifted downwards with respect to the nominal one.
CKIN(5) : (D=1.) lower cutoff on p_T-hat values, in addition to the
    CKIN(3) cut above, for processes which are singular in the
    limit p_T-hat -> 0 (see CKIN(6)).
CKIN(6) : (D=1.) hard 2 -> 2 processes, which do not proceed only
    via an intermediate resonance (i.e. are 2 -> 1 -> 2 processes),
    are classified as singular in the limit p_T-hat -> 0 if either
    or both of the two final state products has a mass m < CKIN(6).
CKIN(7), CKIN(8) : (D=-10.,10.) range of allowed scattering
    subsystem rapidities y* in the CM frame of the event,
    where y* = (1/2) * ln(x_1/x_2).
CKIN(9), CKIN(10) : (D=-10.,10.) range of allowed (true) rapidities
    for the product with largest rapidity in a 2 -> 2 or a
    2 -> 1 -> 2 process, defined in the CM frame of the event,
    i.e. max(y*_3, y*_4).
CKIN(11), CKIN(12) : (D=-10.,10.) range of allowed (true) rapidities
    for the product with smallest rapidity in a 2 -> 2 or a
    2 -> 1 -> 2 process, defined in the CM frame of the event,
    i.e. min(y*_3, y*_4). Consistency thus requires
    CKIN(11) <= CKIN(9) and CKIN(12) <= CKIN(10).
CKIN(13), CKIN(14) : (D=-10.,10.) range of allowed pseudorapidities
    for the product with largest pseudorapidity in a 2 -> 2 or a
    2 -> 1 -> 2 process, defined in the CM frame of the event,
    i.e. max(eta*_3, eta*_4).
CKIN(15), CKIN(16) : (D=-10.,10.) range of allowed pseudorapidities
    for the product with smallest pseudorapidity in a 2 -> 2 or a
    2 -> 1 -> 2 process, defined in the CM frame of the event,
    i.e. min(eta*_3, eta*_4). Consistency thus requires
    CKIN(15) <= CKIN(13) and CKIN(16) <= CKIN(14).
CKIN(17), CKIN(18) : (D=-1.,1.) range of allowed cos(theta*) values
    for the product with largest cos(theta*) value in a 2 -> 2 or a
    2 -> 1 -> 2 process, defined in the CM frame of the event,
    i.e. max(cos(theta*_3),cos(theta*_4)).
CKIN(19), CKIN(20) : (D=-1.,1.) range of allowed cos(theta*) values
    for the product with smallest cos(theta*) value in a 2 -> 2 or a
    2 -> 1 -> 2 process, defined in the CM frame of the event,
    i.e. min(cos(theta*_3),cos(theta*_4)). Consistency thus requires
    CKIN(19) <= CKIN(17) and CKIN(20) <= CKIN(18).
CKIN(21), CKIN(22) : (D=0.,1.) range of allowed x_1 values for
    the parton on side 1 that enters the hard interaction.
CKIN(23), CKIN(24) : (D=0.,1.) range of allowed x_2 values for
    the parton on side 2 that enters the hard interaction.
CKIN(25), CKIN(26) : (D=-1.,1.) range of allowed Feynman-x values,
    where x_F = x_1 - x_2.
CKIN(27), CKIN(28) : (D=-1.,1.) range of allowed cos(theta-hat) values
    in a hard 2 -> 2 scattering, where theta-hat is the scattering
    angle in the rest frame of the hard interaction.
CKIN(31), CKIN(32) : (D=2.,-1.) range of allowed m'-hat values, where
    m'-hat is the mass of the complete three- or four-body final state
    in 2 -> 3 or 2 -> 4 processes (while m-hat, constrained in CKIN(1)
    and CKIN(2), here corresponds to the one- or two-body central
    system). If CKIN(32) < 0., the upper limit is inactive.
CKIN(35), CKIN(36) : (D=0.,-1.) range of allowed abs(t-hat) (= -that)
    values in 2 -> 2 processes. Note that for deep inelastic scattering
    this is nothing but the Q2 scale, in the limit that initial and
    final state radiation is neglected. If CKIN(36) < 0., the upper
    limit is inactive.
CKIN(37), CKIN(38) : (D=0.,-1.) range of allowed abs(u-hat) (= -uhat)
    values in 2 -> 2 processes. If CKIN(38) < 0., the upper limit is
    inactive.
CKIN(41) - CKIN(44) : (D=12.,-1.,12.,-1.) range of allowed mass values
    of the two (or one) resonances produced in a 'true' 2 -> 2
    process, i.e. one not (only) proceeding through a single s-channel
    resonance (2 -> 1 -> 2). (These are the ones listed as 2 -> 2 in the
    table in section 2.2.) Only particles with a width above PARP(41)
    are considered as bona fide resonances and tested against the CKIN
    limits; particles with a smaller width are put on mass-shell
    without applying any cuts. The exact interpretation of the CKIN
    variables depends on the flavours of the produced two resonances.
    For two resonances like Z + W (produced from f + f' -> Z + W),
    which are not identical and which are not each other's
    antiparticles, one has
    CKIN(41) < m1 < CKIN(42), and
    CKIN(43) < m2 < CKIN(44),
    where m1 and m2 are the actually generated masses of the two
    resonances, and 1 and 2 are defined by the order in which they
    given in the production process specification, see section 2.2.
    For two resonances like Z + Z, which are identical, or W+ + W-,
    which are each other's antiparticles, one instead has
    CKIN(41) < min(m1,m2) < CKIN(42), and
    CKIN(43) < max(m1,m2) < CKIN(44).
    In addition, whatever limits are set on CKIN(1) and, in particular,
    CKIN(2) obviously affects the masses actually selected.
    Note 1: If MSTP(42) = 0, so that no mass smearing is allowed,
        the CKIN values have not effect (the same as for particles
        with too narrow a width).
    Note 2: If CKIN(42) < CKIN(41) or CKIN(44) < CKIN(43) it means
        that the CKIN(42) or CKIN(44) limit is inactive.
    Note 3: If limits are active and the resonances are identical, it
        is up to the user to ensure that CKIN(41) <= CKIN(43) and
        CKIN(42) <= CKIN(44).
    Note 4: For identical resonances, it is not possible to preselect
        which of the resonances is the lighter one; if e.g. one Z is
        to decay to leptons and the other to quarks, there is no
        mechanism to guarantee that the lepton pair has a smaller mass
        than the quark one.
    Note 5: The CKIN values are applied to all relevant 2 -> 2
        processes equally, which may not be what one desires if several
        processes are generated simultaneously. Some caution is
        therefore urged in the use of the CKIN(41) - CKIN(44) values.
        Also in other respects, users are recommended to take proper
        care - if a Z is only allowed to decay into b + bbar, e.g.,
        setting its mass range to be 2 - 8 GeV is obviously not a
        good idea.
    Note 6: In principle, the machinery should work for any 2 -> 2
        process with resonances in the final state, but so far it has
        only been checked for processes 22 - 26, so also from this
        point some caution is urged.
CKIN(45) - CKIN(48) : (D=12.,-1.,12.,-1.) range of allowed mass values
    of the two (or one) secondary resonances produced in 2 -> 1 -> 2
    process (like g + g -> H0 -> Z0 + Z0) or even a 2 -> 2 -> 4 (or 3)
    process (like q + qbar -> Z0 + H0 -> Z0 + W+ + W-). Note that these
    CKIN values only affect the secondary resonances; the primary ones
    are constrained by CKIN(1), CKIN(2) and CKIN(41) - CKIN(44).
    (indirectly, of course, the choice of primary resonance masses
    affects the allowed mass range for the secondary ones).
    What is considered to be a resonance is defined by PARP(41);
    particles with a width smaller than this are automatically put on
    mass-shell. The description closely parallels the one given for
    CKIN(41) - CKIN(44). Thus, for two resonances which are not
    identical or each other's antiparticles, one has
    CKIN(45) < m1 < CKIN(46), and
    CKIN(47) < m2 < CKIN(48),
    where m1 and m2 are the actually generated masses of the two
    resonances, and 1 and 2 are defined by the order in which they
    given in the decay channel specification in the program (see e.g.
    output from PYSTAT(2) or LULIST(12)). For two resonances which
    are identical or each other's antiparticles, one instead has
    CKIN(45) < min(m1,m2) < CKIN(46), and
    CKIN(47) < max(m1,m2) < CKIN(48).
    Notes 1 - 5: as for CKIN(41) - CKIN(44), with trivial modifications.
    Note 6: Setting limits on secondary resonance masses is possible
        in any of the channels of the allowed types (see above).
        However, so far only H -> Z0 + Z0 and H -> W+ + W- have been
        fully implemented, such that an arbitrary mass range below
        the naive mass threshold may be picked. For other possible
        resonances, any restrictions made on the allowed mass range
        are not reflected in the cross-section; and further it is not
        recommendable to pick mass windows that makes an on-mass-shell
        decay impossible. These limitations will be relaxed in future
        versions.
CKIN(51) - CKIN(56) : (D=0.,-1.,0.,-1.,0.,-1.) range of allowed
    transverse momenta in a true 2 -> 3 process. Currently two
    different alternatives are around. For subprocess 131, the p_T
    of the first product (the Z) is set by CKIN(3) and CKIN(4),
    while for the quark and antiquark p_T:s one has
    CKIN(51) < min(p_T_q,p_T_q~) < CKIN(52)
    CKIN(53) < max(p_T_q,p_T_q~) < CKIN(54)
    Negative CKIN(52) and CKIN(54) values means that the corresponding
    limits are inactive.
    For subprocesses 121-124, and their H'0 and A0 equivalents (173,
    174,178,179,181,182,186,187), CKIN(51) - CKIN(54) again corresponds
    to p_T ranges for scattered partons, but in order of appearance,
    i.e. CKIN(51) - CKIN(52) for the parton scattered off the beam and
    CKIN(53) - CKIN(54) for the one scattered off the target. CKIN(55)
    and CKIN(56) here sets p_T limits for the third product, the H0,
    i.e. the CKIN(3) and CKIN(4) values have no effect for this
    process. Since the p_T of the Higgs is not one of the primary
    variables selected, any constraints here may mean reduced Monte
    Carlo efficiency, while for these processes CKIN(51) - CKIN(54)
    are 'hardwired'  and therefore do not cost anything.

______________________________________________________________________

2.5.  The General Switches and Parameters

In addition to the event information described in section 2.6, the
PYPARS commonblock contains the status code and parameters which
regulate the performance of the program. All of them are provided with
sensible default values, so that a novice user can neglect them, and
only gradually explore the full range of possibilities.

COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
Purpose: to give access to status code and parameters which regulate the
    performance of the program. If the default values, below denoted by
    (D=...), are not satisfactory, they must in general be changed
    before the PYINIT call. Exceptions, i.e. variables which can be
    changed for each new event, are denoted by (C).
MSTP(1) : (D=3) maximum number of generations. Automatically set <= 4.
MSTP(2) : (D=1) calculation of alpha_strong at hard interaction, in the
   routine ULALPS.
   = 0 : alpha_strong is fixed at value PARU(111).
   = 1 : first order running alpha_strong.
   = 2 : second order running alpha_strong.
MSTP(3) : (D=2) selection of Lambda value in alpha_strong for
    MSTP(2) >= 1.
    = 1 : Lambda is given by PARP(1) for hard interactions, by PARP(61)
        for spacelike showers and by PARJ(81) for timelike ones. This
        Lambda is assumed valid for 5 flavours; for the hard interaction
        the number of flavours assumed can be changed by MSTU(112).
    = 2 : Lambda value is chosen according to the structure function
        parametrizations, i.e. Lambda = 0.20 GeV for EHLQ1, = 0.29 GeV
        for EHLQ2, = 0.20 GeV for DO1, = 0.40 GeV for DO2, = 0.187 GeV
        for MT1, = 0.212 GeV for MT2, = 0.191 GeV for MT3, = 0.155 GeV
        for MT4, = 0.22 GeV for GRV1, = 0.16 GeV for GRV2, etc. for
        structure functions in the PDFLIB and PAKPDF libraries (cf.
        (MSTP(51), MSTP(52)). The choice is always based on the proton
        structure function set selected, i.e. is unaffected by pion and
        photon structure function selection. All the Lambda values above
        are assumed to refer to 4 flavours, and MSTU(112) is set
        accordingly. This Lambda value is used both for the hard
        scattering and the initial and final state radiation. The
        ambiguity in the choice of Q^2 argument still remains (see
        MSTP(32), MSTP(64) and MSTJ(44)). This Lambda value is used
        also for MSTP(57) = 0, but the sensible choice here would be
        to use MSTP(2) = 0 and have no initial or final state radiation.
MSTP(4) : (D=0) treatment of the Higgs sector, predominantly the
    neutral one.
    = 0 : the H is given the standard model Higgs couplings, while
        H'0 and A0 couplings should be set by the user in PARU(171) -
        PARU(175) and PARU(181) - PARU(185), respectively.
    = 1 : the user should set couplings for all three Higgses, for
        the H0 in PARU(161) - PARU(165), and for the H'0 and A0 as
        above.
    = 2 : the mass of H0 in PMAS(25,1) and the tan(beta) value in
        PARU(141) are used to derive H'0, A0 and H+- masses, and H0,
        H'0, A0 and H+- couplings, using the relations of the minimal
        supersymmetric extension to the standard model at Born level
        [GHK90]. Existing mass and couplings are overwritten by the
        derived values. See section 2.3.9 for discussion in parameter
        constraints.
    = 3: as =2, but using relations at the one loop level. This option
        is not yet implemented.
MSTP(5) : (D=0) presence of anomalous couplings in processes.
    = 0 : absent.
    >= 1 : present, wherever implemented. See section 2.3.11 for further
        details.
MSTP(6) : (D=0) usage of the fourth generation fermions to simulate
    other fermions.
    = 0 : none, i.e. can be used as a standard forurht generation.
    = 1 : excited fermions, as present in compositeness scenarios;
        see section 2.3.12.
MSTP(11) : (D=0) use of electron structure function in e+e- and ep
    interactions.
    = 0 : no, i.e. electron carries whole beam energy.
    = 1 : yes, i.e. electron only carries fraction of beam energy in
        agreement with next-to-leading electron structure function
        in [Kle89] (p. 34, eqs. (58), (60)), thereby including the
        effects of initial state bremsstrahlung.
    = 2 : include the effects of beamstrahlung in terms of a simple,
        approximate effective electron structure function, of the
        type D(x) = N * x^a (1-x)^b, where a = PARP(11), b = PARP(12),
        and N a normalization factor so as to make the integral of
        D(x) equal to one. The mean x value is <x> = (1+a)/(2+a+b).
    = 3: include, approximately, the effects of both beamstrahlung and
        bremsstrahlung, in terms of a parametrization as in = 2, but
        with b = PARP(12) + 0.5 * beta, where
        beta = (2.*alpha_em/pi) * (ln(Q^2/m_e^2) - 1)
        is the traditional factor associated with bremsstrahlung.
MSTP(12) : (D=0) use of e- ('sea', i.e. from e -> gamma -> e), e+,
        quark and gluon structure functions inside an electron.
    = 0 : off.
    = 1 : on, provided that MSTP(11) >= 1. The e- and e+ structure
        functions are given by [Che75], divided by two (reduction
        due to Q2 ordering). The other structure functions
        are obtained by numerical convolution of the photon content
        inside an electron (as given by the bremsstrahlung spectrum of
        MSTP(11) = 1) with the quark and gluon content inside a photon
        (as given by the Drees and Grassie structure functions).
        Required numerical precision is set by PARP(14).
        Since the need for numerical integration makes this option
        somewhat more time-consuming than ordinary structure function
        evaluation, one should only use it when studying processes
        where it is needed.
MSTP(13) : (D=1) choice of Q^2 range over which electrons are assumed
    to radiate photons; affects normalization of e- (sea), e+, photon,
    quark and gluon structure functions inside electron.
    = 1 : range set by Q^2 argument of structure function call, i.e.
        by Q^2 scale of hard interaction. Therefore structure functions
        are proportional to ln(Q^2/m_e^2). This is normally most
        appropriate for e+e- annihilation.
    = 2 : range set by user determined Q_max^2, given in PARP(13).
        Structure functions are assumed proportional to
        ln((Q_max^2/m_e^2)*(1-x)/x^2). This is normally most appropriate
        for photoproduction, where the electron is supposed to go
        undetected, i.e. scatter less than Q_max^2.
MSTP(14) : (D=1) structure of incoming photon beam or target (does not
    affect photon inside electron, only photons appearing as argument
    in the PYINIT call).
    = 0 : photon is assumed pointlike, i.e. can only interact in
        processes which contain an incoming photon, such as ISUB =
        33 - 36.
    = 1 : photon is assumed resolved; i.e. can only interact through
        its constituent quarks and gluons, in processes such as
        ISUB = 11 - 31.
    Note: currently it is not possible to mix the two kinds of
        photon-induced processes in one run.
MSTP(21) : (D=1) nature of fermion-fermion scatterings simulated in
    process 10 by t-channel exchange.
    = 0 : all off.
    = 1 : full mixture of gamma/Z neutral current and W charged current.
    = 2 : gamma neutral current only.
    = 3 : Z neutral current only.
    = 4 : gamma/Z neutral current only.
    = 5 : W charged current only.
MSTP(22) : (D=0) special override of normal Q^2 definition used for
    maximum of parton shower evolution. This option only affects
    processes 10 and 83 (Deep Inelastic Scattering) and only in
    lepton-hadron events.
    = 0 : use the scale as given in MSTP(32).
    = 1 : use the DIS Q^2 scale, i.e. -that.
    = 2 : use the DIS W^2 scale, i.e. (1-x)/x * (-that).
    = 3 : use the DIS Q*W scale, i.e. sqrt((1-x)/x) * (-that).
    Note: in all of these alternatives, a multiplicative factor is
       introduced by PARP(67) and PARP(71).
MSTP(23) : (D=0) for Deep Inelastic Scattering processes (10 and 83)
    this option allows the x and Q^2 of the original hard scattering
    to be retained by the final electron.
    = 0 : no correction procedure, i.e. x and Q^2 of the scattered
        electron differ from the originally generated x and Q^2.
    = 1 : post facto correction, i.e. the change of electron momentum
        by initial and final QCD radiation, primordial k_T and beam
        remnant treatment is corrected for by a shuffling of momentum
        between the electron and hadron side in the final state.
    Note: the correction procedure will fail for a fraction of the
        events, which thus are thrown (and new ones generated in their
        place). The correction option is not unambiguous, and should
        not be taken too seriously. For very small Q^2 values, the x
        is not exactly preserved even after this procedure.
MSTP(31) : (D=1) parametrization of total and elastic cross-sections,
    nuclear slope parameter B  and curvature C [Blo85].
    = 1 : Block-Cahn fit 1 for cross-section, fit 1 for slope parameter.
    = 2 : Block-Cahn fit 2 for cross-section, fit 1 for slope parameter.
    = 3 : Block-Cahn fit 3 for cross-section, fit 1 for slope parameter.
    = 4 : Block-Cahn fit 6 for cross-section, fit 2 for slope parameter.
    = 5 : Block-Cahn fit 8 for cross-section, fit 2 for slope parameter.
    Note: sets 1-3 for cross-section and set 1 for slope parameter are
        fits excluding recent measurements from Spp~S, whereas sets 4-5
        for cross-section and set 2 for slope parameter are fits
        including the Spp~S measurements.
MSTP(32) : (D=2) Q^2-definition in hard scattering for 2 -> 2 processes;
    for resonance production Q^2 is always chosen to be m_res^2, where
    m_res is the mass of the resonance.
    = 1 : Q^2 = 2*shat*that*uhat/(shat^2 + that^2 + uhat^2).
    = 2 : Q^2 = (m_T1^2 + m_T2^2)/2.
    = 3 : Q^2 = min(-that, -uhat).
    = 4 : Q^2 = shat.
    = 5 : Q^2 = -that.
MSTP(33) : (D=0) inclusion of K factors in hard cross-sections for
    parton-parton interactions (i.e. for incoming quarks and gluons).
    = 0 : none, i.e. K = 1.
    = 1 : a common K factor is used, as stored in PARP(31).
    = 2 : separate factors are used for ordinary (PARP(31)) and colour
        annihilation graphs (PARP(32)).
    = 3 : A K factor is introduced by a shift in the alpha_strong Q^2
        argument, alpha_strong = alpha_strong(PARP(33)*Q^2).
MSTP(34) : (D=0) use of interference term in matrix elements for QCD
    processes.
    = 0 : excluded (i.e. string-inspired matrix elements).
    = 1 : included (i.e. conventional QCD matrix elements).
    Note: for the option MSTP(34)=1, i.e., interference terms included,
        these are divided between the different possible colour
        configurations according to the pole structure of the
        (string-inspired) matrix elements for the different colour
        configurations.
MSTP(35) : (D=0) threshold behaviour for heavy flavour production, i.e.
    ISUB = 81, 82, 84, 85, and also for Z and Z' decays. The nonstandard
    options are mainly intended for top, but can be used, with less
    theoretical reliability, also for charm and bottom (for Z and Z'
    only top and heavier flavours are affected). See [Fad89] for
    details.
    = 0 : naive lowest order matrix element behaviour.
    = 1 : enhancement or suppression close to threshold, according to
        the colour structure of the process, see [Fad89] for the
        expressions actually used. The alpha_strong value appearing
        in the threshold factor (which is not the same as the
        alpha_strong of the lowest order 2 -> 2 process) is taken to
        be fixed at the value given in PARP(35). The threshold factor
        used in an event is stored in PARI(81).
    = 2 : as =1, but the alpha_strong value appearing in the threshold
        factor is taken to be running, with argument
        Q^2 = m_Q sqrt( (m-hat - 2m_Q)^2 + Gamma_Q^2).
        Here m_Q is the nominal heavy quark mass, Gamma_Q is the width
        of the heavy quark mass distribution, and m-hat is the invariant
        mass of the heavy quark pair. The Gamma_Q value has to be stored
        by the user in PARP(35). The regularization of alpha_strong at
        low Q^2 is given by MSTP(36).
MSTP(36) : (D=2) regularization of alpha_strong in the limit Q^2 -> 0
    for the threshold factor obtainable in the MSTP(35) = 2 option;
    see MSTU(115) for a list of the possibilities.
MSTP(37) : (D=1) inclusion of running quark masses in Higgs production
    (q + qb -> H0) and decay (H0 -> q + qb) couplings.
    = 0 : not included, i.e. fix quark masses are used according to
        the values in the PMAS array.
    = 1 : included, with running starting from the value given in
        the PMAS array, at a Q (i.e. Higgs mass) scale given by
        PARP(37) times the quark mass itself.
MSTP(38) : (D=5) handling of quark loop masses in the box graphs
    g + g -> gamma + gamma and g + g -> g + gamma.
    = 0 : the program will for each flavour automatically choose
        the massless approximation for light quarks and the full
        massive formulae for heavy quarks, with the dividing line
        between light and heavy quarks dependent on the actual
        shat scale.
    >= 1 : the program will use the massless approximation throughout,
        assuming the presence of MSTP(38) effectively massless quark
        species (however, at most 8). Normally one would use = 5 for
        the inclusion of all quarks up to bottom, and = 6 to include
        top as well.
    WARNING: for =0, numerical instabilities may arise for scattering
        at small angles. Users are therefore recommended in this case
        to set CKIN(27) and CKIN(28) so as to exclude the range of
        scattering angles that are not of interest anyway.
MSTP(41) : (D=1) master switch for all resonance decays (Z0, W+-, H0,
    Z'0, W'+-, H'0, A0, H+-, LQ, R).
    = 0 : off.
    = 1 : on.
    Note: also for MSTP(41) = 1 it is possible to switch off the decays
        of specific resonances by using the MDCY(KC,1) switches in
        JETSET. However, since the MDCY values are overwritten in the
        PYINIT call, individual resonances should be switched off after
        the PYINIT call.
    WARNING: leptoquark decays must not be switched off if one later
        on intends to let it decay (with LUEXEC); see 2.3.10.
MSTP(42) : (D=1) mass treatment in 2 -> 2 processes, where the final
    state resonances have finite width (see PARP(41)).
    (Does not apply for the production of a single s-channel resonance,
    where the mass appears explicitly in the cross-section of the
    process, and thus is always selected with width.)
    = 0 : particles are put on mass-shell.
    = 1 : mass generated according to a Breit-Wigner.
MSTP(43) : (D=3) treatment of Z0/gamma* interference in matrix elements.
    So far only implemented in subprocesses 1 and 22; in other processes
    what is called a Z0 is really a Z0 only, without the gamma* piece.
    = 1 : only gamma* included.
    = 2 : only Z0 included.
    = 3 : complete Z0/gamma* structure (with interference) included.
MSTP(44) : (D=7) treatment of Z'0/Z0/gamma* interference in matrix
    elements.
    = 1 : only gamma* included.
    = 2 : only Z0 included.
    = 3 : only Z'0 included.
    = 4 : only Z0/gamma* (with interference) included.
    = 5 : only Z'0/gamma* (with interference) included.
    = 6 : only Z'0/Z0 (with interference) included.
    = 7 : complete Z'0/Z0/gamma* structure (with interference) included.
MSTP(45) : (D=2) treatment of WW -> WW structure (ISUB = 77).
    = 1 : only W+W+ -> W+W+ and W-W- -> W-W- included.
    = 2 : only W+W- -> W+W- included.
    = 3 : all charge combinations WW -> WW included.
MSTP(46) : (D=1) treatment of VV -> V'V' structures (ISUB = 71 - 77).
    = 0 : only s-channel Higgs exchange included (where existing). With
        this option, subprocesses 71 - 72 and 76 - 77 will essentially
        be equivalent to subprocesses 5 and 8, respectively, with the
        proper decay channels (i.e. only Z0Z0 or W+W-) set via MDME.
        The description obtained for subprocesses 5 and 8 in this case
        is more sophisticated, however; see section 2.3.6.
    = 1 : all graphs contributing to VV -> V'V' processes are included.
    = 2 : only graphs not involving Higgs exchange (either in s, t or u
        channel) are included; this option then gives the naive
        behaviour if no Higgs exists, including  unphysical unitarity
        violations at high energies.
    = 3 : the strongly interacting Higgs-like model of Dobado, Herrero
        and Terron [DHT90] with Pade unitarization. Note that to use
        this option it is necessary to set the Higgs mass to a large
        number like 20 TeV (i.e. PMAS(25,1)=20000.). The parameter nu
        is stored in PARP(44), but should not have to be changed.
    = 4 : as =3, but with K-matrix unitarization.
    = 5 : the strongly interacting QCD-like model of Dobado, Herrero
        and Terron [DHT90] with Pade unitarization. The parameter nu
        is stored in PARP(44), but should not have to be changed.
        The effective techni-rho mass in this model is stored in
        PARP(45); by default it is 2054 GeV, which is the expected
        value for three technicolours, based on scaling up the ordinary
        rho mass appropriately.
    = 6 : as =5, but with K-matrix unitarization. While PARP(45) still
        is a parameter of the model, this type of unitarization does
        not give rise to a resonance at a mass of PARP(45).
MSTP(47) : (D=1) (C) angular orientation of decay products of resonances
    (Z, W, H, Z', W', R), either when produced singly or in pairs (also
    from a H0 decay), or in combination with a single quark, gluon or
    photon.
    = 0 : independent decay of each resonance, isotropic in CM frame of
        the resonance.
    = 1 : correlated decay angular distributions according to proper
        matrix elements, to the extent these are known.
MSTP(51) : (D=1) choice of proton structure functions; see also
    MSTP(52).
    = 1 : EHLQ set 1 (1986 updated version).
    = 2 : EHLQ set 2 (1986 updated version).
    = 3 : Duke-Owens set 1.
    = 4 : Duke-Owens set 2.
    = 5 : Morfin-Tung set 1.
    = 6 : Morfin-Tung set 2.
    = 7 : Morfin-Tung set 3.
    = 8 : Morfin-Tung set 4.
    = 9 : Gluck-Reya-Vogt LO set.
    = 10 : Gluck-Reya-Vogt HO set.
    Note: since all parametrizations have some region of applicability,
        the structure functions are assumed frozed below the lowest Q^2
        and above the highest Q^2 covered by the parametrizations. The
        extrapolation to low x is covered by PARP(51).
MSTP(52) : (D=1) choice of proton structure function library.
    = 1 : the internal PYTHIA one, with structure functions according
        to MSTP(51) above.
    = 2 : the PDFLIB one [Plo91], with the PDFLIB MODE (set) number
        to be given in MSTP(51).
    = 3 : the PAKPDF one [Cha91], with the PAKPDF 100*IPAR+ISET number
        to be given in MSTP(51).
    NOTE: to make use of options 2 and 3, it is necessary to enable
        the calls to PDFLIB and PAKPDF, respectively, in subroutines
        PYINIT and PYSTFU. What to do is described inline in the
        Fortran file; look for lines with C! in the first two columns.
    WARNING: For external structure function libraries, PYTHIA does not
        check whether MSTP(51) corresponds to a valid code, neither if
        special x and Q2 restrictions exist for a given set, such that
        crazy values could be returned. This puts an extra
        responsibility on the user.
MSTP(53) : (D=1) choice of pion structure functions; see also
    MSTP(54).
    = 1 : Owens set 1.
    = 2 : Owens set 2.
MSTP(54) : (D=1) choice of pion structure function library.
    = 1 : the internal PYTHIA one, with structure functions according
        to MSTP(53) above.
    = 2 : the PDFLIB one [Plo91], with the PDFLIB MODE (set) number
        to be given in MSTP(53).
    NOTE: to make use of option 2, it is necessary to enable the calls
        to PDFLIB in subroutines PYINIT and PYSTFU. What to do is
        described inline in the Fortran file; look for lines with
        C! in the first two columns.
    WARNING: For external structure function libraries, PYTHIA does not
        check whether MSTP(53) corresponds to a valid code, neither if
        special x and Q2 restrictions exist for a given set, such that
        crazy values could be returned. This puts an extra
        responsibility on the user.
MSTP(55) : (D=1) choice of photon structure functions; see also
    MSTP(52). Currently there are no alternatives available for
    this switch.
    = 1 : Drees-Grassie.
MSTP(56) : (D=1) choice of photon structure function library.
    = 1 : the internal PYTHIA one, with structure functions according
        to MSTP(53) above.
    = 3 : the PHOPDF one [Cha91a], with the PHOPDF 100*IPAR+ISET number
        to be given in MSTP(55).
    NOTE: to make use of option 3, it is necessary to enable the calls
        to PHOPDF, in subroutines PYINIT, PYSTFU and PYSTEL. What to
        do is described inline in the Fortran file; look for lines with
        C! in the first two columns.
    WARNING: For external structure function libraries, PYTHIA does not
        check whether MSTP(53) corresponds to a valid code, neither if
        special x and Q2 restrictions exist for a given set, such that
        crazy values could be returned. This puts an extra
        responsibility on the user.
MSTP(57) : (D=1) choice of Q^2 dependence in structure functions.
    This option only applies to the proton, pion and photon structure
    functions that come with PYTHIA, and not to externally linked ones.
    = 0 : structure functions are evaluated at nominal lower cutoff
        value Q_0^2, i.e. are made Q^2-independent.
    = 1 : the parametrized Q^2 dependence is used.
MSTP(58) : (D=min(6,2*MSTP(1))) maximum number of quark flavours used
    in structure functions, and thus also for initial state spacelike
    showers. If some distributions (notably t) are absent in the
    parametrization selected in MSTP(51), these are obviously
    automatically excluded.
MSTP(61) : (D=1) (C) master switch for initial state QCD and QED
    radiation.
    = 0 : off.
    = 1 : on.
MSTP(62) : (D=2) (C) level of coherence imposed on the spacelike parton
    shower evolution.
    = 1 : none, i.e. neither Q^2 values nor angles need be ordered.
    = 2 : Q^2 values at branches are strictly ordered, increasing
        towards the hard interaction.
    = 3 : Q^2 values and opening angles of emitted (on-mass-shell or
        timelike) partons are both strictly ordered, increasing
        towards the hard interaction.
MSTP(63) : (D=2) (C) structure of associated timelike showers, i.e.
    showers initiated by emission off the incoming spacelike partons.
    = 0 : no associated showers are allowed, i.e. emitted partons are
        put on mass-shell.
    = 1 : a shower may evolve, with maximum allowed timelike virtuality
        set by the phase space only.
    = 2 : a shower may evolve, with maximum allowed timelike virtuality
        set by phase space or by PARP(71) times the Q^2 value of the
        spacelike parton created in the same vertex, whichever is the
        stronger constraint.
    = 3 : a shower may evolve, with maximum allowed timelike virtuality
        set by phase space or by the requirement of angular ordering,
        whichever is the stronger constraint. (Not yet implemented!)
MSTP(64) : (D=2) (C) choice of alpha_strong and Q^2 scale in spacelike
    parton showers.
    = 0 : alpha-strong is taken to be fix at the value PARU(111).
    = 1 : first order running alpha_strong with argument PARP(63)*Q^2.
    = 2 : first order running alpha_strong with argument
        PARP(64)*k_T^2 = PARP(64)*(1-z)*Q^2.
MSTP(65) : (D=1) (C) treatment of soft gluon emission in spacelike
    parton shower evolution.
    = 0 : soft gluons are entirely neglected.
    = 1 : soft gluon emission is resummed and included together
        with the hard radiation as an effective z shift.
MSTP(71) : (D=1) (C) master switch for final state QCD and QED
    radiation.
    = 0 : off.
    = 1 : on.
    Note: additional switches (e.g. for conventional/coherent showers)
        are available in JETSET, see MSTJ(41) - MSTJ(45), PARJ(81) and
        PARJ(82). In particular, to allow photon emission off quarks
        and leptons, it is necessary to put MSTJ(41)=2.
MSTP(81) : (D=1) master switch for multiple interactions.
    = 0 : off.
    = 1 : on.
MSTP(82) : (D=1) structure of multiple interactions. For QCD processes,
    used down to p_T values below p_Tmin, it also affects the choice
    of structure for the one hard/semihard interaction.
    = 0 : simple two-string model without any hard interactions.
    = 1 : multiple interactions assuming the same probability in all
        events, with an abrupt p_Tmin cutoff at PARP(81).
    = 2 : multiple interactions assuming the same probability in all
        events, with a continuous turnoff of the cross-section at
        p_T0 = PARP(82).
    = 3 : multiple interactions assuming a varying impact parameter
        and a hadronic matter overlap consistent with a Gaussian matter
        distribution, with a continuous turnoff of the cross-section at
        p_T0 = PARP(82).
    = 4 : multiple interactions assuming a varying impact parameter
        and a hadronic matter overlap consistent with a double Gaussian
        matter distribution given by PARP(83) and PARP(84), with a
        continuous turnoff of the cross-section at p_T0 = PARP(82).
    Note 1: For MSTP(82) >= 2 and CKIN(3) > PARP(82), cross-sections
        given with PYSTAT(1) may be somewhat too large, since (for
        reasons of efficiency) the probability factor that the hard
        interaction is indeed the hardest in the event is not
        included in the cross-sections. It is included in the event
        selection, however, so the events generated are correctly
        distributed. For CKIN(3) values a couple of times larger than
        PARP(82) this ceases to be a problem.
    Note 2: The PARP(81) and, in particular, PARP(82) values are
        sensitive to the choice of structure functions, Lambda_QCD,
        etc., in the sense that a change in the latter variables
        leads to a net change in the multiple interaction rate, which
        has to be compensated by a retuning of PARP(81) or PARP(82)
        if one wants to keep the net multiple interaction structure the
        same. The default PARP(81) value is consistent with the other
        default values give, i.e. EHLQ set 1 structure functions etc.
        When options MSTP(82) = 2 - 4 are used, the default PARP(82)
        value is to be used in conjunction with MSTP(2) = 2 and
        MSTP(33) = 3. These switches should be set by the user.
MSTP(83) : (D=100) number of Monte Carlo generated phase space points
        per bin (whereof there are 20) in the initialization (in PYINMU)
        of multiple interactions for MSTP(82) >= 2.
MSTP(91) : (D=1) (C) primordial k_T distribution in hadron. See MSTP(93)
    for photon.
    = 0 : no primordial k_T.
    = 1 : Gaussian, width given in PARP(91), upper cutoff in PARP(93).
    = 2 : exponential, width given in PARP(92), upper cutoff in
        PARP(93).
MSTP(92) : (D=4) (C) energy partitioning in hadron remnant. The energy
    fraction chi taken by one of the two objects, with conventions as
    described for PARP(94) - PARP(98), is chosen according to the
    different distributions below. Here c_min = 2*m_quark/sqrt(s) =
    0.6 GeV/sqrt(s).
    = 1 : 1 for meson, 2*(1-chi) for baryon, i.e. simple counting rules.
    = 2 : (k+1)*(1-chi)^k, with k as given in PARP(94) - PARP(97).
    = 3 : as =2 for remnant splitting into hadron plus jet, but
        proportional to (1-chi)^k/(chi^2+c_min^2)^0.25 for remnant
        splitting into two jets, with k given by PARP(94) or PARP(96).
    = 4 : as =2 for remnant splitting into hadron plus jet, but
        proportional to (1-chi)^k/(chi^2+c_min^2)^0.5 for remnant
        splitting into two jets, with k given by PARP(94) or PARP(96).
    = 5 : as =2 for remnant splitting into hadron plus jet, but
        proportional to (1-chi)^k/(chi^2+c_min^2)^(PARP(98)/2) for
        remnant splitting into two jets, with k given by PARP(94) or
        PARP(96).
MSTP(93) : (D=1) (C) primordial k_T distribution in photon, either it
    is one of the incoming particles or inside an electron.
    = 0 : no primordial k_T.
    = 1 : Gaussian, width given in PARP(99), upper cutoff in PARP(100).
    = 2 : exponential, width given in PARP(99), upper cutoff in
        PARP(100).
    = 3 : powerlike of the type d(k_T^2)/(PARP(99)^2+k_T^2)^2, with
        upper k_T cutoff in PARP(100).
    = 4 : powerlike of the type d(k_T^2)/(PARP(99)^2+k_T^2), with
        upper k_T cutoff in PARP(100).
    = 5 : powerlike of the type d(k_T^2)/(PARP(99)^2+k_T^2), with
        upper k_T cutoff given by the p_T of the hard process or by
        PARP(100), whichever is smaller.
    NOTE: for options 1 and 2 the PARP(100) value is of minor
        importance, once PARP(100) >> PARP(99). However, options 3 and
        4 correspond to distributions with infinite <k_T^2> if the
        k_T spectrum is not cut off, and therefore the PARP(100) value
        is as important for the overall distribution as is PARP(99).
MSTP(101) : (D=1) (C) structure of diffractive system.
    = 1 : forward moving diquark + interacting quark.
    = 2 : forward moving diquark + quark joined via interacting gluon
          ('hairpin' configuration).
MSTP(111) : (D=1) (C) master switch for fragmentation and decay, as
    obtained with a LUEXEC call.
    = 0 : off.
    = 1 : on.
    = -1 : only choose kinematical variables for hard scattering,
        i.e. no jets are defined. This is useful e.g. to calculate
        cross-sections (by Monte Carlo integration) without wanting
        to simulate events; information obtained with PYSTAT(1)
        will be correct.
MSTP(112) : (D=1) (C) cuts on partonic events; only affects an
    exceedingly tiny fraction of events.
    = 0 : no cuts (can be used only with independent fragmentation,
        at least in principle).
    = 1 : string cuts (as normally required for fragmentation).
MSTP(113) : (D=1) (C) recalculation of energies of partons from their
    momenta and masses, to be done immediately before fragmentation, to
    compensate in parts for some numerical problems appearing at high
    energies.
    = 0 : not performed.
    = 1 : performed.
MSTP(121) : (D=0) calculation of kinematics selection coefficients and
    differential cross-section maxima for subprocesses included
    (by user or default).
    = 0 : not known; to be calculated at initialization.
    = 1 : not known; to be calculated at initialization; however, the
        maximum value then obtained is to be multiplied by PARP(121)
        (this may be useful if a violation factor has been observed
        in a previous run of the same kind).
    = 2 : known; kinematics selection coefficients stored by user in
        COEF(ISUB,J) (J = 1 - 20) in /PYINT2/ and maximum of the
        corresponding differential cross-section times Jacobians in
        XSEC(ISUB,1) in /PYINT5/. This is to be done for each included
        subprocess ISUB before initialization, with the sum of all
        XSEC(ISUB,1) values, except for ISUB = 95, stored in XSEC(0,1).
MSTP(122) : (D=1) initialization and differential cross-section
    maximization printout (see also MSTP(127)).
    = 0 : none.
    = 1 : short message.
    = 2 : detailed message, including full maximization.
MSTP(123) : (D=2) reaction to violation of maximum differential cross-
    section.
    = 0 : stop generation, print message.
    = 1 : continue generation, print message for each subsequently
          larger violation.
    = 2 : as 1, but also increase value of maximum.
MSTP(124) : (D=1) (C) frame for presentation of event.
    = 1 : as specified in PYINIT.
    = 2 : CM frame of incoming particles.
MSTP(125) : (D=1) (C) documentation of partonic process, see section
    'The event record' for details.
    = 0 : only list ultimate string/particle configuration.
    = 1 : additionally list short summary of the hard process.
    = 2 : list complete documentation of intermediate steps of parton
        shower evolution.
MSTP(126) : (D=20) number of lines in the beginning of event record
    that are reserved for event history information; see section 2.7.
    This value should never be reduced, but may be increased at a
    later date if more complicated processes are included.
MSTP(127) : (D=1) writing of header (version number and last date of
    change) on output file.
    = 0 : not done.
    = 1 : header is written at first PYINIT call, at which time
        MSTP(127) is set =0.
MSTP(128) : (D=0) storing of copy of resonance decay products in
    the documentation section of the event record, and mother
    pointer (K(I,3)) relation of the actual resonance decay
    products (stored in the main section of the event record)
    to the documentation copy.
    = 0 : products are stored also in the documentation section,
        and each product stored in the main section points back
        to the corresponding entry in the documentation section.
    = 1 : products are stored also in the documentation section,
        but the products stored in the main section point back to
        the decaying resonance copy in the main section.
    = 2 : products are not stored in the documentation section,
        the products stored in the main section point back to the
        the decaying resonance copy in the main section.
MSTP(129) : (D=10) for the maximization of 2 -> 3 processes
    (ISET(ISUB)=5) each phase space point in tau, y* and tau' is
    tested MSTP(129) times in the other dimensions (at randomly
    selected points) to determine the effective maximum in the
    (tau, y*, tau') point.
MSTP(131) : (D=0) master switch for pileup events, i.e. several
    independent hadron-hadron interactions generated in the same
    bunch-bunch crossing, with the events following one after the
    other in the event record.
    = 0 : off, i.e. only one event is generated at a time.
    = 1 : on, i.e. several events are allowed in the same event record.
        Information on the processes generated may be found in MSTI(41)
        - MSTI(50).
MSTP(132) : (D=4) the processes that are switched on for pileup events.
    The first event may be set up completely arbitrarily, using the
    switches in the PYSUBS commonblock, while all the subsequent events
    have to be of one of the 'inclusive' processes which dominate the
    cross-section, according to the options below. It is thus not
    possible to generate two rare events in the pileup option.
    = 1 : low-p_T processes (ISUB = 95) only. The low-pT model actually
        used, both in the hard event and in the pileup events, is the
        one set by MSTP(81) etc. This means that implicitly also
        high-pT jets can be generated in the pileup events.
    = 2 : low-p_T + double diffractive processes (ISUB = 95 and 93).
    = 3 : low-p_T + double diffractive + single diffractive processes
        (ISUB = 95, 93 and 92).
    = 4 : low-p_T + double diffractive + single diffractive + elastic
        processes, together corresponding to the full hadron-hadron
        cross-section (ISUB = 95, 93, 92 and 91).
MSTP(133) : (D=0) multiplicity distribution of pileup events.
    = 0 : selected by user, before each PYEVNT call, by giving the
        MSTP(134) value.
    = 1 : a Poissonian multiplicity distribution in the total number
        of pileup events. This is the relevant distribution if the
        switches set for the first event in PYSUBS give the same
        subprocesses as are implied by MSTP(132). In that case the
        mean number of events per beam crossing is nbar =
        (sum of cross-section for allowed processes) * PARP(131).
        Since bunch crossing which do not give any events at all
        (probability exp(-nbar)) are not simulated, the actual average
        number per PYEVNT call is <n> = nbar/(1-exp(-nbar)).
    = 2 : a biased distribution, as is relevant when one of the
        events to be generated is assumed to belong to an event class
        with a cross-section much smaller than the total hadronic
        cross-section. If sigma_hard is the cross-section for this
        rare process (or the sum of the cross-sections of several rare
        processes) and sigma_soft the cross-section for the processes
        allowed by MSTP(132), then define nbar = sigma_soft * PARP(131)
        and f = sigma_hard/sigma_soft. The probability that a bunch
        crossing will give i events is then
        Prob(i) = exp(-nbar) * nbar^i/i! * f * i
        i.e. the naive Poissonian is suppressed by a factor f since
        one of the events will be rare rather than frequent, but
        enhanced by a factor i since any of the i events may be the
        rare one. Only beam crossings which give at least one event
        of the required rare type are simulated, and the distribution
        above normalized accordingly.
    Note: for practical reasons, it is required that nbar < 120, i.e.
        that an average beam crossing does not contain more than 120
        pileup events. The multiplicity distribution is truncated
        above 200, or when the probability for a multiplicity has fallen
        below 10^(-6), whichever occurs sooner. Also low multiplicities
        with probabilities below 10^(-6) are truncated. See also
        PARI(91) - PARI(93).
MSTP(134) : (D=1) a user selected multiplicity, i.e. total number
    of pileup events, to be generated in the next PYEVNT call.
    May be reset for each new event, but must be in the range
    1 <= MSTP(134) <= 200.
MSTP(141) : (D=0) calling of PYKCUT in event generation chain, for
    inclusion of user-specified cuts.
    = 0 : not called.
    = 1 : called.
MSTP(142) : (D=0) calling of PYEVWT in the event generation chain,
    either to give weighted events or to modify standard cross-sections.
    See PYEVWT description in section 2.1 for further details.
    = 0 : not called.
    = 1 : called; the distribution of events among subprocesses and
        in kinematics variables is modified by the factor WTXS, set
        by the user in the PYEVWT call, but events come with a
        compensating weight PARI(10)=1./WTXS, such that total
        cross-sections are unchanged.
    = 2 : called; the cross-section itself is modified by the factor
        WTXS, set by the user in the PYEVWT call.
MSTP(181) : (R) PYTHIA version number.
MSTP(182) : (R) PYTHIA subversion number.
MSTP(183) : (R) last year of change for PYTHIA.
MSTP(184) : (R) last month of change for PYTHIA.
MSTP(185) : (R) last day of change for PYTHIA.

PARP(1) : (D=0.25 GeV) nominal Lambda-QCD used in running alpha-strong
    for hard scattering (see MSTP(3)).
PARP(2) : (D=10. GeV) lowest CM energy for the event-as-a-whole that
    the program will accept to simulate.
PARP(11), PARP(12) : (D=0.,-0.67) coefficients of a beamstrahlung-
    induced effective electron structure function, see MSTP(11) = 2.
    To ensure proper normalization, it is necessary to have
    PARP(11) and PARP(12) both > -1; in practice, values smaller
    than roughly -0.95 are not recommended.
PARP(13) : (D=25. GeV^2) Q_max^2 scale, to be set by user for defining
    maximum scale allowed for photoproduction when using the option
    MSTP(13)=2.
PARP(14) : (D=0.01) in the numerical integration of quark and gluon
    structure functions inside an electron, the successive halvings
    of evaluation point spacing is interrupted when two values agree
    in relative size, abs(new-old)/(new+old), to better than PARP(14).
    There are hardwired lower and upper limits of 2 and 8 halvings,
    respectively.
PARP(31) : (D=1.5) common K factor multiplying the differential
    cross-section for hard parton-parton processes when MSTP(33) = 1
    or 2, with the exception of colour annihilation graphs in the
    latter case.
PARP(32) : (D=2.0) special K-factor multiplying the differential
    cross-section in hard colour annihilation graphs, including
    resonance production, when MSTP(33) = 2.
PARP(33) : (D=0.075) this factor is used to multiply the ordinary Q^2
    scale in alpha_strong at the hard interaction for MSTP(33) = 3.
    The effective K-factor thus obtained is in accordance with the
    results in [Ell86].
PARP(35) : (D=0.20) fix alpha_strong value used in heavy flavour
    threshold factor when MSTP(35) = 1.
PARP(36) : (D=0. GeV) the width Gamma_Q for the heavy flavour studied
    in processes ISUB = 81 or 82; to be used for the threshold factor
    when MSTP(35) = 2.
PARP(37) : (D=2.) for MSTP(37)=1 this regulates the point at which the
    naive fix quark mass in Higgs couplings is assumed defined;
    specifically the running quark mass is assumed to coincide with
    fix one at PARP(37) times the fix quark mass, i.e.
    m_running(PARP(37)*m_fix) = m_fix.
PARP(38) : (D=0.70 GeV^3) the squared wave function at the origin,
    abs(R(0))^2, of the J/Psi wave function. Used for process 86.
    See ref. Glo88.
PARP(39) : (D=0.006 GeV^3) the squared derivative of the wave function
    at the origin, abs(R'(0))^2/m^2, of the chi_c wave functions.
    Used for the processes 87, 88 and 89. See ref. Glo88.
PARP(41) : (D=0.020 GeV) in the process of generating mass for
    resonances, and optionally to force that mass to be in a given
    range, only resonances with a total width in excess of PARP(41)
    are generated according to a Breit-Wigner shape (if allowed by
    MSTP(42)), while more narrow resonances are put on mass-shell.
PARP(42) : (D=2. GeV) minimum mass of resonances assumed allowed when
    evaluating total width of H0 to Z0 + Z0 or W+ + W- for cases
    when the H0 is so light that (at least) one Z/W is forced to be
    off mass-shell. Also generally used as safety check on minimum
    mass of resonance.
PARP(43) : (D=0.10) precision parameter used in numerical integration
    of width into channel with at least one daughter off mass-shell.
PARP(44) : (D=1000.) the nu papameter of the stronly interacting Z/W
    model of Dobado, Herrero and Terron [DHT90].
PARP(45) : (D=2054. GeV) the effective techni-rho mass parameter of the
    strongly interacting model of Dobado, Herrero and Terron [DHT90];
    see MSTP(46)=5. On physical grounds it should not be chosen smaller
    than about 1 TeV or larger than the about default value.
PARP(51) : (D=1.) if structure functions for light flavours have to be
    extrapolated to x values lower than covered by the parametrizations,
    an x^(-PARP(51)) behaviour is assumed in that region. This option
    only applies for the EHLQ and GRV proton structure functions that
    are internal to PYTHIA.
PARP(61) : (D=0.25 GeV) (C) Lambda value used in spacelike parton
    shower (see MSTP(64)). This value may be overwritten, see
    MSTP(3).
PARP(62) : (D=1. GeV) (C) effective cutoff Q or k_T value (see
    MSTP(64)), below which spacelike parton showers are not evolved.
PARP(63) : (D=0.25) (C) in spacelike shower evolution the virtuality
    Q^2 of a parton is multiplied by PARP(63) for use as a scale in
    alpha_strong and structure functions when MSTP(64) = 1.
PARP(64) : (D=1.) C) in spacelike parton shower evolution the transverse
    momentum-squared evolution scale k_T^2 is multiplied by PARP(64)
    for use as a scale in alpha_strong and structure functions when
    MSTP(64) = 2.
PARP(65) : (D=2. GeV) (C) effective minimum energy (in CM frame) of
    timelike or on-shell parton emitted in spacelike shower; see also
    PARP(66).
PARP(66) : (D=0.001) (C) effective lower cutoff on 1-z in spacelike
    showers, in addition to the cut implied by PARP(65).
PARP(67) : (D=4.) (C) the Q^2 scale of the hard scattering (see
    MSTP(32)) is multiplied by PARP(67) to define the maximum parton
    virtuality allowed in spacelike showers. This does not apply to
    s-channel resonances, where the maximum virtuality is set by m^2.
PARP(68) : (D=1E-6) lower Q^2 cutoff for QED spacelike showers.
PARP(71) : (D=4.) (C) the Q^2 scale of the hard scattering (see
    MSTP(32)) is multiplied by PARP(71) to define the maximum parton
    virtuality allowed in timelike showers. This does not apply to
    s-channel resonances, where the maximum virtuality is set by m^2.
PARP(81) : (D=1.45 GeV/c) effective minimum transverse momentum p^T_min
    for multiple interactions with MSTP(82) = 1.
    Note: The default value is changed compared to the one used up till
        version 5.3 (1.60 GeV/c), not because of any physics change,
        but simply because the treatment of alpha_strong and Lambda
        at flavour thresholds has been improved, leading to a smaller
        alpha_strong at small p_T values, which has to be compensated.
PARP(82) : (D=1.70 GeV) regularization scale p^T_0 of the transverse
    momentum spectrum for multiple interactions with MSTP(82) >= 2.
    Note: The default value is changed compared to the one used up till
        version 5.3 (2.00 GeV/c), not because of any physics change,
        but simply because the treatment of alpha_strong and Lambda
        at flavour thresholds has been improved, leading to a smaller
        alpha_strong at small p_T values, which has to be compensated.
PARP(83), PARP(84) : (D=0.5, 0.2) parameters of an assumed double
    Gaussian matter distribution inside the colliding hadrons for
    MSTP(82) = 4, of the form
    (1-PARP(83))*exp(-r^2) + (PARP(83)/PARP(84)^3)*exp(-r^2/PARP(84)^2)
    i.e. with a core of radius PARP(84) of the main radius and
    containing a fraction PARP(83) of the total hadronic matter.
PARP(85) : (D=0.33) probability that an additional interaction in the
    multiple interaction formalism gives two gluons, with colour
    connections to 'nearest neighbours' in momentum space.
PARP(86) : (D=0.66) probability that an additional interaction in the
    multiple interaction formalism gives two gluons, either as
    described in PARP(85) or as a closed gluon loop. Remaining fraction
    is supposed to consist of quark-antiquark pairs.
PARP(87), PARP(88) : (D=0.7, 0.5) in order to account for an assumed
    dominance of valence quarks at low transverse momentum scales, a
    probability is introduced that a gg-scattering according to naive
    cross-section is replaced by a qq~ one; this is used only for
    MSTP(82) >= 2. The probability is parametrized as
    prob = PARP(87)*(1.-(p_T^2/(p_T^2+PARP(88)*(PARP(82))^2))^2).
PARP(91) : (D=0.44 GeV/c) (C) width of Gaussian primordial k_T
    distribution inside hadron for MSTP(91) = 1.
PARP(92) : (D=0.44 GeV/c) (C) width of exponential primordial k_T
    distribution inside hadron for MSTP(91) = 2.
PARP(93) : (D=2. GeV/c) (C) upper cutoff for primordial k_T
    distribution inside hadron.
PARP(94) : (D=1.) (C) for MSTP(92) >= 2 this gives the value of the
    parameter k for the case when a pion remnant is split into two
    fragments (which is which is chosen at random).
PARP(95) : (D=0.) (C) for MSTP(92) >= 2 this gives the value of the
    parameter k for the case when a pion remnant is split into a meson
    and a spectator fragment jet, with chi giving the energy fraction
    taken by the meson.
PARP(96) : (D=3.) (C) for MSTP(92) >= 2 this gives the value of the
    parameter k for the case when a nucleon remnant is split into a
    diquark and a quark fragment, with chi giving the energy fraction
    taken by the quark jet.
PARP(97) : (D=1.) (C) for MSTP(92) >= 2 this gives the value of the
    parameter k for the case when a nucleon remnant is split into a
    baryon and a quark jet or a meson and a diquark jet, with chi giving
    the energy fraction taken by the quark jet or meson, respectively.
PARP(98) : (D=0.75) (C) for MSTP(92) = 5 this gives the power of an
    assumed basic 1/chi^PARP(98) behaviour in the splitting
    distribution.
PARP(99) : (D=0.44 GeV/c) (C) width of primordial k_T distribution
    inside photon; exact meaning depends on MSTP(93) value chosen.
PARP(100) : (D=2. GeV/c) (C) upper cutoff for primordial k_T
    distribution inside photon.
PARP(101) : (D=-0.02 GeV^2) minimum value of t-hat in (diffractive and)
    elastic scattering.
PARP(111) : (D=2. GeV) used to define the minimum invariant mass of
    the remnant hadronic system (i.e. when interacting partons have been
    taken away), together with original hadron masses and extra parton
    masses.
PARP(121) : (D=1.) the maxima obtained at initial maximization are
    multiplied by this factor if PARP(121)=1; typically PARP(121) would
    be given as the product of the violation factors observed (i.e the
    ratio of final maximum value to the initial maximum value) for the
    given process(es).
PARP(122) : (D=0.4) fraction of total probability that is shared
    democratically between the COEF coefficients open for the given
    variable, with remaining fraction distributed according to the
    optimization results of PYMAXI.
PARP(131) : (D=0.01 mb^(-1)) in the pileup events scenario (see
    MSTP(131) - MSTP(133)), PARP(131) gives the assumed luminosity per
    bunch-bunch crossing, i.e. if a subprocess has a cross-section
    sigma, the average number of events of this type per bunch-bunch
    crossing is nbar = sigma * PARP(131). PARP(131) may be obtained
    by dividing the integrated luminosity over a given time (1 s, say)
    by the number of bunch-bunch crossings that this corresponds to.
    Since the program will not generate more than 200 pileup events,
    the initialization procedure will crash if nbar is above 120.
______________________________________________________________________

2.6.  General Event Information

When an event is generated with PYEVNT, some information on this event
is stored in the MSTI and PARI arrays of the LUDATA commonblock (often
copied directly from the internal MINT and VINT variables). Further
information is stored in the complete event record; see section 2.7.

Part of the information is only relevant for some subprocesses; by
default everything irrelevant is set 0. Kindly note that, like the CKIN
constraints described in section 2.4, kinematical variables normally
(where it is not explicitly stated otherwise) refer to the naive hard
scattering, before initial and final state radiation effects have been
included.

COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
Purpose: to provide information on latest event generated or, in a few
    cases, on the accumulated statistics during the run.

MSTI(1) : specifies the general type of subprocess that has occured,
    according to the ISUB code given in section 2.2.
MSTI(2) : whenever MSTI(1) (together with MSTI(15) and MSTI(16)) are
    not enough to specify the type of process uniquely, MSTI(2)
    provides an ordering of the different possibilities. This is
    particularly relevant for the different colour flow topologies
    possible in QCD 2 -> 2 processes. With i = MSTI(15), j = MSTI(16)
    and k = MSTI(2), the QCD possibilities are, in the classification
    scheme of [Ben84]:
    ISUB = 11, i = j, q_i q_i -> q_i q_i;
        k = 1 : colour configuration A.
        k = 2 : colour configuration B.
    ISUB = 11, i /= j, q_i q_j -> q_i q_j;
        k = 1 : only possibility.
    ISUB = 12, q_i q~_i -> q_l q~_l;
        k = 1 : only possibility.
    ISUB = 13, q_i q~_i -> g g;
        k = 1 : colour configuration A.
        k = 2 : colour configuration B.
    ISUB = 28, q_i g -> q_i g;
        k = 1 : colour configuration A.
        k = 2 : colour configuration B.
    ISUB = 53, g g -> q_l q~_l;
        k = 1 : colour configuration A.
        k = 2 : colour configuration B.
    ISUB = 68, g g -> g g;
        k = 1 : colour configuration A.
        k = 2 : colour configuration B.
        k = 3 : colour configuration C.
    ISUB = 83, f q -> f' Q (by t-channel W exchange; does not
        distingusih colour flows but result of user selection);
        k = 1 : heavy flavour Q is produced on side 1.
        k = 2 : heavy flavour Q is produced on side 2.
MSTI(3) : number of partons produced in the hard interactions, i.e.
    the number n of the 2 -> n matrix elements used; is sometimes 3 or
    4 when a basic 2 -> 1 or 2 -> 2 process has been folded with
    two 1 -> 2 initial branchings (like q q' -> q" q'" H0).
MSTI(4) : number of documentation lines in beginning of common block
    LUJETS that are given with K(I,1) = 21; 0 for MSTP(125) = 0.
MSTI(5) : number of events generated to date in current run.
MSTI(6) : current frame of event, cf. MSTP(124).
MSTI(7), MSTI(8) : line number for documentation of outgoing partons/
    particles from hard scattering for 2 -> 2 or 2 -> 1 -> 2 processes
    (else = 0).
MSTI(10) : is 1 if cross-section maximum was violated in current
    event, and 0 if not.
MSTI(11) : KF flavour code for beam (side 1) particle.
MSTI(12) : KF flavour code for target (side 2) particle.
MSTI(13), MSTI(14) : KF flavour codes for side 1 and side 2 initial
    state shower initiators.
MSTI(15), MSTI(16) : KF flavour codes for side 1 and side 2 incoming
    partons to the hard interaction.
MSTI(17), MSTI(18) : flag to signal if particle on side 1 or side 2
    has been scattered diffractively; 0 if no, 1 if yes.
MSTI(21) - MSTI(24) : KF flavour codes for outgoing partons from the
    hard interaction. The number of positions actually used is
    process-dependent, see MSTI(3); trailing positions not used are
    set = 0.
MSTI(25), MSTI(26) : KF flavour codes of the products in the decay
    of a single s-channel resonance formed in the hard interaction.
    Are thus only used when MSTI(3) = 1 and the resonance is allowed
    to decay.
MSTI(31) : number of hard or semihard scatterings that occured in
    current event in the multiple interaction scenario; is = 0 for a
    low-p_T event.
MSTI(41) : the number of pileup events generated in latest PYEVNT
    call (including the first, 'hard' event).
MSTI(42) - MSTI(50) : ISUB codes for the events 2 - 10 generated in
    the pileup events scenario. The first event ISUB code is stored
    in MSTI(1). If MSTI(41) is less than 10, only as many positions are
    filled as there are pileup events. If MSTI(41) is above 10, some
    ISUB codes will not appear anywhere.

PARI(1) : total integrated cross-section for the processes under study,
    in mb. This number is obtained as a by-product of the selection of
    hard process kinematics, and is thus known with better accuracy
    when more event have been generated. The value stored here is based
    on all events up till the latest one generated.
PARI(2) : is the ratio PARI(1)/MSTI(5), i.e. the ratio of total
    integrated cross-section and number of events generated.
    Histograms generated with unit weight for events have to be
    multiplied by this factor, at the end of the run, to convert
    results to mb.
    For MSTP(142)=1, MSTI(5) is replaced by the sum of PARI(10)
    values. Histograms are then filled with weight PARI(10) for
    each event and multiplied by PARI(2) at the end.
PARI(9) : is weight WTXS returned from PYEVWT call when MSTP(142) >= 1,
    else is 1.
PARI(10) : is compensating weight 1./WTXS that should be associated
    to events when MSTP(142) = 1, else is 1.
PARI(11) : E_CM, i.e. total CM energy.
PARI(12) : s, i.e. total CM energy-squared.
PARI(13) : m-hat = sqrt(s-hat), i.e. mass of the hard scattering
    subsystem.
PARI(14) : s-hat of the hard subprocess (2 -> 2 or 2 -> 1).
PARI(15) : t-hat of the hard subprocess (2 -> 2 or 2 -> 1 -> 2).
PARI(16) : u-hat of the hard subprocess (2 -> 2 or 2 -> 1 -> 2).
PARI(17) : p_T-hat of the hard subprocess (2 -> 2 or 2 -> 1 -> 2),
    evaluated in the rest frame of the hard interaction.
PARI(18) : p_T-hat^2 of the hard subprocess; see PARI(17).
PARI(19) : m'-hat, the mass of the complete three- or four-body final
    state in 2 -> 3 or 2 -> 4 processes (while m-hat, given in PARI(13),
    here corresponds to the one- or two-body central system).
    Kinematically m-hat <= m'-hat <= E_CM.
PARI(20) : s'-hat = m'-hat^2; see PARI(19).
PARI(21) : Q of the hard subprocess. The exact definition is
    process-dependent, see MSTP(32).
PARI(22) : Q^2 of the hard subprocess; see PARI(21).
PARI(23) : Q of the outer hard scattering subprocess. Agrees with
    PARI(21) for a 2 -> 1 or 2 -> 2 process. For a 2 -> 3 or 2 -> 4
    W/Z fusion process, it is set by the W/Z mass scale, and for
    subprocesses 121 and 122 by the heavy quark mass.
PARI(24) : Q^2 of the outer hard scattering subprocess; see PARI(23).
PARI(25) : Q scale used as maximum virtuality in parton showers.
    Is equal to PARI(23), except for Deep Inelastic Scattering
    processes when MSTP(22) > 0.
PARI(26) : Q^2 scale in parton showers; see PARI(25).
PARI(31), PARI(32) : the momentum fractions x of the initial
    state parton shower initiators on side 1 and 2, respectively.
PARI(33), PARI(34) : the momentum fractions x taken by the partons
   at the hard interaction, as used e.g. in the structure functions.
PARI(35) : Feynman-x, x_F = x_1 - x_2 = PARI(33) - PARI(34).
PARI(36) : tau = s-hat/s = x_1*x_2 = PARI(33)*PARI(34).
PARI(37) : y* = (1/2) * log(x_1/x_2), i.e. rapidity of the hard
   interaction subsystem in the CM frame of the event-as-a-whole.
PARI(38)  : tau' = s'-hat/s = PARI(20)/PARI(12).
PARI(39), PARI(40) : the primordial k_T values selected in the two
    beam remnants.
PARI(41) : cos(theta-hat), where theta-hat is the scattering angle of
   a 2 -> 2 (or 2 -> 1 -> 2) interaction, defined in the rest frame
   of the hard scattering subsystem.
PARI(42) : x_T, i.e. scaled transverse momentum of the hard
   scattering subprocess, x_T = 2 * p_T-hat/E_CM.
PARI(43), PARI(44) : x_L3 and x_L4, i.e. longitudinal momentum
   fractions of the two scattered partons, in the range -1 < x_L < 1,
   in the CM frame of the event-as-a-whole.
PARI(45), PARI(46) : x_3 and x_4, i.e. scaled energy fractions of the
   two scattered partons, in the CM frame of the event-as-a-whole.
PARI(47), PARI(48) : y*_3 and y*_4, i.e. rapidities of the two
   scattered partons in the CM frame of the event-as-a-whole.
PARI(49), PARI(50) : eta*_3 and eta*_4, i.e. pseudorapidities of the
   two scattered partons in the CM frame of the event-as-a-whole.
PARI(51), PARI(52) : cos(theta*_3) and cos(theta*_4), i.e. cosines of
   the polar angles of the two scattered partons in the CM frame of
   the event-as-a-whole.
PARI(53), PARI(54) : theta*_3 and theta*_4, i.e. polar angles of the
   two scattered partons, defined in the range 0 < theta* < pi, in
   the CM frame of the event-as-a-whole.
PARI(55), PARI(56) : azimuthal angles phi*_3 and phi*_4 of the two
   scattered partons, defined in the range -pi < phi* < pi, in the
   CM frame of the event-as-a-whole.
PARI(61) : multiple interaction enhancement factor for current event.
    A large value corresponds to a central collision and a small value
    to a peripheral one.
PARI(65) : sum of the transverse momenta of partons generated at
    the hardest interaction of the event, excluding initial and
    final state radiation, i.e. 2 * PARI(17).
PARI(66) : sum of the transverse momenta of all partons generated at
    the hardest interaction, including initial and final state
    radiation, resonance decay products, and primordial k_T.
PARI(67) : sum of transverse momenta of partons generated at hard
    interactions, excluding the hardest one (see PARI(65)), and also
    excluding initial and final state radiation. Is nonvanishing only
    in the multiple interaction scenario.
PARI(68) : sum of transverse momenta of all partons generated at hard
    interactions, excluding the hardest one (see PARI(66)), but
    including initial and final state radiation. Is nonvanishing only
    in the multiple interaction scenario.
PARI(69) : sum of transverse momenta of all partons generated in hard
    interactions (PARI(66) + PARI(68)) and, additionally, of all beam
    remnant partons.
PARI(71), PARI(72) : sum of the momentum fractions x taken by initial
    state parton shower initiators on side 1 and and side 2, excluding
    those of the hardest interaction. Is nonvanishing only in the
    multiple interaction scenario.
PARI(73), PARI(74) : sum of the momentum fractions x taken by the
    partons at the hard interaction on side 1 and side 2, excluding
    those of the hardest interaction. Is nonvanishing only in the
    multiple interaction scenario.
PARI(75), PARI(76) : the x value of a photon that branches into
    quarks or gluons, i.e. x at interface between initial state QED
    and QCD cascades.
PARI(77), PARI(78) : the chi values selected for beam remnants that
    are split into two objects, describing how the energy is shared
    (see MSTP(92) etc.); is 0. if no splitting is needed.
PARI(81) : size of the threshold factor (enhancement or suppression)
    in the latest event with heavy flavour production; see MSTP(35).
PARI(91) : average multiplicity nbar of pileup events, defined as
    (cross-section for processes allowed by MSTP(132)) * PARP(131).
    Only relevant for MSTP(133) = 1 or 2.
PARI(92) : average multiplicity of pileup events as actually
    simulated, i.e. with multiplicity = 0 events removed and the
    high-end tail truncated. Only relevant for MSTP(133) = 1 or 2.
PARI(93) : for MSTP(133) = 1 it is the probability that a beam crossing
    will produce a pileup event at all, i.e. that there will be at
    least one hadron-hadron interaction; for MSTP(133) = 2 the
    probability that a beam crossing will produce a pileup event with
    one hadron-hadron interaction of the desired rare type.
______________________________________________________________________

2.7.  The Event Record

When an event is generated by a PYEVNT call, it is stored in the
common block LUJETS. Here each parton and particle is represented by
one line of information, giving a status code, parton/particle code,
line of mother, lines of first and last daughter (or colour
flow information in hard interaction and parton showers),
momentum, energy, mass, production vertex and time, and invariant
lifetime (if unstable). The LUJETS common block is used extensively
both by the PYTHIA and the JETSET routines; indeed it provides the
bridge that allows the general utility routines in JETSET to be used
also for PYTHIA events. A detailed description of LUJETS is found in
the JETSET 7.3 manual, a companion file to the program. The PYTHIA
event listing begins (optionally) with a few lines of event summary,
specific to the hard process simulated and thus not described in the
JETSET manual. Therefore we here give a brief summary on the general
structure of the LUJETS common block, followed by details specific to
PYTHIA.

COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)

N : number of entries, i.e. number of lines used in the K, P and V
    matrices, in the current event.

K(I,1) : status code KS for the parton/particle stored in the line.
    = 0 : empty line.
    = 1 - 10 : undecayed particle or unfragmented parton.
    = 11 - 20 : decayed particle or fragmented or showered parton.
    = 21 : documentation lines for compressed event summary (see below).
K(I,2) : KF flavour code for partons and particles, following the 1988
    Particle Data Group numbering conventions [PDG88]. Some of the most
    frequently used ones are
       1  d            11  e-           21  g            31
       2  u            12  nu_e         22  gamma        32  Z'0
       3  s            13  mu-          23  Z0           33
       4  c            14  nu_mu        24  W+           34  W'+
       5  b            15  tau-         25  H0           35  H'0
       6  t            16  nu_tau       26               36  A0
       7  l            17  chi-         27               37  H+
       8  h            18  nu_chi       28               38
       9               19               29               39  LQ
      10               20               30               40  R0
    1103  dd_1       2103  ud_1       3101  sd_0       3101  su_0
    2101  ud_0       2203  uu_1       3103  sd_1       3103  su_1
     211  pi+         213  rho+       2112  n          1114  Delta-
     311  K0          313  K*0        2212  p          2114  Delta0
     321  K+          323  K*+        3122  Lambda0    2214  Delta+
     411  D+          413  D*+        3112  Sigma-     2224  Delta++
     421  D0          423  D*0        3212  Sigma0     3114  Sigma*-
     431  D_s+        433  D*_s+      3222  Sigma+     3214  Sigma*0
     130  K_L0                        3312  Xi-        3224  Sigma*+
     310  K_S0                        3322  Xi0        3334  Omega-
    A negative KF code, where existing, always corresponds to the
    antiparticle of the one listed above.
    To designate diffractive states, three non-standard codes are made
    available by JETSET.
     210  pi_diffr+      2110  n_diffr       2210  p_diffr
K(I,3) : line number of parent parton or jet, where known, else 0.
K(I,4) : normally line number of the first daughter; for K(I,1) = 3, 13
    or 14 instead special colour flow information.
K(I,5) : normally line number of the last daughter; for K(I,1) = 3, 13
    or 14 instead special colour flow information.

P(I,1) : p_x, momentum in the x direction (in GeV/c).
P(I,2) : p_y, momentum in the y direction (in GeV/c).
P(I,3) : p_z, momentum in the z direction (in GeV/c).
P(I,4) : E, energy (in GeV).
P(I,5) : m, mass (in GeV/c^2). Spacelike partons are
    given with P(I,5) = -sqrt(-m^2).

V(I,1) : x position of production vertex (in mm).
V(I,2) : y position of production vertex (in mm).
V(I,3) : z position of production vertex (in mm).
V(I,4) : time of production (in mm/c = 3.33*10^-12 s).
V(I,5) : proper lifetime of the particle/parton (in mm/c = 3.33*10^-12
    s); is 0 for a stable entry.

In most instances, only the actual partons and particles produced are of
interest. For MSTP(125) = 0, the event record starts off with the parton
configuration existing after hard interaction, initial and final state
radiation, multiple interactions and beam remnants have been considered.
The partons are arranged in colour singlet clusters, ordered as required
for string fragmentation. Also photons and leptons produced as part of
the hard interaction (e.g. from q q~ -> g gamma or u u~ -> Z0 -> e+ e-)
appear in this part of the event record. These original entries
appear with pointer K(I,3) = 0, whereas the products of the subsequent
fragmentation and decay have K(I,3) numbers pointing back to the line
of the parent.

The standard documentation, obtained with MSTP(125) = 1, includes a few
lines in the beginning of the event record, which contain a brief
summary of the process that has taken place. The number of lines used
depends on the nature of the hard process, and is stored in MSTI(4) for
the current event. These lines all have K(I,1) = 21. For all
processes, lines 1 and 2 give the two incoming hadrons. When listed with
LULIST, these two lines will be separated from subsequent lines by a
sequence of ====== signs, to improve readability. For diffractive and
elastic events, the two outgoing states in lines 3 and 4 completes the
list. Otherwise, lines 3 and 4 contain the two partons that initiate
the two initial state parton showers, and 5 and 6 the endproducts of
these showers, i.e. the partons that enter the hard interaction. With
initial state radiation switched off, lines 3 and 5 and lines 4 and 6
coincide. For a simple 2 -> 2 hard scattering, lines 7 and 8 give the
two outgoing partons/particles from the hard  interaction, before any
final state radiation. For 2 -> 2 processes proceeding via an
intermediate resonance such as Z0/gamma*, W+/-, H0 or R, the resonance
is found in 7 and the two outgoing partons/particles in 8 and 9. In some
cases one or both of these may be a resonance in its own right, so that
further pairs of lines are added for subsequent decays. If the decay of
a given resonance has been switched off, then no decay products are
listed either in this initial summary or in the subsequent ordinary
listing. Whenever partons are listed, they are assumed on mass shell for
simplicity. The fact that effective masses may be generated by initial
and final state radiation is taken into account in the actual parton
configuration that is allowed to fragment, however. A special case is
provided by W+ W- or Z0 Z0 fusion to a H0. Then the virtual W:s or Z:s
are shown in lines 7 and 8, the H0 in line 9 and the two recoiling
quarks (that emitted the bosons) in 10 and 11, followed by the Higgs
decay products. Since the W:s and Z:s are spacelike, what is actually
listed as the mass for them is -sqrt(-m^2). The listing of the event
documentation closes with another line made up of ====== signs.

A few examples may help clarify the picture. For a single diffractive
event p + p~ -> p_diffr + p~, the event record will start with
 I K(I,1)   K(I,2) K(I,3)   comment
 1     21     2212      0   incoming p
 2     21    -2212      0   incoming p~
=========================   not part of record; appears in listings
 3     21       27      1   outgoing p_diffr
 4     21    -2212      2   outgoing p~
=========================

The typical QCD 2 -> 2 process would be
 I K(I,1)   K(I,2) K(I,3)   comment
 1     21     2212      0   incoming p
 2     21    -2212      0   incoming p~
=========================
 3     21        2      1   u picked from incoming p
 4     21       -1      2   d~ picked from incoming p~
 5     21       21      3   u evolved to g at hard scattering
 6     21       -1      4   still d~ at hard scattering
 7     21       21      0   outgoing g from hard scattering
 8     21       -1      0   outgoing d~ from hard scattering
=========================
Note that, where well defined, the K(I,3) code does contain information
on which side the different partons come from, e.g. above the gluon in
line 5 points back to the u in line 3, which points back to the proton
in line 1. In the example above it would have been possible to
associate the scattered g in line 7 with the incoming one in line 5, but
this is not possible in the general case, consider e.g. g g -> g g.

As a final example, W+ W- fusion to a H0 might look like
 I K(I,1)   K(I,2) K(I,3)   comment
 1     21     2212      0   first incoming p
 2     21     2212      0   second incoming p
=========================
 3     21        2      1   u picked from first p
 4     21       21      2   g picked from second p
 5     21        2      3   still u after initial state radiation
 6     21       -4      4   g evolved to c~
 7     21       24      5   spacelike W+ emitted by u quark
 8     21      -24      6   spacelike W- emitted by c~ quark
 9     21       25      0   Higgs produced by W+ W- fusion
10     21        1      5   u turned into d by emission of W+
11     21       -3      6   c~ turned into s~ by emission of W-
12     21       23      9   first Z0 coming from decay of H0
13     21       23      9   second Z0 coming from decay of H0
14     21       12     12   electron neutrino from first Z0 decay
15     21      -12     12   electron antineutrino from first Z0 decay
16     21        5     13   b quark from second Z0 decay
17     21       -5     13   b~ quark from second Z0 decay
=========================

After these lines with initial information, the event record looks the
same as for MSTP(125) = 0, i.e. first comes the parton configuration
to be fragmented and, after another separator line ====== in the output
(but not the event record), the products of subsequent fragmentation
and decay chains. The K(I,3) pointers for the partons, as well as
leptons and photons produced in the hard interaction, are now pointing
towards the documentation lines above, however. In particular, beam
remnants point to 1 or 2, depending on which side they belong
to, and partons emitted in the initial state parton showers point to 3
or 4. In the second example above, the partons produced by final state
radiation will be pointing back to 7 and 8; as usual, it should be
remembered that a specific assignment to 7 or 8 need not be unique. For
the third example, final state radiation partons will come both from
partons 10 and 11 and from partons 16 and 17, and additionally there
will be a neutrino-antineutrino pair pointing to 14 and 15. The extra
pairs of partons that are generated by multiple interactions do not
point back to anything, i.e. they have K(I,3) = 0.

There exists a third documentation option, MSTP(125) = 2. Here the
history of initial and final state parton branchings may be traced,
including all details on colour flow. This information has not been
optimized for user-friendliness, and can not be recommended for
general usage. With this option, the initial documentation lines
are the same. They are followed by blank lines, K(I,1) = 0, up to
line 20 (can be changed in MSTP(126)). From line 21 and onwards
each parton with K(I,1) = 3, 13 or 14 appear with special colour flow
information in the K(I,4) and K(I,5) positions; see the JETSET 7.3
manual. For an ordinary 2 -> 2 scattering, the two incoming partons
at the hard scattering are stored in line 21 and 22, and the two
outgoing in 23 and 24. The colour flow between these partons has
to be chosen according to the proper relative probabilities in
cases when many alternatives are possible, see [Ben84]. If there
is initial state radiation, the two partons in line 21 and 22 are
copied down to line 25 and 26, from which the initial state showers are
reconstructed backwards step by step. The branching history may be read
by noting that, for a branching a -> b + c, the K(I,3) codes of b and c
point towards the line number of a. Since the showers are reconstructed
backwards, this actually means that parton b would appear in the listing
before parton a and c, and hence have a pointer to a position below
itself in the list. Associated timelike partons c may initiate timelike
showers, as may the partons of the hard scattering. Again a showering
parton or pair of partons will be copied down towards the end of the
list and allowed to undergo successive branchings c -> d + e, with d and
e pointing towards c. The mass of timelike partons is properly stored in
P(I,5), for spacelike partons instead -sqrt(-m~2) is stored. After this
section containing all the branchings comes the final parton
configuration, properly arranged in colour, followed by all
subsequent fragmentation and decay products, as usual.
______________________________________________________________________

2.8.  Other Routines and Commonblocks

The subroutines and commonblocks that a user will come in direct
contact with have already been described. A number of other routines
and commonblocks exist, and are here briefly listed for the sake of
completeness.

SUBROUTINE PYINKI
Purpose: to initialize the kinematics given by the two incoming
    hadrons/leptons.

SUBROUTINE PYINRE
Purpose: to initialize the widths and effective widths of resonances.

SUBROUTINE PYXTOT
Purpose: to give the parametrized total, double diffractive, single
    diffractive and elastic cross-sections for different energies and
    colliding hadrons.

SUBROUTINE PYMAXI
Purpose: to find optimal coefficients COEF for the selection of
    kinematical variables, and to find the related maxima for
    the differential cross-section times Jacobian factors,
    for each of the subprocesses included.

SUBROUTINE PYPILE
Purpose: to determine the number of pileup events, i.e. events
    appearing in the same beam-beam crossing.

SUBROUTINE PYRAND
Purpose: to generate the quantities characterizing a hard scattering
    on the parton level, according to the relevant matrix elements.

SUBROUTINE PYSCAT
Purpose: to find outgoing flavours and to set up the kinematics and
    colour flow of the hard scattering.

SUBROUTINE PYSSPA(IP1,IP2)
Purpose: to generate the spacelike showers of the initial state
    radiation.

SUBROUTINE PYRESD
Purpose: to allow Z0, W+/-, H0, Z'0, H+/- and R0 resonances to decay,
    including chains of successive decays and parton showers.

SUBROUTINE PYMULT
Purpose: to generate semihard interactions according to the
    multiple interaction formalism.

SUBROUTINE PYREMN(IPU1,IPU2)
Purpose: to add on target remnants and include primordial k_T.

SUBROUTINE PYDIFF
Purpose: to handle diffractive and elastic scattering events.

SUBROUTINE PYDOCU
Purpose: to compute cross-sections of processes, based on current
    Monte Carlo statistics, and to store event information in the
    MSTI and PARI arrays.

SUBROUTINE PYWIDT
Purpose: to calculate widths and effective widths of resonances.

SUBROUTINE PYKLIM
Purpose: to calculate allowed kinematical limits.

SUBROUTINE PYKMAP
Purpose: to calculate the value of a kinematical variable when this
    is selected according to one of the simple pieces.

Purpose: to give the differential cross-section (multiplied by the
    relevant Jacobians) for a given subprocess and kinematical
    setup.

SUBROUTINE PYSTFU(KF,X,Q2,XPQ)
SUBROUTINE PYSTEL(X,Q2,XPEL)
SUBROUTINE PYSTGA(X,Q2,XPGA)
SUBROUTINE PYSTPI(X,Q2,XPPI)
SUBROUTINE PYSTPR(X,Q2,XPPR)
Purpose: to give gluon, quark and antiquark structure functions for
    given x and Q^2 values for e+-, gamma, p, p~, n, n~, pi+-, and
    hyperons.

SUBROUTINE PYSPLI
Purpose: to give hadron remnant or remnants left when the reacting
    parton is kicked out.

FUNCTION PYGAMM(X)
Purpose: to give the value of the ordinary gamma function Gamma(x)
    (used in some structure function parametrizations).

SUBROUTINE PYWAUX(IAUX,EPS,WRE,WIM)
SUBROUTINE PYI3AU(EPS,RAT,Y3RE,Y3IM)
Purpose: to calculate the value of some auxiliary functions appearing
    in the cross-section expressions in PYSIGH.

FUNCTION PYSPEN(XREIN,XIMIN,IREIM)
Purpose: to calculate the real and imaginary part of the Spence
    function.

BLOCK DATA PYDATA
Purpose: to give sensible default values to all status codes and
    parameters.

COMMON/PYINT1/MINT(400),VINT(400)
Purpose: to collect a host of integer and real valued variables used
    internally in the program during the initialization and/or event
    generation stage. These variables must not be changed by the user.

MINT(1) : specifies the general type of subprocess that has occured,
    according to the ISUB code given in section 2.2.
MINT(2) : whenever MINT(1) (together with MINT(15) and MINT(16)) are
    not enough to specify the type of process uniquely, MINT(2)
    provides an ordering of the different possibilities, see
    MSTI(2).
MINT(3) : number of partons produced in the hard interactions, i.e.
    the number n of the 2 -> n matrix elements used; is sometimes 3 or
    4 when a basic 2 -> 1 or 2 -> 2 process has been folded with
    two 1 -> 2 initial branchings (like q q' -> q" q'" H0).
MINT(4) : number of documentation lines in beginning of common block
    LUJETS that are given with K(I,1) = 21; 0 for MSTP(125) = 0.
MINT(5) : number of events generated to date in current run.
MINT(6) : current frame of event, cf. MSTP(124).
MINT(7), MINT(8) : line number for documentation of outgoing partons/
    particles from hard scattering for 2 -> 2 or 2 -> 1 -> 2 processes
    (else = 0).
MINT(10) : is 1 if cross-section maximum was violated in current
    event, and 0 if not.
MINT(11) : KF flavour code for beam (side 1) particle.
MINT(12) : KF flavour code for target (side 2) particle.
MINT(13), MINT(14) : KF flavour codes for side 1 and side 2 initial
    state shower initiators.
MINT(15), MINT(16) : KF flavour codes for side 1 and side 2 incoming
    partons to the hard interaction.
MINT(17), MINT(18) : flag to signal if particle on side 1 or side 2
    has been scattered diffractively; 0 if no, 1 if yes.
MINT(19), MINT(20) : flag to signal initial state structure with
    parton inside photon inside electron on side 1 or side 2;
    0 if no, 1 if yes.
MINT(21) - MINT(24) : KF flavour codes for outgoing partons from the
    hard interaction. The number of positions actually used is
    process-dependent, see MINT(3); trailing positions not used are
    set = 0.
MINT(25), MINT(26) : KF flavour codes of the products in the decay
    of a single s-channel resonance formed in the hard interaction.
    Are thus only used when MINT(3) = 1 and the resonance is allowed
    to decay.
MINT(31) : number of hard or semihard scatterings that occured in
    current event in the multiple interaction scenario; is = 0 for a
    low-p_T event.
MINT(35) : in a true 2 -> 3 process, where one particle is a resonance
    with decay channel selected already before the PYRESD call, the
    decay channel number (in the /LUDAT3/ numbering) is stored here.
MINT(41), MINT(42) : type of incoming beam or target particle;
    1 for lepton and 2 for hadron. A photon counts as a lepton if it
    is not resolved (MSTP(14)=0) and as a hadron if it is resolved
    (MSTP(14)=1).
MINT(43) : combination of incoming beam and target particles. A photon
    counts as a hadron.
    = 1 : lepton on lepton.
    = 2 : lepton on hadron.
    = 3 : hadron on lepton.
    = 4 : hadron on hadron.
MINT(44) : as MINT(43), but a photon counts as a lepton.
MINT(45), MINT(46) : structure of incoming beam and target particles.
    = 1 : no internal structure, i.e. electron or photon which carries
        full beam energy.
    = 2 : defined with structure functions which are not peaked at
        x = 1, i.e. hadrons and a resolved photon.
    = 3 : defined with structure functions which are peaked at x = 1,
        i.e. the electron.
MINT(47) : combination of incoming beam and target particle structure
    function types.
    = 1 : no structure function either for beam or target.
    = 2 : structure functions for target but not for beam.
    = 3 : structure functions for beam but not for target.
    = 4 : structure functions both for beam and target, but not
        both peaked at x = 1.
    = 5 : structure functions both for beam and target, with both
        peaked at x = 1.
MINT(48) : total number of subprocesses switched on.
MINT(49) : number of subprocesses that are switched on, apart from
    elastic scattering and single, double and central diffractive.
MINT(51) : internal flag that event failed cuts.
    = 0 : no problem.
    = 1 : event failed; new one to be generated.
MINT(52) : internal counter for number of lines used (in /LUJETS/)
    before multiple interactions are considered.
MINT(53) : internal counter for number of lines used (in /LUJETS/)
    before beam remnants are considered.
MINT(55) : the heaviest new flavour switched on for QCD processes,
    specifically the flavour to be generated for ISUB = 81, 82 or
    83 or 84.
MINT(56) : the heaviest new flavour switched on for QED processes,
    specifically for ISUB = 85. Note that, unlike MINT(55), the
    heaviest flavour may here be a lepton, and that heavy means the
    one with largest KF code
MINT(61) : internal switch for the mode of operation of resonance
    width calculations in PYWIDT for gamma*/Z0 or gamma*/Z0/Z'0.
    = 0 : without reference to initial state flavours.
    = 1 : with reference to given initial state flavours.
    = 2 : for given final state flavours.
MINT(62) : internal switch for use at initialization of H width.
    = 0 : use widths into Z + Z* or W + W* calculated before.
    = 1 : evaluate widths into Z + Z* or W + W* for current Higgs mass.
MINT(65) : internal switch to indicate initialization without
    specified reaction.
    = 0 : normal initialization.
    = 1 : initialization with argument 'none' in PYINIT call.
MINT(71) : switch whether current process is singular for p_T -> 0
    or not.
    = 0 : non-singular process, i.e. proceeding via an s-channel
        resonance or with both products having a mass above CKIN(6).
    = 1 : singular process.
MINT(72) : number of s-channel resonances which may contribute to
    the cross-section.
MINT(73) : KF code of first s-channel resonance; 0 if there is none.
MINT(74) : KF code of second s-channel resonance; 0 if there is none.
MINT(81) : number of pileup events selected.
MINT(82) : sequence number of currently considered pileup event.
MINT(83) : number of lines in the event record already filled by
    previously considered pileup events.
MINT(84) : MINT(83) + MSTP(126), i.e. number of lines already filled
    by previously considered events plus number of lines to be kept
    free for event documentation.

VINT(1) : ECM, CM energy.
VINT(2) : s (=ECM^2) mass-square of complete system.
VINT(3) : mass of beam particle.
VINT(4) : mass of target particle.
VINT(5) : momentum of beam (and target) particle in CM frame.
VINT(6) - VINT(10) : theta, phi and beta for rotation and boost
    from CM frame to user-specified frame.
VINT(11) : tau_min.
VINT(12) : y*_min.
VINT(13) : cos(theta-hat)_min for cos(theta-hat) <= 0.
VINT(14) : cos(theta-hat)_min for cos(theta-hat) >= 0.
VINT(15) : x_T^2_min.
VINT(16) : tau'_min.
VINT(21) : tau.
VINT(22) : y*.
VINT(23) : cos(theta-hat).
VINT(24) : phi* (azimuthal angle).
VINT(25) : x_T^2.
VINT(26) : tau'.
VINT(31) : tau_max.
VINT(32) : y*_max.
VINT(33) : cos(theta-hat)_max for cos(theta-hat) <= 0.
VINT(34) : cos(theta-hat)_max for cos(theta-hat) >= 0.
VINT(35) : x_T^2_max.
VINT(36) : tau'_max.
VINT(41), VINT(42) : the momentum fractions x taken by the partons
   at the hard interaction, as used e.g. in the structure functions.
VINT(43) : m-hat = sqrt(s-hat), mass of hard scattering subsystem.
VINT(44) : s-hat of the hard subprocess (2 -> 2 or 2 -> 1).
VINT(45) : t-hat of the hard subprocess (2 -> 2 or 2 -> 1 -> 2).
VINT(46) : u-hat of the hard subprocess (2 -> 2 or 2 -> 1 -> 2).
VINT(47) : p_T-hat of the hard subprocess (2 -> 2 or 2 -> 1 -> 2),
    i.e. transverse momentum evaluated in the rest frame of the
    scattering.
VINT(48) : p_T-hat^2 of the hard subprocess; see VINT(47).
VINT(49) : m'-hat, the mass of the complete three- or four-body final
    state in 2 -> 3 or 2 -> 4 processes.
VINT(50) : s'-hat = m'-hat^2; see VINT(49).
VINT(51) : Q of the hard subprocess. The exact definition is
    process-dependent, see MSTP(32).
VINT(52) : Q^2 of the hard subprocess; see VINT(51).
VINT(53) : Q of the outer hard scattering subprocess. Agrees with
    VINT(51) for a 2 -> 1 or 2 -> 2 process. For a 2 -> 3 or 2 -> 4
    W/Z fusion process, it is set by the W/Z mass scale, and for
    subprocesses 121 and 122 by the heavy quark mass.
VINT(54) : Q^2 of the outer hard scattering subprocess; see VINT(53).
VINT(55) : Q scale used as maximum virtuality in parton showers.
    Is equal to VINT(53), except for Deep Inelastic Scattering
    processes when MSTP(22) > 0.
VINT(56) : Q^2 scale in parton showers; see VINT(55).
VINT(57) : alpha_em value of hard process.
VINT(58) : alpha_strong value of hard process.
VINT(59) : sin(theta-hat) (cf. VINT(23)); used for improved numerical
    precision in elastic and diffractive scattering.
VINT(61), VINT(62) : nominal m^2 values, i.e. without initial state
    radiation effects, for the two partons entering the hard
    interaction.
VINT(63), VINT(64) : nominal m^2 values, i.e. without final state
    radiation effects, for the two (or one) partons/particles leaving
    the hard interaction.
VINT(65) : p-hat_init, i.e. common nominal absolute momentum of the two
    partons entering the hard interaction, in their rest frame.
VINT(66) : p-hat_fin, i.e. common nominal absolute momentum of the two
    partons leaving the hard interaction, in their rest frame.
VINT(71) : p_T_min of process, i.e. CKIN(3) or CKIN(5), depending on
    which is larger, and whether the process is singular in p_T -> 0
    or not.
VINT(73) : tau = m^2/s value of first resonance, if any; see MINT(73).
VINT(74) : m*Gamma/s value of first resonance, if any; see MINT(73).
VINT(75) : tau = m^2/s value of second resonance, if any; see MINT(74).
VINT(76) : m*Gamma/s value of second resonance, if any; see MINT(74).
VINT(80) : correction factor (evaluated in PYOFSH) for the cross-section
    of resonances produced in 2 -> 2 processes, if only some mass range
    of the full Breit-Wigner shape is allowed by user-set mass cuts
    (CKIN(2), CKIN(45) - CKIN(48)).
VINT(81) - VINT(84) : the cos(theta) and phi variables of a true 2 -> 3
    process, where one product is a resonance, effectively giving
    2 -> 4. The first two are cos(theta) and phi for the resonance
    decay, the other two ditto for the effective system formed by the
    other two particles.
VINT(85), VINT(86) : transverse momenta in a true 2 -> 3 process; one
    is stored in VINT(47) (that of the Z in g g -> Z Q Q), while the
    smaller of the two others is stored in VINT(85) and the larger
    in VINT(86).
VINT(91), VINT(92) : gives a dimensionless suppression factor, to take
    into account reduction in cross-section due to the allowed channels
    for a W+W+ or W-W- pair, respectively, in the same sense as
    WIDS(24,1) gives it for a W+W- pair.
VINT(98) : is sum of VINT(100) values for current run.
VINT(99) : is weight WTXS returned from PYEVWT call when MSTP(142) >= 1,
    else is 1.
VINT(100) : is compensating weight 1./WTXS that should be associated
    to events when MSTP(142) = 1, else is 1.
VINT(101) : total cross-section.
VINT(102) : elastic cross-section.
VINT(103) : single diffractive cross-section.
VINT(104) : double diffractive cross-section.
VINT(105) : central diffractive cross-section.
VINT(106) : total non-diffractive, inelastic cross-section.
VINT(108) : ratio of maximum differential cross-section observed
    to maximum differential cross-section assumed in generation;
    cf. MSTP(123).
VINT(109) : ratio of minimal (negative!) cross-section observed
    to maximum differential cross-section assumed in generation;
    could only become negative if cross-sections are incorrectly
    included.
VINT(111) - VINT(116) : for MINT(61) = 1 gives kinematical factors
    for the different pieces contributing to gamma*/Z0 or gamma*/Z0/Z'0
    production, for MINT(61) = 2 gives sum of final state weights
    for the same; coefficients are given in the order pure gamma*,
    gamma*-Z0 interference, gamma*-Z'0 interference, pure Z0, Z0-Z'0
    interference and pure Z'0.
VINT(117) : width of Z0; needed in gamma*/Z0/Z'0 production.
VINT(121) : nuclear slope parameter B in t-hat distribution for
    (diffractive and) elastic scattering.
VINT(122) : curvature parameter C in t-hat distribution for (diffractive
    and) elastic scattering.
VINT(131) : total cross-section (in mb) for subprocesses allowed in
    the pileup events scenario according to the MSTP(132) value.
VINT(132) : nbar = VINT(131) * PARP(131), cf. PARI(91).
VINT(133) : <n> = (sum_i Prob_i * i)/(sum_i Prob_i) as actually
    simulated, i.e. 1 <= i <= 100 (or smaller), cf. PARI(92).
VINT(134) : is exp(-nbar) * sum_i nbar^i/i! for MSTP(133) = 1 and
    exp(-nbar) * sum_i nbar^i/(i-1)! for MSTP(133) = 2, cf. PARI(93).
VINT(138) : size of the threshold factor (enhancement or suppression)
    in the latest event with heavy flavour production; see MSTP(35).
VINT(141), VINT(142) : x values for the parton shower initiators of the
    hardest interaction; used to find what is left for multiple
    interactions.
VINT(143), VINT(144) : 1 - (sum of x-values) for all scatterings; used
    for rescaling each new x-value in the multiple interaction structure
    function evaluation.
VINT(145) : estimate of total parton-parton cross-section for multiple
    interactions; used for MSTP(82) >= 2.
VINT(146) : common correction factor f_c in the multiple interaction
    probability; used for MSTP(82) >= 2.
VINT(147) : average hadronic matter overlap; used for MSTP(82) >= 2.
VINT(148) : enhancement factor for current event in multiple interaction
    probability, defined as the actual overlap divided by the average
    one; used for MSTP(82) >= 2.
VINT(149) : x_T^2 cutoff or turnoff for multiple interactions.
    For MSTP(82) <= 1 it is 4*p_Tmin^2/s, for MSTP(82) >= 2 it is
    4*p_T0^2/s.
VINT(150) : probability to keep given event in multiple interaction
    scenario, as given by the "Sudakov" form factor.
VINT(151), VINT(152) : sum of x values for all multiple interaction
    partons.
VINT(153) : current differential cross-section value obtained from
    PYSIGH; used in multiple interactions only.
VINT(155), VINT(156) : the x value of a photon that branches into
    quarks or gluons, i.e. x at interface between initial state QED
    and QCD cascades.
VINT(157), VINT(158) : the primordial k_T values selected in the two
    beam remnants.
VINT(159), VINT(160) : the chi values selected for beam remnants that
    are split into two objects, describing how the energy is shared
    (see MSTP(92) etc.); is 0. if no splitting is needed.
VINT(161) - VINT(200) : sum of Cabibbo-Kobayashi-Maskawa matrix
    element-squared that a given flavour is allowed to couple to.
    Results are stored in format VINT(180+KF) for quark and lepton
    flavours and antiflavours (which need not be the same; see
    MDME(IDC,2). For leptons, these factors are normally unity.
VINT(201) - VINT(220) : additional variables needed in phase space
    selection for 2 -> 3 processes with ISET(ISUB) = 5. Below indices
    1, 2 and 3 refer to scattered partons 1, 2 and 3, except that
    q variables are q_1 + q_2 -> q_1' q_2' q_3' for four-momenta.
    All kinematical variables refer to the internal kinematics of
    the 3-body final state - the kinematics of the system as a whole
    is described by tau' and y*, and the mass distribution of
    particle 3 (a resonance) by tau.
VINT(201) : m_1.
VINT(202) : p_T1^2.
VINT(203) : phi_1.
VINT(204) : M_1 (mass of propagator).
VINT(205) : weight for the p_T1^2 choice.
VINT(206) : m_2.
VINT(207) : p_T2^2.
VINT(208) : phi_2.
VINT(209) : M_2 (mass of propagator).
VINT(210) : weight for the p_T2^2 choice.
VINT(211) : y_3.
VINT(212) : y_3max.
VINT(213) : epsilon = +-1; choice between two mirror solutions 1 <-> 2.
VINT(214) : weight associated to epsilon-choice.
VINT(215) : t_1 = (q_1 - q_1')^2.
VINT(216) : t_2 = (q_2 - q_2')^2.
VINT(217) : q_1*q_2'.
VINT(218) : q_2*q_1'.
VINT(219) : q_1'*q_2'.
VINT(220) : sqrt((m_T12^2-m_T1^2-m_T2^2)^2-4.*m_T1^2*m_T2^2).

COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
Purpose: to store information necessary for efficient generation of the
    different subprocesses, specifically type of generation scheme and
    coefficients of the Jacobian. Also to store cumulative
    cross-sections needed for multiple interaction generation for
    MSTP(82) >= 2. These variables must not be changed by the user.

ISET(ISUB) : gives the type of kinematical variable selection scheme
    used for subprocess ISUB.
    = 0 : elastic, diffractive and low-p_T processes.
    = 1 : 2 -> 1 processes (irrespective of subsequent decays).
    = 2 : 2 -> 2 processes (i.e. the bulk of processes).
    = 3 : 2 -> 3 processes (like q q' -> q" q'" H0).
    = 4 : 2 -> 4 processes (like q q' -> q" q'" W+ W-).
    = 5 : 'true' 2 -> 3 processes, one method.
    = 6 : 'true' 2 -> 3 processes, another method;
        currently only g + g -> Z + Q + Q~.
    = 9 : 2 -> 2 in multiple interactions (pT as kinematics variable).
    = -1 : legitimate process which has not yet been implemented.
    = -2 : ISUB is an undefined process code.

KFPR(ISUB,J) : give the KF flavour codes for the products produced in
    subprocess ISUB. If there is only one product, the J = 2 position
    is left blank. Also, quarks and leptons assumed massless in the
    matrix elements are denoted by 0. The main application is thus to
    identify resonances produced (Z0, W+, H0, etc.).

COEF(ISUB,J) : factors used in the Jacobians in order to speed up the
    selection of kinematical variables. More precisely, the shape of
    the cross-section is given as the sum of terms with different
    behaviour, where the integral over the allowed phase space is
    unity for each term. COEF gives the relative strength of these
    terms, normalized so that the sum of coefficients for each variable
    used is unity. Note that which coefficients are indeed used is
    process-dependent.
    ISUB : standard subprocess code.
    J = 1 : tau selected according 1/tau.
    J = 2 : tau selected according to 1/tau^2.
    J = 3 : tau selected according to 1/(tau*(tau+tau_r)), where
        tau_r = m^2/s is tau value of resonance; only used for
        resonance production.
    J = 4 : tau selected according to Breit-Wigner of form
        1/((tau-tau_r)^2+gam_r^2), where tau_r = m^2/s is tau value
        of resonance and gam_r = m*Gamma/s is its scaled mass times
        width; only used for resonance production.
    J = 5 : tau selected according to 1/(tau*(tau+tau'_r)), where
        tau'_r = m'^2/s is tau value of second resonance; only
        used for simultaneous production of two resonances.
    J = 6 : tau selected according to second Breit-Wigner of form
        1/((tau-tau'_r)^2+gam'_r^2), where tau'_r = m'^2/s is tau value
        of second resonance and gam'_r = m'*Gamma'/s is its scaled
        mass times width; is used only for simultaneous production
        of two resonances, like gamma*/Z0/Z'0.
    J = 7 : tau selected according to 1/(1-tau); only used when both
        structure functions are peaked at x = 1.
    J = 8 : y* selected according to y* - y*min;
    J = 9 : y* selected according to y*max - y*;
    J = 10 : y* selected according to 1/cosh(y*).
    J = 11 : y* selected according to 1/(1-exp(y*-y*max)); only used
        when beam structure function is peaked close to x = 1.
    J = 12 : y* selected according to 1/(1-exp(y*min-y*)); only used
        when target structure function is peaked close to x = 1.
    J = 13 : z = cos(theta-hat) selected evenly between limits.
    J = 14 : z = cos(theta-hat) selected according to 1/(a-z),
        where a = 1 + 2*m_3^2*m_4^2/shat^2, m_3 and m_4 being the
        masses of the two final state particles.
    J = 15 : z = cos(theta-hat) selected according to 1/(a+z),
        with a as above.
    J = 16 : z = cos(theta-hat) selected according to 1/(a-z)^2,
        with a as above.
    J = 17 : z = cos(theta-hat) selected according to 1/(a+z)^2,
        with a as above.
    J = 18 : tau' selected according to 1/tau'.
    J = 19 : tau' selected according to (1 - tau/tau')^3/tau'^2.
    J = 20 : tau' selected according to 1/(1-tau'); only used when
        both structure functions are peaked close to x = 1.

ICOL : contains information on different colour flow topologies in
    hard 2 -> 2 processes.

COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
Purpose: to store information on structure functions, subprocess cross-
    sections and different final state relative weights. These variables
    must not be changed by the user.

XSFX : current values of structure functions (multiplied by x) on
    beam and target side.

ISIG(ICHN,1) : incoming parton/particle on the beam side to the hard
    interaction for allowed channel no. ICHN. The number of channels
    filled with relevant information is given by NCHN, one of the
    arguments returned in a PYSIGH call. Thus only 1 <= ICHN <= NCHN
    is filled with relevant information.
ISIG(ICHN,2) : incoming parton/particle on the target side to the hard
    interaction for allowed channel no. ICHN. See also comment above.
ISIG(ICHN,3) : colour flow type for allowed channel no. ICHN; see
    MSTI(2) list. See also comment above. For 'subprocess' 96 uniquely,
    ISIG(ICHN,3) is also used to translate information on what is the
    correct subprocess number (11, 12, 13, 28, 53 or 68); this is
    used for reassigning subprocess 96 to either of these.

SIGH(ICHN) : evaluated differential cross-section for allowed channel
    no. ICHN, i.e. hard matrix element value times structure
    functions, for current kinematical setup (in addition, Jacobian
    factors are included in the figures, as used to speed up
    generation). See also comment for ISIG(ICHN,1).

COMMON/PYINT4/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3)
Purpose: to store partial and effective decay widths for the different
    resonances. These variables must not be changed by the user.

WIDP(KF,J) : gives partial decay widths of resonances into different
    channels (in GeV), given that all physically allowed final
    states are included.
    KF : standard KF code for resonance considered.
    J : enumerates the different decay channels possible for
        resonance KF, as stored in the JETSET LUDAT3 commonblock,
        with the first channel in J = 1, etc.

WIDE(KF,J) : gives effective decay widths of resonances inte different
    channels (in GeV), given the decay modes actually left open in
    the current run. The on/off status of decay modes is set by the
    MDME switches in JETSET; see section 2.9.
    KF : standard KF code for resonance considered.
    J : enumerates the different decay channels possible for
        resonance KF, as stored in the JETSET LUDAT3 commonblock,
        with the first channel in J = 1, etc.

WIDS(KF,J) : gives a dimensionless suppression factor, which is defined
    as the ratio of the total width of channels switched on to the total
    width of all possible channels (replace width by width-squared for
    a pair of resonances). The on/off status of channels is set by the
    MDME switches in JETSET; see section 2.9. The information in WIDS
    is used e.g. in cross-section calculations.
    KF : standard KF code for resonance considered.
    J = 1 : suppression when a pair of resonances of type KF are
        produced together. When an antiparticle exists, the
        particle-antiparticle pair (like W+ W-) is the relevant
        combination, else the particle-particle one (like Z0 Z0).
    J = 2 : suppression for a particle of type KF when produced
        on its own, or together with a particle of another type.
    J = 3 : suppression for an antiparticle of type KF when produced
        on its own, or together with a particle of another type.

COMMON/PYINT5/NGEN(0:200,3),XSEC(0:200,3)
Purpose: to store information necessary for cross-section calculation
    and differential cross-section maximum violation. These variables
    must not be changed by the user.

NGEN(ISUB,1) : gives the number of times that the differential
    cross-section (times Jacobian factors) has been evaluates for
    subprocess ISUB, with NGEN(0,1) the sum of these.
NGEN(ISUB,2) : gives the number of times that a kinematical setup
    for subproces ISUB is accepted in the generation procedure,
    with NGEN(0,2) the sum of these.
NGEN(ISUB,3) : gives the number of times an event of subprocess type
    ISUB is generated, with NGEN(0,3) the sum of these. Usually
    NGEN(ISUB,3) = NGEN(ISUB,2), i.e. an accepted kinematical
    configuration normally can be used to produce an event.

XSEC(ISUB,1) : estimated maximum differential cross-section (times the
    Jacobian factors used to speed up the generation process) for
    the different subprocesses in use, with XSEC(0,1) the sum of these
    (except low-p_T, i.e. ISUB = 95).
XSEC(ISUB,2) : gives the sum of differential cross-sections (times
    Jacobian factors) for the NGEN(ISUB,1) phase space points
    evaluated so far.
XSEC(ISUB,3) : gives the estimated integrated cross-section for
    subprocess ISUB, based on the statistics accumulated so far,
    with XSEC(0,3) the estimated total cross-section for all
    subprocesses included (all in mb). This is exactly the
    information obtainable by a PYSTAT(1) call.

COMMON/PYINT6/PROC(0:200)
Purpose: to store character strings for the different possible
    subprocesses; used when printing tables.

PROC(ISUB) : name for the different subprocesses, according to ISUB
    code. PROC(0) denotes all processes.
_________________

Finally, in addition a number of routines and commonblocks with names
beginning with RK come with the program. These contain the matrix
element evaluation for the process g + g -> Z + q + q~, based on a
program of Ronald Kleiss, with only minor modifications.
______________________________________________________________________

2.9.  The JETSET routines

Since the fragmentation and decay of an original parton configuration is
handled by the JETSET routines, using the LUJETS common block described
above, all of the non-e+e- routines and parameters are at the disposal
of the PYTHIA user. This section is intended as a brief reminder of some
of the most useful ones, in particular for hadron physics applications.
Further details may be obtained by consulting the JETSET manual.

The subroutine LULIST can be used to list an event; for output that fits
onto a normal terminal screen CALL LULIST(1) is recommended, whereas
CALL LULIST(2) gives a slightly more comprehensive listing.  The event
record nominally contains not only the stable final particles, but also
the original partons and intermediate resonances. With CALL LUEDIT(1)
all partons/particles that have fragmented/decayed are removed from the
event record, and N is updated accordingly. A CALL LUEDIT(2) will remove
neutrinos as well, and CALL LUEDIT(3) will only leave charged, stable
particles. Some of the information indirectly stored in LUJETS can be
accessed more easily by using the PLU and KLU functions, thus e.g.
y = PLU(I,17) gives the true rapidity and eta = PLU(I,19) the
pseudorapidity of particle I.

With MSTP(111) = 0, fragmentation and decays are switched off in
the PYEVNT call. If one desires to fragment the event at some later
opportunity, this can be achieved by a CALL LUEXEC. The action of
the program during this call, which is done automatically from
PYEVNT for MSTP(111) = 1, can be controlled by the switches and
parameters in the LUDAT1 common block:
      COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
Thus, MSTJ(21) = 0 will switch off particle decays but keep jet
fragmentation on, MSTJ(1) can be set to give independent
fragmentation rather than the standard string fragmentation, etc.
Contrary to most PYTHIA variables, the user can freely change the
values during the course of the run, since no initialization is
involved.

The same commonblock also contains a number of conversion factors and
coupling constants. Thus MSTU(101) can be set to used fix or running
alpha_em, PARU(101) to set the alpha_em value, PARU(102) to set
sin^2(theta_W). A host of couplings for particles beyond the
standard model may be found in PARU(121) onwards. These couplings can
not be changed in midrun, since they are used in the initialization to
estimate cross-section maxima.

The common block LUDAT2 contains several pieces of information on
particle and flavour properties:
      COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
The particle masses may be the most interesting variables here.
Since these masses are directly taken over by the PYTHIA routines,
any mass changes should be made before the PYINIT call. Of
particular interest are the following:
PMAS(6,1) : (D=120. GeV/c^2) mass of top quark.
PMAS(7,1) : (D=120. GeV/c^2) mass of l quark.
PMAS(8,1) : (D=200. GeV/c^2) mass of h quark.
PMAS(17,1) : (D=100. GeV/c^2) mass of chi lepton.
PMAS(23,1) : (D=91.2 GeV/c^2) Z0 nominal mass.
PMAS(24,1) : (D=80.0 GeV/c^2) W+/- nominal mass.
PMAS(25,1) : (D=50. GeV/c^2) H0 nominal mass.
PMAS(32,1) : (D=500. GeV/c^2) Z'0 nominal mass.
PMAS(34,1) : (D=500. GeV/c^2) W'+/- nominal mass.
PMAS(35,1) : (D=300. GeV/c^2) H'0 nominal mass.
PMAS(36,1) : (D=300. GeV/c^2) A0 nominal mass.
PMAS(37,1) : (D=300. GeV/c^2) H+/- nominal mass.
PMAS(39,1) : (D=200. GeV/c^2) LQ nominal mass.
PMAS(40,1) : (D=5000. GeV/c^2) R nominal mass.
For the resonances above (particles 23 - 40), the resonance width is
calculated at PYINIT initialization, taking into account the channels
assumed possible for resonance decays (for H0 this will depend
very much on the mass assumed), so the nominal widths stored in
PMAS(KF,2) are actually overwritten at this point.

Note that the first index in the KCHG and PMAS arrays is not the
standard KF code described in section 2.7, which may be both positive
and negative, and extends to very large values. Instead a compressed
code KC is used, which only runs between 1 and 500, i.e. which contains
no distinction between particles and antiparticles. The KF and KC
codes agree between 1 and 40, i.e. for quarks, leptons and gauge bosons,
while meson and baryon codes have been compressed. Thus, several
hadrons containing heavy quarks (some c and b baryons, most t, l
and h hadrons) have been lumped together, with masses and charges
instead derived from the quark content (and, for the mass, the spin).
The user is not supposed to learn the KC code, but instead use the
LUCOMP function to translate from KF to KC code, i.e. KC = LUCOMP(KF).
PMAS(LUCOMP(663),1) = 110. changes the toponium mass to 110 GeV, to
give one example. LUCOMP will return KC = 0 for any particle code
unknown in the program.

The common block LUDAT3 contains all particle decay data:
      COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
It can be used to switch off the decay of some particles selectively,
and even to switch on and off separate decay modes for a given
particle. These switches are used extensively within PYTHIA, both
to select allowed final state flavours and to calculate the reduction
in cross-section due to these restrictions. MDCY(KC,1) is an on/off
switch (= 0 off, = 1 on) for the decay of compressed particle code KC.
The decay of Z0, both in PYTHIA and in JETSET, is thus switched off
by putting MDCY(23,1) = 0. In order to keep the event record from
becoming uncomfortably large, switching off pi0 decays by putting
MDCY(LUCOMP(111),1) = 0 is probably the single most effective action.

All the decay channels defined in the program are stored sequentially
in MDME, BRAT and KFDP. The entry point into this table and the number
of channels stored for a given particle is given by the last two
entries of the MDCY array. In order to obtain a list of particle data,
including all decay channels, do CALL LULIST(12) (and take out a
hardcopy for future reference). For resonances only, the same
information is available with CALL PYSTAT(2). For a given decay channel
IDC, the switch in MDME(IDC,1) is of special interest. It can be used to
switch on or off that particular channel, or leave it selectively open.
Effective branching ratios are automatically recalculated among the
decay channels left open. If a particle is allowed to decay by the
MDCY(KC,1) value, at least one channel must be left open by the user.
Since the MDME(IDC,1) code is of special interest in PYTHIA for
selecting flavours in single gauge boson or gauge boson pair decays,
the JETSET manual description on this point is reproduced in extenso:
= -1 : this is a non-standard model decay mode, which by default is
    assumed not to exist. It is therefore not simulated; neither does
    it contribute to the total width of the particle. Normally, this
    option is used for decays involving fourth generation or H+-
    particles.
= 0 : channel is switched off. It is not simulated, but still
    contributes to the total width of the particle.
= 1 : channel is switched on.
= 2 : channel is switched on for a particle but off for an
    antiparticle. It is also on for a particle its own antiparticle,
    i.e. here it means the same as =1.
= 3 : channel is switched on for an antiparticle but off for a
    particle. It is off for a particle its own antiparticle.
= 4 : in the production of a pair of equal or charge conjugate
    resonances in PYTHIA, say H0 -> W+ W-, either one of the
    resonances is allowed to decay according to this group of
    channels, but not both. If the two particles of the pair
    are different, the channel is on.
    Within JETSET, this option only means that the channel is
    switched off.
= 5 : as =4, but an independent group of channels, such that in
    a pair of equal or charge conjugate resonances the decay of
    either resonance may be specified independently. If the two
    particles in the pair are different, the channel is off.
    Within JETSET, this option only means that the channel is
    switched off.
WARNING: the two values -1 and 0 may look similar, but in fact
    are quite different. In neither case the channel so set is
    generated, but in the latter case the channel still contributes
    to the total width of a resonance, and thus affects both
    simulated line shape and the generated cross-section. Thus the
    value 0 is appropriate to a channel we assume exists, even if
    we are not currently simulating it, while -1 should be used for
    channels we believe do not exist. In particular, users are warned
    unwittingly to set fourth generation channels 0 (rather than -1),
    since by now the support for a fourth generation is small.
All the options above may be freely mixed. The difference,
for those cases where both make sense, between using values
2 and 3 and using 4 and 5 is that the latter automatically
include charge conjugate states, e.g. H0 -> W+ W- ->
e+ nue d ubar or dbar u e- nuebar, but the former only one
of them. In calculations of the joint branching ratio, this
makes a factor 2 difference.

One should note that, whereas the MDME code can be used to set
allowed decay channels both for resonances (Z0, W+/-, H0, Z'0,
W'+/-, H+/-, and R) and for ordinary particles (pi0, K, D, B,
top mesons, etc.), it is only for the resonances that cross-sections
are modified to take into account the allowed channels. The reason
is obvious: resonances are only produced as part of the hard
interaction, or as decay products of other resonances, according
to fairly strict rules, whereas ordinary particles can be produced
anywhere in the event (charm and bottom quark pairs may appear
inside initial and final state showers, e.g.) in a rather haphazard
fashion. It is therefore not possible to generate unbiased events
in a consistent manner for ordinary particles. (One might imagine
doing things properly for the t quark, since it is not likely to be
produced in showers, but this has not been done so far.)

For QCD processes, like the pair production of a new heavy quark
(ISUB = 81 and 82), imagined 'decay modes' are stored for the
gluon among the ordinary ones. Branching ratios here have no meaning,
but the MDME(IDC,1) switch can be used to enable or disable the
production of different new quarks. In the current version of the
massive matrix element treatment, only one heavy flavour can be
generated in a given one; if the user has left several flavours
switched on, only the heaviest is actually used. The same goes for
ISUB = 84, where the photon 'decay modes' can be selected. Similarly,
the 'decay modes' stored for quarks indicate which branchings of the
type q -> q' W are allowed, as used e.g. in W+ W- -> H0 production.

It is also possible to include completely new decay modes, or even a
completely new particle. This is described in the JETSET manual.

The event analysis routines can be called on to evaluate the event
stored in the LUJETS common block. Particularly useful is the cluster
finding routine LUCELL, which assumes a (modifiable) grid of cells in
pseudorapidity and azimuthal angle and carries out a cluster finding
algorithm a la the UA1 one. With a CALL LUCELL(NJET) the number of jets
found is given by NJET and the position and ET of those jets is stored
after the event proper, in lines N+1 through N+NJET.

The routine LUGIVE can be used to feed in commonblock values both to
JETSET and to PYTHIA. Since LUGIVE checks on array boundaries and
also writes on output the variables which have been changes, this may
at times be a better way to set up the values to be used by the program.
A call might look like
      CALL LUGIVE('MSEL=0;MSUB(142)=1;PMAS(34,1)=300.;CKIN(1)=200.')
i.e. several commands can be put in one if separated by semicolon
(but it is of course also possible to do one command at a time and
call LUGIVE as often as desired). Note that for people who want data
card driven programs, this facilitates the construction of a driver
routine for PYTHIA.
______________________________________________________________________

2.10. On cross-sections

This section includes two related topics. Thew first is how events are
generated according to the proper matrix elements, the second how
the final cross-sections (available e.g. with PYSTAT) are calculated.

In PYRAND, the variables used in the generation of hard scattering
2 -> n processes are tau (= s-hat/s, i.e. mass-squared of
scattering subsystem scaled by total CM energy-squared), y*
(rapidity of scattering subsystem in the CM frame of the event),
z (= cos(theta-hat), i.e. cosine of scattering angle in CM frame
of the hard interaction; this applies to 2 -> 2 and 2 -> 4
processes only), and tau' (= s'/s, i.e. mass-squared of complete
three- or four-body final state scaled by total CM energy-squared;
this applies to 2 -> 3 and 2 -> 4 processes only). These variables
are generated according to the distributions h1(tau)/tau, h2(y*),
h3(z), and h4(tau')/tau', where

h1(tau)  = c0 + I0/I1*c1*1/tau + I0/I2*c2*1/(tau + tau_R) +
         + I0/I3*c3*tau/((s*tau - m_R^2)^2 + (m_R*Gamma_R)^2) +
         + I0/I4*c4*1/(tau + tau_R') + I0/I5*c5 *
         * tau/((s*tau - m_R'^2)^2 + (m_R'*Gamma_R')^2)

h2(y*)   = I0/I1*c1*(y* - y*_min) + I0/I1*c2*(y*_max - y*) +
         + I0/I3*c3*1/cosh(y*)

h3(z)    = c0 + I0/I1*c1*1/(a - z) + I0/I2*c1*1/(a + z) +
         + I0/I3*c3*1/(a - z)^2 + I0/I4*c4*1/(a + z)^2

h4(tau') = c0 + I0/I1*c1*(1 - tau/tau')^3/tau',

where subscripts R and R' in h1 denote values (of tau, mass and
width) pertaining to a resonance R or R', and a in h3 denotes
1 + 2*(m_3*m_4/s-hat)^2 (= 1 for massless products).

In each of the above cases, In denotes the integral over the
quantity multiplying coefficient cn; thus, e.g.,

h1: I0 = Int dtau/tau = ln(tau_max/tau_min)

h3: I0 = Int dz = z-_max - z-_min + z+_max - z+_min

h4: I0 = Int dtau'/tau' = ln(tau'_max/tau'_min)

(for h2, where no c0 occurs, I0 is defined to be Int dy* =
= y*_max - y*_min), and the coefficients cn are normalized, such
that Sum(n) cn = 1. (The symbols Int and Sum are used to denote
integral and sum over, respectively.)

For other types of processes, e.g. elastic and diffractive
scattering, other variables are used in the generation, but the
structure of the generation procedure is similar.

After this primary generation, the variables are then weighted
with the relevant matrix elements, etc., according to standard
Monte Carlo techniques.

The process-dependent weights are coded in the routine PYSIGH, as
follows. The cross-sections for 2 -> 1 and 2 -> 2 processes can be
written, respectively,

sigma(2 -> 1) = Int dx_1 Int dx_2 Sum(ijk) f_i f_j *
              * dsigma-hat(ijk)

and

sigma(2 -> 2) = Int dx_1 Int dx_2 Int dt-hat Sum(ijk) f_i f_j *
              * dsigma-hat(ijk)/dt-hat,

where Int and Sum denote respectively integral and sum over, f_i
the structure function for parton i, and dsigma-hat the
differential cross-section for the hard scattering. Converting to
the variables tau, y*, and z= cos(theta-hat), the cross-sections
can then be rewritten as

sigma(2 -> 1) = pi/s * Int dtau h1(tau)/tau Int dy* h2(y*) *
              * 1/(tau*h1(tau)) * 1/h2(y*) * Sum(ijk) F_i F_j *
              * s-hat/pi*dsigma-hat(ijk) =

              = pi/s * I(tau) * I(y*) * 1/(tau*h1) * 1/h2 *
              * Sum(ijk) F_i F_j * s-hat/pi*dsigma-hat(ijk)

and

sigma(2 -> 2) = pi/s * Int dtau h1(tau)/tau Int dy* h2(y*) *
              * Int dz h3(z) * 1/(tau*h1(tau)) * 1/h2(y*) *
              * 1/2*beta34/h3(z) * Sum(ijk) F_i F_j *
              * s-hat^2/pi*dsigma-hat(ijk)/dt-hat =

              = pi/s * I(tau) * I(y*) * I(z) * 1/(tau*h1) *
              * 1/h2 * 1/2*beta34/h3 * Sum(ijk) F_1 F_j *
              * s-hat^2/pi*dsigma-hat(ijk)/dt-hat,

where F_i denotes x*f_i, beta34 = lambda(1,r_3,r_4) with lambda
the usual lambda-function and r_i = m_i^2/s-hat, and

    I(tau) = Int dtau h1(tau)/tau

    I(y*) = Int dy* h2(y*)

    I(z) = Int dz h3(z),

and the h's are distributions used in the actual generation of
the kinematical variables tau, y* and z (see above), such that

    Int dtau h1(tau)/tau = ln(tau_max/tau_min)

    Int dy* h2(y*) = y*_max - y*_min

    Int dz h3(z) = z-_max - z-_min + z+_max - z+_min.

For each subprocess in PYSIGH, what is coded under the relevant
ISUB is thus the dimensionless quantity

    s-hat/pi*dsigma-hat(ijk)                (2 -> 1 processes)

    s-hat^2/pi*dsigma-hat(ijk)/dt-hat     (2 -> 2 processes)

which is then multiplied by the structure functions (or, rather,
x*f) and a common factor, containing pi/s, the conversion factor
from GeV^-2 to mb, I(tau), I(y*) (and, for 2 -> 2 processes,
I(z)) and 1/(tau*h1)), and 1/h2 (and, for 2 -> 2 processes, also
1/2*beta34/h3). It is the latter product which is returned in a
PYSIGH call; we will denote this SIGH(tau,y*,z) in the following.

For 2 -> 2 processes with identical final state particles, the
symmetrization factor of 1/2 is explicitly included at the very
end, when the result is stored in the SIGH array. In the final
cross-section, a factor of 2 is retrieved because of integration
over the full phase space (rather than only half of it).

For 2 -> 1 processes with dsigma-hat given in the zero-width
approximation, the delta function in tau has been replaced by
the (modified) Breit-Wigner

    1/pi*s*H_R/((s*tau - m_R^2)^2 - H_R^2),

where

    H_R = s*tau/m_R*Gamma_R,

with m_R and Gamma_R the mass and the width of the resonance R,
respectively. The Breit-Wigner thus contains an shat-dependent
width which improves the description of the line shape as compared
to the case with widths calculated at nominal resonance mass. In
the 2 -> 2 processes 71 - 77, where the Higgs occurs in s-, t- and
u-channel exchanges, the width is calculated at nominal Higgs mass,
however.

For effective 2 -> 3, 4 processes (e.g. Z0 Z0 -> H0), the common
factor further includes the integration over the extra variable
tau', i.e. a factor I(tau') * f(tau')/h4(tau'), where

    f(tau') = (1 + tau/tau') ln(tau'/tau) - 2 (1 - tau/tau')

and

    I(tau') = Int dtau' h4(tau')/tau' = ln(tau'_max/tau'_min).

A user wanting to include new processes (experts only!), will
have to insert then the function s-hat/pi*dsigma-hat(ijk) or
s-hat^2/pi*dsigma-hat(ijk)/dt-hat under the appropriate ISUB
heading (in addition to changes in e.g. SUBROUTINE PYSCAT).

Let us now describe how these differential cross-sections are used
to derive the total cross-sections (within the user-defined cuts)
for the allowed processes. We will first consider high-p_T event
generation only, and later cover low-p_T physics.

At initialization of a 2 -> 2 process, the coefficients cn appearing
in the expressions h1(tau), h2(y*) and h3(z) are first optimized
so as to fit the actual variation of the cross-section as well as can
be achieved, i.e. to see to it that SIGH(tau,y*,z) is a slowly varying
function of its arguments. Thereafter the maximum of SIGH(tau,y*,z)
is searched for. This gives a preliminary cross-section estimate,
which is printed among the initialization information, and which
represents an overestimate of the actual cross-section. It is this
maximum which is used below to find the relative weight for selecting
a processes. For 2 -> 1, 2-> 3 and 2 -> 4 processes, the same
philosophy is used, with trivial changes in the number of variables
that have to be considered.

In particular, note that SIGH is selected so as to allow a direct
comparison between the different classes of events. More precisely,
the SIGH value returned by PYSIGH is the total integrated cross-
section a process would have inside the simulated phase space volume,
if the actual variation of the cross-section were perfectly described
by h1(tau)/tau * h2(y*) * h3(z) (or equivalent), with overall
normalization given by the differential cross-section value in the
point selected.

In order to generate an event, first a process is chosen, thereafter
a set of kinematical variables for this process is selected, as
described above. The probability for accepting this set is given by
the actual quantity SIGH(tau,y*,z), divided by the maximum of this
quantity, as found earlier. In case of rejection, both process and
kinematical variables are to be chosen anew.

For each individual subprocess, there are counters giving the number
of points tried and the sum of the SIGH values in these points. At
the end of each PYEVNT call, the ratio of these two numbers is
calculated, thus providing a Monte Carlo estimate of the actual
cross-section for each process. If only very few events have been
generated, the statistical fluctuations are large. Therefore it is
recommended to call PYSTAT only after all events have been generated,
although a call in principle could be made at any time. In PYSTAT
the cross-sections will be listed together with the number of phase
space points tried and the number of events actually generated.

The cross-sections for each subprocess are calculated taking into
account which initial and final states are allowed, by the KFIN array
in the common block PYSUBS, and the MDME(IDC,1) values in the common
block LUDAT3, respectively. By default, everything up to and including
top is allowed, while it is assumed that no fourth generation exists.
If some initial states are switched off, the cross-sections will be
reduced accordingly. For several processes, like q + g -> q + g,
the flavours in the final state are the ones of the initial state,
and for these MDME has no effect. In processes like
g + g -> q + qb or q + qb -> Z0/gamma* -> mu+ + mu-, on the other hand,
the quark and lepton flavours allowed in the final state can be selected
using MDME, and cross-sections are only based on those channels that
are actually allowed. We emphasize that, for resonances like Z0 or W+/-,
the total widths are assumed given as the sum over three (or,
optionally, four) full generations of quarks and leptons. By only
allowing a few channels, we do thus in no way alter the basic properties
of the resonances, but only restrict the generation of events to those
particular channels we happen to be interested in, with the cross-
sections reduced accordingly. Particles like H0, Z' and W' can either
decay directly into a q-q~ or l-l~ pair, or into a (gauge) boson pair,
which then can decay further. Therefore the production cross-section
given by the program is reduced, not only by switching off some decay
channels of the primary resonance, but also by switching off subsequent
Z, W or H+ decay modes (and the branching ratios of the primary
resonance into the boson pair is reduced accordingly). Even if the
decay of a resonance is switched off completely, by using the MDCY
array in the LUDAT3 commonblock, cross-sections are still based on
the decay channels allowed by the MDME values in the same commonblock.
Note that what is said here about resonances does not hold for
ordinary particles or even the top quark, where a suppression of some
decay channels are never taken into account in PYTHIA cross-sections,
but have to be compensated for by the user himself.

The cross-section values printed by PYSTAT are stored in XSEC(ISUB,3)
in common block PYINT5. In particular, XSEC(0,3) = PARI(1) gives the
total cross-section for all subprocesses studied in current run. Events
are generated with unit weight. At the end of a run, histogram contents
etc. can be converted from number of events to differential cross-
sections by multiplying by XSEC(0,3) and dividing by the total number of
events, the latter stored in MSTI(5). This ratio is given in PARI(2).
For histograms one should also divide by the bin width.

What has been said so far applies when only hard interactions are
studied. The total hadronic cross-section, as well as the elastic,
single diffractive and double diffractive cross-sections are taken from
parametrizations available in the literature. For the total and elastic
cross-sections this means the fits of Block-Cahn [Blo85], for
diffractive events the ansatz of Goulianos [Gou83] is used.
The cross-section for non-diffractive low-pT events is then
obtained by subtraction. When the multiple interaction formalism is used
to generate these latter events, some of them will actually contain a
semihard interaction, and be listed in PYSTAT as such, with only the
events without any semihard interactions being classified as low-p_T
(ISUB = 95). This matter of classification does not remove the
constraint that the total cross-section is fixed, however. Neither can
the KFIN switches be set to suppress some of the flavours of the
initial hadrons (the inclusion of heavy flavours or not can still be
set using MSTP(54), however).
______________________________________________________________________

2.11. Examples

The program is built as a slave system, i.e. the user supplies the main
program, which calls on the PYTHIA and JETSET routines to perform
specific tasks and then resumes control.

A typical program for analysis of collider events at 540 GeV CM energy
with a minimum p_T of 10 GeV/c at the hard scattering (because of
initial state radiation, fragmentation effects, etc., the actual
p_T-cutoff will be smeared around this value) might look like
      COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
      COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
      COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
      ...                   ! set all common block variables that
      ...                   ! did not have desired default values
      CKIN(3)=10.           ! lower p_T cutoff
      CALL PYINIT('CMS','p','pbar',540.)         ! initialize
      ...                   ! initialize analysis statistics
      DO 100 IEVENT=1,1000                ! loop over events
      CALL PYEVNT                         ! generate event
      IF(IEVENT.EQ.1) CALL LULIST(1)      ! list first event
      ...                   ! insert desired analysis chain for
      ...                   ! each event
  100 CONTINUE
      CALL PYSTAT(1)        ! print cross-sections
      ...                   ! user output
      END
______________________________________________________________________

2.12. How to include external processes in PYTHIA

Despite a large repertory of processes in PYTHIA, clearly the number
of missing ones is even larger, and with time this discrepancy is
likely to increase. There are several reasons why it is not
practicable to imagine a PYTHIA which has 'everything'. One is
the amount of time it takes to implement a process for the single
PYTHIA author, compared to the rate of new cross-section results
produced by the rather larger matrix element calculations community.
Another is the length of currently produced matrix element expressions,
which would make the program very bulky. A third argument is that,
whereas the phase space of 2 -> 1 and 2 -> 2 processes can be set up
once and for all according to a reasonably flexible machinery,
processes with more final state particles are less easy to generate.
To achieve a reasonable efficiency, it is necessary to tailor the
phase space selection procedure to the dynamics of the given process,
and to the desired experimental cuts.

We have therefore added the possibility to include 'external' processes
into PYTHIA. In this section we will describe how it is possible to
specify the partonic state of some hard scattering process in an
interface commonblock. PYTHIA will read this commonblock, and add
initial and final state showers, beam remnants and underlying events,
fragmentation and decays, to build up an event in as much detail as an
ordinary PYTHIA one. You may also use PYTHIA to mix events of different
kinds, and to keep track of cross-section statistics. You have to
provide the matrix elements, the phase space generator, and the storage
of event information in the commonblock.

First a minor comment, however. Some processes may be seen as just
trivial modifications of already existing ones. For instance, you might
want to add some additional term, corresponding to contact interactions,
to the matrix elements of a PYTHIA 2 -> 2 process. In that case it is
not necessary to go through the machinery below, but instead you can use
the PYEVWT routine to introduce an additional weight for the event,
defined as the ratio of the modified to the unmodified differential
cross-section. If you use the option MSTP(142)=2, this weight is
considered as part of the 'true' cross-section of the process, and
the generation is changed accordingly.

The more generic facility for including an external process is a bit
more complicated, and involves two routines and one commonblock. All
names contain 'UP', which is short for User Process.

If you want to include a new process, first you have to pick an unused
subprocess number ISUB, see table in section 2.2. For instance,
currently the numbers 191 - 200 are unused, so this might be a logical
place to put a new process. This number and the 'title' of the process
(plus SIGMAX, to be described below) have to be given in to PYTHIA in
a subroutine call
      CALL PYUPIN(ISUB,TITLE,SIGMAX)
before the call to PYINIT. The TITLE can be any character string up to
28 characters long, e.g.
      CALL PYUPIN(191,'g + g -> t + tbar + gamma',SIGMAX)
The call to PYUPIN tells the program that a process ISUB exists, but not
that you want to generate it. This is done, as with normal processes,
by setting MSUB(ISUB)=1 before the PYINIT call.

Once the event generation chain has been started and PYEVNT is called to
generate an event, this routine may in its turn call the routine
PYUPEV, which is the routine you must supply, in which the next
event is selected. (A dummy copy of PYUPEV has been included at the end
of PYTHIA; depending on machine you may have to comment out this copy
when you link your own.) The call arguments are
      CALL PYUPEV(ISUB,SIGEV)
where ISUB is given by PYEVNT, while SIGEV is to be calculated (see
below) and returned to PYEVNT. If there is only one user-defined process
then the ISUB input is superfluous, else it is necessary to branch to
the relevant process.

The SIGEV variable is supposed to give the differential cross-section
ot the current event times the phase space volume within which events
are generated, expressed in millibarn. This means that, in the limit
that many events are generated, the average value of SIGEV gives the
total cross-section of the simulated process. The SIGMAX value, handed
to PYTHIA in the PYUPIN call, is assumed to be the maximum value that
SIGEV will reach. Events will be accepted with a probability
SIGEV/SIGMAX, i.e. the acceptance/rejection of events according to
differential cross-section is done by PYEVNT, not by the user.
This means that the events that come out in the end all have unit
weight, i.e. the user does not have to worry about events with
different weights. It also allows several subprocesses to  be
generated together, in the proper mixture.

Of course, the tricky part is that the differential cross-section
usually is strongly peaked in a few regions of the phase-space, such
that the average probability to accept an event, <SIGEV>/SIGMAX is
small. It may then be necessary to find a suitable set of transformed
phase space coordinates, for which the correspondingly transformed
differential cross-section is more well-behaved.

To avoid unclarities, here a more formal version of the two paragraphs
above. Call dX the differential phase space, e.g. for a 2 -> 2 process
dX = dx_1 dx_2 d(t-hat), where x_1 and x_2 are the momentum fractions
carried by the two incoming partons and t-hat the Mandelstam variable
of the scattering. Call d(sigma)/dX the differential cross-section
of the process, e.g. for 2 -> 2: d(sigma)/dX =
sum_ij f_i(x_1,Q2) * f_j(x_2,Q2) * d(sigma-hat_ij)/d(t-hat),
i.e. the product of structure functions and hard scattering matrix
elements, summed over all allowed incoming flavours i and j.
The physical cross-section that one then wants to generate is
sigma = integral d(sigma)/dX * dX, where the integral is over the
allowed phase space volume. The event generation procedure consists
in selecting an X uniformly in dX and then evaluating the weight
d(sigma)/dX in this point. SIGEV is now simply
SIGEV = d(sigma)/dX * integral dX, i.e. the differential cross-section
times the considered volume of phase space. Clearly, when averaged over
many events, SIGEV will correctly estimate the desired cross-section. If
SIGEV fluctuates too much, one may try to transform to new variables X',
where events now are picked accordingly to dX' and
SIGEV = d(sigma)/dX' * integral dX'.

If you do not know how big SIGMAX is, you can put it to some very small
value (but bigger than zero, however) and do an exploratory run. When
the program encounters events with SIGEV > SIGMAX, a warning message is
printed which gives the new SIGMAX, that the program will use from then
on. Hopefully such maximum violations only appear at the beginning of
the run, and later it stabilizes to a level which then can be used as
input for a second, correct run.

If you want to do the event rejection yourself, simply put SIGEV
equal to SIGMAX. In that case events will not be rejected by PYTHIA
(except if there is something else wrong with them). If SIGMAX is the
correct total cross-section of the process, event mixing with other
processes will still work fine. You could also decide not to reject
any events, but to use weighted events. In that case you can only have
one ISUB switched on in a run, since the program will not know how to
mix different kinds of events, and you can not use PYTHIA to do
cross-section statistics for you. Therefore you could e.g. put
SIGMAX = SIGEV = 1, and use a commonblock to transfer event weight
and other information from your PYUPEV routine to the main program.

In addition to the SIGEV value returned for each event, it is also
necessary to return the event itself. This is done via the commonblock
      COMMON/PYUPPR/NUP,KUP(20,7),PUP(20,5),NFUP,IFUP(10,2),Q2UP(0:10)

The first part closely parallels the standard event record in the
LUJETS commonblock (see the JETSET manual or section 2.7),
although with a few simplifications. The number NUP gives the number of
particles involved in the process, where a particle may be a quark,
a lepton, a gauge boson, or anything else. The first two are simply
the two incoming particles that initiate the hard scattering, while the
remaining NUP-2 are the outgoing particles from the hard process.
For each particle I, with 1 <= I <= NUP, the following information is
stored:
KUP(I,1) : is = 1 normally. However, if you put it = 2 that signifies
    intermediate states that are not to be treated by PYTHIA, but are
    included only to make the event record more easy to read.
KUP(I,2) : is the flavour code of a particle, i.e. the two incoming
    partons for I=1 and 2, and the outgoing particles for I >= 3.
    The flavour codes are the standard PDG ones, as used elsewhere in
    the program.
KUP(I,3) : may be used to indicate the position of a mother. Such
    information may again make the record more readable, but is not
    really needed, and so one may well put all KUP(I,3)=0.
KUP(I,4) : for a final state parton which carries colour, KUP(I,4)
    gives the position of the parton from which the colour comes;
    else it must be 0.
KUP(I,5) : for a final state parton which carries anticolour, KUP(I,5)
    gives the position of the parton from which the anticolour comes;
    else it must be 0.
KUP(I,6) : for an initial state parton which carries colour, KUP(I,6)
    gives the position of the parton to which the colour goes;
    else it must be 0.
KUP(I,7) : for an initial state parton which carries anticolour,
    KUP(I,7) gives the position of the parton to which the anticolour
    goes; else it must be 0.
PUP(I,1) : p_x, i.e. x momentum.
PUP(I,2) : p_y, i.e. y momentum.
PUP(I,3) : p_z, i.e. z momentum.
PUP(I,4) : E, i.e. energy.
PUP(I,5) : m, i.e. mass.
After this brief summary, we proceed with more details and examples.

To illustrate the issue of documentation in KUP(I,1) and KUP(I,3),
consider the case of W+ production and decay to u + dbar, maybe as
part of a more complex process. The final state particles clearly
are u and dbar, so the W+ need not be given at all, but it would
make the event listing more easy to read. Therefore one should add
the W, but with KUP(I,1)=2. (If the W would have been added with
KUP(I,1)=1, it would later have been treated by PYTHIA/JETSET, which
means it would have been allowed to decay once more.) If the W+ is in
line 3, the u in 4 and the dbar in 5, one could further put KUP(4,3)=3
and KUP(5,3)=3 to indicate that the u and dbar in lines 4 and 5 come
from the W+ stored in line 3.

The colour flow information for coloured particles (quarks, gluons,
leptoquarks, ...) is needed to set up parton showers and
fragmentation properly. Sometimes many different colour flows are
possible for one and the same process, as discussed elsewhere. It is
up to you whether or not you will include all possible colour flows
in the appropriate mixture, but at least you must pick some
representative colour configuration. Consider e.g. the case of
g(1) + g(2) -> q(3) + qbar(4), where the numbers give the position
in the array. It is clear the q must get its colour from either of
the two gluons, which means there are (at least) two possibilities.
Picking it to come from gluon 1, one would thus write KUP(3,4)=1,
to be read 'the colour of parton 3 comes from parton 1'. By implication
therefore also KUP(1,6)=3, i.e. 'the colour of parton 1 goes to parton
3', i.e. colour flow is bookkept doubly. The anticolour now must
flow from parton 2 to parton 4, i.e. KUP(2,7)=4 and KUP(4,5)=2. This
completely specifies the colours of the q and qbar, but not of the
two gluons. In fact, one colour in the initial state 'annihilates'
between the g(1) and g(2), i.e. the anticolour of gluon 1 and the
colour of gluon 2 match, which may be expressed by KUP(1,7)=2 and
KUP(2,6)=1. In other words colour/anticolour of an initial state
parton may either go to a final state parton or to another initial
state parton. Correspondingly, the colour/anticolour of a final
state parton may come either from an initial state parton or from
another final state parton. An example of the latter possibility is
W decays, or generically the decay of any colour singlet particle.
(Thus a third colour flow above is represented by
g + g -> H0 -> q + qbar, where no colour passes through the Higgs,
and therefore colour flows between the two gluons and, separately,
between the q and qbar.)

Storing of momenta should be straightforward, but a few comments must
be made. Even if you in the PYINIT call ask to have events generated
in a fixed target or a user-specified frame, at intermediate stages
PYTHIA will still work in the cm frame of the two incoming beam
particles, with the first beam moving in the +z direction and the
second in the -z one. This cm frame must also be used when giving
the momenta of the process. In addition, the two incoming partons in
lines 1 and 2 are assumed massless. Therefore the initial state
partons are characterized only by the two energies P(1,4) and P(2,4),
with P(1,3)=P(1,4), P(2,3)=-P(2,4), and everything else zero. In the
final state, energies, momenta and masses are free, but must add up to
give the same four-momentum as in the initial state. All momenta are
given in GeV, with speed of light c = 1.

The second part of the PYUPPR commonblock is used to regulate the
initial and final state showering, as follows:
Q2UP(0) : Q2 scale of initial state showers.
NFUP : number of parton pairs that undergo final state showers.
IFUP(IF,1), IFUP(IF,2) : positions of the two partons of a final state
    showering pair, where the index IF runs between 1 and NFUP.
Q2UP(IF) : the Q2 scale of the final state shower between parton
    pair IF above.

If you do not want any showering at all, you can put MSTP(61)=0 and
MSTP(71)=0, and then you do not have to give the quantities above.
In general the scale choices Q2 are not unique, which means some
guesswork is involved. Since the showers add additional partonic
activity at mass-squared scales below the Q2 choices above, the
Q2UP should be of the order of the phase space cutoffs, so as to provide
a reasonably smooth joining between partonic activity from matrix
elements and that from showers. There are a few cases where choices
are rather easy. In the decay of any s-channel colour neutral state,
such as a W, the Q2 scale of final state showers is just set by the
mass-squared of the resonance. For initial state radiation, Q2UP(0)
should be about the same as the Q2 scale used for the evaluation of
structure functions for the hard process, up to some factor of order
unity. (One frequent choice for this factor would be 4, if your
structure function scale is something like squared transverse momentum,
simply because m2 is of order 4 pT2.)

The 'parton' shower evolution actually also can include photon emission
off quarks and leptons, if the shower switches are properly set. It is
not possible to define only one particle in the arrays above, since
then it would not be possible to conserve energy and momentum in the
shower. You can very well have a pair where only one of the two can
branch, however. For instance, in a g + gamma final state, only
the gluon can shower, but the photon can lose energy to the gamma
in such a way that the gluon branchings becomes possible.

Currently, it is not possible to do showering where three or more final
state particles are involved at the same time. This may be added at a
later stage. It is therefore necessary to subdivide suitably into pairs,
and maybe leave some (especially colour neutral) particles unshowered.

You are free to make use of whatever tools you want in your PYUPEV
routine, and normally there would be no contact with the rest of
PYTHIA, except as described above. However, you are free to make use
of some of the tools already available. One attractive possibility
is to use PYSTFU for structure function evaluation; that way the
standard PYTHIA switches could be used to switch between different
parametrizations. Other possible tools could be RLU for random number
generation, ULALPS for alpha_strong evaluation, ULALEM for evaluation
of a running alpha_em, and maybe a few more.

**********************************************************************

      References

Alt89   G. Altarelli, B. Mele, M. Ruiz-Altaba, Z. Physik C45 (1989) 109

Bai83   R. Baier, R. Ruckl, Z. Physik 19 (1983) 251

Bar90   T.L. Barklow, SLAC-PUB-5364 (1990)

Bau90   U. Baur, M. Spira, P. M. Zerwas, Phys. Rev. D42 (1990) 815

Ben84   H.-U. Bengtsson, Computer Phys. Comm. 31 (1984) 323

Ben84a  H.-U. Bengtsson, G. Ingelman, LU TP 84-3, Ref.TH.3820-CERN

Ben85   H.-U. Bengtsson, G. Ingelman, Computer Phys. Comm. 34 (1985) 251

Ben85a  H.-U. Bengtsson, W.-S. Hou, A. Soni, D.H. Stork,
        Phys. Rev. Lett. 55 (1985) 2762

Ben87   H.-U. Bengtsson, T. Sjostrand, Computer Phys. Comm. 46 (1987) 43

Ber84   E. L. Berger, E. Braaten, R. D. Field, Nucl. Phys.
        B239 (1984) 52

Blo85   M. M. Block, R. N. Cahn, Rev. Mod. Phys. 57 (1985) 563
        M. M. Block, R. N. Cahn, in Physics Simulations at High Energy,
        eds. V. Barger, T. Gottschalk, F. Halzen (World Scientific,
        Singapore, 1987), p. 89

BLS87   M. C. Bento, C. H. Llewellyn Smith, Nucl. Phys. B289 (1987) 36

Cah84   R.N. Cahn, S. Dawson, Phys. Lett. 136B (1984) 196
        R.N. Cahn, Nucl. Phys. B255 (1985) 341
        G. Altarelli, B. Mele, F. Pitolli, Nucl. Phys. B287 (1987) 205

Cha85   M. Chanowitz, M.K. Gaillard, Nucl. Phys. B261 (1985) 379

Cha91   K. Charchula, DESY 91-093 (1991)

Cha91a  K. Charchula, ZEUS note 91-73 (1991)

Che75   M.-S. Chen, P. Zerwas, Phys. Rev. D12 (1975) 187;
        P. Zerwas, private communication (1991)

Coc90   D. Cocolicchio, F. Feruglio, G.L. Fogli, J. Terron,
        CERN-TH.5909/90
        F. Feruglio, private communication

Com77   B. L. Combridge, J. Kripfganz, J. Ranft, Phys. Lett. 70B (1977)
        234

Con71   V. Constantini, B. de Tollis, G. Pistoni, Nuovo Cim. 2A (1971)
        733

Cut78   R. Cutler, D. Sivers, Phys. Rev. D17 (1978) 196

DFL88   M. Diemoz, F. Ferroni, E. Longo, G. Martinelli, Z. Physik C39
        (1988) 21

DHT90   A. Dobado, M. J. Herrero, J. Terron, CERN-TH.5670/90 and
        FTUAM/08-90, CERN-TH.5813/90

Dic88   D. A. Dicus, S. S. D. Willenbrock, Phys. Rev. D37 (1988) 1801

Dre85   M. Drees, K. Grassie, Z. Physik C28 (1985) 451

Dre89   M. Drees, J. Ellis, D. Zeppenfeld, Phys. Lett. B223 (1989) 454

Duk82   D.W. Duke, J.F. Owens, Phys. Rev. D26 (1982) 1600

Dun86   M. J. Duncan, G. L. Kane, W. W. Repko, Nucl. Phys. B272 (1986)
        517

EHL84   E. Eichten, I. Hinchliffe, K. Lane, C. Quigg, Rev. Mod. Phys.
        56 (1984) 579; 58 (1985) 1065

Ell86   R. K. Ellis, J. C. Sexton, Nucl. Phys. B269 (1986) 445

Ell88   R. K. Ellis, I. Hinchliffe, M. Soldate, J. J. van der Bij,
        Nucl. Phys. B297 (1988) 221

Fad89   V. Fadin, V. Khoze, JETP Lett. 46 (1987) 417, Yad. Fiz 48
        (1988) 487, in Proceedings of the 24th Winter School of the
        LNPI (Leningrad, 1989), vol. I, p. 3
        V. Fadin, V, Khoze, T. Sjostrand, Z. Phys. C48 (1990) 613

Fon81   M. Fontannaz, B. Pire, D. Schiff, Z. Phys. C11 (1981) 211

Gab86   E. Gabrielli, Mod. Phys. Lett. A1 (1986) 465

Gas87   R. Gastmans, W. Troost, T.T. Wu, Phys. Lett. B184 (1987) 257

GHK90   J.F. Gunion, H.E. Haber, G. Kane, S. Dawson, The Higgs Hunter's
        Guide (Addison-Wesley, 1990)
        A. Djouadi, private communication

Glo88   E.W.N. Glover, A.D. Martin, W.J. Stirling, Z. Physik C38
        (1988) 473

Gou83   K. Goulianos, Phys. Rep. 101 (1983) 169

GRV89   M. Gluck, E. Reya, A. Vogt, DO-TH 89/20 (extended version)

Gun86   J.F. Gunion, Z. Kunszt, Phys. Rev. D33 (1986) 665,
        errata as private communication from the authors

Hag90   K. Hagiwara, H. Iwasaki, A. Miyamoto, H. Murayama,
        D. Zeppenfeld, Durham preprint DTP/90/82 etc. (1990)

Hal78   F. Halzen, D. M. Scott, Phys. Rev. D18 (1978) 3378

Hew88   J.L. Hewett, S. Pakvasa, Phys. Rev. D37 (1988) 3165
        and private communication from the authors

Ing87   G. Ingelman et al., in Proceedings of the HERA Workshop,
        ed. R.D. Peccei, Vol. 1, p. 3

Kat83   M. Katuya, Phys. Lett. 124B (1983) 421

Kle89   R. Kleiss et al., in Z Physics at LEP 1, eds. G. Altarelli,
        R. Kleiss, C. Verzegnassi, CERN 89-08, Vol. 3, p.1

Kle90   R. Kleiss, private communication (1990)

Kun84   Z. Kunszt, Nucl. Phys. B247 (1984) 339

Kun87   Z. Kunszt, CERN 87-07, Vol. 1, p. 123, and private communication

Lan90   K. Lane, private communication (1991)

Mor90   J.G. Morfin, W.-K. Tung, Fermilab-Pub-90/74 and IIT-PHY-90/11

PDG88   Particle Data Group, G. P. Yost et al.,
        Phys. Lett. 204B (1988) 1

Plo91   H. Plothow-Besch, in Proceedings of the 3rd Workshop on Detector
        and Event Simulation in High Energy Physics, Amsterdam,
        8-12 April 1991

Sam90   M.A. Samuel, G. Li, N. Sinha, R. Sinha, M.K. Sundaresan,
        Carleton preprint (1990)

Sjo87   T. Sjostrand, M. van Zijl, Phys. Rev. D36 (1987) 2019

Sjo90   T. Sjostrand, JETSET 7.3 program and manual (unpublished,
        obtainable on request)

Ter90   J. Terron, private communication (paper in preparation)

Tun87   W.-K. Tung, in Physics Simulations at High Energy,
        eds. V. Barger, T. Gottschalk, F. Halzen (World Scientific,
        Singapore, 1987), p. 601
        W.-K. Tung, private communication

Wud86   J. Wudka, Phys. Lett. 167B (1986) 337

Zer90   P. Zerwas, private communication

**********************************************************************