Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyresd.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyresd.f
1 
2 
3 C*********************************************************************
4 
5 C...PYRESD
6 C...Allows resonances to decay (including parton showers for hadronic
7 C...channels).
8 
9  SUBROUTINE pyresd(IRES)
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...Parameter statement for maximum size of showers.
19  parameter(maxnur=1000)
20 C...Commonblocks.
21  common/pypart/npart,npartd,ipart(maxnur),ptpart(maxnur)
22  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
23  common/pyctag/nct,mct(4000,2)
24  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
25  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
26  common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
27  common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
28  common/pypars/mstp(200),parp(200),msti(200),pari(200)
29  common/pyint1/mint(400),vint(400)
30  common/pyint2/iset(500),kfpr(500,2),coef(500,20),icol(40,4,2)
31  common/pyint4/mwid(500),wids(500,5)
32  SAVE /pypart/,/pyjets/,/pyctag/,/pydat1/,/pydat2/,/pydat3/,
33  &/pysubs/,/pypars/,/pyint1/,/pyint2/,/pyint4/
34 C...Local arrays and complex and character variables.
35  dimension iref(50,8),kdcy(3),kfl1(3),kfl2(3),kfl3(3),keql(3),
36  &kcqm(3),kcq1(3),kcq2(3),kcq3(3),nsd(3),pmmn(3),ilin(6),
37  &hgz(3,3),coup(6,4),corl(2,2,2),pk(6,4),pkk(6,6),cthe(3),
38  &phi(3),wdtp(0:400),wdte(0:400,0:5),dpmo(5),xm(5),vdcy(4),
39  &itjunc(3),ctm2(3)
40  COMPLEX fgk,ha(6,6),hc(6,6)
41  REAL tir,uir
42  CHARACTER code*9,mass*9
43 
44 C...The F, Xi and Xj functions of Gunion and Kunszt
45 C...(Phys. Rev. D33, 665, plus errata from the authors).
46  fgk(i1,i2,i3,i4,i5,i6)=4.*ha(i1,i3)*hc(i2,i6)*(ha(i1,i5)*
47  &hc(i1,i4)+ha(i3,i5)*hc(i3,i4))
48  digk(dt,du)=-4d0*d34*d56+dt*(3d0*dt+4d0*du)+dt**2*(dt*du/
49  &(d34*d56)-2d0*(1d0/d34+1d0/d56)*(dt+du)+2d0*(d34/d56+d56/d34))
50  djgk(dt,du)=8d0*(d34+d56)**2-8d0*(d34+d56)*(dt+du)-6d0*dt*du-
51  &2d0*dt*du*(dt*du/(d34*d56)-2d0*(1d0/d34+1d0/d56)*(dt+du)+
52  &2d0*(d34/d56+d56/d34))
53 
54 C...Some general constants.
55  xw=paru(102)
56  xwv=xw
57  IF(mstp(8).GE.2) xw=1d0-(pmas(24,1)/pmas(23,1))**2
58  xw1=1d0-xw
59  sqmz=pmas(23,1)**2
60 
61  gmmz=pmas(23,1)*pmas(23,2)
62  sqmw=pmas(24,1)**2
63  gmmw=pmas(24,1)*pmas(24,2)
64  sh=vint(44)
65 
66 C...Boost and rotate to rest frame of incoming partons,
67 C...to get proper amount of smearing of decay angles.
68  ibst=0
69  IF(ires.EQ.0) THEN
70  ibst=1
71  etotin=p(mint(84)+1,4)+p(mint(84)+2,4)
72  bexin=(p(mint(84)+1,1)+p(mint(84)+2,1))/etotin
73  beyin=(p(mint(84)+1,2)+p(mint(84)+2,2))/etotin
74  bezin=(p(mint(84)+1,3)+p(mint(84)+2,3))/etotin
75  CALL pyrobo(mint(83)+7,n,0d0,0d0,-bexin,-beyin,-bezin)
76  phiin=pyangl(p(mint(84)+1,1),p(mint(84)+1,2))
77  CALL pyrobo(mint(83)+7,n,0d0,-phiin,0d0,0d0,0d0)
78  thein=pyangl(p(mint(84)+1,3),p(mint(84)+1,1))
79  CALL pyrobo(mint(83)+7,n,-thein,0d0,0d0,0d0,0d0)
80  ENDIF
81 
82 C...Reset original resonance configuration.
83  DO 100 jt=1,8
84  iref(1,jt)=0
85  100 CONTINUE
86 
87 C...Define initial one, two or three objects for subprocess.
88  ihdec=0
89  IF(ires.EQ.0) THEN
90  isub=mint(1)
91  IF(iset(isub).EQ.1.OR.iset(isub).EQ.3) THEN
92  iref(1,1)=mint(84)+2+iset(isub)
93  iref(1,4)=mint(83)+6+iset(isub)
94  jtmax=1
95  ELSEIF(iset(isub).EQ.2.OR.iset(isub).EQ.4) THEN
96  iref(1,1)=mint(84)+1+iset(isub)
97  iref(1,2)=mint(84)+2+iset(isub)
98  iref(1,4)=mint(83)+5+iset(isub)
99  iref(1,5)=mint(83)+6+iset(isub)
100  jtmax=2
101  ELSEIF(iset(isub).EQ.5) THEN
102  iref(1,1)=mint(84)+3
103  iref(1,2)=mint(84)+4
104  iref(1,3)=mint(84)+5
105  iref(1,4)=mint(83)+7
106  iref(1,5)=mint(83)+8
107  iref(1,6)=mint(83)+9
108  jtmax=3
109  ENDIF
110 
111 C...Define original resonance for odd cases.
112  ELSE
113  isub=0
114  IF(k(ires,2).EQ.25.OR.k(ires,2).EQ.35.OR.k(ires,2).EQ.36)
115  & ihdec=1
116  IF(ihdec.EQ.1) isub=3
117  iref(1,1)=ires
118  iref(1,4)=k(ires,3)
119  irestm=ires
120  IF(iref(1,4).GT.mint(84)) THEN
121  110 itmpmo=iref(1,4)
122  IF(k(itmpmo,2).EQ.94) THEN
123  iref(1,4)=k(itmpmo,3)+(irestm-itmpmo-1)
124  IF(k(iref(1,4),3).LE.mint(84)) iref(1,4)=k(iref(1,4),3)
125  ELSEIF(k(itmpmo,2).EQ.k(ires,2)) THEN
126  irestm=itmpmo
127 C...Explicitly check that reference particle exists, otherwise stop recursion
128  IF(itmpmo.GT.0.AND.k(itmpmo,3).GT.0) THEN
129  iref(1,4)=k(itmpmo,3)
130  goto 110
131  ENDIF
132  ENDIF
133  ENDIF
134  IF(iref(1,4).GT.mint(84)) THEN
135  ematch=1d10
136  iref14=iref(1,4)
137  DO 120 ii=mint(83)+7,mint(83)+mint(4)
138  IF(k(ii,2).EQ.k(ires,2).AND.abs(p(ii,4)-p(iref14,4)).LT.
139  & ematch) THEN
140  iref(1,4)=ii
141  ematch=abs(p(ii,4)-p(iref14,4))
142  ENDIF
143  120 CONTINUE
144  ENDIF
145  jtmax=1
146  ENDIF
147 
148 C...Check if initial resonance has been moved (in resonance + jet).
149  DO 140 jt=1,3
150  IF(iref(1,jt).GT.0) THEN
151  IF(k(iref(1,jt),1).GT.10) THEN
152  kfa=iabs(k(iref(1,jt),2))
153  IF(kfa.GE.6.AND.kchg(pycomp(kfa),2).NE.0) THEN
154  kda1=mod(k(iref(1,jt),4),mstu(5))
155  kda2=mod(k(iref(1,jt),5),mstu(5))
156  IF(kda1.GT.iref(1,jt).AND.kda1.LE.n) THEN
157  IF(k(kda1,2).EQ.21) kda1=k(kda1,5)/mstu(5)
158  ENDIF
159  IF(kda2.GT.iref(1,jt).AND.kda2.LE.n) THEN
160  IF(k(kda2,2).EQ.21) kda2=k(kda2,4)/mstu(5)
161  ENDIF
162  DO 130 i=iref(1,jt)+1,n
163  IF(k(i,2).EQ.k(iref(1,jt),2).AND.(i.EQ.kda1.OR.
164  & i.EQ.kda2)) THEN
165  iref(1,jt)=i
166  kda1=mod(k(iref(1,jt),4),mstu(5))
167  kda2=mod(k(iref(1,jt),5),mstu(5))
168  IF(kda1.GT.iref(1,jt).AND.kda1.LE.n) THEN
169  IF(k(kda1,2).EQ.21) kda1=k(kda1,5)/mstu(5)
170  ENDIF
171  IF(kda2.GT.iref(1,jt).AND.kda2.LE.n) THEN
172  IF(k(kda2,2).EQ.21) kda2=k(kda2,4)/mstu(5)
173  ENDIF
174  ENDIF
175  130 CONTINUE
176  ELSE
177  kda=mod(k(iref(1,jt),4),mstu(5))
178  IF(mwid(pycomp(kfa)).NE.0.AND.kda.GT.1) iref(1,jt)=kda
179  ENDIF
180  ENDIF
181  ENDIF
182  140 CONTINUE
183 
184 C...Set decay vertex for initial resonances
185  DO 160 jt=1,jtmax
186  DO 150 i=1,4
187  v(iref(1,jt),i)=0d0
188  150 CONTINUE
189  160 CONTINUE
190 
191 C...Loop over decay history.
192  np=1
193  ip=0
194  170 ip=ip+1
195  ninh=0
196  jtmax=2
197  IF(iref(ip,2).EQ.0) jtmax=1
198  IF(iref(ip,3).NE.0) jtmax=3
199  it4=0
200  nsav=n
201 
202 C...Check for Higgs which appears as decay product of user-process.
203  IF(isub.EQ.0) THEN
204  ihdec=0
205  IF(iref(ip,7).EQ.25.OR.iref(ip,7).EQ.35.OR.iref(ip,7)
206  & .EQ.36) ihdec=1
207  IF(ihdec.EQ.1) isub=3
208  ENDIF
209 
210 C...Start treatment of one, two or three resonances in parallel.
211  180 n=nsav
212  DO 340 jt=1,jtmax
213  id=iref(ip,jt)
214  kdcy(jt)=0
215  kfl1(jt)=0
216  kfl2(jt)=0
217  kfl3(jt)=0
218  keql(jt)=0
219  nsd(jt)=id
220  itjunc(jt)=0
221 
222 C...Check whether particle can/is allowed to decay.
223  IF(id.EQ.0) goto 330
224  kfa=iabs(k(id,2))
225  kca=pycomp(kfa)
226  IF(mwid(kca).EQ.0) goto 330
227  IF(k(id,1).GT.10.OR.mdcy(kca,1).EQ.0) goto 330
228  IF(kfa.EQ.6.OR.kfa.EQ.7.OR.kfa.EQ.8.OR.kfa.EQ.17.OR.
229  & kfa.EQ.18) it4=it4+1
230  k(id,4)=mstu(5)*(k(id,4)/mstu(5))
231  k(id,5)=mstu(5)*(k(id,5)/mstu(5))
232 
233 C...Choose lifetime and determine decay vertex.
234  IF(k(id,1).EQ.5) THEN
235  v(id,5)=0d0
236  ELSEIF(k(id,1).NE.4) THEN
237  v(id,5)=-pmas(kca,4)*log(pyr(0))
238  ENDIF
239  DO 190 j=1,4
240  vdcy(j)=v(id,j)+v(id,5)*p(id,j)/p(id,5)
241  190 CONTINUE
242 
243 C...Determine whether decay allowed or not.
244  mout=0
245  IF(mstj(22).EQ.2) THEN
246  IF(pmas(kca,4).GT.parj(71)) mout=1
247  ELSEIF(mstj(22).EQ.3) THEN
248  IF(vdcy(1)**2+vdcy(2)**2+vdcy(3)**2.GT.parj(72)**2) mout=1
249  ELSEIF(mstj(22).EQ.4) THEN
250  IF(vdcy(1)**2+vdcy(2)**2.GT.parj(73)**2) mout=1
251  IF(abs(vdcy(3)).GT.parj(74)) mout=1
252  ENDIF
253  IF(mout.EQ.1.AND.k(id,1).NE.5) THEN
254  k(id,1)=4
255  goto 330
256  ENDIF
257 
258 C...Info for selection of decay channel: sign, pairings.
259  IF(kchg(kca,3).EQ.0) THEN
260  ipm=2
261  ELSE
262  ipm=(5-isign(1,k(id,2)))/2
263  ENDIF
264  kfb=0
265  IF(jtmax.EQ.2) THEN
266  kfb=iabs(k(iref(ip,3-jt),2))
267  ELSEIF(jtmax.EQ.3) THEN
268  jt2=jt+1-3*(jt/3)
269  kfb=iabs(k(iref(ip,jt2),2))
270  IF(kfb.NE.kfa) THEN
271  jt2=jt+2-3*((jt+1)/3)
272  kfb=iabs(k(iref(ip,jt2),2))
273  ENDIF
274  ENDIF
275 
276 C...Select decay channel.
277  IF(isub.EQ.1.OR.isub.EQ.15.OR.isub.EQ.19.OR.isub.EQ.22.OR.
278  & isub.EQ.30.OR.isub.EQ.35.OR.isub.EQ.141) mint(61)=1
279  CALL pywidt(kfa,p(id,5)**2,wdtp,wdte)
280  wdte0s=wdte(0,1)+wdte(0,ipm)+wdte(0,4)
281  IF(kfb.EQ.kfa) wdte0s=wdte0s+wdte(0,5)
282  IF(wdte0s.LE.0d0) goto 330
283  rkfl=wdte0s*pyr(0)
284  idl=0
285  200 idl=idl+1
286  idc=idl+mdcy(kca,2)-1
287  rkfl=rkfl-(wdte(idl,1)+wdte(idl,ipm)+wdte(idl,4))
288  IF(kfb.EQ.kfa) rkfl=rkfl-wdte(idl,5)
289  IF(idl.LT.mdcy(kca,3).AND.rkfl.GT.0d0) goto 200
290 
291 C...Read out flavours and colour charges of decay channel chosen.
292  kcqm(jt)=kchg(kca,2)*isign(1,k(id,2))
293  IF(kcqm(jt).EQ.-2) kcqm(jt)=2
294  kfl1(jt)=kfdp(idc,1)*isign(1,k(id,2))
295  kfc1a=pycomp(iabs(kfl1(jt)))
296  IF(kchg(kfc1a,3).EQ.0) kfl1(jt)=iabs(kfl1(jt))
297  kcq1(jt)=kchg(kfc1a,2)*isign(1,kfl1(jt))
298  IF(kcq1(jt).EQ.-2) kcq1(jt)=2
299  kfl2(jt)=kfdp(idc,2)*isign(1,k(id,2))
300  kfc2a=pycomp(iabs(kfl2(jt)))
301  IF(kchg(kfc2a,3).EQ.0) kfl2(jt)=iabs(kfl2(jt))
302  kcq2(jt)=kchg(kfc2a,2)*isign(1,kfl2(jt))
303  IF(kcq2(jt).EQ.-2) kcq2(jt)=2
304  kfl3(jt)=kfdp(idc,3)*isign(1,k(id,2))
305  kcq3(jt)=0
306  IF(kfl3(jt).NE.0) THEN
307  kfc3a=pycomp(iabs(kfl3(jt)))
308  IF(kchg(kfc3a,3).EQ.0) kfl3(jt)=iabs(kfl3(jt))
309  kcq3(jt)=kchg(kfc3a,2)*isign(1,kfl3(jt))
310  IF(kcq3(jt).EQ.-2) kcq3(jt)=2
311  ENDIF
312 
313 C...Set/save further info on channel.
314  kdcy(jt)=1
315  IF(kfb.EQ.kfa) keql(jt)=mdme(idc,1)
316  nsd(jt)=n
317  hgz(jt,1)=vint(111)
318  hgz(jt,2)=vint(112)
319  hgz(jt,3)=vint(114)
320  jtz=jt
321 
322 C...Select masses; to begin with assume resonances narrow.
323  DO 220 i=1,3
324  p(n+i,5)=0d0
325  pmmn(i)=0d0
326  IF(i.EQ.1) THEN
327  kflw=iabs(kfl1(jt))
328  kcw=kfc1a
329  ELSEIF(i.EQ.2) THEN
330  kflw=iabs(kfl2(jt))
331  kcw=kfc2a
332  ELSEIF(i.EQ.3) THEN
333  IF(kfl3(jt).EQ.0) goto 220
334  kflw=iabs(kfl3(jt))
335  kcw=kfc3a
336  ENDIF
337  p(n+i,5)=pmas(kcw,1)
338 CMRENNA++
339 C...This prevents SUSY/t particles from becoming too light.
340  IF(kflw/ksusy1.EQ.1.OR.kflw/ksusy1.EQ.2) THEN
341  pmmn(i)=pmas(kcw,1)
342  DO 210 idc=mdcy(kcw,2),mdcy(kcw,2)+mdcy(kcw,3)-1
343  IF(mdme(idc,1).GT.0.AND.brat(idc).GT.1e-4) THEN
344  pmsum=pmas(pycomp(kfdp(idc,1)),1)+
345  & pmas(pycomp(kfdp(idc,2)),1)
346  IF(kfdp(idc,3).NE.0) pmsum=pmsum+
347  & pmas(pycomp(kfdp(idc,3)),1)
348  pmmn(i)=min(pmmn(i),pmsum)
349  ENDIF
350  210 CONTINUE
351 CMRENNA--
352  ELSEIF(kflw.EQ.6) THEN
353  pmmn(i)=pmas(24,1)+pmas(5,1)
354  ENDIF
355  220 CONTINUE
356 
357 C...Check which two out of three are widest.
358  iwid1=1
359  iwid2=2
360  pwid1=pmas(kfc1a,2)
361  pwid2=pmas(kfc2a,2)
362  kflw1=iabs(kfl1(jt))
363  kflw2=iabs(kfl2(jt))
364  IF(kfl3(jt).NE.0) THEN
365  pwid3=pmas(kfc3a,2)
366  IF(pwid3.GT.pwid1.AND.pwid2.GE.pwid1) THEN
367  iwid1=3
368  pwid1=pwid3
369  kflw1=iabs(kfl3(jt))
370  ELSEIF(pwid3.GT.pwid2) THEN
371  iwid2=3
372  pwid2=pwid3
373  kflw2=iabs(kfl3(jt))
374  ENDIF
375  ENDIF
376 
377 C...If all narrow then only check that masses consistent.
378  IF(mstp(42).LE.0.OR.(pwid1.LT.parp(41).AND.
379  & pwid2.LT.parp(41))) THEN
380 CMRENNA++
381 C....Handle near degeneracy cases.
382  IF(kfa/ksusy1.EQ.1.OR.kfa/ksusy1.EQ.2) THEN
383  IF(p(n+1,5)+p(n+2,5)+p(n+3,5).GT.p(id,5)) THEN
384  p(n+1,5)=p(id,5)-p(n+2,5)-0.5d0
385  IF(p(n+1,5).LT.0d0) p(n+1,5)=0d0
386  ENDIF
387  ENDIF
388 CMRENNA--
389  IF(p(n+1,5)+p(n+2,5)+p(n+3,5).GT.p(id,5)) THEN
390  CALL pyerrm(13,'(PYRESD:) daughter masses too large')
391  mint(51)=1
392  goto 720
393  ELSEIF(p(n+1,5)+p(n+2,5)+p(n+3,5)+parj(64).GT.p(id,5)) THEN
394  CALL pyerrm(3,'(PYRESD:) daughter masses too large')
395  mint(51)=1
396  goto 720
397  ENDIF
398 
399 C...For three wide resonances select narrower of three
400 C...according to BW decoupled from rest.
401  ELSE
402  pmtot=p(id,5)
403  IF(kfl3(jt).NE.0) THEN
404  iwid3=6-iwid1-iwid2
405  kflw3=iabs(kfl1(jt))+iabs(kfl2(jt))+iabs(kfl3(jt))-
406  & kflw1-kflw2
407  loop=0
408  230 loop=loop+1
409  p(n+iwid3,5)=pymass(kflw3)
410  IF(loop.LE.10.AND. p(n+iwid3,5).LE.pmmn(iwid3)) goto 230
411  pmtot=pmtot-p(n+iwid3,5)
412  ENDIF
413 C...Select other two correlated within remaining phase space.
414  IF(ip.EQ.1) THEN
415  ckin45=ckin(45)
416  ckin47=ckin(47)
417  ckin(45)=max(pmmn(iwid1),ckin(45))
418  ckin(47)=max(pmmn(iwid2),ckin(47))
419  CALL pyofsh(2,kfa,kflw1,kflw2,pmtot,p(n+iwid1,5),
420  & p(n+iwid2,5))
421  ckin(45)=ckin45
422  ckin(47)=ckin47
423  ELSE
424  ckin(49)=pmmn(iwid1)
425  ckin(50)=pmmn(iwid2)
426  CALL pyofsh(5,kfa,kflw1,kflw2,pmtot,p(n+iwid1,5),
427  & p(n+iwid2,5))
428  ckin(49)=0d0
429  ckin(50)=0d0
430  ENDIF
431  IF(mint(51).EQ.1) goto 720
432  ENDIF
433 
434 C...Begin fill decay products, with colour flow for coloured objects.
435  mstu10=mstu(10)
436  mstu(10)=1
437  mstu(19)=1
438 
439 CMRENNA++
440 C...1) Three-body decays of SUSY particles (plus special case top).
441  IF(kfl3(jt).NE.0) THEN
442  DO 250 i=n+1,n+3
443  DO 240 j=1,5
444  k(i,j)=0
445  v(i,j)=0d0
446  240 CONTINUE
447  mct(i,1)=0
448  mct(i,2)=0
449  250 CONTINUE
450  k(n+1,1)=1
451  k(n+1,2)=kfl1(jt)
452  k(n+2,1)=1
453  k(n+2,2)=kfl2(jt)
454  k(n+3,1)=1
455  k(n+3,2)=kfl3(jt)
456  idin=id
457  CALL pytbdy(idin)
458 
459 C...Set colour flow for t -> W + b + Z.
460  IF(kfa.EQ.6) THEN
461  k(n+2,1)=3
462  isid=4
463  IF(kcqm(jt).EQ.-1) isid=5
464  idau=n+2
465  k(id,isid)=k(id,isid)+idau
466  k(idau,isid)=mstu(5)*id
467 
468 C...Set colour flow in three-body decays - programmed as special cases.
469 
470  ELSEIF(kfc2a.LE.6) THEN
471  k(n+2,1)=3
472  k(n+3,1)=3
473  isid=4
474  IF(kfl2(jt).LT.0) isid=5
475  k(n+2,isid)=mstu(5)*(n+3)
476  k(n+3,9-isid)=mstu(5)*(n+2)
477 C...PS++: Bugfix 16 MAR 2006 for 3-body squark decays (e.g. via SLHA)
478  ELSEIF(kfa.GT.ksusy1.AND.mod(kfa,ksusy1).LT.10
479  & .AND.kfl3(jt).NE.0) THEN
480  kqsuma=iabs(kcq1(jt))+iabs(kcq2(jt))+iabs(kcq3(jt))
481 C...3-body decays of squarks to colour singlets plus one quark
482  IF (kqsuma.EQ.1) THEN
483 C...Find quark
484  iq=0
485  IF (kcq1(jt).NE.0) iq=1
486  IF (kcq2(jt).NE.0) iq=2
487  IF (kcq3(jt).NE.0) iq=3
488  isid=4
489  IF (k(n+iq,2).LT.0) isid=5
490  k(n+iq,1)=3
491  k(id,isid)=k(id,isid)+(n+iq)
492  k(n+iq,isid)=mstu(5)*id
493  ENDIF
494 C...PS--
495  ENDIF
496  IF(kfl1(jt).EQ.ksusy1+21) THEN
497  k(n+1,1)=3
498  k(n+2,1)=3
499  k(n+3,1)=3
500  isid=4
501  IF(kfl2(jt).LT.0) isid=5
502  k(n+1,isid)=mstu(5)*(n+2)
503  k(n+1,9-isid)=mstu(5)*(n+3)
504  k(n+2,isid)=mstu(5)*(n+1)
505  k(n+3,9-isid)=mstu(5)*(n+1)
506  ENDIF
507  IF(kfa.EQ.ksusy1+21) THEN
508  k(n+2,1)=3
509  k(n+3,1)=3
510  isid=4
511  IF(kfl2(jt).LT.0) isid=5
512  k(id,isid)=k(id,isid)+(n+2)
513  k(id,9-isid)=k(id,9-isid)+(n+3)
514  k(n+2,isid)=mstu(5)*id
515  k(n+3,9-isid)=mstu(5)*id
516  ENDIF
517  nsav=n
518  n=n+3
519  n=nsav
520 CMRENNA--
521 
522  IF(kfa.GE.ksusy1+22.AND.kfa.LE.ksusy1+37.AND.
523  & iabs(kcq2(jt)).EQ.1) THEN
524  k(n+2,1)=3
525  k(n+3,1)=3
526  isid=4
527  IF(kfl2(jt).LT.0) isid=5
528  k(n+2,isid)=mstu(5)*(n+3)
529  k(n+3,9-isid)=mstu(5)*(n+2)
530  ENDIF
531 
532 C...Set colour flow in three-body decays with baryon number violation.
533 C...Neutralino and chargino decays first.
534  kcqsum=kcq1(jt)+kcq2(jt)+kcq3(jt)
535  IF(kcqm(jt).EQ.0.AND.iabs(kcqsum).EQ.3) THEN
536  itjunc(jt)=(1+(1-kcq1(jt))/2)
537  k(n+4,4)=itjunc(jt)*mstu(5)
538 C...Insert junction to keep track of colours.
539  IF(kcq1(jt).NE.0) k(n+1,1)=3
540  IF(kcq2(jt).NE.0) k(n+2,1)=3
541  IF(kcq3(jt).NE.0) k(n+3,1)=3
542 C...Set special junction codes:
543  k(n+4,1)=42
544  k(n+4,2)=88
545 
546 C...Order decay products by invariant mass. (will be used in PYSTRF).
547  pm12=p(n+1,4)*p(n+2,4)-p(n+1,1)*p(n+2,1)-p(n+1,2)*p(n+2,2)-
548  & p(n+1,3)*p(n+2,3)
549  pm13=p(n+1,4)*p(n+3,4)-p(n+1,1)*p(n+3,1)-p(n+1,2)*p(n+3,2)-
550  & p(n+1,3)*p(n+3,3)
551  pm23=p(n+2,4)*p(n+3,4)-p(n+2,1)*p(n+3,1)-p(n+2,2)*p(n+3,2)-
552  & p(n+2,3)*p(n+3,3)
553  IF(pm12.LT.pm13.AND.pm12.LT.pm23) THEN
554  k(n+4,4)=n+3+k(n+4,4)
555  k(n+4,5)=n+1+mstu(5)*(n+2)
556  ELSEIF(pm13.LT.pm23) THEN
557  k(n+4,4)=n+2+k(n+4,4)
558  k(n+4,5)=n+1+mstu(5)*(n+3)
559  ELSE
560  k(n+4,4)=n+1+k(n+4,4)
561  k(n+4,5)=n+2+mstu(5)*(n+3)
562  ENDIF
563  DO 260 j=1,5
564  p(n+4,j)=0d0
565  v(n+4,j)=0d0
566  260 CONTINUE
567 C...Connect daughters to junction.
568  DO 270 ii=n+1,n+3
569  k(ii,4)=0
570  k(ii,5)=0
571  k(ii,itjunc(jt)+3)=mstu(5)*(n+4)
572  270 CONTINUE
573 C...Particle counter should be stepped up one extra for junction.
574  n=n+1
575 
576 C...Gluino decays.
577  ELSEIF (kcqm(jt).EQ.2.AND.iabs(kcqsum).EQ.3) THEN
578  itjunc(jt)=(5+(1-kcq1(jt))/2)
579  k(n+4,4)=itjunc(jt)*mstu(5)
580 C...Insert junction to keep track of colours.
581  IF(kcq1(jt).NE.0) k(n+1,1)=3
582  IF(kcq2(jt).NE.0) k(n+2,1)=3
583  IF(kcq3(jt).NE.0) k(n+3,1)=3
584  k(n+4,1)=42
585  k(n+4,2)=88
586  DO 280 j=1,5
587  p(n+4,j)=0d0
588  v(n+4,j)=0d0
589  280 CONTINUE
590  ctmsum=0d0
591  DO 290 ii=n+1,n+3
592  k(ii,4)=0
593  k(ii,5)=0
594 C...Start by connecting all daughters to junction.
595  k(ii,itjunc(jt)-1)=mstu(5)*(n+4)
596 C...Only consider colour topologies with off shell resonances.
597  rmq1=pmas(pycomp(k(ii,2)),1)
598  rmres=pmas(pycomp(ksusy1+iabs(k(ii,2))),1)
599  rmglu=pmas(pycomp(ksusy1+21),1)
600  IF (rmglu-rmq1.LT.rmres) THEN
601 C...Calculate propagators for each colour topology.
602  rm2q23=rmglu**2+rmq1**2-2d0*(p(ii,4)*p(id,4)+p(ii,1)
603  & *p(id,1)+p(ii,2)*p(id,2)+p(ii,3)*p(id,3))
604  ctm2(ii-n)=1d0/(rm2q23-rmres**2)**2
605  ELSE
606  ctm2(ii-n)=0d0
607  ENDIF
608  ctmsum=ctmsum+ctm2(ii-n)
609  290 CONTINUE
610  ctmsum=pyr(0)*ctmsum
611 C...Select colour topology J, with most off shell least likely.
612  j=0
613  300 j=j+1
614  ctmsum=ctmsum-ctm2(j)
615  IF (ctmsum.GT.0d0) goto 300
616 C...The lucky winner gets its colour (anti-colour) directly from gluino.
617  k(n+j,itjunc(jt)-1)=mstu(5)*id
618  k(id,itjunc(jt)-1)=n+j+(k(id,itjunc(jt)-1)/mstu(5))*mstu(5)
619 C...The other gluino colour is connected to junction
620  k(id,10-itjunc(jt))=n+4+(k(id,10-itjunc(jt))/mstu(5))*
621  & mstu(5)
622  k(n+4,4)=k(n+4,4)+id
623 C...Lastly, connect junction to remaining daughters.
624  k(n+4,5)=n+1+mod(j,3)+mstu(5)*(n+1+mod(j+1,3))
625 C...Particle counter should be stepped up one extra for junction.
626  n=n+1
627  ENDIF
628 
629 C...Update particle counter.
630  n=n+3
631 
632 C...2) Everything else two-body decay.
633  ELSE
634  CALL py2ent(n+1,kfl1(jt),kfl2(jt),p(id,5))
635  mct(n-1,1)=0
636  mct(n-1,2)=0
637  mct(n,1)=0
638  mct(n,2)=0
639 C...First set colour flow as if mother colour singlet.
640  IF(kcq1(jt).NE.0) THEN
641  k(n-1,1)=3
642  IF(kcq1(jt).NE.-1) k(n-1,4)=mstu(5)*n
643  IF(kcq1(jt).NE.1) k(n-1,5)=mstu(5)*n
644  ENDIF
645  IF(kcq2(jt).NE.0) THEN
646  k(n,1)=3
647  IF(kcq2(jt).NE.-1) k(n,4)=mstu(5)*(n-1)
648  IF(kcq2(jt).NE.1) k(n,5)=mstu(5)*(n-1)
649  ENDIF
650 C...Then redirect colour flow if mother (anti)triplet.
651  IF(kcqm(jt).EQ.0) THEN
652  ELSEIF(kcqm(jt).NE.2) THEN
653  isid=4
654  IF(kcqm(jt).EQ.-1) isid=5
655  idau=n-1
656  IF(kcq1(jt).EQ.0.OR.kcq2(jt).EQ.2) idau=n
657  k(id,isid)=k(id,isid)+idau
658  k(idau,isid)=mstu(5)*id
659 C...Then redirect colour flow if mother octet.
660  ELSEIF(kcq1(jt).EQ.0.OR.kcq2(jt).EQ.0) THEN
661  idau=n-1
662  IF(kcq1(jt).EQ.0) idau=n
663  k(id,4)=k(id,4)+idau
664  k(id,5)=k(id,5)+idau
665  k(idau,4)=mstu(5)*id
666  k(idau,5)=mstu(5)*id
667  ELSE
668  isid=4
669  IF(kcq1(jt).EQ.-1) isid=5
670  IF(kcq1(jt).EQ.2) isid=int(4.5d0+pyr(0))
671  k(id,isid)=k(id,isid)+(n-1)
672  k(id,9-isid)=k(id,9-isid)+n
673  k(n-1,isid)=mstu(5)*id
674  k(n,9-isid)=mstu(5)*id
675  ENDIF
676 
677 C...Insert junction
678  IF(iabs(kcq1(jt)+kcq2(jt)-kcqm(jt)).EQ.3) THEN
679  n=n+1
680 C...~q* mother: type 3 junction. ~q mother: type 4.
681  itjunc(jt)=(7+kcqm(jt))/2
682 C...Specify junction KF and set colour flow from junction
683  k(n,1)=42
684  k(n,2)=88
685  k(n,3)=id
686 C...Junction type encoded together with mother:
687  k(n,4)=id+itjunc(jt)*mstu(5)
688  k(n,5)=n-1+mstu(5)*(n-2)
689 C...Zero P and V for junction (V filled later)
690  DO 310 j=1,5
691  p(n,j)=0d0
692  v(n,j)=0d0
693  310 CONTINUE
694 C...Set colour flow from mother to junction
695  k(id,8-itjunc(jt))= n + mstu(5)*(k(id,8-itjunc(jt))/mstu(5))
696 C...Set colour flow from daughters to junction
697  DO 320 ii=n-2,n-1
698  k(ii,4) = 0
699  k(ii,5) = 0
700 C...(Anti-)colour mother is junction.
701  k(ii,1+itjunc(jt)) = mstu(5)*(n)
702  320 CONTINUE
703  ENDIF
704  ENDIF
705 
706 C...End loop over resonances for daughter flavour and mass selection.
707  mstu(10)=mstu10
708  330 IF(mwid(kca).NE.0.AND.(kfl1(jt).EQ.0.OR.kfl3(jt).NE.0))
709  & ninh=ninh+1
710  IF(ires.GT.0.AND.mwid(kca).NE.0.AND.mdcy(kca,1).NE.0.AND.
711  & kfl1(jt).EQ.0) THEN
712  WRITE(code,'(I9)') k(id,2)
713  WRITE(mass,'(F9.3)') p(id,5)
714  CALL pyerrm(3,'(PYRESD:) Failed to decay particle'//
715  & code//' with mass'//mass)
716  mint(51)=1
717  goto 720
718  ENDIF
719  340 CONTINUE
720 
721 C...Check for allowed combinations. Skip if no decays.
722  IF(jtmax.EQ.1) THEN
723  IF(kdcy(1).EQ.0) goto 710
724  ELSEIF(jtmax.EQ.2) THEN
725  IF(kdcy(1).EQ.0.AND.kdcy(2).EQ.0) goto 710
726  IF(keql(1).EQ.4.AND.keql(2).EQ.4) goto 180
727  IF(keql(1).EQ.5.AND.keql(2).EQ.5) goto 180
728  ELSEIF(jtmax.EQ.3) THEN
729  IF(kdcy(1).EQ.0.AND.kdcy(2).EQ.0.AND.kdcy(3).EQ.0) goto 710
730  IF(keql(1).EQ.4.AND.keql(2).EQ.4) goto 180
731  IF(keql(1).EQ.4.AND.keql(3).EQ.4) goto 180
732  IF(keql(2).EQ.4.AND.keql(3).EQ.4) goto 180
733  IF(keql(1).EQ.5.AND.keql(2).EQ.5) goto 180
734  IF(keql(1).EQ.5.AND.keql(3).EQ.5) goto 180
735  IF(keql(2).EQ.5.AND.keql(3).EQ.5) goto 180
736  ENDIF
737 
738 C...Special case: matrix element option for Z0 decay to quarks.
739  IF(mstp(48).EQ.1.AND.isub.EQ.1.AND.jtmax.EQ.1.AND.
740  &iabs(mint(11)).EQ.11.AND.iabs(kfl1(1)).LE.5) THEN
741 
742 C...Check consistency of MSTJ options set.
743  IF(mstj(109).EQ.2.AND.mstj(110).NE.1) THEN
744  CALL pyerrm(6,
745  & '(PYRESD:) MSTJ(109) value requires MSTJ(110) = 1')
746  mstj(110)=1
747  ENDIF
748  IF(mstj(109).EQ.2.AND.mstj(111).NE.0) THEN
749  CALL pyerrm(6,
750  & '(PYRESD:) MSTJ(109) value requires MSTJ(111) = 0')
751 
752  mstj(111)=0
753  ENDIF
754 
755 C...Select alpha_strong behaviour.
756  mst111=mstu(111)
757  par112=paru(112)
758  mstu(111)=mstj(108)
759  IF(mstj(108).EQ.2.AND.(mstj(101).EQ.0.OR.mstj(101).EQ.1))
760  & mstu(111)=1
761  paru(112)=parj(121)
762  IF(mstu(111).EQ.2) paru(112)=parj(122)
763 
764 C...Find axial fraction in total cross section for scalar gluon model.
765  parj(171)=0d0
766  IF((iabs(mstj(101)).EQ.1.AND.mstj(109).EQ.1).OR.
767  & (mstj(101).EQ.5.AND.mstj(49).EQ.1)) THEN
768  poll=1d0-parj(131)*parj(132)
769  sff=1d0/(16d0*xw*xw1)
770  sfw=p(id,5)**4/((p(id,5)**2-parj(123)**2)**2+
771  & (parj(123)*parj(124))**2)
772  sfi=sfw*(1d0-(parj(123)/p(id,5))**2)
773  ve=4d0*xw-1d0
774  hf1i=sfi*sff*(ve*poll+parj(132)-parj(131))
775  hf1w=sfw*sff**2*((ve**2+1d0)*poll+2d0*ve*
776  & (parj(132)-parj(131)))
777  kflc=iabs(kfl1(1))
778  pmq=pymass(kflc)
779  qf=kchg(kflc,1)/3d0
780  vq=1d0
781  IF(mod(mstj(103),2).EQ.1) vq=sqrt(max(0d0,
782  & 1d0-(2d0*pmq/p(id,5))**2))
783  vf=sign(1d0,qf)-4d0*qf*xw
784  rfv=0.5d0*vq*(3d0-vq**2)*(qf**2*poll-2d0*qf*vf*hf1i+
785  & vf**2*hf1w)+vq**3*hf1w
786  IF(rfv.GT.0d0) parj(171)=min(1d0,vq**3*hf1w/rfv)
787  ENDIF
788 
789 C...Choice of jet configuration.
790  CALL pyxjet(p(id,5),njet,cut)
791  kflc=iabs(kfl1(1))
792  kfln=21
793  IF(njet.EQ.4) THEN
794  CALL pyx4jt(njet,cut,kflc,p(id,5),kfln,x1,x2,x4,x12,x14)
795  ELSEIF(njet.EQ.3) THEN
796  CALL pyx3jt(njet,cut,kflc,p(id,5),x1,x3)
797  ELSE
798  mstj(120)=1
799  ENDIF
800 
801 C...Fill jet configuration; return if incorrect kinematics.
802  nc=n-2
803  IF(njet.EQ.2.AND.mstj(101).NE.5) THEN
804  CALL py2ent(nc+1,kflc,-kflc,p(id,5))
805  ELSEIF(njet.EQ.2) THEN
806  CALL py2ent(-(nc+1),kflc,-kflc,p(id,5))
807  ELSEIF(njet.EQ.3) THEN
808  CALL py3ent(nc+1,kflc,21,-kflc,p(id,5),x1,x3)
809  ELSEIF(kfln.EQ.21) THEN
810  CALL py4ent(nc+1,kflc,kfln,kfln,-kflc,p(id,5),x1,x2,x4,
811  & x12,x14)
812  ELSE
813  CALL py4ent(nc+1,kflc,-kfln,kfln,-kflc,p(id,5),x1,x2,x4,
814  & x12,x14)
815  ENDIF
816  IF(mstu(24).NE.0) THEN
817  mint(51)=1
818  mstu(111)=mst111
819  paru(112)=par112
820  goto 720
821  ENDIF
822 
823 C...Angular orientation according to matrix element.
824  IF(mstj(106).EQ.1) THEN
825  CALL pyxdif(nc,njet,kflc,p(id,5),chiz,thez,phiz)
826  IF(mint(11).LT.0) thez=paru(1)-thez
827  cthe(1)=cos(thez)
828  CALL pyrobo(nc+1,n,0d0,chiz,0d0,0d0,0d0)
829  CALL pyrobo(nc+1,n,thez,phiz,0d0,0d0,0d0)
830  ENDIF
831 
832 C...Boost partons to Z0 rest frame.
833  CALL pyrobo(nc+1,n,0d0,0d0,p(id,1)/p(id,4),
834  & p(id,2)/p(id,4),p(id,3)/p(id,4))
835 
836 C...Mark decayed resonance and add documentation lines,
837  k(id,1)=k(id,1)+10
838  idoc=mint(83)+mint(4)
839  DO 360 i=nc+1,n
840  i1=mint(83)+mint(4)+1
841  k(i,3)=i1
842  IF(mstp(128).GE.1) k(i,3)=id
843  IF(mstp(128).LE.1.AND.mint(4).LT.mstp(126)) THEN
844  mint(4)=mint(4)+1
845  k(i1,1)=21
846  k(i1,2)=k(i,2)
847  k(i1,3)=iref(ip,4)
848  DO 350 j=1,5
849  p(i1,j)=p(i,j)
850  350 CONTINUE
851  ENDIF
852  360 CONTINUE
853 
854 C...Generate parton shower.
855  IF(mstj(101).EQ.5.AND.mint(35).LE.1) THEN
856  CALL pyshow(n-1,n,p(id,5))
857  ELSEIF(mstj(101).EQ.5.AND.mint(35).GE.2) THEN
858  npart=2
859  ipart(1)=n-1
860  ipart(2)=n
861  ptpart(1)=0.5d0*p(id,5)
862  ptpart(2)=ptpart(1)
863  nct=nct+1
864  IF(k(n-1,2).GT.0) THEN
865  mct(n-1,1)=nct
866  mct(n,2)=nct
867  ELSE
868  mct(n-1,2)=nct
869  mct(n,1)=nct
870  ENDIF
871  CALL pyptfs(2,0.5d0*p(id,5),0d0,ptgen)
872  ENDIF
873 
874 C... End special case for Z0: skip ahead.
875  mstu(111)=mst111
876  paru(112)=par112
877  goto 700
878  ENDIF
879 
880 C...Order incoming partons and outgoing resonances.
881  IF(jtmax.EQ.2.AND.isub.NE.0.AND.mstp(47).GE.1.AND.
882  &ninh.EQ.0) THEN
883  ilin(1)=mint(84)+1
884  IF(k(mint(84)+1,2).GT.0) ilin(1)=mint(84)+2
885  IF(k(ilin(1),2).EQ.21.OR.k(ilin(1),2).EQ.22)
886  & ilin(1)=2*mint(84)+3-ilin(1)
887  ilin(2)=2*mint(84)+3-ilin(1)
888  imin=1
889  IF(iref(ip,7).EQ.25.OR.iref(ip,7).EQ.35.OR.iref(ip,7)
890  & .EQ.36) imin=3
891  imax=2
892  iord=1
893  IF(k(iref(ip,1),2).EQ.23) iord=2
894  IF(k(iref(ip,1),2).EQ.24.AND.k(iref(ip,2),2).EQ.-24) iord=2
895  iakipd=iabs(k(iref(ip,iord),2))
896  IF(iakipd.EQ.25.OR.iakipd.EQ.35.OR.iakipd.EQ.36) iord=3-iord
897  IF(kdcy(iord).EQ.0) iord=3-iord
898 
899 C...Order decay products of resonances.
900  DO 370 jt=iord,3-iord,3-2*iord
901  IF(kdcy(jt).EQ.0) THEN
902  ilin(imax+1)=nsd(jt)
903  imax=imax+1
904  ELSEIF(k(nsd(jt)+1,2).GT.0) THEN
905  ilin(imax+1)=n+2*jt-1
906  ilin(imax+2)=n+2*jt
907  imax=imax+2
908  k(n+2*jt-1,2)=k(nsd(jt)+1,2)
909  k(n+2*jt,2)=k(nsd(jt)+2,2)
910  ELSE
911  ilin(imax+1)=n+2*jt
912 
913  ilin(imax+2)=n+2*jt-1
914  imax=imax+2
915  k(n+2*jt-1,2)=k(nsd(jt)+1,2)
916  k(n+2*jt,2)=k(nsd(jt)+2,2)
917  ENDIF
918  370 CONTINUE
919 
920 C...Find charge, isospin, left- and righthanded couplings.
921  DO 390 i=imin,imax
922  DO 380 j=1,4
923  coup(i,j)=0d0
924  380 CONTINUE
925  kfa=iabs(k(ilin(i),2))
926  IF(kfa.EQ.0.OR.kfa.GT.20) goto 390
927  coup(i,1)=kchg(kfa,1)/3d0
928  coup(i,2)=(-1)**mod(kfa,2)
929  coup(i,4)=-2d0*coup(i,1)*xwv
930  coup(i,3)=coup(i,2)+coup(i,4)
931  390 CONTINUE
932 
933 C...Full propagator dependence and flavour correlations for 2 gamma*/Z.
934  IF(isub.EQ.22) THEN
935  DO 420 i=3,5,2
936  i1=iord
937  IF(i.EQ.5) i1=3-iord
938  DO 410 j1=1,2
939  DO 400 j2=1,2
940  corl(i/2,j1,j2)=coup(1,1)**2*hgz(i1,1)*coup(i,1)**2/
941  & 16d0+coup(1,1)*coup(1,j1+2)*hgz(i1,2)*coup(i,1)*
942  & coup(i,j2+2)/4d0+coup(1,j1+2)**2*hgz(i1,3)*
943  & coup(i,j2+2)**2
944  400 CONTINUE
945  410 CONTINUE
946  420 CONTINUE
947  cowt12=(corl(1,1,1)+corl(1,1,2))*(corl(2,1,1)+corl(2,1,2))+
948  & (corl(1,2,1)+corl(1,2,2))*(corl(2,2,1)+corl(2,2,2))
949  comx12=(corl(1,1,1)+corl(1,1,2)+corl(1,2,1)+corl(1,2,2))*
950  & (corl(2,1,1)+corl(2,1,2)+corl(2,2,1)+corl(2,2,2))
951 
952  IF(cowt12.LT.pyr(0)*comx12) goto 180
953  ENDIF
954  ENDIF
955 
956 C...Select angular orientation type - Z'/W' only.
957  mzpwp=0
958  IF(isub.EQ.141) THEN
959  IF(pyr(0).LT.paru(130)) mzpwp=1
960  IF(ip.EQ.2) THEN
961  IF(iabs(k(iref(2,1),2)).EQ.37) mzpwp=2
962  iakir=iabs(k(iref(2,2),2))
963  IF(iakir.EQ.25.OR.iakir.EQ.35.OR.iakir.EQ.36) mzpwp=2
964  IF(iakir.LE.20) mzpwp=2
965  ENDIF
966  IF(ip.GE.3) mzpwp=2
967  ELSEIF(isub.EQ.142) THEN
968  IF(pyr(0).LT.paru(136)) mzpwp=1
969  IF(ip.EQ.2) THEN
970  iakir=iabs(k(iref(2,2),2))
971  IF(iakir.EQ.25.OR.iakir.EQ.35.OR.iakir.EQ.36) mzpwp=2
972  IF(iakir.LE.20) mzpwp=2
973  ENDIF
974  IF(ip.GE.3) mzpwp=2
975  ENDIF
976 
977 C...Select random angles (begin of weighting procedure).
978  430 DO 440 jt=1,jtmax
979  IF(kdcy(jt).EQ.0) goto 440
980  IF(jtmax.EQ.1.AND.isub.NE.0.AND.ihdec.EQ.0) THEN
981  cthe(jt)=vint(13)+(vint(33)-vint(13)+vint(34)-vint(14))*pyr(0)
982  IF(cthe(jt).GT.vint(33)) cthe(jt)=cthe(jt)+vint(14)-vint(33)
983  phi(jt)=vint(24)
984  ELSE
985  cthe(jt)=2d0*pyr(0)-1d0
986  phi(jt)=paru(2)*pyr(0)
987  ENDIF
988  440 CONTINUE
989 
990  IF(jtmax.EQ.2.AND.mstp(47).GE.1.AND.ninh.EQ.0) THEN
991 C...Construct massless four-vectors.
992  DO 460 i=n+1,n+4
993  k(i,1)=1
994  DO 450 j=1,5
995  p(i,j)=0d0
996  v(i,j)=0d0
997  450 CONTINUE
998  460 CONTINUE
999  DO 470 jt=1,jtmax
1000  IF(kdcy(jt).EQ.0) goto 470
1001  id=iref(ip,jt)
1002  p(n+2*jt-1,3)=0.5d0*p(id,5)
1003  p(n+2*jt-1,4)=0.5d0*p(id,5)
1004  p(n+2*jt,3)=-0.5d0*p(id,5)
1005  p(n+2*jt,4)=0.5d0*p(id,5)
1006  CALL pyrobo(n+2*jt-1,n+2*jt,acos(cthe(jt)),phi(jt),
1007  & p(id,1)/p(id,4),p(id,2)/p(id,4),p(id,3)/p(id,4))
1008  470 CONTINUE
1009 
1010 C...Store incoming and outgoing momenta, with random rotation to
1011 C...avoid accidental zeroes in HA expressions.
1012  IF(isub.NE.0) THEN
1013  DO 490 i=imin,imax
1014  k(n+4+i,1)=1
1015  p(n+4+i,4)=sqrt(p(ilin(i),1)**2+p(ilin(i),2)**2+
1016  & p(ilin(i),3)**2+p(ilin(i),5)**2)
1017  p(n+4+i,5)=p(ilin(i),5)
1018  DO 480 j=1,3
1019  p(n+4+i,j)=p(ilin(i),j)
1020  480 CONTINUE
1021  490 CONTINUE
1022  500 therr=acos(2d0*pyr(0)-1d0)
1023  phirr=paru(2)*pyr(0)
1024  CALL pyrobo(n+4+imin,n+4+imax,therr,phirr,0d0,0d0,0d0)
1025  DO 520 i=imin,imax
1026  IF(p(n+4+i,1)**2+p(n+4+i,2)**2.LT.1d-4*(p(n+4+i,1)**2+
1027  & p(n+4+i,2)**2+p(n+4+i,3)**2)) goto 500
1028  DO 510 j=1,4
1029  pk(i,j)=p(n+4+i,j)
1030  510 CONTINUE
1031  520 CONTINUE
1032  ENDIF
1033 
1034 C...Calculate internal products.
1035  IF(isub.EQ.22.OR.isub.EQ.23.OR.isub.EQ.25.OR.isub.EQ.141.OR.
1036  & isub.EQ.142) THEN
1037  DO 540 i1=imin,imax-1
1038  DO 530 i2=i1+1,imax
1039  ha(i1,i2)=sngl(sqrt((pk(i1,4)-pk(i1,3))*(pk(i2,4)+
1040  & pk(i2,3))/(1d-20+pk(i1,1)**2+pk(i1,2)**2)))*
1041  & cmplx(sngl(pk(i1,1)),sngl(pk(i1,2)))-
1042  & sngl(sqrt((pk(i1,4)+pk(i1,3))*(pk(i2,4)-pk(i2,3))/
1043  & (1d-20+pk(i2,1)**2+pk(i2,2)**2)))*
1044  & cmplx(sngl(pk(i2,1)),sngl(pk(i2,2)))
1045  hc(i1,i2)=conjg(ha(i1,i2))
1046  IF(i1.LE.2) ha(i1,i2)=cmplx(0.,1.)*ha(i1,i2)
1047  IF(i1.LE.2) hc(i1,i2)=cmplx(0.,1.)*hc(i1,i2)
1048  ha(i2,i1)=-ha(i1,i2)
1049  hc(i2,i1)=-hc(i1,i2)
1050  530 CONTINUE
1051  540 CONTINUE
1052  ENDIF
1053 
1054 C...Calculate four-products.
1055  IF(isub.NE.0) THEN
1056  DO 560 i=1,2
1057  DO 550 j=1,4
1058  pk(i,j)=-pk(i,j)
1059  550 CONTINUE
1060  560 CONTINUE
1061  DO 580 i1=imin,imax-1
1062  DO 570 i2=i1+1,imax
1063  pkk(i1,i2)=2d0*(pk(i1,4)*pk(i2,4)-pk(i1,1)*pk(i2,1)-
1064  & pk(i1,2)*pk(i2,2)-pk(i1,3)*pk(i2,3))
1065  pkk(i2,i1)=pkk(i1,i2)
1066  570 CONTINUE
1067  580 CONTINUE
1068  ENDIF
1069  ENDIF
1070 
1071  kfagm=iabs(iref(ip,7))
1072  IF(mstp(47).LE.0.OR.ninh.NE.0) THEN
1073 C...Isotropic decay selected by user.
1074  wt=1d0
1075  wtmax=1d0
1076 
1077  ELSEIF(jtmax.EQ.3) THEN
1078 C...Isotropic decay when three mother particles.
1079  wt=1d0
1080  wtmax=1d0
1081 
1082  ELSEIF(it4.GE.1) THEN
1083 C... Isotropic decay t -> b + W etc for 4th generation q and l.
1084  wt=1d0
1085  wtmax=1d0
1086 
1087  ELSEIF(iref(ip,7).EQ.25.OR.iref(ip,7).EQ.35.OR.
1088  & iref(ip,7).EQ.36) THEN
1089 C...Angular weight for h0/A0 -> Z0 + Z0 or W+ + W- -> 4 quarks/leptons.
1090 C...CP-odd case added by Kari Ertresvag Myklevoll.
1091 C...Now also with mixed Higgs CP-states
1092  eta=parp(25)
1093  IF(ip.EQ.1) wtmax=sh**2
1094  IF(ip.GE.2) wtmax=p(iref(ip,8),5)**4
1095  kfa=iabs(k(iref(ip,1),2))
1096  kft=iabs(k(iref(ip,2),2))
1097 
1098  IF((kfa.EQ.kft).AND.(kfa.EQ.23.OR.kfa.EQ.24).AND.
1099  & mstp(25).GE.3) THEN
1100 C...For mixed CP states need epsilon product.
1101  p10=pk(3,4)
1102  p20=pk(4,4)
1103  p30=pk(5,4)
1104  p40=pk(6,4)
1105  p11=pk(3,1)
1106  p21=pk(4,1)
1107  p31=pk(5,1)
1108  p41=pk(6,1)
1109  p12=pk(3,2)
1110  p22=pk(4,2)
1111  p32=pk(5,2)
1112  p42=pk(6,2)
1113  p13=pk(3,3)
1114  p23=pk(4,3)
1115  p33=pk(5,3)
1116  p43=pk(6,3)
1117  epsi=p10*p21*p32*p43-p10*p21*p33*p42-p10*p22*p31*p43+p10*p22*
1118  & p33*p41+p10*p23*p31*p42-p10*p23*p32*p41-p11*p20*p32*p43+p11*
1119  & p20*p33*p42+p11*p22*p30*p43-p11*p22*p33*p40-p11*p23*p30*p42+
1120  & p11*p23*p32*p40+p12*p20*p31*p43-p12*p20*p33*p41-p12*p21*p30*
1121  & p43+p12*p21*p33*p40+p12*p23*p30*p41-p12*p23*p31*p40-p13*p20*
1122  & p31*p42+p13*p20*p32*p41+p13*p21*p30*p42-p13*p21*p32*p40-p13*
1123  & p22*p30*p41+p13*p22*p31*p40
1124 C...For mixed CP states need gauge boson masses.
1125  xma=sqrt(max(0d0,(pk(3,4)+pk(4,4))**2-(pk(3,1)+pk(4,1))**2-
1126  & (pk(3,2)+pk(4,2))**2-(pk(3,3)+pk(4,3))**2))
1127  xmb=sqrt(max(0d0,(pk(5,4)+pk(6,4))**2-(pk(5,1)+pk(6,1))**2-
1128  & (pk(5,2)+pk(6,2))**2-(pk(5,3)+pk(6,3))**2))
1129  xmv=pmas(kfa,1)
1130  ENDIF
1131 
1132 C...Z decay
1133  IF(kfa.EQ.23.AND.(kfa.EQ.kft)) THEN
1134  kflf1a=iabs(kfl1(1))
1135  ef1=kchg(kflf1a,1)/3d0
1136  af1=sign(1d0,ef1+0.1d0)
1137  vf1=af1-4d0*ef1*xwv
1138  kflf2a=iabs(kfl1(2))
1139  ef2=kchg(kflf2a,1)/3d0
1140  af2=sign(1d0,ef2+0.1d0)
1141  vf2=af2-4d0*ef2*xwv
1142  va12as=4d0*vf1*af1*vf2*af2/((vf1**2+af1**2)*(vf2**2+af2**2))
1143  IF((mstp(25).EQ.0.AND.iref(ip,7).NE.36).OR.mstp(25).EQ.1)
1144  & THEN
1145 C...CP-even decay
1146  wt=8d0*(1d0+va12as)*pkk(3,5)*pkk(4,6)+
1147  & 8d0*(1d0-va12as)*pkk(3,6)*pkk(4,5)
1148  ELSEIF(mstp(25).LE.2) THEN
1149 C...CP-odd decay
1150  wt=((pkk(3,5)+pkk(4,6))**2 +(pkk(3,6)+pkk(4,5))**2
1151  & -2*pkk(3,4)*pkk(5,6)
1152  & -2*(pkk(3,5)*pkk(4,6)-pkk(3,6)*pkk(4,5))**2/
1153  & (pkk(3,4)*pkk(5,6))
1154  & +va12as*(pkk(3,5)+pkk(3,6)-pkk(4,5)-pkk(4,6))*
1155  & (pkk(3,5)+pkk(4,5)-pkk(3,6)-pkk(4,6)))/(1+va12as)
1156  ELSE
1157 C...Mixed CP states.
1158  wt=32d0*(0.25d0*((1d0+va12as)*pkk(3,5)*pkk(4,6)
1159  & +(1d0-va12as)*pkk(3,6)*pkk(4,5))
1160  & -0.5d0*eta/xmv**2*epsi*((1d0+va12as)*(pkk(3,5)+pkk(4,6))
1161  & -(1d0-va12as)*(pkk(3,6)+pkk(4,5)))
1162  & +6.25d-2*eta**2/xmv**4*(-2d0*pkk(3,4)**2*pkk(5,6)**2
1163  & -2d0*(pkk(3,5)*pkk(4,6)-pkk(3,6)*pkk(4,5))**2
1164  & +pkk(3,4)*pkk(5,6)
1165  & *((pkk(3,5)+pkk(4,6))**2+(pkk(3,6)+pkk(4,5))**2)
1166  & +va12as*pkk(3,4)*pkk(5,6)
1167  & *(pkk(3,5)+pkk(3,6)-pkk(4,5)-pkk(4,6))
1168  & *(pkk(3,5)-pkk(3,6)+pkk(4,5)-pkk(4,6))))
1169  & /(1d0 +2d0*eta*xma*xmb/xmv**2
1170  & +2d0*(eta*xma*xmb/xmv**2)**2*(1d0+va12as))
1171  ENDIF
1172 
1173 C...W decay
1174  ELSEIF(kfa.EQ.24.AND.(kfa.EQ.kft)) THEN
1175  IF((mstp(25).EQ.0.AND.iref(ip,7).NE.36).OR.mstp(25).EQ.1)
1176  & THEN
1177 C...CP-even decay
1178  wt=16d0*pkk(3,5)*pkk(4,6)
1179  ELSEIF(mstp(25).LE.2) THEN
1180 C...CP-odd decay
1181  wt=0.5d0*((pkk(3,5)+pkk(4,6))**2 +(pkk(3,6)+pkk(4,5))**2
1182  & -2*pkk(3,4)*pkk(5,6)
1183  & -2*(pkk(3,5)*pkk(4,6)-pkk(3,6)*pkk(4,5))**2/
1184  & (pkk(3,4)*pkk(5,6))
1185  & +(pkk(3,5)+pkk(3,6)-pkk(4,5)-pkk(4,6))*
1186  & (pkk(3,5)+pkk(4,5)-pkk(3,6)-pkk(4,6)))
1187  ELSE
1188 C...Mixed CP states.
1189  wt=32d0*(0.25d0*2d0*pkk(3,5)*pkk(4,6)
1190  & -0.5d0*eta/xmv**2*epsi*2d0*(pkk(3,5)+pkk(4,6))
1191  & +6.25d-2*eta**2/xmv**4*(-2d0*pkk(3,4)**2*pkk(5,6)**2
1192  & -2d0*(pkk(3,5)*pkk(4,6)-pkk(3,6)*pkk(4,5))**2
1193  & +pkk(3,4)*pkk(5,6)
1194  & *((pkk(3,5)+pkk(4,6))**2+(pkk(3,6)+pkk(4,5))**2)
1195  & +pkk(3,4)*pkk(5,6)
1196  & *(pkk(3,5)+pkk(3,6)-pkk(4,5)-pkk(4,6))
1197  & *(pkk(3,5)-pkk(3,6)+pkk(4,5)-pkk(4,6))))
1198  & /(1d0 +2d0*eta*xma*xmb/xmv**2
1199  & +(2d0*eta*xma*xmb/xmv**2)**2)
1200  ENDIF
1201 
1202 C...No angular correlations in other Higgs decays.
1203  ELSE
1204  wt=wtmax
1205  ENDIF
1206 
1207  ELSEIF((kfagm.EQ.6.OR.kfagm.EQ.7.OR.kfagm.EQ.8.OR.
1208  & kfagm.EQ.17.OR.kfagm.EQ.18).AND.iabs(k(iref(ip,1),2)).EQ.24)
1209  & THEN
1210 C...Angular correlation in f -> f' + W -> f' + 2 quarks/leptons.
1211  i1=iref(ip,8)
1212  IF(mod(kfagm,2).EQ.0) THEN
1213  i2=n+1
1214  i3=n+2
1215  ELSE
1216  i2=n+2
1217  i3=n+1
1218  ENDIF
1219  i4=iref(ip,2)
1220  wt=(p(i1,4)*p(i2,4)-p(i1,1)*p(i2,1)-p(i1,2)*p(i2,2)-
1221  & p(i1,3)*p(i2,3))*(p(i3,4)*p(i4,4)-p(i3,1)*p(i4,1)-
1222  & p(i3,2)*p(i4,2)-p(i3,3)*p(i4,3))
1223  wtmax=(p(i1,5)**4-p(iref(ip,1),5)**4)/8d0
1224 
1225  ELSEIF(isub.EQ.1) THEN
1226 C...Angular weight for gamma*/Z0 -> 2 quarks/leptons.
1227  ei=kchg(iabs(mint(15)),1)/3d0
1228  ai=sign(1d0,ei+0.1d0)
1229  vi=ai-4d0*ei*xwv
1230  ef=kchg(iabs(kfl1(1)),1)/3d0
1231  af=sign(1d0,ef+0.1d0)
1232 
1233  vf=af-4d0*ef*xwv
1234  rmf=min(1d0,4d0*pmas(iabs(kfl1(1)),1)**2/sh)
1235  wt1=ei**2*vint(111)*ef**2+ei*vi*vint(112)*ef*vf+
1236  & (vi**2+ai**2)*vint(114)*(vf**2+(1d0-rmf)*af**2)
1237  wt2=rmf*(ei**2*vint(111)*ef**2+ei*vi*vint(112)*ef*vf+
1238  & (vi**2+ai**2)*vint(114)*vf**2)
1239  wt3=sqrt(1d0-rmf)*(ei*ai*vint(112)*ef*af+
1240  & 4d0*vi*ai*vint(114)*vf*af)
1241  wt=wt1*(1d0+cthe(1)**2)+wt2*(1d0-cthe(1)**2)+
1242  & 2d0*wt3*cthe(1)*isign(1,mint(15)*kfl1(1))
1243  wtmax=2d0*(wt1+abs(wt3))
1244 
1245  ELSEIF(isub.EQ.2) THEN
1246 C...Angular weight for W+/- -> 2 quarks/leptons.
1247  rm3=pmas(iabs(kfl1(1)),1)**2/sh
1248  rm4=pmas(iabs(kfl2(1)),1)**2/sh
1249  be34=sqrt(max(0d0,(1d0-rm3-rm4)**2-4d0*rm3*rm4))
1250  wt=(1d0+be34*cthe(1)*isign(1,mint(15)*kfl1(1)))**2-(rm3-rm4)**2
1251  wtmax=4d0
1252 
1253  ELSEIF(isub.EQ.15.OR.isub.EQ.19) THEN
1254 C...Angular weight for f + fbar -> gluon/gamma + (gamma*/Z0) ->
1255 C...-> gluon/gamma + 2 quarks/leptons.
1256  clilf=coup(1,1)**2*hgz(jtz,1)*coup(3,1)**2/16d0+
1257  & coup(1,1)*coup(1,3)*hgz(jtz,2)*coup(3,1)*coup(3,3)/4d0+
1258  & coup(1,3)**2*hgz(jtz,3)*coup(3,3)**2
1259  clirf=coup(1,1)**2*hgz(jtz,1)*coup(3,1)**2/16d0+
1260  & coup(1,1)*coup(1,3)*hgz(jtz,2)*coup(3,1)*coup(3,4)/4d0+
1261  & coup(1,3)**2*hgz(jtz,3)*coup(3,4)**2
1262  crilf=coup(1,1)**2*hgz(jtz,1)*coup(3,1)**2/16d0+
1263  & coup(1,1)*coup(1,4)*hgz(jtz,2)*coup(3,1)*coup(3,3)/4d0+
1264  & coup(1,4)**2*hgz(jtz,3)*coup(3,3)**2
1265  crirf=coup(1,1)**2*hgz(jtz,1)*coup(3,1)**2/16d0+
1266  & coup(1,1)*coup(1,4)*hgz(jtz,2)*coup(3,1)*coup(3,4)/4d0+
1267  & coup(1,4)**2*hgz(jtz,3)*coup(3,4)**2
1268  wt=(clilf+crirf)*(pkk(1,3)**2+pkk(2,4)**2)+
1269  & (clirf+crilf)*(pkk(1,4)**2+pkk(2,3)**2)
1270  wtmax=(clilf+clirf+crilf+crirf)*
1271  & ((pkk(1,3)+pkk(1,4))**2+(pkk(2,3)+pkk(2,4))**2)
1272 
1273  ELSEIF(isub.EQ.16.OR.isub.EQ.20) THEN
1274 C...Angular weight for f + fbar' -> gluon/gamma + W+/- ->
1275 C...-> gluon/gamma + 2 quarks/leptons.
1276  wt=pkk(1,3)**2+pkk(2,4)**2
1277  wtmax=(pkk(1,3)+pkk(1,4))**2+(pkk(2,3)+pkk(2,4))**2
1278 
1279  ELSEIF(isub.EQ.22) THEN
1280 C...Angular weight for f + fbar -> Z0 + Z0 -> 4 quarks/leptons.
1281  s34=p(iref(ip,iord),5)**2
1282  s56=p(iref(ip,3-iord),5)**2
1283  ti=pkk(1,3)+pkk(1,4)+s34
1284  ui=pkk(1,5)+pkk(1,6)+s56
1285  tir=REAL(ti)
1286  uir=REAL(ui)
1287  fgk135=abs(fgk(1,2,3,4,5,6)/tir+fgk(1,2,5,6,3,4)/uir)**2
1288  fgk145=abs(fgk(1,2,4,3,5,6)/tir+fgk(1,2,5,6,4,3)/uir)**2
1289  fgk136=abs(fgk(1,2,3,4,6,5)/tir+fgk(1,2,6,5,3,4)/uir)**2
1290  fgk146=abs(fgk(1,2,4,3,6,5)/tir+fgk(1,2,6,5,4,3)/uir)**2
1291  fgk253=abs(fgk(2,1,5,6,3,4)/tir+fgk(2,1,3,4,5,6)/uir)**2
1292  fgk263=abs(fgk(2,1,6,5,3,4)/tir+fgk(2,1,3,4,6,5)/uir)**2
1293  fgk254=abs(fgk(2,1,5,6,4,3)/tir+fgk(2,1,4,3,5,6)/uir)**2
1294  fgk264=abs(fgk(2,1,6,5,4,3)/tir+fgk(2,1,4,3,6,5)/uir)**2
1295 
1296  wt=
1297  & corl(1,1,1)*corl(2,1,1)*fgk135+corl(1,1,2)*corl(2,1,1)*fgk145+
1298  & corl(1,1,1)*corl(2,1,2)*fgk136+corl(1,1,2)*corl(2,1,2)*fgk146+
1299  & corl(1,2,1)*corl(2,2,1)*fgk253+corl(1,2,2)*corl(2,2,1)*fgk263+
1300  & corl(1,2,1)*corl(2,2,2)*fgk254+corl(1,2,2)*corl(2,2,2)*fgk264
1301  wtmax=16d0*((corl(1,1,1)+corl(1,1,2))*(corl(2,1,1)+corl(2,1,2))+
1302  & (corl(1,2,1)+corl(1,2,2))*(corl(2,2,1)+corl(2,2,2)))*s34*s56*
1303  & ((ti**2+ui**2+2d0*sh*(s34+s56))/(ti*ui)-s34*s56*(1d0/ti**2+
1304  & 1d0/ui**2))
1305 
1306  ELSEIF(isub.EQ.23) THEN
1307 C...Angular weight for f + fbar' -> Z0 + W+/- -> 4 quarks/leptons.
1308  d34=p(iref(ip,iord),5)**2
1309  d56=p(iref(ip,3-iord),5)**2
1310  dt=pkk(1,3)+pkk(1,4)+d34
1311  du=pkk(1,5)+pkk(1,6)+d56
1312  facbw=1d0/((sh-sqmw)**2+gmmw**2)
1313  cawz=coup(2,3)/dt-2d0*xw1*coup(1,2)*(sh-sqmw)*facbw
1314  cbwz=coup(1,3)/du+2d0*xw1*coup(1,2)*(sh-sqmw)*facbw
1315  fgk135=abs(REAL(cawz)*fgk(1,2,3,4,5,6)+
1316 
1317  & REAL(cbwz)*fgk(1,2,5,6,3,4))
1318  fgk136=abs(REAL(cawz)*fgk(1,2,3,4,6,5)+
1319  & REAL(cbwz)*fgk(1,2,6,5,3,4))
1320  wt=(coup(5,3)*fgk135)**2+(coup(5,4)*fgk136)**2
1321  wtmax=4d0*d34*d56*(coup(5,3)**2+coup(5,4)**2)*(cawz**2*
1322  & digk(dt,du)+cbwz**2*digk(du,dt)+cawz*cbwz*djgk(dt,du))
1323 
1324  ELSEIF(isub.EQ.24.OR.isub.EQ.171.OR.isub.EQ.176) THEN
1325 C...Angular weight for f + fbar -> Z0 + h0 -> 2 quarks/leptons + h0
1326 C...(or H0, or A0).
1327  wt=((coup(1,3)*coup(3,3))**2+(coup(1,4)*coup(3,4))**2)*
1328  & pkk(1,3)*pkk(2,4)+((coup(1,3)*coup(3,4))**2+(coup(1,4)*
1329  & coup(3,3))**2)*pkk(1,4)*pkk(2,3)
1330  wtmax=(coup(1,3)**2+coup(1,4)**2)*(coup(3,3)**2+coup(3,4)**2)*
1331  & (pkk(1,3)+pkk(1,4))*(pkk(2,3)+pkk(2,4))
1332 
1333  ELSEIF(isub.EQ.25) THEN
1334 C...Angular weight for f + fbar -> W+ + W- -> 4 quarks/leptons.
1335  polr=(1d0+parj(132))*(1d0-parj(131))
1336  poll=(1d0-parj(132))*(1d0+parj(131))
1337  d34=p(iref(ip,iord),5)**2
1338  d56=p(iref(ip,3-iord),5)**2
1339  dt=pkk(1,3)+pkk(1,4)+d34
1340  du=pkk(1,5)+pkk(1,6)+d56
1341  facbw=1d0/((sh-sqmz)**2+sqmz*pmas(23,2)**2)
1342  cdww=(coup(1,3)*sqmz*(sh-sqmz)*facbw+coup(1,2))/sh
1343  caww=cdww+0.5d0*(coup(1,2)+1d0)/dt
1344  cbww=cdww+0.5d0*(coup(1,2)-1d0)/du
1345  ccww=coup(1,4)*sqmz*(sh-sqmz)*facbw/sh
1346  fgk135=abs(REAL(caww)*fgk(1,2,3,4,5,6)-
1347  & REAL(cbww)*fgk(1,2,5,6,3,4))
1348  fgk253=abs(fgk(2,1,5,6,3,4)-fgk(2,1,3,4,5,6))
1349  IF(mstp(50).LE.0) THEN
1350  wt=fgk135**2+(ccww*fgk253)**2
1351  wtmax=4d0*d34*d56*(caww**2*digk(dt,du)+cbww**2*digk(du,dt)-
1352  & caww*cbww*djgk(dt,du)+ccww**2*(digk(dt,du)+digk(du,dt)-
1353  & djgk(dt,du)))
1354  ELSE
1355  wt=poll*fgk135**2+polr*(ccww*fgk253)**2
1356  wtmax=4d0*d34*d56*(poll*(caww**2*digk(dt,du)+
1357  & cbww**2*digk(du,dt)-caww*cbww*djgk(dt,du))+
1358  & polr*ccww**2*(digk(dt,du)+digk(du,dt)-djgk(dt,du)))
1359  ENDIF
1360 
1361  ELSEIF(isub.EQ.26.OR.isub.EQ.172.OR.isub.EQ.177) THEN
1362 C...Angular weight for f + fbar' -> W+/- + h0 -> 2 quarks/leptons + h0
1363 C...(or H0, or A0).
1364  wt=pkk(1,3)*pkk(2,4)
1365  wtmax=(pkk(1,3)+pkk(1,4))*(pkk(2,3)+pkk(2,4))
1366 
1367  ELSEIF(isub.EQ.30.OR.isub.EQ.35) THEN
1368 C...Angular weight for f + g/gamma -> f + (gamma*/Z0)
1369 C...-> f + 2 quarks/leptons.
1370  clilf=coup(1,1)**2*hgz(jtz,1)*coup(3,1)**2/16d0+
1371  & coup(1,1)*coup(1,3)*hgz(jtz,2)*coup(3,1)*coup(3,3)/4d0+
1372  & coup(1,3)**2*hgz(jtz,3)*coup(3,3)**2
1373  clirf=coup(1,1)**2*hgz(jtz,1)*coup(3,1)**2/16d0+
1374  & coup(1,1)*coup(1,3)*hgz(jtz,2)*coup(3,1)*coup(3,4)/4d0+
1375  & coup(1,3)**2*hgz(jtz,3)*coup(3,4)**2
1376  crilf=coup(1,1)**2*hgz(jtz,1)*coup(3,1)**2/16d0+
1377  & coup(1,1)*coup(1,4)*hgz(jtz,2)*coup(3,1)*coup(3,3)/4d0+
1378  & coup(1,4)**2*hgz(jtz,3)*coup(3,3)**2
1379  crirf=coup(1,1)**2*hgz(jtz,1)*coup(3,1)**2/16d0+
1380  & coup(1,1)*coup(1,4)*hgz(jtz,2)*coup(3,1)*coup(3,4)/4d0+
1381  & coup(1,4)**2*hgz(jtz,3)*coup(3,4)**2
1382  IF(k(ilin(1),2).GT.0) wt=(clilf+crirf)*(pkk(1,4)**2+
1383  & pkk(3,5)**2)+(clirf+crilf)*(pkk(1,3)**2+pkk(4,5)**2)
1384  IF(k(ilin(1),2).LT.0) wt=(clilf+crirf)*(pkk(1,3)**2+
1385  & pkk(4,5)**2)+(clirf+crilf)*(pkk(1,4)**2+pkk(3,5)**2)
1386  wtmax=(clilf+clirf+crilf+crirf)*
1387  & ((pkk(1,3)+pkk(1,4))**2+(pkk(3,5)+pkk(4,5))**2)
1388 
1389  ELSEIF(isub.EQ.31.OR.isub.EQ.36) THEN
1390 C...Angular weight for f + g/gamma -> f' + W+/- -> f' + 2 fermions.
1391  IF(k(ilin(1),2).GT.0) wt=pkk(1,4)**2+pkk(3,5)**2
1392  IF(k(ilin(1),2).LT.0) wt=pkk(1,3)**2+pkk(4,5)**2
1393  wtmax=(pkk(1,3)+pkk(1,4))**2+(pkk(3,5)+pkk(4,5))**2
1394 
1395  ELSEIF(isub.EQ.71.OR.isub.EQ.72.OR.isub.EQ.73.OR.isub.EQ.76.OR.
1396  & isub.EQ.77) THEN
1397 C...Angular weight for V_L1 + V_L2 -> V_L3 + V_L4 (V = Z/W).
1398  wt=16d0*pkk(3,5)*pkk(4,6)
1399  wtmax=sh**2
1400 
1401  ELSEIF(isub.EQ.110) THEN
1402 C...Angular weight for f + fbar -> gamma + h0 -> gamma + X is isotropic.
1403  wt=1d0
1404  wtmax=1d0
1405 
1406  ELSEIF(isub.EQ.141) THEN
1407  IF(ip.EQ.1.AND.iabs(kfl1(1)).LT.20) THEN
1408 C...Angular weight for f + fbar -> gamma*/Z0/Z'0 -> 2 quarks/leptons.
1409 C...Couplings of incoming flavour.
1410  kfai=iabs(mint(15))
1411  ei=kchg(kfai,1)/3d0
1412  ai=sign(1d0,ei+0.1d0)
1413  vi=ai-4d0*ei*xwv
1414  kfaic=1
1415  IF(kfai.LE.10.AND.mod(kfai,2).EQ.0) kfaic=2
1416  IF(kfai.GT.10.AND.mod(kfai,2).NE.0) kfaic=3
1417  IF(kfai.GT.10.AND.mod(kfai,2).EQ.0) kfaic=4
1418  IF(kfai.LE.2.OR.kfai.EQ.11.OR.kfai.EQ.12) THEN
1419  vpi=paru(119+2*kfaic)
1420  api=paru(120+2*kfaic)
1421  ELSEIF(kfai.LE.4.OR.kfai.EQ.13.OR.kfai.EQ.14) THEN
1422  vpi=parj(178+2*kfaic)
1423  api=parj(179+2*kfaic)
1424  ELSE
1425  vpi=parj(186+2*kfaic)
1426  api=parj(187+2*kfaic)
1427  ENDIF
1428 C...Couplings of final flavour.
1429  kfaf=iabs(kfl1(1))
1430  ef=kchg(kfaf,1)/3d0
1431  af=sign(1d0,ef+0.1d0)
1432  vf=af-4d0*ef*xwv
1433  kfafc=1
1434  IF(kfaf.LE.10.AND.mod(kfaf,2).EQ.0) kfafc=2
1435  IF(kfaf.GT.10.AND.mod(kfaf,2).NE.0) kfafc=3
1436  IF(kfaf.GT.10.AND.mod(kfaf,2).EQ.0) kfafc=4
1437  IF(kfaf.LE.2.OR.kfaf.EQ.11.OR.kfaf.EQ.12) THEN
1438  vpf=paru(119+2*kfafc)
1439  apf=paru(120+2*kfafc)
1440  ELSEIF(kfaf.LE.4.OR.kfaf.EQ.13.OR.kfaf.EQ.14) THEN
1441  vpf=parj(178+2*kfafc)
1442  apf=parj(179+2*kfafc)
1443  ELSE
1444  vpf=parj(186+2*kfafc)
1445  apf=parj(187+2*kfafc)
1446  ENDIF
1447 C...Asymmetry and weight.
1448  asym=2d0*(ei*ai*vint(112)*ef*af+ei*api*vint(113)*ef*apf+
1449  & 4d0*vi*ai*vint(114)*vf*af+(vi*api+vpi*ai)*vint(115)*
1450  & (vf*apf+vpf*af)+4d0*vpi*api*vint(116)*vpf*apf)/
1451  & (ei**2*vint(111)*ef**2+ei*vi*vint(112)*ef*vf+
1452  & ei*vpi*vint(113)*ef*vpf+(vi**2+ai**2)*vint(114)*
1453  & (vf**2+af**2)+(vi*vpi+ai*api)*vint(115)*(vf*vpf+af*apf)+
1454  & (vpi**2+api**2)*vint(116)*(vpf**2+apf**2))
1455  wt=1d0+asym*cthe(1)*isign(1,mint(15)*kfl1(1))+cthe(1)**2
1456  wtmax=2d0+abs(asym)
1457  ELSEIF(ip.EQ.1.AND.iabs(kfl1(1)).EQ.24) THEN
1458 C...Angular weight for f + fbar -> Z' -> W+ + W-.
1459  rm1=p(nsd(1)+1,5)**2/sh
1460  rm2=p(nsd(1)+2,5)**2/sh
1461  ccos2=-(1d0/16d0)*((1d0-rm1-rm2)**2-4d0*rm1*rm2)*
1462  & (1d0-2d0*rm1-2d0*rm2+rm1**2+rm2**2+10d0*rm1*rm2)
1463  cflat=-ccos2+0.5d0*(rm1+rm2)*(1d0-2d0*rm1-2d0*rm2+
1464  & (rm2-rm1)**2)
1465  wt=cflat+ccos2*cthe(1)**2
1466  wtmax=cflat+max(0d0,ccos2)
1467  ELSEIF(ip.EQ.1.AND.(kfl1(1).EQ.25.OR.kfl1(1).EQ.35.OR.
1468  & iabs(kfl1(1)).EQ.37)) THEN
1469 C...Angular weight for f + fbar -> Z' -> h0 + A0, H0 + A0, H+ + H-.
1470  wt=1d0-cthe(1)**2
1471  wtmax=1d0
1472  ELSEIF(ip.EQ.1.AND.kfl2(1).EQ.25) THEN
1473 C...Angular weight for f + fbar -> Z' -> Z0 + h0.
1474  rm1=p(nsd(1)+1,5)**2/sh
1475  rm2=p(nsd(1)+2,5)**2/sh
1476  flam2=max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2)
1477  wt=1d0+flam2*(1d0-cthe(1)**2)/(8d0*rm1)
1478  wtmax=1d0+flam2/(8d0*rm1)
1479  ELSEIF(mzpwp.EQ.0) THEN
1480 C...Angular weight for f + fbar -> Z' -> W+ + W- -> 4 quarks/leptons
1481 C...(W:s like if intermediate Z).
1482  d34=p(iref(ip,iord),5)**2
1483  d56=p(iref(ip,3-iord),5)**2
1484  dt=pkk(1,3)+pkk(1,4)+d34
1485  du=pkk(1,5)+pkk(1,6)+d56
1486  fgk135=abs(fgk(1,2,3,4,5,6)-fgk(1,2,5,6,3,4))
1487  fgk253=abs(fgk(2,1,5,6,3,4)-fgk(2,1,3,4,5,6))
1488  wt=(coup(1,3)*fgk135)**2+(coup(1,4)*fgk253)**2
1489  wtmax=4d0*d34*d56*(coup(1,3)**2+coup(1,4)**2)*
1490  & (digk(dt,du)+digk(du,dt)-djgk(dt,du))
1491  ELSEIF(mzpwp.EQ.1) THEN
1492 C...Angular weight for f + fbar -> Z' -> W+ + W- -> 4 quarks/leptons
1493 C...(W:s approximately longitudinal, like if intermediate H).
1494  wt=16d0*pkk(3,5)*pkk(4,6)
1495  wtmax=sh**2
1496  ELSE
1497 C...Angular weight for f + fbar -> Z' -> H+ + H-, Z0 + h0, h0 + A0,
1498 C...H0 + A0 -> 4 quarks/leptons, t + tbar -> b + W+ + bbar + W- .
1499  wt=1d0
1500  wtmax=1d0
1501  ENDIF
1502 
1503  ELSEIF(isub.EQ.142) THEN
1504  IF(ip.EQ.1.AND.iabs(kfl1(1)).LT.20) THEN
1505 C...Angular weight for f + fbar' -> W'+/- -> 2 quarks/leptons.
1506  kfai=iabs(mint(15))
1507  kfaic=1
1508  IF(kfai.GT.10) kfaic=2
1509  vi=paru(129+2*kfaic)
1510  ai=paru(130+2*kfaic)
1511  kfaf=iabs(kfl1(1))
1512  kfafc=1
1513  IF(kfaf.GT.10) kfafc=2
1514  vf=paru(129+2*kfafc)
1515  af=paru(130+2*kfafc)
1516  asym=8d0*vi*ai*vf*af/((vi**2+ai**2)*(vf**2+af**2))
1517  wt=1d0+asym*cthe(1)*isign(1,mint(15)*kfl1(1))+cthe(1)**2
1518  wtmax=2d0+abs(asym)
1519  ELSEIF(ip.EQ.1.AND.iabs(kfl2(1)).EQ.23) THEN
1520 C...Angular weight for f + fbar' -> W'+/- -> W+/- + Z0.
1521  rm1=p(nsd(1)+1,5)**2/sh
1522  rm2=p(nsd(1)+2,5)**2/sh
1523  ccos2=-(1d0/16d0)*((1d0-rm1-rm2)**2-4d0*rm1*rm2)*
1524  & (1d0-2d0*rm1-2d0*rm2+rm1**2+rm2**2+10d0*rm1*rm2)
1525  cflat=-ccos2+0.5d0*(rm1+rm2)*(1d0-2d0*rm1-2d0*rm2+
1526  & (rm2-rm1)**2)
1527  wt=cflat+ccos2*cthe(1)**2
1528  wtmax=cflat+max(0d0,ccos2)
1529  ELSEIF(ip.EQ.1.AND.kfl2(1).EQ.25) THEN
1530 C...Angular weight for f + fbar -> W'+/- -> W+/- + h0.
1531  rm1=p(nsd(1)+1,5)**2/sh
1532  rm2=p(nsd(1)+2,5)**2/sh
1533  flam2=max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2)
1534  wt=1d0+flam2*(1d0-cthe(1)**2)/(8d0*rm1)
1535  wtmax=1d0+flam2/(8d0*rm1)
1536  ELSEIF(mzpwp.EQ.0) THEN
1537 C...Angular weight for f + fbar' -> W' -> W + Z0 -> 4 quarks/leptons
1538 C...(W/Z like if intermediate W).
1539  d34=p(iref(ip,iord),5)**2
1540  d56=p(iref(ip,3-iord),5)**2
1541  dt=pkk(1,3)+pkk(1,4)+d34
1542  du=pkk(1,5)+pkk(1,6)+d56
1543  fgk135=abs(fgk(1,2,3,4,5,6)-fgk(1,2,5,6,3,4))
1544  fgk136=abs(fgk(1,2,3,4,6,5)-fgk(1,2,6,5,3,4))
1545  wt=(coup(5,3)*fgk135)**2+(coup(5,4)*fgk136)**2
1546  wtmax=4d0*d34*d56*(coup(5,3)**2+coup(5,4)**2)*
1547  & (digk(dt,du)+digk(du,dt)-djgk(dt,du))
1548  ELSEIF(mzpwp.EQ.1) THEN
1549 C...Angular weight for f + fbar' -> W' -> W + Z0 -> 4 quarks/leptons
1550 C...(W/Z approximately longitudinal, like if intermediate H).
1551  wt=16d0*pkk(3,5)*pkk(4,6)
1552  wtmax=sh**2
1553  ELSE
1554 C...Angular weight for f + fbar -> W' -> W + h0 -> whatever,
1555 C...t + bbar -> t + W + bbar.
1556  wt=1d0
1557  wtmax=1d0
1558  ENDIF
1559 
1560  ELSEIF(isub.EQ.145.OR.isub.EQ.162.OR.isub.EQ.163.OR.isub.EQ.164)
1561  & THEN
1562 C...Isotropic decay of leptoquarks (assumed spin 0).
1563  wt=1d0
1564  wtmax=1d0
1565 
1566  ELSEIF(isub.GE.146.AND.isub.LE.148) THEN
1567 C...Decays of (spin 1/2) q*/e* -> q/e + (g,gamma) or (Z0,W+-).
1568  side=1d0
1569  IF(mint(16).EQ.21.OR.mint(16).EQ.22) side=-1d0
1570  IF(ip.EQ.1.AND.(kfl1(1).EQ.21.OR.kfl1(1).EQ.22)) THEN
1571  wt=1d0+side*cthe(1)
1572  wtmax=2d0
1573  ELSEIF(ip.EQ.1) THEN
1574 
1575  rm1=p(nsd(1)+1,5)**2/sh
1576  wt=1d0+side*cthe(1)*(1d0-0.5d0*rm1)/(1d0+0.5d0*rm1)
1577  wtmax=1d0+(1d0-0.5d0*rm1)/(1d0+0.5d0*rm1)
1578  ELSE
1579 C...W/Z decay assumed isotropic, since not known.
1580  wt=1d0
1581  wtmax=1d0
1582  ENDIF
1583 
1584  ELSEIF(isub.EQ.149) THEN
1585 C...Isotropic decay of techni-eta.
1586  wt=1d0
1587  wtmax=1d0
1588 
1589  ELSEIF(isub.EQ.191) THEN
1590  IF(ip.EQ.1.AND.iabs(kfl1(1)).GT.21) THEN
1591 C...Angular weight for f + fbar -> rho_tc0 -> W+ W-,
1592 C...W+ pi_tc-, pi_tc+ W- or pi_tc+ pi_tc-.
1593  wt=1d0-cthe(1)**2
1594  wtmax=1d0
1595  ELSEIF(ip.EQ.1) THEN
1596 C...Angular weight for f + fbar -> rho_tc0 -> f fbar.
1597  cthesg=cthe(1)*isign(1,mint(15))
1598  xwrht=(1d0-2d0*xw)/(4d0*xw*(1d0-xw))
1599  bwzr=xwrht*sh*(sh-sqmz)/((sh-sqmz)**2+gmmz**2)
1600  bwzi=xwrht*sh*gmmz/((sh-sqmz)**2+gmmz**2)
1601  kfai=iabs(mint(15))
1602  ei=kchg(kfai,1)/3d0
1603  ai=sign(1d0,ei+0.1d0)
1604  vi=ai-4d0*ei*xwv
1605  vali=0.5d0*(vi+ai)
1606  vari=0.5d0*(vi-ai)
1607  alefti=(ei+vali*bwzr)**2+(vali*bwzi)**2
1608  arighi=(ei+vari*bwzr)**2+(vari*bwzi)**2
1609  kfaf=iabs(kfl1(1))
1610  ef=kchg(kfaf,1)/3d0
1611  af=sign(1d0,ef+0.1d0)
1612  vf=af-4d0*ef*xwv
1613  valf=0.5d0*(vf+af)
1614  varf=0.5d0*(vf-af)
1615  aleftf=(ef+valf*bwzr)**2+(valf*bwzi)**2
1616  arighf=(ef+varf*bwzr)**2+(varf*bwzi)**2
1617  asame=alefti*aleftf+arighi*arighf
1618  aflip=alefti*arighf+arighi*aleftf
1619  wt=asame*(1d0+cthesg)**2+aflip*(1d0-cthesg)**2
1620  wtmax=4d0*max(asame,aflip)
1621  ELSE
1622 C...Isotropic decay of W/pi_tc produced in rho_tc decay.
1623  wt=1d0
1624  wtmax=1d0
1625  ENDIF
1626 
1627  ELSEIF(isub.EQ.192) THEN
1628  IF(ip.EQ.1.AND.iabs(kfl1(1)).GT.21) THEN
1629 C...Angular weight for f + fbar' -> rho_tc+ -> W+ Z0,
1630 C...W+ pi_tc0, pi_tc+ Z0 or pi_tc+ pi_tc0.
1631  wt=1d0-cthe(1)**2
1632  wtmax=1d0
1633  ELSEIF(ip.EQ.1) THEN
1634 C...Angular weight for f + fbar' -> rho_tc+ -> f fbar'.
1635  cthesg=cthe(1)*isign(1,mint(15))
1636  wt=(1d0+cthesg)**2
1637  wtmax=4d0
1638  ELSE
1639 C...Isotropic decay of W/Z/pi_tc produced in rho_tc+ decay.
1640  wt=1d0
1641  wtmax=1d0
1642  ENDIF
1643 
1644  ELSEIF(isub.EQ.193) THEN
1645  IF(ip.EQ.1.AND.iabs(kfl1(1)).GT.21) THEN
1646 C...Angular weight for f + fbar -> omega_tc0 ->
1647 C...gamma pi_tc0 or Z0 pi_tc0.
1648  wt=1d0+cthe(1)**2
1649  wtmax=2d0
1650  ELSEIF(ip.EQ.1) THEN
1651 C...Angular weight for f + fbar -> omega_tc0 -> f fbar.
1652  cthesg=cthe(1)*isign(1,mint(15))
1653  bwzr=(0.5d0/(1d0-xw))*sh*(sh-sqmz)/((sh-sqmz)**2+gmmz**2)
1654  bwzi=(0.5d0/(1d0-xw))*sh*gmmz/((sh-sqmz)**2+gmmz**2)
1655  kfai=iabs(mint(15))
1656  ei=kchg(kfai,1)/3d0
1657  ai=sign(1d0,ei+0.1d0)
1658  vi=ai-4d0*ei*xwv
1659  vali=0.5d0*(vi+ai)
1660  vari=0.5d0*(vi-ai)
1661  blefti=(ei-vali*bwzr)**2+(vali*bwzi)**2
1662  brighi=(ei-vari*bwzr)**2+(vari*bwzi)**2
1663  kfaf=iabs(kfl1(1))
1664  ef=kchg(kfaf,1)/3d0
1665  af=sign(1d0,ef+0.1d0)
1666  vf=af-4d0*ef*xwv
1667  valf=0.5d0*(vf+af)
1668  varf=0.5d0*(vf-af)
1669  bleftf=(ef-valf*bwzr)**2+(valf*bwzi)**2
1670  brighf=(ef-varf*bwzr)**2+(varf*bwzi)**2
1671  bsame=blefti*bleftf+brighi*brighf
1672  bflip=blefti*brighf+brighi*bleftf
1673  wt=bsame*(1d0+cthesg)**2+bflip*(1d0-cthesg)**2
1674  wtmax=4d0*max(bsame,bflip)
1675  ELSE
1676 C...Isotropic decay of Z/pi_tc produced in omega_tc decay.
1677  wt=1d0
1678  wtmax=1d0
1679  ENDIF
1680 
1681  ELSEIF(isub.EQ.353) THEN
1682 C...Angular weight for Z_R0 -> 2 quarks/leptons.
1683  ei=kchg(iabs(mint(15)),1)/3d0
1684  ai=sign(1d0,ei+0.1d0)
1685  vi=ai-4d0*ei*xwv
1686  ef=kchg(pycomp(kfl1(1)),1)/3d0
1687  af=sign(1d0,ef+0.1d0)
1688  vf=af-4d0*ef*xwv
1689  rmf=min(1d0,4d0*pmas(pycomp(kfl1(1)),1)**2/sh)
1690  wt1=(vi**2+ai**2)*(vf**2+(1d0-rmf)*af**2)
1691  wt2=rmf*(vi**2+ai**2)*vf**2
1692  wt3=sqrt(1d0-rmf)*4d0*vi*ai*vf*af
1693  wt=wt1*(1d0+cthe(1)**2)+wt2*(1d0-cthe(1)**2)+
1694  & 2d0*wt3*cthe(1)*isign(1,mint(15)*kfl1(1))
1695  wtmax=2d0*(wt1+abs(wt3))
1696 
1697  ELSEIF(isub.EQ.354) THEN
1698 C...Angular weight for W_R+/- -> 2 quarks/leptons.
1699  rm3=pmas(pycomp(kfl1(1)),1)**2/sh
1700  rm4=pmas(pycomp(kfl2(1)),1)**2/sh
1701  be34=sqrt(max(0d0,(1d0-rm3-rm4)**2-4d0*rm3*rm4))
1702  wt=(1d0+be34*cthe(1)*isign(1,mint(15)*kfl1(1)))**2-(rm3-rm4)**2
1703  wtmax=4d0
1704 
1705  ELSEIF(isub.EQ.391) THEN
1706 C...Angular weight for f + fbar -> G* -> f + fbar
1707  IF(ip.EQ.1.AND.iabs(kfl1(1)).LE.18) THEN
1708  wt=1d0-3d0*cthe(1)**2+4d0*cthe(1)**4
1709  wtmax=2d0
1710 C...Angular weight for f + fbar -> G* -> gamma + gamma or g + g
1711 C...implemented by M.-C. Lemaire
1712  ELSEIF(ip.EQ.1.AND.(iabs(kfl1(1)).EQ.21.OR.
1713  & iabs(kfl1(1)).EQ.22)) THEN
1714  wt=1d0-cthe(1)**4
1715  wtmax=1d0
1716 C...Other G* decays not yet implemented angular distributions.
1717  ELSE
1718  wt=1d0
1719  wtmax=1d0
1720  ENDIF
1721 
1722  ELSEIF(isub.EQ.392) THEN
1723 C...Angular weight for g + g -> G* -> f + fbar
1724  IF(ip.EQ.1.AND.iabs(kfl1(1)).LE.18) THEN
1725  wt=1d0-cthe(1)**4
1726  wtmax=1d0
1727 C...Angular weight for g + g -> G* -> gamma +gamma or g + g
1728 C...implemented by M.-C. Lemaire
1729  ELSEIF(ip.EQ.1.AND.(iabs(kfl1(1)).EQ.21.OR.
1730  & iabs(kfl1(1)).EQ.22)) THEN
1731  wt=1d0+6d0*cthe(1)**2+cthe(1)**4
1732  wtmax=8d0
1733 C...Other G* decays not yet implemented angular distributions.
1734  ELSE
1735  wt=1d0
1736  wtmax=1d0
1737  ENDIF
1738 
1739 C...Obtain correct angular distribution by rejection techniques.
1740  ELSE
1741  wt=1d0
1742  wtmax=1d0
1743  ENDIF
1744  IF(wt.LT.pyr(0)*wtmax) goto 430
1745 
1746 C...Construct massive four-vectors using angles chosen.
1747  590 DO 690 jt=1,jtmax
1748  IF(kdcy(jt).EQ.0) goto 690
1749  id=iref(ip,jt)
1750  DO 600 j=1,5
1751  dpmo(j)=p(id,j)
1752  600 CONTINUE
1753  dpmo(4)=sqrt(dpmo(1)**2+dpmo(2)**2+dpmo(3)**2+dpmo(5)**2)
1754 CMRENNA++
1755  IF(kfl3(jt).EQ.0) THEN
1756  CALL pyrobo(nsd(jt)+1,nsd(jt)+2,acos(cthe(jt)),phi(jt),
1757  & dpmo(1)/dpmo(4),dpmo(2)/dpmo(4),dpmo(3)/dpmo(4))
1758  n0=nsd(jt)+2
1759  ELSE
1760  CALL pyrobo(nsd(jt)+1,nsd(jt)+3,acos(cthe(jt)),phi(jt),
1761  & dpmo(1)/dpmo(4),dpmo(2)/dpmo(4),dpmo(3)/dpmo(4))
1762  n0=nsd(jt)+3
1763  ENDIF
1764 
1765  DO 610 j=1,4
1766  vdcy(j)=v(id,j)+v(id,5)*p(id,j)/p(id,5)
1767  610 CONTINUE
1768 C...Fill in position of decay vertex.
1769  DO 630 i=nsd(jt)+1,n0
1770  DO 620 j=1,4
1771  v(i,j)=vdcy(j)
1772  620 CONTINUE
1773  v(i,5)=0d0
1774 
1775  630 CONTINUE
1776 CMRENNA--
1777 
1778 C...Mark decayed resonances; trace history.
1779  k(id,1)=k(id,1)+10
1780  kfa=iabs(k(id,2))
1781  kca=pycomp(kfa)
1782  IF(kcqm(jt).NE.0) THEN
1783 C...Do not kill colour flow through coloured resonance!
1784  ELSE
1785  k(id,4)=nsd(jt)+1
1786  k(id,5)=nsd(jt)+2
1787 C...If 3-body or 2-body with junction:
1788  IF(kfl3(jt).NE.0.OR.itjunc(jt).NE.0) k(id,5)=nsd(jt)+3
1789 C...If 3-body with junction:
1790  IF(itjunc(jt).NE.0.AND.kfl3(jt).NE.0) k(id,5)=nsd(jt)+4
1791  ENDIF
1792 
1793 C...Add documentation lines.
1794  isubrg=max(1,min(500,mint(1)))
1795  IF(ires.EQ.0.OR.iset(isubrg).EQ.11) THEN
1796  idoc=mint(83)+mint(4)
1797 CMRENNA+++
1798  ihi=nsd(jt)+2
1799  IF(kfl3(jt).NE.0) ihi=ihi+1
1800  DO 650 i=nsd(jt)+1,ihi
1801 CMRENNA---
1802  i1=mint(83)+mint(4)+1
1803  k(i,3)=i1
1804  IF(mstp(128).GE.1) k(i,3)=id
1805  IF(mstp(128).LE.1.AND.mint(4).LT.mstp(126)) THEN
1806  mint(4)=mint(4)+1
1807  k(i1,1)=21
1808  k(i1,2)=k(i,2)
1809  k(i1,3)=iref(ip,jt+3)
1810  DO 640 j=1,5
1811  p(i1,j)=p(i,j)
1812  640 CONTINUE
1813  ENDIF
1814  650 CONTINUE
1815  ELSE
1816  k(nsd(jt)+1,3)=id
1817  k(nsd(jt)+2,3)=id
1818 C...If 3-body or 2-body with junction:
1819  IF(kfl3(jt).NE.0.OR.itjunc(jt).GT.0) k(nsd(jt)+3,3)=id
1820 C...If 3-body with junction:
1821  IF(kfl3(jt).NE.0.AND.itjunc(jt).GT.0) k(nsd(jt)+4,3)=id
1822  ENDIF
1823 
1824 C...Do showering of two or three objects.
1825  nshbef=n
1826  IF(mstp(71).GE.1.AND.mint(35).LE.1) THEN
1827  IF(kfl3(jt).EQ.0) THEN
1828  CALL pyshow(nsd(jt)+1,nsd(jt)+2,p(id,5))
1829  ELSE
1830  CALL pyshow(nsd(jt)+1,-3,p(id,5))
1831  ENDIF
1832 
1833 c...For pT-ordered shower need set up first, especially colour tags.
1834 C...(Need to set up colour tags even if MSTP(71) = 0)
1835  ELSEIF(mint(35).GE.2) THEN
1836  npart=2
1837  IF(kfl3(jt).NE.0) npart=3
1838  ipart(1)=nsd(jt)+1
1839  ipart(2)=nsd(jt)+2
1840  ipart(3)=nsd(jt)+3
1841  ptpart(1)=0.5d0*p(id,5)
1842  ptpart(2)=ptpart(1)
1843  ptpart(3)=ptpart(1)
1844  IF(kcq1(jt).EQ.1.OR.kcq1(jt).EQ.2) THEN
1845  mother=k(nsd(jt)+1,4)/mstu(5)
1846  IF(mother.LE.nsd(jt)) THEN
1847  mct(nsd(jt)+1,1)=mct(mother,1)
1848  ELSE
1849  nct=nct+1
1850  mct(nsd(jt)+1,1)=nct
1851  mct(mother,2)=nct
1852  ENDIF
1853  ENDIF
1854  IF(kcq1(jt).EQ.-1.OR.kcq1(jt).EQ.2) THEN
1855  mother=k(nsd(jt)+1,5)/mstu(5)
1856  IF(mother.LE.nsd(jt)) THEN
1857  mct(nsd(jt)+1,2)=mct(mother,2)
1858  ELSE
1859  nct=nct+1
1860  mct(nsd(jt)+1,2)=nct
1861  mct(mother,1)=nct
1862  ENDIF
1863  ENDIF
1864  IF(mct(nsd(jt)+2,1).EQ.0.AND.(kcq2(jt).EQ.1.OR.
1865  & kcq2(jt).EQ.2)) THEN
1866  mother=k(nsd(jt)+2,4)/mstu(5)
1867  IF(mother.LE.nsd(jt)) THEN
1868  mct(nsd(jt)+2,1)=mct(mother,1)
1869  ELSE
1870  nct=nct+1
1871  mct(nsd(jt)+2,1)=nct
1872  mct(mother,2)=nct
1873  ENDIF
1874  ENDIF
1875  IF(mct(nsd(jt)+2,2).EQ.0.AND.(kcq2(jt).EQ.-1.OR.
1876  & kcq2(jt).EQ.2)) THEN
1877  mother=k(nsd(jt)+2,5)/mstu(5)
1878  IF(mother.LE.nsd(jt)) THEN
1879  mct(nsd(jt)+2,2)=mct(mother,2)
1880  ELSE
1881  nct=nct+1
1882  mct(nsd(jt)+2,2)=nct
1883  mct(mother,1)=nct
1884  ENDIF
1885  ENDIF
1886  IF(npart.EQ.3.AND.mct(nsd(jt)+3,1).EQ.0.AND.
1887  & (kcq3(jt).EQ.1.OR. kcq3(jt).EQ.2)) THEN
1888  mother=k(nsd(jt)+3,4)/mstu(5)
1889  mct(nsd(jt)+3,1)=mct(mother,1)
1890  ENDIF
1891  IF(npart.EQ.3.AND.mct(nsd(jt)+3,2).EQ.0.AND.
1892  & (kcq3(jt).EQ.-1.OR.kcq3(jt).EQ.2)) THEN
1893  mother=k(nsd(jt)+3,5)/mstu(5)
1894  mct(nsd(jt)+2,2)=mct(mother,2)
1895  ENDIF
1896  IF (mstp(71).GE.1) CALL pyptfs(2,0.5d0*p(id,5),0d0,ptgen)
1897  ENDIF
1898  nshaft=n
1899  IF(jt.EQ.1) naft1=n
1900 
1901 C...Check if decay products moved by shower.
1902  nsd1=nsd(jt)+1
1903  nsd2=nsd(jt)+2
1904  nsd3=nsd(jt)+3
1905  IF(nshaft.GT.nshbef) THEN
1906  IF(k(nsd1,1).GT.10) THEN
1907  DO 660 i=nshbef+1,nshaft
1908  IF(k(i,1).LT.10.AND.k(i,2).EQ.k(nsd1,2)) nsd1=i
1909  660 CONTINUE
1910  ENDIF
1911  IF(k(nsd2,1).GT.10) THEN
1912  DO 670 i=nshbef+1,nshaft
1913  IF(k(i,1).LT.10.AND.k(i,2).EQ.k(nsd2,2).AND.
1914  & i.NE.nsd1) nsd2=i
1915  670 CONTINUE
1916  ENDIF
1917  IF(kfl3(jt).NE.0.AND.k(nsd3,1).GT.10) THEN
1918  DO 680 i=nshbef+1,nshaft
1919  IF(k(i,1).LT.10.AND.k(i,2).EQ.k(nsd3,2).AND.
1920  & i.NE.nsd1.AND.i.NE.nsd2) nsd3=i
1921  680 CONTINUE
1922  ENDIF
1923  ENDIF
1924 
1925 C...Store decay products for further treatment.
1926  np=np+1
1927  iref(np,1)=nsd1
1928  iref(np,2)=nsd2
1929  iref(np,3)=0
1930  IF(kfl3(jt).NE.0) iref(np,3)=nsd3
1931  iref(np,4)=idoc+1
1932  iref(np,5)=idoc+2
1933  iref(np,6)=0
1934  IF(kfl3(jt).NE.0) iref(np,6)=idoc+3
1935  iref(np,7)=k(iref(ip,jt),2)
1936  iref(np,8)=iref(ip,jt)
1937  690 CONTINUE
1938 
1939 
1940 C...Fill information for 2 -> 1 -> 2.
1941  700 IF(jtmax.EQ.1.AND.kdcy(1).NE.0.AND.isub.NE.0) THEN
1942  mint(7)=mint(83)+6+2*iset(isub)
1943  mint(8)=mint(83)+7+2*iset(isub)
1944  mint(25)=kfl1(1)
1945  mint(26)=kfl2(1)
1946  vint(23)=cthe(1)
1947  rm3=p(n-1,5)**2/sh
1948  rm4=p(n,5)**2/sh
1949  be34=sqrt(max(0d0,(1d0-rm3-rm4)**2-4d0*rm3*rm4))
1950  vint(45)=-0.5d0*sh*(1d0-rm3-rm4-be34*cthe(1))
1951  vint(46)=-0.5d0*sh*(1d0-rm3-rm4+be34*cthe(1))
1952  vint(48)=0.25d0*sh*be34**2*max(0d0,1d0-cthe(1)**2)
1953  vint(47)=sqrt(vint(48))
1954  ENDIF
1955 
1956 C...Possibility of colour rearrangement in W+W- events.
1957  IF((isub.EQ.25.OR.isub.EQ.22).AND.mstp(115).GE.1) THEN
1958  iakf1=iabs(kfl1(1))
1959  iakf2=iabs(kfl1(2))
1960  iakf3=iabs(kfl2(1))
1961  iakf4=iabs(kfl2(2))
1962  IF(min(iakf1,iakf2,iakf3,iakf4).GE.1.AND.
1963  & max(iakf1,iakf2,iakf3,iakf4).LE.5) CALL
1964  & pyreco(iref(1,1),iref(1,2),nsd(1),naft1)
1965  IF(mint(51).NE.0) RETURN
1966  ENDIF
1967 
1968 C...Loop back if needed.
1969  710 IF(ip.LT.np) goto 170
1970 
1971 C...Boost back to standard frame.
1972  720 IF(ibst.EQ.1) CALL pyrobo(mint(83)+7,n,thein,phiin,bexin,beyin,
1973  &bezin)
1974 
1975  RETURN
1976  END