Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pywidt.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pywidt.f
1 
2 C*********************************************************************
3 
4 C...PYWIDT
5 C...Calculates full and partial widths of resonances.
6 
7  SUBROUTINE pywidt(KFLR,SH,WDTP,WDTE)
8 
9 C...Double precision and integer declarations.
10  IMPLICIT DOUBLE PRECISION(a-h, o-z)
11  IMPLICIT INTEGER(i-n)
12  INTEGER pyk,pychge,pycomp
13 C...Parameter statement to help give large particle numbers.
14  parameter(ksusy1=1000000,ksusy2=2000000,ktechn=3000000,
15  &kexcit=4000000,kdimen=5000000)
16 C...Commonblocks.
17  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
18  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
19  common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
20  common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
21  common/pypars/mstp(200),parp(200),msti(200),pari(200)
22  common/pyint1/mint(400),vint(400)
23  common/pyint4/mwid(500),wids(500,5)
24  common/pymssm/imss(0:99),rmss(0:99)
25  common/pyssmt/zmix(4,4),umix(2,2),vmix(2,2),smz(4),smw(2),
26  &sfmix(16,4),zmixi(4,4),umixi(2,2),vmixi(2,2)
27  common/pytcsm/itcm(0:99),rtcm(0:99)
28  SAVE /pydat1/,/pydat2/,/pydat3/,/pysubs/,/pypars/,/pyint1/,
29  &/pyint4/,/pymssm/,/pyssmt/,/pytcsm/
30 C...Local arrays and saved variables.
31  COMPLEX*16 zmixc(4,4),al,bl,ar,br,fl,fr
32  dimension wdtp(0:400),wdte(0:400,0:5),mofsv(3,2),widwsv(3,2),
33  &wid2sv(3,2),wdtpp(0:400),wdtep(0:400,0:5)
34  SAVE mofsv,widwsv,wid2sv
35  DATA mofsv/6*0/,widwsv/6*0d0/,wid2sv/6*0d0/
36 
37 C...Compressed code and sign; mass.
38  kfla=iabs(kflr)
39  kfls=isign(1,kflr)
40  kc=pycomp(kfla)
41  shr=sqrt(sh)
42  pmr=pmas(kc,1)
43 
44 C...Reset width information.
45  DO 110 i=0,mdcy(kc,3)
46  wdtp(i)=0d0
47  DO 100 j=0,5
48  wdte(i,j)=0d0
49  100 CONTINUE
50  110 CONTINUE
51 
52 C...Allow for fudge factor to rescale resonance width.
53  fudge=1d0
54  IF(mstp(110).NE.0.AND.(mwid(kc).EQ.1.OR.mwid(kc).EQ.2.OR.
55  &(mwid(kc).EQ.3.AND.mint(63).EQ.1))) THEN
56  IF(mstp(110).EQ.kfla) THEN
57  fudge=parp(110)
58  ELSEIF(mstp(110).EQ.-1) THEN
59  IF(kfla.NE.6.AND.kfla.NE.23.AND.kfla.NE.24) fudge=parp(110)
60  ELSEIF(mstp(110).EQ.-2) THEN
61  fudge=parp(110)
62  ENDIF
63  ENDIF
64 
65 C...Not to be treated as a resonance: return.
66  IF((mwid(kc).LE.0.OR.mwid(kc).GE.4).AND.kfla.NE.21.AND.
67  &kfla.NE.22) THEN
68  wdtp(0)=1d0
69  wdte(0,0)=1d0
70  mint(61)=0
71  mint(62)=0
72  mint(63)=0
73  RETURN
74 
75 C...Treatment as a resonance based on tabulated branching ratios.
76  ELSEIF(mwid(kc).EQ.2.OR.(mwid(kc).EQ.3.AND.mint(63).EQ.0)) THEN
77 C...Loop over possible decay channels; skip irrelevant ones.
78  DO 120 i=1,mdcy(kc,3)
79  idc=i+mdcy(kc,2)-1
80  IF(mdme(idc,1).LT.0) goto 120
81 
82 C...Read out decay products and nominal masses.
83  kfd1=kfdp(idc,1)
84  kfc1=pycomp(kfd1)
85  IF(kchg(kfc1,3).EQ.1) kfd1=kfls*kfd1
86  pm1=pmas(kfc1,1)
87  kfd2=kfdp(idc,2)
88  kfc2=pycomp(kfd2)
89  IF(kchg(kfc2,3).EQ.1) kfd2=kfls*kfd2
90  pm2=pmas(kfc2,1)
91  kfd3=kfdp(idc,3)
92  pm3=0d0
93  IF(kfd3.NE.0) THEN
94  kfc3=pycomp(kfd3)
95  IF(kchg(kfc3,3).EQ.1) kfd3=kfls*kfd3
96  pm3=pmas(kfc3,1)
97  ENDIF
98 
99 C...Naive partial width and alternative threshold factors.
100  wdtp(i)=pmas(kc,2)*brat(idc)*(shr/pmr)
101  IF(mdme(idc,2).GE.51.AND.mdme(idc,2).LE.53.AND.
102  & pm1+pm2+pm3.GE.shr) THEN
103  wdtp(i)=0d0
104  ELSEIF(mdme(idc,2).EQ.52.AND.kfd3.EQ.0) THEN
105  wdtp(i)=wdtp(i)*sqrt(max(0d0,(sh-pm1**2-pm2**2)**2-
106  & 4d0*pm1**2*pm2**2))/sh
107  ELSEIF(mdme(idc,2).EQ.52) THEN
108  pma=max(pm1,pm2,pm3)
109  pmc=min(pm1,pm2,pm3)
110  pmb=pm1+pm2+pm3-pma-pmc
111  pmbc=pmb+pmc+0.5d0*(shr-pma-pmc-pmc)
112  pman=pma**2/sh
113  pmbn=pmb**2/sh
114  pmcn=pmc**2/sh
115  pmbcn=pmbc**2/sh
116  wdtp(i)=wdtp(i)*sqrt(max(0d0,
117  & ((1d0-pman-pmbcn)**2-4d0*pman*pmbcn)*
118  & ((pmbcn-pmbn-pmcn)**2-4d0*pmbn*pmcn)))*
119  & ((shr-pma)**2-(pmb+pmc)**2)*
120  & (1d0+0.25d0*(pma+pmb+pmc)/shr)/
121  & ((1d0-pmbcn)*pmbcn*sh)
122  ELSEIF(mdme(idc,2).EQ.53.AND.kfd3.EQ.0) THEN
123  wdtp(i)=wdtp(i)*sqrt(
124  & max(0d0,(sh-pm1**2-pm2**2)**2-4d0*pm1**2*pm2**2)/
125  & max(1d-4,(pmr**2-pm1**2-pm2**2)**2-4d0*pm1**2*pm2**2))
126  ELSEIF(mdme(idc,2).EQ.53) THEN
127  pma=max(pm1,pm2,pm3)
128  pmc=min(pm1,pm2,pm3)
129  pmb=pm1+pm2+pm3-pma-pmc
130  pmbc=pmb+pmc+0.5d0*(shr-pma-pmb-pmc)
131  pman=pma**2/sh
132  pmbn=pmb**2/sh
133  pmcn=pmc**2/sh
134  pmbcn=pmbc**2/sh
135  facact=sqrt(max(0d0,
136  & ((1d0-pman-pmbcn)**2-4d0*pman*pmbcn)*
137  & ((pmbcn-pmbn-pmcn)**2-4d0*pmbn*pmcn)))*
138  & ((shr-pma)**2-(pmb+pmc)**2)*
139  & (1d0+0.25d0*(pma+pmb+pmc)/shr)/
140  & ((1d0-pmbcn)*pmbcn*sh)
141  pmbc=pmb+pmc+0.5d0*(pmr-pma-pmb-pmc)
142  pman=pma**2/pmr**2
143  pmbn=pmb**2/pmr**2
144  pmcn=pmc**2/pmr**2
145  pmbcn=pmbc**2/pmr**2
146  facnom=sqrt(max(0d0,
147  & ((1d0-pman-pmbcn)**2-4d0*pman*pmbcn)*
148  & ((pmbcn-pmbn-pmcn)**2-4d0*pmbn*pmcn)))*
149  & ((pmr-pma)**2-(pmb+pmc)**2)*
150  & (1d0+0.25d0*(pma+pmb+pmc)/pmr)/
151  & ((1d0-pmbcn)*pmbcn*pmr**2)
152  wdtp(i)=wdtp(i)*facact/max(1d-6,facnom)
153  ENDIF
154  wdtp(i)=fudge*wdtp(i)
155  wdtp(0)=wdtp(0)+wdtp(i)
156 
157 C...Calculate secondary width (at most two identical/opposite).
158  wid2=1d0
159  IF(mdme(idc,1).GT.0) THEN
160  IF(kfd2.EQ.kfd1) THEN
161  IF(kchg(kfc1,3).EQ.0) THEN
162  wid2=wids(kfc1,1)
163  ELSEIF(kfd1.GT.0) THEN
164  wid2=wids(kfc1,4)
165  ELSE
166  wid2=wids(kfc1,5)
167  ENDIF
168  IF(kfd3.GT.0) THEN
169  wid2=wid2*wids(kfc3,2)
170  ELSEIF(kfd3.LT.0) THEN
171  wid2=wid2*wids(kfc3,3)
172  ENDIF
173  ELSEIF(kfd2.EQ.-kfd1) THEN
174  wid2=wids(kfc1,1)
175  IF(kfd3.GT.0) THEN
176  wid2=wid2*wids(kfc3,2)
177  ELSEIF(kfd3.LT.0) THEN
178  wid2=wid2*wids(kfc3,3)
179  ENDIF
180  ELSEIF(kfd3.EQ.kfd1) THEN
181  IF(kchg(kfc1,3).EQ.0) THEN
182  wid2=wids(kfc1,1)
183  ELSEIF(kfd1.GT.0) THEN
184  wid2=wids(kfc1,4)
185  ELSE
186  wid2=wids(kfc1,5)
187  ENDIF
188  IF(kfd2.GT.0) THEN
189  wid2=wid2*wids(kfc2,2)
190  ELSEIF(kfd2.LT.0) THEN
191  wid2=wid2*wids(kfc2,3)
192  ENDIF
193  ELSEIF(kfd3.EQ.-kfd1) THEN
194  wid2=wids(kfc1,1)
195  IF(kfd2.GT.0) THEN
196  wid2=wid2*wids(kfc2,2)
197  ELSEIF(kfd2.LT.0) THEN
198  wid2=wid2*wids(kfc2,3)
199  ENDIF
200  ELSEIF(kfd3.EQ.kfd2) THEN
201  IF(kchg(kfc2,3).EQ.0) THEN
202  wid2=wids(kfc2,1)
203  ELSEIF(kfd2.GT.0) THEN
204  wid2=wids(kfc2,4)
205  ELSE
206  wid2=wids(kfc2,5)
207  ENDIF
208  IF(kfd1.GT.0) THEN
209  wid2=wid2*wids(kfc1,2)
210  ELSEIF(kfd1.LT.0) THEN
211  wid2=wid2*wids(kfc1,3)
212  ENDIF
213  ELSEIF(kfd3.EQ.-kfd2) THEN
214  wid2=wids(kfc2,1)
215  IF(kfd1.GT.0) THEN
216  wid2=wid2*wids(kfc1,2)
217  ELSEIF(kfd1.LT.0) THEN
218  wid2=wid2*wids(kfc1,3)
219  ENDIF
220  ELSE
221  IF(kfd1.GT.0) THEN
222  wid2=wids(kfc1,2)
223  ELSE
224  wid2=wids(kfc1,3)
225  ENDIF
226  IF(kfd2.GT.0) THEN
227  wid2=wid2*wids(kfc2,2)
228  ELSE
229  wid2=wid2*wids(kfc2,3)
230  ENDIF
231  IF(kfd3.GT.0) THEN
232  wid2=wid2*wids(kfc3,2)
233  ELSEIF(kfd3.LT.0) THEN
234  wid2=wid2*wids(kfc3,3)
235  ENDIF
236  ENDIF
237 
238 C...Store effective widths according to case.
239  wdte(i,mdme(idc,1))=wdtp(i)*wid2
240  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
241  wdte(i,0)=wdte(i,mdme(idc,1))
242  wdte(0,0)=wdte(0,0)+wdte(i,0)
243  ENDIF
244  120 CONTINUE
245 C...Return.
246  mint(61)=0
247  mint(62)=0
248  mint(63)=0
249  RETURN
250  ENDIF
251 
252 C...Here begins detailed dynamical calculation of resonance widths.
253 C...Shared treatment of Higgs states.
254  kfhigg=25
255  ihigg=1
256  IF(kfla.EQ.35.OR.kfla.EQ.36) THEN
257  kfhigg=kfla
258  ihigg=kfla-33
259  ENDIF
260 
261 C...Common electroweak and strong constants.
262  xw=paru(102)
263  xwv=xw
264  IF(mstp(8).GE.2) xw=1d0-(pmas(24,1)/pmas(23,1))**2
265  xw1=1d0-xw
266  aem=pyalem(sh)
267  IF(mstp(8).GE.1) aem=sqrt(2d0)*paru(105)*pmas(24,1)**2*xw/paru(1)
268  as=pyalps(sh)
269  radc=1d0+as/paru(1)
270 
271  IF(kfla.EQ.6) THEN
272 C...t quark.
273  fac=(aem/(16d0*xw))*(sh/pmas(24,1)**2)*shr
274  radct=1d0-2.5d0*as/paru(1)
275  DO 140 i=1,mdcy(kc,3)
276  idc=i+mdcy(kc,2)-1
277  IF(mdme(idc,1).LT.0) goto 140
278  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
279  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
280  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 140
281  wid2=1d0
282  IF(i.GE.4.AND.i.LE.7) THEN
283 C...t -> W + q; including approximate QCD correction factor.
284  wdtp(i)=fac*vckm(3,i-3)*radct*
285  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))*
286  & ((1d0-rm2)**2+(1d0+rm2)*rm1-2d0*rm1**2)
287  IF(kflr.GT.0) THEN
288  wid2=wids(24,2)
289  IF(i.EQ.7) wid2=wid2*wids(7,2)
290  ELSE
291  wid2=wids(24,3)
292  IF(i.EQ.7) wid2=wid2*wids(7,3)
293  ENDIF
294  ELSEIF(i.EQ.9) THEN
295 C...t -> H + b.
296  rm2r=pymrun(kfdp(idc,2),sh)**2/sh
297  wdtp(i)=fac*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))*
298  & ((1d0+rm2-rm1)*(rm2r*paru(141)**2+1d0/paru(141)**2)+
299  & 4d0*sqrt(rm2r*rm2))
300  wid2=wids(37,2)
301  IF(kflr.LT.0) wid2=wids(37,3)
302 CMRENNA++
303  ELSEIF(i.GE.10.AND.i.LE.13.AND.imss(1).NE.0) THEN
304 C...t -> ~t + ~chi_i0, i = 1, 2, 3 or 4.
305  beta=atan(rmss(5))
306  sinb=sin(beta)
307  tanw=sqrt(paru(102)/(1d0-paru(102)))
308  et=kchg(6,1)/3d0
309  t3l=sign(0.5d0,et)
310  kfc1=pycomp(kfdp(idc,1))
311  kfc2=pycomp(kfdp(idc,2))
312  pmnchi=pmas(kfc1,1)
313  pmstop=pmas(kfc2,1)
314  IF(shr.GT.pmnchi+pmstop) THEN
315  iz=i-9
316  DO 130 ik=1,4
317  zmixc(iz,ik)=dcmplx(zmix(iz,ik),zmixi(iz,ik))
318  130 CONTINUE
319  al=shr*dconjg(zmixc(iz,4))/(2.0d0*pmas(24,1)*sinb)
320  ar=-et*zmixc(iz,1)*tanw
321  bl=t3l*(zmixc(iz,2)-zmixc(iz,1)*tanw)-ar
322  br=al
323  fl=sfmix(6,1)*al+sfmix(6,2)*ar
324  fr=sfmix(6,1)*bl+sfmix(6,2)*br
325  pcm=sqrt((sh-(pmnchi+pmstop)**2)*
326  & (sh-(pmnchi-pmstop)**2))/(2d0*shr)
327  wdtp(i)=(0.5d0*pyalem(sh)/paru(102))*pcm*
328  & ((abs(fl)**2+abs(fr)**2)*(sh+pmnchi**2-pmstop**2)+
329  & smz(iz)*4d0*shr*dble(fl*dconjg(fr)))/sh
330  IF(kflr.GT.0) THEN
331  wid2=wids(kfc1,2)*wids(kfc2,2)
332  ELSE
333  wid2=wids(kfc1,2)*wids(kfc2,3)
334  ENDIF
335  ENDIF
336  ELSEIF(i.EQ.14.AND.imss(1).NE.0) THEN
337 C...t -> ~g + ~t
338  kfc1=pycomp(kfdp(idc,1))
339  kfc2=pycomp(kfdp(idc,2))
340  pmnchi=pmas(kfc1,1)
341  pmstop=pmas(kfc2,1)
342  IF(shr.GT.pmnchi+pmstop) THEN
343  rl=sfmix(6,1)
344  rr=-sfmix(6,2)
345  pcm=sqrt((sh-(pmnchi+pmstop)**2)*
346  & (sh-(pmnchi-pmstop)**2))/(2d0*shr)
347  wdtp(i)=4d0/3d0*0.5d0*pyalps(sh)*pcm*((rl**2+rr**2)*
348  & (sh+pmnchi**2-pmstop**2)+pmnchi*4d0*shr*rl*rr)/sh
349  IF(kflr.GT.0) THEN
350  wid2=wids(kfc1,2)*wids(kfc2,2)
351  ELSE
352  wid2=wids(kfc1,2)*wids(kfc2,3)
353  ENDIF
354  ENDIF
355  ELSEIF(i.EQ.15.AND.imss(1).NE.0) THEN
356 C...t -> ~gravitino + ~t
357  xmp2=rmss(29)**2
358  kfc1=pycomp(kfdp(idc,1))
359  xmgr2=pmas(kfc1,1)**2
360  wdtp(i)=sh**2*shr/(96d0*paru(1)*xmp2*xmgr2)*(1d0-rm2)**4
361  kfc2=pycomp(kfdp(idc,2))
362  wid2=wids(kfc2,2)
363  IF(kflr.LT.0) wid2=wids(kfc2,3)
364 CMRENNA--
365  ENDIF
366  wdtp(i)=fudge*wdtp(i)
367  wdtp(0)=wdtp(0)+wdtp(i)
368  IF(mdme(idc,1).GT.0) THEN
369  wdte(i,mdme(idc,1))=wdtp(i)*wid2
370  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
371  wdte(i,0)=wdte(i,mdme(idc,1))
372  wdte(0,0)=wdte(0,0)+wdte(i,0)
373  ENDIF
374  140 CONTINUE
375 
376  ELSEIF(kfla.EQ.7) THEN
377 C...b' quark.
378  fac=(aem/(16d0*xw))*(sh/pmas(24,1)**2)*shr
379  DO 150 i=1,mdcy(kc,3)
380  idc=i+mdcy(kc,2)-1
381  IF(mdme(idc,1).LT.0) goto 150
382  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
383  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
384  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 150
385  wid2=1d0
386  IF(i.GE.4.AND.i.LE.7) THEN
387 C...b' -> W + q.
388  wdtp(i)=fac*vckm(i-3,4)*
389  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))*
390  & ((1d0-rm2)**2+(1d0+rm2)*rm1-2d0*rm1**2)
391  IF(kflr.GT.0) THEN
392  wid2=wids(24,3)
393  IF(i.EQ.6) wid2=wid2*wids(6,2)
394  IF(i.EQ.7) wid2=wid2*wids(8,2)
395  ELSE
396  wid2=wids(24,2)
397  IF(i.EQ.6) wid2=wid2*wids(6,3)
398  IF(i.EQ.7) wid2=wid2*wids(8,3)
399  ENDIF
400  wid2=wids(24,3)
401  IF(kflr.LT.0) wid2=wids(24,2)
402  ELSEIF(i.EQ.9.OR.i.EQ.10) THEN
403 C...b' -> H + q.
404  wdtp(i)=fac*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))*
405  & ((1d0+rm2-rm1)*(paru(141)**2+rm2/paru(141)**2)+4d0*rm2)
406  IF(kflr.GT.0) THEN
407  wid2=wids(37,3)
408  IF(i.EQ.10) wid2=wid2*wids(6,2)
409  ELSE
410  wid2=wids(37,2)
411  IF(i.EQ.10) wid2=wid2*wids(6,3)
412  ENDIF
413  ENDIF
414  wdtp(i)=fudge*wdtp(i)
415  wdtp(0)=wdtp(0)+wdtp(i)
416  IF(mdme(idc,1).GT.0) THEN
417  wdte(i,mdme(idc,1))=wdtp(i)*wid2
418  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
419  wdte(i,0)=wdte(i,mdme(idc,1))
420  wdte(0,0)=wdte(0,0)+wdte(i,0)
421  ENDIF
422  150 CONTINUE
423 
424  ELSEIF(kfla.EQ.8) THEN
425 C...t' quark.
426  fac=(aem/(16d0*xw))*(sh/pmas(24,1)**2)*shr
427  DO 160 i=1,mdcy(kc,3)
428  idc=i+mdcy(kc,2)-1
429  IF(mdme(idc,1).LT.0) goto 160
430  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
431  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
432  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 160
433  wid2=1d0
434  IF(i.GE.4.AND.i.LE.7) THEN
435 C...t' -> W + q.
436  wdtp(i)=fac*vckm(4,i-3)*
437  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))*
438  & ((1d0-rm2)**2+(1d0+rm2)*rm1-2d0*rm1**2)
439  IF(kflr.GT.0) THEN
440  wid2=wids(24,2)
441  IF(i.EQ.7) wid2=wid2*wids(7,2)
442  ELSE
443  wid2=wids(24,3)
444  IF(i.EQ.7) wid2=wid2*wids(7,3)
445  ENDIF
446  ELSEIF(i.EQ.9.OR.i.EQ.10) THEN
447 C...t' -> H + q.
448  wdtp(i)=fac*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))*
449  & ((1d0+rm2-rm1)*(rm2*paru(141)**2+1d0/paru(141)**2)+4d0*rm2)
450  IF(kflr.GT.0) THEN
451  wid2=wids(37,2)
452  IF(i.EQ.10) wid2=wid2*wids(7,2)
453  ELSE
454  wid2=wids(37,3)
455  IF(i.EQ.10) wid2=wid2*wids(7,3)
456  ENDIF
457  ENDIF
458  wdtp(i)=fudge*wdtp(i)
459  wdtp(0)=wdtp(0)+wdtp(i)
460  IF(mdme(idc,1).GT.0) THEN
461  wdte(i,mdme(idc,1))=wdtp(i)*wid2
462  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
463  wdte(i,0)=wdte(i,mdme(idc,1))
464  wdte(0,0)=wdte(0,0)+wdte(i,0)
465  ENDIF
466  160 CONTINUE
467 
468  ELSEIF(kfla.EQ.17) THEN
469 C...tau' lepton.
470  fac=(aem/(16d0*xw))*(sh/pmas(24,1)**2)*shr
471  DO 170 i=1,mdcy(kc,3)
472  idc=i+mdcy(kc,2)-1
473  IF(mdme(idc,1).LT.0) goto 170
474  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
475  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
476  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 170
477  wid2=1d0
478  IF(i.EQ.3) THEN
479 C...tau' -> W + nu'_tau.
480  wdtp(i)=fac*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))*
481  & ((1d0-rm2)**2+(1d0+rm2)*rm1-2d0*rm1**2)
482  IF(kflr.GT.0) THEN
483  wid2=wids(24,3)
484  wid2=wid2*wids(18,2)
485  ELSE
486  wid2=wids(24,2)
487  wid2=wid2*wids(18,3)
488  ENDIF
489  ELSEIF(i.EQ.5) THEN
490 C...tau' -> H + nu'_tau.
491  wdtp(i)=fac*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))*
492  & ((1d0+rm2-rm1)*(paru(141)**2+rm2/paru(141)**2)+4d0*rm2)
493  IF(kflr.GT.0) THEN
494  wid2=wids(37,3)
495  wid2=wid2*wids(18,2)
496  ELSE
497  wid2=wids(37,2)
498  wid2=wid2*wids(18,3)
499  ENDIF
500  ENDIF
501  wdtp(i)=fudge*wdtp(i)
502  wdtp(0)=wdtp(0)+wdtp(i)
503  IF(mdme(idc,1).GT.0) THEN
504  wdte(i,mdme(idc,1))=wdtp(i)*wid2
505  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
506  wdte(i,0)=wdte(i,mdme(idc,1))
507  wdte(0,0)=wdte(0,0)+wdte(i,0)
508  ENDIF
509  170 CONTINUE
510 
511  ELSEIF(kfla.EQ.18) THEN
512 C...nu'_tau neutrino.
513  fac=(aem/(16d0*xw))*(sh/pmas(24,1)**2)*shr
514  DO 180 i=1,mdcy(kc,3)
515  idc=i+mdcy(kc,2)-1
516  IF(mdme(idc,1).LT.0) goto 180
517  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
518  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
519  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 180
520  wid2=1d0
521  IF(i.EQ.2) THEN
522 C...nu'_tau -> W + tau'.
523  wdtp(i)=fac*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))*
524  & ((1d0-rm2)**2+(1d0+rm2)*rm1-2d0*rm1**2)
525  IF(kflr.GT.0) THEN
526  wid2=wids(24,2)
527  wid2=wid2*wids(17,2)
528  ELSE
529  wid2=wids(24,3)
530  wid2=wid2*wids(17,3)
531  ENDIF
532  ELSEIF(i.EQ.3) THEN
533 C...nu'_tau -> H + tau'.
534  wdtp(i)=fac*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))*
535  & ((1d0+rm2-rm1)*(rm2*paru(141)**2+1d0/paru(141)**2)+4d0*rm2)
536  IF(kflr.GT.0) THEN
537  wid2=wids(37,2)
538  wid2=wid2*wids(17,2)
539  ELSE
540  wid2=wids(37,3)
541  wid2=wid2*wids(17,3)
542  ENDIF
543  ENDIF
544  wdtp(i)=fudge*wdtp(i)
545  wdtp(0)=wdtp(0)+wdtp(i)
546  IF(mdme(idc,1).GT.0) THEN
547  wdte(i,mdme(idc,1))=wdtp(i)*wid2
548  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
549  wdte(i,0)=wdte(i,mdme(idc,1))
550  wdte(0,0)=wdte(0,0)+wdte(i,0)
551  ENDIF
552  180 CONTINUE
553 
554  ELSEIF(kfla.EQ.21) THEN
555 C...QCD:
556 C***Note that widths are not given in dimensional quantities here.
557  DO 190 i=1,mdcy(kc,3)
558  idc=i+mdcy(kc,2)-1
559  IF(mdme(idc,1).LT.0) goto 190
560  rm1=pmas(iabs(kfdp(idc,1)),1)**2/sh
561  rm2=pmas(iabs(kfdp(idc,2)),1)**2/sh
562  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 190
563  wid2=1d0
564  IF(i.LE.8) THEN
565 C...QCD -> q + qbar
566  wdtp(i)=(1d0+2d0*rm1)*sqrt(max(0d0,1d0-4d0*rm1))
567  IF(i.EQ.6) wid2=wids(6,1)
568  IF((i.EQ.7.OR.i.EQ.8)) wid2=wids(i,1)
569  ENDIF
570  wdtp(i)=fudge*wdtp(i)
571  wdtp(0)=wdtp(0)+wdtp(i)
572  IF(mdme(idc,1).GT.0) THEN
573  wdte(i,mdme(idc,1))=wdtp(i)*wid2
574  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
575  wdte(i,0)=wdte(i,mdme(idc,1))
576  wdte(0,0)=wdte(0,0)+wdte(i,0)
577  ENDIF
578  190 CONTINUE
579 
580  ELSEIF(kfla.EQ.22) THEN
581 C...QED photon.
582 C***Note that widths are not given in dimensional quantities here.
583  DO 200 i=1,mdcy(kc,3)
584  idc=i+mdcy(kc,2)-1
585  IF(mdme(idc,1).LT.0) goto 200
586  rm1=pmas(iabs(kfdp(idc,1)),1)**2/sh
587  rm2=pmas(iabs(kfdp(idc,2)),1)**2/sh
588  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 200
589  wid2=1d0
590  IF(i.LE.8) THEN
591 C...QED -> q + qbar.
592  ef=kchg(i,1)/3d0
593  fcof=3d0*radc
594  IF(i.GE.6.AND.mstp(35).GE.1) fcof=fcof*pyhfth(sh,sh*rm1,1d0)
595  wdtp(i)=fcof*ef**2*(1d0+2d0*rm1)*sqrt(max(0d0,1d0-4d0*rm1))
596  IF(i.EQ.6) wid2=wids(6,1)
597  IF((i.EQ.7.OR.i.EQ.8)) wid2=wids(i,1)
598  ELSEIF(i.LE.12) THEN
599 C...QED -> l+ + l-.
600  ef=kchg(9+2*(i-8),1)/3d0
601  wdtp(i)=ef**2*(1d0+2d0*rm1)*sqrt(max(0d0,1d0-4d0*rm1))
602  IF(i.EQ.12) wid2=wids(17,1)
603  ENDIF
604  wdtp(i)=fudge*wdtp(i)
605  wdtp(0)=wdtp(0)+wdtp(i)
606  IF(mdme(idc,1).GT.0) THEN
607  wdte(i,mdme(idc,1))=wdtp(i)*wid2
608  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
609  wdte(i,0)=wdte(i,mdme(idc,1))
610  wdte(0,0)=wdte(0,0)+wdte(i,0)
611  ENDIF
612  200 CONTINUE
613 
614  ELSEIF(kfla.EQ.23) THEN
615 C...Z0:
616  icase=1
617  xwc=1d0/(16d0*xw*xw1)
618  fac=(aem*xwc/3d0)*shr
619  210 CONTINUE
620  IF(mint(61).GE.1.AND.icase.EQ.2) THEN
621  vint(111)=0d0
622  vint(112)=0d0
623  vint(114)=0d0
624  ENDIF
625  IF(mint(61).EQ.1.AND.icase.EQ.2) THEN
626  kfi=iabs(mint(15))
627  IF(kfi.GT.20) kfi=iabs(mint(16))
628  ei=kchg(kfi,1)/3d0
629  ai=sign(1d0,ei)
630  vi=ai-4d0*ei*xwv
631  sqmz=pmas(23,1)**2
632  hz=shr*wdtp(0)
633  IF(mstp(43).EQ.1.OR.mstp(43).EQ.3) vint(111)=1d0
634  IF(mstp(43).EQ.3) vint(112)=
635  & 2d0*xwc*sh*(sh-sqmz)/((sh-sqmz)**2+hz**2)
636  IF(mstp(43).EQ.2.OR.mstp(43).EQ.3) vint(114)=
637  & xwc**2*sh**2/((sh-sqmz)**2+hz**2)
638  ENDIF
639  DO 220 i=1,mdcy(kc,3)
640  idc=i+mdcy(kc,2)-1
641  IF(mdme(idc,1).LT.0) goto 220
642  rm1=pmas(iabs(kfdp(idc,1)),1)**2/sh
643  rm2=pmas(iabs(kfdp(idc,2)),1)**2/sh
644  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 220
645  wid2=1d0
646  IF(i.LE.8) THEN
647 C...Z0 -> q + qbar
648  ef=kchg(i,1)/3d0
649  af=sign(1d0,ef+0.1d0)
650  vf=af-4d0*ef*xwv
651  fcof=3d0*radc
652  IF(i.GE.6.AND.mstp(35).GE.1) fcof=fcof*pyhfth(sh,sh*rm1,1d0)
653  IF(i.EQ.6) wid2=wids(6,1)
654  IF((i.EQ.7.OR.i.EQ.8)) wid2=wids(i,1)
655  ELSEIF(i.LE.16) THEN
656 C...Z0 -> l+ + l-, nu + nubar
657  ef=kchg(i+2,1)/3d0
658  af=sign(1d0,ef+0.1d0)
659  vf=af-4d0*ef*xwv
660  fcof=1d0
661  IF((i.EQ.15.OR.i.EQ.16)) wid2=wids(2+i,1)
662  ENDIF
663  be34=sqrt(max(0d0,1d0-4d0*rm1))
664  IF(icase.EQ.1) THEN
665  wdtp(i)=fac*fcof*(vf**2*(1d0+2d0*rm1)+af**2*(1d0-4d0*rm1))*
666  & be34
667  ELSEIF(mint(61).EQ.1.AND.icase.EQ.2) THEN
668  wdtp(i)=fac*fcof*((ei**2*vint(111)*ef**2+ei*vi*vint(112)*
669  & ef*vf+(vi**2+ai**2)*vint(114)*vf**2)*(1d0+2d0*rm1)+
670  & (vi**2+ai**2)*vint(114)*af**2*(1d0-4d0*rm1))*be34
671  ELSEIF(mint(61).EQ.2.AND.icase.EQ.2) THEN
672  fggf=fcof*ef**2*(1d0+2d0*rm1)*be34
673  fgzf=fcof*ef*vf*(1d0+2d0*rm1)*be34
674  fzzf=fcof*(vf**2*(1d0+2d0*rm1)+af**2*(1d0-4d0*rm1))*be34
675  ENDIF
676  IF(icase.EQ.1) wdtp(i)=fudge*wdtp(i)
677  IF(icase.EQ.1) wdtp(0)=wdtp(0)+wdtp(i)
678  IF(mdme(idc,1).GT.0) THEN
679  IF((icase.EQ.1.AND.mint(61).NE.1).OR.
680  & (icase.EQ.2.AND.mint(61).EQ.1)) THEN
681  wdte(i,mdme(idc,1))=wdtp(i)*wid2
682  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+
683  & wdte(i,mdme(idc,1))
684  wdte(i,0)=wdte(i,mdme(idc,1))
685  wdte(0,0)=wdte(0,0)+wdte(i,0)
686  ENDIF
687  IF(mint(61).EQ.2.AND.icase.EQ.2) THEN
688  IF(mstp(43).EQ.1.OR.mstp(43).EQ.3) vint(111)=
689  & vint(111)+fggf*wid2
690  IF(mstp(43).EQ.3) vint(112)=vint(112)+fgzf*wid2
691  IF(mstp(43).EQ.2.OR.mstp(43).EQ.3) vint(114)=
692  & vint(114)+fzzf*wid2
693  ENDIF
694  ENDIF
695  220 CONTINUE
696  IF(mint(61).GE.1) icase=3-icase
697  IF(icase.EQ.2) goto 210
698 
699  ELSEIF(kfla.EQ.24) THEN
700 C...W+/-:
701  fac=(aem/(24d0*xw))*shr
702  DO 230 i=1,mdcy(kc,3)
703  idc=i+mdcy(kc,2)-1
704  IF(mdme(idc,1).LT.0) goto 230
705  rm1=pmas(iabs(kfdp(idc,1)),1)**2/sh
706  rm2=pmas(iabs(kfdp(idc,2)),1)**2/sh
707  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 230
708  wid2=1d0
709  IF(i.LE.16) THEN
710 C...W+/- -> q + qbar'
711  fcof=3d0*radc*vckm((i-1)/4+1,mod(i-1,4)+1)
712  IF(kflr.GT.0) THEN
713  IF(mod(i,4).EQ.3) wid2=wids(6,2)
714  IF(mod(i,4).EQ.0) wid2=wids(8,2)
715  IF(i.GE.13) wid2=wid2*wids(7,3)
716  ELSE
717  IF(mod(i,4).EQ.3) wid2=wids(6,3)
718  IF(mod(i,4).EQ.0) wid2=wids(8,3)
719  IF(i.GE.13) wid2=wid2*wids(7,2)
720  ENDIF
721  ELSEIF(i.LE.20) THEN
722 C...W+/- -> l+/- + nu
723  fcof=1d0
724  IF(kflr.GT.0) THEN
725  IF(i.EQ.20) wid2=wids(17,3)*wids(18,2)
726  ELSE
727  IF(i.EQ.20) wid2=wids(17,2)*wids(18,3)
728  ENDIF
729  ENDIF
730  wdtp(i)=fac*fcof*(2d0-rm1-rm2-(rm1-rm2)**2)*
731  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))
732  wdtp(i)=fudge*wdtp(i)
733  wdtp(0)=wdtp(0)+wdtp(i)
734  IF(mdme(idc,1).GT.0) THEN
735  wdte(i,mdme(idc,1))=wdtp(i)*wid2
736  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
737  wdte(i,0)=wdte(i,mdme(idc,1))
738  wdte(0,0)=wdte(0,0)+wdte(i,0)
739  ENDIF
740  230 CONTINUE
741 
742  ELSEIF(kfla.EQ.25.OR.kfla.EQ.35.OR.kfla.EQ.36) THEN
743 C...h0 (or H0, or A0):
744  shfs=sh
745  fac=(aem/(8d0*xw))*(shfs/pmas(24,1)**2)*shr
746  DO 270 i=1,mdcy(kfhigg,3)
747  idc=i+mdcy(kfhigg,2)-1
748  IF(mdme(idc,1).LT.0) goto 270
749  kfc1=pycomp(kfdp(idc,1))
750  kfc2=pycomp(kfdp(idc,2))
751  rm1=pmas(kfc1,1)**2/sh
752  rm2=pmas(kfc2,1)**2/sh
753  IF(i.NE.16.AND.i.NE.17.AND.sqrt(rm1)+sqrt(rm2).GT.1d0)
754  & goto 270
755  wid2=1d0
756 
757  IF(i.LE.8) THEN
758 C...h0 -> q + qbar
759  wdtp(i)=fac*3d0*(pymrun(kfdp(idc,1),sh)**2/shfs)*
760  & sqrt(max(0d0,1d0-4d0*rm1))*radc
761 C...A0 behaves like beta, ho and H0 like beta**3.
762  IF(ihigg.NE.3) wdtp(i)=wdtp(i)*(1d0-4d0*rm1)
763  IF(mstp(4).GE.1.OR.ihigg.GE.2) THEN
764  IF(mod(i,2).EQ.1) wdtp(i)=wdtp(i)*paru(151+10*ihigg)**2
765  IF(mod(i,2).EQ.0) wdtp(i)=wdtp(i)*paru(152+10*ihigg)**2
766  IF(imss(1).NE.0.AND.kfc1.EQ.5) THEN
767  wdtp(i)=wdtp(i)/(1d0+rmss(41))**2
768  IF(ihigg.NE.3) THEN
769  wdtp(i)=wdtp(i)*(1d0+rmss(41)*paru(152+10*ihigg)/
770  & paru(151+10*ihigg))**2
771  ENDIF
772  ENDIF
773  ENDIF
774  IF(i.EQ.6) wid2=wids(6,1)
775  IF((i.EQ.7.OR.i.EQ.8)) wid2=wids(i,1)
776  ELSEIF(i.LE.12) THEN
777 C...h0 -> l+ + l-
778  wdtp(i)=fac*rm1*sqrt(max(0d0,1d0-4d0*rm1))*(sh/shfs)
779 C...A0 behaves like beta, ho and H0 like beta**3.
780  IF(ihigg.NE.3) wdtp(i)=wdtp(i)*(1d0-4d0*rm1)
781  IF(mstp(4).GE.1.OR.ihigg.GE.2) wdtp(i)=wdtp(i)*
782  & paru(153+10*ihigg)**2
783  IF(i.EQ.12) wid2=wids(17,1)
784 
785  ELSEIF(i.EQ.13) THEN
786 C...h0 -> g + g; quark loop contribution only
787  etare=0d0
788  etaim=0d0
789  DO 240 j=1,2*mstp(1)
790  eps=(2d0*pmas(j,1))**2/sh
791 C...Loop integral; function of eps=4m^2/shat; different for A0.
792  IF(eps.LE.1d0) THEN
793  IF(eps.GT.1d-4) THEN
794  root=sqrt(1d0-eps)
795  rln=log((1d0+root)/(1d0-root))
796  ELSE
797  rln=log(4d0/eps-2d0)
798  ENDIF
799  phire=-0.25d0*(rln**2-paru(1)**2)
800  phiim=0.5d0*paru(1)*rln
801  ELSE
802  phire=(asin(1d0/sqrt(eps)))**2
803  phiim=0d0
804  ENDIF
805  IF(ihigg.LE.2) THEN
806  etarej=-0.5d0*eps*(1d0+(1d0-eps)*phire)
807  etaimj=-0.5d0*eps*(1d0-eps)*phiim
808  ELSE
809  etarej=-0.5d0*eps*phire
810  etaimj=-0.5d0*eps*phiim
811  ENDIF
812 C...Couplings (=1 for standard model Higgs).
813  IF(mstp(4).GE.1.OR.ihigg.GE.2) THEN
814  IF(mod(j,2).EQ.1) THEN
815  etarej=etarej*paru(151+10*ihigg)
816  etaimj=etaimj*paru(151+10*ihigg)
817  ELSE
818  etarej=etarej*paru(152+10*ihigg)
819  etaimj=etaimj*paru(152+10*ihigg)
820  ENDIF
821  ENDIF
822  etare=etare+etarej
823  etaim=etaim+etaimj
824  240 CONTINUE
825  eta2=etare**2+etaim**2
826  wdtp(i)=fac*(as/paru(1))**2*eta2
827 
828  ELSEIF(i.EQ.14) THEN
829 C...h0 -> gamma + gamma; quark, lepton, W+- and H+- loop contributions
830  etare=0d0
831  etaim=0d0
832  jmax=3*mstp(1)+1
833  IF(mstp(4).GE.1.OR.ihigg.GE.2) jmax=jmax+1
834  DO 250 j=1,jmax
835  IF(j.LE.2*mstp(1)) THEN
836  ej=kchg(j,1)/3d0
837  eps=(2d0*pmas(j,1))**2/sh
838  ELSEIF(j.LE.3*mstp(1)) THEN
839  jl=2*(j-2*mstp(1))-1
840  ej=kchg(10+jl,1)/3d0
841  eps=(2d0*pmas(10+jl,1))**2/sh
842  ELSEIF(j.EQ.3*mstp(1)+1) THEN
843  eps=(2d0*pmas(24,1))**2/sh
844  ELSE
845  eps=(2d0*pmas(37,1))**2/sh
846  ENDIF
847 C...Loop integral; function of eps=4m^2/shat.
848  IF(eps.LE.1d0) THEN
849  IF(eps.GT.1d-4) THEN
850  root=sqrt(1d0-eps)
851  rln=log((1d0+root)/(1d0-root))
852  ELSE
853  rln=log(4d0/eps-2d0)
854  ENDIF
855  phire=-0.25d0*(rln**2-paru(1)**2)
856  phiim=0.5d0*paru(1)*rln
857  ELSE
858  phire=(asin(1d0/sqrt(eps)))**2
859  phiim=0d0
860  ENDIF
861  IF(j.LE.3*mstp(1)) THEN
862 C...Fermion loops: loop integral different for A0; charges.
863  IF(ihigg.LE.2) THEN
864  phipre=-0.5d0*eps*(1d0+(1d0-eps)*phire)
865  phipim=-0.5d0*eps*(1d0-eps)*phiim
866  ELSE
867  phipre=-0.5d0*eps*phire
868  phipim=-0.5d0*eps*phiim
869  ENDIF
870  IF(j.LE.2*mstp(1).AND.mod(j,2).EQ.1) THEN
871  ejc=3d0*ej**2
872  ejh=paru(151+10*ihigg)
873  ELSEIF(j.LE.2*mstp(1)) THEN
874  ejc=3d0*ej**2
875  ejh=paru(152+10*ihigg)
876  ELSE
877  ejc=ej**2
878  ejh=paru(153+10*ihigg)
879  ENDIF
880  IF(mstp(4).EQ.0.AND.ihigg.EQ.1) ejh=1d0
881  etarej=ejc*ejh*phipre
882  etaimj=ejc*ejh*phipim
883  ELSEIF(j.EQ.3*mstp(1)+1) THEN
884 C...W loops: loop integral and charges.
885  etarej=0.5d0+0.75d0*eps*(1d0+(2d0-eps)*phire)
886  etaimj=0.75d0*eps*(2d0-eps)*phiim
887  IF(mstp(4).GE.1.OR.ihigg.GE.2) THEN
888  etarej=etarej*paru(155+10*ihigg)
889  etaimj=etaimj*paru(155+10*ihigg)
890  ENDIF
891  ELSE
892 C...Charged H loops: loop integral and charges.
893  fachhh=(pmas(24,1)/pmas(37,1))**2*
894  & paru(158+10*ihigg+2*(ihigg/3))
895  etarej=eps*(1d0-eps*phire)*fachhh
896  etaimj=-eps**2*phiim*fachhh
897  ENDIF
898  etare=etare+etarej
899  etaim=etaim+etaimj
900  250 CONTINUE
901  eta2=etare**2+etaim**2
902  wdtp(i)=fac*(aem/paru(1))**2*0.5d0*eta2
903 
904  ELSEIF(i.EQ.15) THEN
905 C...h0 -> gamma + Z0; quark, lepton, W and H+- loop contributions
906  etare=0d0
907  etaim=0d0
908  jmax=3*mstp(1)+1
909  IF(mstp(4).GE.1.OR.ihigg.GE.2) jmax=jmax+1
910  DO 260 j=1,jmax
911  IF(j.LE.2*mstp(1)) THEN
912  ej=kchg(j,1)/3d0
913  aj=sign(1d0,ej+0.1d0)
914  vj=aj-4d0*ej*xwv
915  eps=(2d0*pmas(j,1))**2/sh
916  epsp=(2d0*pmas(j,1)/pmas(23,1))**2
917  ELSEIF(j.LE.3*mstp(1)) THEN
918  jl=2*(j-2*mstp(1))-1
919  ej=kchg(10+jl,1)/3d0
920  aj=sign(1d0,ej+0.1d0)
921  vj=aj-4d0*ej*xwv
922  eps=(2d0*pmas(10+jl,1))**2/sh
923  epsp=(2d0*pmas(10+jl,1)/pmas(23,1))**2
924  ELSE
925  eps=(2d0*pmas(24,1))**2/sh
926  epsp=(2d0*pmas(24,1)/pmas(23,1))**2
927  ENDIF
928 C...Loop integrals; functions of eps=4m^2/shat and eps'=4m^2/m_Z^2.
929  IF(eps.LE.1d0) THEN
930  root=sqrt(1d0-eps)
931  IF(eps.GT.1d-4) THEN
932  rln=log((1d0+root)/(1d0-root))
933  ELSE
934  rln=log(4d0/eps-2d0)
935  ENDIF
936  phire=-0.25d0*(rln**2-paru(1)**2)
937  phiim=0.5d0*paru(1)*rln
938  psire=0.5d0*root*rln
939  psiim=-0.5d0*root*paru(1)
940  ELSE
941  phire=(asin(1d0/sqrt(eps)))**2
942  phiim=0d0
943  psire=sqrt(eps-1d0)*asin(1d0/sqrt(eps))
944  psiim=0d0
945  ENDIF
946  IF(epsp.LE.1d0) THEN
947  root=sqrt(1d0-epsp)
948  IF(epsp.GT.1d-4) THEN
949  rln=log((1d0+root)/(1d0-root))
950  ELSE
951  rln=log(4d0/epsp-2d0)
952  ENDIF
953  phirep=-0.25d0*(rln**2-paru(1)**2)
954  phiimp=0.5d0*paru(1)*rln
955  psirep=0.5d0*root*rln
956  psiimp=-0.5d0*root*paru(1)
957  ELSE
958  phirep=(asin(1d0/sqrt(epsp)))**2
959  phiimp=0d0
960  psirep=sqrt(epsp-1d0)*asin(1d0/sqrt(epsp))
961  psiimp=0d0
962  ENDIF
963  fxyre=eps*epsp/(8d0*(eps-epsp))*(1d0+eps*epsp/(eps-epsp)*
964  & (phire-phirep)+2d0*eps/(eps-epsp)*(psire-psirep))
965  fxyim=eps**2*epsp/(8d0*(eps-epsp)**2)*
966  & (epsp*(phiim-phiimp)+2d0*(psiim-psiimp))
967  f1re=-eps*epsp/(2d0*(eps-epsp))*(phire-phirep)
968  f1im=-eps*epsp/(2d0*(eps-epsp))*(phiim-phiimp)
969  IF(j.LE.3*mstp(1)) THEN
970 C...Fermion loops: loop integral different for A0; charges.
971  IF(ihigg.EQ.3) fxyre=0d0
972  IF(ihigg.EQ.3) fxyim=0d0
973  IF(j.LE.2*mstp(1).AND.mod(j,2).EQ.1) THEN
974  ejc=-3d0*ej*vj
975  ejh=paru(151+10*ihigg)
976  ELSEIF(j.LE.2*mstp(1)) THEN
977  ejc=-3d0*ej*vj
978  ejh=paru(152+10*ihigg)
979  ELSE
980  ejc=-ej*vj
981  ejh=paru(153+10*ihigg)
982  ENDIF
983  IF(mstp(4).EQ.0.AND.ihigg.EQ.1) ejh=1d0
984  etarej=ejc*ejh*(fxyre-0.25d0*f1re)
985  etaimj=ejc*ejh*(fxyim-0.25d0*f1im)
986  ELSEIF(j.EQ.3*mstp(1)+1) THEN
987 C...W loops: loop integral and charges.
988  heps=(1d0+2d0/eps)*xw/xw1-(5d0+2d0/eps)
989  etarej=-xw1*((3d0-xw/xw1)*f1re+heps*fxyre)
990  etaimj=-xw1*((3d0-xw/xw1)*f1im+heps*fxyim)
991  IF(mstp(4).GE.1.OR.ihigg.GE.2) THEN
992  etarej=etarej*paru(155+10*ihigg)
993  etaimj=etaimj*paru(155+10*ihigg)
994  ENDIF
995  ELSE
996 C...Charged H loops: loop integral and charges.
997  fachhh=(pmas(24,1)/pmas(37,1))**2*(1d0-2d0*xw)*
998  & paru(158+10*ihigg+2*(ihigg/3))
999  etarej=fachhh*fxyre
1000  etaimj=fachhh*fxyim
1001  ENDIF
1002  etare=etare+etarej
1003  etaim=etaim+etaimj
1004  260 CONTINUE
1005  eta2=(etare**2+etaim**2)/(xw*xw1)
1006  wdtp(i)=fac*(aem/paru(1))**2*(1d0-pmas(23,1)**2/sh)**3*eta2
1007  wid2=wids(23,2)
1008 
1009  ELSEIF(i.LE.17) THEN
1010 C...h0 -> Z0 + Z0, W+ + W-
1011  pm1=pmas(iabs(kfdp(idc,1)),1)
1012  pg1=pmas(iabs(kfdp(idc,1)),2)
1013  IF(mint(62).GE.1) THEN
1014  IF(mstp(42).EQ.0.OR.(4d0*(pm1+10d0*pg1)**2.LT.sh.AND.
1015  & ckin(46).LT.ckin(45).AND.ckin(48).LT.ckin(47).AND.
1016  & max(ckin(45),ckin(47)).LT.pm1-10d0*pg1)) THEN
1017  mofsv(ihigg,i-15)=0
1018  widw=(1d0-4d0*rm1+12d0*rm1**2)*sqrt(max(0d0,
1019  & 1d0-4d0*rm1))
1020  wid2=1d0
1021  ELSE
1022  mofsv(ihigg,i-15)=1
1023  rmas=sqrt(max(0d0,sh))
1024  CALL pyofsh(1,kfla,kfdp(idc,1),kfdp(idc,2),rmas,widw,
1025  & wid2)
1026  widwsv(ihigg,i-15)=widw
1027  wid2sv(ihigg,i-15)=wid2
1028  ENDIF
1029  ELSE
1030  IF(mofsv(ihigg,i-15).EQ.0) THEN
1031  widw=(1d0-4d0*rm1+12d0*rm1**2)*sqrt(max(0d0,
1032  & 1d0-4d0*rm1))
1033  wid2=1d0
1034  ELSE
1035  widw=widwsv(ihigg,i-15)
1036  wid2=wid2sv(ihigg,i-15)
1037  ENDIF
1038  ENDIF
1039  wdtp(i)=fac*widw/(2d0*(18-i))
1040  IF(mstp(49).NE.0) wdtp(i)=wdtp(i)*pmas(kfhigg,1)**2/shfs
1041  IF(mstp(4).GE.1.OR.ihigg.GE.2) wdtp(i)=wdtp(i)*
1042  & paru(138+i+10*ihigg)**2
1043  wid2=wid2*wids(7+i,1)
1044 
1045  ELSEIF(i.EQ.18.AND.ihigg.GE.2) THEN
1046 C...H0 -> Z0 + h0, A0-> Z0 + h0
1047  wdtp(i)=fac*0.5d0*sqrt(max(0d0,
1048  & (1d0-rm1-rm2)**2-4d0*rm1*rm2))**3
1049  IF(ihigg.EQ.2) THEN
1050  wdtp(i)=wdtp(i)*paru(179)**2
1051  ELSEIF(ihigg.EQ.3) THEN
1052  wdtp(i)=wdtp(i)*paru(186)**2
1053  ENDIF
1054  wid2=wids(23,2)*wids(25,2)
1055 
1056  ELSEIF(i.EQ.19.AND.ihigg.GE.2) THEN
1057 C...H0 -> h0 + h0, A0-> h0 + h0
1058  wdtp(i)=fac*0.25d0*
1059  & pmas(23,1)**4/sh**2*sqrt(max(0d0,1d0-4d0*rm1))
1060  IF(ihigg.EQ.2) THEN
1061  wdtp(i)=wdtp(i)*paru(176)**2
1062  ELSEIF(ihigg.EQ.3) THEN
1063  wdtp(i)=wdtp(i)*paru(169)**2
1064  ENDIF
1065  wid2=wids(25,1)
1066  ELSEIF((i.EQ.20.OR.i.EQ.21).AND.ihigg.GE.2) THEN
1067 C...H0 -> W+/- + H-/+, A0 -> W+/- + H-/+
1068  wdtp(i)=fac*0.5d0*sqrt(max(0d0,
1069  & (1d0-rm1-rm2)**2-4d0*rm1*rm2))**3
1070  & *paru(195+ihigg)**2
1071  IF(i.EQ.20) THEN
1072  wid2=wids(24,2)*wids(37,3)
1073  ELSEIF(i.EQ.21) THEN
1074  wid2=wids(24,3)*wids(37,2)
1075  ENDIF
1076 
1077  ELSEIF(i.EQ.22.AND.ihigg.EQ.2) THEN
1078 C...H0 -> Z0 + A0.
1079  wdtp(i)=fac*0.5d0*paru(187)**2*sqrt(max(0d0,
1080  & (1d0-rm1-rm2)**2-4d0*rm1*rm2))**3
1081  wid2=wids(36,2)*wids(23,2)
1082 
1083  ELSEIF(i.EQ.23.AND.ihigg.EQ.2) THEN
1084 C...H0 -> h0 + A0.
1085  wdtp(i)=fac*0.5d0*paru(180)**2*
1086  & pmas(23,1)**4/sh**2*sqrt(max(0d0,1d0-4d0*rm1))
1087  wid2=wids(25,2)*wids(36,2)
1088 
1089  ELSEIF(i.EQ.24.AND.ihigg.EQ.2) THEN
1090 C...H0 -> A0 + A0
1091  wdtp(i)=fac*0.25d0*paru(177)**2*
1092  & pmas(23,1)**4/sh**2*sqrt(max(0d0,1d0-4d0*rm1))
1093  wid2=wids(36,1)
1094 
1095 CMRENNA++
1096  ELSE
1097 C...Add in SUSY decays (two-body) by rescaling by phase space factor.
1098  rm10=rm1*sh/pmr**2
1099  rm20=rm2*sh/pmr**2
1100  wfac0=1d0+rm10**2+rm20**2-2d0*(rm10+rm20+rm10*rm20)
1101  wfac=1d0+rm1**2+rm2**2-2d0*(rm1+rm2+rm1*rm2)
1102  IF(wfac.LE.0d0 .OR. wfac0.LE.0d0) THEN
1103  wfac=0d0
1104  ELSE
1105  wfac=wfac/wfac0
1106  ENDIF
1107  wdtp(i)=pmas(kfla,2)*brat(idc)*(shr/pmr)*sqrt(wfac)
1108 CMRENNA--
1109  IF(kfc2.EQ.kfc1) THEN
1110  wid2=wids(kfc1,1)
1111  ELSE
1112  ksgn1=2
1113  IF(kfdp(idc,1).LT.0) ksgn1=3
1114  ksgn2=2
1115  IF(kfdp(idc,2).LT.0) ksgn2=3
1116  wid2=wids(kfc1,ksgn1)*wids(kfc2,ksgn2)
1117  ENDIF
1118  ENDIF
1119  wdtp(i)=fudge*wdtp(i)
1120  wdtp(0)=wdtp(0)+wdtp(i)
1121  IF(mdme(idc,1).GT.0) THEN
1122  wdte(i,mdme(idc,1))=wdtp(i)*wid2
1123  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
1124  wdte(i,0)=wdte(i,mdme(idc,1))
1125  wdte(0,0)=wdte(0,0)+wdte(i,0)
1126  ENDIF
1127  270 CONTINUE
1128 
1129  ELSEIF(kfla.EQ.32) THEN
1130 C...Z'0:
1131  icase=1
1132  xwc=1d0/(16d0*xw*xw1)
1133  fac=(aem*xwc/3d0)*shr
1134  vint(117)=0d0
1135  280 CONTINUE
1136  IF(mint(61).GE.1.AND.icase.EQ.2) THEN
1137  vint(111)=0d0
1138  vint(112)=0d0
1139  vint(113)=0d0
1140  vint(114)=0d0
1141  vint(115)=0d0
1142  vint(116)=0d0
1143  ENDIF
1144  IF(mint(61).EQ.1.AND.icase.EQ.2) THEN
1145  kfai=iabs(mint(15))
1146  ei=kchg(kfai,1)/3d0
1147  ai=sign(1d0,ei+0.1d0)
1148  vi=ai-4d0*ei*xwv
1149  kfaic=1
1150  IF(kfai.LE.10.AND.mod(kfai,2).EQ.0) kfaic=2
1151  IF(kfai.GT.10.AND.mod(kfai,2).NE.0) kfaic=3
1152  IF(kfai.GT.10.AND.mod(kfai,2).EQ.0) kfaic=4
1153  IF(kfai.LE.2.OR.kfai.EQ.11.OR.kfai.EQ.12) THEN
1154  vpi=paru(119+2*kfaic)
1155  api=paru(120+2*kfaic)
1156  ELSEIF(kfai.LE.4.OR.kfai.EQ.13.OR.kfai.EQ.14) THEN
1157  vpi=parj(178+2*kfaic)
1158  api=parj(179+2*kfaic)
1159  ELSE
1160  vpi=parj(186+2*kfaic)
1161  api=parj(187+2*kfaic)
1162  ENDIF
1163  sqmz=pmas(23,1)**2
1164  hz=shr*vint(117)
1165  sqmzp=pmas(32,1)**2
1166  hzp=shr*wdtp(0)
1167  IF(mstp(44).EQ.1.OR.mstp(44).EQ.4.OR.mstp(44).EQ.5.OR.
1168  & mstp(44).EQ.7) vint(111)=1d0
1169  IF(mstp(44).EQ.4.OR.mstp(44).EQ.7) vint(112)=
1170  & 2d0*xwc*sh*(sh-sqmz)/((sh-sqmz)**2+hz**2)
1171  IF(mstp(44).EQ.5.OR.mstp(44).EQ.7) vint(113)=
1172  & 2d0*xwc*sh*(sh-sqmzp)/((sh-sqmzp)**2+hzp**2)
1173  IF(mstp(44).EQ.2.OR.mstp(44).EQ.4.OR.mstp(44).EQ.6.OR.
1174  & mstp(44).EQ.7) vint(114)=xwc**2*sh**2/((sh-sqmz)**2+hz**2)
1175  IF(mstp(44).EQ.6.OR.mstp(44).EQ.7) vint(115)=
1176  & 2d0*xwc**2*sh**2*((sh-sqmz)*(sh-sqmzp)+hz*hzp)/
1177  & (((sh-sqmz)**2+hz**2)*((sh-sqmzp)**2+hzp**2))
1178  IF(mstp(44).EQ.3.OR.mstp(44).EQ.5.OR.mstp(44).EQ.6.OR.
1179  & mstp(44).EQ.7) vint(116)=xwc**2*sh**2/((sh-sqmzp)**2+hzp**2)
1180  ENDIF
1181  DO 290 i=1,mdcy(kc,3)
1182  idc=i+mdcy(kc,2)-1
1183  IF(mdme(idc,1).LT.0) goto 290
1184  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
1185  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
1186  IF(sqrt(rm1)+sqrt(rm2).GT.1d0.OR.mdme(idc,1).LT.0) goto 290
1187  wid2=1d0
1188  IF(i.LE.16) THEN
1189  IF(i.LE.8) THEN
1190 C...Z'0 -> q + qbar
1191  ef=kchg(i,1)/3d0
1192  af=sign(1d0,ef+0.1d0)
1193  vf=af-4d0*ef*xwv
1194  IF(i.LE.2) THEN
1195  vpf=paru(123-2*mod(i,2))
1196  apf=paru(124-2*mod(i,2))
1197  ELSEIF(i.LE.4) THEN
1198  vpf=parj(182-2*mod(i,2))
1199  apf=parj(183-2*mod(i,2))
1200  ELSE
1201  vpf=parj(190-2*mod(i,2))
1202  apf=parj(191-2*mod(i,2))
1203  ENDIF
1204  fcof=3d0*radc
1205  IF(i.GE.6.AND.mstp(35).GE.1) fcof=fcof*
1206  & pyhfth(sh,sh*rm1,1d0)
1207  IF(i.EQ.6) wid2=wids(6,1)
1208  IF((i.EQ.7.OR.i.EQ.8)) wid2=wids(i,1)
1209  ELSEIF(i.LE.16) THEN
1210 C...Z'0 -> l+ + l-, nu + nubar
1211  ef=kchg(i+2,1)/3d0
1212  af=sign(1d0,ef+0.1d0)
1213  vf=af-4d0*ef*xwv
1214  IF(i.LE.10) THEN
1215  vpf=paru(127-2*mod(i,2))
1216  apf=paru(128-2*mod(i,2))
1217  ELSEIF(i.LE.12) THEN
1218  vpf=parj(186-2*mod(i,2))
1219  apf=parj(187-2*mod(i,2))
1220  ELSE
1221  vpf=parj(194-2*mod(i,2))
1222  apf=parj(195-2*mod(i,2))
1223  ENDIF
1224  fcof=1d0
1225  IF((i.EQ.15.OR.i.EQ.16)) wid2=wids(2+i,1)
1226  ENDIF
1227  be34=sqrt(max(0d0,1d0-4d0*rm1))
1228  IF(icase.EQ.1) THEN
1229  wdtpz=fcof*(vf**2*(1d0+2d0*rm1)+af**2*(1d0-4d0*rm1))*be34
1230  wdtp(i)=fac*fcof*(vpf**2*(1d0+2d0*rm1)+
1231  & apf**2*(1d0-4d0*rm1))*be34
1232  ELSEIF(mint(61).EQ.1.AND.icase.EQ.2) THEN
1233  wdtp(i)=fac*fcof*((ei**2*vint(111)*ef**2+ei*vi*vint(112)*
1234  & ef*vf+ei*vpi*vint(113)*ef*vpf+(vi**2+ai**2)*vint(114)*
1235  & vf**2+(vi*vpi+ai*api)*vint(115)*vf*vpf+(vpi**2+api**2)*
1236  & vint(116)*vpf**2)*(1d0+2d0*rm1)+((vi**2+ai**2)*vint(114)*
1237  & af**2+(vi*vpi+ai*api)*vint(115)*af*apf+(vpi**2+api**2)*
1238  & vint(116)*apf**2)*(1d0-4d0*rm1))*be34
1239  ELSEIF(mint(61).EQ.2) THEN
1240  fggf=fcof*ef**2*(1d0+2d0*rm1)*be34
1241  fgzf=fcof*ef*vf*(1d0+2d0*rm1)*be34
1242  fgzpf=fcof*ef*vpf*(1d0+2d0*rm1)*be34
1243  fzzf=fcof*(vf**2*(1d0+2d0*rm1)+af**2*(1d0-4d0*rm1))*be34
1244  fzzpf=fcof*(vf*vpf*(1d0+2d0*rm1)+af*apf*(1d0-4d0*rm1))*
1245  & be34
1246  fzpzpf=fcof*(vpf**2*(1d0+2d0*rm1)+apf**2*(1d0-4d0*rm1))*
1247  & be34
1248  ENDIF
1249  ELSEIF(i.EQ.17) THEN
1250 C...Z'0 -> W+ + W-
1251  wdtpzp=paru(129)**2*xw1**2*
1252  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3*
1253  & (1d0+10d0*rm1+10d0*rm2+rm1**2+rm2**2+10d0*rm1*rm2)
1254  IF(icase.EQ.1) THEN
1255  wdtpz=0d0
1256  wdtp(i)=fac*wdtpzp
1257  ELSEIF(mint(61).EQ.1.AND.icase.EQ.2) THEN
1258  wdtp(i)=fac*(vpi**2+api**2)*vint(116)*wdtpzp
1259  ELSEIF(mint(61).EQ.2) THEN
1260  fggf=0d0
1261  fgzf=0d0
1262  fgzpf=0d0
1263  fzzf=0d0
1264  fzzpf=0d0
1265  fzpzpf=wdtpzp
1266  ENDIF
1267  wid2=wids(24,1)
1268  ELSEIF(i.EQ.18) THEN
1269 C...Z'0 -> H+ + H-
1270  czc=2d0*(1d0-2d0*xw)
1271  be34c=(1d0-4d0*rm1)*sqrt(max(0d0,1d0-4d0*rm1))
1272  IF(icase.EQ.1) THEN
1273  wdtpz=0.25d0*paru(142)**2*czc**2*be34c
1274  wdtp(i)=fac*0.25d0*paru(143)**2*czc**2*be34c
1275  ELSEIF(mint(61).EQ.1.AND.icase.EQ.2) THEN
1276  wdtp(i)=fac*0.25d0*(ei**2*vint(111)+paru(142)*ei*vi*
1277  & vint(112)*czc+paru(143)*ei*vpi*vint(113)*czc+paru(142)**2*
1278  & (vi**2+ai**2)*vint(114)*czc**2+paru(142)*paru(143)*
1279  & (vi*vpi+ai*api)*vint(115)*czc**2+paru(143)**2*
1280  & (vpi**2+api**2)*vint(116)*czc**2)*be34c
1281  ELSEIF(mint(61).EQ.2) THEN
1282  fggf=0.25d0*be34c
1283  fgzf=0.25d0*paru(142)*czc*be34c
1284  fgzpf=0.25d0*paru(143)*czc*be34c
1285  fzzf=0.25d0*paru(142)**2*czc**2*be34c
1286  fzzpf=0.25d0*paru(142)*paru(143)*czc**2*be34c
1287  fzpzpf=0.25d0*paru(143)**2*czc**2*be34c
1288  ENDIF
1289  wid2=wids(37,1)
1290  ELSEIF(i.EQ.19) THEN
1291 C...Z'0 -> Z0 + gamma.
1292  ELSEIF(i.EQ.20) THEN
1293 C...Z'0 -> Z0 + h0
1294  flam=sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))
1295  wdtpzp=paru(145)**2*4d0*abs(1d0-2d0*xw)*
1296  & (3d0*rm1+0.25d0*flam**2)*flam
1297  IF(icase.EQ.1) THEN
1298  wdtpz=0d0
1299  wdtp(i)=fac*wdtpzp
1300  ELSEIF(mint(61).EQ.1.AND.icase.EQ.2) THEN
1301  wdtp(i)=fac*(vpi**2+api**2)*vint(116)*wdtpzp
1302  ELSEIF(mint(61).EQ.2) THEN
1303  fggf=0d0
1304  fgzf=0d0
1305  fgzpf=0d0
1306  fzzf=0d0
1307  fzzpf=0d0
1308  fzpzpf=wdtpzp
1309  ENDIF
1310  wid2=wids(23,2)*wids(25,2)
1311  ELSEIF(i.EQ.21.OR.i.EQ.22) THEN
1312 C...Z' -> h0 + A0 or H0 + A0.
1313  be34c=sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3
1314  IF(i.EQ.21) THEN
1315  czah=paru(186)
1316  czpah=paru(188)
1317  ELSE
1318  czah=paru(187)
1319  czpah=paru(189)
1320  ENDIF
1321  IF(icase.EQ.1) THEN
1322  wdtpz=czah**2*be34c
1323  wdtp(i)=fac*czpah**2*be34c
1324  ELSEIF(mint(61).EQ.1.AND.icase.EQ.2) THEN
1325  wdtp(i)=fac*(czah**2*(vi**2+ai**2)*vint(114)+czah*czpah*
1326  & (vi*vpi+ai*api)*vint(115)+czpah**2*(vpi**2+api**2)*
1327  & vint(116))*be34c
1328  ELSEIF(mint(61).EQ.2) THEN
1329  fggf=0d0
1330  fgzf=0d0
1331  fgzpf=0d0
1332  fzzf=czah**2*be34c
1333  fzzpf=czah*czpah*be34c
1334  fzpzpf=czpah**2*be34c
1335  ENDIF
1336  IF(i.EQ.21) wid2=wids(25,2)*wids(36,2)
1337  IF(i.EQ.22) wid2=wids(35,2)*wids(36,2)
1338  ENDIF
1339  IF(icase.EQ.1) THEN
1340  vint(117)=vint(117)+fac*wdtpz
1341  wdtp(i)=fudge*wdtp(i)
1342  wdtp(0)=wdtp(0)+wdtp(i)
1343  ENDIF
1344  IF(mdme(idc,1).GT.0) THEN
1345  IF((icase.EQ.1.AND.mint(61).NE.1).OR.
1346  & (icase.EQ.2.AND.mint(61).EQ.1)) THEN
1347  wdte(i,mdme(idc,1))=wdtp(i)*wid2
1348  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+
1349  & wdte(i,mdme(idc,1))
1350  wdte(i,0)=wdte(i,mdme(idc,1))
1351  wdte(0,0)=wdte(0,0)+wdte(i,0)
1352  ENDIF
1353  IF(mint(61).EQ.2.AND.icase.EQ.2) THEN
1354  IF(mstp(44).EQ.1.OR.mstp(44).EQ.4.OR.mstp(44).EQ.5.OR.
1355  & mstp(44).EQ.7) vint(111)=vint(111)+fggf*wid2
1356  IF(mstp(44).EQ.4.OR.mstp(44).EQ.7) vint(112)=vint(112)+
1357  & fgzf*wid2
1358  IF(mstp(44).EQ.5.OR.mstp(44).EQ.7) vint(113)=vint(113)+
1359  & fgzpf*wid2
1360  IF(mstp(44).EQ.2.OR.mstp(44).EQ.4.OR.mstp(44).EQ.6.OR.
1361  & mstp(44).EQ.7) vint(114)=vint(114)+fzzf*wid2
1362  IF(mstp(44).EQ.6.OR.mstp(44).EQ.7) vint(115)=vint(115)+
1363  & fzzpf*wid2
1364  IF(mstp(44).EQ.3.OR.mstp(44).EQ.5.OR.mstp(44).EQ.6.OR.
1365  & mstp(44).EQ.7) vint(116)=vint(116)+fzpzpf*wid2
1366  ENDIF
1367  ENDIF
1368  290 CONTINUE
1369  IF(mint(61).GE.1) icase=3-icase
1370  IF(icase.EQ.2) goto 280
1371 
1372  ELSEIF(kfla.EQ.34) THEN
1373 C...W'+/-:
1374  fac=(aem/(24d0*xw))*shr
1375  DO 300 i=1,mdcy(kc,3)
1376  idc=i+mdcy(kc,2)-1
1377  IF(mdme(idc,1).LT.0) goto 300
1378  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
1379  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
1380  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 300
1381  wid2=1d0
1382  IF(i.LE.20) THEN
1383  IF(i.LE.16) THEN
1384 C...W'+/- -> q + qbar'
1385  fcof=3d0*radc*(paru(131)**2+paru(132)**2)*
1386  & vckm((i-1)/4+1,mod(i-1,4)+1)
1387  IF(kflr.GT.0) THEN
1388  IF(mod(i,4).EQ.3) wid2=wids(6,2)
1389  IF(mod(i,4).EQ.0) wid2=wids(8,2)
1390  IF(i.GE.13) wid2=wid2*wids(7,3)
1391  ELSE
1392  IF(mod(i,4).EQ.3) wid2=wids(6,3)
1393  IF(mod(i,4).EQ.0) wid2=wids(8,3)
1394  IF(i.GE.13) wid2=wid2*wids(7,2)
1395  ENDIF
1396  ELSEIF(i.LE.20) THEN
1397 C...W'+/- -> l+/- + nu
1398  fcof=paru(133)**2+paru(134)**2
1399  IF(kflr.GT.0) THEN
1400  IF(i.EQ.20) wid2=wids(17,3)*wids(18,2)
1401  ELSE
1402  IF(i.EQ.20) wid2=wids(17,2)*wids(18,3)
1403  ENDIF
1404  ENDIF
1405  wdtp(i)=fac*fcof*0.5d0*(2d0-rm1-rm2-(rm1-rm2)**2)*
1406  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))
1407  ELSEIF(i.EQ.21) THEN
1408 C...W'+/- -> W+/- + Z0
1409  wdtp(i)=fac*paru(135)**2*0.5d0*xw1*(rm1/rm2)*
1410  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3*
1411  & (1d0+10d0*rm1+10d0*rm2+rm1**2+rm2**2+10d0*rm1*rm2)
1412  IF(kflr.GT.0) wid2=wids(24,2)*wids(23,2)
1413  IF(kflr.LT.0) wid2=wids(24,3)*wids(23,2)
1414  ELSEIF(i.EQ.23) THEN
1415 C...W'+/- -> W+/- + h0
1416  flam=sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))
1417  wdtp(i)=fac*paru(146)**2*2d0*(3d0*rm1+0.25d0*flam**2)*flam
1418  IF(kflr.GT.0) wid2=wids(24,2)*wids(25,2)
1419  IF(kflr.LT.0) wid2=wids(24,3)*wids(25,2)
1420  ENDIF
1421  wdtp(i)=fudge*wdtp(i)
1422  wdtp(0)=wdtp(0)+wdtp(i)
1423  IF(mdme(idc,1).GT.0) THEN
1424  wdte(i,mdme(idc,1))=wdtp(i)*wid2
1425  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
1426  wdte(i,0)=wdte(i,mdme(idc,1))
1427  wdte(0,0)=wdte(0,0)+wdte(i,0)
1428  ENDIF
1429  300 CONTINUE
1430 
1431  ELSEIF(kfla.EQ.37) THEN
1432 C...H+/-:
1433 C IF(MSTP(49).EQ.0) THEN
1434  shfs=sh
1435 C ELSE
1436 C SHFS=PMAS(37,1)**2
1437 C ENDIF
1438  fac=(aem/(8d0*xw))*(shfs/pmas(24,1)**2)*shr
1439  DO 310 i=1,mdcy(kc,3)
1440  idc=i+mdcy(kc,2)-1
1441  IF(mdme(idc,1).LT.0) goto 310
1442  kfc1=pycomp(kfdp(idc,1))
1443  kfc2=pycomp(kfdp(idc,2))
1444  rm1=pmas(kfc1,1)**2/sh
1445  rm2=pmas(kfc2,1)**2/sh
1446  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 310
1447  wid2=1d0
1448  IF(i.LE.4) THEN
1449 C...H+/- -> q + qbar'
1450  rm1r=pymrun(kfdp(idc,1),sh)**2/sh
1451  rm2r=pymrun(kfdp(idc,2),sh)**2/sh
1452  wdtp(i)=fac*3d0*radc*max(0d0,(rm1r*paru(141)**2+
1453  & rm2r/paru(141)**2)*(1d0-rm1r-rm2r)-4d0*rm1r*rm2r)*
1454  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))*(sh/shfs)
1455  IF(kflr.GT.0) THEN
1456  IF(i.EQ.3) wid2=wids(6,2)
1457  IF(i.EQ.4) wid2=wids(7,3)*wids(8,2)
1458  ELSE
1459  IF(i.EQ.3) wid2=wids(6,3)
1460  IF(i.EQ.4) wid2=wids(7,2)*wids(8,3)
1461  ENDIF
1462  ELSEIF(i.LE.8) THEN
1463 C...H+/- -> l+/- + nu
1464  wdtp(i)=fac*((rm1*paru(141)**2+rm2/paru(141)**2)*
1465  & (1d0-rm1-rm2)-4d0*rm1*rm2)*sqrt(max(0d0,
1466  & (1d0-rm1-rm2)**2-4d0*rm1*rm2))*(sh/shfs)
1467  IF(kflr.GT.0) THEN
1468  IF(i.EQ.8) wid2=wids(17,3)*wids(18,2)
1469  ELSE
1470  IF(i.EQ.8) wid2=wids(17,2)*wids(18,3)
1471  ENDIF
1472  ELSEIF(i.EQ.9) THEN
1473 C...H+/- -> W+/- + h0.
1474  wdtp(i)=fac*paru(195)**2*0.5d0*sqrt(max(0d0,
1475  & (1d0-rm1-rm2)**2-4d0*rm1*rm2))**3
1476  IF(kflr.GT.0) wid2=wids(24,2)*wids(25,2)
1477  IF(kflr.LT.0) wid2=wids(24,3)*wids(25,2)
1478 
1479 CMRENNA++
1480  ELSE
1481 C...Add in SUSY decays (two-body) by rescaling by phase space factor.
1482  rm10=rm1*sh/pmr**2
1483  rm20=rm2*sh/pmr**2
1484  wfac0=1d0+rm10**2+rm20**2-2d0*(rm10+rm20+rm10*rm20)
1485  wfac=1d0+rm1**2+rm2**2-2d0*(rm1+rm2+rm1*rm2)
1486  IF(wfac.LE.0d0 .OR. wfac0.LE.0d0) THEN
1487  wfac=0d0
1488  ELSE
1489  wfac=wfac/wfac0
1490  ENDIF
1491  wdtp(i)=pmas(kc,2)*brat(idc)*(shr/pmr)*sqrt(wfac)
1492 CMRENNA--
1493  ksgn1=2
1494  IF(kfls*kfdp(idc,1).LT.0.AND.kchg(kfc1,3).EQ.1) ksgn1=3
1495  ksgn2=2
1496  IF(kfls*kfdp(idc,2).LT.0.AND.kchg(kfc2,3).EQ.1) ksgn2=3
1497  wid2=wids(kfc1,ksgn1)*wids(kfc2,ksgn2)
1498  ENDIF
1499  wdtp(i)=fudge*wdtp(i)
1500  wdtp(0)=wdtp(0)+wdtp(i)
1501  IF(mdme(idc,1).GT.0) THEN
1502  wdte(i,mdme(idc,1))=wdtp(i)*wid2
1503  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
1504  wdte(i,0)=wdte(i,mdme(idc,1))
1505  wdte(0,0)=wdte(0,0)+wdte(i,0)
1506  ENDIF
1507  310 CONTINUE
1508 
1509  ELSEIF(kfla.EQ.41) THEN
1510 C...R:
1511  fac=(aem/(12d0*xw))*shr
1512  DO 320 i=1,mdcy(kc,3)
1513  idc=i+mdcy(kc,2)-1
1514  IF(mdme(idc,1).LT.0) goto 320
1515  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
1516  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
1517  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 320
1518  wid2=1d0
1519  IF(i.LE.6) THEN
1520 C...R -> q + qbar'
1521  fcof=3d0*radc
1522  ELSEIF(i.LE.9) THEN
1523 C...R -> l+ + l'-
1524  fcof=1d0
1525  ENDIF
1526  wdtp(i)=fac*fcof*(2d0-rm1-rm2-(rm1-rm2)**2)*
1527  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))
1528  IF(kflr.GT.0) THEN
1529  IF(i.EQ.4) wid2=wids(6,3)
1530  IF(i.EQ.5) wid2=wids(7,3)
1531  IF(i.EQ.6) wid2=wids(6,2)*wids(8,3)
1532  IF(i.EQ.9) wid2=wids(17,3)
1533  ELSE
1534  IF(i.EQ.4) wid2=wids(6,2)
1535  IF(i.EQ.5) wid2=wids(7,2)
1536  IF(i.EQ.6) wid2=wids(6,3)*wids(8,2)
1537  IF(i.EQ.9) wid2=wids(17,2)
1538  ENDIF
1539  wdtp(i)=fudge*wdtp(i)
1540  wdtp(0)=wdtp(0)+wdtp(i)
1541  IF(mdme(idc,1).GT.0) THEN
1542  wdte(i,mdme(idc,1))=wdtp(i)*wid2
1543  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
1544  wdte(i,0)=wdte(i,mdme(idc,1))
1545  wdte(0,0)=wdte(0,0)+wdte(i,0)
1546  ENDIF
1547  320 CONTINUE
1548 
1549  ELSEIF(kfla.EQ.42) THEN
1550 C...LQ (leptoquark).
1551  fac=(aem/4d0)*paru(151)*shr
1552  DO 330 i=1,mdcy(kc,3)
1553  idc=i+mdcy(kc,2)-1
1554  IF(mdme(idc,1).LT.0) goto 330
1555  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
1556  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
1557  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 330
1558  wdtp(i)=fac*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3
1559  wid2=1d0
1560  ilqq=kfdp(idc,1)*isign(1,kflr)
1561  IF(ilqq.GE.6) wid2=wids(ilqq,2)
1562  IF(ilqq.LE.-6) wid2=wids(-ilqq,3)
1563  ilql=kfdp(idc,2)*isign(1,kflr)
1564  IF(ilql.GE.17) wid2=wid2*wids(ilql,2)
1565  IF(ilql.LE.-17) wid2=wid2*wids(-ilql,3)
1566  wdtp(i)=fudge*wdtp(i)
1567  wdtp(0)=wdtp(0)+wdtp(i)
1568  IF(mdme(idc,1).GT.0) THEN
1569  wdte(i,mdme(idc,1))=wdtp(i)*wid2
1570  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
1571  wdte(i,0)=wdte(i,mdme(idc,1))
1572  wdte(0,0)=wdte(0,0)+wdte(i,0)
1573  ENDIF
1574  330 CONTINUE
1575 
1576  ELSEIF(kfla.EQ.ktechn+111.OR.kfla.EQ.ktechn+221) THEN
1577 C...Techni-pi0 and techni-pi0':
1578  fac=(1d0/(32d0*paru(1)*rtcm(1)**2))*shr
1579  DO 340 i=1,mdcy(kc,3)
1580  idc=i+mdcy(kc,2)-1
1581  IF(mdme(idc,1).LT.0) goto 340
1582  pm1=pmas(pycomp(kfdp(idc,1)),1)
1583  pm2=pmas(pycomp(kfdp(idc,2)),1)
1584  rm1=pm1**2/sh
1585  rm2=pm2**2/sh
1586  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 340
1587  wid2=1d0
1588 C...pi_tc -> g + g
1589  IF(i.EQ.8) THEN
1590  facp=(as/(4d0*paru(1))*itcm(1)/rtcm(1))**2
1591  & /(8d0*paru(1))*sh*shr
1592  IF(kfla.EQ.ktechn+111) THEN
1593  facp=facp*rtcm(9)
1594  ELSE
1595  facp=facp*rtcm(10)
1596  ENDIF
1597  wdtp(i)=facp
1598  ELSE
1599 C...pi_tc -> f + fbar.
1600  fcof=1d0
1601  ika=iabs(kfdp(idc,1))
1602  IF(ika.LT.10) fcof=3d0*radc
1603  hm1=pm1
1604  hm2=pm2
1605  IF(ika.GE.4.AND.ika.LE.6) THEN
1606  fcof=fcof*rtcm(1+ika)**2
1607  hm1=pymrun(kfdp(idc,1),sh)
1608  hm2=pymrun(kfdp(idc,2),sh)
1609  ELSEIF(ika.EQ.15) THEN
1610  fcof=fcof*rtcm(8)**2
1611  ENDIF
1612  wdtp(i)=fac*fcof*(hm1+hm2)**2*
1613  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))
1614  ENDIF
1615  wdtp(i)=fudge*wdtp(i)
1616  wdtp(0)=wdtp(0)+wdtp(i)
1617  IF(mdme(idc,1).GT.0) THEN
1618  wdte(i,mdme(idc,1))=wdtp(i)*wid2
1619  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
1620  wdte(i,0)=wdte(i,mdme(idc,1))
1621  wdte(0,0)=wdte(0,0)+wdte(i,0)
1622  ENDIF
1623  340 CONTINUE
1624 
1625  ELSEIF(kfla.EQ.ktechn+211) THEN
1626 C...pi+_tc
1627  fac=(1d0/(32d0*paru(1)*rtcm(1)**2))*shr
1628  DO 350 i=1,mdcy(kc,3)
1629  idc=i+mdcy(kc,2)-1
1630  IF(mdme(idc,1).LT.0) goto 350
1631  pm1=pmas(pycomp(kfdp(idc,1)),1)
1632  pm2=pmas(pycomp(kfdp(idc,2)),1)
1633  pm3=0d0
1634  IF(i.EQ.5) pm3=pmas(pycomp(kfdp(idc,3)),1)
1635  rm1=pm1**2/sh
1636  rm2=pm2**2/sh
1637  rm3=pm3**2/sh
1638  IF(sqrt(rm1)+sqrt(rm2)+sqrt(rm3).GT.1d0) goto 350
1639  wid2=1d0
1640 C...pi_tc -> f + f'.
1641  fcof=1d0
1642  IF(iabs(kfdp(idc,1)).LT.10) fcof=3d0*radc
1643 C...pi_tc+ -> W b b~
1644  IF(i.EQ.5.AND.shr.LT.pmas(6,1)+pmas(5,1)) THEN
1645  fcof=3d0*radc
1646  xmt2=pmas(6,1)**2/sh
1647  facp=fac/(4d0*paru(1))*fcof*xmt2*rtcm(7)**2
1648  kfc3=pycomp(kfdp(idc,3))
1649  check = sqrt(rm1)+sqrt(rm2)+sqrt(rm3)
1650  check = sqrt(rm1)
1651  t0 = (1d0-check**2)*
1652  & (xmt2*(6d0*xmt2**2+3d0*xmt2*rm1-4d0*rm1**2)-
1653  & (5d0*xmt2**2+2d0*xmt2*rm1-8d0*rm1**2))/(4d0*xmt2**2)
1654  t1 = (1d0-xmt2)*(rm1-xmt2)*((xmt2**2+xmt2*rm1+4d0*rm1**2)
1655  & -3d0*xmt2**2*(xmt2+rm1))/(2d0*xmt2**3)
1656  t3 = rm1**2/xmt2**3*(3d0*xmt2-4d0*rm1+4d0*xmt2*rm1)
1657  wdtp(i)=facp*(t0 + t1*log((xmt2-check**2)/(xmt2-1d0))
1658  & +t3*log(check))
1659  IF(kflr.GT.0) THEN
1660  wid2=wids(24,2)
1661  ELSE
1662  wid2=wids(24,3)
1663  ENDIF
1664  ELSE
1665  fcof=1d0
1666  ika=iabs(kfdp(idc,1))
1667  IF(ika.LT.10) fcof=3d0*radc
1668  hm1=pm1
1669  hm2=pm2
1670  IF(i.GE.1.AND.i.LE.5) THEN
1671  IF(i.LE.2) THEN
1672  fcof=fcof*rtcm(5)**2
1673  ELSEIF(i.LE.4) THEN
1674  fcof=fcof*rtcm(6)**2
1675  ELSEIF(i.EQ.5) THEN
1676  fcof=fcof*rtcm(7)**2
1677  ENDIF
1678  hm1=pymrun(kfdp(idc,1),sh)
1679  hm2=pymrun(kfdp(idc,2),sh)
1680  ELSEIF(i.EQ.8) THEN
1681  fcof=fcof*rtcm(8)**2
1682  ENDIF
1683  wdtp(i)=fac*fcof*(hm1+hm2)**2*
1684  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))
1685  ENDIF
1686  wdtp(i)=fudge*wdtp(i)
1687  wdtp(0)=wdtp(0)+wdtp(i)
1688  IF(mdme(idc,1).GT.0) THEN
1689  wdte(i,mdme(idc,1))=wdtp(i)*wid2
1690  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
1691  wdte(i,0)=wdte(i,mdme(idc,1))
1692  wdte(0,0)=wdte(0,0)+wdte(i,0)
1693  ENDIF
1694  350 CONTINUE
1695 
1696  ELSEIF(kfla.EQ.ktechn+331) THEN
1697 C...Techni-eta.
1698  fac=(sh/parp(46)**2)*shr
1699  DO 360 i=1,mdcy(kc,3)
1700  idc=i+mdcy(kc,2)-1
1701  IF(mdme(idc,1).LT.0) goto 360
1702  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
1703  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
1704  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 360
1705  wid2=1d0
1706  IF(i.LE.2) THEN
1707  wdtp(i)=fac*rm1*sqrt(max(0d0,1d0-4d0*rm1))/(4d0*paru(1))
1708  IF(i.EQ.2) wid2=wids(6,1)
1709  ELSE
1710  wdtp(i)=fac*5d0*as**2/(96d0*paru(1)**3)
1711  ENDIF
1712  wdtp(i)=fudge*wdtp(i)
1713  wdtp(0)=wdtp(0)+wdtp(i)
1714  IF(mdme(idc,1).GT.0) THEN
1715  wdte(i,mdme(idc,1))=wdtp(i)*wid2
1716  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
1717  wdte(i,0)=wdte(i,mdme(idc,1))
1718  wdte(0,0)=wdte(0,0)+wdte(i,0)
1719  ENDIF
1720  360 CONTINUE
1721 
1722  ELSEIF(kfla.EQ.ktechn+113) THEN
1723 C...Techni-rho0:
1724  alprht=2.91d0*(3d0/itcm(1))
1725  fac=(alprht/12d0)*shr
1726  facf=(1d0/6d0)*(aem**2/alprht)*shr
1727  sqmz=pmas(23,1)**2
1728  sqmw=pmas(24,1)**2
1729  shp=sh
1730  CALL pywidx(23,shp,wdtpp,wdtep)
1731  gmmz=shr*wdtpp(0)
1732  xwrht=(1d0-2d0*xw)/(4d0*xw*(1d0-xw))
1733  bwzr=xwrht*sh*(sh-sqmz)/((sh-sqmz)**2+gmmz**2)
1734  bwzi=xwrht*sh*gmmz/((sh-sqmz)**2+gmmz**2)
1735  DO 370 i=1,mdcy(kc,3)
1736  idc=i+mdcy(kc,2)-1
1737  IF(mdme(idc,1).LT.0) goto 370
1738  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
1739  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
1740  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 370
1741  wid2=1d0
1742  IF(i.EQ.1) THEN
1743 C...rho_tc0 -> W+ + W-.
1744  wdtp(i)=fac*rtcm(3)**4*
1745  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3
1746  wid2=wids(24,1)
1747  ELSEIF(i.EQ.2) THEN
1748 C...rho_tc0 -> W+ + pi_tc-.
1749  wdtp(i)=fac*rtcm(3)**2*(1d0-rtcm(3)**2)*
1750  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3+
1751  & aem*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))*
1752  & ((1d0-rm1-rm2)**2-4d0*rm1*rm2 + 6d0*sqmw/sh)*
1753  & (1d0-rtcm(3)**2)/4d0/xw/24d0/rtcm(13)**2*shr**3
1754  wid2=wids(24,2)*wids(pycomp(ktechn+211),3)
1755  ELSEIF(i.EQ.3) THEN
1756 C...rho_tc0 -> pi_tc+ + W-.
1757  wdtp(i)=fac*rtcm(3)**2*(1d0-rtcm(3)**2)*
1758  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3+
1759  & aem*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))*
1760  & ((1d0-rm1-rm2)**2-4d0*rm1*rm2 + 6d0*sqmw/sh)*
1761  & (1d0-rtcm(3)**2)/4d0/xw/24d0/rtcm(13)**2*shr**3
1762  wid2=wids(pycomp(ktechn+211),2)*wids(24,3)
1763  ELSEIF(i.EQ.4) THEN
1764 C...rho_tc0 -> pi_tc+ + pi_tc-.
1765  wdtp(i)=fac*(1d0-rtcm(3)**2)**2*
1766  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3
1767  wid2=wids(pycomp(ktechn+211),1)
1768  ELSEIF(i.EQ.5) THEN
1769 C...rho_tc0 -> gamma + pi_tc0
1770  wdtp(i)=aem*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3*
1771  & (2d0*rtcm(2)-1d0)**2*(1d0-rtcm(3)**2)/24d0/rtcm(12)**2*
1772  & shr**3
1773  wid2=wids(pycomp(ktechn+111),2)
1774  ELSEIF(i.EQ.6) THEN
1775 C...rho_tc0 -> gamma + pi_tc0'
1776  wdtp(i)=aem*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3*
1777  & (1d0-rtcm(4)**2)/24d0/rtcm(12)**2*shr**3
1778  wid2=wids(pycomp(ktechn+221),2)
1779  ELSEIF(i.EQ.7) THEN
1780 C...rho_tc0 -> Z0 + pi_tc0
1781  wdtp(i)=aem*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3*
1782  & (2d0*rtcm(2)-1d0)**2*(1d0-rtcm(3)**2)/24d0/rtcm(12)**2*
1783  & xw/xw1*shr**3
1784  wid2=wids(23,2)*wids(pycomp(ktechn+111),2)
1785  ELSEIF(i.EQ.8) THEN
1786 C...rho_tc0 -> Z0 + pi_tc0'
1787  wdtp(i)=aem*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3*
1788  & (1d0-rtcm(4)**2)/24d0/rtcm(12)**2*(1d0-2d0*xw)**2/4d0/
1789  & xw/xw1*shr**3
1790  wid2=wids(23,2)*wids(pycomp(ktechn+221),2)
1791  ELSE
1792 C...rho_tc0 -> f + fbar.
1793  wid2=1d0
1794  IF(i.LE.16) THEN
1795  ia=i-8
1796  fcof=3d0*radc
1797  IF(ia.GE.6.AND.ia.LE.8) wid2=wids(ia,1)
1798  ELSE
1799  ia=i-6
1800  fcof=1d0
1801  IF(ia.GE.17) wid2=wids(ia,1)
1802  ENDIF
1803  ei=kchg(ia,1)/3d0
1804  ai=sign(1d0,ei+0.1d0)
1805  vi=ai-4d0*ei*xwv
1806  vali=0.5d0*(vi+ai)
1807  vari=0.5d0*(vi-ai)
1808  wdtp(i)=facf*fcof*sqrt(max(0d0,1d0-4d0*rm1))*((1d0-rm1)*
1809  & ((ei+vali*bwzr)**2+(vali*bwzi)**2+
1810  & (ei+vari*bwzr)**2+(vari*bwzi)**2)+6d0*rm1*(
1811  & (ei+vali*bwzr)*(ei+vari*bwzr)+vali*vari*bwzi**2))
1812  ENDIF
1813  wdtp(i)=fudge*wdtp(i)
1814  wdtp(0)=wdtp(0)+wdtp(i)
1815  IF(mdme(idc,1).GT.0) THEN
1816  wdte(i,mdme(idc,1))=wdtp(i)*wid2
1817  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
1818  wdte(i,0)=wdte(i,mdme(idc,1))
1819  wdte(0,0)=wdte(0,0)+wdte(i,0)
1820  ENDIF
1821  370 CONTINUE
1822 
1823  ELSEIF(kfla.EQ.ktechn+213) THEN
1824 C...Techni-rho+/-:
1825  alprht=2.91d0*(3d0/itcm(1))
1826  fac=(alprht/12d0)*shr
1827  sqmz=pmas(23,1)**2
1828  sqmw=pmas(24,1)**2
1829  shp=sh
1830  CALL pywidx(24,shp,wdtpp,wdtep)
1831  gmmw=shr*wdtpp(0)
1832  facf=(1d0/12d0)*(aem**2/alprht)*shr*
1833  & (0.125d0/xw**2)*sh**2/((sh-sqmw)**2+gmmw**2)
1834  DO 380 i=1,mdcy(kc,3)
1835  idc=i+mdcy(kc,2)-1
1836  IF(mdme(idc,1).LT.0) goto 380
1837  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
1838  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
1839  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 380
1840  wid2=1d0
1841  IF(i.EQ.1) THEN
1842 C...rho_tc+ -> W+ + Z0.
1843  wdtp(i)=fac*rtcm(3)**4*
1844  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3
1845  IF(kflr.GT.0) THEN
1846  wid2=wids(24,2)*wids(23,2)
1847  ELSE
1848  wid2=wids(24,3)*wids(23,2)
1849  ENDIF
1850  ELSEIF(i.EQ.2) THEN
1851 C...rho_tc+ -> W+ + pi_tc0.
1852  wdtp(i)=fac*rtcm(3)**2*(1d0-rtcm(3)**2)*
1853  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3+
1854  & aem*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))*
1855  & ((1d0-rm1-rm2)**2-4d0*rm1*rm2 + 6d0*sqmw/sh)*
1856  & (1d0-rtcm(3)**2)/4d0/xw/24d0/rtcm(13)**2*shr**3
1857  IF(kflr.GT.0) THEN
1858  wid2=wids(24,2)*wids(pycomp(ktechn+111),2)
1859  ELSE
1860  wid2=wids(24,3)*wids(pycomp(ktechn+111),2)
1861  ENDIF
1862  ELSEIF(i.EQ.3) THEN
1863 C...rho_tc+ -> pi_tc+ + Z0.
1864  wdtp(i)=fac*rtcm(3)**2*(1d0-rtcm(3)**2)*
1865  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3+
1866  & aem*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))*
1867  & ((1d0-rm1-rm2)**2-4d0*rm1*rm2 + 6d0*sqmz/sh)*
1868  & (1d0-rtcm(3)**2)/4d0/xw/xw1/24d0/rtcm(13)**2*shr**3+
1869  & aem*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3*
1870  & (2d0*rtcm(2)-1d0)**2*(1d0-rtcm(3)**2)/24d0/rtcm(12)**2*
1871  & shr**3*xw/xw1
1872  IF(kflr.GT.0) THEN
1873  wid2=wids(pycomp(ktechn+211),2)*wids(23,2)
1874  ELSE
1875  wid2=wids(pycomp(ktechn+211),3)*wids(23,2)
1876  ENDIF
1877  ELSEIF(i.EQ.4) THEN
1878 C...rho_tc+ -> pi_tc+ + pi_tc0.
1879  wdtp(i)=fac*(1d0-rtcm(3)**2)**2*
1880  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3
1881  IF(kflr.GT.0) THEN
1882  wid2=wids(pycomp(ktechn+211),2)*wids(pycomp(ktechn+111),2)
1883  ELSE
1884  wid2=wids(pycomp(ktechn+211),3)*wids(pycomp(ktechn+111),2)
1885  ENDIF
1886  ELSEIF(i.EQ.5) THEN
1887 C...rho_tc+ -> pi_tc+ + gamma
1888  wdtp(i)=aem*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3*
1889  & (2d0*rtcm(2)-1d0)**2*(1d0-rtcm(3)**2)/24d0/rtcm(12)**2*
1890  & shr**3
1891  IF(kflr.GT.0) THEN
1892  wid2=wids(pycomp(ktechn+211),2)
1893  ELSE
1894  wid2=wids(pycomp(ktechn+211),3)
1895  ENDIF
1896  ELSEIF(i.EQ.6) THEN
1897 C...rho_tc+ -> W+ + pi_tc0'
1898  wdtp(i)=aem*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3*
1899  & (1d0-rtcm(4)**2)/4d0/xw/24d0/rtcm(12)**2*shr**3
1900  IF(kflr.GT.0) THEN
1901  wid2=wids(24,2)*wids(pycomp(ktechn+221),2)
1902  ELSE
1903  wid2=wids(24,3)*wids(pycomp(ktechn+221),2)
1904  ENDIF
1905  ELSE
1906 C...rho_tc+ -> f + fbar'.
1907  ia=i-6
1908  wid2=1d0
1909  IF(ia.LE.16) THEN
1910  fcof=3d0*radc*vckm((ia-1)/4+1,mod(ia-1,4)+1)
1911  IF(kflr.GT.0) THEN
1912  IF(mod(ia,4).EQ.3) wid2=wids(6,2)
1913  IF(mod(ia,4).EQ.0) wid2=wids(8,2)
1914  IF(ia.GE.13) wid2=wid2*wids(7,3)
1915  ELSE
1916  IF(mod(ia,4).EQ.3) wid2=wids(6,3)
1917  IF(mod(ia,4).EQ.0) wid2=wids(8,3)
1918  IF(ia.GE.13) wid2=wid2*wids(7,2)
1919  ENDIF
1920  ELSE
1921  fcof=1d0
1922  IF(kflr.GT.0) THEN
1923  IF(ia.EQ.20) wid2=wids(17,3)*wids(18,2)
1924  ELSE
1925  IF(ia.EQ.20) wid2=wids(17,2)*wids(18,3)
1926  ENDIF
1927  ENDIF
1928  wdtp(i)=facf*fcof*(2d0-rm1-rm2-(rm1-rm2)**2)*
1929  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))
1930  ENDIF
1931  wdtp(i)=fudge*wdtp(i)
1932  wdtp(0)=wdtp(0)+wdtp(i)
1933  IF(mdme(idc,1).GT.0) THEN
1934  wdte(i,mdme(idc,1))=wdtp(i)*wid2
1935  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
1936  wdte(i,0)=wdte(i,mdme(idc,1))
1937  wdte(0,0)=wdte(0,0)+wdte(i,0)
1938  ENDIF
1939  380 CONTINUE
1940 
1941  ELSEIF(kfla.EQ.ktechn+223) THEN
1942 C...Techni-omega:
1943  alprht=2.91d0*(3d0/itcm(1))
1944  fac=(alprht/12d0)*shr
1945  facf=(1d0/6d0)*(aem**2/alprht)*shr*(2d0*rtcm(2)-1d0)**2
1946  sqmz=pmas(23,1)**2
1947  shp=sh
1948  CALL pywidx(23,shp,wdtpp,wdtep)
1949  gmmz=shr*wdtpp(0)
1950  bwzr=(0.5d0/(1d0-xw))*sh*(sh-sqmz)/((sh-sqmz)**2+gmmz**2)
1951  bwzi=-(0.5d0/(1d0-xw))*sh*gmmz/((sh-sqmz)**2+gmmz**2)
1952  DO 390 i=1,mdcy(kc,3)
1953  idc=i+mdcy(kc,2)-1
1954  IF(mdme(idc,1).LT.0) goto 390
1955  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
1956  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
1957  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 390
1958  wid2=1d0
1959  IF(i.EQ.1) THEN
1960 C...omega_tc0 -> gamma + pi_tc0.
1961  wdtp(i)=aem/24d0/rtcm(12)**2*(1d0-rtcm(3)**2)*
1962  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3*shr**3
1963  wid2=wids(pycomp(ktechn+111),2)
1964  ELSEIF(i.EQ.2) THEN
1965 C...omega_tc0 -> Z0 + pi_tc0
1966  wdtp(i)=aem*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3*
1967  & (1d0-rtcm(3)**2)/24d0/rtcm(12)**2*(1d0-2d0*xw)**2/4d0/
1968  & xw/xw1*shr**3
1969  wid2=wids(23,2)*wids(pycomp(ktechn+111),2)
1970  ELSEIF(i.EQ.3) THEN
1971 C...omega_tc0 -> gamma + pi_tc0'
1972  wdtp(i)=aem*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3*
1973  & (2d0*rtcm(2)-1d0)**2*(1d0-rtcm(4)**2)/24d0/rtcm(12)**2*
1974  & shr**3
1975  wid2=wids(pycomp(ktechn+221),2)
1976  ELSEIF(i.EQ.4) THEN
1977 C...omega_tc0 -> Z0 + pi_tc0'
1978  wdtp(i)=aem*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3*
1979  & (2d0*rtcm(2)-1d0)**2*(1d0-rtcm(4)**2)/24d0/rtcm(12)**2*
1980  & xw/xw1*shr**3
1981  wid2=wids(23,2)*wids(pycomp(ktechn+221),2)
1982  ELSEIF(i.EQ.5) THEN
1983 C...omega_tc0 -> W+ + pi_tc-
1984  wdtp(i)=aem*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3*
1985  & (1d0-rtcm(3)**2)/4d0/xw/24d0/rtcm(12)**2*shr**3+
1986  & fac*rtcm(3)**2*(1d0-rtcm(3)**2)*rtcm(11)**2*
1987  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3
1988  wid2=wids(24,2)*wids(pycomp(ktechn+211),3)
1989  ELSEIF(i.EQ.6) THEN
1990 C...omega_tc0 -> pi_tc+ + W-
1991  wdtp(i)=aem*sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3*
1992  & (1d0-rtcm(3)**2)/4d0/xw/24d0/rtcm(12)**2*shr**3+
1993  & fac*rtcm(3)**2*(1d0-rtcm(3)**2)*rtcm(11)**2*
1994  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3
1995  wid2=wids(24,3)*wids(pycomp(ktechn+211),2)
1996  ELSEIF(i.EQ.7) THEN
1997 C...omega_tc0 -> W+ + W-.
1998  wdtp(i)=fac*rtcm(3)**4*rtcm(11)**2*
1999  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3
2000  wid2=wids(24,1)
2001  ELSEIF(i.EQ.8) THEN
2002 C...omega_tc0 -> pi_tc+ + pi_tc-.
2003  wdtp(i)=fac*(1d0-rtcm(3)**2)**2*rtcm(11)**2*
2004  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))**3
2005  wid2=wids(pycomp(ktechn+211),1)
2006  ELSE
2007 C...omega_tc0 -> f + fbar.
2008  wid2=1d0
2009  IF(i.LE.14) THEN
2010  ia=i-8
2011  fcof=3d0*radc
2012  IF(ia.GE.6.AND.ia.LE.8) wid2=wids(ia,1)
2013  ELSE
2014  ia=i-6
2015  fcof=1d0
2016  IF(ia.GE.17) wid2=wids(ia,1)
2017  ENDIF
2018  ei=kchg(ia,1)/3d0
2019  ai=sign(1d0,ei+0.1d0)
2020  vi=ai-4d0*ei*xwv
2021  vali=-0.5d0*(vi+ai)
2022  vari=-0.5d0*(vi-ai)
2023  wdtp(i)=facf*fcof*sqrt(max(0d0,1d0-4d0*rm1))*((1d0-rm1)*
2024  & ((ei+vali*bwzr)**2+(vali*bwzi)**2+
2025  & (ei+vari*bwzr)**2+(vari*bwzi)**2)+6d0*rm1*(
2026  & (ei+vali*bwzr)*(ei+vari*bwzr)+vali*vari*bwzi**2))
2027  ENDIF
2028  wdtp(i)=fudge*wdtp(i)
2029  wdtp(0)=wdtp(0)+wdtp(i)
2030  IF(mdme(idc,1).GT.0) THEN
2031  wdte(i,mdme(idc,1))=wdtp(i)*wid2
2032  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
2033  wdte(i,0)=wdte(i,mdme(idc,1))
2034  wdte(0,0)=wdte(0,0)+wdte(i,0)
2035  ENDIF
2036  390 CONTINUE
2037 
2038 C.....V8 -> quark anti-quark
2039  ELSEIF(kfla.EQ.ktechn+100021) THEN
2040  fac=as/6d0*shr
2041  tant3=rtcm(21)
2042  IF(itcm(2).EQ.0) THEN
2043  imdl=1
2044  ELSEIF(itcm(2).EQ.1) THEN
2045  imdl=2
2046  ENDIF
2047  DO 400 i=1,mdcy(kc,3)
2048  idc=i+mdcy(kc,2)-1
2049  IF(mdme(idc,1).LT.0) goto 400
2050  pm1=pmas(pycomp(kfdp(idc,1)),1)
2051  rm1=pm1**2/sh
2052  IF(rm1.GT.0.25d0) goto 400
2053  wid2=1d0
2054  IF(i.EQ.5.OR.i.EQ.6.OR.imdl.EQ.2) THEN
2055  fmix=1d0/tant3**2
2056  ELSE
2057  fmix=tant3**2
2058  ENDIF
2059  wdtp(i)=fac*(1d0+2d0*rm1)*sqrt(1d0-4d0*rm1)*fmix
2060  IF(i.EQ.6) wid2=wids(6,1)
2061  wdtp(i)=fudge*wdtp(i)
2062  wdtp(0)=wdtp(0)+wdtp(i)
2063  IF(mdme(idc,1).GT.0) THEN
2064  wdte(i,mdme(idc,1))=wdtp(i)*wid2
2065  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
2066  wdte(i,0)=wdte(i,mdme(idc,1))
2067  wdte(0,0)=wdte(0,0)+wdte(i,0)
2068  ENDIF
2069  400 CONTINUE
2070 
2071  ELSEIF(kfla.EQ.ktechn+100111.OR.kfla.EQ.ktechn+200111) THEN
2072  fac=(1d0/(4d0*paru(1)*rtcm(1)**2))*shr
2073  clebf=0d0
2074  DO 410 i=1,mdcy(kc,3)
2075  idc=i+mdcy(kc,2)-1
2076  IF(mdme(idc,1).LT.0) goto 410
2077  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
2078  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
2079  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 410
2080  wid2=1d0
2081 C...pi_tc -> g + g
2082  IF(i.EQ.7) THEN
2083  IF(kfla.EQ.ktechn+100111) THEN
2084  clebg=4d0/3d0
2085  ELSE
2086  clebg=5d0/3d0
2087  ENDIF
2088  facp=(as/(8d0*paru(1))*itcm(1)/rtcm(1))**2
2089  & /(2d0*paru(1))*sh*shr*clebg
2090  wdtp(i)=facp
2091  ELSE
2092 C...pi_tc -> f + fbar.
2093  IF(i.EQ.6) wid2=wids(6,1)
2094  fcof=1d0
2095  ika=iabs(kfdp(idc,1))
2096  IF(ika.LT.10) fcof=3d0*radc
2097  hm1=pymrun(kfdp(idc,1),sh)
2098  wdtp(i)=fac*fcof*hm1**2*clebf*
2099  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))
2100  ENDIF
2101  wdtp(i)=fudge*wdtp(i)
2102  wdtp(0)=wdtp(0)+wdtp(i)
2103  IF(mdme(idc,1).GT.0) THEN
2104  wdte(i,mdme(idc,1))=wdtp(i)*wid2
2105  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
2106  wdte(i,0)=wdte(i,mdme(idc,1))
2107  wdte(0,0)=wdte(0,0)+wdte(i,0)
2108  ENDIF
2109  410 CONTINUE
2110 
2111  ELSEIF(kfla.GE.ktechn+100113.AND.kfla.LE.ktechn+400113) THEN
2112  fac=as/6d0*shr
2113  alprht=2.91d0*(3d0/itcm(1))
2114  tant3=rtcm(21)
2115  sin2t=2d0*tant3/(tant3**2+1d0)
2116  sint3=tant3/sqrt(tant3**2+1d0)
2117  csxpp=rtcm(22)
2118  rm82=rtcm(27)**2
2119  x12=(rtcm(29)*sqrt(1d0-rtcm(29)**2)*cos(rtcm(30))+
2120  & rtcm(31)*sqrt(1d0-rtcm(31)**2)*cos(rtcm(32)))/sqrt(2d0)
2121  x21=(rtcm(29)*sqrt(1d0-rtcm(29)**2)*sin(rtcm(30))+
2122  & rtcm(31)*sqrt(1d0-rtcm(31)**2)*sin(rtcm(32)))/sqrt(2d0)
2123  x11=(.25d0*(rtcm(29)**2+rtcm(31)**2+2d0)-
2124  & sint3**2)*2d0
2125  x22=(.25d0*(2d0-rtcm(29)**2-rtcm(31)**2)-
2126  & sint3**2)*2d0
2127  CALL pywidx(ktechn+100021,sh,wdtpp,wdtep)
2128 
2129  IF(wdtpp(0).GT.rtcm(33)*shr) wdtpp(0)=rtcm(33)*shr
2130  gmv8=shr*wdtpp(0)
2131  rmv8=pmas(pycomp(ktechn+100021),1)
2132  fv8re=sh*(sh-rmv8**2)/((sh-rmv8**2)**2+gmv8**2)
2133  fv8im=sh*gmv8/((sh-rmv8**2)**2+gmv8**2)
2134  IF(itcm(2).EQ.0) THEN
2135  imdl=1
2136  ELSE
2137  imdl=2
2138  ENDIF
2139  DO 420 i=1,mdcy(kc,3)
2140  IF(i.EQ.7.AND.(kfla.EQ.ktechn+200113.OR.
2141  & kfla.EQ.ktechn+300113)) goto 420
2142  idc=i+mdcy(kc,2)-1
2143  IF(mdme(idc,1).LT.0) goto 420
2144  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
2145  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
2146  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 420
2147  wid2=1d0
2148  IF(i.LE.6) THEN
2149  IF(i.EQ.6) wid2=wids(6,1)
2150  xig=1d0
2151  IF(kfla.EQ.ktechn+200113) THEN
2152  xig=0d0
2153  xij=x12
2154  ELSEIF(kfla.EQ.ktechn+300113) THEN
2155  xig=0d0
2156  xij=x21
2157  ELSEIF(kfla.EQ.ktechn+100113) THEN
2158  xij=x11
2159  ELSE
2160  xij=x22
2161  ENDIF
2162  IF(i.EQ.5.OR.i.EQ.6.OR.imdl.EQ.2) THEN
2163  fmix=1d0/tant3/sin2t
2164  ELSE
2165  fmix=-tant3/sin2t
2166  ENDIF
2167  xfac=(xig+fmix*xij*fv8re)**2+(fmix*xij*fv8im)**2
2168  wdtp(i)=fac*(1d0+2d0*rm1)*sqrt(1d0-4d0*rm1)*as/alprht*xfac
2169  ELSEIF(i.EQ.7) THEN
2170  wdtp(i)=shr*as**2/(4d0*alprht)
2171  ELSEIF(kfla.EQ.ktechn+400113.AND.i.LE.9) THEN
2172  psh=shr*(1d0-rm1)/2d0
2173  wdtp(i)=as/9d0*psh**3/rm82
2174  IF(i.EQ.8) THEN
2175  wdtp(i)=2d0*wdtp(i)*csxpp**2
2176  wid2=wids(pycomp(kfdp(idc,1)),2)
2177  ELSE
2178  wdtp(i)=5d0*wdtp(i)
2179  wid2=wids(pycomp(kfdp(idc,1)),2)
2180  ENDIF
2181  ENDIF
2182  wdtp(i)=fudge*wdtp(i)
2183  wdtp(0)=wdtp(0)+wdtp(i)
2184  IF(mdme(idc,1).GT.0) THEN
2185  wdte(i,mdme(idc,1))=wdtp(i)*wid2
2186  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
2187  wdte(i,0)=wdte(i,mdme(idc,1))
2188  wdte(0,0)=wdte(0,0)+wdte(i,0)
2189  ENDIF
2190  420 CONTINUE
2191 
2192  ELSEIF(kfla.EQ.kexcit+1) THEN
2193 C...d* excited quark.
2194  fac=(sh/rtcm(41)**2)*shr
2195  DO 430 i=1,mdcy(kc,3)
2196  idc=i+mdcy(kc,2)-1
2197  IF(mdme(idc,1).LT.0) goto 430
2198  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
2199  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
2200  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 430
2201  wid2=1d0
2202  IF(i.EQ.1) THEN
2203 C...d* -> g + d.
2204  wdtp(i)=fac*as*rtcm(45)**2/3d0
2205  wid2=1d0
2206  ELSEIF(i.EQ.2) THEN
2207 C...d* -> gamma + d.
2208  qf=-rtcm(43)/2d0+rtcm(44)/6d0
2209  wdtp(i)=fac*aem*qf**2/4d0
2210  wid2=1d0
2211  ELSEIF(i.EQ.3) THEN
2212 C...d* -> Z0 + d.
2213  qf=-rtcm(43)*xw1/2d0-rtcm(44)*xw/6d0
2214  wdtp(i)=fac*aem*qf**2/(8d0*xw*xw1)*
2215  & (1d0-rm1)**2*(2d0+rm1)
2216  wid2=wids(23,2)
2217  ELSEIF(i.EQ.4) THEN
2218 C...d* -> W- + u.
2219  wdtp(i)=fac*aem*rtcm(43)**2/(16d0*xw)*
2220  & (1d0-rm1)**2*(2d0+rm1)
2221  IF(kflr.GT.0) wid2=wids(24,3)
2222  IF(kflr.LT.0) wid2=wids(24,2)
2223  ENDIF
2224  wdtp(i)=fudge*wdtp(i)
2225  wdtp(0)=wdtp(0)+wdtp(i)
2226  IF(mdme(idc,1).GT.0) THEN
2227  wdte(i,mdme(idc,1))=wdtp(i)*wid2
2228  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
2229  wdte(i,0)=wdte(i,mdme(idc,1))
2230  wdte(0,0)=wdte(0,0)+wdte(i,0)
2231  ENDIF
2232  430 CONTINUE
2233 
2234  ELSEIF(kfla.EQ.kexcit+2) THEN
2235 C...u* excited quark.
2236  fac=(sh/rtcm(41)**2)*shr
2237  DO 440 i=1,mdcy(kc,3)
2238  idc=i+mdcy(kc,2)-1
2239  IF(mdme(idc,1).LT.0) goto 440
2240  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
2241  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
2242  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 440
2243  wid2=1d0
2244  IF(i.EQ.1) THEN
2245 C...u* -> g + u.
2246  wdtp(i)=fac*as*rtcm(45)**2/3d0
2247  wid2=1d0
2248  ELSEIF(i.EQ.2) THEN
2249 C...u* -> gamma + u.
2250  qf=rtcm(43)/2d0+rtcm(44)/6d0
2251  wdtp(i)=fac*aem*qf**2/4d0
2252  wid2=1d0
2253  ELSEIF(i.EQ.3) THEN
2254 C...u* -> Z0 + u.
2255  qf=rtcm(43)*xw1/2d0-rtcm(44)*xw/6d0
2256  wdtp(i)=fac*aem*qf**2/(8d0*xw*xw1)*
2257  & (1d0-rm1)**2*(2d0+rm1)
2258  wid2=wids(23,2)
2259  ELSEIF(i.EQ.4) THEN
2260 C...u* -> W+ + d.
2261  wdtp(i)=fac*aem*rtcm(43)**2/(16d0*xw)*
2262  & (1d0-rm1)**2*(2d0+rm1)
2263  IF(kflr.GT.0) wid2=wids(24,2)
2264  IF(kflr.LT.0) wid2=wids(24,3)
2265  ENDIF
2266  wdtp(i)=fudge*wdtp(i)
2267  wdtp(0)=wdtp(0)+wdtp(i)
2268  IF(mdme(idc,1).GT.0) THEN
2269  wdte(i,mdme(idc,1))=wdtp(i)*wid2
2270  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
2271  wdte(i,0)=wdte(i,mdme(idc,1))
2272  wdte(0,0)=wdte(0,0)+wdte(i,0)
2273  ENDIF
2274  440 CONTINUE
2275 
2276  ELSEIF(kfla.EQ.kexcit+11) THEN
2277 C...e* excited lepton.
2278  fac=(sh/rtcm(41)**2)*shr
2279  DO 450 i=1,mdcy(kc,3)
2280  idc=i+mdcy(kc,2)-1
2281  IF(mdme(idc,1).LT.0) goto 450
2282  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
2283  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
2284  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 450
2285  wid2=1d0
2286  IF(i.EQ.1) THEN
2287 C...e* -> gamma + e.
2288  qf=-rtcm(43)/2d0-rtcm(44)/2d0
2289  wdtp(i)=fac*aem*qf**2/4d0
2290  wid2=1d0
2291  ELSEIF(i.EQ.2) THEN
2292 C...e* -> Z0 + e.
2293  qf=-rtcm(43)*xw1/2d0+rtcm(44)*xw/2d0
2294  wdtp(i)=fac*aem*qf**2/(8d0*xw*xw1)*
2295  & (1d0-rm1)**2*(2d0+rm1)
2296  wid2=wids(23,2)
2297  ELSEIF(i.EQ.3) THEN
2298 C...e* -> W- + nu.
2299  wdtp(i)=fac*aem*rtcm(43)**2/(16d0*xw)*
2300  & (1d0-rm1)**2*(2d0+rm1)
2301  IF(kflr.GT.0) wid2=wids(24,3)
2302  IF(kflr.LT.0) wid2=wids(24,2)
2303  ENDIF
2304  wdtp(i)=fudge*wdtp(i)
2305  wdtp(0)=wdtp(0)+wdtp(i)
2306  IF(mdme(idc,1).GT.0) THEN
2307  wdte(i,mdme(idc,1))=wdtp(i)*wid2
2308  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
2309  wdte(i,0)=wdte(i,mdme(idc,1))
2310  wdte(0,0)=wdte(0,0)+wdte(i,0)
2311  ENDIF
2312  450 CONTINUE
2313 
2314  ELSEIF(kfla.EQ.kexcit+12) THEN
2315 C...nu*_e excited neutrino.
2316  fac=(sh/rtcm(41)**2)*shr
2317  DO 460 i=1,mdcy(kc,3)
2318  idc=i+mdcy(kc,2)-1
2319  IF(mdme(idc,1).LT.0) goto 460
2320  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
2321  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
2322  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 460
2323  wid2=1d0
2324  IF(i.EQ.1) THEN
2325 C...nu*_e -> Z0 + nu*_e.
2326  qf=rtcm(43)*xw1/2d0+rtcm(44)*xw/2d0
2327  wdtp(i)=fac*aem*qf**2/(8d0*xw*xw1)*
2328  & (1d0-rm1)**2*(2d0+rm1)
2329  wid2=wids(23,2)
2330  ELSEIF(i.EQ.2) THEN
2331 C...nu*_e -> W+ + e.
2332  wdtp(i)=fac*aem*rtcm(43)**2/(16d0*xw)*
2333  & (1d0-rm1)**2*(2d0+rm1)
2334  IF(kflr.GT.0) wid2=wids(24,2)
2335  IF(kflr.LT.0) wid2=wids(24,3)
2336  ENDIF
2337  wdtp(i)=fudge*wdtp(i)
2338  wdtp(0)=wdtp(0)+wdtp(i)
2339  IF(mdme(idc,1).GT.0) THEN
2340  wdte(i,mdme(idc,1))=wdtp(i)*wid2
2341  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
2342  wdte(i,0)=wdte(i,mdme(idc,1))
2343  wdte(0,0)=wdte(0,0)+wdte(i,0)
2344  ENDIF
2345  460 CONTINUE
2346 
2347  ELSEIF(kfla.EQ.kdimen+39) THEN
2348 C...G* (graviton resonance):
2349  fac=(parp(50)**2/paru(1))*shr
2350  DO 470 i=1,mdcy(kc,3)
2351  idc=i+mdcy(kc,2)-1
2352  IF(mdme(idc,1).LT.0) goto 470
2353  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
2354  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
2355  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 470
2356  wid2=1d0
2357  IF(i.LE.8) THEN
2358 C...G* -> q + qbar
2359  fcof=3d0*radc
2360  IF(i.GE.6.AND.mstp(35).GE.1) fcof=fcof*
2361  & pyhfth(sh,sh*rm1,1d0)
2362  wdtp(i)=fac*fcof*sqrt(max(0d0,1d0-4d0*rm1))**3*
2363  & (1d0+8d0*rm1/3d0)/320d0
2364  IF(i.EQ.6) wid2=wids(6,1)
2365  IF(i.EQ.7.OR.i.EQ.8) wid2=wids(i,1)
2366  ELSEIF(i.LE.16) THEN
2367 C...G* -> l+ + l-, nu + nubar
2368  fcof=1d0
2369  wdtp(i)=fac*sqrt(max(0d0,1d0-4d0*rm1))**3*
2370  & (1d0+8d0*rm1/3d0)/320d0
2371  IF(i.EQ.15.OR.i.EQ.16) wid2=wids(2+i,1)
2372  ELSEIF(i.EQ.17) THEN
2373 C...G* -> g + g.
2374  wdtp(i)=fac/20d0
2375  ELSEIF(i.EQ.18) THEN
2376 C...G* -> gamma + gamma.
2377  wdtp(i)=fac/160d0
2378  ELSEIF(i.EQ.19) THEN
2379 C...G* -> Z0 + Z0.
2380  wdtp(i)=fac*sqrt(max(0d0,1d0-4d0*rm1))*(13d0/12d0+
2381  & 14d0*rm1/3d0+4d0*rm1**2)/160d0
2382  wid2=wids(23,1)
2383  ELSEIF(i.EQ.20) THEN
2384 C...G* -> W+ + W-.
2385  wdtp(i)=fac*sqrt(max(0d0,1d0-4d0*rm1))*(13d0/12d0+
2386  & 14d0*rm1/3d0+4d0*rm1**2)/80d0
2387  wid2=wids(24,1)
2388  ENDIF
2389  wdtp(i)=fudge*wdtp(i)
2390  wdtp(0)=wdtp(0)+wdtp(i)
2391  IF(mdme(idc,1).GT.0) THEN
2392  wdte(i,mdme(idc,1))=wdtp(i)*wid2
2393  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
2394  wdte(i,0)=wdte(i,mdme(idc,1))
2395  wdte(0,0)=wdte(0,0)+wdte(i,0)
2396  ENDIF
2397  470 CONTINUE
2398 
2399  ELSEIF(kfla.EQ.9900012.OR.kfla.EQ.9900014.OR.kfla.EQ.9900016) THEN
2400 C...nu_eR, nu_muR, nu_tauR: righthanded Majorana neutrinos.
2401  pmwr=max(1.001d0*shr,pmas(pycomp(9900024),1))
2402  fac=(aem**2/(768d0*paru(1)*xw**2))*shr**5/pmwr**4
2403  DO 480 i=1,mdcy(kc,3)
2404  idc=i+mdcy(kc,2)-1
2405  IF(mdme(idc,1).LT.0) goto 480
2406  pm1=pmas(pycomp(kfdp(idc,1)),1)
2407  pm2=pmas(pycomp(kfdp(idc,2)),1)
2408  pm3=pmas(pycomp(kfdp(idc,3)),1)
2409  IF(pm1+pm2+pm3.GE.shr) goto 480
2410  wid2=1d0
2411  IF(i.LE.9) THEN
2412 C...nu_lR -> l- qbar q'
2413  fcof=3d0*radc*vckm((i-1)/3+1,mod(i-1,3)+1)
2414  IF(mod(i,3).EQ.0) wid2=wids(6,2)
2415  ELSEIF(i.LE.18) THEN
2416 C...nu_lR -> l+ q qbar'
2417  fcof=3d0*radc*vckm((i-10)/3+1,mod(i-10,3)+1)
2418  IF(mod(i-9,3).EQ.0) wid2=wids(6,3)
2419  ELSE
2420 C...nu_lR -> l- l'+ nu_lR' + charge conjugate.
2421  fcof=1d0
2422  wid2=wids(pycomp(kfdp(idc,3)),2)
2423  ENDIF
2424  x=(pm1+pm2+pm3)/shr
2425  fx=1d0-8d0*x**2+8d0*x**6-x**8-24d0*x**4*log(x)
2426  y=(shr/pmwr)**2
2427  fy=(12d0*(1d0-y)*log(1d0-y)+12d0*y-6d0*y**2-2d0*y**3)/y**4
2428  wdtp(i)=fac*fcof*fx*fy
2429  wdtp(i)=fudge*wdtp(i)
2430  wdtp(0)=wdtp(0)+wdtp(i)
2431  IF(mdme(idc,1).GT.0) THEN
2432  wdte(i,mdme(idc,1))=wdtp(i)*wid2
2433  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
2434  wdte(i,0)=wdte(i,mdme(idc,1))
2435  wdte(0,0)=wdte(0,0)+wdte(i,0)
2436  ENDIF
2437  480 CONTINUE
2438 
2439  ELSEIF(kfla.EQ.9900023) THEN
2440 C...Z_R0:
2441  fac=(aem/(48d0*xw*xw1*(1d0-2d0*xw)))*shr
2442  DO 490 i=1,mdcy(kc,3)
2443  idc=i+mdcy(kc,2)-1
2444  IF(mdme(idc,1).LT.0) goto 490
2445  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
2446  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
2447  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 490
2448  wid2=1d0
2449  symmet=1d0
2450  IF(i.LE.6) THEN
2451 C...Z_R0 -> q + qbar
2452  ef=kchg(i,1)/3d0
2453  af=sign(1d0,ef+0.1d0)*(1d0-2d0*xw)
2454  vf=sign(1d0,ef+0.1d0)-4d0*ef*xw
2455  fcof=3d0*radc
2456  IF(i.EQ.6) wid2=wids(6,1)
2457  ELSEIF(i.EQ.7.OR.i.EQ.10.OR.i.EQ.13) THEN
2458 C...Z_R0 -> l+ + l-
2459  af=-(1d0-2d0*xw)
2460  vf=-1d0+4d0*xw
2461  fcof=1d0
2462  ELSEIF(i.EQ.8.OR.i.EQ.11.OR.i.EQ.14) THEN
2463 C...Z0 -> nu_L + nu_Lbar, assumed Majorana.
2464  af=-2d0*xw
2465  vf=0d0
2466  fcof=1d0
2467  symmet=0.5d0
2468  ELSEIF(i.LE.15) THEN
2469 C...Z0 -> nu_R + nu_R, assumed Majorana.
2470  af=2d0*xw1
2471  vf=0d0
2472  fcof=1d0
2473  wid2=wids(pycomp(kfdp(idc,1)),1)
2474  symmet=0.5d0
2475  ENDIF
2476  wdtp(i)=fac*fcof*(vf**2*(1d0+2d0*rm1)+af**2*(1d0-4d0*rm1))*
2477  & sqrt(max(0d0,1d0-4d0*rm1))*symmet
2478  wdtp(i)=fudge*wdtp(i)
2479  wdtp(0)=wdtp(0)+wdtp(i)
2480  IF(mdme(idc,1).GT.0) THEN
2481  wdte(i,mdme(idc,1))=wdtp(i)*wid2
2482  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
2483  wdte(i,0)=wdte(i,mdme(idc,1))
2484  wdte(0,0)=wdte(0,0)+wdte(i,0)
2485  ENDIF
2486  490 CONTINUE
2487 
2488  ELSEIF(kfla.EQ.9900024) THEN
2489 C...W_R+/-:
2490  fac=(aem/(24d0*xw))*shr
2491  DO 500 i=1,mdcy(kc,3)
2492  idc=i+mdcy(kc,2)-1
2493  IF(mdme(idc,1).LT.0) goto 500
2494  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
2495  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
2496  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 500
2497  wid2=1d0
2498  IF(i.LE.9) THEN
2499 C...W_R+/- -> q + qbar'
2500  fcof=3d0*radc*vckm((i-1)/3+1,mod(i-1,3)+1)
2501  IF(kflr.GT.0) THEN
2502  IF(mod(i,3).EQ.0) wid2=wids(6,2)
2503  ELSE
2504  IF(mod(i,3).EQ.0) wid2=wids(6,3)
2505  ENDIF
2506  ELSEIF(i.LE.12) THEN
2507 C...W_R+/- -> l+/- + nu_R
2508  fcof=1d0
2509  ENDIF
2510  wdtp(i)=fac*fcof*(2d0-rm1-rm2-(rm1-rm2)**2)*
2511  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))
2512  wdtp(i)=fudge*wdtp(i)
2513  wdtp(0)=wdtp(0)+wdtp(i)
2514  IF(mdme(idc,1).GT.0) THEN
2515  wdte(i,mdme(idc,1))=wdtp(i)*wid2
2516  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
2517  wdte(i,0)=wdte(i,mdme(idc,1))
2518  wdte(0,0)=wdte(0,0)+wdte(i,0)
2519  ENDIF
2520  500 CONTINUE
2521 
2522  ELSEIF(kfla.EQ.9900041) THEN
2523 C...H_L++/--:
2524  fac=(1d0/(8d0*paru(1)))*shr
2525  DO 510 i=1,mdcy(kc,3)
2526  idc=i+mdcy(kc,2)-1
2527  IF(mdme(idc,1).LT.0) goto 510
2528  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
2529  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
2530  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 510
2531  wid2=1d0
2532  IF(i.LE.6) THEN
2533 C...H_L++/-- -> l+/- + l'+/-
2534  fcof=parp(180+3*((iabs(kfdp(idc,1))-11)/2)+
2535  & (iabs(kfdp(idc,2))-9)/2)**2
2536  IF(kfdp(idc,1).NE.kfdp(idc,2)) fcof=2d0*fcof
2537  ELSEIF(i.EQ.7) THEN
2538 C...H_L++/-- -> W_L+/- + W_L+/-
2539  fcof=0.5d0*parp(190)**4*parp(192)**2/pmas(24,1)**2*
2540  & (3d0*rm1+0.25d0/rm1-1d0)
2541  wid2=wids(24,4+(1-kfls)/2)
2542  ENDIF
2543  wdtp(i)=fac*fcof*
2544  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))
2545  wdtp(i)=fudge*wdtp(i)
2546  wdtp(0)=wdtp(0)+wdtp(i)
2547  IF(mdme(idc,1).GT.0) THEN
2548  wdte(i,mdme(idc,1))=wdtp(i)*wid2
2549  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
2550  wdte(i,0)=wdte(i,mdme(idc,1))
2551  wdte(0,0)=wdte(0,0)+wdte(i,0)
2552  ENDIF
2553  510 CONTINUE
2554 
2555  ELSEIF(kfla.EQ.9900042) THEN
2556 C...H_R++/--:
2557  fac=(1d0/(8d0*paru(1)))*shr
2558  DO 520 i=1,mdcy(kc,3)
2559  idc=i+mdcy(kc,2)-1
2560  IF(mdme(idc,1).LT.0) goto 520
2561  rm1=pmas(pycomp(kfdp(idc,1)),1)**2/sh
2562  rm2=pmas(pycomp(kfdp(idc,2)),1)**2/sh
2563  IF(sqrt(rm1)+sqrt(rm2).GT.1d0) goto 520
2564  wid2=1d0
2565  IF(i.LE.6) THEN
2566 C...H_R++/-- -> l+/- + l'+/-
2567  fcof=parp(180+3*((iabs(kfdp(idc,1))-11)/2)+
2568  & (iabs(kfdp(idc,2))-9)/2)**2
2569  IF(kfdp(idc,1).NE.kfdp(idc,2)) fcof=2d0*fcof
2570  ELSEIF(i.EQ.7) THEN
2571 C...H_R++/-- -> W_R+/- + W_R+/-
2572  fcof=parp(191)**2*(3d0*rm1+0.25d0/rm1-1d0)
2573  wid2=wids(pycomp(9900024),4+(1-kfls)/2)
2574  ENDIF
2575  wdtp(i)=fac*fcof*
2576  & sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))
2577  wdtp(i)=fudge*wdtp(i)
2578  wdtp(0)=wdtp(0)+wdtp(i)
2579  IF(mdme(idc,1).GT.0) THEN
2580  wdte(i,mdme(idc,1))=wdtp(i)*wid2
2581  wdte(0,mdme(idc,1))=wdte(0,mdme(idc,1))+wdte(i,mdme(idc,1))
2582  wdte(i,0)=wdte(i,mdme(idc,1))
2583  wdte(0,0)=wdte(0,0)+wdte(i,0)
2584  ENDIF
2585  520 CONTINUE
2586 
2587  ENDIF
2588  mint(61)=0
2589  mint(62)=0
2590  mint(63)=0
2591  RETURN
2592  END