Added sort method.
[u/mrichter/AliRoot.git] / DPMJET / pythia6115.f
CommitLineData
d30b8254 1C*********************************************************************
2C*********************************************************************
3C* **
4C* March 1997 **
5C* **
6C* The Lund Monte Carlo for Hadronic Processes **
7C* **
8C* PYTHIA version 6.1 **
9C* **
10C* Torbjorn Sjostrand **
11C* Department of Theoretical Physics 2 **
12C* Lund University **
13C* Solvegatan 14A, S-223 62 Lund, Sweden **
14C* phone +46 - 46 - 222 48 16 **
15C* E-mail torbjorn@thep.lu.se **
16C* **
17C* SUSY parts by **
18C* Stephen Mrenna **
19C* Argonne National Laboratory **
20C* 9700 South Cass Avenue, Argonne, IL 60439, USA **
21C* phone + 1 - 630 - 252 - 7615 **
22C* E-mail mrenna@hep.anl.gov **
23C* **
24C* Several parts are written by Hans-Uno Bengtsson **
25C* PYSHOW is written together with Mats Bengtsson **
26C* advanced popcorn baryon production written by Patrik Eden **
27C* CTEQ 3 parton distributions are by the CTEQ collaboration **
28C* GRV 94 parton distributions are by Glueck, Reya and Vogt **
29C* SaS photon parton distributions together with Gerhard Schuler **
30C* g + g and q + qbar -> t + tbar + H code by Zoltan Kunszt **
31C* MSSM Higgs mass calculation code by M. Carena, **
32C* J.R. Espinosa, M. Quiros and C.E.M. Wagner **
33C* PYGAUS adapted from CERN library (K.S. Kolbig) **
34C* **
35C* The latest program version and documentation is found on WWW **
36C* http://www.thep.lu.se/tf2/staff/torbjorn/Pythia.html **
37C* **
38C* Copyright Torbjorn Sjostrand, Lund 1997 **
39C* **
40C*********************************************************************
41C*********************************************************************
42C *
43C List of subprograms in order of appearance, with main purpose *
44C (S = subroutine, F = function, B = block data) *
45C *
46C B PYDATA to contain all default values *
47C S PYTEST to test the proper functioning of the package *
48C S PYHEPC to convert between /PYJETS/ and /HEPEVT/ records *
49C *
50C S PYINIT to administer the initialization procedure *
51C S PYEVNT to administer the generation of an event *
52C S PYSTAT to print cross-section and other information *
53C S PYINRE to initialize treatment of resonances *
54C S PYINBM to read in beam, target and frame choices *
55C S PYINKI to initialize kinematics of incoming particles *
56C S PYINPR to set up the selection of included processes *
57C S PYXTOT to give total, elastic and diffractive cross-sect. *
58C S PYMAXI to find differential cross-section maxima *
59C S PYPILE to select multiplicity of pileup events *
60C S PYSAVE to save alternatives for gamma-p and gamma-gamma *
61C S PYRAND to select subprocess and kinematics for event *
62C S PYSCAT to set up kinematics and colour flow of event *
63C S PYSSPA to simulate initial state spacelike showers *
64C S PYRESD to perform resonance decays *
65C S PYMULT to generate multiple interactions *
66C S PYREMN to add on target remnants *
67C S PYDIFF to set up kinematics for diffractive events *
68C S PYDOCU to compute cross-sections and handle documentation *
69C S PYFRAM to perform boosts between different frames *
70C S PYWIDT to calculate full and partial widths of resonances *
71C S PYOFSH to calculate partial width into off-shell channels *
72C S PYRECO to handle colour reconnection in W+W- events *
73C S PYKLIM to calculate borders of allowed kinematical region *
74C S PYKMAP to construct value of kinematical variable *
75C S PYSIGH to calculate differential cross-sections *
76C S PYPDFU to evaluate parton distributions *
77C S PYPDFL to evaluate parton distributions at low x and Q^2 *
78C S PYPDEL to evaluate electron parton distributions *
79C S PYPDGA to evaluate photon parton distributions (generic) *
80C S PYGGAM to evaluate photon parton distributions (SaS sets) *
81C S PYGVMD to evaluate VMD part of photon parton distributions *
82C S PYGANO to evaluate anomalous part of photon pdf's *
83C S PYGBEH to evaluate Bethe-Heitler part of photon pdf's *
84C S PYGDIR to evaluate direct contribution to photon pdf's *
85C S PYPDPI to evaluate pion parton distributions *
86C S PYPDPR to evaluate proton parton distributions *
87C F PYCTEQ to evaluate the CTEQ 3 proton parton distributions *
88C S PYGRVL to evaluate the GRV 94L pronton parton distributions *
89C S PYGRVM to evaluate the GRV 94M pronton parton distributions *
90C S PYGRVD to evaluate the GRV 94D pronton parton distributions *
91C F PYGRVV auxiliary to the PYGRV* routines *
92C F PYGRVW auxiliary to the PYGRV* routines *
93C F PYGRVS auxiliary to the PYGRV* routines *
94C F PYHFTH to evaluate threshold factor for heavy flavour *
95C S PYSPLI to find flavours left in hadron when one removed *
96C F PYGAMM to evaluate ordinary Gamma function Gamma(x) *
97C S PYWAUX to evaluate auxiliary functions W1(s) and W2(s) *
98C S PYI3AU to evaluate auxiliary function I3(s,t,u,v) *
99C F PYSPEN to evaluate Spence (dilogarithm) function Sp(x) *
100C S PYQQBH to evaluate matrix element for g + g -> Q + Qbar + H *
101C *
102C S PYMSIN to initialize the supersymmetry simulation *
103C S PYAPPS to determine MSSM parameters from SUGRA input *
104C F PYRNMQ to determine running quark masses *
105C F PYRNMT to determine running top mass *
106C S PYTHRG to calculate sfermion third-gen. mass eigenstates *
107C S PYINOM to calculate neutralino/chargino mass eigenstates *
108C F PYRNM3 to determine running M3, gluino mass *
109C S PYEIG4 to calculate eigenvalues and -vectors in 4*4 matrix *
110C S PYHGGM to determine Higgs mass spectrum *
111C S PYSUBH to determine Higgs masses in the MSSM *
112C S PYPOLE to determine Higgs masses in the MSSM *
113C S PYVACU to determine Higgs masses in the MSSM *
114C S PYRGHM auxiliary to PYVACU *
115C S PYGFXX auxiliary to PYRGHM *
116C F PYFINT auxiliary to PYVACU *
117C F PYFISB auxiliary to PYFINT *
118C S PYSFDC to calculate sfermion decay partial widths *
119C S PYGLUI to calculate gluino decay partial widths *
120C S PYTBBN to calculate 3-body decay of gluino to neutralino *
121C S PYTBBC to calculate 3-body decay of gluino to chargino *
122C S PYNJDC to calculate neutralino decay partial widths *
123C S PYCJDC to calculate chargino decay partial widths *
124C F PYXXZ5 auxiliary for neutralino 3-body decay *
125C F PYXXW5 auxiliary for ino charge change 3-body decay *
126C F PYXXGA auxiliary for ino -> ino + gamma decay *
127C F PYX2XG auxiliary for ino -> ino + gauge boson decay *
128C F PYX2XH auxiliary for ino -> ino + Higgs decay *
129C F PYXXZ2 auxiliary for chargino 3-body decay *
130C S PYHEXT to calculate non-SM Higgs decay partial widths *
131C F PYH2XX auxiliary for H -> ino + ino decay *
132C F PYGAUS to perform Gaussian integration *
133C F PYSIMP to perform Simpson integration *
134C F PYLAMF to evaluate the lambda kinematics function *
135C S PYTBDY to perform 3-body decay of gauginos *
136C *
137C S PY1ENT to fill one entry (= parton or particle) *
138C S PY2ENT to fill two entries *
139C S PY3ENT to fill three entries *
140C S PY4ENT to fill four entries *
141C S PYJOIN to connect entries with colour flow information *
142C S PYGIVE to fill (or query) commonblock variables *
143C S PYEXEC to administrate fragmentation and decay chain *
144C S PYPREP to rearrange showered partons along strings *
145C S PYSTRF to do string fragmentation of jet system *
146C S PYINDF to do independent fragmentation of one or many jets *
147C S PYDECY to do the decay of a particle *
148C S PYDCYK to select parton and hadron flavours in decays *
149C S PYKFDI to select parton and hadron flavours in fragm *
150C S PYNMES to select number of popcorn mesons *
151C S PYKFIN to calculate falvour prod. ratios from input params. *
152C S PYPTDI to select transverse momenta in fragm *
153C S PYZDIS to select longitudinal scaling variable in fragm *
154C S PYSHOW to do timelike parton shower evolution *
155C S PYBOEI to include Bose-Einstein effects (crudely) *
156C F PYMASS to give the mass of a particle or parton *
157C S PYNAME to give the name of a particle or parton *
158C F PYCHGE to give three times the electric charge *
159C F PYCOMP to compress standard KF flavour code to internal KC *
160C S PYERRM to write error messages and abort faulty run *
161C F PYALEM to give the alpha_electromagnetic value *
162C F PYALPS to give the alpha_strong value *
163C F PYANGL to give the angle from known x and y components *
164C F PYR to provide a random number generator *
165C S PYRGET to save the state of the random number generator *
166C S PYRSET to set the state of the random number generator *
167C S PYROBO to rotate and/or boost an event *
168C S PYEDIT to remove unwanted entries from record *
169C S PYLIST to list event record or particle data *
170C S PYLOGO to write a logo *
171C S PYUPDA to update particle data *
172C F PYK to provide integer-valued event information *
173C F PYP to provide real-valued event information *
174C S PYSPHE to perform sphericity analysis *
175C S PYTHRU to perform thrust analysis *
176C S PYCLUS to perform three-dimensional cluster analysis *
177C S PYCELL to perform cluster analysis in (eta, phi, E_T) *
178C S PYJMAS to give high and low jet mass of event *
179C S PYFOWO to give Fox-Wolfram moments *
180C S PYTABU to analyze events, with tabular output *
181C *
182C S PYEEVT to administrate the generation of an e+e- event *
183C S PYXTEE to give the total cross-section at given CM energy *
184C S PYRADK to generate initial state photon radiation *
185C S PYXKFL to select flavour of primary qqbar pair *
186C S PYXJET to select (matrix element) jet multiplicity *
187C S PYX3JT to select kinematics of three-jet event *
188C S PYX4JT to select kinematics of four-jet event *
189C S PYXDIF to select angular orientation of event *
190C S PYONIA to perform generation of onium decay to gluons *
191C *
192C S PYBOOK to book a histogram *
193C S PYFILL to fill an entry in a histogram *
194C S PYFACT to multiply histogram contents by a factor *
195C S PYOPER to perform operations between histograms *
196C S PYHIST to print and reset all histograms *
197C S PYPLOT to print a single histogram *
198C S PYNULL to reset contents of a single histogram *
199C S PYDUMP to dump histogram contents onto a file *
200C *
201C S PYKCUT dummy routine for user kinematical cuts *
202C S PYEVWT dummy routine for weighting events *
203C S PYUPIN dummy routine to initialize a user process *
204C S PYUPEV dummy routine to generate a user process event *
205C S PDFSET dummy routine to be removed when using PDFLIB *
206C S STRUCTM dummy routine to be removed when using PDFLIB *
207C S PYTAUD dummy routine for interface to tau decay libraries *
208C S PYTIME dummy routine for giving date and time *
209C *
210C*********************************************************************
211
212C...PYDATA
213C...Default values for switches and parameters,
214C...and particle, decay and process data.
215
216 BLOCK DATA PYDATA
217
218C...Double precision and integer declarations.
219 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
220 INTEGER PYK,PYCHGE,PYCOMP
221C...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
246C...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
312C...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
448C...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
1018C...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
1107C...PYDATR, with initial values for the random number generator.
1108 DATA MRPY/19780503,0,0,97,33,0/
1109
1110C...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
1131C...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
1180C...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
1307C...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
1311C...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
1474C...Cross sections and slope offsets.
1475 DATA SIGT/294*0D0/
1476
1477C...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
1487C...Data for histogramming routines.
1488 DATA IHIST/1000,20000,55,1/
1489 DATA INDX/1000*0/
1490
1491 END
1492
1493C*********************************************************************
1494
1495C...PYTEST
1496C...A simple program (disguised as subroutine) to run at installation
1497C...as a check that the program works as intended.
1498
1499 SUBROUTINE PYTEST(MTEST)
1500
1501C...Double precision and integer declarations.
1502 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1503 INTEGER PYK,PYCHGE,PYCOMP
1504C...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/
1512C...Local arrays.
1513 DIMENSION PSUM(5),PINI(6),PFIN(6)
1514
1515C...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
1532C...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
1537C...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
1554C...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
1565C...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
1578C...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
1590C...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
1614C...Fifty e+e- continuum events with matrix elements.
1615 ELSEIF(IEV.LE.350) THEN
1616 MSTJ(101)=2
1617 CALL PYEEVT(0,40D0)
1618
1619C...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
1626C...Fifty e+e- continuum events with coherent shower.
1627 ELSEIF(IEV.LE.450) THEN
1628 CALL PYEEVT(0,500D0)
1629
1630C...Fifty Upsilon decays to ggg or gammagg with coherent shower.
1631 ELSE
1632 CALL PYONIA(5,9.46D0)
1633 ENDIF
1634
1635C...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
1646C...Check conservation of energy, momentum and charge;
1647C...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
1664C...Check that all KF codes are known ones, and that partons/particles
1665C...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
1681C...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
1689C...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
1698C...Summarize result of run.
1699 IF(MTEST.GE.1) CALL PYTABU(22)
1700
1701C...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
1718C...Second part: complete events of various kinds.
1719C...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
1724C...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
1741C...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
1750C...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
1759C...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
1771C...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
1786C...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
1801C...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
1813C...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
1823C...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
1840C...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
1846C...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
1865C...Check that all KF codes are known ones, and that partons/particles
1866C...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
1882C...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
1895C...List statistics for each process type.
1896 IF(MTEST.GE.1) CALL PYSTAT(1)
1897 230 CONTINUE
1898
1899C...Summarize result of run.
1900 IF(NERR.EQ.0) WRITE(MSTU(11),6500)
1901 IF(NERR.GT.0) WRITE(MSTU(11),6600) NERR
1902
1903C...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
1921C*********************************************************************
1922
1923C...PYHEPC
1924C...Converts PYTHIA event record contents to or from
1925C...the standard event record commonblock.
1926
1927 SUBROUTINE PYHEPC(MCONV)
1928
1929C...Double precision and integer declarations.
1930 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1931 INTEGER PYK,PYCHGE,PYCOMP
1932C...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/
1937C...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
1944C...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
1973C...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
1980C...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
2006C...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
2028C...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
2057C...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
2077C*********************************************************************
2078
2079C...PYINIT
2080C...Initializes the generation procedure; finds maxima of the
2081C...differential cross-sections to be used for weighting.
2082
2083 SUBROUTINE PYINIT(FRAME,BEAM,TARGET,WIN)
2084
2085C...Double precision and integer declarations.
2086 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
2087 INTEGER PYK,PYCHGE,PYCOMP
2088C...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/
2101C...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
2106C...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
2113C...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
2118C...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
2126C...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
2131C...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
2149C...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
2163C...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
2183C...Initialize the SUSY generation: couplings, masses,
2184C...decay modes, branching ratios, and so on.
2185 CALL PYMSIN
2186
2187C...Initialize widths and partial widths for resonances.
2188 CALL PYINRE
2189C...Set Z0 mass and width for e+e- routines.
2190 PARJ(123)=PMAS(23,1)
2191 PARJ(124)=PMAS(23,2)
2192
2193C...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
2200C...For gamma-p or gamma-gamma allow many (3 or 6) alternatives.
2201C...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
2212C...Set up kinematics of process.
2213 CALL PYINKI(0)
2214
2215C...Precalculate flavour selection weights
2216 CALL PYKFIN
2217
2218C...Loop over gamma-p or gamma-gamma alternatives.
2219 DO 160 IGA=1,MINT(121)
2220 MINT(122)=IGA
2221
2222C...Select partonic subprocesses to be included in the simulation.
2223 CALL PYINPR
2224
2225C...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
2248C...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
2256C...Find parametrized total cross-sections.
2257 CALL PYXTOT
2258
2259C...Maxima of differential cross-sections.
2260 IF(MSTP(121).LE.1) CALL PYMAXI
2261
2262C...Initialize possibility of pileup events.
2263 IF(MINT(121).GT.1) MSTP(131)=0
2264 IF(MSTP(131).NE.0) CALL PYPILE(1)
2265
2266C...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
2270C...Save results for gamma-p and gamma-gamma alternatives.
2271 IF(MINT(121).GT.1) CALL PYSAVE(1,IGA)
2272 160 CONTINUE
2273
2274C...Initialization finished.
2275 170 IF(MSTP(122).GE.1) WRITE(MSTU(11),5600)
2276
2277C...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
2294C*********************************************************************
2295
2296C...PYEVNT
2297C...Administers the generation of a high-pT event via calls to
2298C...a number of subroutines.
2299
2300 SUBROUTINE PYEVNT
2301
2302C...Double precision and integer declarations.
2303 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
2304 INTEGER PYK,PYCHGE,PYCOMP
2305C...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/
2317C...Local array.
2318 DIMENSION VTX(4)
2319
2320C...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
2331C...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
2343C...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
2358C...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
2376C...Hard scattering (including low-pT):
2377C...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
2385C...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
2392C...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
2413C...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
2419C...Multiple interactions.
2420 IF(MSTP(81).GE.1.AND.MINT(50).EQ.1) CALL PYMULT(6)
2421 MINT(53)=N
2422
2423C...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
2429C...Diffractive and elastic scattering.
2430 CALL PYDIFF
2431 ENDIF
2432
2433C...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
2447C...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
2456C...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
2495C...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
2506C...Go back to lab frame (needed for vertices, also in fragmentation).
2507 CALL PYFRAM(1)
2508
2509C...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
2522C...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
2535C...Store event information and calculate Monte Carlo estimates of
2536C...subprocess cross-sections.
2537 250 IF(IPILE.EQ.1) CALL PYDOCU
2538
2539C...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
2551C...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
2560C...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
2568C***********************************************************************
2569
2570C...PYSTAT
2571C...Prints out information about cross-sections, decay widths, branching
2572C...ratios, kinematical limits, status codes and parameter values.
2573
2574 SUBROUTINE PYSTAT(MSTAT)
2575
2576C...Double precision and integer declarations.
2577 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
2578 INTEGER PYK,PYCHGE,PYCOMP
2579C...Parameter statement to help give large particle numbers.
2580 PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KEXCIT=4000000)
2581C...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/
2596C...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
2615C...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
2642C...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
2661C...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
2688C...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
2724C...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
2743C...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
2763C...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
2772C...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
2783C...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
2851C*********************************************************************
2852
2853C...PYINRE
2854C...Calculates full and effective widths of gauge bosons, stores
2855C...masses and widths, rescales coefficients to be used for
2856C...resonance production generation.
2857
2858 SUBROUTINE PYINRE
2859
2860C...Double precision and integer declarations.
2861 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
2862 INTEGER PYK,PYCHGE,PYCOMP
2863C...Parameter statement to help give large particle numbers.
2864 PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KEXCIT=4000000)
2865C...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/
2881C...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
2885C...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
2939C...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
2946C...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
2970C...Loop over possible resonances.
2971 DO 180 I=1,NRES
2972 KC=KCORD(I)
2973 KF=KCHG(KC,4)
2974
2975C...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
2987C...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
2998C...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
3005C...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
3032C...Set resonance widths and branching ratios;
3033C...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
3046C...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
3058C...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
3079C...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
3096C...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
3105C...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
3112C*********************************************************************
3113
3114C...PYINBM
3115C...Identifies the two incoming particles and the choice of frame.
3116
3117 SUBROUTINE PYINBM(CHFRAM,CHBEAM,CHTARG,WIN)
3118
3119C...Double precision and integer declarations.
3120 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
3121 INTEGER PYK,PYCHGE,PYCOMP
3122C...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/
3130C...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
3146C...Store initial energy. Default frame.
3147 VINT(290)=WIN
3148 MINT(111)=0
3149
3150C...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
3165C...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
3188 CHIDNT(I)(1:7)='reggeon'
3189 ELSEIF(CHIDNT(I)(1:3).EQ.'pom') THEN
3190 CHIDNT(I)(1:7)='pomeron'
3191 ENDIF
3192 130 CONTINUE
3193
3194C...Identify free initialization.
3195 IF(CHCOM(1)(1:2).EQ.'no') THEN
3196 MINT(65)=1
3197 RETURN
3198 ENDIF
3199
3200C...Identify incoming beam and target particles.
3201 DO 150 I=1,2
3202 DO 140 J=1,29
3203 IF(CHIDNT(I+1).EQ.CHCDE(J)) MINT(10+I)=KCDE(J)
3204 140 CONTINUE
3205 PM(I)=PYMASS(MINT(10+I))
3206 VINT(2+I)=PM(I)
3207 150 CONTINUE
3208 IF(MINT(11).EQ.0) WRITE(MSTU(11),5000) CHBEAM(1:LEN(2))
3209 IF(MINT(12).EQ.0) WRITE(MSTU(11),5100) CHTARG(1:LEN(3))
3210 IF(MINT(11).EQ.0.OR.MINT(12).EQ.0) STOP
3211
3212C...Identify choice of frame and input energies.
3213 CHINIT=' '
3214
3215C...Events defined in the CM frame.
3216 IF(CHCOM(1)(1:2).EQ.'cm') THEN
3217 MINT(111)=1
3218 S=WIN**2
3219 IF(MSTP(122).GE.1) THEN
3220 IF(CHCOM(2)(1:1).NE.'e') THEN
3221 LOFFS=(31-(LEN(2)+LEN(3)))/2
3222 CHINIT(LOFFS+1:76)='PYTHIA will be initialized for a '//
3223 & CHCOM(2)(1:LEN(2))//' on '//CHCOM(3)(1:LEN(3))//
3224 & ' collider'//' '
3225 ELSE
3226 LOFFS=(30-(LEN(2)+LEN(3)))/2
3227 CHINIT(LOFFS+1:76)='PYTHIA will be initialized for an '//
3228 & CHCOM(2)(1:LEN(2))//' on '//CHCOM(3)(1:LEN(3))//
3229 & ' collider'//' '
3230 ENDIF
3231 WRITE(MSTU(11),5200) CHINIT
3232 WRITE(MSTU(11),5300) WIN
3233 ENDIF
3234
3235C...Events defined in fixed target frame.
3236 ELSEIF(CHCOM(1)(1:3).EQ.'fix') THEN
3237 MINT(111)=2
3238 S=PM(1)**2+PM(2)**2+2D0*PM(2)*SQRT(PM(1)**2+WIN**2)
3239 IF(MSTP(122).GE.1) THEN
3240 LOFFS=(29-(LEN(2)+LEN(3)))/2
3241 CHINIT(LOFFS+1:76)='PYTHIA will be initialized for '//
3242 & CHCOM(2)(1:LEN(2))//' on '//CHCOM(3)(1:LEN(3))//
3243 & ' fixed target'//' '
3244 WRITE(MSTU(11),5200) CHINIT
3245 WRITE(MSTU(11),5400) WIN
3246 WRITE(MSTU(11),5500) SQRT(S)
3247 ENDIF
3248
3249C...Frame defined by user three-vectors.
3250 ELSEIF(CHCOM(1)(1:3).EQ.'use') THEN
3251 MINT(111)=3
3252 P(1,5)=PM(1)
3253 P(2,5)=PM(2)
3254 P(1,4)=SQRT(P(1,1)**2+P(1,2)**2+P(1,3)**2+P(1,5)**2)
3255 P(2,4)=SQRT(P(2,1)**2+P(2,2)**2+P(2,3)**2+P(2,5)**2)
3256 S=(P(1,4)+P(2,4))**2-(P(1,1)+P(2,1))**2-(P(1,2)+P(2,2))**2-
3257 & (P(1,3)+P(2,3))**2
3258 IF(MSTP(122).GE.1) THEN
3259 LOFFS=(12-(LEN(2)+LEN(3)))/2
3260 CHINIT(LOFFS+1:76)='PYTHIA will be initialized for '//
3261 & CHCOM(2)(1:LEN(2))//' on '//CHCOM(3)(1:LEN(3))//
3262 & ' user-specified configuration'//' '
3263 WRITE(MSTU(11),5200) CHINIT
3264 WRITE(MSTU(11),5600)
3265 WRITE(MSTU(11),5700) CHCOM(2),P(1,1),P(1,2),P(1,3),P(1,4)
3266 WRITE(MSTU(11),5700) CHCOM(3),P(2,1),P(2,2),P(2,3),P(2,4)
3267 WRITE(MSTU(11),5500) SQRT(MAX(0D0,S))
3268 ENDIF
3269
3270C...Frame defined by user four-vectors.
3271 ELSEIF(CHCOM(1)(1:4).EQ.'four') THEN
3272 MINT(111)=4
3273 PMS1=P(1,4)**2-P(1,1)**2-P(1,2)**2-P(1,3)**2
3274 P(1,5)=SIGN(SQRT(ABS(PMS1)),PMS1)
3275 PMS2=P(2,4)**2-P(2,1)**2-P(2,2)**2-P(2,3)**2
3276 P(2,5)=SIGN(SQRT(ABS(PMS2)),PMS2)
3277 S=(P(1,4)+P(2,4))**2-(P(1,1)+P(2,1))**2-(P(1,2)+P(2,2))**2-
3278 & (P(1,3)+P(2,3))**2
3279 IF(MSTP(122).GE.1) THEN
3280 LOFFS=(12-(LEN(2)+LEN(3)))/2
3281 CHINIT(LOFFS+1:76)='PYTHIA will be initialized for '//
3282 & CHCOM(2)(1:LEN(2))//' on '//CHCOM(3)(1:LEN(3))//
3283 & ' user-specified configuration'//' '
3284 WRITE(MSTU(11),5200) CHINIT
3285 WRITE(MSTU(11),5600)
3286 WRITE(MSTU(11),5700) CHCOM(2),P(1,1),P(1,2),P(1,3),P(1,4)
3287 WRITE(MSTU(11),5700) CHCOM(3),P(2,1),P(2,2),P(2,3),P(2,4)
3288 WRITE(MSTU(11),5500) SQRT(MAX(0D0,S))
3289 ENDIF
3290
3291C...Frame defined by user five-vectors.
3292 ELSEIF(CHCOM(1)(1:4).EQ.'five') THEN
3293 MINT(111)=5
3294 S=(P(1,4)+P(2,4))**2-(P(1,1)+P(2,1))**2-(P(1,2)+P(2,2))**2-
3295 & (P(1,3)+P(2,3))**2
3296 IF(MSTP(122).GE.1) THEN
3297 LOFFS=(12-(LEN(2)+LEN(3)))/2
3298 CHINIT(LOFFS+1:76)='PYTHIA will be initialized for '//
3299 & CHCOM(2)(1:LEN(2))//' on '//CHCOM(3)(1:LEN(3))//
3300 & ' user-specified configuration'//' '
3301 WRITE(MSTU(11),5200) CHINIT
3302 WRITE(MSTU(11),5600)
3303 WRITE(MSTU(11),5700) CHCOM(2),P(1,1),P(1,2),P(1,3),P(1,4)
3304 WRITE(MSTU(11),5700) CHCOM(3),P(2,1),P(2,2),P(2,3),P(2,4)
3305 WRITE(MSTU(11),5500) SQRT(MAX(0D0,S))
3306 ENDIF
3307
3308C...Unknown frame. Error for too low CM energy.
3309 ELSE
3310 WRITE(MSTU(11),5800) CHFRAM(1:LEN(1))
3311 STOP
3312 ENDIF
3313 IF(S.LT.PARP(2)**2) THEN
3314 WRITE(MSTU(11),5900) SQRT(S)
3315 STOP
3316 ENDIF
3317
3318C...Formats for initialization and error information.
3319 5000 FORMAT(1X,'Error: unrecognized beam particle ''',A,'''D0'/
3320 &1X,'Execution stopped!')
3321 5100 FORMAT(1X,'Error: unrecognized target particle ''',A,'''D0'/
3322 &1X,'Execution stopped!')
3323 5200 FORMAT(/1X,78('=')/1X,'I',76X,'I'/1X,'I',A76,'I')
3324 5300 FORMAT(1X,'I',18X,'at',1X,F10.3,1X,'GeV center-of-mass energy',
3325 &19X,'I'/1X,'I',76X,'I'/1X,78('='))
3326 5400 FORMAT(1X,'I',22X,'at',1X,F10.3,1X,'GeV/c lab-momentum',22X,'I')
3327 5500 FORMAT(1X,'I',76X,'I'/1X,'I',11X,'corresponding to',1X,F10.3,1X,
3328 &'GeV center-of-mass energy',12X,'I'/1X,'I',76X,'I'/1X,78('='))
3329 5600 FORMAT(1X,'I',76X,'I'/1X,'I',18X,'px (GeV/c)',3X,'py (GeV/c)',3X,
3330 &'pz (GeV/c)',6X,'E (GeV)',9X,'I')
3331 5700 FORMAT(1X,'I',8X,A8,4(2X,F10.3,1X),8X,'I')
3332 5800 FORMAT(1X,'Error: unrecognized coordinate frame ''',A,'''D0'/
3333 &1X,'Execution stopped!')
3334 5900 FORMAT(1X,'Error: too low CM energy,',F8.3,' GeV for event ',
3335 &'generation.'/1X,'Execution stopped!')
3336
3337 RETURN
3338 END
3339
3340C*********************************************************************
3341
3342C...PYINKI
3343C...Sets up kinematics, including rotations and boosts to/from CM frame.
3344
3345 SUBROUTINE PYINKI(MODKI)
3346
3347C...Double precision and integer declarations.
3348 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
3349 INTEGER PYK,PYCHGE,PYCOMP
3350C...Commonblocks.
3351 COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
3352 COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
3353 COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
3354 COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
3355 COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
3356 COMMON/PYINT1/MINT(400),VINT(400)
3357 SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYSUBS/,/PYPARS/,/PYINT1/
3358
3359C...Set initial flavour state.
3360 N=2
3361 DO 100 I=1,2
3362 K(I,1)=1
3363 K(I,2)=MINT(10+I)
3364 100 CONTINUE
3365
3366C...Reset boost. Do kinematics for various cases.
3367 DO 110 J=6,10
3368 VINT(J)=0D0
3369 110 CONTINUE
3370
3371C...Set up kinematics for events defined in CM frame.
3372 IF(MINT(111).EQ.1) THEN
3373 WIN=VINT(290)
3374 IF(MODKI.EQ.1) WIN=PARP(171)*VINT(290)
3375 S=WIN**2
3376 P(1,5)=VINT(3)
3377 P(2,5)=VINT(4)
3378 P(1,1)=0D0
3379 P(1,2)=0D0
3380 P(2,1)=0D0
3381 P(2,2)=0D0
3382 P(1,3)=SQRT(((S-P(1,5)**2-P(2,5)**2)**2-(2D0*P(1,5)*P(2,5))**2)/
3383 & (4D0*S))
3384 P(2,3)=-P(1,3)
3385 P(1,4)=SQRT(P(1,3)**2+P(1,5)**2)
3386 P(2,4)=SQRT(P(2,3)**2+P(2,5)**2)
3387
3388C...Set up kinematics for fixed target events.
3389 ELSEIF(MINT(111).EQ.2) THEN
3390 WIN=VINT(290)
3391 IF(MODKI.EQ.1) WIN=PARP(171)*VINT(290)
3392 P(1,5)=VINT(3)
3393 P(2,5)=VINT(4)
3394 P(1,1)=0D0
3395 P(1,2)=0D0
3396 P(2,1)=0D0
3397 P(2,2)=0D0
3398 P(1,3)=WIN
3399 P(1,4)=SQRT(P(1,3)**2+P(1,5)**2)
3400 P(2,3)=0D0
3401 P(2,4)=P(2,5)
3402 S=P(1,5)**2+P(2,5)**2+2D0*P(2,4)*P(1,4)
3403 VINT(10)=P(1,3)/(P(1,4)+P(2,4))
3404 CALL PYROBO(0,0,0D0,0D0,0D0,0D0,-VINT(10))
3405
3406C...Set up kinematics for events in user-defined frame.
3407 ELSEIF(MINT(111).EQ.3) THEN
3408 P(1,5)=VINT(3)
3409 P(2,5)=VINT(4)
3410 P(1,4)=SQRT(P(1,1)**2+P(1,2)**2+P(1,3)**2+P(1,5)**2)
3411 P(2,4)=SQRT(P(2,1)**2+P(2,2)**2+P(2,3)**2+P(2,5)**2)
3412 DO 120 J=1,3
3413 VINT(7+J)=(P(1,J)+P(2,J))/(P(1,4)+P(2,4))
3414 120 CONTINUE
3415 CALL PYROBO(0,0,0D0,0D0,-VINT(8),-VINT(9),-VINT(10))
3416 VINT(7)=PYANGL(P(1,1),P(1,2))
3417 CALL PYROBO(0,0,0D0,-VINT(7),0D0,0D0,0D0)
3418 VINT(6)=PYANGL(P(1,3),P(1,1))
3419 CALL PYROBO(0,0,-VINT(6),0D0,0D0,0D0,0D0)
3420 S=P(1,5)**2+P(2,5)**2+2D0*(P(1,4)*P(2,4)-P(1,3)*P(2,3))
3421
3422C...Set up kinematics for events with user-defined four-vectors.
3423 ELSEIF(MINT(111).EQ.4) THEN
3424 PMS1=P(1,4)**2-P(1,1)**2-P(1,2)**2-P(1,3)**2
3425 P(1,5)=SIGN(SQRT(ABS(PMS1)),PMS1)
3426 PMS2=P(2,4)**2-P(2,1)**2-P(2,2)**2-P(2,3)**2
3427 P(2,5)=SIGN(SQRT(ABS(PMS2)),PMS2)
3428 DO 130 J=1,3
3429 VINT(7+J)=(P(1,J)+P(2,J))/(P(1,4)+P(2,4))
3430 130 CONTINUE
3431 CALL PYROBO(0,0,0D0,0D0,-VINT(8),-VINT(9),-VINT(10))
3432 VINT(7)=PYANGL(P(1,1),P(1,2))
3433 CALL PYROBO(0,0,0D0,-VINT(7),0D0,0D0,0D0)
3434 VINT(6)=PYANGL(P(1,3),P(1,1))
3435 CALL PYROBO(0,0,-VINT(6),0D0,0D0,0D0,0D0)
3436 S=(P(1,4)+P(2,4))**2
3437
3438C...Set up kinematics for events with user-defined five-vectors.
3439 ELSEIF(MINT(111).EQ.5) THEN
3440 DO 140 J=1,3
3441 VINT(7+J)=(P(1,J)+P(2,J))/(P(1,4)+P(2,4))
3442 140 CONTINUE
3443 CALL PYROBO(0,0,0D0,0D0,-VINT(8),-VINT(9),-VINT(10))
3444 VINT(7)=PYANGL(P(1,1),P(1,2))
3445 CALL PYROBO(0,0,0D0,-VINT(7),0D0,0D0,0D0)
3446 VINT(6)=PYANGL(P(1,3),P(1,1))
3447 CALL PYROBO(0,0,-VINT(6),0D0,0D0,0D0,0D0)
3448 S=(P(1,4)+P(2,4))**2
3449 ENDIF
3450
3451C...Return or error for too low CM energy.
3452 IF(MODKI.EQ.1.AND.S.LT.PARP(2)**2) THEN
3453 IF(MSTP(172).LE.1) THEN
3454 CALL PYERRM(23,
3455 & '(PYINKI:) too low invariant mass in this event')
3456 ELSE
3457 MSTI(61)=1
3458 RETURN
3459 ENDIF
3460 ENDIF
3461
3462C...Save information on incoming particles.
3463 VINT(1)=SQRT(S)
3464 VINT(2)=S
3465 IF(MINT(111).GE.4) VINT(3)=P(1,5)
3466 IF(MINT(111).GE.4) VINT(4)=P(2,5)
3467 VINT(5)=P(1,3)
3468 IF(MODKI.EQ.0) VINT(289)=S
3469 DO 150 J=1,5
3470 V(1,J)=0D0
3471 V(2,J)=0D0
3472 VINT(290+J)=P(1,J)
3473 VINT(295+J)=P(2,J)
3474 150 CONTINUE
3475
3476C...Store pT cut-off and related constants to be used in generation.
3477 IF(MODKI.EQ.0) VINT(285)=CKIN(3)
3478 IF(MSTP(82).LE.1) THEN
3479 IF(MINT(121).GT.1) PARP(81)=1.30D0+0.15D0*LOG(VINT(1)/200D0)/
3480 & LOG(900D0/200D0)
3481 PTMN=PARP(81)
3482 ELSE
3483 IF(MINT(121).GT.1) PARP(82)=1.25D0+0.15D0*LOG(VINT(1)/200D0)/
3484 & LOG(900D0/200D0)
3485 PTMN=PARP(82)
3486 ENDIF
3487 VINT(149)=4D0*PTMN**2/S
3488
3489 RETURN
3490 END
3491
3492C*********************************************************************
3493
3494C...PYINPR
3495C...Selects partonic subprocesses to be included in the simulation.
3496
3497 SUBROUTINE PYINPR
3498
3499C...Double precision and integer declarations.
3500 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
3501 INTEGER PYK,PYCHGE,PYCOMP
3502C...Commonblocks.
3503 COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
3504 COMMON/PYDAT3/MDCY(500,3),MDME(4000,2),BRAT(4000),KFDP(4000,5)
3505 COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
3506 COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
3507 COMMON/PYINT1/MINT(400),VINT(400)
3508 COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
3509 SAVE /PYDAT1/,/PYDAT3/,/PYSUBS/,/PYPARS/,/PYINT1/,/PYINT2/
3510
3511C...Reset processes to be included.
3512 IF(MSEL.NE.0) THEN
3513 DO 100 I=1,500
3514 MSUB(I)=0
3515 100 CONTINUE
3516 ENDIF
3517
3518C...For e-gamma witn MSTP(14)=10 allow mixture of VMD and anomalous.
3519 IF(MINT(121).EQ.2) THEN
3520 MSUB(10)=1
3521 MINT(123)=MINT(122)+1
3522
3523C...For gamma-p or gamma-gamma with MSTP(14)=10 allow mixture.
3524C...Here also set a few parameters otherwise normally not touched.
3525 ELSEIF(MINT(121).GT.1) THEN
3526
3527C...Parton distributions dampened at small Q2; go to low energies,
3528C...alpha_s <1; no minimum pT cut-off a priori.
3529 MSTP(57)=3
3530 MSTP(85)=0
3531 PARP(2)=2D0
3532 PARU(115)=1D0
3533 CKIN(5)=0.2D0
3534 CKIN(6)=0.2D0
3535
3536C...Define pT cut-off parameters and whether run involves low-pT.
3537 IF(MSTP(82).LE.1) THEN
3538 PTMVMD=1.30D0+0.15D0*LOG(VINT(1)/200D0)/LOG(900D0/200D0)
3539 ELSE
3540 PTMVMD=1.25D0+0.15D0*LOG(VINT(1)/200D0)/LOG(900D0/200D0)
3541 ENDIF
3542 PTMDIR=PARP(15)
3543 PTMANO=PTMVMD
3544 IF(MSTP(15).EQ.5) PTMANO=0.60D0+
3545 & 0.125D0*LOG(1D0+0.10D0*VINT(1))**2
3546 IPTL=1
3547 IF(VINT(285).GT.MAX(PTMVMD,PTMDIR,PTMANO)) IPTL=0
3548 IF(MSEL.EQ.2) IPTL=1
3549
3550C...Set up for p/VMD * VMD.
3551 IF(MINT(122).EQ.1) THEN
3552 MINT(123)=2
3553 MSUB(11)=1
3554 MSUB(12)=1
3555 MSUB(13)=1
3556 MSUB(28)=1
3557 MSUB(53)=1
3558 MSUB(68)=1
3559 IF(IPTL.EQ.1) MSUB(95)=1
3560 IF(MSEL.EQ.2) THEN
3561 MSUB(91)=1
3562 MSUB(92)=1
3563 MSUB(93)=1
3564 MSUB(94)=1
3565 ENDIF
3566 PARP(81)=PTMVMD
3567 PARP(82)=PTMVMD
3568 IF(IPTL.EQ.1) CKIN(3)=0D0
3569
3570C...Set up for p/VMD * direct gamma.
3571 ELSEIF(MINT(122).EQ.2) THEN
3572 MINT(123)=0
3573 IF(MINT(121).EQ.6) MINT(123)=5
3574 MSUB(33)=1
3575 MSUB(54)=1
3576 IF(IPTL.EQ.1) CKIN(3)=PTMDIR
3577
3578C...Set up for p/VMD * anomalous gamma.
3579 ELSEIF(MINT(122).EQ.3) THEN
3580 MINT(123)=3
3581 IF(MINT(121).EQ.6) MINT(123)=7
3582 MSUB(11)=1
3583 MSUB(12)=1
3584 MSUB(13)=1
3585 MSUB(28)=1
3586 MSUB(53)=1
3587 MSUB(68)=1
3588 IF(MSTP(82).GE.2) MSTP(85)=1
3589 IF(IPTL.EQ.1) CKIN(3)=PTMANO
3590
3591C...Set up for direct * direct gamma (switch off leptons).
3592 ELSEIF(MINT(122).EQ.4) THEN
3593 MINT(123)=0
3594 MSUB(58)=1
3595 DO 110 II=MDCY(22,2),MDCY(22,2)+MDCY(22,3)-1
3596 IF(IABS(KFDP(II,1)).GE.10) MDME(II,1)=MIN(0,MDME(II,1))
3597 110 CONTINUE
3598 IF(IPTL.EQ.1) CKIN(3)=PTMDIR
3599
3600C...Set up for direct * anomalous gamma.
3601 ELSEIF(MINT(122).EQ.5) THEN
3602 MINT(123)=6
3603 MSUB(33)=1
3604 MSUB(54)=1
3605 IF(IPTL.EQ.1) CKIN(3)=PTMANO
3606
3607C...Set up for anomalous * anomalous gamma.
3608 ELSEIF(MINT(122).EQ.6) THEN
3609 MINT(123)=3
3610 MSUB(11)=1
3611 MSUB(12)=1
3612 MSUB(13)=1
3613 MSUB(28)=1
3614 MSUB(53)=1
3615 MSUB(68)=1
3616 IF(MSTP(82).GE.2) MSTP(85)=1
3617 IF(IPTL.EQ.1) CKIN(3)=PTMANO
3618 ENDIF
3619
3620C...End of special set up for gamma-p and gamma-gamma.
3621 CKIN(1)=2D0*CKIN(3)
3622 ENDIF
3623
3624C...Flavour information for individual beams.
3625 DO 120 I=1,2
3626 MINT(40+I)=1
3627 IF(MINT(123).GE.1.AND.MINT(10+I).EQ.22) MINT(40+I)=2
3628 IF(IABS(MINT(10+I)).GT.100) MINT(40+I)=2
3629 IF(MINT(10+I).EQ.28.OR.MINT(10+I).EQ.29) MINT(40+I)=2
3630 MINT(44+I)=MINT(40+I)
3631 IF(MSTP(11).GE.1.AND.IABS(MINT(10+I)).EQ.11) MINT(44+I)=3
3632 120 CONTINUE
3633
3634C...If two gammas, whereof one direct, pick the first.
3635 IF(MINT(11).EQ.22.AND.MINT(12).EQ.22) THEN
3636 IF(MINT(123).GE.4.AND.MINT(123).LE.6) THEN
3637 MINT(41)=1
3638 MINT(45)=1
3639 ENDIF
3640 ELSEIF(MINT(11).EQ.22.OR.MINT(12).EQ.22) THEN
3641 IF(MINT(123).GE.4) CALL PYERRM(26,
3642 & '(PYINPR:) unallowed MSTP(14) code for single photon')
3643 ENDIF
3644
3645C...Flavour information on combination of incoming particles.
3646 MINT(43)=2*MINT(41)+MINT(42)-2
3647 MINT(44)=MINT(43)
3648 IF(MINT(123).LE.0) THEN
3649 IF(MINT(11).EQ.22) MINT(43)=MINT(43)+2
3650 IF(MINT(12).EQ.22) MINT(43)=MINT(43)+1
3651 ELSEIF(MINT(123).LE.3) THEN
3652 IF(MINT(11).EQ.22) MINT(44)=MINT(44)-2
3653 IF(MINT(12).EQ.22) MINT(44)=MINT(44)-1
3654 ELSEIF(MINT(11).EQ.22.AND.MINT(12).EQ.22) THEN
3655 MINT(43)=4
3656 MINT(44)=1
3657 ENDIF
3658 MINT(47)=2*MIN(2,MINT(45))+MIN(2,MINT(46))-2
3659 IF(MIN(MINT(45),MINT(46)).EQ.3) MINT(47)=5
3660 MINT(50)=0
3661 IF(MINT(41).EQ.2.AND.MINT(42).EQ.2) MINT(50)=1
3662 IF((MINT(11).EQ.22.OR.MINT(12).EQ.22).AND.MINT(123).GE.3)
3663 &MINT(50)=0
3664 MINT(107)=0
3665 IF(MINT(11).EQ.22) THEN
3666 MINT(107)=MINT(123)
3667 IF(MINT(123).GE.4) MINT(107)=0
3668 IF(MINT(123).EQ.7) MINT(107)=2
3669 ENDIF
3670 MINT(108)=0
3671 IF(MINT(12).EQ.22) THEN
3672 MINT(108)=MINT(123)
3673 IF(MINT(123).GE.4) MINT(108)=MINT(123)-3
3674 IF(MINT(123).EQ.7) MINT(108)=3
3675 ENDIF
3676
3677C...Select default processes according to incoming beams
3678C...(already done for gamma-p and gamma-gamma with MSTP(14)=10).
3679 IF(MINT(121).GT.1) THEN
3680 ELSEIF(MSEL.EQ.1.OR.MSEL.EQ.2) THEN
3681
3682 IF(MINT(43).EQ.1) THEN
3683C...Lepton + lepton -> gamma/Z0 or W.
3684 IF(MINT(11)+MINT(12).EQ.0) MSUB(1)=1
3685 IF(MINT(11)+MINT(12).NE.0) MSUB(2)=1
3686
3687 ELSEIF(MINT(43).LE.3.AND.MINT(123).EQ.0.AND.
3688 & (MINT(11).EQ.22.OR.MINT(12).EQ.22)) THEN
3689C...Unresolved photon + lepton: Compton scattering.
3690 MSUB(34)=1
3691
3692 ELSEIF(MINT(43).LE.3) THEN
3693C...Lepton + hadron: deep inelastic scattering.
3694 MSUB(10)=1