********************************************************************** ********************************************************************** ** ** ** 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 = (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 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 = 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) : = (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, /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 **********************************************************************