12 February 1993 Updates to PYTHIA 5.6 and JETSET 7.3 Physics and Manual Torbjorn Sjostrand Theory Division, CERN CH-1211 Geneva 23 Switzerland -------------------------------------------------- 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.6488/92 described the status as of 20 May 1992. The LaTeX file with this manual will be updated, though less frequently than the programs. (One revised version has appeared, which documents changes made up until 6 September.) 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 after 20 May 1992. 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.3 16 June: LUHEPC modified to avoid problems with event history when pileup events are generated in PYTHIA. 17 June: add option 16 to LUEDIT. With a CALL LUEDIT(16), some missing daughter pointers in K(I,4) amd K(I,5) will be reconstructed from the mother pointers in K(I,3). 17 June: when 3-jet events are generated according to the scalar gluon option of LUEEVT, the full angular distribution (including weak effects) is included according to the formulae in Ref. Lae80. 26 June: a techni-eta is introduced as particle 38, with default mass PMAS(38,1) equal to 350 GeV and the three (dominant) decay modes b + bbar, t + tbar and g + g. So as not to mess up the ordering of existing decay channels, these modes are put in lines 1201 - 1203. The techni-eta has zero spin, is neutral, singlet under weak SU(2), but a colour octet. 10 July: introduce code 110 for a diffractive state of a pi0, a rho0, or a gamma. 29 July: in LUSHOW do not set history pointers in K(I,4) and K(I,5) for photon. 31 July: skip calculation of momentum shift in LUBOEI for pairs already so close that numerical precision problems arise. 2 September: modify expression for maximum weight in tau decay in LUDECY to make generation more efficient. 17 October: introduce extra parameter to determine action of LUCELL calls. Cells with transverse energy below a threshold PARU(58) are not considered, i.e. the energy in those cells is set vanishing. By default PARU(58)=0., i.e. the old behaviour is retained. Also correct mistake in LUCLUS, which gave a wrong value for the minimum distance between two clusters as stored in PARU(63) (does not affect the algorithm as such, only this complementary piece of information). 28 October: introduce possibility for b quark to develop parton shower in a top hadron decay, i.e. T -> W + b + spectator. Is intended for scenarios where the top quark first fragments and only later decays; cf. PYTHIA note of 14 October. This possibility is regulated by switch MSTJ(27), with default value 0 (for backwards compatibility). The possibilities are: MSTJ(27) = 0 : (default) no parton shower, i.e. old behaviour. = 1 : b parton shower; the W is thereafter assumed to decay isotropically. = 2 : b parton shower; the W thereafter decays anisotropically, as in option 0, but with W momentum reduced when the b parton shower initiator is off the mass shell. Differences between options 1 and 2 are minor; in principle 2 is the better but 1 is better matched to what is done when one uses option MSTP(48)=1 and lets top decay before fragmentation. 29 October: correct a minor bug in LUZDIS (alternative fragmentation function f(z) = (1-z)**c could not be reached). 3 November: change the default values for heavy flavour fragmentation functions to correspond to the Peterson parametrizations with values epsilon_c = 0.06, epsilon_b = 0.006 and epsilon_t = 0.00001. This corresponds to a change of PARJ(54) = -0.06, PARJ(55) = -0.006, and PARJ(56) - PARJ(58) = -0.00001. To enable these values, it is recommended to use MSTJ(11)=3. -------------------------------------------------- 3) Changes in PYTHIA 5.6 2 June: modification of the scheme for preservation of (x, Q2) in leptoproduction. As before, this is achieved by a post facto modification of kinematics. The main problem is that modification may often be impossible, due to kinematical constraints. Earlier such events were simply thrown, which meant some bias in the generated distributions. In the new scheme, the modification is done somewhat differently, so that the failure rate has gone down significantly. In addition, in case of failure, up to five new tries are made from the same hard scattering kinematics, i.e. nominal (x, Q2) values, but with new initial and final state shower evolution. The preservation option should therefore be more useful. 22 June: introduce a call to LUEDIT(16) in PYEVNT so as to reconstruct a few lost parts of event history (see JETSET, 17 June). 22 June: introduce parametrizations of the gamma-p total cross-section in PYXTOT. The variable MSTP(30) can be used to select between the two alternatives, as follows. MSTP(30) = 1 : (default) the fit by Landshoff, see LP-HEP proceedings. = 2 : the CERN-HERA fit, see Review of Particle Properties. See further 8 July. 23 June: correct trivial bug which gave too high rate for the decay of supersymmetric Higgs particles to tau+ tau-. 24 June: introduce (simple) description of multiple interactions in gamma-p events, using the same formalism as for p-p events. Thus you can use MSEL=1 to generate a combination of high-pT and low-pT events, where a high-pT event may actually contain several parton-parton interactions. The switches MSTP(81) and so on work exactly as in p-p events. The main limitation is that this only works for a photon beam, i.e. the formalism can not be used for photoproduction events in e-p scattering. (The main reason is that the multiple interactions formalism currently needs to be initialized at a fixed c.m. energy of the hadronic system, while the e-p event generation gives a spectrum of hadronic c.m. energies.) 29 June: implement process 149, g + g -> eta_techni (see also JETSET, 26 June). Three decay modes are included, b + bbar, t + tbar and g + g. If kinematically allowed, the t + tbar decay mode always dominates. The two main free parameters are the techni-eta mass and the decay constant F_pi. The latter appears inversely quadratically in all the partial widths. Since the total cross-section is proportional to the g + g partial width, F_pi also affects the size of the cross-section. F_pi is stored in PARP(46), and has the default value 123 GeV. Further information about the model implemented can be found in T. Appelquist and G. Triantaphyllou, Yale University preprint, in preparation. 3 July: expand process 77 so it also works for like sign longitudinal W scattering, W+ W+ -> W+ W+ and W- W- -> W- W-. 6 July: introduce further option MSTP(22)=4 for parton shower evolution scale in deep inelastic scattering chosen to be Q2 * (1-x) * max(1, ln(1/x)) (see Ingelman, HERA workshop). 6 July: modify simulation of charged Higgs to include the down-type part of quark couplings in process 161, and to include the running of down-type mass in processes 143 and 161 (if requested by MSTP(37)). 7 July: correct an error in the MSTP(13)=2 (non-default) option of PYSTEL, whereby the structure functions of partons inside a photon inside an electron were overestimated by about 10% - 20%. 7 July: add an option 3 to a PYFRAM call, whereby you can boost a leptoproduction event to the hadronic c.m. frame. Mainly intended for deep inelastic scattering, but can also be used in photoproduction. Note that the lepton, and any photons radiated by it, are not removed from the event record; if you only want to study the hadrons, you have to remove it yourself. 8 July: add options MSTP(30)=0 for gamma-p events and MSTP(31)=0 for hadron-hadron ones, wherein you may put your own values for cross-sections, according to the following: PARP(21) - PARP(26) : (D=6*0.) cross-sections and slope parameters that you can set yourself in the options MSTP(30)=0 and MSTP(31)=0, respectively. If these options are used, at least PARP(21) and PARP(25) must be assigned nonvanishing values, or the program is likely to crash. All values are taken to refer to the c.m. energy of the current run, i.e. the energy dependence of parameters need not be specified (unlike the generic parametrizations found in the programs). PARP(21) : total cross-section sigma_tot in mb. PARP(22) : elastic cross-section sigma_el in mb. PARP(23) : single diffractive cross-section sigma_sd in mb. PARP(24) : double diffractive cross-section sigma_dd in mb. PARP(25) : nuclear slope parameter B in GeV^(-2). PARP(26) : curvature parameter C in GeV^(-4). Note added January 1993: when the VDM description of gamma-p events was expanded on 23 September, this section was not expanded accordingly. In particular, with the new description it is also necessary to specify how cross-sections are distributed among rho0, omega and phi VDM states. After the PYINIT call, the cross-sections in VINT(101) - VINT(107) should therefore be split among the new variables VINT(241) - VINT(267) in the relative rho0/omega/phi proportions. In other words, one must select VINT(241) + VINT(251) + VINT(261) = VINT(101), etc. Further, one should put VINT(248) = VINT(258) = VINT(268) = VINT(121), and MINT(101) = 3 if the beam partile (in the PYINIT call) is a photon, while MINT(102) = 3 if the target particle is a photon. 8 July: include elastic and diffractive cross-sections for gamma-p interactions. You can therefore use MSEL=2, as for p-p, to get a mixture of elastic, single diffractive, double diffractive and ordinary non-diffractive events. It is not possible, in the same run, to generate also the direct contributions, where the photon interacts pointlike (see MSTP(14)). A summary of the current status: For gamma-p events, where the photon acts like a hadron in the vector dominance model sense, the default parametrization of the total cross-section is taken from Landshoff. The cross-section for the elastic process gamma p to rho^0 p has been parametrized as sigma_el(s) = sigma_tot(s) * 0.058 * (s / 1 GeV^2)^0.08, based on the numbers predicted by Schuler and Terron. For lack of any other information, the single diffractive cross-section has been assumed equal to the elastic one, and the double diffractive to a third of this. The latter follows from the relation suggested by Goulianos sigma_dd(s) = sigma_sd^2(s) / (3 * sigma_el(s)) . In elastic scattering, the t distribution is always parametrized as d(sigma)/dt = const. * exp(B * t + C * t^2). where B is the slope parameter and C the curvature. For gamma-p, one assumes B = 9.2 GeV^(-2) + 0.5 GeV^(-2) * ln(0.25 GeV^(-2) * s) , which is a simplified version of the form given by Schuler and Terron, while the curvature parameter C has been assumed vanishing. 10 July: replace the existing elastic and diffractive kinematics generation machinery by a new one, as follows: For elastic scattering, t is generated according to the formula given above. If the curvature C is vanishing or negative, then t is allowed to go up to the kinematical limit. For a positive C, i.e. for a distribution which eventually diverges, the distribution is cut off at t = -B/(2C). For diffractive scattering, the distribution in t is the square-root of the distribution for elastic events, i.e. B and C are rescaled by a factor 1/2. For each diffractive system an additional factor dM^2/M^2 is multiplied on to give the mass M of the excited state. Initially the two or three variables t, M_1^2 and M_2^2 are selected independently of each other, within the full allowed range. However, a multiplet of variables which is internally kinematically inconsistent is then rejected and a new one generated in its place. This leads to a slight shift of the average abs(t) value upwards and a non-negligible suppression of the high-mass tails of the diffractive systems. A rho^ formed by gamma to rho^0 in elastic or diffractive scattering is normally longitudinally polarized, and therefore its decay angular distribution in rho^0 to pi^+ pi^- is taken to be proportional to sin^2(theta), where the reference axis is given by the rho^0 direction of motion. The mass spectrum of the diffractive system is assumed to start 0.2 GeV above the mass of the incoming particle, or 0.2 GeV above the rho^ mass for an incoming gamma. This value is stored in PARP(102) and could be changed, but must be at least a pi^0 mass (plus some saftey margin). A light diffractive system, with a mass less than 1 GeV (PARP(103)) above the mass of the incoming particle, is taken to decay isotropically into a two-body state. Single-resonance diffractive states, such as a Delta^+, are therefore not explicily generated, but are assumed described in an average, smeared-out sense. A more massive diffractive system is treated as before, i.e. the 'Pomeron' can either interact with a quark or with a gluon, and thus pull out either a single or a double string (see MSTP(101)). 6 August: modify PYREMN so that check on remnant mass is not made when there is no remnant. 7 August: modify PYREMN so that energy sharing in resolved photon remnant is as for meson; see MSTP(92) etc. 8 August: modify PYSTFU and PYSTEL so that the PDFLIB photon structure functions may be used. To use this feature, put MSTP(56)=2 and select PDFLIB set in MSTP(55). As before, you need to enable lines beginning with C! to allow the proper calls to be made. 12 August: replace PDFLIB calls to STRUCTF by calls to STRUCTM, to allow for a difference between up and down sea quarks. Also allow option MSTP(57)=0 (structure functions evaluated at lower cutoff scale) for PDFLIB. 31 August: modify PYKMAP to protect from numerical precision problem in 2 -> 3 processes. 23 September: introduce new parametrizations for total and partial cross-sections, the former according to Donnachie-Landshoff (CERN-TH.6635/92), the latter according to Schuler-Sjostrand (paper in preparation). The partial cross-sections include elastic, single diffractive on either of the two sides, and double diffractive. The new parametrizations are automatically default for gamma-p events, MSTP(30)=1, while the old parametrizations are retained for the alternative, MSTP(30)=2. For backwards compatibility the default is unchanged for hadron-hadron collisions, and the new option is added as MSTP(31)=6. In principle you are recommended to use this option, since e.g. it reproduces the Tevatron total and elastic cross-sections much better. You should still be warned that some parameters, e.g. for multiple interactions, may need to be retuned somewhat for this alternative. In addition, the treatment of diffractive system is further improved, in two respects. (1) The slope parameter B is evaluated in Pomeron language, which gives predictions that depend on the nature of the incoming particles, on the c.m. energy, and on the mass(es) of the diffractive system(s), separately for elastic, single and double diffractive events. (2) In photoproduction, the photon may be represented either by a rho0, an omega or a phi, in the appropriate proportions. The couplings are inversely proportional to the expressions f_V^2/(4pi), which are stored in PARP(161) - PARP(163) (default values 2.20, 23.6, 18.4) for rho0, omega and phi, respectively. Further, the default value of PARP(102) is changed from 0.2 to 0.28 GeV; this represents the minimum mass allowed for a diffractive system above the mass of the incoming particle. For internal use, a number of MINT and VINT variables are now used, including MINT(92), MINT(101) - MINT(107), VINT(107) and VINT(241) - VINT(271). 25 September: introduce two new options MSTP(14) = 2 and 3 for the structure of a photon. The options now include: MSTP(14) (D=1) : nature of photon beam. = 0 : direct, i.e. without internal structure. = 1 : resolved, i.e. all interactions where a photon structure function is needed. = 2 : resolved, VDM only, i.e. where the photon structure function can be taken to be an appropriately scaled-down version of the rho0 (= pi0) one. = 3 : resolved, anomalous only, i.e. the additional perturbatively calculable part of the photon structure function initiated by branchings gamma -> q + qbar. In principle, options 2 and 3 together correspond to option 1, modulo the fact that existing structure function parametrizations do not give a 1-to-1 correspondence. However, the separate treatments in 2 and 3 allows a more precise modelling e.g. of beam jet structure, than can be obtained with the more rough option 1. This will be described in a paper by Schuler and Sjostrand, in preparation. However, it should be remembered that large uncertainties remain. A full description of photon interactions therefore require separate runs, either two with MSTP(14) = 0 and 1, or three with MSTP(14) = 0, 2 and 3. Kindly note that the current description only allows for one incoming photon, i.e. does not yet cover the gamma-gamma case, where e.g. one photon is direct and another is resolved (according to either of the options). The structure function routine PYSTFU is modified to call the pi+ structure function for MSTP(14)=2, and turn it into a pi0 one which then is scaled down. For MSTP(14)=3 a new routine PYSTAG is called by PYSTFU. This routine (plus help routines PYSTGS and PYDILN), adapted from G. Schuler, contains the anomalous part of the photon structure function. A nominal lower cutoff scale P0 is needed, at which the anomalous structure function is assumed vanishing. This cutoff P0 is stored in PARP(15), and has the default value 0.5 GeV. For charm and bottom flavours the threshold is assumed to appear at Q^2 = PARP(16) * m_q^2. By default PARP(16) = 1. The actual scale P^2, at which the photon is resolved, is logarithmically distributed between the lower cutoff P0^2 and the scale Q^2 of the hard interaction. The current P^2 value is stored in VINT(272). In the initial state radiation machinery, an option MSTP(66) = 1 is included by default, wherein the P^2 scale is used as an additional lower cutoff scale for QCD radiation. The option MSTP(66) = 0 will always do backwards evolution down to the standard cutoff PARP(62). 1 October: introduce routine PYSTFL, which allows a modification of the low-Q^2 and low-x behaviour of the proton structure function. This option is accessed with the switch value MSTP(57)=2. It can be used in conjunction with any internal or PDFLIB proton structure function. Its main effect is to make structure functions vanish in the limit Q^2 -> 0. 5 October: introduce option to dampen the cross-section of an arbitrary hard 2 -> 2 process in the limit pT -> 0, in the same way as is done with QCD processes in options MSTP(82) >= 2 of the multiple interactions scenario. More precisely, the cross-section is multiplied by pT^4/(pT0^2 +pT^2)^2, where pT0 = PARP(82). In addition alpha_S is evaluated at the scale pT0^2 + pT^2 rather than pT^2. This possibility is switched on with MSTP(85) = 1; default is MSTP(85) = 0, i.e. off. 14 October: introduce option to allow top quarks to decay before they fragment. To use this option, put MSTP(48)=1. Additionally you have to change to MDCY(6,1)=1. Default is MSTP(48)=0, MDCY(6,1)=0, i.e. top quarks fragment into top hadrons. Currently only the decay mode t -> b + W is implemented. Note that, in the new mode, the top decays to a W with code +-24, while for top hadron decay the separate code 89 (Wvirt) has been used. This is of importance e.g. if one wishes to select only semileptonic decays. The treatment of top is still not that of the standard resonances, in that cross-section for top production is not reduced when the subsequent W decay is limited to a few modes; you have to multiply with the respective branching rations yourself. Also kindly note that spin correlations are not included, i.e. the t and the W both decay isotropically in their respective rest frames. 28 October: As consequence of JETSET change one can now also include parton shower evolution in PYTHIA. The following summarizes the recent top changes: If a top is lighter than about 120 GeV, then the top quark is sufficiently long-lived that top hadrons T are formed, which subsequently decay. For heavier tops, the decay is faster, and therefore hadrons have no time to be formed. Of course there is a transition region where the issue is not so clear - where things will work out differently in each event. So, in principle, one should try both scenarios to check there are no serious disagreement between them. So far as I can tell, they indeed seem to match well. Earlier, only the first alternative was available. Now you can use MSTP(48) to switch between the two. MSTP(48) : (D=0) flag to be set. = 0: as before; top quark fragments to top hadrons which then decay. Recommended for top mass below 120 GeV. = 1: new; top decays to W + b, and the b subsequently fragments. Recommended for top mass above 120 GeV. You must also change MDCY(6,1)=1 to use this option, see below. Depending on which alternative you choose above, there are other switches that can be set. The first deals with QCD radiation off the b quark. This was negligible when the top was assumed to weigh below 100 GeV, but is more significant for the heavier tops now used. The main effect of such gluon emission is to soften the momentum of the B quarks produced, and add a few extra soft pions. In the option MSTP(48)=0, you can switch on this emission with MSTJ(27) = 1 or 2. The default is MSTJ(27)=0, i.e. no emission, for reasons of backwards compatibility, so you are recommended to change. When MSTP(48) = 1, QCD radiation is automatically included as part of standard, and need not be selected. (It can be switched off, together with other final state radiation, by putting MSTP(71) = 0.) Some theoretical ambiguity exists in the modelling of this gluon emission, but what is offered should be a reasonable first approximation. The t -> b + W is followed by the W decay. In principle this decay carries some spin information, i.e. is not quite isotropic. This information is not available in the MSTP(48) = 1 option i.e. here the W decays isotropically. For MSTP(48) = 0 it is made use of when MSTJ(27) = 2, while MSTJ(27) = 1 implies isotropic W decay. In practice, effects are small, to the best of my knowledge. For MSTP(48) = 0 top hadrons are produced and subsequently decay. Therefore the decay channels of the top are set for the top hadron pseudoparticle 86. For MSTP(48) = 1 no hadrons are formed, and it is instead the t quark itself, with code 6, that decays. Since this was not originally foreseen, you must change to MDCY(6,1)=1. In scheme MSTP(48) = 0 the decay W is in fact the Wvirt with code 89, while the standard W with code 24 is used for MSTP(48) = 1. The separate code 89 is related to the necessity to be more sophisticated if the t mass is below or around the W one, and is therefore now more of a historical curiosity. It implies, however, that you have to select the W decay channels differently in the two scenarios. If you print cross-sections with PYSTAT(1), then remember that these do not include the top and W branching ratios, for the case that you only allow some channels. You therefore have to multiply with these yourself. This holds true in either option. So far, I have not implemented the potential t -> b + H mode. In principle, this is not difficult, but did not get top priority. It is recommended to switch to Peterson fragmentation for the heavy quarks and keep the standard Lund one for the other hadrons. This you do by putting MSTJ(11) = 3. I have included as defaults the recommended epsilon values epsilon_c = 0.06, epsilon_b = 0.006 and epsilon_t = 0.00001 in the new Jetset version, so you need not set these yourself. The c and b values are typical LEP fit results, with QCD emission included. Note that there is some mismatch of Q2 scales between e+e- and pp production of heavy flavours, so the epsilon values above are not holy numbers, but should be a good first guess. 29 October: trace how the position of a t is moved by shower evolution, so that option MSTP(48)=1 takes effect also e.g. in Z0 -> t + tbar. 14 December: fix numerical precision problem in PYDILN. 13 January: correct mass for gamma + p -> phi + p in PYXTOT. Constrain acos evaluation in PYI3AU to avoid math library bug on SUN. 29 January: possibility to run PYTHIA with varying energies from one event to the next. Here comes a summary of the new features: It is now possible to use PYTHIA in a mode where the energy can be varied from one event to the next, without the need to reinitialize with a new PYINIT call. This allows a significant speed-up of execution, although it is not as fast as running at a fixed energy. It can not be used for everything - we will come to the fine print at the end - but should be applicable for most tasks. The master switch to access this possibility is in MSTP(171). By default it is off, so you must set MSTP(171)=1 before initialization. There are two submodes of running, MSTP(172)= 1 or 2. In the former mode, PYTHIA will generate an event at the requested energy. This means that you have to know which energy you want beforehand. In the latter mode, PYTHIA will often return without having generated an event - with flag MSTI(61)=1 to signal that - and you are then requested to give a new energy. The energy spectrum of accepted events will then, in the end, be your naive input spectrum weighted with the cross-section of the processes you study. We will come back to this. The energy can be varied, whichever frame is given in the PYINIT call. In addition two new frames are allowed, to increase the flexibility. The summary is: CALL PYINIT(FRAME,BEAM,TARGET,WIN) FRAME : frame of experiment. = 'CMS' : WIN gives nominal c.m. energy used at initialization. For each event to be generated thereafter, PARP(171) should be filled by the fractional energy of that one, i.e. E_cm = PARP(171)*WIN. Therefore PARP(171) should normally be smaller than unity. = 'FIXT' : WIN gives initialization beam energy, with target at rest. For each event PARP(171) should be filled by the fractional beam energy of that one, i.e. E_beam = PARP(171)*WIN. Therefore PARP(171) should normally be smaller than unity. = 'USER' : give in beam and target three-momenta in P(1,1) - P(1,3) and P(2,1) - P(2,3), respectively. These variables have to be given at initialization and for each new event. The particle masses are assumed to be the nominal ones and the particle energies are calculated by the program. The c.m. energy for events to be generated should normally be smaller than the one used at initialization. = 'FOUR' : as 'USER', except that you should give P(1,1) - P(1,4) and P(2,1) - P(2,4), i.e. also the energies of the two incoming particles. These particles need not be on the mass shell. = 'FIVE' : as 'USER', except that you should give P(1,1) - P(1,5) and P(2,1) - P(2,5), i.e. also the energies and masses of the two incoming particles. These particles need not be on the mass shell. For a particle with m**2 = E**2 - p**2 < 0, one should give P(I,5) = -sqrt(-m**2). Since the mass should be uniquely given by the energy and momentum, the only advantage of this option compared with 'FOUR' is that a small mass may be better known by the user than the program can calculate itself, due to the limited precision. To illustrate the use of the MSTP(172)=2 facility, consider the case of beamstrahlung in e+e- linear colliders. This is just for convenience; what is said here can be translated easily into other situations. Assume that the beam spectrum is given by D(z), where z is the fraction retained by the original e after beamstrahlung. Therefore 0 <= z <=1 and the integral of D(z) is unity. This is not perfectly general; one could imagine branchings e- -> e- + gamma -> e- + e+ + e-, which gives a multiplication in the number of beam particles. This could either be expressed in terms of a D(z) with integral larger than unity or in terms of an increased luminosity. We will assume the latter, and use D(z) properly normalized. Given a nominal s = 4*E_beam**2, the actual s' after beamstrahlung is given by s' = z_1 * z_2 * s. For a process with a cross-section sigma(s) the total cross-section is sigma_tot = Int_0^1 Int_0^1 D(z_1) D(z_2) * * sigma(s' = z_1 * z_2 * s) dz_1 dz_2. The cross-section sigma may in itself be an integral over a number of additional phase space variables. If the maximum of the differential cross-section is known, a correct procedure to generate events is 1) pick z_1 and z_2 according to D(z_1) and D(z_2); 2) pick a set of phase space variables of the process, for the given s' of the event; 3) evaluate sigma(s') and compare with sigma_max; 4) if event is rejected, then return to 1) to generate new variables; 5) else continue generation to give complete event. You as a user are assumed to take care of step 1), and present the resulting kinematics with incoming e+ and e- of varying energy. Thereafter PYTHIA will do steps 2 - 5, and either return an event or put MSTI(61)=1 to signal failure. The maximization procedure of PYTHIA does search in phase space to find sigma_max, but it does not vary the s' energy in this process. Therefore the maximum search in the PYINIT call should be performed where the cross-section is largest. For processes with increasing cross-section as function of energy this means at the largest energy that will ever be encountered, i.e. s' = s in the case above. This is the 'standard' case, but often one encounters other behaviours, where more complicated procedures are needed. One such case would be the process e+ e- -> Z* -> Z + H, which is known to have a cross-section which increases near the threshold but is decreasing asymptotically. If one already knows that the maximum, for a given Higgs mass, appears at 300 GeV, say, then the PYINIT call should be made with that energy, even if subsequently one will be generating events for a 500 GeV collider. In general, it may be necessary to modify the selection of z_1 and z_2 and assign a compensating event weight. For instance, consider a process with a 1/s cross-section. Then the sigma_tot expression above may be rewritten as sigma_tot = Int_0^1 Int_0^1 D(z_1)/z_1 D(z_2)/z_2 * * z_1 z_2 sigma(s' = z_1 * z_2 * s) dz_1 dz_2. The expression z_1 z_2 sigma(s') is now flat in s', i.e. not only can sigma_max be found at a convenient energy such as the maximum one, but additionally the PYTHIA generation efficiency (the likelihood of surviving step 4)) is greatly enhanced. The price to be paid is that z has to be selected according to D(z)/z rather than according to D(z). Note that D(z)/z is not normalized to unity. One therefore needs to define I_D = Int_0^1 D(z)/z and a properly normalized D'(z) = (1/I_D) * D(z)/z. Then sigma_tot = Int_0^1 Int_0^1 D'(z_1) D'(z_2) * * (I_D)**2 z_1 z_2 sigma(s' = z_1 * z_2 * s) dz_1 dz_2. Therefore the proper event weight is WT = (I_D)**2 z_1 z_2. This weight should be stored, for each event, in PARP(173). Additionally you must put MSTP(173)=1 for the program to make use of weights at all. Often D(z) are not known analytically; therefore I_D is also not known beforehand, but may have to be evaluated (by you) during the course of the run. Then you should just use WT = z_1 z_2 in PARP(173) and do the overall normalization yourself in the end. Only the cross-sections are affected by this, the events themselves still are properly distributed in s' and internal phase space. Above it has been assumed tacitly that D(z) -> 0 for z -> 0. If not, D(z)/z is divergent, and it is not possible to define a properly normalized D'(z) = D(z)/z. If the cross-section is truly diverging like 1/s, then a D(z) which is nonvanishing for z -> 0 does imply an infinite total cross-section, whichever way things are considered. In cases like that, it may be necessary to impose a lower cut on z, based on some physics or detector consideration. The most difficult cases are those with a very narrow and high peak, such as the Z0. One could initialize at the energy of maximum cross-section and use D(z) as is, but efficiency might turn out to be very low. One might then be tempted to do more complicated transforms of the kind illustrated above. As a rule it is then convenient to work in the variables tau_z = z_1 * z_2 and y_z = (1/2) ln (z_1/z_2), cf. section 7.2 of the manual. Clearly, the better the behaviour of the cross-section can be modelled in the choice of z_1 and z_2, the better the overall event generation efficieny. Even under the best of circumstances, the efficiency will still be lower than for runs with fix energy. There is also a non-negligible time overhead for using variable energies in the first place, from kinematics reconstruction and (in part) from the phase space selection. One should therefore not use variable energies when not needed, and not use a large range of energies sqrt(s') if in the end only a smaller range is of experimental interest. The following new common block variables therefore exist: MSTP(171) (D=0) possibility of variable energies from one event to the next. = 0 : no; i.e. the energy is fixed at the initialization call. = 1 : yes; i.e. a new energy has to be given for each new event. MSTP(172) (D=2) options for variable energy generation, applicable when MSTP(171)=1. = 1 : an event is generated at the requested energy, i.e. internally a loop is performed over possible event configurations until one is accepted. If the requested c.m. energy of an event is below PARP(2) the run is aborted. Cross-section information can not be trusted with this option, since it depends on how you decided to pick the requested energies. = 2 : only one event configuration is tried. If that is accepted, the event is generated in full. If not, no event is generated, and the status code MSTI(61)=1 is returned. You are then expected to give a new energy, looping until an acceptable event is found. No event is generated if the requested c.m. energy is below PARP(2), instead status code MSTI(61) is set =1 to signal the failure. In principle, cross-sections should come out correctly with this option. MSTP(173) (D=0) possibility for user to give in event weight to compensate for a biased choice of beam spectrum. = 0 : no, i.e. event weight is unity. = 1 : yes; weight to be given for each event in PARP(173). PARP(171) : to be set, event-by-event, when variable energies are allowed, i.e. MSTP(171)=1. If PYINIT is called with FRAME='CMS' (='FIXT'), PARP(171) multiplies the c.m. energy (beam energy) used at initialization. For the options 'USER', 'FOUR' and 'FIVE', PARP(171) is dummy. PARP(173) : event weight to be given by user when MSTP(173)=1. MSTI(61) : status flag set when events are generated. It is only of interest for runs with variable energies, MSTP(171)=1, with the option MSTP(172)=2. = 0 : an event has been generated. = 1 : no event was generated, either because the c.m. energy was too low or because the Monte Carlo phase space point selection machinery rejected the trial point. A new energy is to be picked by the user. Also a few additional internal variables have been introduced: MINT(111) : the frame given in PYINIT call, 0 - 5 for 'NONE', 'CMS', 'FIXT', 'USER', 'FOUR' and 'FIVE', respectively. VINT(289) : squared c.m. energy found in PYINIT call. VINT(290) : the WIN argument of a PYINIT call. VINT(291) - VINT(300) : the two five-vectors of the two incoming particles, as reconstructed in PYINKI or PYEVKI. This facility may be combined with most other aspects of the program. For instance, it is possible to simulate beamstrahlung as above and still include bremsstrahlung with MSTP(11)=1. Further, one may multiply the overall event weight of PARP(173) with a kinematics-dependent weight given by PYEVWT, although it is not recommended (since the chances of making a mistake are also multiplied). However, a few things do NOT work. * It is not possible to use pile-up events, i.e. you must have MSTP(131)=0. * The possibility of giving in your own cross-section optimization coefficients, option MSTP(121) = 2, would require more input than enumerated in the manual, and this option should therefore not be used. You can still use MSTP(121) = 1, however. * The multiple interactions scenario with MSTP(82) >= 2 only works approximately for energies different than the initialization one. If the c.m. energy spread is smaller than a factor 2, say, the approximation should be reasonable, but if the spread is larger one may have to subdivide into subruns of different energy bins. The initialization should be made at the largest energy to be encountered - whenever multiple interactions are possible (i.e. for incoming hadrons and resolved photons) this is where cross-sections are largest anyway, and so is no further constraint. There is no simple possibility to change PARP(82) during the course of the run, i.e. a energy-independent pT0 must be assumed. The default option MSTP(82) = 1 works fine, i.e. does not suffer from the constraints above. If so desired pTmin=PARP(81) can be set differently for each event, as a function of c.m. energy. 11 February: inclusion of J/Psi in the VDM component of photoproduction. It appears on equal footing with rho0, omega and phi when a gamma-p event is generated with MSTP(14)=2. The coupling f_J/Psi**2/4pi is stored in PARP(164). The following points should be noted. * The parametrization of the J/Psi + p total cross section is uncertain. We use the phi-p cross-section as a reference and scale it down by about a factor 10, in agreement with current estimates, but there is an uncertainty by maybe a factor of 2 in this number. What has been measured experimentally (at low energies) is only the elastic gamma + p -> J/Psi + p cross-section, which is proportional to the sigma_tot**2/f_J/Psi**2. In other words, only a combination of sigma_tot and f_J/Psi is measured, which means that a change of sigma_tot is allowed if f_J/Psi is rescaled. * At HERA energies the VDM cross section goes up by about 1.4 mb by the inclusion of this new component. This is so little that it does not affect any of the discussion of photoproduction properties in general, but it gives a non-negligible contribution to charm production in general. If you only want to simulate this component, i.e. switch off the rho0, omega and phi contributions, you can do this by setting PARP(161) - PARP(163) to be very large, e.g. 1E10. * The fraction of elastic J/Psi production is much lower than is the case for rho0. This is expected, based on the smallness of the total cross section (noted above) and the optical theorem. A J/Psi obviously may appear in the final state of elastic and single diffractive events (where the proton becomes a diffractive system). A diffractive J/Psi state normally would produce e.g. D + Dbar mesons, but if the mass of the diffractive state is below this threshold, the program will allow diffractive J/Psi + pi0 states. Therefore also double diffractive events may contain J/Psi. * For the nondiffractive events, one would need knowledge about the J/Psi structure functions to do a good job. For lack of anything else, I assume the same valence spectrum as for pions, which certainly means the c quarks are too soft. Therefore one should be careful not to trust the results too much, in particular not at high transverse momenta. Most are produced in low-pT events, where structure functions are not explicitly used, and where the uncertainty is smaller. * The J/Psi decays isotropically, i.e. no attempt has been made to include the possibility of j/Psi polarization. The manual should be corrected as follows. PARP(164) : (D=11.5) the coupling f_J/Psi**2/4pi used in the VDM model when MSTP(14)=2. VINT(271) - VINT(280) : cross sections for the J/Psi VDM component, parallelling the current description of VINT(241) - VINT(270). VINT(281) : what used to be VINT(271), i.e. ratio between total gamma+X and total pi0+X cross section. VINT(282) : what used to be VINT(272), i.e. scale P**2 at which an anomalous photon is resolved.