Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyscat.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyscat.f
1 
2 C*********************************************************************
3 
4 C...PYSCAT
5 C...Finds outgoing flavours and event type; sets up the kinematics
6 C...and colour flow of the hard scattering
7 
8  SUBROUTINE pyscat
9 
10 C...Double precision and integer declarations
11  IMPLICIT DOUBLE PRECISION(a-h, o-z)
12  IMPLICIT INTEGER(i-n)
13  INTEGER pyk,pychge,pycomp
14 C...Parameter statement to help give large particle numbers.
15  parameter(ksusy1=1000000,ksusy2=2000000,ktechn=3000000,
16  &kexcit=4000000,kdimen=5000000)
17 C...Parameter statement for maximum size of showers.
18  parameter(maxnur=1000)
19 
20 C...User process event common block.
21  INTEGER maxnup
22  parameter(maxnup=500)
23  INTEGER nup,idprup,idup,istup,mothup,icolup
24  DOUBLE PRECISION xwgtup,scalup,aqedup,aqcdup,pup,vtimup,spinup
25  common/hepeup/nup,idprup,xwgtup,scalup,aqedup,aqcdup,idup(maxnup),
26  &istup(maxnup),mothup(2,maxnup),icolup(2,maxnup),pup(5,maxnup),
27  &vtimup(maxnup),spinup(maxnup)
28  SAVE /hepeup/
29 
30 C...Commonblocks.
31  common/pypart/npart,npartd,ipart(maxnur),ptpart(maxnur)
32  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
33  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
34  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
35  common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
36  common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
37  common/pypars/mstp(200),parp(200),msti(200),pari(200)
38  common/pyint1/mint(400),vint(400)
39  common/pyint2/iset(500),kfpr(500,2),coef(500,20),icol(40,4,2)
40  common/pyint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)
41  common/pyint4/mwid(500),wids(500,5)
42  common/pyint5/ngenpd,ngen(0:500,3),xsec(0:500,3)
43  common/pyssmt/zmix(4,4),umix(2,2),vmix(2,2),smz(4),smw(2),
44  &sfmix(16,4),zmixi(4,4),umixi(2,2),vmixi(2,2)
45  common/pytcsm/itcm(0:99),rtcm(0:99)
46  SAVE /pypart/,/pyjets/,/pydat1/,/pydat2/,/pydat3/,/pysubs/,
47  &/pypars/,/pyint1/,/pyint2/,/pyint3/,/pyint4/,/pyint5/,/pyssmt/,
48  &/pytcsm/
49 C...Local arrays and saved variables
50  dimension wdtp(0:400),wdte(0:400,0:5),pmq(2),z(2),cthe(2),
51  &phi(2),kuppo(100),vintsv(41:66),ilab(100)
52  SAVE vintsv
53 
54 C...Read out process
55  isub=mint(1)
56  isubsv=isub
57 
58 C...Restore information for low-pT processes
59  IF(isub.EQ.95.AND.mint(57).GE.1) THEN
60  DO 100 j=41,66
61  100 vint(j)=vintsv(j)
62  ENDIF
63 
64 C...Convert H' or A process into equivalent H one
65  ihigg=1
66  kfhigg=25
67  IF((isub.GE.151.AND.isub.LE.160).OR.(isub.GE.171.AND.
68  &isub.LE.190)) THEN
69  ihigg=2
70  IF(mod(isub-1,10).GE.5) ihigg=3
71  kfhigg=33+ihigg
72  IF(isub.EQ.151.OR.isub.EQ.156) isub=3
73  IF(isub.EQ.152.OR.isub.EQ.157) isub=102
74  IF(isub.EQ.153.OR.isub.EQ.158) isub=103
75  IF(isub.EQ.171.OR.isub.EQ.176) isub=24
76  IF(isub.EQ.172.OR.isub.EQ.177) isub=26
77  IF(isub.EQ.173.OR.isub.EQ.178) isub=123
78  IF(isub.EQ.174.OR.isub.EQ.179) isub=124
79  IF(isub.EQ.181.OR.isub.EQ.186) isub=121
80  IF(isub.EQ.182.OR.isub.EQ.187) isub=122
81  IF(isub.EQ.183.OR.isub.EQ.188) isub=111
82  IF(isub.EQ.184.OR.isub.EQ.189) isub=112
83  IF(isub.EQ.185.OR.isub.EQ.190) isub=113
84  ENDIF
85 
86  IF(isub.EQ.401.OR.isub.EQ.402) kfhigg=kfpr(isub,1)
87 
88 C...Convert bottomonium process into equivalent charmonium ones.
89  IF(isub.GE.461.AND.isub.LE.479) isub=isub-40
90 
91 C...Choice of subprocess, number of documentation lines
92  idoc=6+iset(isub)
93  IF(isub.EQ.95) idoc=8
94  IF(iset(isub).EQ.5) idoc=9
95  IF(iset(isub).EQ.11) idoc=4+nup
96  mint(3)=idoc-6
97  IF(idoc.GE.9.AND.iset(isub).LE.4) idoc=idoc+2
98  mint(4)=idoc
99  ipu1=mint(84)+1
100  ipu2=mint(84)+2
101  ipu3=mint(84)+3
102  ipu4=mint(84)+4
103  ipu5=mint(84)+5
104  ipu6=mint(84)+6
105 
106 C...Reset K, P and V vectors. Store incoming particles
107  DO 120 jt=1,mstp(126)+100
108  i=mint(83)+jt
109  IF(i.GT.mstu(4)) goto 120
110  DO 110 j=1,5
111  k(i,j)=0
112  p(i,j)=0d0
113  v(i,j)=0d0
114  110 CONTINUE
115  120 CONTINUE
116  DO 140 jt=1,2
117  i=mint(83)+jt
118  k(i,1)=21
119  k(i,2)=mint(10+jt)
120  DO 130 j=1,5
121  p(i,j)=vint(285+5*jt+j)
122  130 CONTINUE
123  140 CONTINUE
124  mint(6)=2
125  kfres=0
126 
127 C...Store incoming partons in their CM-frame. Save pdf value.
128  sh=vint(44)
129  shr=sqrt(sh)
130  shp=vint(26)*vint(2)
131  shpr=sqrt(shp)
132  shuser=shr
133  IF(iset(isub).GE.3.AND.iset(isub).LE.5) shuser=shpr
134  DO 150 jt=1,2
135  i=mint(84)+jt
136  k(i,1)=14
137  k(i,2)=mint(14+jt)
138  k(i,3)=mint(83)+2+jt
139  p(i,3)=0.5d0*shuser*(-1d0)**(jt-1)
140  p(i,4)=0.5d0*shuser
141  vint(38+jt)=xsfx(jt,mint(14+jt))
142  150 CONTINUE
143 
144 C...Copy incoming partons to documentation lines
145  DO 170 jt=1,2
146  i1=mint(83)+4+jt
147  i2=mint(84)+jt
148  k(i1,1)=21
149  k(i1,2)=k(i2,2)
150  k(i1,3)=i1-2
151  DO 160 j=1,5
152  p(i1,j)=p(i2,j)
153  160 CONTINUE
154  170 CONTINUE
155 
156 C...Choose new quark/lepton flavour for relevant annihilation graphs
157  IF(isub.EQ.12.OR.isub.EQ.53.OR.isub.EQ.54.OR.isub.EQ.58.OR.
158  &(isub.GE.135.AND.isub.LE.140).OR.isub.EQ.382.OR.isub.EQ.385) THEN
159  iglga=21
160  IF(isub.EQ.58.OR.(isub.GE.137.AND.isub.LE.140)) iglga=22
161  CALL pywidt(iglga,sh,wdtp,wdte)
162  180 rkfl=(wdte(0,1)+wdte(0,2)+wdte(0,4))*pyr(0)
163  DO 190 i=1,mdcy(iglga,3)
164  kflf=kfdp(i+mdcy(iglga,2)-1,1)
165  rkfl=rkfl-(wdte(i,1)+wdte(i,2)+wdte(i,4))
166  IF(rkfl.LE.0d0) goto 200
167  190 CONTINUE
168  200 CONTINUE
169  IF((isub.EQ.53.OR.isub.EQ.385).AND.mint(2).LE.2) THEN
170  IF(kflf.GE.4) goto 180
171  ELSEIF((isub.EQ.53.OR.isub.EQ.385).AND.mint(2).LE.4) THEN
172  kflf=4
173  mint(2)=mint(2)-2
174  ELSEIF(isub.EQ.53.OR.isub.EQ.385) THEN
175  kflf=5
176  mint(2)=mint(2)-4
177  ELSEIF(isub.EQ.382.AND.itcm(5).EQ.1.AND.iabs(mint(15)).LE.2
178  & .AND.iabs(kflf).GE.3) THEN
179  facqqb=vint(58)**2*4d0/9d0*(vint(45)**2+vint(46)**2)/
180  & vint(44)**2
181  faccib=vint(46)**2/rtcm(41)**4
182  IF(facqqb/(facqqb+faccib).LT.pyr(0)) goto 180
183  ELSEIF(isub.EQ.382.AND.itcm(5).EQ.5.AND.mint(2).EQ.2) THEN
184  kflf=5
185  mint(2)=1
186  ELSEIF(isub.EQ.382.AND.itcm(5).EQ.5.AND.mint(2).EQ.1) THEN
187  IF(kflf.EQ.5) goto 180
188  ELSEIF(isub.EQ.54.OR.isub.EQ.135.OR.isub.EQ.136) THEN
189  IF((kchg(pycomp(kflf),1)/2d0)**2.LT.pyr(0)) goto 180
190  ELSEIF(isub.EQ.58.OR.(isub.GE.137.AND.isub.LE.140)) THEN
191  IF((kchg(pycomp(kflf),1)/3d0)**2.LT.pyr(0)) goto 180
192  ENDIF
193  ENDIF
194 
195 C...Final state flavours and colour flow: default values
196  js=1
197  mint(21)=mint(15)
198  mint(22)=mint(16)
199  mint(23)=0
200  mint(24)=0
201  kcc=20
202  kcs=isign(1,mint(15))
203 
204  IF(iset(isub).EQ.11) THEN
205 C...User-defined processes: find products
206  mint(3)=0
207  DO 210 iup=3,nup
208  IF(istup(iup).LT.1.OR.istup(iup).GT.3) THEN
209  ELSEIF(nup.EQ.5.AND.iup.GE.4.AND.mothup(1,4).EQ.3) THEN
210  mint(21+iup)=idup(iup)
211  ELSEIF(istup(iup).EQ.1.AND.(istup(mothup(1,iup)).EQ.2.OR.
212  & istup(mothup(1,iup)).EQ.3).AND.idup(mothup(1,iup)).NE.0) THEN
213  ELSEIF(idup(iup).EQ.0) THEN
214  ELSE
215  mint(3)=mint(3)+1
216  IF(mint(3).LE.6) mint(20+mint(3))=idup(iup)
217  ENDIF
218  210 CONTINUE
219 
220  ELSEIF(isub.LE.10) THEN
221  IF(isub.EQ.1) THEN
222 C...f + fbar -> gamma*/Z0
223  kfres=23
224 
225  ELSEIF(isub.EQ.2) THEN
226 C...f + fbar' -> W+/-
227  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
228  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
229  kfres=isign(24,kch1+kch2)
230 
231  ELSEIF(isub.EQ.3) THEN
232 C...f + fbar -> h0 (or H0, or A0)
233  kfres=kfhigg
234 
235  ELSEIF(isub.EQ.4) THEN
236 C...gamma + W+/- -> W+/-
237 
238  ELSEIF(isub.EQ.5) THEN
239 C...Z0 + Z0 -> h0
240  xh=sh/shp
241  mint(21)=mint(15)
242  mint(22)=mint(16)
243  pmq(1)=pymass(mint(21))
244  pmq(2)=pymass(mint(22))
245  220 jt=int(1.5d0+pyr(0))
246  zmin=2d0*pmq(jt)/shpr
247  zmax=1d0-pmq(3-jt)/shpr-(sh-pmq(jt)**2)/
248  & (shpr*(shpr-pmq(3-jt)))
249  zmax=min(1d0-xh,zmax)
250  z(jt)=zmin+(zmax-zmin)*pyr(0)
251  IF(-1d0+(1d0+xh)/(1d0-z(jt))-xh/(1d0-z(jt))**2.LT.
252  & (1d0-xh)**2/(4d0*xh)*pyr(0)) goto 220
253  sqc1=1d0-4d0*pmq(jt)**2/(z(jt)**2*shp)
254  IF(sqc1.LT.1d-8) goto 220
255  c1=sqrt(sqc1)
256  c2=1d0+2d0*(pmas(23,1)**2-pmq(jt)**2)/(z(jt)*shp)
257  cthe(jt)=(c2-(c2**2-c1**2)/(c2+(2d0*pyr(0)-1d0)*c1))/c1
258  cthe(jt)=min(1d0,max(-1d0,cthe(jt)))
259  z(3-jt)=1d0-xh/(1d0-z(jt))
260  sqc1=1d0-4d0*pmq(3-jt)**2/(z(3-jt)**2*shp)
261  IF(sqc1.LT.1d-8) goto 220
262  c1=sqrt(sqc1)
263  c2=1d0+2d0*(pmas(23,1)**2-pmq(3-jt)**2)/(z(3-jt)*shp)
264  cthe(3-jt)=(c2-(c2**2-c1**2)/(c2+(2d0*pyr(0)-1d0)*c1))/c1
265  cthe(3-jt)=min(1d0,max(-1d0,cthe(3-jt)))
266  phir=paru(2)*pyr(0)
267  cphi=cos(phir)
268  ang=cthe(1)*cthe(2)-sqrt(1d0-cthe(1)**2)*
269  & sqrt(1d0-cthe(2)**2)*cphi
270  z1=2d0-z(jt)
271  z2=ang*sqrt(z(jt)**2-4d0*pmq(jt)**2/shp)
272  z3=1d0-z(jt)-xh+(pmq(1)**2+pmq(2)**2)/shp
273  z(3-jt)=2d0/(z1**2-z2**2)*(z1*z3+z2*sqrt(z3**2-(z1**2-z2**2)*
274  & pmq(3-jt)**2/shp))
275  zmin=2d0*pmq(3-jt)/shpr
276  zmax=1d0-pmq(jt)/shpr-(sh-pmq(3-jt)**2)/(shpr*(shpr-pmq(jt)))
277  zmax=min(1d0-xh,zmax)
278  IF(z(3-jt).LT.zmin.OR.z(3-jt).GT.zmax) goto 220
279  kcc=22
280  kfres=25
281 
282  ELSEIF(isub.EQ.6) THEN
283 C...Z0 + W+/- -> W+/-
284 
285  ELSEIF(isub.EQ.7) THEN
286 C...W+ + W- -> Z0
287 
288  ELSEIF(isub.EQ.8) THEN
289 C...W+ + W- -> h0
290  xh=sh/shp
291  230 DO 260 jt=1,2
292  i=mint(14+jt)
293  ia=iabs(i)
294  IF(ia.LE.10) THEN
295  rvckm=vint(180+i)*pyr(0)
296  DO 240 j=1,mstp(1)
297  ib=2*j-1+mod(ia,2)
298  ipm=(5-isign(1,i))/2
299  idc=j+mdcy(ia,2)+2
300  IF(mdme(idc,1).NE.1.AND.mdme(idc,1).NE.ipm) goto 240
301  mint(20+jt)=isign(ib,i)
302  rvckm=rvckm-vckm((ia+1)/2,(ib+1)/2)
303  IF(rvckm.LE.0d0) goto 250
304  240 CONTINUE
305  ELSE
306  ib=2*((ia+1)/2)-1+mod(ia,2)
307  mint(20+jt)=isign(ib,i)
308  ENDIF
309  250 pmq(jt)=pymass(mint(20+jt))
310  260 CONTINUE
311  jt=int(1.5d0+pyr(0))
312  zmin=2d0*pmq(jt)/shpr
313  zmax=1d0-pmq(3-jt)/shpr-(sh-pmq(jt)**2)/
314  & (shpr*(shpr-pmq(3-jt)))
315  zmax=min(1d0-xh,zmax)
316  IF(zmin.GE.zmax) goto 230
317  z(jt)=zmin+(zmax-zmin)*pyr(0)
318  IF(-1d0+(1d0+xh)/(1d0-z(jt))-xh/(1d0-z(jt))**2.LT.
319  & (1d0-xh)**2/(4d0*xh)*pyr(0)) goto 230
320  sqc1=1d0-4d0*pmq(jt)**2/(z(jt)**2*shp)
321  IF(sqc1.LT.1d-8) goto 230
322  c1=sqrt(sqc1)
323  c2=1d0+2d0*(pmas(24,1)**2-pmq(jt)**2)/(z(jt)*shp)
324  cthe(jt)=(c2-(c2**2-c1**2)/(c2+(2d0*pyr(0)-1d0)*c1))/c1
325  cthe(jt)=min(1d0,max(-1d0,cthe(jt)))
326  z(3-jt)=1d0-xh/(1d0-z(jt))
327  sqc1=1d0-4d0*pmq(3-jt)**2/(z(3-jt)**2*shp)
328  IF(sqc1.LT.1d-8) goto 230
329  c1=sqrt(sqc1)
330  c2=1d0+2d0*(pmas(24,1)**2-pmq(3-jt)**2)/(z(3-jt)*shp)
331  cthe(3-jt)=(c2-(c2**2-c1**2)/(c2+(2d0*pyr(0)-1d0)*c1))/c1
332  cthe(3-jt)=min(1d0,max(-1d0,cthe(3-jt)))
333  phir=paru(2)*pyr(0)
334  cphi=cos(phir)
335  ang=cthe(1)*cthe(2)-sqrt(1d0-cthe(1)**2)*
336  & sqrt(1d0-cthe(2)**2)*cphi
337  z1=2d0-z(jt)
338  z2=ang*sqrt(z(jt)**2-4d0*pmq(jt)**2/shp)
339  z3=1d0-z(jt)-xh+(pmq(1)**2+pmq(2)**2)/shp
340  z(3-jt)=2d0/(z1**2-z2**2)*(z1*z3+z2*sqrt(z3**2-(z1**2-z2**2)*
341  & pmq(3-jt)**2/shp))
342  zmin=2d0*pmq(3-jt)/shpr
343  zmax=1d0-pmq(jt)/shpr-(sh-pmq(3-jt)**2)/(shpr*(shpr-pmq(jt)))
344  zmax=min(1d0-xh,zmax)
345  IF(z(3-jt).LT.zmin.OR.z(3-jt).GT.zmax) goto 230
346  kcc=22
347  kfres=25
348 
349  ELSEIF(isub.EQ.10) THEN
350 C...f + f' -> f + f' (gamma/Z/W exchange); th = (p(f)-p(f))**2
351  IF(mint(2).EQ.1) THEN
352  kcc=22
353  ELSE
354 C...W exchange: need to mix flavours according to CKM matrix
355  DO 280 jt=1,2
356  i=mint(14+jt)
357  ia=iabs(i)
358  IF(ia.LE.10) THEN
359  rvckm=vint(180+i)*pyr(0)
360  DO 270 j=1,mstp(1)
361  ib=2*j-1+mod(ia,2)
362  ipm=(5-isign(1,i))/2
363  idc=j+mdcy(ia,2)+2
364  IF(mdme(idc,1).NE.1.AND.mdme(idc,1).NE.ipm) goto 270
365  mint(20+jt)=isign(ib,i)
366  rvckm=rvckm-vckm((ia+1)/2,(ib+1)/2)
367  IF(rvckm.LE.0d0) goto 280
368  270 CONTINUE
369  ELSE
370  ib=2*((ia+1)/2)-1+mod(ia,2)
371  mint(20+jt)=isign(ib,i)
372  ENDIF
373  280 CONTINUE
374  kcc=22
375  ENDIF
376  ENDIF
377 
378  ELSEIF(isub.LE.20) THEN
379  IF(isub.EQ.11) THEN
380 C...f + f' -> f + f' (g exchange); th = (p(f)-p(f))**2
381  kcc=mint(2)
382  IF(mint(15)*mint(16).LT.0) kcc=kcc+2
383 
384  ELSEIF(isub.EQ.12) THEN
385 C...f + fbar -> f' + fbar'; th = (p(f)-p(f'))**2
386  mint(21)=isign(kflf,mint(15))
387  mint(22)=-mint(21)
388  kcc=4
389 
390  ELSEIF(isub.EQ.13) THEN
391 C...f + fbar -> g + g; th arbitrary
392  mint(21)=21
393  mint(22)=21
394  kcc=mint(2)+4
395 
396  ELSEIF(isub.EQ.14) THEN
397 C...f + fbar -> g + gamma; th arbitrary
398  IF(pyr(0).GT.0.5d0) js=2
399  mint(20+js)=21
400  mint(23-js)=22
401  kcc=17+js
402 
403  ELSEIF(isub.EQ.15) THEN
404 C...f + fbar -> g + Z0; th arbitrary
405  IF(pyr(0).GT.0.5d0) js=2
406  mint(20+js)=21
407  mint(23-js)=23
408  kcc=17+js
409 
410  ELSEIF(isub.EQ.16) THEN
411 C...f + fbar' -> g + W+/-; th = (p(f)-p(W-))**2 or (p(fbar')-p(W+))**2
412  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
413  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
414  IF(mint(15)*(kch1+kch2).LT.0) js=2
415  mint(20+js)=21
416  mint(23-js)=isign(24,kch1+kch2)
417  kcc=17+js
418 
419  ELSEIF(isub.EQ.17) THEN
420 C...f + fbar -> g + h0; th arbitrary
421  IF(pyr(0).GT.0.5d0) js=2
422  mint(20+js)=21
423  mint(23-js)=25
424  kcc=17+js
425 
426  ELSEIF(isub.EQ.18) THEN
427 C...f + fbar -> gamma + gamma; th arbitrary
428  mint(21)=22
429  mint(22)=22
430 
431  ELSEIF(isub.EQ.19) THEN
432 C...f + fbar -> gamma + Z0; th arbitrary
433  IF(pyr(0).GT.0.5d0) js=2
434  mint(20+js)=22
435  mint(23-js)=23
436 
437  ELSEIF(isub.EQ.20) THEN
438 C...f + fbar' -> gamma + W+/-; th = (p(f)-p(W-))**2 or
439 C...(p(fbar')-p(W+))**2
440  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
441  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
442  IF(mint(15)*(kch1+kch2).LT.0) js=2
443  mint(20+js)=22
444  mint(23-js)=isign(24,kch1+kch2)
445  ENDIF
446 
447  ELSEIF(isub.LE.30) THEN
448  IF(isub.EQ.21) THEN
449 C...f + fbar -> gamma + h0; th arbitrary
450  IF(pyr(0).GT.0.5d0) js=2
451  mint(20+js)=22
452  mint(23-js)=25
453 
454  ELSEIF(isub.EQ.22) THEN
455 C...f + fbar -> Z0 + Z0; th arbitrary
456  mint(21)=23
457  mint(22)=23
458 
459  ELSEIF(isub.EQ.23) THEN
460 C...f + fbar' -> Z0 + W+/-; th = (p(f)-p(W-))**2 or (p(fbar')-p(W+))**2
461  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
462  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
463  IF(mint(15)*(kch1+kch2).LT.0) js=2
464  mint(20+js)=23
465  mint(23-js)=isign(24,kch1+kch2)
466 
467  ELSEIF(isub.EQ.24) THEN
468 C...f + fbar -> Z0 + h0 (or H0, or A0); th arbitrary
469  IF(pyr(0).GT.0.5d0) js=2
470  mint(20+js)=23
471  mint(23-js)=kfhigg
472 
473  ELSEIF(isub.EQ.25) THEN
474 C...f + fbar -> W+ + W-; th = (p(f)-p(W-))**2
475  mint(21)=-isign(24,mint(15))
476  mint(22)=-mint(21)
477 
478  ELSEIF(isub.EQ.26) THEN
479 C...f + fbar' -> W+/- + h0 (or H0, or A0);
480 C...th = (p(f)-p(W-))**2 or (p(fbar')-p(W+))**2
481  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
482  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
483  IF(mint(15)*(kch1+kch2).GT.0) js=2
484  mint(20+js)=isign(24,kch1+kch2)
485  mint(23-js)=kfhigg
486 
487  ELSEIF(isub.EQ.27) THEN
488 C...f + fbar -> h0 + h0
489 
490  ELSEIF(isub.EQ.28) THEN
491 C...f + g -> f + g; th = (p(f)-p(f))**2
492  IF(mint(15).EQ.21) js=2
493  kcc=mint(2)+6
494  IF(mint(15).EQ.21) kcc=kcc+2
495  IF(mint(15).NE.21) kcs=isign(1,mint(15))
496  IF(mint(16).NE.21) kcs=isign(1,mint(16))
497 
498  ELSEIF(isub.EQ.29) THEN
499 C...f + g -> f + gamma; th = (p(f)-p(f))**2
500  IF(mint(15).EQ.21) js=2
501  mint(23-js)=22
502  kcc=15+js
503  kcs=isign(1,mint(14+js))
504 
505  ELSEIF(isub.EQ.30) THEN
506 C...f + g -> f + Z0; th = (p(f)-p(f))**2
507  IF(mint(15).EQ.21) js=2
508  mint(23-js)=23
509  kcc=15+js
510  kcs=isign(1,mint(14+js))
511  ENDIF
512 
513  ELSEIF(isub.LE.40) THEN
514  IF(isub.EQ.31) THEN
515 C...f + g -> f' + W+/-; th = (p(f)-p(f'))**2; choose flavour f'
516  IF(mint(15).EQ.21) js=2
517  i=mint(14+js)
518  ia=iabs(i)
519  mint(23-js)=isign(24,kchg(ia,1)*i)
520  rvckm=vint(180+i)*pyr(0)
521  DO 290 j=1,mstp(1)
522  ib=2*j-1+mod(ia,2)
523  ipm=(5-isign(1,i))/2
524  idc=j+mdcy(ia,2)+2
525  IF(mdme(idc,1).NE.1.AND.mdme(idc,1).NE.ipm) goto 290
526  mint(20+js)=isign(ib,i)
527  rvckm=rvckm-vckm((ia+1)/2,(ib+1)/2)
528  IF(rvckm.LE.0d0) goto 300
529  290 CONTINUE
530  300 kcc=15+js
531  kcs=isign(1,mint(14+js))
532 
533  ELSEIF(isub.EQ.32) THEN
534 C...f + g -> f + h0; th = (p(f)-p(f))**2
535  IF(mint(15).EQ.21) js=2
536  mint(23-js)=25
537  kcc=15+js
538  kcs=isign(1,mint(14+js))
539 
540  ELSEIF(isub.EQ.33) THEN
541 C...f + gamma -> f + g; th=(p(f)-p(f))**2
542  IF(mint(15).EQ.22) js=2
543  mint(23-js)=21
544  kcc=24+js
545  kcs=isign(1,mint(14+js))
546 
547  ELSEIF(isub.EQ.34) THEN
548 C...f + gamma -> f + gamma; th=(p(f)-p(f))**2
549  IF(mint(15).EQ.22) js=2
550  kcc=22
551  kcs=isign(1,mint(14+js))
552 
553  ELSEIF(isub.EQ.35) THEN
554 C...f + gamma -> f + Z0; th=(p(f)-p(f))**2
555  IF(mint(15).EQ.22) js=2
556  mint(23-js)=23
557  kcc=22
558 
559  ELSEIF(isub.EQ.36) THEN
560 C...f + gamma -> f' + W+/-; th=(p(f)-p(f'))**2
561  IF(mint(15).EQ.22) js=2
562  i=mint(14+js)
563  ia=iabs(i)
564  mint(23-js)=isign(24,kchg(ia,1)*i)
565  IF(ia.LE.10) THEN
566  rvckm=vint(180+i)*pyr(0)
567  DO 310 j=1,mstp(1)
568  ib=2*j-1+mod(ia,2)
569  ipm=(5-isign(1,i))/2
570  idc=j+mdcy(ia,2)+2
571  IF(mdme(idc,1).NE.1.AND.mdme(idc,1).NE.ipm) goto 310
572  mint(20+js)=isign(ib,i)
573  rvckm=rvckm-vckm((ia+1)/2,(ib+1)/2)
574  IF(rvckm.LE.0d0) goto 320
575  310 CONTINUE
576  ELSE
577  ib=2*((ia+1)/2)-1+mod(ia,2)
578  mint(20+js)=isign(ib,i)
579  ENDIF
580  320 kcc=22
581 
582  ELSEIF(isub.EQ.37) THEN
583 C...f + gamma -> f + h0
584 
585  ELSEIF(isub.EQ.38) THEN
586 C...f + Z0 -> f + g
587 
588  ELSEIF(isub.EQ.39) THEN
589 C...f + Z0 -> f + gamma
590 
591  ELSEIF(isub.EQ.40) THEN
592 C...f + Z0 -> f + Z0
593  ENDIF
594 
595  ELSEIF(isub.LE.50) THEN
596  IF(isub.EQ.41) THEN
597 C...f + Z0 -> f' + W+/-
598 
599  ELSEIF(isub.EQ.42) THEN
600 C...f + Z0 -> f + h0
601 
602  ELSEIF(isub.EQ.43) THEN
603 C...f + W+/- -> f' + g
604 
605  ELSEIF(isub.EQ.44) THEN
606 C...f + W+/- -> f' + gamma
607 
608  ELSEIF(isub.EQ.45) THEN
609 C...f + W+/- -> f' + Z0
610 
611  ELSEIF(isub.EQ.46) THEN
612 C...f + W+/- -> f' + W+/-
613 
614  ELSEIF(isub.EQ.47) THEN
615 C...f + W+/- -> f' + h0
616 
617  ELSEIF(isub.EQ.48) THEN
618 C...f + h0 -> f + g
619 
620  ELSEIF(isub.EQ.49) THEN
621 C...f + h0 -> f + gamma
622 
623  ELSEIF(isub.EQ.50) THEN
624 C...f + h0 -> f + Z0
625  ENDIF
626 
627  ELSEIF(isub.LE.60) THEN
628  IF(isub.EQ.51) THEN
629 C...f + h0 -> f' + W+/-
630 
631  ELSEIF(isub.EQ.52) THEN
632 C...f + h0 -> f + h0
633 
634  ELSEIF(isub.EQ.53) THEN
635 C...g + g -> f + fbar; th arbitrary
636  kcs=(-1)**int(1.5d0+pyr(0))
637  mint(21)=isign(kflf,kcs)
638  mint(22)=-mint(21)
639  kcc=mint(2)+10
640 
641  ELSEIF(isub.EQ.54) THEN
642 C...g + gamma -> f + fbar; th arbitrary
643  kcs=(-1)**int(1.5d0+pyr(0))
644  mint(21)=isign(kflf,kcs)
645  mint(22)=-mint(21)
646  kcc=27
647  IF(mint(16).EQ.21) kcc=28
648 
649  ELSEIF(isub.EQ.55) THEN
650 C...g + Z0 -> f + fbar
651 
652  ELSEIF(isub.EQ.56) THEN
653 C...g + W+/- -> f + fbar'
654 
655  ELSEIF(isub.EQ.57) THEN
656 C...g + h0 -> f + fbar
657 
658  ELSEIF(isub.EQ.58) THEN
659 C...gamma + gamma -> f + fbar; th arbitrary
660  kcs=(-1)**int(1.5d0+pyr(0))
661  mint(21)=isign(kflf,kcs)
662  mint(22)=-mint(21)
663  kcc=21
664 
665  ELSEIF(isub.EQ.59) THEN
666 C...gamma + Z0 -> f + fbar
667 
668  ELSEIF(isub.EQ.60) THEN
669 C...gamma + W+/- -> f + fbar'
670  ENDIF
671 
672  ELSEIF(isub.LE.70) THEN
673  IF(isub.EQ.61) THEN
674 C...gamma + h0 -> f + fbar
675 
676  ELSEIF(isub.EQ.62) THEN
677 C...Z0 + Z0 -> f + fbar
678 
679  ELSEIF(isub.EQ.63) THEN
680 C...Z0 + W+/- -> f + fbar'
681 
682  ELSEIF(isub.EQ.64) THEN
683 C...Z0 + h0 -> f + fbar
684 
685  ELSEIF(isub.EQ.65) THEN
686 C...W+ + W- -> f + fbar
687 
688  ELSEIF(isub.EQ.66) THEN
689 C...W+/- + h0 -> f + fbar'
690 
691  ELSEIF(isub.EQ.67) THEN
692 C...h0 + h0 -> f + fbar
693 
694  ELSEIF(isub.EQ.68) THEN
695 C...g + g -> g + g; th arbitrary
696  kcc=mint(2)+12
697  kcs=(-1)**int(1.5d0+pyr(0))
698 
699  ELSEIF(isub.EQ.69) THEN
700 C...gamma + gamma -> W+ + W-; th arbitrary
701  mint(21)=24
702  mint(22)=-24
703  kcc=21
704 
705  ELSEIF(isub.EQ.70) THEN
706 C...gamma + W+/- -> Z0 + W+/-; th=(p(W)-p(W))**2
707  IF(mint(15).EQ.22) mint(21)=23
708  IF(mint(16).EQ.22) mint(22)=23
709  kcc=21
710  ENDIF
711 
712  ELSEIF(isub.LE.80) THEN
713  IF(isub.EQ.71.OR.isub.EQ.72) THEN
714 C...Z0 + Z0 -> Z0 + Z0; Z0 + Z0 -> W+ + W-
715  xh=sh/shp
716  mint(21)=mint(15)
717  mint(22)=mint(16)
718  pmq(1)=pymass(mint(21))
719  pmq(2)=pymass(mint(22))
720  330 jt=int(1.5d0+pyr(0))
721  zmin=2d0*pmq(jt)/shpr
722  zmax=1d0-pmq(3-jt)/shpr-(sh-pmq(jt)**2)/
723  & (shpr*(shpr-pmq(3-jt)))
724  zmax=min(1d0-xh,zmax)
725  z(jt)=zmin+(zmax-zmin)*pyr(0)
726  IF(-1d0+(1d0+xh)/(1d0-z(jt))-xh/(1d0-z(jt))**2.LT.
727  & (1d0-xh)**2/(4d0*xh)*pyr(0)) goto 330
728  sqc1=1d0-4d0*pmq(jt)**2/(z(jt)**2*shp)
729  IF(sqc1.LT.1d-8) goto 330
730  c1=sqrt(sqc1)
731  c2=1d0+2d0*(pmas(23,1)**2-pmq(jt)**2)/(z(jt)*shp)
732  cthe(jt)=(c2-(c2**2-c1**2)/(c2+(2d0*pyr(0)-1d0)*c1))/c1
733  cthe(jt)=min(1d0,max(-1d0,cthe(jt)))
734  z(3-jt)=1d0-xh/(1d0-z(jt))
735  sqc1=1d0-4d0*pmq(3-jt)**2/(z(3-jt)**2*shp)
736  IF(sqc1.LT.1d-8) goto 330
737  c1=sqrt(sqc1)
738  c2=1d0+2d0*(pmas(23,1)**2-pmq(3-jt)**2)/(z(3-jt)*shp)
739  cthe(3-jt)=(c2-(c2**2-c1**2)/(c2+(2d0*pyr(0)-1d0)*c1))/c1
740  cthe(3-jt)=min(1d0,max(-1d0,cthe(3-jt)))
741  phir=paru(2)*pyr(0)
742  cphi=cos(phir)
743  ang=cthe(1)*cthe(2)-sqrt(1d0-cthe(1)**2)*
744  & sqrt(1d0-cthe(2)**2)*cphi
745  z1=2d0-z(jt)
746  z2=ang*sqrt(z(jt)**2-4d0*pmq(jt)**2/shp)
747  z3=1d0-z(jt)-xh+(pmq(1)**2+pmq(2)**2)/shp
748  z(3-jt)=2d0/(z1**2-z2**2)*(z1*z3+z2*sqrt(z3**2-(z1**2-z2**2)*
749  & pmq(3-jt)**2/shp))
750  zmin=2d0*pmq(3-jt)/shpr
751  zmax=1d0-pmq(jt)/shpr-(sh-pmq(3-jt)**2)/(shpr*(shpr-pmq(jt)))
752  zmax=min(1d0-xh,zmax)
753  IF(z(3-jt).LT.zmin.OR.z(3-jt).GT.zmax) goto 330
754  kcc=22
755 
756  ELSEIF(isub.EQ.73) THEN
757 C...Z0 + W+/- -> Z0 + W+/-
758  js=mint(2)
759  xh=sh/shp
760  340 jt=3-mint(2)
761  i=mint(14+jt)
762  ia=iabs(i)
763  IF(ia.LE.10) THEN
764  rvckm=vint(180+i)*pyr(0)
765  DO 350 j=1,mstp(1)
766  ib=2*j-1+mod(ia,2)
767  ipm=(5-isign(1,i))/2
768  idc=j+mdcy(ia,2)+2
769  IF(mdme(idc,1).NE.1.AND.mdme(idc,1).NE.ipm) goto 350
770  mint(20+jt)=isign(ib,i)
771  rvckm=rvckm-vckm((ia+1)/2,(ib+1)/2)
772  IF(rvckm.LE.0d0) goto 360
773  350 CONTINUE
774  ELSE
775  ib=2*((ia+1)/2)-1+mod(ia,2)
776  mint(20+jt)=isign(ib,i)
777  ENDIF
778  360 pmq(jt)=pymass(mint(20+jt))
779  mint(23-jt)=mint(17-jt)
780  pmq(3-jt)=pymass(mint(23-jt))
781  jt=int(1.5d0+pyr(0))
782  zmin=2d0*pmq(jt)/shpr
783  zmax=1d0-pmq(3-jt)/shpr-(sh-pmq(jt)**2)/
784  & (shpr*(shpr-pmq(3-jt)))
785  zmax=min(1d0-xh,zmax)
786  IF(zmin.GE.zmax) goto 340
787  z(jt)=zmin+(zmax-zmin)*pyr(0)
788  IF(-1d0+(1d0+xh)/(1d0-z(jt))-xh/(1d0-z(jt))**2.LT.
789  & (1d0-xh)**2/(4d0*xh)*pyr(0)) goto 340
790  sqc1=1d0-4d0*pmq(jt)**2/(z(jt)**2*shp)
791  IF(sqc1.LT.1d-8) goto 340
792  c1=sqrt(sqc1)
793  c2=1d0+2d0*(pmas(23,1)**2-pmq(jt)**2)/(z(jt)*shp)
794  cthe(jt)=(c2-(c2**2-c1**2)/(c2+(2d0*pyr(0)-1d0)*c1))/c1
795  cthe(jt)=min(1d0,max(-1d0,cthe(jt)))
796  z(3-jt)=1d0-xh/(1d0-z(jt))
797  sqc1=1d0-4d0*pmq(3-jt)**2/(z(3-jt)**2*shp)
798  IF(sqc1.LT.1d-8) goto 340
799  c1=sqrt(sqc1)
800  c2=1d0+2d0*(pmas(23,1)**2-pmq(3-jt)**2)/(z(3-jt)*shp)
801  cthe(3-jt)=(c2-(c2**2-c1**2)/(c2+(2d0*pyr(0)-1d0)*c1))/c1
802  cthe(3-jt)=min(1d0,max(-1d0,cthe(3-jt)))
803  phir=paru(2)*pyr(0)
804  cphi=cos(phir)
805  ang=cthe(1)*cthe(2)-sqrt(1d0-cthe(1)**2)*
806  & sqrt(1d0-cthe(2)**2)*cphi
807  z1=2d0-z(jt)
808  z2=ang*sqrt(z(jt)**2-4d0*pmq(jt)**2/shp)
809  z3=1d0-z(jt)-xh+(pmq(1)**2+pmq(2)**2)/shp
810  z(3-jt)=2d0/(z1**2-z2**2)*(z1*z3+z2*sqrt(z3**2-(z1**2-z2**2)*
811  & pmq(3-jt)**2/shp))
812  zmin=2d0*pmq(3-jt)/shpr
813  zmax=1d0-pmq(jt)/shpr-(sh-pmq(3-jt)**2)/(shpr*(shpr-pmq(jt)))
814  zmax=min(1d0-xh,zmax)
815  IF(z(3-jt).LT.zmin.OR.z(3-jt).GT.zmax) goto 340
816  kcc=22
817 
818  ELSEIF(isub.EQ.74) THEN
819 C...Z0 + h0 -> Z0 + h0
820 
821  ELSEIF(isub.EQ.75) THEN
822 C...W+ + W- -> gamma + gamma
823 
824  ELSEIF(isub.EQ.76.OR.isub.EQ.77) THEN
825 C...W+ + W- -> Z0 + Z0; W+ + W- -> W+ + W-
826  xh=sh/shp
827  370 DO 400 jt=1,2
828  i=mint(14+jt)
829  ia=iabs(i)
830  IF(ia.LE.10) THEN
831  rvckm=vint(180+i)*pyr(0)
832  DO 380 j=1,mstp(1)
833  ib=2*j-1+mod(ia,2)
834  ipm=(5-isign(1,i))/2
835  idc=j+mdcy(ia,2)+2
836  IF(mdme(idc,1).NE.1.AND.mdme(idc,1).NE.ipm) goto 380
837  mint(20+jt)=isign(ib,i)
838  rvckm=rvckm-vckm((ia+1)/2,(ib+1)/2)
839  IF(rvckm.LE.0d0) goto 390
840  380 CONTINUE
841  ELSE
842  ib=2*((ia+1)/2)-1+mod(ia,2)
843  mint(20+jt)=isign(ib,i)
844  ENDIF
845  390 pmq(jt)=pymass(mint(20+jt))
846  400 CONTINUE
847  jt=int(1.5d0+pyr(0))
848  zmin=2d0*pmq(jt)/shpr
849  zmax=1d0-pmq(3-jt)/shpr-(sh-pmq(jt)**2)/
850  & (shpr*(shpr-pmq(3-jt)))
851  zmax=min(1d0-xh,zmax)
852  IF(zmin.GE.zmax) goto 370
853  z(jt)=zmin+(zmax-zmin)*pyr(0)
854  IF(-1d0+(1d0+xh)/(1d0-z(jt))-xh/(1d0-z(jt))**2.LT.
855  & (1d0-xh)**2/(4d0*xh)*pyr(0)) goto 370
856  sqc1=1d0-4d0*pmq(jt)**2/(z(jt)**2*shp)
857  IF(sqc1.LT.1d-8) goto 370
858  c1=sqrt(sqc1)
859  c2=1d0+2d0*(pmas(24,1)**2-pmq(jt)**2)/(z(jt)*shp)
860  cthe(jt)=(c2-(c2**2-c1**2)/(c2+(2d0*pyr(0)-1d0)*c1))/c1
861  cthe(jt)=min(1d0,max(-1d0,cthe(jt)))
862  z(3-jt)=1d0-xh/(1d0-z(jt))
863  sqc1=1d0-4d0*pmq(3-jt)**2/(z(3-jt)**2*shp)
864  IF(sqc1.LT.1d-8) goto 370
865  c1=sqrt(sqc1)
866  c2=1d0+2d0*(pmas(24,1)**2-pmq(3-jt)**2)/(z(3-jt)*shp)
867  cthe(3-jt)=(c2-(c2**2-c1**2)/(c2+(2d0*pyr(0)-1d0)*c1))/c1
868  cthe(3-jt)=min(1d0,max(-1d0,cthe(3-jt)))
869  phir=paru(2)*pyr(0)
870  cphi=cos(phir)
871  ang=cthe(1)*cthe(2)-sqrt(1d0-cthe(1)**2)*
872  & sqrt(1d0-cthe(2)**2)*cphi
873  z1=2d0-z(jt)
874  z2=ang*sqrt(z(jt)**2-4d0*pmq(jt)**2/shp)
875  z3=1d0-z(jt)-xh+(pmq(1)**2+pmq(2)**2)/shp
876  z(3-jt)=2d0/(z1**2-z2**2)*(z1*z3+z2*sqrt(z3**2-(z1**2-z2**2)*
877  & pmq(3-jt)**2/shp))
878  zmin=2d0*pmq(3-jt)/shpr
879  zmax=1d0-pmq(jt)/shpr-(sh-pmq(3-jt)**2)/(shpr*(shpr-pmq(jt)))
880  zmax=min(1d0-xh,zmax)
881  IF(z(3-jt).LT.zmin.OR.z(3-jt).GT.zmax) goto 370
882  kcc=22
883 
884  ELSEIF(isub.EQ.78) THEN
885 C...W+/- + h0 -> W+/- + h0
886 
887  ELSEIF(isub.EQ.79) THEN
888 C...h0 + h0 -> h0 + h0
889 
890  ELSEIF(isub.EQ.80) THEN
891 C...q + gamma -> q' + pi+/-; th=(p(q)-p(q'))**2
892  IF(mint(15).EQ.22) js=2
893  i=mint(14+js)
894  ia=iabs(i)
895  mint(23-js)=isign(211,kchg(ia,1)*i)
896  ib=3-ia
897  mint(20+js)=isign(ib,i)
898  kcc=22
899  ENDIF
900 
901  ELSEIF(isub.LE.90) THEN
902  IF(isub.EQ.81) THEN
903 C...q + qbar -> Q + Qbar; th = (p(q)-p(Q))**2
904  mint(21)=isign(mint(55),mint(15))
905  mint(22)=-mint(21)
906  kcc=4
907 
908  ELSEIF(isub.EQ.82) THEN
909 C...g + g -> Q + Qbar; th arbitrary
910  kcs=(-1)**int(1.5d0+pyr(0))
911  mint(21)=isign(mint(55),kcs)
912  mint(22)=-mint(21)
913  kcc=mint(2)+10
914 
915  ELSEIF(isub.EQ.83) THEN
916 C...f + q -> f' + Q; th = (p(f) - p(f'))**2
917  kfold=mint(16)
918  IF(mint(2).EQ.2) kfold=mint(15)
919  kfaold=iabs(kfold)
920  IF(kfaold.GT.10) THEN
921  kfanew=kfaold+2*mod(kfaold,2)-1
922  ELSE
923  rckm=vint(180+kfold)*pyr(0)
924  ipm=(5-isign(1,kfold))/2
925  kfanew=-mod(kfaold+1,2)
926  410 kfanew=kfanew+2
927  idc=mdcy(kfaold,2)+(kfanew+1)/2+2
928  IF(mdme(idc,1).EQ.1.OR.mdme(idc,1).EQ.ipm) THEN
929  IF(mod(kfaold,2).EQ.0) rckm=rckm-
930  & vckm(kfaold/2,(kfanew+1)/2)
931  IF(mod(kfaold,2).EQ.1) rckm=rckm-
932  & vckm(kfanew/2,(kfaold+1)/2)
933  ENDIF
934  IF(kfanew.LE.6.AND.rckm.GT.0d0) goto 410
935  ENDIF
936  IF(mint(2).EQ.1) THEN
937  mint(21)=isign(mint(55),mint(15))
938  mint(22)=isign(kfanew,mint(16))
939  ELSE
940  mint(21)=isign(kfanew,mint(15))
941  mint(22)=isign(mint(55),mint(16))
942  js=2
943  ENDIF
944  kcc=22
945 
946  ELSEIF(isub.EQ.84) THEN
947 C...g + gamma -> Q + Qbar; th arbitary
948  kcs=(-1)**int(1.5d0+pyr(0))
949  mint(21)=isign(mint(55),kcs)
950  mint(22)=-mint(21)
951  kcc=27
952  IF(mint(16).EQ.21) kcc=28
953 
954  ELSEIF(isub.EQ.85) THEN
955 C...gamma + gamma -> F + Fbar; th arbitary
956  kcs=(-1)**int(1.5d0+pyr(0))
957  mint(21)=isign(mint(56),kcs)
958  mint(22)=-mint(21)
959  kcc=21
960 
961  ELSEIF(isub.GE.86.AND.isub.LE.89) THEN
962 C...g + g -> (J/Psi, chi_0c, chi_1c or chi_2c) + g
963  mint(21)=kfpr(isub,1)
964  mint(22)=kfpr(isub,2)
965  kcc=24
966  kcs=(-1)**int(1.5d0+pyr(0))
967  ENDIF
968 
969  ELSEIF(isub.LE.100) THEN
970  IF(isub.EQ.95) THEN
971 C...Low-pT ( = energyless g + g -> g + g)
972  kcc=mint(2)+12
973  kcs=(-1)**int(1.5d0+pyr(0))
974 
975  ELSEIF(isub.EQ.96) THEN
976 C...Multiple interactions (should be reassigned to QCD process)
977  ENDIF
978 
979  ELSEIF(isub.LE.110) THEN
980  IF(isub.EQ.101) THEN
981 C...g + g -> gamma*/Z0
982  kcc=21
983  kfres=22
984 
985  ELSEIF(isub.EQ.102) THEN
986 C...g + g -> h0 (or H0, or A0)
987  kcc=21
988  kfres=kfhigg
989 
990  ELSEIF(isub.EQ.103) THEN
991 C...gamma + gamma -> h0 (or H0, or A0)
992  kcc=21
993  kfres=kfhigg
994 
995  ELSEIF(isub.EQ.104.OR.isub.EQ.105) THEN
996 C...g + g -> chi_0c or chi_2c.
997  kcc=21
998  kfres=kfpr(isub,1)
999 
1000  ELSEIF(isub.EQ.106) THEN
1001 C...g + g -> J/Psi + gamma
1002  mint(21)=kfpr(isub,1)
1003  mint(22)=kfpr(isub,2)
1004  kcc=21
1005 
1006  ELSEIF(isub.EQ.107) THEN
1007 C...g + gamma -> J/Psi + g
1008  mint(21)=kfpr(isub,1)
1009  mint(22)=kfpr(isub,2)
1010  kcc=22
1011  IF(mint(16).EQ.22) kcc=33
1012 
1013  ELSEIF(isub.EQ.108) THEN
1014 C...gamma + gamma -> J/Psi + gamma
1015  mint(21)=kfpr(isub,1)
1016  mint(22)=kfpr(isub,2)
1017 
1018  ELSEIF(isub.EQ.110) THEN
1019 C...f + fbar -> gamma + h0; th arbitrary
1020  IF(pyr(0).GT.0.5d0) js=2
1021  mint(20+js)=22
1022  mint(23-js)=kfhigg
1023  ENDIF
1024 
1025  ELSEIF(isub.LE.120) THEN
1026  IF(isub.EQ.111) THEN
1027 C...f + fbar -> g + h0; th arbitrary
1028  IF(pyr(0).GT.0.5d0) js=2
1029  mint(20+js)=21
1030  mint(23-js)=kfhigg
1031  kcc=17+js
1032 
1033  ELSEIF(isub.EQ.112) THEN
1034 C...f + g -> f + h0; th = (p(f) - p(f))**2
1035  IF(mint(15).EQ.21) js=2
1036  mint(23-js)=kfhigg
1037  kcc=15+js
1038  kcs=isign(1,mint(14+js))
1039 
1040  ELSEIF(isub.EQ.113) THEN
1041 C...g + g -> g + h0; th arbitrary
1042  IF(pyr(0).GT.0.5d0) js=2
1043  mint(23-js)=kfhigg
1044  kcc=22+js
1045  kcs=(-1)**int(1.5d0+pyr(0))
1046 
1047  ELSEIF(isub.EQ.114) THEN
1048 C...g + g -> gamma + gamma; th arbitrary
1049  IF(pyr(0).GT.0.5d0) js=2
1050  mint(21)=22
1051  mint(22)=22
1052  kcc=21
1053 
1054  ELSEIF(isub.EQ.115) THEN
1055 C...g + g -> g + gamma; th arbitrary
1056  IF(pyr(0).GT.0.5d0) js=2
1057  mint(23-js)=22
1058  kcc=22+js
1059  kcs=(-1)**int(1.5d0+pyr(0))
1060 
1061  ELSEIF(isub.EQ.116) THEN
1062 C...g + g -> gamma + Z0
1063 
1064  ELSEIF(isub.EQ.117) THEN
1065 C...g + g -> Z0 + Z0
1066 
1067  ELSEIF(isub.EQ.118) THEN
1068 C...g + g -> W+ + W-
1069  ENDIF
1070 
1071  ELSEIF(isub.LE.140) THEN
1072  IF(isub.EQ.121) THEN
1073 C...g + g -> Q + Qbar + h0
1074  kcs=(-1)**int(1.5d0+pyr(0))
1075  mint(21)=isign(kfpr(isubsv,2),kcs)
1076  mint(22)=-mint(21)
1077  kcc=11+int(0.5d0+pyr(0))
1078  kfres=kfhigg
1079 
1080  ELSEIF(isub.EQ.122) THEN
1081 C...q + qbar -> Q + Qbar + h0
1082  mint(21)=isign(kfpr(isubsv,2),mint(15))
1083  mint(22)=-mint(21)
1084  kcc=4
1085  kfres=kfhigg
1086 
1087  ELSEIF(isub.EQ.123) THEN
1088 C...f + f' -> f + f' + h0 (or H0, or A0) (Z0 + Z0 -> h0 as
1089 C...inner process)
1090  kcc=22
1091  kfres=kfhigg
1092 
1093  ELSEIF(isub.EQ.124) THEN
1094 C...f + f' -> f" + f"' + h0 (or H0, or A) (W+ + W- -> h0 as
1095 C...inner process)
1096  DO 430 jt=1,2
1097  i=mint(14+jt)
1098  ia=iabs(i)
1099  IF(ia.LE.10) THEN
1100  rvckm=vint(180+i)*pyr(0)
1101  DO 420 j=1,mstp(1)
1102  ib=2*j-1+mod(ia,2)
1103  ipm=(5-isign(1,i))/2
1104  idc=j+mdcy(ia,2)+2
1105  IF(mdme(idc,1).NE.1.AND.mdme(idc,1).NE.ipm) goto 420
1106  mint(20+jt)=isign(ib,i)
1107  rvckm=rvckm-vckm((ia+1)/2,(ib+1)/2)
1108  IF(rvckm.LE.0d0) goto 430
1109  420 CONTINUE
1110  ELSE
1111  ib=2*((ia+1)/2)-1+mod(ia,2)
1112  mint(20+jt)=isign(ib,i)
1113  ENDIF
1114  430 CONTINUE
1115  kcc=22
1116  kfres=kfhigg
1117 
1118  ELSEIF(isub.EQ.131.OR.isub.EQ.132) THEN
1119 C...f + gamma*_(T,L) -> f + g; th=(p(f)-p(f))**2
1120  IF(mint(15).EQ.22) js=2
1121  mint(23-js)=21
1122  kcc=24+js
1123  kcs=isign(1,mint(14+js))
1124 
1125  ELSEIF(isub.EQ.133.OR.isub.EQ.134) THEN
1126 C...f + gamma*_(T,L) -> f + gamma; th=(p(f)-p(f))**2
1127  IF(mint(15).EQ.22) js=2
1128  kcc=22
1129  kcs=isign(1,mint(14+js))
1130 
1131  ELSEIF(isub.EQ.135.OR.isub.EQ.136) THEN
1132 C...g + gamma*_(T,L) -> f + fbar; th arbitrary
1133  kcs=(-1)**int(1.5d0+pyr(0))
1134  mint(21)=isign(kflf,kcs)
1135  mint(22)=-mint(21)
1136  kcc=27
1137  IF(mint(16).EQ.21) kcc=28
1138 
1139  ELSEIF(isub.GE.137.AND.isub.LE.140) THEN
1140 C...gamma*_(T,L) + gamma*_(T,L) -> f + fbar; th arbitrary
1141  kcs=(-1)**int(1.5d0+pyr(0))
1142  mint(21)=isign(kflf,kcs)
1143  mint(22)=-mint(21)
1144  kcc=21
1145 
1146  ENDIF
1147 
1148  ELSEIF(isub.LE.160) THEN
1149  IF(isub.EQ.141) THEN
1150 C...f + fbar -> gamma*/Z0/Z'0
1151  kfres=32
1152 
1153  ELSEIF(isub.EQ.142) THEN
1154 C...f + fbar' -> W'+/-
1155  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
1156  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
1157  kfres=isign(34,kch1+kch2)
1158 
1159  ELSEIF(isub.EQ.143) THEN
1160 C...f + fbar' -> H+/-
1161  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
1162  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
1163  kfres=isign(37,kch1+kch2)
1164 
1165  ELSEIF(isub.EQ.144) THEN
1166 C...f + fbar' -> R
1167  kfres=isign(41,mint(15)+mint(16))
1168 
1169  ELSEIF(isub.EQ.145) THEN
1170 C...q + l -> LQ (leptoquark)
1171  IF(iabs(mint(16)).LE.8) js=2
1172  kfres=isign(42,mint(14+js))
1173  kcc=28+js
1174  kcs=isign(1,mint(14+js))
1175 
1176  ELSEIF(isub.EQ.146) THEN
1177 C...e + gamma -> e* (excited lepton)
1178  IF(mint(15).EQ.22) js=2
1179  kfres=isign(kfpr(isub,1),mint(14+js))
1180  kcc=22
1181 
1182  ELSEIF(isub.EQ.147.OR.isub.EQ.148) THEN
1183 C...q + g -> q* (excited quark)
1184  IF(mint(15).EQ.21) js=2
1185  kfres=isign(kfpr(isub,1),mint(14+js))
1186  kcc=30+js
1187  kcs=isign(1,mint(14+js))
1188 
1189  ELSEIF(isub.EQ.149) THEN
1190 C...g + g -> eta_tc
1191  kfres=ktechn+331
1192  kcc=23
1193  kcs=(-1)**int(1.5d0+pyr(0))
1194  ENDIF
1195 
1196  ELSEIF(isub.LE.200) THEN
1197  IF(isub.EQ.161) THEN
1198 C...f + g -> f' + H+/-; th = (p(f)-p(f'))**2
1199  IF(mint(15).EQ.21) js=2
1200  i=mint(14+js)
1201  ia=iabs(i)
1202  mint(23-js)=isign(37,kchg(ia,1)*i)
1203  ib=ia+mod(ia,2)-mod(ia+1,2)
1204  mint(20+js)=isign(ib,i)
1205  kcc=15+js
1206  kcs=isign(1,mint(14+js))
1207 
1208  ELSEIF(isub.EQ.162) THEN
1209 C...q + g -> LQ + lbar; LQ=leptoquark; th=(p(q)-p(LQ))^2
1210  IF(mint(15).EQ.21) js=2
1211  mint(20+js)=isign(42,mint(14+js))
1212  kflql=kfdp(mdcy(42,2),2)
1213  mint(23-js)=-isign(kflql,mint(14+js))
1214  kcc=15+js
1215  kcs=isign(1,mint(14+js))
1216 
1217  ELSEIF(isub.EQ.163) THEN
1218 C...g + g -> LQ + LQbar; LQ=leptoquark; th arbitrary
1219  kcs=(-1)**int(1.5d0+pyr(0))
1220  mint(21)=isign(42,kcs)
1221  mint(22)=-mint(21)
1222  kcc=mint(2)+10
1223 
1224  ELSEIF(isub.EQ.164) THEN
1225 C...q + qbar -> LQ + LQbar; LQ=leptoquark; th=(p(q)-p(LQ))**2
1226  mint(21)=isign(42,mint(15))
1227  mint(22)=-mint(21)
1228  kcc=4
1229 
1230  ELSEIF(isub.EQ.165) THEN
1231 C...q + qbar -> l- + l+; th=(p(q)-p(l-))**2
1232  mint(21)=isign(kfpr(isub,1),mint(15))
1233  mint(22)=-mint(21)
1234 
1235  ELSEIF(isub.EQ.166) THEN
1236 C...q + qbar' -> l + nu; th=(p(u)-p(nu))**2 or (p(ubar)-p(nubar))**2
1237  IF(mod(mint(15),2).EQ.0) THEN
1238  mint(21)=isign(kfpr(isub,1)+1,mint(15))
1239  mint(22)=isign(kfpr(isub,1),mint(16))
1240  ELSE
1241  mint(21)=isign(kfpr(isub,1),mint(15))
1242  mint(22)=isign(kfpr(isub,1)+1,mint(16))
1243  ENDIF
1244 
1245  ELSEIF(isub.EQ.167.OR.isub.EQ.168) THEN
1246 C...q + q' -> q" + q* (excited quark)
1247  kfqstr=kfpr(isub,2)
1248  kfqexc=mod(kfqstr,kexcit)
1249  js=mint(2)
1250  mint(20+js)=isign(kfqstr,mint(14+js))
1251  IF(iabs(mint(15)).NE.kfqexc.AND.iabs(mint(16)).NE.kfqexc)
1252  & mint(23-js)=isign(kfqexc,mint(17-js))
1253  kcc=22
1254  js=3-js
1255 
1256  ELSEIF(isub.EQ.169) THEN
1257 C...q + qbar -> e + e* (excited lepton)
1258  kfqstr=kfpr(isub,2)
1259  kfqexc=mod(kfqstr,kexcit)
1260  js=mint(2)
1261  mint(20+js)=isign(kfqstr,mint(14+js))
1262  mint(23-js)=isign(kfqexc,mint(17-js))
1263  js=3-js
1264 
1265  ELSEIF(isub.EQ.191) THEN
1266 C...f + fbar -> rho_tc0.
1267  kfres=ktechn+113
1268 
1269  ELSEIF(isub.EQ.192) THEN
1270 C...f + fbar' -> rho_tc+/-
1271  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
1272  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
1273  kfres=isign(ktechn+213,kch1+kch2)
1274 
1275  ELSEIF(isub.EQ.193) THEN
1276 C...f + fbar -> omega_tc0.
1277  kfres=ktechn+223
1278 
1279  ELSEIF(isub.EQ.194) THEN
1280 C...f + fbar -> f' + fbar' via mixture of s-channel
1281 C...rho_tc and omega_tc; th=(p(f)-p(f'))**2
1282  mint(21)=isign(kfpr(isub,1),mint(15))
1283  mint(22)=-mint(21)
1284 
1285  ELSEIF(isub.EQ.195) THEN
1286 C...f + fbar' -> f'' + fbar''' via s-channel
1287 C...rho_tc+ th=(p(f)-p(f'))**2
1288 C...q + qbar' -> l + nu; th=(p(u)-p(nu))**2 or (p(ubar)-p(nubar))**2
1289  IF(mod(mint(15),2).EQ.0) THEN
1290  mint(21)=isign(kfpr(isub,1)+1,mint(15))
1291  mint(22)=isign(kfpr(isub,1),mint(16))
1292  ELSE
1293  mint(21)=isign(kfpr(isub,1),mint(15))
1294  mint(22)=isign(kfpr(isub,1)+1,mint(16))
1295  ENDIF
1296  ENDIF
1297 
1298 CMRENNA++
1299  ELSEIF(isub.LE.215) THEN
1300  IF(isub.EQ.201) THEN
1301 C...f + fbar -> ~e_L + ~e_Lbar
1302  mint(21)=isign(ksusy1+11,kcs)
1303  mint(22)=-mint(21)
1304 
1305  ELSEIF(isub.EQ.202) THEN
1306 C...f + fbar -> ~e_R + ~e_Rbar
1307  mint(21)=isign(ksusy2+11,kcs)
1308  mint(22)=-mint(21)
1309 
1310  ELSEIF(isub.EQ.203) THEN
1311 C...f + fbar -> ~e_L + ~e_Rbar
1312  IF(mint(15).LT.0) js=2
1313  IF(mint(2).EQ.1) THEN
1314  mint(20+js)=kfpr(isub,1)
1315  mint(23-js)=-kfpr(isub,2)
1316  ELSE
1317  mint(20+js)=-kfpr(isub,1)
1318  mint(23-js)=kfpr(isub,2)
1319  ENDIF
1320 
1321  ELSEIF(isub.EQ.204) THEN
1322 C...f + fbar -> ~mu_L + ~mu_Lbar
1323  mint(21)=isign(ksusy1+13,kcs)
1324  mint(22)=-mint(21)
1325 
1326  ELSEIF(isub.EQ.205) THEN
1327 C...f + fbar -> ~mu_R + ~mu_Rbar
1328  mint(21)=isign(ksusy2+13,kcs)
1329  mint(22)=-mint(21)
1330 
1331  ELSEIF(isub.EQ.206) THEN
1332 C...f + fbar -> ~mu_L + ~mu_Rbar
1333  IF(mint(15).LT.0) js=2
1334  IF(mint(2).EQ.1) THEN
1335  mint(20+js)=kfpr(isub,1)
1336  mint(23-js)=-kfpr(isub,2)
1337  ELSE
1338  mint(20+js)=-kfpr(isub,1)
1339  mint(23-js)=kfpr(isub,2)
1340  ENDIF
1341 
1342  ELSEIF(isub.EQ.207) THEN
1343 C...f + fbar -> ~tau_1 + ~tau_1bar
1344  mint(21)=isign(ksusy1+15,kcs)
1345  mint(22)=-mint(21)
1346 
1347  ELSEIF(isub.EQ.208) THEN
1348 C...f + fbar -> ~tau_2 + ~tau_2bar
1349  mint(21)=isign(ksusy2+15,kcs)
1350  mint(22)=-mint(21)
1351 
1352  ELSEIF(isub.EQ.209) THEN
1353 C...f + fbar -> ~tau_1 + ~tau_2bar
1354  IF(mint(15).LT.0) js=2
1355  IF(mint(2).EQ.1) THEN
1356  mint(20+js)=kfpr(isub,1)
1357  mint(23-js)=-kfpr(isub,2)
1358  ELSE
1359  mint(20+js)=-kfpr(isub,1)
1360  mint(23-js)=kfpr(isub,2)
1361  ENDIF
1362 
1363  ELSEIF(isub.EQ.210) THEN
1364 C...q + qbar' -> ~l_L + ~nulbar; th arbitrary
1365  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
1366  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
1367  mint(21)=-isign(kfpr(isub,1),kch1+kch2)
1368  mint(22)=isign(kfpr(isub,2),kch1+kch2)
1369 
1370  ELSEIF(isub.EQ.211) THEN
1371 C...q + qbar'-> ~tau_1 + ~nutaubar; th arbitrary
1372  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
1373  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
1374  mint(21)=-isign(ksusy1+15,kch1+kch2)
1375  mint(22)=isign(ksusy1+16,kch1+kch2)
1376 
1377  ELSEIF(isub.EQ.212) THEN
1378 C...q + qbar'-> ~tau_2 + ~nutaubar; th arbitrary
1379  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
1380  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
1381  mint(21)=-isign(ksusy2+15,kch1+kch2)
1382  mint(22)=isign(ksusy1+16,kch1+kch2)
1383 
1384  ELSEIF(isub.EQ.213) THEN
1385 C...f + fbar -> ~nul + ~nulbar
1386  mint(21)=isign(kfpr(isub,1),kcs)
1387  mint(22)=-mint(21)
1388 
1389  ELSEIF(isub.EQ.214) THEN
1390 C...f + fbar -> ~nutau + ~nutaubar
1391  mint(21)=isign(ksusy1+16,kcs)
1392  mint(22)=-mint(21)
1393  ENDIF
1394 
1395  ELSEIF(isub.LE.225) THEN
1396  IF(isub.EQ.216) THEN
1397 C...f + fbar -> ~chi01 + ~chi01
1398  mint(21)=ksusy1+22
1399  mint(22)=ksusy1+22
1400 
1401  ELSEIF(isub.EQ.217) THEN
1402 C...f + fbar -> ~chi02 + ~chi02
1403  mint(21)=ksusy1+23
1404  mint(22)=ksusy1+23
1405 
1406  ELSEIF(isub.EQ.218 ) THEN
1407 C...f + fbar -> ~chi03 + ~chi03
1408  mint(21)=ksusy1+25
1409  mint(22)=ksusy1+25
1410 
1411  ELSEIF(isub.EQ.219 ) THEN
1412 C...f + fbar -> ~chi04 + ~chi04
1413  mint(21)=ksusy1+35
1414  mint(22)=ksusy1+35
1415 
1416  ELSEIF(isub.EQ.220 ) THEN
1417 C...f + fbar -> ~chi01 + ~chi02
1418  IF(mint(15).LT.0) js=2
1419 C IF(PYR(0).GT.0.5D0) JS=2
1420  mint(20+js)=ksusy1+22
1421  mint(23-js)=ksusy1+23
1422 
1423  ELSEIF(isub.EQ.221 ) THEN
1424 C...f + fbar -> ~chi01 + ~chi03
1425  IF(mint(15).LT.0) js=2
1426 C IF(PYR(0).GT.0.5D0) JS=2
1427  mint(20+js)=ksusy1+22
1428  mint(23-js)=ksusy1+25
1429 
1430  ELSEIF(isub.EQ.222) THEN
1431 C...f + fbar -> ~chi01 + ~chi04
1432  IF(mint(15).LT.0) js=2
1433 C IF(PYR(0).GT.0.5D0) JS=2
1434  mint(20+js)=ksusy1+22
1435  mint(23-js)=ksusy1+35
1436 
1437  ELSEIF(isub.EQ.223) THEN
1438 C...f + fbar -> ~chi02 + ~chi03
1439  IF(mint(15).LT.0) js=2
1440 C IF(PYR(0).GT.0.5D0) JS=2
1441  mint(20+js)=ksusy1+23
1442  mint(23-js)=ksusy1+25
1443 
1444  ELSEIF(isub.EQ.224) THEN
1445 C...f + fbar -> ~chi02 + ~chi04
1446  IF(mint(15).LT.0) js=2
1447 C IF(PYR(0).GT.0.5D0) JS=2
1448  mint(20+js)=ksusy1+23
1449  mint(23-js)=ksusy1+35
1450 
1451  ELSEIF(isub.EQ.225) THEN
1452 C...f + fbar -> ~chi03 + ~chi04
1453  IF(mint(15).LT.0) js=2
1454 C IF(PYR(0).GT.0.5D0) JS=2
1455  mint(20+js)=ksusy1+25
1456  mint(23-js)=ksusy1+35
1457  ENDIF
1458 
1459  ELSEIF(isub.LE.236) THEN
1460  IF(isub.EQ.226) THEN
1461 C...f + fbar -> ~chi+-1 + ~chi-+1
1462 C...th=(p(q)-p(chi+))**2 or (p(qbar)-p(chi-))**2
1463  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
1464  mint(21)=isign(ksusy1+24,kch1)
1465  mint(22)=-mint(21)
1466 
1467  ELSEIF(isub.EQ.227) THEN
1468 C...f + fbar -> ~chi+-2 + ~chi-+2
1469  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
1470  mint(21)=isign(ksusy1+37,kch1)
1471  mint(22)=-mint(21)
1472 
1473  ELSEIF(isub.EQ.228) THEN
1474 C...f + fbar -> ~chi+-1 + ~chi-+2
1475 C...th=(p(q)-p(chi1+))**2 or th=(p(qbar)-p(chi1-))**2
1476 C...js=1 if pyr<.5, js=2 if pyr>.5
1477 C...if 15=q, 16=qbar and js=1, chi1+ + chi2-, th=(q-chi1+)**2
1478 C...if 15=qbar, 16=q and js=1, chi2- + chi1+, th=(q-chi1+)**2
1479 C...if 15=q, 16=qbar and js=2, chi1- + chi2+, th=(qbar-chi1-)**2
1480 C...if 15=qbar, 16=q and js=2, chi2+ + chi1-, th=(q-chi1-)**2
1481  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
1482  kch2=int(1-kch1)/2
1483  IF(mint(2).EQ.1) THEN
1484  mint(21)= isign(ksusy1+24,kch1)
1485  mint(22)= -isign(ksusy1+37,kch1)
1486 c IF(KCH2.EQ.0) JS=2
1487  ELSE
1488  mint(21)= isign(ksusy1+37,kch1)
1489  mint(22)= -isign(ksusy1+24,kch1)
1490  js=2
1491 c IF(KCH2.EQ.1) JS=2
1492  ENDIF
1493 
1494  ELSEIF(isub.EQ.229) THEN
1495 C...q + qbar' -> ~chi01 + ~chi+-1
1496 C...th=(p(u)-p(chi+))**2 or (p(ubar)-p(chi-))**2
1497  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
1498  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
1499 C...CHECK THIS
1500  IF(mod(mint(15),2).EQ.0) js=2
1501  mint(20+js)=ksusy1+22
1502  mint(23-js)=isign(ksusy1+24,kch1+kch2)
1503 
1504  ELSEIF(isub.EQ.230) THEN
1505 C...q + qbar' -> ~chi02 + ~chi+-1
1506  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
1507  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
1508  IF(mod(mint(15),2).EQ.0) js=2
1509  mint(20+js)=ksusy1+23
1510  mint(23-js)=isign(ksusy1+24,kch1+kch2)
1511 
1512  ELSEIF(isub.EQ.231) THEN
1513 C...q + qbar' -> ~chi03 + ~chi+-1
1514  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
1515  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
1516  IF(mod(mint(15),2).EQ.0) js=2
1517  mint(20+js)=ksusy1+25
1518  mint(23-js)=isign(ksusy1+24,kch1+kch2)
1519 
1520  ELSEIF(isub.EQ.232) THEN
1521 C...q + qbar' -> ~chi04 + ~chi+-1
1522  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
1523  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
1524  IF(mod(mint(15),2).EQ.0) js=2
1525  mint(20+js)=ksusy1+35
1526  mint(23-js)=isign(ksusy1+24,kch1+kch2)
1527 
1528  ELSEIF(isub.EQ.233) THEN
1529 C...q + qbar' -> ~chi01 + ~chi+-2
1530  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
1531  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
1532  IF(mod(mint(15),2).EQ.0) js=2
1533  mint(20+js)=ksusy1+22
1534  mint(23-js)=isign(ksusy1+37,kch1+kch2)
1535 
1536  ELSEIF(isub.EQ.234) THEN
1537 C...q + qbar' -> ~chi02 + ~chi+-2
1538  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
1539  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
1540  IF(mod(mint(15),2).EQ.0) js=2
1541  mint(20+js)=ksusy1+23
1542  mint(23-js)=isign(ksusy1+37,kch1+kch2)
1543 
1544  ELSEIF(isub.EQ.235) THEN
1545 C...q + qbar' -> ~chi03 + ~chi+-2
1546  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
1547  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
1548  IF(mod(mint(15),2).EQ.0) js=2
1549  mint(20+js)=ksusy1+25
1550  mint(23-js)=isign(ksusy1+37,kch1+kch2)
1551 
1552  ELSEIF(isub.EQ.236) THEN
1553 C...q + qbar' -> ~chi04 + ~chi+-2
1554  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
1555  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
1556  IF(mod(mint(15),2).EQ.0) js=2
1557  mint(20+js)=ksusy1+35
1558  mint(23-js)=isign(ksusy1+37,kch1+kch2)
1559  ENDIF
1560 
1561  ELSEIF(isub.LE.245) THEN
1562  IF(isub.EQ.237) THEN
1563 C...q + qbar -> ~chi01 + ~g
1564 C...th arbitrary
1565  IF(pyr(0).GT.0.5d0) js=2
1566  mint(20+js)=ksusy1+21
1567  mint(23-js)=ksusy1+22
1568  kcc=17+js
1569 
1570  ELSEIF(isub.EQ.238) THEN
1571 C...q + qbar -> ~chi02 + ~g
1572 C...th arbitrary
1573  IF(pyr(0).GT.0.5d0) js=2
1574  mint(20+js)=ksusy1+21
1575  mint(23-js)=ksusy1+23
1576  kcc=17+js
1577 
1578  ELSEIF(isub.EQ.239) THEN
1579 C...q + qbar -> ~chi03 + ~g
1580 C...th arbitrary
1581  IF(pyr(0).GT.0.5d0) js=2
1582  mint(20+js)=ksusy1+21
1583  mint(23-js)=ksusy1+25
1584  kcc=17+js
1585 
1586  ELSEIF(isub.EQ.240) THEN
1587 C...q + qbar -> ~chi04 + ~g
1588 C...th arbitrary
1589  IF(pyr(0).GT.0.5d0) js=2
1590  mint(20+js)=ksusy1+21
1591  mint(23-js)=ksusy1+35
1592  kcc=17+js
1593 
1594  ELSEIF(isub.EQ.241) THEN
1595 C...q + qbar' -> ~chi+-1 + ~g
1596 C...if 15=u, 16=dbar, then (kch1+kch2)>0, js=1, chi+
1597 C...if 15=d, 16=ubar, then (kch1+kch2)<0, js=2, chi-
1598 C...if 15=ubar, 16=d, then (kch1+kch2)<0, js=1, chi-
1599 C...if 15=dbar, 16=u, then (kch1+kch2)>0, js=2, chi+
1600 C...th=(p(q)-p(chi+))**2 or (p(qbar')-p(chi-))**2
1601  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
1602  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
1603  js=1
1604  IF(mint(15)*(kch1+kch2).GT.0) js=2
1605  mint(20+js)=ksusy1+21
1606  mint(23-js)=isign(ksusy1+24,kch1+kch2)
1607  kcc=17+js
1608 
1609  ELSEIF(isub.EQ.242) THEN
1610 C...q + qbar' -> ~chi+-2 + ~g
1611 C...if 15=u, 16=dbar, then (kch1+kch2)>0, js=1, chi+
1612 C...if 15=d, 16=ubar, then (kch1+kch2)<0, js=2, chi-
1613 C...if 15=ubar, 16=d, then (kch1+kch2)<0, js=1, chi-
1614 C...if 15=dbar, 16=u, then (kch1+kch2)>0, js=2, chi+
1615 C...th=(p(q)-p(chi+))**2 or (p(qbar')-p(chi-))**2
1616  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
1617  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
1618  js=1
1619  IF(mint(15)*(kch1+kch2).GT.0) js=2
1620  mint(20+js)=ksusy1+21
1621  mint(23-js)=isign(ksusy1+37,kch1+kch2)
1622  kcc=17+js
1623 
1624  ELSEIF(isub.EQ.243) THEN
1625 C...q + qbar -> ~g + ~g ; th arbitrary
1626  mint(21)=ksusy1+21
1627  mint(22)=ksusy1+21
1628  kcc=mint(2)+4
1629 
1630  ELSEIF(isub.EQ.244) THEN
1631 C...g + g -> ~g + ~g ; th arbitrary
1632  kcc=mint(2)+12
1633  kcs=(-1)**int(1.5d0+pyr(0))
1634  mint(21)=ksusy1+21
1635  mint(22)=ksusy1+21
1636  ENDIF
1637 
1638  ELSEIF(isub.LE.260) THEN
1639  IF(isub.EQ.246) THEN
1640 C...qj + g -> ~qj_L + ~chi01
1641  IF(mint(15).EQ.21) js=2
1642  i=mint(14+js)
1643  ia=iabs(i)
1644  mint(20+js)=isign(ksusy1+ia,i)
1645  mint(23-js)=ksusy1+22
1646  kcc=15+js
1647  kcs=isign(1,mint(14+js))
1648 
1649  ELSEIF(isub.EQ.247) THEN
1650 C...qj + g -> ~qj_R + ~chi01
1651  IF(mint(15).EQ.21) js=2
1652  i=mint(14+js)
1653  ia=iabs(i)
1654  mint(20+js)=isign(ksusy2+ia,i)
1655  mint(23-js)=ksusy1+22
1656  kcc=15+js
1657  kcs=isign(1,mint(14+js))
1658 
1659  ELSEIF(isub.EQ.248) THEN
1660 C...qj + g -> ~qj_L + ~chi02
1661  IF(mint(15).EQ.21) js=2
1662  i=mint(14+js)
1663  ia=iabs(i)
1664  mint(20+js)=isign(ksusy1+ia,i)
1665  mint(23-js)=ksusy1+23
1666  kcc=15+js
1667  kcs=isign(1,mint(14+js))
1668 
1669  ELSEIF(isub.EQ.249) THEN
1670 C...qj + g -> ~qj_R + ~chi02
1671  IF(mint(15).EQ.21) js=2
1672  i=mint(14+js)
1673  ia=iabs(i)
1674  mint(20+js)=isign(ksusy2+ia,i)
1675  mint(23-js)=ksusy1+23
1676  kcc=15+js
1677  kcs=isign(1,mint(14+js))
1678 
1679  ELSEIF(isub.EQ.250) THEN
1680 C...qj + g -> ~qj_L + ~chi03
1681  IF(mint(15).EQ.21) js=2
1682  i=mint(14+js)
1683  ia=iabs(i)
1684  mint(20+js)=isign(ksusy1+ia,i)
1685  mint(23-js)=ksusy1+25
1686  kcc=15+js
1687  kcs=isign(1,mint(14+js))
1688 
1689  ELSEIF(isub.EQ.251) THEN
1690 C...qj + g -> ~qj_R + ~chi03
1691  IF(mint(15).EQ.21) js=2
1692  i=mint(14+js)
1693  ia=iabs(i)
1694  mint(20+js)=isign(ksusy2+ia,i)
1695  mint(23-js)=ksusy1+25
1696  kcc=15+js
1697  kcs=isign(1,mint(14+js))
1698 
1699  ELSEIF(isub.EQ.252) THEN
1700 C...qj + g -> ~qj_L + ~chi04
1701  IF(mint(15).EQ.21) js=2
1702  i=mint(14+js)
1703  ia=iabs(i)
1704  mint(20+js)=isign(ksusy1+ia,i)
1705  mint(23-js)=ksusy1+35
1706  kcc=15+js
1707  kcs=isign(1,mint(14+js))
1708 
1709  ELSEIF(isub.EQ.253) THEN
1710 C...qj + g -> ~qj_R + ~chi04
1711  IF(mint(15).EQ.21) js=2
1712  i=mint(14+js)
1713  ia=iabs(i)
1714  mint(20+js)=isign(ksusy2+ia,i)
1715  mint(23-js)=ksusy1+35
1716  kcc=15+js
1717  kcs=isign(1,mint(14+js))
1718 
1719  ELSEIF(isub.EQ.254) THEN
1720 C...qj + g -> ~qk_L + ~chi+-1
1721  IF(mint(15).EQ.21) js=2
1722  i=mint(14+js)
1723  ia=iabs(i)
1724  mint(23-js)=isign(ksusy1+24,kchg(ia,1)*i)
1725  ib=-ia+int((ia+1)/2)*4-1
1726  mint(20+js)=isign(ksusy1+ib,i)
1727  kcc=15+js
1728  kcs=isign(1,mint(14+js))
1729 
1730  ELSEIF(isub.EQ.255) THEN
1731 C...qj + g -> ~qk_L + ~chi+-1
1732  IF(mint(15).EQ.21) js=2
1733  i=mint(14+js)
1734  ia=iabs(i)
1735  mint(23-js)=isign(ksusy1+24,kchg(ia,1)*i)
1736  ib=-ia+int((ia+1)/2)*4-1
1737  mint(20+js)=isign(ksusy2+ib,i)
1738  kcc=15+js
1739  kcs=isign(1,mint(14+js))
1740 
1741  ELSEIF(isub.EQ.256) THEN
1742 C...qj + g -> ~qk_L + ~chi+-2
1743  IF(mint(15).EQ.21) js=2
1744  i=mint(14+js)
1745  ia=iabs(i)
1746  ib=-ia+int((ia+1)/2)*4-1
1747  mint(20+js)=isign(ksusy1+ib,i)
1748  mint(23-js)=isign(ksusy1+37,kchg(ia,1)*i)
1749  kcc=15+js
1750  kcs=isign(1,mint(14+js))
1751 
1752  ELSEIF(isub.EQ.257) THEN
1753 C...qj + g -> ~qk_R + ~chi+-2
1754  IF(mint(15).EQ.21) js=2
1755  i=mint(14+js)
1756  ia=iabs(i)
1757  ib=-ia+int((ia+1)/2)*4-1
1758  mint(20+js)=isign(ksusy2+ib,i)
1759  mint(23-js)=isign(ksusy1+37,kchg(ia,1)*i)
1760  kcc=15+js
1761  kcs=isign(1,mint(14+js))
1762 
1763  ELSEIF(isub.EQ.258) THEN
1764 C...qj + g -> ~qj_L + ~g
1765  IF(mint(15).EQ.21) js=2
1766  i=mint(14+js)
1767  ia=iabs(i)
1768  mint(20+js)=isign(ksusy1+ia,i)
1769  mint(23-js)=ksusy1+21
1770  kcc=mint(2)+6
1771  IF(js.EQ.2) kcc=kcc+2
1772  kcs=isign(1,i)
1773 
1774  ELSEIF(isub.EQ.259) THEN
1775 C...qj + g -> ~qj_R + ~g
1776  IF(mint(15).EQ.21) js=2
1777  i=mint(14+js)
1778  ia=iabs(i)
1779  mint(20+js)=isign(ksusy2+ia,i)
1780  mint(23-js)=ksusy1+21
1781  kcc=mint(2)+6
1782  IF(js.EQ.2) kcc=kcc+2
1783  kcs=isign(1,i)
1784  ENDIF
1785 
1786  ELSEIF(isub.LE.270) THEN
1787  IF(isub.EQ.261) THEN
1788 C...f + fbar -> ~t_1 + ~t_1bar; th = (p(q)-p(sq))**2
1789  isgn=1
1790  IF(mint(43).EQ.1.AND.pyr(0).GT.0.5d0) isgn=-1
1791  mint(21)=isgn*isign(kfpr(isub,1),kcs)
1792  mint(22)=-mint(21)
1793 C...Correct color combination
1794  IF(mint(43).EQ.4) kcc=4
1795 
1796  ELSEIF(isub.EQ.262) THEN
1797 C...f + fbar -> ~t_2 + ~t_2bar; th = (p(q)-p(sq))**2
1798  isgn=1
1799  IF(mint(43).EQ.1.AND.pyr(0).GT.0.5d0) isgn=-1
1800  mint(21)=isgn*isign(kfpr(isub,1),kcs)
1801  mint(22)=-mint(21)
1802 C...Correct color combination
1803  IF(mint(43).EQ.4) kcc=4
1804 
1805  ELSEIF(isub.EQ.263) THEN
1806 C...f + fbar -> ~t_1 + ~t_2bar; th = (p(q)-p(sq))**2
1807  IF((kcs.GT.0.AND.mint(2).EQ.1).OR.
1808  & (kcs.LT.0.AND.mint(2).EQ.2)) THEN
1809  mint(21)=isign(kfpr(isub,1),kcs)
1810  mint(22)=-isign(kfpr(isub,2),kcs)
1811  ELSE
1812  js=2
1813  mint(21)=isign(kfpr(isub,2),kcs)
1814  mint(22)=-isign(kfpr(isub,1),kcs)
1815  ENDIF
1816 C...Correct color combination
1817  IF(mint(43).EQ.4) kcc=4
1818 
1819  ELSEIF(isub.EQ.264) THEN
1820 C...g + g -> ~t_1 + ~t_1bar; th arbitrary
1821  kcs=(-1)**int(1.5d0+pyr(0))
1822  mint(21)=isign(kfpr(isub,1),kcs)
1823  mint(22)=-mint(21)
1824  kcc=mint(2)+10
1825 
1826  ELSEIF(isub.EQ.265) THEN
1827 C...g + g -> ~t_2 + ~t_2bar; th arbitrary
1828  kcs=(-1)**int(1.5d0+pyr(0))
1829  mint(21)=isign(kfpr(isub,1),kcs)
1830  mint(22)=-mint(21)
1831  kcc=mint(2)+10
1832  ENDIF
1833 
1834  ELSEIF(isub.LE.296) THEN
1835  IF(isub.EQ.271.OR.isub.EQ.281.OR.isub.EQ.291) THEN
1836 C...qi + qj -> ~qi_L + ~qj_L
1837  kcc=mint(2)
1838  IF(mint(15)*mint(16).LT.0) kcc=kcc+2
1839  mint(21)=isign(ksusy1+iabs(mint(15)),mint(15))
1840  mint(22)=isign(ksusy1+iabs(mint(16)),mint(16))
1841 
1842  ELSEIF(isub.EQ.272.OR.isub.EQ.282.OR.isub.EQ.292) THEN
1843 C...qi + qj -> ~qi_R + ~qj_R
1844  kcc=mint(2)
1845  IF(mint(15)*mint(16).LT.0) kcc=kcc+2
1846  mint(21)=isign(ksusy2+iabs(mint(15)),mint(15))
1847  mint(22)=isign(ksusy2+iabs(mint(16)),mint(16))
1848 
1849  ELSEIF(isub.EQ.273.OR.isub.EQ.283.OR.isub.EQ.293) THEN
1850 C...qi + qj -> ~qi_L + ~qj_R
1851  mint(21)=isign(kfpr(isub,1),mint(15))
1852  mint(22)=isign(kfpr(isub,2),mint(16))
1853  kcc=mint(2)
1854  IF(mint(15)*mint(16).LT.0) kcc=kcc+2
1855 
1856  ELSEIF(isub.EQ.274.OR.isub.EQ.284) THEN
1857 C...qi + qjbar -> ~qi_L + ~qj_Lbar; th = (p(f)-p(sf'))**2
1858  mint(21)=isign(ksusy1+iabs(mint(15)),mint(15))
1859  mint(22)=isign(ksusy1+iabs(mint(16)),mint(16))
1860  kcc=mint(2)
1861  IF(mint(15)*mint(16).LT.0) kcc=kcc+2
1862 
1863  ELSEIF(isub.EQ.275.OR.isub.EQ.285) THEN
1864 C...qi + qjbar -> ~qi_R + ~qj_Rbar ; th = (p(f)-p(sf'))**2
1865  mint(21)=isign(ksusy2+iabs(mint(15)),mint(15))
1866  mint(22)=isign(ksusy2+iabs(mint(16)),mint(16))
1867  kcc=mint(2)
1868  IF(mint(15)*mint(16).LT.0) kcc=kcc+2
1869 
1870  ELSEIF(isub.EQ.276.OR.isub.EQ.286.OR.isub.EQ.296) THEN
1871 C...qi + qjbar -> ~qi_L + ~qj_Rbar ; th = (p(f)-p(sf'))**2
1872  mint(21)=isign(kfpr(isub,1),mint(15))
1873  mint(22)=isign(kfpr(isub,2),mint(16))
1874  kcc=mint(2)
1875  IF(mint(15)*mint(16).LT.0) kcc=kcc+2
1876 
1877  ELSEIF(isub.EQ.277.OR.isub.EQ.287) THEN
1878 C...f + fbar -> ~qi_L + ~qi_Lbar ; th = (p(q)-p(sq))**2
1879  isgn=1
1880  IF(mint(43).EQ.1.AND.pyr(0).GT.0.5d0) isgn=-1
1881  mint(21)=isgn*isign(kfpr(isub,1),kcs)
1882  mint(22)=-mint(21)
1883  IF(mint(43).EQ.4) kcc=4
1884 
1885  ELSEIF(isub.EQ.278.OR.isub.EQ.288) THEN
1886 C...f + fbar -> ~qi_R + ~qi_Rbar; th = (p(q)-p(sq))**2
1887  isgn=1
1888  IF(mint(43).EQ.1.AND.pyr(0).GT.0.5d0) isgn=-1
1889  mint(21)=isgn*isign(kfpr(isub,1),kcs)
1890  mint(22)=-mint(21)
1891  IF(mint(43).EQ.4) kcc=4
1892 
1893  ELSEIF(isub.EQ.279.OR.isub.EQ.289) THEN
1894 C...g + g -> ~qi_L + ~qi_Lbar ; th arbitrary
1895 C...pure LL + RR
1896  kcs=(-1)**int(1.5d0+pyr(0))
1897  mint(21)=isign(kfpr(isub,1),kcs)
1898  mint(22)=-mint(21)
1899  kcc=mint(2)+10
1900 
1901  ELSEIF(isub.EQ.280.OR.isub.EQ.290) THEN
1902 C...g + g -> ~qi_R + ~qi_Rbar ; th arbitrary
1903  kcs=(-1)**int(1.5d0+pyr(0))
1904  mint(21)=isign(kfpr(isub,1),kcs)
1905  mint(22)=-mint(21)
1906  kcc=mint(2)+10
1907 
1908  ELSEIF(isub.EQ.294) THEN
1909 C...qj + g -> ~qj_L + ~g
1910  IF(mint(15).EQ.21) js=2
1911  i=mint(14+js)
1912  ia=iabs(i)
1913  mint(20+js)=isign(ksusy1+ia,i)
1914  mint(23-js)=ksusy1+21
1915  kcc=mint(2)+6
1916  IF(js.EQ.2) kcc=kcc+2
1917  kcs=isign(1,i)
1918 
1919  ELSEIF(isub.EQ.295) THEN
1920 C...qj + g -> ~qj_R + ~g
1921  IF(mint(15).EQ.21) js=2
1922  i=mint(14+js)
1923  ia=iabs(i)
1924  mint(20+js)=isign(ksusy2+ia,i)
1925  mint(23-js)=ksusy1+21
1926  kcc=mint(2)+6
1927  IF(js.EQ.2) kcc=kcc+2
1928  kcs=isign(1,i)
1929  ENDIF
1930 
1931  ELSEIF(isub.LE.340) THEN
1932 
1933  IF(isub.EQ.297.OR.isub.EQ.298) THEN
1934 C...q + qbar' -> H+ + H0
1935  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
1936  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
1937  IF(mint(15)*(kch1+kch2).GT.0) js=2
1938  mint(20+js)=isign(37,kch1+kch2)
1939  mint(23-js)=kfpr(isub,2)
1940  ELSEIF(isub.EQ.299.OR.isub.EQ.300) THEN
1941 C...f + fbar -> A0 + H0; th arbitrary
1942  IF(pyr(0).GT.0.5d0) js=2
1943  mint(20+js)=kfpr(isub,1)
1944  mint(23-js)=kfpr(isub,2)
1945  ELSEIF(isub.EQ.301) THEN
1946 C...f + fbar -> H+ H-
1947  mint(21)=isign(kfpr(isub,1),kcs)
1948  mint(22)=-mint(21)
1949  ENDIF
1950 CMRENNA--
1951 
1952  ELSEIF(isub.LE.360) THEN
1953 
1954  IF(isub.EQ.341.OR.isub.EQ.342) THEN
1955 C...l + l -> H_L++/--, H_R++/--
1956  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
1957  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
1958  kfres=isign(kfpr(isub,1),kch1+kch2)
1959 
1960  ELSEIF(isub.GE.343.AND.isub.LE.348) THEN
1961 C...l + gamma -> l' + H++/--; th=(p(l)-p(H))**2
1962  IF(mint(15).EQ.22) js=2
1963  mint(20+js)=isign(kfpr(isub,1),-mint(14+js))
1964  mint(23-js)=isign(kfpr(isub,2),-mint(14+js))
1965  kcc=22
1966 
1967  ELSEIF(isub.EQ.349.OR.isub.EQ.350) THEN
1968 C...f + fbar -> H++ + H--; th = (p(f)-p(H--))**2
1969  mint(21)=-isign(kfpr(isub,1),mint(15))
1970  mint(22)=-mint(21)
1971 
1972  ELSEIF(isub.EQ.351.OR.isub.EQ.352) THEN
1973 C...f + f' -> f" + f"' + H++/-- (W+/- + W+/- -> H++/--
1974 C...as inner process).
1975  DO 450 jt=1,2
1976  i=mint(14+jt)
1977  ia=iabs(i)
1978  IF(ia.LE.10) THEN
1979  rvckm=vint(180+i)*pyr(0)
1980  DO 440 j=1,mstp(1)
1981  ib=2*j-1+mod(ia,2)
1982  ipm=(5-isign(1,i))/2
1983  idc=j+mdcy(ia,2)+2
1984  IF(mdme(idc,1).NE.1.AND.mdme(idc,1).NE.ipm) goto 440
1985  mint(20+jt)=isign(ib,i)
1986  rvckm=rvckm-vckm((ia+1)/2,(ib+1)/2)
1987  IF(rvckm.LE.0d0) goto 450
1988  440 CONTINUE
1989  ELSE
1990  ib=2*((ia+1)/2)-1+mod(ia,2)
1991  mint(20+jt)=isign(ib,i)
1992  ENDIF
1993  450 CONTINUE
1994  kcc=22
1995  kfres=isign(kfpr(isub,1),mint(15))
1996  IF(mod(mint(15),2).EQ.1) kfres=-kfres
1997 
1998  ELSEIF(isub.EQ.353) THEN
1999 C...f + fbar -> Z_R0
2000  kfres=kfpr(isub,1)
2001 
2002  ELSEIF(isub.EQ.354) THEN
2003 C...f + fbar' -> W+/-
2004  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
2005  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
2006  kfres=isign(kfpr(isub,1),kch1+kch2)
2007 
2008  ENDIF
2009 
2010  ELSEIF(isub.LE.380) THEN
2011 
2012  IF(isub.LE.363.OR.isub.EQ.368) THEN
2013 C...f + fbar -> charged+ charged- technicolor
2014  ksw=(-1)**int(1.5d0+pyr(0))
2015  mint(21)=isign(kfpr(isub,1),ksw)
2016  mint(22)=-isign(kfpr(isub,2),ksw)
2017 
2018  ELSEIF(isub.LE.367) THEN
2019 C...f + fbar -> neutral neutral technicolor
2020  mint(21)=kfpr(isub,1)
2021  mint(22)=kfpr(isub,2)
2022 
2023  ELSEIF(isub.EQ.374.OR.isub.EQ.375) THEN
2024 C...f + fbar' -> neutral charged technicolor
2025  in=1
2026  ic=2
2027  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
2028  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
2029  IF(mint(15)*(kch1+kch2).LT.0) js=2
2030  mint(23-js)=isign(kfpr(isub,ic),kch1+kch2)
2031  mint(20+js)=kfpr(isub,in)
2032 
2033  ELSEIF(isub.GE.370.AND.isub.LE.377) THEN
2034 C...f + fbar' -> charged neutral technicolor
2035  in=2
2036  ic=1
2037  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
2038  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
2039  IF(mint(15)*(kch1+kch2).GT.0) js=2
2040  mint(20+js)=isign(kfpr(isub,ic),kch1+kch2)
2041  mint(23-js)=kfpr(isub,in)
2042  ENDIF
2043 
2044  ELSEIF(isub.LE.400) THEN
2045  IF(isub.EQ.381) THEN
2046 C...f + f' -> f + f' (g exchange); th = (p(f)-p(f))**2, TC extensions
2047  kcc=mint(2)
2048  IF(mint(15)*mint(16).LT.0) kcc=kcc+2
2049 
2050  ELSEIF(isub.EQ.382) THEN
2051 C...f + fbar -> f' + fbar'; th = (p(f)-p(f'))**2, TC extensions
2052  mint(21)=isign(kflf,mint(15))
2053  mint(22)=-mint(21)
2054  kcc=4
2055 
2056  ELSEIF(isub.EQ.383) THEN
2057 C...f + fbar -> g + g; th arbitrary, TC extensions
2058  mint(21)=21
2059  mint(22)=21
2060  kcc=mint(2)+4
2061 
2062  ELSEIF(isub.EQ.384) THEN
2063 C...f + g -> f + g; th = (p(f)-p(f))**2, TC extensions
2064  IF(mint(15).EQ.21) js=2
2065  kcc=mint(2)+6
2066  IF(mint(15).EQ.21) kcc=kcc+2
2067  IF(mint(15).NE.21) kcs=isign(1,mint(15))
2068  IF(mint(16).NE.21) kcs=isign(1,mint(16))
2069 
2070  ELSEIF(isub.EQ.385) THEN
2071 C...g + g -> f + fbar; th arbitrary, TC extensions
2072  kcs=(-1)**int(1.5d0+pyr(0))
2073  mint(21)=isign(kflf,kcs)
2074  mint(22)=-mint(21)
2075  kcc=mint(2)+10
2076 
2077  ELSEIF(isub.EQ.386) THEN
2078 C...g + g -> g + g; th arbitrary, TC extensions
2079  kcc=mint(2)+12
2080  kcs=(-1)**int(1.5d0+pyr(0))
2081 
2082  ELSEIF(isub.EQ.387) THEN
2083 C...q + qbar -> Q + Qbar; th = (p(q)-p(Q))**2, TC extensions
2084  mint(21)=isign(mint(55),mint(15))
2085  mint(22)=-mint(21)
2086  kcc=4
2087 
2088  ELSEIF(isub.EQ.388) THEN
2089 C...g + g -> Q + Qbar; th arbitrary, TC extensions
2090  kcs=(-1)**int(1.5d0+pyr(0))
2091  mint(21)=isign(mint(55),kcs)
2092  mint(22)=-mint(21)
2093  kcc=mint(2)+10
2094 
2095  ELSEIF(isub.EQ.391) THEN
2096 C...f + fbar -> G*.
2097  kfres=kfpr(isub,1)
2098 
2099  ELSEIF(isub.EQ.392) THEN
2100 C...g + g -> G*.
2101  kcc=21
2102  kfres=kfpr(isub,1)
2103 
2104  ELSEIF(isub.EQ.393) THEN
2105 C...q + qbar -> g + G*; th arbitrary.
2106  IF(pyr(0).GT.0.5d0) js=2
2107  mint(20+js)=kfpr(isub,1)
2108  mint(23-js)=kfpr(isub,2)
2109  kcc=17+js
2110 
2111  ELSEIF(isub.EQ.394) THEN
2112 C...q + g -> q + G*; th = (p(f) - p(f))**2
2113  IF(mint(15).EQ.21) js=2
2114  mint(23-js)=kfpr(isub,2)
2115  kcc=15+js
2116  kcs=isign(1,mint(14+js))
2117 
2118  ELSEIF(isub.EQ.395) THEN
2119 C...g + g -> G* + g; th arbitrary.
2120  IF(pyr(0).GT.0.5d0) js=2
2121  mint(23-js)=kfpr(isub,2)
2122  kcc=22+js
2123  ENDIF
2124 
2125  ELSEIF(isub.LE.420) THEN
2126  IF(isub.EQ.401) THEN
2127 C...g + g -> t + b + H+/-
2128  kcs=(-1)**int(1.5d0+pyr(0))
2129  mint(21)=isign(kfpr(isubsv,2),kcs)
2130  mint(22)=isign(5,-kcs)
2131  kcc=11+int(0.5d0+pyr(0))
2132  kfres=isign(kfhigg,-kcs)
2133 
2134  ELSEIF(isub.EQ.402) THEN
2135 C...q + qbar -> t + b + H+/-
2136  kfl=(-1)**int(1.5d0+pyr(0))
2137  mint(21)=isign(int(6.+.5*kfl),kcs)
2138  mint(22)=isign(int(6.-.5*kfl),-kcs)
2139  kcc=4
2140  kfres=isign(kfhigg,-kfl*kcs)
2141  ENDIF
2142 
2143 C...QUARKONIA+++
2144 C...Additional code by Stefan Wolf
2145  ELSEIF(isub.LE.430) THEN
2146  IF(isub.GE.421.AND.isub.LE.424) THEN
2147 C...g + g -> QQ~[n] + g
2148 C...MINT(21), MINT(22) copied from ISUB.EQ.86-89
2149 C...[g + g -> (J/Psi, chi_0c, chi_1c or chi_2c) + g]
2150 C...KCC and KCS copied from ISUB.EQ.86-89 (for ISUB.EQ.421)
2151 C...[g + g -> (J/Psi, chi_0c, chi_1c or chi_2c) + g]
2152 C...or from ISUB.EQ.68 (for ISUB.NE.421)
2153 C...[g + g -> g + g; th arbitrary]
2154  mint(21)=kfpr(isubsv,1)
2155  mint(22)=kfpr(isubsv,2)
2156  IF(isub.EQ.421) THEN
2157  kcc=24
2158  kcs=(-1)**int(1.5d0+pyr(0))
2159  ELSE
2160  kcc=mint(2)+12
2161  kcs=(-1)**int(1.5d0+pyr(0))
2162  ENDIF
2163 
2164  ELSEIF(isub.GE.425.AND.isub.LE.427) THEN
2165 C...q + g -> q + QQ~[n]
2166 C...MINT(21), MINT(22) "copied" from ISUB.EQ.112
2167 C...[f + g -> f + h0; th = (p(f)-p(f))**2; (q + g -> q + h0 only)]
2168 C...KCC copied from ISUB.EQ.28
2169 C...[f + g -> f + g; th = (p(f)-p(f))**2; (q + g -> q + g only)]
2170  IF(mint(15).EQ.21) js=2
2171  mint(23-js)=kfpr(isubsv,2)
2172  kcc=mint(2)+6
2173  IF(mint(15).EQ.21) kcc=kcc+2
2174  IF(mint(15).NE.21) kcs=isign(1,mint(15))
2175  IF(mint(16).NE.21) kcs=isign(1,mint(16))
2176 
2177  ELSEIF(isub.GE.428.AND.isub.LE.430) THEN
2178 C...q + q~ -> g + QQ~[n]
2179 C...MINT(21), MINT(22) "copied" from ISUB.EQ.111
2180 C...[f + fbar -> g + h0; th arbitrary; (q + qbar -> g + h0 only)]
2181 C...KCC copied from ISUB.EQ.13
2182 C...[f + fbar -> g + g; th arbitrary; (q + qbar -> g + g only)]
2183  IF(pyr(0).GT.0.5) js=2
2184  mint(20+js)=21
2185  mint(23-js)=kfpr(isubsv,2)
2186  kcc=mint(2)+4
2187  ENDIF
2188 
2189  ELSEIF(isub.LE.440) THEN
2190  IF(isub.GE.431.AND.isub.LE.433) THEN
2191 C...g + g -> QQ~[n] + g
2192 C...MINT(21), MINT(22) copied from ISUB.EQ.86-89
2193 C...[g + g -> (J/Psi, chi_0c, chi_1c or chi_2c) + g]
2194 C...KCC and KCS copied from ISUB.EQ.86-89
2195 C...[g + g -> (J/Psi, chi_0c, chi_1c or chi_2c) + g]
2196  mint(21)=kfpr(isubsv,1)
2197  mint(22)=kfpr(isubsv,2)
2198  kcc=24
2199  kcs=(-1)**int(1.5d0+pyr(0))
2200 
2201  ELSEIF(isub.GE.434.AND.isub.LE.436) THEN
2202 C...q + g -> q + QQ~[n]
2203 C...MINT(21), MINT(22) "copied" from ISUB.EQ.112
2204 C...[f + g -> f + h0; th = (p(f)-p(f))**2; (q + g -> q + h0 only)]
2205 C...KCC and KCS copied from ISUB.EQ.112
2206 C...[f + g -> f + h0; th = (p(f)-p(f))**2; (q + g -> q + h0 only)]
2207  IF(mint(15).EQ.21) js=2
2208  mint(23-js)=kfpr(isubsv,2)
2209  kcc=15+js
2210  kcs=isign(1,mint(14+js))
2211 
2212  ELSEIF(isub.GE.437.AND.isub.LE.439) THEN
2213 C...q + q~ -> g + QQ~[n]
2214 C...MINT(21), MINT(22) "copied" from ISUB.EQ.111
2215 C...[f + fbar -> g + h0; th arbitrary; (q + qbar -> g + h0 only)]
2216 C...KCC copied from ISUB.EQ.111
2217 C...[f + fbar -> g + h0; th arbitrary; (q + qbar -> g + h0 only)]
2218  IF(pyr(0).GT.0.5) js=2
2219  mint(20+js)=21
2220  mint(23-js)=kfpr(isubsv,2)
2221  kcc=17+js
2222  ENDIF
2223 C...QUARKONIA---
2224 
2225  ENDIF
2226 
2227  IF(iset(isub).EQ.11) THEN
2228 C...Store documentation for user-defined processes
2229  bezup=(pup(3,1)+pup(3,2))/(pup(4,1)+pup(4,2))
2230  kuppo(1)=mint(83)+5
2231  kuppo(2)=mint(83)+6
2232  i=mint(83)+6
2233  DO 470 iup=3,nup
2234  kuppo(iup)=0
2235  IF(mstp(128).GE.2.AND.mothup(1,iup).GE.3) THEN
2236  idoc=idoc-1
2237  mint(4)=mint(4)-1
2238  goto 470
2239  ENDIF
2240  i=i+1
2241  kuppo(iup)=i
2242  k(i,1)=21
2243  k(i,2)=idup(iup)
2244  IF(idup(iup).EQ.0) k(i,2)=90
2245  k(i,3)=0
2246  IF(mothup(1,iup).GE.3) k(i,3)=kuppo(mothup(1,iup))
2247  k(i,4)=0
2248  k(i,5)=0
2249  DO 460 j=1,5
2250  p(i,j)=pup(j,iup)
2251  460 CONTINUE
2252  v(i,5)=vtimup(iup)
2253  470 CONTINUE
2254  CALL pyrobo(mint(83)+7,mint(83)+4+nup,0d0,vint(24),0d0,0d0,
2255  & -bezup)
2256 
2257 C...Store final state partons for user-defined processes
2258  n=ipu2
2259  DO 490 iup=3,nup
2260  n=n+1
2261  k(n,1)=1
2262  IF(istup(iup).EQ.2.OR.istup(iup).EQ.3) k(n,1)=11
2263  k(n,2)=idup(iup)
2264  IF(idup(iup).EQ.0) k(n,2)=90
2265  IF(mstp(128).LE.0.OR.mothup(1,iup).EQ.0) THEN
2266  k(n,3)=kuppo(iup)
2267  ELSE
2268  k(n,3)=mint(84)+mothup(1,iup)
2269  ENDIF
2270  k(n,4)=0
2271  k(n,5)=0
2272 C...Search for daughters of intermediate colourless particles.
2273  IF(k(n,1).EQ.11.AND.kchg(pycomp(k(n,2)),2).EQ.0) THEN
2274  DO 475 iupdau=iup+1,nup
2275  IF(mothup(1,iupdau).EQ.iup.AND.k(n,4).EQ.0) k(n,4)=
2276  & n+iupdau-iup
2277  IF(mothup(1,iupdau).EQ.iup) k(n,5)=n+iupdau-iup
2278  475 CONTINUE
2279  ENDIF
2280  DO 480 j=1,5
2281  p(n,j)=pup(j,iup)
2282  480 CONTINUE
2283  v(n,5)=vtimup(iup)
2284  490 CONTINUE
2285  CALL pyrobo(ipu3,n,0d0,vint(24),0d0,0d0,-bezup)
2286 
2287 C...Arrange colour flow for user-defined processes
2288  nlbl=0
2289  DO 540 iup1=1,nup
2290  i1=mint(84)+iup1
2291  IF(kchg(pycomp(k(i1,2)),2).EQ.0) goto 540
2292  IF(k(i1,1).EQ.1) k(i1,1)=3
2293  IF(k(i1,1).EQ.11) k(i1,1)=14
2294 C...Find a not yet considered colour/anticolour line.
2295  DO 530 isde1=1,2
2296  IF(icolup(isde1,iup1).EQ.0) goto 530
2297  nmat=0
2298  DO 500 ilbl=1,nlbl
2299  IF(icolup(isde1,iup1).EQ.ilab(ilbl)) nmat=1
2300  500 CONTINUE
2301  IF(nmat.EQ.0) THEN
2302  nlbl=nlbl+1
2303  ilab(nlbl)=icolup(isde1,iup1)
2304 C...Find all others belonging to same line.
2305  i3=i1
2306  i4=0
2307  DO 520 iup2=iup1+1,nup
2308  i2=mint(84)+iup2
2309  DO 510 isde2=1,2
2310  IF(icolup(isde2,iup2).EQ.icolup(isde1,iup1)) THEN
2311  IF(isde2.EQ.isde1) THEN
2312  k(i3,3+isde2)=k(i3,3+isde2)+i2
2313  k(i2,3+isde2)=k(i2,3+isde2)+mstu(5)*i3
2314  i3=i2
2315  ELSEIF(i4.NE.0) THEN
2316  k(i4,3+isde2)=k(i4,3+isde2)+i2
2317  k(i2,3+isde2)=k(i2,3+isde2)+mstu(5)*i4
2318  i4=i2
2319  ELSEIF(iup2.LE.2) THEN
2320  k(i1,3+isde1)=k(i1,3+isde1)+i2
2321  k(i2,3+isde2)=k(i2,3+isde2)+i1
2322  i4=i2
2323  ELSE
2324  k(i1,3+isde1)=k(i1,3+isde1)+mstu(5)*i2
2325  k(i2,3+isde2)=k(i2,3+isde2)+mstu(5)*i1
2326  i4=i2
2327  ENDIF
2328  ENDIF
2329  510 CONTINUE
2330  520 CONTINUE
2331  ENDIF
2332  530 CONTINUE
2333  540 CONTINUE
2334 
2335  ELSEIF(idoc.EQ.7) THEN
2336 C...Resonance not decaying; store kinematics
2337  i=mint(83)+7
2338  k(ipu3,1)=1
2339  k(ipu3,2)=kfres
2340  k(ipu3,3)=i
2341  p(ipu3,4)=shuser
2342  p(ipu3,5)=shuser
2343  k(i,1)=21
2344  k(i,2)=kfres
2345  p(i,4)=shuser
2346  p(i,5)=shuser
2347  n=ipu3
2348  mint(21)=kfres
2349  mint(22)=0
2350 
2351 C...Special cases: colour flow in coloured resonances
2352  kcres=pycomp(kfres)
2353  IF(kchg(kcres,2).NE.0) THEN
2354  k(ipu3,1)=3
2355  DO 550 j=1,2
2356  jc=j
2357  IF(kcs.EQ.-1) jc=3-j
2358  IF(icol(kcc,1,jc).NE.0.AND.k(ipu1,1).EQ.14) k(ipu1,j+3)=
2359  & mint(84)+icol(kcc,1,jc)
2360  IF(icol(kcc,2,jc).NE.0.AND.k(ipu2,1).EQ.14) k(ipu2,j+3)=
2361  & mint(84)+icol(kcc,2,jc)
2362  IF(icol(kcc,3,jc).NE.0.AND.k(ipu3,1).EQ.3) k(ipu3,j+3)=
2363  & mstu(5)*(mint(84)+icol(kcc,3,jc))
2364  550 CONTINUE
2365  ELSE
2366  k(ipu1,4)=ipu2
2367  k(ipu1,5)=ipu2
2368  k(ipu2,4)=ipu1
2369  k(ipu2,5)=ipu1
2370  ENDIF
2371 
2372  ELSEIF(idoc.EQ.8) THEN
2373 C...2 -> 2 processes: store outgoing partons in their CM-frame
2374  DO 560 jt=1,2
2375  i=mint(84)+2+jt
2376  kca=pycomp(mint(20+jt))
2377  k(i,1)=1
2378  IF(kchg(kca,2).NE.0) k(i,1)=3
2379  k(i,2)=mint(20+jt)
2380  k(i,3)=mint(83)+idoc+jt-2
2381  kfaa=iabs(k(i,2))
2382  IF(kfpr(isubsv,1+mod(js+jt,2)).NE.0) THEN
2383  p(i,5)=sqrt(vint(63+mod(js+jt,2)))
2384  ELSE
2385  p(i,5)=pymass(k(i,2))
2386  ENDIF
2387  IF((kfaa.EQ.6.OR.kfaa.EQ.7.OR.kfaa.EQ.8).AND.
2388  & p(i,5).LT.parp(42)) p(i,5)=pymass(k(i,2))
2389  560 CONTINUE
2390  IF(p(ipu3,5)+p(ipu4,5).GE.shr) THEN
2391  kfa1=iabs(mint(21))
2392  kfa2=iabs(mint(22))
2393  IF((kfa1.GT.3.AND.kfa1.NE.21).OR.(kfa2.GT.3.AND.kfa2.NE.21))
2394  & THEN
2395  mint(51)=1
2396  RETURN
2397  ENDIF
2398  p(ipu3,5)=0d0
2399  p(ipu4,5)=0d0
2400  ENDIF
2401  p(ipu3,4)=0.5d0*(shr+(p(ipu3,5)**2-p(ipu4,5)**2)/shr)
2402  p(ipu3,3)=sqrt(max(0d0,p(ipu3,4)**2-p(ipu3,5)**2))
2403  p(ipu4,4)=shr-p(ipu3,4)
2404  p(ipu4,3)=-p(ipu3,3)
2405  n=ipu4
2406  mint(7)=mint(83)+7
2407  mint(8)=mint(83)+8
2408 
2409 C...Rotate outgoing partons using cos(theta)=(th-uh)/lam(sh,sqm3,sqm4)
2410  CALL pyrobo(ipu3,ipu4,acos(vint(23)),vint(24),0d0,0d0,0d0)
2411 
2412  ELSEIF(idoc.EQ.9) THEN
2413 C...2 -> 3 processes: store outgoing partons in their CM frame
2414  DO 570 jt=1,2
2415  i=mint(84)+2+jt
2416  kca=pycomp(mint(20+jt))
2417  k(i,1)=1
2418  IF(kchg(kca,2).NE.0) k(i,1)=3
2419  k(i,2)=mint(20+jt)
2420  k(i,3)=mint(83)+idoc+jt-3
2421  jta=jt
2422 C...t and b in opposide order in event list as compared to
2423 C...matrix element?
2424  IF(isub.EQ.402.AND.iabs(mint(21)).EQ.5) jta=3-jt
2425  IF(iabs(k(i,2)).LE.22) THEN
2426  p(i,5)=pymass(k(i,2))
2427  ELSE
2428  p(i,5)=sqrt(vint(63+mod(js+jta,2)))
2429  ENDIF
2430  pt=sqrt(max(0d0,vint(197+5*jta)-p(i,5)**2+vint(196+5*jta)**2))
2431  p(i,1)=pt*cos(vint(198+5*jta))
2432  p(i,2)=pt*sin(vint(198+5*jta))
2433  570 CONTINUE
2434  k(ipu5,1)=1
2435  k(ipu5,2)=kfres
2436  k(ipu5,3)=mint(83)+idoc
2437  p(ipu5,5)=shr
2438  p(ipu5,1)=-p(ipu3,1)-p(ipu4,1)
2439  p(ipu5,2)=-p(ipu3,2)-p(ipu4,2)
2440  pms1=p(ipu3,5)**2+p(ipu3,1)**2+p(ipu3,2)**2
2441  pms2=p(ipu4,5)**2+p(ipu4,1)**2+p(ipu4,2)**2
2442  pms3=p(ipu5,5)**2+p(ipu5,1)**2+p(ipu5,2)**2
2443  pmt3=sqrt(pms3)
2444  p(ipu5,3)=pmt3*sinh(vint(211))
2445  p(ipu5,4)=pmt3*cosh(vint(211))
2446  pms12=(shpr-p(ipu5,4))**2-p(ipu5,3)**2
2447  sql12=(pms12-pms1-pms2)**2-4d0*pms1*pms2
2448  IF(sql12.LE.0d0) THEN
2449  mint(51)=1
2450  RETURN
2451  ENDIF
2452  p(ipu3,3)=(-p(ipu5,3)*(pms12+pms1-pms2)+
2453  & vint(213)*(shpr-p(ipu5,4))*sqrt(sql12))/(2d0*pms12)
2454  p(ipu4,3)=-p(ipu3,3)-p(ipu5,3)
2455  IF(isub.EQ.402.AND.iabs(mint(21)).EQ.5) THEN
2456 C...t and b in opposide order in event list as compared to
2457 C...matrix element
2458  p(ipu4,3)=(-p(ipu5,3)*(pms12+pms2-pms1)+
2459  & vint(213)*(shpr-p(ipu5,4))*sqrt(sql12))/(2d0*pms12)
2460  p(ipu3,3)=-p(ipu4,3)-p(ipu5,3)
2461  END IF
2462  p(ipu3,4)=sqrt(pms1+p(ipu3,3)**2)
2463  p(ipu4,4)=sqrt(pms2+p(ipu4,3)**2)
2464  mint(23)=kfres
2465  n=ipu5
2466  mint(7)=mint(83)+7
2467  mint(8)=mint(83)+8
2468 
2469  ELSEIF(idoc.EQ.11) THEN
2470 C...Z0 + Z0 -> h0, W+ + W- -> h0: store Higgs and outgoing partons
2471  phi(1)=paru(2)*pyr(0)
2472  phi(2)=phi(1)-phir
2473  DO 580 jt=1,2
2474  i=mint(84)+2+jt
2475  k(i,1)=1
2476  IF(kchg(pycomp(mint(20+jt)),2).NE.0) k(i,1)=3
2477  k(i,2)=mint(20+jt)
2478  k(i,3)=mint(83)+idoc+jt-2
2479  p(i,5)=pymass(k(i,2))
2480  IF(0.5d0*shpr*z(jt).LE.p(i,5)) THEN
2481  mint(51)=1
2482  RETURN
2483  ENDIF
2484  pabs=sqrt(max(0d0,(0.5d0*shpr*z(jt))**2-p(i,5)**2))
2485  ptabs=pabs*sqrt(max(0d0,1d0-cthe(jt)**2))
2486  p(i,1)=ptabs*cos(phi(jt))
2487  p(i,2)=ptabs*sin(phi(jt))
2488  p(i,3)=pabs*cthe(jt)*(-1)**(jt+1)
2489  p(i,4)=0.5d0*shpr*z(jt)
2490  izw=mint(83)+6+jt
2491  k(izw,1)=21
2492  k(izw,2)=23
2493  IF(isub.EQ.8) k(izw,2)=isign(24,pychge(mint(14+jt)))
2494  k(izw,3)=izw-2
2495  p(izw,1)=-p(i,1)
2496  p(izw,2)=-p(i,2)
2497  p(izw,3)=(0.5d0*shpr-pabs*cthe(jt))*(-1)**(jt+1)
2498  p(izw,4)=0.5d0*shpr*(1d0-z(jt))
2499  p(izw,5)=-sqrt(max(0d0,p(izw,3)**2+ptabs**2-p(izw,4)**2))
2500  580 CONTINUE
2501  i=mint(83)+9
2502  k(ipu5,1)=1
2503  k(ipu5,2)=kfres
2504  k(ipu5,3)=i
2505  p(ipu5,5)=shr
2506  p(ipu5,1)=-p(ipu3,1)-p(ipu4,1)
2507  p(ipu5,2)=-p(ipu3,2)-p(ipu4,2)
2508  p(ipu5,3)=-p(ipu3,3)-p(ipu4,3)
2509  p(ipu5,4)=shpr-p(ipu3,4)-p(ipu4,4)
2510  k(i,1)=21
2511  k(i,2)=kfres
2512  DO 590 j=1,5
2513  p(i,j)=p(ipu5,j)
2514  590 CONTINUE
2515  n=ipu5
2516  mint(23)=kfres
2517 
2518  ELSEIF(idoc.EQ.12) THEN
2519 C...Z0 and W+/- scattering: store bosons and outgoing partons
2520  phi(1)=paru(2)*pyr(0)
2521  phi(2)=phi(1)-phir
2522  jtran=int(1.5d0+pyr(0))
2523  DO 600 jt=1,2
2524  i=mint(84)+2+jt
2525  k(i,1)=1
2526  IF(kchg(pycomp(mint(20+jt)),2).NE.0) k(i,1)=3
2527  k(i,2)=mint(20+jt)
2528  k(i,3)=mint(83)+idoc+jt-2
2529  p(i,5)=pymass(k(i,2))
2530  IF(0.5d0*shpr*z(jt).LE.p(i,5)) p(i,5)=0d0
2531  pabs=sqrt(max(0d0,(0.5d0*shpr*z(jt))**2-p(i,5)**2))
2532  ptabs=pabs*sqrt(max(0d0,1d0-cthe(jt)**2))
2533  p(i,1)=ptabs*cos(phi(jt))
2534  p(i,2)=ptabs*sin(phi(jt))
2535  p(i,3)=pabs*cthe(jt)*(-1)**(jt+1)
2536  p(i,4)=0.5d0*shpr*z(jt)
2537  izw=mint(83)+6+jt
2538  k(izw,1)=21
2539  IF(mint(14+jt).EQ.mint(20+jt)) THEN
2540  k(izw,2)=23
2541  ELSE
2542  k(izw,2)=isign(24,pychge(mint(14+jt))-pychge(mint(20+jt)))
2543  ENDIF
2544  k(izw,3)=izw-2
2545  p(izw,1)=-p(i,1)
2546  p(izw,2)=-p(i,2)
2547  p(izw,3)=(0.5d0*shpr-pabs*cthe(jt))*(-1)**(jt+1)
2548  p(izw,4)=0.5d0*shpr*(1d0-z(jt))
2549  p(izw,5)=-sqrt(max(0d0,p(izw,3)**2+ptabs**2-p(izw,4)**2))
2550  ipu=mint(84)+4+jt
2551  k(ipu,1)=3
2552  k(ipu,2)=kfpr(isub,jt)
2553  IF(isub.EQ.72.AND.jt.EQ.jtran) k(ipu,2)=-k(ipu,2)
2554  IF(isub.EQ.73.OR.isub.EQ.77) k(ipu,2)=k(izw,2)
2555  k(ipu,3)=mint(83)+8+jt
2556  IF(iabs(k(ipu,2)).LE.10.OR.k(ipu,2).EQ.21) THEN
2557  p(ipu,5)=pymass(k(ipu,2))
2558  ELSE
2559  p(ipu,5)=sqrt(vint(63+mod(js+jt,2)))
2560  ENDIF
2561  mint(22+jt)=k(ipu,2)
2562  600 CONTINUE
2563 C...Find rotation and boost for hard scattering subsystem
2564  i1=mint(83)+7
2565  i2=mint(83)+8
2566  bexcm=(p(i1,1)+p(i2,1))/(p(i1,4)+p(i2,4))
2567  beycm=(p(i1,2)+p(i2,2))/(p(i1,4)+p(i2,4))
2568  bezcm=(p(i1,3)+p(i2,3))/(p(i1,4)+p(i2,4))
2569  gamcm=(p(i1,4)+p(i2,4))/shr
2570  bepcm=bexcm*p(i1,1)+beycm*p(i1,2)+bezcm*p(i1,3)
2571  px=p(i1,1)+gamcm*(gamcm/(1d0+gamcm)*bepcm-p(i1,4))*bexcm
2572  py=p(i1,2)+gamcm*(gamcm/(1d0+gamcm)*bepcm-p(i1,4))*beycm
2573  pz=p(i1,3)+gamcm*(gamcm/(1d0+gamcm)*bepcm-p(i1,4))*bezcm
2574  thecm=pyangl(pz,sqrt(px**2+py**2))
2575  phicm=pyangl(px,py)
2576 C...Store hard scattering subsystem. Rotate and boost it
2577  sqlam=(sh-p(ipu5,5)**2-p(ipu6,5)**2)**2-4d0*p(ipu5,5)**2*
2578  & p(ipu6,5)**2
2579  pabs=sqrt(max(0d0,sqlam/(4d0*sh)))
2580  cthwz=vint(23)
2581  sthwz=sqrt(max(0d0,1d0-cthwz**2))
2582  phiwz=vint(24)-phicm
2583  p(ipu5,1)=pabs*sthwz*cos(phiwz)
2584  p(ipu5,2)=pabs*sthwz*sin(phiwz)
2585  p(ipu5,3)=pabs*cthwz
2586  p(ipu5,4)=sqrt(pabs**2+p(ipu5,5)**2)
2587  p(ipu6,1)=-p(ipu5,1)
2588  p(ipu6,2)=-p(ipu5,2)
2589  p(ipu6,3)=-p(ipu5,3)
2590  p(ipu6,4)=sqrt(pabs**2+p(ipu6,5)**2)
2591  CALL pyrobo(ipu5,ipu6,thecm,phicm,bexcm,beycm,bezcm)
2592  DO 620 jt=1,2
2593  i1=mint(83)+8+jt
2594  i2=mint(84)+4+jt
2595  k(i1,1)=21
2596  k(i1,2)=k(i2,2)
2597  DO 610 j=1,5
2598  p(i1,j)=p(i2,j)
2599  610 CONTINUE
2600  620 CONTINUE
2601  n=ipu6
2602  mint(7)=mint(83)+9
2603  mint(8)=mint(83)+10
2604  ENDIF
2605 
2606  IF(iset(isub).EQ.11) THEN
2607  ELSEIF(idoc.GE.8) THEN
2608 C...Store colour connection indices
2609  DO 630 j=1,2
2610  jc=j
2611  IF(kcs.EQ.-1) jc=3-j
2612  IF(icol(kcc,1,jc).NE.0.AND.k(ipu1,1).EQ.14) k(ipu1,j+3)=
2613  & k(ipu1,j+3)+mint(84)+icol(kcc,1,jc)
2614  IF(icol(kcc,2,jc).NE.0.AND.k(ipu2,1).EQ.14) k(ipu2,j+3)=
2615  & k(ipu2,j+3)+mint(84)+icol(kcc,2,jc)
2616  IF(icol(kcc,3,jc).NE.0.AND.k(ipu3,1).EQ.3) k(ipu3,j+3)=
2617  & mstu(5)*(mint(84)+icol(kcc,3,jc))
2618  IF(icol(kcc,4,jc).NE.0.AND.k(ipu4,1).EQ.3) k(ipu4,j+3)=
2619  & mstu(5)*(mint(84)+icol(kcc,4,jc))
2620  630 CONTINUE
2621 
2622 C...Copy outgoing partons to documentation lines
2623  imax=2
2624  IF(idoc.EQ.9) imax=3
2625  DO 650 i=1,imax
2626  i1=mint(83)+idoc-imax+i
2627  i2=mint(84)+2+i
2628  k(i1,1)=21
2629  k(i1,2)=k(i2,2)
2630  IF(idoc.LE.9) k(i1,3)=0
2631  IF(idoc.GE.11) k(i1,3)=mint(83)+2+i
2632  DO 640 j=1,5
2633  p(i1,j)=p(i2,j)
2634  640 CONTINUE
2635  650 CONTINUE
2636 
2637  ELSEIF(idoc.EQ.9) THEN
2638 C...Store colour connection indices
2639  DO 660 j=1,2
2640  jc=j
2641  IF(kcs.EQ.-1) jc=3-j
2642  IF(icol(kcc,1,jc).NE.0.AND.k(ipu1,1).EQ.14) k(ipu1,j+3)=
2643  & k(ipu1,j+3)+mint(84)+icol(kcc,1,jc)+
2644  & max(0,min(1,icol(kcc,1,jc)-2))
2645  IF(icol(kcc,2,jc).NE.0.AND.k(ipu2,1).EQ.14) k(ipu2,j+3)=
2646  & k(ipu2,j+3)+mint(84)+icol(kcc,2,jc)+
2647  & max(0,min(1,icol(kcc,2,jc)-2))
2648  IF(icol(kcc,3,jc).NE.0.AND.k(ipu4,1).EQ.3) k(ipu4,j+3)=
2649  & mstu(5)*(mint(84)+icol(kcc,3,jc))
2650  IF(icol(kcc,4,jc).NE.0.AND.k(ipu5,1).EQ.3) k(ipu5,j+3)=
2651  & mstu(5)*(mint(84)+icol(kcc,4,jc))
2652  660 CONTINUE
2653 
2654 C...Copy outgoing partons to documentation lines
2655  DO 680 i=1,3
2656  i1=mint(83)+idoc-3+i
2657  i2=mint(84)+2+i
2658  k(i1,1)=21
2659  k(i1,2)=k(i2,2)
2660  k(i1,3)=0
2661  DO 670 j=1,5
2662  p(i1,j)=p(i2,j)
2663  670 CONTINUE
2664  680 CONTINUE
2665  ENDIF
2666 
2667 C...Copy outgoing partons to list of allowed radiators.
2668  npart=0
2669  IF(mint(35).GE.2.AND.iset(isub).NE.0) THEN
2670  DO 690 i=mint(84)+3,n
2671  npart=npart+1
2672  ipart(npart)=i
2673  ptpart(npart)=sqrt(p(i,5)**2+p(i,1)**2+p(i,2)**2)
2674  690 CONTINUE
2675  ENDIF
2676 
2677 C...Low-pT events: remove gluons used for string drawing purposes
2678  IF(isub.EQ.95) THEN
2679  IF(mint(35).LE.1) THEN
2680  k(ipu3,1)=k(ipu3,1)+10
2681  k(ipu4,1)=k(ipu4,1)+10
2682  ENDIF
2683  DO 700 j=41,66
2684  vintsv(j)=vint(j)
2685  vint(j)=0d0
2686  700 CONTINUE
2687  DO 720 i=mint(83)+5,mint(83)+8
2688  DO 710 j=1,5
2689  p(i,j)=0d0
2690  710 CONTINUE
2691  720 CONTINUE
2692  ENDIF
2693 
2694  RETURN
2695  END