Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pysgex.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pysgex.f
1 
2 C*********************************************************************
3 
4 C...PYSGEX
5 C...Subprocess cross sections for assorted exotic processes,
6 C...including Z'/W'/LQ/R/f*/H++/Z_R/W_R/G*.
7 C...Auxiliary to PYSIGH.
8 
9  SUBROUTINE pysgex(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/pytcsm/itcm(0:99),rtcm(0:99)
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/,/pytcsm/,/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.160) THEN
40  IF(isub.EQ.141) THEN
41 C...f + fbar -> gamma*/Z0/Z'0
42  sqmzp=pmas(32,1)**2
43  mint(61)=2
44  CALL pywidt(32,sh,wdtp,wdte)
45  hp0=aem/3d0*sh
46  hp1=aem/3d0*xwc*sh
47  hp2=hp1
48  hs=shr*vint(117)
49  hsp=shr*wdtp(0)
50  faczp=4d0*comfac*3d0
51  DO 100 i=mmina,mmaxa
52  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 100
53  ei=kchg(iabs(i),1)/3d0
54  ai=sign(1d0,ei)
55  vi=ai-4d0*ei*xwv
56  ia=iabs(i)
57  IF(ia.LT.10) THEN
58  IF(ia.LE.2) THEN
59  vpi=paru(123-2*mod(iabs(i),2))
60  api=paru(124-2*mod(iabs(i),2))
61  ELSEIF(ia.LE.4) THEN
62  vpi=parj(182-2*mod(iabs(i),2))
63  api=parj(183-2*mod(iabs(i),2))
64  ELSE
65  vpi=parj(190-2*mod(iabs(i),2))
66  api=parj(191-2*mod(iabs(i),2))
67  ENDIF
68  ELSE
69  IF(ia.LE.12) THEN
70  vpi=paru(127-2*mod(iabs(i),2))
71  api=paru(128-2*mod(iabs(i),2))
72  ELSEIF(ia.LE.14) THEN
73  vpi=parj(186-2*mod(iabs(i),2))
74  api=parj(187-2*mod(iabs(i),2))
75  ELSE
76  vpi=parj(194-2*mod(iabs(i),2))
77  api=parj(195-2*mod(iabs(i),2))
78  ENDIF
79  ENDIF
80  hi0=hp0
81  IF(iabs(i).LE.10) hi0=hi0*faca/3d0
82  hi1=hp1
83  IF(iabs(i).LE.10) hi1=hi1*faca/3d0
84  hi2=hp2
85  IF(iabs(i).LE.10) hi2=hi2*faca/3d0
86  nchn=nchn+1
87  isig(nchn,1)=i
88  isig(nchn,2)=-i
89  isig(nchn,3)=1
90  sigh(nchn)=faczp*(ei**2/sh2*hi0*hp0*vint(111)+ei*vi*
91  & (1d0-sqmz/sh)/((sh-sqmz)**2+hs**2)*(hi0*hp1+hi1*hp0)*
92  & vint(112)+ei*vpi*(1d0-sqmzp/sh)/((sh-sqmzp)**2+hsp**2)*
93  & (hi0*hp2+hi2*hp0)*vint(113)+(vi**2+ai**2)/
94  & ((sh-sqmz)**2+hs**2)*hi1*hp1*vint(114)+(vi*vpi+ai*api)*
95  & ((sh-sqmz)*(sh-sqmzp)+hs*hsp)/(((sh-sqmz)**2+hs**2)*
96  & ((sh-sqmzp)**2+hsp**2))*(hi1*hp2+hi2*hp1)*vint(115)+
97  & (vpi**2+api**2)/((sh-sqmzp)**2+hsp**2)*hi2*hp2*vint(116))
98  100 CONTINUE
99 
100  ELSEIF(isub.EQ.142) THEN
101 C...f + fbar' -> W'+/-
102  sqmwp=pmas(34,1)**2
103  CALL pywidt(34,sh,wdtp,wdte)
104  hs=shr*wdtp(0)
105  facbw=4d0*comfac/((sh-sqmwp)**2+hs**2)*3d0
106  hp=aem/(24d0*xw)*sh
107  DO 120 i=mmin1,mmax1
108  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 120
109  ia=iabs(i)
110  DO 110 j=mmin2,mmax2
111  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 110
112  ja=iabs(j)
113  IF(i*j.GT.0.OR.mod(ia+ja,2).EQ.0) goto 110
114  IF((ia.LE.10.AND.ja.GT.10).OR.(ia.GT.10.AND.ja.LE.10))
115  & goto 110
116  kchw=(kchg(ia,1)*isign(1,i)+kchg(ja,1)*isign(1,j))/3
117  hi=hp*(paru(133)**2+paru(134)**2)
118  IF(ia.LE.10) hi=hp*(paru(131)**2+paru(132)**2)*
119  & vckm((ia+1)/2,(ja+1)/2)*faca/3d0
120  nchn=nchn+1
121  isig(nchn,1)=i
122  isig(nchn,2)=j
123  isig(nchn,3)=1
124  hf=shr*(wdte(0,1)+wdte(0,(5-kchw)/2)+wdte(0,4))
125  sigh(nchn)=hi*facbw*hf
126  110 CONTINUE
127  120 CONTINUE
128 
129  ELSEIF(isub.EQ.144) THEN
130 C...f + fbar' -> R
131  sqmr=pmas(41,1)**2
132  CALL pywidt(41,sh,wdtp,wdte)
133  hs=shr*wdtp(0)
134  facbw=4d0*comfac/((sh-sqmr)**2+hs**2)*3d0
135  hp=aem/(12d0*xw)*sh
136  DO 140 i=mmin1,mmax1
137  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 140
138  ia=iabs(i)
139  DO 130 j=mmin2,mmax2
140  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 130
141  ja=iabs(j)
142  IF(i*j.GT.0.OR.iabs(ia-ja).NE.2) goto 130
143  hi=hp
144  IF(ia.LE.10) hi=hi*faca/3d0
145  hf=shr*(wdte(0,1)+wdte(0,(10-(i+j))/4)+wdte(0,4))
146  nchn=nchn+1
147  isig(nchn,1)=i
148  isig(nchn,2)=j
149  isig(nchn,3)=1
150  sigh(nchn)=hi*facbw*hf
151  130 CONTINUE
152  140 CONTINUE
153 
154  ELSEIF(isub.EQ.145) THEN
155 C...q + l -> LQ (leptoquark)
156  sqmlq=pmas(42,1)**2
157  CALL pywidt(42,sh,wdtp,wdte)
158  hs=shr*wdtp(0)
159  facbw=4d0*comfac/((sh-sqmlq)**2+hs**2)
160  IF(abs(shr-pmas(42,1)).GT.parp(48)*pmas(42,2)) facbw=0d0
161  hp=aem/4d0*sh
162  kflqq=kfdp(mdcy(42,2),1)
163  kflql=kfdp(mdcy(42,2),2)
164  DO 160 i=mmin1,mmax1
165  IF(kfac(1,i).EQ.0) goto 160
166  ia=iabs(i)
167  IF(ia.NE.kflqq.AND.ia.NE.iabs(kflql)) goto 160
168  DO 150 j=mmin2,mmax2
169  IF(kfac(2,j).EQ.0) goto 150
170  ja=iabs(j)
171  IF(ja.NE.kflqq.AND.ja.NE.iabs(kflql)) goto 150
172  IF(i*j.NE.kflqq*kflql) goto 150
173  IF(ja.EQ.ia) goto 150
174  IF(ia.EQ.kflqq) kchlq=isign(1,i)
175  IF(ja.EQ.kflqq) kchlq=isign(1,j)
176  hi=hp*paru(151)
177  hf=shr*(wdte(0,1)+wdte(0,(5-kchlq)/2)+wdte(0,4))
178  nchn=nchn+1
179  isig(nchn,1)=i
180  isig(nchn,2)=j
181  isig(nchn,3)=1
182  sigh(nchn)=hi*facbw*hf
183  150 CONTINUE
184  160 CONTINUE
185 
186  ELSEIF(isub.EQ.146) THEN
187 C...e + gamma* -> e* (excited lepton)
188  kfqstr=kfpr(isub,1)
189  kcqstr=pycomp(kfqstr)
190  kfqexc=mod(kfqstr,kexcit)
191  CALL pywidt(kfqstr,sh,wdtp,wdte)
192  hs=shr*wdtp(0)
193  facbw=comfac/((sh-pmas(kcqstr,1)**2)**2+hs**2)
194  qf=-rtcm(43)/2d0-rtcm(44)/2d0
195  facbw=facbw*aem*qf**2*sh/rtcm(41)**2
196  IF(abs(shr-pmas(kcqstr,1)).GT.parp(48)*pmas(kcqstr,2))
197  & facbw=0d0
198  hp=sh
199  DO 180 i=-kfqexc,kfqexc,2*kfqexc
200  DO 170 isde=1,2
201  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,22).EQ.0) goto 170
202  IF(isde.EQ.2.AND.kfac(1,22)*kfac(2,i).EQ.0) goto 170
203  hi=hp
204  IF(i.GT.0) hf=shr*(wdte(0,1)+wdte(0,2)+wdte(0,4))
205  IF(i.LT.0) hf=shr*(wdte(0,1)+wdte(0,3)+wdte(0,4))
206  nchn=nchn+1
207  isig(nchn,isde)=i
208  isig(nchn,3-isde)=22
209  isig(nchn,3)=1
210  sigh(nchn)=hi*facbw*hf
211  170 CONTINUE
212  180 CONTINUE
213 
214  ELSEIF(isub.EQ.147.OR.isub.EQ.148) THEN
215 C...d + g -> d* and u + g -> u* (excited quarks)
216  kfqstr=kfpr(isub,1)
217  kcqstr=pycomp(kfqstr)
218  kfqexc=mod(kfqstr,kexcit)
219  CALL pywidt(kfqstr,sh,wdtp,wdte)
220  hs=shr*wdtp(0)
221  facbw=comfac/((sh-pmas(kcqstr,1)**2)**2+hs**2)
222  facbw=facbw*as*rtcm(45)**2*sh/(3d0*rtcm(41)**2)
223  IF(abs(shr-pmas(kcqstr,1)).GT.parp(48)*pmas(kcqstr,2))
224  & facbw=0d0
225  hp=sh
226  DO 200 i=-kfqexc,kfqexc,2*kfqexc
227  DO 190 isde=1,2
228  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 190
229  IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 190
230  hi=hp
231  IF(i.GT.0) hf=shr*(wdte(0,1)+wdte(0,2)+wdte(0,4))
232  IF(i.LT.0) hf=shr*(wdte(0,1)+wdte(0,3)+wdte(0,4))
233  nchn=nchn+1
234  isig(nchn,isde)=i
235  isig(nchn,3-isde)=21
236  isig(nchn,3)=1
237  sigh(nchn)=hi*facbw*hf
238  190 CONTINUE
239  200 CONTINUE
240  ENDIF
241 
242  ELSEIF(isub.LE.190) THEN
243  IF(isub.EQ.162) THEN
244 C...q + g -> LQ + lbar; LQ=leptoquark
245  sqmlq=pmas(42,1)**2
246  faclq=comfac*faca*paru(151)*(as*aem/6d0)*(-th/sh)*
247  & (uh2+sqmlq**2)/(uh-sqmlq)**2
248  kflqq=kfdp(mdcy(42,2),1)
249  DO 220 i=mmina,mmaxa
250  IF(iabs(i).NE.kflqq) goto 220
251  kchlq=isign(1,i)
252  DO 210 isde=1,2
253  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 210
254  IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 210
255  nchn=nchn+1
256  isig(nchn,isde)=i
257  isig(nchn,3-isde)=21
258  isig(nchn,3)=1
259  sigh(nchn)=faclq*wids(42,(5-kchlq)/2)
260  210 CONTINUE
261  220 CONTINUE
262 
263  ELSEIF(isub.EQ.163) THEN
264 C...g + g -> LQ + LQbar; LQ=leptoquark
265  sqmlq=pmas(42,1)**2
266  faclq=comfac*faca*wids(42,1)*(as**2/2d0)*
267  & (7d0/48d0+3d0*(uh-th)**2/(16d0*sh2))*(1d0+2d0*sqmlq*th/
268  & (th-sqmlq)**2+2d0*sqmlq*uh/(uh-sqmlq)**2+4d0*sqmlq**2/
269  & ((th-sqmlq)*(uh-sqmlq)))
270  IF(kfac(1,21)*kfac(2,21).EQ.0) goto 230
271  nchn=nchn+1
272  isig(nchn,1)=21
273  isig(nchn,2)=21
274 C...Since don't know proper colour flow, randomize between alternatives
275  isig(nchn,3)=int(1.5d0+pyr(0))
276  sigh(nchn)=faclq
277  230 CONTINUE
278 
279  ELSEIF(isub.EQ.164) THEN
280 C...q + qbar -> LQ + LQbar; LQ=leptoquark
281  delta=0.25d0*(sqm3-sqm4)**2/sh
282  sqmlq=0.5d0*(sqm3+sqm4)-delta
283  th=th-delta
284  uh=uh-delta
285 C SQMLQ=PMAS(42,1)**2
286  faclqa=comfac*wids(42,1)*(as**2/9d0)*
287  & (sh*(sh-4d0*sqmlq)-(uh-th)**2)/sh2
288  faclqs=comfac*wids(42,1)*((paru(151)**2*aem**2/8d0)*
289  & (-sh*th-(sqmlq-th)**2)/th2+(paru(151)*aem*as/18d0)*
290  & ((sqmlq-th)*(uh-th)+sh*(sqmlq+th))/(sh*th))
291  kflqq=kfdp(mdcy(42,2),1)
292  DO 240 i=mmina,mmaxa
293  IF(i.EQ.0.OR.iabs(i).GT.mstp(58).OR.
294  & kfac(1,i)*kfac(2,-i).EQ.0) goto 240
295  nchn=nchn+1
296  isig(nchn,1)=i
297  isig(nchn,2)=-i
298  isig(nchn,3)=1
299  sigh(nchn)=faclqa
300  IF(iabs(i).EQ.kflqq) sigh(nchn)=faclqa+faclqs
301  240 CONTINUE
302 
303  ELSEIF(isub.EQ.167.OR.isub.EQ.168) THEN
304 C...q + q' -> q" + d* and q + q' -> q" + u* (excited quarks)
305  kfqstr=kfpr(isub,2)
306  kcqstr=pycomp(kfqstr)
307  kfqexc=mod(kfqstr,kexcit)
308  facqsa=comfac*(sh/rtcm(41)**2)**2*(1d0-sqm4/sh)
309  facqsb=comfac*0.25d0*(sh/rtcm(41)**2)**2*(1d0-sqm4/sh)*
310  & (1d0+sqm4/sh)*(1d0+cth)*(1d0+((sh-sqm4)/(sh+sqm4))*cth)
311 C...Propagators: as simulated in PYOFSH and as desired
312  gmmq=pmas(kcqstr,1)*pmas(kcqstr,2)
313  hbw4=gmmq/((sqm4-pmas(kcqstr,1)**2)**2+gmmq**2)
314  CALL pywidt(kfqstr,sqm4,wdtp,wdte)
315  gmmqc=sqrt(sqm4)*wdtp(0)
316  hbw4c=gmmqc/((sqm4-pmas(kcqstr,1)**2)**2+gmmqc**2)
317  facqsa=facqsa*hbw4c/hbw4
318  facqsb=facqsb*hbw4c/hbw4
319 C...Branching ratios.
320  brpos=(wdte(0,1)+wdte(0,2)+wdte(0,4))/wdtp(0)
321  brneg=(wdte(0,1)+wdte(0,3)+wdte(0,4))/wdtp(0)
322  DO 260 i=mmin1,mmax1
323  ia=iabs(i)
324  IF(i.EQ.0.OR.ia.GT.6.OR.kfac(1,i).EQ.0) goto 260
325  DO 250 j=mmin2,mmax2
326  ja=iabs(j)
327  IF(j.EQ.0.OR.ja.GT.6.OR.kfac(2,j).EQ.0) goto 250
328  IF(ia.EQ.kfqexc.AND.i.EQ.j) THEN
329  nchn=nchn+1
330  isig(nchn,1)=i
331  isig(nchn,2)=j
332  isig(nchn,3)=1
333  IF(i.GT.0) sigh(nchn)=(4d0/3d0)*facqsa*brpos
334  IF(i.LT.0) sigh(nchn)=(4d0/3d0)*facqsa*brneg
335  nchn=nchn+1
336  isig(nchn,1)=i
337  isig(nchn,2)=j
338  isig(nchn,3)=2
339  IF(j.GT.0) sigh(nchn)=(4d0/3d0)*facqsa*brpos
340  IF(j.LT.0) sigh(nchn)=(4d0/3d0)*facqsa*brneg
341  ELSEIF((ia.EQ.kfqexc.OR.ja.EQ.kfqexc).AND.i*j.GT.0) THEN
342  nchn=nchn+1
343  isig(nchn,1)=i
344  isig(nchn,2)=j
345  isig(nchn,3)=1
346  IF(ja.EQ.kfqexc) isig(nchn,3)=2
347  IF(isig(nchn,isig(nchn,3)).GT.0) sigh(nchn)=facqsa*brpos
348  IF(isig(nchn,isig(nchn,3)).LT.0) sigh(nchn)=facqsa*brneg
349  ELSEIF(ia.EQ.kfqexc.AND.i.EQ.-j) THEN
350  nchn=nchn+1
351  isig(nchn,1)=i
352  isig(nchn,2)=j
353  isig(nchn,3)=1
354  IF(i.GT.0) sigh(nchn)=(8d0/3d0)*facqsb*brpos
355  IF(i.LT.0) sigh(nchn)=(8d0/3d0)*facqsb*brneg
356  nchn=nchn+1
357  isig(nchn,1)=i
358  isig(nchn,2)=j
359  isig(nchn,3)=2
360  IF(j.GT.0) sigh(nchn)=(8d0/3d0)*facqsb*brpos
361  IF(j.LT.0) sigh(nchn)=(8d0/3d0)*facqsb*brneg
362  ELSEIF(i.EQ.-j) THEN
363  nchn=nchn+1
364  isig(nchn,1)=i
365  isig(nchn,2)=j
366  isig(nchn,3)=1
367  IF(i.GT.0) sigh(nchn)=facqsb*brpos
368  IF(i.LT.0) sigh(nchn)=facqsb*brneg
369  nchn=nchn+1
370  isig(nchn,1)=i
371  isig(nchn,2)=j
372  isig(nchn,3)=2
373  IF(j.GT.0) sigh(nchn)=facqsb*brpos
374  IF(j.LT.0) sigh(nchn)=facqsb*brneg
375  ELSEIF(ia.EQ.kfqexc.OR.ja.EQ.kfqexc) THEN
376  nchn=nchn+1
377  isig(nchn,1)=i
378  isig(nchn,2)=j
379  isig(nchn,3)=1
380  IF(ja.EQ.kfqexc) isig(nchn,3)=2
381  IF(isig(nchn,isig(nchn,3)).GT.0) sigh(nchn)=facqsb*brpos
382  IF(isig(nchn,isig(nchn,3)).LT.0) sigh(nchn)=facqsb*brneg
383  ENDIF
384  250 CONTINUE
385  260 CONTINUE
386 
387  ELSEIF(isub.EQ.169) THEN
388 C...q + qbar -> e + e* (excited lepton)
389  kfqstr=kfpr(isub,2)
390  kcqstr=pycomp(kfqstr)
391  kfqexc=mod(kfqstr,kexcit)
392  facqsb=(comfac/12d0)*(sh/rtcm(41)**2)**2*(1d0-sqm4/sh)*
393  & (1d0+sqm4/sh)*(1d0+cth)*(1d0+((sh-sqm4)/(sh+sqm4))*cth)
394 C...Propagators: as simulated in PYOFSH and as desired
395  gmmq=pmas(kcqstr,1)*pmas(kcqstr,2)
396  hbw4=gmmq/((sqm4-pmas(kcqstr,1)**2)**2+gmmq**2)
397  CALL pywidt(kfqstr,sqm4,wdtp,wdte)
398  gmmqc=sqrt(sqm4)*wdtp(0)
399  hbw4c=gmmqc/((sqm4-pmas(kcqstr,1)**2)**2+gmmqc**2)
400  facqsb=facqsb*hbw4c/hbw4
401 C...Branching ratios.
402  brpos=(wdte(0,1)+wdte(0,2)+wdte(0,4))/wdtp(0)
403  brneg=(wdte(0,1)+wdte(0,3)+wdte(0,4))/wdtp(0)
404  DO 270 i=mmin1,mmax1
405  ia=iabs(i)
406  IF(i.EQ.0.OR.ia.GT.6.OR.kfac(1,i).EQ.0) goto 270
407  j=-i
408  ja=iabs(j)
409  IF(j.EQ.0.OR.ja.GT.6.OR.kfac(2,j).EQ.0) goto 270
410  nchn=nchn+1
411  isig(nchn,1)=i
412  isig(nchn,2)=j
413  isig(nchn,3)=1
414  IF(i.GT.0) sigh(nchn)=facqsb*brpos
415  IF(i.LT.0) sigh(nchn)=facqsb*brneg
416  nchn=nchn+1
417  isig(nchn,1)=i
418  isig(nchn,2)=j
419  isig(nchn,3)=2
420  IF(j.GT.0) sigh(nchn)=facqsb*brpos
421  IF(j.LT.0) sigh(nchn)=facqsb*brneg
422  270 CONTINUE
423  ENDIF
424 
425  ELSEIF(isub.LE.360) THEN
426  IF(isub.EQ.341.OR.isub.EQ.342) THEN
427 C...l + l -> H_L++/-- or H_R++/--.
428  kfres=kfpr(isub,1)
429  kfrec=pycomp(kfres)
430  CALL pywidt(kfres,sh,wdtp,wdte)
431  hs=shr*wdtp(0)
432  facbw=8d0*comfac/((sh-pmas(kfrec,1)**2)**2+hs**2)
433  DO 290 i=mmin1,mmax1
434  ia=iabs(i)
435  IF((ia.NE.11.AND.ia.NE.13.AND.ia.NE.15).OR.kfac(1,i).EQ.0)
436  & goto 290
437  DO 280 j=mmin2,mmax2
438  ja=iabs(j)
439  IF((ja.NE.11.AND.ja.NE.13.AND.ja.NE.15).OR.kfac(2,j).EQ.0)
440  & goto 280
441  IF(i*j.LT.0) goto 280
442  kchh=(kchg(ia,1)*isign(1,i)+kchg(ja,1)*isign(1,j))/3
443  nchn=nchn+1
444  isig(nchn,1)=i
445  isig(nchn,2)=j
446  isig(nchn,3)=1
447  hi=sh*parp(181+3*((ia-11)/2)+(ja-11)/2)**2/(8d0*paru(1))
448  hf=shr*(wdte(0,1)+wdte(0,(5-kchh/2)/2)+wdte(0,4))
449  sigh(nchn)=hi*facbw*hf
450  280 CONTINUE
451  290 CONTINUE
452 
453  ELSEIF(isub.GE.343.AND.isub.LE.348) THEN
454 C...l + gamma -> H_L++/-- l' or l + gamma -> H_R++/-- l'.
455  kfres=kfpr(isub,1)
456  kfrec=pycomp(kfres)
457 C...Propagators: as simulated in PYOFSH and as desired
458  hbw3=pmas(kfrec,1)*pmas(kfrec,2)/((sqm3-pmas(kfrec,1)**2)**2+
459  & (pmas(kfrec,1)*pmas(kfrec,2))**2)
460  CALL pywidt(kfres,sqm3,wdtp,wdte)
461  gmmc=sqrt(sqm3)*wdtp(0)
462  hbw3c=gmmc/((sqm3-pmas(kfrec,1)**2)**2+gmmc**2)
463  fhcc=comfac*aem*hbw3c/hbw3
464  DO 310 i=mmina,mmaxa
465  ia=iabs(i)
466  IF(ia.NE.11.AND.ia.NE.13.AND.ia.NE.15) goto 310
467  sqml=pmas(ia,1)**2
468  j=isign(kfpr(isub,2),-i)
469  kchh=isign(2,kchg(ia,1)*isign(1,i))
470  widsc=(wdte(0,1)+wdte(0,(5-kchh/2)/2)+wdte(0,4))/wdtp(0)
471  smm1=8d0*(sh+th-sqm3)*(sh+th-2d0*sqm3-sqml-sqm4)/
472  & (uh-sqm3)**2
473  smm2=2d0*((2d0*sqm3-3d0*sqml)*sqm4+(sqml-2d0*sqm4)*th-
474  & (th-sqm4)*sh)/(th-sqm4)**2
475  smm3=2d0*((2d0*sqm3-3d0*sqm4+th)*sqml-(2d0*sqml-sqm4+th)*
476  & sh)/(sh-sqml)**2
477  smm12=4d0*((2d0*sqml-sqm4-2d0*sqm3+th)*sh+(th-3d0*sqm3-
478  & 3d0*sqm4)*th+(2d0*sqm3-2d0*sqml+3d0*sqm4)*sqm3)/
479  & ((uh-sqm3)*(th-sqm4))
480  smm13=-4d0*((th+sqml-2d0*sqm4)*th-(sqm3+3d0*sqml-2d0*sqm4)*
481  & sqm3+(sqm3+3d0*sqml+th)*sh-(th-sqm3+sh)**2)/
482  & ((uh-sqm3)*(sh-sqml))
483  smm23=-4d0*((sqml-sqm4+sqm3)*th-sqm3**2+sqm3*(sqml+sqm4)-
484  & 3d0*sqml*sqm4-(sqml-sqm4-sqm3+th)*sh)/
485  & ((sh-sqml)*(th-sqm4))
486  smm=(sh/(sh-sqml))**2*(smm1+smm2+smm3+smm12+smm13+smm23)*
487  & parp(181+3*((ia-11)/2)+(iabs(j)-11)/2)**2/(4d0*paru(1))
488  DO 300 isde=1,2
489  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,22).EQ.0) goto 300
490  IF(isde.EQ.2.AND.kfac(1,22)*kfac(2,i).EQ.0) goto 300
491  nchn=nchn+1
492  isig(nchn,isde)=i
493  isig(nchn,3-isde)=22
494  isig(nchn,3)=0
495  sigh(nchn)=fhcc*smm*widsc
496  300 CONTINUE
497  310 CONTINUE
498 
499  ELSEIF(isub.EQ.349.OR.isub.EQ.350) THEN
500 C...f + fbar -> H_L++ + H_L-- or H_R++ + H_R--
501  kfres=kfpr(isub,1)
502  kfrec=pycomp(kfres)
503  sqmh=pmas(kfrec,1)**2
504  gmmh=pmas(kfrec,1)*pmas(kfrec,2)
505 C...Propagators: H++/-- as simulated in PYOFSH and as desired
506  hbw3=gmmh/((sqm3-sqmh)**2+gmmh**2)
507  CALL pywidt(kfres,sqm3,wdtp,wdte)
508  gmmh3=sqrt(sqm3)*wdtp(0)
509  hbw3c=gmmh3/((sqm3-sqmh)**2+gmmh3**2)
510  hbw4=gmmh/((sqm4-sqmh)**2+gmmh**2)
511  CALL pywidt(kfres,sqm4,wdtp,wdte)
512  gmmh4=sqrt(sqm4)*wdtp(0)
513  hbw4c=gmmh4/((sqm4-sqmh)**2+gmmh4**2)
514 C...Kinematical and coupling functions
515  fachh=comfac*(hbw3c/hbw3)*(hbw4c/hbw4)*(th*uh-sqm3*sqm4)
516  xwhh=(1d0-2d0*xwv)/(8d0*xwv*(1d0-xwv))
517 C...Loop over allowed flavours
518  DO 320 i=mmina,mmaxa
519  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 320
520  ei=kchg(iabs(i),1)/3d0
521  ai=sign(1d0,ei+0.1d0)
522  vi=ai-4d0*ei*xwv
523  fcoi=1d0
524  IF(iabs(i).LE.10) fcoi=faca/3d0
525  IF(isub.EQ.349) THEN
526  hbwz=1d0/((sh-sqmz)**2+gmmz**2)
527  IF(iabs(i).LT.10) THEN
528  dsighh=8d0*aem**2*(ei**2/sh2+
529  & 2d0*ei*vi*xwhh*(sh-sqmz)*hbwz/sh+
530  & (vi**2+ai**2)*xwhh**2*hbwz)
531  ELSE
532  iaoff=181+3*((iabs(i)-11)/2)
533  hsum=(parp(iaoff)**2+parp(iaoff+1)**2+parp(iaoff+2)**2)/
534  & (4d0*paru(1))
535  dsighh=8d0*aem**2*(ei**2/sh2+
536  & 2d0*ei*vi*xwhh*(sh-sqmz)*hbwz/sh+
537  & (vi**2+ai**2)*xwhh**2*hbwz)+
538  & 8d0*aem*(ei*hsum/(sh*th)+
539  & (vi+ai)*xwhh*hsum*(sh-sqmz)*hbwz/th)+
540  & 4d0*hsum**2/th2
541  ENDIF
542  ELSE
543  IF(iabs(i).LT.10) THEN
544  dsighh=8d0*aem**2*ei**2/sh2
545  ELSE
546  iaoff=181+3*((iabs(i)-11)/2)
547  hsum=(parp(iaoff)**2+parp(iaoff+1)**2+parp(iaoff+2)**2)/
548  & (4d0*paru(1))
549  dsighh=8d0*aem**2*ei**2/sh2+8d0*aem*ei*hsum/(sh*th)+
550  & 4d0*hsum**2/th2
551  ENDIF
552  ENDIF
553  nchn=nchn+1
554  isig(nchn,1)=i
555  isig(nchn,2)=-i
556  isig(nchn,3)=1
557  sigh(nchn)=fachh*fcoi*dsighh
558  320 CONTINUE
559 
560  ELSEIF(isub.EQ.351.OR.isub.EQ.352) THEN
561 C...f + f' -> f" + f"' + H++/-- (W+/- + W+/- -> H++/-- as inner process)
562  kfres=kfpr(isub,1)
563  kfrec=pycomp(kfres)
564  sqmh=pmas(kfrec,1)**2
565  IF(isub.EQ.351) facnor=parp(190)**8*parp(192)**2
566  IF(isub.EQ.352) facnor=parp(191)**6*2d0*
567  & pmas(pycomp(9900024),1)**2
568  facww=comfac*facnor*taup*vint(2)*vint(219)
569  facprt=1d0/((vint(204)**2-vint(215))*
570  & (vint(209)**2-vint(216)))
571  facpru=1d0/((vint(204)**2+2d0*vint(217))*
572  & (vint(209)**2+2d0*vint(218)))
573  CALL pywidt(kfres,sh,wdtp,wdte)
574  hs=shr*wdtp(0)
575  facbw=(1d0/paru(1))*vint(2)/((sh-sqmh)**2+hs**2)
576  IF(abs(shr-pmas(kfrec,1)).GT.parp(48)*pmas(kfrec,2))
577  & facbw=0d0
578  DO 340 i=mmin1,mmax1
579  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 340
580  IF(isub.EQ.352.AND.iabs(i).GT.10) goto 340
581  kchwi=(1-2*mod(iabs(i),2))*isign(1,i)
582  DO 330 j=mmin2,mmax2
583  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 330
584  IF(isub.EQ.352.AND.iabs(j).GT.10) goto 330
585  kchwj=(1-2*mod(iabs(j),2))*isign(1,j)
586  kchh=kchwi+kchwj
587  IF(iabs(kchh).NE.2) goto 330
588  faclr=vint(180+i)*vint(180+j)
589  hf=shr*(wdte(0,1)+wdte(0,(5-kchh/2)/2)+wdte(0,4))
590  IF(i.EQ.j.AND.iabs(i).GT.10) THEN
591  facprp=0.5d0*(facprt+facpru)**2
592  ELSE
593  facprp=facprt**2
594  ENDIF
595  nchn=nchn+1
596  isig(nchn,1)=i
597  isig(nchn,2)=j
598  isig(nchn,3)=1
599  sigh(nchn)=faclr*facww*facprp*facbw*hf
600  330 CONTINUE
601  340 CONTINUE
602 
603  ELSEIF(isub.EQ.353) THEN
604 C...f + fbar -> Z_R0
605  sqmzr=pmas(pycomp(kfpr(isub,1)),1)**2
606  CALL pywidt(kfpr(isub,1),sh,wdtp,wdte)
607  hs=shr*wdtp(0)
608  facbw=4d0*comfac/((sh-sqmzr)**2+hs**2)*3d0
609  hf=shr*(wdte(0,1)+wdte(0,2)+wdte(0,4))
610  hp=(aem/(3d0*(1d0-2d0*xw)))*xwc*sh
611  DO 350 i=mmina,mmaxa
612  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 350
613  IF(iabs(i).LE.8) THEN
614  ei=kchg(iabs(i),1)/3d0
615  ai=sign(1d0,ei+0.1d0)*(1d0-2d0*xw)
616  vi=sign(1d0,ei+0.1d0)-4d0*ei*xw
617  ELSE
618  ai=-(1d0-2d0*xw)
619  vi=-1d0+4d0*xw
620  ENDIF
621  hi=hp*(vi**2+ai**2)
622  IF(iabs(i).LE.10) hi=hi*faca/3d0
623  nchn=nchn+1
624  isig(nchn,1)=i
625  isig(nchn,2)=-i
626  isig(nchn,3)=1
627  sigh(nchn)=hi*facbw*hf
628  350 CONTINUE
629 
630  ELSEIF(isub.EQ.354) THEN
631 C...f + fbar' -> W_R+/-
632  sqmwr=pmas(pycomp(kfpr(isub,1)),1)**2
633  CALL pywidt(kfpr(isub,1),sh,wdtp,wdte)
634  hs=shr*wdtp(0)
635  facbw=4d0*comfac/((sh-sqmwr)**2+hs**2)*3d0
636  hp=aem/(24d0*xw)*sh
637  DO 370 i=mmin1,mmax1
638  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 370
639  ia=iabs(i)
640  DO 360 j=mmin2,mmax2
641  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 360
642  ja=iabs(j)
643  IF(i*j.GT.0.OR.mod(ia+ja,2).EQ.0) goto 360
644  IF((ia.LE.10.AND.ja.GT.10).OR.(ia.GT.10.AND.ja.LE.10))
645  & goto 360
646  kchw=(kchg(ia,1)*isign(1,i)+kchg(ja,1)*isign(1,j))/3
647  hi=hp*2d0
648  IF(ia.LE.10) hi=hi*vckm((ia+1)/2,(ja+1)/2)*faca/3d0
649  nchn=nchn+1
650  isig(nchn,1)=i
651  isig(nchn,2)=j
652  isig(nchn,3)=1
653  hf=shr*(wdte(0,1)+wdte(0,(5-kchw)/2)+wdte(0,4))
654  sigh(nchn)=hi*facbw*hf
655  360 CONTINUE
656  370 CONTINUE
657  ENDIF
658 
659  ELSEIF(isub.LE.400) THEN
660  IF(isub.EQ.391) THEN
661 C...f + fbar -> G*.
662  kfgstr=kfpr(isub,1)
663  kcgstr=pycomp(kfgstr)
664  CALL pywidt(kfgstr,sh,wdtp,wdte)
665  hs=shr*wdtp(0)
666  hf=shr*(wdte(0,1)+wdte(0,2)+wdte(0,4))
667  facg=comfac*parp(50)**2/(16d0*paru(1))*sh*hf/
668  & ((sh-pmas(kcgstr,1)**2)**2+hs**2)
669 C...Modify cross section in wings of peak.
670  facg = facg * sh**2 / pmas(kcgstr,1)**4
671  DO 380 i=mmina,mmaxa
672  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 380
673  hi=1d0
674  IF(iabs(i).LE.10) hi=hi*faca/3d0
675  nchn=nchn+1
676  isig(nchn,1)=i
677  isig(nchn,2)=-i
678  isig(nchn,3)=1
679  sigh(nchn)=facg*hi
680  380 CONTINUE
681 
682  ELSEIF(isub.EQ.392) THEN
683 C...g + g -> G*.
684  kfgstr=kfpr(isub,1)
685  kcgstr=pycomp(kfgstr)
686  CALL pywidt(kfgstr,sh,wdtp,wdte)
687  hs=shr*wdtp(0)
688  hf=shr*(wdte(0,1)+wdte(0,2)+wdte(0,4))
689  facg=comfac*parp(50)**2/(32d0*paru(1))*sh*hf/
690  & ((sh-pmas(kcgstr,1)**2)**2+hs**2)
691 C...Modify cross section in wings of peak.
692  facg = facg * sh**2 / pmas(kcgstr,1)**4
693  IF(kfac(1,21)*kfac(2,21).EQ.0) goto 390
694  nchn=nchn+1
695  isig(nchn,1)=21
696  isig(nchn,2)=21
697  isig(nchn,3)=1
698  sigh(nchn)=facg
699  390 CONTINUE
700 
701  ELSEIF(isub.EQ.393) THEN
702 C...q + qbar -> g + G*.
703  kfgstr=kfpr(isub,2)
704  kcgstr=pycomp(kfgstr)
705  facg=comfac*parp(50)**2*as*sh/(72d0*paru(1)*sqm4)*
706  & (4d0*(th2+uh2)/sh2+9d0*(th+uh)/sh+(th2/uh+uh2/th)/sh+
707  & 3d0*(4d0+th/uh+uh/th)+4d0*(sh/uh+sh/th)+
708  & 2d0*sh2/(th*uh))
709 C...Propagators: as simulated in PYOFSH and as desired
710  gmmg=pmas(kcgstr,1)*pmas(kcgstr,2)
711  hbw4=gmmg/((sqm4-pmas(kcgstr,1)**2)**2+gmmg**2)
712  CALL pywidt(kfgstr,sqm4,wdtp,wdte)
713  hs=sqrt(sqm4)*wdtp(0)
714  hf=sqrt(sqm4)*(wdte(0,1)+wdte(0,2)+wdte(0,4))
715  hbw4c=hf/((sqm4-pmas(kcgstr,1)**2)**2+hs**2)
716  facg=facg*hbw4c/hbw4
717  DO 400 i=mmina,mmaxa
718  IF(i.EQ.0.OR.iabs(i).GT.mstp(58).OR.
719  & kfac(1,i)*kfac(2,-i).EQ.0) goto 400
720  nchn=nchn+1
721  isig(nchn,1)=i
722  isig(nchn,2)=-i
723  isig(nchn,3)=1
724  sigh(nchn)=facg
725  400 CONTINUE
726 
727  ELSEIF(isub.EQ.394) THEN
728 C...q + g -> q + G*.
729  kfgstr=kfpr(isub,2)
730  kcgstr=pycomp(kfgstr)
731  facg=-comfac*parp(50)**2*as*sh/(192d0*paru(1)*sqm4)*
732  & (4d0*(sh2+uh2)/(th*sh)+9d0*(sh+uh)/sh+sh/uh+uh2/sh2+
733  & 3d0*th*(4d0+sh/uh+uh/sh)/sh+4d0*th2*(1d0/uh+1d0/sh)/sh+
734  & 2d0*th2*th/(uh*sh2))
735 C...Propagators: as simulated in PYOFSH and as desired
736  gmmg=pmas(kcgstr,1)*pmas(kcgstr,2)
737  hbw4=gmmg/((sqm4-pmas(kcgstr,1)**2)**2+gmmg**2)
738  CALL pywidt(kfgstr,sqm4,wdtp,wdte)
739  hs=sqrt(sqm4)*wdtp(0)
740  hf=sqrt(sqm4)*(wdte(0,1)+wdte(0,2)+wdte(0,4))
741  hbw4c=hf/((sqm4-pmas(kcgstr,1)**2)**2+hs**2)
742  facg=facg*hbw4c/hbw4
743  DO 420 i=mmina,mmaxa
744  IF(i.EQ.0.OR.iabs(i).GT.mstp(58)) goto 420
745  DO 410 isde=1,2
746  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 410
747  IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 410
748  nchn=nchn+1
749  isig(nchn,isde)=i
750  isig(nchn,3-isde)=21
751  isig(nchn,3)=1
752  sigh(nchn)=facg
753  410 CONTINUE
754  420 CONTINUE
755 
756  ELSEIF(isub.EQ.395) THEN
757 C...g + g -> g + G*.
758  kfgstr=kfpr(isub,2)
759  kcgstr=pycomp(kfgstr)
760  facg=comfac*3d0*parp(50)**2*as*sh/(32d0*paru(1)*sqm4)*
761  & ((th2+th*uh+uh2)**2/(sh2*th*uh)+2d0*(th2/uh+uh2/th)/sh+
762  & 3d0*(th/uh+uh/th)+2d0*(sh/uh+sh/th)+sh2/(th*uh))
763 C...Propagators: as simulated in PYOFSH and as desired
764  gmmg=pmas(kcgstr,1)*pmas(kcgstr,2)
765  hbw4=gmmg/((sqm4-pmas(kcgstr,1)**2)**2+gmmg**2)
766  CALL pywidt(kfgstr,sqm4,wdtp,wdte)
767  hs=sqrt(sqm4)*wdtp(0)
768  hf=sqrt(sqm4)*(wdte(0,1)+wdte(0,2)+wdte(0,4))
769  hbw4c=hf/((sqm4-pmas(kcgstr,1)**2)**2+hs**2)
770  facg=facg*hbw4c/hbw4
771  IF(kfac(1,21)*kfac(2,21).NE.0) THEN
772  nchn=nchn+1
773  isig(nchn,1)=21
774  isig(nchn,2)=21
775  isig(nchn,3)=1
776  sigh(nchn)=facg
777  ENDIF
778  ENDIF
779  ENDIF
780 
781  RETURN
782  END