Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pysgwz.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pysgwz.f
1 
2 C*********************************************************************
3 
4 C...PYSGWZ
5 C...Subprocess cross sections for W/Z processes,
6 C...except that longitudinal WW scattering is in Higgs sector.
7 C...Auxiliary to PYSIGH.
8 
9  SUBROUTINE pysgwz(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/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
23  common/pypars/mstp(200),parp(200),msti(200),pari(200)
24  common/pyint1/mint(400),vint(400)
25  common/pyint2/iset(500),kfpr(500,2),coef(500,20),icol(40,4,2)
26  common/pyint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)
27  common/pyint4/mwid(500),wids(500,5)
28  common/pytcsm/itcm(0:99),rtcm(0:99)
29  common/pysgcm/isub,isubsv,mmin1,mmax1,mmin2,mmax2,mmina,mmaxa,
30  &kfac(2,-40:40),comfac,fack,faca,sh,th,uh,sh2,th2,uh2,sqm3,sqm4,
31  &shr,sqpth,taup,be34,cth,x(2),sqmz,sqmw,gmmz,gmmw,
32  &aem,as,xw,xw1,xwc,xwv,poll,polr,polll,polrr
33  SAVE /pydat1/,/pydat2/,/pydat3/,/pysubs/,/pypars/,/pyint1/,
34  &/pyint2/,/pyint3/,/pyint4/,/pytcsm/,/pysgcm/
35 C...Local arrays and complex numbers
36  dimension wdtp(0:400),wdte(0:400,0:5),hgz(6,3),hl3(3),hr3(3),
37  &hl4(3),hr4(3)
38  COMPLEX*16 coulck,coulcp,coulcd,coulcr,coulcs
39 
40 C...Differential cross section expressions.
41 
42  IF(isub.LE.20) THEN
43  IF(isub.EQ.1) THEN
44 C...f + fbar -> gamma*/Z0
45  mint(61)=2
46  CALL pywidt(23,sh,wdtp,wdte)
47  hs=shr*wdtp(0)
48  facz=4d0*comfac*3d0
49  hp0=aem/3d0*sh
50  hp1=aem/3d0*xwc*sh
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  hi0=hp0
57  IF(iabs(i).LE.10) hi0=hi0*faca/3d0
58  hi1=hp1
59  IF(iabs(i).LE.10) hi1=hi1*faca/3d0
60  nchn=nchn+1
61  isig(nchn,1)=i
62  isig(nchn,2)=-i
63  isig(nchn,3)=1
64  sigh(nchn)=facz*(ei**2/sh2*hi0*hp0*vint(111)+
65  & ei*vi*(1d0-sqmz/sh)/((sh-sqmz)**2+hs**2)*
66  & (hi0*hp1+hi1*hp0)*vint(112)+(vi**2+ai**2)/
67  & ((sh-sqmz)**2+hs**2)*hi1*hp1*vint(114))
68  100 CONTINUE
69 
70  ELSEIF(isub.EQ.2) THEN
71 C...f + fbar' -> W+/-
72  CALL pywidt(24,sh,wdtp,wdte)
73  hs=shr*wdtp(0)
74  facbw=4d0*comfac/((sh-sqmw)**2+hs**2)*3d0
75  hp=aem/(24d0*xw)*sh
76  DO 120 i=mmin1,mmax1
77  IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 120
78  ia=iabs(i)
79  DO 110 j=mmin2,mmax2
80  IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 110
81  ja=iabs(j)
82  IF(i*j.GT.0.OR.mod(ia+ja,2).EQ.0) goto 110
83  IF((ia.LE.10.AND.ja.GT.10).OR.(ia.GT.10.AND.ja.LE.10))
84  & goto 110
85  kchw=(kchg(ia,1)*isign(1,i)+kchg(ja,1)*isign(1,j))/3
86  hi=hp*2d0
87  IF(ia.LE.10) hi=hi*vckm((ia+1)/2,(ja+1)/2)*faca/3d0
88  nchn=nchn+1
89  isig(nchn,1)=i
90  isig(nchn,2)=j
91  isig(nchn,3)=1
92  hf=shr*(wdte(0,1)+wdte(0,(5-kchw)/2)+wdte(0,4))
93  sigh(nchn)=hi*facbw*hf
94  110 CONTINUE
95  120 CONTINUE
96 
97  ELSEIF(isub.EQ.15) THEN
98 C...f + fbar -> g + (gamma*/Z0) (q + qbar -> g + (gamma*/Z0) only)
99  faczg=comfac*as*aem*(8d0/9d0)*(th2+uh2+2d0*sqm4*sh)/(th*uh)
100 C...gamma, gamma/Z interference and Z couplings to final fermion pairs
101  hfgg=0d0
102  hfgz=0d0
103  hfzz=0d0
104  radc4=1d0+pyalps(sqm4)/paru(1)
105  DO 130 i=1,min(16,mdcy(23,3))
106  idc=i+mdcy(23,2)-1
107  IF(mdme(idc,1).LT.0) goto 130
108  imdm=0
109  IF(mdme(idc,1).EQ.1.OR.mdme(idc,1).EQ.2.OR.mdme(idc,1).EQ.4)
110  & imdm=1
111  IF(i.LE.8) THEN
112  ef=kchg(i,1)/3d0
113  af=sign(1d0,ef+0.1d0)
114  vf=af-4d0*ef*xwv
115  ELSEIF(i.LE.16) THEN
116  ef=kchg(i+2,1)/3d0
117  af=sign(1d0,ef+0.1d0)
118  vf=af-4d0*ef*xwv
119  ENDIF
120  rm1=pmas(iabs(kfdp(idc,1)),1)**2/sqm4
121  IF(4d0*rm1.LT.1d0) THEN
122  fcof=1d0
123  IF(i.LE.8) fcof=3d0*radc4
124  be34=sqrt(max(0d0,1d0-4d0*rm1))
125  IF(imdm.EQ.1) THEN
126  hfgg=hfgg+fcof*ef**2*(1d0+2d0*rm1)*be34
127  hfgz=hfgz+fcof*ef*vf*(1d0+2d0*rm1)*be34
128  hfzz=hfzz+fcof*(vf**2*(1d0+2d0*rm1)+
129  & af**2*(1d0-4d0*rm1))*be34
130  ENDIF
131  ENDIF
132  130 CONTINUE
133 C...Propagators: as simulated in PYOFSH and as desired
134  hbw4=(1d0/paru(1))*gmmz/((sqm4-sqmz)**2+gmmz**2)
135  mint15=mint(15)
136  mint(15)=1
137  mint(61)=1
138  CALL pywidt(23,sqm4,wdtp,wdte)
139  mint(15)=mint15
140  hfaem=(paru(108)/paru(2))*(2d0/3d0)
141  hfgg=hfgg*hfaem*vint(111)/sqm4
142  hfgz=hfgz*hfaem*vint(112)/sqm4
143  hfzz=hfzz*hfaem*vint(114)/sqm4
144 C...Loop over flavours; consider full gamma/Z structure
145  DO 140 i=mmina,mmaxa
146  IF(i.EQ.0.OR.iabs(i).GT.mstp(58).OR.
147  & kfac(1,i)*kfac(2,-i).EQ.0) goto 140
148  ei=kchg(iabs(i),1)/3d0
149  ai=sign(1d0,ei)
150  vi=ai-4d0*ei*xwv
151  nchn=nchn+1
152  isig(nchn,1)=i
153  isig(nchn,2)=-i
154  isig(nchn,3)=1
155  sigh(nchn)=faczg*(ei**2*hfgg+ei*vi*hfgz+
156  & (vi**2+ai**2)*hfzz)/hbw4
157  140 CONTINUE
158 
159  ELSEIF(isub.EQ.16) THEN
160 C...f + fbar' -> g + W+/- (q + qbar' -> g + W+/- only)
161  facwg=comfac*as*aem/xw*2d0/9d0*(th2+uh2+2d0*sqm4*sh)/(th*uh)
162 C...Propagators: as simulated in PYOFSH and as desired
163  hbw4=gmmw/((sqm4-sqmw)**2+gmmw**2)
164  CALL pywidt(24,sqm4,wdtp,wdte)
165  gmmwc=sqrt(sqm4)*wdtp(0)
166  hbw4c=gmmwc/((sqm4-sqmw)**2+gmmwc**2)
167  facwg=facwg*hbw4c/hbw4
168  DO 160 i=mmin1,mmax1
169  ia=iabs(i)
170  IF(i.EQ.0.OR.ia.GT.10.OR.kfac(1,i).EQ.0) goto 160
171  DO 150 j=mmin2,mmax2
172  ja=iabs(j)
173  IF(j.EQ.0.OR.ja.GT.10.OR.kfac(2,j).EQ.0) goto 150
174  IF(i*j.GT.0.OR.mod(ia+ja,2).EQ.0) goto 150
175  kchw=(kchg(ia,1)*isign(1,i)+kchg(ja,1)*isign(1,j))/3
176  widsc=(wdte(0,1)+wdte(0,(5-kchw)/2)+wdte(0,4))/wdtp(0)
177  fckm=vckm((ia+1)/2,(ja+1)/2)
178  nchn=nchn+1
179  isig(nchn,1)=i
180  isig(nchn,2)=j
181  isig(nchn,3)=1
182  sigh(nchn)=facwg*fckm*widsc
183  150 CONTINUE
184  160 CONTINUE
185 
186  ELSEIF(isub.EQ.19) THEN
187 C...f + fbar -> gamma + (gamma*/Z0)
188  facgz=comfac*2d0*aem**2*(th2+uh2+2d0*sqm4*sh)/(th*uh)
189 C...gamma, gamma/Z interference and Z couplings to final fermion pairs
190  hfgg=0d0
191  hfgz=0d0
192  hfzz=0d0
193  radc4=1d0+pyalps(sqm4)/paru(1)
194  DO 170 i=1,min(16,mdcy(23,3))
195  idc=i+mdcy(23,2)-1
196  IF(mdme(idc,1).LT.0) goto 170
197  imdm=0
198  IF(mdme(idc,1).EQ.1.OR.mdme(idc,1).EQ.2.OR.mdme(idc,1).EQ.4)
199  & imdm=1
200  IF(i.LE.8) THEN
201  ef=kchg(i,1)/3d0
202  af=sign(1d0,ef+0.1d0)
203  vf=af-4d0*ef*xwv
204  ELSEIF(i.LE.16) THEN
205  ef=kchg(i+2,1)/3d0
206  af=sign(1d0,ef+0.1d0)
207  vf=af-4d0*ef*xwv
208  ENDIF
209  rm1=pmas(iabs(kfdp(idc,1)),1)**2/sqm4
210  IF(4d0*rm1.LT.1d0) THEN
211  fcof=1d0
212  IF(i.LE.8) fcof=3d0*radc4
213  be34=sqrt(max(0d0,1d0-4d0*rm1))
214  IF(imdm.EQ.1) THEN
215  hfgg=hfgg+fcof*ef**2*(1d0+2d0*rm1)*be34
216  hfgz=hfgz+fcof*ef*vf*(1d0+2d0*rm1)*be34
217  hfzz=hfzz+fcof*(vf**2*(1d0+2d0*rm1)+
218  & af**2*(1d0-4d0*rm1))*be34
219  ENDIF
220  ENDIF
221  170 CONTINUE
222 C...Propagators: as simulated in PYOFSH and as desired
223  hbw4=(1d0/paru(1))*gmmz/((sqm4-sqmz)**2+gmmz**2)
224  mint15=mint(15)
225  mint(15)=1
226  mint(61)=1
227  CALL pywidt(23,sqm4,wdtp,wdte)
228  mint(15)=mint15
229  hfaem=(paru(108)/paru(2))*(2d0/3d0)
230  hfgg=hfgg*hfaem*vint(111)/sqm4
231  hfgz=hfgz*hfaem*vint(112)/sqm4
232  hfzz=hfzz*hfaem*vint(114)/sqm4
233 C...Loop over flavours; consider full gamma/Z structure
234  DO 180 i=mmina,mmaxa
235  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 180
236  ei=kchg(iabs(i),1)/3d0
237  ai=sign(1d0,ei)
238  vi=ai-4d0*ei*xwv
239  fcoi=1d0
240  IF(iabs(i).LE.10) fcoi=faca/3d0
241  nchn=nchn+1
242  isig(nchn,1)=i
243  isig(nchn,2)=-i
244  isig(nchn,3)=1
245  sigh(nchn)=facgz*fcoi*ei**2*(ei**2*hfgg+ei*vi*hfgz+
246  & (vi**2+ai**2)*hfzz)/hbw4
247  180 CONTINUE
248 
249  ELSEIF(isub.EQ.20) THEN
250 C...f + fbar' -> gamma + W+/-
251  facgw=comfac*0.5d0*aem**2/xw
252 C...Propagators: as simulated in PYOFSH and as desired
253  hbw4=gmmw/((sqm4-sqmw)**2+gmmw**2)
254  CALL pywidt(24,sqm4,wdtp,wdte)
255  gmmwc=sqrt(sqm4)*wdtp(0)
256  hbw4c=gmmwc/((sqm4-sqmw)**2+gmmwc**2)
257  facgw=facgw*hbw4c/hbw4
258 C...Anomalous couplings
259  term1=(th2+uh2+2d0*sqm4*sh)/(th*uh)
260  term2=0d0
261  term3=0d0
262  IF(itcm(5).GE.1.AND.itcm(5).LE.4) THEN
263  term2=rtcm(46)*(th-uh)/(th+uh)
264  term3=0.5d0*rtcm(46)**2*(th*uh+(th2+uh2)*sh/
265  & (4d0*sqmw))/(th+uh)**2
266  ENDIF
267  DO 200 i=mmin1,mmax1
268  ia=iabs(i)
269  IF(i.EQ.0.OR.ia.GT.20.OR.kfac(1,i).EQ.0) goto 200
270  DO 190 j=mmin2,mmax2
271  ja=iabs(j)
272  IF(j.EQ.0.OR.ja.GT.20.OR.kfac(2,j).EQ.0) goto 190
273  IF(i*j.GT.0.OR.mod(ia+ja,2).EQ.0) goto 190
274  IF((ia.LE.10.AND.ja.GT.10).OR.(ia.GT.10.AND.ja.LE.10))
275  & goto 190
276  kchw=(kchg(ia,1)*isign(1,i)+kchg(ja,1)*isign(1,j))/3
277  widsc=(wdte(0,1)+wdte(0,(5-kchw)/2)+wdte(0,4))/wdtp(0)
278  IF(ia.LE.10) THEN
279  facwr=uh/(th+uh)-1d0/3d0
280  fckm=vckm((ia+1)/2,(ja+1)/2)
281  fcoi=faca/3d0
282  ELSE
283  facwr=-th/(th+uh)
284  fckm=1d0
285  fcoi=1d0
286  ENDIF
287  facwk=term1*facwr**2+term2*facwr+term3
288  nchn=nchn+1
289  isig(nchn,1)=i
290  isig(nchn,2)=j
291  isig(nchn,3)=1
292  sigh(nchn)=facgw*facwk*fcoi*fckm*widsc
293  190 CONTINUE
294  200 CONTINUE
295  ENDIF
296 
297  ELSEIF(isub.LE.40) THEN
298  IF(isub.EQ.22) THEN
299 C...f + fbar -> (gamma*/Z0) + (gamma*/Z0)
300 C...Kinematics dependence
301  faczz=comfac*aem**2*((th2+uh2+2d0*(sqm3+sqm4)*sh)/(th*uh)-
302  & sqm3*sqm4*(1d0/th2+1d0/uh2))
303 C...gamma, gamma/Z interference and Z couplings to final fermion pairs
304  DO 220 i=1,6
305  DO 210 j=1,3
306  hgz(i,j)=0d0
307  210 CONTINUE
308  220 CONTINUE
309  radc3=1d0+pyalps(sqm3)/paru(1)
310  radc4=1d0+pyalps(sqm4)/paru(1)
311  DO 230 i=1,min(16,mdcy(23,3))
312  idc=i+mdcy(23,2)-1
313  IF(mdme(idc,1).LT.0) goto 230
314  imdm=0
315  IF(mdme(idc,1).EQ.1.OR.mdme(idc,1).EQ.2) imdm=1
316  IF(mdme(idc,1).EQ.4.OR.mdme(idc,1).EQ.5) imdm=mdme(idc,1)-2
317  IF(i.LE.8) THEN
318  ef=kchg(i,1)/3d0
319  af=sign(1d0,ef+0.1d0)
320  vf=af-4d0*ef*xwv
321  ELSEIF(i.LE.16) THEN
322  ef=kchg(i+2,1)/3d0
323  af=sign(1d0,ef+0.1d0)
324  vf=af-4d0*ef*xwv
325  ENDIF
326  rm1=pmas(iabs(kfdp(idc,1)),1)**2/sqm3
327  IF(4d0*rm1.LT.1d0) THEN
328  fcof=1d0
329  IF(i.LE.8) fcof=3d0*radc3
330  be34=sqrt(max(0d0,1d0-4d0*rm1))
331  IF(imdm.GE.1) THEN
332  hgz(1,imdm)=hgz(1,imdm)+fcof*ef**2*(1d0+2d0*rm1)*be34
333  hgz(2,imdm)=hgz(2,imdm)+fcof*ef*vf*(1d0+2d0*rm1)*be34
334  hgz(3,imdm)=hgz(3,imdm)+fcof*(vf**2*(1d0+2d0*rm1)+
335  & af**2*(1d0-4d0*rm1))*be34
336  ENDIF
337  ENDIF
338  rm1=pmas(iabs(kfdp(idc,1)),1)**2/sqm4
339  IF(4d0*rm1.LT.1d0) THEN
340  fcof=1d0
341  IF(i.LE.8) fcof=3d0*radc4
342  be34=sqrt(max(0d0,1d0-4d0*rm1))
343  IF(imdm.GE.1) THEN
344  hgz(4,imdm)=hgz(4,imdm)+fcof*ef**2*(1d0+2d0*rm1)*be34
345  hgz(5,imdm)=hgz(5,imdm)+fcof*ef*vf*(1d0+2d0*rm1)*be34
346  hgz(6,imdm)=hgz(6,imdm)+fcof*(vf**2*(1d0+2d0*rm1)+
347  & af**2*(1d0-4d0*rm1))*be34
348  ENDIF
349  ENDIF
350  230 CONTINUE
351 C...Propagators: as simulated in PYOFSH and as desired
352  hbw3=(1d0/paru(1))*gmmz/((sqm3-sqmz)**2+gmmz**2)
353  hbw4=(1d0/paru(1))*gmmz/((sqm4-sqmz)**2+gmmz**2)
354  mint15=mint(15)
355  mint(15)=1
356  mint(61)=1
357  CALL pywidt(23,sqm3,wdtp,wdte)
358  mint(15)=mint15
359  hfaem=(paru(108)/paru(2))*(2d0/3d0)
360  DO 240 j=1,3
361  hgz(1,j)=hgz(1,j)*hfaem*vint(111)/sqm3
362  hgz(2,j)=hgz(2,j)*hfaem*vint(112)/sqm3
363  hgz(3,j)=hgz(3,j)*hfaem*vint(114)/sqm3
364  240 CONTINUE
365  mint15=mint(15)
366  mint(15)=1
367  mint(61)=1
368  CALL pywidt(23,sqm4,wdtp,wdte)
369  mint(15)=mint15
370  hfaem=(paru(108)/paru(2))*(2d0/3d0)
371  DO 250 j=1,3
372  hgz(4,j)=hgz(4,j)*hfaem*vint(111)/sqm4
373  hgz(5,j)=hgz(5,j)*hfaem*vint(112)/sqm4
374  hgz(6,j)=hgz(6,j)*hfaem*vint(114)/sqm4
375  250 CONTINUE
376 C...Loop over flavours; separate left- and right-handed couplings
377  DO 270 i=mmina,mmaxa
378  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 270
379  ei=kchg(iabs(i),1)/3d0
380  ai=sign(1d0,ei)
381  vi=ai-4d0*ei*xwv
382  vali=vi-ai
383  vari=vi+ai
384  fcoi=1d0
385  IF(iabs(i).LE.10) fcoi=faca/3d0
386  DO 260 j=1,3
387  hl3(j)=ei**2*hgz(1,j)+ei*vali*hgz(2,j)+vali**2*hgz(3,j)
388  hr3(j)=ei**2*hgz(1,j)+ei*vari*hgz(2,j)+vari**2*hgz(3,j)
389  hl4(j)=ei**2*hgz(4,j)+ei*vali*hgz(5,j)+vali**2*hgz(6,j)
390  hr4(j)=ei**2*hgz(4,j)+ei*vari*hgz(5,j)+vari**2*hgz(6,j)
391  260 CONTINUE
392  faclr=hl3(1)*hl4(1)+hl3(1)*(hl4(2)+hl4(3))+
393  & hl4(1)*(hl3(2)+hl3(3))+hl3(2)*hl4(3)+hl4(2)*hl3(3)+
394  & hr3(1)*hr4(1)+hr3(1)*(hr4(2)+hr4(3))+
395  & hr4(1)*(hr3(2)+hr3(3))+hr3(2)*hr4(3)+hr4(2)*hr3(3)
396  nchn=nchn+1
397  isig(nchn,1)=i
398  isig(nchn,2)=-i
399  isig(nchn,3)=1
400  sigh(nchn)=0.5d0*faczz*fcoi*faclr/(hbw3*hbw4)
401  270 CONTINUE
402 
403  ELSEIF(isub.EQ.23) THEN
404 C...f + fbar' -> Z0 + W+/- (Z0 only, i.e. no gamma* admixture.)
405  faczw=comfac*0.5d0*(aem/xw)**2
406  faczw=faczw*wids(23,2)
407  thuh=max(th*uh-sqm3*sqm4,sh*ckin(3)**2)
408  facbw=1d0/((sh-sqmw)**2+gmmw**2)
409  DO 290 i=mmin1,mmax1
410  ia=iabs(i)
411  IF(i.EQ.0.OR.ia.GT.20.OR.kfac(1,i).EQ.0) goto 290
412  DO 280 j=mmin2,mmax2
413  ja=iabs(j)
414  IF(j.EQ.0.OR.ja.GT.20.OR.kfac(2,j).EQ.0) goto 280
415  IF(i*j.GT.0.OR.mod(ia+ja,2).EQ.0) goto 280
416  IF((ia.LE.10.AND.ja.GT.10).OR.(ia.GT.10.AND.ja.LE.10))
417  & goto 280
418  kchw=(kchg(ia,1)*isign(1,i)+kchg(ja,1)*isign(1,j))/3
419  ei=kchg(ia,1)/3d0
420  ai=sign(1d0,ei+0.1d0)
421  vi=ai-4d0*ei*xwv
422  ej=kchg(ja,1)/3d0
423  aj=sign(1d0,ej+0.1d0)
424  vj=aj-4d0*ej*xwv
425  IF(vi+ai.GT.0) THEN
426  visav=vi
427  aisav=ai
428  vi=vj
429  ai=aj
430  vj=visav
431  aj=aisav
432  ENDIF
433  fckm=1d0
434  IF(ia.LE.10) fckm=vckm((ia+1)/2,(ja+1)/2)
435  fcoi=1d0
436  IF(ia.LE.10) fcoi=faca/3d0
437  nchn=nchn+1
438  isig(nchn,1)=i
439  isig(nchn,2)=j
440  isig(nchn,3)=1
441  sigh(nchn)=faczw*fcoi*fckm*(facbw*((9d0-8d0*xw)/4d0*thuh+
442  & (8d0*xw-6d0)/4d0*sh*(sqm3+sqm4))+(thuh-sh*(sqm3+sqm4))*
443  & (sh-sqmw)*facbw*0.5d0*((vj+aj)/th-(vi+ai)/uh)+
444  & thuh/(16d0*xw1)*((vj+aj)**2/th2+(vi+ai)**2/uh2)+
445  & sh*(sqm3+sqm4)/(8d0*xw1)*(vi+ai)*(vj+aj)/(th*uh))*
446  & wids(24,(5-kchw)/2)
447 C***Protect against slightly negative cross sections. (Reason yet to be
448 C***sorted out. One possibility: addition of width to the W propagator.)
449  sigh(nchn)=max(0d0,sigh(nchn))
450  280 CONTINUE
451  290 CONTINUE
452 
453  ELSEIF(isub.EQ.25) THEN
454 C...f + fbar -> W+ + W-
455 C...Propagators: Z0, W+- as simulated in PYOFSH and as desired
456  gmmzc=gmmz
457  hbwzc=sh**2/((sh-sqmz)**2+gmmzc**2)
458  hbw3=gmmw/((sqm3-sqmw)**2+gmmw**2)
459  CALL pywidt(24,sqm3,wdtp,wdte)
460  gmmw3=sqrt(sqm3)*wdtp(0)
461  hbw3c=gmmw3/((sqm3-sqmw)**2+gmmw3**2)
462  hbw4=gmmw/((sqm4-sqmw)**2+gmmw**2)
463  CALL pywidt(24,sqm4,wdtp,wdte)
464  gmmw4=sqrt(sqm4)*wdtp(0)
465  hbw4c=gmmw4/((sqm4-sqmw)**2+gmmw4**2)
466 C...Kinematical functions
467  thuh=max(th*uh-sqm3*sqm4,sh*ckin(3)**2)
468  thuh34=(2d0*sh*(sqm3+sqm4)+thuh)/(sqm3*sqm4)
469  gs=(((sh-sqm3-sqm4)**2-4d0*sqm3*sqm4)*thuh34+12d0*thuh)/sh2
470  gt=thuh34+4d0*thuh/th2
471  gst=((sh-sqm3-sqm4)*thuh34+4d0*(sh*(sqm3+sqm4)-thuh)/th)/sh
472  gu=thuh34+4d0*thuh/uh2
473  gsu=((sh-sqm3-sqm4)*thuh34+4d0*(sh*(sqm3+sqm4)-thuh)/uh)/sh
474 C...Common factors and couplings
475  facww=comfac*(hbw3c/hbw3)*(hbw4c/hbw4)
476  facww=facww*wids(24,1)
477  cgg=aem**2/2d0
478  cgz=aem**2/(4d0*xw)*hbwzc*(1d0-sqmz/sh)
479  czz=aem**2/(32d0*xw**2)*hbwzc
480  cng=aem**2/(4d0*xw)
481  cnz=aem**2/(16d0*xw**2)*hbwzc*(1d0-sqmz/sh)
482  cnn=aem**2/(16d0*xw**2)
483 C...Coulomb factor for W+W- pair
484  IF(mstp(40).GE.1.AND.mstp(40).LE.3) THEN
485  coule=(sh-4d0*sqmw)/(4d0*pmas(24,1))
486  coulp=max(1d-10,0.5d0*be34*sqrt(sh))
487  IF(coule.LT.100d0*pmas(24,2)) THEN
488  coulp1=sqrt(0.5d0*pmas(24,1)*(sqrt(coule**2+
489  & pmas(24,2)**2)-coule))
490  ELSE
491  coulp1=sqrt(0.5d0*pmas(24,1)*(0.5d0*pmas(24,2)**2/coule))
492  ENDIF
493  IF(coule.GT.-100d0*pmas(24,2)) THEN
494  coulp2=sqrt(0.5d0*pmas(24,1)*(sqrt(coule**2+
495  & pmas(24,2)**2)+coule))
496  ELSE
497  coulp2=sqrt(0.5d0*pmas(24,1)*(0.5d0*pmas(24,2)**2/
498  & abs(coule)))
499  ENDIF
500  IF(mstp(40).EQ.1) THEN
501  couldc=paru(1)-2d0*atan((coulp1**2+coulp2**2-coulp**2)/
502  & max(1d-10,2d0*coulp*coulp1))
503  faccou=1d0+0.5d0*paru(101)*couldc/max(1d-5,be34)
504  ELSEIF(mstp(40).EQ.2) THEN
505  coulck=dcmplx(dble(coulp1),dble(coulp2))
506  coulcp=dcmplx(0d0,dble(coulp))
507  coulcd=(coulck+coulcp)/(coulck-coulcp)
508  coulcr=1d0+dble(paru(101)*sqrt(sh))/
509  & (4d0*coulcp)*log(coulcd)
510  coulcs=dcmplx(0d0,0d0)
511  nstp=100
512  DO 300 istp=1,nstp
513  coulxx=(istp-0.5)/nstp
514  coulcs=coulcs+(1d0/coulxx)*log((1d0+coulxx*coulcd)/
515  & (1d0+coulxx/coulcd))
516  300 CONTINUE
517  coulcr=coulcr+dble(paru(101)**2*sh)/(16d0*coulcp*coulck)*
518  & (coulcs/nstp)
519  faccou=abs(coulcr)**2
520  ELSEIF(mstp(40).EQ.3) THEN
521  couldc=paru(1)-2d0*(1d0-be34)**2*atan((coulp1**2+
522  & coulp2**2-coulp**2)/max(1d-10,2d0*coulp*coulp1))
523  faccou=1d0+0.5d0*paru(101)*couldc/max(1d-5,be34)
524  ENDIF
525  ELSEIF(mstp(40).EQ.4) THEN
526  faccou=1d0+0.5d0*paru(101)*paru(1)/max(1d-5,be34)
527  ELSE
528  faccou=1d0
529  ENDIF
530  vint(95)=faccou
531  facww=facww*faccou
532 C...Loop over allowed flavours
533  DO 310 i=mmina,mmaxa
534  IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 310
535  ei=kchg(iabs(i),1)/3d0
536  ai=sign(1d0,ei+0.1d0)
537  vi=ai-4d0*ei*xwv
538  fcoi=1d0
539  IF(iabs(i).LE.10) fcoi=faca/3d0
540  IF(mstp(50).LE.0.OR.iabs(i).LE.10) THEN
541  IF(ai.LT.0d0) THEN
542  dsigww=(cgg*ei**2+cgz*vi*ei+czz*(vi**2+ai**2))*gs+
543  & (cng*ei+cnz*(vi+ai))*gst+cnn*gt
544  ELSE
545  dsigww=(cgg*ei**2+cgz*vi*ei+czz*(vi**2+ai**2))*gs-
546  & (cng*ei+cnz*(vi+ai))*gsu+cnn*gu
547  ENDIF
548  ELSE
549  xmw02=0.5d0*(sqm3+sqm4)-0.25d0*(sqm3-sqm4)**2/sh
550  bet=sqrt(1d0-4d0*xmw02/sh)
551  gat=1d0/sqrt(1d0-bet**2)
552  sthe2=1d0-cth**2
553  ampzg=bet**3*(16d0+(4d0*bet**2*gat**2+3d0/gat**2)*sthe2)
554  ampnu=bet*(2d0+bet**2*gat**2*sthe2/2d0+
555  & 2d0*bet**2*(1d0-bet**2)*sthe2/(1d0-2d0*bet*cth+bet**2)**2)
556  ampng=bet*((1d0+bet**2)*(4d0+bet**2*gat**2*sthe2)+
557  & 2d0*(1d0-bet**2)*(bet**2*sthe2-2d0*(1d0-bet**2))/
558  & (1d0-2d0*bet*cth+bet**2))
559  propi1=(0.25d0*sqmz/xmw02)*hbwzc*(1d0-sqmz/sh)
560  propi2=(0.25d0*sqmz/xmw02)**2*hbwzc
561  a0=(2d0*(xmw02/sqmz)-(1d0-bet**2)*xw)*poll
562  a1=(2d0*(xmw02/sqmz)**2-2*xmw02/sqmz*(1d0-bet**2)*xw)*poll
563  a2=(1d0-bet**2)**2*xw**2*(polr+poll)/2d0
564  atot=ampnu*poll+(a1+a2)*propi2*ampzg-a0*propi1*ampng
565  atot=atot*cnn/sqmw*sh/bet*2d0
566  dsigww=atot
567  ENDIF
568  nchn=nchn+1
569  isig(nchn,1)=i
570  isig(nchn,2)=-i
571  isig(nchn,3)=1
572  sigh(nchn)=facww*fcoi*dsigww
573  310 CONTINUE
574 
575  ELSEIF(isub.EQ.30) THEN
576 C...f + g -> f + (gamma*/Z0) (q + g -> q + (gamma*/Z0) only)
577  fzq=comfac*faca*as*aem*(1d0/3d0)*(sh2+uh2+2d0*sqm4*th)/
578  & (-sh*uh)
579 C...gamma, gamma/Z interference and Z couplings to final fermion pairs
580  hfgg=0d0
581  hfgz=0d0
582  hfzz=0d0
583  radc4=1d0+pyalps(sqm4)/paru(1)
584  DO 320 i=1,min(16,mdcy(23,3))
585  idc=i+mdcy(23,2)-1
586  IF(mdme(idc,1).LT.0) goto 320
587  imdm=0
588  IF(mdme(idc,1).EQ.1.OR.mdme(idc,1).EQ.2.OR.mdme(idc,1).EQ.4)
589  & imdm=1
590  IF(i.LE.8) THEN
591  ef=kchg(i,1)/3d0
592  af=sign(1d0,ef+0.1d0)
593  vf=af-4d0*ef*xwv
594  ELSEIF(i.LE.16) THEN
595  ef=kchg(i+2,1)/3d0
596  af=sign(1d0,ef+0.1d0)
597  vf=af-4d0*ef*xwv
598  ENDIF
599  rm1=pmas(iabs(kfdp(idc,1)),1)**2/sqm4
600  IF(4d0*rm1.LT.1d0) THEN
601  fcof=1d0
602  IF(i.LE.8) fcof=3d0*radc4
603  be34=sqrt(max(0d0,1d0-4d0*rm1))
604  IF(imdm.EQ.1) THEN
605  hfgg=hfgg+fcof*ef**2*(1d0+2d0*rm1)*be34
606  hfgz=hfgz+fcof*ef*vf*(1d0+2d0*rm1)*be34
607  hfzz=hfzz+fcof*(vf**2*(1d0+2d0*rm1)+
608  & af**2*(1d0-4d0*rm1))*be34
609  ENDIF
610  ENDIF
611  320 CONTINUE
612 C...Propagators: as simulated in PYOFSH and as desired
613  hbw4=(1d0/paru(1))*gmmz/((sqm4-sqmz)**2+gmmz**2)
614  mint15=mint(15)
615  mint(15)=1
616  mint(61)=1
617  CALL pywidt(23,sqm4,wdtp,wdte)
618  mint(15)=mint15
619  hfaem=(paru(108)/paru(2))*(2d0/3d0)
620  hfgg=hfgg*hfaem*vint(111)/sqm4
621  hfgz=hfgz*hfaem*vint(112)/sqm4
622  hfzz=hfzz*hfaem*vint(114)/sqm4
623 C...Loop over flavours; consider full gamma/Z structure
624  DO 340 i=mmina,mmaxa
625  IF(i.EQ.0.OR.iabs(i).GT.mstp(58)) goto 340
626  ei=kchg(iabs(i),1)/3d0
627  ai=sign(1d0,ei)
628  vi=ai-4d0*ei*xwv
629  faczq=fzq*(ei**2*hfgg+ei*vi*hfgz+
630  & (vi**2+ai**2)*hfzz)/hbw4
631  DO 330 isde=1,2
632  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 330
633  IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 330
634  nchn=nchn+1
635  isig(nchn,isde)=i
636  isig(nchn,3-isde)=21
637  isig(nchn,3)=1
638  sigh(nchn)=faczq
639  330 CONTINUE
640  340 CONTINUE
641 
642  ELSEIF(isub.EQ.31) THEN
643 C...f + g -> f' + W+/- (q + g -> q' + W+/- only)
644  facwq=comfac*faca*as*aem/xw*1d0/12d0*
645  & (sh2+uh2+2d0*sqm4*th)/(-sh*uh)
646 C...Propagators: as simulated in PYOFSH and as desired
647  hbw4=gmmw/((sqm4-sqmw)**2+gmmw**2)
648  CALL pywidt(24,sqm4,wdtp,wdte)
649  gmmwc=sqrt(sqm4)*wdtp(0)
650  hbw4c=gmmwc/((sqm4-sqmw)**2+gmmwc**2)
651  facwq=facwq*hbw4c/hbw4
652  DO 360 i=mmina,mmaxa
653  IF(i.EQ.0.OR.iabs(i).GT.mstp(58)) goto 360
654  ia=iabs(i)
655  kchw=isign(1,kchg(ia,1)*isign(1,i))
656  widsc=(wdte(0,1)+wdte(0,(5-kchw)/2)+wdte(0,4))/wdtp(0)
657  DO 350 isde=1,2
658  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 350
659  IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 350
660  nchn=nchn+1
661  isig(nchn,isde)=i
662  isig(nchn,3-isde)=21
663  isig(nchn,3)=1
664  sigh(nchn)=facwq*vint(180+i)*widsc
665  350 CONTINUE
666  360 CONTINUE
667 
668  ELSEIF(isub.EQ.35) THEN
669 C...f + gamma -> f + (gamma*/Z0)
670  IF(mint(15).EQ.22.AND.vint(3).LT.0d0) THEN
671  fzqn=sh2+uh2+2d0*(sqm4-vint(3)**2)*th
672  fzqdtm=vint(3)**2*sqm4-sh*(uh-vint(4)**2)
673  ELSEIF(mint(16).EQ.22.AND.vint(4).LT.0d0) THEN
674  fzqn=sh2+uh2+2d0*(sqm4-vint(4)**2)*th
675  fzqdtm=vint(4)**2*sqm4-sh*(uh-vint(3)**2)
676  ELSE
677  fzqn=sh2+uh2+2d0*sqm4*th
678  fzqdtm=-sh*uh
679  ENDIF
680  fzqn=comfac*2d0*aem**2*max(0d0,fzqn)
681 C...gamma, gamma/Z interference and Z couplings to final fermion pairs
682  hfgg=0d0
683  hfgz=0d0
684  hfzz=0d0
685  radc4=1d0+pyalps(sqm4)/paru(1)
686  DO 370 i=1,min(16,mdcy(23,3))
687  idc=i+mdcy(23,2)-1
688  IF(mdme(idc,1).LT.0) goto 370
689  imdm=0
690  IF(mdme(idc,1).EQ.1.OR.mdme(idc,1).EQ.2.OR.mdme(idc,1).EQ.4)
691  & imdm=1
692  IF(i.LE.8) THEN
693  ef=kchg(i,1)/3d0
694  af=sign(1d0,ef+0.1d0)
695  vf=af-4d0*ef*xwv
696  ELSEIF(i.LE.16) THEN
697  ef=kchg(i+2,1)/3d0
698  af=sign(1d0,ef+0.1d0)
699  vf=af-4d0*ef*xwv
700  ENDIF
701  rm1=pmas(iabs(kfdp(idc,1)),1)**2/sqm4
702  IF(4d0*rm1.LT.1d0) THEN
703  fcof=1d0
704  IF(i.LE.8) fcof=3d0*radc4
705  be34=sqrt(max(0d0,1d0-4d0*rm1))
706  IF(imdm.EQ.1) THEN
707  hfgg=hfgg+fcof*ef**2*(1d0+2d0*rm1)*be34
708  hfgz=hfgz+fcof*ef*vf*(1d0+2d0*rm1)*be34
709  hfzz=hfzz+fcof*(vf**2*(1d0+2d0*rm1)+
710  & af**2*(1d0-4d0*rm1))*be34
711  ENDIF
712  ENDIF
713  370 CONTINUE
714 C...Propagators: as simulated in PYOFSH and as desired
715  hbw4=(1d0/paru(1))*gmmz/((sqm4-sqmz)**2+gmmz**2)
716  mint15=mint(15)
717  mint(15)=1
718  mint(61)=1
719  CALL pywidt(23,sqm4,wdtp,wdte)
720  mint(15)=mint15
721  hfaem=(paru(108)/paru(2))*(2d0/3d0)
722  hfgg=hfgg*hfaem*vint(111)/sqm4
723  hfgz=hfgz*hfaem*vint(112)/sqm4
724  hfzz=hfzz*hfaem*vint(114)/sqm4
725 C...Loop over flavours; consider full gamma/Z structure
726  DO 390 i=mmina,mmaxa
727  IF(i.EQ.0) goto 390
728  ei=kchg(iabs(i),1)/3d0
729  ai=sign(1d0,ei)
730  vi=ai-4d0*ei*xwv
731  faczq=ei**2*(ei**2*hfgg+ei*vi*hfgz+
732  & (vi**2+ai**2)*hfzz)/hbw4
733  fzqd=max(pmas(iabs(i),1)**2*sqm4,fzqdtm)
734  DO 380 isde=1,2
735  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,22).EQ.0) goto 380
736  IF(isde.EQ.2.AND.kfac(1,22)*kfac(2,i).EQ.0) goto 380
737  nchn=nchn+1
738  isig(nchn,isde)=i
739  isig(nchn,3-isde)=22
740  isig(nchn,3)=1
741  sigh(nchn)=faczq*fzqn/fzqd
742  380 CONTINUE
743  390 CONTINUE
744 
745  ELSEIF(isub.EQ.36) THEN
746 C...f + gamma -> f' + W+/-
747  fwq=comfac*aem**2/(2d0*xw)*
748  & (sh2+uh2+2d0*sqm4*th)/(sqpth*sqm4-sh*uh)
749 C...Propagators: as simulated in PYOFSH and as desired
750  hbw4=gmmw/((sqm4-sqmw)**2+gmmw**2)
751  CALL pywidt(24,sqm4,wdtp,wdte)
752  gmmwc=sqrt(sqm4)*wdtp(0)
753  hbw4c=gmmwc/((sqm4-sqmw)**2+gmmwc**2)
754  fwq=fwq*hbw4c/hbw4
755  DO 410 i=mmina,mmaxa
756  IF(i.EQ.0) goto 410
757  ia=iabs(i)
758  eia=abs(kchg(iabs(i),1)/3d0)
759  facwq=fwq*(eia-sh/(sh+uh))**2
760  kchw=isign(1,kchg(ia,1)*isign(1,i))
761  widsc=(wdte(0,1)+wdte(0,(5-kchw)/2)+wdte(0,4))/wdtp(0)
762  DO 400 isde=1,2
763  IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,22).EQ.0) goto 400
764  IF(isde.EQ.2.AND.kfac(1,22)*kfac(2,i).EQ.0) goto 400
765  nchn=nchn+1
766  isig(nchn,isde)=i
767  isig(nchn,3-isde)=22
768  isig(nchn,3)=1
769  sigh(nchn)=facwq*vint(180+i)*widsc
770  400 CONTINUE
771  410 CONTINUE
772  ENDIF
773 
774  ELSEIF(isub.LE.100) THEN
775  IF(isub.EQ.69) THEN
776 C...gamma + gamma -> W+ + W-
777  sqmwe=max(0.5d0*sqmw,sqrt(sqm3*sqm4))
778  fprop=sh2/((sqmwe-th)*(sqmwe-uh))
779  facww=comfac*6d0*aem**2*(1d0-fprop*(4d0/3d0+2d0*sqmwe/sh)+
780  & fprop**2*(2d0/3d0+2d0*(sqmwe/sh)**2))*wids(24,1)
781  IF(kfac(1,22)*kfac(2,22).EQ.0) goto 420
782  nchn=nchn+1
783  isig(nchn,1)=22
784  isig(nchn,2)=22
785  isig(nchn,3)=1
786  sigh(nchn)=facww
787  420 CONTINUE
788 
789  ELSEIF(isub.EQ.70) THEN
790 C...gamma + W+/- -> Z0 + W+/-
791  sqmwe=max(0.5d0*sqmw,sqrt(sqm3*sqm4))
792  fprop=(th-sqmwe)**2/(-sh*(sqmwe-uh))
793  faczw=comfac*6d0*aem**2*(xw1/xw)*
794  & (1d0-fprop*(4d0/3d0+2d0*sqmwe/(th-sqmwe))+
795  & fprop**2*(2d0/3d0+2d0*(sqmwe/(th-sqmwe))**2))*wids(23,2)
796  DO 440 kchw=1,-1,-2
797  DO 430 isde=1,2
798  IF(kfac(isde,22)*kfac(3-isde,24*kchw).EQ.0) goto 430
799  nchn=nchn+1
800  isig(nchn,isde)=22
801  isig(nchn,3-isde)=24*kchw
802  isig(nchn,3)=1
803  sigh(nchn)=faczw*wids(24,(5-kchw)/2)
804  430 CONTINUE
805  440 CONTINUE
806  ENDIF
807  ENDIF
808 
809  RETURN
810  END