Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyhisigh.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyhisigh.f
1 
2 C***********************************************************************
3 
4  SUBROUTINE pyhisigh(NCHN,SIGS)
5 
6 C...Differential matrix elements for all included subprocesses.
7 C...Note that what is coded is (disregarding the COMFAC factor)
8 C...1) for 2 -> 1 processes: s-hat/pi*d(sigma-hat), where,
9 C...when d(sigma-hat) is given in the zero-width limit, the delta
10 C...function in tau is replaced by a Breit-Wigner:
11 C...1/pi*(s*m_res*Gamma_res)/((s*tau-m_res^2)^2+(m_res*Gamma_res)^2);
12 C...2) for 2 -> 2 processes: (s-hat)**2/pi*d(sigma-hat)/d(t-hat);
13 C...i.e., dimensionless quantities. COMFAC contains the factor
14 C...pi/s and the conversion factor from GeV^-2 to mb.
15  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
16  SAVE /ludat1/
17  common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
18  SAVE /ludat2/
19  common/ludat3/mdcy(500,3),mdme(2000,2),brat(2000),kfdp(2000,5)
20  SAVE /ludat3/
21  common/pyhisubs/msel,msub(200),kfin(2,-40:40),ckin(200)
22  SAVE /pyhisubs/
23  common/pyhipars/mstp(200),parp(200),msti(200),pari(200)
24  SAVE /pyhipars/
25  common/pyhiint1/mint(400),vint(400)
26  SAVE /pyhiint1/
27  common/pyhiint2/iset(200),kfpr(200,2),coef(200,20),icol(40,4,2)
28  SAVE /pyhiint2/
29  common/pyhiint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)
30  SAVE /pyhiint3/
31  common/pyhiint4/widp(21:40,0:40),wide(21:40,0:40),wids(21:40,3)
32  SAVE /pyhiint4/
33  common/pyhiint5/ngen(0:200,3),xsec(0:200,3)
34  SAVE /pyhiint5/
35  dimension x(2),xpq(-6:6),kfac(2,-40:40),wdtp(0:40),wdte(0:40,0:5)
36 
37 C...Reset number of channels and cross-section.
38  nchn=0
39  sigs=0.
40 
41 C...Read kinematical variables and limits.
42  isub=mint(1)
43  taumin=vint(11)
44  ystmin=vint(12)
45  ctnmin=vint(13)
46  ctpmin=vint(14)
47  xt2min=vint(15)
48  taupmn=vint(16)
49  tau=vint(21)
50  yst=vint(22)
51  cth=vint(23)
52  xt2=vint(25)
53  taup=vint(26)
54  taumax=vint(31)
55  ystmax=vint(32)
56  ctnmax=vint(33)
57  ctpmax=vint(34)
58  xt2max=vint(35)
59  taupmx=vint(36)
60 
61 C...Derive kinematical quantities.
62  IF(iset(isub).LE.2.OR.iset(isub).EQ.5) THEN
63  x(1)=sqrt(tau)*exp(yst)
64  x(2)=sqrt(tau)*exp(-yst)
65  ELSE
66  x(1)=sqrt(taup)*exp(yst)
67  x(2)=sqrt(taup)*exp(-yst)
68  ENDIF
69  IF(mint(43).EQ.4.AND.iset(isub).GE.1.AND.
70  &(x(1).GT.0.999.OR.x(2).GT.0.999)) RETURN
71  sh=tau*vint(2)
72  sqm3=vint(63)
73  sqm4=vint(64)
74  rm3=sqm3/sh
75  rm4=sqm4/sh
76  be34=sqrt((1.-rm3-rm4)**2-4.*rm3*rm4)
77  rpts=4.*vint(71)**2/sh
78  be34l=sqrt(max(0.,(1.-rm3-rm4)**2-4.*rm3*rm4-rpts))
79  rm34=2.*rm3*rm4
80  rsqm=1.+rm34
81  rthm=(4.*rm3*rm4+rpts)/(1.-rm3-rm4+be34l)
82  th=-0.5*sh*max(rthm,1.-rm3-rm4-be34*cth)
83  uh=-0.5*sh*max(rthm,1.-rm3-rm4+be34*cth)
84  sqpth=0.25*sh*be34**2*(1.-cth**2)
85  sh2=sh**2
86  th2=th**2
87  uh2=uh**2
88 
89 C...Choice of Q2 scale.
90  IF(iset(isub).EQ.1.OR.iset(isub).EQ.3) THEN
91  q2=sh
92  ELSEIF(mod(iset(isub),2).EQ.0.OR.iset(isub).EQ.5) THEN
93  IF(mstp(32).EQ.1) THEN
94  q2=2.*sh*th*uh/(sh**2+th**2+uh**2)
95  ELSEIF(mstp(32).EQ.2) THEN
96  q2=sqpth+0.5*(sqm3+sqm4)
97  ELSEIF(mstp(32).EQ.3) THEN
98  q2=min(-th,-uh)
99  ELSEIF(mstp(32).EQ.4) THEN
100  q2=sh
101  ENDIF
102  IF(iset(isub).EQ.5.AND.mstp(82).GE.2) q2=q2+parp(82)**2
103  ENDIF
104 
105 C...Store derived kinematical quantities.
106  vint(41)=x(1)
107  vint(42)=x(2)
108  vint(44)=sh
109  vint(43)=sqrt(sh)
110  vint(45)=th
111  vint(46)=uh
112  vint(48)=sqpth
113  vint(47)=sqrt(sqpth)
114  vint(50)=taup*vint(2)
115  vint(49)=sqrt(max(0.,vint(50)))
116  vint(52)=q2
117  vint(51)=sqrt(q2)
118 
119 C...Calculate parton structure functions.
120  IF(iset(isub).LE.0) goto 145
121  IF(mint(43).GE.2) THEN
122  q2sf=q2
123  IF(iset(isub).EQ.3.OR.iset(isub).EQ.4) THEN
124  q2sf=pmas(23,1)**2
125  IF(isub.EQ.8.OR.isub.EQ.76.OR.isub.EQ.77) q2sf=pmas(24,1)**2
126  ENDIF
127  DO 100 i=3-mint(41),mint(42)
128  xsf=x(i)
129  IF(iset(isub).EQ.5) xsf=x(i)/vint(142+i)
130  CALL pyhistfu(mint(10+i),xsf,q2sf,xpq,i)
131  DO 100 kfl=-6,6
132  100 xsfx(i,kfl)=xpq(kfl)
133  ENDIF
134 
135 C...Calculate alpha_strong and K-factor.
136  IF(mstp(33).NE.3) as=ulalps(q2)
137  fack=1.
138  faca=1.
139  IF(mstp(33).EQ.1) THEN
140  fack=parp(31)
141  ELSEIF(mstp(33).EQ.2) THEN
142  fack=parp(31)
143  faca=parp(32)/parp(31)
144  ELSEIF(mstp(33).EQ.3) THEN
145  q2as=parp(33)*q2
146  IF(iset(isub).EQ.5.AND.mstp(82).GE.2) q2as=q2as+
147  & paru(112)*parp(82)
148  as=ulalps(q2as)
149  ENDIF
150  radc=1.+as/paru(1)
151 
152 C...Set flags for allowed reacting partons/leptons.
153  DO 130 i=1,2
154  DO 110 j=-40,40
155  110 kfac(i,j)=0
156  IF(mint(40+i).EQ.1) THEN
157  kfac(i,mint(10+i))=1
158  ELSE
159  DO 120 j=-40,40
160  kfac(i,j)=kfin(i,j)
161  IF(abs(j).GT.mstp(54).AND.j.NE.21) kfac(i,j)=0
162  IF(abs(j).LE.6) THEN
163  IF(xsfx(i,j).LT.1.e-10) kfac(i,j)=0
164  ELSEIF(j.EQ.21) THEN
165  IF(xsfx(i,0).LT.1.e-10) kfac(i,21)=0
166  ENDIF
167  120 CONTINUE
168  ENDIF
169  130 CONTINUE
170 
171 C...Lower and upper limit for flavour loops.
172  min1=0
173  max1=0
174  min2=0
175  max2=0
176  DO 140 j=-20,20
177  IF(kfac(1,-j).EQ.1) min1=-j
178  IF(kfac(1,j).EQ.1) max1=j
179  IF(kfac(2,-j).EQ.1) min2=-j
180  IF(kfac(2,j).EQ.1) max2=j
181  140 CONTINUE
182  mina=min(min1,min2)
183  maxa=max(max1,max2)
184 
185 C...Common conversion factors (including Jacobian) for subprocesses.
186  sqmz=pmas(23,1)**2
187  gmmz=pmas(23,1)*pmas(23,2)
188  sqmw=pmas(24,1)**2
189  gmmw=pmas(24,1)*pmas(24,2)
190  sqmh=pmas(25,1)**2
191  gmmh=pmas(25,1)*pmas(25,2)
192  sqmzp=pmas(32,1)**2
193  gmmzp=pmas(32,1)*pmas(32,2)
194  sqmhc=pmas(37,1)**2
195  gmmhc=pmas(37,1)*pmas(37,2)
196  sqmr=pmas(40,1)**2
197  gmmr=pmas(40,1)*pmas(40,2)
198  aem=paru(101)
199  xw=paru(102)
200 
201 C...Phase space integral in tau and y*.
202  comfac=paru(1)*paru(5)/vint(2)
203  IF(mint(43).EQ.4) comfac=comfac*fack
204  IF((mint(43).GE.2.OR.iset(isub).EQ.3.OR.iset(isub).EQ.4).AND.
205  &iset(isub).NE.5) THEN
206  atau0=log(taumax/taumin)
207  atau1=(taumax-taumin)/(taumax*taumin)
208  h1=coef(isub,1)+(atau0/atau1)*coef(isub,2)/tau
209  IF(mint(72).GE.1) THEN
210  taur1=vint(73)
211  gamr1=vint(74)
212  atau2=log(taumax/taumin*(taumin+taur1)/(taumax+taur1))/taur1
213  atau3=(atan((taumax-taur1)/gamr1)-atan((taumin-taur1)/gamr1))/
214  & gamr1
215  h1=h1+(atau0/atau2)*coef(isub,3)/(tau+taur1)+
216  & (atau0/atau3)*coef(isub,4)*tau/((tau-taur1)**2+gamr1**2)
217  ENDIF
218  IF(mint(72).EQ.2) THEN
219  taur2=vint(75)
220  gamr2=vint(76)
221  atau4=log(taumax/taumin*(taumin+taur2)/(taumax+taur2))/taur2
222  atau5=(atan((taumax-taur2)/gamr2)-atan((taumin-taur2)/gamr2))/
223  & gamr2
224  h1=h1+(atau0/atau4)*coef(isub,5)/(tau+taur2)+
225  & (atau0/atau5)*coef(isub,6)*tau/((tau-taur2)**2+gamr2**2)
226  ENDIF
227  comfac=comfac*atau0/(tau*h1)
228  ENDIF
229  IF(mint(43).EQ.4.AND.iset(isub).NE.5) THEN
230  ayst0=ystmax-ystmin
231  ayst1=0.5*(ystmax-ystmin)**2
232  ayst2=ayst1
233  ayst3=2.*(atan(exp(ystmax))-atan(exp(ystmin)))
234  h2=(ayst0/ayst1)*coef(isub,7)*(yst-ystmin)+(ayst0/ayst2)*
235  & coef(isub,8)*(ystmax-yst)+(ayst0/ayst3)*coef(isub,9)/cosh(yst)
236  comfac=comfac*ayst0/h2
237  ENDIF
238 
239 C...2 -> 1 processes: reduction in angular part of phase space integral
240 C...for case of decaying resonance.
241  acth0=ctnmax-ctnmin+ctpmax-ctpmin
242  IF((iset(isub).EQ.1.OR.iset(isub).EQ.3).AND.
243  &mdcy(kfpr(isub,1),1).EQ.1) THEN
244  IF(kfpr(isub,1).EQ.25.OR.kfpr(isub,1).EQ.37) THEN
245  comfac=comfac*0.5*acth0
246  ELSE
247  comfac=comfac*0.125*(3.*acth0+ctnmax**3-ctnmin**3+
248  & ctpmax**3-ctpmin**3)
249  ENDIF
250 
251 C...2 -> 2 processes: angular part of phase space integral.
252  ELSEIF(iset(isub).EQ.2.OR.iset(isub).EQ.4) THEN
253  acth1=log((max(rm34,rsqm-ctnmin)*max(rm34,rsqm-ctpmin))/
254  & (max(rm34,rsqm-ctnmax)*max(rm34,rsqm-ctpmax)))
255  acth2=log((max(rm34,rsqm+ctnmax)*max(rm34,rsqm+ctpmax))/
256  & (max(rm34,rsqm+ctnmin)*max(rm34,rsqm+ctpmin)))
257  acth3=1./max(rm34,rsqm-ctnmax)-1./max(rm34,rsqm-ctnmin)+
258  & 1./max(rm34,rsqm-ctpmax)-1./max(rm34,rsqm-ctpmin)
259  acth4=1./max(rm34,rsqm+ctnmin)-1./max(rm34,rsqm+ctnmax)+
260  & 1./max(rm34,rsqm+ctpmin)-1./max(rm34,rsqm+ctpmax)
261  h3=coef(isub,10)+
262  & (acth0/acth1)*coef(isub,11)/max(rm34,rsqm-cth)+
263  & (acth0/acth2)*coef(isub,12)/max(rm34,rsqm+cth)+
264  & (acth0/acth3)*coef(isub,13)/max(rm34,rsqm-cth)**2+
265  & (acth0/acth4)*coef(isub,14)/max(rm34,rsqm+cth)**2
266  comfac=comfac*acth0*0.5*be34/h3
267  ENDIF
268 
269 C...2 -> 3, 4 processes: phace space integral in tau'.
270  IF(mint(43).GE.2.AND.(iset(isub).EQ.3.OR.iset(isub).EQ.4)) THEN
271  ataup0=log(taupmx/taupmn)
272  ataup1=((1.-tau/taupmx)**4-(1.-tau/taupmn)**4)/(4.*tau)
273  h4=coef(isub,15)+
274  & ataup0/ataup1*coef(isub,16)/taup*(1.-tau/taup)**3
275  IF(1.-tau/taup.GT.1.e-4) THEN
276  fzw=(1.+tau/taup)*log(taup/tau)-2.*(1.-tau/taup)
277  ELSE
278  fzw=1./6.*(1.-tau/taup)**3*tau/taup
279  ENDIF
280  comfac=comfac*ataup0*fzw/h4
281  ENDIF
282 
283 C...Phase space integral for low-pT and multiple interactions.
284  IF(iset(isub).EQ.5) THEN
285  comfac=paru(1)*paru(5)*fack*0.5*vint(2)/sh2
286  atau0=log(2.*(1.+sqrt(1.-xt2))/xt2-1.)
287  atau1=2.*atan(1./xt2-1.)/sqrt(xt2)
288  h1=coef(isub,1)+(atau0/atau1)*coef(isub,2)/sqrt(tau)
289  comfac=comfac*atau0/h1
290  ayst0=ystmax-ystmin
291  ayst1=0.5*(ystmax-ystmin)**2
292  ayst3=2.*(atan(exp(ystmax))-atan(exp(ystmin)))
293  h2=(ayst0/ayst1)*coef(isub,7)*(yst-ystmin)+(ayst0/ayst1)*
294  & coef(isub,8)*(ystmax-yst)+(ayst0/ayst3)*coef(isub,9)/cosh(yst)
295  comfac=comfac*ayst0/h2
296  IF(mstp(82).LE.1) comfac=comfac*xt2**2*(1./vint(149)-1.)
297 C...For MSTP(82)>=2 an additional factor (xT2/(xT2+VINT(149))**2 is
298 C...introduced to make cross-section finite for xT2 -> 0.
299  IF(mstp(82).GE.2) comfac=comfac*xt2**2/(vint(149)*
300  & (1.+vint(149)))
301  ENDIF
302 
303 C...A: 2 -> 1, tree diagrams.
304 
305  145 IF(isub.LE.10) THEN
306  IF(isub.EQ.1) THEN
307 C...f + fb -> gamma*/Z0.
308  mint(61)=2
309  CALL pyhiwidt(23,sqrt(sh),wdtp,wdte)
310  facz=comfac*aem**2*4./3.
311  DO 150 i=mina,maxa
312  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 150
313  ei=kchg(iabs(i),1)/3.
314  ai=sign(1.,ei)
315  vi=ai-4.*ei*xw
316  facf=1.
317  IF(iabs(i).LE.10) facf=faca/3.
318  nchn=nchn+1
319  isig(nchn,1)=i
320  isig(nchn,2)=-i
321  isig(nchn,3)=1
322  sigh(nchn)=facf*facz*(ei**2*vint(111)+ei*vi/(8.*xw*(1.-xw))*
323  & sh*(sh-sqmz)/((sh-sqmz)**2+gmmz**2)*vint(112)+(vi**2+ai**2)/
324  & (16.*xw*(1.-xw))**2*sh2/((sh-sqmz)**2+gmmz**2)*vint(114))
325  150 CONTINUE
326 
327  ELSEIF(isub.EQ.2) THEN
328 C...f + fb' -> W+/-.
329  CALL pyhiwidt(24,sqrt(sh),wdtp,wdte)
330  facw=comfac*(aem/xw)**2*1./24*sh2/((sh-sqmw)**2+gmmw**2)
331  DO 170 i=min1,max1
332  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 170
333  ia=iabs(i)
334  DO 160 j=min2,max2
335  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 160
336  ja=iabs(j)
337  IF(i*j.GT.0.OR.mod(ia+ja,2).EQ.0) goto 160
338  IF((ia.LE.10.AND.ja.GT.10).OR.(ia.GT.10.AND.ja.LE.10)) goto 160
339  kchw=(kchg(ia,1)*isign(1,i)+kchg(ja,1)*isign(1,j))/3
340  facf=1.
341  IF(ia.LE.10) facf=vckm((ia+1)/2,(ja+1)/2)*faca/3.
342  nchn=nchn+1
343  isig(nchn,1)=i
344  isig(nchn,2)=j
345  isig(nchn,3)=1
346  sigh(nchn)=facf*facw*(wdte(0,1)+wdte(0,(5-kchw)/2)+wdte(0,4))
347  160 CONTINUE
348  170 CONTINUE
349 
350  ELSEIF(isub.EQ.3) THEN
351 C...f + fb -> H0.
352  CALL pyhiwidt(25,sqrt(sh),wdtp,wdte)
353  fach=comfac*(aem/xw)**2*1./48.*(sh/sqmw)**2*
354  & sh2/((sh-sqmh)**2+gmmh**2)*(wdte(0,1)+wdte(0,2)+wdte(0,4))
355  DO 180 i=mina,maxa
356  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 180
357  rmq=pmas(iabs(i),1)**2/sh
358  nchn=nchn+1
359  isig(nchn,1)=i
360  isig(nchn,2)=-i
361  isig(nchn,3)=1
362  sigh(nchn)=fach*rmq*sqrt(max(0.,1.-4.*rmq))
363  180 CONTINUE
364 
365  ELSEIF(isub.EQ.4) THEN
366 C...gamma + W+/- -> W+/-.
367 
368  ELSEIF(isub.EQ.5) THEN
369 C...Z0 + Z0 -> H0.
370  CALL pyhiwidt(25,sqrt(sh),wdtp,wdte)
371  fach=comfac*1./(128.*paru(1)**2*16.*(1.-xw)**3)*(aem/xw)**4*
372  & (sh/sqmw)**2*sh2/((sh-sqmh)**2+gmmh**2)*
373  & (wdte(0,1)+wdte(0,2)+wdte(0,4))
374  DO 200 i=min1,max1
375  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 200
376  DO 190 j=min2,max2
377  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 190
378  ei=kchg(iabs(i),1)/3.
379  ai=sign(1.,ei)
380  vi=ai-4.*ei*xw
381  ej=kchg(iabs(j),1)/3.
382  aj=sign(1.,ej)
383  vj=aj-4.*ej*xw
384  nchn=nchn+1
385  isig(nchn,1)=i
386  isig(nchn,2)=j
387  isig(nchn,3)=1
388  sigh(nchn)=fach*(vi**2+ai**2)*(vj**2+aj**2)
389  190 CONTINUE
390  200 CONTINUE
391 
392  ELSEIF(isub.EQ.6) THEN
393 C...Z0 + W+/- -> W+/-.
394 
395  ELSEIF(isub.EQ.7) THEN
396 C...W+ + W- -> Z0.
397 
398  ELSEIF(isub.EQ.8) THEN
399 C...W+ + W- -> H0.
400  CALL pyhiwidt(25,sqrt(sh),wdtp,wdte)
401  fach=comfac*1./(128*paru(1)**2)*(aem/xw)**4*(sh/sqmw)**2*
402  & sh2/((sh-sqmh)**2+gmmh**2)*(wdte(0,1)+wdte(0,2)+wdte(0,4))
403  DO 220 i=min1,max1
404  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 220
405  ei=sign(1.,float(i))*kchg(iabs(i),1)
406  DO 210 j=min2,max2
407  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 210
408  ej=sign(1.,float(j))*kchg(iabs(j),1)
409  IF(ei*ej.GT.0.) goto 210
410  nchn=nchn+1
411  isig(nchn,1)=i
412  isig(nchn,2)=j
413  isig(nchn,3)=1
414  sigh(nchn)=fach*vint(180+i)*vint(180+j)
415  210 CONTINUE
416  220 CONTINUE
417  ENDIF
418 
419 C...B: 2 -> 2, tree diagrams.
420 
421  ELSEIF(isub.LE.20) THEN
422  IF(isub.EQ.11) THEN
423 C...f + f' -> f + f'.
424  facqq1=comfac*as**2*4./9.*(sh2+uh2)/th2
425  facqqb=comfac*as**2*4./9.*((sh2+uh2)/th2*faca-
426  & mstp(34)*2./3.*uh2/(sh*th))
427  facqq2=comfac*as**2*4./9.*((sh2+th2)/uh2-
428  & mstp(34)*2./3.*sh2/(th*uh))
429  DO 240 i=min1,max1
430  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 240
431  DO 230 j=min2,max2
432  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 230
433  nchn=nchn+1
434  isig(nchn,1)=i
435  isig(nchn,2)=j
436  isig(nchn,3)=1
437  sigh(nchn)=facqq1
438  IF(i.EQ.-j) sigh(nchn)=facqqb
439  IF(i.EQ.j) THEN
440  sigh(nchn)=0.5*sigh(nchn)
441  nchn=nchn+1
442  isig(nchn,1)=i
443  isig(nchn,2)=j
444  isig(nchn,3)=2
445  sigh(nchn)=0.5*facqq2
446  ENDIF
447  230 CONTINUE
448  240 CONTINUE
449 
450  ELSEIF(isub.EQ.12) THEN
451 C...f + fb -> f' + fb' (q + qb -> q' + qb' only).
452  CALL pyhiwidt(21,sqrt(sh),wdtp,wdte)
453  facqqb=comfac*as**2*4./9.*(th2+uh2)/sh2*(wdte(0,1)+wdte(0,2)+
454  & wdte(0,3)+wdte(0,4))
455  DO 250 i=mina,maxa
456  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 250
457  nchn=nchn+1
458  isig(nchn,1)=i
459  isig(nchn,2)=-i
460  isig(nchn,3)=1
461  sigh(nchn)=facqqb
462  250 CONTINUE
463 
464  ELSEIF(isub.EQ.13) THEN
465 C...f + fb -> g + g (q + qb -> g + g only).
466  facgg1=comfac*as**2*32./27.*(uh/th-(2.+mstp(34)*1./4.)*uh2/sh2)
467  facgg2=comfac*as**2*32./27.*(th/uh-(2.+mstp(34)*1./4.)*th2/sh2)
468  DO 260 i=mina,maxa
469  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 260
470  nchn=nchn+1
471  isig(nchn,1)=i
472  isig(nchn,2)=-i
473  isig(nchn,3)=1
474  sigh(nchn)=0.5*facgg1
475  nchn=nchn+1
476  isig(nchn,1)=i
477  isig(nchn,2)=-i
478  isig(nchn,3)=2
479  sigh(nchn)=0.5*facgg2
480  260 CONTINUE
481 
482  ELSEIF(isub.EQ.14) THEN
483 C...f + fb -> g + gamma (q + qb -> g + gamma only).
484  facgg=comfac*as*aem*8./9.*(th2+uh2)/(th*uh)
485  DO 270 i=mina,maxa
486  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 270
487  ei=kchg(iabs(i),1)/3.
488  nchn=nchn+1
489  isig(nchn,1)=i
490  isig(nchn,2)=-i
491  isig(nchn,3)=1
492  sigh(nchn)=facgg*ei**2
493  270 CONTINUE
494 
495  ELSEIF(isub.EQ.15) THEN
496 C...f + fb -> g + Z0 (q + qb -> g + Z0 only).
497  faczg=comfac*as*aem/(xw*(1.-xw))*1./18.*
498  & (th2+uh2+2.*sqm4*sh)/(th*uh)
499  faczg=faczg*wids(23,2)
500  DO 280 i=mina,maxa
501  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 280
502  ei=kchg(iabs(i),1)/3.
503  ai=sign(1.,ei)
504  vi=ai-4.*ei*xw
505  nchn=nchn+1
506  isig(nchn,1)=i
507  isig(nchn,2)=-i
508  isig(nchn,3)=1
509  sigh(nchn)=faczg*(vi**2+ai**2)
510  280 CONTINUE
511 
512  ELSEIF(isub.EQ.16) THEN
513 C...f + fb' -> g + W+/- (q + qb' -> g + W+/- only).
514  facwg=comfac*as*aem/xw*2./9.*(th2+uh2+2.*sqm4*sh)/(th*uh)
515  DO 300 i=min1,max1
516  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 300
517  ia=iabs(i)
518  DO 290 j=min2,max2
519  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 290
520  ja=iabs(j)
521  IF(i*j.GT.0.OR.mod(ia+ja,2).EQ.0) goto 290
522  kchw=(kchg(ia,1)*isign(1,i)+kchg(ja,1)*isign(1,j))/3
523  fckm=1.
524  IF(mint(43).EQ.4) fckm=vckm((ia+1)/2,(ja+1)/2)
525  nchn=nchn+1
526  isig(nchn,1)=i
527  isig(nchn,2)=j
528  isig(nchn,3)=1
529  sigh(nchn)=facwg*fckm*wids(24,(5-kchw)/2)
530  290 CONTINUE
531  300 CONTINUE
532 
533  ELSEIF(isub.EQ.17) THEN
534 C...f + fb -> g + H0 (q + qb -> g + H0 only).
535 
536  ELSEIF(isub.EQ.18) THEN
537 C...f + fb -> gamma + gamma.
538  facgg=comfac*faca*aem**2*1./3.*(th2+uh2)/(th*uh)
539  DO 310 i=mina,maxa
540  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 310
541  ei=kchg(iabs(i),1)/3.
542  nchn=nchn+1
543  isig(nchn,1)=i
544  isig(nchn,2)=-i
545  isig(nchn,3)=1
546  sigh(nchn)=facgg*ei**4
547  310 CONTINUE
548 
549  ELSEIF(isub.EQ.19) THEN
550 C...f + fb -> gamma + Z0.
551  facgz=comfac*faca*aem**2/(xw*(1.-xw))*1./24.*
552  & (th2+uh2+2.*sqm4*sh)/(th*uh)
553  facgz=facgz*wids(23,2)
554  DO 320 i=mina,maxa
555  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 320
556  ei=kchg(iabs(i),1)/3.
557  ai=sign(1.,ei)
558  vi=ai-4.*ei*xw
559  nchn=nchn+1
560  isig(nchn,1)=i
561  isig(nchn,2)=-i
562  isig(nchn,3)=1
563  sigh(nchn)=facgz*ei**2*(vi**2+ai**2)
564  320 CONTINUE
565 
566  ELSEIF(isub.EQ.20) THEN
567 C...f + fb' -> gamma + W+/-.
568  facgw=comfac*faca*aem**2/xw*1./6.*
569  & ((2.*uh-th)/(3.*(sh-sqm4)))**2*(th2+uh2+2.*sqm4*sh)/(th*uh)
570  DO 340 i=min1,max1
571  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 340
572  ia=iabs(i)
573  DO 330 j=min2,max2
574  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 330
575  ja=iabs(j)
576  IF(i*j.GT.0.OR.mod(ia+ja,2).EQ.0) goto 330
577  kchw=(kchg(ia,1)*isign(1,i)+kchg(ja,1)*isign(1,j))/3
578  fckm=1.
579  IF(mint(43).EQ.4) fckm=vckm((ia+1)/2,(ja+1)/2)
580  nchn=nchn+1
581  isig(nchn,1)=i
582  isig(nchn,2)=j
583  isig(nchn,3)=1
584  sigh(nchn)=facgw*fckm*wids(24,(5-kchw)/2)
585  330 CONTINUE
586  340 CONTINUE
587  ENDIF
588 
589  ELSEIF(isub.LE.30) THEN
590  IF(isub.EQ.21) THEN
591 C...f + fb -> gamma + H0.
592 
593  ELSEIF(isub.EQ.22) THEN
594 C...f + fb -> Z0 + Z0.
595  faczz=comfac*faca*(aem/(xw*(1.-xw)))**2*1./768.*
596  & (uh/th+th/uh+2.*(sqm3+sqm4)*sh/(th*uh)-
597  & sqm3*sqm4*(1./th2+1./uh2))
598  faczz=faczz*wids(23,1)
599  DO 350 i=mina,maxa
600  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 350
601  ei=kchg(iabs(i),1)/3.
602  ai=sign(1.,ei)
603  vi=ai-4.*ei*xw
604  nchn=nchn+1
605  isig(nchn,1)=i
606  isig(nchn,2)=-i
607  isig(nchn,3)=1
608  sigh(nchn)=faczz*(vi**4+6.*vi**2*ai**2+ai**4)
609  350 CONTINUE
610 
611  ELSEIF(isub.EQ.23) THEN
612 C...f + fb' -> Z0 + W+/-.
613  faczw=comfac*faca*(aem/xw)**2*1./6.
614  faczw=faczw*wids(23,2)
615  thuh=max(th*uh-sqm3*sqm4,sh*ckin(3)**2)
616  DO 370 i=min1,max1
617  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 370
618  ia=iabs(i)
619  DO 360 j=min2,max2
620  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 360
621  ja=iabs(j)
622  IF(i*j.GT.0.OR.mod(ia+ja,2).EQ.0) goto 360
623  kchw=(kchg(ia,1)*isign(1,i)+kchg(ja,1)*isign(1,j))/3
624  ei=kchg(ia,1)/3.
625  ai=sign(1.,ei)
626  vi=ai-4.*ei*xw
627  ej=kchg(ja,1)/3.
628  aj=sign(1.,ej)
629  vj=aj-4.*ej*xw
630  IF(vi+ai.GT.0) THEN
631  visav=vi
632  aisav=ai
633  vi=vj
634  ai=aj
635  vj=visav
636  aj=aisav
637  ENDIF
638  fckm=1.
639  IF(mint(43).EQ.4) fckm=vckm((ia+1)/2,(ja+1)/2)
640  nchn=nchn+1
641  isig(nchn,1)=i
642  isig(nchn,2)=j
643  isig(nchn,3)=1
644  sigh(nchn)=faczw*fckm*(1./(sh-sqmw)**2*
645  & ((9.-8.*xw)/4.*thuh+(8.*xw-6.)/4.*sh*(sqm3+sqm4))+
646  & (thuh-sh*(sqm3+sqm4))/(2.*(sh-sqmw))*((vj+aj)/th-(vi+ai)/uh)+
647  & thuh/(16.*(1.-xw))*((vj+aj)**2/th2+(vi+ai)**2/uh2)+
648  & sh*(sqm3+sqm4)/(8.*(1.-xw))*(vi+ai)*(vj+aj)/(th*uh))*
649  & wids(24,(5-kchw)/2)
650  360 CONTINUE
651  370 CONTINUE
652 
653  ELSEIF(isub.EQ.24) THEN
654 C...f + fb -> Z0 + H0.
655  thuh=max(th*uh-sqm3*sqm4,sh*ckin(3)**2)
656  fachz=comfac*faca*(aem/(xw*(1.-xw)))**2*1./96.*
657  & (thuh+2.*sh*sqmz)/(sh-sqmz)**2
658  fachz=fachz*wids(23,2)*wids(25,2)
659  DO 380 i=mina,maxa
660  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 380
661  ei=kchg(iabs(i),1)/3.
662  ai=sign(1.,ei)
663  vi=ai-4.*ei*xw
664  nchn=nchn+1
665  isig(nchn,1)=i
666  isig(nchn,2)=-i
667  isig(nchn,3)=1
668  sigh(nchn)=fachz*(vi**2+ai**2)
669  380 CONTINUE
670 
671  ELSEIF(isub.EQ.25) THEN
672 C...f + fb -> W+ + W-.
673  facww=comfac*faca*(aem/xw)**2*1./12.
674  facww=facww*wids(24,1)
675  thuh=max(th*uh-sqm3*sqm4,sh*ckin(3)**2)
676  DO 390 i=mina,maxa
677  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 390
678  ei=kchg(iabs(i),1)/3.
679  ai=sign(1.,ei)
680  vi=ai-4.*ei*xw
681  dsigww=thuh/sh2*(3.-(sh-3.*(sqm3+sqm4))/(sh-sqmz)*
682  & (vi+ai)/(2.*ai*(1.-xw))+(sh/(sh-sqmz))**2*
683  & (1.-2.*(sqm3+sqm4)/sh+12.*sqm3*sqm4/sh2)*(vi**2+ai**2)/
684  & (8.*(1.-xw)**2))-2.*sqmz/(sh-sqmz)*(vi+ai)/ai+
685  & sqmz*sh/(sh-sqmz)**2*(1.-2.*(sqm3+sqm4)/sh)*(vi**2+ai**2)/
686  & (2.*(1.-xw))
687  IF(kchg(iabs(i),1).LT.0) THEN
688  dsigww=dsigww+2.*(1.+sqmz/(sh-sqmz)*(vi+ai)/(2.*ai))*
689  & (thuh/(sh*th)-(sqm3+sqm4)/th)+thuh/th2
690  ELSE
691  dsigww=dsigww+2.*(1.+sqmz/(sh-sqmz)*(vi+ai)/(2.*ai))*
692  & (thuh/(sh*uh)-(sqm3+sqm4)/uh)+thuh/uh2
693  ENDIF
694  nchn=nchn+1
695  isig(nchn,1)=i
696  isig(nchn,2)=-i
697  isig(nchn,3)=1
698  sigh(nchn)=facww*dsigww
699  390 CONTINUE
700 
701  ELSEIF(isub.EQ.26) THEN
702 C...f + fb' -> W+/- + H0.
703  thuh=max(th*uh-sqm3*sqm4,sh*ckin(3)**2)
704  fachw=comfac*faca*(aem/xw)**2*1./24.*(thuh+2.*sh*sqmw)/
705  & (sh-sqmw)**2
706  fachw=fachw*wids(25,2)
707  DO 410 i=min1,max1
708  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 410
709  ia=iabs(i)
710  DO 400 j=min2,max2
711  IF(j.EQ.0.OR.kfac(1,j).EQ.0) goto 400
712  ja=iabs(j)
713  IF(i*j.GT.0.OR.mod(ia+ja,2).EQ.0) goto 400
714  kchw=(kchg(ia,1)*isign(1,i)+kchg(ja,1)*isign(1,j))/3
715  fckm=1.
716  IF(mint(43).EQ.4) fckm=vckm((ia+1)/2,(ja+1)/2)
717  nchn=nchn+1
718  isig(nchn,1)=i
719  isig(nchn,2)=j
720  isig(nchn,3)=1
721  sigh(nchn)=fachw*fckm*wids(24,(5-kchw)/2)
722  400 CONTINUE
723  410 CONTINUE
724 
725  ELSEIF(isub.EQ.27) THEN
726 C...f + fb -> H0 + H0.
727 
728  ELSEIF(isub.EQ.28) THEN
729 C...f + g -> f + g (q + g -> q + g only).
730  facqg1=comfac*as**2*4./9.*((2.+mstp(34)*1./4.)*uh2/th2-uh/sh)*
731  & faca
732  facqg2=comfac*as**2*4./9.*((2.+mstp(34)*1./4.)*sh2/th2-sh/uh)
733  DO 430 i=mina,maxa
734  IF(i.EQ.0) goto 430
735  DO 420 isde=1,2
736  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 420
737  IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 420
738  nchn=nchn+1
739  isig(nchn,isde)=i
740  isig(nchn,3-isde)=21
741  isig(nchn,3)=1
742  sigh(nchn)=facqg1
743  nchn=nchn+1
744  isig(nchn,isde)=i
745  isig(nchn,3-isde)=21
746  isig(nchn,3)=2
747  sigh(nchn)=facqg2
748  420 CONTINUE
749  430 CONTINUE
750 
751  ELSEIF(isub.EQ.29) THEN
752 C...f + g -> f + gamma (q + g -> q + gamma only).
753  fgq=comfac*faca*as*aem*1./3.*(sh2+uh2)/(-sh*uh)
754  DO 450 i=mina,maxa
755  IF(i.EQ.0) goto 450
756  ei=kchg(iabs(i),1)/3.
757  facgq=fgq*ei**2
758  DO 440 isde=1,2
759  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 440
760  IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 440
761  nchn=nchn+1
762  isig(nchn,isde)=i
763  isig(nchn,3-isde)=21
764  isig(nchn,3)=1
765  sigh(nchn)=facgq
766  440 CONTINUE
767  450 CONTINUE
768 
769  ELSEIF(isub.EQ.30) THEN
770 C...f + g -> f + Z0 (q + g -> q + Z0 only).
771  fzq=comfac*faca*as*aem/(xw*(1.-xw))*1./48.*
772  & (sh2+uh2+2.*sqm4*th)/(-sh*uh)
773  fzq=fzq*wids(23,2)
774  DO 470 i=mina,maxa
775  IF(i.EQ.0) goto 470
776  ei=kchg(iabs(i),1)/3.
777  ai=sign(1.,ei)
778  vi=ai-4.*ei*xw
779  faczq=fzq*(vi**2+ai**2)
780  DO 460 isde=1,2
781  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 460
782  IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 460
783  nchn=nchn+1
784  isig(nchn,isde)=i
785  isig(nchn,3-isde)=21
786  isig(nchn,3)=1
787  sigh(nchn)=faczq
788  460 CONTINUE
789  470 CONTINUE
790  ENDIF
791 
792  ELSEIF(isub.LE.40) THEN
793  IF(isub.EQ.31) THEN
794 C...f + g -> f' + W+/- (q + g -> q' + W+/- only).
795  facwq=comfac*faca*as*aem/xw*1./12.*
796  & (sh2+uh2+2.*sqm4*th)/(-sh*uh)
797  DO 490 i=mina,maxa
798  IF(i.EQ.0) goto 490
799  ia=iabs(i)
800  kchw=isign(1,kchg(ia,1)*isign(1,i))
801  DO 480 isde=1,2
802  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 480
803  IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 480
804  nchn=nchn+1
805  isig(nchn,isde)=i
806  isig(nchn,3-isde)=21
807  isig(nchn,3)=1
808  sigh(nchn)=facwq*vint(180+i)*wids(24,(5-kchw)/2)
809  480 CONTINUE
810  490 CONTINUE
811 
812  ELSEIF(isub.EQ.32) THEN
813 C...f + g -> f + H0 (q + g -> q + H0 only).
814 
815  ELSEIF(isub.EQ.33) THEN
816 C...f + gamma -> f + g (q + gamma -> q + g only).
817 
818  ELSEIF(isub.EQ.34) THEN
819 C...f + gamma -> f + gamma.
820 
821  ELSEIF(isub.EQ.35) THEN
822 C...f + gamma -> f + Z0.
823 
824  ELSEIF(isub.EQ.36) THEN
825 C...f + gamma -> f' + W+/-.
826 
827  ELSEIF(isub.EQ.37) THEN
828 C...f + gamma -> f + H0.
829 
830  ELSEIF(isub.EQ.38) THEN
831 C...f + Z0 -> f + g (q + Z0 -> q + g only).
832 
833  ELSEIF(isub.EQ.39) THEN
834 C...f + Z0 -> f + gamma.
835 
836  ELSEIF(isub.EQ.40) THEN
837 C...f + Z0 -> f + Z0.
838  ENDIF
839 
840  ELSEIF(isub.LE.50) THEN
841  IF(isub.EQ.41) THEN
842 C...f + Z0 -> f' + W+/-.
843 
844  ELSEIF(isub.EQ.42) THEN
845 C...f + Z0 -> f + H0.
846 
847  ELSEIF(isub.EQ.43) THEN
848 C...f + W+/- -> f' + g (q + W+/- -> q' + g only).
849 
850  ELSEIF(isub.EQ.44) THEN
851 C...f + W+/- -> f' + gamma.
852 
853  ELSEIF(isub.EQ.45) THEN
854 C...f + W+/- -> f' + Z0.
855 
856  ELSEIF(isub.EQ.46) THEN
857 C...f + W+/- -> f' + W+/-.
858 
859  ELSEIF(isub.EQ.47) THEN
860 C...f + W+/- -> f' + H0.
861 
862  ELSEIF(isub.EQ.48) THEN
863 C...f + H0 -> f + g (q + H0 -> q + g only).
864 
865  ELSEIF(isub.EQ.49) THEN
866 C...f + H0 -> f + gamma.
867 
868  ELSEIF(isub.EQ.50) THEN
869 C...f + H0 -> f + Z0.
870  ENDIF
871 
872  ELSEIF(isub.LE.60) THEN
873  IF(isub.EQ.51) THEN
874 C...f + H0 -> f' + W+/-.
875 
876  ELSEIF(isub.EQ.52) THEN
877 C...f + H0 -> f + H0.
878 
879  ELSEIF(isub.EQ.53) THEN
880 C...g + g -> f + fb (g + g -> q + qb only).
881  CALL pyhiwidt(21,sqrt(sh),wdtp,wdte)
882  facqq1=comfac*as**2*1./6.*(uh/th-(2.+mstp(34)*1./4.)*uh2/sh2)*
883  & (wdte(0,1)+wdte(0,2)+wdte(0,3)+wdte(0,4))*faca
884  facqq2=comfac*as**2*1./6.*(th/uh-(2.+mstp(34)*1./4.)*th2/sh2)*
885  & (wdte(0,1)+wdte(0,2)+wdte(0,3)+wdte(0,4))*faca
886  IF(kfac(1,21)*kfac(2,21).EQ.0) goto 500
887  nchn=nchn+1
888  isig(nchn,1)=21
889  isig(nchn,2)=21
890  isig(nchn,3)=1
891  sigh(nchn)=facqq1
892  nchn=nchn+1
893  isig(nchn,1)=21
894  isig(nchn,2)=21
895  isig(nchn,3)=2
896  sigh(nchn)=facqq2
897  500 CONTINUE
898 
899  ELSEIF(isub.EQ.54) THEN
900 C...g + gamma -> f + fb (g + gamma -> q + qb only).
901 
902  ELSEIF(isub.EQ.55) THEN
903 C...g + gamma -> f + fb (g + gamma -> q + qb only).
904 
905  ELSEIF(isub.EQ.56) THEN
906 C...g + gamma -> f + fb (g + gamma -> q + qb only).
907 
908  ELSEIF(isub.EQ.57) THEN
909 C...g + gamma -> f + fb (g + gamma -> q + qb only).
910 
911  ELSEIF(isub.EQ.58) THEN
912 C...gamma + gamma -> f + fb.
913 
914  ELSEIF(isub.EQ.59) THEN
915 C...gamma + Z0 -> f + fb.
916 
917  ELSEIF(isub.EQ.60) THEN
918 C...gamma + W+/- -> f + fb'.
919  ENDIF
920 
921  ELSEIF(isub.LE.70) THEN
922  IF(isub.EQ.61) THEN
923 C...gamma + H0 -> f + fb.
924 
925  ELSEIF(isub.EQ.62) THEN
926 C...Z0 + Z0 -> f + fb.
927 
928  ELSEIF(isub.EQ.63) THEN
929 C...Z0 + W+/- -> f + fb'.
930 
931  ELSEIF(isub.EQ.64) THEN
932 C...Z0 + H0 -> f + fb.
933 
934  ELSEIF(isub.EQ.65) THEN
935 C...W+ + W- -> f + fb.
936 
937  ELSEIF(isub.EQ.66) THEN
938 C...W+/- + H0 -> f + fb'.
939 
940  ELSEIF(isub.EQ.67) THEN
941 C...H0 + H0 -> f + fb.
942 
943  ELSEIF(isub.EQ.68) THEN
944 C...g + g -> g + g.
945  facgg1=comfac*as**2*9./4.*(sh2/th2+2.*sh/th+3.+2.*th/sh+
946  & th2/sh2)*faca
947  facgg2=comfac*as**2*9./4.*(uh2/sh2+2.*uh/sh+3.+2.*sh/uh+
948  & sh2/uh2)*faca
949  facgg3=comfac*as**2*9./4.*(th2/uh2+2.*th/uh+3+2.*uh/th+uh2/th2)
950  IF(kfac(1,21)*kfac(2,21).EQ.0) goto 510
951  nchn=nchn+1
952  isig(nchn,1)=21
953  isig(nchn,2)=21
954  isig(nchn,3)=1
955  sigh(nchn)=0.5*facgg1
956  nchn=nchn+1
957  isig(nchn,1)=21
958  isig(nchn,2)=21
959  isig(nchn,3)=2
960  sigh(nchn)=0.5*facgg2
961  nchn=nchn+1
962  isig(nchn,1)=21
963  isig(nchn,2)=21
964  isig(nchn,3)=3
965  sigh(nchn)=0.5*facgg3
966  510 CONTINUE
967 
968  ELSEIF(isub.EQ.69) THEN
969 C...gamma + gamma -> W+ + W-.
970 
971  ELSEIF(isub.EQ.70) THEN
972 C...gamma + W+/- -> gamma + W+/-.
973  ENDIF
974 
975  ELSEIF(isub.LE.80) THEN
976  IF(isub.EQ.71) THEN
977 C...Z0 + Z0 -> Z0 + Z0.
978  be2=1.-4.*sqmz/sh
979  th=-0.5*sh*be2*(1.-cth)
980  uh=-0.5*sh*be2*(1.+cth)
981  shang=1./(1.-xw)*sqmw/sqmz*(1.+be2)**2
982  ashre=(sh-sqmh)/((sh-sqmh)**2+gmmh**2)*shang
983  ashim=-gmmh/((sh-sqmh)**2+gmmh**2)*shang
984  thang=1./(1.-xw)*sqmw/sqmz*(be2-cth)**2
985  athre=(th-sqmh)/((th-sqmh)**2+gmmh**2)*thang
986  athim=-gmmh/((th-sqmh)**2+gmmh**2)*thang
987  uhang=1./(1.-xw)*sqmw/sqmz*(be2+cth)**2
988  auhre=(uh-sqmh)/((uh-sqmh)**2+gmmh**2)*uhang
989  auhim=-gmmh/((uh-sqmh)**2+gmmh**2)*uhang
990  fach=0.5*comfac*1./(4096.*paru(1)**2*16.*(1.-xw)**2)*
991  & (aem/xw)**4*(sh/sqmw)**2*((ashre+athre+auhre)**2+
992  & (ashim+athim+auhim)**2)*sqmz/sqmw
993  DO 530 i=min1,max1
994  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 530
995  ei=kchg(iabs(i),1)/3.
996  ai=sign(1.,ei)
997  vi=ai-4.*ei*xw
998  avi=ai**2+vi**2
999  DO 520 j=min2,max2
1000  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 520
1001  ej=kchg(iabs(j),1)/3.
1002  aj=sign(1.,ej)
1003  vj=aj-4.*ej*xw
1004  avj=aj**2+vj**2
1005  nchn=nchn+1
1006  isig(nchn,1)=i
1007  isig(nchn,2)=j
1008  isig(nchn,3)=1
1009  sigh(nchn)=fach*avi*avj
1010  520 CONTINUE
1011  530 CONTINUE
1012 
1013  ELSEIF(isub.EQ.72) THEN
1014 C...Z0 + Z0 -> W+ + W-.
1015  be2=sqrt((1.-4.*sqmw/sh)*(1.-4.*sqmz/sh))
1016  cth2=cth**2
1017  th=-0.5*sh*(1.-2.*(sqmw+sqmz)/sh-be2*cth)
1018  uh=-0.5*sh*(1.-2.*(sqmw+sqmz)/sh+be2*cth)
1019  shang=4.*sqrt(sqmw/(sqmz*(1.-xw)))*(1.-2.*sqmw/sh)*
1020  & (1.-2.*sqmz/sh)
1021  ashre=(sh-sqmh)/((sh-sqmh)**2+gmmh**2)*shang
1022  ashim=-gmmh/((sh-sqmh)**2+gmmh**2)*shang
1023  atwre=(1.-xw)/sqmz*sh/(th-sqmw)*((cth-be2)**2*(3./2.+be2/2.*cth-
1024  & (sqmw+sqmz)/sh+(sqmw-sqmz)**2/(sh*sqmw))+4.*((sqmw+sqmz)/sh*
1025  & (1.-3.*cth2)+8.*sqmw*sqmz/sh2*(2.*cth2-1.)+
1026  & 4.*(sqmw**2+sqmz**2)/sh2*cth2+2.*(sqmw+sqmz)/sh*be2*cth))
1027  atwim=0.
1028  auwre=(1.-xw)/sqmz*sh/(uh-sqmw)*((cth+be2)**2*(3./2.-be2/2.*cth-
1029  & (sqmw+sqmz)/sh+(sqmw-sqmz)**2/(sh*sqmw))+4.*((sqmw+sqmz)/sh*
1030  & (1.-3.*cth2)+8.*sqmw*sqmz/sh2*(2.*cth2-1.)+
1031  & 4.*(sqmw**2+sqmz**2)/sh2*cth2-2.*(sqmw+sqmz)/sh*be2*cth))
1032  auwim=0.
1033  a4re=2.*(1.-xw)/sqmz*(3.-cth2-4.*(sqmw+sqmz)/sh)
1034  a4im=0.
1035  fach=comfac*1./(4096.*paru(1)**2*16.*(1.-xw)**2)*(aem/xw)**4*
1036  & (sh/sqmw)**2*((ashre+atwre+auwre+a4re)**2+
1037  & (ashim+atwim+auwim+a4im)**2)*sqmz/sqmw
1038  DO 550 i=min1,max1
1039  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 550
1040  ei=kchg(iabs(i),1)/3.
1041  ai=sign(1.,ei)
1042  vi=ai-4.*ei*xw
1043  avi=ai**2+vi**2
1044  DO 540 j=min2,max2
1045  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 540
1046  ej=kchg(iabs(j),1)/3.
1047  aj=sign(1.,ej)
1048  vj=aj-4.*ej*xw
1049  avj=aj**2+vj**2
1050  nchn=nchn+1
1051  isig(nchn,1)=i
1052  isig(nchn,2)=j
1053  isig(nchn,3)=1
1054  sigh(nchn)=fach*avi*avj
1055  540 CONTINUE
1056  550 CONTINUE
1057 
1058  ELSEIF(isub.EQ.73) THEN
1059 C...Z0 + W+/- -> Z0 + W+/-.
1060  be2=1.-2.*(sqmz+sqmw)/sh+((sqmz-sqmw)/sh)**2
1061  ep1=1.+(sqmz-sqmw)/sh
1062  ep2=1.-(sqmz-sqmw)/sh
1063  th=-0.5*sh*be2*(1.-cth)
1064  uh=(sqmz-sqmw)**2/sh-0.5*sh*be2*(1.+cth)
1065  thang=sqrt(sqmw/(sqmz*(1.-xw)))*(be2-ep1*cth)*(be2-ep2*cth)
1066  athre=(th-sqmh)/((th-sqmh)**2+gmmh**2)*thang
1067  athim=-gmmh/((th-sqmh)**2+gmmh**2)*thang
1068  aswre=(1.-xw)/sqmz*sh/(sh-sqmw)*(-be2*(ep1+ep2)**4*cth+
1069  & 1./4.*(be2+ep1*ep2)**2*((ep1-ep2)**2-4.*be2*cth)+
1070  & 2.*be2*(be2+ep1*ep2)*(ep1+ep2)**2*cth-
1071  & 1./16.*sh/sqmw*(ep1**2-ep2**2)**2*(be2+ep1*ep2)**2)
1072  aswim=0.
1073  auwre=(1.-xw)/sqmz*sh/(uh-sqmw)*(-be2*(ep2+ep1*cth)*
1074  & (ep1+ep2*cth)*(be2+ep1*ep2)+be2*(ep2+ep1*cth)*
1075  & (be2+ep1*ep2*cth)*(2.*ep2-ep2*cth+ep1)-be2*(ep2+ep1*cth)**2*
1076  & (be2-ep2**2*cth)-1./8.*(be2+ep1*ep2*cth)**2*((ep1+ep2)**2+
1077  & 2.*be2*(1.-cth))+1./32.*sh/sqmw*(be2+ep1*ep2*cth)**2*
1078  & (ep1**2-ep2**2)**2-be2*(ep1+ep2*cth)*(ep2+ep1*cth)*
1079  & (be2+ep1*ep2)+be2*(ep1+ep2*cth)*(be2+ep1*ep2*cth)*
1080  & (2.*ep1-ep1*cth+ep2)-be2*(ep1+ep2*cth)**2*(be2-ep1**2*cth)-
1081  & 1./8.*(be2+ep1*ep2*cth)**2*((ep1+ep2)**2+2.*be2*(1.-cth))+
1082  & 1./32.*sh/sqmw*(be2+ep1*ep2*cth)**2*(ep1**2-ep2**2)**2)
1083  auwim=0.
1084  a4re=(1.-xw)/sqmz*(ep1**2*ep2**2*(cth**2-1.)-
1085  & 2.*be2*(ep1**2+ep2**2+ep1*ep2)*cth-2.*be2*ep1*ep2)
1086  a4im=0.
1087  fach=comfac*1./(4096.*paru(1)**2*4.*(1.-xw))*(aem/xw)**4*
1088  & (sh/sqmw)**2*((athre+aswre+auwre+a4re)**2+
1089  & (athim+aswim+auwim+a4im)**2)*sqrt(sqmz/sqmw)
1090  DO 570 i=min1,max1
1091  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 570
1092  ei=kchg(iabs(i),1)/3.
1093  ai=sign(1.,ei)
1094  vi=ai-4.*ei*xw
1095  avi=ai**2+vi**2
1096  DO 560 j=min2,max2
1097  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 560
1098  ej=kchg(iabs(j),1)/3.
1099  aj=sign(1.,ej)
1100  vj=ai-4.*ej*xw
1101  avj=aj**2+vj**2
1102  nchn=nchn+1
1103  isig(nchn,1)=i
1104  isig(nchn,2)=j
1105  isig(nchn,3)=1
1106  sigh(nchn)=fach*(avi*vint(180+j)+vint(180+i)*avj)
1107  560 CONTINUE
1108  570 CONTINUE
1109 
1110  ELSEIF(isub.EQ.75) THEN
1111 C...W+ + W- -> gamma + gamma.
1112 
1113  ELSEIF(isub.EQ.76) THEN
1114 C...W+ + W- -> Z0 + Z0.
1115  be2=sqrt((1.-4.*sqmw/sh)*(1.-4.*sqmz/sh))
1116  cth2=cth**2
1117  th=-0.5*sh*(1.-2.*(sqmw+sqmz)/sh-be2*cth)
1118  uh=-0.5*sh*(1.-2.*(sqmw+sqmz)/sh+be2*cth)
1119  shang=4.*sqrt(sqmw/(sqmz*(1.-xw)))*(1.-2.*sqmw/sh)*
1120  & (1.-2.*sqmz/sh)
1121  ashre=(sh-sqmh)/((sh-sqmh)**2+gmmh**2)*shang
1122  ashim=-gmmh/((sh-sqmh)**2+gmmh**2)*shang
1123  atwre=(1.-xw)/sqmz*sh/(th-sqmw)*((cth-be2)**2*(3./2.+be2/2.*cth-
1124  & (sqmw+sqmz)/sh+(sqmw-sqmz)**2/(sh*sqmw))+4.*((sqmw+sqmz)/sh*
1125  & (1.-3.*cth2)+8.*sqmw*sqmz/sh2*(2.*cth2-1.)+
1126  & 4.*(sqmw**2+sqmz**2)/sh2*cth2+2.*(sqmw+sqmz)/sh*be2*cth))
1127  atwim=0.
1128  auwre=(1.-xw)/sqmz*sh/(uh-sqmw)*((cth+be2)**2*(3./2.-be2/2.*cth-
1129  & (sqmw+sqmz)/sh+(sqmw-sqmz)**2/(sh*sqmw))+4.*((sqmw+sqmz)/sh*
1130  & (1.-3.*cth2)+8.*sqmw*sqmz/sh2*(2.*cth2-1.)+
1131  & 4.*(sqmw**2+sqmz**2)/sh2*cth2-2.*(sqmw+sqmz)/sh*be2*cth))
1132  auwim=0.
1133  a4re=2.*(1.-xw)/sqmz*(3.-cth2-4.*(sqmw+sqmz)/sh)
1134  a4im=0.
1135  fach=0.5*comfac*1./(4096.*paru(1)**2)*(aem/xw)**4*(sh/sqmw)**2*
1136  & ((ashre+atwre+auwre+a4re)**2+(ashim+atwim+auwim+a4im)**2)
1137  DO 590 i=min1,max1
1138  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 590
1139  ei=sign(1.,float(i))*kchg(iabs(i),1)
1140  DO 580 j=min2,max2
1141  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 580
1142  ej=sign(1.,float(j))*kchg(iabs(j),1)
1143  IF(ei*ej.GT.0.) goto 580
1144  nchn=nchn+1
1145  isig(nchn,1)=i
1146  isig(nchn,2)=j
1147  isig(nchn,3)=1
1148  sigh(nchn)=fach*vint(180+i)*vint(180+j)
1149  580 CONTINUE
1150  590 CONTINUE
1151 
1152  ELSEIF(isub.EQ.77) THEN
1153 C...W+/- + W+/- -> W+/- + W+/-.
1154  be2=1.-4.*sqmw/sh
1155  be4=be2**2
1156  cth2=cth**2
1157  cth3=cth**3
1158  th=-0.5*sh*be2*(1.-cth)
1159  uh=-0.5*sh*be2*(1.+cth)
1160  shang=(1.+be2)**2
1161  ashre=(sh-sqmh)/((sh-sqmh)**2+gmmh**2)*shang
1162  ashim=-gmmh/((sh-sqmh)**2+gmmh**2)*shang
1163  thang=(be2-cth)**2
1164  athre=(th-sqmh)/((th-sqmh)**2+gmmh**2)*thang
1165  athim=-gmmh/((th-sqmh)**2+gmmh**2)*thang
1166  sgzang=1./sqmw*be2*(3.-be2)**2*cth
1167  asgre=xw*sgzang
1168  asgim=0.
1169  aszre=(1.-xw)*sh/(sh-sqmz)*sgzang
1170  aszim=0.
1171  tgzang=1./sqmw*(be2*(4.-2.*be2+be4)+be2*(4.-10.*be2+be4)*cth+
1172  & (2.-11.*be2+10.*be4)*cth2+be2*cth3)
1173  atgre=0.5*xw*sh/th*tgzang
1174  atgim=0.
1175  atzre=0.5*(1.-xw)*sh/(th-sqmz)*tgzang
1176  atzim=0.
1177  a4re=1./sqmw*(1.+2.*be2-6.*be2*cth-cth2)
1178  a4im=0.
1179  fach=comfac*1./(4096.*paru(1)**2)*(aem/xw)**4*(sh/sqmw)**2*
1180  & ((ashre+athre+asgre+aszre+atgre+atzre+a4re)**2+
1181  & (ashim+athim+asgim+aszim+atgim+atzim+a4im)**2)
1182  DO 610 i=min1,max1
1183  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 610
1184  ei=sign(1.,float(i))*kchg(iabs(i),1)
1185  DO 600 j=min2,max2
1186  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 600
1187  ej=sign(1.,float(j))*kchg(iabs(j),1)
1188  IF(ei*ej.GT.0.) goto 600
1189  nchn=nchn+1
1190  isig(nchn,1)=i
1191  isig(nchn,2)=j
1192  isig(nchn,3)=1
1193  sigh(nchn)=fach*vint(180+i)*vint(180+j)
1194  600 CONTINUE
1195  610 CONTINUE
1196 
1197  ELSEIF(isub.EQ.78) THEN
1198 C...W+/- + H0 -> W+/- + H0.
1199 
1200  ELSEIF(isub.EQ.79) THEN
1201 C...H0 + H0 -> H0 + H0.
1202 
1203  ENDIF
1204 
1205 C...C: 2 -> 2, tree diagrams with masses.
1206 
1207  ELSEIF(isub.LE.90) THEN
1208  IF(isub.EQ.81) THEN
1209 C...q + qb -> Q + QB.
1210  facqqb=comfac*as**2*4./9.*(((th-sqm3)**2+
1211  & (uh-sqm3)**2)/sh2+2.*sqm3/sh)
1212  IF(mstp(35).GE.1) THEN
1213  IF(mstp(35).EQ.1) THEN
1214  alssg=parp(35)
1215  ELSE
1216  mst115=mstu(115)
1217  mstu(115)=mstp(36)
1218  q2bn=sqrt(sqm3*((sqrt(sh)-2.*sqrt(sqm3))**2+parp(36)**2))
1219  alssg=ulalps(q2bn)
1220  mstu(115)=mst115
1221  ENDIF
1222  xrepu=paru(1)*alssg/(6.*sqrt(max(1e-20,1.-4.*sqm3/sh)))
1223  frepu=xrepu/(exp(min(100.,xrepu))-1.)
1224  pari(81)=frepu
1225  facqqb=facqqb*frepu
1226  ENDIF
1227  DO 620 i=mina,maxa
1228  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 620
1229  nchn=nchn+1
1230  isig(nchn,1)=i
1231  isig(nchn,2)=-i
1232  isig(nchn,3)=1
1233  sigh(nchn)=facqqb
1234  620 CONTINUE
1235 
1236  ELSEIF(isub.EQ.82) THEN
1237 C...g + g -> Q + QB.
1238  facqq1=comfac*faca*as**2*1./6.*((uh-sqm3)/(th-sqm3)-
1239  & 2.*(uh-sqm3)**2/sh2+4.*sqm3/sh*(th*uh-sqm3**2)/(th-sqm3)**2)
1240  facqq2=comfac*faca*as**2*1./6.*((th-sqm3)/(uh-sqm3)-
1241  & 2.*(th-sqm3)**2/sh2+4.*sqm3/sh*(th*uh-sqm3**2)/(uh-sqm3)**2)
1242  IF(mstp(35).GE.1) THEN
1243  IF(mstp(35).EQ.1) THEN
1244  alssg=parp(35)
1245  ELSE
1246  mst115=mstu(115)
1247  mstu(115)=mstp(36)
1248  q2bn=sqrt(sqm3*((sqrt(sh)-2.*sqrt(sqm3))**2+parp(36)**2))
1249  alssg=ulalps(q2bn)
1250  mstu(115)=mst115
1251  ENDIF
1252  xattr=4.*paru(1)*alssg/(3.*sqrt(max(1e-20,1.-4.*sqm3/sh)))
1253  fattr=xattr/(1.-exp(-min(100.,xattr)))
1254  xrepu=paru(1)*alssg/(6.*sqrt(max(1e-20,1.-4.*sqm3/sh)))
1255  frepu=xrepu/(exp(min(100.,xrepu))-1.)
1256  fatre=(2.*fattr+5.*frepu)/7.
1257  pari(81)=fatre
1258  facqq1=facqq1*fatre
1259  facqq2=facqq2*fatre
1260  ENDIF
1261  IF(kfac(1,21)*kfac(2,21).EQ.0) goto 630
1262  nchn=nchn+1
1263  isig(nchn,1)=21
1264  isig(nchn,2)=21
1265  isig(nchn,3)=1
1266  sigh(nchn)=facqq1
1267  nchn=nchn+1
1268  isig(nchn,1)=21
1269  isig(nchn,2)=21
1270  isig(nchn,3)=2
1271  sigh(nchn)=facqq2
1272  630 CONTINUE
1273 
1274  ENDIF
1275 
1276 C...D: Mimimum bias processes.
1277 
1278  ELSEIF(isub.LE.100) THEN
1279  IF(isub.EQ.91) THEN
1280 C...Elastic scattering.
1281  sigs=xsec(isub,1)
1282 
1283  ELSEIF(isub.EQ.92) THEN
1284 C...Single diffractive scattering.
1285  sigs=xsec(isub,1)
1286 
1287  ELSEIF(isub.EQ.93) THEN
1288 C...Double diffractive scattering.
1289  sigs=xsec(isub,1)
1290 
1291  ELSEIF(isub.EQ.94) THEN
1292 C...Central diffractive scattering.
1293  sigs=xsec(isub,1)
1294 
1295  ELSEIF(isub.EQ.95) THEN
1296 C...Low-pT scattering.
1297  sigs=xsec(isub,1)
1298 
1299  ELSEIF(isub.EQ.96) THEN
1300 C...Multiple interactions: sum of QCD processes.
1301  CALL pyhiwidt(21,sqrt(sh),wdtp,wdte)
1302 
1303 C...q + q' -> q + q'.
1304  facqq1=comfac*as**2*4./9.*(sh2+uh2)/th2
1305  facqqb=comfac*as**2*4./9.*((sh2+uh2)/th2*faca-
1306  & mstp(34)*2./3.*uh2/(sh*th))
1307  facqq2=comfac*as**2*4./9.*((sh2+th2)/uh2-
1308  & mstp(34)*2./3.*sh2/(th*uh))
1309  DO 650 i=-3,3
1310  IF(i.EQ.0) goto 650
1311  DO 640 j=-3,3
1312  IF(j.EQ.0) goto 640
1313  nchn=nchn+1
1314  isig(nchn,1)=i
1315  isig(nchn,2)=j
1316  isig(nchn,3)=111
1317  sigh(nchn)=facqq1
1318  IF(i.EQ.-j) sigh(nchn)=facqqb
1319  IF(i.EQ.j) THEN
1320  sigh(nchn)=0.5*sigh(nchn)
1321  nchn=nchn+1
1322  isig(nchn,1)=i
1323  isig(nchn,2)=j
1324  isig(nchn,3)=112
1325  sigh(nchn)=0.5*facqq2
1326  ENDIF
1327  640 CONTINUE
1328  650 CONTINUE
1329 
1330 C...q + qb -> q' + qb' or g + g.
1331  facqqb=comfac*as**2*4./9.*(th2+uh2)/sh2*(wdte(0,1)+wdte(0,2)+
1332  & wdte(0,3)+wdte(0,4))
1333  facgg1=comfac*as**2*32./27.*(uh/th-(2.+mstp(34)*1./4.)*uh2/sh2)
1334  facgg2=comfac*as**2*32./27.*(th/uh-(2.+mstp(34)*1./4.)*th2/sh2)
1335  DO 660 i=-3,3
1336  IF(i.EQ.0) goto 660
1337  nchn=nchn+1
1338  isig(nchn,1)=i
1339  isig(nchn,2)=-i
1340  isig(nchn,3)=121
1341  sigh(nchn)=facqqb
1342  nchn=nchn+1
1343  isig(nchn,1)=i
1344  isig(nchn,2)=-i
1345  isig(nchn,3)=131
1346  sigh(nchn)=0.5*facgg1
1347  nchn=nchn+1
1348  isig(nchn,1)=i
1349  isig(nchn,2)=-i
1350  isig(nchn,3)=132
1351  sigh(nchn)=0.5*facgg2
1352  660 CONTINUE
1353 
1354 C...q + g -> q + g.
1355  facqg1=comfac*as**2*4./9.*((2.+mstp(34)*1./4.)*uh2/th2-uh/sh)*
1356  & faca
1357  facqg2=comfac*as**2*4./9.*((2.+mstp(34)*1./4.)*sh2/th2-sh/uh)
1358  DO 680 i=-3,3
1359  IF(i.EQ.0) goto 680
1360  DO 670 isde=1,2
1361  nchn=nchn+1
1362  isig(nchn,isde)=i
1363  isig(nchn,3-isde)=21
1364  isig(nchn,3)=281
1365  sigh(nchn)=facqg1
1366  nchn=nchn+1
1367  isig(nchn,isde)=i
1368  isig(nchn,3-isde)=21
1369  isig(nchn,3)=282
1370  sigh(nchn)=facqg2
1371  670 CONTINUE
1372  680 CONTINUE
1373 
1374 C...g + g -> q + qb or g + g.
1375  facqq1=comfac*as**2*1./6.*(uh/th-(2.+mstp(34)*1./4.)*uh2/sh2)*
1376  & (wdte(0,1)+wdte(0,2)+wdte(0,3)+wdte(0,4))*faca
1377  facqq2=comfac*as**2*1./6.*(th/uh-(2.+mstp(34)*1./4.)*th2/sh2)*
1378  & (wdte(0,1)+wdte(0,2)+wdte(0,3)+wdte(0,4))*faca
1379  facgg1=comfac*as**2*9./4.*(sh2/th2+2.*sh/th+3.+2.*th/sh+
1380  & th2/sh2)*faca
1381  facgg2=comfac*as**2*9./4.*(uh2/sh2+2.*uh/sh+3.+2.*sh/uh+
1382  & sh2/uh2)*faca
1383  facgg3=comfac*as**2*9./4.*(th2/uh2+2.*th/uh+3+2.*uh/th+uh2/th2)
1384  nchn=nchn+1
1385  isig(nchn,1)=21
1386  isig(nchn,2)=21
1387  isig(nchn,3)=531
1388  sigh(nchn)=facqq1
1389  nchn=nchn+1
1390  isig(nchn,1)=21
1391  isig(nchn,2)=21
1392  isig(nchn,3)=532
1393  sigh(nchn)=facqq2
1394  nchn=nchn+1
1395  isig(nchn,1)=21
1396  isig(nchn,2)=21
1397  isig(nchn,3)=681
1398  sigh(nchn)=0.5*facgg1
1399  nchn=nchn+1
1400  isig(nchn,1)=21
1401  isig(nchn,2)=21
1402  isig(nchn,3)=682
1403  sigh(nchn)=0.5*facgg2
1404  nchn=nchn+1
1405  isig(nchn,1)=21
1406  isig(nchn,2)=21
1407  isig(nchn,3)=683
1408  sigh(nchn)=0.5*facgg3
1409  ENDIF
1410 
1411 C...E: 2 -> 1, loop diagrams.
1412 
1413  ELSEIF(isub.LE.110) THEN
1414  IF(isub.EQ.101) THEN
1415 C...g + g -> gamma*/Z0.
1416 
1417  ELSEIF(isub.EQ.102) THEN
1418 C...g + g -> H0.
1419  CALL pyhiwidt(25,sqrt(sh),wdtp,wdte)
1420  etare=0.
1421  etaim=0.
1422  DO 690 i=1,2*mstp(1)
1423  eps=4.*pmas(i,1)**2/sh
1424  IF(eps.LE.1.) THEN
1425  IF(eps.GT.1.e-4) THEN
1426  root=sqrt(1.-eps)
1427  rln=log((1.+root)/(1.-root))
1428  ELSE
1429  rln=log(4./eps-2.)
1430  ENDIF
1431  phire=0.25*(rln**2-paru(1)**2)
1432  phiim=0.5*paru(1)*rln
1433  ELSE
1434  phire=-(asin(1./sqrt(eps)))**2
1435  phiim=0.
1436  ENDIF
1437  etare=etare+0.5*eps*(1.+(eps-1.)*phire)
1438  etaim=etaim+0.5*eps*(eps-1.)*phiim
1439  690 CONTINUE
1440  eta2=etare**2+etaim**2
1441  fach=comfac*faca*(as/paru(1)*aem/xw)**2*1./512.*
1442  & (sh/sqmw)**2*eta2*sh2/((sh-sqmh)**2+gmmh**2)*
1443  & (wdte(0,1)+wdte(0,2)+wdte(0,4))
1444  IF(kfac(1,21)*kfac(2,21).EQ.0) goto 700
1445  nchn=nchn+1
1446  isig(nchn,1)=21
1447  isig(nchn,2)=21
1448  isig(nchn,3)=1
1449  sigh(nchn)=fach
1450  700 CONTINUE
1451 
1452  ENDIF
1453 
1454 C...F: 2 -> 2, box diagrams.
1455 
1456  ELSEIF(isub.LE.120) THEN
1457  IF(isub.EQ.111) THEN
1458 C...f + fb -> g + H0 (q + qb -> g + H0 only).
1459  a5stur=0.
1460  a5stui=0.
1461  DO 710 i=1,2*mstp(1)
1462  sqmq=pmas(i,1)**2
1463  epss=4.*sqmq/sh
1464  epsh=4.*sqmq/sqmh
1465  a5stur=a5stur+sqmq/sqmh*(4.+4.*sh/(th+uh)*(pyhiw1au(epss,1)-
1466  & pyhiw1au(epsh,1))+(1.-4.*sqmq/(th+uh))*(pyhiw2au(epss,1)-
1467  & pyhiw2au(epsh,1)))
1468  a5stui=a5stui+sqmq/sqmh*(4.*sh/(th+uh)*(pyhiw1au(epss,2)-
1469  & pyhiw1au(epsh,2))+(1.-4.*sqmq/(th+uh))*(pyhiw2au(epss,2)-
1470  & pyhiw2au(epsh,2)))
1471  710 CONTINUE
1472  facgh=comfac*faca/(144.*paru(1)**2)*aem/xw*as**3*sqmh/sqmw*
1473  & sqmh/sh*(uh**2+th**2)/(uh+th)**2*(a5stur**2+a5stui**2)
1474  facgh=facgh*wids(25,2)
1475  DO 720 i=mina,maxa
1476  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 720
1477  nchn=nchn+1
1478  isig(nchn,1)=i
1479  isig(nchn,2)=-i
1480  isig(nchn,3)=1
1481  sigh(nchn)=facgh
1482  720 CONTINUE
1483 
1484  ELSEIF(isub.EQ.112) THEN
1485 C...f + g -> f + H0 (q + g -> q + H0 only).
1486  a5tsur=0.
1487  a5tsui=0.
1488  DO 730 i=1,2*mstp(1)
1489  sqmq=pmas(i,1)**2
1490  epst=4.*sqmq/th
1491  epsh=4.*sqmq/sqmh
1492  a5tsur=a5tsur+sqmq/sqmh*(4.+4.*th/(sh+uh)*(pyhiw1au(epst,1)-
1493  & pyhiw1au(epsh,1))+(1.-4.*sqmq/(sh+uh))*(pyhiw2au(epst,1)-
1494  & pyhiw2au(epsh,1)))
1495  a5tsui=a5tsui+sqmq/sqmh*(4.*th/(sh+uh)*(pyhiw1au(epst,2)-
1496  & pyhiw1au(epsh,2))+(1.-4.*sqmq/(sh+uh))*(pyhiw2au(epst,2)-
1497  & pyhiw2au(epsh,2)))
1498  730 CONTINUE
1499  facqh=comfac*faca/(384.*paru(1)**2)*aem/xw*as**3*sqmh/sqmw*
1500  & sqmh/(-th)*(uh**2+sh**2)/(uh+sh)**2*(a5tsur**2+a5tsui**2)
1501  facqh=facqh*wids(25,2)
1502  DO 750 i=mina,maxa
1503  IF(i.EQ.0) goto 750
1504  DO 740 isde=1,2
1505  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 740
1506  IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 740
1507  nchn=nchn+1
1508  isig(nchn,isde)=i
1509  isig(nchn,3-isde)=21
1510  isig(nchn,3)=1
1511  sigh(nchn)=facqh
1512  740 CONTINUE
1513  750 CONTINUE
1514 
1515  ELSEIF(isub.EQ.113) THEN
1516 C...g + g -> g + H0.
1517  a2stur=0.
1518  a2stui=0.
1519  a2ustr=0.
1520  a2usti=0.
1521  a2tusr=0.
1522  a2tusi=0.
1523  a4stur=0.
1524  a4stui=0.
1525  DO 760 i=6,2*mstp(1)
1526 C'''Only t-quarks yet included
1527  sqmq=pmas(i,1)**2
1528  epss=4.*sqmq/sh
1529  epst=4.*sqmq/th
1530  epsu=4.*sqmq/uh
1531  epsh=4.*sqmq/sqmh
1532  IF(epsh.LT.1.e-6) goto 760
1533  bestu=0.5*(1.+sqrt(1.+epss*th/uh))
1534  beust=0.5*(1.+sqrt(1.+epsu*sh/th))
1535  betus=0.5*(1.+sqrt(1.+epst*uh/sh))
1536  beuts=bestu
1537  betsu=beust
1538  besut=betus
1539  w3stur=pyhii3au(bestu,epsh,1)-pyhii3au(bestu,epss,1)-
1540  & pyhii3au(bestu,epsu,1)
1541  w3stui=pyhii3au(bestu,epsh,2)-pyhii3au(bestu,epss,2)-
1542  & pyhii3au(bestu,epsu,2)
1543  w3sutr=pyhii3au(besut,epsh,1)-pyhii3au(besut,epss,1)-
1544  & pyhii3au(besut,epst,1)
1545  w3suti=pyhii3au(besut,epsh,2)-pyhii3au(besut,epss,2)-
1546  & pyhii3au(besut,epst,2)
1547  w3tsur=pyhii3au(betsu,epsh,1)-pyhii3au(betsu,epst,1)-
1548  & pyhii3au(betsu,epsu,1)
1549  w3tsui=pyhii3au(betsu,epsh,2)-pyhii3au(betsu,epst,2)-
1550  & pyhii3au(betsu,epsu,2)
1551  w3tusr=pyhii3au(betus,epsh,1)-pyhii3au(betus,epst,1)-
1552  & pyhii3au(betus,epss,1)
1553  w3tusi=pyhii3au(betus,epsh,2)-pyhii3au(betus,epst,2)-
1554  & pyhii3au(betus,epss,2)
1555  w3ustr=pyhii3au(beust,epsh,1)-pyhii3au(beust,epsu,1)-
1556  & pyhii3au(beust,epst,1)
1557  w3usti=pyhii3au(beust,epsh,2)-pyhii3au(beust,epsu,2)-
1558  & pyhii3au(beust,epst,2)
1559  w3utsr=pyhii3au(beuts,epsh,1)-pyhii3au(beuts,epsu,1)-
1560  & pyhii3au(beuts,epss,1)
1561  w3utsi=pyhii3au(beuts,epsh,2)-pyhii3au(beuts,epsu,2)-
1562  & pyhii3au(beuts,epss,2)
1563  b2stur=sqmq/sqmh**2*(sh*(uh-sh)/(sh+uh)+2.*th*uh*(uh+2.*sh)/
1564  & (sh+uh)**2*(pyhiw1au(epst,1)-pyhiw1au(epsh,1))+(sqmq-sh/4.)*
1565  & (0.5*pyhiw2au(epss,1)+0.5*pyhiw2au(epsh,1)-pyhiw2au(epst,1)+
1566  & w3stur)+
1567  & sh**2*(2.*sqmq/(sh+uh)**2-0.5/(sh+uh))*(pyhiw2au(epst,1)-
1568  & pyhiw2au(epsh,1))+0.5*th*uh/sh*(pyhiw2au(epsh,1)-
1569  & 2.*pyhiw2au(epst,1))+
1570  & 0.125*(sh-12.*sqmq-4.*th*uh/sh)*w3tsur)
1571  b2stui=sqmq/sqmh**2*(2.*th*uh*(uh+2.*sh)/(sh+uh)**2*
1572  & (pyhiw1au(epst,2)-pyhiw1au(epsh,2))+(sqmq-sh/4.)*
1573  & (0.5*pyhiw2au(epss,2)+0.5*pyhiw2au(epsh,2)-pyhiw2au(epst,2)+
1574  & w3stui)+
1575  & sh**2*(2.*sqmq/(sh+uh)**2-0.5/(sh+uh))*(pyhiw2au(epst,2)-
1576  & pyhiw2au(epsh,2))+0.5*th*uh/sh*(pyhiw2au(epsh,2)-
1577  & 2.*pyhiw2au(epst,2))+
1578  & 0.125*(sh-12.*sqmq-4.*th*uh/sh)*w3tsui)
1579  b2sutr=sqmq/sqmh**2*(sh*(th-sh)/(sh+th)+2.*uh*th*(th+2.*sh)/
1580  & (sh+th)**2*(pyhiw1au(epsu,1)-pyhiw1au(epsh,1))+(sqmq-sh/4.)*
1581  & (0.5*pyhiw2au(epss,1)+0.5*pyhiw2au(epsh,1)-pyhiw2au(epsu,1)+
1582  & w3sutr)+
1583  & sh**2*(2.*sqmq/(sh+th)**2-0.5/(sh+th))*(pyhiw2au(epsu,1)-
1584  & pyhiw2au(epsh,1))+0.5*uh*th/sh*(pyhiw2au(epsh,1)-
1585  & 2.*pyhiw2au(epsu,1))+
1586  & 0.125*(sh-12.*sqmq-4.*uh*th/sh)*w3ustr)
1587  b2suti=sqmq/sqmh**2*(2.*uh*th*(th+2.*sh)/(sh+th)**2*
1588  & (pyhiw1au(epsu,2)-pyhiw1au(epsh,2))+(sqmq-sh/4.)*
1589  & (0.5*pyhiw2au(epss,2)+0.5*pyhiw2au(epsh,2)-pyhiw2au(epsu,2)+
1590  & w3suti)+
1591  & sh**2*(2.*sqmq/(sh+th)**2-0.5/(sh+th))*(pyhiw2au(epsu,2)-
1592  & pyhiw2au(epsh,2))+0.5*uh*th/sh*(pyhiw2au(epsh,2)-
1593  & 2.*pyhiw2au(epsu,2))+
1594  & 0.125*(sh-12.*sqmq-4.*uh*th/sh)*w3usti)
1595  b2tsur=sqmq/sqmh**2*(th*(uh-th)/(th+uh)+2.*sh*uh*(uh+2.*th)/
1596  & (th+uh)**2*(pyhiw1au(epss,1)-pyhiw1au(epsh,1))+(sqmq-th/4.)*
1597  & (0.5*pyhiw2au(epst,1)+0.5*pyhiw2au(epsh,1)-pyhiw2au(epss,1)+
1598  & w3tsur)+
1599  & th**2*(2.*sqmq/(th+uh)**2-0.5/(th+uh))*(pyhiw2au(epss,1)-
1600  & pyhiw2au(epsh,1))+0.5*sh*uh/th*(pyhiw2au(epsh,1)-
1601  & 2.*pyhiw2au(epss,1))+
1602  & 0.125*(th-12.*sqmq-4.*sh*uh/th)*w3stur)
1603  b2tsui=sqmq/sqmh**2*(2.*sh*uh*(uh+2.*th)/(th+uh)**2*
1604  & (pyhiw1au(epss,2)-pyhiw1au(epsh,2))+(sqmq-th/4.)*
1605  & (0.5*pyhiw2au(epst,2)+0.5*pyhiw2au(epsh,2)-pyhiw2au(epss,2)+
1606  & w3tsui)+
1607  & th**2*(2.*sqmq/(th+uh)**2-0.5/(th+uh))*(pyhiw2au(epss,2)-
1608  & pyhiw2au(epsh,2))+0.5*sh*uh/th*(pyhiw2au(epsh,2)-
1609  & 2.*pyhiw2au(epss,2))+
1610  & 0.125*(th-12.*sqmq-4.*sh*uh/th)*w3stui)
1611  b2tusr=sqmq/sqmh**2*(th*(sh-th)/(th+sh)+2.*uh*sh*(sh+2.*th)/
1612  & (th+sh)**2*(pyhiw1au(epsu,1)-pyhiw1au(epsh,1))+(sqmq-th/4.)*
1613  & (0.5*pyhiw2au(epst,1)+0.5*pyhiw2au(epsh,1)-pyhiw2au(epsu,1)+
1614  & w3tusr)+
1615  & th**2*(2.*sqmq/(th+sh)**2-0.5/(th+sh))*(pyhiw2au(epsu,1)-
1616  & pyhiw2au(epsh,1))+0.5*uh*sh/th*(pyhiw2au(epsh,1)-
1617  & 2.*pyhiw2au(epsu,1))+
1618  & 0.125*(th-12.*sqmq-4.*uh*sh/th)*w3utsr)
1619  b2tusi=sqmq/sqmh**2*(2.*uh*sh*(sh+2.*th)/(th+sh)**2*
1620  & (pyhiw1au(epsu,2)-pyhiw1au(epsh,2))+(sqmq-th/4.)*
1621  & (0.5*pyhiw2au(epst,2)+0.5*pyhiw2au(epsh,2)-pyhiw2au(epsu,2)+
1622  & w3tusi)+
1623  & th**2*(2.*sqmq/(th+sh)**2-0.5/(th+sh))*(pyhiw2au(epsu,2)-
1624  & pyhiw2au(epsh,2))+0.5*uh*sh/th*(pyhiw2au(epsh,2)-
1625  & 2.*pyhiw2au(epsu,2))+
1626  & 0.125*(th-12.*sqmq-4.*uh*sh/th)*w3utsi)
1627  b2ustr=sqmq/sqmh**2*(uh*(th-uh)/(uh+th)+2.*sh*th*(th+2.*uh)/
1628  & (uh+th)**2*(pyhiw1au(epss,1)-pyhiw1au(epsh,1))+(sqmq-uh/4.)*
1629  & (0.5*pyhiw2au(epsu,1)+0.5*pyhiw2au(epsh,1)-pyhiw2au(epss,1)+
1630  & w3ustr)+
1631  & uh**2*(2.*sqmq/(uh+th)**2-0.5/(uh+th))*(pyhiw2au(epss,1)-
1632  & pyhiw2au(epsh,1))+0.5*sh*th/uh*(pyhiw2au(epsh,1)-
1633  & 2.*pyhiw2au(epss,1))+
1634  & 0.125*(uh-12.*sqmq-4.*sh*th/uh)*w3sutr)
1635  b2usti=sqmq/sqmh**2*(2.*sh*th*(th+2.*uh)/(uh+th)**2*
1636  & (pyhiw1au(epss,2)-pyhiw1au(epsh,2))+(sqmq-uh/4.)*
1637  & (0.5*pyhiw2au(epsu,2)+0.5*pyhiw2au(epsh,2)-pyhiw2au(epss,2)+
1638  & w3usti)+
1639  & uh**2*(2.*sqmq/(uh+th)**2-0.5/(uh+th))*(pyhiw2au(epss,2)-
1640  & pyhiw2au(epsh,2))+0.5*sh*th/uh*(pyhiw2au(epsh,2)-
1641  & 2.*pyhiw2au(epss,2))+
1642  & 0.125*(uh-12.*sqmq-4.*sh*th/uh)*w3suti)
1643  b2utsr=sqmq/sqmh**2*(uh*(sh-uh)/(uh+sh)+2.*th*sh*(sh+2.*uh)/
1644  & (uh+sh)**2*(pyhiw1au(epst,1)-pyhiw1au(epsh,1))+(sqmq-uh/4.)*
1645  & (0.5*pyhiw2au(epsu,1)+0.5*pyhiw2au(epsh,1)-pyhiw2au(epst,1)+
1646  & w3utsr)+
1647  & uh**2*(2.*sqmq/(uh+sh)**2-0.5/(uh+sh))*(pyhiw2au(epst,1)-
1648  & pyhiw2au(epsh,1))+0.5*th*sh/uh*(pyhiw2au(epsh,1)-
1649  & 2.*pyhiw2au(epst,1))+
1650  & 0.125*(uh-12.*sqmq-4.*th*sh/uh)*w3tusr)
1651  b2utsi=sqmq/sqmh**2*(2.*th*sh*(sh+2.*uh)/(uh+sh)**2*
1652  & (pyhiw1au(epst,2)-pyhiw1au(epsh,2))+(sqmq-uh/4.)*
1653  & (0.5*pyhiw2au(epsu,2)+0.5*pyhiw2au(epsh,2)-pyhiw2au(epst,2)+
1654  & w3utsi)+
1655  & uh**2*(2.*sqmq/(uh+sh)**2-0.5/(uh+sh))*(pyhiw2au(epst,2)-
1656  & pyhiw2au(epsh,2))+0.5*th*sh/uh*(pyhiw2au(epsh,2)-
1657  & 2.*pyhiw2au(epst,2))+
1658  & 0.125*(uh-12.*sqmq-4.*th*sh/uh)*w3tusi)
1659  b4stur=sqmq/sqmh*(-2./3.+(sqmq/sqmh-1./4.)*(pyhiw2au(epss,1)-
1660  & pyhiw2au(epsh,1)+w3stur))
1661  b4stui=sqmq/sqmh*(sqmq/sqmh-1./4.)*(pyhiw2au(epss,2)-
1662  & pyhiw2au(epsh,2)+w3stui)
1663  b4tusr=sqmq/sqmh*(-2./3.+(sqmq/sqmh-1./4.)*(pyhiw2au(epst,1)-
1664  & pyhiw2au(epsh,1)+w3tusr))
1665  b4tusi=sqmq/sqmh*(sqmq/sqmh-1./4.)*(pyhiw2au(epst,2)-
1666  & pyhiw2au(epsh,2)+w3tusi)
1667  b4ustr=sqmq/sqmh*(-2./3.+(sqmq/sqmh-1./4.)*(pyhiw2au(epsu,1)-
1668  & pyhiw2au(epsh,1)+w3ustr))
1669  b4usti=sqmq/sqmh*(sqmq/sqmh-1./4.)*(pyhiw2au(epsu,2)-
1670  & pyhiw2au(epsh,2)+w3usti)
1671  a2stur=a2stur+b2stur+b2sutr
1672  a2stui=a2stui+b2stui+b2suti
1673  a2ustr=a2ustr+b2ustr+b2utsr
1674  a2usti=a2usti+b2usti+b2utsi
1675  a2tusr=a2tusr+b2tusr+b2tsur
1676  a2tusi=a2tusi+b2tusi+b2tsui
1677  a4stur=a4stur+b4stur+b4ustr+b4tusr
1678  a4stui=a4stui+b4stui+b4usti+b4tusi
1679  760 CONTINUE
1680  facgh=comfac*faca*3./(128.*paru(1)**2)*aem/xw*as**3*
1681  & sqmh/sqmw*sqmh**3/(sh*th*uh)*(a2stur**2+a2stui**2+a2ustr**2+
1682  & a2usti**2+a2tusr**2+a2tusi**2+a4stur**2+a4stui**2)
1683  facgh=facgh*wids(25,2)
1684  IF(kfac(1,21)*kfac(2,21).EQ.0) goto 770
1685  nchn=nchn+1
1686  isig(nchn,1)=21
1687  isig(nchn,2)=21
1688  isig(nchn,3)=1
1689  sigh(nchn)=facgh
1690  770 CONTINUE
1691 
1692  ELSEIF(isub.EQ.114) THEN
1693 C...g + g -> gamma + gamma.
1694  asre=0.
1695  asim=0.
1696  DO 780 i=1,2*mstp(1)
1697  ei=kchg(iabs(i),1)/3.
1698  sqmq=pmas(i,1)**2
1699  epss=4.*sqmq/sh
1700  epst=4.*sqmq/th
1701  epsu=4.*sqmq/uh
1702  IF(epss+abs(epst)+abs(epsu).LT.3.e-6) THEN
1703  a0stur=1.+(th-uh)/sh*log(th/uh)+0.5*(th2+uh2)/sh2*
1704  & (log(th/uh)**2+paru(1)**2)
1705  a0stui=0.
1706  a0tsur=1.+(sh-uh)/th*log(-sh/uh)+0.5*(sh2+uh2)/th2*
1707  & log(-sh/uh)**2
1708  a0tsui=-paru(1)*((sh-uh)/th+(sh2+uh2)/th2*log(-sh/uh))
1709  a0utsr=1.+(th-sh)/uh*log(-th/sh)+0.5*(th2+sh2)/uh2*
1710  & log(-th/sh)**2
1711  a0utsi=paru(1)*((th-sh)/uh+(th2+sh2)/uh2*log(-th/sh))
1712  a1stur=-1.
1713  a1stui=0.
1714  a2stur=-1.
1715  a2stui=0.
1716  ELSE
1717  bestu=0.5*(1.+sqrt(1.+epss*th/uh))
1718  beust=0.5*(1.+sqrt(1.+epsu*sh/th))
1719  betus=0.5*(1.+sqrt(1.+epst*uh/sh))
1720  beuts=bestu
1721  betsu=beust
1722  besut=betus
1723  a0stur=1.+(1.+2.*th/sh)*pyhiw1au(epst,1)+(1.+2.*uh/sh)*
1724  & pyhiw1au(epsu,1)+0.5*((th2+uh2)/sh2-epss)*(pyhiw2au(epst,1)+
1725  & pyhiw2au(epsu,1))-
1726  & 0.25*epst*(1.-0.5*epss)*(pyhii3au(besut,epss,1)+
1727  & pyhii3au(besut,epst,1))-0.25*epsu*(1.-0.5*epss)*
1728  & (pyhii3au(bestu,epss,1)+pyhii3au(bestu,epsu,1))+
1729  & 0.25*(-2.*(th2+uh2)/sh2+4.*epss+epst+epsu+0.5*epst*epsu)*
1730  & (pyhii3au(betsu,epst,1)+pyhii3au(betsu,epsu,1))
1731  a0stui=(1.+2.*th/sh)*pyhiw1au(epst,2)+(1.+2.*uh/sh)*
1732  & pyhiw1au(epsu,2)+0.5*((th2+uh2)/sh2-epss)*(pyhiw2au(epst,2)+
1733  & pyhiw2au(epsu,2))-
1734  & 0.25*epst*(1.-0.5*epss)*(pyhii3au(besut,epss,2)+
1735  & pyhii3au(besut,epst,2))-0.25*epsu*(1.-0.5*epss)*
1736  & (pyhii3au(bestu,epss,2)+pyhii3au(bestu,epsu,2))+
1737  & 0.25*(-2.*(th2+uh2)/sh2+4.*epss+epst+epsu+0.5*epst*epsu)*
1738  & (pyhii3au(betsu,epst,2)+pyhii3au(betsu,epsu,2))
1739  a0tsur=1.+(1.+2.*sh/th)*pyhiw1au(epss,1)+(1.+2.*uh/th)*
1740  & pyhiw1au(epsu,1)+0.5*((sh2+uh2)/th2-epst)*(pyhiw2au(epss,1)+
1741  & pyhiw2au(epsu,1))-
1742  & 0.25*epss*(1.-0.5*epst)*(pyhii3au(betus,epst,1)+
1743  & pyhii3au(betus,epss,1))-0.25*epsu*(1.-0.5*epst)*
1744  & (pyhii3au(betsu,epst,1)+pyhii3au(betsu,epsu,1))+
1745  & 0.25*(-2.*(sh2+uh2)/th2+4.*epst+epss+epsu+0.5*epss*epsu)*
1746  & (pyhii3au(bestu,epss,1)+pyhii3au(bestu,epsu,1))
1747  a0tsui=(1.+2.*sh/th)*pyhiw1au(epss,2)+(1.+2.*uh/th)*
1748  & pyhiw1au(epsu,2)+0.5*((sh2+uh2)/th2-epst)*(pyhiw2au(epss,2)+
1749  & pyhiw2au(epsu,2))-
1750  & 0.25*epss*(1.-0.5*epst)*(pyhii3au(betus,epst,2)+
1751  & pyhii3au(betus,epss,2))-0.25*epsu*(1.-0.5*epst)*
1752  & (pyhii3au(betsu,epst,2)+pyhii3au(betsu,epsu,2))+
1753  & 0.25*(-2.*(sh2+uh2)/th2+4.*epst+epss+epsu+0.5*epss*epsu)*
1754  & (pyhii3au(bestu,epss,2)+pyhii3au(bestu,epsu,2))
1755  a0utsr=1.+(1.+2.*th/uh)*pyhiw1au(epst,1)+(1.+2.*sh/uh)*
1756  & pyhiw1au(epss,1)+0.5*((th2+sh2)/uh2-epsu)*(pyhiw2au(epst,1)+
1757  & pyhiw2au(epss,1))-
1758  & 0.25*epst*(1.-0.5*epsu)*(pyhii3au(beust,epsu,1)+
1759  & pyhii3au(beust,epst,1))-0.25*epss*(1.-0.5*epsu)*
1760  & (pyhii3au(beuts,epsu,1)+pyhii3au(beuts,epss,1))+
1761  & 0.25*(-2.*(th2+sh2)/uh2+4.*epsu+epst+epss+0.5*epst*epss)*
1762  & (pyhii3au(betus,epst,1)+pyhii3au(betus,epss,1))
1763  a0utsi=(1.+2.*th/uh)*pyhiw1au(epst,2)+(1.+2.*sh/uh)*
1764  & pyhiw1au(epss,2)+0.5*((th2+sh2)/uh2-epsu)*(pyhiw2au(epst,2)+
1765  & pyhiw2au(epss,2))-
1766  & 0.25*epst*(1.-0.5*epsu)*(pyhii3au(beust,epsu,2)+
1767  & pyhii3au(beust,epst,2))-0.25*epss*(1.-0.5*epsu)*
1768  & (pyhii3au(beuts,epsu,2)+pyhii3au(beuts,epss,2))+
1769  & 0.25*(-2.*(th2+sh2)/uh2+4.*epsu+epst+epss+0.5*epst*epss)*
1770  & (pyhii3au(betus,epst,2)+pyhii3au(betus,epss,2))
1771  a1stur=-1.-0.25*(epss+epst+epsu)*(pyhiw2au(epss,1)+
1772  & pyhiw2au(epst,1)+pyhiw2au(epsu,1))+0.25*(epsu+0.5*epss*epst)*
1773  & (pyhii3au(besut,epss,1)+pyhii3au(besut,epst,1))+
1774  & 0.25*(epst+0.5*epss*epsu)*(pyhii3au(bestu,epss,1)+
1775  & pyhii3au(bestu,epsu,1))+0.25*(epss+0.5*epst*epsu)*
1776  & (pyhii3au(betsu,epst,1)+pyhii3au(betsu,epsu,1))
1777  a1stui=-0.25*(epss+epst+epsu)*(pyhiw2au(epss,2)+
1778  & pyhiw2au(epst,2)+
1779  & pyhiw2au(epsu,2))+0.25*(epsu+0.5*epss*epst)*
1780  & (pyhii3au(besut,epss,2)+pyhii3au(besut,epst,2))+
1781  & 0.25*(epst+0.5*epss*epsu)*(pyhii3au(bestu,epss,2)+
1782  & pyhii3au(bestu,epsu,2))+0.25*(epss+0.5*epst*epsu)*
1783  & (pyhii3au(betsu,epst,2)+pyhii3au(betsu,epsu,2))
1784  a2stur=-1.+0.125*epss*epst*(pyhii3au(besut,epss,1)+
1785  & pyhii3au(besut,epst,1))+
1786  & 0.125*epss*epsu*(pyhii3au(bestu,epss,1)+
1787  & pyhii3au(bestu,epsu,1))+
1788  & 0.125*epst*epsu*(pyhii3au(betsu,epst,1)+
1789  & pyhii3au(betsu,epsu,1))
1790  a2stui=0.125*epss*epst*(pyhii3au(besut,epss,2)+
1791  & pyhii3au(besut,epst,2))+
1792  & 0.125*epss*epsu*(pyhii3au(bestu,epss,2)+
1793  & pyhii3au(bestu,epsu,2))+
1794  & 0.125*epst*epsu*(pyhii3au(betsu,epst,2)+
1795  & pyhii3au(betsu,epsu,2))
1796  ENDIF
1797  asre=asre+ei**2*(a0stur+a0tsur+a0utsr+4.*a1stur+a2stur)
1798  asim=asim+ei**2*(a0stui+a0tsui+a0utsi+4.*a1stui+a2stui)
1799  780 CONTINUE
1800  facgg=comfac*faca/(8.*paru(1)**2)*as**2*aem**2*(asre**2+asim**2)
1801  IF(kfac(1,21)*kfac(2,21).EQ.0) goto 790
1802  nchn=nchn+1
1803  isig(nchn,1)=21
1804  isig(nchn,2)=21
1805  isig(nchn,3)=1
1806  sigh(nchn)=facgg
1807  790 CONTINUE
1808 
1809  ELSEIF(isub.EQ.115) THEN
1810 C...g + g -> gamma + Z0.
1811 
1812  ELSEIF(isub.EQ.116) THEN
1813 C...g + g -> Z0 + Z0.
1814 
1815  ELSEIF(isub.EQ.117) THEN
1816 C...g + g -> W+ + W-.
1817 
1818  ENDIF
1819 
1820 C...G: 2 -> 3, tree diagrams.
1821 
1822  ELSEIF(isub.LE.140) THEN
1823  IF(isub.EQ.121) THEN
1824 C...g + g -> f + fb + H0.
1825 
1826  ENDIF
1827 
1828 C...H: 2 -> 1, tree diagrams, non-standard model processes.
1829 
1830  ELSEIF(isub.LE.160) THEN
1831  IF(isub.EQ.141) THEN
1832 C...f + fb -> gamma*/Z0/Z'0.
1833  mint(61)=2
1834  CALL pyhiwidt(32,sqrt(sh),wdtp,wdte)
1835  faczp=comfac*aem**2*4./9.
1836  DO 800 i=mina,maxa
1837  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 800
1838  ei=kchg(iabs(i),1)/3.
1839  ai=sign(1.,ei)
1840  vi=ai-4.*ei*xw
1841  api=sign(1.,ei)
1842  vpi=api-4.*ei*xw
1843  nchn=nchn+1
1844  isig(nchn,1)=i
1845  isig(nchn,2)=-i
1846  isig(nchn,3)=1
1847  sigh(nchn)=faczp*(ei**2*vint(111)+ei*vi/(8.*xw*(1.-xw))*
1848  & sh*(sh-sqmz)/((sh-sqmz)**2+gmmz**2)*vint(112)+ei*vpi/(8.*xw*
1849  & (1.-xw))*sh*(sh-sqmzp)/((sh-sqmzp)**2+gmmzp**2)*vint(113)+
1850  & (vi**2+ai**2)/(16.*xw*(1.-xw))**2*sh2/((sh-sqmz)**2+gmmz**2)*
1851  & vint(114)+2.*(vi*vpi+ai*api)/(16.*xw*(1.-xw))**2*sh2*
1852  & ((sh-sqmz)*(sh-sqmzp)+gmmz*gmmzp)/(((sh-sqmz)**2+gmmz**2)*
1853  & ((sh-sqmzp)**2+gmmzp**2))*vint(115)+(vpi**2+api**2)/
1854  & (16.*xw*(1.-xw))**2*sh2/((sh-sqmzp)**2+gmmzp**2)*vint(116))
1855  800 CONTINUE
1856 
1857  ELSEIF(isub.EQ.142) THEN
1858 C...f + fb' -> H+/-.
1859  CALL pyhiwidt(37,sqrt(sh),wdtp,wdte)
1860  fhc=comfac*(aem/xw)**2*1./48.*(sh/sqmw)**2*sh2/
1861  & ((sh-sqmhc)**2+gmmhc**2)
1862 C'''No construction yet for leptons
1863  DO 840 i=1,mstp(54)/2
1864  il=2*i-1
1865  iu=2*i
1866  rmql=pmas(il,1)**2/sh
1867  rmqu=pmas(iu,1)**2/sh
1868  fachc=fhc*((rmql*paru(121)+rmqu/paru(121))*(1.-rmql-rmqu)-
1869  & 4.*rmql*rmqu)/sqrt(max(0.,(1.-rmql-rmqu)**2-4.*rmql*rmqu))
1870  IF(kfac(1,il)*kfac(2,-iu).EQ.0) goto 810
1871  kchhc=(kchg(il,1)-kchg(iu,1))/3
1872  nchn=nchn+1
1873  isig(nchn,1)=il
1874  isig(nchn,2)=-iu
1875  isig(nchn,3)=1
1876  sigh(nchn)=fachc*(wdte(0,1)+wdte(0,(5-kchhc)/2)+wdte(0,4))
1877  810 IF(kfac(1,-il)*kfac(2,iu).EQ.0) goto 820
1878  kchhc=(-kchg(il,1)+kchg(iu,1))/3
1879  nchn=nchn+1
1880  isig(nchn,1)=-il
1881  isig(nchn,2)=iu
1882  isig(nchn,3)=1
1883  sigh(nchn)=fachc*(wdte(0,1)+wdte(0,(5-kchhc)/2)+wdte(0,4))
1884  820 IF(kfac(1,iu)*kfac(2,-il).EQ.0) goto 830
1885  kchhc=(kchg(iu,1)-kchg(il,1))/3
1886  nchn=nchn+1
1887  isig(nchn,1)=iu
1888  isig(nchn,2)=-il
1889  isig(nchn,3)=1
1890  sigh(nchn)=fachc*(wdte(0,1)+wdte(0,(5-kchhc)/2)+wdte(0,4))
1891  830 IF(kfac(1,-iu)*kfac(2,il).EQ.0) goto 840
1892  kchhc=(-kchg(iu,1)+kchg(il,1))/3
1893  nchn=nchn+1
1894  isig(nchn,1)=-iu
1895  isig(nchn,2)=il
1896  isig(nchn,3)=1
1897  sigh(nchn)=fachc*(wdte(0,1)+wdte(0,(5-kchhc)/2)+wdte(0,4))
1898  840 CONTINUE
1899 
1900  ELSEIF(isub.EQ.143) THEN
1901 C...f + fb -> R.
1902  CALL pyhiwidt(40,sqrt(sh),wdtp,wdte)
1903  facr=comfac*(aem/xw)**2*1./9.*sh2/((sh-sqmr)**2+gmmr**2)
1904  DO 860 i=min1,max1
1905  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 860
1906  ia=iabs(i)
1907  DO 850 j=min2,max2
1908  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 850
1909  ja=iabs(j)
1910  IF(i*j.GT.0.OR.iabs(ia-ja).NE.2) goto 850
1911  nchn=nchn+1
1912  isig(nchn,1)=i
1913  isig(nchn,2)=j
1914  isig(nchn,3)=1
1915  sigh(nchn)=facr*(wdte(0,1)+wdte(0,(10-(i+j))/4)+wdte(0,4))
1916  850 CONTINUE
1917  860 CONTINUE
1918 
1919  ENDIF
1920 
1921 C...I: 2 -> 2, tree diagrams, non-standard model processes.
1922 
1923  ELSE
1924  IF(isub.EQ.161) THEN
1925 C...f + g -> f' + H+/- (q + g -> q' + H+/- only).
1926  fhcq=comfac*faca*as*aem/xw*1./24
1927  DO 900 i=1,mstp(54)
1928  iu=i+mod(i,2)
1929  sqmq=pmas(iu,1)**2
1930  fachcq=fhcq/paru(121)*sqmq/sqmw*(sh/(sqmq-uh)+
1931  & 2.*sqmq*(sqmhc-uh)/(sqmq-uh)**2+(sqmq-uh)/sh+
1932  & 2.*sqmq/(sqmq-uh)+2.*(sqmhc-uh)/(sqmq-uh)*(sqmhc-sqmq-sh)/sh)
1933  IF(kfac(1,-i)*kfac(2,21).EQ.0) goto 870
1934  kchhc=isign(1,-kchg(i,1))
1935  nchn=nchn+1
1936  isig(nchn,1)=-i
1937  isig(nchn,2)=21
1938  isig(nchn,3)=1
1939  sigh(nchn)=fachcq*(wdte(0,1)+wdte(0,(5-kchhc)/2)+wdte(0,4))
1940  870 IF(kfac(1,i)*kfac(2,21).EQ.0) goto 880
1941  kchhc=isign(1,kchg(i,1))
1942  nchn=nchn+1
1943  isig(nchn,1)=i
1944  isig(nchn,2)=21
1945  isig(nchn,3)=1
1946  sigh(nchn)=fachcq*(wdte(0,1)+wdte(0,(5-kchhc)/2)+wdte(0,4))
1947  880 IF(kfac(1,21)*kfac(2,-i).EQ.0) goto 890
1948  kchhc=isign(1,-kchg(i,1))
1949  nchn=nchn+1
1950  isig(nchn,1)=21
1951  isig(nchn,2)=-i
1952  isig(nchn,3)=1
1953  sigh(nchn)=fachcq*(wdte(0,1)+wdte(0,(5-kchhc)/2)+wdte(0,4))
1954  890 IF(kfac(1,21)*kfac(2,i).EQ.0) goto 900
1955  kchhc=isign(1,kchg(i,1))
1956  nchn=nchn+1
1957  isig(nchn,1)=21
1958  isig(nchn,2)=i
1959  isig(nchn,3)=1
1960  sigh(nchn)=fachcq*(wdte(0,1)+wdte(0,(5-kchhc)/2)+wdte(0,4))
1961  900 CONTINUE
1962 
1963  ENDIF
1964  ENDIF
1965 
1966 C...Multiply with structure functions.
1967  IF(isub.LE.90.OR.isub.GE.96) THEN
1968  DO 910 ichn=1,nchn
1969  IF(mint(41).EQ.2) THEN
1970  kfl1=isig(ichn,1)
1971  IF(kfl1.EQ.21) kfl1=0
1972  sigh(ichn)=sigh(ichn)*xsfx(1,kfl1)
1973  ENDIF
1974  IF(mint(42).EQ.2) THEN
1975  kfl2=isig(ichn,2)
1976  IF(kfl2.EQ.21) kfl2=0
1977  sigh(ichn)=sigh(ichn)*xsfx(2,kfl2)
1978  ENDIF
1979  910 sigs=sigs+sigh(ichn)
1980  ENDIF
1981 
1982  RETURN
1983  END