CMake: removing qpythia from the depedencies
[u/mrichter/AliRoot.git] / PYTHIA6 / pythia_doc.update
1                   ****************************
2                   *                          *
3                   *    PYTHIA version 6.1    *
4                   *                          *
5                   ****************************
6
7                   (Last updated 17 August 2000)
8
9 The new PYTHIA version is a logical continuation of previous versions.
10 Therefore a user should not feel lost. However, many details have been 
11 changed. The major changes (so far) are:
12 - The supersymmetric process machinery of SPYTHIA has been included.
13 - PYTHIA and JETSET have been merged.
14 - All real variables are declared in double precision.
15 - The internal mapping of particle codes has changed.
16 - Extended capabilities to handle reactions of virtual photons.
17 - Baryon production according to advanced popcorn scheme (new option as 
18   of version 6.110, with some consequences also for default behaviour).
19
20 Below follows a more extensive list of main changes, performed to move 
21 from Pythia 5.7 and Jetset 7.4 to Pythia 6.1. Eventually this file
22 will be complemented by a completely updated manual. However, based
23 on the information here and some common sense it should be possible
24 to use the program already now, if you are familiar with previous 
25 versions. 
26
27 -----------------------------------------------------------------------
28
29 PYTHIA/JETSET CODE MERGING
30
31 * The PYTHIA and JETSET routines have been joined into one single file.
32   - LUDATA and PYDATA have been joined to a single BLOCK DATA.
33   - LUTEST and PYTEST have been joined to a single test program.
34   - SAVE statements are common for former JETSET and PYTHIA routines.
35   - version information (for the title page) is based on MSTP(181-185) 
36     while MSTP(186) and MSTU(181-186) are not used any longer. 
37
38 * All JETSET routines and commonblocks have been renamed to begin 
39   with PY. 
40   - In most cases just by letting LU -> PY or UL -> PY.
41   - Special cases: RLU -> PYR (also PYRGET, PYRSET), KLU -> PYK,
42     PLU -> PYP. Also commonblock and internal variables of the form
43     *RLU* are replaced by *RPY* (where * represents "wildcard" 
44     characters).
45   - To declare integer functions, a line INTEGER PYK,PYCHGE,PYCOMP 
46     has been added at the beginning of all routines.
47   - To avoid a name clash, LUXTOT becomes PYXTEE.
48   - Some comment lines in code have been modified; also the name
49     JETSET becomes PYTHIA.
50
51 * Persons who have code that relies on the /LUJETS/ single precision
52   commonblock could easily write a translation routine to copy the
53   /PYJETS/ double precision information to /LUJETS/. In fact, only
54   the LUGIVE and LULOGO routines of JETSET have access to some PYTHIA
55   commonblocks, and therefore these are the only ones that need to
56   be modified if one, for some reason, would like to combine the 
57   new PYTHIA with the old JETSET routines. Similarly, it would only
58   require minor changes in the PYHEPC routine code to allow the
59   /HEPEVT/ commonblock to be in single precision, as before.  
60
61 -----------------------------------------------------------------------
62
63 DOUBLE PRECISION
64
65 * Conversion from single to double precision. 
66   - All real constants by explicitly exchanging E for D where present
67     and else adding D0. 
68   - All real variables with an IMPLICIT DOUBLE PRECISION(A-H, O-Z)
69     at beginning of all routines.
70   - Commonblocks with an odd number of integers before real variables
71     have been padded with a dummy integer variable or reordered:
72       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
73       COMMON/PYSUBS/MSEL,MSELPD,MSUB(200),KFIN(2,-40:40),CKIN(200)
74       COMMON/PYUPPR/NUP,KUP(20,7),NFUP,IFUP(10,2),PUP(20,5),Q2UP(0:10)
75       COMMON/PYINT5/NGENPD,NGEN(0:200,3),XSEC(0:200,3)
76     (In the HARD PROCESSES section below is described further changes;
77     specifically MSUB, NGEN and XSEC are expanded to 500 processes.)
78   - Obsolete conversions with DBLE(), SNGL() etc. have been removed.
79   - pi, 2pi and GeV <-> fm conversion given with more decimals.
80   - Some blanks have been removed and lines reformatted or split when 
81     lines have become too long.
82   - PYUPDA(3,..) writes values with D0 added.
83   - Commonblock /PYINT9/ with differential cross section sum in double
84     precision is superfluous and has been removed. 
85   - Note that in the Fortran 77 standard COMPLEX cannot be defined
86     as double precision. COMPLEX is used only very sparingly, and only
87     in the PYRESD, PYRAND and PYSIGH routines. Some small pieces of code 
88     therefore still use single precision. If you have a compiler
89     option for automatic promotion of single to double you are
90     welcome to use it to handle also these parts, but otherwise 
91     they should not be harmful.
92
93 -----------------------------------------------------------------------
94
95 PARTICLE CODES AND DATA
96
97 * Some particles have been renamed:
98    7 from l      to b' 
99    8 from h      to t' 
100   17 from chi    to tau' 
101   18 from nu_chi to nu'_tau 
102   25 from H      to h
103   35 from H'     to H
104   Furthermore, in all names where a tilde was previously used to
105   indicate an antiparticle, the previous alternative 'bar' is now
106   used throughout, both in particle names and process titles 
107   (to avoid confusion with supersymmetry, see below).   
108
109 * Some particle codes have been changed according to
110   the "LEP 2 standard" (see LEP 2 workshop proceedings):
111   psi'     from 30443 to  100443
112   Upsilon' from 30553 to  100553  
113   d*       from     7 to 4000001
114   u*       from     8 to 4000002
115   e*-      from    17 to 4000011
116   nu*_e    from    18 to 4000012
117   The codes 7, 8, 17 and 18 are now exclusively used for fourth
118   generation fermions. The switch MSTP(6) is then superfluous.
119   PYRESD and other code pieces have been rewritten to take into
120   account the change.  
121
122 * Supersymmetric particle codes have been introduced according to
123   the "LEP 2 standard" (see LEP 2 workshop proceedings):
124   1000001 ~d_L         2000001 ~d_R         1000021 ~g
125   1000002 ~u_L         2000002 ~u_R         1000022 ~chi_10
126   1000003 ~s_L         2000003 ~s_R         1000023 ~chi_20
127   1000004 ~c_L         2000004 ~c_R         1000024 ~chi_1+
128   1000005 ~b_1         2000005 ~b_2         1000025 ~chi_30
129   1000006 ~t_1         2000006 ~t_2         1000035 ~chi_40
130   1000011 ~e_L-        2000011 ~e_R-        1000037 ~chi_2+
131   1000012 ~nu_eL       2000012 ~nu_eR       1000039 ~G  
132   1000013 ~mu_L-       2000013 ~mu_R-
133   1000014 ~nu_muL      2000014 ~nu_muR
134   1000015 ~tau_1-      2000015 ~tau_2-
135   1000016 ~nu_tauL     2000016 ~nu_tauR
136   In the third generation the left and right states are assumed
137   to mix to nontrivial mass eigenstates, while mixing is not included
138   in the first two. Note that all sparticle names begin with a tilde.
139   Default masses are arbitrary and branching ratios not set at all. 
140   This is taken care of at initialization if IMSS(1) is positive 
141   (see below).  
142
143 * A hint on large particle numbers: if you want to avoid mistyping
144   the number of zeros, it may pay off to define a line like
145       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KEXCIT=4000000)
146   at the beginning of your program and then refer to particles as
147   KSUSY1+1 = ~d_L and so on. This then also agrees with the internal
148   notation.
149
150 * A number of technicolour particle codes have been added:
151   51 pi_tech0          54 rho_tech0   
152   52 pi_tech+          55 rho_tech+
153   53 pi'_tech0         56 omega_tech0
154
155 * Some new particle codes for doubly charged Higgs production in
156   left-right-symmetric scenarios.
157   61 H_L++             64 nu_Re 
158   62 H_R++             65 nu_Rmu        
159   63 W_R+              66 nu_Rtau
160   The indices _L and _R indicate belonging to left or right SU(2)
161   gauge group. 
162
163 * Top and fourth generation hadrons are gone. Henceforth the t, b' and
164   t' quarks are always assumed to decay before they would have time to
165   hadronize.
166   - All t, b' and t' hadron codes are unknown to the program.
167     (In a pinch, such a hadron could be represented e.g. by a string 
168     with a top quark and an antiquark or diquark, with string mass
169     equated to the expected hadron mass.)
170   - The MSTP(48) and MSTP(49) switches are removed; decay treatment is 
171     as the old default option 2. 
172   - Extra code has been inserted in PYEVNT and PYEXEC to decay any
173     leftover resonances (including these quarks) before fragmentation 
174     routines are called. 
175   - The particles codes 86 - 89, previously used for generic t, b' and
176     t' hadron decays are gone. 
177   - The decay channel of particle 17 to 89 is removed, and PYRESD
178     changed accordingly. 
179   - The PYLIST(11) listing is changed. 
180   - Also the PYLIST(12) listing is changed. The MSTU(14) switch is gone;
181     restrictions on which codes are listed can only be applied with
182     MSTU(1) and MSTU(2). 
183   - In the decay description, matrix element codes 45 and 46 are now
184     superfluous. MSTJ(25) and MSTJ(27) are no longer used, and 
185     MSTJ(23) is less used.
186   - PYTEST is reduced so it does not involve these hadrons.
187
188 * There is a new scheme to relate the standard KF codes with the 
189   compressed KC codes. 
190   - We remind that KF codes essentially follow the PDG standard for 
191     particle numbering, and with the introduction of SUSY now range up 
192     to seven-digit codes (plus a sign). They therefore cannot be used to 
193     directly access information in particle data tables. The compressed 
194     KC codes range between 1 and 500, and give the index to the KCHG,
195     PMAS, MDCY, CHAF and MWID arrays.
196   - Each KF code known to the program is now one-to-one associated 
197     with a KC code; the only doublevaluedness left is that both 
198     the particle KF and its antiparticle -KF (if existing) is mapped
199     to the same KC. This specifically means that all charm and bottom
200     hadrons and all diquarks now are separately defined.
201   - Whereas KF codes below 100 still obey KC=KF, the mapping of codes
202     above 100 is completely changed. The mapping is no longer hardcoded 
203     in PYCOMP, but defined by the fourth component of the KCHG array 
204     (see below). Therefore it can be changed or expanded during the 
205     course of a run, either by PYUPDA calls or by direct user 
206     intervention.
207  
208 * The KCHG array in /PYDAT2/ has been expanded with a fourth component,
209   where KCHG(KC,4)=KF. It thus gives the inverse mapping from KC codes
210   to KF ones (see above).  
211
212 * The PYCOMP code has been completely rewritten. It gives the direct
213   mapping from the KF codes to the KC ones (see above). 
214   - Internally the PYCOMP uses a binary search in a table, with KF codes 
215     arranged in increasing order, based on the KCHG(KC,4) array. This 
216     table is constructed the first time PYCOMP is called, at which time 
217     MSTU(20) is set to 1. In case of a user change of the KCHG(KC,4) 
218     array one should reset MSTU(20)=0 to force a reinitialization at the
219     next PYCOMP call (this is automatically done in PYUPDA calls).
220     To speed up execution, the latest (KF,KC) pair is kept in memory
221     and checked before the standard binary search.
222   - Code has been changed thoughout the program to be compatible with
223     this new mode of PYCOMP operation.
224
225 * Particle data is now stored and read out for each particle separately.
226   - PYMASS uses tables of charm and bottom hadron masses rather than
227     mass formulae.
228   - The array CHAF(500) has been expanded to CHAF(500,2), where the
229     first component gives the particle name and the second the
230     antiparticle one (where existing). 
231   - PYNAME accesses ready-constructed names rather than constructs
232     the names from scratch. 
233   - PYCHGE accesses charges directly from the KCHG(KC,1) array. 
234
235 * PYUPDA has been changed.
236   - Option 1 writes out a table of all particle codes defined.
237     For each particle is given its KF code, particle and antiparticle
238     names in CHAF , the three first KCHG components, the PMAS components,
239     MWID (see below) and the first MDCY component. The information on
240     decay channels is unchanged, but the format expanded.
241   - Option 2 reads in the table written in the,form described above,
242     and replaces all existing data, including the KF<->KC mapping, 
243     with the new ones.
244   - Option 3 reads in a table, like option 2, but uses it as a 
245     complement to rather than a replacement of existing data.
246     The input file should therefore only contain new particles and
247     particles with changed data. New particles are added to the
248     bottom of the KC and decay channel tables. Changed particles
249     retain their KC codes and hence the position of particle data, but 
250     their old decay channels are removed, this space is recuperated,
251     and new decay channels are added at the end.
252   - Option 4 corresponds to the old option 3, i.e. writes existing
253     data to DATA statements for inclusion in the default program 
254     code.  
255
256 * The maximum number of decay channels has been expanded from 2000
257   to 4000; this affects the arrays MDME, BRAT and KFDP in PYDAT3,
258   and MSTU(7).
259
260 * PYLIST, PYSTAT and PYUPDA are changed to allow for the larger 
261   particle codes that may now appear.
262
263 -----------------------------------------------------------------------
264
265 RESONANCE DECAYS
266
267 * The dimensions of the WDTP and WDTE return arrays of PYWIDT have 
268   been expanded from a maximum of 40 to 200 decay channels.
269
270 * PYWIDT has been modified so that it returns total and partial widths
271   in units of GeV. Previously most widths were given in dimensionless
272   units, with an extra multiplicative factor added elsewhere, e.g. in
273   PYINRE or PYSIGH. Therefore also these routines are modified. Also
274   VINT(117) is now in dimensions of GeV.
275
276 * Commonblock PYINT4 is completely reorganized as
277       COMMON/PYINT4/MWID(500),WIDS(500,5)  
278   The WIDP and WIDE arrays were essentially only used by PYSTAT(2)
279   and have been eliminated.
280   Where before the resonances could only be found in the range
281   21:40, in the new description any compressed code KC between
282   1 and 500 can be used to represent a resonance.
283
284   MWID(KC) gives the character of particle with compressed code KC,
285       mainly as used in PYWIDT to calculate widths of resonances
286       (not necessarily at the nominal mass).
287       = 0 : an ordinary particle; not to be treated as resonance.
288       = 1 : a resonance for which the partial and total widths
289           (and hence branching ratios) are dynamically calculated 
290           in PYWIDT calls; i.e. special code has to exist for each 
291           such particle. The effects of allowed/unallowed secondary
292           decays are included, both in the relative composition
293           of decays and in the process cross section.
294       = 2 : The total width is taken to be the one stored in PMAS(KC,2)
295           and the relative branching ratios the ones in BRAT(IDC) for 
296           decay channels IDC. There is then no need for any special
297           code in PYWIDT to handle a resonance. During the run,
298           the stored PMAS(KC,2) and BRAT(IDC) values are used to 
299           calculate the total and partial widths of the decay channels. 
300           Some extra information and assumptions are then used.
301           Firstly, the stored BRAT values are assumed to be the full 
302           branching ratios, including all possible channels and
303           all secondary decays. The actual relative branching fractions 
304           are modified to take into account that the simulation of some
305           channels may be switched off (even selectively for a particle
306           and an antiparticle), as given by MDME(IDC,1), and that
307           some secondary channels may not be allowed, as expressed by
308           the WIDS factors. This also goes into process cross sections.
309           Secondly, it is assumed that all widths scale like sqrt(shat)/m, 
310           the ratio of the actual to the nominal mass. A further nontrivial 
311           change as a function of the actual mass can be set for each 
312           channel by the MDME(IDC,2) value, see below.
313       = 3 : a hybrid version of options 1 and 2 above. At initialization
314           the PYWIDT code is used to calculate PMAS(KC,2) and BRAT(IDC)
315           at the nominal mass of the resonance. Special code must then 
316           exist in PYWIDT for the particle. The PMAS(KC,2) and BRAT(IDC) 
317           values overwrite the default ones. In the subsequent generation 
318           of events, the simpler scheme of option 2 is used, thus saving
319           some execution time. 
320       Note: the Z and Z' cannot be used with options 2 and 3, since the
321           more complicated interference structure implemented for those
322           particles is only handled correctly for option 1. 
323
324   WIDS(KC,J) : gives the suppression factor to cross sections caused 
325       by the closing of some secondary decays, as calculated in PYWIDT.
326       Is built up recursively from the lightest particle to the 
327       heaviest one at initialization, with the exception that W and Z 
328       are done already from the beginning (since these often are
329       forced off the mass shell). WIDS can go wrong in case you 
330       have perverse situations where the branching ratios vary
331       rapidly as a function of energy, across the resonance shape.
332       This then influences process cross sections.  
333       The J components store information according to
334       = 1 : a (matched) resonance-antiresonance pair or two identical
335             resonances (e.g. W+W- or Z0Z0).
336       = 2 : a single resonance (e.g. W+ or Z0).
337       = 3 : a single antiresonance (e.g. W-).
338       = 4 : a (matched) resonance-resonance pair, for particle which
339           has a nonidentical antiparticle (e.g. W+W+).
340       = 5 : a (matched) antiresonance-antiresonance pair (e.g. W-W-).
341
342 * The MDME(IDC,2) matrix element codes for a specific decay channel have
343   been expanded with further values that can be used for decay channels
344   treated by the PYRESD/PYWIDT decay machinery. These codes have no
345   meaning in the framework of ordinary particle decays in PYDECY.
346       = 50 : (default behaviour, also obtained for any other code value
347           apart from the ones listed below) do not include any special
348           threshold factors. That is, a decay channel is left open even 
349           if the sum of daughter nominal masses is above the mother 
350           actual mass, which is possible if at least one of the daughters 
351           can be pushed off the mass shell.
352       = 51 : a step threshold, i.e. a channel is switched off when
353           the sum of daughter nominal masses is above the mother actual
354           mass.
355       = 52 : a beta-factor threshold, i.e. 
356           sqrt( (1-m1**2/m**2-m2**2/m**2)**2 - 4*m1**2*m2**2/m**4),
357           assuming that the values stored in PMAS(KC,2) and BRAT(IDC)
358           did not include any threshold effects at all.      
359       = 53 : as =52, but assuming that PMAS(KC,2) and BRAT(IDC) did
360           include the threshold effects, so that the weight should be
361           beta(at the actual mass)/beta(at the nominal mass).
362       = 54 - 59 : free for further options.        
363
364 * VINT(91), VINT(92) are obsolete and replaced by WIDS(24,4), WIDS(24,5).
365
366 * The decay angles in H -> Z0 Z0 -> 4 fermions were previously selected
367   in the same way as for H -> W+ W- -> 4 fermions. Now the correct angular
368   correlations are included also for this case. Reference:
369   O. Linossier and R. Zitoun, internal ATLAS note and private communication. 
370
371 * PYRESD now takes an argument IRES. The standard call from PYEVNT, 
372   for the hard process, has IRES=0, and then finds resonances to be
373   treated based on the subprocess number ISUB. In case of a nonzero
374   IRES only the resonance in position IRES of the event record is 
375   considered. This is used by PYEVNT and PYEXEC to decay leftover
376   resonances. (Example: a b -> W + t branching may give a t quark as
377   beam remnant.)
378
379 * Now also three decay products can be handled by PYRESD.
380
381 * CKIN(49) and CKIN(50) have been introduced to allow minimum mass
382   limits to be passed from PYRESD to PYOFSH. They are used for
383   tertiary and higher resonances, i.e. those not controlled by 
384   CKIN(41)-CKIN(48). They need not be touched by the user.   
385
386 *  An approximate 1 - 2.5 alpha_s/pi QCD correction factor has been
387    introduced for the width of the top decay t -> b + W.
388
389 * New default behaviour of the Higgs resonance shape.
390   MSTP(49) : (D=1) assumed variation of the Higgs width as a function
391       of the actual mass mhat = sqrt(shat) and the nominal mass m_H. 
392       = 0 : the width is proportional to mhat**3; thus the high-mass
393           tail of the Breit-Wigner is enhanced.
394       = 1 : the width is proportional to m_H**2 * mhat. For a fixed
395           Higgs mass m_H this means a width variation across the 
396           Breit-Wigner more in accord with other resonances (such as
397           the Z0).  This alternative gives more emphasis to the
398           low-mass tail, where the parton distributions are peaked
399           (for hadron colliders). This option is favoured by 
400           resummation studies [M. Seymour, Phys. Lett. B354 (1995)
401           409].
402       Note : this switch does not affect processes 71 - 77, where a 
403           fixed Higgs width is used in order to control cancellation
404           of divergences.   
405
406 -----------------------------------------------------------------------
407
408 HARD PROCESSES
409
410 * The maximum number of processes has been expanded from 200 to 500;
411   this affects MSUB, ISET, KFPR, COEF, NGEN, XSEC and PROC in 
412   several commonblocks.
413
414 * SUSY processes have been introduced according to the SPYTHIA program.
415   - Look in the publication
416         SPYTHIA: A Supersymmetric Extension of PYTHIA 5.7
417         S. Mrenna, Computer Physics Commun. 101 (1997) 232
418         (hep-ph/9609360)
419     for a description of the physics that has been implemented.
420   - The list of new processes and process numbers is according to 
421     tables 2 and 3 in the SPYTHIA manual. Also the MSEL values
422     of table 4 can be used.
423   - Switches and free parameters that can be used to select a wide
424     variety of SUSY scenarios are accessed in 
425         COMMON/PYMSSM/IMSS(0:99),RMSS(0:99)   
426     according to the description in SPYTHIA manual section 3.1.
427   - The notation for sparticles follows the LEP 2 standard outlined 
428     above, and thus disagrees with the one in the SPYTHIA manual.
429   - The supersymmetric code has largely been taken over unchanged
430     from SPYTHIA, but a number of minor changes and bug fixes have 
431     been introduced. As examples, the sparticle mass selection has 
432     been improved, as has the selection of showering parton system.
433     The strategy for the selection of slepton and squark flavour, in
434     processes with several flavours allowed, has also been changed. 
435   - Some routine and commonblock names have been changed, but none
436     of the major ones listed in the SPYTHIA manual. The dependence
437     on CERN library routines has been eliminated.
438   - A major administrative change is that it is now possible to set
439     allowed decay channels of sparticles using the MDME array,
440     as for ordinary resonances, and have this reflected in the
441     process cross sections.
442   - One difference between the SUSY simulation and the other parts of 
443     the program is that it is not beforehand known which sparticles
444     may be stable. Normally this would mean either the chi_1 or the
445     gravitino, but in principle also other sparticles could be
446     stable. The ones found to be stable have their MWID(KC) and 
447     MDCY(KC,1) values set zero at initialization. If several
448     PYINIT calls are made in the same run, with different SUSY
449     parameters, the ones set zero above are not necessarily set
450     back to nonzero values (the exception is chi_1), since the
451     original values are not saved anywhere. This may then have to
452     be done by hand, or else some particles that ought to decay will
453     not do that. 
454   - Bottom squark production is now treated separately as for
455     the top squark.  However, there are more processes because bottom
456     is in the PDF.  The new processes are:
457     281   b q -> ~b_1 ~q_L    (q not b)
458     282   b q -> ~b_2 ~q_R
459     283   b q -> ~b_1 ~q_R + ~b_2 ~q_L
460     284   b qbar -> ~b_1 ~q_Lbar
461     285   b qbar -> ~b_2 ~q_Rbar
462     286   b qbar -> ~b_1 ~q_Rbar + ~b_2 ~q_Lbar
463     287   q qbar -> ~b_1 ~b_1bar
464     288                2    2
465     289   g g -> ~b_1 ~b_1bar
466     290             2    2
467     291   b b -> ~b_1 ~b_1
468     292             2    2
469     293             1    2
470     294   b g -> ~b_1 ~g
471     295             2
472     296   b bbar -> ~b_1 ~b_2bar + ~b_1bar ~b_2
473     MSEL = 45 has been added specifically to switch on these processes.
474   - New parameter 
475     IMMS(5) : (D=0) allows the user to set the stop, sbottom, and stau
476     masses and mixings by hand. 
477         = 0 : no, the program calculates itself.
478         = 1 : Yes, calculate from given input. In that case, 
479             RMMS(10) = lightest stop, RMSS(12) = heaviest stop,
480             RMSS(11) = lightest sbottom, RMSS(13) = lightest stau, 
481             RMSS(14) = heaviest stau, and RMSS(26,27,28) are the 
482             (1,1) elements of the (2x2) mixing matrix for sbottom, 
483             stop, and stau.
484
485 * Higgs pair production now added as explicit processes. (Since before
486     some processes exist as Z' decay modes, where the Z' part can be 
487     switched off to simulate the expected behaviour within the MSSM.)
488     297   q qbar' -> H+/- h0
489     298   q qbar' -> H+/- H0
490     299   q qbar -> A0 h0
491     300   q qbar -> A0 H0
492     301   q qbar -> H+ H-
493
494 * A new machinery has been introduced to generate the spectrum of
495   transverse and longitudinal photons in a lepton beam, and to
496   convolute that with the appropriate matrix elements, including
497   the virtuality of the photons, see C. Friberg and T. Sjostrand, 
498   Eur. Phys. J. C 13 (2000) 151.  
499   - In order to obtain it, the PYINIT beam or target code should
500     be given in the form 'gamma/lepton', where lepton can be either
501     of e-, e+, mu-, mu+, tau- or tau+. Thus,
502     for HERA : BEAM,TARGET = 'gamma/e-','p'
503     for LEP  :             = 'gamma/e-','gamma/e+'
504     Kinematics information in the PYINIT call should refer to the
505     full energy available, with the program itself generating the 
506     fraction given to the photon(s).
507   - The documentation section at the beginning of the event record
508     has been expanded to reflect the new layer of administration. 
509     Positions 1 and 2 contain the original beam particles, e.g. 
510     e and p (or e+ and e-). In position 3 (and 4 for e+e-) 
511     is (are) the scattered outgoing lepton(s). Thereafter comes 
512     the normal documentation, but starting from the photon rather
513     than a lepton. For ep, this means 4 and 5 are the gamma* and p,
514     6 and 7 the shower initiators, 8 and 9 the incoming partons to
515     the hard interaction, and 10 and 11 the outgoing ones. Thus the
516     documentation is 3 lines longer (4 for e+e-) than normally.
517   - A number of new CKIN cuts have been introduced to restrict
518     the range of kinematics for the photons generated off the 
519     lepton beams. In each quartet of numbers, the first two corresponds
520     to the range allowed on incoming side 1 (beam) and the last two
521     to side 2 (target). The cuts are only applicable for a lepton
522     beam. Note that the x and Q2 (P2) variables are the basis
523     for the generation, and so can be restricted with no loss of
524     efficiency. For leptoproduction the W is uniquely given by the
525     one x value of the problem, so here also W cuts are fully efficient.
526     The other cuts may imply a slowdown of the program, but not as much 
527     as if equivalent cuts only are introduced after events are fully 
528     generated.
529     CKIN(61) - CKIN(64) : (D=0.0001,0.99,0.0001,0.99) allowed range for
530         the energy fractions x that the photon take of the respective 
531         incoming lepton energy. These fractions are defined in the
532         cm frame of the collision, and differ from energy fractions 
533         as defined in another frame. (Watch out at HERA!) In order to 
534         avoid some technical problems, absolute lower and upper limits 
535         are set internally at 0.0001 and 0.9999.
536     CKIN(65) - CKIN(68) : (D=0.,-1.,0.,-1. GeV^2) allowed range for the
537         spacelike virtuality of the photon, conventionally called either 
538         Q2 or P2, depending on process. A negative number means that the
539         upper limit is inactive, i.e. purely given by kinematics. A nonzero
540         lower limit is implicitly given by kinematics constraints. 
541     CKIN(69) - CKIN(72) : (D=0.,-1.,0.,-1.) allowed range of the
542         scattering angle theta of the lepton, defined in the cm frame
543         of the event. (Watch out at HERA!) A negative number means that 
544         the upper limit is inactive, i.e. equal to pi.
545     CKIN(73) - CKIN(76) : (D=0.0001,0.99,0.0001,0.99) allowed range for 
546         the lightcone fraction y that the photon take of the respective 
547         incoming lepton energy. The lightcone is defined by the 
548         four-momentum of the lepton or hadron on the other side of the
549         event (and thus deviates from true lightcone fraction by mass
550         effects that normally are negligible). The y value is related to 
551         the x and Q2 (P2) values by y = x + Q2/s if mass terms are
552         neglected.   
553     CKIN(77), CKIN(78) : (D=2.,-1. GeV) allowed range for W, i.e. either
554         the photon-hadron or photon-photon invariant mass. A negative 
555         number means that the upper limit is inactive. 
556   - This machinery cannot be combined with the variable-energy option 
557     obtainable for MSTP(171)=1. The reason is that a variable-energy
558     machinery is now used internally for the gamma-hadron or gamma-gamma
559     subsystem, with some information saved at initialization for the full
560     energy.     
561   - Internally, some new variables are used:
562     MINT(141), MINT(142) : KF code for incoming lepton beam or target 
563         particles, while MINT(11) and MINT(12) is then the photon code.
564         A nonzero value is the main check whether the photon emission
565         machinery should be called at all.    
566     MINT(143) : the number of tries before a successful kinematics 
567         configuration is found in PYGAGA. Used for the cross section
568         updating in PYRAND.
569     VINT(301) : cm energy for the full collision, while VINT(1) 
570         gives the gamma-hadron or gamma-gamma subsystem energy.
571     VINT(302) : full squared cm energy, while VINT(2) gives the subsystem
572         squared energy.
573     VINT(303), VINT(304) : mass of beam or target lepton, while VINT(3) 
574         or VINT(4) give the mass of a photon  emitted off it.
575     VINT(305), VINT(306) : x values, i.e. respective photon energy 
576         fractions of the incoming lepton in the cm frame of the event.
577     VINT(307), VINT(308) : Q2 or P2, virtuality of the respective photon
578         (thus the square of VINT(3), VINT(4)).    
579     VINT(309), VINT(310) : y values, i.e. respective photon lightcone 
580         energy fraction of the lepton. 
581     VINT(311), VINT(312) : theta, scattering angle of the respective 
582         lepton in the cm frame of the event.   
583     VINT(313), VINT(314) : phi, azimuthal angle of the respective 
584         scattered lepton in the cm frame of the event. 
585     VINT(319) : photon flux factor in PYGAGA for current event.
586     VINT(320) : photon flux factor in PYGAGA at initialization.  
587   - Some of these values are also saved in the MSTI and PARI arrays at
588     the end of the event generation. (In the case of pileup events,
589     values stored here refer to the first event, while the MINT/VINT
590     ones are for the latest one, as usual.)
591     MSTI(71), MSTI(72) : KF code for incoming lepton beam or target 
592         particles, while MSTI(11) and MSTI(12) is then the photon code.
593     PARI(101) : cm energy for the full collision, while PARI(11) 
594         gives the gamma-hadron or gamma-gamma subsystem energy.
595     PARI(102) : full squared cm energy, while PARI(12) gives the subsystem
596         squared energy.
597     PARI(103), PARI(104) : x values, i.e. respective photon energy 
598         fractions of the incoming lepton in the cm frame of the event.
599     PARI(105), PARI(106) : Q2 or P2, virtuality of the respective photon
600         (thus the square of PARI(3), PARI(4)).    
601     PARI(107), PARI(108) : y values, i.e. respective photon lightcone 
602         energy fraction of the lepton. 
603     PARI(109), PARI(110) : theta, scattering angle of the respective 
604         lepton in the cm frame of the event.   
605     PARI(111), PARI(112) : phi, azimuthal angle of the respective 
606         scattered lepton in the cm frame of the event.   
607   - A new routine has been added for internal use:
608     SUBROUTINE PYGAGA(IGA)
609     IGA = 1 : call at initialization to set up x and Q2 limits etc.
610         = 2 : call at maximization step to give estimate of maximal 
611             photon flux factor.
612         = 3 : call at the beginning of the event generation to select 
613             the kinematics of the photon emission and to give the
614             flux factor. 
615         = 4 : call at the end of the event generation to set up the
616             full kinematics of the photon emission.
617   - Since there are currently no processes associated with resolved
618     longitudinal photons, the effect of these can be approximated by
619     some nonzero MSTP(17) and PARP(165). (Additionally, PARP(167) or
620     PARP(168) may need to be set.) 
621     MSTP(17) : (D=4) possibility of a extrafactor for resolved processes, 
622         to approximately take into accound the effects of longitudinal 
623         photons. Given on the form
624         R = 1 + PARP(165) * r(Q^2,mu^2) * f_L(y,Q^2)/f_T(y,Q^2).
625         Here the 1 represents the basic transverse contribution,
626         PARP(165) is some arbitrary overall factor, and f_L/f_T 
627         the (known) ratio of longitudinal to transverse photon 
628         flux factors. The arbitrary function r depends on the photon  
629         virtuality Q^2 and the hard scale mu^2 of the process.
630         = 0 : No contribution, i.e. r=0.
631         = 1 : r = 4 * mu^2 * Q^2 / (mu^2 + Q^2)^2.
632         = 2 : r = 4 * Q^2 / (mu^2 + Q^2).
633         = 3 : r = 4 * Q^2 / (m_{rho}^2 + Q^2). 
634         = 4 : r = 4 * m_V^2 * Q^2 / (m_V^2 + Q^2)^2.
635         = 5 : r = 4 * Q^2 / (m_V^2 + Q^2).
636         In options 4 and 5 m_V is the vector meson mass for VMD 
637         and 2 * k_T for GVMD states. Since there is no mu dependence
638         for these options (as well as for =3) they also affect
639         minimum-bias cross sections, where mu would be vanishing. 
640         Currently the rho mass is used also in options 4 and 5, for 
641         simplicity.
642         NOTE: For a photon given by the gamma/e option in the PYINIT call, 
643             the y spectrum is dynamically generated and y is thus known
644             from event to event. For a photon beam in the PYINIT call,
645             y is unknown from the onset, and has to be provided by the 
646             user if any longitudinal factor is to be included. So long
647             as these values, in PARP(167) and PARP(168), are at their
648             default values, 0, it is assumed they have not been set and
649             thus the MSTP(17) and PARP(165) values are inactive.
650     PARP(165) : (D=0.5)  a simple multiplicative factor applied to the 
651         cross section for the transverse resolved photons, see above
652         in MSTP(17). No preferred value, but typically one could use 
653         PARP(165)=1 as main contrast to the no-effect =0, with the 
654         default arbitrarily chosen in the middle. 
655     PARP(167), PARP(168): (D=2*0) the longitudinal energy fraction 
656         y of an incoming photon, side 1 or 2, used in the R expression 
657         to evaluate f_L(y,Q^2)/f_T(y,Q^2). Need not be supplied when 
658         a photon spectrum is generated inside a lepton beam, but only 
659         when a photon is directly given as argument in the PYINIT call.
660     VINT(315), VINT(316): internal storage of the R factor above, for 
661         each of the two sides.
662     PARI(113), PARI(114); values of the R factors above, for each of 
663         the two sides.
664
665 * New total cross sections have been introduced into PYXTOT for 
666   Generalized Vector Meson Dominance (GVMD) states, and both VMD and 
667   GVMD parameterizations have been extended also to include virtual
668   photons. Further details in C. Friberg and T. Sjostrand, in 
669   preparation.
670   - The GVMD states are seen as a continuous spectrum, characterized
671     by the k_T scale of the gamma -> q + qbar branching, with
672     k_T stretching between k_0 and p_Tmin(W^2). The rate of fluctuation
673     into such states is given by perturbative QED, while the
674     hadronic cross section for a given state is assumed to obey geometric
675     scaling, i.e. fall off like k_rho^2/k_T^2 relative to a VMD state
676     for a real photon, where k_rho is a reference scale.
677   - The jet cross sections for these GVMD states are associated with 
678     the anomalous part of the photon structure function, just like
679     the homogeneous part is associated with the VMD states. 
680   - GVMD state also have "elastic" and diffractive cross sections 
681     obtained by the same scaling of VMD cross sections as indicated 
682     above for the total cross section. The mass selection of the 
683     GVMD state is according to dm^2/(m^2+Q^2)^2 between limits
684     2 k_0 < m < 2 p_Tmin(W^2), i.e. the mass is associated with
685     2 k_T of the state. See VINT(69), VINT(79) below. A GVMD state 
686     is bookkept as a diffractive state in event listing, even when
687     it scatters "elastically", since the subsequent hadronization
688     descriptions are very similar.
689   - Whether or not minimum bias events are simulated depends on the
690     CKIN(3) value. For a low CKIN(3), CKIN(3) < p_Tmin(W_init^2),
691     like the default value CKIN(3) = 0, low-pT physics is switched 
692     on together with jet production, with the latter properly 
693     eikonalized to be lower than the total one. For a high CKIN(3), 
694     CKIN(3) > p_Tmin(W_init^2), only jet production is included.
695     This is just like for hadron-hadron collisions, except that the
696     initialization energy scale W_init is selected in the allowed
697     W range rather than to be the full CM energy. When MSEL=2, also
698     elastic and diffractive events are simulated.
699   - Multiple interactions become possible in both the VMD and GVMD
700     sector, with the average number of interactions given by the 
701     ratio of the jet to the total cross section. Currently only
702     the simpler default scenario MSTP(82)=1 is implemented, i.e.
703     the more sophisticated variable-impact-parameter ones need further
704     physics studies and model development.
705   - For a virtual photon of virtuality Q^2, the total cross section is 
706     reduced by a dipole factor (m_rho^2/(m_rho^2 + Q^2))^2 for a VMD
707     state and by (4 k_T^2/(4 k_T^2 + Q^2))^2 for a GVMD one. That
708     is, the "mass" of a GVMD state is taken to be 2 k_T.Properly 
709     each VMD state should have own mass, but so far this has not been 
710     implemented. This would mainly be of relevance for J/psi, where 
711     however also other complications enter. 
712   - gamma* gamma* cross sections are obtained by simple multiplicative
713     factors as above, one for each photon, relative to rho rho events
714     (and other vector mesons).
715   - The primordial kT selection is described in the section on MSTP(66).
716     For clarity, we point out that elastic and diffractive events are
717     characterized by the mass of the diffractive states but without
718     any primordial kT, while jet production involves a primordial kT
719     but no mass selection. Both are thus not used at the same time,
720     but implicitly they are associated as m = 2 k_T.
721   - New or modified commonblock variables:
722     MSTP(15) : (D=0) modified default, to give same pT cutoff procedure 
723         as for VMD jet cross sections.
724     PARP(15) : (D=0.5 GeV) k_0 scale where GVMD k_T spectrum begins.
725     MINT(50) : now set = 1 also for anomalous states, to indicate that
726         total cross sections are defined for them.
727     VINT(67), VINT(68) : the mass of a VMD state; for GVMD photons
728         the VMD state with the equivalent flavour content.
729     VINT(69), VINT(70) : the actual mass of a VMD or GVMD state;
730         agrees with the above for VMD but is selected as a larger 
731         number for GVMD. Required for elastic and diffractive events.
732     VINT(63), VINT(64) : the squared (!) mass of the outgoing states;
733         for elastic events equal to VINT(69)^2 and VINT(70)^2 and for
734         diffractive events above that.    
735     VINT(154) : current p_Tmin(W^2) value; see section on UNDERLYING 
736         EVENTS for details. 
737     VINT(149) : the scaled value 4 p_Tmin(W^2)^2/W^2; denominator 
738         changed from s and therefore needs to be recalculated for each 
739         new event (like VINT(154)).
740     PARP(18) : (D=0.4 GeV) scale k_rho, such that the cross sections 
741         of a GVMD state of scale k_T is suppressed by a factor 
742         k_rho^2/k_T^2 relative to those of a VMD state. Order should be
743         m_rho/2, with some finetuning to fit data.
744     VINT(317) : dipole suppression factor in PYXTOT for current event.
745     VINT(318) : dipole suppression factor in PYXTOT at initialization. 
746     MSTP(66) : (D=5) see separate note below. 
747   - The MSTP(84) and MSTP(85) switches have been made obsolete by these
748     changes and no longer exist.    
749
750 * An additional suppression of resolved (VMD or GVMD) cross sections is 
751   introduced to compensate for an overlap with DIS processes in the 
752   region of intermediate Q^2 and rather small W^2.
753   - MSTP(20) : (D=3) suppression of resolved cross sections.
754         = 0 : no; used as is.
755         > 0 : yes, by a factor (W^2/(W^2 + Q_1^2 + Q_2^2))^MSTP(20).
756               (where Q_i^2 = 0 for an incoming hadron).
757   - The suppression factor is joined with the dipole suppression       
758     stored in VINT(317) and VINT(318).
759
760 * New processes have been introduced for incoming virtual (spacelike)
761   photons, as obtained e.g. in ep and e+e- collisions. These are
762   thus extensions of processes previously encoded for real photons.
763   - 131   f_i + gamma*_T -> f_i + g           (cf. proc 33)
764     132   f_i + gamma*_L -> f_i + g
765     133   f_i + gamma*_T -> f_i + gamma       (cf. proc 34)
766     134   f_i + gamma*_L -> f_i + gamma
767     135   g + gamma*_T -> f_i + fbar_i        (cf. proc 54)
768     136   g + gamma*_L -> f_i + fbar_i
769     137   gamma*_T + gamma*_T -> f_i + fbar_i (cf. proc 58)
770     138   gamma*_T + gamma*_L -> f_i + fbar_i 
771     139   gamma*_L + gamma*_T -> f_i + fbar_i 
772     140   gamma*_L + gamma*_L -> f_i + fbar_i 
773   - Here indices _T and _L represent transverse and longitudinal
774     photons, respectively. In the limit of vanishing virtuality,
775     the _T photon cross section approaches that for a real photon,
776     while the _L one vanishes.
777   - The virtuality of the photon or photons can be stored in P(1,5) 
778     and P(2,5), respectively, provided PYINIT is called with the 
779     'FIVE' option. A spacelike photon of virtuality Q**2 (or P**2,
780     depending on notational convention followed) would thus have
781     P(i,5) = -Q (or -P). The virtuality could be varied from one 
782     event to the next, but then it is convenient to initialize
783     for the lowest virtuality likely to be encountered.
784   - In several of the standard MSEL options, processes selected for 
785     real photons have been replaced by the corresponding processes 
786     for virtual ones. 
787
788 * Direct processes in the range of k_T values stretching between k_0 and 
789   p_Tmin(W^2) are, by an eikonalization process, associated with the 
790   low-pT part of the GVMD states above. Further details in C. Friberg 
791   and T. Sjostrand, in preparation.
792   - As a consequence, the minimum pT for direct processes should be 
793     increased from k_0 to p_Tmin(W^2). 
794   - New variable:
795     MSTP(18) : (D=3) choice of pTmin for direct processes:
796         = 1 : same as for VMD and GVMD states, as explained above.. 
797         = 2 : pTmin is chosen to be PARP(15), i.e. the old behaviour.
798             In this case, also parton distributions, jet cross sections
799             and alpha_strong values were dampened for small pT.
800         = 3 : as =1, but if the Q scale of the virtual photon is 
801             above the VMD/GVMD p_Tmin(W^2), pTmin is chosen equal to Q.
802             This is part of the strategy to mix in DIS processes at 
803             pT below Q, e.g. in MSTP(14)=30.
804
805 * New process 99 for DIS scattering, by photon exchange only. Thus, in 
806   this sense less powerful than process 10, but allows the use of the 
807   same photon flux machinery as for other gamma*-p and gamma*-gamma*
808   processes, and thus offers a unified description in the region of
809   intermediate Q2 values.
810   - Notice that it counts as a "total cross section" process, in the
811     sense that the hard subprocess in itself contains no high-pT
812     scale. Therefore, it will be switched off in event class mixes 
813     such as MSTP(14)=30 if CKIN(3) is above pTmin(W^2) and MSEL 
814     is not 2.  
815   - 99    f_i + gamma* -> f_i.
816   - Since the standard 2 -> 1 kinematics machinery is not relevant for
817     this process - shat = 0 - a new code ISET(ISUB)=8 is introduced
818     for the kinematics selection machinery in PYRAND, and a new routine
819     PYDISG for setting up the kinematics, beam remnants and showers.
820   - New variable to select DIS cross section.
821     MSTP(19) : (D=4) choice of partonic cross section in process 99.
822         = 0 : QPM answer 4 pi^2 alpha_em/Q^2 * 
823             \sum_q e_q^2 (x q(x,Q^2) + x qbar(x,Q^2))   
824             (with parton distributions frozen below the lowest Q
825             allowed in the parameterization). Note that this answer
826             is divergent for Q^2 -> 0 and thus violates gauge
827             invariance.
828         = 1 : QPM answer is modified by a factor Q^2/(Q^2 + m_rho^2)
829             to provide a finite cross section in the Q^2 -> 0 limit.
830             A minimal regularization recipe.
831         = 2 : QPM answer is modified by a factor Q^4/(Q^2 + m_rho^2)^2
832             to provide a vanishing cross section in the Q^2 -> 0 limit.
833             Appropriate if one assumes that the normal photoproduction 
834             description gives the total cross section for Q^2 = 0,
835             without any DIS contribution.
836         = 3 : as = 2, but additionally suppression by a parameterized
837             factor f(W^2,Q^2) (different for gamma*-p and gamma*-gamma*)
838             that avoids doublecounting the direct-process region where
839             p_T > Q. Shower evolution for DIS events is then also
840             restricted to be at scales below Q, whereas evolution all
841             the way up to W is allowed in the other options above.
842         = 4 : as = 3, but additionally include factor 1/(1-x) for
843             conversion from F_2 to sigma. This is formally required, 
844             but is only relevant for small W2 and therefore often 
845             neglected.
846   - MINT(107),MINT(108) = 4 denotes DIS photon on respective side.
847     MINT(123) = 8 denotes DIS*VMD/p or vice verse, = 9 DIS*anomalous   
848         or vice versa.
849     In MINT(41)-MINT(46), a DIS photon is treated same way as a direct 
850         one. 
851   - Many of the normal kinematical variables for 2 -> 2 processes are 
852     not defined for this process. The pT in PARI(17) is explicitly set 
853     =0, but some others may well contain irrelevant junk.        
854
855 * MSTP(14) is extended with new possibilities to select the nature 
856   of incoming virtual photons. The reason is that the existing 
857   options specify e.g. direct * VMD, summing over the possibilities
858   of which photon is direct and which anomalous. This is allowed
859   when the situation is symmetric, i.e. for two incoming real photons,
860   but not if one is virtual. Some of  the new options agree with 
861   previous ones, but are included to allow a more consistent pattern.
862   MSTP(14): (D=30) structure of incoming photon beam or target.
863       = 11 : direct * direct (see note 4).
864       = 12 : direct * VMD (i.e. first photon direct, second VMD). 
865       = 13 : direct * anomalous.
866       = 14 : VMD * direct.  
867       = 15 : VMD * VMD.  
868       = 16 : VMD * anomalous.  
869       = 17 : anomalous * direct.  
870       = 18 : anomalous * VDM.  
871       = 19 : anomalous * anomalous.
872       = 20 : a mixture of the nine above components, in the same 
873           spirit as =10 provides a mixture for real gammas (or
874           a virtual gamma on a hadron). For gamma-hadron, this 
875           option coincides with =10.
876       = 21 : direct * direct (see note 4).
877       = 22 : direct * resolved. 
878       = 23 : resolved * direct.
879       = 24 : resolved * resolved.     
880       = 25 : a mixture of the four above components, offering a 
881           simpler alternative to =20 in cases where the parton
882           distributions of the photon have not been split into VMD 
883           and anomalous components. For gamma-hadron, only two
884           components need be mixed.
885       = 26 : DIS * VMD/p.
886       = 27 : DIS * anomalous.
887       = 28 : VMD/p * DIS.
888       = 29 : anomalous * DIS.
889       = 30 : a mixture of all the 4 (for gamma*-p) or 13 (for
890           gamma*-gamma*) that are available, is as = 20 with the
891           DIS processes 26-29 mixed in.    
892   Note 1: The MSTP(14) options apply for a photon defined by a 'gamma'
893       or 'gamma/lepton' beam in the PYINIT call, but not to those 
894       photons implicitly obtained in a 'lepton' beam with the
895       MSTP(12)=1 option. This latter approach to resolved photons is
896       more primitive and is no longer recommended. 
897   Note 2: these new options are not needed and therefore not defined
898       for e-p collisions. The recommended 'best' values thus are
899       MSTP(14)=30, which also is the new default value.
900   Note 3: as a consequence of the appearance of new event classes, 
901       the MINT(122) and MSTI(9) code is not the same for gamma* gamma*
902       events as for gamma p, gamma* p or gamma gamma ones. 
903       Instead the code is 3*(icode_1 - 1) + icode_2, where icode is
904       1 for direct, 2 for VMD and 3 for anomalous/GVMD and indices 
905       refer to the two incoming photons. For gamma* p code 4 is DIS,
906       and for gamma* gamma* codes 10-13 corresponds to the MSTP(14)
907       codes 26-29. As before, MINT(122) and MSTI(9) are only defined
908       when several processes are to be mixed, not when generating one       
909       at a time. Also the MINT(123) code is modified (not shown here).
910   Note 4: The direct * direct event class excludes lepton pair 
911       production when run with the default MSEL=1 option (or MSEL=2), 
912       in order not to confuse users. You can obtain lepton pairs as well,
913       e.g. by running with MSEL=0 and switching on the desired processes
914       by hand.  
915
916 * MSTP(16) is new variable to select momentum variable of 
917   e -> gamma branching.
918   MSTP(16) (D=1) choice of definition of the fractional momentum 
919   taken by a photon radiated off a lepton. Enters in the flux
920   factor for the photon rate, and thereby in cross sections.
921   = 0 ; x, i.e. energy fraction in the rest frame of the event.
922   = 1 ; y, i.e. lightcone fraction.
923
924 * MSTP(32) : (D=8) has been expanded with new options for the choice 
925     of Q2 scale, specifically intended for processes with incoming
926     virtual photons. The new options are ordered from a "minimal"
927     dependence on the virtualities to a "maximal" one, based on
928     reasonable kinematics considerations. The old default value
929     MSTP(32)=2 forms the starting point, with no dependence at 
930     all, and the new default is some intermediate choice. 
931     Notation is that P1**2 and P2**2 are the virtualities of the 
932     two incoming particles, pT the transverse momentum of the 
933     scattering process, and m3 and m4 the masses of the two 
934     outgoing partons. For a direct photon, P**2 is the photon
935     virtuality and x=1. For a resolved photon, P**2 still refers
936     to the photon, rather than the unknown virtuality of the 
937     reacting parton in the photon, and x is the momentum fraction
938     taken by this parton.
939     = 6 : Q2 = (1 + x1*P1**2/shat + x2*P2**2/shat)*
940         (pT**2 + m3**2/2 + m4**2/2).
941     = 7 : Q2 = (1 + P1**2/shat + P2**2/shat)*
942         (pT**2 + m3**2/2 + m4**2/2).
943     = 8 : Q2 = pT**2 + (P1**2 + P2**2 +m3**2 + m4**2)/2.  
944     = 9 : Q2 = pT**2 + P1**2 + P2**2 +m3**2 + m4**2. 
945     = 10 : s (the full energy-squared of the process).
946     Note: options 6 and 7 are motivated by assuming that one
947         wants a scale that interpolates between that for small 
948         that and uhat for small uhat, such as 
949         Q2 = - that*uhat/(that+uhat). When kinematics for 
950         the 2 -> 2 process is constructed as if an incoming
951         photon is massless when it is not, it gives rise to a
952         mismatch factor 1 + P**2/shat (neglecting the other 
953         masses) in this Q2 definition, which is then what is
954         used in option 7 (with the neglect of some small 
955         cross-terms when both photons are virtual). When a
956         virtual photon is resolved, the virtuality of the
957         incoming parton can be anything from x*P**2 and upwards.
958         So option 6 uses the smallest kinematically possible 
959         value, while 7 is more representative of the typical 
960         scale. Option 8 and 9 are more handwaving extensions of
961         the default option, with 9 specially constructed to 
962         ensure that the Q2 scale is always bigger than P**2.  
963
964 * MSTP(66) has been expanded with new default option for the
965     selection of lower parton-shower cut-off (and primordial kT). 
966     MSTP(66) : (D=5) choice of lower cut-off for initial-state QCD
967         radiation in VMD or anomalous photoproduction events.
968         = 0 : the lower Q2 cut-off is the standard one in PARP(62)^2.
969         = 1 : for anomalous photons, the lower Q2 cut-off  is the 
970             larger of PARP(62)^2 and VINT(283) or VINT(284),
971             where the latter is the virtuality scale for the 
972             gamma -> q qbar vertex on the appropriate side of 
973             the event. The VINT values are selected logarithmically
974             even between PARP(15)^2 and the Q2 scale of the
975             parton distributions of the hard process.    
976         = 2 : extended option of the above, intended for virtual
977             photons. For VMD photons, the lower Q2 cut-off is the
978             larger of PARP(62)^2 and the P^2_{int} scale of the
979             SaS parton distributions. For anomalous photons,
980             the lower cut-off is chosen as for =1, but the
981             VINT(283) and VINT(284) are here selected logarithmically
982             even between P^2_{int} and the Q2 scale of the
983             parton distributions of the hard process.  
984         = 3 : simplified option, default in versions 6.143 -  6.147. 
985             The k_T of the anomalous/GVMD component is distributed
986             like 1/k_T^2 between k_0 and p_Tmin(W^2). Apart from
987             the change of the upper limit, this option works just 
988             like = 1.
989         = 4 : a stronger damping at large k_T, like 
990             dk_T^2/(k_T^2 + Q^2/4)^2 with 
991             k_0 < k_T < p_Tmin(W^2). Apart from this, 
992             it works like = 1.
993         = 5 : a k_T generated as in =4 is added vectorially with a 
994             standard Gaussian k_T generated like for VMD states.      
995             Ensures that GVMD has typical k_T's above those of VMD,
996             in spite of the large primordial k_T's implied by hadronic
997             physics. (Probably attributable to a lack of soft QCD
998             radiation in parton showers.)
999    
1000 * New processes
1001   - 146   e + gamma -> e*
1002     169   q + qbar -> e + e*
1003   - similar to existing processes 147,148 or 167,168 for q*.
1004
1005 * Several new processes for technicolour production.
1006   NOTE: as of version 6.126 changes/additions appear according
1007   to the next section.
1008   - 191   f_i + fbar_i -> rho_techni0
1009     192   f_i + fbar_j -> rho_techni+-
1010     193   f_i + fbar_i -> omega_techni0
1011     194   f_i + fbar_i -> f_k + fbar_k 
1012   - The first three processes are based on s-channel production of
1013     the respective resonance. All decay modes implemented can be
1014     simulated separately or in combination, in the standard fashion. 
1015     These include pairs of fermions, of gauge bosons, of technipions, 
1016     and of mixtures gauge bosons + technipions. 
1017   - Process 194 includes full interference between rho_techni0 and
1018     omega_techni0. It can only be used for one final-state flavour
1019     at a time. This flavour is set in KFPR(194,1).
1020   - The physics parameters of the technicolour scenario are:
1021     PARP(140) : (D=0.0) multiplicative fudge factor, entering 
1022         quadratically in the width for pi_tech+ -> W+ b bbar. 
1023     PARP(141) : (D=0.33333) sin(chi), sinus of mixing angle between 
1024         gague bosons and technipions in the decay of technirhos;
1025         if 0 the decay is entirely to technipions and if 1 entirely
1026         to gauge bosons.  
1027     PARP(142) : (D=82 GeV) F_T, decay constant of the technipion
1028         states; the technipion widths are proportional to 1/F_T^2.
1029     PARP(143) : (D=1.0) Q_U, charge of the up-type technifermions;
1030         the down-type ones have Q_D = Q_U - 1 and thus do not
1031         require a separate parameter.     
1032     PARP(144) : (D=4.0) N_TC, the number of technicolours, that
1033         enters in several cross sections and decay rates. 
1034     PARP(145) : (D=200 GeV) M_T, mass parameter for the decay
1035         omega_techni0 -> gamma + pi_techni0; the partial width
1036         is proportional to 1/M_T^2.
1037     PARP(146) - PARP(148) : (D=1.0, 1.0, 1.0) multiplicative fudge 
1038         factors, entering quadratically in the widths of technipions 
1039         to a fermion pair. The three numbers are for pi_tech0,
1040         pi_tech+ and pi_tech'0, respectively.
1041     PARP(149) - PARP(150) : (D=1.0, 0.0) multiplicative fudge factors, 
1042         entering linearly in the widths of technipions to a gluon pair. 
1043         The two numbers are for pi_tech0 and pi_tech'0, respectively. 
1044   - The main references are
1045     E. Eichten and K. Lane, Phys. Lett. B388 (1996) 803
1046     E. Eichten, K. Lane and J. Womersley, in preparation
1047
1048 * Starting with version 6.126, the simulation of the production and 
1049   decays of technicolor particles has been substantially upgraded.
1050   - The processes 149, 191, 192, and 193 are to be considered obsolete, 
1051     and are temporarily retained to allow cross checking with the new 
1052     processes.
1053   - Process 194 has been changed to more accurately represent the 
1054     mixing between the photon, Z, techni_rho0, and techni_omega 
1055     particles in the Drell-Yan process.  Process 195 is the analogous 
1056     process including W and techni_rho+/- mixing.  By default, the final 
1057     state fermions are e+ e- and e+/- nu_e, respectively.  These can be 
1058     changed through the parameters KFPR(194,1) and KFPR(195,1), 
1059     respectively (the KFPR value should represent a charged fermion).
1060   - The full set of recommended processes are:
1061     Drell--Yan (ETC == Extended TechniColor)
1062     194   f+fbar -> f'+fbar' (ETC)      
1063     195   f+fbar' -> f"+fbar"' (ETC)
1064     techni_rho0/omega
1065     361   f + fbar -> W_L+ W_L-         
1066     362   f + fbar -> W_L+/- pi_T-/+    
1067     363   f + fbar -> pi_T+ pi_T-       
1068     364   f + fbar -> gamma pi_T0       
1069     365   f + fbar -> gamma pi_T0'      
1070     366   f + fbar -> Z0 pi_T0          
1071     367   f + fbar -> Z0 pi_T0'         
1072     368   f + fbar -> W+/- pi_T-/+      
1073     charged techni_rho
1074     370   f + fbar' -> W_L+/- Z_L0      
1075     371   f + fbar' -> W_L+/- pi_T0     
1076     372   f + fbar' -> pi_T+/- Z_L0     
1077     373   f + fbar' -> pi_T+/- pi_T0    
1078     374   f + fbar' -> gamma pi_T+/-    
1079     375   f + fbar' -> Z0 pi_T+/-       
1080     376   f + fbar' -> W+/- pi_T0       
1081     377   f + fbar' -> W+/- pi_T0'      
1082   - All of the processes from 361 to 377 can be accessed at once
1083     by setting MSEL=50.
1084   - The production and decay rates depend on several "Straw Man"
1085     technicolor parameters:
1086     Techniparticle masses
1087     PMAS(51,1) : (D=110.0 GeV)  neutral techni_pi mass
1088     PMAS(52,1) : (D=110.0 GeV)  charged techni_pi mass
1089     PMAS(53,1) : (D=110.0 GeV)  neutral techni_pi' mass
1090     PMAS(54,1) : (D=210.0 GeV)  neutral techni_rho mass
1091     PMAS(55,1) : (D=210.0 GeV)  charged techni_rho mass
1092     PMAS(56,1) : (D=210.0 GeV)  techni_omega mass
1093     Note:  the rho and omega masses are not pole masses
1094     Lagrangian parameters
1095     PARP(141) :  (D= 0.33333)  sine of chi, the mixing angle between
1096         technipion interaction and mass eigenstates
1097     PARP(142) :  (D=82.0000 GeV)   F_T, the technipion decay constant
1098     PARP(143) :  (D= 1.3333)   Q_U, charge of up-type technifermion; 
1099         the down-type technifermion has a charge Q_D=Q_U-1   
1100     PARP(144) :  (D= 4.0000)   N_TC, number of technicolors; fixes the
1101         relative values of g_em and g_etc
1102     PARP(145) :  (D= 1.0000)   C_c, coefficient of the technipion decays
1103         to charm; appears squared in the decay rate  
1104     PARP(146) :  (D= 1.0000)   C_b, coefficient of the technipion decays
1105         to bottom; appears squared in the decay rate  
1106     PARP(147) :  (D= 0.0182)   C_t, coefficient of the technipion decays
1107         to top, estimated to be m_b/m_t; appears squared in the decay rate  
1108     PARP(148) :  (D= 1.0000)   C_tau, coefficient of the technipion decays
1109         to tau; appears squared in the decay rate  
1110     PARP(149) :  (D=0.00000)   C_pi, coefficient of technipion decays
1111         to gg
1112     PARP(150) :  (D=1.33333)   C_pi', coefficient of technipion' decays
1113         to gg
1114     ****Note the switch from PARP to PARJ****
1115     PARJ(172) :  (D=200.000 GeV) M_V, vector mass parameter for technivector
1116         decays to transverse gauge bosons and technipions
1117     PARJ(173) :  (D=200.000 GeV) M_A, axial mass parameter for technivector
1118         decays to transverse gauge bosons and technipions
1119     PARJ(174) :  (D=0.33300)   sine of chi', the mixing angle between
1120         the technipion' interaction and mass eigenstates 
1121     PARJ(175) :  (D=0.05000)   isospin violating technirho/techniomega
1122         mixing amplitude
1123   - As a final comment, it is worth mentioning that the decays products
1124     of the W and Z bosons are distributed according to phase space, 
1125     regardless of their designation as W_L/Z_L or transverse gauge bosons.  
1126     The exact meaning of longitudinal or transverse polarizations in this 
1127     case requires more thought.
1128   -  References:
1129      K. Lane, hep-ph/9903369
1130      K. Lane, hep-ph/9903372
1131
1132 * Several new processes for doubly charged Higgs production in
1133   left-right-symmetric models, with an additional righthanded SU(2)
1134   gauge group.
1135   - 341   l + l -> H_L++/--
1136     342   l + l -> H_R++/-- 
1137     343   l + gamma -> H_L++/-- + e-/+ 
1138     344   l + gamma -> H_R++/-- + e-/+ 
1139     345   l + gamma -> H_L++/-- + mu-/+ 
1140     346   l + gamma -> H_R++/-- + mu-/+ 
1141     347   l + gamma -> H_L++/-- + tau-/+
1142     348   l + gamma -> H_R++/-- + tau-/+
1143     349   f + fbar -> H_L++ + H_L--   
1144     350   f + fbar -> H_R++ + H_R-- 
1145     351   f_i + f_j -> f_k + f_l + H_L++/--
1146     352   f_i + f_j -> f_k + f_l + H_R++/--
1147   - Default model masses are
1148     code   name     mass (GeV)
1149      61   H_L++      200       
1150      62   H_R++      200             
1151      63   W_R+       750    
1152      64   nu_Re      750
1153      65   nu_Rmu     750   
1154      66   nu_Rtau    750
1155   - Main decays implemented are
1156      H_L++ -> l_i+ l_j+    (i, j generation index)
1157            -> W_L+ W_L+
1158      H_R++ -> l_i+ l_j+    
1159            -> W_R+ W_R+
1160      W_R+  -> q_i qbar_j   (assuming standard CKM matrix)
1161            -> l_i+ nu_Ri   (if kinematically allowed)
1162   - The physics parameters of this scenario are 
1163     PARP(181) - PARP(189) : (D = 0.1, 0.01, 0.01, 0.01, 0.1, 0.01,
1164         0.01, 0.01, 0.3) Yukawa couplings of leptons to H++, assumed 
1165         same for H_L++ and H_R++. Is a symmetric 3*3 array, where 
1166         PARP(177+3*i+j) gives the coupling to a lepton pair with 
1167         generation indices i and j. Thus the default matrix is 
1168         dominated by the diagonal elements and especially by the 
1169         tau-tau one.
1170     PARP(190) : (D=0.64) g_L = e/sin(theta_W).
1171     PARP(191) : (D=0.64) g_R, assumed same as g_L.
1172     PARP(192) : (D=5 GeV) v_L vacuum expectation value of the 
1173         left-triplet. The corresponding v_R is assumed given by
1174         v_R = sqrt(2) M_W_R / g_R and is not stored explicitly. 
1175   - The main references are
1176     K. Huitu, J. Maalampi, A. Pietil\"a and M. Raidal, 
1177         Nucl. Phys. B487 (1997) 27 and private communication;
1178     G. Barenboim, K. Huitu, J. Maalampi and M. Raidal, 
1179         Phys. Lett. B394 (1997) 132.    
1180
1181 * Two new processes for chi_c production.
1182   - 104   g + g -> chi_0c
1183     105   g + g -> chi_2c
1184   - These are the lowest-order equivalents of processes 87 and 89.
1185     Note that g + g -> chi_1c is forbidden, and so not included as
1186     a match to process 88.
1187   - Reference Bai83.
1188
1189 * Three new processes for J/psi production.
1190   - 106   g + g -> J/psi + gamma
1191     107   g + gamma -> J/psi + g
1192     108   gamma + gamma -> J/psi + gamma
1193   - All of these are closely related to the existing process 86,
1194     g + g -> J/psi + g; only the colour- and charge-related 
1195     prefactors differ in the matrix element expressions.
1196   - References: 
1197     106: M. Drees and C.S. Kim, Z. Phys. C53 (1991) 673
1198     107: E.L. Berger and D. Jones, Phys. Rev. D23 (1981) 1521
1199     108: H. Jung, private communication;
1200          H. Kharraziha, private communication
1201
1202 * Process 1 has been modified so that masses are included in the
1203   expression for the decay polar angle distribution.
1204
1205 * Processes 15, 19, 22, 30 and 35 have been corrected for an inconsistent
1206   use of width definition in the Breit-Wigner shape, which affected
1207   the low-mass part of the gamma*/Z0 spectrum.
1208
1209 * The previous process 131, g + g -> Z + b + bbar, has been removed, 
1210   for reasons of inefficient phase space generation. This implies that
1211   - all the  RK... routines are gone
1212   - all code related to ISET(ISUB)=6 is removed
1213   - variables MINT(35) and VINT(81-84) are unused
1214   - cuts CKIN(51-56) are restricted in usage
1215   - options 5 and 6 of PYOFSH are removed (with old option 7 moved to 5)
1216
1217 * Processes 147, 148, 167 and 168 have been modified to take into
1218   account changes in d*, u*, e* and nu*_e codes.
1219
1220 * New parameter and new default behaviour of the program.
1221   MSTP(9) : (D=0) inclusion of top (and fourth generation quarks) as
1222       allowed remnant flavours q' in processes that involve q -> q' + W
1223       branchings and where the matrix elements have been calculated 
1224       under the assumption that q' is massless.
1225       = 0 : no.
1226       = 1 : yes, but it is possible, as before, to switch off individual
1227           channels by the setting of MDME switches. Mass effects are
1228           taken into account, in a crude fashion, by rejecting events 
1229           where kinematics becomes inconsistent when the q' mass is 
1230           included.
1231
1232 * Many processes proceed via an s-channel resonance, described by a
1233   Breit-Wigner. In some instances this description is not really valid
1234   far away from the resonance position, e.g. because interference with
1235   other graphs should then be included. The wings of the Breit-Wigner
1236   are therefore routinely cut out in most processes, though not all.
1237   This cut has been modified and can now be set by the user.
1238   PARP(48) : (D=50.) the Breit-Wigner factor in the cross section is
1239       set to vanish for masses that deviate from the nominal one by
1240       more than PARP(48) times the nominal resonance width (i.e. the
1241       width evaluated at the nominal mass). Is used in most processes 
1242       with a single s-channel resonance, but there are some exceptions,
1243       notably gamma/Z0 and W+-. 
1244
1245 * The PYSTAT routine has been expanded with a new option 6. 
1246   A CALL PYSTAT(6) will produce a list of all subprocesses implemented
1247   in the program.
1248
1249 * New parameter.
1250   PARP(104) : (D=0.8 GeV) minimum energy above threshold for the
1251     evaluation of total, elastic and diffractive cross sections.
1252     Below this lower cut, the cross section is made to vanish.
1253
1254 -----------------------------------------------------------------------
1255
1256 THE E+E- ROUTINES
1257
1258 * The e+e- routines PYEEVT and PYONIA (formerly LUEEVT and LUONIA)
1259   have been kept in this version, but may disappear in the future.
1260   The functionality of PYEEVT is obtained with PYTHIA subprocess 1 and
1261   that of PYONIA by the decay of Upsilon in PYDECY. Some differences
1262   exist between the respective alternatives. 
1263   - The PYEEVT flavour selection and resonance shape handling is not as 
1264     good as the subprocess 1 one. 
1265   - The PYEEVT initial-state photon radiation is based on exact first 
1266     order rather than exponentiated structure functions, and is inferior 
1267     in terms of total photon energy radiated but may be better for 
1268     high-angle photons (here subprecess 19 can also be used, however).
1269   - Further examples could be given. 
1270   
1271 * Formerly, the main difference was that LUEEVT also had a number of
1272   matrix-element options, in addition to the default parton-shower one. 
1273   These options are now available also from PYTHIA subprocess 1,
1274   as follows.
1275   - MSTP(48) : (D=0) switch for the treatment of gamma*/Z0 decay for
1276         process 1 in e+e- events.
1277         = 0 : normal PYTHIA machinery.
1278         = 1 : if the decay of the Z0 is to either of the five lighter
1279             quarks, d, u, s, c or b, the special treatment of Z0
1280             decay is accessed, including matrix element options.
1281   - This option is based on the machinery of the PYEEVT and associated 
1282     routines when it comes to the description of QCD multijet structure 
1283     and the angular orientation of jets, but relies on the normal 
1284     PYEVNT machinery for everything else: cross section calculation, 
1285     initial state photon radiation, flavour composition of decays 
1286     (i.e. information on channels allowed), etc. 
1287   - The initial state has to be e+e; forward-backward asymmetries would    
1288     not come out right for quark-annihilation production of the gamma*/Z0
1289     and therefore the machinery defaults to the standard one in such
1290     cases.
1291   - You can set the behaviour for the MSTP(48) option using the normal
1292     matrix element related switches. This especially means MSTJ(101) for
1293     the selection of first- or second-order matrix elements (=1 and =2,
1294     respectively). Further selectivity is obtained with the switches
1295     and parameters MSTJ(102) (for the angular orientation part only), 
1296     MSTJ(103) (except the production threshold factor part), MSTJ(106),
1297     MSTJ(108) - MSTJ(111), PARJ(121), PARJ(122), PARJ(125) - PARJ(129).
1298     Information can be read from MSTJ(120), MSTJ(121), PARJ(150), 
1299     PARJ(152) - PARJ(156), PARJ(168), PARJ(169), PARJ(171).
1300   - The other e+e- switches and parameters should not be touched. In most 
1301     cases they are simply not accessed, since the related part is handled 
1302     by the PYEVNT machinery instead. In other cases they could give
1303     incorrect or misleading results. Beam polarization as set by
1304     PARJ(131) - PARJ(134), for instance, is only included for the 
1305     angular orientation, but is missing for the cross section information.
1306     PARJ(123) and PARJ(124) for the Z0 mass and width are set in the
1307     PYINIT call based on the input mass and calculated widths.
1308   - The cross section calculation is unaffected by the matrix element 
1309     machinery. Thus also for negative MSTJ(101) values, where only specific 
1310     jet multiplicities are generated, the PYSTAT cross section is the
1311     full one.  
1312             
1313 -----------------------------------------------------------------------
1314
1315 PARTON DISTRIBUTIONS
1316
1317 * The "structure function" expression is replaced by "parton distribution".
1318   Hence routines PYST** are renamed PYPD**.
1319
1320 * A number of old proton distributions have been removed, and newer ones
1321   inserted. The current list of distributions available in MSTP(51) is 
1322   - MSTP(51) : (D=4)
1323         = 1 : CTEQ 3L (leading order). 
1324         = 2 : CTEQ 3M (MSbar).
1325         = 3 : CTEQ 3D (DIS).
1326         = 4 : GRV 94L (leading order).
1327         = 5 : GRV 94M (MSbar).
1328         = 6 : GRV 94D (DIS).
1329         = 7 : CTEQ 5L (leading order). 
1330         = 8 : CTEQ 5M1 (MSbar; slightly updated version of CTEQ 5M).
1331         = 11: GRV 92L (leading order).
1332         = 12: EHLQ 1 (leading order).
1333         = 13: EHLQ 1 (leading order).
1334         = 14: DO 1 (leading order).
1335         = 15: DO 2 (leading order).
1336   - References: CTEQ Collaboration, H.L. Lai et al., Phys. Rev.
1337     D51 (1995) 4763, hep-ph/9903282; 
1338     M. Gluck, E. Reya and A. Vogt, Z. Phys. C67 (1995) 433 
1339   - New routines: PYCTEQ, PYGRVL, PYGRVM, PYGRVD, PYGRVV, PYGRVW, PYGRVS
1340     PYCT5L, PYCT5M and PYPDPO.
1341   - Note 1: distributions 11-15 are obsolete and should not be used for 
1342     current physics studies. They are only implemented to have some sets 
1343     in common between Pythia 5 and 6, for crosschecks. 
1344   - Note 2: distribution 16 is an undocumented toy model with all parton
1345     distributions like 1/x; see code for details.
1346
1347 * The SaS parton distributions (accessed with MSTP(55) = 5 - 12) have 
1348   been upgraded from version 1 to version 2 of the SaSgam library (as
1349   before, with routines and commonblocks renamed). This is no change for 
1350   real photons, as have been studied so far, but allows in the future 
1351   several new alternatives to extend the distributions to virtual photons.
1352   - MSTP(60) : (D=7) extension of the SaS real-photon distributions to
1353         off-shell photons, especially for the anomalous component. For 
1354         details, see G.A. Schuler and T. Sjostrand, Phys. Lett. B376 (1996) 
1355         193.
1356          = 1 : dipole dampening by integration; very time-consuming.
1357          = 2 : P_0^2 = max( Q_0^2, P^2 )
1358          = 3 : P'_0^2 = Q_0^2 + P^2.
1359          = 4 : P_{eff} that preserves momentum sum. 
1360          = 5 : P_{int} that preserves momentum and average evolution range.
1361          = 6 : P_{eff}, matched to P_0 in P2 -> Q2 limit.
1362          = 7 : P_{int}, matched to P_0 in P2 -> Q2 limit.
1363   - The PYGGAM argument list is expanded with one further input parameter IP2
1364         SUBROUTINE PYGGAM(ISET,X,Q2,P2,IP2,F2GM,XPDFGM) 
1365     where MSTP(60) is used as IP2 value in internal calls.
1366   - PYGVMD and PYGANO has VXPGA(-6:6) as new last argument. The array contains
1367     the valence part ofd the distributions at return.  
1368   - COMMON/PYINT9/VXPVMD(-6:6),VXPANL(-6:6),VXPANH(-6:6),VXPDGM(-6:6) 
1369     Purpose: to give the valence parts of the photon parton distributions
1370         (x-weighted, as usual) when the PYGGAM routine is called. Companion to 
1371         /PYINT8/, which gives the total parton distributions.
1372         VXPVMD(KFL) : valence distributions of the VMD part; matches
1373             XPVMD in /PYINT8/.
1374         VXPANL(KFL) : valence distributions of the anomalous part of light 
1375             quarks; matches XPANL in /PYINT8/.
1376         VXPANH(KFL) : valence distributions of the anomalous part of heavy
1377             quarks; matches XPANH in /PYINT8/.
1378         VXPDGM(KFL) : gives the sum of valence distributions parts; 
1379             matches XPDFGM in the PYGGAM call.
1380     Note 1: the Bethe-Heitler and direct contributions in XPBEH(KFL) and
1381         XPDIR(KFL) in /PYINT8/ are pure valence-like, andtherefore are not
1382         duplicated here.        
1383     Note 2: the sea parts of the distributions can be obtained by taking the
1384         appropriate differences between the total distributions and the
1385         valence distributions.
1386
1387 * The default set for the parton distributions of the pion has been changed:
1388   MSTP(53) (D=3) choice of pion parton distribution set; default is
1389       GRV LO (updated version).  
1390
1391 * New argument KFA for PYPDEL routine, which now also can handle muons and 
1392   taus.
1393             
1394 -----------------------------------------------------------------------
1395
1396 PARTON SHOWERS
1397
1398 * PYSHOW allows the photon emission cutoff parameter to be set 
1399   separately for quarks and leptons. The former function remains with
1400   PARJ(83), while the latter introduces the new parameter PARJ(90),
1401   with default 0.0001 GeV. Thus photon emission off leptons becomes
1402   more realistic, covering a larger part of the phase space. Since 
1403   the lepton mass is not explicitly included in the shower formalism,
1404   the emission rate is still not well reproduced (underestimated!) for
1405   lepton-photon invariant masses smaller than roughly twice the 
1406   lepton mass itself. 
1407
1408 * PYSHOW has been modified and expanded with new options related to 
1409   mass effects in the shower.
1410   - First, the emission of gluons off primary quarks in gamma*/Z0 decays
1411     has been modified. Specifically, the matrix-element correction factor 
1412     obtained for MSTJ(47)=2 or 3 (default) has been modified, to better
1413     take into account how the shower populates the phase space for 
1414     massive quarks (which is used as denominator if the corrective 
1415     weight). This increases the amount of gluon radiation, so that e.g.
1416     the amount of b bbar g three-jet events at LEP1 (within some  
1417     reasonable 3-jet region) goes up by about 5%. Light quarks are not
1418     affected.
1419   - Second, the description of g -> q qbar branchings has been expanded
1420     with several new options, in order to explore a larger range of 
1421     uncertainty in predictions. 
1422     MSTJ(42) (D=2) coherence level in shower.
1423         = 3 : in the definition of the angle in a g -> q qbar 
1424             branchings, the naive massless expression is reduced by a 
1425             factor sqrt(1 - 4 m_q^2/m_g^2), which can be motivated
1426             by a corresponding actual reduction in the p_T by mass 
1427             effects. The requirement of angular ordering then kills
1428             fewer potential g -> q qbar branchings, i.e. the rate of 
1429             such comes up. The q -> q g and g -> g g branchings are not
1430             changed from =2. This option is fully within the range of
1431             uncertainty that exists.
1432         = 4 : no angular ordering requirement conditions at all are 
1433             imposed on g -> q qbar branchings, while angular ordering
1434             still is required for q -> q g and g -> gg. This is an
1435             unrealistic extreme, and results obtained with it should 
1436             not be overstressed. However, for some studies it is of 
1437             interest. For instance, it not only gives a much higher
1438             rate of charm and bottom production in showers, but also 
1439             affects the kinematical distributions of such pairs.
1440     MSTJ(44) (D=2) alpha-strong argument in shower.
1441         = 3 : While pT^2 is used as alpha_strong argument in q -> q g 
1442             and g -> gg branchings, as in =2, instead m^2/4 is used as 
1443             argument for g -> q qbar ones. The argument is that the 
1444             soft-gluon resummation results suggesting the pT^2 scale
1445             in the former processes is not valid for the latter one,
1446             so that any multiple of the mass of the branching parton
1447             is a perfectly valid alternative. The m^2/4 ones then gives
1448             continuity with pT^2 for z=1/2. Furthermore, with this 
1449             choice,it is no longer necessary to have the requirement 
1450             of a minimum pT in branchings, else required in order to
1451             avoid having alpha_strong blow up. Therefore, in this option,
1452             that cut has been removed for g -> gg branchings. 
1453             Specifically, when combined with MSTJ(42)=4, it is possible 
1454             to reproduce the naive 1 + cos^2(theta) angular distribution
1455             of g -> gg branchings, which is not possible in any other
1456             approach. (However, as noted above, it may give too high
1457             a charm and bottom production rate in showers.)    
1458
1459 * PYSSPA is improved by new scheme for merging matrix-element and 
1460   parton-shower descriptions of initial-state radiation in the
1461   production of a single gauge-boson resonance, i.e. Z/gamma*, W,
1462   Z', W' and R. Also for other processes the possibility of changing
1463   the maximal scale of shower radiation is introduced, although here
1464   without the complete matrix-element correction machinery. The scheme
1465   is based on the study by
1466   G. Miu and T. Sjostrand, LU TP 98-30 and hep-ph/9812455;
1467   see also G. Miu, LU TP 98-9, hep-ph/9804317.
1468   As a consequence, the MSTP(68) variable takes on a new meaning. 
1469   MSTP(68) : (D=1) choice of maximum virtuality scale and matrix-element
1470       matching scheme for initial-state radiation.
1471       = 0 : maximum shower virtuality is the same as the Q2 choice for
1472           the parton distributions, see MSTP(32). (Except that the
1473           multiplicative extra factor PARP(34) is absent and instead
1474           PARP(67) can be used for this purpose. No matrix-element 
1475           correction. 
1476       = 1 : as =0 for most processes, but new scheme for processes 1, 2,
1477           141, 142 and 144, i.e. single s-channel colourless gauge boson
1478           production: gamma*/Z0, W+-, Z', W'+- and R. Here the maximum
1479           scale of shower evolution is s, the total squared energy.
1480           The nearest branching on either side of the hard scattering,
1481           of the types q -> q + g, f -> f + gamma, g -> q + qbar or
1482           gamma -> f + fbar, are corrected by the ratio of the 
1483           first-order matrix-element weight to the parton-shower one,
1484           so as to obtain an improved description. See the references
1485           above for a detailed description. Note that the improvements 
1486           apply both for incoming hadron and lepton beams.
1487       = 2 : the maximum scale for initial-state shower evolution is always 
1488           selected to be s, except for the 2 -> 2 QCD processes 11, 12, 
1489           13, 28, 53 and 68. Based on the experience in the references
1490           above, there is reason to assume that this does give an 
1491           improved qualitative description of the high-pT tail, although
1492           the quantitative agreement is currently beyond our control.
1493           No matrix-element corrections, even for the processes in =1.
1494           The QCD exception is to avoid the doublecounting issues that
1495           could easily arise here.
1496       = -1 : as =0, except that there is no requirement on uhat being
1497           negative, see point 2 below. 
1498   Note that the default behaviour of PYSSPA is changed, in several 
1499   respects. 
1500   1) For the processes listed in MSTP(68)=1, as a consequence of the new 
1501   default MSP(68) value (the old one was =0). This clearly is very
1502   important for the high-pT tail and hence implies a significant
1503   improvement here.  
1504   2) A new cut is imposed on the combination of z and Q2 values
1505   in a branching, 
1506     uhat = Q2 - shat_old * (1-z)/z = Q2 - shat_new * (1-z) < 0,
1507   where the association with the uhat variable is relevant if the branching
1508   is reinterpreted in terms of a 2 -> 2 scattering. Usually such a 
1509   requirement comes out of the kinematics, and therefore is imposed 
1510   eventually anyway. The corner of emissions that do not respect this 
1511   requirement is that where the Q2 value of the spacelike emitting
1512   parton is little changed and the z value of the branching is close 
1513   to unity. (That is, such branchings are kinematically allowed, but 
1514   since the mapping to matrix-element variables would assume the first 
1515   parton to have Q2=0, this mapping gives an unphysical uhat, and hence 
1516   no possibility to impose a matrix-element correction factor.) The 
1517   correct behaviour in this region is beyond leading-log preditivity.
1518   It is mainly important for the hardest emission, i.e. with largest Q2.
1519   The effect of this change is to reduce the total amount of emission
1520   by a non-negligible amount when no matrix-element correction is applied.
1521   (Witness results with the special option MSTP(68)=-1.) For matrix-element 
1522   corrections to be applied, this requirement must be used for the hardest 
1523   branching, and then whether it is used or not for the softer ones is 
1524   less relevant.
1525   3) The PARP(65) parameter for minimum gluon energy emitted in a 
1526   spacelike showers is modified by an extra factor roughly corresponding 
1527   to the 1/gamma factor for the boost to the hard subprocess frame. 
1528   Earlier, when a subsystem was strongly boosted, i.e. at large rapidities 
1529   in the cm frame of the collision, the PARP(65) requirement became quite 
1530   stringent on the low-energy incoming side. Therefore much radiation 
1531   could be cut out. Since this point gives changes of opposite sign to 
1532   point 2, the net result of both changes gives a small net result.
1533   4) The angular-ordering requirement is based on ordering p_T/p rather 
1534   than p_T/p_L, i.e. replacing tan(theta) by sin(theta). Earlier the 
1535   starting value tan(theta)_max = 10 could actually be violated by some 
1536   bona fide emissions for strongly boosted subsystems, where one side has 
1537   small p_L. Therefore some emissions were incorrectly removed when 
1538   MSTP(62)=3, i.e. default. 
1539   5) For incoming muon (and tau) beams the kinematical variables are
1540   better selected to represent the differences in lepton mass. 
1541
1542 * PARP(67): (D=1) new default value to better represent matching between
1543   hard-process scale and initial-statee-radiation scale in QCD processes.
1544
1545 * New variable MSTP(69), replacing some of the functionality previously
1546   provided by MSTP(68) but removed with the change in PYSSPA with 
1547   Pythia 6.119:
1548   MSTP(69) (D=0) possibility to change Q2 scale for parton distributions
1549   from the MSTP(32) choice, especially for e+e-.
1550   = 0 : use MSTP(32) scale.
1551   = 1 : in lepton-lepton collisions, the QED lepton-inside-lepton 
1552       parton distributions are evaluated with s, the full squared CM
1553       energy, as scale.
1554   = 2 : s is used as parton distribution scale also in other 
1555       processes. 
1556    
1557 -----------------------------------------------------------------------
1558
1559 UNDERLYING EVENTS
1560
1561 * The assumption of a (more or less) energy-independent pTmin or pT0
1562   lower cut-off of jet production in multiple interactions was developed
1563   in the days when parton distributions were normally assumed flat at
1564   small x (i.e. x*f_i(x,Q2) -> constant for x -> 0 at small Q2). In view
1565   of the HERA data this is no longer a valid assumption, and parton
1566   distributions have evolved to reflect this. The consequence is that
1567   the jet rate above some fixed small pTmin is increasing much faster
1568   than originally assumed. If unchecked, it leads to a much too fast
1569   increase in the multiple interaction rate, based on comparisons with
1570   data or with physics intuition. While we have no understanding of the
1571   detailed physics mechanisms, it appears sensible to introduce an
1572   explicit energy-dependence of pTmin and pT0 that closely matches
1573   this small-x dependence or, alternatively, the increase of the total
1574   cross section by the Pomeron term. Therefore the form used internally
1575   is now
1576       pTmin = PARP(81) * (ECM/PARP(89))**PARP(90), alternatively   
1577       pT0   = PARP(82) * (ECM/PARP(89))**PARP(90),
1578   with ECM the current center of mass energy. Two new parameters are
1579   introduced, PARP(89) and PARP(90), and the default values of PARP(81) 
1580   and PARP(82) are changed.
1581   PARP(81) : (D=1.9 GeV) effective mimimum transverse momentum pTmin
1582       for multiple interactions with MSTP(82) = 1, at the reference 
1583       energy scale PARP(89), with the degree of energy rescaling 
1584       given by PARP(90).
1585   PARP(82) : (D=2.1 GeV) regularization scale pT0 of the transverse 
1586       momentum spectrum for multiple interactions with MSTP(82) >= 2,
1587       at the reference energy scale PARP(89), with the degree of energy 
1588       rescaling given by PARP(90). (Current default based on MSTP(82)=4
1589       option, without any change of MSTP(2) or MSTP(33).)
1590   PARP(89) : (D=1000 GeV) reference energy scale, at which PARP(81) and
1591       PARP(82) give the pTmin and pT0 values directly. Has no physical
1592       meaning in itself, but is used for convenience only. (A form 
1593       pTmin = PARP(81) * ECM**PARP(90) would have been equally possible,
1594       but then with a less transparent meaning of PARP(81).) For studies  
1595       of the pTmin dependence at some specific energy it is convenient to
1596       choose PARP(89) equal to this energy.
1597   PARP(90) : (D=0.16) power of the energy-rescaling term of the pTmin and
1598       pT0 parameters, which are assumed proportional to ECM**PARP(90).
1599       The default value is inspired by the rise of the total cross
1600       section by the Pomeron term, s**epsilon = ECM**(2*epsilon) =
1601       ECM**(2*0.08), which is not inconsistent with the small-x behaviour.
1602       It is also reasonably consistent with the energy-dependence
1603       implied with a comparison with the UA5 multiplicity distributions
1604       at 200 and 900 GeV. PARP(90) = 0 is an allowed value, i.e. it is 
1605       possible to have energy-independent parameters.
1606    
1607 -----------------------------------------------------------------------
1608
1609 HADRONIZATION
1610
1611 * Additions and changes to the baryon production models: 
1612   see separate section below.
1613
1614 * The possibility of colour rearrangement has been introduced for 
1615   subprocess 25, W+W- pair production. 
1616   - For a description of the basic physics ideas and the scenarios 
1617     implemented, see
1618     T. Sjostrand and V.A. Khoze, Z. Phys. C62 (1994) 281.
1619     Available is also an alternative (the GH one)  loosely based on
1620     G. Gustafson and J. Hakkinen, Z. Phys. C64 (1994) 659. 
1621   - Only events where both W's decay hadronically (and not to top) 
1622     are affected. At most one reconnection is allowed per event.
1623   - The code is based on the one available since long for
1624     the PYTHIA 5.7 program as a freestanding extension. The former 
1625     version provided further information on where in an event a 
1626     reconnection occurs, but was less well integrated in the PYTHIA 
1627     framework than is the current implementation.
1628   - A new subroutine
1629     SUBROUTINE PYRECO(IW1,IW2,NSD1,NAFT1)
1630     has been added with the different scenarios. It is called from
1631     PYRESD where appropriate.
1632   - MSTP(115) : (D=0) (C) choice of colour rearrangement scenario.
1633         = 0 : no reconnection.
1634         = 1 : scenario I, reconnection inspired by a type I 
1635             superconductor, with the reconnection probability related 
1636             to overlap volume in space and time between the W+ and W- 
1637             strings. Related parameters are found in PARP(115) - 
1638             PARP(119), with PARP(117) of special interest.
1639         = 2 : scenario II, reconnection inspired by a type II 
1640             superconductor, with reconnection possible when two string 
1641             cores cross. Related parameter in PARP(115).
1642         = 3 : scenario II', as model II but with the additional 
1643             requirement that a reconnection will only occur if the 
1644             total string length is reduced by it.
1645         = 5 : the GH scenario, where the reconnection can occur that
1646             reduces the total string length (Lambda measure) most. 
1647             PARP(120) gives the fraction of such event where a 
1648             reconnection is actually made; since almost all events
1649             could allow a reconnection that would reduce the string
1650             length, PARP(120) is almost the same as the reconnection
1651             probability.
1652         = 11 : the intermediate scenario, where a reconnection is
1653             made at the "origin" of events, based on the subdivision
1654             of all radiation of a q-qbar system as coming either from 
1655             the q or the qbar. PARP(120) gives the assumed probability
1656             that a reconnection will occur. A somewhat simpleminded
1657             model, but not quite unrealistic.   
1658         = 12 : the instantaneous scenario, where a reconnection is 
1659             allowed to occur before the parton showers, and showering
1660             is performed inside the reconnected systems with maximum
1661             virtuality set by the mass of the reconnected systems.
1662             PARP(120) gives the assumed probability that a reconnection 
1663             will occur. Is completely unrealistic, but useful as an
1664             extreme example with very large effects. 
1665   - PARP(115) : (D=1.5 fm) (C) the average fragmentation time of a 
1666         string, giving the exponential suppression that a reconnection 
1667         cannot occur if strings decayed before crossing. Is implicitly 
1668         fixed by the string constant and the fragmentation function 
1669         parameters, and so a significant change is not recommended.
1670   - PARP(116) : (D=0.5 fm) (C) width of the type I string, giving the
1671         radius of the Gaussian distribution in x and y separately.
1672   - PARP(117) : (D=0.6) (C) k_I, the main free parameter in the 
1673         reconnection probability for scenario I; the probability is 
1674         given by PARP(117) times the overlap volume, up to saturation 
1675         effects.
1676   - PARP(118), PARP(119) : (D=2.5,2.0) (C) f_r and f_t, respectively,
1677         used in the Monte Carlo sampling of the phase space volume 
1678         in scenario I. There is no real reason to change these numbers. 
1679   - PARP(120) : (D=1.0) (D) (C) fraction of events in the GH, intermediate
1680         and instantaneous scenarios where a reconnection is allowed to
1681         occur. For the GH one a further suppression of the reconnection
1682         rate occurs from the requirement of reduced string length in a
1683         reconnection.    
1684   - MSTI(32) : information on whether a reconnection occured in the
1685         current event; is 0 normally but 1 in case of reconnection. 
1686   - MINT(32) : information on whether a reconnection occured in the
1687         current event; is 0 normally but 1 in case of reconnection. 
1688  
1689 * The Bose-Einstein description is expanded with several new options. 
1690     The earlier global method to ensure energy conservation
1691     (MSTJ(54)=0) has been replaced by a local one, i.e. the default 
1692     behaviour haschanged.
1693   - For a description of the basic physics ideas and the scenarios 
1694     implemented, see
1695     L. Lonnblad and T. Sjostrand, Eur. Phys. J. C2 (1998) 165.
1696   - Some new switches and parameters have been added.
1697   - MSTJ(53) (D=0) In e+e- -> W+W-, apply BE algorithm 
1698         = 0 : on all pion pairs. 
1699         = 1 : only on pairs were both pions come from the same W. 
1700         = 2 : only on pairs were the pions come from different Ws.
1701         = -1 : on all pairs except unequal pions coming from 
1702             different Ws. 
1703         = -2 : when calculating balancing shifts for pions from same W,
1704             only consider pairs from this W. 
1705   - MSTJ(54) (D=2) Alternative local energy compensation. (Notation
1706         in brackets refer to the one used in the above paper.)
1707         = 0 : global energy compensation (BE_0). 
1708         = 1 : compensate with identical pairs by negative BE 
1709             enhancement with a third of the radius (BE_3).
1710         = 2 : ditto, but with the compensation constrained to vanish
1711             at Q=0, by an addional 1-exp(-Q2*R2/4) factor (BE_32). 
1712         = -1 : compensate with pair giving the smallest invariant mass
1713             (BE_m). 
1714         = -2 : compansate with pair giving the smallest string length
1715             (BE_lambda). 
1716   - MSTJ(55) (D=0) Calculation of difference vector. 
1717         = 0 : in the lab frame. 
1718         = 1 : in the CMS of the given pair. 
1719   - MSTJ(56) (D=0) In e+e- -> W+W-, include distance between W's. 
1720         = 0 : radius is the same for all pairs. 
1721         = 1 : radius for pairs from different W's is R+deltaR_WW. 
1722             (When considering W pairs with an energy well above threshold,
1723             this should give more realistic results.)
1724   - MSTJ(57) (D=1) Penalty for shifting particles with close-by 
1725         identical neighbors in local energy compensation, MSTJ(54) < 0. 
1726         = 0 : no penalty. 
1727         = 1 : penalty. 
1728   - PARJ(95) :(R) Set to the energy imbalance after the BE algorithm, 
1729         before rescaling of momenta. 
1730   - PARJ(96) : (R) Set to the alpha needed to retain energy-momentum 
1731         conservation in each event for relevant models. 
1732              
1733 * Particle data has in part been updated to the PDG 1996 edition
1734   (Particle Data Group, R.M. Barnett et al.,Phys. Rev. D54 (1996) 1.
1735   Not yet updated are the weakly decaying charm and bottom hadron 
1736   branching ratios - here new information is added continuously, but 
1737   still without a coherent picture. (In the PDG attempts at constrained 
1738   fits, i.e. the numbers relevent for generators, 24% of D0 decays
1739   are left unexplained, 36% of D+, 82% of D_s, 93% of Lambda_c,
1740   and no attempt at all is made in the bottom sector.)
1741  
1742 * A new decay channel may be selected in PYDECY after 200 failures.
1743    
1744 -----------------------------------------------------------------------
1745
1746 BARYON PRODUCTION MODELS
1747
1748 * New advanced scheme for baryon production with the popcorn mechanism,
1749   plus some minor changes in the default older popcorn scheme.
1750   - New code written by Patrik Eden, patrik@thep.lu.se.
1751   - For a description of the new popcorn scheme, see
1752     P. Eden and G. Gustafson, Z. Phys. C75 (1997) 41.
1753
1754 * The default baryon production option, MSTJ(12)=2, is not changed in 
1755   any significant way. Advance warning is given, however, that the
1756   default may be changed in future versions, at least to MSTJ(12)=3
1757   and possibly to MSTJ(12)=5.  
1758
1759 * Three new features for the baryon production are introduced.
1760   - Improved treatment of SU(6) symmetry requirements.
1761     + If q -> B + (qq)bar is SU(6)-rejected it may now change to 
1762       q -> M + q'.
1763     + If qq -> M + qq', SU(6) symmetry is included in the weights for qq'; 
1764       qq is kept with unit probability.
1765     + As before, qq is kept and only q is reselected when qq -> B + qbar 
1766       is SU(6)-rejected.
1767     + As before, the joining qq + q -> B (when joining the two string 
1768       sides) suffers no SU(6) suppression.
1769     The arguments for this procedure are presented below. It should not 
1770     be regarded as a new model, rather a more correct implementation of 
1771     the old. However, in order to enable the user to see the effects of 
1772     the SU(6) weighting separately, both procedures are available as 
1773     different options. 
1774   - Suppression of diquark vertices with small Gamma values. This is 
1775     based on a study of the production dynamics of the three quarks that
1776     form a baryon. The main experimental consequence is a suppression 
1777     of the baryon production rate at large momentum fraction.
1778   - New flavour algorithm for baryons and popcorn (also using the 
1779     small-Gamma suppression). While the old popcorn alternative allowed
1780     at most one meson to be produced in between the baryon and the 
1781     antibaryon, the new model allows an arbitrary number. The new 
1782     flavour model makes explicit use of the popcorn suppression 
1783     exp(-2*m_T*M_T/k), where m_T is the transverse mass of the quark 
1784     creating the colour fluctuation, M_T is the total invariant 
1785     transverse mass of the popcorn meson system, and k is the string 
1786     tension constant. Thus two parameters, representing the mean 
1787     2*m_T/k for light quarks and s-quarks, respectively, governs both 
1788     diquark and popcorn meson production. A corresponding parameter is 
1789     introduced for the fragmentation of strings that contain diquarks 
1790     already from the beginning, i.e. baryon remnants.
1791
1792 * The arguments for the the new flavour SU(6) rules are as follows. 
1793   - In case of rejection due to SU(6) when q -> B + (qq)bar, one again 
1794     chooses between a diquark or a quark. If choosing diquark, a new one 
1795     is selected and tested, etc. In earlier versions of JETSET and PYTHIA, 
1796     the algorithm was instead to always produce a new diquark if the 
1797     previous one had been rejected. This leads to a slightly faster 
1798     algorithm and a better interpretation of the input parameter for the 
1799     diquark-to-quark production rate. However, the probability that a 
1800     quark will produce a baryon and a antidiquark is then flavour 
1801     independent, which is not in agreement with the model. With 
1802     JETSET 7.4 default values, this leads e.g. to an enhancement of the 
1803     Omega- relative to primary proton production with approximately a 
1804     factor 1.2.
1805   - When selecting flavours according to the popcorn model for 
1806     qq -> M + qq', the quark coming from the accepted qq is kept, and the 
1807     other member of qq', as well as the spin of qq', is chosen with weights 
1808     taking SU(6) symmetry into account. Thus the flavour of qq is not 
1809     influenced by SU(6) factors for qq', but the popcorn meson is.
1810   - When a diquark has been fitted into a symmetrical three-particle 
1811     state, it should not suffer any further SU(6) suppressions. Thus the 
1812     accompanying antidiquark should "survive" with unit probability. When 
1813     producing a quark to go with a previously produced diquark, this is 
1814     achieved by testing the configuration against the proper SU(6) factor, 
1815     and in case of rejection keep the diquark and pick a new quark, which 
1816     then is tested, etc.
1817   - There is no obvious corresponding algorithm available when a quark 
1818     from one side and a diquark from the other are joined to form the 
1819     last hadron of the string. In this case the quark is a member of a 
1820     pair, in which the antiquark already has formed a specific hadron. 
1821     Thus the quark flavour cannot be reselected. One could consider the 
1822     SU(6) rejection as a major joining failure, and restart the 
1823     fragmentation of the original string, but then the the already accepted 
1824     diquark DOES suffer extra SU(6) suppression. In the program the joining 
1825     of a quark and a diquark is always accepted.
1826
1827 * While the default behaviour of the older diquark and popcorn 
1828   scenarios is essentially unchanged, the new implementation of baryon
1829   production does imply some minor differences also here. 
1830   - The "brute force" suppression of rank 1 baryons by a paramter PARJ(19), 
1831     is no longer a special option (previously MSTJ(12)=3). For backward  
1832     compatibility, it is however not removed. Instead it is in fact always  
1833     on, but is effectively off by keeping PARJ(19)=1, its default value.
1834   - New, but fairly equivalent treatment of baryon production in closed 
1835     strings. (Calling the probability for diquark production x, the 
1836     probability for baryon production is changed from x at 1st vertex 
1837     and (1-x)x at 2nd, to 0 at 1st and x+(1-x)x at 2nd.)
1838   - New treatment of random flavour selection - the 'rndmflav' decay 
1839     product - in hadron decays. If the production of daughters fails, 
1840     a new rndmflav is now selected. Previously the same one was used
1841     until successful. (This is changed for the following reason: if 
1842     rndmflav is a diquark, at least one BB~ pair is produced, which 
1843     makes it more difficult to fulfill energy conservation, especially 
1844     if the decaying hadron is light.)
1845   - New treatment of a low-mass closed string = cluster -> 2 hadrons.
1846     (If splitting the cluster by a diquark, the old model approximation 
1847     of only one popcorn meson means that only one member of the 
1848     diquark-antidiquark pair should be allowed to split to a popcorn 
1849     meson. This is accounted for when splitting larger closed strings
1850     in PYSTRF, and when selecting rndmflav's in PYDECY. However, 
1851     it was previously not done in PYPREP.)
1852
1853 * How to use the baryon production options.
1854   - Use of the old diquark and popcorn models, MSTJ(12) = 1 and 2, is 
1855     essentially unchanged. Note, however, that PARJ(19) is available
1856     for an ad-hoc suppression of first-rank baryon production.
1857   - Use of the old popcorn model with new SU(6) weighting:
1858     + Set MSTJ(12)=3.
1859     + Increase PARJ(1) by approximately a factor 1.2 to retain about the 
1860       same effective baryon production rate as in MSTJ(12)=2.
1861     + Note: the new SU(6) weighting e.g. implies that the total 
1862       production rate of charm and bottom baryons is reduced.
1863   - Use of the old flavour model with new SU(6) treatment and modified 
1864     fragmentation function for diquark vertices (which softens baryon 
1865     spectra):
1866     + Set MSTJ(12)=4.
1867     + Increase PARJ(1) by about a factor 1.7 and PARJ(5) by about a 
1868       factor 1.2 to restore the baryon and popcorn rates of the 
1869       MSTJ(12)=2 default.
1870   - Use of the new flavour model (automatically with modified diquark 
1871     fragmentation function.)
1872     + Set MSTJ(12)=5.
1873     + Increase PARJ(1) by approximately a factor 2.
1874     + Change PARJ(18) from 1 to approx. 0.19.
1875     + Instead of PARJ(3-7), tune PARJ(8-10,18). (Here PARJ(10) is used 
1876       only in collisions having remnants of baryon beam particles.)
1877     + Note: the proposed parameter values are based on a global fit to
1878       all baryon production rates. This e.g. means that the proton rate 
1879       is lower than in the MSTJ(12)=2 option, with current data 
1880       somewhere in between. The PARJ(1) value would have to be about
1881       3 times higher in MSTJ(12)=5 than in =2 to have the same total
1882       baryon production rate (=proton+neutron), but then other baryon
1883       rates would not match at all. 
1884   - The new options MSTJ(12)=4 and =5 (and, to some extent, =3) soften 
1885     baryon spectra in such a way that PARJ(45) (delta-a for diquarks in 
1886     the Lund symmetric fragmentation function) is available for a retune. 
1887     It affects i.e. baryon-antibaryon rapidity correlations and the 
1888     baryon excess over antibaryons in quark jets.
1889
1890 * Changes in and additions to the commonblocks.
1891   MSTU(121-125) : Internal flags and counters; only MSTU(123) may be
1892       touched by user.
1893       MSTU(121) : Popcorn meson counter.
1894       MSTU(122) : Points at the proper diquark production weights, to 
1895           distinguish between ordinary popcorn and rank 0 diquark 
1896           systems. Only needed if MSTJ(12)=5.
1897       MSTU(123) : Initalization flag. If MSTU(123) is 0 in a PYKFDI call, 
1898           PYKFIN is called and MSTU(123) set to 1. Would need to be 
1899           reset by the user if flavour parameters are changed in the 
1900           middle of a run. 
1901       MSTU(124) : First parton flavour in decay call, stored to easily 
1902           find random flavour partner in a popcorn system.
1903       MSTU(125) : Maximum number of popcorn mesons allowed in decay flavour 
1904           generation. If a larger popcorn system passes the fake string 
1905           suppressions, the error KF=0 is returned and the flavour 
1906           generation for the decay is restarted.
1907   MSTU(131-140) : Store of popcorn meson flavour codes in decay algorithm.
1908       Purely internal.
1909   MSTJ(12) : (D=2) Main switch for choice of baryon production model.
1910       Suppression of rank 1 baryons by a parameter PARJ(19) is no longer 
1911       governed by the MSTJ(12) switch, but instead turned on by setting 
1912       PARJ(19)<1.
1913       Three new options are available: 
1914       = 3 : as =2, but with improved SU(6) treatment.
1915       = 4 : as =3, but also suppressing diquark vertices with low Gamma 
1916             values.
1917       = 5 : Revised popcorn model. Independent of PARJ(3-7). Depending 
1918             on PARJ(8-10). Including the same kind of suppression as =4.
1919   PARJ(8), PARJ(9) : (D=0.6,1.2 GeV^-1) The new popcorn parameters B_u 
1920       and dB = B_s - B_u. Used to suppress popcorn mesons of total 
1921       invariant mass M_T by exp(-B_q*M_T).
1922   PARJ(10) : (D=0.6 GeV^-1) Corresponding parameter for suppression of 
1923       leading rank mesons of transverse mass M_T in the fragmentation of 
1924       diquark jets.
1925   PARF(131-187) : Different diquark and popcorn weights, calculated in 
1926       PYKFIN.
1927   PARF(191) : (D=0.2 GeV) Non-constituent mass of ud_0 diquark, which has 
1928       a slight influence on the weights in the new algorithm.
1929   PARF(192) : (D=0.5) Gamma suppression parameter. The suppression factor 
1930       is 1 - PARF(192)**Gamma, with Gamma in GeV^2.
1931   PARF(193,194) : Store of some parameters used by the present popcorn 
1932       system.
1933   PARF(201-1400) : Weights for every possible popcorn meson construction 
1934        in the MSTJ(12)=5 option. Thus MSTJ(12)=5 is forbidden to be 
1935        combined with MSTJ(15)=1.
1936
1937 * In summary, all commonblock variables are completely internal, except 
1938   MSTU(123), MSTJ(12), PARJ(8-10) and PARF(191-192).
1939   - PARF(191-192) should not need to be changed.
1940   - MSTU(123) should be 0 when starting, and reset to 0 whenever changing 
1941     a switch or parameter which influences flavour weights.
1942   - With MSTJ(12)=4, PARJ(5) may need to increase.
1943   - With MSTJ(12)=5, a preliminary tune suggests 
1944     PARJ(8,9,10) = 0.6, 1.2, 0.6, PARJ(1)=0.20 and PARJ(18)=0.19.
1945
1946 * Three new subroutines are added, but are only needed for internal use.
1947   SUBROUTINE PYKFIN : Precalculates a set of diquark and popcorn weights.
1948       Called by PYKFDI if MSTU(123)=0. Sets MSTU(123) to 1.    
1949   SUBROUTINE PYNMES(KFDIQ) : If KFDIQ=0, it generates the number of 
1950       popcorn mesons and stores some relevant parameters. If KFDIQ not 0 
1951       it generates number of leading rank mesons in the fragmentation of 
1952       a diquark string with original diquark KFDIQ. Called by PYKFDI.
1953   SUBROUTINE PYDCYK(KFL1,KFL2,KFL3,KF) : Handles flavour production in 
1954       the decay of unstable particles and small string clusters. Is 
1955       essentially the same as PYKFDI, but takes into acount the effects 
1956       of string dynamics on flavour selection in the MSTJ(12)>3 options. 
1957       KFL1,KFL2,KFL3 and KF are the same as for PYKFDI. Called by PYDECY 
1958       and PYPREP.
1959
1960 * The complete list of subprogram changes is as follows.
1961   PYCOMP : Taking internal popcorn flags on diquarks into account.
1962   PYDECY, PYMASS : No longer checking diquarks for popcorn flags 
1963       before calling PYCOMP.
1964   PYKFDI : Quite differently formulated, but equivalent algorithm 
1965       introduced. Improved treatment of SU(6) symmetry requirements.
1966       New flavour algorithm, based on advanced popcorn model. 
1967   PYNMES : New function. Selects number of popcorns mesons in a 
1968       popcorn system. (Also used in the reformulated algorithm of 
1969       the old model, when it always returns 0 or 1 popcorn meson.)
1970   PYKFIN : New subroutine. Precalculates a large set of flavour 
1971       production weights from the input parameters.
1972   PYSTRF : The rank 1 baryon suppression no longer depends on any switch, 
1973       but merely on the suppression parameter. Default is no suppression.
1974       New option, suppressing diquark vertices at small Gamma, 
1975       introduced. New (but corresponding) treatment of baryon production 
1976       at first and second vertex of closed string. Suppression factors of 
1977       popcorn meson system due to its transverse mass in new flavour 
1978       algorithm introduced. Junction strings forbidden to be combined 
1979       with new popcorn options.
1980   PYINDF : The rank 1 baryon suppression no longer depends on any switch, 
1981       but merely on the suppression parameter. Default is no suppression.
1982       Warning message if trying to combine with new popcorn options.
1983       (No new options implemented.)
1984   PYDCYK : New subroutine, handles flavour selection in new popcorn model 
1985       for the case of cluster and hadron decays, where no dynamical 
1986       string variables are present. Generalized to take care of old 
1987       flavour models as well.
1988   PYDECY : Uses PYDCYK instead of PYKFDI to a large extent. Reselects 
1989       random flavour which failed.
1990   PYPREP : Uses PYDCYK instead of PYKFDI in cluster decays. This implies 
1991       a better treatment of closed string clusters, where previously both
1992       a random flavour diquark and its antidiquark partner was tested for 
1993       popcorn.
1994   PYDATA : PARJ(8-10) given default values for new flavour algorithm.
1995       Old model kept as default in MSTJ(12), PARJ(1) and PARJ(18).
1996       PARF(131-194,201-1400) and MSTU(121-140) used internally.
1997
1998 * Internally the diquark codes have been extended to store the necessary
1999   further popcorn information. As before, an initially existing diquark 
2000   has a code of the type 1000*q_a + 100*q_b + 2s+1, where q_a > q_b. 
2001   Diquarks created in the fragmentation process now have the longer code
2002   10000*q_c + 1000*q_a + 100*q_b + 2s+1, i.e. one further digit is set.
2003   Here q_c is the curtain quark, i.e. the flavour of the quark-antiquark
2004   pair that is shared between the baryon and the antibaryon, either
2005   q_a or q_b. The non-curtain quark, the other of q_a and q_b, may have 
2006   its antiquark partner in a popcorn meson. In case there are no popcorn 
2007   mesons this information is not needed, but is still set at random to be 
2008   either of q_a and q_b. The extended code is used internally in PYSTRF 
2009   and PYDECY and in some routines called by them, but is not visible in 
2010   any event listings.  
2011
2012 -----------------------------------------------------------------------
2013
2014 INTERFACES TO OTHER GENERATORS
2015
2016 * In e+e- annihilation events, a convenient classification of electroweak
2017   physics is by the number of fermions in the final state. Two fermions
2018   from Z0 decay is LEP1 physics, four fermions can come e.g. from W+W-
2019   or Z0Z0 events at LEP2, and at higher energies six fermions are produced
2020   by three-gauge-boson production or top-antitop. Often interference terms
2021   are non-negligible, requiring much more complex matrix-element expressions 
2022   than are normally provided in PYTHIA. Dedicated electroweak generators 
2023   often exist, however, and the task is therefore to interface them to
2024   the generic parton showering and hadronization machinery available in
2025   PYTHIA. In the LEP2 workshop (I.G. Knowles et al., in Physics at LEP2,
2026   CERN 96-01, eds. G.Altarelli, T. Sjostrand and F. Zwirner, p. 103) one
2027   possible strategy was outline to allow reasonably standardized 
2028   interfaces between the electroweak and the QCD generators. The LU4FRM
2029   routine was provided for the key four-fermion case. This routine is now
2030   included here, in slightly modified form, together with two new siblings 
2031   for two and six fermions. The former is trivial and included mainly for
2032   completeness, while the latter is rather more delicate.   
2033   - CALL PY2FRM(IRAD,ITAU,ICOM)
2034     Purpose: to allow a parton shower to develop and partons to hadronize
2035         from a two-fermion starting point. The initial list is supposed to
2036         be ordered such that the fermion precedes the antifermion. In 
2037         addition, an arbitrary number of photons may be included, e.g. from
2038         initial-state radiation; these will not be affected by the operation
2039         and can be put anywhere. The scale for QCD (and QED) radiation is 
2040         automatically set to be the mass of the fermion-antifermion pair. 
2041         (It is thus not suited for Bhabha scattering.)
2042     IRAD : final-state QED radiation.
2043         = 0 : no final-state photon radiation, only QCD showers.
2044         = 1 : photon radiation inside each final fermion pair, also leptons,
2045             in addition to the QCD one for quarks.  
2046     ITAU : handling of tau lepton decay (where PYTHIA does not include
2047         spin effects, although some generators provide the helicity
2048         information that would allow a more sophisticated modelling).
2049         = 0 : taus are considered stable.
2050         = 1 : taus are allowed to decay.
2051     ICOM : place where information about the event (flavours, momenta etc.)
2052         is stored at input and output.
2053         = 0 : in the HEPEVT commonblock (meaning that information is 
2054             automatically translated to PYJETS before treatment and back 
2055             afterwards).
2056         = 1 : in the PYJETS commonblock. All fermions and photons can e.g.
2057             be given with status code K(I,1)=1, flavour code in K(I,2)
2058             and five-momentum (momentum, energy, mass) in P(I,J). The
2059             V vector and remaining components in the K one are best put
2060             to zero. Also remember to set the total number of entries N.
2061   - CALL PY4FRM(ATOTSQ,A1SQ,A2SQ,ISTRAT,IRAD,ITAU,ICOM)
2062     Purpose: to allow a parton shower to develop and partons to hadronize
2063         from a four-fermion starting point. The initial list of fermions 
2064         is supposed to be ordered in the sequence fermion (1) - 
2065         antifermion (2) - fermion (3) - antifermion (4). The flavour pairs 
2066         should be arranged so that, if possible, the first two could come 
2067         from a W+ and the second two from a W-; else each pair should have 
2068         flavours consistent with a Z0. In addition, an arbitrary number of 
2069         photons may be included, e.g. from initial-state radiation; these 
2070         will not be affected by the operation and can be put anywhere. 
2071         Since the colour flow need not be unique, three real and one 
2072         integer numbers are providing further input. Once the colour 
2073         pairing is determined, the scale for QCD (and QED) radiation is 
2074         automatically set to be the mass of the fermion-antifermion pair. 
2075         (This is the relevant choice for normal fermion pair production 
2076         from resonance decay, but is not suited e.g. for 2-gamma processes 
2077         dominated by small-t propagators.) The pairing is also meaningful 
2078         for QED radiation, in the sense that a four-lepton final state is 
2079         subdivided into two radiating subsystems in the same way. Only if 
2080         the event consists of one lepton pair and one quark pair is the 
2081         information superfluous.
2082     ATOTSQ : total squared amplitude for the event, irrespective of
2083         colour flow.
2084     A1SQ : squared amplitude for the configuration with fermions 1 + 2 and
2085         3 + 4 as the two colour singlets.
2086     A2SQ : squared amplitude for the configuration with fermions 1 + 4 and
2087        3 + 2 as the two colour singlets.
2088     ISTRAT : the choice of strategy to select either of the two possible 
2089         colour configurations. Here 0 is supposed to represent a reasonable 
2090         compromize, while 1 and 2 are selected so as to give the largest 
2091         reasonable spread one could imagine.
2092         = 0 : pick configurations according to relative probabilities 
2093             A1SQ : A2SQ.
2094         = 1 : assign the interference contribution to maximize the 1 + 2 
2095             and 3 + 4 pairing of fermions.
2096         = 2 : assign the interference contribution to maximize the 1 + 4 
2097             and 3 + 2 pairing of fermions.
2098     IRAD : final-state QED radiation.
2099         = 0 : no final-state photon radiation, only QCD showers.
2100         = 1 : photon radiation inside each final fermion pair, also leptons,
2101             in addition to the QCD one for quarks.  
2102     ITAU : handling of tau lepton decay (where PYTHIA does not include
2103         spin effects, although some generators provide the helicity
2104         information that would allow a more sophisticated modelling).
2105         = 0 : taus are considered stable.
2106         = 1 : taus are allowed to decay.
2107     ICOM : place where information about the event (flavours, momenta etc.)
2108         is stored at input and output.
2109         = 0 : in the HEPEVT commonblock (meaning that information is 
2110             automatically translated to PYJETS before treatment and back 
2111             afterwards).
2112         = 1 : in the PYJETS commonblock. All fermions and photons can e.g.
2113             be given with status code K(I,1)=1, flavour code in K(I,2)
2114             and five-momentum (momentum, energy, mass) in P(I,J). The
2115             V vector and remaining components in the K one are best put
2116             to zero. Also remember to set the total number of entries N.
2117   - CALL PY6FRM(P12,P13,P21,P23,P31,P32,PTOP,IRAD,ITAU,ICOM)         
2118     Purpose: to allow a parton shower to develop and partons to hadronize
2119         from a six-fermion starting point. The initial list of fermions is 
2120         supposed to be ordered in the sequence fermion (1) - antifermion (2) - 
2121         fermion (3) - antifermion (4) - fermion (5) - antifermion (6). The 
2122         flavour pairs should be arranged so that, if possible, the first two 
2123         could come from a Z0, the middle two from a W+ and the last two from 
2124         a W-; else each pair should have flavours consistent with a Z0.
2125         Specifically, this means that in a t-tbar event, the t decay products
2126         would be found in 1 (b) and 3 and 4 (from the W+ decay) and the tbar
2127         ones in 2 (bbar) and 5 and 6 (from the W- decay). In addition, an 
2128         arbitrary number of photons may be included, e.g. from initial-state 
2129         radiation; these will not be affected by the operation and can be put 
2130         anywhere. Since the colour flow need not be unique, further input is
2131         needed to specify this. The number of possible interference 
2132         contributions being much larger than for the four-fermion case, we 
2133         have not tried to implement different strategies. Instead six 
2134         probabilities may be input for the different pairings, that the user 
2135         e.g. could pick at the six possible squared amplitudes, or according 
2136         to some more complicated scheme for how to handle the interference 
2137         terms. The treatment of cascades must be quite different for top 
2138         events and the rest. For a normal three-boson event, each fermion 
2139         pair would form one radiating system, with scale set equal to the 
2140         fermion-antifermion invariant mass. (This is the relevant choice for 
2141         normal fermion pair production from resonance decay, but is not 
2142         suited e.g. for 2-gamma processes dominated by small-t propagators.) 
2143         In the top case, on the other hand, the b (bbar) would be radiating 
2144         with a recoil taken by the W+ (W-) in such a way that the t (tbar) 
2145         mass is preserved, while the W dipoles would radiate as normal. 
2146         Therefore the user need also supply a probability for the event to 
2147         be a top one, again e.g. based on some squared amplitude.  
2148     P12, P13, P21, P23, P31, P32 : relative probabilities for the six possible
2149         pairings of fermions with antifermions. The first (second) digit tells 
2150         which antifermion the first (second) fermion is paired with, with the 
2151         third pairing given by elimination. Thus e.g. P23 means the first 
2152         fermion is paired with the second antifermion, the second fermion 
2153         with the third antifermion and the third fermion with the first 
2154         antifermion. Pairings are only possible between quarks and leptons 
2155         separately. The sum of probabilities for allowed pairings is 
2156         automatically normalized to unity. 
2157     PTOP : the probability that the configuration is a top one; a number 
2158         between 0 and one. In this case, it is important that the order 
2159         described above is respected, with the b and bbar coming first. 
2160         No colour ambiguity exists if the top interpretation is selected, 
2161         so then the P12 - P32 numbers are not used. (One could imagine 
2162         colour reconnection at later stages of the process, e.g. between 
2163         the two W's. However, we are then no longer speaking of ambiguities 
2164         related to the hard process itself but rather to the possibility of 
2165         subsquent reconnection, e.g. at the nonperturbative level. This is 
2166         an interesting topic in itself, but not the one addressed here, 
2167         where the colour assignment is used for the full cascade evolution.) 
2168     IRAD : final-state QED radiation.
2169         = 0 : no final-state photon radiation, only QCD showers.
2170         = 1 : photon radiation inside each final fermion pair, also leptons,
2171             in addition to the QCD one for quarks.  
2172     ITAU : handling of tau lepton decay (where PYTHIA does not include
2173         spin effects, although some generators provide the helicity
2174         information that would allow a more sophisticated modelling).
2175         = 0 : taus are considered stable.
2176         = 1 : taus are allowed to decay.
2177     ICOM : place where information about the event (flavours, momenta etc.)
2178         is stored at input and output.
2179         = 0 : in the HEPEVT commonblock (meaning that information is 
2180             automatically translated to PYJETS before treatment and back 
2181             afterwards).
2182         = 1 : in the PYJETS commonblock. All fermions and photons can e.g.
2183             be given with status code K(I,1)=1, flavour code in K(I,2)
2184             and five-momentum (momentum, energy, mass) in P(I,J). The
2185             V vector and remaining components in the K one are best put
2186             to zero. Also remember to set the total number of entries N.
2187
2188 * The above routines are not set up to handle QCD four-jet events, i.e.
2189   events of the types q qbar g g and q qbar q' qbar' (with q' qbar' coming 
2190   from a gluon branching). Such events are generated in normal parton
2191   showers, but not necessarily at the right rate (a problem that may be 
2192   especially interesting for massive quarks like b). Therefore one would 
2193   like to start a QCD parton shower from a given four-parton configuration.
2194   Already some time ago, a machinery was developed to handle this kind of
2195   occurences, see J. Andre and T. Sjostrand, Phys. Rev. D57 (1998) 5767. 
2196   This approach has now been adapted to Pythia 6.1, in a somewhat modified 
2197   form. The main change is that, in the original work, the colour flow was 
2198   picked in a separate first step (not discussed in the publication, since 
2199   it is part of the standard Jetset 4-parton configuration machinery), 
2200   which reduces the number of allowed q qbar g g parton-shower histories. 
2201   In the current implementation, more geared towards completely external 
2202   generators, no colour flow assumptions are made, meaning a few more 
2203   possible shower histories to pick between. Another change is that mass 
2204   effects are better respected by the z definition. In its structure, the 
2205   new code is rather different from the original Jetset 7.4 based one. 
2206   The code contains one new user routime, PY4JET, two new auxiliary ones,
2207   PY4JTW and PY4JTS, and significant additions to the PYSHOW showering
2208   routine.
2209   - CALL PY4JET(PMAX,IRAD,ICOM)
2210     Purpose: to allow a parton shower to develop and partons to hadronize
2211         from a q qbar g g or q qbar q' qbar' original configuration. The
2212         partons should be ordered exactly as indicated above, with the
2213         primary q qbar pair first and thereafter the two gluons or the
2214         secondary q' qbar' pair. (Strictly speaking, the definition of 
2215         primary and secondary fermion pair is ambiguous. In practice,
2216         however, differences in topological variables like the pair mass
2217         should make it feasible to have some sensible criterion on an event 
2218         by event basis.) Within each pair, fermion should precede antifermion. 
2219         In addition, an arbitrary number of photons may be included, e.g. from
2220         initial-state radiation; these will not be affected by the operation
2221         and can be put anywhere. The program will select a possible 
2222         parton shower history from the given parton configuration, and then 
2223         continue the shower from there on. The history selected is displayed
2224         in lines Nold+1 to Nold+6, where Nold is the N value before the 
2225         routine is called. Here the masses and energies of intermediate 
2226         partons are clearly displayed. The lines Nold+7 and Nold+8 contain
2227         the equivalent on-mass-shell parton pair from which the shower is
2228         started.  
2229     PMAX : the maximum mass scale (in GeV) from which the shower is started
2230         in those branches that are not already fixed by the matrix-element 
2231         history. If PMAX is set zero (actually below PARJ(82), the shower 
2232         cutoff scale), the shower starting scale is instead set to be equal 
2233         to the smallest mass of the virtual partons in the reconstructed
2234         shower history. A fixed PMAX can thus be used to obtain a reasonably
2235         exclusive set of four-jet events (to that PMAX scale), with little 
2236         five-jet contamination, while the PMAX=0 option gives a more
2237         inclusive interpretation, with five- or more-jet events possible.
2238         Note that the shower is based on evolution in mass, meaning the cut
2239         is really one of mass, not of pT, and that it may therefore be
2240         advantageous to set up the matrix elements cuts accordingly if one
2241         wishes to mix different event classes. This is not a requirement, 
2242         however.
2243     IRAD : final-state QED radiation.
2244         = 0 : no final-state photon radiation, only QCD showers.
2245         = 1 : photon radiation is allowed in the QCD shower (but currently
2246             a photon cannot be one of the four original partons).
2247     ICOM : place where information about the event (flavours, momenta etc.)
2248         is stored at input and output.
2249         = 0 : in the HEPEVT commonblock (meaning that information is 
2250             automatically translated to PYJETS before treatment and back 
2251             afterwards).
2252         = 1 : in the PYJETS commonblock. All fermions and photons can e.g.
2253             be given with status code K(I,1)=1, flavour code in K(I,2)
2254             and five-momentum (momentum, energy, mass) in P(I,J). The
2255             V vector and remaining components in the K one are best put
2256             to zero. Also remember to set the total number of entries N.
2257
2258 -----------------------------------------------------------------------
2259
2260 HISTOGRAMS
2261
2262 * The GBOOK package was written in 1979, at a time when HBOOK was not
2263   available in Fortran 77. It has been used since as a small and simple
2264   histogramming program. For this new version of PYTHIA the program has
2265   been updated to run together with PYTHIA in double precision. Only the
2266   one-dimensional histogram part has been retained, and subroutine names
2267   have been changed to fit PYTHIA conventions. These modified routines 
2268   are now distributed together with PYTHIA. They would not be used for
2269   final graphics, but may be handy for simple checks. 
2270
2271 * Basic principles.
2272   - There is a maximum of 1000 histograms at the disposal of the user, 
2273     numbered in the range 1 to 1000. Before a histogram can be filled, 
2274     space must be reserved (booked) for it, and histogram information 
2275     provided. 
2276   - Histogram contents are stored in a commonblock of dimension 20000,
2277     in the order they are booked. Each booked histogram requires NX+28 
2278     numbers, where NX is the number of x bins and the 28 include limits, 
2279     under/overflow and the title. If you run out of space, the program 
2280     can be recompiled with larger dimensions. 
2281   - Histograms can be manipulated with a few routines.
2282   - Histogram output is "line printer" style, i.e. no graphics. 
2283
2284 * CALL PYBOOK(ID,TITLE,NX,XL,XU)
2285   Purpose: to book a one-dimensional histogram.
2286   ID : histogram number, integer between 1 and 1000.
2287   TITLE : histogram title, at most 60 characters.
2288   NX : number of bins (in x direction) in histogram, integer between
2289       1 and 100.
2290   XL, XU : lower and upper bound, respectively, on the x range
2291       covered by the histogram.
2292
2293 * CALL PYFILL(ID,X,W)
2294   Purpose: to fill a one-dimensional histogram.
2295   ID : histogram number.
2296   X : x coordinate of point.
2297   W : weight to be added in this point.
2298
2299 * CALL PYFACT(ID,F)
2300   Purpose: to rescale the contents of a histogram.
2301   ID : histogram number.
2302   F : rescaling factor, i.e. a factor that all bin contents (including
2303       overflow etc.) are multiplied by.
2304   Remark: a typical rescaling factor could be f = 
2305       1/(bin size * number of events) = NX/(XU-XL) * 1/(number of events).
2306
2307 * CALL PYOPER(ID1,OPER,ID2,ID3,F1,F2)
2308   Purpose: this is a general purpose routine for editing one or several
2309       histograms, which all are assumed to have the same number of
2310       bins. Operations are carried out bin by bin, including overflow
2311       bins etc.
2312   OPER: gives the type of operation to be carried out, a one-character
2313       string or a CHARACTER*1 variable.
2314       = '+', '-', '*', '/': add, subract, multiply or divide the
2315           contents in ID1 and ID2 and put the result in ID3. F1 and F2, 
2316           if not 1D0, give factors by which the ID1 and ID2 bin contents 
2317           are multiplied before the indicated operation. (Division with a
2318           vanishing bin content will give 0.)
2319       = 'A', 'S', 'L': for 'S' the square root of the content in ID1
2320           is taken (result 0 for negative bin contents) and for 'L' the
2321           10-logarithm is taken (a nonpositive bin content is before that
2322           replaced by 0.8 times the smallest positive bin content).
2323           Thereafter, in all three cases, the content is multiplied by F1
2324           and added with F2, and the result is placed in ID3. Thus ID2
2325           is dummy in these cases.
2326       = 'M': intended for statistical analysis, bin-by-bin mean and
2327           standard deviation of a variable, assuming that ID1 contains
2328           accumulated weights, ID2 accumulated weight*variable and
2329           ID3 accumulated weight*variable-squared. Afterwards ID2 will
2330           contain the mean values (=ID2/ID1) and ID3 the standard
2331           deviations (=sqrt(ID3/ID1-(ID2/ID1)**2)). In the end, F1
2332           multiplies ID1 (for normalization purposes), while F2 is dummy.
2333     ID1, ID2, ID3 : histogram numbers, used as described above.
2334     F1, F2 : factors or offsets, used as described above.
2335
2336 * CALL PYHIST
2337   Purpose: to print all histograms that have been filled, and
2338       thereafter reset their bin contents to 0.
2339
2340 * CALL PYPLOT(ID)
2341   Purpose: to print out a single histogram.
2342   ID : histogram to be printed.
2343
2344 * CALL PYNULL(ID)
2345   Purpose: to reset all bin contents, including overflow etc., to 0.
2346   ID : histogram to be reset.
2347
2348 * CALL PYDUMP(MDUMP,LFN,NHI,IHI)
2349   Purpose: to dump the contents of existing histograms on an external 
2350       file, from which they could be read in to another program.
2351   MDUMP : the action to be taken.
2352        = 1 : dump histograms, each with the first line giving histogram 
2353            number and title, the second the number of x bins and lower 
2354            and upper limit, the third the total number of entries and
2355            under-, inside- and overflow, and subsequent ones the bin 
2356            contents grouped five per line. If NHI=0 all existing 
2357            histograms are dumped and IHI is dummy, else the NHI 
2358            histograms with numbers IHI(1) through IHI(NHI) are dumped.
2359        = 2 : read in histograms dumped with MDUMP=1 and book and 
2360            fill histograms according to this information. (With
2361            modest modifications this option could instead be used
2362            to write the info to HBOOK/HPLOT format, or whatever.) 
2363            NHI and IHI are dummy.
2364        = 3 : dump histogram contents in column style, where the
2365            first column contains the x values (average of respective
2366            bin) of the first histogram, and subsequent columns the
2367            histogram contents. All histograms dumped this way must
2368            have the same number of x bins, but it is not checked whether
2369            the x range is also the same. If NHI=0 all existing histograms 
2370            are dumped and IHI is dummy, else the NHI histograms with 
2371            numbers IHI(1) through IHI(NHI) are dumped. A file
2372            written this way can be read e.g. by GNUPLOT.                
2373   LFN : the file number to which the contents should be written. 
2374       You must see to it that this file is properly opened for write
2375       (since the definition of file names is machine dependent). 
2376   NHI : number of histograms to be dumped; if 0 then all existing
2377       histograms are dumped.
2378   IHI : array containing histogram numbers in the first NHI positions 
2379       for NHI nonzero.        
2380
2381 * COMMON/PYBINS/IHIST(4),INDX(1000),BIN(20000)
2382   Purpose: to contain all information on histograms.
2383   IHIST(1) : (D=1000) maximum allowed histogram number, i.e. dimension
2384       of the INDX array. 
2385   IHIST(2) : (D=20000) size of histogram storage, i.e. dimension of 
2386       the BIN array.
2387   IHIST(3) : (D=55) maximum number of lines per page assumed for 
2388       printing histograms. 18 lines are reserved for title, 
2389       bin contents and statistics, while the rest can be used for the
2390       histogram proper.
2391   IHIST(4) : internal counter for space usage in the BIN array.
2392   INDX : gives the initial address in BIN for each histogram.
2393       If this array is expanded, also IHIST(1) should be changed.
2394   BIN : gives bin contents and some further histogram information for 
2395       the booked histograms.  If this array is expanded, also IHIST(2) 
2396       should be changed.
2397
2398 -----------------------------------------------------------------------
2399
2400 MISCELLANEOUS
2401
2402 * Improved clarity of code and comments.
2403   - The contents of DO loops are indented two steps.
2404   - The header info given for each subroutine has been moved and modified.
2405   - Title page with PYLOGO has been modified.
2406
2407 * LUDBRB has been removed. The new PYROBO always requires two
2408   integer arguments to give range of action, followed by the angles
2409   and the boost vector. The integer arguments can be picked 0 to indicate 
2410   standard range (1-N). 
2411
2412 * MSTP(126) is now by default 50, giving the number of documentation
2413   lines at the beginning of the record.
2414
2415 * PYGIVE has been updated with new commonblock variables and changed
2416   array dimensions. 
2417
2418 * The random number generator PYR now works in double precision,
2419   i.e. 48 bits are set. The Marsaglia-Zaman algorithm is used, as before,
2420   with a minor extension at the initialization stage. 
2421
2422 * PYCLUS has been expanded with new options 5 and 6, which do the
2423   Durham algorithm as option 3 and 4 do the JADE one.
2424
2425 * A new subroutine PYTIME has been added to give the date and time,
2426   for use in PYLOGO and elsewhere. Since Fortran 77 does not contain 
2427   a standard way of obtaining this information, the routine is dummy,
2428   to be replaced by the user. The output is given in an integer array
2429   ITIME(6), with components year, month, day, hour, minute and second.
2430   If there should be no such information available on a system, it is
2431   acceptable to put all the numbers above to 0.
2432
2433 * Extra check in PYSCAT for low remnant energies (mainly for heavy
2434   quarks).
2435
2436 * A new function PYMRUN to allow running (Q2-dependent) masses.
2437   - PM = PYMRUN(KF,Q2) 
2438     Purpose: to give running masses of d, u, s, c and b quarks. For all other 
2439         particles, the PYMASS function is called by PYMRUN to give the normal 
2440         mass. 
2441     KF : flavour code.
2442     Q2 : the scale at which the mass is evaluated. 
2443     Note: The nominal values, valid at a reference scale 
2444         Q2ref = max((PARP(37)*nominalmass)^2 , 4*Lambda^2),
2445         are stored in PARF(91)-PARF(95). 
2446   - PARF(91) - PARF(95) : (D  = 0.0099, 0.0056, 0.199, 1.35, 4.5 GeV) default
2447         nominal masses, used to give the running masses. (Note change of b 
2448         quark mass from the 5 GeV previously used.) 
2449   - The result is that, for the d, u, s, c and b quarks, there are now 
2450     three different sets of masses  in use in the program.
2451     PMAS(KF,1) : the "on-shell" masses used to set up the kinematics of
2452         partonic state produced in an event.
2453     PARF(100+KF) : constituent masses, used in the fragmentation description,
2454         recommended not to change.
2455     PARF(90+KF) : the current algebra style masses, used as input for running
2456         masses in Higgs physics. 
2457     For diquarks, only the first two exist, and for the others only the first 
2458     one.    
2459
2460 * For the HEPEVT common, NMXHEP is 4000 rather than 2000 and real variables  
2461   are DOUBLE PRECISION, to conform with the LEP 2 workshop agreement.
2462
2463 * Some bug fixes.
2464
2465 -----------------------------------------------------------------------
2466
2467 CHANGES FROM BASELINE VERSION
2468
2469 6.100 : 4 March 1997 - baseline.
2470
2471 6.101 : 17 March 1997
2472   - PYRECO: DETER(I,J,K) -> DETER(I,J,L) to avoid problems with some
2473     compilers.
2474   - PYDUMP: bug END=180 -> END=170.
2475   - PYWIDT: calculation of beta threshold factor reorganized to avoid
2476     overflow at high energies and to remove an inconsistency.
2477   - PYSTAT: option 2 changed to allow listing of third decay product
2478     in some channels.
2479   - PYTIME: alternative timing suited for GNU LINUX libU77.
2480   - PYRAND: information on where in phase space a maximum has been
2481     violated has been reduced (MSTP(122)=0 : not at all; =1 : only 
2482     when error (i.e. not for warnings); =2 : always).  
2483
2484 6.102 : 22 April 1997
2485   - PYMASS: the special options for MSTJ(93) nonzero, used especially
2486     in the fragmentation process, have been corrected. This corrects 
2487     an error in the translation from JETSET 7.4 to PYTHIA 6.1. The 
2488     error has somewhat suppressed the amount of baryon production 
2489     relative to JETSET 7.4, but effects are not drastic.
2490   - PYMULT: the comparison XT2.LE.0.01D0*VINT(149) has been changed to
2491     0.01001 to avoid possibility of infinite loop.
2492   - PYSIGH: further check for process 145 that IA not equal to JA
2493     (purely preventive; not known to have caused any problems).
2494
2495 6.103 : 23 May 1997
2496   - PYSIGH, PYVACU, PYHGGM: some updates/corrections of the SUSY 
2497     generation. 
2498   - PYUPIN: allow external process numbers up to 500.
2499
2500 6.104 : 30 June 1997 
2501   - Three new processes for J/psi production: 106 - 108, see above,
2502     in the section on `hard processes'.
2503   - PYRESD: a major bug in the angular distribution of process 1,
2504     caused by a missing factor of 2 in the WTMAX expression. This
2505     leads to an essentially flat distribution in cos(theta).
2506
2507 6.110 : 10 October 1997
2508   - Modified code for baryon production. The default behaviour is
2509     essentially unchanged, while an advanced popcorn scheme has been
2510     added as a further option. Also some intermediate new options
2511     are implemented. The physics aspects are described above, in 
2512     the section on `baryon production models'. 
2513   - The Breit-Wigner evaluation in process 35 corrected in the same
2514     way as has already been implemented for the other Z-production
2515     processes (but apparently overlooked here).
2516   - Restore bug fix to process 145, erroneously not carried over to 
2517     version 6.104. 
2518   - In the fixed-alpha_s option MSTU(111)=0 the Lambda=PARU(117) is 
2519     set so that the first-order running alpha_s agrees with the 
2520     desired fixed alpha_s for the Q2 value used. Of no consequence 
2521     except as extra safety.
2522   - Error message if PYFILL is used with an unbooked histogram number.
2523   - Further line added to output/input for PYDUMP options 1 and 2,
2524     giving information on the total number of entries and under-, 
2525     inside- and overflow. 
2526
2527 6.111 : 27 October 1997
2528   - Forgotten values for XLO and XHI inserted in PYFINT routine.
2529   - Change of sign convention for RMSS(16) in PYAPPS routine.
2530
2531 6.112 : 30 October 1997
2532   - PYRESD has been modified to cope with the decay t -> W + b + Z
2533     (note order of decay products), by including the necessary 
2534     colour flow option and by setting angular weight according to
2535     isotropic decay of the W and Z. The program does not calculate
2536     the partial width to this potential channel. 
2537
2538 6.113 : 11 November 1997
2539   - PYEIG4 has been expanded to cover a missed ambiguity in the solution
2540     of a fourth-degree equation. This ambiguity could, for some parameter 
2541     values, give the wrong mass eigenstates in the neutralino sector. 
2542
2543 6.114 : 19 November 1997
2544   - GOTO jump into IF...ENDIF block removed from PYSTRF.
2545   - Underscore replaced by W in some PYKFIN variable names.
2546
2547 6.115 : 27 January 1998
2548   - In the intermediate scenario of colour reconnection, MSTP(115)=11,
2549     the QCD radiation has been reduced until now by an untentional 
2550     application of the colour interference machinery. This is now 
2551     solved by having MSTJ(50)=0 during the shower call.
2552   - A factor 1/SH has been missing in the width expression for
2553     t -> stop + neutralino, thus giving too large partial width.
2554   - Two errors in PYRESD corrected for the case the routine is called
2555     from outside the standard PYINIT/PYEVNT machinery, i.e. without having
2556     a subprocess number defined. The first ensures isotropic decay angles, 
2557     the second correct history pointers in K(I,3).
2558   - D-format changed to E-format in PYDUMP(3), to be consistent with
2559     GNUPLOT input conventions.
2560   - Further check on allowed histogram numbers in PYFILL, PYFACT,
2561     PYOPER, PYPLOT and PYNULL.
2562   - Removed redundant/erroneous check on MSTU(183) in PYLOGO.
2563   - MSTP(48) default changed from 2 to 0 as intended. (Should not
2564     have mattered anywhere.)
2565    
2566 6.116 : 8 July 1998
2567   - Initial-state radiation for a muon beam is now allowed (and also for a
2568     tau beam). The radiation machinery is as for an electron, with a
2569     trivial replacement of the electron mass. To distinguish the e/mu/tau 
2570     cases, the PYPDEL routine has KFA as a further argument and PYPDFU is
2571     modified accordingly.
2572   - Ten new processes, 131 - 140, for reactions involving virtual photons. 
2573     The (square root with appropriate sign of the) photon virtualities can 
2574     be set in P(1,5) and P(2,5) when PYINIT is called with the 'FIVE' option.
2575   - Two new processes, 104 and 105, for chi_c production.
2576   - PYPDFU and other routines are modified to allow virtual photons. A dummy 
2577     copy of STRUCTP (the PDFLIB routine for virtual photons) is included in 
2578     case PDFLIB is not linked.
2579   - New variable VINT(120) coincides with VINT(3) or VINT(4), depending on 
2580     which side of the event is considered. Is used to bring information on the
2581     user-defined virtuality of a photon beam to the parton distributions 
2582     of the photon. 
2583   - GRV 92L parton distribution is reinserted, for crosschecks with 
2584     Pythia 5.7. Affects PYPDPR, PYPDFU and PYINIT.
2585   - The technipi partial width to quarks corrected down by factor 3
2586     (avoiding doublecounting of colour factor).
2587   - A fudge factor PARP(146) has been introduced for the technipi partial
2588     width to a fermion pair. 
2589   - Address of the Pythia webpage is updated.
2590   - In PYSHOW an IF-test has been broken into two nested ones to avoid 
2591     testing on meaningless condition.
2592   - PYMAXI has been modified to handle the case when a user-defined 
2593     process is implemented but switched off (calculation of XSEC(0,1)).
2594   - Protection against square root of negative number in PYTHRG.
2595     
2596 6.117 : 19 August 1998
2597   - New options 11 - 25 for MSTP(14) to mix alternatives for virtual photons.
2598   - PYCLUS and PYCELL modified to ensure that N is unchanged and MSTU(3)=0
2599     when NJET is negative (to signal failure of the algorithm). 
2600
2601 6.118 : 13 September 1998
2602   - Bottom squark production is now treated separately, as for the top 
2603     squark.  However, there are more processes because bottom is in the 
2604     PDF. The new processes 281 - 296 are listed in the Hard Subprocess 
2605     section above.
2606   - Displaced vertices are now produced for resonances.  This can be
2607     particularly important for delayed decays of SUSY particles to
2608     gravitinos, e.g. ~chi0_2 -> ~gravitino + photon.
2609   - The angular distribution in chargino pair production has been
2610     reversed (i.e. that <-> uhat)  for some charge combinations.
2611   - The width for ~g -> ~squark + quark has been fixed.  The sign of 
2612     a squark mixing angle was reversed.
2613   - PYMSIN modified so that several (SUSY parameter) initializations can 
2614     be done in a single run without setting up conflicting information.
2615   - Some bugs in the technicolor decay widths have been fixed, and some
2616     new options are now available, see PARP(146) - PARP(151).
2617   - New option IMSS(5)=1 added.
2618   - New Higgs pair production processes 297-301. A few of these are
2619     already available as Z' decays, where the Z' part can be killed,
2620     but this provides a more direct implementation.
2621   - Expanded top decays to include gravitino stop and gluino stop.
2622     Added entries for virtual chargino decays of stop that might
2623     be important for light stop and light staus:
2624          ~t_1 -> ~nu_tauL tau+   b       
2625          ~t_1 -> ~tau_1+  nu_tau b       
2626     Also added entries for the neutralino:
2627          ~chi_10 -> c dbar e-      
2628                  -> d sbar nu_e    
2629     The latter two would be R-parity violating decays.
2630     The status of these decays modes is -1, and they have not
2631     been tested.
2632   - The branching ratios are zeroed out before refilling when
2633     initializing SUSY decays. 
2634
2635 6.119 : 25 September 1998
2636   - Machinery introduced to allow photon inside lepton beam.
2637     See further description above, section on hard processes. 
2638   - Extended Bose-Einstein treatment, with many new options for
2639     W pair studies, see above on hadronization. Default behaviour 
2640     changed.
2641   - PYSSPA modified so that the PARP(65) parameter for minimum gluon 
2642     energy emitted in spacelike showers is modifed by an extra factor
2643     roughly corresponding to the 1/gamma factor for the boost to the
2644     hard subprocess frame. Earlier, when a subsystem was strongly 
2645     boosted, i.e. at large rapidities in the cm frame of the collision,
2646     the PARP(65) requirement became quite stringent on the low-energy
2647     incoming side. Therefore much radiation could be cut out.
2648   - PYSSPA modified so that the angular ordering requirement is based 
2649     on ordering p_T/p rather than p_T/p_L, i.e. replacing tan(theta)
2650     by sin(theta). Earlier the starting value tan(theta)_max = 10
2651     could actually be violated by some bona fide emissions for strongly
2652     boosted subsystems, where one side has small p_L. Therefore some
2653     emissions were incorrectly removed when MSTP(62)=3, i.e. default.
2654   - PYSSPA now sets relevant mass for QED emission to be mu or tau one
2655     rather than e one for such incoming beams. (For a collider between
2656     two different lepton species, the more massive one is used as
2657     reference.)
2658   - PYSHOW modified, so that photon emission off a lepton is governed
2659     by the PARJ(90) parameter rather than PARJ(83) (see PARTON SHOWERS).
2660   - PYUPDA corrected for bug in calculation of phase space available
2661     in decay (generated unnecessary warnings).
2662   - PYRESD modified to avoid calculation of undefined four-products
2663     when called for an odd resonance (i.e. one not part of the 
2664     standard PYTHIA machinery, e.g. filled with PY1ENT). 
2665   - In pair production of heavy flavour (top) in processes 81,82, 84 
2666     and 85, earlier only one of the masses was used in the matrix element,
2667     under the assumption that the two were identical. Since we do not
2668     have expressions involving the two separately, we now use an average
2669     value constructed so that the beta kinematics factor is the same
2670     for both having the average as for each having its correct value.
2671   - Move technicolour parameter PARP(151) to PARP(140) to avoid clash.
2672   - Effects of secondary widths included if leptoquark decays to top
2673     (or fourth-generation fermions).
2674
2675 6.120 : 1 October 1998
2676   - The pTmin and pT0 cutoff parameters of the multiple interactions
2677     scenario(s) are now made explicitly energy-dependent (see 
2678     MISCELLANEOUS).
2679   - MINT(45), MINT(46) set correctly to allow photon radiation off a
2680     muon beam. Also some other minor bugs corrected for muon beams.
2681     Note, however, that the MSTP(12)=1 option to obtain e.g. electrons
2682     inside photons inside electrons does not work for muons. 
2683   - PYSSPA modified so lower Q2 cutoff for QED radiation off lepton 
2684     is always at least twice the mass-squared, in addition to the 
2685     cutoff provided by PARP(68).
2686   - W2 limits in CKIN(39) and CKIN(40) not checked if process 10 is 
2687     called for two lepton beams.
2688   - Labels cleaned up. 
2689  
2690 6.121 : 15 October 1998
2691   - New routines PY2FRM, PY4FRM and PY6FRM added as generic interfaces 
2692     to two-, four- and six-fermion generators, see MISCELLANEOUS.
2693   - The MSTP(14) switch has been expanded so that MSTP(14)=20 and =25
2694     works also for gamma-hadron, not only for gamma-gamma. These two 
2695     values would therefore be the two main alternatives for users. 
2696     The default has been changed to MSTP(14)=20.
2697   - The MSTP(32) parameter for choice of Q2 scale has been expanded 
2698     with new options intended for virtual incoming photons.
2699   - New function PYMRUN(KF,Q2) gives running (MSbar) mass of d, u, s, 
2700     c and b quarks. For all other KF, the PYMASS function is called by
2701     PYMRUN to give the normal mass. PYWIDT and PYSIGH has been modified 
2702     for Higgs (and some technicolour) widths and production processes
2703     to call PYMRUN rather than to implement the running inline. The
2704     code for the running is identical, so the difference is that now
2705     the PMAS(KC,1) masses can be set to the "on-shell" values expected
2706     rather than the MSbar ones. The nominal b quark mass has been reduced 
2707     from 5 to 4.5 GeV, affecting some Higgs branching ratios. 
2708     The technipi rate to leptons has been somewhat changed.
2709     See also MISCELLANEOUS.
2710   - Correct minor bug in partial width of top to gravitino + stop. 
2711   - Reimplement PARP(146)-(148) in code (had been lost).
2712   - Minor correction in the initialization printout. 
2713  
2714 6.122 : 4 January 1999
2715   - New matrix-element correction scheme for initial-state radiation,
2716     especially relevant for the production of a single s-channel
2717     resonance. This allows much better description e.g. of the pT 
2718     properties of W and Z produced in hadron colliders. See
2719     MISCELLANEOUS for further details.
2720   - Change in PYRESD so that, when Z' or W' decays to (one or two) 
2721     top quarks, these are allowed to decay isotropically. (Previously
2722     the matrix element for Z' -> W+ W- -> 4 fermions was erroneously
2723     called.)
2724   - Minor change in PYSHOW to catch one case where K(I,1) values 
2725     can become incorrectly set if the routine is called for a lepton
2726     pair of very low mass (roughly below 1 GeV). 
2727   - The lepton-inside-lepton parton distribution is changed. Previously
2728     f_e^e(x) was normal for x < 1 - 10^-4, scaled up for 
2729     1 - 10^-4 < x < 1 - 10^-6 and 0 for x > 1 - 10^-6, where the
2730     rescaling was arranged so as to give the correct integral of 
2731     f_e^e(x) from 0 to 1. Now the border at 1 - 10^-4 has been moved
2732     to 1 - 10^-7 and the one at 1 - 10^-6 to 1 - 10^-10. This way any 
2733     irregularities in the line shape has been pushed into a much narrower 
2734     region; of some interest e.g. for a muon collider.
2735   - Angular distribution included in decay of W in process 36,
2736     gamma + f -> W + f', by analogy with process 31.
2737   - DATA PARU split in two statements to avoid the 19-continuation-lines
2738     limit.
2739   - Extra safety check in PYREMN to avoid division by zero if
2740     chi = 0 or 1.
2741   - Matrix-element code MDME(IDC,2)=32 restored for h0, H0, A0, H+- ->
2742     q qbar (set 0 in recent versions). This code is irrelevant when 
2743     resonance decays is performed in PYRESD, as is almost always the
2744     case.
2745   - PYMSIN modified so that MWID(KC) and MDCY(KC,1) values are saved
2746     and restored for the lightest supersymmetric particle. Is relevent
2747     where a single run contains several PYINIT calls for different
2748     SUSY parameter sets, and hence different LSP's: it switches back on 
2749     the decays of a particle that was LSP but no longer is it.
2750   - PYLOGO, PYTIME and PYHIST slightly modified for year 
2751     2000-compatible output.
2752   - New option MSTP(39)=5, where the Q2 scale of the gg, qqbar -> QQbarH 
2753     processes is set equal to the squared nominal Higgs mass 
2754     (cf. MSTP(39)=3 is the actual Higgs mass, i.e. fluctuating between
2755     events).
2756   - Introduce a line
2757           IMPLICIT INTEGER(I-N)
2758     in routines. This helps avoid a bug in the SGI Fortran compiler.
2759
2760 6.123 : 2 February 1999
2761   - New process machinery for doubly charged Higgs production in a 
2762     left-right-symmetric scenario. Includes new particles and new hard
2763     subprocesses; see these subsections.
2764   - Introduce missing shat factor in cross section for process 140.
2765   - Correct logic of photon virtuality choice in processes 131 - 136,
2766     which gave erroneous results for the direct*resolved cases of
2767     gamma*gamma* events.
2768   - Explicit DOUBLE PRECISION declaration for EXTERNAL functions and
2769     some DATA statements moved after all declarations to avoid problems
2770     on some compilers.
2771   - The pT^2 fluctuation margin allowed for independent fragmentation 
2772     in PYTEST increased. 
2773  
2774 6.124 : 7 February 1999
2775   - The effects of longitudinal resolved photons can be approximated 
2776     by a multiplicative factor to the transverse resolved cross sections,
2777     see PARP(165) in the hard processes section.
2778   - Possibility to choose between e -> gamma splitting variable 
2779     being energy fraction x or lightcome fraction y, see MSTP(16).
2780   - Cross sections for direct photon processes 137-140 corrected by a
2781     factor shat/lambda, usually very close to unity, to better describe
2782     phase space relations. 
2783   - A few bug corrections in the new popcorn scenario (see section 
2784     above on baryon production models). Especially, one bug also came 
2785     to affect the default baryon production scenario, and could in
2786     some cases result in charge and baryon number nonconservation in
2787     the beam remnant fragmentation process (PYREMN).
2788   - PYKFIN extensively rewritten. Mostly cosmetics, but also 
2789     1) For MSTJ(12)=5, a factor 2 was misplaced for ud_1 and uu_1 
2790     diquark production in the process (rank 0 qq) -> ... M + B + ...
2791     2) In the old algorithm the diquark SU(6) survival factor was not
2792     considered when generating the final diquark of a popcorn
2793     system. In Pythia 6.110, this factor was introduced for the new
2794     options MSTJ(12)>2, but unintentionally also for MSTJ(12)=2. For
2795     backward compatibility, the diquark SU(6) survival factors are now
2796     set to 1 if MSTJ(12)<3.
2797   - IN PYRAND the VINT(25) = x_T^2 calculation was incorrect for a 
2798     user-defined process; normalization now changed from VINT(1) to
2799     VINT(2). Will have given too high a starting x_T for multiple
2800     interactions.
2801   - The well-known but harmless rho0 -> eta gamma and a_2 -> eta' pi
2802     possibilities of looping in PYDECY no longer cause a warning 
2803     message.
2804   - In process 23 the cross section in PYSIGH is explicitly ensured to
2805     be non-negative. This is likely a problem of the far-out wings of 
2806     the Breit-Wigners, which the cross section is not set up to handle. 
2807  
2808 6.125 : 21 February 1999
2809   - PYSTRF corrected for a bug in the choice of the string region 
2810     which defines the longitudinal directions of the final two hadrons.
2811     In principle the bug is severe, but in practice its consequences 
2812     are limited, since in many events the final string region is 
2813     uniquely defined so that the choice is irrelevent, and since,
2814     even when there is a choice, the procedure would work whichever 
2815     of the two possible regions is selected.
2816   - PYSSPA treatment of QED showers corrected, in three respect. First,
2817     lower x cutoff (XEE) changed from 1D-7 to 1D-10 to match changes
2818     in lepton-in-lepton distributions of 6.122. Second, the 
2819     matrix-element matching can made also for QED processes.
2820     Third, a scattered lepton does not (occasionally) get K(I,1)=3.
2821   - New (default) option for lower parton-shower cutoff (and `primordial 
2822     kT') in resolved photons, see MSTP(66) above.
2823   - New parameter and default behaviour for multiple interactions in the 
2824     VMD component of virtual photons, see MSTP(84) above.
2825   - For non-QCD processes, a photon is now assumed unresolved when 
2826     MSTP(14)=10, 20 or 25. (In principle, both the resolved and direct
2827     possibilities ought to be explored, but this mixing is not currently 
2828     implemented, so picking direct at least will explore one of the two 
2829     main alternatives rather than none.) Also a minor change in PYMAXI 
2830     to correct the calculation of number of points in y* for a photon 
2831     beam.
2832   - New option MSTP(32)=10 : Q2 scale is equal to CM energy. No special
2833     reason except as extreme contrast (or not so extreme, for many e+e-
2834     processes).
2835
2836 6.126 : 26 March 1999
2837   - The simulation of the production and decays of technicolor particles
2838     has been substantially upgraded. The processes 149, 191, 192, and 193 
2839     are to be considered obsolete, and are temporarily retained to allow 
2840     cross checking with the new processes. Process 194 has been changed 
2841     to more accurately represent the mixing between the photon, Z, 
2842     techni_rho0, and techni_omega particles in the Drell-Yan process.  
2843     Process 195 is the analogous process including W and techni_rho+/- 
2844     mixing. Processes 361 - 377 are completely new. For a listing of
2845     processes and parameters, see the description in the Hard Processes
2846     section.   
2847   - The possibility of flavor--dependent Z0' couplings has been considered.  
2848     The previous set of parameters PARU(121)--PARU(128) now affect only 
2849     the first generation of fermions.  As before, these parameters represent 
2850     the V and A couplings for the d-quark, u-quark, electron, and nu_e, 
2851     respectively.  The parameters PARJ(180)--PARJ(187) and 
2852     PARJ(188)--PARJ(195) represent the V and A couplings of the s-quark, 
2853     c-quark, muon, nu_mu and b-quark, t-quark, tau, and nu_tau, 
2854     respectively.  The default value for all parameters are the standard 
2855     model values.
2856   - PYMSIN : improve zeroing of branching ratios when several parameter
2857     sets are considered in the same run.
2858
2859 6.127 : 18 May 1999
2860   - In process 226, for chargino pair production, the sign in the
2861     u quark inteference term in the cross section is changed.  
2862   - In PYRESD, the HGZ array of relative Z/gamma weights in processes
2863     15, 19, 30 and 35 was not always stored and read out with the same
2864     index, leading to a potentially incorrect angular distribution in 
2865     the Z decay, specifically concerning forward-backward asymmetries. 
2866   - In process 25, W pair production, the contribution from Z exchange
2867     to the cross section is now evaluated with a fixed width for the Z 
2868     in the propagator, in PYSIGH. That is, the GMMZ = mass * width
2869     in the denominator of the progagator is used with the nominal mass
2870     and width of the Z. Previously the actual mass and the running width
2871     were used, which gave rise to divering cross sections, by imperfect
2872     gauge cancellation, at large energies. 
2873   - In process 140, for longitudinal photon interactions, the cross section
2874     corrected for an erroneous factor shat too much.
2875   - In photon physics, the setting of MINT(15), MINT(16), VINT(307) and
2876     VINT(308) have been corrected for some cases, affecting PYRAND,
2877     PYGAGA, PYMAXI and PYINPR.
2878   - The normalizing cross section used for multiple interactions in
2879     photon collisions is scaled by a factor m^4/(m^2+Q^2)^2 for virtual
2880     photons, rather than only the square root of that.
2881   - New variable MSTP(69), replacing some of the functionality previously
2882     provided by MSTP(68) but removed with the change in PYSSPA with 
2883     Pythia 6.119:
2884     MSTP(69) : (D=0) possibility to change Q2 scale for parton distributions
2885     from the MSTP(32) choice, especially for e+e-.
2886     = 0 : use MSTP(32) scale.
2887     = 1 : in lepton-lepton collisions, the QED lepton-inside-lepton 
2888         parton distributions are evaluated with s, the full squared CM
2889         energy, as scale.
2890     = 2 : s is used as parton distribution scale also in other 
2891         processes. 
2892   - Insert WID2=1 in a few more places in PYWIDT, to avoid it being 
2893     undefined.
2894   - THE, PHI, CHI -> THEZ, PHIZ, CHIZ in the special e+e- -> Z option
2895     of PYRESD, to avoid a name clash.
2896   - Insert extra demand when storing THE2T in PYSSPA, for consistency
2897     (to avoid storing an undefined variable).
2898   - Replace SQMW*PMAS(24,2)**2 by GMMW**2 in PYSIGH. 
2899   - Remove unused SR2 in PYSIGH.
2900   - Further examples of PYTIME solutions on some machines have been added.
2901
2902 6.128 : 3 June 1999
2903   - Introduce new options for MINT(47):
2904     = 6 : parton distribution is peaked at x=1 for target and no 
2905         distribution at all for beam.
2906     = 7 : parton distribution is peaked at x=1 for beam and no 
2907         distribution at all for target.
2908     This prompts modifications in several routines, especially
2909     PYSIGH and PYKMAP, with modified checks and phase space factors,
2910     and also e.g. the possibility of having a 1/(1-tau) term in
2911     the selecion procedure of e-gamma collisions.  
2912   - Put MINT(15 or 16) = 22 in PYGAGA for a photon whenever MSTP(14)=0.
2913   - Modify the cross section for process 35 to depend on the 
2914     virtuality of the incoming photon in e gamma* -> e Z0.
2915     (Exact form ambiguous, but hopefully sensible choice.)
2916   - In PYOFSH the modification of the allowed resonace mass range by
2917     CKIN values is modified when the other particle is not a resonance.
2918
2919 6.129 : 9 July 1999
2920   - Correct severe bug in colour reconnection with PYRECO, whereby the 
2921     W+ and W- decay vertices were swapped if the resonaces were given 
2922     in the order W- W+. This is the case when PYINIT is called with the 
2923     beam arguments 'e-','e+' rather than 'e+','e-'. In scenarios I, II 
2924     and II', the rearrangement rate is then overestimated.
2925   - Allow to switch on colour reconnection also for subprocess 22, Z0Z0 
2926     pair production, analogously with rearrangement in W+W- events.
2927     Note, however, that the Z0 decay vertex position is calculated 
2928     without any regard to the gamma* component of the cross section. 
2929     Thus, the description in scenarios I, II and II' would not be 
2930     sensible e.g. for a light-mass gamma*gamma* pair.
2931   - The width of the A0 higgs particle to a fermion pair is corrected
2932     to be like beta (=velocity), rather than like beta**3.
2933   - Default decay status for Higgs modes to supersymmetric particles
2934     changed from -1 to 1 in MDME(IDC,1) in PYDATA.
2935   - New options to take into account the effects of resolved longitudinal
2936     photons, see above MSTP(17), PARP(165), PARP(167) and PARP(168)
2937     (PYSIGH routine).
2938   - When calling virtual-photon PDF's via PDFLIB, ensure that P^2 < Q^2.
2939     For GRS, also ensure the specific cuts in that parameterization,
2940     specifically P^2 < Q^2/5.   
2941
2942 6.130 : 6 September 1999
2943   - pi0 decay was unintentionally switched off by default from version
2944     6.126 onwards, but is now again allowed to decay.
2945   - A major bug has been corrected in the PYBOEI routine for Bose-Einstein
2946     corrections. It does not affect the BE_0, BE_3 and BE_32 (default) 
2947     options, but only the BE_m and BE_lambda alternatives, MSTJ(54)=-1
2948     and -2 (when MSTJ(57)=1, which is default). The weight used to 
2949     define the most likely particles to carry the energy/momentum
2950     compensation of BE pairs contained a sign error in an exponential,
2951     which meant that not always the intended particles were selected.
2952     This affects several of the distributions obtained with these 
2953     algorithms in the past. The predicted average W mass shift is among
2954     the quantities changed, but stays of the same order.  
2955   - PYSHOW has been modified and expanded with new options, see the
2956     PARTON SHOWERS section above for further details. 
2957     First, the emission of gluons off primary quarks in gamma*/Z0 decays
2958     has been modified. This increases the amount of gluon radiation off
2959     heavier quarks like b's (by about 5% at LEP1), while light quarks 
2960     are not affected.
2961     Second, the description of g -> q qbar branchings has been expanded
2962     with several new options, MSTJ(42)=3 and 4 and MSTJ(44)=3, in order 
2963     to explore a larger range of uncertainty in predictions. 
2964   - PY6FRM has been improved, so that if an event is classified as a
2965     ttbar one, the t pair is allowed to radiate gluons before the top 
2966     decay. Radiation off the b's and in the W decays is there as earlier.
2967
2968 6.131 : 13 September 1999
2969   - New routine PY4JET introduced (with auxiliary routines PY4JTW and
2970     PY4JTS, and significant additions to PYSHOW) to provide interface 
2971     from a four-jet QCD generator to a parton-shower evolution. See
2972     further description in section on INTERFACES TO OTHER GENERATORS. 
2973
2974 6.132 : 23 September 1999
2975   - Default (pseudo)rapidity limits in CKIN(9)-CKIN(16) changed from 
2976     +-10 to +-40, since former still can imply an unwanted pTmin cut.
2977   - Also PYKLIM modified to avoid erroneous rejections at extreme
2978     (pseudo)rapidities.
2979   - When MINT(15)=1 is set before a PYWIDT call, the original value
2980     is restored afterwards.
2981   - Incoming/outgoing lepton mass included in cross section for 
2982     process 35.
2983
2984 6.133 : 29 September 1999
2985   - Correct the calculation of that and uhat for process 35 when the
2986     incoming photon is virtual, so that masses are assigned assuming
2987     that incoming parton number 2 is the photon (for the internal 
2988     numbering).
2989   - Introduce a missing factor of 1/2 in the cross section for processes
2990     351 and 352, H++/H-- production, and change the rules for when
2991     t^ and u^ contributions should be symmetrized so it is only done
2992     for identical leptons.  
2993
2994 6.134 : 10 October 1999
2995   - New internally available parton distribution parameterizations for 
2996     the proton, CTEQ 5L and CTEQ 5M1. These are obtained with 
2997     MSTP(51)=7 and =8, respectively. Internally: changes in PYPDPR and 
2998     PYINIT, and additions of new routines PYCT5L and PYCT5M. All is based
2999     on code written by Jon Pumplin, with minor modifications to fit the
3000     PYTHIA framework.
3001
3002 6.135 : 3 November 1999
3003   - Major modifications in routine PYPREP, intended e.g. to improve
3004     modelling of charm/bottom production from small-mass parton systems,
3005     "clusters". It is now possible to select between a few alternative
3006     descriptions of how energy/momentum is shuffled when a cluster 
3007     collapses to a single particle, and to have anisotropic decay when
3008     a cluster gives two particles. See E. Norrbin and T. Sjostrand,
3009     Phys. Lett. B442 (1998) 407 and in preparation.  
3010     MSTJ(16) : (D=2) mode of cluster treatment.
3011         = 0 : old scheme. Cluster decays (to two hadrons) are isotropic.
3012             In cluster&n