Moved to $(FLUPRO) directory.
[u/mrichter/AliRoot.git] / DPMJET / pythia6115.f
1 C*********************************************************************
2 C*********************************************************************
3 C*                                                                  **
4 C*                                                    March 1997    **
5 C*                                                                  **
6 C*           The Lund Monte Carlo for Hadronic Processes            **
7 C*                                                                  **
8 C*                        PYTHIA version 6.1                        **
9 C*                                                                  **
10 C*                        Torbjorn Sjostrand                        **
11 C*                Department of Theoretical Physics 2               **
12 C*                         Lund University                          **
13 C*               Solvegatan 14A, S-223 62 Lund, Sweden              **
14 C*                    phone +46 - 46 - 222 48 16                    **
15 C*                    E-mail torbjorn@thep.lu.se                    **
16 C*                                                                  **
17 C*                          SUSY parts by                           **
18 C*                         Stephen Mrenna                           **
19 C*                    Argonne National Laboratory                   **
20 C*          9700 South Cass Avenue, Argonne, IL 60439, USA          **
21 C*                   phone + 1 - 630 - 252 - 7615                   **
22 C*                    E-mail mrenna@hep.anl.gov                     **
23 C*                                                                  **
24 C*         Several parts are written by Hans-Uno Bengtsson          **
25 C*          PYSHOW is written together with Mats Bengtsson          **
26 C*     advanced popcorn baryon production written by Patrik Eden    **
27 C*     CTEQ 3 parton distributions are by the CTEQ collaboration    **
28 C*      GRV 94 parton distributions are by Glueck, Reya and Vogt    **
29 C*   SaS photon parton distributions together with Gerhard Schuler  **
30 C*     g + g and q + qbar -> t + tbar + H code by Zoltan Kunszt     **
31 C*         MSSM Higgs mass calculation code by M. Carena,           **
32 C*           J.R. Espinosa, M. Quiros and C.E.M. Wagner             **
33 C*         PYGAUS adapted from CERN library (K.S. Kolbig)           **
34 C*                                                                  **
35 C*   The latest program version and documentation is found on WWW   **
36 C*       http://www.thep.lu.se/tf2/staff/torbjorn/Pythia.html       **
37 C*                                                                  **
38 C*              Copyright Torbjorn Sjostrand, Lund 1997             **
39 C*                                                                  **
40 C*********************************************************************
41 C*********************************************************************
42 C                                                                    *
43 C  List of subprograms in order of appearance, with main purpose     *
44 C  (S = subroutine, F = function, B = block data)                    *
45 C                                                                    *
46 C  B   PYDATA   to contain all default values                        *
47 C  S   PYTEST   to test the proper functioning of the package        *
48 C  S   PYHEPC   to convert between /PYJETS/ and /HEPEVT/ records     *
49 C                                                                    *
50 C  S   PYINIT   to administer the initialization procedure           *
51 C  S   PYEVNT   to administer the generation of an event             *
52 C  S   PYSTAT   to print cross-section and other information         *
53 C  S   PYINRE   to initialize treatment of resonances                *
54 C  S   PYINBM   to read in beam, target and frame choices            *
55 C  S   PYINKI   to initialize kinematics of incoming particles       *
56 C  S   PYINPR   to set up the selection of included processes        *
57 C  S   PYXTOT   to give total, elastic and diffractive cross-sect.   *
58 C  S   PYMAXI   to find differential cross-section maxima            *
59 C  S   PYPILE   to select multiplicity of pileup events              *
60 C  S   PYSAVE   to save alternatives for gamma-p and gamma-gamma     *
61 C  S   PYRAND   to select subprocess and kinematics for event        *
62 C  S   PYSCAT   to set up kinematics and colour flow of event        *
63 C  S   PYSSPA   to simulate initial state spacelike showers          *
64 C  S   PYRESD   to perform resonance decays                          *
65 C  S   PYMULT   to generate multiple interactions                    *
66 C  S   PYREMN   to add on target remnants                            *
67 C  S   PYDIFF   to set up kinematics for diffractive events          *
68 C  S   PYDOCU   to compute cross-sections and handle documentation   *
69 C  S   PYFRAM   to perform boosts between different frames           *
70 C  S   PYWIDT   to calculate full and partial widths of resonances   *
71 C  S   PYOFSH   to calculate partial width into off-shell channels   *
72 C  S   PYRECO   to handle colour reconnection in W+W- events         *
73 C  S   PYKLIM   to calculate borders of allowed kinematical region   *
74 C  S   PYKMAP   to construct value of kinematical variable           *
75 C  S   PYSIGH   to calculate differential cross-sections             *
76 C  S   PYPDFU   to evaluate parton distributions                     *
77 C  S   PYPDFL   to evaluate parton distributions at low x and Q^2    *
78 C  S   PYPDEL   to evaluate electron parton distributions            *
79 C  S   PYPDGA   to evaluate photon parton distributions (generic)    *
80 C  S   PYGGAM   to evaluate photon parton distributions (SaS sets)   *
81 C  S   PYGVMD   to evaluate VMD part of photon parton distributions  *
82 C  S   PYGANO   to evaluate anomalous part of photon pdf's           *
83 C  S   PYGBEH   to evaluate Bethe-Heitler part of photon pdf's       *
84 C  S   PYGDIR   to evaluate direct contribution to photon pdf's      *
85 C  S   PYPDPI   to evaluate pion parton distributions                *
86 C  S   PYPDPR   to evaluate proton parton distributions              *
87 C  F   PYCTEQ   to evaluate the CTEQ 3 proton parton distributions   *
88 C  S   PYGRVL   to evaluate the GRV 94L pronton parton distributions *
89 C  S   PYGRVM   to evaluate the GRV 94M pronton parton distributions *
90 C  S   PYGRVD   to evaluate the GRV 94D pronton parton distributions *
91 C  F   PYGRVV   auxiliary to the PYGRV* routines                     *
92 C  F   PYGRVW   auxiliary to the PYGRV* routines                     *
93 C  F   PYGRVS   auxiliary to the PYGRV* routines                     *
94 C  F   PYHFTH   to evaluate threshold factor for heavy flavour       *
95 C  S   PYSPLI   to find flavours left in hadron when one removed     *
96 C  F   PYGAMM   to evaluate ordinary Gamma function Gamma(x)         *
97 C  S   PYWAUX   to evaluate auxiliary functions W1(s) and W2(s)      *
98 C  S   PYI3AU   to evaluate auxiliary function I3(s,t,u,v)           *
99 C  F   PYSPEN   to evaluate Spence (dilogarithm) function Sp(x)      *
100 C  S   PYQQBH   to evaluate matrix element for g + g -> Q + Qbar + H *
101 C                                                                    *
102 C  S   PYMSIN   to initialize the supersymmetry simulation           *
103 C  S   PYAPPS   to determine MSSM parameters from SUGRA input        *
104 C  F   PYRNMQ   to determine running quark masses                    *
105 C  F   PYRNMT   to determine running top mass                        *
106 C  S   PYTHRG   to calculate sfermion third-gen. mass eigenstates    *
107 C  S   PYINOM   to calculate neutralino/chargino mass eigenstates    *
108 C  F   PYRNM3   to determine running M3, gluino mass                 *
109 C  S   PYEIG4   to calculate eigenvalues and -vectors in 4*4 matrix  *
110 C  S   PYHGGM   to determine Higgs mass spectrum                     *
111 C  S   PYSUBH   to determine Higgs masses in the MSSM                *
112 C  S   PYPOLE   to determine Higgs masses in the MSSM                *
113 C  S   PYVACU   to determine Higgs masses in the MSSM                *
114 C  S   PYRGHM   auxiliary to PYVACU                                  *
115 C  S   PYGFXX   auxiliary to PYRGHM                                  *
116 C  F   PYFINT   auxiliary to PYVACU                                  *
117 C  F   PYFISB   auxiliary to PYFINT                                  *
118 C  S   PYSFDC   to calculate sfermion decay partial widths           *
119 C  S   PYGLUI   to calculate gluino decay partial widths             *
120 C  S   PYTBBN   to calculate 3-body decay of gluino to neutralino    *
121 C  S   PYTBBC   to calculate 3-body decay of gluino to chargino      *
122 C  S   PYNJDC   to calculate neutralino decay partial widths         *
123 C  S   PYCJDC   to calculate chargino decay partial widths           *
124 C  F   PYXXZ5   auxiliary for neutralino 3-body decay                *
125 C  F   PYXXW5   auxiliary for ino charge change 3-body decay         *
126 C  F   PYXXGA   auxiliary for ino -> ino + gamma decay               *
127 C  F   PYX2XG   auxiliary for ino -> ino + gauge boson decay         *
128 C  F   PYX2XH   auxiliary for ino -> ino + Higgs decay               *
129 C  F   PYXXZ2   auxiliary for chargino 3-body decay                  *
130 C  S   PYHEXT   to calculate non-SM Higgs decay partial widths       *
131 C  F   PYH2XX   auxiliary for H -> ino + ino decay                   *
132 C  F   PYGAUS   to perform Gaussian integration                      *
133 C  F   PYSIMP   to perform Simpson integration                       *
134 C  F   PYLAMF   to evaluate the lambda kinematics function           *
135 C  S   PYTBDY   to perform 3-body decay of gauginos                  *
136 C                                                                    *
137 C  S   PY1ENT   to fill one entry (= parton or particle)             *
138 C  S   PY2ENT   to fill two entries                                  *
139 C  S   PY3ENT   to fill three entries                                *
140 C  S   PY4ENT   to fill four entries                                 *
141 C  S   PYJOIN   to connect entries with colour flow information      *
142 C  S   PYGIVE   to fill (or query) commonblock variables             *
143 C  S   PYEXEC   to administrate fragmentation and decay chain        *
144 C  S   PYPREP   to rearrange showered partons along strings          *
145 C  S   PYSTRF   to do string fragmentation of jet system             *
146 C  S   PYINDF   to do independent fragmentation of one or many jets  *
147 C  S   PYDECY   to do the decay of a particle                        *
148 C  S   PYDCYK   to select parton and hadron flavours in decays       *
149 C  S   PYKFDI   to select parton and hadron flavours in fragm        *
150 C  S   PYNMES   to select number of popcorn mesons                   *
151 C  S   PYKFIN   to calculate falvour prod. ratios from input params. *
152 C  S   PYPTDI   to select transverse momenta in fragm                *
153 C  S   PYZDIS   to select longitudinal scaling variable in fragm     *
154 C  S   PYSHOW   to do timelike parton shower evolution               *
155 C  S   PYBOEI   to include Bose-Einstein effects (crudely)           *
156 C  F   PYMASS   to give the mass of a particle or parton             *
157 C  S   PYNAME   to give the name of a particle or parton             *
158 C  F   PYCHGE   to give three times the electric charge              *
159 C  F   PYCOMP   to compress standard KF flavour code to internal KC  *
160 C  S   PYERRM   to write error messages and abort faulty run         *
161 C  F   PYALEM   to give the alpha_electromagnetic value              *
162 C  F   PYALPS   to give the alpha_strong value                       *
163 C  F   PYANGL   to give the angle from known x and y components      *
164 C  F   PYR      to provide a random number generator                 *
165 C  S   PYRGET   to save the state of the random number generator     *
166 C  S   PYRSET   to set the state of the random number generator      *
167 C  S   PYROBO   to rotate and/or boost an event                      *
168 C  S   PYEDIT   to remove unwanted entries from record               *
169 C  S   PYLIST   to list event record or particle data                *
170 C  S   PYLOGO   to write a logo                                      *
171 C  S   PYUPDA   to update particle data                              *
172 C  F   PYK      to provide integer-valued event information          *
173 C  F   PYP      to provide real-valued event information             *
174 C  S   PYSPHE   to perform sphericity analysis                       *
175 C  S   PYTHRU   to perform thrust analysis                           *
176 C  S   PYCLUS   to perform three-dimensional cluster analysis        *
177 C  S   PYCELL   to perform cluster analysis in (eta, phi, E_T)       *
178 C  S   PYJMAS   to give high and low jet mass of event               *
179 C  S   PYFOWO   to give Fox-Wolfram moments                          *
180 C  S   PYTABU   to analyze events, with tabular output               *
181 C                                                                    *
182 C  S   PYEEVT   to administrate the generation of an e+e- event      *
183 C  S   PYXTEE   to give the total cross-section at given CM energy   *
184 C  S   PYRADK   to generate initial state photon radiation           *
185 C  S   PYXKFL   to select flavour of primary qqbar pair              *
186 C  S   PYXJET   to select (matrix element) jet multiplicity          *
187 C  S   PYX3JT   to select kinematics of three-jet event              *
188 C  S   PYX4JT   to select kinematics of four-jet event               *
189 C  S   PYXDIF   to select angular orientation of event               *
190 C  S   PYONIA   to perform generation of onium decay to gluons       *
191 C                                                                    *
192 C  S   PYBOOK   to book a histogram                                  *
193 C  S   PYFILL   to fill an entry in a histogram                      *
194 C  S   PYFACT   to multiply histogram contents by a factor           *
195 C  S   PYOPER   to perform operations between histograms             *
196 C  S   PYHIST   to print and reset all histograms                    *
197 C  S   PYPLOT   to print a single histogram                          *
198 C  S   PYNULL   to reset contents of a single histogram              *
199 C  S   PYDUMP   to dump histogram contents onto a file               *
200 C                                                                    *
201 C  S   PYKCUT   dummy routine for user kinematical cuts              *
202 C  S   PYEVWT   dummy routine for weighting events                   *
203 C  S   PYUPIN   dummy routine to initialize a user process           *
204 C  S   PYUPEV   dummy routine to generate a user process event       *
205 C  S   PDFSET   dummy routine to be removed when using PDFLIB        *
206 C  S   STRUCTM  dummy routine to be removed when using PDFLIB        *
207 C  S   PYTAUD   dummy routine for interface to tau decay libraries   *
208 C  S   PYTIME   dummy routine for giving date and time               *
209 C                                                                    *
210 C*********************************************************************
211
212 C...PYDATA
213 C...Default values for switches and parameters,
214 C...and particle, decay and process data.
215
216       BLOCK DATA PYDATA
217
218 C...Double precision and integer declarations.
219       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
220       INTEGER PYK,PYCHGE,PYCOMP
221 C...Commonblocks.
222       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
223       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
224       COMMON/PYDAT3/MDCY(500,3),MDME(4000,2),BRAT(4000),KFDP(4000,5)
225       COMMON/PYDAT4/CHAF(500,2)
226       CHARACTER CHAF*16
227       COMMON/PYDATR/MRPY(6),RRPY(100)
228       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
229       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
230       COMMON/PYINT1/MINT(400),VINT(400)
231       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
232       COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
233       COMMON/PYINT4/MWID(500),WIDS(500,5)
234       COMMON/PYINT5/NGENPD,NGEN(0:500,3),XSEC(0:500,3)
235       COMMON/PYINT6/PROC(0:500)
236       CHARACTER PROC*28
237       COMMON/PYINT7/SIGT(0:6,0:6,0:5)
238       COMMON/PYMSSM/IMSS(0:99),RMSS(0:99)
239       COMMON/PYSSMT/ZMIX(4,4),UMIX(2,2),VMIX(2,2),SMZ(4),SMW(2),
240      &SFMIX(16,4)
241       COMMON/PYBINS/IHIST(4),INDX(1000),BIN(20000)
242       SAVE /PYDAT1/,/PYDAT2/,/PYDAT3/,/PYDAT4/,/PYDATR/,/PYSUBS/,
243      &/PYPARS/,/PYINT1/,/PYINT2/,/PYINT3/,/PYINT4/,/PYINT5/,
244      &/PYINT6/,/PYINT7/,/PYMSSM/,/PYSSMT/,/PYBINS/
245
246 C...PYDAT1, containing status codes and most parameters.
247       DATA MSTU/
248      &   0,    0,    0, 4000,10000,  500, 4000,    0,    0,    2,
249      1   6,    1,    1,    0,    1,    1,    0,    0,    0,    0,
250      2   2,   10,    0,    0,    1,   10,    0,    0,    0,    0,
251      3   0,    0,    0,    0,    0,    0,    0,    0,    0,    0,
252      4   2,    2,    1,    4,    2,    1,    1,    0,    0,    0,
253      5  25,   24,    0,    1,    0,    0,    0,    0,    0,    0,
254      6   0,    0,    0,    0,    0,    0,    0,    0,    0,    0,
255      7  30*0,
256      1   1,    0,    0,    0,    0,    0,    0,    0,    0,    0,
257      2   1,    5,    3,    5,    0,    0,    0,    0,    0,    0,
258      &  80*0/
259       DATA PARU/
260      &  3.141592653589793D0, 6.283185307179586D0,
261      &  0.197327D0, 5.06773D0, 0.389380D0, 2.56819D0,  4*0D0,
262      1  0.001D0, 0.09D0, 0.01D0, 0D0, 0D0, 0D0, 0D0, 0D0, 0D0, 0D0,
263      2  0D0,   0D0,   0D0,   0D0,  0D0,  0D0,  0D0,  0D0,  0D0,  0D0,
264      3  0D0,   0D0,   0D0,   0D0,  0D0,  0D0,  0D0,  0D0,  0D0,  0D0,
265      4  2.0D0,  1.0D0, 0.25D0,  2.5D0, 0.05D0,
266      4  0D0,   0D0, 0.0001D0, 0D0,   0D0,
267      5  2.5D0,1.5D0,7.0D0,1.0D0,0.5D0,2.0D0,3.2D0, 0D0, 0D0, 0D0,
268      6  40*0D0,
269      &  0.00729735D0, 0.232D0, 0.007764D0, 1.0D0, 1.16639D-5,
270      &  0D0, 0D0, 0D0, 0D0,  0D0,
271      1  0.20D0, 0.25D0, 1.0D0, 4.0D0, 10D0, 0D0, 0D0,  0D0, 0D0, 0D0,
272      2 -0.693D0, -1.0D0, 0.387D0, 1.0D0, -0.08D0,
273      2 -1.0D0,  1.0D0,  1.0D0,  1.0D0,  0D0,
274      3  1.0D0,-1.0D0, 1.0D0,-1.0D0, 1.0D0,  0D0,  0D0, 0D0, 0D0, 0D0,
275      4  5.0D0, 1.0D0, 1.0D0,  0D0, 1.0D0, 1.0D0,  0D0, 0D0, 0D0, 0D0,
276      5  1.0D0, 0D0, 0D0, 0D0, 1000D0, 1.0D0, 1.0D0, 1.0D0, 1.0D0,0D0,
277      6  1.0D0, 1.0D0, 1.0D0, 1.0D0, 1.0D0,  0D0,  0D0, 0D0, 0D0, 0D0,
278      7  1.0D0, 1.0D0, 1.0D0, 1.0D0, 1.0D0, 1.0D0, 1.0D0, 0D0,0D0,0D0,
279      8  1.0D0, 1.0D0, 1.0D0, 0.0D0, 0.0D0, 1.0D0, 1.0D0, 0D0,0D0,0D0,
280      9  0D0,  0D0,  0D0,  0D0, 1.0D0,  0D0,  0D0, 0D0, 0D0, 0D0/
281       DATA MSTJ/
282      &  1,    3,    0,    0,    0,    0,    0,    0,    0,    0,
283      1  4,    2,    0,    1,    0,    0,    0,    0,    0,    0,
284      2  2,    1,    1,    2,    1,    2,    2,    0,    0,    0,
285      3  0,    0,    0,    0,    0,    0,    0,    0,    0,    0,
286      4  2,    2,    4,    2,    5,    3,    3,    0,    0,    3,
287      5  0,    3,    0,    0,    0,    0,    0,    0,    0,    0,
288      6  40*0,
289      &  5,    2,    7,    5,    1,    1,    0,    2,    0,    2,
290      1  0,    0,    0,    0,    1,    1,    0,    0,    0,    0,
291      2  80*0/
292       DATA PARJ/
293      &  0.10D0, 0.30D0, 0.40D0, 0.05D0, 0.50D0,
294      &  0.50D0, 0.50D0,   0.6D0,   1.2D0,   0.6D0,
295      1  0.50D0,0.60D0,0.75D0, 0D0, 0D0, 0D0, 0D0, 1.0D0, 1.0D0, 0D0,
296      2  0.36D0, 1.0D0,0.01D0, 2.0D0,1.0D0,0.4D0, 0D0, 0D0, 0D0, 0D0,
297      3  0.10D0, 1.0D0, 0.8D0, 1.5D0,0D0,2.0D0,0.2D0,2.5D0,0.6D0,0D0,
298      4  0.3D0, 0.58D0, 0.5D0, 0.9D0,0.5D0,1.0D0,1.0D0,1.0D0,0D0,0D0,
299      5  0.77D0, 0.77D0, 0.77D0, -0.05D0, -0.005D0,
300      5 -0.00001D0, -0.00001D0, -0.00001D0, 1.0D0, 0D0,
301      6  4.5D0, 0.7D0, 0D0,0.003D0, 0.5D0, 0.5D0, 0D0, 0D0, 0D0, 0D0,
302      7  10D0, 1000D0, 100D0, 1000D0, 0D0, 0.7D0,10D0, 0D0, 0D0, 0D0,
303      8  0.29D0, 1.0D0, 1.0D0,  0D0,  10D0, 10D0, 0D0, 0D0, 0D0, 0D0,
304      9  0.02D0, 1.0D0, 0.2D0,  0D0,  0D0,  0D0,  0D0, 0D0, 0D0, 0D0,
305      &  0D0,  0D0,  0D0,  0D0,   0D0,   0D0,  0D0,  0D0,  0D0,  0D0,
306      1  0D0,  0D0,  0D0,  0D0,   0D0,   0D0,  0D0,  0D0,  0D0,  0D0,
307      2  1.0D0, 0.25D0,91.187D0,2.489D0, 0.01D0,
308      2  2.0D0,  1.0D0, 0.25D0,0.002D0,   0D0,
309      3  0D0, 0D0, 0D0, 0D0, 0.01D0, 0.99D0, 0D0, 0D0,  0.2D0,   0D0,
310      4  60*0D0/
311
312 C...PYDAT2, with particle data and flavour treatment parameters.
313       DATA (KCHG(I,1),I=   1, 500)/-1,2,-1,2,-1,2,-1,2,2*0,-3,0,-3,0,
314      &-3,0,-3,6*0,3,9*0,3,2*0,3,0,-1,12*0,3,2*0,3,28*0,2,-1,20*0,4*3,
315      &8*0,3*3,4*0,3*3,3*0,3*3,7*0,3*3,3*0,3*3,3*0,-2,-3,2*1,3*0,4,3*3,
316      &6,2*-2,2*-3,0,2*1,2*0,2*3,-2,2*-3,2*0,-3,2*1,2*0,3,0,2*4,2*3,2*6,
317      &3,2*1,2*0,2*3,2*0,4,2*3,2*6,2*3,6,2*-2,2*-3,0,-3,0,2*1,2*0,2*3,0,
318      &3,2*-2,2*-3,2*0,2*-3,0,2*1,2*0,2*3,2*0,2*3,-2,2*-3,2*0,2*-3,2*0,
319      &-3,2*0,2*3,4*0,2*3,2*0,2*3,2*0,2*3,4*0,2*3,2*0,2*3,3*0,3,2*0,3,0,
320      &3,0,3,2*0,3,0,3,3*0,-1,2,-1,2,-1,2,-3,0,-3,0,-3,4*0,3,2*0,3,0,-1,
321      &2,-1,2,-1,2,-3,0,-3,0,-3,0,-1,2,-3,164*0/
322       DATA (KCHG(I,2),I=   1, 500)/8*1,12*0,2,16*0,2,1,113*0,-1,0,2*-1,
323      &3*0,-1,4*0,2*-1,3*0,2*-1,4*0,-1,5*0,2*-1,4*0,2*-1,5*0,2*-1,6*0,
324      &-1,7*0,2*-1,5*0,2*-1,6*0,2*-1,7*0,2*-1,8*0,-1,56*0,6*1,6*0,2,7*0,
325      &6*1,6*0,2*1,165*0/
326       DATA (KCHG(I,3),I=   1, 500)/8*1,2*0,8*1,5*0,1,9*0,1,2*0,1,0,2*1,
327      &11*0,1,2*0,1,26*0,1,0,2*1,20*0,4*1,5*0,6*1,4*0,9*1,4*0,12*1,3*0,
328      &102*1,2*0,2*1,2*0,4*1,2*0,6*1,2*0,8*1,3*0,1,0,2*1,0,3*1,0,4*1,
329      &3*0,12*1,3*0,1,2*0,1,0,16*1,163*0/
330       DATA (KCHG(I,4),I=   1, 293)/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,
331      &16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,
332      &37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,
333      &58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,
334      &79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,
335      &100,110,111,113,115,130,210,211,213,215,220,221,223,225,310,311,
336      &313,315,321,323,325,330,331,333,335,411,413,415,421,423,425,431,
337      &433,435,440,441,443,445,511,513,515,521,523,525,531,533,535,541,
338      &543,545,551,553,555,1103,1114,2101,2103,2110,2112,2114,2203,2210,
339      &2212,2214,2224,3101,3103,3112,3114,3122,3201,3203,3212,3214,3222,
340      &3224,3303,3312,3314,3322,3324,3334,4101,4103,4112,4114,4122,4132,
341      &4201,4203,4212,4214,4222,4224,4232,4301,4303,4312,4314,4322,4324,
342      &4332,4334,4403,4412,4414,4422,4424,4432,4434,4444,5101,5103,5112,
343      &5114,5122,5132,5142,5201,5203,5212,5214,5222,5224,5232,5242,5301,
344      &5303,5312,5314,5322,5324,5332,5334,5342,5401,5403,5412,5414,5422,
345      &5424,5432,5434,5442,5444,5503,5512,5514,5522,5524,5532,5534,5542,
346      &5544,5554,10111,10113,10211,10213,10221,10223,10311,10313,10321,
347      &10323,10331,10333,10411,10413,10421,10423,10431,10433,10441,
348      &10443,10511,10513,10521,10523,10531,10533,10541,10543,10551,
349      &10553,20113,20213,20223,20313,20323,20333,20413,20423,20433/
350       DATA (KCHG(I,4),I= 294, 500)/20443,20513,20523,20533,20543,20553,
351      &100443,100553,1000001,1000002,1000003,1000004,1000005,1000006,
352      &1000011,1000012,1000013,1000014,1000015,1000016,1000021,1000022,
353      &1000023,1000024,1000025,1000035,1000037,1000039,2000001,2000002,
354      &2000003,2000004,2000005,2000006,2000011,2000012,2000013,2000014,
355      &2000015,2000016,4000001,4000002,4000011,4000012,163*0/
356       DATA (PMAS(I,1),I=   1, 214)/0.0099D0,0.0056D0,0.199D0,1.35D0,
357      &5D0,175D0,2*400D0,2*0D0,0.00051D0,0D0,0.10566D0,0D0,1.777D0,0D0,
358      &400D0,5*0D0,91.187D0,80.33D0,80D0,6*0D0,500D0,900D0,500D0,
359      &3*300D0,350D0,200D0,5000D0,10*0D0,3*100D0,3*200D0,26*0D0,1D0,2D0,
360      &5D0,16*0D0,0.13498D0,0.7685D0,1.318D0,0.49767D0,0D0,0.13957D0,
361      &0.7669D0,1.318D0,0D0,0.54745D0,0.78194D0,1.275D0,2*0.49767D0,
362      &0.8961D0,1.432D0,0.4936D0,0.8916D0,1.425D0,0D0,0.95777D0,
363      &1.0194D0,1.525D0,1.8693D0,2.01D0,2.46D0,1.8645D0,2.0067D0,2.46D0,
364      &1.9685D0,2.1124D0,2.5735D0,0D0,2.9798D0,3.09688D0,3.5562D0,
365      &5.2792D0,5.3248D0,5.83D0,5.2789D0,5.3248D0,5.83D0,5.3693D0,
366      &5.4163D0,6.07D0,6.594D0,6.602D0,7.35D0,9.4D0,9.4603D0,9.9132D0,
367      &0.77133D0,1.234D0,0.57933D0,0.77133D0,0D0,0.93957D0,1.233D0,
368      &0.77133D0,0D0,0.93827D0,1.232D0,1.231D0,0.80473D0,0.92953D0,
369      &1.19744D0,1.3872D0,1.11568D0,0.80473D0,0.92953D0,1.19255D0,
370      &1.3837D0,1.18937D0,1.3828D0,1.09361D0,1.3213D0,1.535D0,1.3149D0,
371      &1.5318D0,1.67245D0,1.96908D0,2.00808D0,2.4521D0,2.5D0,2.2849D0,
372      &2.4703D0,1.96908D0,2.00808D0,2.4535D0,2.5D0,2.4529D0,2.5D0,
373      &2.4656D0,2.15432D0,2.17967D0,2.55D0,2.63D0,2.55D0,2.63D0,2.704D0,
374      &2.8D0,3.27531D0,3.59798D0,3.65648D0,3.59798D0,3.65648D0,
375      &3.78663D0,3.82466D0,4.91594D0,5.38897D0,5.40145D0,5.8D0,5.81D0/
376       DATA (PMAS(I,1),I= 215, 500)/5.641D0,5.84D0,7.00575D0,5.38897D0,
377      &5.40145D0,5.8D0,5.81D0,5.8D0,5.81D0,5.84D0,7.00575D0,5.56725D0,
378      &5.57536D0,5.96D0,5.97D0,5.96D0,5.97D0,6.12D0,6.13D0,7.19099D0,
379      &6.67143D0,6.67397D0,7.03724D0,7.0485D0,7.03724D0,7.0485D0,
380      &7.21101D0,7.219D0,8.30945D0,8.31325D0,10.07354D0,10.42272D0,
381      &10.44144D0,10.42272D0,10.44144D0,10.60209D0,10.61426D0,
382      &11.70767D0,11.71147D0,15.11061D0,0.9835D0,1.231D0,0.9835D0,
383      &1.231D0,1D0,1.17D0,1.429D0,1.29D0,1.429D0,1.29D0,2*1.4D0,2.272D0,
384      &2.424D0,2.272D0,2.424D0,2.5D0,2.536D0,3.4151D0,3.46D0,5.68D0,
385      &5.73D0,5.68D0,5.73D0,5.92D0,5.97D0,7.25D0,7.3D0,9.8598D0,9.875D0,
386      &2*1.23D0,1.282D0,2*1.402D0,1.427D0,2*2.372D0,2.56D0,3.5106D0,
387      &2*5.78D0,6.02D0,7.3D0,9.8919D0,3.686D0,10.0233D0,32*500D0,
388      &4*400D0,163*0D0/
389       DATA (PMAS(I,2),I=   1, 500)/5*0D0,1.4D0,16*0D0,2.47833D0,
390      &2.069D0,0.00295D0,6*0D0,14.67788D0,0D0,16.79392D0,8.45231D0,
391      &4.93534D0,5.80468D0,19.1898D0,0.39162D0,417.35283D0,62*0D0,
392      &0.151D0,0.107D0,3*0D0,0.149D0,0.107D0,2*0D0,0.00843D0,0.185D0,
393      &2*0D0,0.0505D0,0.109D0,0D0,0.0498D0,0.098D0,0D0,0.0002D0,
394      &0.00443D0,0.076D0,2*0D0,0.023D0,2*0D0,0.023D0,2*0D0,0.015D0,0D0,
395      &0.0013D0,0D0,0.002D0,2*0D0,0.02D0,2*0D0,0.02D0,2*0D0,0.02D0,
396      &2*0D0,0.02D0,4*0D0,0.12D0,4*0D0,0.12D0,3*0D0,2*0.12D0,3*0D0,
397      &0.0394D0,4*0D0,0.036D0,0D0,0.0358D0,2*0D0,0.0099D0,0D0,0.0091D0,
398      &74*0D0,0.06D0,0.142D0,0.06D0,0.142D0,0D0,0.36D0,0.287D0,0.09D0,
399      &0.287D0,0.09D0,0.25D0,0.08D0,0.05D0,0.02D0,0.05D0,0.02D0,0.05D0,
400      &0D0,0.014D0,0.01D0,8*0.05D0,0D0,0.01D0,2*0.4D0,0.025D0,2*0.174D0,
401      &0.053D0,3*0.05D0,0.0009D0,4*0.05D0,3*0D0,19*1D0,0D0,7*1D0,0D0,
402      &1D0,0D0,1D0,0D0,2.60511D0,2.60839D0,0.42904D0,0.41921D0,163*0D0/
403       DATA (PMAS(I,3),I=   1, 500)/5*0D0,14D0,16*0D0,24.78326D0,
404      &20.69D0,0.02954D0,6*0D0,146.77876D0,0D0,167.93924D0,84.52308D0,
405      &49.35344D0,58.04675D0,191.89803D0,3.91624D0,4173.5283D0,62*0D0,
406      &0.4D0,0.25D0,3*0D0,0.4D0,0.25D0,2*0D0,0.1D0,0.17D0,2*0D0,0.2D0,
407      &0.12D0,0D0,0.2D0,0.12D0,0D0,0.002D0,0.015D0,0.2D0,2*0D0,0.12D0,
408      &2*0D0,0.12D0,2*0D0,0.05D0,0D0,0.005D0,0D0,0.01D0,2*0D0,0.05D0,
409      &2*0D0,0.05D0,2*0D0,0.05D0,2*0D0,0.05D0,4*0D0,0.14D0,4*0D0,0.14D0,
410      &3*0D0,2*0.14D0,3*0D0,0.04D0,4*0D0,0.035D0,0D0,0.035D0,2*0D0,
411      &0.05D0,0D0,0.05D0,74*0D0,0.05D0,0.25D0,0.05D0,0.25D0,0D0,0.2D0,
412      &0.4D0,0.005D0,0.4D0,0.01D0,0.35D0,0.001D0,0.1D0,0.08D0,0.1D0,
413      &0.08D0,0.1D0,0D0,0.05D0,0.02D0,6*0.1D0,0.05D0,0.1D0,0D0,0.02D0,
414      &2*0.3D0,0.05D0,2*0.3D0,0.02D0,2*0.1D0,0.03D0,0.001D0,4*0.1D0,
415      &3*0D0,19*10D0,0.00001D0,7*10D0,0.00001D0,10D0,0.00001D0,10D0,
416      &0.00001D0,26.05109D0,26.08388D0,4.29043D0,4.19206D0,163*0D0/
417       DATA (PMAS(I,4),I=   1, 500)/12*0D0,658654D0,0D0,0.0872D0,68*0D0,
418      &0.1D0,0.387D0,16*0D0,0.00003D0,2*0D0,15500D0,0D0,7804.5D0,6*0D0,
419      &26.762D0,3*0D0,3709D0,6*0D0,0.317D0,2*0D0,0.1244D0,2*0D0,0.14D0,
420      &6*0D0,0.468D0,2*0D0,0.462D0,2*0D0,0.483D0,2*0D0,0.15D0,19*0D0,
421      &44.34D0,0D0,78.88D0,4*0D0,23.96D0,2*0D0,49.1D0,0D0,87.1D0,0D0,
422      &24.6D0,4*0D0,0.0618D0,0.029D0,6*0D0,0.106D0,6*0D0,0.019D0,2*0D0,
423      &7*0.1D0,4*0D0,0.342D0,2*0.387D0,6*0D0,2*0.387D0,6*0D0,0.387D0,
424      &0D0,0.387D0,2*0D0,8*0.387D0,0D0,9*0.387D0,83*0D0,163*0D0/
425       DATA PARF/
426      &  0.5D0,0.25D0, 0.5D0,0.25D0, 1D0, 0.5D0,  0D0,  0D0,  0D0, 0D0,
427      1  0.5D0,  0D0, 0.5D0,  0D0,  1D0,  1D0,  0D0,  0D0,  0D0, 0D0,
428      2  0.5D0,  0D0, 0.5D0,  0D0,  1D0,  1D0,  0D0,  0D0,  0D0, 0D0,
429      3  0.5D0,  0D0, 0.5D0,  0D0,  1D0,  1D0,  0D0,  0D0,  0D0, 0D0,
430      4  0.5D0,  0D0, 0.5D0,  0D0,  1D0,  1D0,  0D0,  0D0,  0D0, 0D0,
431      5  0.5D0,  0D0, 0.5D0,  0D0,  1D0,  1D0,  0D0,  0D0,  0D0, 0D0,
432      6  0.75D0, 0.5D0, 0D0,0.1667D0,0.0833D0,0.1667D0,0D0,0D0,0D0, 0D0,
433      7  0D0,  0D0,  1D0,0.3333D0,0.6667D0,0.3333D0,0D0,0D0,0D0, 0D0,
434      8  0D0,  0D0,  0D0,  0D0,  0D0,  0D0,  0D0,  0D0,  0D0, 0D0,
435      9  0D0,  0D0,  0D0,  0D0,  0D0,  0D0,  0D0,  0D0,  0D0, 0D0,
436      & 0.325D0,0.325D0,0.5D0,1.6D0, 5.0D0,  0D0,  0D0,  0D0,  0D0, 0D0,
437      1 0D0,0.11D0,0.16D0,0.048D0,0.50D0,0.45D0,0.55D0,0.60D0,0D0,0D0,
438      2 0.2D0, 0.1D0,  0D0,  0D0,  0D0,  0D0,  0D0,  0D0,  0D0, 0D0,
439      3 60*0D0,
440      4 0.2D0,  0.5D0,  8*0D0,
441      5 1800*0D0/
442       DATA ((VCKM(I,J),J=1,4),I=1,4)/
443      &  0.95113D0,  0.04884D0,  0.00003D0,  0.00000D0,
444      &  0.04884D0,  0.94940D0,  0.00176D0,  0.00000D0,
445      &  0.00003D0,  0.00176D0,  0.99821D0,  0.00000D0,
446      &  0.00000D0,  0.00000D0,  0.00000D0,  1.00000D0/
447
448 C...PYDAT3, with particle decay parameters and data.
449       DATA (MDCY(I,1),I=   1, 500)/5*0,3*1,6*0,1,0,1,5*0,3*1,6*0,1,0,
450      &7*1,10*0,2*1,0,3*1,26*0,3*1,16*0,3*1,3*0,2*1,0,7*1,0,2*1,0,12*1,
451      &0,18*1,0,1,4*0,1,3*0,2*1,2*0,3*1,2*0,4*1,0,5*1,2*0,4*1,2*0,5*1,
452      &2*0,6*1,0,7*1,2*0,5*1,2*0,6*1,2*0,7*1,2*0,8*1,0,75*1,0,7*1,0,1,0,
453      &1,0,4*1,163*0/
454       DATA (MDCY(I,2),I=   1, 500)/1,9,17,25,33,41,54,64,2*0,74,78,80,
455      &85,87,141,143,148,2*0,151,160,172,188,208,6*0,287,0,309,332,414,
456      &494,521,524,525,10*0,534,539,0,544,564,588,26*0,606,607,611,16*0,
457      &620,622,627,636,0,645,647,649,0,656,664,670,679,681,683,686,696,
458      &702,705,0,716,722,733,739,802,805,813,874,876,884,917,919,0,923,
459      &924,927,929,965,966,974,1010,1011,1019,1058,1059,1063,1094,1095,
460      &1099,1100,1109,0,1111,4*0,1112,3*0,1115,1118,2*0,1119,1121,1124,
461      &2*0,1128,1129,1132,1135,0,1138,1143,1145,1148,1150,2*0,1154,1155,
462      &1156,1232,2*0,1236,1237,1238,1239,1240,2*0,1244,1245,1247,1248,
463      &1250,1254,0,1255,1259,1263,1267,1271,1275,1279,2*0,1283,1284,
464      &1285,1302,1311,2*0,1320,1321,1322,1323,1324,1333,2*0,1342,1343,
465      &1344,1345,1346,1355,1356,2*0,1365,1374,1383,1392,1401,1410,1419,
466      &1428,0,1437,1446,1455,1464,1473,1482,1491,1500,1509,1518,1519,
467      &1520,1521,1522,1527,1530,1532,1537,1539,1544,1551,1555,1557,1559,
468      &1561,1563,1565,1567,1569,1570,1572,1574,1576,1578,1580,1582,1584,
469      &1586,1588,1589,1591,1593,1607,1609,1611,1615,1617,1619,1621,1623,
470      &1625,1627,1629,1631,1633,1644,1658,1670,1682,1694,1706,1718,1731,
471      &1742,1753,1764,1775,1786,1797,1858,1863,1965,2021,2139,2273,0,
472      &2344,2360,2376,2392,2408,2424,2440,0,2455,0,2470,0,2485,2489,
473      &2493,2496,163*0/
474       DATA (MDCY(I,3),I=   1, 500)/5*8,13,2*10,2*0,4,2,5,2,54,2,5,3,
475      &2*0,9,12,16,20,79,6*0,22,0,23,82,80,27,3,1,9,10*0,2*5,0,20,24,18,
476      &26*0,1,4,9,16*0,2,5,2*9,0,2*2,7,0,8,6,9,2*2,3,10,6,3,11,0,6,11,6,
477      &63,3,8,61,2,8,33,2,4,0,1,3,2,36,1,8,36,1,8,39,1,4,31,1,4,1,9,2,0,
478      &1,4*0,3,3*0,3,1,2*0,2,3,4,2*0,1,3*3,0,5,2,3,2,4,2*0,2*1,76,4,2*0,
479      &4*1,4,2*0,1,2,1,2,4,1,0,7*4,2*0,2*1,17,2*9,2*0,4*1,2*9,2*0,4*1,9,
480      &1,9,2*0,8*9,0,9*9,4*1,5,3,2,5,2,5,7,4,7*2,1,9*2,1,2*2,14,2*2,4,
481      &9*2,11,14,5*12,13,6*11,61,5,102,56,118,134,71,0,6*16,15,0,15,0,
482      &15,0,2*4,3,2,163*0/
483       DATA (MDME(I,1),I=   1,4000)/6*1,-1,7*1,-1,7*1,-1,7*1,-1,7*1,-1,
484      &7*1,-1,1,-1,12*1,2*-1,8*1,2*-1,73*1,-1,2*1,-1,6*1,2*-1,7*1,2*-1,
485      &3*1,-1,6*1,2*-1,6*1,2*-1,3*1,-1,3*1,-1,3*1,5*-1,3*1,-1,85*1,2*-1,
486      &6*1,8*-1,3*1,-1,3*1,-1,3*1,5*-1,3*1,4*-1,197*1,2*-1,2*1,-1,20*1,
487      &2*-1,6*1,2*-1,7*1,-1,3*1,-1,3*1,5*-1,3*1,-1,1,-1,6*1,2*-1,6*1,
488      &2*-1,1892*1,1503*0/
489       DATA (MDME(I,2),I=   1,4000)/43*102,4*0,102,0,4*53,3*102,4*0,102,
490      &2*0,3*102,4*0,102,2*0,6*102,42,6*102,2*42,2*0,8*41,2*0,36*41,
491      &8*102,0,102,0,102,2*0,21*102,8*32,8*0,16*32,21*0,62*53,8*32,14*0,
492      &16*32,27*0,62*53,18*0,62*53,9*0,18*53,3*32,0,6*32,3*0,2*32,3*0,
493      &2*32,7*0,8*32,12*0,16*32,6*0,8*32,8*0,12,2*42,2*11,9*42,0,2,3,
494      &15*0,4*42,5*0,3,12*0,2,3*0,1,0,3,16*0,2*3,15*0,2*42,2*3,18*0,2*3,
495      &3*0,1,11*0,22*42,41*0,2*3,9*0,16*42,45*0,3,10*0,10*42,20*0,2*13,
496      &6*0,12,2*0,12,0,12,14*42,16*0,48,3*13,2*42,9*0,14*42,16*0,48,
497      &3*13,2*42,9*0,14*42,19*0,48,3*13,2*42,6*0,2*11,28*42,5*0,32,3*0,
498      &4*32,2*4,0,32,45*0,14*42,52*0,10*13,2*42,2*11,4*0,2*42,2*11,6*0,
499      &2*42,2*11,0,2*42,2*11,2*42,2*11,2*42,2*11,2*42,2*11,2*42,2*11,
500      &2*42,2*11,2*42,2*11,2*0,3*42,8*0,48,3*13,20*42,4*0,18*42,4*0,
501      &9*42,0,162*42,50*0,2*12,17*0,2*32,33*0,12,9*0,32,2*0,12,11*0,
502      &4*32,2*4,5*0,828*53,1515*0/
503       DATA (BRAT(I)  ,I=   1, 418)/43*0D0,0.00003D0,0.00177D0,0.9982D0,
504      &33*0D0,1D0,6*0D0,0.1783D0,0.1735D0,0.1131D0,0.2494D0,0.003D0,
505      &0.09D0,0.0027D0,0.01D0,0.0014D0,0.0012D0,2*0.00025D0,0.0071D0,
506      &0.012D0,0.0004D0,0.00075D0,0.00006D0,2*0.00078D0,0.0034D0,0.08D0,
507      &0.011D0,0.0191D0,0.00006D0,0.005D0,0.0133D0,0.0067D0,0.0005D0,
508      &0.0035D0,0.0006D0,0.0015D0,0.00021D0,0.0002D0,0.00075D0,0.0001D0,
509      &0.0002D0,0.0011D0,3*0.0002D0,0.00022D0,0.0004D0,0.0001D0,
510      &2*0.00205D0,2*0.00069D0,0.00025D0,0.00051D0,0.00025D0,35*0D0,
511      &0.15403D0,0.11945D0,0.15402D0,0.11931D0,0.15215D0,3*0D0,
512      &0.03357D0,0.0668D0,0.03357D0,0.0668D0,0.0335D0,0.0668D0,2*0D0,
513      &0.32139D0,0.0165D0,2*0D0,0.0165D0,0.32067D0,2*0D0,0.00001D0,
514      &0.00059D0,6*0D0,2*0.10814D0,0.10806D0,3*0D0,0.00031D0,0.04438D0,
515      &0.88031D0,4*0D0,0.0002D0,0.05531D0,0D0,0.01838D0,0.00071D0,0D0,
516      &0.00009D0,0.00032D0,62*0D0,0.14449D0,0.11223D0,0.14449D0,
517      &0.11223D0,0.14443D0,0.05782D0,2*0D0,0.03172D0,0.06305D0,
518      &0.03172D0,0.06305D0,0.03172D0,0.06305D0,8*0D0,0.24928D0,0.0128D0,
519      &0.00001D0,0D0,0.0128D0,0.24882D0,0.00039D0,0D0,0.00001D0,
520      &0.00046D0,0.22153D0,5*0D0,2*0.08464D0,0.08463D0,7*0D0,0.00005D0,
521      &0.00097D0,5*0D0,0.00007D0,0D0,0.00049D0,0.00001D0,0.00006D0,
522      &0.30591D0,0.68863D0,0D0,0.0038D0,66*0D0,0.00008D0,0.00167D0/
523       DATA (BRAT(I)  ,I= 419, 722)/5*0D0,0.00013D0,0D0,0.00294D0,
524      &0.00001D0,3*0D0,0.99517D0,63*0D0,0.00002D0,0.07231D0,2*0D0,
525      &0.00001D0,0.00269D0,0D0,0.92497D0,18*0D0,0.0024D0,0.99483D0,
526      &0.00278D0,1D0,3*0.21511D0,0.21478D0,2*0D0,2*0.06995D0,2*0D0,1D0,
527      &3*0D0,0.95D0,0.05D0,3*0D0,4*0.25D0,16*0D0,4*0.25D0,20*0D0,1D0,
528      &17*0D0,1D0,2*0.08D0,0.76D0,0.08D0,2*0.105D0,0.04D0,0.5D0,0.08D0,
529      &0.14D0,0.01D0,0.015D0,0.005D0,0.988D0,0.012D0,0.998739D0,
530      &0.00079D0,0.00038D0,0.000046D0,0.000045D0,2*0.34725D0,0.144D0,
531      &0.104D0,0.0245D0,2*0.01225D0,0.0028D0,0.0057D0,0.2112D0,0.1256D0,
532      &2*0.1939D0,2*0.1359D0,0.002D0,0.001D0,0.0006D0,0.999877D0,
533      &0.000123D0,0.99955D0,0.00045D0,2*0.34725D0,0.144D0,0.104D0,
534      &0.049D0,0.0028D0,0.0057D0,0.3923D0,0.321D0,0.2317D0,0.0478D0,
535      &0.0049D0,0.0013D0,0.0003D0,0.0007D0,0.89D0,0.08693D0,0.0221D0,
536      &0.00083D0,2*0.00007D0,0.564D0,0.282D0,0.072D0,0.028D0,0.023D0,
537      &2*0.0115D0,0.005D0,0.003D0,0.6861D0,0.3139D0,2*0.5D0,0.665D0,
538      &0.333D0,0.002D0,0.333D0,0.166D0,0.168D0,0.084D0,0.087D0,0.043D0,
539      &0.059D0,2*0.029D0,0.002D0,0.6352D0,0.2116D0,0.0559D0,0.0173D0,
540      &0.0482D0,0.0318D0,0.666D0,0.333D0,0.001D0,0.332D0,0.166D0,
541      &0.168D0,0.084D0,0.086D0,0.043D0,0.059D0,2*0.029D0,2*0.002D0,
542      &0.437D0,0.208D0,0.302D0,0.0302D0,0.0212D0,0.0016D0,0.48947D0/
543       DATA (BRAT(I)  ,I= 723, 897)/0.34D0,3*0.043D0,0.027D0,0.0126D0,
544      &0.0013D0,0.0003D0,0.00025D0,0.00008D0,0.444D0,2*0.222D0,0.104D0,
545      &2*0.004D0,0.07D0,0.065D0,2*0.005D0,2*0.011D0,5*0.001D0,0.07D0,
546      &0.065D0,2*0.005D0,2*0.011D0,5*0.001D0,0.026D0,0.019D0,0.066D0,
547      &0.041D0,0.045D0,0.076D0,0.0073D0,2*0.0047D0,0.026D0,0.001D0,
548      &0.0006D0,0.0066D0,0.005D0,2*0.003D0,2*0.0006D0,2*0.001D0,0.006D0,
549      &0.005D0,0.012D0,0.0057D0,0.067D0,0.008D0,0.0022D0,0.027D0,
550      &0.004D0,0.019D0,0.012D0,0.002D0,0.009D0,0.0218D0,0.001D0,0.022D0,
551      &0.087D0,0.001D0,0.0019D0,0.0015D0,0.0028D0,0.683D0,0.306D0,
552      &0.011D0,0.3D0,0.15D0,0.16D0,0.08D0,0.13D0,0.06D0,0.08D0,0.04D0,
553      &0.034D0,0.027D0,2*0.002D0,2*0.004D0,2*0.002D0,0.034D0,0.027D0,
554      &2*0.002D0,2*0.004D0,2*0.002D0,0.0365D0,0.045D0,0.073D0,0.062D0,
555      &3*0.021D0,0.0061D0,0.015D0,0.025D0,0.0088D0,0.074D0,0.0109D0,
556      &0.0041D0,0.002D0,0.0035D0,0.0011D0,0.001D0,0.0027D0,2*0.0016D0,
557      &0.0018D0,0.011D0,0.0063D0,0.0052D0,0.018D0,0.016D0,0.0034D0,
558      &0.0036D0,0.0009D0,0.0006D0,0.015D0,0.0923D0,0.018D0,0.022D0,
559      &0.0077D0,0.009D0,0.0075D0,0.024D0,0.0085D0,0.067D0,0.0511D0,
560      &0.017D0,0.0004D0,0.0028D0,0.619D0,0.381D0,0.3D0,0.15D0,0.16D0,
561      &0.08D0,0.13D0,0.06D0,0.08D0,0.04D0,0.01D0,2*0.02D0,0.03D0,
562      &2*0.005D0,2*0.02D0,0.03D0,2*0.005D0,0.015D0,0.037D0,0.028D0/
563       DATA (BRAT(I)  ,I= 898,1063)/0.079D0,0.095D0,0.052D0,0.0078D0,
564      &4*0.001D0,0.028D0,0.033D0,0.026D0,0.05D0,0.01D0,4*0.005D0,0.25D0,
565      &0.0952D0,0.94D0,0.06D0,2*0.4D0,2*0.1D0,1D0,0.0602D0,0.0601D0,
566      &0.8797D0,0.135D0,0.865D0,0.02D0,0.055D0,2*0.005D0,0.008D0,
567      &0.012D0,0.02D0,0.055D0,2*0.005D0,0.008D0,0.012D0,0.01D0,0.03D0,
568      &0.0035D0,0.011D0,0.0055D0,0.0042D0,0.009D0,0.018D0,0.015D0,
569      &0.0185D0,0.0135D0,0.025D0,0.0004D0,0.0007D0,0.0008D0,0.0014D0,
570      &0.0019D0,0.0025D0,0.4291D0,0.08D0,0.07D0,0.02D0,0.015D0,0.005D0,
571      &1D0,0.3D0,0.15D0,0.16D0,0.08D0,0.13D0,0.06D0,0.08D0,0.04D0,
572      &0.02D0,0.055D0,2*0.005D0,0.008D0,0.012D0,0.02D0,0.055D0,
573      &2*0.005D0,0.008D0,0.012D0,0.01D0,0.03D0,0.0035D0,0.011D0,
574      &0.0055D0,0.0042D0,0.009D0,0.018D0,0.015D0,0.0185D0,0.0135D0,
575      &0.025D0,0.0004D0,0.0007D0,0.0008D0,0.0014D0,0.0019D0,0.0025D0,
576      &0.4291D0,0.08D0,0.07D0,0.02D0,0.015D0,0.005D0,1D0,0.3D0,0.15D0,
577      &0.16D0,0.08D0,0.13D0,0.06D0,0.08D0,0.04D0,0.02D0,0.055D0,
578      &2*0.005D0,0.008D0,0.012D0,0.02D0,0.055D0,2*0.005D0,0.008D0,
579      &0.012D0,0.01D0,0.03D0,0.0035D0,0.011D0,0.0055D0,0.0042D0,0.009D0,
580      &0.018D0,0.015D0,0.0185D0,0.0135D0,0.025D0,2*0.0002D0,0.0007D0,
581      &2*0.0004D0,0.0014D0,0.001D0,0.0009D0,0.0025D0,0.4291D0,0.08D0,
582      &0.07D0,0.02D0,0.015D0,0.005D0,1D0,2*0.3D0,2*0.2D0,0.047D0/
583       DATA (BRAT(I)  ,I=1064,1254)/0.122D0,0.006D0,0.012D0,0.035D0,
584      &0.012D0,0.035D0,0.003D0,0.007D0,0.15D0,0.037D0,0.008D0,0.002D0,
585      &0.05D0,0.015D0,0.003D0,0.001D0,0.014D0,0.042D0,0.014D0,0.042D0,
586      &0.24D0,0.065D0,0.012D0,0.003D0,0.001D0,0.002D0,0.001D0,0.002D0,
587      &0.014D0,0.003D0,1D0,2*0.3D0,2*0.2D0,1D0,0.0252D0,0.0248D0,
588      &0.0267D0,0.015D0,0.045D0,0.015D0,0.045D0,0.7743D0,0.029D0,0.22D0,
589      &0.78D0,1D0,0.331D0,0.663D0,0.006D0,0.663D0,0.331D0,0.006D0,1D0,
590      &0.999D0,0.001D0,0.88D0,2*0.06D0,0.639D0,0.358D0,0.002D0,0.001D0,
591      &1D0,0.88D0,2*0.06D0,0.516D0,0.483D0,0.001D0,0.88D0,2*0.06D0,
592      &0.9988D0,0.0001D0,0.0006D0,0.0004D0,0.0001D0,0.667D0,0.333D0,
593      &0.9954D0,0.0011D0,0.0035D0,0.333D0,0.667D0,0.676D0,0.234D0,
594      &0.085D0,0.005D0,2*1D0,0.018D0,2*0.005D0,0.003D0,0.002D0,
595      &2*0.006D0,0.018D0,2*0.005D0,0.003D0,0.002D0,2*0.006D0,0.0066D0,
596      &0.025D0,0.016D0,0.0088D0,2*0.005D0,0.0058D0,0.005D0,0.0055D0,
597      &4*0.004D0,2*0.002D0,2*0.004D0,0.003D0,0.002D0,2*0.003D0,
598      &3*0.002D0,2*0.001D0,0.002D0,2*0.001D0,2*0.002D0,0.0013D0,
599      &0.0018D0,5*0.001D0,4*0.003D0,2*0.005D0,2*0.002D0,2*0.001D0,
600      &2*0.002D0,2*0.001D0,0.2432D0,0.057D0,2*0.035D0,0.15D0,2*0.075D0,
601      &0.03D0,2*0.015D0,2*0.08D0,0.76D0,0.08D0,4*1D0,2*0.08D0,0.76D0,
602      &0.08D0,1D0,2*0.5D0,1D0,2*0.5D0,2*0.08D0,0.76D0,0.08D0,1D0/
603       DATA (BRAT(I)  ,I=1255,1447)/2*0.08D0,0.76D0,3*0.08D0,0.76D0,
604      &3*0.08D0,0.76D0,3*0.08D0,0.76D0,3*0.08D0,0.76D0,3*0.08D0,0.76D0,
605      &3*0.08D0,0.76D0,0.08D0,2*1D0,2*0.105D0,0.04D0,0.0077D0,0.02D0,
606      &0.0235D0,0.0285D0,0.0435D0,0.0011D0,0.0022D0,0.0044D0,0.4291D0,
607      &0.08D0,0.07D0,0.02D0,0.015D0,0.005D0,2*0.105D0,0.04D0,0.5D0,
608      &0.08D0,0.14D0,0.01D0,0.015D0,0.005D0,2*0.105D0,0.04D0,0.5D0,
609      &0.08D0,0.14D0,0.01D0,0.015D0,0.005D0,4*1D0,2*0.105D0,0.04D0,
610      &0.5D0,0.08D0,0.14D0,0.01D0,0.015D0,0.005D0,2*0.105D0,0.04D0,
611      &0.5D0,0.08D0,0.14D0,0.01D0,0.015D0,0.005D0,4*1D0,2*0.105D0,
612      &0.04D0,0.5D0,0.08D0,0.14D0,0.01D0,0.015D0,0.005D0,1D0,2*0.105D0,
613      &0.04D0,0.5D0,0.08D0,0.14D0,0.01D0,0.015D0,0.005D0,2*0.105D0,
614      &0.04D0,0.5D0,0.08D0,0.14D0,0.01D0,0.015D0,0.005D0,2*0.105D0,
615      &0.04D0,0.5D0,0.08D0,0.14D0,0.01D0,0.015D0,0.005D0,2*0.105D0,
616      &0.04D0,0.5D0,0.08D0,0.14D0,0.01D0,0.015D0,0.005D0,2*0.105D0,
617      &0.04D0,0.5D0,0.08D0,0.14D0,0.01D0,0.015D0,0.005D0,2*0.105D0,
618      &0.04D0,0.5D0,0.08D0,0.14D0,0.01D0,0.015D0,0.005D0,2*0.105D0,
619      &0.04D0,0.5D0,0.08D0,0.14D0,0.01D0,0.015D0,0.005D0,2*0.105D0,
620      &0.04D0,0.5D0,0.08D0,0.14D0,0.01D0,0.015D0,0.005D0,2*0.105D0,
621      &0.04D0,0.5D0,0.08D0,0.14D0,0.01D0,0.015D0,0.005D0,2*0.105D0,
622      &0.04D0,0.5D0,0.08D0,0.14D0,0.01D0,0.015D0,0.005D0,2*0.105D0/
623       DATA (BRAT(I)  ,I=1448,1648)/0.04D0,0.5D0,0.08D0,0.14D0,0.01D0,
624      &0.015D0,0.005D0,2*0.105D0,0.04D0,0.5D0,0.08D0,0.14D0,0.01D0,
625      &0.015D0,0.005D0,2*0.105D0,0.04D0,0.5D0,0.08D0,0.14D0,0.01D0,
626      &0.015D0,0.005D0,2*0.105D0,0.04D0,0.5D0,0.08D0,0.14D0,0.01D0,
627      &0.015D0,0.005D0,2*0.105D0,0.04D0,0.5D0,0.08D0,0.14D0,0.01D0,
628      &0.015D0,0.005D0,2*0.105D0,0.04D0,0.5D0,0.08D0,0.14D0,0.01D0,
629      &0.015D0,0.005D0,2*0.105D0,0.04D0,0.5D0,0.08D0,0.14D0,0.01D0,
630      &0.015D0,0.005D0,2*0.105D0,0.04D0,0.5D0,0.08D0,0.14D0,0.01D0,
631      &0.015D0,0.005D0,4*1D0,0.52D0,0.26D0,0.11D0,2*0.055D0,0.333D0,
632      &0.334D0,0.333D0,0.667D0,0.333D0,0.28D0,0.14D0,0.313D0,0.157D0,
633      &0.11D0,0.667D0,0.333D0,0.28D0,0.14D0,0.313D0,0.157D0,0.11D0,
634      &0.36D0,0.18D0,0.03D0,2*0.015D0,2*0.2D0,4*0.25D0,0.667D0,0.333D0,
635      &0.667D0,0.333D0,0.667D0,0.333D0,0.667D0,0.333D0,4*0.5D0,0.007D0,
636      &0.993D0,1D0,0.667D0,0.333D0,0.667D0,0.333D0,0.667D0,0.333D0,
637      &0.667D0,0.333D0,8*0.5D0,0.02D0,0.98D0,1D0,4*0.5D0,3*0.146D0,
638      &3*0.05D0,0.15D0,2*0.05D0,4*0.024D0,0.066D0,0.667D0,0.333D0,
639      &0.667D0,0.333D0,4*0.25D0,0.667D0,0.333D0,0.667D0,0.333D0,2*0.5D0,
640      &0.273D0,0.727D0,0.667D0,0.333D0,0.667D0,0.333D0,4*0.5D0,0.35D0,
641      &0.65D0,2*0.0083D0,0.1866D0,0.324D0,0.184D0,0.027D0,0.001D0,
642      &0.093D0,0.087D0,0.078D0,0.0028D0,3*0.014D0,0.008D0,0.024D0/
643       DATA (BRAT(I)  ,I=1649,4000)/0.008D0,0.024D0,0.425D0,0.02D0,
644      &0.185D0,0.088D0,0.043D0,0.067D0,0.066D0,827*0D0,0.8516D0,
645      &0.00539D0,0.04483D0,0.09819D0,0.85053D0,0.02152D0,0.02989D0,
646      &0.09806D0,0.29439D0,0.10943D0,0.59618D0,0.38983D0,0.61017D0,
647      &1503*0D0/
648       DATA (KFDP(I,1),I=   1, 375)/21,22,23,4*-24,25,21,22,23,4*24,25,
649      &21,22,23,4*-24,25,21,22,23,4*24,25,21,22,23,4*-24,25,21,22,23,
650      &4*24,25,37,1000022,1000023,1000025,1000035,21,22,23,4*-24,25,
651      &2*-37,21,22,23,4*24,25,2*37,22,23,-24,25,23,24,-12,22,23,-24,25,
652      &23,24,-12,-14,48*16,22,23,-24,25,23,24,22,23,-24,25,-37,23,24,37,
653      &1,2,3,4,5,6,7,8,21,1,2,3,4,5,6,7,8,11,13,15,17,1,2,3,4,5,6,7,8,
654      &11,12,13,14,15,16,17,18,4*-1,4*-3,4*-5,4*-7,-11,-13,-15,-17,1,2,
655      &3,4,5,6,7,8,11,13,15,17,21,2*22,23,24,1000022,2*1000023,
656      &3*1000025,4*1000035,2*1000024,2*1000037,1000001,2000001,1000001,
657      &-1000001,1000002,2000002,1000002,-1000002,1000003,2000003,
658      &1000003,-1000003,1000004,2000004,1000004,-1000004,1000005,
659      &2000005,1000005,-1000005,1000006,2000006,1000006,-1000006,
660      &1000011,2000011,1000011,-1000011,1000012,2000012,1000012,
661      &-1000012,1000013,2000013,1000013,-1000013,1000014,2000014,
662      &1000014,-1000014,1000015,2000015,1000015,-1000015,1000016,
663      &2000016,1000016,-1000016,1,2,3,4,5,6,7,8,11,12,13,14,15,16,17,18,
664      &24,37,2*23,25,35,4*-1,4*-3,4*-5,4*-7,-11,-13,-15,-17,3*24,1,2,3,
665      &4,5,6,7,8,11,13,15,17,21,2*22,23,24,23,25,36,1000022,2*1000023,
666      &3*1000025,4*1000035,2*1000024,2*1000037,1000001,2000001,1000001,
667      &-1000001,1000002,2000002,1000002,-1000002,1000003,2000003/
668       DATA (KFDP(I,1),I= 376, 606)/1000003,-1000003,1000004,2000004,
669      &1000004,-1000004,1000005,2000005,1000005,-1000005,1000006,
670      &2000006,1000006,-1000006,1000011,2000011,1000011,-1000011,
671      &1000012,2000012,1000012,-1000012,1000013,2000013,1000013,
672      &-1000013,1000014,2000014,1000014,-1000014,1000015,2000015,
673      &1000015,-1000015,1000016,2000016,1000016,-1000016,1,2,3,4,5,6,7,
674      &8,11,13,15,17,21,2*22,23,24,23,1000022,2*1000023,3*1000025,
675      &4*1000035,2*1000024,2*1000037,1000001,2000001,1000001,-1000001,
676      &1000002,2000002,1000002,-1000002,1000003,2000003,1000003,
677      &-1000003,1000004,2000004,1000004,-1000004,1000005,2000005,
678      &1000005,-1000005,1000006,2000006,1000006,-1000006,1000011,
679      &2000011,1000011,-1000011,1000012,2000012,1000012,-1000012,
680      &1000013,2000013,1000013,-1000013,1000014,2000014,1000014,
681      &-1000014,1000015,2000015,1000015,-1000015,1000016,2000016,
682      &1000016,-1000016,-1,-3,-5,-7,-11,-13,-15,-17,24,2*1000022,
683      &2*1000023,2*1000025,2*1000035,1000006,2000006,1000006,2000006,
684      &-1000001,-1000003,-1000011,-1000013,-1000015,-2000015,5,6,21,2,1,
685      &2,3,4,5,6,11,13,15,4,5,11,13,15,2*4,-11,-13,-15,2*24,2*52,1,2,3,
686      &4,5,6,7,8,11,12,13,14,15,16,17,18,2*24,2*52,4*-1,4*-3,4*-5,4*-7,
687      &-11,-13,-15,-17,22,23,1,2,3,4,5,6,7,8,11,12,13,14,15,16,17,18,82/
688       DATA (KFDP(I,1),I= 607,1001)/-11,-13,2*2,-12,-14,-16,2*-2,2*-4,
689      &-2,-4,2*22,211,111,221,13,11,213,-213,221,223,321,130,310,111,
690      &331,111,211,-12,12,-14,14,211,111,22,-13,-11,2*211,213,113,221,
691      &223,321,211,331,22,111,211,2*22,211,22,111,211,22,211,221,111,11,
692      &211,111,2*211,321,130,310,221,111,211,111,130,310,321,2*311,321,
693      &311,323,313,323,313,321,3*311,-13,3*211,12,14,311,2*321,311,321,
694      &313,323,313,323,311,4*321,211,111,3*22,111,321,130,-213,113,213,
695      &211,22,111,11,13,211,321,130,310,221,211,111,11*-11,11*-13,-311,
696      &-313,-311,-313,-20313,2*-311,-313,-311,-313,2*111,2*221,2*331,
697      &2*113,2*223,2*333,-311,-313,2*-321,211,-311,-321,333,-311,-313,
698      &-321,211,2*-321,2*-311,-321,211,113,421,2*411,421,411,423,413,
699      &423,413,421,411,8*-11,8*-13,-321,-323,-321,-323,-311,2*-313,-311,
700      &-313,2*-311,-321,-10323,-321,-323,-321,-311,2*-313,211,111,333,
701      &3*-321,-311,-313,-321,-313,310,333,211,2*-321,-311,-313,-311,211,
702      &-321,3*-311,211,113,321,2*421,411,421,413,423,413,423,411,421,
703      &-15,5*-11,5*-13,221,331,333,221,331,333,10221,211,213,211,213,
704      &321,323,321,323,2212,221,331,333,221,2*2,2*431,421,411,423,413,
705      &82,11,13,82,443,82,6*12,6*14,2*16,3*-411,3*-413,2*-411,2*-413,
706      &2*441,2*443,2*20443,2*2,2*4,2,4,511,521,511,523,513,523,513,521,
707      &511,6*12,6*14,2*16,3*-421,3*-423,2*-421,2*-423,2*441,2*443/
708       DATA (KFDP(I,1),I=1002,1428)/2*20443,2*2,2*4,2,4,521,511,521,513,
709      &523,513,523,511,521,6*12,6*14,2*16,3*-431,3*-433,2*-431,2*-433,
710      &3*441,3*443,3*20443,2*2,2*4,2,4,531,521,511,523,513,16,2*4,2*12,
711      &2*14,2*16,4*2,4*4,2*-11,2*-13,2*-1,2*-3,2*-11,2*-13,2*-1,541,511,
712      &521,513,523,21,11,13,15,1,2,3,4,21,22,553,21,2112,2212,2*2112,
713      &2212,2112,2*2212,2112,-12,3122,3212,3112,2212,2*2112,-12,2*3122,
714      &3222,3112,2212,2112,2212,3122,3222,3212,3122,3112,-12,-14,-12,
715      &3322,3312,2*3122,3212,3322,3312,3122,3322,3312,-12,2*4122,7*-11,
716      &7*-13,2*2224,2*2212,2*2214,2*3122,2*3212,2*3214,5*3222,4*3224,
717      &2*3322,3324,2*2224,7*2212,5*2214,2*2112,2*2114,2*3122,2*3212,
718      &2*3214,2*3222,2*3224,4*2,3,2*2,1,2*2,-11,-13,2*2,4*4122,-11,-13,
719      &2*2,3*4132,3*4232,-11,-13,2*2,4332,-11,-13,2*2,-11,-13,2*2,-11,
720      &-13,2*2,-11,-13,2*2,-11,-13,2*2,-11,-13,2*2,-11,-13,2*2,2*5122,
721      &-12,-14,-16,5*4122,441,443,20443,2*-2,2*-4,-2,-4,-12,-14,-16,
722      &2*-2,2*-4,-2,-4,-12,-14,-16,2*-2,2*-4,-2,-4,4*5122,-12,-14,-16,
723      &2*-2,2*-4,-2,-4,-12,-14,-16,2*-2,2*-4,-2,-4,2*5132,2*5232,-12,
724      &-14,-16,2*-2,2*-4,-2,-4,5332,-12,-14,-16,2*-2,2*-4,-2,-4,-12,-14,
725      &-16,2*-2,2*-4,-2,-4,-12,-14,-16,2*-2,2*-4,-2,-4,-12,-14,-16,2*-2,
726      &2*-4,-2,-4,-12,-14,-16,2*-2,2*-4,-2,-4,-12,-14,-16,2*-2,2*-4,-2,
727      &-4,-12,-14,-16,2*-2,2*-4,-2,-4,-12,-14,-16,2*-2,2*-4,-2,-4,-12/
728       DATA (KFDP(I,1),I=1429,1710)/-14,-16,2*-2,2*-4,-2,-4,-12,-14,-16,
729      &2*-2,2*-4,-2,-4,-12,-14,-16,2*-2,2*-4,-2,-4,-12,-14,-16,2*-2,
730      &2*-4,-2,-4,-12,-14,-16,2*-2,2*-4,-2,-4,-12,-14,-16,2*-2,2*-4,-2,
731      &-4,-12,-14,-16,2*-2,2*-4,-2,-4,-12,-14,-16,2*-2,2*-4,-2,-4,-12,
732      &-14,-16,2*-2,2*-4,-2,-4,-12,-14,-16,2*-2,2*-4,-2,-4,221,223,221,
733      &223,211,111,321,130,310,213,113,-213,321,311,321,311,323,313,
734      &2*311,321,311,321,313,323,321,211,111,321,130,310,2*211,313,-313,
735      &323,-323,421,411,423,413,411,421,413,423,411,421,423,413,443,
736      &2*82,521,511,523,513,511,521,513,523,521,511,523,513,511,521,513,
737      &523,553,2*21,213,-213,113,213,10211,10111,-10211,2*221,213,2*113,
738      &-213,2*321,2*311,113,323,2*313,323,313,-313,323,-323,423,2*413,
739      &2*423,413,443,82,523,2*513,2*523,2*513,523,553,21,11,13,82,4*443,
740      &10441,20443,445,441,11,13,15,1,2,3,4,21,22,2*553,10551,20553,555,
741      &1000039,-1000024,-1000037,1000022,1000023,1000025,1000035,
742      &1000002,2000002,1000002,2000002,1000021,1000039,1000024,1000037,
743      &1000022,1000023,1000025,1000035,1000001,2000001,1000001,2000001,
744      &1000021,1000039,-1000024,-1000037,1000022,1000023,1000025,
745      &1000035,1000004,2000004,1000004,2000004,1000021,1000039,1000024,
746      &1000037,1000022,1000023,1000025,1000035,1000003,2000003,1000003,
747      &2000003,1000021,1000039,-1000024,-1000037,1000022,1000023/
748       DATA (KFDP(I,1),I=1711,1900)/1000025,1000035,1000006,2000006,
749      &1000006,2000006,1000021,1000039,1000024,1000037,1000022,1000023,
750      &1000025,1000035,1000005,2000005,1000005,2000005,1000021,1000022,
751      &1000039,-1000024,-1000037,1000022,1000023,1000025,1000035,
752      &1000012,2000012,1000012,2000012,1000039,1000024,1000037,1000022,
753      &1000023,1000025,1000035,1000011,2000011,1000011,2000011,1000039,
754      &-1000024,-1000037,1000022,1000023,1000025,1000035,1000014,
755      &2000014,1000014,2000014,1000039,1000024,1000037,1000022,1000023,
756      &1000025,1000035,1000013,2000013,1000013,2000013,1000039,-1000024,
757      &-1000037,1000022,1000023,1000025,1000035,1000016,2000016,1000016,
758      &2000016,1000039,1000024,1000037,1000022,1000023,1000025,1000035,
759      &1000015,2000015,1000015,2000015,1000039,1000001,-1000001,2000001,
760      &-2000001,1000002,-1000002,2000002,-2000002,1000003,-1000003,
761      &2000003,-2000003,1000004,-1000004,2000004,-2000004,1000005,
762      &-1000005,2000005,-2000005,1000006,-1000006,2000006,-2000006,
763      &6*1000022,6*1000023,6*1000025,6*1000035,1000024,-1000024,1000024,
764      &-1000024,1000024,-1000024,1000037,-1000037,1000037,-1000037,
765      &1000037,-1000037,10*1000039,16*1000022,1000024,-1000024,1000024,
766      &-1000024,1000024,-1000024,1000024,-1000024,1000024,-1000024,
767      &1000024,-1000024,1000037,-1000037,1000037,-1000037,1000037/
768       DATA (KFDP(I,1),I=1901,2095)/-1000037,1000037,-1000037,1000037,
769      &-1000037,1000037,-1000037,1000024,-1000024,1000037,-1000037,
770      &1000001,-1000001,2000001,-2000001,1000002,-1000002,2000002,
771      &-2000002,1000003,-1000003,2000003,-2000003,1000004,-1000004,
772      &2000004,-2000004,1000005,-1000005,2000005,-2000005,1000006,
773      &-1000006,2000006,-2000006,1000011,-1000011,2000011,-2000011,
774      &1000012,-1000012,2000012,-2000012,1000013,-1000013,2000013,
775      &-2000013,1000014,-1000014,2000014,-2000014,1000015,-1000015,
776      &2000015,-2000015,1000016,-1000016,2000016,-2000016,5*1000021,
777      &2*1000039,6*1000022,6*1000023,6*1000025,6*1000035,1000022,
778      &1000023,1000025,1000035,1000002,2000002,-1000001,-2000001,
779      &1000004,2000004,-1000003,-2000003,1000006,2000006,-1000005,
780      &-2000005,1000012,2000012,-1000011,-2000011,1000014,2000014,
781      &-1000013,-2000013,1000016,2000016,-1000015,-2000015,2*1000021,
782      &5*1000039,16*1000022,16*1000023,1000024,-1000024,1000024,
783      &-1000024,1000024,-1000024,1000024,-1000024,1000024,-1000024,
784      &1000024,-1000024,1000037,-1000037,1000037,-1000037,1000037,
785      &-1000037,1000037,-1000037,1000037,-1000037,1000037,-1000037,
786      &1000024,-1000024,1000037,-1000037,1000001,-1000001,2000001,
787      &-2000001,1000002,-1000002,2000002,-2000002,1000003,-1000003/
788       DATA (KFDP(I,1),I=2096,2323)/2000003,-2000003,1000004,-1000004,
789      &2000004,-2000004,1000005,-1000005,2000005,-2000005,1000006,
790      &-1000006,2000006,-2000006,1000011,-1000011,2000011,-2000011,
791      &1000012,-1000012,2000012,-2000012,1000013,-1000013,2000013,
792      &-2000013,1000014,-1000014,2000014,-2000014,1000015,-1000015,
793      &2000015,-2000015,1000016,-1000016,2000016,-2000016,5*1000021,
794      &5*1000039,16*1000022,16*1000023,16*1000025,1000024,-1000024,
795      &1000024,-1000024,1000024,-1000024,1000024,-1000024,1000024,
796      &-1000024,1000024,-1000024,1000037,-1000037,1000037,-1000037,
797      &1000037,-1000037,1000037,-1000037,1000037,-1000037,1000037,
798      &-1000037,1000024,-1000024,1000037,-1000037,1000001,-1000001,
799      &2000001,-2000001,1000002,-1000002,2000002,-2000002,1000003,
800      &-1000003,2000003,-2000003,1000004,-1000004,2000004,-2000004,
801      &1000005,-1000005,2000005,-2000005,1000006,-1000006,2000006,
802      &-2000006,1000011,-1000011,2000011,-2000011,1000012,-1000012,
803      &2000012,-2000012,1000013,-1000013,2000013,-2000013,1000014,
804      &-1000014,2000014,-2000014,1000015,-1000015,2000015,-2000015,
805      &1000016,-1000016,2000016,-2000016,5*1000021,2*1000039,15*1000024,
806      &6*1000022,6*1000023,6*1000025,6*1000035,1000022,1000023,1000025,
807      &1000035,1000002,2000002,-1000001,-2000001,1000004,2000004/
808       DATA (KFDP(I,1),I=2324,4000)/-1000003,-2000003,1000006,2000006,
809      &-1000005,-2000005,1000012,2000012,-1000011,-2000011,1000014,
810      &2000014,-1000013,-2000013,1000016,2000016,-1000015,-2000015,
811      &2*1000021,1000039,-1000024,-1000037,1000022,1000023,1000025,
812      &1000035,4*1000001,1000002,2000002,1000002,2000002,1000021,
813      &1000039,1000024,1000037,1000022,1000023,1000025,1000035,
814      &4*1000002,1000001,2000001,1000001,2000001,1000021,1000039,
815      &-1000024,-1000037,1000022,1000023,1000025,1000035,4*1000003,
816      &1000004,2000004,1000004,2000004,1000021,1000039,1000024,1000037,
817      &1000022,1000023,1000025,1000035,4*1000004,1000003,2000003,
818      &1000003,2000003,1000021,1000039,-1000024,-1000037,1000022,
819      &1000023,1000025,1000035,4*1000005,1000006,2000006,1000006,
820      &2000006,1000021,1000039,1000024,1000037,1000022,1000023,1000025,
821      &1000035,4*1000006,1000005,2000005,1000005,2000005,1000021,
822      &1000039,-1000024,-1000037,1000022,1000023,1000025,1000035,
823      &4*1000011,1000012,2000012,1000012,2000012,1000039,-1000024,
824      &-1000037,1000022,1000023,1000025,1000035,4*1000013,1000014,
825      &2000014,1000014,2000014,1000039,-1000024,-1000037,1000022,
826      &1000023,1000025,1000035,4*1000015,1000016,2000016,1000016,
827      &2000016,21,22,23,-24,21,22,23,24,22,23,-24,23,24,1503*0/
828       DATA (KFDP(I,2),I=   1, 337)/3*1,2,4,6,8,1,3*2,1,3,5,7,2,3*3,2,4,
829      &6,8,3,3*4,1,3,5,7,4,3*5,2,4,6,8,5,3*6,1,3,5,7,6,5,4*1000006,3*7,
830      &2,4,6,8,7,4,6,3*8,1,3,5,7,8,5,7,2*11,12,11,12,2*11,2*13,14,13,14,
831      &13,11,13,-211,-213,-211,-213,-211,-213,-211,-213,2*-211,-321,
832      &-323,-321,2*-323,3*-321,4*-211,-213,-211,-213,-211,-213,-211,
833      &-213,-211,-213,3*-211,-213,4*-211,-323,-321,2*-211,2*-321,3*-211,
834      &2*15,16,15,16,15,2*17,18,17,2*18,2*17,-1,-2,-3,-4,-5,-6,-7,-8,21,
835      &-1,-2,-3,-4,-5,-6,-7,-8,-11,-13,-15,-17,-1,-2,-3,-4,-5,-6,-7,-8,
836      &-11,-12,-13,-14,-15,-16,-17,-18,2,4,6,8,2,4,6,8,2,4,6,8,2,4,6,8,
837      &12,14,16,18,-1,-2,-3,-4,-5,-6,-7,-8,-11,-13,-15,-17,21,22,2*23,
838      &-24,2*1000022,1000023,1000022,1000023,1000025,1000022,1000023,
839      &1000025,1000035,-1000024,-1000037,-1000024,-1000037,-1000001,
840      &2*-2000001,2000001,-1000002,2*-2000002,2000002,-1000003,
841      &2*-2000003,2000003,-1000004,2*-2000004,2000004,-1000005,
842      &2*-2000005,2000005,-1000006,2*-2000006,2000006,-1000011,
843      &2*-2000011,2000011,-1000012,2*-2000012,2000012,-1000013,
844      &2*-2000013,2000013,-1000014,2*-2000014,2000014,-1000015,
845      &2*-2000015,2000015,-1000016,2*-2000016,2000016,-1,-2,-3,-4,-5,-6,
846      &-7,-8,-11,-12,-13,-14,-15,-16,-17,-18,-24,-37,22,25,2*36,2,4,6,8,
847      &2,4,6,8,2,4,6,8,2,4,6,8,12,14,16,18,23,22,25,-1,-2,-3,-4,-5,-6/
848       DATA (KFDP(I,2),I= 338, 524)/-7,-8,-11,-13,-15,-17,21,22,2*23,
849      &-24,2*25,36,2*1000022,1000023,1000022,1000023,1000025,1000022,
850      &1000023,1000025,1000035,-1000024,-1000037,-1000024,-1000037,
851      &-1000001,2*-2000001,2000001,-1000002,2*-2000002,2000002,-1000003,
852      &2*-2000003,2000003,-1000004,2*-2000004,2000004,-1000005,
853      &2*-2000005,2000005,-1000006,2*-2000006,2000006,-1000011,
854      &2*-2000011,2000011,-1000012,2*-2000012,2000012,-1000013,
855      &2*-2000013,2000013,-1000014,2*-2000014,2000014,-1000015,
856      &2*-2000015,2000015,-1000016,2*-2000016,2000016,-1,-2,-3,-4,-5,-6,
857      &-7,-8,-11,-13,-15,-17,21,22,2*23,-24,25,2*1000022,1000023,
858      &1000022,1000023,1000025,1000022,1000023,1000025,1000035,-1000024,
859      &-1000037,-1000024,-1000037,-1000001,2*-2000001,2000001,-1000002,
860      &2*-2000002,2000002,-1000003,2*-2000003,2000003,-1000004,
861      &2*-2000004,2000004,-1000005,2*-2000005,2000005,-1000006,
862      &2*-2000006,2000006,-1000011,2*-2000011,2000011,-1000012,
863      &2*-2000012,2000012,-1000013,2*-2000013,2000013,-1000014,
864      &2*-2000014,2000014,-1000015,2*-2000015,2000015,-1000016,
865      &2*-2000016,2000016,2,4,6,8,12,14,16,18,25,1000024,1000037,
866      &1000024,1000037,1000024,1000037,1000024,1000037,2*-1000005,
867      &2*-2000005,1000002,1000004,1000012,1000014,2*1000016,-5,-6,21,11/
868       DATA (KFDP(I,2),I= 525, 940)/-3,-4,-5,-6,-7,-8,-13,-15,-17,-4,-5,
869      &-11,-13,-15,-5,-3,12,14,16,-24,-52,-24,-52,-1,-2,-3,-4,-5,-6,-7,
870      &-8,-11,-12,-13,-14,-15,-16,-17,-18,23,51,23,51,2,4,6,8,2,4,6,8,2,
871      &4,6,8,2,4,6,8,12,14,16,18,2*51,-1,-2,-3,-4,-5,-6,-7,-8,-11,-12,
872      &-13,-14,-15,-16,-17,-18,-82,12,14,-1,-3,11,13,15,1,4,3,4,1,3,22,
873      &11,-211,2*22,-13,-11,-211,211,111,211,-321,130,310,22,2*111,-211,
874      &11,-11,13,-13,-211,111,22,14,12,111,22,111,3*211,-311,22,211,22,
875      &111,-211,211,11,-211,13,22,-211,111,-211,22,111,-11,-211,111,
876      &2*-211,-321,130,310,221,111,-211,111,2*0,-211,111,22,-211,111,
877      &-211,111,-211,211,-213,113,223,221,14,111,211,111,-11,-13,211,
878      &111,22,211,111,211,111,2*211,213,113,223,221,22,-211,111,113,223,
879      &22,111,-321,310,211,111,2*-211,221,22,-11,-13,-211,-321,130,310,
880      &221,-211,111,11*12,11*14,2*211,2*213,211,20213,2*321,2*323,211,
881      &213,211,213,211,213,211,213,211,213,211,213,3*211,213,211,2*321,
882      &8*211,2*113,3*211,111,22,211,111,211,111,4*211,8*12,8*14,2*211,
883      &2*213,2*111,221,2*113,223,333,20213,211,2*321,323,2*311,313,-211,
884      &111,113,2*211,321,2*211,311,321,310,211,-211,4*211,321,4*211,113,
885      &2*211,-321,111,22,-211,111,-211,111,-211,211,-211,211,16,5*12,
886      &5*14,3*211,3*213,211,2*111,2*113,2*-311,2*-313,-2112,3*321,323,
887      &2*-1,22,111,321,311,321,311,-82,-11,-13,-82,22,-82,6*-11,6*-13/
888       DATA (KFDP(I,2),I= 941,1318)/2*-15,211,213,20213,211,213,20213,
889      &431,433,431,433,311,313,311,313,311,313,-1,-4,-3,-4,-1,-3,22,
890      &-211,111,-211,111,-211,211,-211,211,6*-11,6*-13,2*-15,211,213,
891      &20213,211,213,20213,431,433,431,433,321,323,321,323,321,323,-1,
892      &-4,-3,-4,-1,-3,22,211,111,211,111,4*211,6*-11,6*-13,2*-15,211,
893      &213,20213,211,213,20213,431,433,431,433,221,331,333,221,331,333,
894      &221,331,333,-1,-4,-3,-4,-1,-3,22,-321,-311,-321,-311,-15,-3,-1,
895      &2*-11,2*-13,2*-15,-1,-4,-3,-4,-3,-4,-1,-4,2*12,2*14,2,3,2,3,2*12,
896      &2*14,2,1,22,411,421,411,421,21,-11,-13,-15,-1,-2,-3,-4,2*21,22,
897      &21,2*-211,111,22,111,211,22,211,-211,11,2*-211,111,-211,111,22,
898      &11,22,111,-211,211,111,211,22,211,111,211,-211,22,11,13,11,-211,
899      &2*111,2*22,111,211,-321,-211,111,11,2*-211,7*12,7*14,-321,-323,
900      &-311,-313,-311,-313,211,213,211,213,211,213,111,221,331,113,223,
901      &111,221,113,223,321,323,321,-211,-213,111,221,331,113,223,333,
902      &10221,111,221,331,113,223,211,213,211,213,321,323,321,323,321,
903      &323,311,313,311,313,2*-1,-3,-1,2203,3201,3203,2203,2101,2103,12,
904      &14,-1,-3,2*111,2*211,12,14,-1,-3,22,111,2*22,111,22,12,14,-1,-3,
905      &22,12,14,-1,-3,12,14,-1,-3,12,14,-1,-3,12,14,-1,-3,12,14,-1,-3,
906      &12,14,-1,-3,12,14,-1,-3,2*-211,11,13,15,-211,-213,-20213,-431,
907      &-433,3*3122,1,4,3,4,1,3,11,13,15,1,4,3,4,1,3,11,13,15,1,4,3,4,1/
908       DATA (KFDP(I,2),I=1319,1774)/3,2*111,2*211,11,13,15,1,4,3,4,1,3,
909      &11,13,15,1,4,3,4,1,3,4*22,11,13,15,1,4,3,4,1,3,22,11,13,15,1,4,3,
910      &4,1,3,11,13,15,1,4,3,4,1,3,11,13,15,1,4,3,4,1,3,11,13,15,1,4,3,4,
911      &1,3,11,13,15,1,4,3,4,1,3,11,13,15,1,4,3,4,1,3,11,13,15,1,4,3,4,1,
912      &3,11,13,15,1,4,3,4,1,3,11,13,15,1,4,3,4,1,3,11,13,15,1,4,3,4,1,3,
913      &11,13,15,1,4,3,4,1,3,11,13,15,1,4,3,4,1,3,11,13,15,1,4,3,4,1,3,
914      &11,13,15,1,4,3,4,1,3,11,13,15,1,4,3,4,1,3,11,13,15,1,4,3,4,1,3,
915      &11,13,15,1,4,3,4,1,3,11,13,15,1,4,3,4,1,3,2*111,2*211,-211,111,
916      &-321,130,310,-211,111,211,-211,111,-213,113,-211,111,223,211,111,
917      &213,113,211,111,223,-211,111,-321,130,310,2*-211,-311,311,-321,
918      &321,211,111,211,111,-211,111,-211,111,311,2*321,311,22,2*-82,
919      &-211,111,-211,111,211,111,211,111,-321,-311,-321,-311,411,421,
920      &411,421,22,2*21,-211,2*211,111,-211,111,2*211,111,-211,211,111,
921      &211,-321,2*-311,-321,22,-211,111,211,111,-311,311,-321,321,211,
922      &111,-211,111,321,311,22,-82,-211,111,211,111,-321,-311,411,421,
923      &22,21,-11,-13,-82,211,111,221,111,4*22,-11,-13,-15,-1,-2,-3,-4,
924      &2*21,211,111,3*22,1,2*2,4*1,2*-24,2*-37,1,2,2*1,4*2,2*24,2*37,2,
925      &3,2*4,4*3,2*-24,2*-37,3,4,2*3,4*4,2*24,2*37,4,5,2*6,4*5,2*-24,
926      &2*-37,5,6,2*5,4*6,2*24,2*37,6,4,11,2*12,4*11,2*-24,2*-37,12,2*11,
927      &4*12,2*24,2*37,13,2*14,4*13,2*-24,2*-37,14,2*13,4*14,2*24,2*37/
928       DATA (KFDP(I,2),I=1775,2218)/15,2*16,4*15,2*-24,2*-37,16,2*15,
929      &4*16,2*24,2*37,21,-1,1,-1,1,-2,2,-2,2,-3,3,-3,3,-4,4,-4,4,-5,5,
930      &-5,5,-6,6,-6,6,1,3,5,2,4,6,1,3,5,2,4,6,1,3,5,2,4,6,1,3,5,2,4,6,1,
931      &-1,3,-3,5,-5,1,-1,3,-3,5,-5,22,23,25,35,36,22,23,25,35,36,22,23,
932      &11,13,15,12,14,16,1,3,5,2,4,25,35,36,-24,24,11,-11,13,-13,15,-15,
933      &1,-1,3,-3,-24,24,11,-11,13,-13,15,-15,1,-1,3,-3,-37,37,-37,37,-1,
934      &1,-1,1,-2,2,-2,2,-3,3,-3,3,-4,4,-4,4,-5,5,-5,5,-6,6,-6,6,-11,11,
935      &-11,11,-12,12,-12,12,-13,13,-13,13,-14,14,-14,14,-15,15,-15,15,
936      &-16,16,-16,16,1,3,5,2,4,24,37,24,-11,-13,-15,-1,-3,24,-11,-13,
937      &-15,-1,-3,24,-11,-13,-15,-1,-3,24,-11,-13,-15,-1,-3,4*37,2*-1,
938      &2*2,2*-3,2*4,2*-5,2*6,2*-11,2*12,2*-13,2*14,2*-15,2*16,-1,-3,22,
939      &23,25,35,36,22,23,11,13,15,12,14,16,1,3,5,2,4,25,35,36,22,23,11,
940      &13,15,12,14,16,1,3,5,2,4,25,35,36,-24,24,11,-11,13,-13,15,-15,1,
941      &-1,3,-3,-24,24,11,-11,13,-13,15,-15,1,-1,3,-3,-37,37,-37,37,-1,1,
942      &-1,1,-2,2,-2,2,-3,3,-3,3,-4,4,-4,4,-5,5,-5,5,-6,6,-6,6,-11,11,
943      &-11,11,-12,12,-12,12,-13,13,-13,13,-14,14,-14,14,-15,15,-15,15,
944      &-16,16,-16,16,1,3,5,2,4,22,23,25,35,36,22,23,11,13,15,12,14,16,1,
945      &3,5,2,4,25,35,36,22,23,11,13,15,12,14,16,1,3,5,2,4,25,35,36,22,
946      &23,11,13,15,12,14,16,1,3,5,2,4,25,35,36,-24,24,11,-11,13,-13,15,
947      &-15,1,-1,3,-3,-24,24,11,-11,13,-13,15,-15,1,-1,3,-3,-37,37,-37/
948       DATA (KFDP(I,2),I=2219,4000)/37,-1,1,-1,1,-2,2,-2,2,-3,3,-3,3,-4,
949      &4,-4,4,-5,5,-5,5,-6,6,-6,6,-11,11,-11,11,-12,12,-12,12,-13,13,
950      &-13,13,-14,14,-14,14,-15,15,-15,15,-16,16,-16,16,1,3,5,2,4,24,37,
951      &23,11,13,15,12,14,16,1,3,5,2,4,25,35,36,24,-11,-13,-15,-1,-3,24,
952      &-11,-13,-15,-1,-3,24,-11,-13,-15,-1,-3,24,-11,-13,-15,-1,-3,4*37,
953      &2*-1,2*2,2*-3,2*4,2*-5,2*6,2*-11,2*12,2*-13,2*14,2*-15,2*16,-1,
954      &-3,1,2*2,4*1,23,25,35,36,2*-24,2*-37,1,2,2*1,4*2,23,25,35,36,
955      &2*24,2*37,2,3,2*4,4*3,23,25,35,36,2*-24,2*-37,3,4,2*3,4*4,23,25,
956      &35,36,2*24,2*37,4,5,2*6,4*5,23,25,35,36,2*-24,2*-37,5,6,2*5,4*6,
957      &23,25,35,36,2*24,2*37,6,11,2*12,4*11,23,25,35,36,2*-24,2*-37,13,
958      &2*14,4*13,23,25,35,36,2*-24,2*-37,15,2*16,4*15,23,25,35,36,2*-24,
959      &2*-37,3*1,4*2,1,2*11,2*12,11,1503*0/
960       DATA (KFDP(I,3),I=   1,1087)/79*0,14,6*0,2*16,2*0,6*111,310,130,
961      &2*0,3*111,310,130,321,113,211,223,221,2*113,2*211,2*223,2*221,
962      &2*113,221,2*113,2*213,-213,113,2*111,310,130,310,130,2*310,130,
963      &470*0,4*3,4*4,1,4,3,2*2,0,-11,8*0,-211,5*0,2*111,211,-211,211,
964      &-211,10*0,111,4*0,2*111,-211,-11,11,-13,22,111,3*0,22,3*0,111,
965      &211,4*0,111,11*0,111,-211,6*0,-211,3*111,7*0,111,-211,5*0,2*221,
966      &3*0,111,5*0,111,11*0,-311,-313,-311,-321,-313,-323,111,221,331,
967      &113,223,-311,-313,-311,-321,-313,-323,111,221,331,113,223,22*0,
968      &111,113,2*211,-211,-311,211,111,3*211,-211,7*211,7*0,111,-211,
969      &111,-211,-321,-323,-311,-321,-313,-323,-211,-213,-321,-323,-311,
970      &-321,-313,-323,-211,-213,22*0,111,113,-311,2*-211,211,-211,310,
971      &-211,2*111,211,2*-211,-321,-211,2*211,-211,111,-211,2*211,6*0,
972      &111,-211,111,-211,0,221,331,333,321,311,221,331,333,321,311,20*0,
973      &3,13*0,-411,-413,-10413,-10411,-20413,-415,-411,-413,-10413,
974      &-10411,-20413,-415,-411,-413,16*0,-4,-1,-4,-3,2*-2,5*0,111,-211,
975      &111,-211,-421,-423,-10423,-10421,-20423,-425,-421,-423,-10423,
976      &-10421,-20423,-425,-421,-423,16*0,-4,-1,-4,-3,2*-2,5*0,111,-211,
977      &111,-211,-431,-433,-10433,-10431,-20433,-435,-431,-433,-10433,
978      &-10431,-20433,-435,-431,-433,19*0,-4,-1,-4,-3,2*-2,8*0,441,443,
979      &441,443,441,443,-4,-1,-4,-3,-4,-3,-4,-1,531,533,531,533,3,2,3,2/
980       DATA (KFDP(I,3),I=1088,2186)/511,513,511,513,1,2,13*0,2*21,11*0,
981      &2112,6*0,2212,12*0,2*3122,3212,10*0,3322,2*0,3122,3212,3214,2112,
982      &2114,2212,2112,3122,3212,3214,2112,2114,2212,2112,52*0,3*3,1,6*0,
983      &4*3,4*0,4*3,6*0,4*3,0,28*3,2*0,3*4122,8*0,4,1,4,3,2*2,4*4,1,4,3,
984      &2*2,4*4,1,4,3,2*2,4*0,4*4,1,4,3,2*2,4*4,1,4,3,2*2,4*0,4*4,1,4,3,
985      &2*2,0,4*4,1,4,3,2*2,4*4,1,4,3,2*2,4*4,1,4,3,2*2,4*4,1,4,3,2*2,
986      &4*4,1,4,3,2*2,4*4,1,4,3,2*2,4*4,1,4,3,2*2,4*4,1,4,3,2*2,4*4,1,4,
987      &3,2*2,4*4,1,4,3,2*2,4*4,1,4,3,2*2,4*4,1,4,3,2*2,4*4,1,4,3,2*2,
988      &4*4,1,4,3,2*2,4*4,1,4,3,2*2,4*4,1,4,3,2*2,4*4,1,4,3,2*2,4*4,1,4,
989      &3,2*2,31*0,211,111,45*0,-211,2*111,-211,3*111,-211,111,211,30*0,
990      &-211,111,13*0,2*21,-211,111,167*0,-1,-3,-5,-2,-4,-6,-1,-3,-5,-2,
991      &-4,-6,-1,-3,-5,-2,-4,-6,-1,-3,-5,-2,-4,-6,-2,2,-4,4,-6,6,-2,2,-4,
992      &4,-6,6,12*0,-11,-13,-15,-12,-14,-16,-1,-3,-5,-2,-4,5*0,-12,12,
993      &-14,14,-16,16,-2,2,-4,4,2*0,-12,12,-14,14,-16,16,-2,2,-4,4,52*0,
994      &-1,-3,-5,-2,-4,3*0,12,14,16,2,4,0,12,14,16,2,4,0,12,14,16,2,4,0,
995      &12,14,16,2,4,28*0,2,4,7*0,-11,-13,-15,-12,-14,-16,-1,-3,-5,-2,-4,
996      &5*0,-11,-13,-15,-12,-14,-16,-1,-3,-5,-2,-4,5*0,-12,12,-14,14,-16,
997      &16,-2,2,-4,4,2*0,-12,12,-14,14,-16,16,-2,2,-4,4,52*0,-1,-3,-5,-2,
998      &-4,7*0,-11,-13,-15,-12,-14,-16,-1,-3,-5,-2,-4,5*0,-11,-13,-15,
999      &-12,-14,-16,-1,-3,-5,-2,-4,5*0,-11,-13,-15,-12,-14,-16,-1,-3,-5/
1000       DATA (KFDP(I,3),I=2187,4000)/-2,-4,5*0,-12,12,-14,14,-16,16,-2,2,
1001      &-4,4,2*0,-12,12,-14,14,-16,16,-2,2,-4,4,52*0,-1,-3,-5,-2,-4,3*0,
1002      &-11,-13,-15,-12,-14,-16,-1,-3,-5,-2,-4,4*0,12,14,16,2,4,0,12,14,
1003      &16,2,4,0,12,14,16,2,4,0,12,14,16,2,4,28*0,2,4,1657*0/
1004       DATA (KFDP(I,4),I=   1,4000)/92*0,4*111,6*0,111,2*0,-211,0,-211,
1005      &3*0,111,2*-211,0,111,0,2*111,113,221,2*111,-213,-211,211,113,
1006      &6*111,310,2*130,470*0,13*81,41*0,-11,10*0,111,-211,4*0,111,62*0,
1007      &111,211,111,211,7*0,111,211,111,211,35*0,2*-211,2*111,211,111,
1008      &-211,2*211,2*-211,13*0,-211,111,-211,111,4*0,-211,111,-211,111,
1009      &34*0,111,-211,3*111,3*-211,2*111,3*-211,14*0,-321,-311,3*0,-321,
1010      &-311,20*0,-3,43*0,6*1,39*0,6*2,42*0,6*3,14*0,8*4,4*0,4*-5,4*0,
1011      &2*-5,67*0,-211,111,5*0,-211,111,52*0,2101,2103,2*2101,6*0,4*81,
1012      &4*0,4*81,6*0,4*81,0,28*81,13*0,6*2101,18*81,4*0,18*81,4*0,9*81,0,
1013      &162*81,31*0,-211,111,2450*0/
1014       DATA (KFDP(I,5),I=   1,4000)/94*0,2*111,17*0,111,7*0,2*111,0,
1015      &3*111,0,111,665*0,-211,2*111,-211,111,-211,111,65*0,111,-211,
1016      &3*111,-211,111,3127*0/
1017
1018 C...PYDAT4, with particle names (character strings).
1019       DATA (CHAF(I,1),I=   1, 190)/'d','u','s','c','b','t','b''','t''',
1020      &2*' ','e-','nu_e','mu-','nu_mu','tau-','nu_tau','tau''-',
1021      &'nu''_tau',2*' ','g','gamma','Z0','W+','h0',2*' ','reggeon',
1022      &'pomeron',2*' ','Z''0','Z"0','W''+','H0','A0','H+','eta_tech0',
1023      &'LQ_ue','R0',10*' ','pi_tech0','pi_tech+','pi''_tech0',
1024      &'rho_tech0','rho_tech+','omega_tech',24*' ','specflav',
1025      &'rndmflav','phasespa','c-hadron','b-hadron',5*' ','cluster',
1026      &'string','indep.','CMshower','SPHEaxis','THRUaxis','CLUSjet',
1027      &'CELLjet','table',' ','rho_diff0','pi0','rho0','a_20','K_L0',
1028      &'pi_diffr+','pi+','rho+','a_2+','omega_di','eta','omega','f_2',
1029      &'K_S0','K0','K*0','K*_20','K+','K*+','K*_2+','phi_diff','eta''',
1030      &'phi','f''_2','D+','D*+','D*_2+','D0','D*0','D*_20','D_s+',
1031      &'D*_s+','D*_2s+','J/psi_di','eta_c','J/psi','chi_2c','B0','B*0',
1032      &'B*_20','B+','B*+','B*_2+','B_s0','B*_s0','B*_2s0','B_c+',
1033      &'B*_c+','B*_2c+','eta_b','Upsilon','chi_2b','dd_1','Delta-',
1034      &'ud_0','ud_1','n_diffr0','n0','Delta0','uu_1','p_diffr+','p+',
1035      &'Delta+','Delta++','sd_0','sd_1','Sigma-','Sigma*-','Lambda0',
1036      &'su_0','su_1','Sigma0','Sigma*0','Sigma+','Sigma*+','ss_1','Xi-',
1037      &'Xi*-','Xi0','Xi*0','Omega-','cd_0','cd_1','Sigma_c0',
1038      &'Sigma*_c0','Lambda_c+','Xi_c0','cu_0','cu_1','Sigma_c+'/
1039       DATA (CHAF(I,1),I= 191, 317)/'Sigma*_c+','Sigma_c++',
1040      &'Sigma*_c++','Xi_c+','cs_0','cs_1','Xi''_c0','Xi*_c0','Xi''_c+',
1041      &'Xi*_c+','Omega_c0','Omega*_c0','cc_1','Xi_cc+','Xi*_cc+',
1042      &'Xi_cc++','Xi*_cc++','Omega_cc+','Omega*_cc+','Omega*_ccc++',
1043      &'bd_0','bd_1','Sigma_b-','Sigma*_b-','Lambda_b0','Xi_b-',
1044      &'Xi_bc0','bu_0','bu_1','Sigma_b0','Sigma*_b0','Sigma_b+',
1045      &'Sigma*_b+','Xi_b0','Xi_bc+','bs_0','bs_1','Xi''_b-','Xi*_b-',
1046      &'Xi''_b0','Xi*_b0','Omega_b-','Omega*_b-','Omega_bc0','bc_0',
1047      &'bc_1','Xi''_bc0','Xi*_bc0','Xi''_bc+','Xi*_bc+','Omega''_bc0',
1048      &'Omega*_bc0','Omega_bcc+','Omega*_bcc+','bb_1','Xi_bb-',
1049      &'Xi*_bb-','Xi_bb0','Xi*_bb0','Omega_bb-','Omega*_bb-',
1050      &'Omega_bbc0','Omega*_bbc0','Omega*_bbb-','a_00','b_10','a_0+',
1051      &'b_1+','f_0','h_1','K*_00','K_10','K*_0+','K_1+','f''_0','h''_1',
1052      &'D*_0+','D_1+','D*_00','D_10','D*_0s+','D_1s+','chi_0c','h_1c',
1053      &'B*_00','B_10','B*_0+','B_1+','B*_0s0','B_1s0','B*_0c+','B_1c+',
1054      &'chi_0b','h_1b','a_10','a_1+','f_1','K*_10','K*_1+','f''_1',
1055      &'D*_1+','D*_10','D*_1s+','chi_1c','B*_10','B*_1+','B*_1s0',
1056      &'B*_1c+','chi_1b','psi''','Upsilon''','~d_L','~u_L','~s_L',
1057      &'~c_L','~b_1','~t_1','~e_L-','~nu_eL','~mu_L-','~nu_muL',
1058      &'~tau_1-','~nu_tauL','~g','~chi_10','~chi_20','~chi_1+'/
1059       DATA (CHAF(I,1),I= 318, 500)/'~chi_30','~chi_40','~chi_2+',
1060      &'~gravitino','~d_R','~u_R','~s_R','~c_R','~b_2','~t_2','~e_R-',
1061      &'~nu_eR','~mu_R-','~nu_muR','~tau_2-','~nu_tauR','d*','u*','e*-',
1062      &'nu*_e0',163*' '/
1063       DATA (CHAF(I,2),I=   1, 206)/'dbar','ubar','sbar','cbar','bbar',
1064      &'tbar','b''bar','t''bar',2*' ','e+','nu_ebar','mu+','nu_mubar',
1065      &'tau+','nu_taubar','tau''+','nu''_taubar',5*' ','W-',9*' ',
1066      &'W''-',2*' ','H-',' ','LQ_uebar','Rbar0',11*' ','pi_tech-',2*' ',
1067      &'rho_tech-',26*' ','rndmflavbar',' ','c-hadronbar','b-hadronbar',
1068      &20*' ','pi_diffr-','pi-','rho-','a_2-',5*' ','Kbar0','K*bar0',
1069      &'K*_2bar0','K-','K*-','K*_2-',4*' ','D-','D*-','D*_2-','Dbar0',
1070      &'D*bar0','D*_2bar0','D_s-','D*_s-','D*_2s-',4*' ','Bbar0',
1071      &'B*bar0','B*_2bar0','B-','B*-','B*_2-','B_sbar0','B*_sbar0',
1072      &'B*_2sbar0','B_c-','B*_c-','B*_2c-',3*' ','dd_1bar','Deltabar+',
1073      &'ud_0bar','ud_1bar','n_diffrbar0','nbar0','Deltabar0','uu_1bar',
1074      &'p_diffrbar-','pbar-','Deltabar-','Deltabar--','sd_0bar',
1075      &'sd_1bar','Sigmabar+','Sigma*bar+','Lambdabar0','su_0bar',
1076      &'su_1bar','Sigmabar0','Sigma*bar0','Sigmabar-','Sigma*bar-',
1077      &'ss_1bar','Xibar+','Xi*bar+','Xibar0','Xi*bar0','Omegabar+',
1078      &'cd_0bar','cd_1bar','Sigma_cbar0','Sigma*_cbar0','Lambda_cbar-',
1079      &'Xi_cbar0','cu_0bar','cu_1bar','Sigma_cbar-','Sigma*_cbar-',
1080      &'Sigma_cbar--','Sigma*_cbar--','Xi_cbar-','cs_0bar','cs_1bar',
1081      &'Xi''_cbar0','Xi*_cbar0','Xi''_cbar-','Xi*_cbar-','Omega_cbar0',
1082      &'Omega*_cbar0','cc_1bar','Xi_ccbar-','Xi*_ccbar-','Xi_ccbar--'/
1083       DATA (CHAF(I,2),I= 207, 324)/'Xi*_ccbar--','Omega_ccbar-',
1084      &'Omega*_ccbar-','Omega*_cccbar-','bd_0bar','bd_1bar',
1085      &'Sigma_bbar+','Sigma*_bbar+','Lambda_bbar0','Xi_bbar+',
1086      &'Xi_bcbar0','bu_0bar','bu_1bar','Sigma_bbar0','Sigma*_bbar0',
1087      &'Sigma_bbar-','Sigma*_bbar-','Xi_bbar0','Xi_bcbar-','bs_0bar',
1088      &'bs_1bar','Xi''_bbar+','Xi*_bbar+','Xi''_bbar0','Xi*_bbar0',
1089      &'Omega_bbar+','Omega*_bbar+','Omega_bcbar0','bc_0bar','bc_1bar',
1090      &'Xi''_bcbar0','Xi*_bcbar0','Xi''_bcbar-','Xi*_bcbar-',
1091      &'Omega''_bcba','Omega*_bcbar0','Omega_bccbar-','Omega*_bccbar-',
1092      &'bb_1bar','Xi_bbbar+','Xi*_bbbar+','Xi_bbbar0','Xi*_bbbar0',
1093      &'Omega_bbbar+','Omega*_bbbar+','Omega_bbcbar0','Omega*_bbcbar0',
1094      &'Omega*_bbbbar+',2*' ','a_0-','b_1-',2*' ','K*_0bar0','K_1bar0',
1095      &'K*_0-','K_1-',2*' ','D*_0-','D_1-','D*_0bar0','D_1bar0',
1096      &'D*_0s-','D_1s-',2*' ','B*_0bar0','B_1bar0','B*_0-','B_1-',
1097      &'B*_0sbar0','B_1sbar0','B*_0c-','B_1c-',3*' ','a_1-',' ',
1098      &'K*_1bar0','K*_1-',' ','D*_1-','D*_1bar0','D*_1s-',' ',
1099      &'B*_1bar0','B*_1-','B*_1sbar0','B*_1c-',3*' ','~d_Lbar',
1100      &'~u_Lbar','~s_Lbar','~c_Lbar','~b_1bar','~t_1bar','~e_L+',
1101      &'~nu_eLbar','~mu_L+','~nu_muLbar','~tau_1+','~nu_tauLbar',3*' ',
1102      &'~chi_1-',2*' ','~chi_2-',' ','~d_Rbar','~u_Rbar','~s_Rbar'/
1103       DATA (CHAF(I,2),I= 325, 500)/'~c_Rbar','~b_2bar','~t_2bar',
1104      &'~e_R+','~nu_eRbar','~mu_R+','~nu_muRbar','~tau_2+',
1105      &'~nu_tauRbar','d*bar','u*bar','e*bar+','nu*_ebar0',163*' '/
1106
1107 C...PYDATR, with initial values for the random number generator.
1108       DATA MRPY/19780503,0,0,97,33,0/
1109
1110 C...Default values for allowed processes and kinematics constraints.
1111       DATA MSEL/1/
1112       DATA MSUB/500*0/
1113       DATA ((KFIN(I,J),J=-40,40),I=1,2)/16*0,4*1,4*0,6*1,5*0,5*1,0,
1114      &5*1,5*0,6*1,4*0,4*1,16*0,16*0,4*1,4*0,6*1,5*0,5*1,0,5*1,5*0,
1115      &6*1,4*0,4*1,16*0/
1116       DATA CKIN/
1117      &  2.0D0, -1.0D0,  0.0D0, -1.0D0,  1.0D0,
1118      &  1.0D0,  -10D0,   10D0,  -10D0,   10D0,
1119      1  -10D0,   10D0,  -10D0,   10D0,  -10D0,
1120      1   10D0, -1.0D0,  1.0D0, -1.0D0,  1.0D0,
1121      2  0.0D0,  1.0D0,  0.0D0,  1.0D0, -1.0D0,
1122      2  1.0D0, -1.0D0,  1.0D0,    0D0,    0D0,
1123      3  2.0D0, -1.0D0,    0D0,    0D0,  0.0D0,
1124      3 -1.0D0,  0.0D0, -1.0D0,  4.0D0, -1.0D0,
1125      4 12.0D0, -1.0D0, 12.0D0, -1.0D0, 12.0D0,
1126      4 -1.0D0, 12.0D0, -1.0D0,    0D0,    0D0,
1127      5  0.0D0, -1.0D0,  0.0D0, -1.0D0,  0.0D0,
1128      5 -1.0D0,    0D0,    0D0,    0D0,    0D0,
1129      6  140*0D0/
1130
1131 C...Default values for main switches and parameters. Reset information.
1132       DATA (MSTP(I),I=1,100)/
1133      &  3,    1,    2,    0,    0,    0,    0,    0,    0,    0,
1134      1  1,    0,    1,    0,    5,    0,    0,    0,    0,    0,
1135      2  1,    0,    1,    0,    0,    0,    0,    0,    0,    1,
1136      3  1,    2,    0,    1,    0,    2,    1,    5,    2,    0,
1137      4  1,    1,    3,    7,    3,    1,    1,    0,    1,    0,
1138      5  4,    1,    3,    1,    5,    1,    1,    6,    1,    7,
1139      6  1,    3,    2,    2,    1,    1,    2,    0,    0,    0,
1140      7  1,    0,    0,    0,    0,    0,    0,    0,    0,    0,
1141      8  1,    1,  100,    0,    0,    0,    0,    0,    0,    0,
1142      9  1,    4,    1,    2,    0,    0,    0,    0,    0,    0/
1143       DATA (MSTP(I),I=101,200)/
1144      &  3,    1,    0,    0,    0,    0,    0,    0,    0,    0,
1145      1  1,    1,    1,    0,    0,    0,    0,    0,    0,    0,
1146      2  0,    1,    2,    1,    1,   50,    0,    0,   10,    0,
1147      3  0,    4,    0,    1,    0,    0,    0,    0,    0,    0,
1148      4  0,    0,    0,    0,    0,    0,    0,    0,    0,    0,
1149      5  0,    0,    0,    0,    0,    0,    0,    0,    0,    0,
1150      6  0,    0,    0,    0,    0,    0,    0,    0,    0,    0,
1151      7  0,    2,    0,    0,    0,    0,    0,    0,    0,    0,
1152      8  6,  115, 1998,   01,   27,    0,    0,    0,    0,    0,
1153      9  0,    0,    0,    0,    0,    0,    0,    0,    0,    0/
1154       DATA (PARP(I),I=1,100)/
1155      &  0.25D0,  10D0, 8*0D0,
1156      1  0D0,   0D0,  1.0D0, 0.01D0,  0.6D0,  1.0D0,  1.0D0, 3*0D0,
1157      2  10*0D0,
1158      3  1.5D0,2.0D0,0.075D0,1.0D0,0.2D0,0D0,2.0D0,0.70D0,0.006D0,0D0,
1159      4  0.02D0,2.0D0,0.10D0,1000D0,2054D0, 123D0, 246D0, 50D0, 2*0D0,
1160      5  1.0D0, 9*0D0,
1161      6  0.25D0, 1.0D0,0.25D0, 1.0D0, 2.0D0,1D-3, 4.0D0,1D-3,2*0D0,
1162      7  4.0D0, 0.25D0, 8*0D0,
1163      8  1.40D0,1.55D0,0.5D0, 0.2D0,0.33D0,0.66D0, 0.7D0, 0.5D0,2*0D0,
1164      9  0.44D0,0.20D0,2.0D0,1.0D0,0D0,3.0D0,1.0D0,0.75D0,0.44D0,2.0D0/
1165       DATA (PARP(I),I=101,200)/
1166      &  0.5D0, 0.28D0,  1.0D0, 0.8D0, 6*0D0,
1167      1  2.0D0, 3*0D0, 1.5D0, 0.5D0, 0.6D0, 2.5D0, 2.0D0, 1.0D0,
1168      2  1.0D0,  0.4D0, 8*0D0,
1169      3  0.01D0, 9*0D0,
1170      4  0.33333D0, 82D0, 1D0, 4D0, 200D0, 5*0D0,
1171      5  0D0,   0D0,   0D0,   0D0, 6*0D0,
1172      6  2.20D0, 23.6D0, 18.4D0, 11.5D0, 6*0D0,
1173      7  0D0,   0D0,   0D0,  1.0D0, 6*0D0,
1174      8  20*0D0/
1175       DATA MSTI/200*0/
1176       DATA PARI/200*0D0/
1177       DATA MINT/400*0/
1178       DATA VINT/400*0D0/
1179
1180 C...Constants for the generation of the various processes.
1181       DATA (ISET(I),I=1,100)/
1182      &  1,    1,    1,   -1,    3,   -1,   -1,    3,   -2,    2,
1183      1  2,    2,    2,    2,    2,    2,   -1,    2,    2,    2,
1184      2 -1,    2,    2,    2,    2,    2,   -1,    2,    2,    2,
1185      3  2,   -1,    2,    2,    2,    2,   -1,   -1,   -1,   -1,
1186      4 -1,   -1,   -1,   -1,   -1,   -1,   -1,   -1,   -1,   -1,
1187      5 -1,   -1,    2,    2,   -1,   -1,   -1,    2,   -1,   -1,
1188      6 -1,   -1,   -1,   -1,   -1,   -1,   -1,    2,    2,    2,
1189      7  4,    4,    4,   -1,   -1,    4,    4,   -1,   -1,    2,
1190      8  2,    2,    2,    2,    2,    2,    2,    2,    2,   -2,
1191      9  0,    0,    0,    0,    0,    9,   -2,   -2,   -2,   -2/
1192       DATA (ISET(I),I=101,200)/
1193      & -1,    1,    1,   -2,   -2,    2,    2,    2,   -2,    2,
1194      1  2,    2,    2,    2,    2,   -1,   -1,   -1,   -2,   -2,
1195      2  5,    5,    5,    5,   -2,   -2,   -2,   -2,   -2,   -2,
1196      3 -1,   -2,   -2,   -2,   -2,   -2,   -2,   -2,   -2,   -2,
1197      4  1,    1,    1,    1,    1,   -2,    1,    1,    1,   -2,
1198      5  1,    1,    1,   -2,   -2,    1,    1,    1,   -2,   -2,
1199      6  2,    2,    2,    2,    2,    2,    2,    2,   -2,   -2,
1200      7  2,    2,    5,    5,   -2,    2,    2,    5,    5,   -2,
1201      8  5,    5,   -2,   -2,   -2,    5,    5,   -2,   -2,   -2,
1202      9  1,    1,    1,    2,   -2,   -2,   -2,   -2,   -2,   -2/
1203       DATA (ISET(I),I=201,300)/
1204      &  2,    2,    2,    2,    2,    2,    2,    2,    2,    2,
1205      1  2,    2,    2,    2,   -2,    2,    2,    2,    2,    2,
1206      2  2,    2,    2,    2,    2,    2,    2,    2,    2,    2,
1207      3  2,    2,    2,    2,    2,    2,    2,    2,    2,    2,
1208      4  2,    2,    2,    2,   -1,    2,    2,    2,    2,    2,
1209      5  2,    2,    2,    2,   -1,    2,   -1,    2,    2,   -2,
1210      6  2,    2,    2,    2,    2,   -1,   -1,   -1,   -1,   -1,
1211      7  2,    2,    2,    2,    2,    2,    2,    2,    2,    2,
1212      8 -2,   -2,   -2,   -2,   -2,   -2,   -2,   -2,   -2,   -2,
1213      9 -2,   -2,   -2,   -2,   -2,   -2,   -2,   -2,   -2,   -2/
1214       DATA (ISET(I),I=301,500)/200*-2/
1215       DATA ((KFPR(I,J),J=1,2),I=1,50)/
1216      &  23,    0,   24,    0,   25,    0,   24,    0,   25,    0,
1217      &  24,    0,   23,    0,   25,    0,    0,    0,    0,    0,
1218      1   0,    0,    0,    0,   21,   21,   21,   22,   21,   23,
1219      1  21,   24,   21,   25,   22,   22,   22,   23,   22,   24,
1220      2  22,   25,   23,   23,   23,   24,   23,   25,   24,   24,
1221      2  24,   25,   25,   25,    0,   21,    0,   22,    0,   23,
1222      3   0,   24,    0,   25,    0,   21,    0,   22,    0,   23,
1223      3   0,   24,    0,   25,    0,   21,    0,   22,    0,   23,
1224      4   0,   24,    0,   25,    0,   21,    0,   22,    0,   23,
1225      4   0,   24,    0,   25,    0,   21,    0,   22,    0,   23/
1226       DATA ((KFPR(I,J),J=1,2),I=51,100)/
1227      5   0,   24,    0,   25,    0,    0,    0,    0,    0,    0,
1228      5   0,    0,    0,    0,    0,    0,    0,    0,    0,    0,
1229      6   0,    0,    0,    0,    0,    0,    0,    0,    0,    0,
1230      6   0,    0,    0,    0,   21,   21,   24,   24,   23,   24,
1231      7  23,   23,   24,   24,   23,   24,   23,   25,   22,   22,
1232      7  23,   23,   24,   24,   24,   25,   25,   25,    0,  211,
1233      8   0,    0,    0,    0,    0,    0,    0,    0,    0,    0,
1234      8 443,   21,10441,   21,20443,   21,  445,   21,    0,    0,
1235      9   0,    0,    0,    0,    0,    0,    0,    0,    0,    0,
1236      9   0,    0,    0,    0,    0,    0,    0,    0,    0,    0/
1237       DATA ((KFPR(I,J),J=1,2),I=101,150)/
1238      &  23,    0,   25,    0,   25,    0,    0,    0,    0,    0,
1239      & 443,   22,  443,   21,  443,   22,    0,    0,   22,   25,
1240      1  21,   25,    0,   25,   21,   25,   22,   22,   21,   22,
1241      1  22,   23,   23,   23,   24,   24,    0,    0,    0,    0,
1242      2  25,    6,   25,    6,   25,    0,   25,    0,    0,    0,
1243      2   0,    0,    0,    0,    0,    0,    0,    0,    0,    0,
1244      3  23,    5,    0,    0,    0,    0,    0,    0,    0,    0,
1245      3   0,    0,    0,    0,    0,    0,    0,    0,    0,    0,
1246      4  32,    0,   34,    0,   37,    0,   40,    0,   39,    0,
1247      4   0,    0, 4000001, 0, 4000002, 0,   38,    0,    0,    0/
1248       DATA ((KFPR(I,J),J=1,2),I=151,200)/
1249      5  35,    0,   35,    0,   35,    0,    0,    0,    0,    0,
1250      5  36,    0,   36,    0,   36,    0,    0,    0,    0,    0,
1251      6   6,   37,   39,    0,   39,   39,   39,   39,   11,    0,
1252      6  11,    0, 0, 4000001, 0, 4000002,    0,    0,    0,    0,
1253      7  23,   35,   24,   35,   35,    0,   35,    0,    0,    0,
1254      7  23,   36,   24,   36,   36,    0,   36,    0,    0,    0,
1255      8  35,    6,   35,    6,    0,    0,    0,    0,    0,    0,
1256      8  36,    6,   36,    6,    0,    0,    0,    0,    0,    0,
1257      9  54,    0,   55,    0,   56,    0,   11,    0,    0,    0,
1258      9   0,    0,    0,    0,    0,    0,    0,    0,    0,    0/
1259       DATA ((KFPR(I,J),J=1,2),I=201,240)/
1260      &  1000011,   1000011,   2000011,   2000011,   1000011,
1261      &  2000011,   1000013,   1000013,   2000013,   2000013,
1262      &  1000013,   2000013,   1000015,   1000015,   2000015,
1263      &  2000015,   1000015,   2000015,   1000011,   1000012,
1264      1  1000015,   1000016,   2000015,   1000016,   1000012,
1265      1  1000012,   1000016,   1000016,         0,         0,
1266      1  1000022,   1000022,   1000023,   1000023,   1000025,
1267      1  1000025,   1000035,   1000035,   1000022,   1000023,
1268      2  1000022,   1000025,   1000022,   1000035,   1000023,
1269      2  1000025,   1000023,   1000035,   1000025,   1000035,
1270      2  1000024,   1000024,   1000037,   1000037,   1000024,
1271      2  1000037,   1000022,   1000024,   1000023,   1000024,
1272      3  1000025,   1000024,   1000035,   1000024,   1000022,
1273      3  1000037,   1000023,   1000037,   1000025,   1000037,
1274      3  1000035,   1000037,   1000021,   1000022,   1000021,
1275      3  1000023,   1000021,   1000025,   1000021,   1000035/
1276       DATA ((KFPR(I,J),J=1,2),I=241,280)/
1277      4  1000021,   1000024,   1000021,   1000037,   1000021,
1278      4  1000021,   1000021,   1000021,         0,         0,
1279      4  1000002,   1000022,   2000002,   1000022,   1000002,
1280      4  1000023,   2000002,   1000023,   1000002,   1000025,
1281      5  2000002,   1000025,   1000002,   1000035,   2000002,
1282      5  1000035,   1000001,   1000024,   2000005,   1000024,
1283      5  1000001,   1000037,   2000005,   1000037,   1000002,
1284      5  1000021,   2000002,   1000021,         0,         0,
1285      6  1000006,   1000006,   2000006,   2000006,   1000006,
1286      6  2000006,   1000006,   1000006,   2000006,   2000006,
1287      6        0,         0,         0,         0,         0,
1288      6        0,         0,         0,         0,         0,
1289      7  1000002,   1000002,   2000002,   2000002,   1000002,
1290      7  2000002,   1000002,   1000002,   2000002,   2000002,
1291      7  1000002,   2000002,   1000002,   1000002,   2000002,
1292      7  2000002,   1000002,   1000002,   2000002,   2000002/
1293       DATA ((KFPR(I,J),J=1,2),I=281,500)/440*0/
1294       DATA COEF/10000*0D0/
1295       DATA (((ICOL(I,J,K),K=1,2),J=1,4),I=1,40)/
1296      &4,0,3,0,2,0,1,0,3,0,4,0,1,0,2,0,2,0,0,1,4,0,0,3,3,0,0,4,1,0,0,2,
1297      &3,0,0,4,1,4,3,2,4,0,0,3,4,2,1,3,2,0,4,1,4,0,2,3,4,0,3,4,2,0,1,2,
1298      &3,2,1,0,1,4,3,0,4,3,3,0,2,1,1,0,3,2,1,4,1,0,0,2,2,4,3,1,2,0,0,1,
1299      &3,2,1,4,1,4,3,2,4,2,1,3,4,2,1,3,3,4,4,3,1,2,2,1,2,0,3,1,2,0,0,0,
1300      &4,2,1,0,0,0,1,0,3,0,0,3,1,2,0,0,4,0,0,4,0,0,1,2,2,0,0,1,4,4,3,3,
1301      &2,2,1,1,4,4,3,3,3,3,4,4,1,1,2,2,3,2,1,3,1,2,0,0,4,2,1,4,0,0,1,2,
1302      &4,0,0,0,4,0,1,3,0,0,3,0,2,4,3,0,3,4,0,0,1,0,0,1,0,0,3,4,2,0,0,2,
1303      &3,0,0,0,1,0,0,0,0,0,3,0,2,0,0,0,2,0,3,1,2,0,0,0,3,2,1,0,1,0,0,0,
1304      &4,4,3,3,2,2,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1305      &0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/
1306
1307 C...Treatment of resonances.
1308       DATA (MWID(I)  ,I=   1, 500)/5*0,3*1,8*0,1,5*0,3*1,6*0,1,0,7*1,
1309      &10*0,2*1,0,3*1,245*0,19*2,0,7*2,0,2,0,2,0,4*1,163*0/
1310
1311 C...Character constants: name of processes.
1312       DATA PROC(0)/                    'All included subprocesses   '/
1313       DATA (PROC(I),I=1,20)/
1314      &'f + fbar -> gamma*/Z0       ',  'f + fbar'' -> W+/-           ',
1315      &'f + fbar -> h0              ',  'gamma + W+/- -> W+/-        ',
1316      &'Z0 + Z0 -> h0               ',  'Z0 + W+/- -> W+/-           ',
1317      &'                            ',  'W+ + W- -> h0               ',
1318      &'                            ',  'f + f'' -> f + f'' (QFD)      ',
1319      1'f + f'' -> f + f'' (QCD)      ','f + fbar -> f'' + fbar''      ',
1320      1'f + fbar -> g + g           ',  'f + fbar -> g + gamma       ',
1321      1'f + fbar -> g + Z0          ',  'f + fbar'' -> g + W+/-       ',
1322      1'f + fbar -> g + h0          ',  'f + fbar -> gamma + gamma   ',
1323      1'f + fbar -> gamma + Z0      ',  'f + fbar'' -> gamma + W+/-   '/
1324       DATA (PROC(I),I=21,40)/
1325      2'f + fbar -> gamma + h0      ',  'f + fbar -> Z0 + Z0         ',
1326      2'f + fbar'' -> Z0 + W+/-      ', 'f + fbar -> Z0 + h0         ',
1327      2'f + fbar -> W+ + W-         ',  'f + fbar'' -> W+/- + h0      ',
1328      2'f + fbar -> h0 + h0         ',  'f + g -> f + g              ',
1329      2'f + g -> f + gamma          ',  'f + g -> f + Z0             ',
1330      3'f + g -> f'' + W+/-          ', 'f + g -> f + h0             ',
1331      3'f + gamma -> f + g          ',  'f + gamma -> f + gamma      ',
1332      3'f + gamma -> f + Z0         ',  'f + gamma -> f'' + W+/-      ',
1333      3'f + gamma -> f + h0         ',  'f + Z0 -> f + g             ',
1334      3'f + Z0 -> f + gamma         ',  'f + Z0 -> f + Z0            '/
1335       DATA (PROC(I),I=41,60)/
1336      4'f + Z0 -> f'' + W+/-         ', 'f + Z0 -> f + h0            ',
1337      4'f + W+/- -> f'' + g          ', 'f + W+/- -> f'' + gamma      ',
1338      4'f + W+/- -> f'' + Z0         ', 'f + W+/- -> f'' + W+/-       ',
1339      4'f + W+/- -> f'' + h0         ', 'f + h0 -> f + g             ',
1340      4'f + h0 -> f + gamma         ',  'f + h0 -> f + Z0            ',
1341      5'f + h0 -> f'' + W+/-         ', 'f + h0 -> f + h0            ',
1342      5'g + g -> f + fbar           ',  'g + gamma -> f + fbar       ',
1343      5'g + Z0 -> f + fbar          ',  'g + W+/- -> f + fbar''       ',
1344      5'g + h0 -> f + fbar          ',  'gamma + gamma -> f + fbar   ',
1345      5'gamma + Z0 -> f + fbar      ',  'gamma + W+/- -> f + fbar''   '/
1346       DATA (PROC(I),I=61,80)/
1347      6'gamma + h0 -> f + fbar      ',  'Z0 + Z0 -> f + fbar         ',
1348      6'Z0 + W+/- -> f + fbar''      ', 'Z0 + h0 -> f + fbar         ',
1349      6'W+ + W- -> f + fbar         ',  'W+/- + h0 -> f + fbar''      ',
1350      6'h0 + h0 -> f + fbar         ',  'g + g -> g + g              ',
1351      6'gamma + gamma -> W+ + W-    ',  'gamma + W+/- -> Z0 + W+/-   ',
1352      7'Z0 + Z0 -> Z0 + Z0          ',  'Z0 + Z0 -> W+ + W-          ',
1353      7'Z0 + W+/- -> Z0 + W+/-      ',  'Z0 + Z0 -> Z0 + h0          ',
1354      7'W+ + W- -> gamma + gamma    ',  'W+ + W- -> Z0 + Z0          ',
1355      7'W+/- + W+/- -> W+/- + W+/-  ',  'W+/- + h0 -> W+/- + h0      ',
1356      7'h0 + h0 -> h0 + h0          ',  'q + gamma -> q'' + pi+/-     '/
1357       DATA (PROC(I),I=81,100)/
1358      8'q + qbar -> Q + Qbar, mass  ',  'g + g -> Q + Qbar, massive  ',
1359      8'f + q -> f'' + Q, massive    ', 'g + gamma -> Q + Qbar, mass ',
1360      8'gamma + gamma -> F + Fbar, m',  'g + g -> J/Psi + g          ',
1361      8'g + g -> chi_0c + g         ',  'g + g -> chi_1c + g         ',
1362      8'g + g -> chi_2c + g         ',  '                            ',
1363      9'Elastic scattering          ',  'Single diffractive (XB)     ',
1364      9'Single diffractive (AX)     ',  'Double  diffractive         ',
1365      9'Low-pT scattering           ',  'Semihard QCD 2 -> 2         ',
1366      9'                            ',  '                            ',
1367      9'                            ',  '                            '/
1368       DATA (PROC(I),I=101,120)/
1369      &'g + g -> gamma*/Z0          ',  'g + g -> h0                 ',
1370      &'gamma + gamma -> h0         ',  '                            ',
1371      &'                            ',  'g + g -> J/Psi + gamma      ',
1372      &'gamma + g -> J/Psi + g      ',  'gamma+gamma -> J/Psi + gamma',
1373      &'                            ',  'f + fbar -> gamma + h0      ',
1374      1'f + fbar -> g + h0          ',  'q + g -> q + h0             ',
1375      1'g + g -> g + h0             ',  'g + g -> gamma + gamma      ',
1376      1'g + g -> g + gamma          ',  'g + g -> gamma + Z0         ',
1377      1'g + g -> Z0 + Z0            ',  'g + g -> W+ + W-            ',
1378      1'                            ',  '                            '/
1379       DATA (PROC(I),I=121,140)/
1380      2'g + g -> Q + Qbar + h0      ',  'q + qbar -> Q + Qbar + h0   ',
1381      2'f + f'' -> f + f'' + h0       ',
1382      2'f + f'' -> f" + f"'' + h0     ',
1383      2'                            ',  '                            ',
1384      2'                            ',  '                            ',
1385      2'                            ',  '                            ',
1386      3'g + g -> Z0 + q + qbar      ',  '                            ',
1387      3'                            ',  '                            ',
1388      3'                            ',  '                            ',
1389      3'                            ',  '                            ',
1390      3'                            ',  '                            '/
1391       DATA (PROC(I),I=141,160)/
1392      4'f + fbar -> gamma*/Z0/Z''0   ', 'f + fbar'' -> W''+/-          ',
1393      4'f + fbar'' -> H+/-           ', 'f + fbar'' -> R              ',
1394      4'q + l -> LQ                 ',  '                            ',
1395      4'd + g -> d*                 ',  'u + g -> u*                 ',
1396      4'g + g -> eta_techni         ',  '                            ',
1397      5'f + fbar -> H0              ',  'g + g -> H0                 ',
1398      5'gamma + gamma -> H0         ',  '                            ',
1399      5'                            ',  'f + fbar -> A0              ',
1400      5'g + g -> A0                 ',  'gamma + gamma -> A0         ',
1401      5'                            ',  '                            '/
1402       DATA (PROC(I),I=161,180)/
1403      6'f + g -> f'' + H+/-          ', 'q + g -> LQ + lbar          ',
1404      6'g + g -> LQ + LQbar         ',  'q + qbar -> LQ + LQbar      ',
1405      6'f + fbar -> f'' + fbar'' (g/Z)',
1406      6'f +fbar'' -> f" + fbar"'' (W) ',
1407      6'q + q'' -> q" + d*           ',  'q + q'' -> q" + u*           ',
1408      6'                            ',  '                            ',
1409      7'f + fbar -> Z0 + H0         ', 'f + fbar'' -> W+/- + H0      ',
1410      7'f + f'' -> f + f'' + H0       ',
1411      7'f + f'' -> f" + f"'' + H0     ',
1412      7'                            ',  'f + fbar -> Z0 + A0         ',
1413      7'f + fbar'' -> W+/- + A0      ',
1414      7'f + f'' -> f + f'' + A0       ',
1415      7'f + f'' -> f" + f"'' + A0     ',
1416      7'                            '/
1417       DATA (PROC(I),I=181,200)/
1418      8'g + g -> Q + Qbar + H0      ',  'q + qbar -> Q + Qbar + H0   ',
1419      8'                            ',  '                            ',
1420      8'                            ',  'g + g -> Q + Qbar + A0      ',
1421      8'q + qbar -> Q + Qbar + A0   ',  '                            ',
1422      8'                            ',  '                            ',
1423      9'f + fbar -> rho_tech0       ',  'f + f'' -> rho_tech+/-       ',
1424      9'f + fbar -> omega_tech0     ',  'f+fbar -> f''+fbar'' (technic)',
1425      9'                            ',  '                            ',
1426      9'                            ',  '                            ',
1427      9'                            ',  '                            '/
1428       DATA (PROC(I),I=201,220)/
1429      &'f + fbar -> ~e_L + ~e_Lbar  ',  'f + fbar -> ~e_R + ~e_Rbar  ',
1430      &'f + fbar -> ~e_R + ~e_Lbar  ',  'f + fbar -> ~mu_L + ~mu_Lbar',
1431      &'f + fbar -> ~mu_R + ~mu_Rbar',  'f + fbar -> ~mu_L + ~mu_Rbar',
1432      &'f+fbar -> ~tau_1 + ~tau_1bar',  'f+fbar -> ~tau_2 + ~tau_2bar',
1433      &'f+fbar -> ~tau_1 + ~tau_2bar',  'q + qbar'' -> ~l_L + ~nulbar ',
1434      1'q+qbar''-> ~tau_1 + ~nutaubar', 'q+qbar''-> ~tau_2 + ~nutaubar',
1435      1'f + fbar -> ~nul + ~nulbar  ',  'f+fbar -> ~nutau + ~nutaubar',
1436      1'                            ',  'f + fbar -> ~chi1 + ~chi1   ',
1437      1'f + fbar -> ~chi2 + ~chi2   ',  'f + fbar -> ~chi3 + ~chi3   ',
1438      1'f + fbar -> ~chi4 + ~chi4   ',  'f + fbar -> ~chi1 + ~chi2   '/
1439       DATA (PROC(I),I=221,240)/
1440      2'f + fbar -> ~chi1 + ~chi3   ',  'f + fbar -> ~chi1 + ~chi4   ',
1441      2'f + fbar -> ~chi2 + ~chi3   ',  'f + fbar -> ~chi2 + ~chi4   ',
1442      2'f + fbar -> ~chi3 + ~chi4   ',  'f+fbar -> ~chi+-1 + ~chi-+1 ',
1443      2'f+fbar -> ~chi+-2 + ~chi-+2 ',  'f+fbar -> ~chi+-1 + ~chi-+2 ',
1444      2'q + qbar'' -> ~chi1 + ~chi+-1', 'q + qbar'' -> ~chi2 + ~chi+-1',
1445      3'q + qbar'' -> ~chi3 + ~chi+-1', 'q + qbar'' -> ~chi4 + ~chi+-1',
1446      3'q + qbar'' -> ~chi1 + ~chi+-2', 'q + qbar'' -> ~chi2 + ~chi+-2',
1447      3'q + qbar'' -> ~chi3 + ~chi+-2', 'q + qbar'' -> ~chi4 + ~chi+-2',
1448      3'q + qbar -> ~chi1 + ~g      ',  'q + qbar -> ~chi2 + ~g      ',
1449      3'q + qbar -> ~chi3 + ~g      ',  'q + qbar -> ~chi4 + ~g      '/
1450       DATA (PROC(I),I=241,260)/
1451      4'q + qbar'' -> ~chi+-1 + ~g   ', 'q + qbar'' -> ~chi+-2 + ~g  ',
1452      4'q + qbar -> ~g + ~g         ',  'g + g -> ~g + ~g            ',
1453      4'                            ',  'qj + g -> ~qj_L + ~chi1     ',
1454      4'qj + g -> ~qj_R + ~chi1     ',  'qj + g -> ~qj_L + ~chi2     ',
1455      4'qj + g -> ~qj_R + ~chi2     ',  'qj + g -> ~qj_L + ~chi3     ',
1456      5'qj + g -> ~qj_R + ~chi3     ',  'qj + g -> ~qj_L + ~chi4     ',
1457      5'qj + g -> ~qj_R + ~chi4     ',  'qj + g -> ~qk_L + ~chi+-1   ',
1458      5'qj + g -> ~qk_R + ~chi+-1   ',  'qj + g -> ~qk_L + ~chi+-2   ',
1459      5'qj + g -> ~qk_R + ~chi+-2   ',  'qj + g -> ~qj_L + ~g        ',
1460      5'qj + g -> ~qj_R + ~g        ',  '                            '/
1461       DATA (PROC(I),I=261,280)/
1462      6'f + fbar -> ~t_1 + ~t_1bar  ',  'f + fbar -> ~t_2 + ~t_2bar  ',
1463      6'f + fbar -> ~t_1 + ~t_2bar  ',  'g + g -> ~t_1 + ~t_1bar     ',
1464      6'g + g -> ~t_2 + ~t_2bar     ',  '                            ',
1465      6'                            ',  '                            ',
1466      6'                            ',  '                            ',
1467      7'qi + qj -> ~qi_L + ~qj_L    ',  'qi + qj -> ~qi_R + ~qj_R    ',
1468      7'qi + qj -> ~qi_L + ~qj_R    ',  'qi+qjbar -> ~qi_L + ~qj_Lbar',
1469      7'qi+qjbar -> ~qi_R + ~qj_Rbar',  'qi+qjbar -> ~qi_L + ~qj_Rbar',
1470      7'f + fbar -> ~qi_L + ~qi_Lbar',  'f + fbar -> ~qi_R + ~qi_Rbar',
1471      7'g + g -> ~qi_L + ~qi_Lbar   ',  'g + g -> ~qi_R + ~qi_Rbar   '/
1472       DATA (PROC(I),I=281,500)/220*'                            '/
1473
1474 C...Cross sections and slope offsets.
1475       DATA SIGT/294*0D0/
1476
1477 C...Supersymmetry switches and parameters.
1478       DATA IMSS/0,
1479      &  0,  0,  0,  1,  0,  0,  0,  1,  0,  0,
1480      1  89*0/
1481       DATA RMSS/0D0,
1482      &  80D0,160D0,500D0,800D0,2D0,250D0,200D0,800D0,700D0,800D0,
1483      1  700D0,500D0,250D0,200D0,800D0,400D0,0D0,0.1D0,850D0,0.041D0,
1484      2   1D0,800D0,1D4,1D4,1D4,0D0,0D0,24D17,2*0D0,
1485      3  69*0D0/
1486
1487 C...Data for histogramming routines.
1488       DATA IHIST/1000,20000,55,1/
1489       DATA INDX/1000*0/
1490
1491       END
1492
1493 C*********************************************************************
1494
1495 C...PYTEST
1496 C...A simple program (disguised as subroutine) to run at installation
1497 C...as a check that the program works as intended.
1498
1499       SUBROUTINE PYTEST(MTEST)
1500
1501 C...Double precision and integer declarations.
1502       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1503       INTEGER PYK,PYCHGE,PYCOMP
1504 C...Commonblocks.
1505       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
1506       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
1507       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
1508       COMMON/PYDAT3/MDCY(500,3),MDME(4000,2),BRAT(4000),KFDP(4000,5)
1509       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
1510       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
1511       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYDAT3/,/PYSUBS/,/PYPARS/
1512 C...Local arrays.
1513       DIMENSION PSUM(5),PINI(6),PFIN(6)
1514
1515 C...Save defaults for values that are changed.
1516       MSTJ1=MSTJ(1)
1517       MSTJ3=MSTJ(3)
1518       MSTJ11=MSTJ(11)
1519       MSTJ42=MSTJ(42)
1520       MSTJ43=MSTJ(43)
1521       MSTJ44=MSTJ(44)
1522       PARJ17=PARJ(17)
1523       PARJ22=PARJ(22)
1524       PARJ43=PARJ(43)
1525       PARJ54=PARJ(54)
1526       MST101=MSTJ(101)
1527       MST104=MSTJ(104)
1528       MST105=MSTJ(105)
1529       MST107=MSTJ(107)
1530       MST116=MSTJ(116)
1531
1532 C...First part: loop over simple events to be generated.
1533       IF(MTEST.GE.1) CALL PYTABU(20)
1534       NERR=0
1535       DO 180 IEV=1,500
1536
1537 C...Reset parameter values. Switch on some nonstandard features.
1538         MSTJ(1)=1
1539         MSTJ(3)=0
1540         MSTJ(11)=1
1541         MSTJ(42)=2
1542         MSTJ(43)=4
1543         MSTJ(44)=2
1544         PARJ(17)=0.1D0
1545         PARJ(22)=1.5D0
1546         PARJ(43)=1D0
1547         PARJ(54)=-0.05D0
1548         MSTJ(101)=5
1549         MSTJ(104)=5
1550         MSTJ(105)=0
1551         MSTJ(107)=1
1552         IF(IEV.EQ.301.OR.IEV.EQ.351.OR.IEV.EQ.401) MSTJ(116)=3
1553
1554 C...Ten events each for some single jets configurations.
1555         IF(IEV.LE.50) THEN
1556           ITY=(IEV+9)/10
1557           MSTJ(3)=-1
1558           IF(ITY.EQ.3.OR.ITY.EQ.4) MSTJ(11)=2
1559           IF(ITY.EQ.1) CALL PY1ENT(1,1,15D0,0D0,0D0)
1560           IF(ITY.EQ.2) CALL PY1ENT(1,3101,15D0,0D0,0D0)
1561           IF(ITY.EQ.3) CALL PY1ENT(1,-2203,15D0,0D0,0D0)
1562           IF(ITY.EQ.4) CALL PY1ENT(1,-4,30D0,0D0,0D0)
1563           IF(ITY.EQ.5) CALL PY1ENT(1,21,15D0,0D0,0D0)
1564
1565 C...Ten events each for some simple jet systems; string fragmentation.
1566         ELSEIF(IEV.LE.130) THEN
1567           ITY=(IEV-41)/10
1568           IF(ITY.EQ.1) CALL PY2ENT(1,1,-1,40D0)
1569           IF(ITY.EQ.2) CALL PY2ENT(1,4,-4,30D0)
1570           IF(ITY.EQ.3) CALL PY2ENT(1,2,2103,100D0)
1571           IF(ITY.EQ.4) CALL PY2ENT(1,21,21,40D0)
1572           IF(ITY.EQ.5) CALL PY3ENT(1,2101,21,-3203,30D0,0.6D0,0.8D0)
1573           IF(ITY.EQ.6) CALL PY3ENT(1,5,21,-5,40D0,0.9D0,0.8D0)
1574           IF(ITY.EQ.7) CALL PY3ENT(1,21,21,21,60D0,0.7D0,0.5D0)
1575           IF(ITY.EQ.8) CALL PY4ENT(1,2,21,21,-2,40D0,
1576      &    0.4D0,0.64D0,0.6D0,0.12D0,0.2D0)
1577
1578 C...Seventy events with independent fragmentation and momentum cons.
1579         ELSEIF(IEV.LE.200) THEN
1580           ITY=1+(IEV-131)/16
1581           MSTJ(2)=1+MOD(IEV-131,4)
1582           MSTJ(3)=1+MOD((IEV-131)/4,4)
1583           IF(ITY.EQ.1) CALL PY2ENT(1,4,-5,40D0)
1584           IF(ITY.EQ.2) CALL PY3ENT(1,3,21,-3,40D0,0.9D0,0.4D0)
1585           IF(ITY.EQ.3) CALL PY4ENT(1,2,21,21,-2,40D0,
1586      &    0.4D0,0.64D0,0.6D0,0.12D0,0.2D0)
1587           IF(ITY.GE.4) CALL PY4ENT(1,2,-3,3,-2,40D0,
1588      &    0.4D0,0.64D0,0.6D0,0.12D0,0.2D0)
1589
1590 C...A hundred events with random jets (check invariant mass).
1591         ELSEIF(IEV.LE.300) THEN
1592   100     DO 110 J=1,5
1593             PSUM(J)=0D0
1594   110     CONTINUE
1595           NJET=2D0+6D0*PYR(0)
1596           DO 130 I=1,NJET
1597             KFL=21
1598             IF(I.EQ.1) KFL=INT(1D0+4D0*PYR(0))
1599             IF(I.EQ.NJET) KFL=-INT(1D0+4D0*PYR(0))
1600             EJET=5D0+20D0*PYR(0)
1601             THETA=ACOS(2D0*PYR(0)-1D0)
1602             PHI=6.2832D0*PYR(0)
1603             IF(I.LT.NJET) CALL PY1ENT(-I,KFL,EJET,THETA,PHI)
1604             IF(I.EQ.NJET) CALL PY1ENT(I,KFL,EJET,THETA,PHI)
1605             IF(I.EQ.1.OR.I.EQ.NJET) MSTJ(93)=1
1606             IF(I.EQ.1.OR.I.EQ.NJET) PSUM(5)=PSUM(5)+PYMASS(KFL)
1607             DO 120 J=1,4
1608               PSUM(J)=PSUM(J)+P(I,J)
1609   120       CONTINUE
1610   130     CONTINUE
1611           IF(PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2.LT.
1612      &    (PSUM(5)+PARJ(32))**2) GOTO 100
1613
1614 C...Fifty e+e- continuum events with matrix elements.
1615         ELSEIF(IEV.LE.350) THEN
1616           MSTJ(101)=2
1617           CALL PYEEVT(0,40D0)
1618
1619 C...Fifty e+e- continuum event with varying shower options.
1620         ELSEIF(IEV.LE.400) THEN
1621           MSTJ(42)=1+MOD(IEV,2)
1622           MSTJ(43)=1+MOD(IEV/2,4)
1623           MSTJ(44)=MOD(IEV/8,3)
1624           CALL PYEEVT(0,90D0)
1625
1626 C...Fifty e+e- continuum events with coherent shower.
1627         ELSEIF(IEV.LE.450) THEN
1628           CALL PYEEVT(0,500D0)
1629
1630 C...Fifty Upsilon decays to ggg or gammagg with coherent shower.
1631         ELSE
1632           CALL PYONIA(5,9.46D0)
1633         ENDIF
1634
1635 C...Generate event. Find total momentum, energy and charge.
1636         DO 140 J=1,4
1637           PINI(J)=PYP(0,J)
1638   140   CONTINUE
1639         PINI(6)=PYP(0,6)
1640         CALL PYEXEC
1641         DO 150 J=1,4
1642           PFIN(J)=PYP(0,J)
1643   150   CONTINUE
1644         PFIN(6)=PYP(0,6)
1645
1646 C...Check conservation of energy, momentum and charge;
1647 C...usually exact, but only approximate for single jets.
1648         MERR=0
1649         IF(IEV.LE.50) THEN
1650           IF((PFIN(1)-PINI(1))**2+(PFIN(2)-PINI(2))**2.GE.4D0)
1651      &    MERR=MERR+1
1652           EPZREM=PINI(4)+PINI(3)-PFIN(4)-PFIN(3)
1653           IF(EPZREM.LT.0D0.OR.EPZREM.GT.2D0*PARJ(31)) MERR=MERR+1
1654           IF(ABS(PFIN(6)-PINI(6)).GT.2.1D0) MERR=MERR+1
1655         ELSE
1656           DO 160 J=1,4
1657             IF(ABS(PFIN(J)-PINI(J)).GT.0.0001D0*PINI(4)) MERR=MERR+1
1658   160     CONTINUE
1659           IF(ABS(PFIN(6)-PINI(6)).GT.0.1D0) MERR=MERR+1
1660         ENDIF
1661         IF(MERR.NE.0) WRITE(MSTU(11),5000) (PINI(J),J=1,4),PINI(6),
1662      &  (PFIN(J),J=1,4),PFIN(6)
1663
1664 C...Check that all KF codes are known ones, and that partons/particles
1665 C...satisfy energy-momentum-mass relation. Store particle statistics.
1666         DO 170 I=1,N
1667           IF(K(I,1).GT.20) GOTO 170
1668           IF(PYCOMP(K(I,2)).EQ.0) THEN
1669             WRITE(MSTU(11),5100) I
1670             MERR=MERR+1
1671           ENDIF
1672           PD=P(I,4)**2-P(I,1)**2-P(I,2)**2-P(I,3)**2-P(I,5)**2
1673           IF(ABS(PD).GT.MAX(0.1D0,0.001D0*P(I,4)**2).OR.P(I,4).LT.0D0)
1674      &    THEN
1675             WRITE(MSTU(11),5200) I
1676             MERR=MERR+1
1677           ENDIF
1678   170   CONTINUE
1679         IF(MTEST.GE.1) CALL PYTABU(21)
1680
1681 C...List all erroneous events and some normal ones.
1682         IF(MERR.NE.0.OR.MSTU(24).NE.0.OR.MSTU(28).NE.0) THEN
1683           IF(MERR.GE.1) WRITE(MSTU(11),6400)
1684           CALL PYLIST(2)
1685         ELSEIF(MTEST.GE.1.AND.MOD(IEV-5,100).EQ.0) THEN
1686           CALL PYLIST(1)
1687         ENDIF
1688
1689 C...Stop execution if too many errors.
1690         IF(MERR.NE.0) NERR=NERR+1
1691         IF(NERR.GE.10) THEN
1692           WRITE(MSTU(11),6300)
1693           CALL PYLIST(1)
1694           STOP
1695         ENDIF
1696   180 CONTINUE
1697
1698 C...Summarize result of run.
1699       IF(MTEST.GE.1) CALL PYTABU(22)
1700
1701 C...Reset commonblock variables changed during run.
1702       MSTJ(1)=MSTJ1
1703       MSTJ(3)=MSTJ3
1704       MSTJ(11)=MSTJ11
1705       MSTJ(42)=MSTJ42
1706       MSTJ(43)=MSTJ43
1707       MSTJ(44)=MSTJ44
1708       PARJ(17)=PARJ17
1709       PARJ(22)=PARJ22
1710       PARJ(43)=PARJ43
1711       PARJ(54)=PARJ54
1712       MSTJ(101)=MST101
1713       MSTJ(104)=MST104
1714       MSTJ(105)=MST105
1715       MSTJ(107)=MST107
1716       MSTJ(116)=MST116
1717
1718 C...Second part: complete events of various kinds.
1719 C...Common initial values. Loop over initiating conditions.
1720       MSTP(122)=MAX(0,MIN(2,MTEST))
1721       MDCY(PYCOMP(111),1)=0
1722       DO 230 IPROC=1,8
1723
1724 C...Reset process type, kinematics cuts, and the flags used.
1725         MSEL=0
1726         DO 190 ISUB=1,500
1727           MSUB(ISUB)=0
1728   190   CONTINUE
1729         CKIN(1)=2D0
1730         CKIN(3)=0D0
1731         MSTP(2)=1
1732         MSTP(11)=0
1733         MSTP(33)=0
1734         MSTP(81)=1
1735         MSTP(82)=1
1736         MSTP(111)=1
1737         MSTP(131)=0
1738         MSTP(133)=0
1739         PARP(131)=0.01D0
1740
1741 C...Prompt photon production at fixed target.
1742         IF(IPROC.EQ.1) THEN
1743           PZSUM=300D0
1744           PESUM=SQRT(PZSUM**2+PYMASS(211)**2)+PYMASS(2212)
1745           PQSUM=2D0
1746           MSEL=10
1747           CKIN(3)=5D0
1748           CALL PYINIT('FIXT','pi+','p',PZSUM)
1749
1750 C...QCD processes at ISR energies.
1751         ELSEIF(IPROC.EQ.2) THEN
1752           PESUM=63D0
1753           PZSUM=0D0
1754           PQSUM=2D0
1755           MSEL=1
1756           CKIN(3)=5D0
1757           CALL PYINIT('CMS','p','p',PESUM)
1758
1759 C...W production + multiple interactions at CERN Collider.
1760         ELSEIF(IPROC.EQ.3) THEN
1761           PESUM=630D0
1762           PZSUM=0D0
1763           PQSUM=0D0
1764           MSEL=12
1765           CKIN(1)=20D0
1766           MSTP(82)=4
1767           MSTP(2)=2
1768           MSTP(33)=3
1769           CALL PYINIT('CMS','p','pbar',PESUM)
1770
1771 C...W/Z gauge boson pairs + pileup events at the Tevatron.
1772         ELSEIF(IPROC.EQ.4) THEN
1773           PESUM=1800D0
1774           PZSUM=0D0
1775           PQSUM=0D0
1776           MSUB(22)=1
1777           MSUB(23)=1
1778           MSUB(25)=1
1779           CKIN(1)=200D0
1780           MSTP(111)=0
1781           MSTP(131)=1
1782           MSTP(133)=2
1783           PARP(131)=0.04D0
1784           CALL PYINIT('CMS','p','pbar',PESUM)
1785
1786 C...Higgs production at LHC.
1787         ELSEIF(IPROC.EQ.5) THEN
1788           PESUM=15400D0
1789           PZSUM=0D0
1790           PQSUM=2D0
1791           MSUB(3)=1
1792           MSUB(102)=1
1793           MSUB(123)=1
1794           MSUB(124)=1
1795           PMAS(25,1)=300D0
1796           CKIN(1)=200D0
1797           MSTP(81)=0
1798           MSTP(111)=0
1799           CALL PYINIT('CMS','p','p',PESUM)
1800
1801 C...Z' production at SSC.
1802         ELSEIF(IPROC.EQ.6) THEN
1803           PESUM=40000D0
1804           PZSUM=0D0
1805           PQSUM=2D0
1806           MSEL=21
1807           PMAS(32,1)=600D0
1808           CKIN(1)=400D0
1809           MSTP(81)=0
1810           MSTP(111)=0
1811           CALL PYINIT('CMS','p','p',PESUM)
1812
1813 C...W pair production at 1 TeV e+e- collider.
1814         ELSEIF(IPROC.EQ.7) THEN
1815           PESUM=1000D0
1816           PZSUM=0D0
1817           PQSUM=0D0
1818           MSUB(25)=1
1819           MSUB(69)=1
1820           MSTP(11)=1
1821           CALL PYINIT('CMS','e+','e-',PESUM)
1822
1823 C...Deep inelastic scattering at a LEP+LHC ep collider.
1824         ELSEIF(IPROC.EQ.8) THEN
1825           P(1,1)=0D0
1826           P(1,2)=0D0
1827           P(1,3)=8000D0
1828           P(2,1)=0D0
1829           P(2,2)=0D0
1830           P(2,3)=-80D0
1831           PESUM=8080D0
1832           PZSUM=7920D0
1833           PQSUM=0D0
1834           MSUB(10)=1
1835           CKIN(3)=50D0
1836           MSTP(111)=0
1837           CALL PYINIT('USER','p','e-',PESUM)
1838         ENDIF
1839
1840 C...Generate 20 events of each required type.
1841         DO 220 IEV=1,20
1842           CALL PYEVNT
1843           PESUMM=PESUM
1844           IF(IPROC.EQ.4) PESUMM=MSTI(41)*PESUM
1845
1846 C...Check conservation of energy/momentum/flavour.
1847           PINI(1)=0D0
1848           PINI(2)=0D0
1849           PINI(3)=PZSUM
1850           PINI(4)=PESUMM
1851           PINI(6)=PQSUM
1852           DO 200 J=1,4
1853             PFIN(J)=PYP(0,J)
1854   200     CONTINUE
1855           PFIN(6)=PYP(0,6)
1856           MERR=0
1857           DEVE=ABS(PFIN(4)-PINI(4))+ABS(PFIN(3)-PINI(3))
1858           DEVT=ABS(PFIN(1)-PINI(1))+ABS(PFIN(2)-PINI(2))
1859           DEVQ=ABS(PFIN(6)-PINI(6))
1860           IF(DEVE.GT.2D-3*PESUM.OR.DEVT.GT.MAX(0.01D0,1D-4*PESUM).OR.
1861      &    DEVQ.GT.0.1D0) MERR=1
1862           IF(MERR.NE.0) WRITE(MSTU(11),5000) (PINI(J),J=1,4),PINI(6),
1863      &    (PFIN(J),J=1,4),PFIN(6)
1864
1865 C...Check that all KF codes are known ones, and that partons/particles
1866 C...satisfy energy-momentum-mass relation.
1867           DO 210 I=1,N
1868             IF(K(I,1).GT.20) GOTO 210
1869             IF(PYCOMP(K(I,2)).EQ.0) THEN
1870               WRITE(MSTU(11),5100) I
1871               MERR=MERR+1
1872             ENDIF
1873             PD=P(I,4)**2-P(I,1)**2-P(I,2)**2-P(I,3)**2-P(I,5)**2*
1874      &      SIGN(1D0,P(I,5))
1875             IF(ABS(PD).GT.MAX(0.1D0,0.002D0*P(I,4)**2,0.002D0*P(I,5)**2)
1876      &      .OR.(P(I,5).GE.0D0.AND.P(I,4).LT.0D0)) THEN
1877               WRITE(MSTU(11),5200) I
1878               MERR=MERR+1
1879             ENDIF
1880   210     CONTINUE
1881
1882 C...Listing of erroneous events, and first event of each type.
1883           IF(MERR.GE.1) NERR=NERR+1
1884           IF(NERR.GE.10) THEN
1885             WRITE(MSTU(11),6300)
1886             CALL PYLIST(1)
1887             STOP
1888           ENDIF
1889           IF(MTEST.GE.1.AND.(MERR.GE.1.OR.IEV.EQ.1)) THEN
1890             IF(MERR.GE.1) WRITE(MSTU(11),6400)
1891             CALL PYLIST(1)
1892           ENDIF
1893   220   CONTINUE
1894
1895 C...List statistics for each process type.
1896         IF(MTEST.GE.1) CALL PYSTAT(1)
1897   230 CONTINUE
1898
1899 C...Summarize result of run.
1900       IF(NERR.EQ.0) WRITE(MSTU(11),6500)
1901       IF(NERR.GT.0) WRITE(MSTU(11),6600) NERR
1902
1903 C...Format statements for output.
1904  5000 FORMAT(/' Momentum, energy and/or charge were not conserved ',
1905      &'in following event'/' sum of',9X,'px',11X,'py',11X,'pz',11X,
1906      &'E',8X,'charge'/' before',2X,4(1X,F12.5),1X,F8.2/' after',3X,
1907      &4(1X,F12.5),1X,F8.2)
1908  5100 FORMAT(/5X,'Entry no.',I4,' in following event not known code')
1909  5200 FORMAT(/5X,'Entry no.',I4,' in following event has faulty ',
1910      &'kinematics')
1911  6300 FORMAT(/5X,'This is the tenth error experienced! Something is ',
1912      &'wrong.'/5X,'Execution will be stopped after listing of event.')
1913  6400 FORMAT(5X,'Faulty event follows:')
1914  6500 FORMAT(//5X,'End result of PYTEST: no errors detected.')
1915  6600 FORMAT(//5X,'End result of PYTEST:',I2,' errors detected.'/
1916      &5X,'This should not have happened!')
1917
1918       RETURN
1919       END
1920
1921 C*********************************************************************
1922
1923 C...PYHEPC
1924 C...Converts PYTHIA event record contents to or from
1925 C...the standard event record commonblock.
1926
1927       SUBROUTINE PYHEPC(MCONV)
1928
1929 C...Double precision and integer declarations.
1930       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1931       INTEGER PYK,PYCHGE,PYCOMP
1932 C...Commonblocks.
1933       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
1934       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
1935       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
1936       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/
1937 C...HEPEVT commonblock.
1938       PARAMETER (NMXHEP=4000)
1939       COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
1940      &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
1941       DOUBLE PRECISION PHEP,VHEP
1942       SAVE /HEPEVT/
1943
1944 C...Conversion from PYTHIA to standard, the easy part.
1945       IF(MCONV.EQ.1) THEN
1946         NEVHEP=0
1947         IF(N.GT.NMXHEP) CALL PYERRM(8,
1948      &  '(PYHEPC:) no more space in /HEPEVT/')
1949         NHEP=MIN(N,NMXHEP)
1950         DO 140 I=1,NHEP
1951           ISTHEP(I)=0
1952           IF(K(I,1).GE.1.AND.K(I,1).LE.10) ISTHEP(I)=1
1953           IF(K(I,1).GE.11.AND.K(I,1).LE.20) ISTHEP(I)=2
1954           IF(K(I,1).GE.21.AND.K(I,1).LE.30) ISTHEP(I)=3
1955           IF(K(I,1).GE.31.AND.K(I,1).LE.100) ISTHEP(I)=K(I,1)
1956           IDHEP(I)=K(I,2)
1957           JMOHEP(1,I)=K(I,3)
1958           JMOHEP(2,I)=0
1959           IF(K(I,1).NE.3.AND.K(I,1).NE.13.AND.K(I,1).NE.14) THEN
1960             JDAHEP(1,I)=K(I,4)
1961             JDAHEP(2,I)=K(I,5)
1962           ELSE
1963             JDAHEP(1,I)=0
1964             JDAHEP(2,I)=0
1965           ENDIF
1966           DO 100 J=1,5
1967             PHEP(J,I)=P(I,J)
1968   100     CONTINUE
1969           DO 110 J=1,4
1970             VHEP(J,I)=V(I,J)
1971   110     CONTINUE
1972
1973 C...Check if new event (from pileup).
1974           IF(I.EQ.1) THEN
1975             INEW=1
1976           ELSE
1977             IF(K(I,1).EQ.21.AND.K(I-1,1).NE.21) INEW=I
1978           ENDIF
1979
1980 C...Fill in missing mother information.
1981           IF(I.GE.INEW+2.AND.K(I,1).EQ.21.AND.K(I,3).EQ.0) THEN
1982             IMO1=I-2
1983             IF(I.GE.INEW+3.AND.K(I-1,1).EQ.21.AND.K(I-1,3).EQ.0)
1984      &      IMO1=IMO1-1
1985             JMOHEP(1,I)=IMO1
1986             JMOHEP(2,I)=IMO1+1
1987           ELSEIF(K(I,2).GE.91.AND.K(I,2).LE.93) THEN
1988             I1=K(I,3)-1
1989   120       I1=I1+1
1990             IF(I1.GE.I) CALL PYERRM(8,
1991      &      '(PYHEPC:) translation of inconsistent event history')
1992             IF(I1.LT.I.AND.K(I1,1).NE.1.AND.K(I1,1).NE.11) GOTO 120
1993             KC=PYCOMP(K(I1,2))
1994             IF(I1.LT.I.AND.KC.EQ.0) GOTO 120
1995             IF(I1.LT.I.AND.KCHG(KC,2).EQ.0) GOTO 120
1996             JMOHEP(2,I)=I1
1997           ELSEIF(K(I,2).EQ.94) THEN
1998             NJET=2
1999             IF(NHEP.GE.I+3.AND.K(I+3,3).LE.I) NJET=3
2000             IF(NHEP.GE.I+4.AND.K(I+4,3).LE.I) NJET=4
2001             JMOHEP(2,I)=MOD(K(I+NJET,4)/MSTU(5),MSTU(5))
2002             IF(JMOHEP(2,I).EQ.JMOHEP(1,I)) JMOHEP(2,I)=
2003      &      MOD(K(I+1,4)/MSTU(5),MSTU(5))
2004           ENDIF
2005
2006 C...Fill in missing daughter information.
2007           IF(K(I,2).EQ.94.AND.MSTU(16).NE.2) THEN
2008             DO 130 I1=JDAHEP(1,I),JDAHEP(2,I)
2009               I2=MOD(K(I1,4)/MSTU(5),MSTU(5))
2010               JDAHEP(1,I2)=I
2011   130       CONTINUE
2012           ENDIF
2013           IF(K(I,2).GE.91.AND.K(I,2).LE.94) GOTO 140
2014           I1=JMOHEP(1,I)
2015           IF(I1.LE.0.OR.I1.GT.NHEP) GOTO 140
2016           IF(K(I1,1).NE.13.AND.K(I1,1).NE.14) GOTO 140
2017           IF(JDAHEP(1,I1).EQ.0) THEN
2018             JDAHEP(1,I1)=I
2019           ELSE
2020             JDAHEP(2,I1)=I
2021           ENDIF
2022   140   CONTINUE
2023         DO 150 I=1,NHEP
2024           IF(K(I,1).NE.13.AND.K(I,1).NE.14) GOTO 150
2025           IF(JDAHEP(2,I).EQ.0) JDAHEP(2,I)=JDAHEP(1,I)
2026   150   CONTINUE
2027
2028 C...Conversion from standard to PYTHIA, the easy part.
2029       ELSE
2030         IF(NHEP.GT.MSTU(4)) CALL PYERRM(8,
2031      &  '(PYHEPC:) no more space in /PYJETS/')
2032         N=MIN(NHEP,MSTU(4))
2033         NKQ=0
2034         KQSUM=0
2035         DO 180 I=1,N
2036           K(I,1)=0
2037           IF(ISTHEP(I).EQ.1) K(I,1)=1
2038           IF(ISTHEP(I).EQ.2) K(I,1)=11
2039           IF(ISTHEP(I).EQ.3) K(I,1)=21
2040           K(I,2)=IDHEP(I)
2041           K(I,3)=JMOHEP(1,I)
2042           K(I,4)=JDAHEP(1,I)
2043           K(I,5)=JDAHEP(2,I)
2044           DO 160 J=1,5
2045             P(I,J)=PHEP(J,I)
2046   160     CONTINUE
2047           DO 170 J=1,4
2048             V(I,J)=VHEP(J,I)
2049   170     CONTINUE
2050           V(I,5)=0D0
2051           IF(ISTHEP(I).EQ.2.AND.PHEP(4,I).GT.PHEP(5,I)) THEN
2052             I1=JDAHEP(1,I)
2053             IF(I1.GT.0.AND.I1.LE.NHEP) V(I,5)=(VHEP(4,I1)-VHEP(4,I))*
2054      &      PHEP(5,I)/PHEP(4,I)
2055           ENDIF
2056
2057 C...Fill in missing information on colour connection in jet systems.
2058           IF(ISTHEP(I).EQ.1) THEN
2059             KC=PYCOMP(K(I,2))
2060             KQ=0
2061             IF(KC.NE.0) KQ=KCHG(KC,2)*ISIGN(1,K(I,2))
2062             IF(KQ.NE.0) NKQ=NKQ+1
2063             IF(KQ.NE.2) KQSUM=KQSUM+KQ
2064             IF(KQ.NE.0.AND.KQSUM.NE.0) THEN
2065               K(I,1)=2
2066             ELSEIF(KQ.EQ.2.AND.I.LT.N) THEN
2067               IF(K(I+1,2).EQ.21) K(I,1)=2
2068             ENDIF
2069           ENDIF
2070   180   CONTINUE
2071         IF(NKQ.EQ.1.OR.KQSUM.NE.0) CALL PYERRM(8,
2072      &  '(PYHEPC:) input parton configuration not colour singlet')
2073       ENDIF
2074
2075       END
2076
2077 C*********************************************************************
2078
2079 C...PYINIT
2080 C...Initializes the generation procedure; finds maxima of the
2081 C...differential cross-sections to be used for weighting.
2082
2083       SUBROUTINE PYINIT(FRAME,BEAM,TARGET,WIN)
2084
2085 C...Double precision and integer declarations.
2086       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
2087       INTEGER PYK,PYCHGE,PYCOMP
2088 C...Commonblocks.
2089       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
2090       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
2091       COMMON/PYDAT3/MDCY(500,3),MDME(4000,2),BRAT(4000),KFDP(4000,5)
2092       COMMON/PYDAT4/CHAF(500,2)
2093       CHARACTER CHAF*16
2094       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
2095       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
2096       COMMON/PYINT1/MINT(400),VINT(400)
2097       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
2098       COMMON/PYINT5/NGENPD,NGEN(0:500,3),XSEC(0:500,3)
2099       SAVE /PYDAT1/,/PYDAT2/,/PYDAT3/,/PYDAT4/,/PYSUBS/,/PYPARS/,
2100      &/PYINT1/,/PYINT2/,/PYINT5/
2101 C...Local arrays and character variables.
2102       DIMENSION ALAMIN(20),NFIN(20)
2103       CHARACTER*(*) FRAME,BEAM,TARGET
2104       CHARACTER CHFRAM*8,CHBEAM*8,CHTARG*8,CHLH(2)*6
2105
2106 C...Interface to PDFLIB.
2107       COMMON/W50512/QCDL4,QCDL5
2108       SAVE /W50512/
2109       DOUBLE PRECISION VALUE(20),QCDL4,QCDL5
2110       CHARACTER*20 PARM(20)
2111       DATA VALUE/20*0D0/,PARM/20*' '/
2112
2113 C...Data:Lambda and n_f values for parton distributions; months.
2114       DATA ALAMIN/0.177D0,0.239D0,0.247D0,0.2322D0,0.248D0,0.248D0,
2115      &14*0.2D0/,NFIN/20*4/
2116       DATA CHLH/'lepton','hadron'/
2117
2118 C...Reset MINT and VINT arrays. Write headers.
2119       DO 100 J=1,400
2120         MINT(J)=0
2121         VINT(J)=0D0
2122   100 CONTINUE
2123       IF(MSTU(12).GE.1) CALL PYLIST(0)
2124       IF(MSTP(122).GE.1) WRITE(MSTU(11),5100)
2125
2126 C...Maximum 4 generations; set maximum number of allowed flavours.
2127       MSTP(1)=MIN(4,MSTP(1))
2128       MSTU(114)=MIN(MSTU(114),2*MSTP(1))
2129       MSTP(58)=MIN(MSTP(58),2*MSTP(1))
2130
2131 C...Sum up Cabibbo-Kobayashi-Maskawa factors for each quark/lepton.
2132       DO 120 I=-20,20
2133         VINT(180+I)=0D0
2134         IA=IABS(I)
2135         IF(IA.GE.1.AND.IA.LE.2*MSTP(1)) THEN
2136           DO 110 J=1,MSTP(1)
2137             IB=2*J-1+MOD(IA,2)
2138             IF(IB.GE.6.AND.MSTP(9).EQ.0) GOTO 110
2139             IPM=(5-ISIGN(1,I))/2
2140             IDC=J+MDCY(IA,2)+2
2141             IF(MDME(IDC,1).EQ.1.OR.MDME(IDC,1).EQ.IPM) VINT(180+I)=
2142      &      VINT(180+I)+VCKM((IA+1)/2,(IB+1)/2)
2143   110     CONTINUE
2144         ELSEIF(IA.GE.11.AND.IA.LE.10+2*MSTP(1)) THEN
2145           VINT(180+I)=1D0
2146         ENDIF
2147   120 CONTINUE
2148
2149 C...Initialize parton distributions: PDFLIB.
2150       IF(MSTP(52).EQ.2) THEN
2151         PARM(1)='NPTYPE'
2152         VALUE(1)=1
2153         PARM(2)='NGROUP'
2154         VALUE(2)=MSTP(51)/1000
2155         PARM(3)='NSET'
2156         VALUE(3)=MOD(MSTP(51),1000)
2157         PARM(4)='TMAS'
2158         VALUE(4)=PMAS(6,1)
2159         CALL PDFSET(PARM,VALUE)
2160         MINT(93)=1000000+MSTP(51)
2161       ENDIF
2162
2163 C...Choose Lambda value to use in alpha-strong.
2164       MSTU(111)=MSTP(2)
2165       IF(MSTP(3).GE.2) THEN
2166         ALAM=0.2D0
2167         NF=4
2168         IF(MSTP(52).EQ.1.AND.MSTP(51).GE.1.AND.MSTP(51).LE.10) THEN
2169           ALAM=ALAMIN(MSTP(51))
2170           NF=NFIN(MSTP(51))
2171         ELSEIF(MSTP(52).EQ.2) THEN
2172           ALAM=QCDL4
2173           NF=4
2174         ENDIF
2175         PARP(1)=ALAM
2176         PARP(61)=ALAM
2177         PARP(72)=ALAM
2178         PARU(112)=ALAM
2179         MSTU(112)=NF
2180         IF(MSTP(3).EQ.3) PARJ(81)=ALAM
2181       ENDIF
2182
2183 C...Initialize the SUSY generation: couplings, masses,
2184 C...decay modes, branching ratios, and so on.
2185       CALL PYMSIN
2186
2187 C...Initialize widths and partial widths for resonances.
2188       CALL PYINRE
2189 C...Set Z0 mass and width for e+e- routines.
2190       PARJ(123)=PMAS(23,1)
2191       PARJ(124)=PMAS(23,2)
2192
2193 C...Identify beam and target particles and frame of process.
2194       CHFRAM=FRAME//' '
2195       CHBEAM=BEAM//' '
2196       CHTARG=TARGET//' '
2197       CALL PYINBM(CHFRAM,CHBEAM,CHTARG,WIN)
2198       IF(MINT(65).EQ.1) GOTO 170
2199
2200 C...For gamma-p or gamma-gamma allow many (3 or 6) alternatives.
2201 C...For e-gamma allow 2 alternatives.
2202       MINT(121)=1
2203       MINT(123)=MSTP(14)
2204       IF(MSTP(14).EQ.10.AND.(MSEL.EQ.1.OR.MSEL.EQ.2)) THEN
2205         IF((MINT(11).EQ.22.OR.MINT(12).EQ.22).AND.
2206      &  (IABS(MINT(11)).GE.28.OR.IABS(MINT(12)).GE.28)) MINT(121)=3
2207         IF(MINT(11).EQ.22.AND.MINT(12).EQ.22) MINT(121)=6
2208         IF((MINT(11).EQ.22.OR.MINT(12).EQ.22).AND.
2209      &  (IABS(MINT(11)).EQ.11.OR.IABS(MINT(12)).EQ.11)) MINT(121)=2
2210       ENDIF
2211
2212 C...Set up kinematics of process.
2213       CALL PYINKI(0)
2214
2215 C...Precalculate flavour selection weights
2216       CALL PYKFIN
2217
2218 C...Loop over gamma-p or gamma-gamma alternatives.
2219       DO 160 IGA=1,MINT(121)
2220         MINT(122)=IGA
2221
2222 C...Select partonic subprocesses to be included in the simulation.
2223         CALL PYINPR
2224
2225 C...Count number of subprocesses on.
2226         MINT(48)=0
2227         DO 130 ISUB=1,500
2228           IF(MINT(50).EQ.0.AND.ISUB.GE.91.AND.ISUB.LE.96.AND.
2229      &    MSUB(ISUB).EQ.1) THEN
2230             WRITE(MSTU(11),5200) ISUB,CHLH(MINT(41)),CHLH(MINT(42))
2231             STOP
2232           ELSEIF(MSUB(ISUB).EQ.1.AND.ISET(ISUB).EQ.-1) THEN
2233             WRITE(MSTU(11),5300) ISUB
2234             STOP
2235           ELSEIF(MSUB(ISUB).EQ.1.AND.ISET(ISUB).LE.-2) THEN
2236             WRITE(MSTU(11),5400) ISUB
2237             STOP
2238           ELSEIF(MSUB(ISUB).EQ.1) THEN
2239             MINT(48)=MINT(48)+1
2240           ENDIF
2241   130   CONTINUE
2242         IF(MINT(48).EQ.0) THEN
2243           WRITE(MSTU(11),5500)
2244           STOP
2245         ENDIF
2246         MINT(49)=MINT(48)-MSUB(91)-MSUB(92)-MSUB(93)-MSUB(94)
2247
2248 C...Reset variables for cross-section calculation.
2249         DO 150 I=0,500
2250           DO 140 J=1,3
2251             NGEN(I,J)=0
2252             XSEC(I,J)=0D0
2253   140     CONTINUE
2254   150   CONTINUE
2255
2256 C...Find parametrized total cross-sections.
2257         CALL PYXTOT
2258
2259 C...Maxima of differential cross-sections.
2260         IF(MSTP(121).LE.1) CALL PYMAXI
2261
2262 C...Initialize possibility of pileup events.
2263         IF(MINT(121).GT.1) MSTP(131)=0
2264         IF(MSTP(131).NE.0) CALL PYPILE(1)
2265
2266 C...Initialize multiple interactions with variable impact parameter.
2267         IF(MINT(50).EQ.1.AND.(MINT(49).NE.0.OR.MSTP(131).NE.0).AND.
2268      &  MSTP(82).GE.2) CALL PYMULT(1)
2269
2270 C...Save results for gamma-p and gamma-gamma alternatives.
2271         IF(MINT(121).GT.1) CALL PYSAVE(1,IGA)
2272   160 CONTINUE
2273
2274 C...Initialization finished.
2275   170 IF(MSTP(122).GE.1) WRITE(MSTU(11),5600)
2276
2277 C...Formats for initialization information.
2278  5100 FORMAT('1',18('*'),1X,'PYINIT: initialization of PYTHIA ',
2279      &'routines',1X,17('*'))
2280  5200 FORMAT(1X,'Error: process number ',I3,' not meaningful for ',A6,
2281      &'-',A6,' interactions.'/1X,'Execution stopped!')
2282  5300 FORMAT(1X,'Error: requested subprocess',I4,' not implemented.'/
2283      &1X,'Execution stopped!')
2284  5400 FORMAT(1X,'Error: requested subprocess',I4,' not existing.'/
2285      &1X,'Execution stopped!')
2286  5500 FORMAT(1X,'Error: no subprocess switched on.'/
2287      &1X,'Execution stopped.')
2288  5600 FORMAT(/1X,22('*'),1X,'PYINIT: initialization completed',1X,
2289      &22('*'))
2290
2291       RETURN
2292       END
2293
2294 C*********************************************************************
2295
2296 C...PYEVNT
2297 C...Administers the generation of a high-pT event via calls to
2298 C...a number of subroutines.
2299
2300       SUBROUTINE PYEVNT
2301
2302 C...Double precision and integer declarations.
2303       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
2304       INTEGER PYK,PYCHGE,PYCOMP
2305 C...Commonblocks.
2306       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
2307       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
2308       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
2309       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
2310       COMMON/PYINT1/MINT(400),VINT(400)
2311       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
2312       COMMON/PYINT4/MWID(500),WIDS(500,5)
2313       COMMON/PYINT5/NGENPD,NGEN(0:500,3),XSEC(0:500,3)
2314       COMMON/PYUPPR/NUP,KUP(20,7),NFUP,IFUP(10,2),PUP(20,5),Q2UP(0:10)
2315       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYPARS/,/PYINT1/,/PYINT2/,
2316      &/PYINT4/,/PYINT5/,/PYUPPR/
2317 C...Local array.
2318       DIMENSION VTX(4)
2319
2320 C...Initial values for some counters.
2321       N=0
2322       MINT(5)=MINT(5)+1
2323       MINT(7)=0
2324       MINT(8)=0
2325       MINT(83)=0
2326       MINT(84)=MSTP(126)
2327       MSTU(24)=0
2328       MSTU70=0
2329       MSTJ14=MSTJ(14)
2330
2331 C...If variable energies: redo incoming kinematics and cross-section.
2332       MSTI(61)=0
2333       IF(MSTP(171).EQ.1) THEN
2334         CALL PYINKI(1)
2335         IF(MSTI(61).EQ.1) THEN
2336           MINT(5)=MINT(5)-1
2337           RETURN
2338         ENDIF
2339         IF(MINT(121).GT.1) CALL PYSAVE(3,1)
2340         CALL PYXTOT
2341       ENDIF
2342
2343 C...Loop over number of pileup events; check space left.
2344       IF(MSTP(131).LE.0) THEN
2345         NPILE=1
2346       ELSE
2347         CALL PYPILE(2)
2348         NPILE=MINT(81)
2349       ENDIF
2350       DO 260 IPILE=1,NPILE
2351         IF(MINT(84)+100.GE.MSTU(4)) THEN
2352           CALL PYERRM(11,
2353      &    '(PYEVNT:) no more space in PYJETS for pileup events')
2354           IF(MSTU(21).GE.1) GOTO 270
2355         ENDIF
2356         MINT(82)=IPILE
2357
2358 C...Generate variables of hard scattering.
2359         MINT(51)=0
2360         MSTI(52)=0
2361   100   CONTINUE
2362         IF(MINT(51).NE.0.OR.MSTU(24).NE.0) MSTI(52)=MSTI(52)+1
2363         MINT(31)=0
2364         MINT(51)=0
2365         MINT(57)=0
2366         CALL PYRAND
2367         IF(MSTI(61).EQ.1) THEN
2368           MINT(5)=MINT(5)-1
2369           RETURN
2370         ENDIF
2371         IF(MINT(51).EQ.2) RETURN
2372         ISUB=MINT(1)
2373         IF(MSTP(111).EQ.-1) GOTO 250
2374
2375         IF(ISUB.LE.90.OR.ISUB.GE.95) THEN
2376 C...Hard scattering (including low-pT):
2377 C...reconstruct kinematics and colour flow of hard scattering.
2378   110     MINT(51)=0
2379           CALL PYSCAT
2380           IF(MINT(51).EQ.1) GOTO 100
2381           IPU1=MINT(84)+1
2382           IPU2=MINT(84)+2
2383           IF(ISUB.EQ.95) GOTO 130
2384
2385 C...Showering of initial state partons (optional).
2386           ALAMSV=PARJ(81)
2387           PARJ(81)=PARP(72)
2388           IF(MSTP(61).GE.1.AND.MINT(47).GE.2) CALL PYSSPA(IPU1,IPU2)
2389           PARJ(81)=ALAMSV
2390           IF(MINT(51).EQ.1) GOTO 100
2391
2392 C...Showering of final state partons (optional).
2393           ALAMSV=PARJ(81)
2394           PARJ(81)=PARP(72)
2395           IF(MSTP(71).GE.1.AND.ISET(ISUB).GE.2.AND.ISET(ISUB).LE.10)
2396      &    THEN
2397             IPU3=MINT(84)+3
2398             IPU4=MINT(84)+4
2399             IF(ISET(ISUB).EQ.5) IPU4=-3
2400             QMAX=VINT(55)
2401             IF(ISET(ISUB).EQ.2) QMAX=SQRT(PARP(71))*VINT(55)
2402             CALL PYSHOW(IPU3,IPU4,QMAX)
2403           ELSEIF(MSTP(71).GE.1.AND.ISET(ISUB).EQ.11.AND.NFUP.GE.1) THEN
2404             DO 120 IUP=1,NFUP
2405               IPU3=IFUP(IUP,1)+MINT(84)
2406               IPU4=IFUP(IUP,2)+MINT(84)
2407               QMAX=SQRT(MAX(0D0,Q2UP(IUP)))
2408               CALL PYSHOW(IPU3,IPU4,QMAX)
2409   120       CONTINUE
2410           ENDIF
2411           PARJ(81)=ALAMSV
2412
2413 C...Decay of final state resonances.
2414           MINT(32)=0
2415           IF(MSTP(41).GE.1.AND.ISET(ISUB).LE.10) CALL PYRESD(0)
2416           IF(MINT(51).EQ.1) GOTO 100
2417           MINT(52)=N
2418
2419 C...Multiple interactions.
2420           IF(MSTP(81).GE.1.AND.MINT(50).EQ.1) CALL PYMULT(6)
2421           MINT(53)=N
2422
2423 C...Hadron remnants and primordial kT.
2424   130     CALL PYREMN(IPU1,IPU2)
2425           IF(MINT(51).EQ.1.AND.MINT(57).GE.1.AND.MINT(57).LE.5) GOTO 110
2426           IF(MINT(51).EQ.1) GOTO 100
2427
2428         ELSE
2429 C...Diffractive and elastic scattering.
2430           CALL PYDIFF
2431         ENDIF
2432
2433 C...Check that no odd resonance left undecayed.
2434         IF(MSTP(111).GE.1) THEN
2435           NFIX=N
2436           DO 140 I=MINT(84)+1,NFIX
2437             IF(K(I,1).GE.1.AND.K(I,1).LE.10.AND.K(I,2).NE.21.AND.
2438      &      K(I,2).NE.22) THEN
2439               IF(MWID(PYCOMP(K(I,2))).NE.0) THEN
2440                 CALL PYRESD(I)
2441                 IF(MINT(51).EQ.1) GOTO 100
2442               ENDIF
2443             ENDIF
2444   140     CONTINUE
2445         ENDIF
2446
2447 C...Recalculate energies from momenta and masses (if desired).
2448         IF(MSTP(113).GE.1) THEN
2449           DO 150 I=MINT(83)+1,N
2450             IF(K(I,1).GT.0.AND.K(I,1).LE.10) P(I,4)=SQRT(P(I,1)**2+
2451      &      P(I,2)**2+P(I,3)**2+P(I,5)**2)
2452   150     CONTINUE
2453           NRECAL=N
2454         ENDIF
2455
2456 C...Rearrange partons along strings, check invariant mass cuts.
2457         MSTU(28)=0
2458         IF(MSTP(111).LE.0) MSTJ(14)=-1
2459         CALL PYPREP(MINT(84)+1)
2460         MSTJ(14)=MSTJ14
2461         IF(MSTP(112).EQ.1.AND.MSTU(28).EQ.3) GOTO 100
2462         IF(MSTP(125).EQ.0.OR.MSTP(125).EQ.1) THEN
2463           DO 180 I=MINT(84)+1,N
2464             IF(K(I,2).EQ.94) THEN
2465               DO 170 I1=I+1,MIN(N,I+3)
2466                 IF(K(I1,3).EQ.I) THEN
2467                   K(I1,3)=MOD(K(I1,4)/MSTU(5),MSTU(5))
2468                   IF(K(I1,3).EQ.0) THEN
2469                     DO 160 II=MINT(84)+1,I-1
2470                         IF(K(II,2).EQ.K(I1,2)) THEN
2471                           IF(MOD(K(II,4),MSTU(5)).EQ.I1.OR.
2472      &                    MOD(K(II,5),MSTU(5)).EQ.I1) K(I1,3)=II
2473                         ENDIF
2474   160               CONTINUE
2475                     IF(K(I+1,3).EQ.0) K(I+1,3)=K(I,3)
2476                   ENDIF
2477                 ENDIF
2478   170         CONTINUE
2479             ENDIF
2480   180     CONTINUE
2481           CALL PYEDIT(12)
2482           CALL PYEDIT(14)
2483           IF(MSTP(125).EQ.0) CALL PYEDIT(15)
2484           IF(MSTP(125).EQ.0) MINT(4)=0
2485           DO 200 I=MINT(83)+1,N
2486             IF(K(I,1).EQ.11.AND.K(I,4).EQ.0.AND.K(I,5).EQ.0) THEN
2487               DO 190 I1=I+1,N
2488                 IF(K(I1,3).EQ.I.AND.K(I,4).EQ.0) K(I,4)=I1
2489                 IF(K(I1,3).EQ.I) K(I,5)=I1
2490   190         CONTINUE
2491             ENDIF
2492   200     CONTINUE
2493         ENDIF
2494
2495 C...Introduce separators between sections in PYLIST event listing.
2496         IF(IPILE.EQ.1.AND.MSTP(125).LE.0) THEN
2497           MSTU70=1
2498           MSTU(71)=N
2499         ELSEIF(IPILE.EQ.1) THEN
2500           MSTU70=3
2501           MSTU(71)=2
2502           MSTU(72)=MINT(4)
2503           MSTU(73)=N
2504         ENDIF
2505
2506 C...Go back to lab frame (needed for vertices, also in fragmentation).
2507         CALL PYFRAM(1)
2508
2509 C...Set nonvanishing production vertex (optional).
2510         IF(MSTP(151).EQ.1) THEN
2511           DO 210 J=1,4
2512             VTX(J)=PARP(150+J)*SQRT(-2D0*LOG(MAX(1D-10,PYR(0))))*
2513      &      SIN(PARU(2)*PYR(0))
2514   210     CONTINUE
2515           DO 230 I=MINT(83)+1,N
2516             DO 220 J=1,4
2517               V(I,J)=V(I,J)+VTX(J)
2518   220       CONTINUE
2519   230     CONTINUE
2520         ENDIF
2521
2522 C...Perform hadronization (if desired).
2523         IF(MSTP(111).GE.1) THEN
2524           CALL PYEXEC
2525           IF(MSTU(24).NE.0) GOTO 100
2526         ENDIF
2527         IF(MSTP(113).GE.1) THEN
2528           DO 240 I=NRECAL,N
2529             IF(P(I,5).GT.0D0) P(I,4)=SQRT(P(I,1)**2+
2530      &      P(I,2)**2+P(I,3)**2+P(I,5)**2)
2531   240     CONTINUE
2532         ENDIF
2533         IF(MSTP(125).EQ.0.OR.MSTP(125).EQ.1) CALL PYEDIT(14)
2534
2535 C...Store event information and calculate Monte Carlo estimates of
2536 C...subprocess cross-sections.
2537   250   IF(IPILE.EQ.1) CALL PYDOCU
2538
2539 C...Set counters for current pileup event and loop to next one.
2540         MSTI(41)=IPILE
2541         IF(IPILE.GE.2.AND.IPILE.LE.10) MSTI(40+IPILE)=ISUB
2542         IF(MSTU70.LT.10) THEN
2543           MSTU70=MSTU70+1
2544           MSTU(70+MSTU70)=N
2545         ENDIF
2546         MINT(83)=N
2547         MINT(84)=N+MSTP(126)
2548         IF(IPILE.LT.NPILE) CALL PYFRAM(2)
2549   260 CONTINUE
2550
2551 C...Generic information on pileup events. Reconstruct missing history.
2552       IF(MSTP(131).EQ.1.AND.MSTP(133).GE.1) THEN
2553         PARI(91)=VINT(132)
2554         PARI(92)=VINT(133)
2555         PARI(93)=VINT(134)
2556         IF(MSTP(133).GE.2) PARI(93)=PARI(93)*XSEC(0,3)/VINT(131)
2557       ENDIF
2558       CALL PYEDIT(16)
2559
2560 C...Transform to the desired coordinate frame.
2561   270 CALL PYFRAM(MSTP(124))
2562       MSTU(70)=MSTU70
2563       PARU(21)=VINT(1)
2564
2565       RETURN
2566       END
2567
2568 C***********************************************************************
2569
2570 C...PYSTAT
2571 C...Prints out information about cross-sections, decay widths, branching
2572 C...ratios, kinematical limits, status codes and parameter values.
2573
2574       SUBROUTINE PYSTAT(MSTAT)
2575
2576 C...Double precision and integer declarations.
2577       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
2578       INTEGER PYK,PYCHGE,PYCOMP
2579 C...Parameter statement to help give large particle numbers.
2580       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KEXCIT=4000000)
2581 C...Commonblocks.
2582       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
2583       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
2584       COMMON/PYDAT3/MDCY(500,3),MDME(4000,2),BRAT(4000),KFDP(4000,5)
2585       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
2586       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
2587       COMMON/PYINT1/MINT(400),VINT(400)
2588       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
2589       COMMON/PYINT4/MWID(500),WIDS(500,5)
2590       COMMON/PYINT5/NGENPD,NGEN(0:500,3),XSEC(0:500,3)
2591       COMMON/PYINT6/PROC(0:500)
2592       CHARACTER PROC*28
2593       COMMON/PYMSSM/IMSS(0:99),RMSS(0:99)
2594       SAVE /PYDAT1/,/PYDAT2/,/PYDAT3/,/PYSUBS/,/PYPARS/,/PYINT1/,
2595      &/PYINT2/,/PYINT4/,/PYINT5/,/PYINT6/,/PYMSSM/
2596 C...Local arrays, character variables and data.
2597       DIMENSION WDTP(0:200),WDTE(0:200,0:5)
2598       CHARACTER PROGA(6)*28,CHAU*16,CHKF*16,CHD1*16,CHD2*16,CHD3*16,
2599      &CHIN(2)*12,STATE(-1:5)*4,CHKIN(21)*18,DISGA(2)*28
2600       DATA PROGA/
2601      &'VMD/hadron * VMD            ','VMD/hadron * direct         ',
2602      &'VMD/hadron * anomalous      ','direct * direct             ',
2603      &'direct * anomalous          ','anomalous * anomalous       '/
2604       DATA DISGA/'e * VMD','e * anomalous'/
2605       DATA STATE/'----','off ','on  ','on/+','on/-','on/1','on/2'/,
2606      &CHKIN/' m_hard (GeV/c^2) ',' p_T_hard (GeV/c) ',
2607      &'m_finite (GeV/c^2)','   y*_subsystem   ','     y*_large     ',
2608      &'     y*_small     ','    eta*_large    ','    eta*_small    ',
2609      &'cos(theta*)_large ','cos(theta*)_small ','       x_1        ',
2610      &'       x_2        ','       x_F        ',' cos(theta_hard)  ',
2611      &'m''_hard (GeV/c^2) ','       tau        ','        y*        ',
2612      &'cos(theta_hard^-) ','cos(theta_hard^+) ','      x_T^2       ',
2613      &'       tau''       '/
2614
2615 C...Cross-sections.
2616       IF(MSTAT.LE.1) THEN
2617         IF(MINT(121).GT.1) CALL PYSAVE(5,0)
2618         WRITE(MSTU(11),5000)
2619         WRITE(MSTU(11),5100)
2620         WRITE(MSTU(11),5200) 0,PROC(0),NGEN(0,3),NGEN(0,1),XSEC(0,3)
2621         DO 100 I=1,500
2622           IF(MSUB(I).NE.1) GOTO 100
2623           WRITE(MSTU(11),5200) I,PROC(I),NGEN(I,3),NGEN(I,1),XSEC(I,3)
2624   100   CONTINUE
2625         IF(MINT(121).GT.1) THEN
2626           WRITE(MSTU(11),5300)
2627           DO 110 IGA=1,MINT(121)
2628             CALL PYSAVE(3,IGA)
2629             IF(MINT(121).EQ.2) THEN
2630               WRITE(MSTU(11),5200) IGA,DISGA(IGA),NGEN(0,3),NGEN(0,1),
2631      &        XSEC(0,3)
2632             ELSE
2633               WRITE(MSTU(11),5200) IGA,PROGA(IGA),NGEN(0,3),NGEN(0,1),
2634      &        XSEC(0,3)
2635             ENDIF
2636   110     CONTINUE
2637           CALL PYSAVE(5,0)
2638         ENDIF
2639         WRITE(MSTU(11),5400) 1D0-DBLE(NGEN(0,3))/
2640      &  MAX(1D0,DBLE(NGEN(0,2)))
2641
2642 C...Decay widths and branching ratios.
2643       ELSEIF(MSTAT.EQ.2) THEN
2644         WRITE(MSTU(11),5500)
2645         WRITE(MSTU(11),5600)
2646         DO 140 KC=1,500
2647           KF=KCHG(KC,4)
2648           CALL PYNAME(KF,CHKF)
2649           IOFF=0
2650           IF(KC.LE.22) THEN
2651             IF(KC.GT.2*MSTP(1).AND.KC.LE.10) GOTO 140
2652             IF(KC.GT.10+2*MSTP(1).AND.KC.LE.20) GOTO 140
2653             IF(KC.LE.5.OR.(KC.GE.11.AND.KC.LE.16)) IOFF=1
2654             IF(KC.EQ.18.AND.PMAS(18,1).LT.1D0) IOFF=1
2655             IF(KC.EQ.21.OR.KC.EQ.22) IOFF=1
2656           ELSE
2657             IF(MWID(KC).LE.0) GOTO 140
2658             IF(IMSS(1).LE.0.AND.(KF/KSUSY1.EQ.1.OR.
2659      &      KF/KSUSY1.EQ.2)) GOTO 140
2660           ENDIF
2661 C...Off-shell branchings.
2662           IF(IOFF.EQ.1) THEN
2663             NGP=0
2664             IF(KC.LE.20) NGP=(MOD(KC,10)+1)/2
2665             IF(NGP.LE.MSTP(1)) WRITE(MSTU(11),5700) KF,CHKF(1:10),
2666      &      PMAS(KC,1),0D0,0D0,STATE(MDCY(KC,1)),0D0
2667             DO 120 J=1,MDCY(KC,3)
2668               IDC=J+MDCY(KC,2)-1
2669               NGP1=0
2670               IF(IABS(KFDP(IDC,1)).LE.20) NGP1=
2671      &        (MOD(IABS(KFDP(IDC,1)),10)+1)/2
2672               NGP2=0
2673               IF(IABS(KFDP(IDC,2)).LE.20) NGP2=
2674      &        (MOD(IABS(KFDP(IDC,2)),10)+1)/2
2675               CALL PYNAME(KFDP(IDC,1),CHD1)
2676               CALL PYNAME(KFDP(IDC,2),CHD2)
2677               IF(KFDP(IDC,3).EQ.0) THEN
2678                 IF(MDME(IDC,2).EQ.102.AND.NGP1.LE.MSTP(1).AND.
2679      &          NGP2.LE.MSTP(1)) WRITE(MSTU(11),5800) IDC,CHD1(1:10),
2680      &          CHD2(1:10),0D0,0D0,STATE(MDME(IDC,1)),0D0
2681               ELSE
2682                 CALL PYNAME(KFDP(IDC,3),CHD3)
2683                 IF(MDME(IDC,2).EQ.102.AND.NGP1.LE.MSTP(1).AND.
2684      &          NGP2.LE.MSTP(1)) WRITE(MSTU(11),5900) IDC,CHD1(1:10),
2685      &          CHD2(1:10),CHD3(1:10),0D0,0D0,STATE(MDME(IDC,1)),0D0
2686               ENDIF
2687   120       CONTINUE
2688 C...On-shell decays.
2689           ELSE
2690             CALL PYWIDT(KF,PMAS(KC,1)**2,WDTP,WDTE)
2691             BRFIN=1D0
2692             IF(WDTE(0,0).LE.0D0) BRFIN=0D0
2693             WRITE(MSTU(11),5700) KF,CHKF(1:10),PMAS(KC,1),WDTP(0),1D0,
2694      &      STATE(MDCY(KC,1)),BRFIN
2695             DO 130 J=1,MDCY(KC,3)
2696               IDC=J+MDCY(KC,2)-1
2697               NGP1=0
2698               IF(IABS(KFDP(IDC,1)).LE.20) NGP1=
2699      &        (MOD(IABS(KFDP(IDC,1)),10)+1)/2
2700               NGP2=0
2701               IF(IABS(KFDP(IDC,2)).LE.20) NGP2=
2702      &        (MOD(IABS(KFDP(IDC,2)),10)+1)/2
2703               BRFIN=0D0
2704               IF(WDTE(0,0).GT.0D0) BRFIN=WDTE(J,0)/WDTE(0,0)
2705               CALL PYNAME(KFDP(IDC,1),CHD1)
2706               CALL PYNAME(KFDP(IDC,2),CHD2)
2707               IF(KFDP(IDC,3).EQ.0) THEN
2708                 IF(NGP1.LE.MSTP(1).AND.NGP2.LE.MSTP(1))
2709      &          WRITE(MSTU(11),5800) IDC,CHD1(1:10),
2710      &          CHD2(1:10),WDTP(J),WDTP(J)/WDTP(0),
2711      &          STATE(MDME(IDC,1)),BRFIN
2712               ELSE
2713                 CALL PYNAME(KFDP(IDC,3),CHD3)
2714                 IF(NGP1.LE.MSTP(1).AND.NGP2.LE.MSTP(1))
2715      &          WRITE(MSTU(11),5900) IDC,CHD1(1:10),
2716      &          CHD2(1:10),CHD3(1:10),WDTP(J),WDTP(J)/WDTP(0),
2717      &          STATE(MDME(IDC,1)),BRFIN
2718               ENDIF
2719   130       CONTINUE
2720           ENDIF
2721   140   CONTINUE
2722         WRITE(MSTU(11),6000)
2723
2724 C...Allowed incoming partons/particles at hard interaction.
2725       ELSEIF(MSTAT.EQ.3) THEN
2726         WRITE(MSTU(11),6100)
2727         CALL PYNAME(MINT(11),CHAU)
2728         CHIN(1)=CHAU(1:12)
2729         CALL PYNAME(MINT(12),CHAU)
2730         CHIN(2)=CHAU(1:12)
2731         WRITE(MSTU(11),6200) CHIN(1),CHIN(2)
2732         DO 150 I=-20,22
2733           IF(I.EQ.0) GOTO 150
2734           IA=IABS(I)
2735           IF(IA.GT.MSTP(58).AND.IA.LE.10) GOTO 150
2736           IF(IA.GT.10+2*MSTP(1).AND.IA.LE.20) GOTO 150
2737           CALL PYNAME(I,CHAU)
2738           WRITE(MSTU(11),6300) CHAU,STATE(KFIN(1,I)),CHAU,
2739      &    STATE(KFIN(2,I))
2740   150   CONTINUE
2741         WRITE(MSTU(11),6400)
2742
2743 C...User-defined limits on kinematical variables.
2744       ELSEIF(MSTAT.EQ.4) THEN
2745         WRITE(MSTU(11),6500)
2746         WRITE(MSTU(11),6600)
2747         SHRMAX=CKIN(2)
2748         IF(SHRMAX.LT.0D0) SHRMAX=VINT(1)
2749         WRITE(MSTU(11),6700) CKIN(1),CHKIN(1),SHRMAX
2750         PTHMIN=MAX(CKIN(3),CKIN(5))
2751         PTHMAX=CKIN(4)
2752         IF(PTHMAX.LT.0D0) PTHMAX=0.5D0*SHRMAX
2753         WRITE(MSTU(11),6800) CKIN(3),PTHMIN,CHKIN(2),PTHMAX
2754         WRITE(MSTU(11),6900) CHKIN(3),CKIN(6)
2755         DO 160 I=4,14
2756           WRITE(MSTU(11),6700) CKIN(2*I-1),CHKIN(I),CKIN(2*I)
2757   160   CONTINUE
2758         SPRMAX=CKIN(32)
2759         IF(SPRMAX.LT.0D0) SPRMAX=VINT(1)
2760         WRITE(MSTU(11),6700) CKIN(31),CHKIN(15),SPRMAX
2761         WRITE(MSTU(11),7000)
2762
2763 C...Status codes and parameter values.
2764       ELSEIF(MSTAT.EQ.5) THEN
2765         WRITE(MSTU(11),7100)
2766         WRITE(MSTU(11),7200)
2767         DO 170 I=1,100
2768           WRITE(MSTU(11),7300) I,MSTP(I),PARP(I),100+I,MSTP(100+I),
2769      &    PARP(100+I)
2770   170   CONTINUE
2771
2772 C...List of all processes implemented in the program.
2773       ELSEIF(MSTAT.EQ.6) THEN
2774         WRITE(MSTU(11),7400)
2775         WRITE(MSTU(11),7500)
2776         DO 180 I=1,500
2777           IF(ISET(I).LT.0) GOTO 180
2778           WRITE(MSTU(11),7600) I,PROC(I),ISET(I),KFPR(I,1),KFPR(I,2)
2779   180   CONTINUE
2780         WRITE(MSTU(11),7700)
2781       ENDIF
2782
2783 C...Formats for printouts.
2784  5000 FORMAT('1',9('*'),1X,'PYSTAT:  Statistics on Number of ',
2785      &'Events and Cross-sections',1X,9('*'))
2786  5100 FORMAT(/1X,78('=')/1X,'I',34X,'I',28X,'I',12X,'I'/1X,'I',12X,
2787      &'Subprocess',12X,'I',6X,'Number of points',6X,'I',4X,'Sigma',3X,
2788      &'I'/1X,'I',34X,'I',28X,'I',12X,'I'/1X,'I',34('-'),'I',28('-'),
2789      &'I',4X,'(mb)',4X,'I'/1X,'I',34X,'I',28X,'I',12X,'I'/1X,'I',1X,
2790      &'N:o',1X,'Type',25X,'I',4X,'Generated',9X,'Tried',1X,'I',12X,
2791      &'I'/1X,'I',34X,'I',28X,'I',12X,'I'/1X,78('=')/1X,'I',34X,'I',28X,
2792      &'I',12X,'I')
2793  5200 FORMAT(1X,'I',1X,I3,1X,A28,1X,'I',1X,I12,1X,I13,1X,'I',1X,1P,
2794      &D10.3,1X,'I')
2795  5300 FORMAT(1X,'I',34X,'I',28X,'I',12X,'I'/1X,78('=')/
2796      &1X,'I',34X,'I',28X,'I',12X,'I')
2797  5400 FORMAT(1X,'I',34X,'I',28X,'I',12X,'I'/1X,78('=')//
2798      &1X,'********* Fraction of events that fail fragmentation ',
2799      &'cuts =',1X,F8.5,' *********'/)
2800  5500 FORMAT('1',27('*'),1X,'PYSTAT:  Decay Widths and Branching ',
2801      &'Ratios',1X,27('*'))
2802  5600 FORMAT(/1X,98('=')/1X,'I',49X,'I',13X,'I',12X,'I',6X,'I',12X,'I'/
2803      &1X,'I',5X,'Mother  -->  Branching/Decay Channel',8X,'I',1X,
2804      &'Width (GeV)',1X,'I',7X,'B.R.',1X,'I',1X,'Stat',1X,'I',2X,
2805      &'Eff. B.R.',1X,'I'/1X,'I',49X,'I',13X,'I',12X,'I',6X,'I',12X,'I'/
2806      &1X,98('='))
2807  5700 FORMAT(1X,'I',49X,'I',13X,'I',12X,'I',6X,'I',12X,'I'/1X,'I',1X,
2808      &I8,2X,A10,3X,'(m =',F10.3,')',2X,'-->',5X,'I',2X,1P,D10.3,0P,1X,
2809      &'I',1X,1P,D10.3,0P,1X,'I',1X,A4,1X,'I',1X,1P,D10.3,0P,1X,'I')
2810  5800 FORMAT(1X,'I',1X,I8,2X,A10,1X,'+',1X,A10,15X,'I',2X,
2811      &1P,D10.3,0P,1X,'I',1X,1P,D10.3,0P,1X,'I',1X,A4,1X,'I',1X,
2812      &1P,D10.3,0P,1X,'I')
2813  5900 FORMAT(1X,'I',1X,I8,2X,A10,1X,'+',1X,A10,1X,'+',1X,A10,2X,'I',2X,
2814      &1P,D10.3,0P,1X,'I',1X,1P,D10.3,0P,1X,'I',1X,A4,1X,'I',1X,
2815      &1P,D10.3,0P,1X,'I')
2816  6000 FORMAT(1X,'I',49X,'I',13X,'I',12X,'I',6X,'I',12X,'I'/1X,98('='))
2817  6100 FORMAT('1',7('*'),1X,'PYSTAT: Allowed Incoming Partons/',
2818      &'Particles at Hard Interaction',1X,7('*'))
2819  6200 FORMAT(/1X,78('=')/1X,'I',38X,'I',37X,'I'/1X,'I',1X,
2820      &'Beam particle:',1X,A12,10X,'I',1X,'Target particle:',1X,A12,7X,
2821      &'I'/1X,'I',38X,'I',37X,'I'/1X,'I',1X,'Content',6X,'State',19X,
2822      &'I',1X,'Content',6X,'State',18X,'I'/1X,'I',38X,'I',37X,'I'/1X,
2823      &78('=')/1X,'I',38X,'I',37X,'I')
2824  6300 FORMAT(1X,'I',1X,A9,5X,A4,19X,'I',1X,A9,5X,A4,18X,'I')
2825  6400 FORMAT(1X,'I',38X,'I',37X,'I'/1X,78('='))
2826  6500 FORMAT('1',12('*'),1X,'PYSTAT: User-Defined Limits on ',
2827      &'Kinematical Variables',1X,12('*'))
2828  6600 FORMAT(/1X,78('=')/1X,'I',76X,'I')
2829  6700 FORMAT(1X,'I',16X,1P,D10.3,0P,1X,'<',1X,A,1X,'<',1X,1P,D10.3,0P,
2830      &16X,'I')
2831  6800 FORMAT(1X,'I',3X,1P,D10.3,0P,1X,'(',1P,D10.3,0P,')',1X,'<',1X,A,
2832      &1X,'<',1X,1P,D10.3,0P,16X,'I')
2833  6900 FORMAT(1X,'I',29X,A,1X,'=',1X,1P,D10.3,0P,16X,'I')
2834  7000 FORMAT(1X,'I',76X,'I'/1X,78('='))
2835  7100 FORMAT('1',12('*'),1X,'PYSTAT: Summary of Status Codes and ',
2836      &'Parameter Values',1X,12('*'))
2837  7200 FORMAT(/3X,'I',4X,'MSTP(I)',9X,'PARP(I)',20X,'I',4X,'MSTP(I)',9X,
2838      &'PARP(I)'/)
2839  7300 FORMAT(1X,I3,5X,I6,6X,1P,D10.3,0P,18X,I3,5X,I6,6X,1P,D10.3)
2840  7400 FORMAT('1',13('*'),1X,'PYSTAT: List of implemented processes',
2841      &1X,13('*'))
2842  7500 FORMAT(/1X,65('=')/1X,'I',34X,'I',28X,'I'/1X,'I',12X,
2843      &'Subprocess',12X,'I',1X,'ISET',2X,'KFPR(I,1)',2X,'KFPR(I,2)',1X,
2844      &'I'/1X,'I',34X,'I',28X,'I'/1X,65('=')/1X,'I',34X,'I',28X,'I')
2845  7600 FORMAT(1X,'I',1X,I3,1X,A28,1X,'I',1X,I4,1X,I10,1X,I10,1X,'I')
2846  7700 FORMAT(1X,'I',34X,'I',28X,'I'/1X,65('='))
2847
2848       RETURN
2849       END
2850
2851 C*********************************************************************
2852
2853 C...PYINRE
2854 C...Calculates full and effective widths of gauge bosons, stores
2855 C...masses and widths, rescales coefficients to be used for
2856 C...resonance production generation.
2857
2858       SUBROUTINE PYINRE
2859
2860 C...Double precision and integer declarations.
2861       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
2862       INTEGER PYK,PYCHGE,PYCOMP
2863 C...Parameter statement to help give large particle numbers.
2864       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KEXCIT=4000000)
2865 C...Commonblocks.
2866       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
2867       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
2868       COMMON/PYDAT3/MDCY(500,3),MDME(4000,2),BRAT(4000),KFDP(4000,5)
2869       COMMON/PYDAT4/CHAF(500,2)
2870       CHARACTER CHAF*16
2871       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
2872       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
2873       COMMON/PYINT1/MINT(400),VINT(400)
2874       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
2875       COMMON/PYINT4/MWID(500),WIDS(500,5)
2876       COMMON/PYINT6/PROC(0:500)
2877       CHARACTER PROC*28
2878       COMMON/PYMSSM/IMSS(0:99),RMSS(0:99)
2879       SAVE /PYDAT1/,/PYDAT2/,/PYDAT3/,/PYDAT4/,/PYSUBS/,/PYPARS/,
2880      &/PYINT1/,/PYINT2/,/PYINT4/,/PYINT6/,/PYMSSM/
2881 C...Local arrays and data.
2882       DIMENSION WDTP(0:200),WDTE(0:200,0:5),WDTPM(0:200),
2883      &WDTEM(0:200,0:5),KCORD(500),PMORD(500)
2884
2885 C...Born level couplings in MSSM Higgs doublet sector.
2886       XW=PARU(102)
2887       XWV=XW
2888       IF(MSTP(8).GE.2) XW=1D0-(PMAS(24,1)/PMAS(23,1))**2
2889       XW1=1D0-XW
2890       IF(MSTP(4).EQ.2) THEN
2891         TANBE=PARU(141)
2892         RATBE=((1D0-TANBE**2)/(1D0+TANBE**2))**2
2893         SQMZ=PMAS(23,1)**2
2894         SQMW=PMAS(24,1)**2
2895         SQMH=PMAS(25,1)**2
2896         SQMA=SQMH*(SQMZ-SQMH)/(SQMZ*RATBE-SQMH)
2897         SQMHP=0.5D0*(SQMA+SQMZ+SQRT((SQMA+SQMZ)**2-4D0*SQMA*SQMZ*RATBE))
2898         SQMHC=SQMA+SQMW
2899         IF(SQMH.GE.SQMZ.OR.MIN(SQMA,SQMHP,SQMHC).LE.0D0) THEN
2900           WRITE(MSTU(11),5000)
2901           STOP
2902         ENDIF
2903         PMAS(35,1)=SQRT(SQMHP)
2904         PMAS(36,1)=SQRT(SQMA)
2905         PMAS(37,1)=SQRT(SQMHC)
2906         ALSU=0.5D0*ATAN(2D0*TANBE*(SQMA+SQMZ)/((1D0-TANBE**2)*
2907      &  (SQMA-SQMZ)))
2908         BESU=ATAN(TANBE)
2909         PARU(142)=1D0
2910         PARU(143)=1D0
2911         PARU(161)=-SIN(ALSU)/COS(BESU)
2912         PARU(162)=COS(ALSU)/SIN(BESU)
2913         PARU(163)=PARU(161)
2914         PARU(164)=SIN(BESU-ALSU)
2915         PARU(165)=PARU(164)
2916         PARU(168)=SIN(BESU-ALSU)+0.5D0*COS(2D0*BESU)*SIN(BESU+ALSU)/XW
2917         PARU(171)=COS(ALSU)/COS(BESU)
2918         PARU(172)=SIN(ALSU)/SIN(BESU)
2919         PARU(173)=PARU(171)
2920         PARU(174)=COS(BESU-ALSU)
2921         PARU(175)=PARU(174)
2922         PARU(176)=COS(2D0*ALSU)*COS(BESU+ALSU)-2D0*SIN(2D0*ALSU)*
2923      &  SIN(BESU+ALSU)
2924         PARU(177)=COS(2D0*BESU)*COS(BESU+ALSU)
2925         PARU(178)=COS(BESU-ALSU)-0.5D0*COS(2D0*BESU)*COS(BESU+ALSU)/XW
2926         PARU(181)=TANBE
2927         PARU(182)=1D0/TANBE
2928         PARU(183)=PARU(181)
2929         PARU(184)=0D0
2930         PARU(185)=PARU(184)
2931         PARU(186)=COS(BESU-ALSU)
2932         PARU(187)=SIN(BESU-ALSU)
2933         PARU(188)=PARU(186)
2934         PARU(189)=PARU(187)
2935         PARU(190)=0D0
2936         PARU(195)=COS(BESU-ALSU)
2937       ENDIF
2938
2939 C...Reset effective widths of gauge bosons.
2940       DO 110 I=1,500
2941         DO 100 J=1,5
2942           WIDS(I,J)=1D0
2943   100   CONTINUE
2944   110 CONTINUE
2945
2946 C...Order resonances by increasing mass (except Z0 and W+/-).
2947       NRES=0
2948       DO 140 KC=1,500
2949         KF=KCHG(KC,4)
2950         IF(KF.EQ.0) GOTO 140
2951         IF(MWID(KC).EQ.0) GOTO 140
2952         IF(KC.EQ.7.OR.KC.EQ.8.OR.KC.EQ.17.OR.KC.EQ.18) THEN
2953           IF(MSTP(1).LE.3) GOTO 140
2954         ENDIF
2955         IF(KF/KSUSY1.EQ.1.OR.KF/KSUSY1.EQ.2) THEN
2956           IF(IMSS(1).LE.0) GOTO 140
2957         ENDIF
2958         NRES=NRES+1
2959         PMRES=PMAS(KC,1)
2960         IF(KC.EQ.23.OR.KC.EQ.24) PMRES=0D0
2961         DO 120 I1=NRES-1,1,-1
2962           IF(PMRES.GE.PMORD(I1)) GOTO 130
2963           KCORD(I1+1)=KCORD(I1)
2964           PMORD(I1+1)=PMORD(I1)
2965   120   CONTINUE
2966   130   KCORD(I1+1)=KC
2967         PMORD(I1+1)=PMRES
2968   140 CONTINUE
2969
2970 C...Loop over possible resonances.
2971       DO 180 I=1,NRES
2972         KC=KCORD(I)
2973         KF=KCHG(KC,4)
2974
2975 C...Check that no fourth generation channels on by mistake.
2976         IF(MSTP(1).LE.3) THEN
2977           DO 150 J=1,MDCY(KC,3)
2978             IDC=J+MDCY(KC,2)-1
2979             KFA1=IABS(KFDP(IDC,1))
2980             KFA2=IABS(KFDP(IDC,2))
2981             IF(KFA1.EQ.7.OR.KFA1.EQ.8.OR.KFA1.EQ.17.OR.KFA1.EQ.18.OR.
2982      &      KFA2.EQ.7.OR.KFA2.EQ.8.OR.KFA2.EQ.17.OR.KFA2.EQ.18)
2983      &      MDME(IDC,1)=-1
2984   150     CONTINUE
2985         ENDIF
2986
2987 C...Check that no supersymmetric channels on by mistake.
2988         IF(IMSS(1).LE.0) THEN
2989           DO 160 J=1,MDCY(KC,3)
2990             IDC=J+MDCY(KC,2)-1
2991             KFA1S=IABS(KFDP(IDC,1))/KSUSY1
2992             KFA2S=IABS(KFDP(IDC,2))/KSUSY1
2993             IF(KFA1S.EQ.1.OR.KFA1S.EQ.2.OR.KFA2S.EQ.1.OR.KFA2S.EQ.2)
2994      &      MDME(IDC,1)=-1
2995   160     CONTINUE
2996         ENDIF
2997
2998 C...Find mass and evaluate width.
2999         PMR=PMAS(KC,1)
3000         IF(KF.EQ.25.OR.KF.EQ.35.OR.KF.EQ.36) MINT(62)=1
3001         IF(MWID(KC).EQ.3) MINT(63)=1
3002         CALL PYWIDT(KF,PMR**2,WDTP,WDTE)
3003         MINT(51)=0
3004
3005 C...Evaluate suppression factors due to non-simulated channels.
3006         IF(KCHG(KC,3).EQ.0) THEN
3007           WIDS(KC,1)=((WDTE(0,1)+WDTE(0,2))**2+
3008      &    2D0*(WDTE(0,1)+WDTE(0,2))*(WDTE(0,4)+WDTE(0,5))+
3009      &    2D0*WDTE(0,4)*WDTE(0,5))/WDTP(0)**2
3010           WIDS(KC,2)=(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))/WDTP(0)
3011           WIDS(KC,3)=0D0
3012           WIDS(KC,4)=0D0
3013           WIDS(KC,5)=0D0
3014         ELSE
3015           IF(MWID(KC).EQ.3) MINT(63)=1
3016           CALL PYWIDT(-KF,PMR**2,WDTPM,WDTEM)
3017           MINT(51)=0
3018           WIDS(KC,1)=((WDTE(0,1)+WDTE(0,2))*(WDTEM(0,1)+WDTEM(0,3))+
3019      &    (WDTE(0,1)+WDTE(0,2))*(WDTEM(0,4)+WDTEM(0,5))+
3020      &    (WDTE(0,4)+WDTE(0,5))*(WDTEM(0,1)+WDTEM(0,3))+
3021      &    WDTE(0,4)*WDTEM(0,5)+WDTE(0,5)*WDTEM(0,4))/WDTP(0)**2
3022           WIDS(KC,2)=(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))/WDTP(0)
3023           WIDS(KC,3)=(WDTEM(0,1)+WDTEM(0,3)+WDTEM(0,4))/WDTP(0)
3024           WIDS(KC,4)=((WDTE(0,1)+WDTE(0,2))**2+
3025      &    2D0*(WDTE(0,1)+WDTE(0,2))*(WDTE(0,4)+WDTE(0,5))+
3026      &    2D0*WDTE(0,4)*WDTE(0,5))/WDTP(0)**2
3027           WIDS(KC,5)=((WDTEM(0,1)+WDTEM(0,3))**2+
3028      &    2D0*(WDTEM(0,1)+WDTEM(0,3))*(WDTEM(0,4)+WDTEM(0,5))+
3029      &    2D0*WDTEM(0,4)*WDTEM(0,5))/WDTP(0)**2
3030         ENDIF
3031
3032 C...Set resonance widths and branching ratios;
3033 C...also on/off switch for decays.
3034         IF(MWID(KC).EQ.1.OR.MWID(KC).EQ.3) THEN
3035           PMAS(KC,2)=WDTP(0)
3036           PMAS(KC,3)=MIN(0.9D0*PMAS(KC,1),10D0*PMAS(KC,2))
3037           MDCY(KC,1)=MSTP(41)
3038           DO 170 J=1,MDCY(KC,3)
3039             IDC=J+MDCY(KC,2)-1
3040             BRAT(IDC)=0D0
3041             IF(WDTP(0).GT.0D0) BRAT(IDC)=WDTP(J)/WDTP(0)
3042   170     CONTINUE
3043         ENDIF
3044   180 CONTINUE
3045
3046 C...Flavours of leptoquark: redefine charge and name.
3047       KFLQQ=KFDP(MDCY(39,2),1)
3048       KFLQL=KFDP(MDCY(39,2),2)
3049       KCHG(39,1)=KCHG(PYCOMP(KFLQQ),1)*ISIGN(1,KFLQQ)+
3050      &KCHG(PYCOMP(KFLQL),1)*ISIGN(1,KFLQL)
3051       LL=1
3052       IF(IABS(KFLQL).EQ.13) LL=2
3053       IF(IABS(KFLQL).EQ.15) LL=3
3054       CHAF(39,1)='LQ_'//CHAF(IABS(KFLQQ),1)(1:1)//
3055      &CHAF(IABS(KFLQL),1)(1:LL)//' '
3056       CHAF(39,2)=CHAF(39,2)(1:4+LL)//'bar '
3057
3058 C...Special cases in treatment of gamma*/Z0: redefine process name.
3059       IF(MSTP(43).EQ.1) THEN
3060         PROC(1)='f + fbar -> gamma*'
3061         PROC(15)='f + fbar -> g + gamma*'
3062         PROC(19)='f + fbar -> gamma + gamma*'
3063         PROC(30)='f + g -> f + gamma*'
3064         PROC(35)='f + gamma -> f + gamma*'
3065       ELSEIF(MSTP(43).EQ.2) THEN
3066         PROC(1)='f + fbar -> Z0'
3067         PROC(15)='f + fbar -> g + Z0'
3068         PROC(19)='f + fbar -> gamma + Z0'
3069         PROC(30)='f + g -> f + Z0'
3070         PROC(35)='f + gamma -> f + Z0'
3071       ELSEIF(MSTP(43).EQ.3) THEN
3072         PROC(1)='f + fbar -> gamma*/Z0'
3073         PROC(15)='f + fbar -> g + gamma*/Z0'
3074         PROC(19)='f + fbar -> gamma + gamma*/Z0'
3075         PROC(30)='f + g -> f + gamma*/Z0'
3076         PROC(35)='f + gamma -> f + gamma*/Z0'
3077       ENDIF
3078
3079 C...Special cases in treatment of gamma*/Z0/Z'0: redefine process name.
3080       IF(MSTP(44).EQ.1) THEN
3081         PROC(141)='f + fbar -> gamma*'
3082       ELSEIF(MSTP(44).EQ.2) THEN
3083         PROC(141)='f + fbar -> Z0'
3084       ELSEIF(MSTP(44).EQ.3) THEN
3085         PROC(141)='f + fbar -> Z''0'
3086       ELSEIF(MSTP(44).EQ.4) THEN
3087         PROC(141)='f + fbar -> gamma*/Z0'
3088       ELSEIF(MSTP(44).EQ.5) THEN
3089         PROC(141)='f + fbar -> gamma*/Z''0'
3090       ELSEIF(MSTP(44).EQ.6) THEN
3091         PROC(141)='f + fbar -> Z0/Z''0'
3092       ELSEIF(MSTP(44).EQ.7) THEN
3093         PROC(141)='f + fbar -> gamma*/Z0/Z''0'
3094       ENDIF
3095
3096 C...Special cases in treatment of WW -> WW: redefine process name.
3097       IF(MSTP(45).EQ.1) THEN
3098         PROC(77)='W+ + W+ -> W+ + W+'
3099       ELSEIF(MSTP(45).EQ.2) THEN
3100         PROC(77)='W+ + W- -> W+ + W-'
3101       ELSEIF(MSTP(45).EQ.3) THEN
3102         PROC(77)='W+/- + W+/- -> W+/- + W+/-'
3103       ENDIF
3104
3105 C...Format for error information.
3106  5000 FORMAT(1X,'Error: unphysical input tan^2(beta) and m_H ',
3107      &'combination'/1X,'Execution stopped!')
3108
3109       RETURN
3110       END
3111
3112 C*********************************************************************
3113
3114 C...PYINBM
3115 C...Identifies the two incoming particles and the choice of frame.
3116
3117        SUBROUTINE PYINBM(CHFRAM,CHBEAM,CHTARG,WIN)
3118
3119 C...Double precision and integer declarations.
3120       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
3121       INTEGER PYK,PYCHGE,PYCOMP
3122 C...Commonblocks.
3123       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
3124       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
3125       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
3126       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
3127       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
3128       COMMON/PYINT1/MINT(400),VINT(400)
3129       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYSUBS/,/PYPARS/,/PYINT1/
3130 C...Local arrays, character variables and data.
3131       CHARACTER CHFRAM*8,CHBEAM*8,CHTARG*8,CHCOM(3)*8,CHALP(2)*26,
3132      &CHIDNT(3)*8,CHTEMP*8,CHCDE(29)*8,CHINIT*76
3133       DIMENSION LEN(3),KCDE(29),PM(2)
3134       DATA CHALP/'abcdefghijklmnopqrstuvwxyz',
3135      &'ABCDEFGHIJKLMNOPQRSTUVWXYZ'/
3136       DATA CHCDE/'e-      ','e+      ','nu_e    ','nu_ebar ',
3137      &'mu-     ','mu+     ','nu_mu   ','nu_mubar','tau-    ',
3138      &'tau+    ','nu_tau  ','nu_tauba','pi+     ','pi-     ',
3139      &'n0      ','nbar0   ','p+      ','pbar-   ','gamma   ',
3140      &'lambda0 ','sigma-  ','sigma0  ','sigma+  ','xi-     ',
3141      &'xi0     ','omega-  ','pi0     ','reggeon ','pomeron '/
3142       DATA KCDE/11,-11,12,-12,13,-13,14,-14,15,-15,16,-16,
3143      &211,-211,2112,-2112,2212,-2212,22,3122,3112,3212,3222,
3144      &3312,3322,3334,111,28,29/
3145
3146 C...Store initial energy. Default frame.
3147       VINT(290)=WIN
3148       MINT(111)=0
3149
3150 C...Convert character variables to lowercase and find their length.
3151       CHCOM(1)=CHFRAM
3152       CHCOM(2)=CHBEAM
3153       CHCOM(3)=CHTARG
3154       DO 130 I=1,3
3155         LEN(I)=8
3156         DO 110 LL=8,1,-1
3157           IF(LEN(I).EQ.LL.AND.CHCOM(I)(LL:LL).EQ.' ') LEN(I)=LL-1
3158           DO 100 LA=1,26
3159             IF(CHCOM(I)(LL:LL).EQ.CHALP(2)(LA:LA)) CHCOM(I)(LL:LL)=
3160      &      CHALP(1)(LA:LA)
3161   100     CONTINUE
3162   110   CONTINUE
3163         CHIDNT(I)=CHCOM(I)
3164
3165 C...Fix up bar, underscore and charge in particle name (if needed).
3166         DO 120 LL=1,6
3167           IF(CHIDNT(I)(LL:LL).EQ.'~') THEN
3168             CHTEMP=CHIDNT(I)
3169             CHIDNT(I)=CHTEMP(1:LL-1)//'bar'//CHTEMP(LL+1:6)//'  '
3170           ENDIF
3171   120   CONTINUE
3172         IF(CHIDNT(I)(7:7).EQ.'~') CHIDNT(I)(7:8)='ba'
3173         IF(CHIDNT(I)(1:2).EQ.'nu'.AND.CHIDNT(I)(3:3).NE.'_') THEN
3174           CHTEMP=CHIDNT(I)
3175           CHIDNT(I)='nu_'//CHTEMP(3:7)
3176         ELSEIF(CHIDNT(I)(1:2).EQ.'n ') THEN
3177           CHIDNT(I)(1:3)='n0 '
3178         ELSEIF(CHIDNT(I)(1:4).EQ.'nbar') THEN
3179           CHIDNT(I)(1:5)='nbar0'
3180         ELSEIF(CHIDNT(I)(1:2).EQ.'p ') THEN
3181           CHIDNT(I)(1:3)='p+ '
3182         ELSEIF(CHIDNT(I)(1:4).EQ.'pbar'.OR.
3183      &    CHIDNT(I)(1:2).EQ.'p-') THEN
3184           CHIDNT(I)(1:5)='pbar-'
3185         ELSEIF(CHIDNT(I)(1:6).EQ.'lambda') THEN
3186           CHIDNT(I)(7:7)='0'
3187         ELSEIF(CHIDNT(I)(1:3).EQ.'reg') THEN