Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pysghg.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pysghg.f
1 
2 C*********************************************************************
3 
4 C...PYSGHG
5 C...Subprocess cross sections for Higgs processes,
6 C...except Higgs pairs in PYSGSU, but including WW scattering.
7 C...Auxiliary to PYSIGH.
8 
9  SUBROUTINE pysghg(NCHN,SIGS)
10 
11 C...Double precision and integer declarations
12  IMPLICIT DOUBLE PRECISION(a-h, o-z)
13  IMPLICIT INTEGER(i-n)
14  INTEGER pyk,pychge,pycomp
15 C...Parameter statement to help give large particle numbers.
16  parameter(ksusy1=1000000,ksusy2=2000000,ktechn=3000000,
17  &kexcit=4000000,kdimen=5000000)
18 C...Commonblocks
19  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
20  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
21  common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
22  common/pypars/mstp(200),parp(200),msti(200),pari(200)
23  common/pyint1/mint(400),vint(400)
24  common/pyint2/iset(500),kfpr(500,2),coef(500,20),icol(40,4,2)
25  common/pyint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)
26  common/pyint4/mwid(500),wids(500,5)
27  common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
28  common/pymssm/imss(0:99),rmss(0:99)
29  common/pysgcm/isub,isubsv,mmin1,mmax1,mmin2,mmax2,mmina,mmaxa,
30  &kfac(2,-40:40),comfac,fack,faca,sh,th,uh,sh2,th2,uh2,sqm3,sqm4,
31  &shr,sqpth,taup,be34,cth,x(2),sqmz,sqmw,gmmz,gmmw,
32  &aem,as,xw,xw1,xwc,xwv,poll,polr,polll,polrr
33  SAVE /pydat1/,/pydat2/,/pydat3/,/pypars/,/pyint1/,/pyint2/,
34  &/pyint3/,/pyint4/,/pysubs/,/pymssm/,/pysgcm/
35 C...Local arrays and complex variables
36  dimension wdtp(0:400),wdte(0:400,0:5)
37  COMPLEX*16 a004,a204,a114,a00u,a20u,a11u
38  COMPLEX*16 cigtot,ciztot,f0alp,f1alp,f2alp,f0bet,f1bet,f2bet,fif
39 
40 C...Convert H or A process into equivalent h one
41  ihigg=1
42  kfhigg=25
43  IF(isub.EQ.401.OR.isub.EQ.402) THEN
44  kfhigg=kfpr(isub,1)
45  END IF
46  IF((isub.GE.151.AND.isub.LE.160).OR.(isub.GE.171.AND.
47  &isub.LE.190)) THEN
48  ihigg=2
49  IF(mod(isub-1,10).GE.5) ihigg=3
50  kfhigg=33+ihigg
51  IF(isub.EQ.151.OR.isub.EQ.156) isub=3
52  IF(isub.EQ.152.OR.isub.EQ.157) isub=102
53  IF(isub.EQ.153.OR.isub.EQ.158) isub=103
54  IF(isub.EQ.171.OR.isub.EQ.176) isub=24
55  IF(isub.EQ.172.OR.isub.EQ.177) isub=26
56  IF(isub.EQ.173.OR.isub.EQ.178) isub=123
57  IF(isub.EQ.174.OR.isub.EQ.179) isub=124
58  IF(isub.EQ.181.OR.isub.EQ.186) isub=121
59  IF(isub.EQ.182.OR.isub.EQ.187) isub=122
60  IF(isub.EQ.183.OR.isub.EQ.188) isub=111
61  IF(isub.EQ.184.OR.isub.EQ.189) isub=112
62  IF(isub.EQ.185.OR.isub.EQ.190) isub=113
63  ENDIF
64  sqmh=pmas(kfhigg,1)**2
65  gmmh=pmas(kfhigg,1)*pmas(kfhigg,2)
66 
67 C...Strongly interacting Z_L/W_L model of Dobado, Herrero, Terron
68  IF((mstp(46).GE.3.AND.mstp(46).LE.6).AND.(isub.EQ.71.OR.isub.EQ.
69  &72.OR.isub.EQ.73.OR.isub.EQ.76.OR.isub.EQ.77)) THEN
70 C...Calculate M_R and N_R functions for Higgs-like and QCD-like models
71  IF(mstp(46).LE.4) THEN
72  hdtlh=log(pmas(25,1)/parp(44))
73  hdtmr=(4.5d0*paru(1)/sqrt(3d0)-74d0/9d0)/8d0+hdtlh/12d0
74  hdtnr=-1d0/18d0+hdtlh/6d0
75  ELSE
76  hdtnm=0.125d0*(1d0/(288d0*paru(1)**2)+(parp(47)/parp(45))**2)
77  hdtlq=log(parp(45)/parp(44))
78  hdtmr=-(4d0*paru(1))**2*0.5d0*hdtnm+hdtlq/12d0
79  hdtnr=(4d0*paru(1))**2*hdtnm+hdtlq/6d0
80  ENDIF
81 
82 C...Calculate lowest and next-to-lowest order partial wave amplitudes
83  hdtv=1d0/(16d0*paru(1)*parp(47)**2)
84  a00l=dble(hdtv*sh)
85  a20l=-0.5d0*a00l
86  a11l=a00l/6d0
87  hdtls=log(sh/parp(44)**2)
88  a004=dble((hdtv*sh)**2/(4d0*paru(1)))*
89  & cmplx(dble((176d0*hdtmr+112d0*hdtnr)/3d0+11d0/27d0-
90  & (50d0/9d0)*hdtls),dble(4d0*paru(1)))
91  a204=dble((hdtv*sh)**2/(4d0*paru(1)))*
92  & cmplx(dble(32d0*(hdtmr+2d0*hdtnr)/3d0+25d0/54d0-
93  & (20d0/9d0)*hdtls),dble(paru(1)))
94  a114=dble((hdtv*sh)**2/(6d0*paru(1)))*
95  & cmplx(dble(4d0*(-2d0*hdtmr+hdtnr)-1d0/18d0),dble(paru(1)/6d0))
96 
97 C...Unitarize partial wave amplitudes with Pade or K-matrix method
98  IF(mstp(46).EQ.3.OR.mstp(46).EQ.5) THEN
99  a00u=a00l/(1d0-a004/a00l)
100  a20u=a20l/(1d0-a204/a20l)
101  a11u=a11l/(1d0-a114/a11l)
102  ELSE
103  a00u=(a00l+dble(a004))/(1d0-dcmplx(0.d0,a00l+dble(a004)))
104  a20u=(a20l+dble(a204))/(1d0-dcmplx(0.d0,a20l+dble(a204)))
105  a11u=(a11l+dble(a114))/(1d0-dcmplx(0.d0,a11l+dble(a114)))
106  ENDIF
107  ENDIF
108 
109 C...Differential cross section expressions.
110 
111  IF(isub.LE.60) THEN
112  IF(isub.EQ.3) THEN
113 C...f + fbar -> h0 (or H0, or A0)
114  CALL pywidt(kfhigg,sh,wdtp,wdte)
115  hs=shr*wdtp(0)
116  facbw=4d0*comfac/((sh-sqmh)**2+hs**2)
117  IF(abs(shr-pmas(kfhigg,1)).GT.parp(48)*pmas(kfhigg,2))
118  & facbw=0d0
119  hp=aem/(8d0*xw)*sh/sqmw*sh
120  hf=shr*(wdte(0,1)+wdte(0,2)+wdte(0,4))
121  DO 100 i=mmina,mmaxa
122  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 100
123  ia=iabs(i)
124  rmq=pymrun(ia,sh)**2/sh
125  hi=hp*rmq
126  IF(ia.LE.10) hi=hp*rmq*faca/3d0
127  IF(mstp(4).GE.1.OR.ihigg.GE.2) THEN
128  ikfi=1
129  IF(ia.LE.10.AND.mod(ia,2).EQ.0) ikfi=2
130  IF(ia.GT.10) ikfi=3
131  hi=hi*paru(150+10*ihigg+ikfi)**2
132  IF(imss(1).NE.0.AND.ia.EQ.5) THEN
133  hi=hi/(1d0+rmss(41))**2
134  IF(ihigg.NE.3) THEN
135  hi=hi*(1d0+rmss(41)*paru(152+10*ihigg)/
136  & paru(151+10*ihigg))**2
137  ENDIF
138  ENDIF
139  ENDIF
140  nchn=nchn+1
141  isig(nchn,1)=i
142  isig(nchn,2)=-i
143  isig(nchn,3)=1
144  sigh(nchn)=hi*facbw*hf
145  100 CONTINUE
146 
147  ELSEIF(isub.EQ.5) THEN
148 C...Z0 + Z0 -> h0
149  CALL pywidt(25,sh,wdtp,wdte)
150  hs=shr*wdtp(0)
151  facbw=4d0*comfac/((sh-sqmh)**2+hs**2)
152  IF(abs(shr-pmas(25,1)).GT.parp(48)*pmas(25,2)) facbw=0d0
153  hp=aem/(8d0*xw)*sh/sqmw*sh
154  hf=shr*(wdte(0,1)+wdte(0,2)+wdte(0,4))
155  hi=hp/4d0
156  faci=8d0/(paru(1)**2*xw1)*(aem*xwc)**2
157  DO 120 i=mmin1,mmax1
158  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 120
159  DO 110 j=mmin2,mmax2
160  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 110
161  ei=kchg(iabs(i),1)/3d0
162  ai=sign(1d0,ei)
163  vi=ai-4d0*ei*xwv
164  ej=kchg(iabs(j),1)/3d0
165  aj=sign(1d0,ej)
166  vj=aj-4d0*ej*xwv
167  nchn=nchn+1
168  isig(nchn,1)=i
169  isig(nchn,2)=j
170  isig(nchn,3)=1
171  sigh(nchn)=faci*(vi**2+ai**2)*(vj**2+aj**2)*hi*facbw*hf
172  110 CONTINUE
173  120 CONTINUE
174 
175  ELSEIF(isub.EQ.8) THEN
176 C...W+ + W- -> h0
177  CALL pywidt(25,sh,wdtp,wdte)
178  hs=shr*wdtp(0)
179  facbw=4d0*comfac/((sh-sqmh)**2+hs**2)
180  IF(abs(shr-pmas(25,1)).GT.parp(48)*pmas(25,2)) facbw=0d0
181  hp=aem/(8d0*xw)*sh/sqmw*sh
182  hf=shr*(wdte(0,1)+wdte(0,2)+wdte(0,4))
183  hi=hp/2d0
184  faci=1d0/(4d0*paru(1)**2)*(aem/xw)**2
185  DO 140 i=mmin1,mmax1
186  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 140
187  ei=sign(1d0,dble(i))*kchg(iabs(i),1)
188  DO 130 j=mmin2,mmax2
189  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 130
190  ej=sign(1d0,dble(j))*kchg(iabs(j),1)
191  IF(ei*ej.GT.0d0) goto 130
192  nchn=nchn+1
193  isig(nchn,1)=i
194  isig(nchn,2)=j
195  isig(nchn,3)=1
196  sigh(nchn)=faci*vint(180+i)*vint(180+j)*hi*facbw*hf
197  130 CONTINUE
198  140 CONTINUE
199 
200  ELSEIF(isub.EQ.24) THEN
201 C...f + fbar -> Z0 + h0 (or H0, or A0)
202 C...Propagators: Z0, h0 as simulated in PYOFSH and as desired
203  hbw3=gmmz/((sqm3-sqmz)**2+gmmz**2)
204  CALL pywidt(23,sqm3,wdtp,wdte)
205  gmmz3=sqrt(sqm3)*wdtp(0)
206  hbw3c=gmmz3/((sqm3-sqmz)**2+gmmz3**2)
207  hbw4=gmmh/((sqm4-sqmh)**2+gmmh**2)
208  CALL pywidt(kfhigg,sqm4,wdtp,wdte)
209  gmmh4=sqrt(sqm4)*wdtp(0)
210  hbw4c=gmmh4/((sqm4-sqmh)**2+gmmh4**2)
211  thuh=max(th*uh-sqm3*sqm4,sh*ckin(3)**2)
212  fachz=comfac*(hbw3c/hbw3)*(hbw4c/hbw4)*8d0*(aem*xwc)**2*
213  & (thuh+2d0*sh*sqm3)/((sh-sqmz)**2+gmmz**2)
214  fachz=fachz*wids(23,2)*wids(kfhigg,2)
215  IF(mstp(4).GE.1.OR.ihigg.GE.2) fachz=fachz*
216  & paru(154+10*ihigg)**2
217  DO 150 i=mmina,mmaxa
218  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 150
219  ei=kchg(iabs(i),1)/3d0
220  ai=sign(1d0,ei)
221  vi=ai-4d0*ei*xwv
222  fcoi=1d0
223  IF(iabs(i).LE.10) fcoi=faca/3d0
224  nchn=nchn+1
225  isig(nchn,1)=i
226  isig(nchn,2)=-i
227  isig(nchn,3)=1
228  sigh(nchn)=fachz*fcoi*(vi**2+ai**2)
229  150 CONTINUE
230 
231  ELSEIF(isub.EQ.26) THEN
232 C...f + fbar' -> W+/- + h0 (or H0, or A0)
233 C...Propagators: W+-, h0 as simulated in PYOFSH and as desired
234  hbw3=gmmw/((sqm3-sqmw)**2+gmmw**2)
235  CALL pywidt(24,sqm3,wdtp,wdte)
236  gmmw3=sqrt(sqm3)*wdtp(0)
237  hbw3c=gmmw3/((sqm3-sqmw)**2+gmmw3**2)
238  hbw4=gmmh/((sqm4-sqmh)**2+gmmh**2)
239  CALL pywidt(kfhigg,sqm4,wdtp,wdte)
240  gmmh4=sqrt(sqm4)*wdtp(0)
241  hbw4c=gmmh4/((sqm4-sqmh)**2+gmmh4**2)
242  thuh=max(th*uh-sqm3*sqm4,sh*ckin(3)**2)
243  fachw=comfac*0.125d0*(aem/xw)**2*(thuh+2d0*sh*sqm3)/
244  & ((sh-sqmw)**2+gmmw**2)*(hbw3c/hbw3)*(hbw4c/hbw4)
245  fachw=fachw*wids(kfhigg,2)
246  IF(mstp(4).GE.1.OR.ihigg.GE.2) fachw=fachw*
247  & paru(155+10*ihigg)**2
248  DO 170 i=mmin1,mmax1
249  ia=iabs(i)
250  IF(i.EQ.0.OR.ia.GT.20.OR.kfac(1,i).EQ.0) goto 170
251  DO 160 j=mmin2,mmax2
252  ja=iabs(j)
253  IF(j.EQ.0.OR.ja.GT.20.OR.kfac(1,j).EQ.0) goto 160
254  IF(i*j.GT.0.OR.mod(ia+ja,2).EQ.0) goto 160
255  IF((ia.LE.10.AND.ja.GT.10).OR.(ia.GT.10.AND.ja.LE.10))
256  & goto 160
257  kchw=(kchg(ia,1)*isign(1,i)+kchg(ja,1)*isign(1,j))/3
258  fckm=1d0
259  IF(ia.LE.10) fckm=vckm((ia+1)/2,(ja+1)/2)
260  fcoi=1d0
261  IF(ia.LE.10) fcoi=faca/3d0
262  nchn=nchn+1
263  isig(nchn,1)=i
264  isig(nchn,2)=j
265  isig(nchn,3)=1
266  sigh(nchn)=fachw*fcoi*fckm*wids(24,(5-kchw)/2)
267  160 CONTINUE
268  170 CONTINUE
269 
270  ELSEIF(isub.EQ.32) THEN
271 C...f + g -> f + h0 (q + g -> q + h0 only)
272  fhcq=comfac*faca*as*aem/xw*1d0/24d0
273 C...H propagator: as simulated in PYOFSH and as desired
274  sqmhc=pmas(25,1)**2
275  gmmhc=pmas(25,1)*pmas(25,2)
276  hbw4=gmmhc/((sqm4-sqmhc)**2+gmmhc**2)
277  CALL pywidt(25,sqm4,wdtp,wdte)
278  gmmhcc=sqrt(sqm4)*wdtp(0)
279  hbw4c=gmmhcc/((sqm4-sqmhc)**2+gmmhcc**2)
280  fhcq=fhcq*hbw4c/hbw4
281  DO 190 i=mmina,mmaxa
282  ia=iabs(i)
283  IF(ia.NE.5) goto 190
284  sqml=pymrun(ia,sh)**2
285  sqmq=pmas(ia,1)**2
286  fachcq=fhcq*sqml/sqmw*
287  & (sh/(sqmq-uh)+2d0*sqmq*(sqm4-uh)/(sqmq-uh)**2+(sqmq-uh)/sh-
288  & 2d0*sqmq/(sqmq-uh)+2d0*(sqm4-uh)/(sqmq-uh)*
289  & (sqm4-sqmq-sh)/sh)
290  DO 180 isde=1,2
291  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 180
292  IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 180
293  nchn=nchn+1
294  isig(nchn,isde)=i
295  isig(nchn,3-isde)=21
296  isig(nchn,3)=1
297  sigh(nchn)=fachcq*wids(25,2)
298  180 CONTINUE
299  190 CONTINUE
300  ENDIF
301 
302  ELSEIF(isub.LE.80) THEN
303  IF(isub.EQ.71) THEN
304 C...Z0 + Z0 -> Z0 + Z0
305  IF(sh.LE.4.01d0*sqmz) goto 220
306 
307  IF(mstp(46).LE.2) THEN
308 C...Exact scattering ME:s for on-mass-shell gauge bosons
309  be2=1d0-4d0*sqmz/sh
310  th=-0.5d0*sh*be2*(1d0-cth)
311  uh=-0.5d0*sh*be2*(1d0+cth)
312  IF(max(th,uh).GT.-1d0) goto 220
313  shang=1d0/xw1*sqmw/sqmz*(1d0+be2)**2
314  ashre=(sh-sqmh)/((sh-sqmh)**2+gmmh**2)*shang
315  ashim=-gmmh/((sh-sqmh)**2+gmmh**2)*shang
316  thang=1d0/xw1*sqmw/sqmz*(be2-cth)**2
317  athre=(th-sqmh)/((th-sqmh)**2+gmmh**2)*thang
318  athim=-gmmh/((th-sqmh)**2+gmmh**2)*thang
319  uhang=1d0/xw1*sqmw/sqmz*(be2+cth)**2
320  auhre=(uh-sqmh)/((uh-sqmh)**2+gmmh**2)*uhang
321  auhim=-gmmh/((uh-sqmh)**2+gmmh**2)*uhang
322  faczz=comfac*1d0/(4096d0*paru(1)**2*16d0*xw1**2)*
323  & (aem/xw)**4*(sh/sqmw)**2*(sqmz/sqmw)*sh2
324  IF(mstp(46).LE.0) faczz=faczz*(ashre**2+ashim**2)
325  IF(mstp(46).EQ.1) faczz=faczz*((ashre+athre+auhre)**2+
326  & (ashim+athim+auhim)**2)
327  IF(mstp(46).EQ.2) faczz=0d0
328 
329  ELSE
330 C...Strongly interacting Z_L/W_L model of Dobado, Herrero, Terron
331  faczz=comfac*(aem/(16d0*paru(1)*xw*xw1))**2*(64d0/9d0)*
332  & abs(a00u+2d0*a20u)**2
333  ENDIF
334  faczz=faczz*wids(23,1)
335 
336  DO 210 i=mmin1,mmax1
337  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 210
338  ei=kchg(iabs(i),1)/3d0
339  ai=sign(1d0,ei)
340  vi=ai-4d0*ei*xwv
341  avi=ai**2+vi**2
342  DO 200 j=mmin2,mmax2
343  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 200
344  ej=kchg(iabs(j),1)/3d0
345  aj=sign(1d0,ej)
346  vj=aj-4d0*ej*xwv
347  avj=aj**2+vj**2
348  nchn=nchn+1
349  isig(nchn,1)=i
350  isig(nchn,2)=j
351  isig(nchn,3)=1
352  sigh(nchn)=0.5d0*faczz*avi*avj
353  200 CONTINUE
354  210 CONTINUE
355  220 CONTINUE
356 
357  ELSEIF(isub.EQ.72) THEN
358 C...Z0 + Z0 -> W+ + W-
359  IF(sh.LE.4.01d0*sqmz) goto 250
360 
361  IF(mstp(46).LE.2) THEN
362 C...Exact scattering ME:s for on-mass-shell gauge bosons
363  be2=sqrt((1d0-4d0*sqmw/sh)*(1d0-4d0*sqmz/sh))
364  cth2=cth**2
365  th=-0.5d0*sh*(1d0-2d0*(sqmw+sqmz)/sh-be2*cth)
366  uh=-0.5d0*sh*(1d0-2d0*(sqmw+sqmz)/sh+be2*cth)
367  IF(max(th,uh).GT.-1d0) goto 250
368  shang=4d0*sqrt(sqmw/(sqmz*xw1))*(1d0-2d0*sqmw/sh)*
369  & (1d0-2d0*sqmz/sh)
370  ashre=(sh-sqmh)/((sh-sqmh)**2+gmmh**2)*shang
371  ashim=-gmmh/((sh-sqmh)**2+gmmh**2)*shang
372  atwre=xw1/sqmz*sh/(th-sqmw)*((cth-be2)**2*(3d0/2d0+be2/2d0*
373  & cth-(sqmw+sqmz)/sh+(sqmw-sqmz)**2/(sh*sqmw))+4d0*
374  & ((sqmw+sqmz)/sh*(1d0-3d0*cth2)+8d0*sqmw*sqmz/sh2*
375  & (2d0*cth2-1d0)+4d0*(sqmw**2+sqmz**2)/sh2*cth2+
376  & 2d0*(sqmw+sqmz)/sh*be2*cth))
377  atwim=0d0
378  auwre=xw1/sqmz*sh/(uh-sqmw)*((cth+be2)**2*(3d0/2d0-be2/2d0*
379  & cth-(sqmw+sqmz)/sh+(sqmw-sqmz)**2/(sh*sqmw))+4d0*
380  & ((sqmw+sqmz)/sh*(1d0-3d0*cth2)+8d0*sqmw*sqmz/sh2*
381  & (2d0*cth2-1d0)+4d0*(sqmw**2+sqmz**2)/sh2*cth2-
382  & 2d0*(sqmw+sqmz)/sh*be2*cth))
383  auwim=0d0
384  a4re=2d0*xw1/sqmz*(3d0-cth2-4d0*(sqmw+sqmz)/sh)
385  a4im=0d0
386  facww=comfac*1d0/(4096d0*paru(1)**2*16d0*xw1**2)*
387  & (aem/xw)**4*(sh/sqmw)**2*(sqmz/sqmw)*sh2
388  IF(mstp(46).LE.0) facww=facww*(ashre**2+ashim**2)
389  IF(mstp(46).EQ.1) facww=facww*((ashre+atwre+auwre+a4re)**2+
390  & (ashim+atwim+auwim+a4im)**2)
391  IF(mstp(46).EQ.2) facww=facww*((atwre+auwre+a4re)**2+
392  & (atwim+auwim+a4im)**2)
393 
394  ELSE
395 C...Strongly interacting Z_L/W_L model of Dobado, Herrero, Terron
396  facww=comfac*(aem/(16d0*paru(1)*xw*xw1))**2*(64d0/9d0)*
397  & abs(a00u-a20u)**2
398  ENDIF
399  facww=facww*wids(24,1)
400 
401  DO 240 i=mmin1,mmax1
402  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 240
403  ei=kchg(iabs(i),1)/3d0
404  ai=sign(1d0,ei)
405  vi=ai-4d0*ei*xwv
406  avi=ai**2+vi**2
407  DO 230 j=mmin2,mmax2
408  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 230
409  ej=kchg(iabs(j),1)/3d0
410  aj=sign(1d0,ej)
411  vj=aj-4d0*ej*xwv
412  avj=aj**2+vj**2
413  nchn=nchn+1
414  isig(nchn,1)=i
415  isig(nchn,2)=j
416  isig(nchn,3)=1
417  sigh(nchn)=facww*avi*avj
418  230 CONTINUE
419  240 CONTINUE
420  250 CONTINUE
421 
422  ELSEIF(isub.EQ.73) THEN
423 C...Z0 + W+/- -> Z0 + W+/-
424  IF(sh.LE.2d0*sqmz+2d0*sqmw) goto 280
425 
426  IF(mstp(46).LE.2) THEN
427 C...Exact scattering ME:s for on-mass-shell gauge bosons
428  be2=1d0-2d0*(sqmz+sqmw)/sh+((sqmz-sqmw)/sh)**2
429  ep1=1d0-(sqmz-sqmw)/sh
430  ep2=1d0+(sqmz-sqmw)/sh
431  th=-0.5d0*sh*be2*(1d0-cth)
432  uh=(sqmz-sqmw)**2/sh-0.5d0*sh*be2*(1d0+cth)
433  IF(max(th,uh).GT.-1d0) goto 280
434  thang=(be2-ep1*cth)*(be2-ep2*cth)
435  athre=(th-sqmh)/((th-sqmh)**2+gmmh**2)*thang
436  athim=-gmmh/((th-sqmh)**2+gmmh**2)*thang
437  aswre=-xw1/sqmz*sh/(sh-sqmw)*(-be2*(ep1+ep2)**4*cth+
438  & 1d0/4d0*(be2+ep1*ep2)**2*((ep1-ep2)**2-4d0*be2*cth)+
439  & 2d0*be2*(be2+ep1*ep2)*(ep1+ep2)**2*cth-
440  & 1d0/16d0*sh/sqmw*(ep1**2-ep2**2)**2*(be2+ep1*ep2)**2)
441  aswim=0d0
442  auwre=xw1/sqmz*sh/(uh-sqmw)*(-be2*(ep2+ep1*cth)*
443  & (ep1+ep2*cth)*(be2+ep1*ep2)+be2*(ep2+ep1*cth)*
444  & (be2+ep1*ep2*cth)*(2d0*ep2-ep2*cth+ep1)-
445  & be2*(ep2+ep1*cth)**2*(be2-ep2**2*cth)-1d0/8d0*
446  & (be2+ep1*ep2*cth)**2*((ep1+ep2)**2+2d0*be2*(1d0-cth))+
447  & 1d0/32d0*sh/sqmw*(be2+ep1*ep2*cth)**2*
448  & (ep1**2-ep2**2)**2-be2*(ep1+ep2*cth)*(ep2+ep1*cth)*
449  & (be2+ep1*ep2)+be2*(ep1+ep2*cth)*(be2+ep1*ep2*cth)*
450  & (2d0*ep1-ep1*cth+ep2)-be2*(ep1+ep2*cth)**2*
451  & (be2-ep1**2*cth)-1d0/8d0*(be2+ep1*ep2*cth)**2*
452  & ((ep1+ep2)**2+2d0*be2*(1d0-cth))+1d0/32d0*sh/sqmw*
453  & (be2+ep1*ep2*cth)**2*(ep1**2-ep2**2)**2)
454  auwim=0d0
455  a4re=xw1/sqmz*(ep1**2*ep2**2*(cth**2-1d0)-
456  & 2d0*be2*(ep1**2+ep2**2+ep1*ep2)*cth-2d0*be2*ep1*ep2)
457  a4im=0d0
458  faczw=comfac*1d0/(4096d0*paru(1)**2*4d0*xw1)*(aem/xw)**4*
459  & (sh/sqmw)**2*sqrt(sqmz/sqmw)*sh2
460  IF(mstp(46).LE.0) faczw=0d0
461  IF(mstp(46).EQ.1) faczw=faczw*((athre+aswre+auwre+a4re)**2+
462  & (athim+aswim+auwim+a4im)**2)
463  IF(mstp(46).EQ.2) faczw=faczw*((aswre+auwre+a4re)**2+
464  & (aswim+auwim+a4im)**2)
465 
466  ELSE
467 C...Strongly interacting Z_L/W_L model of Dobado, Herrero, Terron
468  faczw=comfac*aem**2/(64d0*paru(1)**2*xw**2*xw1)*16d0*
469  & abs(a20u+3d0*a11u*dble(cth))**2
470  ENDIF
471  faczw=faczw*wids(23,2)
472 
473  DO 270 i=mmin1,mmax1
474  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 270
475  ei=kchg(iabs(i),1)/3d0
476  ai=sign(1d0,ei)
477  vi=ai-4d0*ei*xwv
478  avi=ai**2+vi**2
479  kchwi=isign(1,kchg(iabs(i),1)*isign(1,i))
480  DO 260 j=mmin2,mmax2
481  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 260
482  ej=kchg(iabs(j),1)/3d0
483  aj=sign(1d0,ej)
484  vj=ai-4d0*ej*xwv
485  avj=aj**2+vj**2
486  kchwj=isign(1,kchg(iabs(j),1)*isign(1,j))
487  nchn=nchn+1
488  isig(nchn,1)=i
489  isig(nchn,2)=j
490  isig(nchn,3)=1
491  sigh(nchn)=faczw*avi*vint(180+j)*wids(24,(5-kchwj)/2)
492  nchn=nchn+1
493  isig(nchn,1)=i
494  isig(nchn,2)=j
495  isig(nchn,3)=2
496  sigh(nchn)=faczw*vint(180+i)*wids(24,(5-kchwi)/2)*avj
497  260 CONTINUE
498  270 CONTINUE
499  280 CONTINUE
500 
501  ELSEIF(isub.EQ.75) THEN
502 C...W+ + W- -> gamma + gamma
503 
504  ELSEIF(isub.EQ.76) THEN
505 C...W+ + W- -> Z0 + Z0
506  IF(sh.LE.4.01d0*sqmz) goto 310
507 
508  IF(mstp(46).LE.2) THEN
509 C...Exact scattering ME:s for on-mass-shell gauge bosons
510  be2=sqrt((1d0-4d0*sqmw/sh)*(1d0-4d0*sqmz/sh))
511  cth2=cth**2
512  th=-0.5d0*sh*(1d0-2d0*(sqmw+sqmz)/sh-be2*cth)
513  uh=-0.5d0*sh*(1d0-2d0*(sqmw+sqmz)/sh+be2*cth)
514  IF(max(th,uh).GT.-1d0) goto 310
515  shang=4d0*sqrt(sqmw/(sqmz*xw1))*(1d0-2d0*sqmw/sh)*
516  & (1d0-2d0*sqmz/sh)
517  ashre=(sh-sqmh)/((sh-sqmh)**2+gmmh**2)*shang
518  ashim=-gmmh/((sh-sqmh)**2+gmmh**2)*shang
519  atwre=xw1/sqmz*sh/(th-sqmw)*((cth-be2)**2*(3d0/2d0+be2/2d0*
520  & cth-(sqmw+sqmz)/sh+(sqmw-sqmz)**2/(sh*sqmw))+4d0*
521  & ((sqmw+sqmz)/sh*(1d0-3d0*cth2)+8d0*sqmw*sqmz/sh2*
522  & (2d0*cth2-1d0)+4d0*(sqmw**2+sqmz**2)/sh2*cth2+
523  & 2d0*(sqmw+sqmz)/sh*be2*cth))
524  atwim=0d0
525  auwre=xw1/sqmz*sh/(uh-sqmw)*((cth+be2)**2*(3d0/2d0-be2/2d0*
526  & cth-(sqmw+sqmz)/sh+(sqmw-sqmz)**2/(sh*sqmw))+4d0*
527  & ((sqmw+sqmz)/sh*(1d0-3d0*cth2)+8d0*sqmw*sqmz/sh2*
528  & (2d0*cth2-1d0)+4d0*(sqmw**2+sqmz**2)/sh2*cth2-
529  & 2d0*(sqmw+sqmz)/sh*be2*cth))
530  auwim=0d0
531  a4re=2d0*xw1/sqmz*(3d0-cth2-4d0*(sqmw+sqmz)/sh)
532  a4im=0d0
533  faczz=comfac*1d0/(4096d0*paru(1)**2)*(aem/xw)**4*
534  & (sh/sqmw)**2*sh2
535  IF(mstp(46).LE.0) faczz=faczz*(ashre**2+ashim**2)
536  IF(mstp(46).EQ.1) faczz=faczz*((ashre+atwre+auwre+a4re)**2+
537  & (ashim+atwim+auwim+a4im)**2)
538  IF(mstp(46).EQ.2) faczz=faczz*((atwre+auwre+a4re)**2+
539  & (atwim+auwim+a4im)**2)
540 
541  ELSE
542 C...Strongly interacting Z_L/W_L model of Dobado, Herrero, Terron
543  faczz=comfac*(aem/(4d0*paru(1)*xw))**2*(64d0/9d0)*
544  & abs(a00u-a20u)**2
545  ENDIF
546  faczz=faczz*wids(23,1)
547 
548  DO 300 i=mmin1,mmax1
549  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 300
550  ei=sign(1d0,dble(i))*kchg(iabs(i),1)
551  DO 290 j=mmin2,mmax2
552  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 290
553  ej=sign(1d0,dble(j))*kchg(iabs(j),1)
554  IF(ei*ej.GT.0d0) goto 290
555  nchn=nchn+1
556  isig(nchn,1)=i
557  isig(nchn,2)=j
558  isig(nchn,3)=1
559  sigh(nchn)=0.5d0*faczz*vint(180+i)*vint(180+j)
560  290 CONTINUE
561  300 CONTINUE
562  310 CONTINUE
563 
564  ELSEIF(isub.EQ.77) THEN
565 C...W+/- + W+/- -> W+/- + W+/-
566  IF(sh.LE.4.01d0*sqmw) goto 340
567 
568  IF(mstp(46).LE.2) THEN
569 C...Exact scattering ME:s for on-mass-shell gauge bosons
570  be2=1d0-4d0*sqmw/sh
571  be4=be2**2
572  cth2=cth**2
573  cth3=cth**3
574  th=-0.5d0*sh*be2*(1d0-cth)
575  uh=-0.5d0*sh*be2*(1d0+cth)
576  IF(max(th,uh).GT.-1d0) goto 340
577  shang=(1d0+be2)**2
578  ashre=(sh-sqmh)/((sh-sqmh)**2+gmmh**2)*shang
579  ashim=-gmmh/((sh-sqmh)**2+gmmh**2)*shang
580  thang=(be2-cth)**2
581  athre=(th-sqmh)/((th-sqmh)**2+gmmh**2)*thang
582  athim=-gmmh/((th-sqmh)**2+gmmh**2)*thang
583  uhang=(be2+cth)**2
584  auhre=(uh-sqmh)/((uh-sqmh)**2+gmmh**2)*uhang
585  auhim=-gmmh/((uh-sqmh)**2+gmmh**2)*uhang
586  sgzang=1d0/sqmw*be2*(3d0-be2)**2*cth
587  asgre=xw*sgzang
588  asgim=0d0
589  aszre=xw1*sh/(sh-sqmz)*sgzang
590  aszim=0d0
591  tgzang=1d0/sqmw*(be2*(4d0-2d0*be2+be4)+be2*(4d0-10d0*be2+
592  & be4)*cth+(2d0-11d0*be2+10d0*be4)*cth2+be2*cth3)
593  atgre=0.5d0*xw*sh/th*tgzang
594  atgim=0d0
595  atzre=0.5d0*xw1*sh/(th-sqmz)*tgzang
596  atzim=0d0
597  ugzang=1d0/sqmw*(be2*(4d0-2d0*be2+be4)-be2*(4d0-10d0*be2+
598  & be4)*cth+(2d0-11d0*be2+10d0*be4)*cth2-be2*cth3)
599  augre=0.5d0*xw*sh/uh*ugzang
600  augim=0d0
601  auzre=0.5d0*xw1*sh/(uh-sqmz)*ugzang
602  auzim=0d0
603  a4are=1d0/sqmw*(1d0+2d0*be2-6d0*be2*cth-cth2)
604  a4aim=0d0
605  a4sre=2d0/sqmw*(1d0+2d0*be2-cth2)
606  a4sim=0d0
607  fww=comfac*1d0/(4096d0*paru(1)**2)*(aem/xw)**4*
608  & (sh/sqmw)**2*sh2
609  IF(mstp(46).LE.0) THEN
610  awware=ashre
611  awwaim=ashim
612  awwsre=0d0
613  awwsim=0d0
614  ELSEIF(mstp(46).EQ.1) THEN
615  awware=ashre+athre+asgre+aszre+atgre+atzre+a4are
616  awwaim=ashim+athim+asgim+aszim+atgim+atzim+a4aim
617  awwsre=-athre-auhre+atgre+atzre+augre+auzre+a4sre
618  awwsim=-athim-auhim+atgim+atzim+augim+auzim+a4sim
619  ELSE
620  awware=asgre+aszre+atgre+atzre+a4are
621  awwaim=asgim+aszim+atgim+atzim+a4aim
622  awwsre=atgre+atzre+augre+auzre+a4sre
623  awwsim=atgim+atzim+augim+auzim+a4sim
624  ENDIF
625  awwa2=awware**2+awwaim**2
626  awws2=awwsre**2+awwsim**2
627 
628  ELSE
629 C...Strongly interacting Z_L/W_L model of Dobado, Herrero, Terron
630  fwwa=comfac*(aem/(4d0*paru(1)*xw))**2*(64d0/9d0)*
631  & abs(a00u+0.5d0*a20u+4.5d0*a11u*dble(cth))**2
632  fwws=comfac*(aem/(4d0*paru(1)*xw))**2*64d0*abs(a20u)**2
633  ENDIF
634 
635  DO 330 i=mmin1,mmax1
636  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 330
637  ei=sign(1d0,dble(i))*kchg(iabs(i),1)
638  DO 320 j=mmin2,mmax2
639  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 320
640  ej=sign(1d0,dble(j))*kchg(iabs(j),1)
641  IF(ei*ej.LT.0d0) THEN
642 C...W+W-
643  IF(mstp(45).EQ.1) goto 320
644  IF(mstp(46).LE.2) facww=fww*awwa2*wids(24,1)
645  IF(mstp(46).GE.3) facww=fwwa*wids(24,1)
646  ELSE
647 C...W+W+/W-W-
648  IF(mstp(45).EQ.2) goto 320
649  IF(mstp(46).LE.2) facww=fww*awws2
650  IF(mstp(46).GE.3) facww=fwws
651  IF(ei.GT.0d0) facww=facww*wids(24,4)
652  IF(ei.LT.0d0) facww=facww*wids(24,5)
653  ENDIF
654  nchn=nchn+1
655  isig(nchn,1)=i
656  isig(nchn,2)=j
657  isig(nchn,3)=1
658  sigh(nchn)=facww*vint(180+i)*vint(180+j)
659  IF(ei*ej.GT.0d0) sigh(nchn)=0.5d0*sigh(nchn)
660  320 CONTINUE
661  330 CONTINUE
662  340 CONTINUE
663  ENDIF
664 
665  ELSEIF(isub.LE.120) THEN
666  IF(isub.EQ.102) THEN
667 C...g + g -> h0 (or H0, or A0)
668  CALL pywidt(kfhigg,sh,wdtp,wdte)
669  wdtp13=0d0
670  DO 345 idc=mdcy(kfhigg,2),mdcy(kfhigg,2)+mdcy(kfhigg,3)-1
671  IF(kfdp(idc,1).EQ.21.AND.kfdp(idc,2).EQ.21.AND.
672  & kfdp(idc,3).EQ.0) wdtp13=pmas(kfhigg,2)*brat(idc)
673  345 CONTINUE
674  IF(wdtp13.EQ.0d0) CALL pyerrm(26,
675  & '(PYSGHG:) did not find Higgs -> g g channel')
676  hs=shr*wdtp(0)
677  hf=shr*(wdte(0,1)+wdte(0,2)+wdte(0,4))
678  facbw=4d0*comfac/((sh-sqmh)**2+hs**2)
679  IF(abs(shr-pmas(kfhigg,1)).GT.parp(48)*pmas(kfhigg,2))
680  & facbw=0d0
681  hi=shr*wdtp13/32d0
682  IF(kfac(1,21)*kfac(2,21).EQ.0) goto 350
683  nchn=nchn+1
684  isig(nchn,1)=21
685  isig(nchn,2)=21
686  isig(nchn,3)=1
687  sigh(nchn)=hi*facbw*hf
688  350 CONTINUE
689 
690  ELSEIF(isub.EQ.103) THEN
691 C...gamma + gamma -> h0 (or H0, or A0)
692  CALL pywidt(kfhigg,sh,wdtp,wdte)
693  wdtp14=0d0
694  DO 355 idc=mdcy(kfhigg,2),mdcy(kfhigg,2)+mdcy(kfhigg,3)-1
695  IF(kfdp(idc,1).EQ.22.AND.kfdp(idc,2).EQ.22.AND.
696  & kfdp(idc,3).EQ.0) wdtp14=pmas(kfhigg,2)*brat(idc)
697  355 CONTINUE
698  IF(wdtp14.EQ.0d0) CALL pyerrm(26,
699  & '(PYSGHG:) did not find Higgs -> gamma gamma channel')
700  hs=shr*wdtp(0)
701  hf=shr*(wdte(0,1)+wdte(0,2)+wdte(0,4))
702  facbw=4d0*comfac/((sh-sqmh)**2+hs**2)
703  IF(abs(shr-pmas(kfhigg,1)).GT.parp(48)*pmas(kfhigg,2))
704  & facbw=0d0
705  hi=shr*wdtp14*2d0
706  IF(kfac(1,22)*kfac(2,22).EQ.0) goto 360
707  nchn=nchn+1
708  isig(nchn,1)=22
709  isig(nchn,2)=22
710  isig(nchn,3)=1
711  sigh(nchn)=hi*facbw*hf
712  360 CONTINUE
713 
714  ELSEIF(isub.EQ.110) THEN
715 C...f + fbar -> gamma + h0
716  thuh=max(th*uh,sh*ckin(3)**2)
717  fachg=comfac*(3d0*aem**4)/(2d0*paru(1)**2*xw*sqmw)*sh*thuh
718  fachg=fachg*wids(kfhigg,2)
719 C...Calculate loop contributions for intermediate gamma* and Z0
720  cigtot=dcmplx(0d0,0d0)
721  ciztot=dcmplx(0d0,0d0)
722  jmax=3*mstp(1)+1
723  DO 370 j=1,jmax
724  IF(j.LE.2*mstp(1)) THEN
725  fnc=1d0
726  ej=kchg(j,1)/3d0
727  aj=sign(1d0,ej+0.1d0)
728  vj=aj-4d0*ej*xwv
729  balp=sqm4/(2d0*pmas(j,1))**2
730  bbet=sh/(2d0*pmas(j,1))**2
731  ELSEIF(j.LE.3*mstp(1)) THEN
732  fnc=3d0
733  jl=2*(j-2*mstp(1))-1
734  ej=kchg(10+jl,1)/3d0
735  aj=sign(1d0,ej+0.1d0)
736  vj=aj-4d0*ej*xwv
737  balp=sqm4/(2d0*pmas(10+jl,1))**2
738  bbet=sh/(2d0*pmas(10+jl,1))**2
739  ELSE
740  balp=sqm4/(2d0*pmas(24,1))**2
741  bbet=sh/(2d0*pmas(24,1))**2
742  ENDIF
743  babi=1d0/(balp-bbet)
744  IF(balp.LT.1d0) THEN
745  f0alp=dcmplx(dble(asin(sqrt(balp))),0d0)
746  f1alp=f0alp**2
747  ELSE
748  f0alp=dcmplx(dble(log(sqrt(balp)+sqrt(balp-1d0))),
749  & -dble(0.5d0*paru(1)))
750  f1alp=-f0alp**2
751  ENDIF
752  f2alp=dble(sqrt(abs(balp-1d0)/balp))*f0alp
753  IF(bbet.LT.1d0) THEN
754  f0bet=dcmplx(dble(asin(sqrt(bbet))),0d0)
755  f1bet=f0bet**2
756  ELSE
757  f0bet=dcmplx(dble(log(sqrt(bbet)+sqrt(bbet-1d0))),
758  & -dble(0.5d0*paru(1)))
759  f1bet=-f0bet**2
760  ENDIF
761  f2bet=dble(sqrt(abs(bbet-1d0)/bbet))*f0bet
762  IF(j.LE.3*mstp(1)) THEN
763  fif=dble(0.5d0*babi)+dble(babi**2)*(dble(0.5d0*(1d0-balp+
764  & bbet))*(f1bet-f1alp)+dble(bbet)*(f2bet-f2alp))
765  cigtot=cigtot+dble(fnc*ej**2)*fif
766  ciztot=ciztot+dble(fnc*ej*vj)*fif
767  ELSE
768  txw=xw/xw1
769  cigtot=cigtot-0.5*(dble(babi*(1.5d0+balp))+dble(babi**2)*
770  & (dble(1.5d0-3d0*balp+4d0*bbet)*(f1bet-f1alp)+
771  & dble(bbet*(2d0*balp+3d0))*(f2bet-f2alp)))
772  ciztot=ciztot-dble(0.5d0*babi*xw1)*(dble(5d0-txw+2d0*balp*
773  & (1d0-txw))*(1d0+dble(2d0*babi*bbet)*(f2bet-f2alp))+
774  & dble(babi*(4d0*bbet*(3d0-txw)-(2d0*balp-1d0)*(5d0-txw)))*
775  & (f1bet-f1alp))
776  ENDIF
777  370 CONTINUE
778  cigtot=cigtot/dble(sh)
779  ciztot=ciztot*dble(xwc)/dcmplx(dble(sh-sqmz),dble(gmmz))
780 C...Loop over initial flavours
781  DO 380 i=mmina,mmaxa
782  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 380
783  ei=kchg(iabs(i),1)/3d0
784  ai=sign(1d0,ei)
785  vi=ai-4d0*ei*xwv
786  fcoi=1d0
787  IF(iabs(i).LE.10) fcoi=faca/3d0
788  nchn=nchn+1
789  isig(nchn,1)=i
790  isig(nchn,2)=-i
791  isig(nchn,3)=1
792  sigh(nchn)=fachg*fcoi*(abs(dble(ei)*cigtot+dble(vi)*
793  & ciztot)**2+ai**2*abs(ciztot)**2)
794  380 CONTINUE
795 
796  ELSEIF(isub.EQ.111) THEN
797 C...f + fbar -> g + h0 (q + qbar -> g + h0 only)
798  IF(mstp(38).NE.0) THEN
799 C...Simple case: only do gg <-> h exactly.
800  CALL pywidt(kfhigg,sqm4,wdtp,wdte)
801  wdtp13=0d0
802  DO 385 idc=mdcy(kfhigg,2),mdcy(kfhigg,2)+mdcy(kfhigg,3)-1
803  IF(kfdp(idc,1).EQ.21.AND.kfdp(idc,2).EQ.21.AND.
804  & kfdp(idc,3).EQ.0) wdtp13=pmas(kfhigg,2)*brat(idc)
805  385 CONTINUE
806  IF(wdtp13.EQ.0d0) CALL pyerrm(26,
807  & '(PYSGHG:) did not find Higgs -> g g channel')
808  facgh=comfac*faca*(2d0/9d0)*as*(wdtp13/sqrt(sqm4))*
809  & (th**2+uh**2)/(sh*sqm4)
810 C...Propagators: as simulated in PYOFSH and as desired
811  hbw4=gmmh/((sqm4-sqmh)**2+gmmh**2)
812  gmmhc=sqrt(sqm4)*wdtp(0)
813  hbw4c=sqrt(sqm4)*(wdte(0,1)+wdte(0,2)+wdte(0,4))/
814  & ((sqm4-sqmh)**2+gmmhc**2)
815  facgh=facgh*hbw4c/hbw4
816  ELSE
817 C...Messy case: do full loop integrals
818  a5stur=0d0
819  a5stui=0d0
820  DO 390 i=1,2*mstp(1)
821  sqmq=pmas(i,1)**2
822  epss=4d0*sqmq/sh
823  epsh=4d0*sqmq/sqmh
824  CALL pywaux(1,epss,w1sr,w1si)
825  CALL pywaux(1,epsh,w1hr,w1hi)
826  CALL pywaux(2,epss,w2sr,w2si)
827  CALL pywaux(2,epsh,w2hr,w2hi)
828  a5stur=a5stur+epsh*(1d0+sh/(th+uh)*(w1sr-w1hr)+
829  & (0.25d0-sqmq/(th+uh))*(w2sr-w2hr))
830  a5stui=a5stui+epsh*(sh/(th+uh)*(w1si-w1hi)+
831  & (0.25d0-sqmq/(th+uh))*(w2si-w2hi))
832  390 CONTINUE
833  facgh=comfac*faca/(144d0*paru(1)**2)*aem/xw*as**3*sqmh/sqmw*
834  & sqmh/sh*(uh**2+th**2)/(uh+th)**2*(a5stur**2+a5stui**2)
835  facgh=facgh*wids(25,2)
836  ENDIF
837  DO 400 i=mmina,mmaxa
838  IF(i.EQ.0.OR.iabs(i).GT.mstp(58).OR.
839  & kfac(1,i)*kfac(2,-i).EQ.0) goto 400
840  nchn=nchn+1
841  isig(nchn,1)=i
842  isig(nchn,2)=-i
843  isig(nchn,3)=1
844  sigh(nchn)=facgh
845  400 CONTINUE
846 
847  ELSEIF(isub.EQ.112) THEN
848 C...f + g -> f + h0 (q + g -> q + h0 only)
849  IF(mstp(38).NE.0) THEN
850 C...Simple case: only do gg <-> h exactly.
851  CALL pywidt(kfhigg,sqm4,wdtp,wdte)
852  wdtp13=0d0
853  DO 405 idc=mdcy(kfhigg,2),mdcy(kfhigg,2)+mdcy(kfhigg,3)-1
854  IF(kfdp(idc,1).EQ.21.AND.kfdp(idc,2).EQ.21.AND.
855  & kfdp(idc,3).EQ.0) wdtp13=pmas(kfhigg,2)*brat(idc)
856  405 CONTINUE
857  IF(wdtp13.EQ.0d0) CALL pyerrm(26,
858  & '(PYSGHG:) did not find Higgs -> g g channel')
859  facqh=comfac*faca*(1d0/12d0)*as*(wdtp13/sqrt(sqm4))*
860  & (sh**2+uh**2)/(-th*sqm4)
861 C...Propagators: as simulated in PYOFSH and as desired
862  hbw4=gmmh/((sqm4-sqmh)**2+gmmh**2)
863  gmmhc=sqrt(sqm4)*wdtp(0)
864  hbw4c=sqrt(sqm4)*(wdte(0,1)+wdte(0,2)+wdte(0,4))/
865  & ((sqm4-sqmh)**2+gmmhc**2)
866  facqh=facqh*hbw4c/hbw4
867  ELSE
868 C...Messy case: do full loop integrals
869  a5tsur=0d0
870  a5tsui=0d0
871  DO 410 i=1,2*mstp(1)
872  sqmq=pmas(i,1)**2
873  epst=4d0*sqmq/th
874  epsh=4d0*sqmq/sqmh
875  CALL pywaux(1,epst,w1tr,w1ti)
876  CALL pywaux(1,epsh,w1hr,w1hi)
877  CALL pywaux(2,epst,w2tr,w2ti)
878  CALL pywaux(2,epsh,w2hr,w2hi)
879  a5tsur=a5tsur+epsh*(1d0+th/(sh+uh)*(w1tr-w1hr)+
880  & (0.25d0-sqmq/(sh+uh))*(w2tr-w2hr))
881  a5tsui=a5tsui+epsh*(th/(sh+uh)*(w1ti-w1hi)+
882  & (0.25d0-sqmq/(sh+uh))*(w2ti-w2hi))
883  410 CONTINUE
884  facqh=comfac*faca/(384d0*paru(1)**2)*aem/xw*as**3*sqmh/sqmw*
885  & sqmh/(-th)*(uh**2+sh**2)/(uh+sh)**2*(a5tsur**2+a5tsui**2)
886  facqh=facqh*wids(25,2)
887  ENDIF
888  DO 430 i=mmina,mmaxa
889  IF(i.EQ.0.OR.iabs(i).GT.mstp(58)) goto 430
890  DO 420 isde=1,2
891  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 420
892  IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 420
893  nchn=nchn+1
894  isig(nchn,isde)=i
895  isig(nchn,3-isde)=21
896  isig(nchn,3)=1
897  sigh(nchn)=facqh
898  420 CONTINUE
899  430 CONTINUE
900 
901  ELSEIF(isub.EQ.113) THEN
902 C...g + g -> g + h0
903  IF(mstp(38).NE.0) THEN
904 C...Simple case: only do gg <-> h exactly.
905  CALL pywidt(kfhigg,sqm4,wdtp,wdte)
906  wdtp13=0d0
907  DO 435 idc=mdcy(kfhigg,2),mdcy(kfhigg,2)+mdcy(kfhigg,3)-1
908  IF(kfdp(idc,1).EQ.21.AND.kfdp(idc,2).EQ.21.AND.
909  & kfdp(idc,3).EQ.0) wdtp13=pmas(kfhigg,2)*brat(idc)
910  435 CONTINUE
911  IF(wdtp13.EQ.0d0) CALL pyerrm(26,
912  & '(PYSGHG:) did not find Higgs -> g g channel')
913  facgh=comfac*faca*(3d0/16d0)*as*(wdtp13/sqrt(sqm4))*
914  & (sh**4+th**4+uh**4+sqm4**4)/(sh*th*uh*sqm4)
915 C...Propagators: as simulated in PYOFSH and as desired
916  hbw4=gmmh/((sqm4-sqmh)**2+gmmh**2)
917  gmmhc=sqrt(sqm4)*wdtp(0)
918  hbw4c=sqrt(sqm4)*(wdte(0,1)+wdte(0,2)+wdte(0,4))/
919  & ((sqm4-sqmh)**2+gmmhc**2)
920  facgh=facgh*hbw4c/hbw4
921  ELSE
922 C...Messy case: do full loop integrals
923  a2stur=0d0
924  a2stui=0d0
925  a2ustr=0d0
926  a2usti=0d0
927  a2tusr=0d0
928  a2tusi=0d0
929  a4stur=0d0
930  a4stui=0d0
931  DO 440 i=1,2*mstp(1)
932  sqmq=pmas(i,1)**2
933  epss=4d0*sqmq/sh
934  epst=4d0*sqmq/th
935  epsu=4d0*sqmq/uh
936  epsh=4d0*sqmq/sqmh
937  IF(epsh.LT.1d-6) goto 440
938  CALL pywaux(1,epss,w1sr,w1si)
939  CALL pywaux(1,epst,w1tr,w1ti)
940  CALL pywaux(1,epsu,w1ur,w1ui)
941  CALL pywaux(1,epsh,w1hr,w1hi)
942  CALL pywaux(2,epss,w2sr,w2si)
943  CALL pywaux(2,epst,w2tr,w2ti)
944  CALL pywaux(2,epsu,w2ur,w2ui)
945  CALL pywaux(2,epsh,w2hr,w2hi)
946  CALL pyi3au(epss,th/uh,y3stur,y3stui)
947  CALL pyi3au(epss,uh/th,y3sutr,y3suti)
948  CALL pyi3au(epst,sh/uh,y3tsur,y3tsui)
949  CALL pyi3au(epst,uh/sh,y3tusr,y3tusi)
950  CALL pyi3au(epsu,sh/th,y3ustr,y3usti)
951  CALL pyi3au(epsu,th/sh,y3utsr,y3utsi)
952  CALL pyi3au(epsh,sqmh/sh*th/uh,yhstur,yhstui)
953  CALL pyi3au(epsh,sqmh/sh*uh/th,yhsutr,yhsuti)
954  CALL pyi3au(epsh,sqmh/th*sh/uh,yhtsur,yhtsui)
955  CALL pyi3au(epsh,sqmh/th*uh/sh,yhtusr,yhtusi)
956  CALL pyi3au(epsh,sqmh/uh*sh/th,yhustr,yhusti)
957  CALL pyi3au(epsh,sqmh/uh*th/sh,yhutsr,yhutsi)
958  w3stur=yhstur-y3stur-y3utsr
959  w3stui=yhstui-y3stui-y3utsi
960  w3sutr=yhsutr-y3sutr-y3tusr
961  w3suti=yhsuti-y3suti-y3tusi
962  w3tsur=yhtsur-y3tsur-y3ustr
963  w3tsui=yhtsui-y3tsui-y3usti
964  w3tusr=yhtusr-y3tusr-y3sutr
965  w3tusi=yhtusi-y3tusi-y3suti
966  w3ustr=yhustr-y3ustr-y3tsur
967  w3usti=yhusti-y3usti-y3tsui
968  w3utsr=yhutsr-y3utsr-y3stur
969  w3utsi=yhutsi-y3utsi-y3stui
970  b2stur=sqmq/sqmh**2*(sh*(uh-sh)/(sh+uh)+2d0*th*uh*
971  & (uh+2d0*sh)/(sh+uh)**2*(w1tr-w1hr)+(sqmq-sh/4d0)*
972  & (0.5d0*w2sr+0.5d0*w2hr-w2tr+w3stur)+sh2*(2d0*sqmq/
973  & (sh+uh)**2-0.5d0/(sh+uh))*(w2tr-w2hr)+0.5d0*th*uh/sh*
974  & (w2hr-2d0*w2tr)+0.125d0*(sh-12d0*sqmq-4d0*th*uh/sh)*w3tsur)
975  b2stui=sqmq/sqmh**2*(2d0*th*uh*(uh+2d0*sh)/(sh+uh)**2*
976  & (w1ti-w1hi)+(sqmq-sh/4d0)*(0.5d0*w2si+0.5d0*w2hi-w2ti+
977  & w3stui)+sh2*(2d0*sqmq/(sh+uh)**2-0.5d0/(sh+uh))*
978  & (w2ti-w2hi)+0.5d0*th*uh/sh*(w2hi-2d0*w2ti)+0.125d0*
979  & (sh-12d0*sqmq-4d0*th*uh/sh)*w3tsui)
980  b2sutr=sqmq/sqmh**2*(sh*(th-sh)/(sh+th)+2d0*uh*th*
981  & (th+2d0*sh)/(sh+th)**2*(w1ur-w1hr)+(sqmq-sh/4d0)*
982  & (0.5d0*w2sr+0.5d0*w2hr-w2ur+w3sutr)+sh2*(2d0*sqmq/
983  & (sh+th)**2-0.5d0/(sh+th))*(w2ur-w2hr)+0.5d0*uh*th/sh*
984  & (w2hr-2d0*w2ur)+0.125d0*(sh-12d0*sqmq-4d0*uh*th/sh)*w3ustr)
985  b2suti=sqmq/sqmh**2*(2d0*uh*th*(th+2d0*sh)/(sh+th)**2*
986  & (w1ui-w1hi)+(sqmq-sh/4d0)*(0.5d0*w2si+0.5d0*w2hi-w2ui+
987  & w3suti)+sh2*(2d0*sqmq/(sh+th)**2-0.5d0/(sh+th))*
988  & (w2ui-w2hi)+0.5d0*uh*th/sh*(w2hi-2d0*w2ui)+0.125d0*
989  & (sh-12d0*sqmq-4d0*uh*th/sh)*w3usti)
990  b2tsur=sqmq/sqmh**2*(th*(uh-th)/(th+uh)+2d0*sh*uh*
991  & (uh+2d0*th)/(th+uh)**2*(w1sr-w1hr)+(sqmq-th/4d0)*
992  & (0.5d0*w2tr+0.5d0*w2hr-w2sr+w3tsur)+th2*(2d0*sqmq/
993  & (th+uh)**2-0.5d0/(th+uh))*(w2sr-w2hr)+0.5d0*sh*uh/th*
994  & (w2hr-2d0*w2sr)+0.125d0*(th-12d0*sqmq-4d0*sh*uh/th)*w3stur)
995  b2tsui=sqmq/sqmh**2*(2d0*sh*uh*(uh+2d0*th)/(th+uh)**2*
996  & (w1si-w1hi)+(sqmq-th/4d0)*(0.5d0*w2ti+0.5d0*w2hi-w2si+
997  & w3tsui)+th2*(2d0*sqmq/(th+uh)**2-0.5d0/(th+uh))*
998  & (w2si-w2hi)+0.5d0*sh*uh/th*(w2hi-2d0*w2si)+0.125d0*
999  & (th-12d0*sqmq-4d0*sh*uh/th)*w3stui)
1000  b2tusr=sqmq/sqmh**2*(th*(sh-th)/(th+sh)+2d0*uh*sh*
1001  & (sh+2d0*th)/(th+sh)**2*(w1ur-w1hr)+(sqmq-th/4d0)*
1002  & (0.5d0*w2tr+0.5d0*w2hr-w2ur+w3tusr)+th2*(2d0*sqmq/
1003  & (th+sh)**2-0.5d0/(th+sh))*(w2ur-w2hr)+0.5d0*uh*sh/th*
1004  & (w2hr-2d0*w2ur)+0.125d0*(th-12d0*sqmq-4d0*uh*sh/th)*w3utsr)
1005  b2tusi=sqmq/sqmh**2*(2d0*uh*sh*(sh+2d0*th)/(th+sh)**2*
1006  & (w1ui-w1hi)+(sqmq-th/4d0)*(0.5d0*w2ti+0.5d0*w2hi-w2ui+
1007  & w3tusi)+th2*(2d0*sqmq/(th+sh)**2-0.5d0/(th+sh))*
1008  & (w2ui-w2hi)+0.5d0*uh*sh/th*(w2hi-2d0*w2ui)+0.125d0*
1009  & (th-12d0*sqmq-4d0*uh*sh/th)*w3utsi)
1010  b2ustr=sqmq/sqmh**2*(uh*(th-uh)/(uh+th)+2d0*sh*th*
1011  & (th+2d0*uh)/(uh+th)**2*(w1sr-w1hr)+(sqmq-uh/4d0)*
1012  & (0.5d0*w2ur+0.5d0*w2hr-w2sr+w3ustr)+uh2*(2d0*sqmq/
1013  & (uh+th)**2-0.5d0/(uh+th))*(w2sr-w2hr)+0.5d0*sh*th/uh*
1014  & (w2hr-2d0*w2sr)+0.125d0*(uh-12d0*sqmq-4d0*sh*th/uh)*w3sutr)
1015  b2usti=sqmq/sqmh**2*(2d0*sh*th*(th+2d0*uh)/(uh+th)**2*
1016  & (w1si-w1hi)+(sqmq-uh/4d0)*(0.5d0*w2ui+0.5d0*w2hi-w2si+
1017  & w3usti)+uh2*(2d0*sqmq/(uh+th)**2-0.5d0/(uh+th))*
1018  & (w2si-w2hi)+0.5d0*sh*th/uh*(w2hi-2d0*w2si)+0.125d0*
1019  & (uh-12d0*sqmq-4d0*sh*th/uh)*w3suti)
1020  b2utsr=sqmq/sqmh**2*(uh*(sh-uh)/(uh+sh)+2d0*th*sh*
1021  & (sh+2d0*uh)/(uh+sh)**2*(w1tr-w1hr)+(sqmq-uh/4d0)*
1022  & (0.5d0*w2ur+0.5d0*w2hr-w2tr+w3utsr)+uh2*(2d0*sqmq/
1023  & (uh+sh)**2-0.5d0/(uh+sh))*(w2tr-w2hr)+0.5d0*th*sh/uh*
1024  & (w2hr-2d0*w2tr)+0.125d0*(uh-12d0*sqmq-4d0*th*sh/uh)*w3tusr)
1025  b2utsi=sqmq/sqmh**2*(2d0*th*sh*(sh+2d0*uh)/(uh+sh)**2*
1026  & (w1ti-w1hi)+(sqmq-uh/4d0)*(0.5d0*w2ui+0.5d0*w2hi-w2ti+
1027  & w3utsi)+uh2*(2d0*sqmq/(uh+sh)**2-0.5d0/(uh+sh))*
1028  & (w2ti-w2hi)+0.5d0*th*sh/uh*(w2hi-2d0*w2ti)+0.125d0*
1029  & (uh-12d0*sqmq-4d0*th*sh/uh)*w3tusi)
1030  b4stur=0.25d0*epsh*(-2d0/3d0+0.25d0*(epsh-1d0)*
1031  & (w2sr-w2hr+w3stur))
1032  b4stui=0.25d0*epsh*0.25d0*(epsh-1d0)*(w2si-w2hi+w3stui)
1033  b4tusr=0.25d0*epsh*(-2d0/3d0+0.25d0*(epsh-1d0)*
1034  & (w2tr-w2hr+w3tusr))
1035  b4tusi=0.25d0*epsh*0.25d0*(epsh-1d0)*(w2ti-w2hi+w3tusi)
1036  b4ustr=0.25d0*epsh*(-2d0/3d0+0.25d0*(epsh-1d0)*
1037  & (w2ur-w2hr+w3ustr))
1038  b4usti=0.25d0*epsh*0.25d0*(epsh-1d0)*(w2ui-w2hi+w3usti)
1039  a2stur=a2stur+b2stur+b2sutr
1040  a2stui=a2stui+b2stui+b2suti
1041  a2ustr=a2ustr+b2ustr+b2utsr
1042  a2usti=a2usti+b2usti+b2utsi
1043  a2tusr=a2tusr+b2tusr+b2tsur
1044  a2tusi=a2tusi+b2tusi+b2tsui
1045  a4stur=a4stur+b4stur+b4ustr+b4tusr
1046  a4stui=a4stui+b4stui+b4usti+b4tusi
1047  440 CONTINUE
1048  facgh=comfac*faca*3d0/(128d0*paru(1)**2)*aem/xw*as**3*
1049  & sqmh/sqmw*sqmh**3/(sh*th*uh)*(a2stur**2+a2stui**2+a2ustr**2+
1050  & a2usti**2+a2tusr**2+a2tusi**2+a4stur**2+a4stui**2)
1051  facgh=facgh*wids(25,2)
1052  ENDIF
1053  IF(kfac(1,21)*kfac(2,21).EQ.0) goto 450
1054  nchn=nchn+1
1055  isig(nchn,1)=21
1056  isig(nchn,2)=21
1057  isig(nchn,3)=1
1058  sigh(nchn)=facgh
1059  450 CONTINUE
1060  ENDIF
1061 
1062  ELSEIF(isub.LE.170) THEN
1063  IF(isub.EQ.121) THEN
1064 C...g + g -> Q + Qbar + h0
1065  IF(kfac(1,21)*kfac(2,21).EQ.0) goto 460
1066  ia=kfpr(isubsv,2)
1067  pmf=pymrun(ia,sh)
1068  facqqh=comfac*(4d0*paru(1)*aem/xw)*(4d0*paru(1)*as)**2*
1069  & (0.5d0*pmf/pmas(24,1))**2
1070  wid2=1d0
1071  IF(ia.EQ.6.OR.ia.EQ.7.OR.ia.EQ.8) wid2=wids(ia,1)
1072  facqqh=facqqh*wid2
1073  IF(mstp(4).GE.1.OR.ihigg.GE.2) THEN
1074  ikfi=1
1075  IF(ia.LE.10.AND.mod(ia,2).EQ.0) ikfi=2
1076  IF(ia.GT.10) ikfi=3
1077  facqqh=facqqh*paru(150+10*ihigg+ikfi)**2
1078  IF(imss(1).NE.0.AND.ia.EQ.5) THEN
1079  facqqh=facqqh/(1d0+rmss(41))**2
1080  IF(ihigg.NE.3) THEN
1081  facqqh=facqqh*(1d0+rmss(41)*paru(152+10*ihigg)/
1082  & paru(151+10*ihigg))**2
1083  ENDIF
1084  ENDIF
1085  ENDIF
1086  CALL pyqqbh(wtqqbh)
1087  CALL pywidt(kfhigg,sh,wdtp,wdte)
1088  hs=shr*wdtp(0)
1089  hf=shr*(wdte(0,1)+wdte(0,2)+wdte(0,4))
1090  facbw=(1d0/paru(1))*vint(2)*hf/((sh-sqmh)**2+hs**2)
1091  IF(abs(shr-pmas(kfhigg,1)).GT.parp(48)*pmas(kfhigg,2))
1092  & facbw=0d0
1093  nchn=nchn+1
1094  isig(nchn,1)=21
1095  isig(nchn,2)=21
1096  isig(nchn,3)=1
1097  sigh(nchn)=facqqh*wtqqbh*facbw
1098  460 CONTINUE
1099 
1100  ELSEIF(isub.EQ.122) THEN
1101 C...q + qbar -> Q + Qbar + h0
1102  ia=kfpr(isubsv,2)
1103  pmf=pymrun(ia,sh)
1104  facqqh=comfac*(4d0*paru(1)*aem/xw)*(4d0*paru(1)*as)**2*
1105  & (0.5d0*pmf/pmas(24,1))**2
1106  wid2=1d0
1107  IF(ia.EQ.6.OR.ia.EQ.7.OR.ia.EQ.8) wid2=wids(ia,1)
1108  facqqh=facqqh*wid2
1109  IF(mstp(4).GE.1.OR.ihigg.GE.2) THEN
1110  ikfi=1
1111  IF(ia.LE.10.AND.mod(ia,2).EQ.0) ikfi=2
1112  IF(ia.GT.10) ikfi=3
1113  facqqh=facqqh*paru(150+10*ihigg+ikfi)**2
1114  IF(imss(1).NE.0.AND.ia.EQ.5) THEN
1115  facqqh=facqqh/(1d0+rmss(41))**2
1116  IF(ihigg.NE.3) THEN
1117  facqqh=facqqh*(1d0+rmss(41)*paru(152+10*ihigg)/
1118  & paru(151+10*ihigg))**2
1119  ENDIF
1120  ENDIF
1121  ENDIF
1122  CALL pyqqbh(wtqqbh)
1123  CALL pywidt(kfhigg,sh,wdtp,wdte)
1124  hs=shr*wdtp(0)
1125  hf=shr*(wdte(0,1)+wdte(0,2)+wdte(0,4))
1126  facbw=(1d0/paru(1))*vint(2)*hf/((sh-sqmh)**2+hs**2)
1127  IF(abs(shr-pmas(kfhigg,1)).GT.parp(48)*pmas(kfhigg,2))
1128  & facbw=0d0
1129  DO 470 i=mmina,mmaxa
1130  IF(i.EQ.0.OR.iabs(i).GT.mstp(58).OR.
1131  & kfac(1,i)*kfac(2,-i).EQ.0) goto 470
1132  nchn=nchn+1
1133  isig(nchn,1)=i
1134  isig(nchn,2)=-i
1135  isig(nchn,3)=1
1136  sigh(nchn)=facqqh*wtqqbh*facbw
1137  470 CONTINUE
1138 
1139  ELSEIF(isub.EQ.123) THEN
1140 C...f + f' -> f + f' + h0 (or H0, or A0) (Z0 + Z0 -> h0 as
1141 C...inner process)
1142  facnor=comfac*(4d0*paru(1)*aem/(xw*xw1))**3*sqmz/32d0
1143  IF(mstp(4).GE.1.OR.ihigg.GE.2) facnor=facnor*
1144  & paru(154+10*ihigg)**2
1145  facprp=1d0/((vint(215)-vint(204)**2)*
1146  & (vint(216)-vint(209)**2))**2
1147  faczz1=facnor*facprp*(0.5d0*taup*vint(2))*vint(219)
1148  faczz2=facnor*facprp*vint(217)*vint(218)
1149  CALL pywidt(kfhigg,sh,wdtp,wdte)
1150  hs=shr*wdtp(0)
1151  hf=shr*(wdte(0,1)+wdte(0,2)+wdte(0,4))
1152  facbw=(1d0/paru(1))*vint(2)*hf/((sh-sqmh)**2+hs**2)
1153  IF(abs(shr-pmas(kfhigg,1)).GT.parp(48)*pmas(kfhigg,2))
1154  & facbw=0d0
1155  DO 490 i=mmin1,mmax1
1156  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 490
1157  ia=iabs(i)
1158  DO 480 j=mmin2,mmax2
1159  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 480
1160  ja=iabs(j)
1161  ei=kchg(ia,1)*isign(1,i)/3d0
1162  ai=sign(1d0,kchg(ia,1)+0.5d0)*isign(1,i)
1163  vi=ai-4d0*ei*xwv
1164  ej=kchg(ja,1)*isign(1,j)/3d0
1165  aj=sign(1d0,kchg(ja,1)+0.5d0)*isign(1,j)
1166  vj=aj-4d0*ej*xwv
1167  faclr1=(vi**2+ai**2)*(vj**2+aj**2)+4d0*vi*ai*vj*aj
1168  faclr2=(vi**2+ai**2)*(vj**2+aj**2)-4d0*vi*ai*vj*aj
1169  nchn=nchn+1
1170  isig(nchn,1)=i
1171  isig(nchn,2)=j
1172  isig(nchn,3)=1
1173  sigh(nchn)=(faclr1*faczz1+faclr2*faczz2)*facbw
1174  480 CONTINUE
1175  490 CONTINUE
1176 
1177  ELSEIF(isub.EQ.124) THEN
1178 C...f + f' -> f" + f"' + h0 (or H0, or A0) (W+ + W- -> h0 as
1179 C...inner process)
1180  facnor=comfac*(4d0*paru(1)*aem/xw)**3*sqmw
1181  IF(mstp(4).GE.1.OR.ihigg.GE.2) facnor=facnor*
1182  & paru(155+10*ihigg)**2
1183  facprp=1d0/((vint(215)-vint(204)**2)*
1184  & (vint(216)-vint(209)**2))**2
1185  facww=facnor*facprp*(0.5d0*taup*vint(2))*vint(219)
1186  CALL pywidt(kfhigg,sh,wdtp,wdte)
1187  hs=shr*wdtp(0)
1188  hf=shr*(wdte(0,1)+wdte(0,2)+wdte(0,4))
1189  facbw=(1d0/paru(1))*vint(2)*hf/((sh-sqmh)**2+hs**2)
1190  IF(abs(shr-pmas(kfhigg,1)).GT.parp(48)*pmas(kfhigg,2))
1191  & facbw=0d0
1192  DO 510 i=mmin1,mmax1
1193  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 510
1194  ei=sign(1d0,dble(i))*kchg(iabs(i),1)
1195  DO 500 j=mmin2,mmax2
1196  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 500
1197  ej=sign(1d0,dble(j))*kchg(iabs(j),1)
1198  IF(ei*ej.GT.0d0) goto 500
1199  faclr=vint(180+i)*vint(180+j)
1200  nchn=nchn+1
1201  isig(nchn,1)=i
1202  isig(nchn,2)=j
1203  isig(nchn,3)=1
1204  sigh(nchn)=faclr*facww*facbw
1205  500 CONTINUE
1206  510 CONTINUE
1207 
1208  ELSEIF(isub.EQ.143) THEN
1209 C...f + fbar' -> H+/-
1210  sqmhc=pmas(37,1)**2
1211  CALL pywidt(37,sh,wdtp,wdte)
1212  hs=shr*wdtp(0)
1213  facbw=4d0*comfac/((sh-sqmhc)**2+hs**2)
1214  hp=aem/(8d0*xw)*sh/sqmw*sh
1215  DO 530 i=mmin1,mmax1
1216  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 530
1217  ia=iabs(i)
1218  im=(mod(ia,10)+1)/2
1219  DO 520 j=mmin2,mmax2
1220  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 520
1221  ja=iabs(j)
1222  jm=(mod(ja,10)+1)/2
1223  IF(i*j.GT.0.OR.ia.EQ.ja.OR.im.NE.jm) goto 520
1224  IF((ia.LE.10.AND.ja.GT.10).OR.(ia.GT.10.AND.ja.LE.10))
1225  & goto 520
1226  IF(mod(ia,2).EQ.0) THEN
1227  iu=ia
1228  il=ja
1229  ELSE
1230  iu=ja
1231  il=ia
1232  ENDIF
1233  rml=pymrun(il,sh)**2/sh
1234  rmu=pymrun(iu,sh)**2/sh
1235  hi=hp*(rml*paru(141)**2+rmu/paru(141)**2)
1236  IF(ia.LE.10) hi=hi*faca/3d0
1237  kchhc=(kchg(ia,1)*isign(1,i)+kchg(ja,1)*isign(1,j))/3
1238  hf=shr*(wdte(0,1)+wdte(0,(5-kchhc)/2)+wdte(0,4))
1239  nchn=nchn+1
1240  isig(nchn,1)=i
1241  isig(nchn,2)=j
1242  isig(nchn,3)=1
1243  sigh(nchn)=hi*facbw*hf
1244  520 CONTINUE
1245  530 CONTINUE
1246 
1247  ELSEIF(isub.EQ.161) THEN
1248 C...f + g -> f' + H+/- (b + g -> t + H+/- only)
1249 C...(choice of only b and t to avoid kinematics problems)
1250  fhcq=comfac*faca*as*aem/xw*1d0/24
1251 C...H propagator: as simulated in PYOFSH and as desired
1252  sqmhc=pmas(37,1)**2
1253  gmmhc=pmas(37,1)*pmas(37,2)
1254  hbw4=gmmhc/((sqm4-sqmhc)**2+gmmhc**2)
1255  CALL pywidt(37,sqm4,wdtp,wdte)
1256  gmmhcc=sqrt(sqm4)*wdtp(0)
1257  hbw4c=gmmhcc/((sqm4-sqmhc)**2+gmmhcc**2)
1258  fhcq=fhcq*hbw4c/hbw4
1259  q2rm=sh
1260  IF(mstp(32).EQ.12) q2rm=parp(194)
1261  DO 550 i=mmina,mmaxa
1262  ia=iabs(i)
1263  IF(ia.NE.5) goto 550
1264  sqml=pymrun(ia,q2rm)**2
1265  iua=ia+mod(ia,2)
1266  sqmq=pymrun(iua,q2rm)**2
1267  fachcq=fhcq*(sqml*paru(141)**2+sqmq/paru(141)**2)/sqmw*
1268  & (sh/(sqmq-uh)+2d0*sqmq*(sqmhc-uh)/(sqmq-uh)**2+(sqmq-uh)/sh-
1269  & 2d0*sqmq/(sqmq-uh)+2d0*(sqmhc-uh)/(sqmq-uh)*
1270  & (sqmhc-sqmq-sh)/sh)
1271  kchhc=isign(1,kchg(ia,1)*isign(1,i))
1272  DO 540 isde=1,2
1273  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 540
1274  IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 540
1275  nchn=nchn+1
1276  isig(nchn,isde)=i
1277  isig(nchn,3-isde)=21
1278  isig(nchn,3)=1
1279  sigh(nchn)=fachcq*wids(37,(5-kchhc)/2)
1280  IF(iua.EQ.6) sigh(nchn)=sigh(nchn)*wids(6,(5+kchhc)/2)
1281  540 CONTINUE
1282  550 CONTINUE
1283  ENDIF
1284 
1285  ELSEIF(isub.LE.402) THEN
1286  IF(isub.EQ.401) THEN
1287 C... g + g -> t + bbar + H-
1288  IF(kfac(1,21)*kfac(2,21).EQ.0) goto 560
1289  ia=kfpr(isubsv,2)
1290  CALL pystbh(wttbh)
1291  CALL pywidt(kfhigg,sh,wdtp,wdte)
1292  hs=shr*wdtp(0)
1293  facbw=(1d0/paru(1))*vint(2)*hs/((sh-sqmh)**2+hs**2)
1294  IF(abs(shr-pmas(kfhigg,1)).GT.parp(48)*pmas(kfhigg,2))
1295  & facbw=0d0
1296  nchn=nchn+1
1297  isig(nchn,1)=21
1298  isig(nchn,2)=21
1299  isig(nchn,3)=1
1300  sigh(nchn)=2d0*comfac*wttbh*facbw
1301 c Since we don't know yet if H+ or H-, assume H+
1302 c when calculating suppression due to closed channels.
1303  sigh(nchn)=sigh(nchn)*wids(37,2)*wids(6,3)
1304  IF(abs(wids(37,2)-wids(37,3))
1305  & .GE.1d-6*(wids(37,2)+wids(37,3)).OR.
1306  & abs(wids(6,2)-wids(6,3))
1307  & .GE.1d-6*(wids(6,2)+wids(6,3))) THEN
1308  WRITE(*,*)'Error: Process 401 cannot handle different'
1309  WRITE(*,*)'decays for H+ and H- or t and tbar.'
1310  WRITE(*,*)'Execution stopped.'
1311  CALL pystop(108)
1312  END IF
1313  560 CONTINUE
1314 
1315  ELSEIF(isub.EQ.402) THEN
1316 C... q + qbar -> t + bbar + H-
1317  ia=kfpr(isubsv,2)
1318  CALL pystbh(wttbh)
1319  CALL pywidt(kfhigg,sh,wdtp,wdte)
1320  hs=shr*wdtp(0)
1321  facbw=(1d0/paru(1))*vint(2)*hs/((sh-sqmh)**2+hs**2)
1322  IF(abs(shr-pmas(kfhigg,1)).GT.parp(48)*pmas(kfhigg,2))
1323  & facbw=0d0
1324  DO 570 i=mmina,mmaxa
1325  IF(i.EQ.0.OR.iabs(i).GT.mstp(58).OR.
1326  & kfac(1,i)*kfac(2,-i).EQ.0) goto 570
1327  nchn=nchn+1
1328  isig(nchn,1)=i
1329  isig(nchn,2)=-i
1330  isig(nchn,3)=1
1331  sigh(nchn)=2d0*comfac*wttbh*facbw
1332 c Since we don't know yet if H+ or H-, assume H+
1333 c when calculating suppression due to closed channels.
1334  sigh(nchn)=sigh(nchn)*wids(37,2)*wids(6,3)
1335  IF(abs(wids(37,2)-wids(37,3))/(wids(37,2)+wids(37,3))
1336  & .GE.1d-6.OR.
1337  & abs(wids(6,2)-wids(6,3))/(wids(6,2)+wids(6,3))
1338  & .GE.1d-6) THEN
1339  WRITE(*,*)'Error: Process 402 cannot handle different'
1340  WRITE(*,*)'decays for H+ and H- or t and tbar.'
1341  WRITE(*,*)'Execution stopped.'
1342  CALL pystop(108)
1343  END IF
1344  570 CONTINUE
1345  ENDIF
1346  ENDIF
1347 
1348  RETURN
1349  END