Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pysgqc.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pysgqc.f
1 
2 C*********************************************************************
3 
4 C...PYSGQC
5 C...Subprocess cross sections for QCD processes,
6 C...including photons.
7 C...Auxiliary to PYSIGH.
8 
9  SUBROUTINE pysgqc(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/pyint7/sigt(0:6,0:6,0:5)
28  common/pysgcm/isub,isubsv,mmin1,mmax1,mmin2,mmax2,mmina,mmaxa,
29  &kfac(2,-40:40),comfac,fack,faca,sh,th,uh,sh2,th2,uh2,sqm3,sqm4,
30  &shr,sqpth,taup,be34,cth,x(2),sqmz,sqmw,gmmz,gmmw,
31  &aem,as,xw,xw1,xwc,xwv,poll,polr,polll,polrr
32  SAVE /pydat1/,/pydat2/,/pydat3/,/pypars/,/pyint1/,/pyint2/,
33  &/pyint3/,/pyint4/,/pyint7/,/pysgcm/
34 C...Local arrays
35  dimension wdtp(0:400),wdte(0:400,0:5)
36 
37 C...Differential cross section expressions.
38 
39  IF(isub.LE.20) THEN
40  IF(isub.EQ.10) THEN
41 C...f + f' -> f + f' (gamma/Z/W exchange)
42  facggf=comfac*aem**2*2d0*(sh2+uh2)/th2
43  facgzf=comfac*aem**2*xwc*4d0*sh2/(th*(th-sqmz))
44  faczzf=comfac*(aem*xwc)**2*2d0*sh2/(th-sqmz)**2
45  facwwf=comfac*(0.5d0*aem/xw)**2*sh2/(th-sqmw)**2
46  DO 110 i=mmin1,mmax1
47  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 110
48  ia=iabs(i)
49  DO 100 j=mmin2,mmax2
50  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 100
51  ja=iabs(j)
52 C...Electroweak couplings
53  ei=kchg(ia,1)*isign(1,i)/3d0
54  ai=sign(1d0,kchg(ia,1)+0.5d0)*isign(1,i)
55  vi=ai-4d0*ei*xwv
56  ej=kchg(ja,1)*isign(1,j)/3d0
57  aj=sign(1d0,kchg(ja,1)+0.5d0)*isign(1,j)
58  vj=aj-4d0*ej*xwv
59  epsij=isign(1,i*j)
60 C...gamma/Z exchange, only gamma exchange, or only Z exchange
61  IF(mstp(21).GE.1.AND.mstp(21).LE.4) THEN
62  IF(mstp(21).EQ.1.OR.mstp(21).EQ.4) THEN
63  facncf=facggf*ei**2*ej**2+facgzf*ei*ej*
64  & (vi*vj*(1d0+uh2/sh2)+ai*aj*epsij*(1d0-uh2/sh2))+
65  & faczzf*((vi**2+ai**2)*(vj**2+aj**2)*(1d0+uh2/sh2)+
66  & 4d0*vi*vj*ai*aj*epsij*(1d0-uh2/sh2))
67  ELSEIF(mstp(21).EQ.2) THEN
68  facncf=facggf*ei**2*ej**2
69  ELSE
70  facncf=faczzf*((vi**2+ai**2)*(vj**2+aj**2)*
71  & (1d0+uh2/sh2)+4d0*vi*vj*ai*aj*epsij*(1d0-uh2/sh2))
72  ENDIF
73 C...Extrafactor 2 for only one incoming neutrino spin state.
74  IF(ia.GT.10.AND.mod(ia,2).EQ.0) facncf=2d0*facncf
75  IF(ja.GT.10.AND.mod(ja,2).EQ.0) facncf=2d0*facncf
76  nchn=nchn+1
77  isig(nchn,1)=i
78  isig(nchn,2)=j
79  isig(nchn,3)=1
80  sigh(nchn)=facncf
81  ENDIF
82 C...W exchange
83  IF((mstp(21).EQ.1.OR.mstp(21).EQ.5).AND.ai*aj.LT.0d0) THEN
84  facccf=facwwf*vint(180+i)*vint(180+j)
85  IF(epsij.LT.0d0) facccf=facccf*uh2/sh2
86  IF(ia.GT.10.AND.mod(ia,2).EQ.0) facccf=2d0*facccf
87  IF(ja.GT.10.AND.mod(ja,2).EQ.0) facccf=2d0*facccf
88  nchn=nchn+1
89  isig(nchn,1)=i
90  isig(nchn,2)=j
91  isig(nchn,3)=2
92  sigh(nchn)=facccf
93  ENDIF
94  100 CONTINUE
95  110 CONTINUE
96 
97  ELSEIF(isub.EQ.11) THEN
98 C...f + f' -> f + f' (g exchange)
99  facqq1=comfac*as**2*4d0/9d0*(sh2+uh2)/th2
100  facqqb=comfac*as**2*4d0/9d0*((sh2+uh2)/th2*faca-
101  & mstp(34)*2d0/3d0*uh2/(sh*th))
102  facqq2=comfac*as**2*4d0/9d0*((sh2+th2)/uh2-
103  & mstp(34)*2d0/3d0*sh2/(th*uh))
104  DO 130 i=mmin1,mmax1
105  ia=iabs(i)
106  IF(i.EQ.0.OR.ia.GT.mstp(58).OR.kfac(1,i).EQ.0) goto 130
107  DO 120 j=mmin2,mmax2
108  ja=iabs(j)
109  IF(j.EQ.0.OR.ja.GT.mstp(58).OR.kfac(2,j).EQ.0) goto 120
110  nchn=nchn+1
111  isig(nchn,1)=i
112  isig(nchn,2)=j
113  isig(nchn,3)=1
114  sigh(nchn)=facqq1
115  IF(i.EQ.-j) sigh(nchn)=facqqb
116  IF(i.EQ.j) THEN
117  sigh(nchn)=0.5d0*sigh(nchn)
118  nchn=nchn+1
119  isig(nchn,1)=i
120  isig(nchn,2)=j
121  isig(nchn,3)=2
122  sigh(nchn)=0.5d0*facqq2
123  ENDIF
124  120 CONTINUE
125  130 CONTINUE
126 
127  ELSEIF(isub.EQ.12) THEN
128 C...f + fbar -> f' + fbar' (q + qbar -> q' + qbar' only)
129  CALL pywidt(21,sh,wdtp,wdte)
130  facqqb=comfac*as**2*4d0/9d0*(th2+uh2)/sh2*
131  & (wdte(0,1)+wdte(0,2)+wdte(0,4))
132  DO 140 i=mmina,mmaxa
133  IF(i.EQ.0.OR.iabs(i).GT.mstp(58).OR.
134  & kfac(1,i)*kfac(2,-i).EQ.0) goto 140
135  nchn=nchn+1
136  isig(nchn,1)=i
137  isig(nchn,2)=-i
138  isig(nchn,3)=1
139  sigh(nchn)=facqqb
140  140 CONTINUE
141 
142  ELSEIF(isub.EQ.13) THEN
143 C...f + fbar -> g + g (q + qbar -> g + g only)
144  facgg1=comfac*as**2*32d0/27d0*(uh/th-(2d0+mstp(34)*1d0/4d0)*
145  & uh2/sh2)
146  facgg2=comfac*as**2*32d0/27d0*(th/uh-(2d0+mstp(34)*1d0/4d0)*
147  & th2/sh2)
148  DO 150 i=mmina,mmaxa
149  IF(i.EQ.0.OR.iabs(i).GT.mstp(58).OR.
150  & kfac(1,i)*kfac(2,-i).EQ.0) goto 150
151  nchn=nchn+1
152  isig(nchn,1)=i
153  isig(nchn,2)=-i
154  isig(nchn,3)=1
155  sigh(nchn)=0.5d0*facgg1
156  nchn=nchn+1
157  isig(nchn,1)=i
158  isig(nchn,2)=-i
159  isig(nchn,3)=2
160  sigh(nchn)=0.5d0*facgg2
161  150 CONTINUE
162 
163  ELSEIF(isub.EQ.14) THEN
164 C...f + fbar -> g + gamma (q + qbar -> g + gamma only)
165  facgg=comfac*as*aem*8d0/9d0*(th2+uh2)/(th*uh)
166  DO 160 i=mmina,mmaxa
167  IF(i.EQ.0.OR.iabs(i).GT.mstp(58).OR.
168  & kfac(1,i)*kfac(2,-i).EQ.0) goto 160
169  ei=kchg(iabs(i),1)/3d0
170  nchn=nchn+1
171  isig(nchn,1)=i
172  isig(nchn,2)=-i
173  isig(nchn,3)=1
174  sigh(nchn)=facgg*ei**2
175  160 CONTINUE
176 
177  ELSEIF(isub.EQ.18) THEN
178 C...f + fbar -> gamma + gamma
179  facgg=comfac*aem**2*2d0*(th2+uh2)/(th*uh)
180  DO 170 i=mmina,mmaxa
181  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 170
182  ei=kchg(iabs(i),1)/3d0
183  fcoi=1d0
184  IF(iabs(i).LE.10) fcoi=faca/3d0
185  nchn=nchn+1
186  isig(nchn,1)=i
187  isig(nchn,2)=-i
188  isig(nchn,3)=1
189  sigh(nchn)=0.5d0*facgg*fcoi*ei**4
190  170 CONTINUE
191  ENDIF
192 
193  ELSEIF(isub.LE.40) THEN
194  IF(isub.EQ.28) THEN
195 C...f + g -> f + g (q + g -> q + g only)
196  facqg1=comfac*as**2*4d0/9d0*((2d0+mstp(34)*1d0/4d0)*uh2/th2-
197  & uh/sh)*faca
198  facqg2=comfac*as**2*4d0/9d0*((2d0+mstp(34)*1d0/4d0)*sh2/th2-
199  & sh/uh)
200  DO 190 i=mmina,mmaxa
201  IF(i.EQ.0.OR.iabs(i).GT.10) goto 190
202  DO 180 isde=1,2
203  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 180
204  IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 180
205  nchn=nchn+1
206  isig(nchn,isde)=i
207  isig(nchn,3-isde)=21
208  isig(nchn,3)=1
209  sigh(nchn)=facqg1
210  nchn=nchn+1
211  isig(nchn,isde)=i
212  isig(nchn,3-isde)=21
213  isig(nchn,3)=2
214  sigh(nchn)=facqg2
215  180 CONTINUE
216  190 CONTINUE
217 
218  ELSEIF(isub.EQ.29) THEN
219 C...f + g -> f + gamma (q + g -> q + gamma only)
220  fgq=comfac*faca*as*aem*1d0/3d0*(sh2+uh2)/(-sh*uh)
221  DO 210 i=mmina,mmaxa
222  IF(i.EQ.0.OR.iabs(i).GT.mstp(58)) goto 210
223  ei=kchg(iabs(i),1)/3d0
224  facgq=fgq*ei**2
225  DO 200 isde=1,2
226  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 200
227  IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 200
228  nchn=nchn+1
229  isig(nchn,isde)=i
230  isig(nchn,3-isde)=21
231  isig(nchn,3)=1
232  sigh(nchn)=facgq
233  200 CONTINUE
234  210 CONTINUE
235 
236  ELSEIF(isub.EQ.33) THEN
237 C...f + gamma -> f + g (q + gamma -> q + g only)
238  fgq=comfac*as*aem*8d0/3d0*(sh2+uh2)/(-sh*uh)
239  DO 230 i=mmina,mmaxa
240  IF(i.EQ.0.OR.iabs(i).GT.mstp(58)) goto 230
241  ei=kchg(iabs(i),1)/3d0
242  facgq=fgq*ei**2
243  DO 220 isde=1,2
244  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,22).EQ.0) goto 220
245  IF(isde.EQ.2.AND.kfac(1,22)*kfac(2,i).EQ.0) goto 220
246  nchn=nchn+1
247  isig(nchn,isde)=i
248  isig(nchn,3-isde)=22
249  isig(nchn,3)=1
250  sigh(nchn)=facgq
251  220 CONTINUE
252  230 CONTINUE
253 
254  ELSEIF(isub.EQ.34) THEN
255 C...f + gamma -> f + gamma
256  fgq=comfac*aem**2*2d0*(sh2+uh2)/(-sh*uh)
257  DO 250 i=mmina,mmaxa
258  IF(i.EQ.0) goto 250
259  ei=kchg(iabs(i),1)/3d0
260  facgq=fgq*ei**4
261  DO 240 isde=1,2
262  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,22).EQ.0) goto 240
263  IF(isde.EQ.2.AND.kfac(1,22)*kfac(2,i).EQ.0) goto 240
264  nchn=nchn+1
265  isig(nchn,isde)=i
266  isig(nchn,3-isde)=22
267  isig(nchn,3)=1
268  sigh(nchn)=facgq
269  240 CONTINUE
270  250 CONTINUE
271  ENDIF
272 
273  ELSEIF(isub.LE.80) THEN
274  IF(isub.EQ.53) THEN
275 C...g + g -> f + fbar (g + g -> q + qbar only)
276  IF(kfac(1,21)*kfac(2,21).EQ.0) goto 270
277  idc0=mdcy(21,2)-1
278 C...Begin by d, u, s flavours.
279  flavwt=0d0
280  IF(mdme(idc0+1,1).GE.1) flavwt=flavwt+
281  & sqrt(max(0d0,1d0-4d0*pmas(1,1)**2/sh))
282  IF(mdme(idc0+2,1).GE.1) flavwt=flavwt+
283  & sqrt(max(0d0,1d0-4d0*pmas(2,1)**2/sh))
284  IF(mdme(idc0+3,1).GE.1) flavwt=flavwt+
285  & sqrt(max(0d0,1d0-4d0*pmas(3,1)**2/sh))
286  facqq1=comfac*as**2*1d0/6d0*(uh/th-(2d0+mstp(34)*1d0/4d0)*
287  & uh2/sh2)*flavwt*faca
288  facqq2=comfac*as**2*1d0/6d0*(th/uh-(2d0+mstp(34)*1d0/4d0)*
289  & th2/sh2)*flavwt*faca
290  nchn=nchn+1
291  isig(nchn,1)=21
292  isig(nchn,2)=21
293  isig(nchn,3)=1
294  sigh(nchn)=facqq1
295  nchn=nchn+1
296  isig(nchn,1)=21
297  isig(nchn,2)=21
298  isig(nchn,3)=2
299  sigh(nchn)=facqq2
300 C...Next c and b flavours: modified that and uhat for fixed
301 C...cos(theta-hat).
302  DO 260 ifl=4,5
303  sqmavg=pmas(ifl,1)**2
304  IF(mdme(idc0+ifl,1).GE.1.AND.sh.GT.4.04d0*sqmavg) THEN
305  be34=sqrt(1d0-4d0*sqmavg/sh)
306  thq=-0.5d0*sh*(1d0-be34*cth)
307  uhq=-0.5d0*sh*(1d0+be34*cth)
308  thuhq=thq*uhq-sqmavg*sh
309  IF(mstp(34).EQ.0) THEN
310  facqq1=uhq/thq-2d0*uhq**2/sh2+4d0*(sqmavg/sh)*thuhq/thq**2
311  facqq2=thq/uhq-2d0*thq**2/sh2+4d0*(sqmavg/sh)*thuhq/uhq**2
312  ELSE
313  facqq1=uhq/thq-2.25d0*uhq**2/sh2+4.5d0*(sqmavg/sh)*thuhq/
314  & thq**2+0.5d0*sqmavg*(thq+sqmavg)/thq**2-sqmavg**2/(sh*thq)
315  facqq2=thq/uhq-2.25d0*thq**2/sh2+4.5d0*(sqmavg/sh)*thuhq/
316  & uhq**2+0.5d0*sqmavg*(uhq+sqmavg)/uhq**2-sqmavg**2/(sh*uhq)
317  ENDIF
318  facqq1=comfac*faca*as**2*(1d0/6d0)*facqq1*be34
319  facqq2=comfac*faca*as**2*(1d0/6d0)*facqq2*be34
320  nchn=nchn+1
321  isig(nchn,1)=21
322  isig(nchn,2)=21
323  isig(nchn,3)=1+2*(ifl-3)
324  sigh(nchn)=facqq1
325  nchn=nchn+1
326  isig(nchn,1)=21
327  isig(nchn,2)=21
328  isig(nchn,3)=2+2*(ifl-3)
329  sigh(nchn)=facqq2
330  ENDIF
331  260 CONTINUE
332  270 CONTINUE
333 
334  ELSEIF(isub.EQ.54) THEN
335 C...g + gamma -> f + fbar (g + gamma -> q + qbar only)
336  CALL pywidt(21,sh,wdtp,wdte)
337  wdtesu=0d0
338  DO 280 i=1,min(8,mdcy(21,3))
339  ef=kchg(i,1)/3d0
340  wdtesu=wdtesu+ef**2*(wdte(i,1)+wdte(i,2)+wdte(i,3)+
341  & wdte(i,4))
342  280 CONTINUE
343  facqq=comfac*aem*as*wdtesu*(th2+uh2)/(th*uh)
344  IF(kfac(1,21)*kfac(2,22).NE.0) THEN
345  nchn=nchn+1
346  isig(nchn,1)=21
347  isig(nchn,2)=22
348  isig(nchn,3)=1
349  sigh(nchn)=facqq
350  ENDIF
351  IF(kfac(1,22)*kfac(2,21).NE.0) THEN
352  nchn=nchn+1
353  isig(nchn,1)=22
354  isig(nchn,2)=21
355  isig(nchn,3)=1
356  sigh(nchn)=facqq
357  ENDIF
358 
359  ELSEIF(isub.EQ.58) THEN
360 C...gamma + gamma -> f + fbar
361  CALL pywidt(22,sh,wdtp,wdte)
362  wdtesu=0d0
363  DO 290 i=1,min(12,mdcy(22,3))
364  IF(i.LE.8) ef= kchg(i,1)/3d0
365  IF(i.GE.9) ef= kchg(9+2*(i-8),1)/3d0
366  wdtesu=wdtesu+ef**2*(wdte(i,1)+wdte(i,2)+wdte(i,3)+
367  & wdte(i,4))
368  290 CONTINUE
369  facff=comfac*aem**2*wdtesu*2d0*(th2+uh2)/(th*uh)
370  IF(kfac(1,22)*kfac(2,22).NE.0) THEN
371  nchn=nchn+1
372  isig(nchn,1)=22
373  isig(nchn,2)=22
374  isig(nchn,3)=1
375  sigh(nchn)=facff
376  ENDIF
377 
378  ELSEIF(isub.EQ.68) THEN
379 C...g + g -> g + g
380  IF(kfac(1,21)*kfac(2,21).EQ.0) goto 300
381  facgg1=comfac*as**2*9d0/4d0*(sh2/th2+2d0*sh/th+3d0+2d0*th/sh+
382  & th2/sh2)*faca
383  facgg2=comfac*as**2*9d0/4d0*(uh2/sh2+2d0*uh/sh+3d0+2d0*sh/uh+
384  & sh2/uh2)*faca
385  facgg3=comfac*as**2*9d0/4d0*(th2/uh2+2d0*th/uh+3d0+2d0*uh/th+
386  & uh2/th2)
387  nchn=nchn+1
388  isig(nchn,1)=21
389  isig(nchn,2)=21
390  isig(nchn,3)=1
391  sigh(nchn)=0.5d0*facgg1
392  nchn=nchn+1
393  isig(nchn,1)=21
394  isig(nchn,2)=21
395  isig(nchn,3)=2
396  sigh(nchn)=0.5d0*facgg2
397  nchn=nchn+1
398  isig(nchn,1)=21
399  isig(nchn,2)=21
400  isig(nchn,3)=3
401  sigh(nchn)=0.5d0*facgg3
402  300 CONTINUE
403 
404  ELSEIF(isub.EQ.80) THEN
405 C...q + gamma -> q' + pi+/-
406  fqpi=comfac*(2d0*aem/9d0)*(-sh/th)*(1d0/sh2+1d0/th2)
407  assh=pyalps(max(0.5d0,0.5d0*sh))
408  q2fpsh=0.55d0/log(max(2d0,2d0*sh))
409  delsh=uh*sqrt(assh*q2fpsh)
410  asuh=pyalps(max(0.5d0,-0.5d0*uh))
411  q2fpuh=0.55d0/log(max(2d0,-2d0*uh))
412  deluh=sh*sqrt(asuh*q2fpuh)
413  DO 320 i=max(-2,mmina),min(2,mmaxa)
414  IF(i.EQ.0) goto 320
415  ei=kchg(iabs(i),1)/3d0
416  ej=sign(1d0-abs(ei),ei)
417  DO 310 isde=1,2
418  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,22).EQ.0) goto 310
419  IF(isde.EQ.2.AND.kfac(1,22)*kfac(2,i).EQ.0) goto 310
420  nchn=nchn+1
421  isig(nchn,isde)=i
422  isig(nchn,3-isde)=22
423  isig(nchn,3)=1
424  sigh(nchn)=fqpi*(ei*delsh+ej*deluh)**2
425  310 CONTINUE
426  320 CONTINUE
427  ENDIF
428 
429  ELSEIF(isub.LE.100) THEN
430  IF(isub.EQ.91) THEN
431 C...Elastic scattering
432  sigs=vint(315)*vint(316)*sigt(0,0,1)
433 
434  ELSEIF(isub.EQ.92) THEN
435 C...Single diffractive scattering (first side, i.e. XB)
436  sigs=vint(315)*vint(316)*sigt(0,0,2)
437 
438  ELSEIF(isub.EQ.93) THEN
439 C...Single diffractive scattering (second side, i.e. AX)
440  sigs=vint(315)*vint(316)*sigt(0,0,3)
441 
442  ELSEIF(isub.EQ.94) THEN
443 C...Double diffractive scattering
444  sigs=vint(315)*vint(316)*sigt(0,0,4)
445 
446  ELSEIF(isub.EQ.95) THEN
447 C...Low-pT scattering
448  sigs=vint(315)*vint(316)*sigt(0,0,5)
449 
450  ELSEIF(isub.EQ.96) THEN
451 C...Multiple interactions: sum of QCD processes
452  CALL pywidt(21,sh,wdtp,wdte)
453 
454 C...q + q' -> q + q'
455  facqq1=comfac*as**2*4d0/9d0*(sh2+uh2)/th2
456  facqqb=comfac*as**2*4d0/9d0*((sh2+uh2)/th2*faca-
457  & mstp(34)*2d0/3d0*uh2/(sh*th))
458  facqq2=comfac*as**2*4d0/9d0*(sh2+th2)/uh2
459  facqqi=-comfac*as**2*4d0/9d0*mstp(34)*2d0/3d0*sh2/(th*uh)
460  ratqqi=(facqq1+facqq2+facqqi)/(facqq1+facqq2)
461  DO 340 i=-5,5
462  IF(i.EQ.0) goto 340
463  DO 330 j=-5,5
464  IF(j.EQ.0) goto 330
465  nchn=nchn+1
466  isig(nchn,1)=i
467  isig(nchn,2)=j
468  isig(nchn,3)=111
469  sigh(nchn)=facqq1
470  IF(i.EQ.-j) sigh(nchn)=facqqb
471  IF(i.EQ.j) THEN
472  sigh(nchn)=0.5d0*facqq1*ratqqi
473  nchn=nchn+1
474  isig(nchn,1)=i
475  isig(nchn,2)=j
476  isig(nchn,3)=112
477  sigh(nchn)=0.5d0*facqq2*ratqqi
478  ENDIF
479  330 CONTINUE
480  340 CONTINUE
481 
482 C...q + qbar -> q' + qbar' or g + g
483  facqqb=comfac*as**2*4d0/9d0*(th2+uh2)/sh2*
484  & (wdte(0,1)+wdte(0,2)+wdte(0,3)+wdte(0,4))
485  facgg1=comfac*as**2*32d0/27d0*(uh/th-(2d0+mstp(34)*1d0/4d0)*
486  & uh2/sh2)
487  facgg2=comfac*as**2*32d0/27d0*(th/uh-(2d0+mstp(34)*1d0/4d0)*
488  & th2/sh2)
489  DO 350 i=-5,5
490  IF(i.EQ.0) goto 350
491  nchn=nchn+1
492  isig(nchn,1)=i
493  isig(nchn,2)=-i
494  isig(nchn,3)=121
495  sigh(nchn)=facqqb
496  nchn=nchn+1
497  isig(nchn,1)=i
498  isig(nchn,2)=-i
499  isig(nchn,3)=131
500  sigh(nchn)=0.5d0*facgg1
501  nchn=nchn+1
502  isig(nchn,1)=i
503  isig(nchn,2)=-i
504  isig(nchn,3)=132
505  sigh(nchn)=0.5d0*facgg2
506  350 CONTINUE
507 
508 C...q + g -> q + g
509  facqg1=comfac*as**2*4d0/9d0*((2d0+mstp(34)*1d0/4d0)*uh2/th2-
510  & uh/sh)*faca
511  facqg2=comfac*as**2*4d0/9d0*((2d0+mstp(34)*1d0/4d0)*sh2/th2-
512  & sh/uh)
513  DO 370 i=-5,5
514  IF(i.EQ.0) goto 370
515  DO 360 isde=1,2
516  nchn=nchn+1
517  isig(nchn,isde)=i
518  isig(nchn,3-isde)=21
519  isig(nchn,3)=281
520  sigh(nchn)=facqg1
521  nchn=nchn+1
522  isig(nchn,isde)=i
523  isig(nchn,3-isde)=21
524  isig(nchn,3)=282
525  sigh(nchn)=facqg2
526  360 CONTINUE
527  370 CONTINUE
528 
529 C...g + g -> q + qbar (only d, u, s)
530  idc0=mdcy(21,2)-1
531  flavwt=0d0
532  IF(mdme(idc0+1,1).GE.1) flavwt=flavwt+
533  & sqrt(max(0d0,1d0-4d0*pmas(1,1)**2/sh))
534  IF(mdme(idc0+2,1).GE.1) flavwt=flavwt+
535  & sqrt(max(0d0,1d0-4d0*pmas(2,1)**2/sh))
536  IF(mdme(idc0+3,1).GE.1) flavwt=flavwt+
537  & sqrt(max(0d0,1d0-4d0*pmas(3,1)**2/sh))
538  facqq1=comfac*as**2*1d0/6d0*(uh/th-(2d0+mstp(34)*1d0/4d0)*
539  & uh2/sh2)*flavwt*faca
540  facqq2=comfac*as**2*1d0/6d0*(th/uh-(2d0+mstp(34)*1d0/4d0)*
541  & th2/sh2)*flavwt*faca
542  nchn=nchn+1
543  isig(nchn,1)=21
544  isig(nchn,2)=21
545  isig(nchn,3)=531
546  sigh(nchn)=facqq1
547  nchn=nchn+1
548  isig(nchn,1)=21
549  isig(nchn,2)=21
550  isig(nchn,3)=532
551  sigh(nchn)=facqq2
552 
553 C...g + g -> c + cbar, b + bbar: modified that/uhat for fixed
554 C...cos(theta-hat)
555  DO 380 ifl=4,5
556  sqmavg=pmas(ifl,1)**2
557  IF(mdme(idc0+ifl,1).GE.1.AND.sh.GT.4.04d0*sqmavg) THEN
558  be34=sqrt(1d0-4d0*sqmavg/sh)
559  thq=-0.5d0*sh*(1d0-be34*cth)
560  uhq=-0.5d0*sh*(1d0+be34*cth)
561  thuhq=thq*uhq-sqmavg*sh
562  IF(mstp(34).EQ.0) THEN
563  facqq1=uhq/thq-2d0*uhq**2/sh2+4d0*(sqmavg/sh)*thuhq/thq**2
564  facqq2=thq/uhq-2d0*thq**2/sh2+4d0*(sqmavg/sh)*thuhq/uhq**2
565  ELSE
566  facqq1=uhq/thq-2.25d0*uhq**2/sh2+4.5d0*(sqmavg/sh)*thuhq/
567  & thq**2+0.5d0*sqmavg*(thq+sqmavg)/thq**2-sqmavg**2/(sh*thq)
568  facqq2=thq/uhq-2.25d0*thq**2/sh2+4.5d0*(sqmavg/sh)*thuhq/
569  & uhq**2+0.5d0*sqmavg*(uhq+sqmavg)/uhq**2-sqmavg**2/(sh*uhq)
570  ENDIF
571  facqq1=comfac*faca*as**2*(1d0/6d0)*facqq1*be34
572  facqq2=comfac*faca*as**2*(1d0/6d0)*facqq2*be34
573  nchn=nchn+1
574  isig(nchn,1)=21
575  isig(nchn,2)=21
576  isig(nchn,3)=531+2*(ifl-3)
577  sigh(nchn)=facqq1
578  nchn=nchn+1
579  isig(nchn,1)=21
580  isig(nchn,2)=21
581  isig(nchn,3)=532+2*(ifl-3)
582  sigh(nchn)=facqq2
583  ENDIF
584  380 CONTINUE
585 
586 C...g + g -> g + g
587  facgg1=comfac*as**2*9d0/4d0*(sh2/th2+2d0*sh/th+3d0+
588  & 2d0*th/sh+th2/sh2)*faca
589  facgg2=comfac*as**2*9d0/4d0*(uh2/sh2+2d0*uh/sh+3d0+
590  & 2d0*sh/uh+sh2/uh2)*faca
591  facgg3=comfac*as**2*9d0/4d0*(th2/uh2+2d0*th/uh+3+
592  & 2d0*uh/th+uh2/th2)
593  nchn=nchn+1
594  isig(nchn,1)=21
595  isig(nchn,2)=21
596  isig(nchn,3)=681
597  sigh(nchn)=0.5d0*facgg1
598  nchn=nchn+1
599  isig(nchn,1)=21
600  isig(nchn,2)=21
601  isig(nchn,3)=682
602  sigh(nchn)=0.5d0*facgg2
603  nchn=nchn+1
604  isig(nchn,1)=21
605  isig(nchn,2)=21
606  isig(nchn,3)=683
607  sigh(nchn)=0.5d0*facgg3
608 
609  ELSEIF(isub.EQ.99) THEN
610 C...f + gamma* -> f.
611  IF(mint(107).EQ.4) THEN
612  q2ga=vint(307)
613  p2ga=vint(308)
614  isde=2
615  ELSE
616  q2ga=vint(308)
617  p2ga=vint(307)
618  isde=1
619  ENDIF
620  comfac=paru(5)*4d0*paru(1)**2*paru(101)*vint(315)*vint(316)
621  pm2rho=pmas(pycomp(113),1)**2
622  IF(mstp(19).EQ.0) THEN
623  comfac=comfac/q2ga
624 C... To use MSTP(19).EQ.1 (less Q2 suppression) with the right factor (1-x)^-1
625  ELSEIF(mstp(19).EQ.1) THEN
626  comfac=comfac/(q2ga+pm2rho)
627  w2ga=vint(2)
628  IF(mint(11).EQ.22.AND.mint(12).EQ.22) THEN
629  xga=q2ga/(w2ga+vint(307)+vint(308))
630  ELSE
631  xga=q2ga/(w2ga+q2ga-pmas(pycomp(mint(10+isde)),1)**2)
632  ENDIF
633  comfac=comfac/max(1d-2,1d0-xga)
634 
635  ELSEIF(mstp(19).EQ.2) THEN
636  comfac=comfac*q2ga/(q2ga+pm2rho)**2
637  ELSE
638  comfac=comfac*q2ga/(q2ga+pm2rho)**2
639  w2ga=vint(2)
640  IF(mint(11).EQ.22.AND.mint(12).EQ.22) THEN
641  rdrds=4.1d-3*w2ga**2.167d0/((q2ga+0.15d0*w2ga)**2*
642  & q2ga**0.75d0)*(1d0+0.11d0*q2ga*p2ga/(1d0+0.02d0*p2ga**2))
643  xga=q2ga/(w2ga+vint(307)+vint(308))
644  ELSE
645  rdrds=1.5d-4*w2ga**2.167d0/((q2ga+0.041d0*w2ga)**2*
646  & q2ga**0.57d0)
647  xga=q2ga/(w2ga+q2ga-pmas(pycomp(mint(10+isde)),1)**2)
648  ENDIF
649  comfac=comfac*exp(-max(1d-10,rdrds))
650  IF(mstp(19).EQ.4) comfac=comfac/max(1d-2,1d0-xga)
651  ENDIF
652  DO 390 i=mmina,mmaxa
653  IF(i.EQ.0.OR.kfac(isde,i).EQ.0) goto 390
654  IF(iabs(i).LT.10.AND.iabs(i).GT.mstp(58)) goto 390
655  ei=kchg(iabs(i),1)/3d0
656  nchn=nchn+1
657  isig(nchn,isde)=i
658  isig(nchn,3-isde)=22
659  isig(nchn,3)=1
660  sigh(nchn)=comfac*ei**2
661  390 CONTINUE
662  ENDIF
663 
664  ELSE
665  IF(isub.EQ.114.OR.isub.EQ.115) THEN
666 C...g + g -> gamma + gamma or g + g -> g + gamma
667  a0stur=0d0
668  a0stui=0d0
669  a0tsur=0d0
670  a0tsui=0d0
671  a0utsr=0d0
672  a0utsi=0d0
673  a1stur=0d0
674  a1stui=0d0
675  a2stur=0d0
676  a2stui=0d0
677  alst=log(-sh/th)
678  alsu=log(-sh/uh)
679  altu=log(th/uh)
680  imax=2*mstp(1)
681  IF(mstp(38).GE.1.AND.mstp(38).LE.8) imax=mstp(38)
682  DO 400 i=1,imax
683  ei=kchg(iabs(i),1)/3d0
684  eiwt=ei**2
685  IF(isub.EQ.115) eiwt=ei
686  sqmq=pmas(i,1)**2
687  epss=4d0*sqmq/sh
688  epst=4d0*sqmq/th
689  epsu=4d0*sqmq/uh
690  IF((mstp(38).GE.1.AND.mstp(38).LE.8).OR.epss.LT.1d-4) THEN
691  b0stur=1d0+(th-uh)/sh*altu+0.5d0*(th2+uh2)/sh2*(altu**2+
692  & paru(1)**2)
693  b0stui=0d0
694  b0tsur=1d0+(sh-uh)/th*alsu+0.5d0*(sh2+uh2)/th2*alsu**2
695  b0tsui=-paru(1)*((sh-uh)/th+(sh2+uh2)/th2*alsu)
696  b0utsr=1d0+(sh-th)/uh*alst+0.5d0*(sh2+th2)/uh2*alst**2
697  b0utsi=-paru(1)*((sh-th)/uh+(sh2+th2)/uh2*alst)
698  b1stur=-1d0
699  b1stui=0d0
700  b2stur=-1d0
701  b2stui=0d0
702  ELSE
703  CALL pywaux(1,epss,w1sr,w1si)
704  CALL pywaux(1,epst,w1tr,w1ti)
705  CALL pywaux(1,epsu,w1ur,w1ui)
706  CALL pywaux(2,epss,w2sr,w2si)
707  CALL pywaux(2,epst,w2tr,w2ti)
708  CALL pywaux(2,epsu,w2ur,w2ui)
709  CALL pyi3au(epss,th/uh,y3stur,y3stui)
710  CALL pyi3au(epss,uh/th,y3sutr,y3suti)
711  CALL pyi3au(epst,sh/uh,y3tsur,y3tsui)
712  CALL pyi3au(epst,uh/sh,y3tusr,y3tusi)
713  CALL pyi3au(epsu,sh/th,y3ustr,y3usti)
714  CALL pyi3au(epsu,th/sh,y3utsr,y3utsi)
715  b0stur=1d0+(1d0+2d0*th/sh)*w1tr+(1d0+2d0*uh/sh)*w1ur+
716  & 0.5d0*((th2+uh2)/sh2-epss)*(w2tr+w2ur)-
717  & 0.25d0*epst*(1d0-0.5d0*epss)*(y3sutr+y3tusr)-
718  & 0.25d0*epsu*(1d0-0.5d0*epss)*(y3stur+y3utsr)+
719  & 0.25d0*(-2d0*(th2+uh2)/sh2+4d0*epss+epst+epsu+
720  & 0.5d0*epst*epsu)*(y3tsur+y3ustr)
721  b0stui=(1d0+2d0*th/sh)*w1ti+(1d0+2d0*uh/sh)*w1ui+
722  & 0.5d0*((th2+uh2)/sh2-epss)*(w2ti+w2ui)-
723  & 0.25d0*epst*(1d0-0.5d0*epss)*(y3suti+y3tusi)-
724  & 0.25d0*epsu*(1d0-0.5d0*epss)*(y3stui+y3utsi)+
725  & 0.25d0*(-2d0*(th2+uh2)/sh2+4d0*epss+epst+epsu+
726  & 0.5d0*epst*epsu)*(y3tsui+y3usti)
727  b0tsur=1d0+(1d0+2d0*sh/th)*w1sr+(1d0+2d0*uh/th)*w1ur+
728  & 0.5d0*((sh2+uh2)/th2-epst)*(w2sr+w2ur)-
729  & 0.25d0*epss*(1d0-0.5d0*epst)*(y3tusr+y3sutr)-
730  & 0.25d0*epsu*(1d0-0.5d0*epst)*(y3tsur+y3ustr)+
731  & 0.25d0*(-2d0*(sh2+uh2)/th2+4d0*epst+epss+epsu+
732  & 0.5d0*epss*epsu)*(y3stur+y3utsr)
733  b0tsui=(1d0+2d0*sh/th)*w1si+(1d0+2d0*uh/th)*w1ui+
734  & 0.5d0*((sh2+uh2)/th2-epst)*(w2si+w2ui)-
735  & 0.25d0*epss*(1d0-0.5d0*epst)*(y3tusi+y3suti)-
736  & 0.25d0*epsu*(1d0-0.5d0*epst)*(y3tsui+y3usti)+
737  & 0.25d0*(-2d0*(sh2+uh2)/th2+4d0*epst+epss+epsu+
738  & 0.5d0*epss*epsu)*(y3stui+y3utsi)
739  b0utsr=1d0+(1d0+2d0*th/uh)*w1tr+(1d0+2d0*sh/uh)*w1sr+
740  & 0.5d0*((th2+sh2)/uh2-epsu)*(w2tr+w2sr)-
741  & 0.25d0*epst*(1d0-0.5d0*epsu)*(y3ustr+y3tsur)-
742  & 0.25d0*epss*(1d0-0.5d0*epsu)*(y3utsr+y3stur)+
743  & 0.25d0*(-2d0*(th2+sh2)/uh2+4d0*epsu+epst+epss+
744  & 0.5d0*epst*epss)*(y3tusr+y3sutr)
745  b0utsi=(1d0+2d0*th/uh)*w1ti+(1d0+2d0*sh/uh)*w1si+
746  & 0.5d0*((th2+sh2)/uh2-epsu)*(w2ti+w2si)-
747  & 0.25d0*epst*(1d0-0.5d0*epsu)*(y3usti+y3tsui)-
748  & 0.25d0*epss*(1d0-0.5d0*epsu)*(y3utsi+y3stui)+
749  & 0.25d0*(-2d0*(th2+sh2)/uh2+4d0*epsu+epst+epss+
750  & 0.5d0*epst*epss)*(y3tusi+y3suti)
751  b1stur=-1d0-0.25d0*(epss+epst+epsu)*(w2sr+w2tr+w2ur)+
752  & 0.25d0*(epsu+0.5d0*epss*epst)*(y3sutr+y3tusr)+
753  & 0.25d0*(epst+0.5d0*epss*epsu)*(y3stur+y3utsr)+
754  & 0.25d0*(epss+0.5d0*epst*epsu)*(y3tsur+y3ustr)
755  b1stui=-0.25d0*(epss+epst+epsu)*(w2si+w2ti+w2ui)+
756  & 0.25d0*(epsu+0.5d0*epss*epst)*(y3suti+y3tusi)+
757  & 0.25d0*(epst+0.5d0*epss*epsu)*(y3stui+y3utsi)+
758  & 0.25d0*(epss+0.5d0*epst*epsu)*(y3tsui+y3usti)
759  b2stur=-1d0+0.125d0*epss*epst*(y3sutr+y3tusr)+
760  & 0.125d0*epss*epsu*(y3stur+y3utsr)+
761  & 0.125d0*epst*epsu*(y3tsur+y3ustr)
762  b2stui=0.125d0*epss*epst*(y3suti+y3tusi)+
763  & 0.125d0*epss*epsu*(y3stui+y3utsi)+
764  & 0.125d0*epst*epsu*(y3tsui+y3usti)
765  ENDIF
766  a0stur=a0stur+eiwt*b0stur
767  a0stui=a0stui+eiwt*b0stui
768  a0tsur=a0tsur+eiwt*b0tsur
769  a0tsui=a0tsui+eiwt*b0tsui
770  a0utsr=a0utsr+eiwt*b0utsr
771  a0utsi=a0utsi+eiwt*b0utsi
772  a1stur=a1stur+eiwt*b1stur
773  a1stui=a1stui+eiwt*b1stui
774  a2stur=a2stur+eiwt*b2stur
775  a2stui=a2stui+eiwt*b2stui
776  400 CONTINUE
777  asqsum=a0stur**2+a0stui**2+a0tsur**2+a0tsui**2+a0utsr**2+
778  & a0utsi**2+4d0*a1stur**2+4d0*a1stui**2+a2stur**2+a2stui**2
779  facgg=comfac*faca/(16d0*paru(1)**2)*as**2*aem**2*asqsum
780  facgp=comfac*faca*5d0/(192d0*paru(1)**2)*as**3*aem*asqsum
781  IF(kfac(1,21)*kfac(2,21).EQ.0) goto 410
782  nchn=nchn+1
783  isig(nchn,1)=21
784  isig(nchn,2)=21
785  isig(nchn,3)=1
786  IF(isub.EQ.114) sigh(nchn)=0.5d0*facgg
787  IF(isub.EQ.115) sigh(nchn)=facgp
788  410 CONTINUE
789 
790  ELSEIF(isub.EQ.131.OR.isub.EQ.132) THEN
791 C...f + gamma*_(T,L) -> f + g (q + gamma*_(T,L) -> q + g only)
792  ph=0d0
793  IF(mint(15).EQ.22.AND.mint(107).EQ.0.AND.vint(3).LT.0d0)
794  & ph=vint(3)**2
795  IF(mint(16).EQ.22.AND.mint(108).EQ.0.AND.vint(4).LT.0d0)
796  & ph=vint(4)**2
797  IF(isub.EQ.131) THEN
798  fgq=comfac*as*aem*8d0/3d0*sh**2/(sh+ph)**2*
799  & ((sh2+uh2-2d0*ph*th)/(-sh*uh)-2d0*ph*th/(sh+ph)**2)
800  ELSE
801  fgq=comfac*as*aem*8d0/3d0*sh**2/(sh+ph)**4*(-4d0*ph*th)
802  ENDIF
803  DO 430 i=mmina,mmaxa
804  IF(i.EQ.0.OR.iabs(i).GT.mstp(58)) goto 430
805  ei=kchg(iabs(i),1)/3d0
806  facgq=fgq*ei**2
807  DO 420 isde=1,2
808  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,22).EQ.0) goto 420
809  IF(isde.EQ.2.AND.kfac(1,22)*kfac(2,i).EQ.0) goto 420
810  nchn=nchn+1
811  isig(nchn,isde)=i
812  isig(nchn,3-isde)=22
813  isig(nchn,3)=1
814  sigh(nchn)=facgq
815  420 CONTINUE
816  430 CONTINUE
817 
818  ELSEIF(isub.EQ.133.OR.isub.EQ.134) THEN
819 C...f + gamma*_(T,L) -> f + gamma
820  ph=0d0
821  IF(mint(15).EQ.22.AND.mint(107).EQ.0.AND.vint(3).LT.0d0)
822  & ph=vint(3)**2
823  IF(mint(16).EQ.22.AND.mint(108).EQ.0.AND.vint(4).LT.0d0)
824  & ph=vint(4)**2
825  IF(isub.EQ.133) THEN
826  fgq=comfac*aem**2*2d0*sh**2/(sh+ph)**2*
827  & ((sh2+uh2-2d0*ph*th)/(-sh*uh)-2d0*ph*th/(sh+ph)**2)
828  ELSE
829  fgq=comfac*aem**2*2d0*sh**2/(sh+ph)**4*(-4d0*ph*th)
830  ENDIF
831  DO 450 i=mmina,mmaxa
832  IF(i.EQ.0) goto 450
833  ei=kchg(iabs(i),1)/3d0
834  facgq=fgq*ei**4
835  DO 440 isde=1,2
836  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,22).EQ.0) goto 440
837  IF(isde.EQ.2.AND.kfac(1,22)*kfac(2,i).EQ.0) goto 440
838  nchn=nchn+1
839  isig(nchn,isde)=i
840  isig(nchn,3-isde)=22
841  isig(nchn,3)=1
842  sigh(nchn)=facgq
843  440 CONTINUE
844  450 CONTINUE
845 
846  ELSEIF(isub.EQ.135.OR.isub.EQ.136) THEN
847 C...g + gamma*_(T,L) -> f + fbar (g + gamma*_(T,L) -> q + qbar only)
848  ph=0d0
849  IF(mint(15).EQ.22.AND.mint(107).EQ.0.AND.vint(3).LT.0d0)
850  & ph=vint(3)**2
851  IF(mint(16).EQ.22.AND.mint(108).EQ.0.AND.vint(4).LT.0d0)
852  & ph=vint(4)**2
853  CALL pywidt(21,sh,wdtp,wdte)
854  wdtesu=0d0
855  DO 460 i=1,min(8,mdcy(21,3))
856  ef=kchg(i,1)/3d0
857  wdtesu=wdtesu+ef**2*(wdte(i,1)+wdte(i,2)+wdte(i,3)+
858  & wdte(i,4))
859  460 CONTINUE
860  IF(isub.EQ.135) THEN
861  facqq=comfac*aem*as*wdtesu*sh**2/(sh+ph)**2*
862  & ((th2+uh2-2d0*ph*sh)/(th*uh)+4d0*ph*sh/(sh+ph)**2)
863  ELSE
864  facqq=comfac*aem*as*wdtesu*sh**2/(sh+ph)**4*8d0*ph*sh
865  ENDIF
866  IF(kfac(1,21)*kfac(2,22).NE.0) THEN
867  nchn=nchn+1
868  isig(nchn,1)=21
869  isig(nchn,2)=22
870  isig(nchn,3)=1
871  sigh(nchn)=facqq
872  ENDIF
873  IF(kfac(1,22)*kfac(2,21).NE.0) THEN
874  nchn=nchn+1
875  isig(nchn,1)=22
876  isig(nchn,2)=21
877  isig(nchn,3)=1
878  sigh(nchn)=facqq
879  ENDIF
880 
881  ELSEIF(isub.GE.137.AND.isub.LE.140) THEN
882 C...gamma*_(T,L) + gamma*_(T,L) -> f + fbar
883  ph1=0d0
884  IF(vint(3).LT.0d0) ph1=vint(3)**2
885  ph2=0d0
886  IF(vint(4).LT.0d0) ph2=vint(4)**2
887  CALL pywidt(22,sh,wdtp,wdte)
888  wdtesu=0d0
889  DO 470 i=1,min(12,mdcy(22,3))
890  IF(i.LE.8) ef= kchg(i,1)/3d0
891  IF(i.GE.9) ef= kchg(9+2*(i-8),1)/3d0
892  wdtesu=wdtesu+ef**2*(wdte(i,1)+wdte(i,2)+wdte(i,3)+
893  & wdte(i,4))
894  470 CONTINUE
895  dlamb2=(th+uh)**2-4d0*ph1*ph2
896  IF(isub.EQ.137) THEN
897  fparam=-sh*(th+uh)/dlamb2
898  facff=comfac*aem**2*wdtesu*2d0*sh2/(dlamb2*th2*uh2)*
899  & (th*uh-ph1*ph2)*((th2+uh2)*(1d0-2d0*fparam*(1d0-fparam))-
900  & 2d0*ph1*ph2*fparam**2)
901  ELSEIF(isub.EQ.138) THEN
902  facff=comfac*aem**2*wdtesu*4d0*sh2*sh/(dlamb2**2*th2*uh2)*
903  & ph2*(4d0*(th*uh-ph1*ph2)*(th*uh+ph1*sh*(th-uh)**2/dlamb2)+
904  & 2d0*ph1**2*(th-uh)**2)
905  ELSEIF(isub.EQ.139) THEN
906  facff=comfac*aem**2*wdtesu*4d0*sh2*sh/(dlamb2**2*th2*uh2)*
907  & ph1*(4d0*(th*uh-ph1*ph2)*(th*uh+ph2*sh*(th-uh)**2/dlamb2)+
908  & 2d0*ph2**2*(th-uh)**2)
909  ELSE
910  facff=comfac*aem**2*wdtesu*32d0*sh2**2/(dlamb2**3*th2*uh2)*
911  & ph1*ph2*(th*uh-ph1*ph2)*(th-uh)**2
912  ENDIF
913  IF(kfac(1,22)*kfac(2,22).NE.0) THEN
914  nchn=nchn+1
915  isig(nchn,1)=22
916  isig(nchn,2)=22
917  isig(nchn,3)=1
918  sigh(nchn)=facff
919  ENDIF
920 
921  ENDIF
922  ENDIF
923 
924  RETURN
925  END