PYTHIA/pythia/pyresd.F

   1
   2 C*********************************************************************
   3
   4       SUBROUTINE PYRESD
   5
   6 C...Allows resonances to decay (including parton showers for hadronic
   7 C...channels).
   8       IMPLICIT DOUBLE PRECISION(D)
   9       COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
  10       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
  11       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
  12       COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
  13       COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
  14       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
  15       COMMON/PYINT1/MINT(400),VINT(400)
  16       COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
  17       COMMON/PYINT4/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3)
  18       SAVE /LUJETS/,/LUDAT1/,/LUDAT2/,/LUDAT3/
  19       SAVE /PYSUBS/,/PYPARS/,/PYINT1/,/PYINT2/,/PYINT4/
  20       DIMENSION IREF(20,8),KDCY(3),KFL1(3),KFL2(3),KEQL(3),NSD(3),
  21      &ILIN(6),HGZ(3,3),COUP(6,4),CORL(2,2,2),PK(6,4),PKK(6,6),
  22      &CTHE(3),PHI(3),WDTP(0:40),WDTE(0:40,0:5),DBEZQQ(3),DPMO(5)
  23       COMPLEX FGK,HA(6,6),HC(6,6)
  24
  25 C...The F, Xi and Xj functions of Gunion and Kunszt
  26 C...(Phys. Rev. D33, 665, plus errata from the authors).
  27       FGK(I1,I2,I3,I4,I5,I6)=4.*HA(I1,I3)*HC(I2,I6)*(HA(I1,I5)*
  28      &HC(I1,I4)+HA(I3,I5)*HC(I3,I4))
  29       DIGK(DT,DU)=-4.*D34*D56+DT*(3.*DT+4.*DU)+DT**2*(DT*DU/(D34*D56)-
  30      &2.*(1./D34+1./D56)*(DT+DU)+2.*(D34/D56+D56/D34))
  31       DJGK(DT,DU)=8.*(D34+D56)**2-8.*(D34+D56)*(DT+DU)-6.*DT*DU-
  32      &2.*DT*DU*(DT*DU/(D34*D56)-2.*(1./D34+1./D56)*(DT+DU)+
  33      &2.*(D34/D56+D56/D34))
  34
  35 C...Some general constants.
  36       XW=PARU(102)
  37       XWV=XW
  38       IF(MSTP(8).GE.2) XW=1.-(PMAS(24,1)/PMAS(23,1))**2
  39       XW1=1.-XW
  40       SQMZ=PMAS(23,1)**2
  41       SQMW=PMAS(24,1)**2
  42       SH=VINT(44)
  43
  44 C...Define initial one, two or three objects.
  45       ISUB=MINT(1)
  46       DO 100 JT=1,8
  47       IREF(1,JT)=0
  48   100 CONTINUE
  49       IF(ISET(ISUB).EQ.1.OR.ISET(ISUB).EQ.3) THEN
  50         IREF(1,1)=MINT(84)+2+ISET(ISUB)
  51         IREF(1,4)=MINT(83)+6+ISET(ISUB)
  52       ELSEIF(ISET(ISUB).EQ.2.OR.ISET(ISUB).EQ.4) THEN
  53         IREF(1,1)=MINT(84)+1+ISET(ISUB)
  54         IREF(1,2)=MINT(84)+2+ISET(ISUB)
  55         IREF(1,4)=MINT(83)+5+ISET(ISUB)
  56         IREF(1,5)=MINT(83)+6+ISET(ISUB)
  57       ELSEIF(ISET(ISUB).EQ.5) THEN
  58         IREF(1,1)=MINT(84)+3
  59         IREF(1,2)=MINT(84)+4
  60         IREF(1,3)=MINT(84)+5
  61         IREF(1,4)=MINT(83)+7
  62         IREF(1,5)=MINT(83)+8
  63         IREF(1,6)=MINT(83)+9
  64       ELSEIF(ISET(ISUB).EQ.6) THEN
  65         IREF(1,1)=MINT(84)+4
  66         IREF(1,2)=MINT(84)+5
  67         IREF(1,3)=MINT(84)+3
  68         IREF(1,4)=MINT(83)+8
  69         IREF(1,5)=MINT(83)+9
  70         IREF(1,6)=MINT(83)+7
  71       ENDIF
  72
  73 C...Check if initial resonance has been moved (in resonance + jet).
  74       DO 120 JT=1,3
  75       IF(IREF(1,JT).GT.0) THEN
  76         IF(K(IREF(1,JT),1).GT.10) THEN
  77           KFA=IABS(K(IREF(1,JT),2))
  78           IF(KFA.EQ.6.OR.KFA.EQ.7.OR.KFA.EQ.8.OR.KFA.EQ.39) THEN
  79             DO 110 I=IREF(1,JT)+1,N
  80             IF(K(I,1).LE.10.AND.K(I,2).EQ.K(IREF(1,JT),2)) IREF(1,JT)=I
  81   110       CONTINUE
  82           ELSE
  83             KDA=MOD(K(IREF(1,JT),4),MSTU(4))
  84             IF(KFA.GE.23.AND.KFA.LE.40.AND.KDA.GT.1) IREF(1,JT)=KDA
  85           ENDIF
  86         ENDIF
  87       ENDIF
  88   120 CONTINUE
  89
  90 C...Loop over decay history.
  91       NP=1
  92       IP=0
  93   130 IP=IP+1
  94       NINH=0
  95       JTMAX=2
  96       IF(IP.EQ.1.AND.IREF(1,2).EQ.0) JTMAX=1
  97       IF(IP.EQ.1.AND.IREF(1,3).NE.0) JTMAX=3
  98       ITLH=0
  99       NSAV=N
 100
 101 C...Start treatment of one or two resonances in parallel.
 102   140 N=NSAV
 103       DO 170 JT=1,JTMAX
 104       ID=IREF(IP,JT)
 105       KDCY(JT)=0
 106       KFL1(JT)=0
 107       KFL2(JT)=0
 108       KEQL(JT)=0
 109       NSD(JT)=ID
 110       IF(ID.EQ.0) GOTO 160
 111       KFA=IABS(K(ID,2))
 112       IF((KFA.LT.23.OR.KFA.GT.40).AND.KFA.NE.6.AND.KFA.NE.7.AND.
 113      &KFA.NE.8.AND.KFA.NE.17.AND.KFA.NE.18) GOTO 160
 114       IF(MSTP(48).LE.0.AND.KFA.EQ.6) GOTO 160
 115       IF(MSTP(6).NE.1.AND.MSTP(49).LE.0.AND.(KFA.EQ.7.OR.KFA.EQ.8.OR.
 116      &KFA.EQ.17.OR.KFA.EQ.18)) GOTO 160
 117       IF(K(ID,1).GT.10.OR.MDCY(KFA,1).EQ.0) GOTO 160
 118       IF(KFA.EQ.6.OR.(MSTP(6).NE.1.AND.(KFA.EQ.7.OR.KFA.EQ.8.OR.
 119      &KFA.EQ.17.OR.KFA.EQ.18))) ITLH=ITLH+1
 120       K(ID,4)=MSTU(5)*(K(ID,4)/MSTU(5))
 121       K(ID,5)=MSTU(5)*(K(ID,5)/MSTU(5))
 122
 123 C...Select decay channel.
 124       KFB=0
 125       IF(ISET(ISUB).NE.6.OR.JT.NE.3) THEN
 126         IF(ISUB.EQ.1.OR.ISUB.EQ.15.OR.ISUB.EQ.19.OR.ISUB.EQ.22.OR.
 127      &  ISUB.EQ.30.OR.ISUB.EQ.35.OR.ISUB.EQ.141) MINT(61)=1
 128         CALL PYWIDT(KFA,P(ID,5)**2,WDTP,WDTE)
 129         IF(KCHG(KFA,3).EQ.0) THEN
 130           IPM=2
 131         ELSE
 132           IPM=(5-ISIGN(1,K(ID,2)))/2
 133         ENDIF
 134         IF(JTMAX.GE.2.AND.JT.LE.2) KFB=IABS(K(IREF(IP,3-JT),2))
 135         WDTE0S=WDTE(0,1)+WDTE(0,IPM)+WDTE(0,4)
 136         IF(KFB.EQ.KFA) WDTE0S=WDTE0S+WDTE(0,5)
 137         IF(WDTE0S.LE.0.) THEN
 138           IF(KFA.EQ.6.OR.KFA.EQ.7.OR.KFA.EQ.8.OR.KFA.EQ.17.OR.
 139      &    KFA.EQ.18) THEN
 140             MINT(51)=1
 141             RETURN
 142           ELSE
 143             GOTO 160
 144           ENDIF
 145         ENDIF
 146         RKFL=WDTE0S*RLU(0)
 147         IDL=0
 148   150   IDL=IDL+1
 149         IDC=IDL+MDCY(KFA,2)-1
 150         RKFL=RKFL-(WDTE(IDL,1)+WDTE(IDL,IPM)+WDTE(IDL,4))
 151         IF(KFB.EQ.KFA) RKFL=RKFL-WDTE(IDL,5)
 152         IF(IDL.LT.MDCY(KFA,3).AND.RKFL.GT.0.) GOTO 150
 153       ELSE
 154         IDC=MINT(35)
 155       ENDIF
 156
 157 C...Read out and classify decay channel chosen.
 158       KFL1(JT)=KFDP(IDC,1)*ISIGN(1,K(ID,2))
 159       KFC1A=IABS(KFL1(JT))
 160       IF(KFC1A.GT.100) KFC1A=LUCOMP(KFC1A)
 161       IF(KCHG(KFC1A,3).EQ.0) KFL1(JT)=IABS(KFL1(JT))
 162       KFL2(JT)=KFDP(IDC,2)*ISIGN(1,K(ID,2))
 163       KFC2A=IABS(KFL2(JT))
 164       IF(KFC2A.GT.100) KFC2A=LUCOMP(KFC2A)
 165       IF(KCHG(KFC2A,3).EQ.0) KFL2(JT)=IABS(KFL2(JT))
 166       KDCY(JT)=2
 167       IF(IABS(KFL1(JT)).LE.10.OR.KFL1(JT).EQ.21) KDCY(JT)=1
 168       IF(IABS(KFL1(JT)).GE.23.AND.IABS(KFL1(JT)).LE.40) KDCY(JT)=3
 169       IF(KFB.EQ.KFA) KEQL(JT)=MDME(IDC,1)
 170       NSD(JT)=N
 171       HGZ(JT,1)=VINT(111)
 172       HGZ(JT,2)=VINT(112)
 173       HGZ(JT,3)=VINT(114)
 174
 175 C...Select masses and check that mass sum not too large.
 176       IF(MSTP(42).LE.0.OR.(PMAS(KFC1A,2).LT.PARP(41).AND.
 177      &PMAS(KFC2A,2).LT.PARP(41))) THEN
 178         P(N+1,5)=PMAS(KFC1A,1)
 179         P(N+2,5)=PMAS(KFC2A,1)
 180         IF(P(N+1,5)+P(N+2,5)+PARJ(64).GT.P(ID,5)) THEN
 181           CALL LUERRM(13,'(PYRESD:) daughter masses too large')
 182           MINT(51)=1
 183           RETURN
 184         ENDIF
 185       ELSEIF(IP.EQ.1) THEN
 186         CALL PYOFSH(2,KFA,KFL1(JT),KFL2(JT),P(ID,5),P(N+1,5),P(N+2,5))
 187         IF(MINT(51).EQ.1) RETURN
 188       ELSE
 189         CALL PYOFSH(7,KFA,KFL1(JT),KFL2(JT),P(ID,5),P(N+1,5),P(N+2,5))
 190         IF(MINT(51).EQ.1) RETURN
 191       ENDIF
 192
 193 C...Fill decay products, prepared for parton showers for quarks.
 194 C...Special cases, done by hand, for techni-eta, t, leptoquark and q*.
 195       MSTU(10)=1
 196       IF(KFA.EQ.38.OR.KFA.EQ.39.OR.((MSTP(6).EQ.1.OR.MSTP(49).GE.1).AND.
 197      &(KFA.EQ.7.OR.KFA.EQ.8)).OR.KFA.EQ.6) THEN
 198         MSTU(19)=1
 199         CALL LU2ENT(N+1,KFL1(JT),KFL2(JT),P(ID,5))
 200         ISID=4
 201         IF(K(ID,2).LT.0) ISID=5
 202         IF(KFA.EQ.38) THEN
 203           IF(KFC1A.EQ.21.AND.RLU(0).GT.0.5) ISID=9-ISID
 204           K(N-1,1)=3
 205           K(N,1)=3
 206           K(ID,ISID)=K(ID,ISID)+(N-1)
 207           K(ID,9-ISID)=K(ID,9-ISID)+N
 208           K(N-1,ISID)=MSTU(5)*ID
 209           K(N-1,9-ISID)=MSTU(5)*N
 210           K(N,ISID)=MSTU(5)*(N-1)
 211           K(N,9-ISID)=MSTU(5)*ID
 212         ELSEIF(KFA.EQ.6.OR.(MSTP(6).NE.1.AND.(KFA.EQ.7.OR.KFA.EQ.8)))
 213      &  THEN
 214           K(N-1,1)=1
 215           K(N,1)=3
 216           K(ID,ISID)=K(ID,ISID)+N
 217           K(N,ISID)=MSTU(5)*ID
 218         ELSEIF(KFA.EQ.39) THEN
 219           K(N-1,1)=3
 220           K(N,1)=1
 221           K(ID,ISID)=K(ID,ISID)+(N-1)
 222           K(N-1,ISID)=MSTU(5)*ID
 223         ELSEIF(KFL1(JT).NE.21) THEN
 224           K(N-1,1)=1
 225           K(N,1)=3
 226           K(ID,ISID)=K(ID,ISID)+N
 227           K(N,ISID)=MSTU(5)*ID
 228         ELSE
 229           K(N-1,1)=3
 230           K(N,1)=3
 231           K(ID,ISID)=K(ID,ISID)+(N-1)
 232           K(N-1,ISID)=MSTU(5)*ID
 233           K(N-1,9-ISID)=MSTU(5)*N
 234           K(N,ISID)=MSTU(5)*(N-1)
 235         ENDIF
 236       ELSEIF(KDCY(JT).EQ.1) THEN
 237         CALL LU2ENT(-(N+1),KFL1(JT),KFL2(JT),P(ID,5))
 238       ELSE
 239         CALL LU2ENT(N+1,KFL1(JT),KFL2(JT),P(ID,5))
 240       ENDIF
 241       MSTU(10)=2
 242   160 IF(KFA.GE.23.AND.KFA.LE.40.AND.KFL1(JT).EQ.0) NINH=NINH+1
 243   170 CONTINUE
 244
 245 C...Check for allowed combinations. Skip if no decays.
 246       IF(JTMAX.GE.2) THEN
 247         IF(KEQL(1).EQ.4.AND.KEQL(2).EQ.4) GOTO 140
 248         IF(KEQL(1).EQ.5.AND.KEQL(2).EQ.5) GOTO 140
 249       ENDIF
 250       IF(JTMAX.EQ.1.AND.KDCY(1).EQ.0) GOTO 480
 251       IF(JTMAX.EQ.2.AND.KDCY(1).EQ.0.AND.KDCY(2).EQ.0) GOTO 480
 252       IF(JTMAX.EQ.3.AND.KDCY(1).EQ.0.AND.KDCY(2).EQ.0.AND.
 253      &KDCY(3).EQ.0) GOTO 480
 254
 255 C...Order incoming partons and outgoing resonances.
 256       IF(JTMAX.EQ.2.AND.MSTP(47).GE.1.AND.NINH.EQ.0) THEN
 257         ILIN(1)=MINT(84)+1
 258         IF(K(MINT(84)+1,2).GT.0) ILIN(1)=MINT(84)+2
 259         IF(K(ILIN(1),2).EQ.21) ILIN(1)=2*MINT(84)+3-ILIN(1)
 260         ILIN(2)=2*MINT(84)+3-ILIN(1)
 261         IMIN=1
 262         IF(IREF(IP,7).EQ.25.OR.IREF(IP,7).EQ.35.OR.IREF(IP,7)
 263      &  .EQ.36) IMIN=3
 264         IMAX=2
 265         IORD=1
 266         IF(K(IREF(IP,1),2).EQ.23) IORD=2
 267         IF(K(IREF(IP,1),2).EQ.24.AND.K(IREF(IP,2),2).EQ.-24) IORD=2
 268         IAKIPD=IABS(K(IREF(IP,IORD),2))
 269         IF(IAKIPD.EQ.25.OR.IAKIPD.EQ.35.OR.IAKIPD.EQ.36) IORD=3-IORD
 270         IF(KDCY(IORD).EQ.0) IORD=3-IORD
 271
 272 C...Order decay products of resonances.
 273         DO 180 JT=IORD,3-IORD,3-2*IORD
 274         IF(KDCY(JT).EQ.0) THEN
 275           ILIN(IMAX+1)=NSD(JT)
 276           IMAX=IMAX+1
 277         ELSEIF(K(NSD(JT)+1,2).GT.0) THEN
 278           ILIN(IMAX+1)=N+2*JT-1
 279           ILIN(IMAX+2)=N+2*JT
 280           IMAX=IMAX+2
 281           K(N+2*JT-1,2)=K(NSD(JT)+1,2)
 282           K(N+2*JT,2)=K(NSD(JT)+2,2)
 283         ELSE
 284           ILIN(IMAX+1)=N+2*JT
 285           ILIN(IMAX+2)=N+2*JT-1
 286           IMAX=IMAX+2
 287           K(N+2*JT-1,2)=K(NSD(JT)+1,2)
 288           K(N+2*JT,2)=K(NSD(JT)+2,2)
 289         ENDIF
 290   180   CONTINUE
 291
 292 C...Find charge, isospin, left- and righthanded couplings.
 293         DO 200 I=IMIN,IMAX
 294         DO 190 J=1,4
 295         COUP(I,J)=0.
 296   190   CONTINUE
 297         KFA=IABS(K(ILIN(I),2))
 298         IF(KFA.EQ.0.OR.KFA.GT.20) GOTO 200
 299         COUP(I,1)=KCHG(KFA,1)/3.
 300         COUP(I,2)=(-1)**MOD(KFA,2)
 301         COUP(I,4)=-2.*COUP(I,1)*XWV
 302         COUP(I,3)=COUP(I,2)+COUP(I,4)
 303   200   CONTINUE
 304
 305 C...Full propagator dependence and flavour correlations for 2 gamma*/Z.
 306         IF(ISUB.EQ.22) THEN
 307           DO 230 I=3,5,2
 308           I1=IORD
 309           IF(I.EQ.5) I1=3-IORD
 310           DO 220 J1=1,2
 311           DO 210 J2=1,2
 312           CORL(I/2,J1,J2)=COUP(1,1)**2*HGZ(I1,1)*COUP(I,1)**2/16.+
 313      &    COUP(1,1)*COUP(1,J1+2)*HGZ(I1,2)*COUP(I,1)*COUP(I,J2+2)/4.+
 314      &    COUP(1,J1+2)**2*HGZ(I1,3)*COUP(I,J2+2)**2
 315   210     CONTINUE
 316   220     CONTINUE
 317   230     CONTINUE
 318           COWT12=(CORL(1,1,1)+CORL(1,1,2))*(CORL(2,1,1)+CORL(2,1,2))+
 319      &    (CORL(1,2,1)+CORL(1,2,2))*(CORL(2,2,1)+CORL(2,2,2))
 320           COMX12=(CORL(1,1,1)+CORL(1,1,2)+CORL(1,2,1)+CORL(1,2,2))*
 321      &    (CORL(2,1,1)+CORL(2,1,2)+CORL(2,2,1)+CORL(2,2,2))
 322           IF(COWT12.LT.RLU(0)*COMX12) GOTO 140
 323         ENDIF
 324       ENDIF
 325
 326 C...Select angular orientation type - Z'/W' only.
 327       MZPWP=0
 328       IF(ISUB.EQ.141) THEN
 329         IF(RLU(0).LT.PARU(130)) MZPWP=1
 330         IF(IP.EQ.2) THEN
 331           IF(IABS(K(IREF(2,1),2)).EQ.37) MZPWP=2
 332           IAKIR=IABS(K(IREF(2,2),2))
 333           IF(IAKIR.EQ.25.OR.IAKIR.EQ.35.OR.IAKIR.EQ.36) MZPWP=2
 334         ENDIF
 335         IF(IP.GE.3) MZPWP=2
 336       ELSEIF(ISUB.EQ.142) THEN
 337         IF(RLU(0).LT.PARU(136)) MZPWP=1
 338         IF(IP.EQ.2) THEN
 339           IAKIR=IABS(K(IREF(2,2),2))
 340           IF(IAKIR.EQ.25.OR.IAKIR.EQ.35.OR.IAKIR.EQ.36) MZPWP=2
 341         ENDIF
 342         IF(IP.GE.3) MZPWP=2
 343       ENDIF
 344
 345 C...Select random angles (begin of weighting procedure).
 346   240 DO 250 JT=1,JTMAX
 347       IF(KDCY(JT).EQ.0) GOTO 250
 348       IF(ISET(ISUB).EQ.6.AND.JT.EQ.3) GOTO 250
 349       IF(JTMAX.EQ.1) THEN
 350         CTHE(JT)=VINT(13)+(VINT(33)-VINT(13)+VINT(34)-VINT(14))*RLU(0)
 351         IF(CTHE(JT).GT.VINT(33)) CTHE(JT)=CTHE(JT)+VINT(14)-VINT(33)
 352         PHI(JT)=VINT(24)
 353       ELSE
 354         CTHE(JT)=2.*RLU(0)-1.
 355         PHI(JT)=PARU(2)*RLU(0)
 356       ENDIF
 357   250 CONTINUE
 358
 359       IF(JTMAX.EQ.2.AND.MSTP(47).GE.1.AND.NINH.EQ.0) THEN
 360 C...Construct massless four-vectors.
 361         DO 270 I=N+1,N+4
 362         K(I,1)=1
 363         DO 260 J=1,5
 364         P(I,J)=0.
 365         V(I,J)=0.
 366   260   CONTINUE
 367   270   CONTINUE
 368         DO 280 JT=1,JTMAX
 369         IF(KDCY(JT).EQ.0) GOTO 280
 370         ID=IREF(IP,JT)
 371         P(N+2*JT-1,3)=0.5*P(ID,5)
 372         P(N+2*JT-1,4)=0.5*P(ID,5)
 373         P(N+2*JT,3)=-0.5*P(ID,5)
 374         P(N+2*JT,4)=0.5*P(ID,5)
 375         CALL LUDBRB(N+2*JT-1,N+2*JT,ACOS(CTHE(JT)),PHI(JT),DBLE(P(ID,1)/
 376      &  P(ID,4)),DBLE(P(ID,2)/P(ID,4)),DBLE(P(ID,3)/P(ID,4)))
 377   280   CONTINUE
 378
 379 C...Store incoming and outgoing momenta, with random rotation to
 380 C...avoid accidental zeroes in HA expressions.
 381         DO 300 I=1,IMAX
 382         K(N+4+I,1)=1
 383         P(N+4+I,4)=SQRT(P(ILIN(I),1)**2+P(ILIN(I),2)**2+P(ILIN(I),3)**2+
 384      &  P(ILIN(I),5)**2)
 385         P(N+4+I,5)=P(ILIN(I),5)
 386         DO 290 J=1,3
 387         P(N+4+I,J)=P(ILIN(I),J)
 388   290   CONTINUE
 389   300   CONTINUE
 390   310   THERR=ACOS(2.*RLU(0)-1.)
 391         PHIRR=PARU(2)*RLU(0)
 392         CALL LUDBRB(N+5,N+4+IMAX,THERR,PHIRR,0D0,0D0,0D0)
 393         DO 330 I=1,IMAX
 394         IF(P(N+4+I,1)**2+P(N+4+I,2)**2.LT.1E-4*P(N+4+I,4)**2) GOTO 310
 395         DO 320 J=1,4
 396         PK(I,J)=P(N+4+I,J)
 397   320   CONTINUE
 398   330   CONTINUE
 399
 400 C...Calculate internal products.
 401         IF(ISUB.EQ.22.OR.ISUB.EQ.23.OR.ISUB.EQ.25.OR.ISUB.EQ.141.OR.
 402      &  ISUB.EQ.142) THEN
 403           DO 350 I1=IMIN,IMAX-1
 404           DO 340 I2=I1+1,IMAX
 405           HA(I1,I2)=SQRT((PK(I1,4)-PK(I1,3))*(PK(I2,4)+PK(I2,3))/
 406      &    (1E-20+PK(I1,1)**2+PK(I1,2)**2))*CMPLX(PK(I1,1),PK(I1,2))-
 407      &    SQRT((PK(I1,4)+PK(I1,3))*(PK(I2,4)-PK(I2,3))/
 408      &    (1E-20+PK(I2,1)**2+PK(I2,2)**2))*CMPLX(PK(I2,1),PK(I2,2))
 409           HC(I1,I2)=CONJG(HA(I1,I2))
 410           IF(I1.LE.2) HA(I1,I2)=CMPLX(0.,1.)*HA(I1,I2)
 411           IF(I1.LE.2) HC(I1,I2)=CMPLX(0.,1.)*HC(I1,I2)
 412           HA(I2,I1)=-HA(I1,I2)
 413           HC(I2,I1)=-HC(I1,I2)
 414   340     CONTINUE
 415   350     CONTINUE
 416         ENDIF
 417         DO 370 I=1,2
 418         DO 360 J=1,4
 419         PK(I,J)=-PK(I,J)
 420   360   CONTINUE
 421   370   CONTINUE
 422         DO 390 I1=IMIN,IMAX-1
 423         DO 380 I2=I1+1,IMAX
 424         PKK(I1,I2)=2.*(PK(I1,4)*PK(I2,4)-PK(I1,1)*PK(I2,1)-
 425      &  PK(I1,2)*PK(I2,2)-PK(I1,3)*PK(I2,3))
 426         PKK(I2,I1)=PKK(I1,I2)
 427   380   CONTINUE
 428   390   CONTINUE
 429       ENDIF
 430
 431       KFAGM=IABS(IREF(IP,7))
 432       IF(MSTP(47).LE.0.OR.NINH.NE.0) THEN
 433 C...Isotropic decay selected by user.
 434         WT=1.
 435         WTMAX=1.
 436
 437       ELSEIF(ITLH.GE.1) THEN
 438 C... Isotropic decay t -> b + W etc for 4th generation q and l.
 439         WT=1.
 440         WTMAX=1.
 441
 442       ELSEIF(IREF(IP,7).EQ.25.OR.IREF(IP,7).EQ.35.OR.
 443      &IREF(IP,7).EQ.36) THEN
 444 C...Angular weight for H0 -> Z0 + Z0 or W+ + W- -> 4 quarks/leptons.
 445         WT=16.*PKK(3,5)*PKK(4,6)
 446         IF(IP.EQ.1) WTMAX=SH**2
 447         IF(IP.GE.2) WTMAX=P(IREF(IP,8),5)**4
 448         KFA=IABS(K(IREF(IP,1),2))
 449         IF(KFA.NE.23.AND.KFA.NE.24) WT=WTMAX
 450
 451       ELSEIF((KFAGM.EQ.6.OR.(MSTP(6).NE.1.AND.(KFAGM.EQ.7.OR.
 452      &KFAGM.EQ.8.OR.KFAGM.EQ.17.OR.KFAGM.EQ.18))).AND.
 453      &IABS(K(IREF(IP,1),2)).EQ.24) THEN
 454 C...Angular correlation in f -> f' + W -> f' + 2 quarks/leptons.
 455         I1=IREF(IP,8)
 456         IF(MOD(KFAGM,2).EQ.0) THEN
 457           I2=N+1
 458           I3=N+2
 459         ELSE
 460           I2=N+2
 461           I3=N+1
 462         ENDIF
 463         I4=IREF(IP,2)
 464         WT=(P(I1,4)*P(I2,4)-P(I1,1)*P(I2,1)-P(I1,2)*P(I2,2)-
 465      &  P(I1,3)*P(I2,3))*(P(I3,4)*P(I4,4)-P(I3,1)*P(I4,1)-
 466      &  P(I3,2)*P(I4,2)-P(I3,3)*P(I4,3))
 467         WTMAX=(P(I1,5)**4-P(IREF(IP,1),5)**4)/8.
 468         IF(KFAGM.EQ.6.AND.MSTP(48).LE.1) WT=WTMAX
 469         IF(KFAGM.NE.6.AND.MSTP(49).LE.1) WT=WTMAX
 470
 471       ELSEIF(ISUB.EQ.1) THEN
 472 C...Angular weight for gamma*/Z0 -> 2 quarks/leptons.
 473         EI=KCHG(IABS(MINT(15)),1)/3.
 474         AI=SIGN(1.,EI+0.1)
 475         VI=AI-4.*EI*XWV
 476         EF=KCHG(IABS(KFL1(1)),1)/3.
 477         AF=SIGN(1.,EF+0.1)
 478         VF=AF-4.*EF*XWV
 479         ASYM=2.*(EI*AI*VINT(112)*EF*AF+4.*VI*AI*VINT(114)*VF*AF)/
 480      &  (EI**2*VINT(111)*EF**2+EI*VI*VINT(112)*EF*VF+
 481      &  (VI**2+AI**2)*VINT(114)*(VF**2+AF**2))
 482         WT=1.+ASYM*CTHE(1)*ISIGN(1,MINT(15)*KFL1(1))+CTHE(1)**2
 483         WTMAX=2.+ABS(ASYM)
 484
 485       ELSEIF(ISUB.EQ.2) THEN
 486 C...Angular weight for W+/- -> 2 quarks/leptons.
 487         WT=(1.+CTHE(1)*ISIGN(1,MINT(15)*KFL1(1)))**2
 488         WTMAX=4.
 489
 490       ELSEIF(ISUB.EQ.15.OR.ISUB.EQ.19) THEN
 491 C...Angular weight for f + f~ -> gluon/gamma + (gamma*/Z0) ->
 492 C...-> gluon/gamma + 2 quarks/leptons.
 493         CLILF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
 494      &  COUP(1,1)*COUP(1,3)*HGZ(2,2)*COUP(3,1)*COUP(3,3)/4.+
 495      &  COUP(1,3)**2*HGZ(2,3)*COUP(3,3)**2
 496         CLIRF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
 497      &  COUP(1,1)*COUP(1,3)*HGZ(2,2)*COUP(3,1)*COUP(3,4)/4.+
 498      &  COUP(1,3)**2*HGZ(2,3)*COUP(3,4)**2
 499         CRILF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
 500      &  COUP(1,1)*COUP(1,4)*HGZ(2,2)*COUP(3,1)*COUP(3,3)/4.+
 501      &  COUP(1,4)**2*HGZ(2,3)*COUP(3,3)**2
 502         CRIRF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
 503      &  COUP(1,1)*COUP(1,4)*HGZ(2,2)*COUP(3,1)*COUP(3,4)/4.+
 504      &  COUP(1,4)**2*HGZ(2,3)*COUP(3,4)**2
 505         WT=(CLILF+CRIRF)*(PKK(1,3)**2+PKK(2,4)**2)+
 506      &  (CLIRF+CRILF)*(PKK(1,4)**2+PKK(2,3)**2)
 507         WTMAX=(CLILF+CLIRF+CRILF+CRIRF)*
 508      &  ((PKK(1,3)+PKK(1,4))**2+(PKK(2,3)+PKK(2,4))**2)
 509
 510       ELSEIF(ISUB.EQ.16.OR.ISUB.EQ.20) THEN
 511 C...Angular weight for f + f~' -> gluon/gamma + W+/- ->
 512 C...-> gluon/gamma + 2 quarks/leptons.
 513         WT=PKK(1,3)**2+PKK(2,4)**2
 514         WTMAX=(PKK(1,3)+PKK(1,4))**2+(PKK(2,3)+PKK(2,4))**2
 515
 516       ELSEIF(ISUB.EQ.22) THEN
 517 C...Angular weight for f + f~ -> Z0 + Z0 -> 4 quarks/leptons.
 518         S34=P(IREF(IP,IORD),5)**2
 519         S56=P(IREF(IP,3-IORD),5)**2
 520         TI=PKK(1,3)+PKK(1,4)+S34
 521         UI=PKK(1,5)+PKK(1,6)+S56
 522         FGK135=ABS(FGK(1,2,3,4,5,6)/TI+FGK(1,2,5,6,3,4)/UI)**2
 523         FGK145=ABS(FGK(1,2,4,3,5,6)/TI+FGK(1,2,5,6,4,3)/UI)**2
 524         FGK136=ABS(FGK(1,2,3,4,6,5)/TI+FGK(1,2,6,5,3,4)/UI)**2
 525         FGK146=ABS(FGK(1,2,4,3,6,5)/TI+FGK(1,2,6,5,4,3)/UI)**2
 526         FGK253=ABS(FGK(2,1,5,6,3,4)/TI+FGK(2,1,3,4,5,6)/UI)**2
 527         FGK263=ABS(FGK(2,1,6,5,3,4)/TI+FGK(2,1,3,4,6,5)/UI)**2
 528         FGK254=ABS(FGK(2,1,5,6,4,3)/TI+FGK(2,1,4,3,5,6)/UI)**2
 529         FGK264=ABS(FGK(2,1,6,5,4,3)/TI+FGK(2,1,4,3,6,5)/UI)**2
 530         WT=
 531      &  CORL(1,1,1)*CORL(2,1,1)*FGK135+CORL(1,1,2)*CORL(2,1,1)*FGK145+
 532      &  CORL(1,1,1)*CORL(2,1,2)*FGK136+CORL(1,1,2)*CORL(2,1,2)*FGK146+
 533      &  CORL(1,2,1)*CORL(2,2,1)*FGK253+CORL(1,2,2)*CORL(2,2,1)*FGK263+
 534      &  CORL(1,2,1)*CORL(2,2,2)*FGK254+CORL(1,2,2)*CORL(2,2,2)*FGK264
 535         WTMAX=16.*((CORL(1,1,1)+CORL(1,1,2))*(CORL(2,1,1)+CORL(2,1,2))+
 536      &  (CORL(1,2,1)+CORL(1,2,2))*(CORL(2,2,1)+CORL(2,2,2)))*S34*S56*
 537      &  ((TI**2+UI**2+2.*SH*(S34+S56))/(TI*UI)-S34*S56*(1./TI**2+
 538      &  1./UI**2))
 539
 540       ELSEIF(ISUB.EQ.23) THEN
 541 C...Angular weight for f + f~' -> Z0 + W+/- -> 4 quarks/leptons.
 542         D34=P(IREF(IP,IORD),5)**2
 543         D56=P(IREF(IP,3-IORD),5)**2
 544         DT=PKK(1,3)+PKK(1,4)+D34
 545         DU=PKK(1,5)+PKK(1,6)+D56
 546         FACBW=1./((SH-SQMW)**2+SQMW*PMAS(24,2)**2)
 547         CAWZ=COUP(2,3)/SNGL(DT)-2.*XW1*COUP(1,2)*(SH-SQMW)*FACBW
 548         CBWZ=COUP(1,3)/SNGL(DU)+2.*XW1*COUP(1,2)*(SH-SQMW)*FACBW
 549         FGK135=ABS(CAWZ*FGK(1,2,3,4,5,6)+CBWZ*FGK(1,2,5,6,3,4))
 550         FGK136=ABS(CAWZ*FGK(1,2,3,4,6,5)+CBWZ*FGK(1,2,6,5,3,4))
 551         WT=(COUP(5,3)*FGK135)**2+(COUP(5,4)*FGK136)**2
 552         WTMAX=4.*D34*D56*(COUP(5,3)**2+COUP(5,4)**2)*(CAWZ**2*
 553      &  DIGK(DT,DU)+CBWZ**2*DIGK(DU,DT)+CAWZ*CBWZ*DJGK(DT,DU))
 554
 555       ELSEIF(ISUB.EQ.24.OR.ISUB.EQ.171.OR.ISUB.EQ.176) THEN
 556 C...Angular weight for f + f~ -> Z0 + H0 -> 2 quarks/leptons + H0
 557 C...(or H'0, or A0).
 558         WT=((COUP(1,3)*COUP(3,3))**2+(COUP(1,4)*COUP(3,4))**2)*
 559      &  PKK(1,3)*PKK(2,4)+((COUP(1,3)*COUP(3,4))**2+(COUP(1,4)*
 560      &  COUP(3,3))**2)*PKK(1,4)*PKK(2,3)
 561         WTMAX=(COUP(1,3)**2+COUP(1,4)**2)*(COUP(3,3)**2+COUP(3,4)**2)*
 562      &  (PKK(1,3)+PKK(1,4))*(PKK(2,3)+PKK(2,4))
 563
 564       ELSEIF(ISUB.EQ.25) THEN
 565 C...Angular weight for f + f~ -> W+ + W- -> 4 quarks/leptons.
 566         D34=P(IREF(IP,IORD),5)**2
 567         D56=P(IREF(IP,3-IORD),5)**2
 568         DT=PKK(1,3)+PKK(1,4)+D34
 569         DU=PKK(1,5)+PKK(1,6)+D56
 570         FACBW=1./((SH-SQMZ)**2+SQMZ*PMAS(23,2)**2)
 571         CDWW=(COUP(1,3)*SQMZ*(SH-SQMZ)*FACBW+COUP(1,2))/SH
 572         CAWW=CDWW+0.5*(COUP(1,2)+1.)/SNGL(DT)
 573         CBWW=CDWW+0.5*(COUP(1,2)-1.)/SNGL(DU)
 574         CCWW=COUP(1,4)*SQMZ*(SH-SQMZ)*FACBW/SH
 575         FGK135=ABS(CAWW*FGK(1,2,3,4,5,6)-CBWW*FGK(1,2,5,6,3,4))
 576         FGK253=ABS(FGK(2,1,5,6,3,4)-FGK(2,1,3,4,5,6))
 577         WT=FGK135**2+(CCWW*FGK253)**2
 578         WTMAX=4.*D34*D56*(CAWW**2*DIGK(DT,DU)+CBWW**2*DIGK(DU,DT)-CAWW*
 579      &  CBWW*DJGK(DT,DU)+CCWW**2*(DIGK(DT,DU)+DIGK(DU,DT)-DJGK(DT,DU)))
 580
 581       ELSEIF(ISUB.EQ.26.OR.ISUB.EQ.172.OR.ISUB.EQ.177) THEN
 582 C...Angular weight for f + f~' -> W+/- + H0 -> 2 quarks/leptons + H0
 583 C...(or H'0, or A0).
 584         WT=PKK(1,3)*PKK(2,4)
 585         WTMAX=(PKK(1,3)+PKK(1,4))*(PKK(2,3)+PKK(2,4))
 586
 587       ELSEIF(ISUB.EQ.30.OR.ISUB.EQ.35) THEN
 588 C...Angular weight for f + g/gamma -> f + (gamma*/Z0)
 589 C...-> f + 2 quarks/leptons.
 590         CLILF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
 591      &  COUP(1,1)*COUP(1,3)*HGZ(2,2)*COUP(3,1)*COUP(3,3)/4.+
 592      &  COUP(1,3)**2*HGZ(2,3)*COUP(3,3)**2
 593         CLIRF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
 594      &  COUP(1,1)*COUP(1,3)*HGZ(2,2)*COUP(3,1)*COUP(3,4)/4.+
 595      &  COUP(1,3)**2*HGZ(2,3)*COUP(3,4)**2
 596         CRILF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
 597      &  COUP(1,1)*COUP(1,4)*HGZ(2,2)*COUP(3,1)*COUP(3,3)/4.+
 598      &  COUP(1,4)**2*HGZ(2,3)*COUP(3,3)**2
 599         CRIRF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
 600      &  COUP(1,1)*COUP(1,4)*HGZ(2,2)*COUP(3,1)*COUP(3,4)/4.+
 601      &  COUP(1,4)**2*HGZ(2,3)*COUP(3,4)**2
 602         IF(K(ILIN(1),2).GT.0) WT=(CLILF+CRIRF)*(PKK(1,4)**2+
 603      &  PKK(3,5)**2)+(CLIRF+CRILF)*(PKK(1,3)**2+PKK(4,5)**2)
 604         IF(K(ILIN(1),2).LT.0) WT=(CLILF+CRIRF)*(PKK(1,3)**2+
 605      &  PKK(4,5)**2)+(CLIRF+CRILF)*(PKK(1,4)**2+PKK(3,5)**2)
 606         WTMAX=(CLILF+CLIRF+CRILF+CRIRF)*
 607      &  ((PKK(1,3)+PKK(1,4))**2+(PKK(3,5)+PKK(4,5))**2)
 608
 609       ELSEIF(ISUB.EQ.31) THEN
 610 C...Angular weight for f + g -> f' + W+/- -> f' + 2 quarks/leptons.
 611         IF(K(ILIN(1),2).GT.0) WT=PKK(1,4)**2+PKK(3,5)**2
 612         IF(K(ILIN(1),2).LT.0) WT=PKK(1,3)**2+PKK(4,5)**2
 613         WTMAX=(PKK(1,3)+PKK(1,4))**2+(PKK(3,5)+PKK(4,5))**2
 614
 615       ELSEIF(ISUB.EQ.71.OR.ISUB.EQ.72.OR.ISUB.EQ.73.OR.ISUB.EQ.76.OR.
 616      &ISUB.EQ.77) THEN
 617 C...Angular weight for V_L1 + V_L2 -> V_L3 + V_L4 (V = Z/W).
 618         WT=16.*PKK(3,5)*PKK(4,6)
 619         WTMAX=SH**2
 620
 621       ELSEIF(ISUB.EQ.110) THEN
 622 C...Angular weight for f + f~ -> gamma + H0 -> gamma + X is isotropic.
 623         WT=1.
 624         WTMAX=1.
 625
 626       ELSEIF(ISUB.EQ.141) THEN
 627         IF(IP.EQ.1.AND.(KDCY(1).EQ.1.OR.KDCY(1).EQ.2)) THEN
 628 C...Angular weight for f + f~ -> gamma*/Z0/Z'0 -> 2 quarks/leptons.
 629 C...Couplings of incoming flavour.
 630           KFAI=IABS(MINT(15))
 631           EI=KCHG(KFAI,1)/3.
 632           AI=SIGN(1.,EI+0.1)
 633           VI=AI-4.*EI*XWV
 634           KFAIC=1
 635           IF(KFAI.LE.10.AND.MOD(KFAI,2).EQ.0) KFAIC=2
 636           IF(KFAI.GT.10.AND.MOD(KFAI,2).NE.0) KFAIC=3
 637           IF(KFAI.GT.10.AND.MOD(KFAI,2).EQ.0) KFAIC=4
 638           VPI=PARU(119+2*KFAIC)
 639           API=PARU(120+2*KFAIC)
 640 C...Couplings of final flavour.
 641           KFAF=IABS(KFL1(1))
 642           EF=KCHG(KFAF,1)/3.
 643           AF=SIGN(1.,EF+0.1)
 644           VF=AF-4.*EF*XWV
 645           KFAFC=1
 646           IF(KFAF.LE.10.AND.MOD(KFAF,2).EQ.0) KFAFC=2
 647           IF(KFAF.GT.10.AND.MOD(KFAF,2).NE.0) KFAFC=3
 648           IF(KFAF.GT.10.AND.MOD(KFAF,2).EQ.0) KFAFC=4
 649           VPF=PARU(119+2*KFAFC)
 650           APF=PARU(120+2*KFAFC)
 651 C...Asymmetry and weight.
 652           ASYM=2.*(EI*AI*VINT(112)*EF*AF+EI*API*VINT(113)*EF*APF+
 653      &    4.*VI*AI*VINT(114)*VF*AF+(VI*API+VPI*AI)*VINT(115)*
 654      &    (VF*APF+VPF*AF)+4.*VPI*API*VINT(116)*VPF*APF)/
 655      &    (EI**2*VINT(111)*EF**2+EI*VI*VINT(112)*EF*VF+
 656      &    EI*VPI*VINT(113)*EF*VPF+(VI**2+AI**2)*VINT(114)*
 657      &    (VF**2+AF**2)+(VI*VPI+AI*API)*VINT(115)*(VF*VPF+AF*APF)+
 658      &    (VPI**2+API**2)*VINT(116)*(VPF**2+APF**2))
 659           WT=1.+ASYM*CTHE(1)*ISIGN(1,MINT(15)*KFL1(1))+CTHE(1)**2
 660           WTMAX=2.+ABS(ASYM)
 661         ELSEIF(IP.EQ.1.AND.IABS(KFL1(1)).EQ.24) THEN
 662 C...Angular weight for f + f~ -> Z' -> W+ + W-.
 663           RM1=P(NSD(1)+1,5)**2/SH
 664           RM2=P(NSD(1)+2,5)**2/SH
 665           CCOS2=-(1./16.)*((1.-RM1-RM2)**2-4.*RM1*RM2)*
 666      &    (1.-2.*RM1-2.*RM2+RM1**2+RM2**2+10.*RM1*RM2)
 667           CFLAT=-CCOS2+0.5*(RM1+RM2)*(1.-2.*RM1-2.*RM2+(RM2-RM1)**2)
 668           WT=CFLAT+CCOS2*CTHE(1)**2
 669           WTMAX=CFLAT+MAX(0.,CCOS2)
 670         ELSEIF(IP.EQ.1.AND.(KFL1(1).EQ.25.OR.KFL1(1).EQ.35.OR.
 671      &  IABS(KFL1(1)).EQ.37)) THEN
 672 C...Angular weight for f + f~ -> Z' -> H0 + A0, H'0 + A0, H+ + H-.
 673           WT=1.-CTHE(1)**2
 674           WTMAX=1.
 675         ELSEIF(IP.EQ.1.AND.KDCY(1).EQ.3) THEN
 676 C...Angular weight for f + f~ -> Z' -> Z0 + H0.
 677           RM1=P(NSD(1)+1,5)**2/SH
 678           RM2=P(NSD(1)+2,5)**2/SH
 679           FLAM2=MAX(0.,(1.-RM1-RM2)**2-4.*RM1*RM2)
 680           WT=1.+FLAM2*(1.-CTHE(1)**2)/(8.*RM1)
 681           WTMAX=1.+FLAM2/(8.*RM1)
 682         ELSEIF(MZPWP.EQ.0) THEN
 683 C...Angular weight for f + f~ -> Z' -> W+ + W- -> 4 quarks/leptons
 684 C...(W:s like if intermediate Z).
 685           D34=P(IREF(IP,IORD),5)**2
 686           D56=P(IREF(IP,3-IORD),5)**2
 687           DT=PKK(1,3)+PKK(1,4)+D34
 688           DU=PKK(1,5)+PKK(1,6)+D56
 689           FGK135=ABS(FGK(1,2,3,4,5,6)-FGK(1,2,5,6,3,4))
 690           FGK253=ABS(FGK(2,1,5,6,3,4)-FGK(2,1,3,4,5,6))
 691           WT=(COUP(1,3)*FGK135)**2+(COUP(1,4)*FGK253)**2
 692           WTMAX=4.*D34*D56*(COUP(1,3)**2+COUP(1,4)**2)*
 693      &    (DIGK(DT,DU)+DIGK(DU,DT)-DJGK(DT,DU))
 694         ELSEIF(MZPWP.EQ.1) THEN
 695 C...Angular weight for f + f~ -> Z' -> W+ + W- -> 4 quarks/leptons
 696 C...(W:s approximately longitudinal, like if intermediate H).
 697           WT=16.*PKK(3,5)*PKK(4,6)
 698           WTMAX=SH**2
 699         ELSE
 700 C...Angular weight for f + f~ -> Z' -> H+ + H-, Z0 + H0, H0 + A0,
 701 C...H'0 + A0 -> 4 quarks/leptons.
 702           WT=1.
 703           WTMAX=1.
 704         ENDIF
 705
 706       ELSEIF(ISUB.EQ.142) THEN
 707         IF(IP.EQ.1.AND.(KDCY(1).EQ.1.OR.KDCY(1).EQ.2)) THEN
 708 C...Angular weight for f + f~' -> W'+/- -> 2 quarks/leptons.
 709           KFAI=IABS(MINT(15))
 710           KFAIC=1
 711           IF(KFAI.GT.10) KFAIC=2
 712           VI=PARU(129+2*KFAIC)
 713           AI=PARU(130+2*KFAIC)
 714           KFAF=IABS(KFL1(1))
 715           KFAFC=1
 716           IF(KFAF.GT.10) KFAFC=2
 717           VF=PARU(129+2*KFAFC)
 718           AF=PARU(130+2*KFAFC)
 719           ASYM=8.*VI*AI*VF*AF/((VI**2+AI**2)*(VF**2+AF**2))
 720           WT=1.+ASYM*CTHE(1)*ISIGN(1,MINT(15)*KFL1(1))+CTHE(1)**2
 721           WTMAX=2.+ABS(ASYM)
 722         ELSEIF(IP.EQ.1.AND.IABS(KFL2(1)).EQ.23) THEN
 723 C...Angular weight for f + f~' -> W'+/- -> W+/- + Z0.
 724           RM1=P(NSD(1)+1,5)**2/SH
 725           RM2=P(NSD(1)+2,5)**2/SH
 726           CCOS2=-(1./16.)*((1.-RM1-RM2)**2-4.*RM1*RM2)*
 727      &    (1.-2.*RM1-2.*RM2+RM1**2+RM2**2+10.*RM1*RM2)
 728           CFLAT=-CCOS2+0.5*(RM1+RM2)*(1.-2.*RM1-2.*RM2+(RM2-RM1)**2)
 729           WT=CFLAT+CCOS2*CTHE(1)**2
 730           WTMAX=CFLAT+MAX(0.,CCOS2)
 731         ELSEIF(IP.EQ.1.AND.KDCY(1).EQ.3) THEN
 732 C...Angular weight for f + f~ -> W'+/- -> W+/- + H0.
 733           RM1=P(NSD(1)+1,5)**2/SH
 734           RM2=P(NSD(1)+2,5)**2/SH
 735           FLAM2=MAX(0.,(1.-RM1-RM2)**2-4.*RM1*RM2)
 736           WT=1.+FLAM2*(1.-CTHE(1)**2)/(8.*RM1)
 737           WTMAX=1.+FLAM2/(8.*RM1)
 738         ELSEIF(MZPWP.EQ.0) THEN
 739 C...Angular weight for f + f~' -> W' -> W + Z0 -> 4 quarks/leptons
 740 C...(W/Z like if intermediate W).
 741           D34=P(IREF(IP,IORD),5)**2
 742           D56=P(IREF(IP,3-IORD),5)**2
 743           DT=PKK(1,3)+PKK(1,4)+D34
 744           DU=PKK(1,5)+PKK(1,6)+D56
 745           FGK135=ABS(FGK(1,2,3,4,5,6)-FGK(1,2,5,6,3,4))
 746           FGK136=ABS(FGK(1,2,3,4,6,5)-FGK(1,2,6,5,3,4))
 747           WT=(COUP(5,3)*FGK135)**2+(COUP(5,4)*FGK136)**2
 748           WTMAX=4.*D34*D56*(COUP(5,3)**2+COUP(5,4)**2)*
 749      &    (DIGK(DT,DU)+DIGK(DU,DT)-DJGK(DT,DU))
 750         ELSEIF(MZPWP.EQ.1) THEN
 751 C...Angular weight for f + f~' -> W' -> W + Z0 -> 4 quarks/leptons
 752 C...(W/Z approximately longitudinal, like if intermediate H).
 753           WT=16.*PKK(3,5)*PKK(4,6)
 754           WTMAX=SH**2
 755         ELSE
 756 C...Angular weight for f + f~ -> W' -> W + H0 -> whatever.
 757           WT=1.
 758           WTMAX=1.
 759         ENDIF
 760
 761       ELSEIF(ISUB.EQ.145.OR.ISUB.EQ.162.OR.ISUB.EQ.163.OR.ISUB.EQ.164)
 762      &THEN
 763 C...Isotropic decay of leptoquarks (assumed spin 0).
 764         WT=1.
 765         WTMAX=1.
 766
 767       ELSEIF(ISUB.EQ.147.OR.ISUB.EQ.148) THEN
 768 C...Decays of (spin 1/2) q* -> q + (g,gamma) or (Z0,W+-).
 769         SIDE=1.
 770         IF(MINT(16).EQ.21) SIDE=-1.
 771         IF(IP.EQ.1.AND.(KFL1(1).EQ.21.OR.KFL1(1).EQ.22)) THEN
 772           WT=1.+SIDE*CTHE(1)
 773           WTMAX=2.
 774         ELSEIF(IP.EQ.1) THEN
 775           RM1=P(NSD(1)+1,5)**2/SH
 776           WT=1.+SIDE*CTHE(1)*(1.-0.5*RM1)/(1.+0.5*RM1)
 777           WTMAX=1.+(1.-0.5*RM1)/(1.+0.5*RM1)
 778         ELSE
 779 C...W/Z decay assumed isotropic, since not known.
 780           WT=1.
 781           WTMAX=1.
 782         ENDIF
 783
 784       ELSEIF(ISUB.EQ.149) THEN
 785 C...Isotropic decay of techni-eta.
 786         WT=1.
 787         WTMAX=1.
 788
 789 C...Obtain correct angular distribution by rejection techniques.
 790       ELSE
 791         WT=1.
 792         WTMAX=1.
 793       ENDIF
 794       IF(WT.LT.RLU(0)*WTMAX) GOTO 240
 795
 796 C...Construct massive four-vectors using angles chosen. Mark decayed
 797 C...resonances, add documentation lines. Shower evolution.
 798   400 DO 470 JT=1,JTMAX
 799       IF(KDCY(JT).EQ.0) GOTO 470
 800       ID=IREF(IP,JT)
 801       IF(ISET(ISUB).NE.6.OR.JT.NE.3) THEN
 802         DO 410 J=1,5
 803         DPMO(J)=P(ID,J)
 804   410   CONTINUE
 805         DPMO(4)=SQRT(DPMO(1)**2+DPMO(2)**2+DPMO(3)**2+DPMO(5)**2)
 806         CALL LUDBRB(NSD(JT)+1,NSD(JT)+2,ACOS(CTHE(JT)),PHI(JT),
 807      &  DPMO(1)/DPMO(4),DPMO(2)/DPMO(4),DPMO(3)/DPMO(4))
 808       ELSE
 809 C...Z + q + q~ : angles already known, in rest frame of system.
 810         DO 420 J=1,3
 811         DBEZQQ(J)=(P(ID,J)+P(ID+1,J)+P(ID+2,J))/(P(ID,4)+P(ID+1,4)+
 812      &  P(ID+2,4))
 813   420   CONTINUE
 814         K(N+1,1)=1
 815         DO 430 J=1,5
 816         P(N+1,J)=P(ID,J)
 817   430   CONTINUE
 818         CALL LUDBRB(N+1,N+1,0.,0.,-DBEZQQ(1),-DBEZQQ(2),-DBEZQQ(3))
 819         PHIZQQ=ULANGL(P(N+1,1),P(N+1,2))
 820         THEZQQ=ULANGL(P(N+1,3),SQRT(P(N+1,1)**2+P(N+1,2)**2))
 821         CALL LUDBRB(NSD(JT)+1,NSD(JT)+2,ACOS(VINT(81)),VINT(82),
 822      &  0D0,0D0,DBLE(SQRT(P(N+1,1)**2+P(N+1,2)**2+P(N+1,3)**2)/
 823      &  P(N+1,4)))
 824         CALL LUDBRB(NSD(JT)+1,NSD(JT)+2,THEZQQ,PHIZQQ,DBEZQQ(1),
 825      &  DBEZQQ(2),DBEZQQ(3))
 826       ENDIF
 827       K(ID,1)=K(ID,1)+10
 828       KFA=IABS(K(ID,2))
 829       IF(KFA.EQ.38.OR.KFA.EQ.39.OR.((MSTP(6).EQ.1.OR.MSTP(49).GE.1).AND.
 830      &(KFA.EQ.7.OR.KFA.EQ.8)).OR.(MSTP(48).GE.1.AND.KFA.EQ.6)) THEN
 831 C...Do not kill colour flow through techni-eta, t, leptoquark or q*!
 832       ELSE
 833         K(ID,4)=NSD(JT)+1
 834         K(ID,5)=NSD(JT)+2
 835       ENDIF
 836       IDOC=MINT(83)+MINT(4)
 837       DO 450 I=NSD(JT)+1,NSD(JT)+2
 838       I1=MINT(83)+MINT(4)+1
 839       K(I,3)=I1
 840       IF(MSTP(128).GE.1) K(I,3)=ID
 841       IF(MSTP(128).LE.1.AND.MINT(4).LT.MSTP(126)) THEN
 842         MINT(4)=MINT(4)+1
 843         K(I1,1)=21
 844         K(I1,2)=K(I,2)
 845         K(I1,3)=IREF(IP,JT+3)
 846         DO 440 J=1,5
 847         P(I1,J)=P(I,J)
 848   440   CONTINUE
 849       ENDIF
 850   450 CONTINUE
 851 C...Shower - top currently special case.
 852       NSHBEF=N
 853       IF(MSTP(71).GE.1.AND.(KDCY(JT).LE.2.OR.KFA.EQ.6.OR.KFA.EQ.7.OR.
 854      &KFA.EQ.8)) CALL LUSHOW(NSD(JT)+1,NSD(JT)+2,P(ID,5))
 855       NSHAFT=N
 856
 857 C...Check if new resonances were produced.
 858       KNSDA1=IABS(K(NSD(JT)+1,2))
 859       KNSDA2=IABS(K(NSD(JT)+2,2))
 860       IF(KNSDA1.EQ.6.OR.KNSDA2.EQ.6.OR.KNSDA1.EQ.7.OR.KNSDA2.EQ.7.OR.
 861      &KNSDA1.EQ.8.OR.KNSDA2.EQ.8) THEN
 862         NSD1=0
 863         NSD2=0
 864         DO 460 I=NSD(JT)+1,N
 865         IF(K(I,1).LT.10.AND.K(I,2).EQ.K(NSD(JT)+1,2)) NSD1=I
 866         IF(K(I,1).LT.10.AND.K(I,2).EQ.K(NSD(JT)+2,2)) NSD2=I
 867   460   CONTINUE
 868         IF(NSD1.NE.0.AND.NSD2.NE.0) THEN
 869           NP=NP+1
 870           IREF(NP,1)=NSD1
 871           IREF(NP,2)=NSD2
 872           IREF(NP,3)=0
 873           IREF(NP,4)=IDOC+1
 874           IREF(NP,5)=IDOC+2
 875           IREF(NP,6)=0
 876           IREF(NP,7)=K(IREF(IP,JT),2)
 877           IREF(NP,8)=IREF(IP,JT)
 878         ENDIF
 879       ELSEIF(KDCY(JT).EQ.3.OR.KFA.EQ.6.OR.KFA.EQ.7.OR.KFA.EQ.8) THEN
 880         NP=NP+1
 881         IREF(NP,1)=NSD(JT)+1
 882         IREF(NP,2)=NSD(JT)+2
 883         IF(NSHAFT-NSHBEF.GT.0) THEN
 884           IREF(NP,1)=NSHBEF+2
 885           IREF(NP,2)=NSHBEF+3
 886         ENDIF
 887         IREF(NP,3)=0
 888         IREF(NP,4)=IDOC+1
 889         IREF(NP,5)=IDOC+2
 890         IREF(NP,6)=0
 891         IREF(NP,7)=K(IREF(IP,JT),2)
 892         IREF(NP,8)=IREF(IP,JT)
 893       ENDIF
 894   470 CONTINUE
 895
 896 C...Fill information for 2 -> 1 -> 2. Loop back if needed.
 897       IF(JTMAX.EQ.1.AND.KDCY(1).NE.0) THEN
 898         MINT(7)=MINT(83)+6+2*ISET(ISUB)
 899         MINT(8)=MINT(83)+7+2*ISET(ISUB)
 900         MINT(25)=KFL1(1)
 901         MINT(26)=KFL2(1)
 902         VINT(23)=CTHE(1)
 903         RM3=P(N-1,5)**2/SH
 904         RM4=P(N,5)**2/SH
 905         BE34=SQRT(MAX(0.,(1.-RM3-RM4)**2-4.*RM3*RM4))
 906         VINT(45)=-0.5*SH*(1.-RM3-RM4-BE34*CTHE(1))
 907         VINT(46)=-0.5*SH*(1.-RM3-RM4+BE34*CTHE(1))
 908         VINT(48)=0.25*SH*BE34**2*MAX(0.,1.-CTHE(1)**2)
 909         VINT(47)=SQRT(VINT(48))
 910       ENDIF
 911   480 IF(IP.LT.NP) GOTO 130
 912
 913       RETURN
 914       END