Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pysghf.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pysghf.f
1 
2 C*********************************************************************
3 
4 C...PYSGHF
5 C...Subprocess cross sections for heavy flavour production,
6 C...open and closed.
7 C...Auxiliary to PYSIGH.
8 
9  SUBROUTINE pysghf(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/pypars/mstp(200),parp(200),msti(200),pari(200)
22  common/pyint1/mint(400),vint(400)
23  common/pyint2/iset(500),kfpr(500,2),coef(500,20),icol(40,4,2)
24  common/pyint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)
25  common/pyint4/mwid(500),wids(500,5)
26  common/pysgcm/isub,isubsv,mmin1,mmax1,mmin2,mmax2,mmina,mmaxa,
27  &kfac(2,-40:40),comfac,fack,faca,sh,th,uh,sh2,th2,uh2,sqm3,sqm4,
28  &shr,sqpth,taup,be34,cth,x(2),sqmz,sqmw,gmmz,gmmw,
29  &aem,as,xw,xw1,xwc,xwv,poll,polr,polll,polrr
30  SAVE /pydat1/,/pydat2/,/pypars/,/pyint1/,/pyint2/,/pyint3/,
31  &/pyint4/,/pysgcm/
32 C...Local arrays
33  dimension wdtp(0:400),wdte(0:400,0:5)
34 
35 C...Determine where are charmonium/bottomonium wave function parameters.
36  ionium=140
37  IF(isub.GE.461.AND.isub.LE.479) ionium=145
38 
39 C...Convert bottomonium process into equivalent charmonium ones.
40  IF(isub.GE.461.AND.isub.LE.479) isub=isub-40
41 
42 C...Differential cross section expressions.
43 
44  IF(isub.LE.100) THEN
45  IF(isub.EQ.81) THEN
46 C...q + qbar -> Q + Qbar
47  sqmavg=0.5d0*(sqm3+sqm4)-0.25d0*(sqm3-sqm4)**2/sh
48  thq=-0.5d0*sh*(1d0-be34*cth)
49  uhq=-0.5d0*sh*(1d0+be34*cth)
50  facqqb=comfac*as**2*4d0/9d0*((thq**2+uhq**2)/sh2+
51  & 2d0*sqmavg/sh)
52  IF(mstp(35).GE.1) facqqb=facqqb*pyhfth(sh,sqmavg,0d0)
53  wid2=1d0
54  IF(mint(55).EQ.6) wid2=wids(6,1)
55  IF(mint(55).EQ.7.OR.mint(55).EQ.8) wid2=wids(mint(55),1)
56  facqqb=facqqb*wid2
57  DO 100 i=mmina,mmaxa
58  IF(i.EQ.0.OR.iabs(i).GT.mstp(58).OR.
59  & kfac(1,i)*kfac(2,-i).EQ.0) goto 100
60  nchn=nchn+1
61  isig(nchn,1)=i
62  isig(nchn,2)=-i
63  isig(nchn,3)=1
64  sigh(nchn)=facqqb
65  100 CONTINUE
66 
67  ELSEIF(isub.EQ.82) THEN
68 C...g + g -> Q + Qbar
69  sqmavg=0.5d0*(sqm3+sqm4)-0.25d0*(sqm3-sqm4)**2/sh
70  thq=-0.5d0*sh*(1d0-be34*cth)
71  uhq=-0.5d0*sh*(1d0+be34*cth)
72  thuhq=thq*uhq-sqmavg*sh
73  IF(mstp(34).EQ.0) THEN
74  facqq1=uhq/thq-2d0*uhq**2/sh2+4d0*(sqmavg/sh)*thuhq/thq**2
75  facqq2=thq/uhq-2d0*thq**2/sh2+4d0*(sqmavg/sh)*thuhq/uhq**2
76  ELSE
77  facqq1=uhq/thq-2.25d0*uhq**2/sh2+4.5d0*(sqmavg/sh)*thuhq/
78  & thq**2+0.5d0*sqmavg*(thq+sqmavg)/thq**2-sqmavg**2/(sh*thq)
79  facqq2=thq/uhq-2.25d0*thq**2/sh2+4.5d0*(sqmavg/sh)*thuhq/
80  & uhq**2+0.5d0*sqmavg*(uhq+sqmavg)/uhq**2-sqmavg**2/(sh*uhq)
81  ENDIF
82  facqq1=comfac*faca*as**2*(1d0/6d0)*facqq1
83  facqq2=comfac*faca*as**2*(1d0/6d0)*facqq2
84  IF(mstp(35).GE.1) THEN
85  fatre=pyhfth(sh,sqmavg,2d0/7d0)
86  facqq1=facqq1*fatre
87  facqq2=facqq2*fatre
88  ENDIF
89  wid2=1d0
90  IF(mint(55).EQ.6) wid2=wids(6,1)
91  IF(mint(55).EQ.7.OR.mint(55).EQ.8) wid2=wids(mint(55),1)
92  facqq1=facqq1*wid2
93  facqq2=facqq2*wid2
94  IF(kfac(1,21)*kfac(2,21).EQ.0) goto 110
95  nchn=nchn+1
96  isig(nchn,1)=21
97  isig(nchn,2)=21
98  isig(nchn,3)=1
99  sigh(nchn)=facqq1
100  nchn=nchn+1
101  isig(nchn,1)=21
102  isig(nchn,2)=21
103  isig(nchn,3)=2
104  sigh(nchn)=facqq2
105  110 CONTINUE
106 
107  ELSEIF(isub.EQ.83) THEN
108 C...f + q -> f' + Q
109  facqqs=comfac*(0.5d0*aem/xw)**2*sh*(sh-sqm3)/(sqmw-th)**2
110  facqqu=comfac*(0.5d0*aem/xw)**2*uh*(uh-sqm3)/(sqmw-th)**2
111  DO 130 i=mmin1,mmax1
112  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 130
113  DO 120 j=mmin2,mmax2
114  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 120
115  IF(i*j.GT.0.AND.mod(iabs(i+j),2).EQ.0) goto 120
116  IF(i*j.LT.0.AND.mod(iabs(i+j),2).EQ.1) goto 120
117  IF(iabs(i).LT.mint(55).AND.mod(iabs(i+mint(55)),2).EQ.1)
118  & THEN
119  nchn=nchn+1
120  isig(nchn,1)=i
121  isig(nchn,2)=j
122  isig(nchn,3)=1
123  IF(mod(mint(55),2).EQ.0) facckm=vckm(mint(55)/2,
124  & (iabs(i)+1)/2)*vint(180+j)
125  IF(mod(mint(55),2).EQ.1) facckm=vckm(iabs(i)/2,
126  & (mint(55)+1)/2)*vint(180+j)
127  wid2=1d0
128  IF(i.GT.0) THEN
129  IF(mint(55).EQ.6) wid2=wids(6,2)
130  IF(mint(55).EQ.7.OR.mint(55).EQ.8) wid2=
131  & wids(mint(55),2)
132  ELSE
133  IF(mint(55).EQ.6) wid2=wids(6,3)
134  IF(mint(55).EQ.7.OR.mint(55).EQ.8) wid2=
135  & wids(mint(55),3)
136  ENDIF
137  IF(i*j.GT.0) sigh(nchn)=facqqs*facckm*wid2
138  IF(i*j.LT.0) sigh(nchn)=facqqu*facckm*wid2
139  ENDIF
140  IF(iabs(j).LT.mint(55).AND.mod(iabs(j+mint(55)),2).EQ.1)
141  & THEN
142  nchn=nchn+1
143  isig(nchn,1)=i
144  isig(nchn,2)=j
145  isig(nchn,3)=2
146  IF(mod(mint(55),2).EQ.0) facckm=vckm(mint(55)/2,
147  & (iabs(j)+1)/2)*vint(180+i)
148  IF(mod(mint(55),2).EQ.1) facckm=vckm(iabs(j)/2,
149  & (mint(55)+1)/2)*vint(180+i)
150  IF(j.GT.0) THEN
151  IF(mint(55).EQ.6) wid2=wids(6,2)
152  IF(mint(55).EQ.7.OR.mint(55).EQ.8) wid2=
153  & wids(mint(55),2)
154  ELSE
155  IF(mint(55).EQ.6) wid2=wids(6,3)
156  IF(mint(55).EQ.7.OR.mint(55).EQ.8) wid2=
157  & wids(mint(55),3)
158  ENDIF
159  IF(i*j.GT.0) sigh(nchn)=facqqs*facckm*wid2
160  IF(i*j.LT.0) sigh(nchn)=facqqu*facckm*wid2
161  ENDIF
162  120 CONTINUE
163  130 CONTINUE
164 
165  ELSEIF(isub.EQ.84) THEN
166 C...g + gamma -> Q + Qbar
167  sqmavg=0.5d0*(sqm3+sqm4)-0.25d0*(sqm3-sqm4)**2/sh
168  thq=-0.5d0*sh*(1d0-be34*cth)
169  uhq=-0.5d0*sh*(1d0+be34*cth)
170  facqq=comfac*as*aem*(kchg(iabs(mint(55)),1)/3d0)**2*
171  & (thq**2+uhq**2+4d0*sqmavg*sh*(1d0-sqmavg*sh/(thq*uhq)))/
172  & (thq*uhq)
173  IF(mstp(35).GE.1) facqq=facqq*pyhfth(sh,sqmavg,0d0)
174  wid2=1d0
175  IF(mint(55).EQ.6) wid2=wids(6,1)
176  IF(mint(55).EQ.7.OR.mint(55).EQ.8) wid2=wids(mint(55),1)
177  facqq=facqq*wid2
178  IF(kfac(1,21)*kfac(2,22).NE.0) THEN
179  nchn=nchn+1
180  isig(nchn,1)=21
181  isig(nchn,2)=22
182  isig(nchn,3)=1
183  sigh(nchn)=facqq
184  ENDIF
185  IF(kfac(1,22)*kfac(2,21).NE.0) THEN
186  nchn=nchn+1
187  isig(nchn,1)=22
188  isig(nchn,2)=21
189  isig(nchn,3)=1
190  sigh(nchn)=facqq
191  ENDIF
192 
193  ELSEIF(isub.EQ.85) THEN
194 C...gamma + gamma -> F + Fbar (heavy fermion, quark or lepton)
195  sqmavg=0.5d0*(sqm3+sqm4)-0.25d0*(sqm3-sqm4)**2/sh
196  thq=-0.5d0*sh*(1d0-be34*cth)
197  uhq=-0.5d0*sh*(1d0+be34*cth)
198  facff=comfac*aem**2*(kchg(iabs(mint(56)),1)/3d0)**4*2d0*
199  & ((1d0-parj(131)*parj(132))*(thq*uhq-sqmavg*sh)*
200  & (uhq**2+thq**2+2d0*sqmavg*sh)+(1d0+parj(131)*parj(132))*
201  & sqmavg*sh**2*(sh-2d0*sqmavg))/(thq*uhq)**2
202  IF(iabs(mint(56)).LT.10) facff=3d0*facff
203  IF(iabs(mint(56)).LT.10.AND.mstp(35).GE.1)
204  & facff=facff*pyhfth(sh,sqmavg,1d0)
205  wid2=1d0
206  IF(mint(56).EQ.6) wid2=wids(6,1)
207  IF(mint(56).EQ.7.OR.mint(56).EQ.8) wid2=wids(mint(56),1)
208  IF(mint(56).EQ.17) wid2=wids(17,1)
209  facff=facff*wid2
210  IF(kfac(1,22)*kfac(2,22).NE.0) THEN
211  nchn=nchn+1
212  isig(nchn,1)=22
213  isig(nchn,2)=22
214  isig(nchn,3)=1
215  sigh(nchn)=facff
216  ENDIF
217 
218  ELSEIF(isub.EQ.86) THEN
219 C...g + g -> J/Psi + g
220  facqqg=comfac*as**3*(5d0/9d0)*parp(38)*sqrt(sqm3)*
221  & (((sh*(sh-sqm3))**2+(th*(th-sqm3))**2+(uh*(uh-sqm3))**2)/
222  & ((th-sqm3)*(uh-sqm3))**2)/(sh-sqm3)**2
223  IF(kfac(1,21)*kfac(2,21).NE.0) THEN
224  nchn=nchn+1
225  isig(nchn,1)=21
226  isig(nchn,2)=21
227  isig(nchn,3)=1
228  sigh(nchn)=facqqg
229  ENDIF
230 
231  ELSEIF(isub.EQ.87) THEN
232 C...g + g -> chi_0c + g
233  pgtw=(sh*th+th*uh+uh*sh)/sh2
234  qgtw=(sh*th*uh)/sh**3
235  rgtw=sqm3/sh
236  facqqg=comfac*as**3*4d0*(parp(39)/sqrt(sqm3))*(1d0/sh)*
237  & (9d0*rgtw**2*pgtw**4*(rgtw**4-2d0*rgtw**2*pgtw+pgtw**2)-
238  & 6d0*rgtw*pgtw**3*qgtw*(2d0*rgtw**4-5d0*rgtw**2*pgtw+pgtw**2)-
239  & pgtw**2*qgtw**2*(rgtw**4+2d0*rgtw**2*pgtw-pgtw**2)+
240  & 2d0*rgtw*pgtw*qgtw**3*(rgtw**2-pgtw)+6d0*rgtw**2*qgtw**4)/
241  & (qgtw*(qgtw-rgtw*pgtw)**4)
242  IF(kfac(1,21)*kfac(2,21).NE.0) THEN
243  nchn=nchn+1
244  isig(nchn,1)=21
245  isig(nchn,2)=21
246  isig(nchn,3)=1
247  sigh(nchn)=facqqg
248  ENDIF
249 
250  ELSEIF(isub.EQ.88) THEN
251 C...g + g -> chi_1c + g
252  pgtw=(sh*th+th*uh+uh*sh)/sh2
253  qgtw=(sh*th*uh)/sh**3
254  rgtw=sqm3/sh
255  facqqg=comfac*as**3*12d0*(parp(39)/sqrt(sqm3))*(1d0/sh)*
256  & pgtw**2*(rgtw*pgtw**2*(rgtw**2-4d0*pgtw)+2d0*qgtw*(-rgtw**4+
257  & 5d0*rgtw**2*pgtw+pgtw**2)-15d0*rgtw*qgtw**2)/
258  & (qgtw-rgtw*pgtw)**4
259  IF(kfac(1,21)*kfac(2,21).NE.0) THEN
260  nchn=nchn+1
261  isig(nchn,1)=21
262  isig(nchn,2)=21
263  isig(nchn,3)=1
264  sigh(nchn)=facqqg
265  ENDIF
266 
267  ELSEIF(isub.EQ.89) THEN
268 C...g + g -> chi_2c + g
269  pgtw=(sh*th+th*uh+uh*sh)/sh2
270  qgtw=(sh*th*uh)/sh**3
271  rgtw=sqm3/sh
272  facqqg=comfac*as**3*4d0*(parp(39)/sqrt(sqm3))*(1d0/sh)*
273  & (12d0*rgtw**2*pgtw**4*(rgtw**4-2d0*rgtw**2*pgtw+pgtw**2)-
274  & 3d0*rgtw*pgtw**3*qgtw*(8d0*rgtw**4-rgtw**2*pgtw+4d0*pgtw**2)+
275  & 2d0*pgtw**2*qgtw**2*(-7d0*rgtw**4+43d0*rgtw**2*pgtw+pgtw**2)+
276  & rgtw*pgtw*qgtw**3*(16d0*rgtw**2-61d0*pgtw)+12d0*rgtw**2*
277  & qgtw**4)/(qgtw*(qgtw-rgtw*pgtw)**4)
278  IF(kfac(1,21)*kfac(2,21).NE.0) THEN
279  nchn=nchn+1
280  isig(nchn,1)=21
281  isig(nchn,2)=21
282  isig(nchn,3)=1
283  sigh(nchn)=facqqg
284  ENDIF
285  ENDIF
286 
287  ELSEIF(isub.LE.200) THEN
288  IF(isub.EQ.104) THEN
289 C...g + g -> chi_c0.
290  kc=pycomp(10441)
291  facbw=comfac*12d0*as**2*parp(39)*pmas(kc,2)/
292  & ((sh-pmas(kc,1)**2)**2+(pmas(kc,1)*pmas(kc,2))**2)
293  IF(abs(sqrt(sh)-pmas(kc,1)).GT.50d0*pmas(kc,2)) facbw=0d0
294  IF(kfac(1,21)*kfac(2,21).NE.0) THEN
295  nchn=nchn+1
296  isig(nchn,1)=21
297  isig(nchn,2)=21
298  isig(nchn,3)=1
299  sigh(nchn)=facbw
300  ENDIF
301 
302  ELSEIF(isub.EQ.105) THEN
303 C...g + g -> chi_c2.
304  kc=pycomp(445)
305  facbw=comfac*16d0*as**2*parp(39)*pmas(kc,2)/
306  & ((sh-pmas(kc,1)**2)**2+(pmas(kc,1)*pmas(kc,2))**2)
307  IF(abs(sqrt(sh)-pmas(kc,1)).GT.50d0*pmas(kc,2)) facbw=0d0
308  IF(kfac(1,21)*kfac(2,21).NE.0) THEN
309  nchn=nchn+1
310  isig(nchn,1)=21
311  isig(nchn,2)=21
312  isig(nchn,3)=1
313  sigh(nchn)=facbw
314  ENDIF
315 
316  ELSEIF(isub.EQ.106) THEN
317 C...g + g -> J/Psi + gamma.
318  eq=kchg(mod(kfpr(isub,1)/10,10),1)/3d0
319  facqqg=comfac*aem*eq**2*as**2*(4d0/3d0)*parp(38)*sqrt(sqm3)*
320  & (((sh*(sh-sqm3))**2+(th*(th-sqm3))**2+(uh*(uh-sqm3))**2)/
321  & ((th-sqm3)*(uh-sqm3))**2)/(sh-sqm3)**2
322  IF(kfac(1,21)*kfac(2,21).NE.0) THEN
323  nchn=nchn+1
324  isig(nchn,1)=21
325  isig(nchn,2)=21
326  isig(nchn,3)=1
327  sigh(nchn)=facqqg
328  ENDIF
329 
330  ELSEIF(isub.EQ.107) THEN
331 C...g + gamma -> J/Psi + g.
332  eq=kchg(mod(kfpr(isub,1)/10,10),1)/3d0
333  facqqg=comfac*aem*eq**2*as**2*(32d0/3d0)*parp(38)*sqrt(sqm3)*
334  & (((sh*(sh-sqm3))**2+(th*(th-sqm3))**2+(uh*(uh-sqm3))**2)/
335  & ((th-sqm3)*(uh-sqm3))**2)/(sh-sqm3)**2
336  IF(kfac(1,21)*kfac(2,22).NE.0) THEN
337  nchn=nchn+1
338  isig(nchn,1)=21
339  isig(nchn,2)=22
340  isig(nchn,3)=1
341  sigh(nchn)=facqqg
342  ENDIF
343  IF(kfac(1,22)*kfac(2,21).NE.0) THEN
344  nchn=nchn+1
345  isig(nchn,1)=22
346  isig(nchn,2)=21
347  isig(nchn,3)=1
348  sigh(nchn)=facqqg
349  ENDIF
350 
351  ELSEIF(isub.EQ.108) THEN
352 C...gamma + gamma -> J/Psi + gamma.
353  eq=kchg(mod(kfpr(isub,1)/10,10),1)/3d0
354  facqqg=comfac*aem**3*eq**6*384d0*parp(38)*sqrt(sqm3)*
355  & (((sh*(sh-sqm3))**2+(th*(th-sqm3))**2+(uh*(uh-sqm3))**2)/
356  & ((th-sqm3)*(uh-sqm3))**2)/(sh-sqm3)**2
357  IF(kfac(1,22)*kfac(2,22).NE.0) THEN
358  nchn=nchn+1
359  isig(nchn,1)=22
360  isig(nchn,2)=22
361  isig(nchn,3)=1
362  sigh(nchn)=facqqg
363  ENDIF
364  ENDIF
365 
366 C...QUARKONIA+++
367 C...Additional code by Stefan Wolf
368  ELSE
369 
370 C...Common code for quarkonium production.
371  shth=sh+th
372  thuh=th+uh
373  uhsh=uh+sh
374  shth2=shth**2
375  thuh2=thuh**2
376  uhsh2=uhsh**2
377  IF ( (isub.GE.421.AND.isub.LE.424).OR.
378  & (isub.GE.431.AND.isub.LE.433)) THEN
379  sqmqq=sqm3
380  ELSEIF((isub.GE.425.AND.isub.LE.430).OR.
381  & (isub.GE.434.AND.isub.LE.439)) THEN
382  sqmqq=sqm4
383  ENDIF
384  sqmqqr=sqrt(sqmqq)
385  IF(mstp(145).EQ.1) THEN
386  IF ( (isub.GE.421.AND.isub.LE.427).OR.
387  & (isub.GE.431.AND.isub.LE.436)) THEN
388  aq=uhsh/(2d0*x(1)) + shth/(2d0*x(2))
389  bq=uhsh/(2d0*x(1)) - shth/(2d0*x(2))
390  atilk1=x(1)*vint(2)/2d0-uhsh/(2d0*sqmqq)*aq
391  atilk2=x(2)*vint(2)/2d0-shth/(2d0*sqmqq)*aq
392  btilk1=-x(1)*vint(2)/2d0-uhsh/(2d0*sqmqq)*bq
393  btilk2=x(2)*vint(2)/2d0-shth/(2d0*sqmqq)*bq
394  ELSEIF( (isub.GE.428.AND.isub.LE.430).OR.
395  & isub.GE.437) THEN
396  aq=shth/(2d0*x(1)) + uhsh/(2d0*x(2))
397  bq=shth/(2d0*x(1)) - uhsh/(2d0*x(2))
398  atilk1=x(1)*vint(2)/2d0-shth/(2d0*sqmqq)*aq
399  atilk2=x(2)*vint(2)/2d0-uhsh/(2d0*sqmqq)*aq
400  btilk1=-x(1)*vint(2)/2d0-shth/(2d0*sqmqq)*bq
401  btilk2=x(2)*vint(2)/2d0-uhsh/(2d0*sqmqq)*bq
402  ENDIF
403  aq2=aq**2
404  bq2=bq**2
405  smqq2=sqmqq*vint(2)
406 C...Polarisation frames
407  IF(mstp(146).EQ.1) THEN
408 C...Recoil frame
409  polh1=sqrt(aq2-smqq2)
410  polh2=sqrt(vint(2)*(aq2-bq2-smqq2))
411  az=-sqmqqr/polh1
412  bz=0d0
413  ax=aq*bq/(polh1*polh2)
414  bx=-polh1/polh2
415  ELSEIF(mstp(146).EQ.2) THEN
416 C...Gottfried Jackson frame
417  polh1=aq+bq
418  polh2=polh1*sqrt(vint(2)*(aq2-bq2-smqq2))
419  az=sqmqqr/polh1
420  bz=az
421  ax=-(bq2+aq*bq+smqq2)/polh2
422  bx=(aq2+aq*bq-smqq2)/polh2
423  ELSEIF(mstp(146).EQ.3) THEN
424 C...Target frame
425  polh1=aq-bq
426  polh2=polh1*sqrt(vint(2)*(aq2-bq2-smqq2))
427  az=-sqmqqr/polh1
428  bz=-az
429  ax=-(bq2-aq*bq+smqq2)/polh2
430  bx=-(aq2-aq*bq-smqq2)/polh2
431  ELSEIF(mstp(146).EQ.4) THEN
432 C...Collins Soper frame
433  polh1=aq2-bq2
434  polh2=sqrt(vint(2)*polh1)
435  az=-bq/polh2
436  bz=aq/polh2
437  ax=-sqmqqr*aq/sqrt(polh1*(polh1-smqq2))
438  bx=sqmqqr*bq/sqrt(polh1*(polh1-smqq2))
439  ENDIF
440 C...Contract EL1(lam) EL2(lam') with K1 and K2 (initial parton momenta)
441  el1k10=az*atilk1+bz*btilk1
442  el1k20=az*atilk2+bz*btilk2
443  el2k10=el1k10
444  el2k20=el1k20
445  el1k11=1d0/sqrt(2d0)*(ax*atilk1+bx*btilk1)
446  el1k21=1d0/sqrt(2d0)*(ax*atilk2+bx*btilk2)
447  el2k11=el1k11
448  el2k21=el1k21
449  ENDIF
450 
451  IF(isub.EQ.421) THEN
452 C...g + g -> QQ~[3S11] + g
453  IF(mstp(145).EQ.0) THEN
454 * FACQQG=COMFAC*PARU(1)*AS**3*(10D0/81D0)*SQMQQR*
455 * & (SH2*THUH2+TH2*UHSH2+UH2*SHTH2)/(SHTH2*THUH2*UHSH2)
456  facqqg=comfac*paru(1)*as**3*(10d0/81d0)*sqmqqr*
457  & (sh2*thuh2+th2*uhsh2+uh2*shth2)/shth2/thuh2/uhsh2
458 * FACQQG=COMFAC*PARU(1)*AS**3*(10D0/81D0)*SQMQQR*
459 * & (SH2/(SHTH2*UHSH2)+TH2/(SHTH2*THUH2)+UH2/(THUH2*UHSH2))
460  ELSE
461  ff=-paru(1)*as**3*(10d0/81d0)*sqmqqr/thuh2/shth2/uhsh2
462  aa=(shth2*uh2+uhsh2*th2+thuh2*sh2)/2d0
463  bb=2d0*(sh2+th2)
464  cc=2d0*(sh2+uh2)
465  dd=2d0*sh2
466  IF(mstp(147).EQ.0) THEN
467  facqqg=-aa+sqmqq*(bb*el1k10*el2k10+cc*el1k20*el2k20
468  & +dd*(el1k10*el2k20+el1k20*el2k10))
469  ELSEIF(mstp(147).EQ.1) THEN
470  facqqg=2d0*(-aa+sqmqq*(bb*el1k11*el2k11+cc*el1k21*el2k21
471  & +dd*(el1k11*el2k21+el1k21*el2k11)))
472  ELSEIF(mstp(147).EQ.3) THEN
473  facqqg=-aa+sqmqq*(bb*el1k10*el2k10+cc*el1k20*el2k20
474  & +dd*(el1k10*el2k20+el1k20*el2k10))
475  ELSEIF(mstp(147).EQ.4) THEN
476  facqqg=-aa+sqmqq*(bb*el1k11*el2k11+cc*el1k21*el2k21
477  & +dd*(el1k11*el2k21+el1k21*el2k11))
478  ELSEIF(mstp(147).EQ.5) THEN
479  facqqg=sqmqq*(bb*el1k11*el2k10+cc*el1k21*el2k20
480  & +dd*(el1k11*el2k20+el1k21*el2k10))
481  ELSEIF(mstp(147).EQ.6) THEN
482  facqqg=sqmqq*(bb*el1k11*el2k11+cc*el1k21*el2k21
483  & +dd*(el1k11*el2k21+el1k21*el2k11))
484  ENDIF
485  facqqg=comfac*ff*facqqg
486  ENDIF
487  IF(kfac(1,21)*kfac(2,21).NE.0) THEN
488  nchn=nchn+1
489  isig(nchn,1)=21
490  isig(nchn,2)=21
491  isig(nchn,3)=1
492  sigh(nchn)=facqqg*parp(ionium+1)
493  ENDIF
494 
495  ELSEIF(isub.EQ.422) THEN
496 C...g + g -> QQ~[3S18] + g
497  IF(mstp(145).EQ.0) THEN
498  facqqg=-comfac*paru(1)*as**3*(1d0/72d0)*
499  & (16d0*sqmqq**2-27d0*(shth2+thuh2+uhsh2))/
500  & (sqmqq*sqmqqr)*
501  & ((sh2*thuh2+th2*uhsh2+uh2*shth2)/shth2/thuh2/uhsh2)
502  ELSE
503  ff=paru(1)*as**3*(16d0*sqmqq**2-27d0*(shth2+thuh2+uhsh2))/
504  & (72d0*sqmqq*sqmqqr*shth2*thuh2*uhsh2)
505  aa=(shth2*uh2+uhsh2*th2+thuh2*sh2)/2d0
506  bb=2d0*(sh2+th2)
507  cc=2d0*(sh2+uh2)
508  dd=2d0*sh2
509  IF(mstp(147).EQ.0) THEN
510  facqqg=-aa+sqmqq*(bb*el1k10*el2k10+cc*el1k20*el2k20
511  & +dd*(el1k10*el2k20+el1k20*el2k10))
512  ELSEIF(mstp(147).EQ.1) THEN
513  facqqg=2d0*(-aa+sqmqq*(bb*el1k11*el2k11+cc*el1k21*el2k21
514  & +dd*(el1k11*el2k21+el1k21*el2k11)))
515  ELSEIF(mstp(147).EQ.3) THEN
516  facqqg=-aa+sqmqq*(bb*el1k10*el2k10+cc*el1k20*el2k20
517  & +dd*(el1k10*el2k20+el1k20*el2k10))
518  ELSEIF(mstp(147).EQ.4) THEN
519  facqqg=-aa+sqmqq*(bb*el1k11*el2k11+cc*el1k21*el2k21
520  & +dd*(el1k11*el2k21+el1k21*el2k11))
521  ELSEIF(mstp(147).EQ.5) THEN
522  facqqg=sqmqq*(bb*el1k11*el2k10+cc*el1k21*el2k20
523  & +dd*(el1k11*el2k20+el1k21*el2k10))
524  ELSEIF(mstp(147).EQ.6) THEN
525  facqqg=sqmqq*(bb*el1k11*el2k11+cc*el1k21*el2k21
526  & +dd*(el1k11*el2k21+el1k21*el2k11))
527  ENDIF
528  facqqg=comfac*ff*facqqg
529  ENDIF
530 C...Split total contribution into different colour flows just like
531 C...in g g -> g g (recalculate kinematics for massless partons).
532  thp=-0.5d0*sh*(1d0-cth)
533  uhp=-0.5d0*sh*(1d0+cth)
534  facgg1=(sh/thp)**2+2d0*sh/thp+3d0+2d0*thp/sh+(thp/sh)**2
535  facgg2=(uhp/sh)**2+2d0*uhp/sh+3d0+2d0*sh/uhp+(sh/uhp)**2
536  facgg3=(thp/uhp)**2+2d0*thp/uhp+3d0+2d0*uhp/thp+(uhp/thp)**2
537  facggs=facgg1+facgg2+facgg3
538  IF(kfac(1,21)*kfac(2,21).NE.0) THEN
539  nchn=nchn+1
540  isig(nchn,1)=21
541  isig(nchn,2)=21
542  isig(nchn,3)=1
543  sigh(nchn)=facqqg*parp(ionium+2)*facgg1/facggs
544  nchn=nchn+1
545  isig(nchn,1)=21
546  isig(nchn,2)=21
547  isig(nchn,3)=2
548  sigh(nchn)=facqqg*parp(ionium+2)*facgg2/facggs
549  nchn=nchn+1
550  isig(nchn,1)=21
551  isig(nchn,2)=21
552  isig(nchn,3)=3
553  sigh(nchn)=facqqg*parp(ionium+2)*facgg3/facggs
554  ENDIF
555 
556  ELSEIF(isub.EQ.423) THEN
557 C...g + g -> QQ~[1S08] + g
558  IF(mstp(145).EQ.0) THEN
559 * FACQQG=COMFAC*PARU(1)*AS**3*(5D0/16D0)*
560 * & (SHTH2*UH2+THUH2*SH2+UHSH2*TH2)/(SQMQQR*SH*TH*UH)*
561 * & (12D0*SQMQQ*SH*TH*UH+SHTH2**2+THUH2**2+UHSH2**2)/
562 * & (SHTH2*THUH2*UHSH2)
563  facqqg=comfac*paru(1)*as**3*(5d0/16d0)*sqmqqr*
564  & (uh2/(thuh2*uhsh2)+sh2/(shth2*uhsh2)+
565  & th2/(shth2*thuh2))*
566  & (12d0+(shth2**2+thuh2**2+uhsh2**2)/(sqmqq*sh*th*uh))
567  ELSE
568  fa=paru(1)*as**3*(5d0/48d0)*sqmqqr*
569  & (uh2/(thuh2*uhsh2)+sh2/(shth2*uhsh2)+
570  & th2/(shth2*thuh2))*
571  & (12d0+(shth2**2+thuh2**2+uhsh2**2)/(sqmqq*sh*th*uh))
572  IF(mstp(147).EQ.0) THEN
573  facqqg=comfac*fa
574  ELSEIF(mstp(147).EQ.1) THEN
575  facqqg=comfac*2d0*fa
576  ELSEIF(mstp(147).EQ.3) THEN
577  facqqg=comfac*fa
578  ELSEIF(mstp(147).EQ.4) THEN
579  facqqg=comfac*fa
580  ELSEIF(mstp(147).EQ.5) THEN
581  facqqg=0d0
582  ELSEIF(mstp(147).EQ.6) THEN
583  facqqg=0d0
584  ENDIF
585  ENDIF
586 C...Split total contribution into different colour flows just like
587 C...in g g -> g g (recalculate kinematics for massless partons).
588  thp=-0.5d0*sh*(1d0-cth)
589  uhp=-0.5d0*sh*(1d0+cth)
590  facgg1=(sh/thp)**2+2d0*sh/thp+3d0+2d0*thp/sh+(thp/sh)**2
591  facgg2=(uhp/sh)**2+2d0*uhp/sh+3d0+2d0*sh/uhp+(sh/uhp)**2
592  facgg3=(thp/uhp)**2+2d0*thp/uhp+3d0+2d0*uhp/thp+(uhp/thp)**2
593  facggs=facgg1+facgg2+facgg3
594  IF(kfac(1,21)*kfac(2,21).NE.0) THEN
595  nchn=nchn+1
596  isig(nchn,1)=21
597  isig(nchn,2)=21
598  isig(nchn,3)=1
599  sigh(nchn)=facqqg*parp(ionium+3)*facgg1/facggs
600  nchn=nchn+1
601  isig(nchn,1)=21
602  isig(nchn,2)=21
603  isig(nchn,3)=2
604  sigh(nchn)=facqqg*parp(ionium+3)*facgg2/facggs
605  nchn=nchn+1
606  isig(nchn,1)=21
607  isig(nchn,2)=21
608  isig(nchn,3)=3
609  sigh(nchn)=facqqg*parp(ionium+3)*facgg3/facggs
610  ENDIF
611 
612  ELSEIF(isub.EQ.424) THEN
613 C...g + g -> QQ~[3PJ8] + g
614  poly=sh2+sh*th+th2
615  IF(mstp(145).EQ.0) THEN
616  facqqg=comfac*5d0*paru(1)*as**3*(3d0*sh*th*shth*poly**4
617  & -sqmqq*poly**2*(7d0*sh**6+36d0*sh**5*th+45d0*sh**4*th2
618  & +28d0*sh**3*th**3+45d0*sh2*th**4+36d0*sh*th**5
619  & +7d0*th**6)
620  & +sqmqq**2*shth*(35d0*sh**8+169d0*sh**7*th
621  & +299d0*sh**6*th2+401d0*sh**5*th**3+418d0*sh**4*th**4
622  & +401d0*sh**3*th**5+299d0*sh2*th**6+169d0*sh*th**7
623  & +35d0*th**8)
624  & -sqmqq**3*(84d0*sh**8+432d0*sh**7*th+905d0*sh**6*th2
625  & +1287d0*sh**5*th**3+1436d0*sh**4*th**4
626  & +1287d0*sh**3*th**5+905d0*sh2*th**6+432d0*sh*th**7
627  & +84d0*th**8)
628  & +sqmqq**4*shth*(126d0*sh**6+451d0*sh**5*th
629  & +677d0*sh**4*th2+836d0*sh**3*th**3+677d0*sh2*th**4
630  & +451d0*sh*th**5+126d0*th**6)
631  & -3d0*sqmqq**5*(42d0*sh**6+171d0*sh**5*th
632  & +304d0*sh**4*th2+362d0*sh**3*th**3+304d0*sh2*th**4
633  & +171d0*sh*th**5+42d0*th**6)
634  & +2d0*sqmqq**6*shth*(42d0*sh**4+106d0*sh**3*th
635  & +119d0*sh2*th2+106d0*sh*th**3+42d0*th**4)
636  & -sqmqq**7*(35d0*sh**4+99d0*sh**3*th+120d0*sh2*th2
637  & +99d0*sh*th**3+35d0*th**4)
638  & +7d0*sqmqq**8*shth*poly)/
639  & (sh*th*uh*sqmqqr*sqmqq*
640  & shth*shth2*thuh*thuh2*uhsh*uhsh2)
641  ELSE
642  ff=-5d0*paru(1)*as**3/(sh2*th2*uh2
643  & *sqmqqr*sqmqq*shth*shth2*thuh*thuh2*uhsh*uhsh2)
644  aa=sh*th*uh*(sh*th*shth*poly**4
645  & -sqmqq*shth2*poly**2*
646  & (sh**4+6d0*sh**3*th-6d0*sh2*th2+6d0*sh*th**3+th**4)
647  & +sqmqq**2*shth*(5d0*sh**8+35d0*sh**7*th+49d0*sh**6*th2
648  & +57d0*sh**5*th**3+46d0*sh**4*th**4+57d0*sh**3*th**5
649  & +49d0*sh2*th**6+35d0*sh*th**7+5d0*th**8)
650  & -sqmqq**3*(16d0*sh**8+104d0*sh**7*th+215d0*sh**6*th2
651  & +291d0*sh**5*th**3+316d0*sh**4*th**4+291d0*sh**3*th**5
652  & +215d0*sh2*th**6+104d0*sh*th**7+16d0*th**8)
653  & +sqmqq**4*shth*(34d0*sh**6+145d0*sh**5*th
654  & +211d0*sh**4*th2+262d0*sh**3*th**3+211d0*sh2*th**4
655  & +145d0*sh*th**5+34d0*th**6)
656  & -sqmqq**5*(44d0*sh**6+193d0*sh**5*th+346d0*sh**4*th2
657  & +410d0*sh**3*th**3+346d0*sh2*th**4+193d0*sh*th**5
658  & +44d0*th**6)
659  & +2d0*sqmqq**6*shth*(17d0*sh**4+45d0*sh**3*th
660  & +49d0*sh2*th2+45d0*sh*th**3+17d0*th**4)
661  & -sqmqq**7*(3d0*sh2+2d0*sh*th+3d0*th2)
662  & *(5d0*sh2+11d0*sh*th+5d0*th2)
663  & +3d0*sqmqq**8*shth*poly)
664  bb=4d0*shth2*poly**3
665  & *(sh**4+sh**3*th-sh2*th2+sh*th**3+th**4)
666  & -sqmqq*shth*(20d0*sh**10+84d0*sh**9*th+166d0*sh**8*th2
667  & +231d0*sh**7*th**3+250d0*sh**6*th**4+250d0*sh**5*th**5
668  & +250d0*sh**4*th**6+231d0*sh**3*th**7+166d0*sh2*th**8
669  & +84d0*sh*th**9+20d0*th**10)
670  & +sqmqq**2*shth2*(40d0*sh**8+86d0*sh**7*th
671  & +66d0*sh**6*th2+67d0*sh**5*th**3+6d0*sh**4*th**4
672  & +67d0*sh**3*th**5+66d0*sh2*th**6+86d0*sh*th**7
673  & +40d0*th**8)
674  & -sqmqq**3*shth*(40d0*sh**8+57d0*sh**7*th
675  & -110d0*sh**6*th2-263d0*sh**5*th**3-384d0*sh**4*th**4
676  & -263d0*sh**3*th**5-110d0*sh2*th**6+57d0*sh*th**7
677  & +40d0*th**8)
678  & +sqmqq**4*(20d0*sh**8-33d0*sh**7*th-368d0*sh**6*th2
679  & -751d0*sh**5*th**3-920d0*sh**4*th**4-751d0*sh**3*th**5
680  & -368d0*sh2*th**6-33d0*sh*th**7+20d0*th**8)
681  & -sqmqq**5*shth*(4d0*sh**6-81d0*sh**5*th-242d0*sh**4*th2
682  & -250d0*sh**3*th**3-242d0*sh2*th**4-81d0*sh*th**5
683  & +4d0*th**6)
684  & -sqmqq**6*sh*th*(41d0*sh**4+120d0*sh**3*th
685  & +142d0*sh2*th2+120d0*sh*th**3+41d0*th**4)
686  & +8d0*sqmqq**7*sh*th*shth*poly
687  cc=4d0*th2*poly**3
688  & *(-sh**4-2d0*sh**3*th+2d0*sh2*th2+3d0*sh*th**3+th**4)
689  & -sqmqq*th2*(-20d0*sh**9-56d0*sh**8*th-24d0*sh**7*th2
690  & +147d0*sh**6*th**3+409d0*sh**5*th**4+599d0*sh**4*th**5
691  & +571d0*sh**3*th**6+370d0*sh2*th**7+148d0*sh*th**8
692  & +28d0*th**9)
693  & +sqmqq**2*(4d0*sh**10+20d0*sh**9*th-16d0*sh**8*th2
694  & -48d0*sh**7*th**3+150d0*sh**6*th**4+611d0*sh**5*th**5
695  & +1060d0*sh**4*th**6+1155d0*sh**3*th**7+854d0*sh2*th**8
696  & +394d0*sh*th**9+84d0*th**10)
697  & -sqmqq**3*shth*(20d0*sh**8+68d0*sh**7*th-20d0*sh**6*th2
698  & +32d0*sh**5*th**3+286d0*sh**4*th**4+577d0*sh**3*th**5
699  & +618d0*sh2*th**6+443d0*sh*th**7+140d0*th**8)
700  & +sqmqq**4*(40d0*sh**8+152d0*sh**7*th+94d0*sh**6*th2
701  & +38d0*sh**5*th**3+290d0*sh**4*th**4+631d0*sh**3*th**5
702  & +738d0*sh2*th**6+513d0*sh*th**7+140d0*th**8)
703  & -sqmqq**5*(40d0*sh**7+129d0*sh**6*th+53d0*sh**5*th2
704  & +7d0*sh**4*th**3+129d0*sh**3*th**4+264d0*sh2*th**5
705  & +266d0*sh*th**6+84d0*th**7)
706  & +sqmqq**6*(20d0*sh**6+55d0*sh**5*th+2d0*sh**4*th2
707  & -15d0*sh**3*th**3+30d0*sh2*th**4+76d0*sh*th**5
708  & +28d0*th**6)
709  & -sqmqq**7*shth*(4d0*sh**4+7d0*sh**3*th-14d0*sh2*th2
710  & +7d0*sh*th**3+4*th**4)
711  & +sqmqq**8*sh*(sh-th)**2*th
712  dd=2d0*th2*shth2*poly**3
713  & *(-sh2+2*sh*th+2*th2)
714  & +sqmqq*(4d0*sh**11+22d0*sh**10*th+70d0*sh**9*th2
715  & +115d0*sh**8*th**3+71d0*sh**7*th**4-119d0*sh**6*th**5
716  & -381d0*sh**5*th**6-552d0*sh**4*th**7-512d0*sh**3*th**8
717  & -320d0*sh2*th**9-126d0*sh*th**10-24d0*th**11)
718  & -sqmqq**2*shth*(20d0*sh**9+84d0*sh**8*th
719  & +212d0*sh**7*th2+247d0*sh**6*th**3+105d0*sh**5*th**4
720  & -178d0*sh**4*th**5-380d0*sh**3*th**6-364d0*sh2*th**7
721  & -210d0*sh*th**8-60d0*th**9)
722  & +sqmqq**3*shth*(40d0*sh**8+159d0*sh**7*th
723  & +374d0*sh**6*th2+404d0*sh**5*th**3+192d0*sh**4*th**4
724  & -141d0*sh**3*th**5-264d0*sh2*th**6-216d0*sh*th**7
725  & -80d0*th**8)
726  & -sqmqq**4*(40d0*sh**8+197d0*sh**7*th+506d0*sh**6*th2
727  & +672d0*sh**5*th**3+460d0*sh**4*th**4+79d0*sh**3*th**5
728  & -138d0*sh2*th**6-164d0*sh*th**7-60d0*th**8)
729  & +sqmqq**5*(20d0*sh**7+107d0*sh**6*th+267d0*sh**5*th2
730  & +307d0*sh**4*th**3+185d0*sh**3*th**4+56d0*sh2*th**5
731  & -30d0*sh*th**6-24d0*th**7)
732  & -sqmqq**6*(4d0*sh**6+31d0*sh**5*th+74d0*sh**4*th2
733  & +71d0*sh**3*th**3+46d0*sh2*th**4+10d0*sh*th**5
734  & -4d0*th**6)
735  & +4d0*sqmqq**7*sh*th*shth*poly
736  IF(mstp(147).EQ.0) THEN
737  facqqg=-aa+sqmqq*(bb*el1k10*el2k10+cc*el1k20*el2k20
738  & +dd*(el1k10*el2k20+el1k20*el2k10))
739  ELSEIF(mstp(147).EQ.1) THEN
740  facqqg=2d0*(-aa+sqmqq*(bb*el1k11*el2k11+cc*el1k21*el2k21
741  & +dd*(el1k11*el2k21+el1k21*el2k11)))
742  ELSEIF(mstp(147).EQ.3) THEN
743  facqqg=-aa+sqmqq*(bb*el1k10*el2k10+cc*el1k20*el2k20
744  & +dd*(el1k10*el2k20+el1k20*el2k10))
745  ELSEIF(mstp(147).EQ.4) THEN
746  facqqg=-aa+sqmqq*(bb*el1k11*el2k11+cc*el1k21*el2k21
747  & +dd*(el1k11*el2k21+el1k21*el2k11))
748  ELSEIF(mstp(147).EQ.5) THEN
749  facqqg=sqmqq*(bb*el1k11*el2k10+cc*el1k21*el2k20
750  & +dd*(el1k11*el2k20+el1k21*el2k10))
751  ELSEIF(mstp(147).EQ.6) THEN
752  facqqg=sqmqq*(bb*el1k11*el2k11+cc*el1k21*el2k21
753  & +dd*(el1k11*el2k21+el1k21*el2k11))
754  ENDIF
755  facqqg=comfac*ff*facqqg
756  ENDIF
757 C...Split total contribution into different colour flows just like
758 C...in g g -> g g (recalculate kinematics for massless partons).
759  thp=-0.5d0*sh*(1d0-cth)
760  uhp=-0.5d0*sh*(1d0+cth)
761  facgg1=(sh/thp)**2+2d0*sh/thp+3d0+2d0*thp/sh+(thp/sh)**2
762  facgg2=(uhp/sh)**2+2d0*uhp/sh+3d0+2d0*sh/uhp+(sh/uhp)**2
763  facgg3=(thp/uhp)**2+2d0*thp/uhp+3d0+2d0*uhp/thp+(uhp/thp)**2
764  facggs=facgg1+facgg2+facgg3
765  IF(kfac(1,21)*kfac(2,21).NE.0) THEN
766  nchn=nchn+1
767  isig(nchn,1)=21
768  isig(nchn,2)=21
769  isig(nchn,3)=1
770  sigh(nchn)=facqqg*parp(ionium+4)*facgg1/facggs
771  nchn=nchn+1
772  isig(nchn,1)=21
773  isig(nchn,2)=21
774  isig(nchn,3)=2
775  sigh(nchn)=facqqg*parp(ionium+4)*facgg2/facggs
776  nchn=nchn+1
777  isig(nchn,1)=21
778  isig(nchn,2)=21
779  isig(nchn,3)=3
780  sigh(nchn)=facqqg*parp(ionium+4)*facgg3/facggs
781  ENDIF
782 
783  ELSEIF(isub.EQ.425) THEN
784 C...q + g -> q + QQ~[3S18]
785  IF(mstp(145).EQ.0) THEN
786  facqqg=-comfac*paru(1)*as**3*(1d0/27d0)*
787  & (4d0*(sh2+uh2)-sh*uh)*(shth2+thuh2)/
788  & (sqmqq*sqmqqr*sh*uh*uhsh2)
789  ELSE
790  ff=paru(1)*as**3*(4d0*(sh2+uh2)-sh*uh)/
791  & (54d0*sqmqq*sqmqqr*sh*uh*uhsh2)
792  aa=shth2+thuh2
793  bb=4d0
794  cc=8d0
795  dd=4d0
796  IF(mstp(147).EQ.0) THEN
797  facqqg=-aa+sqmqq*(bb*el1k10*el2k10+cc*el1k20*el2k20
798  & +dd*(el1k10*el2k20+el1k20*el2k10))
799  ELSEIF(mstp(147).EQ.1) THEN
800  facqqg=2d0*(-aa+sqmqq*(bb*el1k11*el2k11+cc*el1k21*el2k21
801  & +dd*(el1k11*el2k21+el1k21*el2k11)))
802  ELSEIF(mstp(147).EQ.3) THEN
803  facqqg=-aa+sqmqq*(bb*el1k10*el2k10+cc*el1k20*el2k20
804  & +dd*(el1k10*el2k20+el1k20*el2k10))
805  ELSEIF(mstp(147).EQ.4) THEN
806  facqqg=-aa+sqmqq*(bb*el1k11*el2k11+cc*el1k21*el2k21
807  & +dd*(el1k11*el2k21+el1k21*el2k11))
808  ELSEIF(mstp(147).EQ.5) THEN
809  facqqg=sqmqq*(bb*el1k11*el2k10+cc*el1k21*el2k20
810  & +dd*(el1k11*el2k20+el1k21*el2k10))
811  ELSEIF(mstp(147).EQ.6) THEN
812  facqqg=sqmqq*(bb*el1k11*el2k11+cc*el1k21*el2k21
813  & +dd*(el1k11*el2k21+el1k21*el2k11))
814  ENDIF
815  facqqg=comfac*ff*facqqg
816  ENDIF
817 C...Split total contribution into different colour flows just like
818 C...in ISUB.EQ.28 [f + g -> f + g (q + g -> q + g only)]
819 C...(recalculate kinematics for massless partons).
820  thp=-0.5d0*sh*(1d0-cth)
821  uhp=-0.5d0*sh*(1d0+cth)
822  facqg1=9d0/4d0*(uhp/thp)**2-uhp/sh
823  facqg2=9d0/4d0*(sh/thp)**2-sh/uhp
824  facqgs=facqg1+facqg2
825  DO 2442 i=mmina,mmaxa
826  IF(i.EQ.0.OR.iabs(i).GT.mstp(58)) goto 2442
827  DO 2441 isde=1,2
828  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 2441
829  IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 2441
830  nchn=nchn+1
831  isig(nchn,isde)=i
832  isig(nchn,3-isde)=21
833  isig(nchn,3)=1
834  sigh(nchn)=facqqg*parp(ionium+2)*facqg1/facqgs
835  nchn=nchn+1
836  isig(nchn,isde)=i
837  isig(nchn,3-isde)=21
838  isig(nchn,3)=2
839  sigh(nchn)=facqqg*parp(ionium+2)*facqg2/facqgs
840  2441 CONTINUE
841  2442 CONTINUE
842 
843  ELSEIF(isub.EQ.426) THEN
844 C...q + g -> q + QQ~[1S08]
845  IF(mstp(145).EQ.0) THEN
846  facqqg=-comfac*paru(1)*as**3*(5d0/18d0)*
847  & (sh2+uh2)/(sqmqqr*th*uhsh2)
848  ELSE
849  fa=-paru(1)*as**3*(5d0/54d0)*(sh2+uh2)/(sqmqqr*th*uhsh2)
850  IF(mstp(147).EQ.0) THEN
851  facqqg=comfac*fa
852  ELSEIF(mstp(147).EQ.1) THEN
853  facqqg=comfac*2d0*fa
854  ELSEIF(mstp(147).EQ.3) THEN
855  facqqg=comfac*fa
856  ELSEIF(mstp(147).EQ.4) THEN
857  facqqg=comfac*fa
858  ELSEIF(mstp(147).EQ.5) THEN
859  facqqg=0d0
860  ELSEIF(mstp(147).EQ.6) THEN
861  facqqg=0d0
862  ENDIF
863  ENDIF
864 C...Split total contribution into different colour flows just like
865 C...in ISUB.EQ.28 [f + g -> f + g (q + g -> q + g only)]
866 C...(recalculate kinematics for massless partons).
867  thp=-0.5d0*sh*(1d0-cth)
868  uhp=-0.5d0*sh*(1d0+cth)
869  facqg1=9d0/4d0*(uhp/thp)**2-uhp/sh
870  facqg2=9d0/4d0*(sh/thp)**2-sh/uhp
871  facqgs=facqg1+facqg2
872  DO 2444 i=mmina,mmaxa
873  IF(i.EQ.0.OR.iabs(i).GT.mstp(58)) goto 2444
874  DO 2443 isde=1,2
875  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 2443
876  IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 2443
877  nchn=nchn+1
878  isig(nchn,isde)=i
879  isig(nchn,3-isde)=21
880  isig(nchn,3)=1
881  sigh(nchn)=facqqg*parp(ionium+3)*facqg1/facqgs
882  nchn=nchn+1
883  isig(nchn,isde)=i
884  isig(nchn,3-isde)=21
885  isig(nchn,3)=2
886  sigh(nchn)=facqqg*parp(ionium+3)*facqg2/facqgs
887  2443 CONTINUE
888  2444 CONTINUE
889 
890  ELSEIF(isub.EQ.427) THEN
891 C...q + g -> q + QQ~[3PJ8]
892  IF(mstp(145).EQ.0) THEN
893  facqqg=-comfac*paru(1)*as**3*(10d0/9d0)*
894  & ((7d0*uhsh+8d0*th)*(sh2+uh2)
895  & +4d0*th*(2d0*sqmqq**2-shth2-thuh2))/
896  & (sqmqq*sqmqqr*th*uhsh2*uhsh)
897  ELSE
898  ff=10d0*paru(1)*as**3/
899  & (9d0*sqmqq*sqmqqr*th2*uhsh2*uhsh)
900  aa=th*uhsh*(2d0*sqmqq**2+shth2+thuh2)
901  bb=8d0*(shth2+th*uh)
902  cc=8d0*uhsh*(shth+thuh)
903  dd=4d0*(2d0*sqmqq*sh+th*uhsh)
904  IF(mstp(147).EQ.0) THEN
905  facqqg=-aa+sqmqq*(bb*el1k10*el2k10+cc*el1k20*el2k20
906  & +dd*(el1k10*el2k20+el1k20*el2k10))
907  ELSEIF(mstp(147).EQ.1) THEN
908  facqqg=2d0*(-aa+sqmqq*(bb*el1k11*el2k11+cc*el1k21*el2k21
909  & +dd*(el1k11*el2k21+el1k21*el2k11)))
910  ELSEIF(mstp(147).EQ.3) THEN
911  facqqg=-aa+sqmqq*(bb*el1k10*el2k10+cc*el1k20*el2k20
912  & +dd*(el1k10*el2k20+el1k20*el2k10))
913  ELSEIF(mstp(147).EQ.4) THEN
914  facqqg=-aa+sqmqq*(bb*el1k11*el2k11+cc*el1k21*el2k21
915  & +dd*(el1k11*el2k21+el1k21*el2k11))
916  ELSEIF(mstp(147).EQ.5) THEN
917  facqqg=sqmqq*(bb*el1k11*el2k10+cc*el1k21*el2k20
918  & +dd*(el1k11*el2k20+el1k21*el2k10))
919  ELSEIF(mstp(147).EQ.6) THEN
920  facqqg=sqmqq*(bb*el1k11*el2k11+cc*el1k21*el2k21
921  & +dd*(el1k11*el2k21+el1k21*el2k11))
922  ENDIF
923  facqqg=comfac*ff*facqqg
924  ENDIF
925 C...Split total contribution into different colour flows just like
926 C...in ISUB.EQ.28 [f + g -> f + g (q + g -> q + g only)]
927 C...(recalculate kinematics for massless partons).
928  thp=-0.5d0*sh*(1d0-cth)
929  uhp=-0.5d0*sh*(1d0+cth)
930  facqg1=9d0/4d0*(uhp/thp)**2-uhp/sh
931  facqg2=9d0/4d0*(sh/thp)**2-sh/uhp
932  facqgs=facqg1+facqg2
933  DO 2446 i=mmina,mmaxa
934  IF(i.EQ.0.OR.iabs(i).GT.mstp(58)) goto 2446
935  DO 2445 isde=1,2
936  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 2445
937  IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 2445
938  nchn=nchn+1
939  isig(nchn,isde)=i
940  isig(nchn,3-isde)=21
941  isig(nchn,3)=1
942  sigh(nchn)=facqqg*parp(ionium+4)*facqg1/facqgs
943  nchn=nchn+1
944  isig(nchn,isde)=i
945  isig(nchn,3-isde)=21
946  isig(nchn,3)=2
947  sigh(nchn)=facqqg*parp(ionium+4)*facqg2/facqgs
948  2445 CONTINUE
949  2446 CONTINUE
950 
951  ELSEIF(isub.EQ.428) THEN
952 C...q + q~ -> g + QQ~[3S18]
953  IF(mstp(145).EQ.0) THEN
954  facqqg=comfac*paru(1)*as**3*(8d0/81d0)*
955  & (4d0*(th2+uh2)-th*uh)*(shth2+uhsh2)/
956  & (sqmqq*sqmqqr*th*uh*thuh2)
957  ELSE
958  ff=-4d0*paru(1)*as**3*(4d0*(th2+uh2)-th*uh)/
959  & (81d0*sqmqq*sqmqqr*th*uh*thuh2)
960  aa=shth2+uhsh2
961  bb=4d0
962  cc=4d0
963  dd=0d0
964  IF(mstp(147).EQ.0) THEN
965  facqqg=-aa+sqmqq*(bb*el1k10*el2k10+cc*el1k20*el2k20
966  & +dd*(el1k10*el2k20+el1k20*el2k10))
967  ELSEIF(mstp(147).EQ.1) THEN
968  facqqg=2d0*(-aa+sqmqq*(bb*el1k11*el2k11+cc*el1k21*el2k21
969  & +dd*(el1k11*el2k21+el1k21*el2k11)))
970  ELSEIF(mstp(147).EQ.3) THEN
971  facqqg=-aa+sqmqq*(bb*el1k10*el2k10+cc*el1k20*el2k20
972  & +dd*(el1k10*el2k20+el1k20*el2k10))
973  ELSEIF(mstp(147).EQ.4) THEN
974  facqqg=-aa+sqmqq*(bb*el1k11*el2k11+cc*el1k21*el2k21
975  & +dd*(el1k11*el2k21+el1k21*el2k11))
976  ELSEIF(mstp(147).EQ.5) THEN
977  facqqg=sqmqq*(bb*el1k11*el2k10+cc*el1k21*el2k20
978  & +dd*(el1k11*el2k20+el1k21*el2k10))
979  ELSEIF(mstp(147).EQ.6) THEN
980  facqqg=sqmqq*(bb*el1k11*el2k11+cc*el1k21*el2k21
981  & +dd*(el1k11*el2k21+el1k21*el2k11))
982  ENDIF
983  facqqg=comfac*ff*facqqg
984  ENDIF
985 C...Split total contribution into different colour flows just like
986 C...in ISUB.EQ.13 [f + fbar -> g + g (q + qbar -> g + g only)]
987 C...(recalculate kinematics for massless partons).
988  thp=-0.5d0*sh*(1d0-cth)
989  uhp=-0.5d0*sh*(1d0+cth)
990  facgg1=uh/th-9d0/4d0*uh2/sh2
991  facgg2=th/uh-9d0/4d0*th2/sh2
992  facggs=facgg1+facgg2
993  DO 2447 i=mmina,mmaxa
994  IF(i.EQ.0.OR.iabs(i).GT.mstp(58).OR.
995  & kfac(1,i)*kfac(2,-i).EQ.0) goto 2447
996  nchn=nchn+1
997  isig(nchn,1)=i
998  isig(nchn,2)=-i
999  isig(nchn,3)=1
1000  sigh(nchn)=facqqg*parp(ionium+2)*facgg1/facggs
1001  nchn=nchn+1
1002  isig(nchn,1)=i
1003  isig(nchn,2)=-i
1004  isig(nchn,3)=2
1005  sigh(nchn)=facqqg*parp(ionium+2)*facgg2/facggs
1006  2447 CONTINUE
1007 
1008  ELSEIF(isub.EQ.429) THEN
1009 C...q + q~ -> g + QQ~[1S08]
1010  IF(mstp(145).EQ.0) THEN
1011  facqqg=comfac*paru(1)*as**3*(20d0/27d0)*
1012  & (th2+uh2)/(sqmqqr*sh*thuh2)
1013  ELSE
1014  fa=paru(1)*as**3*(20d0/81d0)*(th2+uh2)/(sqmqqr*sh*thuh2)
1015  IF(mstp(147).EQ.0) THEN
1016  facqqg=comfac*fa
1017  ELSEIF(mstp(147).EQ.1) THEN
1018  facqqg=comfac*2d0*fa
1019  ELSEIF(mstp(147).EQ.3) THEN
1020  facqqg=comfac*fa
1021  ELSEIF(mstp(147).EQ.4) THEN
1022  facqqg=comfac*fa
1023  ELSEIF(mstp(147).EQ.5) THEN
1024  facqqg=0d0
1025  ELSEIF(mstp(147).EQ.6) THEN
1026  facqqg=0d0
1027  ENDIF
1028  ENDIF
1029 C...Split total contribution into different colour flows just like
1030 C...in ISUB.EQ.13 [f + fbar -> g + g (q + qbar -> g + g only)]
1031 C...(recalculate kinematics for massless partons).
1032  thp=-0.5d0*sh*(1d0-cth)
1033  uhp=-0.5d0*sh*(1d0+cth)
1034  facgg1=uh/th-9d0/4d0*uh2/sh2
1035  facgg2=th/uh-9d0/4d0*th2/sh2
1036  facggs=facgg1+facgg2
1037  DO 2448 i=mmina,mmaxa
1038  IF(i.EQ.0.OR.iabs(i).GT.mstp(58).OR.
1039  & kfac(1,i)*kfac(2,-i).EQ.0) goto 2448
1040  nchn=nchn+1
1041  isig(nchn,1)=i
1042  isig(nchn,2)=-i
1043  isig(nchn,3)=1
1044  sigh(nchn)=facqqg*parp(ionium+3)*facgg1/facggs
1045  nchn=nchn+1
1046  isig(nchn,1)=i
1047  isig(nchn,2)=-i
1048  isig(nchn,3)=2
1049  sigh(nchn)=facqqg*parp(ionium+3)*facgg2/facggs
1050  2448 CONTINUE
1051 
1052  ELSEIF(isub.EQ.430) THEN
1053 C...q + q~ -> g + QQ~[3PJ8]
1054  IF(mstp(145).EQ.0) THEN
1055  facqqg=comfac*paru(1)*as**3*(80d0/27d0)*
1056  & ((7d0*thuh+8d0*sh)*(th2+uh2)
1057  & +4d0*sh*(2d0*sqmqq**2-shth2-uhsh2))/
1058  & (sqmqq*sqmqqr*sh*thuh2*thuh)
1059  ELSE
1060  ff=-80d0*paru(1)*as**3/(27d0*sqmqq*sqmqqr*sh2*thuh2*thuh)
1061  aa=sh*thuh*(2d0*sqmqq**2+shth2+uhsh2)
1062  bb=8d0*(uhsh2+sh*th)
1063  cc=8d0*(shth2+sh*uh)
1064  dd=4d0*(shth2+uhsh2+sh*sqmqq-sqmqq**2)
1065  IF(mstp(147).EQ.0) THEN
1066  facqqg=-aa+sqmqq*(bb*el1k10*el2k10+cc*el1k20*el2k20
1067  & +dd*(el1k10*el2k20+el1k20*el2k10))
1068  ELSEIF(mstp(147).EQ.1) THEN
1069  facqqg=2d0*(-aa+sqmqq*(bb*el1k11*el2k11+cc*el1k21*el2k21
1070  & +dd*(el1k11*el2k21+el1k21*el2k11)))
1071  ELSEIF(mstp(147).EQ.3) THEN
1072  facqqg=-aa+sqmqq*(bb*el1k10*el2k10+cc*el1k20*el2k20
1073  & +dd*(el1k10*el2k20+el1k20*el2k10))
1074  ELSEIF(mstp(147).EQ.4) THEN
1075  facqqg=-aa+sqmqq*(bb*el1k11*el2k11+cc*el1k21*el2k21
1076  & +dd*(el1k11*el2k21+el1k21*el2k11))
1077  ELSEIF(mstp(147).EQ.5) THEN
1078  facqqg=sqmqq*(bb*el1k11*el2k10+cc*el1k21*el2k20
1079  & +dd*(el1k11*el2k20+el1k21*el2k10))
1080  ELSEIF(mstp(147).EQ.6) THEN
1081  facqqg=sqmqq*(bb*el1k11*el2k11+cc*el1k21*el2k21
1082  & +dd*(el1k11*el2k21+el1k21*el2k11))
1083  ENDIF
1084  facqqg=comfac*ff*facqqg
1085  ENDIF
1086 C...Split total contribution into different colour flows just like
1087 C...in ISUB.EQ.13 [f + fbar -> g + g (q + qbar -> g + g only)]
1088 C...(recalculate kinematics for massless partons).
1089  thp=-0.5d0*sh*(1d0-cth)
1090  uhp=-0.5d0*sh*(1d0+cth)
1091  facgg1=uh/th-9d0/4d0*uh2/sh2
1092  facgg2=th/uh-9d0/4d0*th2/sh2
1093  facggs=facgg1+facgg2
1094  DO 2449 i=mmina,mmaxa
1095  IF(i.EQ.0.OR.iabs(i).GT.mstp(58).OR.
1096  & kfac(1,i)*kfac(2,-i).EQ.0) goto 2449
1097  nchn=nchn+1
1098  isig(nchn,1)=i
1099  isig(nchn,2)=-i
1100  isig(nchn,3)=1
1101  sigh(nchn)=facqqg*parp(ionium+4)*facgg1/facggs
1102  nchn=nchn+1
1103  isig(nchn,1)=i
1104  isig(nchn,2)=-i
1105  isig(nchn,3)=2
1106  sigh(nchn)=facqqg*parp(ionium+4)*facgg2/facggs
1107  2449 CONTINUE
1108 
1109  ELSEIF(isub.EQ.431) THEN
1110 C...g + g -> QQ~[3P01] + g
1111  pgtw=(sh*th+th*uh+uh*sh)/sh2
1112  qgtw=(sh*th*uh)/sh**3
1113  rgtw=sqmqq/sh
1114  IF(mstp(145).EQ.0) THEN
1115  facqqg=comfac*paru(1)*as**3*8d0/(9d0*sqmqqr*sh)*
1116  & (9d0*rgtw**2*pgtw**4*
1117  & (rgtw**4-2d0*rgtw**2*pgtw+pgtw**2)
1118  & -6d0*rgtw*pgtw**3*qgtw*
1119  & (2d0*rgtw**4-5d0*rgtw**2*pgtw+pgtw**2)
1120  & -pgtw**2*qgtw**2*(rgtw**4+2d0*rgtw**2*pgtw-pgtw**2)
1121  & +2d0*rgtw*pgtw*qgtw**3*(rgtw**2-pgtw)
1122  & +6d0*rgtw**2*qgtw**4)/(qgtw*(qgtw-rgtw*pgtw)**4)
1123  ELSE
1124  fc1=paru(1)*as**3*8d0/(27d0*sqmqqr*sh)*
1125  & (9d0*rgtw**2*pgtw**4*
1126  & (rgtw**4-2d0*rgtw**2*pgtw+pgtw**2)
1127  & -6d0*rgtw*pgtw**3*qgtw*
1128  & (2d0*rgtw**4-5d0*rgtw**2*pgtw+pgtw**2)
1129  & -pgtw**2*qgtw**2*(rgtw**4+2d0*rgtw**2*pgtw-pgtw**2)
1130  & +2d0*rgtw*pgtw*qgtw**3*(rgtw**2-pgtw)
1131  & +6d0*rgtw**2*qgtw**4)/(qgtw*(qgtw-rgtw*pgtw)**4)
1132  IF(mstp(147).EQ.0) THEN
1133  facqqg=comfac*fc1
1134  ELSEIF(mstp(147).EQ.1) THEN
1135  facqqg=comfac*2d0*fc1
1136  ELSEIF(mstp(147).EQ.3) THEN
1137  facqqg=comfac*fc1
1138  ELSEIF(mstp(147).EQ.4) THEN
1139  facqqg=comfac*fc1
1140  ELSEIF(mstp(147).EQ.5) THEN
1141  facqqg=0d0
1142  ELSEIF(mstp(147).EQ.6) THEN
1143  facqqg=0d0
1144  ENDIF
1145  ENDIF
1146  IF(kfac(1,21)*kfac(2,21).NE.0) THEN
1147  nchn=nchn+1
1148  isig(nchn,1)=21
1149  isig(nchn,2)=21
1150  isig(nchn,3)=1
1151  sigh(nchn)=facqqg*parp(ionium+5)
1152  ENDIF
1153 
1154  ELSEIF(isub.EQ.432) THEN
1155 C...g + g -> QQ~[3P11] + g
1156  pgtw=(sh*th+th*uh+uh*sh)/sh2
1157  qgtw=(sh*th*uh)/sh**3
1158  rgtw=sqmqq/sh
1159  IF(mstp(145).EQ.0) THEN
1160  facqqg=comfac*paru(1)*as**3*8d0/(3d0*sqmqqr*sh)*
1161  & pgtw**2*(rgtw*pgtw**2*(rgtw**2-4d0*pgtw)
1162  & +2d0*qgtw*(-rgtw**4+5d0*rgtw**2*pgtw+pgtw**2)
1163  & -15d0*rgtw*qgtw**2)/(qgtw-rgtw*pgtw)**4
1164  ELSE
1165  ff=4d0/3d0*paru(1)*as**3*sqmqqr/shth2**2/thuh2**2/uhsh2**2
1166  c1=(4d0*pgtw**5+23d0*pgtw**2*qgtw**2
1167  & +(-14d0*pgtw**3*qgtw+3d0*qgtw**3)*rgtw
1168  & -(pgtw**4+2d0*pgtw*qgtw**2)*rgtw**2
1169  & +3d0*pgtw**2*qgtw*rgtw**3)*sh2**5
1170  c2=2d0*shth2*(sh2*thuh*(sh*thuh*(sh-th)*(sh-uh)
1171  & -th*uh*(th-uh)**2)+sh2**2*(th-uh)*(th2+uh2-sh*thuh)
1172  & *(pgtw**2-qgtw*(sh+2d0*uh)/sh))
1173  c3=2d0*uhsh2*(sh2*thuh*(sh*thuh*(sh-th)*(sh-uh)
1174  & -th*uh*(th-uh)**2)-sh2**2*(th-uh)*(th2+uh2-sh*thuh)
1175  & *(pgtw**2-qgtw*(sh+2d0*th)/sh))
1176  c4=-4d0*thuh*(th-uh)**2*
1177  & (th**3*uh**3+sh2**2*(2d0*th+uh)*(th+2d0*uh)
1178  & -sh2*th*uh*(th2+uh2))
1179  & +4d0*thuh2*(sh**3*(sh2**2+th2**2+uh2**2)
1180  & -sh*th*uh*(sh2**2+th*uh*(th2-3d0*th*uh+uh2)
1181  & +sh2*(5d0*thuh2-17d0*th*uh)))
1182  IF(mstp(147).EQ.0) THEN
1183  facqqg=-c1+c2*el1k10*el2k10+c3*el1k20*el2k20
1184  & +c4*(el1k10*el2k20+el1k20*el2k10)/2d0
1185  ELSEIF(mstp(147).EQ.1) THEN
1186  facqqg=2d0*(-c1+c2*el1k11*el2k11+c3*el1k21*el2k21
1187  & +c4*(el1k11*el2k21+el1k21*el2k11)/2d0)
1188  ELSEIF(mstp(147).EQ.3) THEN
1189  facqqg=-c1+c2*el1k10*el2k10+c3*el1k20*el2k20
1190  & +c4*(el1k10*el2k20+el1k20*el2k10)/2d0
1191  ELSEIF(mstp(147).EQ.4) THEN
1192  facqqg=-c1+c2*el1k11*el2k11+c3*el1k21*el2k21
1193  & +c4*(el1k11*el2k21+el1k21*el2k11)/2d0
1194  ELSEIF(mstp(147).EQ.5) THEN
1195  facqqg=c2*el1k11*el2k10+c3*el1k21*el2k20
1196  & +c4*(el1k11*el2k20+el1k21*el2k10)/2d0
1197  ELSEIF(mstp(147).EQ.6) THEN
1198  facqqg=c2*el1k11*el2k11+c3*el1k21*el2k21
1199  & +c4*(el1k11*el2k21+el1k21*el2k11)/2d0
1200  ENDIF
1201  facqqg=comfac*ff*facqqg
1202  ENDIF
1203  IF(kfac(1,21)*kfac(2,21).NE.0) THEN
1204  nchn=nchn+1
1205  isig(nchn,1)=21
1206  isig(nchn,2)=21
1207  isig(nchn,3)=1
1208  sigh(nchn)=facqqg*parp(ionium+5)
1209  ENDIF
1210 
1211  ELSEIF(isub.EQ.433) THEN
1212 C...g + g -> QQ~[3P21] + g
1213  pgtw=(sh*th+th*uh+uh*sh)/sh2
1214  qgtw=(sh*th*uh)/sh**3
1215  rgtw=sqmqq/sh
1216  IF(mstp(145).EQ.0) THEN
1217  facqqg=comfac*paru(1)*as**3*8d0/(9d0*sqmqqr*sh)*
1218  & (12d0*rgtw**2*pgtw**4*
1219  & (rgtw**4-2d0*rgtw**2*pgtw+pgtw**2)
1220  & -3d0*rgtw*pgtw**3*qgtw*
1221  & (8d0*rgtw**4-rgtw**2*pgtw+4d0*pgtw**2)
1222  & +2d0*pgtw**2*qgtw**2*
1223  & (-7d0*rgtw**4+43d0*rgtw**2*pgtw+pgtw**2)
1224  & +rgtw*pgtw*qgtw**3*(16d0*rgtw**2-61d0*pgtw)
1225  & +12d0*rgtw**2*qgtw**4)/(qgtw*(qgtw-rgtw*pgtw)**4)
1226  ELSE
1227  ff=(16d0*paru(1)*as**3*sqmqq*sqmqqr)/
1228  & (3d0*sh2*th2*uh2*shth2**2*thuh2**2*uhsh2**2)
1229  c1=pgtw**2*qgtw*(pgtw*rgtw-qgtw)**2*(rgtw**2-2d0*pgtw)
1230  & *sh*sh2**7
1231  c2=2d0*shth2*(-sh2**3*th2**3-sh**5*th**5*uh*shth
1232  & +sh2**2*th2**2*uh2*(8d0*shth2-5d0*sh*th)
1233  & +sh**3*th**3*uh**3*shth*(17d0*shth2-2d0*sh*th)
1234  & +sh2*th2*uh2**2*(105d0*sh2*th2+64d0*sh*th*(sh2+th2)
1235  & +10d0*(sh2**2+th2**2))
1236  & +sh2*th2*uh**5*shth*(32d0*shth2+7d0*sh*th)
1237  & -uh2**3*(sh2**3-87d0*sh**3*th**3+th2**3
1238  & -45d0*sh2*th2*(sh2+th2)-5d0*sh*th*(sh2**2+th2**2))
1239  & +sh*th*uh**7*shth*(7d0*shth2+12d0*sh*th)
1240  & +4d0*sh*th*uh2**4*shth2)
1241  c3=2d0*uhsh2*(-sh2**3*uh2**3-sh**5*uh**5*th*uhsh
1242  & +sh2**2*uh2**2*th2*(8d0*uhsh2-5d0*sh*uh)
1243  & +sh**3*uh**3*th**3*uhsh*(17d0*uhsh2-2d0*sh*uh)
1244  & +sh2*uh2*th2**2*(105d0*sh2*uh2+64d0*sh*uh*(sh2+uh2)
1245  & +10d0*(sh2**2+uh2**2))
1246  & +sh2*uh2*th**5*uhsh*(32d0*uhsh2+7d0*sh*uh)
1247  & -th2**3*(sh2**3-87d0*sh**3*uh**3+uh2**3
1248  & -45d0*sh2*uh2*(sh2+uh2)-5d0*sh*uh*(sh2**2+uh2**2))
1249  & +sh*uh*th**7*uhsh*(7d0*uhsh2+12d0*sh*uh)
1250  & +4d0*sh*uh*th2**4*uhsh2)
1251  c4=-2d0*shth*uhsh*(-2d0*th2**3*uh2**3
1252  & -sh**5*th2*uh2*thuh*(5d0*th+3d0*uh)*(3d0*th+5d0*uh)
1253  & +sh2**3*(2d0*th+uh)*(th+2d0*uh)*(th2-uh2)**2
1254  & -sh*th2**2*uh2**2*thuh*(5d0*thuh2-4d0*th*uh)
1255  & -sh2*th**3*uh**3*thuh2*(13d0*thuh2-16d0*th*uh)
1256  & -sh**3*th2*uh2*(92d0*th2*uh2*thuh
1257  & +53d0*th*uh*(th**3+uh**3)+11d0*(th**5+uh**5))
1258  & -sh2**2*th*uh*(114d0*th**3*uh**3
1259  & +83d0*th2*uh2*(th2+uh2)+28d0*th*uh*(th2**2+uh2**2)
1260  & +3d0*(th2**3+uh2**3)))
1261  c5=4d0*sh*th*uh2*shth2*(2d0*sh*th+sh*uh+th*uh)**2
1262  & *(2d0*uh*sqmqq**2+shth*(sh*th-uh2))
1263  c6=4d0*sh*uh*th2*uhsh2*(2d0*sh*uh+sh*th+th*uh)**2
1264  & *(2d0*th*sqmqq**2+uhsh*(sh*uh-th2))
1265  c7=4d0*sh*th*uh2*shth*(sh2**2*th**3*(11d0*sh+16d0*th)
1266  & +sh**3*th2*uh*(31d0*sh2+83d0*sh*th+61d0*th2)
1267  & +sh2*th*uh2*(19d0*sh**3+110d0*sh2*th+156d0*sh*th2+
1268  & 82d0*th**3)
1269  & +sh*th*uh**3*(43d0*sh**3+132d0*sh2*th+124d0*sh*th2
1270  & +45d0*th**3)
1271  & +th*uh2**2*(37d0*sh**3+68d0*sh2*th+43d0*sh*th2+
1272  & 8d0*th**3)
1273  & +th*uh**5*(11d0*sh2+13d0*sh*th+5d0*th2)
1274  & +sh**3*uh**3*(3d0*uhsh2-2d0*sh*uh)
1275  & +th**5*uhsh*(5d0*uhsh2+2d0*sh*uh))
1276  c8=4d0*sh*uh*th2*uhsh*(sh2**2*uh**3*(11d0*sh+16d0*uh)
1277  & +sh**3*uh2*th*(31d0*sh2+83d0*sh*uh+61d0*uh2)
1278  & +sh2*uh*th2*(19d0*sh**3+110d0*sh2*uh+156d0*sh*uh2+
1279  & 82d0*uh**3)
1280  & +sh*uh*th**3*(43d0*sh**3+132d0*sh2*uh+124d0*sh*uh2
1281  & +45d0*uh**3)
1282  & +uh*th2**2*(37d0*sh**3+68d0*sh2*uh+43d0*sh*uh2+
1283  & 8d0*uh**3)
1284  & +uh*th**5*(11d0*sh2+13d0*sh*uh+5d0*uh2)
1285  & +sh**3*th**3*(3d0*shth2-2d0*sh*th)
1286  & +uh**5*shth*(5d0*shth2+2d0*sh*th))
1287  c9=4d0*shth*uhsh*(2d0*th**5*uh**5*thuh
1288  & +4d0*sh*th2**2*uh2**2*thuh2
1289  & -sh2*th**3*uh**3*thuh*(th2+uh2)
1290  & -2d0*sh**3*th2*uh2*(thuh2**2+2d0*th*uh*thuh2-th2*uh2)
1291  & +sh2**2*th*uh*thuh*(-th*uh*thuh2+3d0*(th2**2+uh2**2))
1292  & +sh**5*(4d0*th2*uh2*(thuh2-th*uh)
1293  & +5d0*th*uh*(th2**2+uh2**2)+2d0*(th2**3+uh2**3)))
1294  c0=-4d0*(2d0*th2**3*uh2**3*sqmqq
1295  & -sh2*th2**2*uh2**2*thuh*(19d0*thuh2-4d0*th*uh)
1296  & -sh**3*th**3*uh**3*thuh2*(32d0*thuh2+29d0*th*uh)
1297  & -sh2**2*th2*uh2*thuh*(264d0*th2*uh2
1298  & +136d0*th*uh*(th2+uh2)+15d0*(th2**2+uh2**2))
1299  & +sh**5*th*uh*(-428d0*th**3*uh**3
1300  & -256d0*th2*uh2*(th2+uh2)-43d0*th*uh*(th2**2+uh2**2)
1301  & +2d0*(th2**3+uh2**3))
1302  & +sh**7*(-46d0*th**3*uh**3-21d0*th2*uh2*(th2+uh2)
1303  & +2d0*th*uh*(th2**2+uh2**2)+2d0*(th2**3+uh2**3))
1304  & +sh2**3*thuh*(-134*th**3*uh**3-53d0*th2*uh2*(th2+uh2)
1305  & +4d0*th*uh*(th2**2+uh2**2)+2d0*(th2**3+uh2**3)))
1306  IF(mstp(147).EQ.0) THEN
1307  facqqg=1d0/3d0*(c1*3d0
1308  & -c2*(2d0*el1k10*el2k10+el1k11*el2k11)
1309  & -c3*(2d0*el1k20*el2k20+el1k21*el2k21)
1310  & -c4*(2d0*el1k10*el2k20+el1k11*el2k21)
1311  & +c5*2d0*(el1k10*el2k10-el1k11*el2k11)**2
1312  & +c6*2d0*(el1k20*el2k20-el1k21*el2k21)**2
1313  & +c7*2d0*(el1k10*el2k10-el1k11*el2k11)
1314  & *(el1k10*el2k20-el1k11*el2k21)
1315  & +c8*2d0*(el1k20*el2k20-el1k21*el2k21)
1316  & *(el1k10*el2k20-el1k11*el2k21)
1317  & +c9*2d0*(el1k10*el2k10-el1k11*el2k11)
1318  & *(el1k20*el2k20-el1k21*el2k21)
1319  & +c0*2d0*(el1k10*el2k20-el1k11*el2k21)**2)
1320  ELSEIF(mstp(147).EQ.1) THEN
1321  facqqg=c1*2d0
1322  & -c2*(el1k10*el2k10+el1k11*el2k11)
1323  & -c3*(el1k20*el2k20+el1k21*el2k21)
1324  & -c4*(el1k10*el2k20+el1k11*el2k21)
1325  & +c5*4d0*el1k10*el2k10*el1k11*el2k11
1326  & +c6*4d0*el1k20*el2k20*el1k21*el2k21
1327  & +c7*2d0*(el1k10*el2k10*el1k11*el2k21
1328  & +el1k10*el2k20*el1k11*el2k11)
1329  & +c8*2d0*(el1k20*el2k20*el1k11*el2k21
1330  & +el1k10*el2k20*el1k21*el2k21)
1331  & +c9*4d0*el1k10*el2k20*el1k11*el2k21
1332  & +c0*(el1k10*el2k10*el1k21*el2k21
1333  & +2d0*el1k10*el2k20*el1k11*el2k21
1334  & +el1k20*el2k20*el1k11*el2k11)
1335  ELSEIF(mstp(147).EQ.2) THEN
1336  facqqg=2d0*(c1
1337  & -c2*el1k11*el2k11
1338  & -c3*el1k21*el2k21
1339  & -c4*el1k11*el2k21
1340  & +c5*(el1k11*el2k11)**2
1341  & +c6*(el1k21*el2k21)**2
1342  & +c7*el1k11*el2k11*el1k11*el2k21
1343  & +c8*el1k21*el2k21*el1k11*el2k21
1344  & +(c9+c0)*(el1k11*el2k21)**2)
1345  ENDIF
1346  facqqg=comfac*ff*facqqg
1347  ENDIF
1348  IF(kfac(1,21)*kfac(2,21).NE.0) THEN
1349  nchn=nchn+1
1350  isig(nchn,1)=21
1351  isig(nchn,2)=21
1352  isig(nchn,3)=1
1353  sigh(nchn)=facqqg*parp(ionium+5)
1354  ENDIF
1355 
1356  ELSEIF(isub.EQ.434) THEN
1357 C...q + g -> q + QQ~[3P01]
1358  IF(mstp(145).EQ.0) THEN
1359  facqqg=-comfac*paru(1)*as**3*(16d0/81d0)*
1360  & (th-3d0*sqmqq)**2*(sh2+uh2)/(sqmqqr*th*uhsh2**2)
1361  ELSE
1362  fa=-paru(1)*as**3*(16d0/243d0)*
1363  & (th-3d0*sqmqq)**2*(sh2+uh2)/(sqmqqr*th*uhsh2**2)
1364  IF(mstp(147).EQ.0) THEN
1365  facqqg=comfac*fa
1366  ELSEIF(mstp(147).EQ.1) THEN
1367  facqqg=comfac*2d0*fa
1368  ELSEIF(mstp(147).EQ.3) THEN
1369  facqqg=comfac*fa
1370  ELSEIF(mstp(147).EQ.4) THEN
1371  facqqg=comfac*fa
1372  ELSEIF(mstp(147).EQ.5) THEN
1373  facqqg=0d0
1374  ELSEIF(mstp(147).EQ.6) THEN
1375  facqqg=0d0
1376  ENDIF
1377  ENDIF
1378  DO 2452 i=mmina,mmaxa
1379  IF(i.EQ.0.OR.iabs(i).GT.mstp(58)) goto 2452
1380  DO 2451 isde=1,2
1381  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 2451
1382  IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 2451
1383  nchn=nchn+1
1384  isig(nchn,isde)=i
1385  isig(nchn,3-isde)=21
1386  isig(nchn,3)=1
1387  sigh(nchn)=facqqg*parp(ionium+5)
1388  2451 CONTINUE
1389  2452 CONTINUE
1390 
1391  ELSEIF(isub.EQ.435) THEN
1392 C...q + g -> q + QQ~[3P11]
1393  IF(mstp(145).EQ.0) THEN
1394  facqqg=-comfac*paru(1)*as**3*(32d0/27d0)*
1395  & (4d0*sqmqq*sh*uh+th*(sh2+uh2))/(sqmqqr*uhsh2**2)
1396  ELSE
1397  ff=(64d0*paru(1)*as**3*sqmqqr)/(27d0*uhsh2**2)
1398  c1=sh*uh
1399  c2=2d0*sh
1400  c3=0d0
1401  c4=2d0*(sh-uh)
1402  IF(mstp(147).EQ.0) THEN
1403  facqqg=-c1+c2*el1k10*el2k10+c3*el1k20*el2k20
1404  & +c4*(el1k10*el2k20+el1k20*el2k10)/2d0
1405  ELSEIF(mstp(147).EQ.1) THEN
1406  facqqg=2d0*(-c1+c2*el1k11*el2k11+c3*el1k21*el2k21
1407  & +c4*(el1k11*el2k21+el1k21*el2k11)/2d0)
1408  ELSEIF(mstp(147).EQ.3) THEN
1409  facqqg=-c1+c2*el1k10*el2k10+c3*el1k20*el2k20
1410  & +c4*(el1k10*el2k20+el1k20*el2k10)/2d0
1411  ELSEIF(mstp(147).EQ.4) THEN
1412  facqqg=-c1+c2*el1k11*el2k11+c3*el1k21*el2k21
1413  & +c4*(el1k11*el2k21+el1k21*el2k11)/2d0
1414  ELSEIF(mstp(147).EQ.5) THEN
1415  facqqg=c2*el1k11*el2k10+c3*el1k21*el2k20
1416  & +c4*(el1k11*el2k20+el1k21*el2k10)/2d0
1417  ELSEIF(mstp(147).EQ.6) THEN
1418  facqqg=c2*el1k11*el2k11+c3*el1k21*el2k21
1419  & +c4*(el1k11*el2k21+el1k21*el2k11)/2d0
1420  ENDIF
1421  facqqg=comfac*ff*facqqg
1422  ENDIF
1423  DO 2454 i=mmina,mmaxa
1424  IF(i.EQ.0.OR.iabs(i).GT.mstp(58)) goto 2454
1425  DO 2453 isde=1,2
1426  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 2453
1427  IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 2453
1428  nchn=nchn+1
1429  isig(nchn,isde)=i
1430  isig(nchn,3-isde)=21
1431  isig(nchn,3)=1
1432  sigh(nchn)=facqqg*parp(ionium+5)
1433  2453 CONTINUE
1434  2454 CONTINUE
1435 
1436  ELSEIF(isub.EQ.436) THEN
1437 C...q + g -> q + QQ~[3P21]
1438  IF(mstp(145).EQ.0) THEN
1439  facqqg=-comfac*paru(1)*as**3*(32d0/81d0)*
1440  & ((6d0*sqmqq**2+th2)*uhsh2
1441  & -2d0*sh*uh*(th2+6d0*sqmqq*uhsh))/
1442  & (sqmqqr*th*uhsh2**2)
1443  ELSE
1444  ff=-(32d0*paru(1)*as**3*sqmqq*sqmqqr)/(27d0*th2*uhsh2**2)
1445  c1=th*uhsh2
1446  c2=4d0*(sh2+th2+2d0*th*uhsh)
1447  c3=4d0*uhsh2
1448  c4=8d0*sh*uhsh
1449  c5=8d0*th
1450  c6=0d0
1451  c7=16d0*th
1452  c8=0d0
1453  c9=-16d0*uhsh
1454  c0=16d0*sqmqq
1455  IF(mstp(147).EQ.0) THEN
1456  facqqg=1d0/3d0*(c1*3d0
1457  & -c2*(2d0*el1k10*el2k10+el1k11*el2k11)
1458  & -c3*(2d0*el1k20*el2k20+el1k21*el2k21)
1459  & -c4*(2d0*el1k10*el2k20+el1k11*el2k21)
1460  & +c5*2d0*(el1k10*el2k10-el1k11*el2k11)**2
1461  & +c6*2d0*(el1k20*el2k20-el1k21*el2k21)**2
1462  & +c7*2d0*(el1k10*el2k10-el1k11*el2k11)
1463  & *(el1k10*el2k20-el1k11*el2k21)
1464  & +c8*2d0*(el1k20*el2k20-el1k21*el2k21)
1465  & *(el1k10*el2k20-el1k11*el2k21)
1466  & +c9*2d0*(el1k10*el2k10-el1k11*el2k11)
1467  & *(el1k20*el2k20-el1k21*el2k21)
1468  & +c0*2d0*(el1k10*el2k20-el1k11*el2k21)**2)
1469  ELSEIF(mstp(147).EQ.1) THEN
1470  facqqg=c1*2d0
1471  & -c2*(el1k10*el2k10+el1k11*el2k11)
1472  & -c3*(el1k20*el2k20+el1k21*el2k21)
1473  & -c4*(el1k10*el2k20+el1k11*el2k21)
1474  & +c5*4d0*el1k10*el2k10*el1k11*el2k11
1475  & +c6*4d0*el1k20*el2k20*el1k21*el2k21
1476  & +c7*2d0*(el1k10*el2k10*el1k11*el2k21
1477  & +el1k10*el2k20*el1k11*el2k11)
1478  & +c8*2d0*(el1k20*el2k20*el1k11*el2k21
1479  & +el1k10*el2k20*el1k21*el2k21)
1480  & +c9*4d0*el1k10*el2k20*el1k11*el2k21
1481  & +c0*(el1k10*el2k10*el1k21*el2k21
1482  & +2d0*el1k10*el2k20*el1k11*el2k21
1483  & +el1k20*el2k20*el1k11*el2k11)
1484  ELSEIF(mstp(147).EQ.2) THEN
1485  facqqg=2d0*(c1
1486  & -c2*el1k11*el2k11
1487  & -c3*el1k21*el2k21
1488  & -c4*el1k11*el2k21
1489  & +c5*(el1k11*el2k11)**2
1490  & +c6*(el1k21*el2k21)**2
1491  & +c7*el1k11*el2k11*el1k11*el2k21
1492  & +c8*el1k21*el2k21*el1k11*el2k21
1493  & +(c9+c0)*(el1k11*el2k21)**2)
1494  ENDIF
1495  facqqg=comfac*ff*facqqg
1496  ENDIF
1497  DO 2456 i=mmina,mmaxa
1498  IF(i.EQ.0.OR.iabs(i).GT.mstp(58)) goto 2456
1499  DO 2455 isde=1,2
1500  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 2455
1501  IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 2455
1502  nchn=nchn+1
1503  isig(nchn,isde)=i
1504  isig(nchn,3-isde)=21
1505  isig(nchn,3)=1
1506  sigh(nchn)=facqqg*parp(ionium+5)
1507  2455 CONTINUE
1508  2456 CONTINUE
1509 
1510  ELSEIF(isub.EQ.437) THEN
1511 C...q + q~ -> g + QQ~[3P01]
1512  IF(mstp(145).EQ.0) THEN
1513  facqqg=comfac*paru(1)*as**3*(128d0/243d0)*
1514  & (sh-3d0*sqmqq)**2*(th2+uh2)/(sqmqqr*sh*thuh2**2)
1515  ELSE
1516  fa=paru(1)*as**3*(128d0/729d0)*
1517  & (sh-3d0*sqmqq)**2*(th2+uh2)/(sqmqqr*sh*thuh2**2)
1518  IF(mstp(147).EQ.0) THEN
1519  facqqg=comfac*fa
1520  ELSEIF(mstp(147).EQ.1) THEN
1521  facqqg=comfac*2d0*fa
1522  ELSEIF(mstp(147).EQ.3) THEN
1523  facqqg=comfac*fa
1524  ELSEIF(mstp(147).EQ.4) THEN
1525  facqqg=comfac*fa
1526  ELSEIF(mstp(147).EQ.5) THEN
1527  facqqg=0d0
1528  ELSEIF(mstp(147).EQ.6) THEN
1529  facqqg=0d0
1530  ENDIF
1531  ENDIF
1532  DO 2457 i=mmina,mmaxa
1533  IF(i.EQ.0.OR.iabs(i).GT.mstp(58).OR.
1534  & kfac(1,i)*kfac(2,-i).EQ.0) goto 2457
1535  nchn=nchn+1
1536  isig(nchn,1)=i
1537  isig(nchn,2)=-i
1538  isig(nchn,3)=1
1539  sigh(nchn)=facqqg*parp(ionium+5)
1540  2457 CONTINUE
1541 
1542  ELSEIF(isub.EQ.438) THEN
1543 C...q + q~ -> g + QQ~[3P11]
1544  IF(mstp(145).EQ.0) THEN
1545  facqqg=comfac*paru(1)*as**3*256d0/81d0*
1546  & (4d0*sqmqq*th*uh+sh*(th2+uh2))/(sqmqqr*thuh2**2)
1547  ELSE
1548  ff=-(512d0*paru(1)*as**3*sqmqqr)/(81d0*thuh2**2)
1549  c1=th*uh
1550  c2=2d0*uh
1551  c3=2d0*th
1552  c4=2d0*thuh
1553  IF(mstp(147).EQ.0) THEN
1554  facqqg=-c1+c2*el1k10*el2k10+c3*el1k20*el2k20
1555  & +c4*(el1k10*el2k20+el1k20*el2k10)/2d0
1556  ELSEIF(mstp(147).EQ.1) THEN
1557  facqqg=2d0*(-c1+c2*el1k11*el2k11+c3*el1k21*el2k21
1558  & +c4*(el1k11*el2k21+el1k21*el2k11)/2d0)
1559  ELSEIF(mstp(147).EQ.3) THEN
1560  facqqg=-c1+c2*el1k10*el2k10+c3*el1k20*el2k20
1561  & +c4*(el1k10*el2k20+el1k20*el2k10)/2d0
1562  ELSEIF(mstp(147).EQ.4) THEN
1563  facqqg=-c1+c2*el1k11*el2k11+c3*el1k21*el2k21
1564  & +c4*(el1k11*el2k21+el1k21*el2k11)/2d0
1565  ELSEIF(mstp(147).EQ.5) THEN
1566  facqqg=c2*el1k11*el2k10+c3*el1k21*el2k20
1567  & +c4*(el1k11*el2k20+el1k21*el2k10)/2d0
1568  ELSEIF(mstp(147).EQ.6) THEN
1569  facqqg=c2*el1k11*el2k11+c3*el1k21*el2k21
1570  & +c4*(el1k11*el2k21+el1k21*el2k11)/2d0
1571  ENDIF
1572  facqqg=comfac*ff*facqqg
1573  ENDIF
1574  DO 2458 i=mmina,mmaxa
1575  IF(i.EQ.0.OR.iabs(i).GT.mstp(58).OR.
1576  & kfac(1,i)*kfac(2,-i).EQ.0) goto 2458
1577  nchn=nchn+1
1578  isig(nchn,1)=i
1579  isig(nchn,2)=-i
1580  isig(nchn,3)=1
1581  sigh(nchn)=facqqg*parp(ionium+5)
1582  2458 CONTINUE
1583 
1584  ELSEIF(isub.EQ.439) THEN
1585 C...q + q~ -> g + QQ~[3P21]
1586  IF(mstp(145).EQ.0) THEN
1587  facqqg=comfac*paru(1)*as**3*(256d0/243d0)*
1588  & ((6d0*sqmqq**2+sh2)*thuh2
1589  & -2d0*th*uh*(sh2+6d0*sqmqq*thuh))/
1590  & (sqmqqr*sh*thuh2**2)
1591  ELSE
1592  ff=(256d0*paru(1)*as**3*sqmqq*sqmqqr)/(81d0*sh2*thuh2**2)
1593  c1=sh*thuh2
1594  c2=4d0*(sh2+uh2+2d0*sh*thuh)
1595  c3=4d0*(sh2+th2+2d0*sh*thuh)
1596  c4=8d0*(sh2-th*uh+2d0*sh*thuh)
1597  c5=8d0*sh
1598  c6=c5
1599  c7=16d0*sh
1600  c8=c7
1601  c9=-16d0*thuh
1602  c0=16d0*sqmqq
1603  IF(mstp(147).EQ.0) THEN
1604  facqqg=1d0/3d0*(c1*3d0
1605  & -c2*(2d0*el1k10*el2k10+el1k11*el2k11)
1606  & -c3*(2d0*el1k20*el2k20+el1k21*el2k21)
1607  & -c4*(2d0*el1k10*el2k20+el1k11*el2k21)
1608  & +c5*2d0*(el1k10*el2k10-el1k11*el2k11)**2
1609  & +c6*2d0*(el1k20*el2k20-el1k21*el2k21)**2
1610  & +c7*2d0*(el1k10*el2k10-el1k11*el2k11)
1611  & *(el1k10*el2k20-el1k11*el2k21)
1612  & +c8*2d0*(el1k20*el2k20-el1k21*el2k21)
1613  & *(el1k10*el2k20-el1k11*el2k21)
1614  & +c9*2d0*(el1k10*el2k10-el1k11*el2k11)
1615  & *(el1k20*el2k20-el1k21*el2k21)
1616  & +c0*2d0*(el1k10*el2k20-el1k11*el2k21)**2)
1617  ELSEIF(mstp(147).EQ.1) THEN
1618  facqqg=c1*2d0
1619  & -c2*(el1k10*el2k10+el1k11*el2k11)
1620  & -c3*(el1k20*el2k20+el1k21*el2k21)
1621  & -c4*(el1k10*el2k20+el1k11*el2k21)
1622  & +c5*4d0*el1k10*el2k10*el1k11*el2k11
1623  & +c6*4d0*el1k20*el2k20*el1k21*el2k21
1624  & +c7*2d0*(el1k10*el2k10*el1k11*el2k21
1625  & +el1k10*el2k20*el1k11*el2k11)
1626  & +c8*2d0*(el1k20*el2k20*el1k11*el2k21
1627  & +el1k10*el2k20*el1k21*el2k21)
1628  & +c9*4d0*el1k10*el2k20*el1k11*el2k21
1629  & +c0*(el1k10*el2k10*el1k21*el2k21
1630  & +2d0*el1k10*el2k20*el1k11*el2k21
1631  & +el1k20*el2k20*el1k11*el2k11)
1632  ELSEIF(mstp(147).EQ.2) THEN
1633  facqqg=2d0*(c1
1634  & -c2*el1k11*el2k11
1635  & -c3*el1k21*el2k21
1636  & -c4*el1k11*el2k21
1637  & +c5*(el1k11*el2k11)**2
1638  & +c6*(el1k21*el2k21)**2
1639  & +c7*el1k11*el2k11*el1k11*el2k21
1640  & +c8*el1k21*el2k21*el1k11*el2k21
1641  & +(c9+c0)*(el1k11*el2k21)**2)
1642  ENDIF
1643  facqqg=comfac*ff*facqqg
1644  ENDIF
1645  DO 2459 i=mmina,mmaxa
1646  IF(i.EQ.0.OR.iabs(i).GT.mstp(58).OR.
1647  & kfac(1,i)*kfac(2,-i).EQ.0) goto 2459
1648  nchn=nchn+1
1649  isig(nchn,1)=i
1650  isig(nchn,2)=-i
1651  isig(nchn,3)=1
1652  sigh(nchn)=facqqg*parp(ionium+5)
1653  2459 CONTINUE
1654  ENDIF
1655 C...QUARKONIA---
1656 
1657  ENDIF
1658 
1659  RETURN
1660  END