Possibility to have different binaries in the same tree introduced
[u/mrichter/AliRoot.git] / PYTHIA / doc / jetset74.lis
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1996/03/08 17:32:28  mclareni
6 * jetset74
7 *
8 *
9 * This directory was created from jetset74.car patch update
10
11
12
13                                                  23 August 1995
14  
15  
16                  Updates to
17  
18           PYTHIA 5.7 and JETSET 7.4
19              Physics and Manual
20  
21              Torbjorn Sjostrand
22      Department of theoretical physics 2
23              University of Lund
24                Solvegatan 14A
25                S-223 62 Lund
26                    Sweden
27  
28 -----------------------------------------------------------------
29  
30 1) Introduction.
31  
32 Since the PYTHIA/JETSET programs are being updated frequently,
33 it is also important to keep the documentation up to date.
34 The big manual that appeared as CERN-TH.7112/93 described the
35 status as of 15 December 1993. The LaTeX file with this manual 
36 will be updated, though less frequently than the programs. 
37 Further, it is not very economical to have to get a new copy 
38 of a long manual just to see if any interesting new features 
39 have been added recently. Therefore I have here collected the 
40 main changes that have taken place since the beginning of 1994. 
41 The changes are indexed by sub-subversion number and date. This 
42 file will be updated as regularly as the programs themselves. 
43 It is not intended to be as complete as the ordinary manual, 
44 but should be sufficient for the intended purpose.
45  
46 -----------------------------------------------------------------
47  
48 2) Changes in JETSET 7.4
49  
50 -----------
51  
52 00, 13 December 1993: baseline version.
53  
54 -----------
55  
56 01, 11 February 1994:
57  
58 LUZDIS has been changed to protect against an overflow in an exponent
59 (harmless physicswise, but enough to crash the program on some 
60 machines).
61  
62 -----------
63  
64 02, 7 April 1994:
65  
66 A possibility has been introduced into LUSHOW to suppress either hard
67 or soft radiation, in the connection of radiation off rapidly decaying
68 objects. The algorithm used is not exact, but still gives some 
69 impression of potential effects. The switch ought to have appeared at
70 the end of the current list of shower switches, but because of lack of 
71 space it appears immediately before.
72  
73 MSTJ(40) : (D=0) possibility to suppress the branching probability for 
74 a branching q -> q + g (or q -> q + gamma) of a quark produced in the
75 decay of an unstable particle with width Gamma, where this width has to
76 be specified by the user in PARJ(89). Can be changed for each new
77 LUSHOW call.  
78 = 0 : no suppression, i.e. the standard parton-shower machinery.
79 = 1 : suppress radiation by a factor chi(omega) = 
80     Gamma**2 / (Gamma**2 + omega**2), where omega is the energy of the
81     gluon (or photon) in the rest frame of the radiating dipole.
82     Essentially this means that hard radiation with omega > Gamma
83     is removed. 
84 = 2 : suppress radiation by a factor 1 - chi(omega) = 
85     omega**2 / (Gamma**2 + omega**2), where omega is the energy of the
86     gluon (or photon) in the rest frame of the radiating dipole.
87     Essentially this means that soft radiation with omega < Gamma
88     is removed. 
89  
90 PARJ(89) : (D=0. GeV) the width of the unstable particle studied for the
91     MSTJ(40) > 0 options; to be set by the user (separately for each
92     LUSHOW call, if need be). 
93  
94 -----
95  
96 A generic interface has been included to an external tau decay library.
97 This should allow the handling of tau polarization, which is not done
98 by JETSET. To use this facility you have to set the switch MSTJ(28),
99 include your own interface routine LUTAUD and see to it that the dummy
100 routine LUTAUD in JETSET is not linked. The dummy routine is there 
101 only to avoid unresolved external references when no user-supplied 
102 interface is linked.
103  
104 MSTJ(28) : (D=0) call to an external tau decay library.
105 = 0 : not done, i.e. the internal LUDECY treatment is used.
106 = 1 : done whenever the tau mother particle species can be identified,
107     else the internal LUDECY treatment is used. Normally the mother 
108     particle should always be identified, but it is possible for 
109     a user to remove event history information or to add extra taus 
110     directly to the event record, and then the mother is not known.
111 = 2 : always done. 
112  
113 CALL LUTAUD(ITAU,IORIG,KFORIG,NDECAY)
114 Purpose: to act as an interface between the generic decay routine
115     LUDECY and a user-supplied tau decay library. The latter library 
116     would normally know how to handle polarized taus, given the tau
117     polarization, so one task of the interface routine is to construct
118     the tau polarization/helicity from the information available.
119     Input to the routine is provided in the first three arguments, 
120     while the last argument and some event record information have 
121     to be set before return.
122 ITAU : line number in the event record where the tau is stored. 
123     The four-momentum of this tau has first been boosted back to the
124     rest frame of the decaying mother and thereafter rotated to move 
125     out along the +z axis. This choice of frame should help the 
126     calculation of the helicity configuration. After the LUTAUD call
127     the tau (and its decay products) will be rotated and boosted back.
128     However, seemingly, the event record does not conserve momentum
129     at this intermediate stage.  
130 IORIG : line number where the mother particle to the tau is stored.
131     Is 0 if the mother is not stored. This does not have to mean the
132     mother is unknown, e.g. in semileptonic B decay the mother is a
133     W+-, and its momentum can be obtained by adding the tau and
134     nu_tau momentum, but there is no line in the event record. 
135     When several copies of the mother is stored (e.g. one in the 
136     documentation section of the event record and one in the main
137     section), IORIG points to the last. If a branchings like 
138     tau -> tau + gamma occurs, the 'grandmother' is given, i.e. the
139     mother of the direct tau before branching.
140 KFORIG : flavour code for the mother particle. Is 0 if the mother
141     is unknown. The mother would typically be a resonance such as 
142     Z0 (23), W+- (+-24), H0 (25), or H+- (+-37).
143     Often the helicity choice would be clear just by the knowledge
144     of this mother species, e.g. W+- vs. H+-. However, sometimes
145     further complications may exist. For instance, the KF code 23
146     represents a mixture of gamma* and Z0; a knowledge of the mother
147     mass (in P(IORIG,5)) would here be required to make the choice
148     of helicities. Further, a W and Z may either be (predominantly)
149     transverse or longitudinal, depending on the production process      
150     under study.
151 NDECAY : the number of decay products of the tau; to be given by the
152     user. You must also store the KF flavour codes of those decay
153     products in the positions K(I,2), N+1 <= I <= N+NDECAY, of the 
154     event record. The corresponding five-momentum (momentum, energy
155     and mass) should be stored in the associated P(I,J) positions,
156     1 <= J <= 5. The four-momenta are expected to add up to the
157     four-momentum of the tau in position ITAU. You should not change
158     the N value or any of the other K or V values (neither for the 
159     tau nor for its decay products) since this is automatically done 
160     in LUDECY.  
161  
162 -----
163  
164 In a few places, a dot has been moved from the end of one line to the
165 beginning of the next continuation line, or the other way around, 
166 to keep together tokens such as .EQ. or .AND., since some debuggers
167 may otherwise complain. 
168  
169 A source of (harmless) division by zero in LUSHOW has been removed.
170  
171 -----------
172  
173 03, 15 July 1994:
174  
175 The LUBOEI routine has been changed to avoid an unintentional gap
176 in the limits of the very first bin. 
177  
178 Further, leptons and photons which are unrelated to the system 
179 feeling the Bose-Einstein effects do not have their energies and
180 momenta changed in the global rescaling step. (Example: W+W- events,
181 where one W decays leptonically; before these lepton momenta could be 
182 slightly changed, but now not.)
183  
184 -----
185  
186 The option LUEDIT(16) (used e.g. from PYEVNT) has been improved with
187 a more extensive search for missing daughter pointers.
188  
189 -----
190  
191 The KLU(I,16) procedure for finding rank has been rewritten to work
192 in the current JETSET version, which it did not before. However, note
193 that it will only work for MSTU(16)=2. As a general comment, the
194 options 14 - 17 of KLU were written at a time when possible event 
195 histories were less complex, and can not be guaranteed always to work 
196 today.
197  
198 -----------
199  
200 04, 25 August 1994:
201  
202 LUSHOW has been corrected, so that if t,  l or h quarks (or d* or u* 
203 quarks masked as l or h) are given with masses that vary from event
204 to event (a Breit-Wigner shape, e.g.), the current mass rather than the
205 nominal mass is used to define the cut-off scales of parton shower
206 evolution.
207  
208 -----
209  
210 LULOGO has been modified to take into account that a new PYTHIA/JETSET
211 description has been published in
212 T. Sjostrand, Computer Phys. Commun. 82 (1994) 74
213 and is from now on the standard reference to these two programs.   
214  
215 -----------
216  
217 05, 27 January 1995:
218  
219 LUCELL has been corrected, in that in the option with smearing of
220 energy rather than transverse energy, the conversion factor between 
221 the two was applied in the wrong direction.
222  
223 -----
224  
225 LUSHOW has been corrected in one place where the PMTH array was 
226 addressed with the wrong order of the indices. This affected quark 
227 mass corrections in the matching to the three-jet matrix elements.
228  
229 -----
230  
231 An additional check has been included in LUBOEI that there are at
232 least two particles involved in the Bose-Einstein effects. (No
233 problem except in some bizarre situations.) 
234  
235 -----------
236  
237 06, 20 February 1995:
238  
239 A new option has been added for the behaviour of the running
240 alpha_em(Q2) in ULALEM. This is not added as a true physics scenario, 
241 but only to produce results with a given, fixed value for the hard 
242 events, while still keeping the conventional value in the Q2=0 limit.   
243 MSTU(101) = 2 : if Q2 is less than PARU(104) then alpha_em is
244     assigned the value PARU(101) (=1/137), while for Q2 above
245     PARU(104) the fixed value PARU(103) (=1/128.8) is used.
246 PARU(103) : (D=0.007764=1/128.8) alpha_em used for hard processes
247     in the option MSTU(101)=2.
248 PARU(104) : (D=1 GeV^2) dividing line for 'low' and 'high'
249     Q2 values in the MSTU(101)=2 option of ULALEM.
250  
251 Additionally, the G_F constant has been added to the parameter list.
252 PARU(105) : (D=1.16639E-5 GeV^-2) G_F, the Fermi constant of weak
253     interactions.
254  
255 -----
256  
257 The LULOGO routine has been updated to reflect my change of
258 affiliation.
259  
260 -----------
261  
262 07, 21 June 1995:
263  
264 Header and LULOGO have been updated with respect to phone number
265 and WWW access.
266  
267 -----
268  
269 The PHEP and VHEP variables in the /HEPEVT/ common block are now 
270 assumed to be in DOUBLE PRECISION, in accord with the proposed
271 LEP 2 workshop addendum to the standard.
272  
273 -----
274  
275 In LUTEST a missing decimal point on the energy check has been
276 reinstated (0001 -> 0.0001).
277  
278 -----
279  
280 In LUINDF the expression PR/(Z*W) has been protected against vanishing
281 denominator.
282  
283 -----------
284  
285 08, 23 August 1995:
286  
287 Check against division by zero in LUSHOW.
288  
289 -----------------------------------------------------------------
290  
291 3) Changes in PYTHIA 5.7
292  
293 -----------
294  
295 00, 13 December 1993: baseline version.
296  
297 -----------
298  
299 01, 27 January 1994:
300  
301 The machinery to handle gamma-gamma interactions is expanded.
302 In particular, several new options have been added to MSTP(14).
303 The updated description of this variable reads as follows.
304  
305 MSTP(14) : (D=0) structure of incoming photon beam or target
306     (does not affect photon inside electron, only photons appearing 
307     as argument in the PYINIT call).
308   = 0 : a photon is assumed to be point-like (a  direct photon), 
309       i.e. can only interact in processes which explicitly contain 
310       the incoming photon, such as f_i gamma -> f_i g for gamma-p 
311       interactions. In gamma-gamma interactions both photons are
312       direct, i.e the main process is gamma gamma -> f_i fbar_i.
313   = 1 : a photon is assumed to be resolved, i.e. can only interact
314       through its constituent quarks and gluons, giving either high-pT
315       parton-parton scatterings or low-pT events. Hard processes are 
316       calculated with the use of the full photon parton distributions. 
317       In gamma-gamma interactions both photons are resolved.
318   = 2 : a photon is assumed resolved, but only the VMD piece is 
319       included in the parton distributions, which therefore mainly 
320       are scaled-down versions of the rho0/pi0 ones. Both high-pT
321       parton-parton scatterings and low-pT events are allowed. In 
322       gamma-gamma interactions both photons are VMD-like.
323   = 3 : a photon is assumed resolved, but only the anomalous piece of 
324       the photon parton distributions is included. Only high-pT 
325       parton-parton scatterings are allowed. In gamma-gamma 
326       interactions both photons are anomalous.
327   = 4 : in gamma-gamma interactions one photon is direct and the other 
328       resolved. A typical process is thus f_i gamma -> f_i g. Hard 
329       processes are calculated with the use of the full photon 
330       parton distributions for the resolved photon. Both possibilities 
331       of which photon is direct are included, in event topologies and 
332       in cross sections. This option cannot be used in configurations 
333       with only one incoming photon. 
334   = 5 : in gamma-gamma interactions one photon is direct and the other 
335       VMD-like. Both possibilities of which photon is direct are 
336       included, in event topologies and in cross sections. This option 
337       cannot be used in configurations with only one incoming photon. 
338   = 6 : in gamma-gamma interactions one photon is direct and the other 
339       anomalous. Both possibilities of which photon is direct are 
340       included, in event topologies and in cross sections. This option 
341       cannot be used in configurations with only one incoming photon. 
342   = 7 : in gamma-gamma interactions one photon is VMD-like and the other 
343       anomalous. Only high-pT parton-parton scatterings are allowed.
344       Both possibilities of which photon is VMD-like are included,
345       in event topologies and in cross sections. This option cannot be 
346       used in configurations with only one incoming photon. 
347   Note: a complete description requires separate runs for the components
348       above, i.e. it is not possible to mix them in a single run. Our
349       best understanding of gamma-p interactions [Sch93,Sch93a] is to
350       have three separate components, 0 + 2 + 3. A simpler alternative
351       is based on two only, 0 + 1. Our best understanding of gamma-gamma
352       interactions [in preparation] requires six separate components, 
353       0 + 2 + 3 + 5 + 6 + 7. A simpler alternative is based on three 
354       only, 0 + 1 + 4.
355  
356 In addition, one new option has been introduced and a few internal
357 variables modified.
358  
359 MSTP(59) : (D=0) possibility to modify the Q2 scale used in the
360       anomalous parton distributions of the photon, as used in the
361       options MSTP(14) = 3, 6 and 7.
362   = 0 : no change of Q2 scale compared to what is normally used.
363   = 1 : the input Q2 scaled is divided by PARP(59)**2 to define
364       the Q2 scale used as argument for the anomalous parton 
365       distributions.
366  
367 PARP(59) : (D=1.) rescaling factor used for the Q2 argument of the 
368       anomalous parton distributions of the photon, see MSTP(59).
369  
370 MINT(105) : is MINT(103) or MINT(104), depending on which side
371     of the event currently is being studied.
372  
373 MINT(107), MINT(108) : if either or both of the two incoming particles 
374     is a photon, then the respective value gives the nature assumed for
375     that photon. The code follows the one used for MSTP(14):
376   = 0 : direct photon.
377   = 1 : resolved photon.
378   = 2 : VMD-like photon.
379   = 3 : anomalous photon.
380  
381 MINT(109) : is either MINT(107) or MINT(108), depending on which side
382     of the event currently is being studied.
383  
384 VINT(282) : no longer used.
385  
386 VINT(283), VINT(284): virtuality scale at which an anomalous photon 
387     on the beam or target side of the event is being resolved. More
388     precisely, it gives the p_T^2 pf the gamma -> q qbar vertex.
389  
390 -----
391  
392 A number of bugs have also been corrected:
393 * Jet + low-pT event generation could give incorrect cross section 
394   information with PYSTAT(1) at low energies. The event generation
395   itself is correct. (The error was introduced when variable energies 
396   became allowed.)
397 * Introduce rejection of top events where top mass (in the tails of the
398   Breit-Wigner distribution) is too low to allow decays t -> W + b.
399 * Plus a few minor bugs, probably harmless.   
400  
401 -----------
402  
403 02, 13 February 1994:
404  
405 The interface to PDFLIB has been modified to reflect that 'TMAS' should
406 no longer be set except in first PDFSET call. (Else a huge amount of 
407 irrelevant warning messages are generated by PDFLIB.)
408  
409 -----
410  
411 The STOP statement in a few dummy routines has been modifed to avoid 
412 irrelevant compilation warning messages on IBM mainframes. 
413  
414 -----
415  
416 A few labels have been renumbered. 
417  
418 -----------
419  
420 03, 22 February 1994:
421  
422 Removal of a bug in PYRESD, which could give (under some specific 
423 conditions) errors in the colour flow. 
424  
425 -----------
426  
427 04, 7 April 1994:
428  
429 Process 11 has been corrected, for the part that concerns anomalous
430 couplings (contact interactions) in the q + q' -> q + q' process.
431 The error was present in the expression for u + dbar -> u + dbar
432 and obvious permutations, while u + d -> u + d, u + ubar -> u + ubar
433 and the others were correct. Thanks to J.-J. Dugne, M. Perrottet and 
434 K. Lane for communications on this point.
435  
436 -----
437  
438 The option MSTP(23)=1 for post-facto (x,Q^2) conservation in deep
439 inelastic scattering can give infinite loops when applied to process
440 83, in particular if one asks for the production of a top. (Remember
441 that the standard DIS kinematics only is defined for massless quarks.)
442 Therefore the switch MSTP(23) has been modifed as follows:
443 MSTP(23) : (D=1) (x, Q^2) correction level in DIS.
444     = 0 : no correction procedure applied.
445     = 1 : correction applied for process 10, but not for process 83.
446     = 2 : correction applied both for process 10 and 83. This latter
447         option could still work fine for charm and bottom, if the
448         energy is sufficient.
449  
450 -----
451  
452 PYRESD is modified to ensure isotropic angular distributions in the 
453 decays of the top or a fourth generation particle, i.e. in t -> b + W+. 
454 This may not be the correct distribution but, unless explicit knowledge
455 exists for a given process, this should always be the default. 
456  
457 -----
458  
459 In processes 16, 20, 31 and 36 the W propagator has been modified to 
460 include s-dependent widths in the Breit-Wigner shape. The most notable
461 effect is a suppression of the low-mass tail of the W mass spectrum.
462  
463 -----
464  
465 When PDFLIB is used, PDFSET is now only called whenever a different
466 structure function is requested. For pp events therefore only one call
467 is made, while gamma-p interactions still involves a call to PDFSET 
468 for each STRUCTM one, since gamma and p structure functions have to be 
469 called alternatingly. To this end, MINT(93) is reset to 
470 1000000 * Nptype + 1000 * Ngroup + Nset after each PDFSET call.
471  
472 -----
473  
474 In a few places, a dot has been moved from the end of one line to the
475 beginning of the next continuation line, or the other way around, 
476 to keep together tokens such as .EQ. or .AND., since some debuggers
477 may otherwise complain. Also some other purely cosmetics changes 
478 for the same reason.  
479  
480 A number of minor errors have been corrected.
481  
482 -----------
483  
484 05, 15 July 1994:
485  
486 A new option has been introduced, MSTP(14)=10, whereby it is possible 
487 to obtain a mixture of the various allowed photon components. For
488 gamma-hadron collisions, this means a mixture of VMD, direct and 
489 anomalous events, for gamma-gamma collisions a mixture of VMD*VMD, 
490 VMD*direct, VMD*anomalous, direct*direct, direct*anomalous and
491 anomalous*anomalous. The mixture is properly given according to 
492 the relative cross sections. 
493  
494 Note that this introduces a completely new layer of administration in
495 PYTHIA. For instance, a subprocess such as q + g -> q + g is allowed
496 in the VMD*VMD, VMD*anomalous and anomalous*anomalous classes, but
497 appear with different sets of parton distributions and with different
498 pT cut-offs. In order to handle this, various information is initialized
499 separately for each event class, and subsequently saved and restored 
500 as the generation switches back and forth between the event classes. 
501 This introduces some limitations on what you may and may not do.
502  
503 First of all, the MSTP(14) switch is only applicable for incoming photon
504 beams, i.e. when 'gamma' is the argument in the PYINIT call. A
505 convolution with the bremsstrahlung photon spectrum in an electron beam
506 may come one day, but not in the immediate future.
507  
508 Secondly, the machinery has only been set up to generate standard
509 QCD physics, specifically either 'minimum bias' one or high-pT jets.
510 For minimum bias, you are not allowed to use the CKIN variables at all.
511 This is not a major limitation, since it is in the spirit of minimum
512 bias physics not to impose any contraints on allowed jet production. 
513 (If you still do, these cuts will be ineffective for the VMD processes 
514 but take effect for the other ones, giving inconsistencies.) Further, 
515 some variables are internally recalculated and reset: CKIN(1), CKIN(3), 
516 CKIN(5), CKIN(6), MSTP(57), MSTP(85), PARP(2), PARP(81), PARP(82), 
517 PARU(115) and MDME(22,J). These can not be modified without changing
518 PYINPR and recompiling the program. The minimum bias physics option 
519 is obtained by default; by switching from MSEL=1 to MSEL=2 also the 
520 elastic and diffractive components of the VMD part are included. 
521 High-pT jet production is obtained by setting the CKIN(3) cut-off 
522 larger than the (energy-dependent) cut-off scales for the VMD, direct 
523 and  anomalous components; typically this means at least 3 GeV. For 
524 lower input CKIN(3) the program will automatically switch back to 
525 minimum bias physics.
526  
527 Finally, pileup events are not at all allowed.  
528  
529 Here is a survey of common block variables affected:
530  
531 MSTP(14) (D=0) strucure of incoming photon beam or target;
532     see description above for PYTHIA 5.701.
533     = 10 : new option where the VMD, direct and anomalous components 
534         are automatically mixed, as described above. Works equally well
535         for gamma-p and gamma-gamma.
536  
537 MSTI(9) : event class used in current event.
538     = 1 : VMD (for gamma-p) or VMD*VMD (for gamma-gamma).
539     = 2 : direct (for gamma-p) or VMD*direct (for gamma-gamma).
540     = 3 : anomalous (for gamma-p) or VMD*anomalous (for gamma-gamma).
541     = 4 : direct*direct (for gamma-gamma).
542     = 5 : direct*anomalous (for gamma-gamma).
543     = 6 : anomalous*anomalous (for gamma-gamma).
544  
545 MINT(121) : number of separate event classes to initialize and mix.
546     = 1 : the normal value.
547     = 3 : for a gamma-hadron interaction when MSTP(14)=10.
548     = 6 : for a gamma-gamma interaction when MSTP(14)=10.
549  
550 MINT(122) : event class used in current event. Code as explained for
551     MSTI(9).
552  
553 MINT(123) : event class used in the current event, with the same list 
554     of possibilities as MSTP(14), except that MSTP(14) = 1, 4 or 10 
555     do not appear.  
556  
557 VINT(285) : the CKIN(3) value provided by the user at initialization;
558     subsequently CKIN(3) may be overwritten (for MSTP(14)=10) but 
559     VINT(285) stays.
560  
561 In addition, the structure of the initialization has been partly
562 reorganized. The routine PYEVKI has been removed, new routines
563 PYINBM, PYINPR and PYSAVE created, and some material has been moved 
564 to or from PYINIT, PYINRE and PYINKI.
565  
566 SUBROUTINE PYINBM : to read in and identify beam and target particles
567     and frame as given in the PYINIT call (used to be done in PYINKI).
568  
569 SUBROUTINE PYINKI(MODKI) : to set up event kinematics, either at
570     initialization (MODKI=0) or for each separate event when varying
571     kinematics (MODKI=1). (The latter task used to be done in PYEVKI.)
572  
573 SUBROUTINE PYINPR : to set up the partonic subprocesses selected with
574     MSEL and, for gamma-p and gamma-gamma, MSTP(14).
575  
576 SUBROUTINE PYSAVE : saves and restores parameters and cross section
577     values between the 3 gamma-p and the 6 gamma-gamma alternatives
578     of MSTP(14)=10. Also makes a random choice for each new event
579     between the allowed alternatives.
580  
581 Among other changes, note that PYSTAT(1) now has been extended so
582 that the subdivision into the various gamma-p and gamma-gamma classes
583 is shown.     
584  
585 -----
586  
587 Further changes of particular relevance for gamma-p and gamma-gamma,
588 but independent of the major revisions above:
589  
590 MSTP(59) and PARP(59) have been removed. Instead the following options
591 are available:
592  
593 MSTP(15) : (D=5) possibility to modify the nature of the anomalous
594 photon component, in particular with respect to the scale choices and
595 cut-offs of hard processes.
596     = 0 : none, i.e. the same treatment as for the VMD component.
597     = 1 : evaluate the anomalous structure functions at a scale
598         Q2/PARP(17)^2.
599     = 2 : as =1, but instead of PARP(17) use PARP(81)/PARP(15) or
600         PARP(82)/PARP(15), depending on MSTP(82) value.
601     = 3 : evaluate anomalous structure function as
602         f^(anom)(x, Q2, p_0^2) - f^(anom)(x, Q2, r^2*Q2)
603         with r = PARP(17).
604     = 4 : as =3, but instead of PARP(17) use PARP(81)/PARP(15) or
605         PARP(82)/PARP(15), depending on MSTP(82) value.
606     = 5 : use larger pTmin for anomalous component than for VMD one,
607         but otherwise no difference.
608  
609 PARP(17) : (D=1) rescaling factor used as described for MSTP(15).
610  
611 MSTP(51) : new option added.
612     = 11 : the GRV p LO parametrization.
613  
614 MSTP(53) : new option added.
615     = 3 : the GRV pi LO parametrization.  
616  
617 MSTP(56) : new option added.
618     = 3 : when the anomalous photon structure function is requested,
619         the homogeneous solution is provided, evolved from a starting
620         value PARP(15) to the requested Q scale. The homogeneous 
621         solution is normalized so that the net momentum is unity,
622         i.e. any factors of alpha_em/2pi and charge have been left out. 
623         The flavour of the original q is given in MSTP(55) (1, 2, 3, 4 
624         or 5 for d, u, s, c or b); the value 0 gives a mixture 
625         according to squared charge, with the exception that c and b 
626         are only allowed above the respective mass threshold (Q > m_q).
627         The four-flavour Lambda value is assumed given in PARP(1);
628         it is automatically recalculated for 3 or 5 flavours at 
629         thresholds. This option is not intended for standard event 
630         generation, but is useful for some theoretical studies.
631   
632 -----
633  
634 Option MSTP(92)=5 for beam remnant treatment erroneously missed some
635 statements which now have been inserted. 
636  
637 Further, new options have been added for the splitting of momentum 
638 between two beam remnants. MSTP(92) keeps its current role for the 
639 production of diquark or quark jets. However, for the splitting into 
640 a hadron plus a quark/diquark jet, MSTP(94) should now be used.
641  
642 MSTP(94) : (D=2) (C) energy partitioning in hadron or resolved photon
643 remnant, when this is split into a hadron plus a remainder-jet. The 
644 energy fraction chi is taken by one of the two objects, with 
645 conventions as described below or for PARP(95) and PARP(97).
646     = 1 : 1 for meson or resolved photon, 2(1-chi) for baryon, i.e.
647         simple counting rules.
648     = 2 : (k+1)*(1-chi)**k, with k as given in PARP(95) or PARP(97).
649     = 3 : the chi of the hadron is selected according to the normal
650         fragmentation function used for the hadron in jet fragmentation,
651         see MSTJ(11). The possibility of a changed fragmentation 
652         function shape in diquark fragmentation (see PARJ(45)) is not 
653         included.  
654     = 4 : as =3, but the shape is changed as allowed in diquark
655         fragmentation (see PARJ(45)); this change is here also allowed 
656         for meson production.    
657  
658 -----
659  
660 In PYDIFF the recoiling gluon energy is calculated in a numerically more
661 stable fashion. 
662  
663 -----
664  
665 A counter has been added to PYSSPA to avoid infinite loops in the
666 angular ordering constraint due to interference with the final state
667 colour charges.
668  
669 -----------
670  
671 06, 25 August 1994:
672  
673 New processes 167 and 168 have been added for the contact interaction
674 production of d* or u* excited quarks
675 167   q q' -> q" d*
676 168   q q' -> q" u*
677 where the different allowed quark and antiquark combinations are given
678 according to eqs. (15) - (19) in U. Baur, M. Spira and P.M. Zerwas,
679 Phys. Rev. D42 (1990) 815. The d* and u* are defined in the same
680 way as for processes 147 and 148. Thus one needs to put MSTP(6)=1 to
681 use l (7) and h (8) for representing the d* and u*. The couplings of
682 the allowed decay channels are given by PARU(157) - PARU(159), and
683 the Lambda scale parameter by PARU(155).
684  
685 At the same time, some minor changes has been introduced in the code
686 for processes 147 and 148, for uniformity.
687  
688 -----
689  
690 Option MSTP(57)=3 now also allows a dampening of pi+- parton 
691 distributions.
692  
693 -----
694  
695 A few minor errors have been corrected.
696  
697 -----------
698  
699 07, 20 October 1994:
700  
701 A major bug discovered in processes 121 and 122 (and thus also affecting
702 181, 182, 186 and 187), g g or q qbar -> Q Qbar H: the kinematics was 
703 incorrectly handed on to the Kunszt matrix elements. This affected the
704 default option Q = t, but effects were especially dramatic when the
705 alternative Q = b was used. 
706  
707 The choice of appropriate Q2 scale for structure functions introduces
708 a further uncertainty in cross sections for the processes above. So long
709 as only t quarks are considered, the t mass is a reasonable choice, but
710 for the Q = b alternative this is presumably too low. Therefore new 
711 options have been introduced as below, with the default behaviour 
712 changed (the old one is obtainable with MSTP(39)=1).
713  
714 MSTP(39) : (D=2) choice of Q2 scale for structure functions and initial
715 state parton showers in processes g g or q qbar -> Q Qbar H.
716     = 1 : m_Q**2.
717     = 2 : max(mT_Q**2 , mT_Qbar**2) = 
718         m_Q**2 + max(pT_Q**2 , pT_Qbar**2).
719     = 3 : m_H**2.
720     = 4 : shat = (p_H + p_Q + p_Qbar)**2.
721  
722 -----
723  
724 Another important bug corrected in the calculation of the reduction of
725 t+tbar cross section when decay modes are forced. This occured when both
726 t and tbar produced a W, and W+ and W- decay modes were set differently.
727  
728 -----------
729  
730 08, 25 October 1994:
731  
732 A few further places changed to make processes 181, 182, 186 and 187 
733 work (see version 5.707 above).
734  
735 -----------
736  
737 09, 26 October 1994:
738  
739 The matrix element for f + fbar -> W+ + W- has been replaced, using the
740 formulae of 
741 D. Bardin, M. Bilenky, D. Lehner, A. Olchevski and T. Riemann,
742 CERN-TH.7295/94,
743 but with the dependence on the t-hat variable not integrated out
744 (D. Bardin, private communication).
745 This avoids some problems encountered in the old expressions when 
746 one or both W's were far off the mass shell.
747  
748 -----
749  
750 Change in calls to PDFLIB, so that the input Q is always at least the 
751 QMIN of the respective set.
752  
753 -----
754  
755 Extra protection against infinite loops in PYSSPA.
756  
757 -----------
758  
759 10, 27 January 1995:
760  
761 The dimensions of the HGZ array in PYRESD has been expanded to avoid
762 accidental writing outside the bounds.
763  
764 -----
765  
766 VINT(41) - VINT(66) are saved and restored in PYSCAT, for use in low-pT
767 events, when beam remnant treatment has failed (with nonzero MINT(57)).
768  
769 -----
770  
771 The routine PYSTGH has been replaced by the routine PYSTHG. This 
772 contains an improved parametrization of the homogeneous evolution 
773 of an anomalous photon from some given initial scale. The argument
774 NF of the PYSTGH routine has been removed; now Lambda is always
775 automatically converted to the relevant NF-flavour value from its
776 4-flavour one, at flavour thresholds. 
777  
778 -----------
779  
780 11, 20 February 1995:
781  
782 New possibilities have been added to switch between electroweak
783 couplings being expressed in terms of a running alpha_em(Q2)  or 
784 in terms of a fixed Fermi constant G_F. This affects both decay widths
785 and process cross sections, in the routines PYINRE, PYRESD, PYWIDT and 
786 PYSIGH. There are three main options, with default agreeing with the
787 old standard.
788  
789 MSTP(8) : (D=0) choice of electroweak parameters to use in decay
790     widths of resonances (W, Z, H, ...) and cross sections (production
791     of W's, Z's, H's, ...). 
792     = 0 : everything is expressed in terms of a running alpha_em(Q2) 
793         and a fixed sin^2(theta_W), i.e. G_F is nowhere used.
794     = 1 : a replacement is made according to 
795         alpha_em(Q2) -> sqrt(2) G_F m_W^2 sin^2(theta_W) / pi
796         in all widths and cross sections. If G_F and m_Z are considered 
797         as given, this means that sin^2(theta_W) and m_W are the only 
798         free electroweak parameter.
799     = 2 : a replacement is made as for =1, but additionally 
800         sin^2(theta_W) is constrained by the relation 
801         sin^2(theta_W) = 1 - m_W^2/m_Z^2
802         This means that m_W remains as a free parameter, but that the
803         sin^(theta_W) value in PARU(102) is never used, EXCEPT in
804         the vector couplings in the combination
805         v = a - 4 sin^2(theta_W) e.
806         This degree of freedom enters e.g. for forward-backward
807         asymmetries in Z^0 decays.
808 Note : This option does not affect the emission of real photons in the 
809     initial and final state, where alpha_em is always used. However, 
810     it does affect also purely electromagnetic hard processes, such as 
811     q + qbar -> gamma + gamma.
812  
813 -----
814  
815 The option MSTP(37)=1, with running quark masses in couplings to Higgs
816 bosons, only works when alpha_s is allowed to run (so one can define
817 a Lambda value). Therefore a check has been introduced in PYWIDT and
818 PYSIGH that the option MSTP(37)=1 is only executed if additionally
819 MSTP(2) is 1 or bigger.
820  
821 -----
822  
823 Some non-physics changes have been made in the RKBBV and STRUCTM codes
824 so as to avoid some (in principle harmless) compiler warnings.
825  
826 -----------
827  
828 12, 15 March 1995:
829  
830 A serious error has been corrected in the MSTP(173)=1 option, i.e. when
831 the program is run with user-defined weights that should compensate for 
832 a biased choice of variable beam energies. This both affected the 
833 relative admixture of low- and high-pT events and the total cross 
834 section obtained by Monte Carlo integration. (PYRAND changed.)
835  
836 In order to improve the flexibility and efficiency of the variable-energy
837 option, the user should now set PARP(174) before the PYINIT call, and
838 thereafter not change it. This allows PARP(173) weights of arbitrary
839 size. (PYRAND and PYMAXI changed.)
840 PARP(174) : (D=1.) maximum event weight that will be encountered in
841     PARP(173) during the course of a run with MSTP(173)=1; to be used 
842     to optimize the efficiency of the event generation. It is always 
843     allowed to use a larger bound than the true one, but with a 
844     corresponding loss in efficiency.
845  
846 MSTI(5) (and MINT(5)) are now changed so they count the number of 
847 successfully generated events, rather than the number of tries made.
848 This change only affects runs with variable energies, MSTP(171)=1 and 
849 MSTP(172)=2, where MSTI(61)=1 signals that a user-provided energy
850 has been rejected in the weighting. This change also affects PARI(2),
851 which becomes the cross section per fully generated event. (PYEVNT
852 changed.)
853  
854 -----
855  
856 The option MSTP(14)=10 has now been extended so that it also works for
857 deep inelastic sacattering of an electron off a (real) photon, i.e.
858 process ISUB = 10. What is obtained is a mixture of the photon acting
859 as a vector meson and it acting as an anomalous state. This should 
860 therefore be the sum of what can be obtained with MSTP(14)=2 and =3.
861 It is distinct from MSTP(14)=1 in that different sets are used for
862 the parton distributions - in MSTP(14)=1 all the contributions to the
863 photon distributions are lumped together, while they are split in 
864 VMD and anomalous parts for MSTP(14)=10. Also the beam remnant treatment
865 is different, with a simple Gaussian distribution (at least by default)
866 for MSTP(14)=1 and the VMD part of MSTP(14)=10, but a powerlike 
867 distribution d(kT^2)/kT^2 between PARP(15) and Q for the anomalous
868 part of MSTP(14)=10. (PYINIT, PYINPR and PYSTAT changed.)
869  
870 To access this option for e and gamma as incoming beams, it is only 
871 necessary to set MSTP(14)=10 and keep MSEL at its default value.
872 Unlike the corresponding option for gamma-p and gamma-gamma, no
873 cuts are overwritten, i.e. it is still the responsability of the user
874 to set these appropriately. Those especially appropriate for DIS usage
875 are CKIN(21)-CKIN(22) or CKIN(23)-CKIN(24) for the x range (former or
876 latter depending on which side is the incoming real photon), 
877 and CKIN(35)-CKIN(36) for the Q2 range. A further new option has been 
878 added (in PYKLIM):
879 CKIN(39), CKIN(40) : (D=4., -1. GeV^2) the W2 range allowed in DIS
880     processes, i.e. subprocess number 10. If CKIN(40) < 0., the upper
881     limit is inactive. Here W2 is defined in terms of W2 = Q2 * (1-x)/x.
882     This formula is not quite correct, in that (i) it neglects the 
883     target mass (for a proton), and (ii) it neglects initial-state 
884     photon radiation off the incoming electron. It should be good 
885     enough for loose cuts, however. 
886  
887 A warning about the usage of PDFLIB for photons. So long as MSTP(14)=1, 
888 i.e. the photon is not split up, PDFLIB is accessed by MSTP(56)=2 and
889 MSTP(55) = parton distribution set, as described in the manual. 
890 However, when the VMD and anomalous pieces are split, the VMD part
891 is based on a rescaling of pion distributions by VMD factors
892 (except for the SaS sets, that already come with a separate VMD piece).
893 Therefore, to access PDFLIB for MSTP(14)=10, it is not correct to set 
894 MSTP(56)=2 and a photon distribution in MSTP(55). Instead, one should
895 put MSTP(56)=2, MSTP(54)=2 and a pion distribution code in MSTP(53),
896 while MSTP(55) has no function. The anomalous part is still based on
897 the SaS parametrization, with PARP(15) as main free parameter.   
898  
899 -----
900  
901 A change has been made in PYREMN to reduce the possibility of infinite 
902 loops. 
903  
904 -----------
905  
906 13, 22 March 1995:
907  
908 The SaS parton distributions of the photons are now available.
909 For details on these sets, see 
910 G.A. Schuler and T. Sjostrand,
911 "Low- and high-mass components of the photon distribution functions",
912 CERN-TH/95-62 and LU TP 95-6.
913 There are four new sets. These differ in that two use a Q0=0.6 GeV and
914 two a Q0=2 GeV, and in that two use the DIS and two the MSbar conventions 
915 for the dominant non-leading contributions. (However, the fits are 
916 formally still leading-order, in that not all next-to-leading
917 contributions have been included.) New default is the SaS 1D set. 
918 Furthermore, for the definition of F_2^gamma, additional terms appear 
919 that do not form part of the parton distributions itself. To partly 
920 take this into account, an additional doubling of the possibilities 
921 has been included. These possibilites can be accesed with MSTP(55):
922  
923 MSTP(55) : (D=5) choice of parton-distribution set of the photon;
924     see also MSTP(56).
925     = 1  : Drees-Grassie.
926     = 5  : SaS 1D (in DIS scheme, with Q0=0.6 GeV).
927     = 6  : SaS 1M (in MSbar scheme, with Q0=0.6 GeV).
928     = 7  : SaS 2D (in DIS scheme, with Q0=2 GeV).
929     = 8  : SaS 2M (in MSbar scheme, with Q0=2 GeV).
930     = 9  : SaS 1D (in DIS scheme, with Q0=0.6 GeV).
931     = 10 : SaS 1M (in MSbar scheme, with Q0=0.6 GeV).
932     = 11 : SaS 2D (in DIS scheme, with Q0=2 GeV).
933     = 12 : SaS 2M (in MSbar scheme, with Q0=2 GeV).
934     Note 1 : sets 5 - 8 use the parton distributions of the respective
935         set, and nothing else. These are appropriate for most 
936         applications, e.g. jet production in gamma-p and gamma-gamma
937         collisions. Sets 9 - 12 instead are appropriate for
938         gamma*-gamma processes, i.e. DIS scattering on a photon, 
939         as measured in F_2^gamma. Here the anomalous contribution
940         for c and b quarks are handled by the Bethe-Heitler formulae,
941         and the direct term is artificially lumped with the anomalous
942         one, so that the event simulation more closely agrees with what
943         will be experimentally observed in these processes. The agreement
944         with the F_2^gamma parametrization is still not perfect, e.g.
945         in the treatment of heavy flavours close to threshold. 
946     Note 2 : Sets 5 - 12 contain both VMD pieces and anomalous pieces,
947         separately parametrized. Therefore the respective piece is 
948         automatically called, whatever MSTP(14) value is used to select 
949         only a part of allowed photon interactions. For other sets 
950         (set 1 above or PDFLIB sets), usually there is no corresponding 
951         subdivision. Then an option like MSTP(14)=2 (VMD part of photon 
952         only) is based on a rescaling of the pion distributions, while
953         MSTP(14)=3 gives the SaS anomalous parametrization.     
954     Note 3 : Formally speaking, the k0 (or p0) cut-off in PARP(15) 
955         need not be set in any relation to the Q0 cut-off scales used by 
956         the various parametrizations. Indeed, due to the familiar scale
957         choice ambiguity problem, there could well be some offset between
958         the two. However, unless you know what you are doing, it is 
959         strongly recommended that you let the two agree, i.e. set
960         PARP(15)=0.6 for the SaS 1 sets and =2 for the SaS 2 sets.
961  
962 PARP(15) : (D=0.6 GeV) default value changed for k0 cut-off for 
963     separation between direct, VMD and anomalous photons; see Note 3 
964     for MSTP(55) above.       
965         
966 The generic routine PYSTFU has been rewritten to handle the interfacing.
967 The old routines PYSTAG, PYSTGS, PYDILN and PYSTHG have been removed.
968 Instead the routines of the SaSgam library have been inserted. In order
969 to avoid any clashes, the routines SAS*** have been renamed PYG***.
970 Thus new routines are PYGGAM, PYGVMD, PYGANO, PYGBEH and PYGDIR.
971 The common block SASCOM is renamed PYINT8. If you want to use the
972 parton distributions for standalone purposes, you are encouraged to
973 use the original SaSgam routines rather than going the way via the
974 Pythia adaptations. 
975  
976       COMMON/PYINT8/XPVMD(-6:6),XPANL(-6:6),XPANH(-6:6),XPBEH(-6:6),
977      &XPDIR(-6:6)
978 Purpose: to store the various components of the parton distributions
979     when the PYGGAM routine is called.
980 XPVMD(KFL) : gives distributions of the VMD part (rho, omega and
981     phi).      
982 XPANL(KFL) : gives distributions of the anomalous part of light quarks
983     (d, u and s).
984 XPANH(KFL) : gives distributions of the anomalous part of heavy quarks
985     (c and b).
986 XPBEH(KFL) : gives Bethe-Heitler distributions of heavy quarks (c and b).
987     This provides an alternative to XPANH, i.e. both should not be used
988     at the same time.
989 XPDIR(KFL) : gives direct correction to the production of light quarks
990     (d, u and s). This term is nonvanishing only in the MSbar scheme,
991     and is applicable for F_2^gamma rather than for the parton
992     distributions themselves.
993  
994 -----
995  
996 PYDOCU has been corrected so that PARI(2) refers to the full cross 
997 section for gamma-p and gamma-gamma processes, rather than that of 
998 the latest subprocess considered.
999  
1000 An additional check has been inserted into PYREMN. 
1001  
1002 -----------
1003  
1004 14, 23 March 1995:
1005  
1006 Some minor modifications to PYSTFU and PYGGAM in the wake of the
1007 changes of the previous version.
1008  
1009 -----------
1010  
1011 15, 24 April 1995:
1012  
1013 An unfortunate choice of default values has been corrected:
1014 the old MSTP(3)=2 value implied that Lambda_QCD was entirely based on 
1015 the Lambda value of the proton structure function; also e.g. for e+e-
1016 annihilation events. Thus the Lambda in PARJ(81) was overwritten,
1017 i.e. did not keep the value required by standard phenomenology, which
1018 typically gave too narrow jets. (While switching to MSTP(3)=1 it worked
1019 fine.) In the modified option MSTP(3)=2 this has been corrected, to 
1020 better agree with user expectations. Change affects PYINIT and
1021 PYRESD. (See further version 16 for additional changes.)
1022  
1023 MSTP(3) : (D=2) choice of Lambda_QCD values.
1024     = 1 : separate for hard scattering, initial showers and final 
1025         showers, as before. Additionally separate for resonance
1026         decays, given in PARP(3).
1027     = 2 : most Lambda values are set to that of the proton structure
1028         function used, except for the Lambda used in the decay of
1029         a resonance (as treated in PYRESD). There the PARP(3) value
1030         is used, with default as in e+e-.
1031     = 3 : all Lambda values are set to that of the proton structure
1032         function used, as was the case for =2 before.
1033  
1034 PARP(3) : (D=0.29 GeV) the Lambda value used in timelike parton showers 
1035     in the decay of a resonance (in PYRESD).
1036  
1037 -----
1038  
1039 The form for PTMANO, the pTmin for anomalous processes, as used in 
1040 PYINPR when processes are mixed for gamma-p or gamma-gamma events,
1041 has been updated to match (as well as can be expected) the SaS 1D 
1042 photon distributions.
1043  
1044 -----------
1045  
1046 16, 30 June 1995:
1047  
1048 The strategy for the changes to MSTP(3) in version 15 above have been
1049 modified for better transparency. The parameter PARP(3) has been removed,
1050 and instead PARP(72) has been introduced. Now PARJ(81) is used for
1051 resonance decays (including e.g. Z0 decay, from which it is determined),
1052 and PARP(72) for other timelike showers. PARJ(81) is not overwritten
1053 for MSTP(3) = 2, but only for = 3. Changes affect PYINIT, PYEVNT and
1054 PYRESD.
1055  
1056 PARP(72) : (D=0.25 GeV) the Lambda value used in timelike parton showers 
1057     except in the decay of a resonance.
1058  
1059 -----
1060  
1061 A new multiplicative factor has been introduced for the hard scattering 
1062 in PYSIGH.
1063  
1064 PARP(34) : (D=1.) the Q2 scale defined by MSTP(32) is multiplied by PARP(34)
1065     when it is used as argument for structure functions and alpha_s at the
1066     hard interaction. It does not affect alpha_s when MSTP(33)=3, nor
1067     does it change the Q2 argument of parton showers.
1068  
1069 -----
1070  
1071 PYREMN has been corrected for occasional too large boost factors.
1072  
1073 An error in PYSIGH for process 148 has been corrected.
1074  
1075 The MSTP(62)=1 option of PYSSPA is modified to avoid division by zero.
1076  
1077 Header has been updated with WWW-information.
1078  
1079 -----------
1080  
1081 17, 23 August 1995:
1082  
1083 MIN1, MIN2, MAX1, MAX2, MINA and MAXA in PYSIGH have had an extra M 
1084 prefixed to avoid confusion with Fortran functions.
1085  
1086 Protect against MDCY(0,1) being accessed in PYSIGH.
1087  
1088 Protect against THB=0 in PYRAND.
1089  
1090 Protect against YSTMAX-YSTMIN = 0 in PYSIGH.
1091  
1092 Check for moved leptoquark at beginning of PYRESD just like for
1093 other particles with colour.    
1094  
1095 -----------------------------------------------------------------