14 September 1995 Updates to PYTHIA 5.7 and JETSET 7.4 Physics and Manual Torbjorn Sjostrand Department of theoretical physics 2 University of Lund Solvegatan 14A S-223 62 Lund Sweden ----------------------------------------------------------------- 1) Introduction. Since the PYTHIA/JETSET programs are being updated frequently, it is also important to keep the documentation up to date. The big manual that appeared as CERN-TH.7112/93 described the status as of 15 December 1993. The LaTeX file with this manual will be updated, though less frequently than the programs. Further, it is not very economical to have to get a new copy of a long manual just to see if any interesting new features have been added recently. Therefore I have here collected the main changes that have taken place since the beginning of 1994. The changes are indexed by sub-subversion number and date. This file will be updated as regularly as the programs themselves. It is not intended to be as complete as the ordinary manual, but should be sufficient for the intended purpose. ----------------------------------------------------------------- 2) Changes in JETSET 7.4 ----------- 00, 13 December 1993: baseline version. ----------- 01, 11 February 1994: LUZDIS has been changed to protect against an overflow in an exponent (harmless physicswise, but enough to crash the program on some machines). ----------- 02, 7 April 1994: A possibility has been introduced into LUSHOW to suppress either hard or soft radiation, in the connection of radiation off rapidly decaying objects. The algorithm used is not exact, but still gives some impression of potential effects. The switch ought to have appeared at the end of the current list of shower switches, but because of lack of space it appears immediately before. MSTJ(40) : (D=0) possibility to suppress the branching probability for a branching q -> q + g (or q -> q + gamma) of a quark produced in the decay of an unstable particle with width Gamma, where this width has to be specified by the user in PARJ(89). Can be changed for each new LUSHOW call. = 0 : no suppression, i.e. the standard parton-shower machinery. = 1 : suppress radiation by a factor chi(omega) = Gamma**2 / (Gamma**2 + omega**2), where omega is the energy of the gluon (or photon) in the rest frame of the radiating dipole. Essentially this means that hard radiation with omega > Gamma is removed. = 2 : suppress radiation by a factor 1 - chi(omega) = omega**2 / (Gamma**2 + omega**2), where omega is the energy of the gluon (or photon) in the rest frame of the radiating dipole. Essentially this means that soft radiation with omega < Gamma is removed. PARJ(89) : (D=0. GeV) the width of the unstable particle studied for the MSTJ(40) > 0 options; to be set by the user (separately for each LUSHOW call, if need be). ----- A generic interface has been included to an external tau decay library. This should allow the handling of tau polarization, which is not done by JETSET. To use this facility you have to set the switch MSTJ(28), include your own interface routine LUTAUD and see to it that the dummy routine LUTAUD in JETSET is not linked. The dummy routine is there only to avoid unresolved external references when no user-supplied interface is linked. MSTJ(28) : (D=0) call to an external tau decay library. = 0 : not done, i.e. the internal LUDECY treatment is used. = 1 : done whenever the tau mother particle species can be identified, else the internal LUDECY treatment is used. Normally the mother particle should always be identified, but it is possible for a user to remove event history information or to add extra taus directly to the event record, and then the mother is not known. = 2 : always done. CALL LUTAUD(ITAU,IORIG,KFORIG,NDECAY) Purpose: to act as an interface between the generic decay routine LUDECY and a user-supplied tau decay library. The latter library would normally know how to handle polarized taus, given the tau polarization, so one task of the interface routine is to construct the tau polarization/helicity from the information available. Input to the routine is provided in the first three arguments, while the last argument and some event record information have to be set before return. ITAU : line number in the event record where the tau is stored. The four-momentum of this tau has first been boosted back to the rest frame of the decaying mother and thereafter rotated to move out along the +z axis. This choice of frame should help the calculation of the helicity configuration. After the LUTAUD call the tau (and its decay products) will be rotated and boosted back. However, seemingly, the event record does not conserve momentum at this intermediate stage. IORIG : line number where the mother particle to the tau is stored. Is 0 if the mother is not stored. This does not have to mean the mother is unknown, e.g. in semileptonic B decay the mother is a W+-, and its momentum can be obtained by adding the tau and nu_tau momentum, but there is no line in the event record. When several copies of the mother is stored (e.g. one in the documentation section of the event record and one in the main section), IORIG points to the last. If a branchings like tau -> tau + gamma occurs, the 'grandmother' is given, i.e. the mother of the direct tau before branching. KFORIG : flavour code for the mother particle. Is 0 if the mother is unknown. The mother would typically be a resonance such as Z0 (23), W+- (+-24), H0 (25), or H+- (+-37). Often the helicity choice would be clear just by the knowledge of this mother species, e.g. W+- vs. H+-. However, sometimes further complications may exist. For instance, the KF code 23 represents a mixture of gamma* and Z0; a knowledge of the mother mass (in P(IORIG,5)) would here be required to make the choice of helicities. Further, a W and Z may either be (predominantly) transverse or longitudinal, depending on the production process under study. NDECAY : the number of decay products of the tau; to be given by the user. You must also store the KF flavour codes of those decay products in the positions K(I,2), N+1 <= I <= N+NDECAY, of the event record. The corresponding five-momentum (momentum, energy and mass) should be stored in the associated P(I,J) positions, 1 <= J <= 5. The four-momenta are expected to add up to the four-momentum of the tau in position ITAU. You should not change the N value or any of the other K or V values (neither for the tau nor for its decay products) since this is automatically done in LUDECY. ----- In a few places, a dot has been moved from the end of one line to the beginning of the next continuation line, or the other way around, to keep together tokens such as .EQ. or .AND., since some debuggers may otherwise complain. A source of (harmless) division by zero in LUSHOW has been removed. ----------- 03, 15 July 1994: The LUBOEI routine has been changed to avoid an unintentional gap in the limits of the very first bin. Further, leptons and photons which are unrelated to the system feeling the Bose-Einstein effects do not have their energies and momenta changed in the global rescaling step. (Example: W+W- events, where one W decays leptonically; before these lepton momenta could be slightly changed, but now not.) ----- The option LUEDIT(16) (used e.g. from PYEVNT) has been improved with a more extensive search for missing daughter pointers. ----- The KLU(I,16) procedure for finding rank has been rewritten to work in the current JETSET version, which it did not before. However, note that it will only work for MSTU(16)=2. As a general comment, the options 14 - 17 of KLU were written at a time when possible event histories were less complex, and can not be guaranteed always to work today. ----------- 04, 25 August 1994: LUSHOW has been corrected, so that if t, l or h quarks (or d* or u* quarks masked as l or h) are given with masses that vary from event to event (a Breit-Wigner shape, e.g.), the current mass rather than the nominal mass is used to define the cut-off scales of parton shower evolution. ----- LULOGO has been modified to take into account that a new PYTHIA/JETSET description has been published in T. Sjostrand, Computer Phys. Commun. 82 (1994) 74 and is from now on the standard reference to these two programs. ----------- 05, 27 January 1995: LUCELL has been corrected, in that in the option with smearing of energy rather than transverse energy, the conversion factor between the two was applied in the wrong direction. ----- LUSHOW has been corrected in one place where the PMTH array was addressed with the wrong order of the indices. This affected quark mass corrections in the matching to the three-jet matrix elements. ----- An additional check has been included in LUBOEI that there are at least two particles involved in the Bose-Einstein effects. (No problem except in some bizarre situations.) ----------- 06, 20 February 1995: A new option has been added for the behaviour of the running alpha_em(Q2) in ULALEM. This is not added as a true physics scenario, but only to produce results with a given, fixed value for the hard events, while still keeping the conventional value in the Q2=0 limit. MSTU(101) = 2 : if Q2 is less than PARU(104) then alpha_em is assigned the value PARU(101) (=1/137), while for Q2 above PARU(104) the fixed value PARU(103) (=1/128.8) is used. PARU(103) : (D=0.007764=1/128.8) alpha_em used for hard processes in the option MSTU(101)=2. PARU(104) : (D=1 GeV^2) dividing line for 'low' and 'high' Q2 values in the MSTU(101)=2 option of ULALEM. Additionally, the G_F constant has been added to the parameter list. PARU(105) : (D=1.16639E-5 GeV^-2) G_F, the Fermi constant of weak interactions. ----- The LULOGO routine has been updated to reflect my change of affiliation. ----------- 07, 21 June 1995: Header and LULOGO have been updated with respect to phone number and WWW access. ----- The PHEP and VHEP variables in the /HEPEVT/ common block are now assumed to be in DOUBLE PRECISION, in accord with the proposed LEP 2 workshop addendum to the standard. ----- In LUTEST a missing decimal point on the energy check has been reinstated (0001 -> 0.0001). ----- In LUINDF the expression PR/(Z*W) has been protected against vanishing denominator. ----------- 08, 23 August 1995: Check against division by zero in LUSHOW. ----------------------------------------------------------------- 3) Changes in PYTHIA 5.7 ----------- 00, 13 December 1993: baseline version. ----------- 01, 27 January 1994: The machinery to handle gamma-gamma interactions is expanded. In particular, several new options have been added to MSTP(14). The updated description of this variable reads as follows. MSTP(14) : (D=0) structure of incoming photon beam or target (does not affect photon inside electron, only photons appearing as argument in the PYINIT call). = 0 : a photon is assumed to be point-like (a direct photon), i.e. can only interact in processes which explicitly contain the incoming photon, such as f_i gamma -> f_i g for gamma-p interactions. In gamma-gamma interactions both photons are direct, i.e the main process is gamma gamma -> f_i fbar_i. = 1 : a photon is assumed to be resolved, i.e. can only interact through its constituent quarks and gluons, giving either high-pT parton-parton scatterings or low-pT events. Hard processes are calculated with the use of the full photon parton distributions. In gamma-gamma interactions both photons are resolved. = 2 : a photon is assumed resolved, but only the VMD piece is included in the parton distributions, which therefore mainly are scaled-down versions of the rho0/pi0 ones. Both high-pT parton-parton scatterings and low-pT events are allowed. In gamma-gamma interactions both photons are VMD-like. = 3 : a photon is assumed resolved, but only the anomalous piece of the photon parton distributions is included. Only high-pT parton-parton scatterings are allowed. In gamma-gamma interactions both photons are anomalous. = 4 : in gamma-gamma interactions one photon is direct and the other resolved. A typical process is thus f_i gamma -> f_i g. Hard processes are calculated with the use of the full photon parton distributions for the resolved photon. Both possibilities of which photon is direct are included, in event topologies and in cross sections. This option cannot be used in configurations with only one incoming photon. = 5 : in gamma-gamma interactions one photon is direct and the other VMD-like. Both possibilities of which photon is direct are included, in event topologies and in cross sections. This option cannot be used in configurations with only one incoming photon. = 6 : in gamma-gamma interactions one photon is direct and the other anomalous. Both possibilities of which photon is direct are included, in event topologies and in cross sections. This option cannot be used in configurations with only one incoming photon. = 7 : in gamma-gamma interactions one photon is VMD-like and the other anomalous. Only high-pT parton-parton scatterings are allowed. Both possibilities of which photon is VMD-like are included, in event topologies and in cross sections. This option cannot be used in configurations with only one incoming photon. Note: a complete description requires separate runs for the components above, i.e. it is not possible to mix them in a single run. Our best understanding of gamma-p interactions [Sch93,Sch93a] is to have three separate components, 0 + 2 + 3. A simpler alternative is based on two only, 0 + 1. Our best understanding of gamma-gamma interactions [in preparation] requires six separate components, 0 + 2 + 3 + 5 + 6 + 7. A simpler alternative is based on three only, 0 + 1 + 4. In addition, one new option has been introduced and a few internal variables modified. MSTP(59) : (D=0) possibility to modify the Q2 scale used in the anomalous parton distributions of the photon, as used in the options MSTP(14) = 3, 6 and 7. = 0 : no change of Q2 scale compared to what is normally used. = 1 : the input Q2 scaled is divided by PARP(59)**2 to define the Q2 scale used as argument for the anomalous parton distributions. PARP(59) : (D=1.) rescaling factor used for the Q2 argument of the anomalous parton distributions of the photon, see MSTP(59). MINT(105) : is MINT(103) or MINT(104), depending on which side of the event currently is being studied. MINT(107), MINT(108) : if either or both of the two incoming particles is a photon, then the respective value gives the nature assumed for that photon. The code follows the one used for MSTP(14): = 0 : direct photon. = 1 : resolved photon. = 2 : VMD-like photon. = 3 : anomalous photon. MINT(109) : is either MINT(107) or MINT(108), depending on which side of the event currently is being studied. VINT(282) : no longer used. VINT(283), VINT(284): virtuality scale at which an anomalous photon on the beam or target side of the event is being resolved. More precisely, it gives the p_T^2 pf the gamma -> q qbar vertex. ----- A number of bugs have also been corrected: * Jet + low-pT event generation could give incorrect cross section information with PYSTAT(1) at low energies. The event generation itself is correct. (The error was introduced when variable energies became allowed.) * Introduce rejection of top events where top mass (in the tails of the Breit-Wigner distribution) is too low to allow decays t -> W + b. * Plus a few minor bugs, probably harmless. ----------- 02, 13 February 1994: The interface to PDFLIB has been modified to reflect that 'TMAS' should no longer be set except in first PDFSET call. (Else a huge amount of irrelevant warning messages are generated by PDFLIB.) ----- The STOP statement in a few dummy routines has been modifed to avoid irrelevant compilation warning messages on IBM mainframes. ----- A few labels have been renumbered. ----------- 03, 22 February 1994: Removal of a bug in PYRESD, which could give (under some specific conditions) errors in the colour flow. ----------- 04, 7 April 1994: Process 11 has been corrected, for the part that concerns anomalous couplings (contact interactions) in the q + q' -> q + q' process. The error was present in the expression for u + dbar -> u + dbar and obvious permutations, while u + d -> u + d, u + ubar -> u + ubar and the others were correct. Thanks to J.-J. Dugne, M. Perrottet and K. Lane for communications on this point. ----- The option MSTP(23)=1 for post-facto (x,Q^2) conservation in deep inelastic scattering can give infinite loops when applied to process 83, in particular if one asks for the production of a top. (Remember that the standard DIS kinematics only is defined for massless quarks.) Therefore the switch MSTP(23) has been modifed as follows: MSTP(23) : (D=1) (x, Q^2) correction level in DIS. = 0 : no correction procedure applied. = 1 : correction applied for process 10, but not for process 83. = 2 : correction applied both for process 10 and 83. This latter option could still work fine for charm and bottom, if the energy is sufficient. ----- PYRESD is modified to ensure isotropic angular distributions in the decays of the top or a fourth generation particle, i.e. in t -> b + W+. This may not be the correct distribution but, unless explicit knowledge exists for a given process, this should always be the default. ----- In processes 16, 20, 31 and 36 the W propagator has been modified to include s-dependent widths in the Breit-Wigner shape. The most notable effect is a suppression of the low-mass tail of the W mass spectrum. ----- When PDFLIB is used, PDFSET is now only called whenever a different structure function is requested. For pp events therefore only one call is made, while gamma-p interactions still involves a call to PDFSET for each STRUCTM one, since gamma and p structure functions have to be called alternatingly. To this end, MINT(93) is reset to 1000000 * Nptype + 1000 * Ngroup + Nset after each PDFSET call. ----- In a few places, a dot has been moved from the end of one line to the beginning of the next continuation line, or the other way around, to keep together tokens such as .EQ. or .AND., since some debuggers may otherwise complain. Also some other purely cosmetics changes for the same reason. A number of minor errors have been corrected. ----------- 05, 15 July 1994: A new option has been introduced, MSTP(14)=10, whereby it is possible to obtain a mixture of the various allowed photon components. For gamma-hadron collisions, this means a mixture of VMD, direct and anomalous events, for gamma-gamma collisions a mixture of VMD*VMD, VMD*direct, VMD*anomalous, direct*direct, direct*anomalous and anomalous*anomalous. The mixture is properly given according to the relative cross sections. Note that this introduces a completely new layer of administration in PYTHIA. For instance, a subprocess such as q + g -> q + g is allowed in the VMD*VMD, VMD*anomalous and anomalous*anomalous classes, but appear with different sets of parton distributions and with different pT cut-offs. In order to handle this, various information is initialized separately for each event class, and subsequently saved and restored as the generation switches back and forth between the event classes. This introduces some limitations on what you may and may not do. First of all, the MSTP(14) switch is only applicable for incoming photon beams, i.e. when 'gamma' is the argument in the PYINIT call. A convolution with the bremsstrahlung photon spectrum in an electron beam may come one day, but not in the immediate future. Secondly, the machinery has only been set up to generate standard QCD physics, specifically either 'minimum bias' one or high-pT jets. For minimum bias, you are not allowed to use the CKIN variables at all. This is not a major limitation, since it is in the spirit of minimum bias physics not to impose any contraints on allowed jet production. (If you still do, these cuts will be ineffective for the VMD processes but take effect for the other ones, giving inconsistencies.) Further, some variables are internally recalculated and reset: CKIN(1), CKIN(3), CKIN(5), CKIN(6), MSTP(57), MSTP(85), PARP(2), PARP(81), PARP(82), PARU(115) and MDME(22,J). These can not be modified without changing PYINPR and recompiling the program. The minimum bias physics option is obtained by default; by switching from MSEL=1 to MSEL=2 also the elastic and diffractive components of the VMD part are included. High-pT jet production is obtained by setting the CKIN(3) cut-off larger than the (energy-dependent) cut-off scales for the VMD, direct and anomalous components; typically this means at least 3 GeV. For lower input CKIN(3) the program will automatically switch back to minimum bias physics. Finally, pileup events are not at all allowed. Here is a survey of common block variables affected: MSTP(14) (D=0) strucure of incoming photon beam or target; see description above for PYTHIA 5.701. = 10 : new option where the VMD, direct and anomalous components are automatically mixed, as described above. Works equally well for gamma-p and gamma-gamma. MSTI(9) : event class used in current event. = 1 : VMD (for gamma-p) or VMD*VMD (for gamma-gamma). = 2 : direct (for gamma-p) or VMD*direct (for gamma-gamma). = 3 : anomalous (for gamma-p) or VMD*anomalous (for gamma-gamma). = 4 : direct*direct (for gamma-gamma). = 5 : direct*anomalous (for gamma-gamma). = 6 : anomalous*anomalous (for gamma-gamma). MINT(121) : number of separate event classes to initialize and mix. = 1 : the normal value. = 3 : for a gamma-hadron interaction when MSTP(14)=10. = 6 : for a gamma-gamma interaction when MSTP(14)=10. MINT(122) : event class used in current event. Code as explained for MSTI(9). MINT(123) : event class used in the current event, with the same list of possibilities as MSTP(14), except that MSTP(14) = 1, 4 or 10 do not appear. VINT(285) : the CKIN(3) value provided by the user at initialization; subsequently CKIN(3) may be overwritten (for MSTP(14)=10) but VINT(285) stays. In addition, the structure of the initialization has been partly reorganized. The routine PYEVKI has been removed, new routines PYINBM, PYINPR and PYSAVE created, and some material has been moved to or from PYINIT, PYINRE and PYINKI. SUBROUTINE PYINBM : to read in and identify beam and target particles and frame as given in the PYINIT call (used to be done in PYINKI). SUBROUTINE PYINKI(MODKI) : to set up event kinematics, either at initialization (MODKI=0) or for each separate event when varying kinematics (MODKI=1). (The latter task used to be done in PYEVKI.) SUBROUTINE PYINPR : to set up the partonic subprocesses selected with MSEL and, for gamma-p and gamma-gamma, MSTP(14). SUBROUTINE PYSAVE : saves and restores parameters and cross section values between the 3 gamma-p and the 6 gamma-gamma alternatives of MSTP(14)=10. Also makes a random choice for each new event between the allowed alternatives. Among other changes, note that PYSTAT(1) now has been extended so that the subdivision into the various gamma-p and gamma-gamma classes is shown. ----- Further changes of particular relevance for gamma-p and gamma-gamma, but independent of the major revisions above: MSTP(59) and PARP(59) have been removed. Instead the following options are available: MSTP(15) : (D=5) possibility to modify the nature of the anomalous photon component, in particular with respect to the scale choices and cut-offs of hard processes. = 0 : none, i.e. the same treatment as for the VMD component. = 1 : evaluate the anomalous structure functions at a scale Q2/PARP(17)^2. = 2 : as =1, but instead of PARP(17) use PARP(81)/PARP(15) or PARP(82)/PARP(15), depending on MSTP(82) value. = 3 : evaluate anomalous structure function as f^(anom)(x, Q2, p_0^2) - f^(anom)(x, Q2, r^2*Q2) with r = PARP(17). = 4 : as =3, but instead of PARP(17) use PARP(81)/PARP(15) or PARP(82)/PARP(15), depending on MSTP(82) value. = 5 : use larger pTmin for anomalous component than for VMD one, but otherwise no difference. PARP(17) : (D=1) rescaling factor used as described for MSTP(15). MSTP(51) : new option added. = 11 : the GRV p LO parametrization. MSTP(53) : new option added. = 3 : the GRV pi LO parametrization. MSTP(56) : new option added. = 3 : when the anomalous photon structure function is requested, the homogeneous solution is provided, evolved from a starting value PARP(15) to the requested Q scale. The homogeneous solution is normalized so that the net momentum is unity, i.e. any factors of alpha_em/2pi and charge have been left out. The flavour of the original q is given in MSTP(55) (1, 2, 3, 4 or 5 for d, u, s, c or b); the value 0 gives a mixture according to squared charge, with the exception that c and b are only allowed above the respective mass threshold (Q > m_q). The four-flavour Lambda value is assumed given in PARP(1); it is automatically recalculated for 3 or 5 flavours at thresholds. This option is not intended for standard event generation, but is useful for some theoretical studies. ----- Option MSTP(92)=5 for beam remnant treatment erroneously missed some statements which now have been inserted. Further, new options have been added for the splitting of momentum between two beam remnants. MSTP(92) keeps its current role for the production of diquark or quark jets. However, for the splitting into a hadron plus a quark/diquark jet, MSTP(94) should now be used. MSTP(94) : (D=2) (C) energy partitioning in hadron or resolved photon remnant, when this is split into a hadron plus a remainder-jet. The energy fraction chi is taken by one of the two objects, with conventions as described below or for PARP(95) and PARP(97). = 1 : 1 for meson or resolved photon, 2(1-chi) for baryon, i.e. simple counting rules. = 2 : (k+1)*(1-chi)**k, with k as given in PARP(95) or PARP(97). = 3 : the chi of the hadron is selected according to the normal fragmentation function used for the hadron in jet fragmentation, see MSTJ(11). The possibility of a changed fragmentation function shape in diquark fragmentation (see PARJ(45)) is not included. = 4 : as =3, but the shape is changed as allowed in diquark fragmentation (see PARJ(45)); this change is here also allowed for meson production. ----- In PYDIFF the recoiling gluon energy is calculated in a numerically more stable fashion. ----- A counter has been added to PYSSPA to avoid infinite loops in the angular ordering constraint due to interference with the final state colour charges. ----------- 06, 25 August 1994: New processes 167 and 168 have been added for the contact interaction production of d* or u* excited quarks 167 q q' -> q" d* 168 q q' -> q" u* where the different allowed quark and antiquark combinations are given according to eqs. (15) - (19) in U. Baur, M. Spira and P.M. Zerwas, Phys. Rev. D42 (1990) 815. The d* and u* are defined in the same way as for processes 147 and 148. Thus one needs to put MSTP(6)=1 to use l (7) and h (8) for representing the d* and u*. The couplings of the allowed decay channels are given by PARU(157) - PARU(159), and the Lambda scale parameter by PARU(155). At the same time, some minor changes has been introduced in the code for processes 147 and 148, for uniformity. ----- Option MSTP(57)=3 now also allows a dampening of pi+- parton distributions. ----- A few minor errors have been corrected. ----------- 07, 20 October 1994: A major bug discovered in processes 121 and 122 (and thus also affecting 181, 182, 186 and 187), g g or q qbar -> Q Qbar H: the kinematics was incorrectly handed on to the Kunszt matrix elements. This affected the default option Q = t, but effects were especially dramatic when the alternative Q = b was used. The choice of appropriate Q2 scale for structure functions introduces a further uncertainty in cross sections for the processes above. So long as only t quarks are considered, the t mass is a reasonable choice, but for the Q = b alternative this is presumably too low. Therefore new options have been introduced as below, with the default behaviour changed (the old one is obtainable with MSTP(39)=1). MSTP(39) : (D=2) choice of Q2 scale for structure functions and initial state parton showers in processes g g or q qbar -> Q Qbar H. = 1 : m_Q**2. = 2 : max(mT_Q**2 , mT_Qbar**2) = m_Q**2 + max(pT_Q**2 , pT_Qbar**2). = 3 : m_H**2. = 4 : shat = (p_H + p_Q + p_Qbar)**2. ----- Another important bug corrected in the calculation of the reduction of t+tbar cross section when decay modes are forced. This occured when both t and tbar produced a W, and W+ and W- decay modes were set differently. ----------- 08, 25 October 1994: A few further places changed to make processes 181, 182, 186 and 187 work (see version 5.707 above). ----------- 09, 26 October 1994: The matrix element for f + fbar -> W+ + W- has been replaced, using the formulae of D. Bardin, M. Bilenky, D. Lehner, A. Olchevski and T. Riemann, CERN-TH.7295/94, but with the dependence on the t-hat variable not integrated out (D. Bardin, private communication). This avoids some problems encountered in the old expressions when one or both W's were far off the mass shell. ----- Change in calls to PDFLIB, so that the input Q is always at least the QMIN of the respective set. ----- Extra protection against infinite loops in PYSSPA. ----------- 10, 27 January 1995: The dimensions of the HGZ array in PYRESD has been expanded to avoid accidental writing outside the bounds. ----- VINT(41) - VINT(66) are saved and restored in PYSCAT, for use in low-pT events, when beam remnant treatment has failed (with nonzero MINT(57)). ----- The routine PYSTGH has been replaced by the routine PYSTHG. This contains an improved parametrization of the homogeneous evolution of an anomalous photon from some given initial scale. The argument NF of the PYSTGH routine has been removed; now Lambda is always automatically converted to the relevant NF-flavour value from its 4-flavour one, at flavour thresholds. ----------- 11, 20 February 1995: New possibilities have been added to switch between electroweak couplings being expressed in terms of a running alpha_em(Q2) or in terms of a fixed Fermi constant G_F. This affects both decay widths and process cross sections, in the routines PYINRE, PYRESD, PYWIDT and PYSIGH. There are three main options, with default agreeing with the old standard. MSTP(8) : (D=0) choice of electroweak parameters to use in decay widths of resonances (W, Z, H, ...) and cross sections (production of W's, Z's, H's, ...). = 0 : everything is expressed in terms of a running alpha_em(Q2) and a fixed sin^2(theta_W), i.e. G_F is nowhere used. = 1 : a replacement is made according to alpha_em(Q2) -> sqrt(2) G_F m_W^2 sin^2(theta_W) / pi in all widths and cross sections. If G_F and m_Z are considered as given, this means that sin^2(theta_W) and m_W are the only free electroweak parameter. = 2 : a replacement is made as for =1, but additionally sin^2(theta_W) is constrained by the relation sin^2(theta_W) = 1 - m_W^2/m_Z^2 This means that m_W remains as a free parameter, but that the sin^(theta_W) value in PARU(102) is never used, EXCEPT in the vector couplings in the combination v = a - 4 sin^2(theta_W) e. This degree of freedom enters e.g. for forward-backward asymmetries in Z^0 decays. Note : This option does not affect the emission of real photons in the initial and final state, where alpha_em is always used. However, it does affect also purely electromagnetic hard processes, such as q + qbar -> gamma + gamma. ----- The option MSTP(37)=1, with running quark masses in couplings to Higgs bosons, only works when alpha_s is allowed to run (so one can define a Lambda value). Therefore a check has been introduced in PYWIDT and PYSIGH that the option MSTP(37)=1 is only executed if additionally MSTP(2) is 1 or bigger. ----- Some non-physics changes have been made in the RKBBV and STRUCTM codes so as to avoid some (in principle harmless) compiler warnings. ----------- 12, 15 March 1995: A serious error has been corrected in the MSTP(173)=1 option, i.e. when the program is run with user-defined weights that should compensate for a biased choice of variable beam energies. This both affected the relative admixture of low- and high-pT events and the total cross section obtained by Monte Carlo integration. (PYRAND changed.) In order to improve the flexibility and efficiency of the variable-energy option, the user should now set PARP(174) before the PYINIT call, and thereafter not change it. This allows PARP(173) weights of arbitrary size. (PYRAND and PYMAXI changed.) PARP(174) : (D=1.) maximum event weight that will be encountered in PARP(173) during the course of a run with MSTP(173)=1; to be used to optimize the efficiency of the event generation. It is always allowed to use a larger bound than the true one, but with a corresponding loss in efficiency. MSTI(5) (and MINT(5)) are now changed so they count the number of successfully generated events, rather than the number of tries made. This change only affects runs with variable energies, MSTP(171)=1 and MSTP(172)=2, where MSTI(61)=1 signals that a user-provided energy has been rejected in the weighting. This change also affects PARI(2), which becomes the cross section per fully generated event. (PYEVNT changed.) ----- The option MSTP(14)=10 has now been extended so that it also works for deep inelastic sacattering of an electron off a (real) photon, i.e. process ISUB = 10. What is obtained is a mixture of the photon acting as a vector meson and it acting as an anomalous state. This should therefore be the sum of what can be obtained with MSTP(14)=2 and =3. It is distinct from MSTP(14)=1 in that different sets are used for the parton distributions - in MSTP(14)=1 all the contributions to the photon distributions are lumped together, while they are split in VMD and anomalous parts for MSTP(14)=10. Also the beam remnant treatment is different, with a simple Gaussian distribution (at least by default) for MSTP(14)=1 and the VMD part of MSTP(14)=10, but a powerlike distribution d(kT^2)/kT^2 between PARP(15) and Q for the anomalous part of MSTP(14)=10. (PYINIT, PYINPR and PYSTAT changed.) To access this option for e and gamma as incoming beams, it is only necessary to set MSTP(14)=10 and keep MSEL at its default value. Unlike the corresponding option for gamma-p and gamma-gamma, no cuts are overwritten, i.e. it is still the responsability of the user to set these appropriately. Those especially appropriate for DIS usage are CKIN(21)-CKIN(22) or CKIN(23)-CKIN(24) for the x range (former or latter depending on which side is the incoming real photon), and CKIN(35)-CKIN(36) for the Q2 range. A further new option has been added (in PYKLIM): CKIN(39), CKIN(40) : (D=4., -1. GeV^2) the W2 range allowed in DIS processes, i.e. subprocess number 10. If CKIN(40) < 0., the upper limit is inactive. Here W2 is defined in terms of W2 = Q2 * (1-x)/x. This formula is not quite correct, in that (i) it neglects the target mass (for a proton), and (ii) it neglects initial-state photon radiation off the incoming electron. It should be good enough for loose cuts, however. A warning about the usage of PDFLIB for photons. So long as MSTP(14)=1, i.e. the photon is not split up, PDFLIB is accessed by MSTP(56)=2 and MSTP(55) = parton distribution set, as described in the manual. However, when the VMD and anomalous pieces are split, the VMD part is based on a rescaling of pion distributions by VMD factors (except for the SaS sets, that already come with a separate VMD piece). Therefore, to access PDFLIB for MSTP(14)=10, it is not correct to set MSTP(56)=2 and a photon distribution in MSTP(55). Instead, one should put MSTP(56)=2, MSTP(54)=2 and a pion distribution code in MSTP(53), while MSTP(55) has no function. The anomalous part is still based on the SaS parametrization, with PARP(15) as main free parameter. ----- A change has been made in PYREMN to reduce the possibility of infinite loops. ----------- 13, 22 March 1995: The SaS parton distributions of the photons are now available. For details on these sets, see G.A. Schuler and T. Sjostrand, "Low- and high-mass components of the photon distribution functions", CERN-TH/95-62 and LU TP 95-6. There are four new sets. These differ in that two use a Q0=0.6 GeV and two a Q0=2 GeV, and in that two use the DIS and two the MSbar conventions for the dominant non-leading contributions. (However, the fits are formally still leading-order, in that not all next-to-leading contributions have been included.) New default is the SaS 1D set. Furthermore, for the definition of F_2^gamma, additional terms appear that do not form part of the parton distributions itself. To partly take this into account, an additional doubling of the possibilities has been included. These possibilites can be accesed with MSTP(55): MSTP(55) : (D=5) choice of parton-distribution set of the photon; see also MSTP(56). = 1 : Drees-Grassie. = 5 : SaS 1D (in DIS scheme, with Q0=0.6 GeV). = 6 : SaS 1M (in MSbar scheme, with Q0=0.6 GeV). = 7 : SaS 2D (in DIS scheme, with Q0=2 GeV). = 8 : SaS 2M (in MSbar scheme, with Q0=2 GeV). = 9 : SaS 1D (in DIS scheme, with Q0=0.6 GeV). = 10 : SaS 1M (in MSbar scheme, with Q0=0.6 GeV). = 11 : SaS 2D (in DIS scheme, with Q0=2 GeV). = 12 : SaS 2M (in MSbar scheme, with Q0=2 GeV). Note 1 : sets 5 - 8 use the parton distributions of the respective set, and nothing else. These are appropriate for most applications, e.g. jet production in gamma-p and gamma-gamma collisions. Sets 9 - 12 instead are appropriate for gamma*-gamma processes, i.e. DIS scattering on a photon, as measured in F_2^gamma. Here the anomalous contribution for c and b quarks are handled by the Bethe-Heitler formulae, and the direct term is artificially lumped with the anomalous one, so that the event simulation more closely agrees with what will be experimentally observed in these processes. The agreement with the F_2^gamma parametrization is still not perfect, e.g. in the treatment of heavy flavours close to threshold. Note 2 : Sets 5 - 12 contain both VMD pieces and anomalous pieces, separately parametrized. Therefore the respective piece is automatically called, whatever MSTP(14) value is used to select only a part of allowed photon interactions. For other sets (set 1 above or PDFLIB sets), usually there is no corresponding subdivision. Then an option like MSTP(14)=2 (VMD part of photon only) is based on a rescaling of the pion distributions, while MSTP(14)=3 gives the SaS anomalous parametrization. Note 3 : Formally speaking, the k0 (or p0) cut-off in PARP(15) need not be set in any relation to the Q0 cut-off scales used by the various parametrizations. Indeed, due to the familiar scale choice ambiguity problem, there could well be some offset between the two. However, unless you know what you are doing, it is strongly recommended that you let the two agree, i.e. set PARP(15)=0.6 for the SaS 1 sets and =2 for the SaS 2 sets. PARP(15) : (D=0.6 GeV) default value changed for k0 cut-off for separation between direct, VMD and anomalous photons; see Note 3 for MSTP(55) above. The generic routine PYSTFU has been rewritten to handle the interfacing. The old routines PYSTAG, PYSTGS, PYDILN and PYSTHG have been removed. Instead the routines of the SaSgam library have been inserted. In order to avoid any clashes, the routines SAS*** have been renamed PYG***. Thus new routines are PYGGAM, PYGVMD, PYGANO, PYGBEH and PYGDIR. The common block SASCOM is renamed PYINT8. If you want to use the parton distributions for standalone purposes, you are encouraged to use the original SaSgam routines rather than going the way via the Pythia adaptations. COMMON/PYINT8/XPVMD(-6:6),XPANL(-6:6),XPANH(-6:6),XPBEH(-6:6), &XPDIR(-6:6) Purpose: to store the various components of the parton distributions when the PYGGAM routine is called. XPVMD(KFL) : gives distributions of the VMD part (rho, omega and phi). XPANL(KFL) : gives distributions of the anomalous part of light quarks (d, u and s). XPANH(KFL) : gives distributions of the anomalous part of heavy quarks (c and b). XPBEH(KFL) : gives Bethe-Heitler distributions of heavy quarks (c and b). This provides an alternative to XPANH, i.e. both should not be used at the same time. XPDIR(KFL) : gives direct correction to the production of light quarks (d, u and s). This term is nonvanishing only in the MSbar scheme, and is applicable for F_2^gamma rather than for the parton distributions themselves. ----- PYDOCU has been corrected so that PARI(2) refers to the full cross section for gamma-p and gamma-gamma processes, rather than that of the latest subprocess considered. An additional check has been inserted into PYREMN. ----------- 14, 23 March 1995: Some minor modifications to PYSTFU and PYGGAM in the wake of the changes of the previous version. ----------- 15, 24 April 1995: An unfortunate choice of default values has been corrected: the old MSTP(3)=2 value implied that Lambda_QCD was entirely based on the Lambda value of the proton structure function; also e.g. for e+e- annihilation events. Thus the Lambda in PARJ(81) was overwritten, i.e. did not keep the value required by standard phenomenology, which typically gave too narrow jets. (While switching to MSTP(3)=1 it worked fine.) In the modified option MSTP(3)=2 this has been corrected, to better agree with user expectations. Change affects PYINIT and PYRESD. (See further version 16 for additional changes.) MSTP(3) : (D=2) choice of Lambda_QCD values. = 1 : separate for hard scattering, initial showers and final showers, as before. Additionally separate for resonance decays, given in PARP(3). = 2 : most Lambda values are set to that of the proton structure function used, except for the Lambda used in the decay of a resonance (as treated in PYRESD). There the PARP(3) value is used, with default as in e+e-. = 3 : all Lambda values are set to that of the proton structure function used, as was the case for =2 before. PARP(3) : (D=0.29 GeV) the Lambda value used in timelike parton showers in the decay of a resonance (in PYRESD). ----- The form for PTMANO, the pTmin for anomalous processes, as used in PYINPR when processes are mixed for gamma-p or gamma-gamma events, has been updated to match (as well as can be expected) the SaS 1D photon distributions. ----------- 16, 30 June 1995: The strategy for the changes to MSTP(3) in version 15 above have been modified for better transparency. The parameter PARP(3) has been removed, and instead PARP(72) has been introduced. Now PARJ(81) is used for resonance decays (including e.g. Z0 decay, from which it is determined), and PARP(72) for other timelike showers. PARJ(81) is not overwritten for MSTP(3) = 2, but only for = 3. Changes affect PYINIT, PYEVNT and PYRESD. PARP(72) : (D=0.25 GeV) the Lambda value used in timelike parton showers except in the decay of a resonance. ----- A new multiplicative factor has been introduced for the hard scattering in PYSIGH. PARP(34) : (D=1.) the Q2 scale defined by MSTP(32) is multiplied by PARP(34) when it is used as argument for structure functions and alpha_s at the hard interaction. It does not affect alpha_s when MSTP(33)=3, nor does it change the Q2 argument of parton showers. ----- PYREMN has been corrected for occasional too large boost factors. An error in PYSIGH for process 148 has been corrected. The MSTP(62)=1 option of PYSSPA is modified to avoid division by zero. Header has been updated with WWW-information. ----------- 17, 23 August 1995: MIN1, MIN2, MAX1, MAX2, MINA and MAXA in PYSIGH have had an extra M prefixed to avoid confusion with Fortran functions. Protect against MDCY(0,1) being accessed in PYSIGH. Protect against THB=0 in PYRAND. Protect against YSTMAX-YSTMIN = 0 in PYSIGH. Check for moved leptoquark at beginning of PYRESD just like for other particles with colour. ----------- 18, 14 September 1995: The protection above against YSTMAX-YSTMIN = 0 in PYSIGH turned out to be flawed: when used with electron-inside-electron structure functions it cut out part of the allowed phase space and thus gave too low a cross section for several e+e- processes. This is now corrected. -----------------------------------------------------------------