pythiaeRHIC/pythia/pyresd.f

0001
0002
0003 C*********************************************************************
0004
0005 C...PYRESD
0006 C...Allows resonances to decay (including parton showers for hadronic
0007 C...channels).
0008
0009       SUBROUTINE PYRESD(IRES)
0010
0011 C...Double precision and integer declarations.
0012       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0013       IMPLICIT INTEGER(I-N)
0014       INTEGER PYK,PYCHGE,PYCOMP
0015 C...Parameter statement to help give large particle numbers.
0016       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000,
0017      &KEXCIT=4000000,KDIMEN=5000000)
0018 C...Parameter statement for maximum size of showers.
0019       PARAMETER (MAXNUR=1000)
0020 C...Commonblocks.
0021       COMMON/PYPART/NPART,NPARTD,IPART(MAXNUR),PTPART(MAXNUR)
0022       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0023       COMMON/PYCTAG/NCT,MCT(4000,2)
0024       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0025       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0026       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
0027       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
0028       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0029       COMMON/PYINT1/MINT(400),VINT(400)
0030       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
0031       COMMON/PYINT4/MWID(500),WIDS(500,5)
0032       SAVE /PYPART/,/PYJETS/,/PYCTAG/,/PYDAT1/,/PYDAT2/,/PYDAT3/,
0033      &/PYSUBS/,/PYPARS/,/PYINT1/,/PYINT2/,/PYINT4/
0034 C...Local arrays and complex and character variables.
0035       DIMENSION IREF(50,8),KDCY(3),KFL1(3),KFL2(3),KFL3(3),KEQL(3),
0036      &KCQM(3),KCQ1(3),KCQ2(3),KCQ3(3),NSD(3),PMMN(3),ILIN(6),
0037      &HGZ(3,3),COUP(6,4),CORL(2,2,2),PK(6,4),PKK(6,6),CTHE(3),
0038      &PHI(3),WDTP(0:400),WDTE(0:400,0:5),DPMO(5),XM(5),VDCY(4),
0039      &ITJUNC(3),CTM2(3)
0040       COMPLEX FGK,HA(6,6),HC(6,6)
0041       REAL TIR,UIR
0042       CHARACTER CODE*9,MASS*9
0043
0044 C...The F, Xi and Xj functions of Gunion and Kunszt
0045 C...(Phys. Rev. D33, 665, plus errata from the authors).
0046       FGK(I1,I2,I3,I4,I5,I6)=4.*HA(I1,I3)*HC(I2,I6)*(HA(I1,I5)*
0047      &HC(I1,I4)+HA(I3,I5)*HC(I3,I4))
0048       DIGK(DT,DU)=-4D0*D34*D56+DT*(3D0*DT+4D0*DU)+DT**2*(DT*DU/
0049      &(D34*D56)-2D0*(1D0/D34+1D0/D56)*(DT+DU)+2D0*(D34/D56+D56/D34))
0050       DJGK(DT,DU)=8D0*(D34+D56)**2-8D0*(D34+D56)*(DT+DU)-6D0*DT*DU-
0051      &2D0*DT*DU*(DT*DU/(D34*D56)-2D0*(1D0/D34+1D0/D56)*(DT+DU)+
0052      &2D0*(D34/D56+D56/D34))
0053
0054 C...Some general constants.
0055       XW=PARU(102)
0056       XWV=XW
0057       IF(MSTP(8).GE.2) XW=1D0-(PMAS(24,1)/PMAS(23,1))**2
0058       XW1=1D0-XW
0059       SQMZ=PMAS(23,1)**2
0060
0061       GMMZ=PMAS(23,1)*PMAS(23,2)
0062       SQMW=PMAS(24,1)**2
0063       GMMW=PMAS(24,1)*PMAS(24,2)
0064       SH=VINT(44)
0065
0066 C...Boost and rotate to rest frame of incoming partons,
0067 C...to get proper amount of smearing of decay angles.
0068       IBST=0
0069       IF(IRES.EQ.0) THEN
0070         IBST=1
0071         ETOTIN=P(MINT(84)+1,4)+P(MINT(84)+2,4)
0072         BEXIN=(P(MINT(84)+1,1)+P(MINT(84)+2,1))/ETOTIN
0073         BEYIN=(P(MINT(84)+1,2)+P(MINT(84)+2,2))/ETOTIN
0074         BEZIN=(P(MINT(84)+1,3)+P(MINT(84)+2,3))/ETOTIN
0075         CALL PYROBO(MINT(83)+7,N,0D0,0D0,-BEXIN,-BEYIN,-BEZIN)
0076         PHIIN=PYANGL(P(MINT(84)+1,1),P(MINT(84)+1,2))
0077         CALL PYROBO(MINT(83)+7,N,0D0,-PHIIN,0D0,0D0,0D0)
0078         THEIN=PYANGL(P(MINT(84)+1,3),P(MINT(84)+1,1))
0079         CALL PYROBO(MINT(83)+7,N,-THEIN,0D0,0D0,0D0,0D0)
0080       ENDIF
0081
0082 C...Reset original resonance configuration.
0083       DO 100 JT=1,8
0084         IREF(1,JT)=0
0085   100 CONTINUE
0086
0087 C...Define initial one, two or three objects for subprocess.
0088       IHDEC=0
0089       IF(IRES.EQ.0) THEN
0090         ISUB=MINT(1)
0091         IF(ISET(ISUB).EQ.1.OR.ISET(ISUB).EQ.3) THEN
0092           IREF(1,1)=MINT(84)+2+ISET(ISUB)
0093           IREF(1,4)=MINT(83)+6+ISET(ISUB)
0094           JTMAX=1
0095         ELSEIF(ISET(ISUB).EQ.2.OR.ISET(ISUB).EQ.4) THEN
0096           IREF(1,1)=MINT(84)+1+ISET(ISUB)
0097           IREF(1,2)=MINT(84)+2+ISET(ISUB)
0098           IREF(1,4)=MINT(83)+5+ISET(ISUB)
0099           IREF(1,5)=MINT(83)+6+ISET(ISUB)
0100           JTMAX=2
0101         ELSEIF(ISET(ISUB).EQ.5) THEN
0102           IREF(1,1)=MINT(84)+3
0103           IREF(1,2)=MINT(84)+4
0104           IREF(1,3)=MINT(84)+5
0105           IREF(1,4)=MINT(83)+7
0106           IREF(1,5)=MINT(83)+8
0107           IREF(1,6)=MINT(83)+9
0108           JTMAX=3
0109         ENDIF
0110
0111 C...Define original resonance for odd cases.
0112       ELSE
0113         ISUB=0
0114         IF(K(IRES,2).EQ.25.OR.K(IRES,2).EQ.35.OR.K(IRES,2).EQ.36)
0115      &  IHDEC=1
0116         IF(IHDEC.EQ.1) ISUB=3
0117         IREF(1,1)=IRES
0118         IREF(1,4)=K(IRES,3)
0119         IRESTM=IRES
0120         IF(IREF(1,4).GT.MINT(84)) THEN
0121   110     ITMPMO=IREF(1,4)
0122           IF(K(ITMPMO,2).EQ.94) THEN
0123             IREF(1,4)=K(ITMPMO,3)+(IRESTM-ITMPMO-1)
0124             IF(K(IREF(1,4),3).LE.MINT(84)) IREF(1,4)=K(IREF(1,4),3)
0125           ELSEIF(K(ITMPMO,2).EQ.K(IRES,2)) THEN
0126             IRESTM=ITMPMO
0127 C...Explicitly check that reference particle exists, otherwise stop recursion
0128             IF(ITMPMO.GT.0.AND.K(ITMPMO,3).GT.0) THEN
0129               IREF(1,4)=K(ITMPMO,3)
0130               GOTO 110
0131             ENDIF
0132           ENDIF
0133         ENDIF
0134         IF(IREF(1,4).GT.MINT(84)) THEN
0135           EMATCH=1D10
0136           IREF14=IREF(1,4)
0137           DO 120 II=MINT(83)+7,MINT(83)+MINT(4)
0138             IF(K(II,2).EQ.K(IRES,2).AND.ABS(P(II,4)-P(IREF14,4)).LT.
0139      &      EMATCH) THEN
0140               IREF(1,4)=II
0141               EMATCH=ABS(P(II,4)-P(IREF14,4))
0142             ENDIF
0143   120     CONTINUE
0144         ENDIF
0145         JTMAX=1
0146       ENDIF
0147
0148 C...Check if initial resonance has been moved (in resonance + jet).
0149       DO 140 JT=1,3
0150         IF(IREF(1,JT).GT.0) THEN
0151           IF(K(IREF(1,JT),1).GT.10) THEN
0152             KFA=IABS(K(IREF(1,JT),2))
0153             IF(KFA.GE.6.AND.KCHG(PYCOMP(KFA),2).NE.0) THEN
0154               KDA1=MOD(K(IREF(1,JT),4),MSTU(5))
0155               KDA2=MOD(K(IREF(1,JT),5),MSTU(5))
0156               IF(KDA1.GT.IREF(1,JT).AND.KDA1.LE.N) THEN
0157                 IF(K(KDA1,2).EQ.21) KDA1=K(KDA1,5)/MSTU(5)
0158               ENDIF
0159               IF(KDA2.GT.IREF(1,JT).AND.KDA2.LE.N) THEN
0160                 IF(K(KDA2,2).EQ.21) KDA2=K(KDA2,4)/MSTU(5)
0161               ENDIF
0162               DO 130 I=IREF(1,JT)+1,N
0163                 IF(K(I,2).EQ.K(IREF(1,JT),2).AND.(I.EQ.KDA1.OR.
0164      &          I.EQ.KDA2)) THEN
0165                   IREF(1,JT)=I
0166                   KDA1=MOD(K(IREF(1,JT),4),MSTU(5))
0167                   KDA2=MOD(K(IREF(1,JT),5),MSTU(5))
0168                   IF(KDA1.GT.IREF(1,JT).AND.KDA1.LE.N) THEN
0169                     IF(K(KDA1,2).EQ.21) KDA1=K(KDA1,5)/MSTU(5)
0170                   ENDIF
0171                   IF(KDA2.GT.IREF(1,JT).AND.KDA2.LE.N) THEN
0172                     IF(K(KDA2,2).EQ.21) KDA2=K(KDA2,4)/MSTU(5)
0173                   ENDIF
0174                 ENDIF
0175   130         CONTINUE
0176             ELSE
0177               KDA=MOD(K(IREF(1,JT),4),MSTU(5))
0178               IF(MWID(PYCOMP(KFA)).NE.0.AND.KDA.GT.1) IREF(1,JT)=KDA
0179             ENDIF
0180           ENDIF
0181         ENDIF
0182   140 CONTINUE
0183
0184 C...Set decay vertex for initial resonances
0185       DO 160 JT=1,JTMAX
0186         DO 150 I=1,4
0187           V(IREF(1,JT),I)=0D0
0188   150   CONTINUE
0189   160 CONTINUE
0190
0191 C...Loop over decay history.
0192       NP=1
0193       IP=0
0194   170 IP=IP+1
0195       NINH=0
0196       JTMAX=2
0197       IF(IREF(IP,2).EQ.0) JTMAX=1
0198       IF(IREF(IP,3).NE.0) JTMAX=3
0199       IT4=0
0200       NSAV=N
0201
0202 C...Check for Higgs which appears as decay product of user-process.
0203       IF(ISUB.EQ.0) THEN
0204         IHDEC=0
0205         IF(IREF(IP,7).EQ.25.OR.IREF(IP,7).EQ.35.OR.IREF(IP,7)
0206      &  .EQ.36) IHDEC=1
0207         IF(IHDEC.EQ.1) ISUB=3
0208       ENDIF
0209
0210 C...Start treatment of one, two or three resonances in parallel.
0211   180 N=NSAV
0212       DO 340 JT=1,JTMAX
0213         ID=IREF(IP,JT)
0214         KDCY(JT)=0
0215         KFL1(JT)=0
0216         KFL2(JT)=0
0217         KFL3(JT)=0
0218         KEQL(JT)=0
0219         NSD(JT)=ID
0220         ITJUNC(JT)=0
0221
0222 C...Check whether particle can/is allowed to decay.
0223         IF(ID.EQ.0) GOTO 330
0224         KFA=IABS(K(ID,2))
0225         KCA=PYCOMP(KFA)
0226         IF(MWID(KCA).EQ.0) GOTO 330
0227         IF(K(ID,1).GT.10.OR.MDCY(KCA,1).EQ.0) GOTO 330
0228         IF(KFA.EQ.6.OR.KFA.EQ.7.OR.KFA.EQ.8.OR.KFA.EQ.17.OR.
0229      &  KFA.EQ.18) IT4=IT4+1
0230         K(ID,4)=MSTU(5)*(K(ID,4)/MSTU(5))
0231         K(ID,5)=MSTU(5)*(K(ID,5)/MSTU(5))
0232
0233 C...Choose lifetime and determine decay vertex.
0234         IF(K(ID,1).EQ.5) THEN
0235           V(ID,5)=0D0
0236         ELSEIF(K(ID,1).NE.4) THEN
0237           V(ID,5)=-PMAS(KCA,4)*LOG(PYR(0))
0238         ENDIF
0239         DO 190 J=1,4
0240           VDCY(J)=V(ID,J)+V(ID,5)*P(ID,J)/P(ID,5)
0241   190   CONTINUE
0242
0243 C...Determine whether decay allowed or not.
0244         MOUT=0
0245         IF(MSTJ(22).EQ.2) THEN
0246           IF(PMAS(KCA,4).GT.PARJ(71)) MOUT=1
0247         ELSEIF(MSTJ(22).EQ.3) THEN
0248           IF(VDCY(1)**2+VDCY(2)**2+VDCY(3)**2.GT.PARJ(72)**2) MOUT=1
0249         ELSEIF(MSTJ(22).EQ.4) THEN
0250           IF(VDCY(1)**2+VDCY(2)**2.GT.PARJ(73)**2) MOUT=1
0251           IF(ABS(VDCY(3)).GT.PARJ(74)) MOUT=1
0252         ENDIF
0253         IF(MOUT.EQ.1.AND.K(ID,1).NE.5) THEN
0254           K(ID,1)=4
0255           GOTO 330
0256         ENDIF
0257
0258 C...Info for selection of decay channel: sign, pairings.
0259         IF(KCHG(KCA,3).EQ.0) THEN
0260           IPM=2
0261         ELSE
0262           IPM=(5-ISIGN(1,K(ID,2)))/2
0263         ENDIF
0264         KFB=0
0265         IF(JTMAX.EQ.2) THEN
0266           KFB=IABS(K(IREF(IP,3-JT),2))
0267         ELSEIF(JTMAX.EQ.3) THEN
0268           JT2=JT+1-3*(JT/3)
0269           KFB=IABS(K(IREF(IP,JT2),2))
0270           IF(KFB.NE.KFA) THEN
0271             JT2=JT+2-3*((JT+1)/3)
0272             KFB=IABS(K(IREF(IP,JT2),2))
0273           ENDIF
0274         ENDIF
0275
0276 C...Select decay channel.
0277         IF(ISUB.EQ.1.OR.ISUB.EQ.15.OR.ISUB.EQ.19.OR.ISUB.EQ.22.OR.
0278      &  ISUB.EQ.30.OR.ISUB.EQ.35.OR.ISUB.EQ.141) MINT(61)=1
0279         CALL PYWIDT(KFA,P(ID,5)**2,WDTP,WDTE)
0280         WDTE0S=WDTE(0,1)+WDTE(0,IPM)+WDTE(0,4)
0281         IF(KFB.EQ.KFA) WDTE0S=WDTE0S+WDTE(0,5)
0282         IF(WDTE0S.LE.0D0) GOTO 330
0283         RKFL=WDTE0S*PYR(0)
0284         IDL=0
0285   200   IDL=IDL+1
0286         IDC=IDL+MDCY(KCA,2)-1
0287         RKFL=RKFL-(WDTE(IDL,1)+WDTE(IDL,IPM)+WDTE(IDL,4))
0288         IF(KFB.EQ.KFA) RKFL=RKFL-WDTE(IDL,5)
0289         IF(IDL.LT.MDCY(KCA,3).AND.RKFL.GT.0D0) GOTO 200
0290
0291 C...Read out flavours and colour charges of decay channel chosen.
0292         KCQM(JT)=KCHG(KCA,2)*ISIGN(1,K(ID,2))
0293         IF(KCQM(JT).EQ.-2) KCQM(JT)=2
0294         KFL1(JT)=KFDP(IDC,1)*ISIGN(1,K(ID,2))
0295         KFC1A=PYCOMP(IABS(KFL1(JT)))
0296         IF(KCHG(KFC1A,3).EQ.0) KFL1(JT)=IABS(KFL1(JT))
0297         KCQ1(JT)=KCHG(KFC1A,2)*ISIGN(1,KFL1(JT))
0298         IF(KCQ1(JT).EQ.-2) KCQ1(JT)=2
0299         KFL2(JT)=KFDP(IDC,2)*ISIGN(1,K(ID,2))
0300         KFC2A=PYCOMP(IABS(KFL2(JT)))
0301         IF(KCHG(KFC2A,3).EQ.0) KFL2(JT)=IABS(KFL2(JT))
0302         KCQ2(JT)=KCHG(KFC2A,2)*ISIGN(1,KFL2(JT))
0303         IF(KCQ2(JT).EQ.-2) KCQ2(JT)=2
0304         KFL3(JT)=KFDP(IDC,3)*ISIGN(1,K(ID,2))
0305         KCQ3(JT)=0
0306         IF(KFL3(JT).NE.0) THEN
0307           KFC3A=PYCOMP(IABS(KFL3(JT)))
0308           IF(KCHG(KFC3A,3).EQ.0) KFL3(JT)=IABS(KFL3(JT))
0309           KCQ3(JT)=KCHG(KFC3A,2)*ISIGN(1,KFL3(JT))
0310           IF(KCQ3(JT).EQ.-2) KCQ3(JT)=2
0311         ENDIF
0312
0313 C...Set/save further info on channel.
0314         KDCY(JT)=1
0315         IF(KFB.EQ.KFA) KEQL(JT)=MDME(IDC,1)
0316         NSD(JT)=N
0317         HGZ(JT,1)=VINT(111)
0318         HGZ(JT,2)=VINT(112)
0319         HGZ(JT,3)=VINT(114)
0320         JTZ=JT
0321
0322 C...Select masses; to begin with assume resonances narrow.
0323         DO 220 I=1,3
0324           P(N+I,5)=0D0
0325           PMMN(I)=0D0
0326           IF(I.EQ.1) THEN
0327             KFLW=IABS(KFL1(JT))
0328             KCW=KFC1A
0329           ELSEIF(I.EQ.2) THEN
0330             KFLW=IABS(KFL2(JT))
0331             KCW=KFC2A
0332           ELSEIF(I.EQ.3) THEN
0333             IF(KFL3(JT).EQ.0) GOTO 220
0334             KFLW=IABS(KFL3(JT))
0335             KCW=KFC3A
0336           ENDIF
0337           P(N+I,5)=PMAS(KCW,1)
0338 CMRENNA++
0339 C...This prevents SUSY/t particles from becoming too light.
0340           IF(KFLW/KSUSY1.EQ.1.OR.KFLW/KSUSY1.EQ.2) THEN
0341             PMMN(I)=PMAS(KCW,1)
0342             DO 210 IDC=MDCY(KCW,2),MDCY(KCW,2)+MDCY(KCW,3)-1
0343               IF(MDME(IDC,1).GT.0.AND.BRAT(IDC).GT.1E-4) THEN
0344                 PMSUM=PMAS(PYCOMP(KFDP(IDC,1)),1)+
0345      &          PMAS(PYCOMP(KFDP(IDC,2)),1)
0346                 IF(KFDP(IDC,3).NE.0) PMSUM=PMSUM+
0347      &          PMAS(PYCOMP(KFDP(IDC,3)),1)
0348                 PMMN(I)=MIN(PMMN(I),PMSUM)
0349               ENDIF
0350   210       CONTINUE
0351 CMRENNA--
0352           ELSEIF(KFLW.EQ.6) THEN
0353             PMMN(I)=PMAS(24,1)+PMAS(5,1)
0354           ENDIF
0355   220   CONTINUE
0356
0357 C...Check which two out of three are widest.
0358         IWID1=1
0359         IWID2=2
0360         PWID1=PMAS(KFC1A,2)
0361         PWID2=PMAS(KFC2A,2)
0362         KFLW1=IABS(KFL1(JT))
0363         KFLW2=IABS(KFL2(JT))
0364         IF(KFL3(JT).NE.0) THEN
0365           PWID3=PMAS(KFC3A,2)
0366           IF(PWID3.GT.PWID1.AND.PWID2.GE.PWID1) THEN
0367             IWID1=3
0368             PWID1=PWID3
0369             KFLW1=IABS(KFL3(JT))
0370           ELSEIF(PWID3.GT.PWID2) THEN
0371             IWID2=3
0372             PWID2=PWID3
0373             KFLW2=IABS(KFL3(JT))
0374           ENDIF
0375         ENDIF
0376
0377 C...If all narrow then only check that masses consistent.
0378         IF(MSTP(42).LE.0.OR.(PWID1.LT.PARP(41).AND.
0379      &  PWID2.LT.PARP(41))) THEN
0380 CMRENNA++
0381 C....Handle near degeneracy cases.
0382           IF(KFA/KSUSY1.EQ.1.OR.KFA/KSUSY1.EQ.2) THEN
0383             IF(P(N+1,5)+P(N+2,5)+P(N+3,5).GT.P(ID,5)) THEN
0384               P(N+1,5)=P(ID,5)-P(N+2,5)-0.5D0
0385               IF(P(N+1,5).LT.0D0) P(N+1,5)=0D0
0386             ENDIF
0387           ENDIF
0388 CMRENNA--
0389           IF(P(N+1,5)+P(N+2,5)+P(N+3,5).GT.P(ID,5)) THEN
0390             CALL PYERRM(13,'(PYRESD:) daughter masses too large')
0391             MINT(51)=1
0392             GOTO 720
0393           ELSEIF(P(N+1,5)+P(N+2,5)+P(N+3,5)+PARJ(64).GT.P(ID,5)) THEN
0394             CALL PYERRM(3,'(PYRESD:) daughter masses too large')
0395             MINT(51)=1
0396             GOTO 720
0397           ENDIF
0398
0399 C...For three wide resonances select narrower of three
0400 C...according to BW decoupled from rest.
0401         ELSE
0402           PMTOT=P(ID,5)
0403           IF(KFL3(JT).NE.0) THEN
0404             IWID3=6-IWID1-IWID2
0405             KFLW3=IABS(KFL1(JT))+IABS(KFL2(JT))+IABS(KFL3(JT))-
0406      &      KFLW1-KFLW2
0407             LOOP=0
0408   230       LOOP=LOOP+1
0409             P(N+IWID3,5)=PYMASS(KFLW3)
0410             IF(LOOP.LE.10.AND. P(N+IWID3,5).LE.PMMN(IWID3)) GOTO 230
0411             PMTOT=PMTOT-P(N+IWID3,5)
0412           ENDIF
0413 C...Select other two correlated within remaining phase space.
0414           IF(IP.EQ.1) THEN
0415             CKIN45=CKIN(45)
0416             CKIN47=CKIN(47)
0417             CKIN(45)=MAX(PMMN(IWID1),CKIN(45))
0418             CKIN(47)=MAX(PMMN(IWID2),CKIN(47))
0419             CALL PYOFSH(2,KFA,KFLW1,KFLW2,PMTOT,P(N+IWID1,5),
0420      &      P(N+IWID2,5))
0421             CKIN(45)=CKIN45
0422             CKIN(47)=CKIN47
0423           ELSE
0424             CKIN(49)=PMMN(IWID1)
0425             CKIN(50)=PMMN(IWID2)
0426             CALL PYOFSH(5,KFA,KFLW1,KFLW2,PMTOT,P(N+IWID1,5),
0427      &      P(N+IWID2,5))
0428             CKIN(49)=0D0
0429             CKIN(50)=0D0
0430           ENDIF
0431           IF(MINT(51).EQ.1) GOTO 720
0432         ENDIF
0433
0434 C...Begin fill decay products, with colour flow for coloured objects.
0435         MSTU10=MSTU(10)
0436         MSTU(10)=1
0437         MSTU(19)=1
0438
0439 CMRENNA++
0440 C...1) Three-body decays of SUSY particles (plus special case top).
0441         IF(KFL3(JT).NE.0) THEN
0442           DO 250 I=N+1,N+3
0443             DO 240 J=1,5
0444               K(I,J)=0
0445               V(I,J)=0D0
0446   240       CONTINUE
0447             MCT(I,1)=0
0448             MCT(I,2)=0
0449   250     CONTINUE
0450           K(N+1,1)=1
0451           K(N+1,2)=KFL1(JT)
0452           K(N+2,1)=1
0453           K(N+2,2)=KFL2(JT)
0454           K(N+3,1)=1
0455           K(N+3,2)=KFL3(JT)
0456           IDIN=ID
0457           CALL PYTBDY(IDIN)
0458
0459 C...Set colour flow for t -> W + b + Z.
0460           IF(KFA.EQ.6) THEN
0461             K(N+2,1)=3
0462             ISID=4
0463             IF(KCQM(JT).EQ.-1) ISID=5
0464             IDAU=N+2
0465             K(ID,ISID)=K(ID,ISID)+IDAU
0466             K(IDAU,ISID)=MSTU(5)*ID
0467
0468 C...Set colour flow in three-body decays - programmed as special cases.
0469
0470           ELSEIF(KFC2A.LE.6) THEN
0471             K(N+2,1)=3
0472             K(N+3,1)=3
0473             ISID=4
0474             IF(KFL2(JT).LT.0) ISID=5
0475             K(N+2,ISID)=MSTU(5)*(N+3)
0476             K(N+3,9-ISID)=MSTU(5)*(N+2)
0477 C...PS++: Bugfix 16 MAR 2006 for 3-body squark decays (e.g. via SLHA)
0478           ELSEIF(KFA.GT.KSUSY1.AND.MOD(KFA,KSUSY1).LT.10
0479      &          .AND.KFL3(JT).NE.0) THEN
0480             KQSUMA=IABS(KCQ1(JT))+IABS(KCQ2(JT))+IABS(KCQ3(JT))
0481 C...3-body decays of squarks to colour singlets plus one quark
0482             IF (KQSUMA.EQ.1) THEN
0483 C...Find quark
0484               IQ=0
0485               IF (KCQ1(JT).NE.0) IQ=1
0486               IF (KCQ2(JT).NE.0) IQ=2
0487               IF (KCQ3(JT).NE.0) IQ=3
0488               ISID=4
0489               IF (K(N+IQ,2).LT.0) ISID=5
0490               K(N+IQ,1)=3
0491               K(ID,ISID)=K(ID,ISID)+(N+IQ)
0492               K(N+IQ,ISID)=MSTU(5)*ID
0493             ENDIF
0494 C...PS--
0495           ENDIF
0496           IF(KFL1(JT).EQ.KSUSY1+21) THEN
0497             K(N+1,1)=3
0498             K(N+2,1)=3
0499             K(N+3,1)=3
0500             ISID=4
0501             IF(KFL2(JT).LT.0) ISID=5
0502             K(N+1,ISID)=MSTU(5)*(N+2)
0503             K(N+1,9-ISID)=MSTU(5)*(N+3)
0504             K(N+2,ISID)=MSTU(5)*(N+1)
0505             K(N+3,9-ISID)=MSTU(5)*(N+1)
0506           ENDIF
0507           IF(KFA.EQ.KSUSY1+21) THEN
0508             K(N+2,1)=3
0509             K(N+3,1)=3
0510             ISID=4
0511             IF(KFL2(JT).LT.0) ISID=5
0512             K(ID,ISID)=K(ID,ISID)+(N+2)
0513             K(ID,9-ISID)=K(ID,9-ISID)+(N+3)
0514             K(N+2,ISID)=MSTU(5)*ID
0515             K(N+3,9-ISID)=MSTU(5)*ID
0516           ENDIF
0517           NSAV=N
0518           N=N+3
0519           N=NSAV
0520 CMRENNA--
0521
0522           IF(KFA.GE.KSUSY1+22.AND.KFA.LE.KSUSY1+37.AND.
0523      &    IABS(KCQ2(JT)).EQ.1) THEN
0524             K(N+2,1)=3
0525             K(N+3,1)=3
0526             ISID=4
0527             IF(KFL2(JT).LT.0) ISID=5
0528             K(N+2,ISID)=MSTU(5)*(N+3)
0529             K(N+3,9-ISID)=MSTU(5)*(N+2)
0530           ENDIF
0531
0532 C...Set colour flow in three-body decays with baryon number violation.
0533 C...Neutralino and chargino decays first.
0534           KCQSUM=KCQ1(JT)+KCQ2(JT)+KCQ3(JT)
0535           IF(KCQM(JT).EQ.0.AND.IABS(KCQSUM).EQ.3) THEN
0536             ITJUNC(JT)=(1+(1-KCQ1(JT))/2)
0537             K(N+4,4)=ITJUNC(JT)*MSTU(5)
0538 C...Insert junction to keep track of colours.
0539             IF(KCQ1(JT).NE.0) K(N+1,1)=3
0540             IF(KCQ2(JT).NE.0) K(N+2,1)=3
0541             IF(KCQ3(JT).NE.0) K(N+3,1)=3
0542 C...Set special junction codes:
0543             K(N+4,1)=42
0544             K(N+4,2)=88
0545
0546 C...Order decay products by invariant mass. (will be used in PYSTRF).
0547             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)-
0548      &      P(N+1,3)*P(N+2,3)
0549             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)-
0550      &      P(N+1,3)*P(N+3,3)
0551             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)-
0552      &      P(N+2,3)*P(N+3,3)
0553             IF(PM12.LT.PM13.AND.PM12.LT.PM23) THEN
0554               K(N+4,4)=N+3+K(N+4,4)
0555               K(N+4,5)=N+1+MSTU(5)*(N+2)
0556             ELSEIF(PM13.LT.PM23) THEN
0557               K(N+4,4)=N+2+K(N+4,4)
0558               K(N+4,5)=N+1+MSTU(5)*(N+3)
0559             ELSE
0560               K(N+4,4)=N+1+K(N+4,4)
0561               K(N+4,5)=N+2+MSTU(5)*(N+3)
0562             ENDIF
0563             DO 260 J=1,5
0564               P(N+4,J)=0D0
0565               V(N+4,J)=0D0
0566   260       CONTINUE
0567 C...Connect daughters to junction.
0568             DO 270 II=N+1,N+3
0569               K(II,4)=0
0570               K(II,5)=0
0571               K(II,ITJUNC(JT)+3)=MSTU(5)*(N+4)
0572   270       CONTINUE
0573 C...Particle counter should be stepped up one extra for junction.
0574             N=N+1
0575
0576 C...Gluino decays.
0577           ELSEIF (KCQM(JT).EQ.2.AND.IABS(KCQSUM).EQ.3) THEN
0578             ITJUNC(JT)=(5+(1-KCQ1(JT))/2)
0579             K(N+4,4)=ITJUNC(JT)*MSTU(5)
0580 C...Insert junction to keep track of colours.
0581             IF(KCQ1(JT).NE.0) K(N+1,1)=3
0582             IF(KCQ2(JT).NE.0) K(N+2,1)=3
0583             IF(KCQ3(JT).NE.0) K(N+3,1)=3
0584             K(N+4,1)=42
0585             K(N+4,2)=88
0586             DO 280 J=1,5
0587               P(N+4,J)=0D0
0588               V(N+4,J)=0D0
0589   280       CONTINUE
0590             CTMSUM=0D0
0591             DO 290 II=N+1,N+3
0592               K(II,4)=0
0593               K(II,5)=0
0594 C...Start by connecting all daughters to junction.
0595               K(II,ITJUNC(JT)-1)=MSTU(5)*(N+4)
0596 C...Only consider colour topologies with off shell resonances.
0597               RMQ1=PMAS(PYCOMP(K(II,2)),1)
0598               RMRES=PMAS(PYCOMP(KSUSY1+IABS(K(II,2))),1)
0599               RMGLU=PMAS(PYCOMP(KSUSY1+21),1)
0600               IF (RMGLU-RMQ1.LT.RMRES) THEN
0601 C...Calculate propagators for each colour topology.
0602                 RM2Q23=RMGLU**2+RMQ1**2-2D0*(P(II,4)*P(ID,4)+P(II,1)
0603      &               *P(ID,1)+P(II,2)*P(ID,2)+P(II,3)*P(ID,3))
0604                 CTM2(II-N)=1D0/(RM2Q23-RMRES**2)**2
0605               ELSE
0606                 CTM2(II-N)=0D0
0607               ENDIF
0608               CTMSUM=CTMSUM+CTM2(II-N)
0609   290       CONTINUE
0610             CTMSUM=PYR(0)*CTMSUM
0611 C...Select colour topology J, with most off shell least likely.
0612             J=0
0613   300       J=J+1
0614             CTMSUM=CTMSUM-CTM2(J)
0615             IF (CTMSUM.GT.0D0) GOTO 300
0616 C...The lucky winner gets its colour (anti-colour) directly from gluino.
0617             K(N+J,ITJUNC(JT)-1)=MSTU(5)*ID
0618             K(ID,ITJUNC(JT)-1)=N+J+(K(ID,ITJUNC(JT)-1)/MSTU(5))*MSTU(5)
0619 C...The other gluino colour is connected to junction
0620             K(ID,10-ITJUNC(JT))=N+4+(K(ID,10-ITJUNC(JT))/MSTU(5))*
0621      &      MSTU(5)
0622             K(N+4,4)=K(N+4,4)+ID
0623 C...Lastly, connect junction to remaining daughters.
0624             K(N+4,5)=N+1+MOD(J,3)+MSTU(5)*(N+1+MOD(J+1,3))
0625 C...Particle counter should be stepped up one extra for junction.
0626             N=N+1
0627          ENDIF
0628
0629 C...Update particle counter.
0630           N=N+3
0631
0632 C...2) Everything else two-body decay.
0633         ELSE
0634           CALL PY2ENT(N+1,KFL1(JT),KFL2(JT),P(ID,5))
0635           MCT(N-1,1)=0
0636           MCT(N-1,2)=0
0637           MCT(N,1)=0
0638           MCT(N,2)=0
0639 C...First set colour flow as if mother colour singlet.
0640           IF(KCQ1(JT).NE.0) THEN
0641             K(N-1,1)=3
0642             IF(KCQ1(JT).NE.-1) K(N-1,4)=MSTU(5)*N
0643             IF(KCQ1(JT).NE.1) K(N-1,5)=MSTU(5)*N
0644           ENDIF
0645           IF(KCQ2(JT).NE.0) THEN
0646             K(N,1)=3
0647             IF(KCQ2(JT).NE.-1) K(N,4)=MSTU(5)*(N-1)
0648             IF(KCQ2(JT).NE.1) K(N,5)=MSTU(5)*(N-1)
0649           ENDIF
0650 C...Then redirect colour flow if mother (anti)triplet.
0651           IF(KCQM(JT).EQ.0) THEN
0652           ELSEIF(KCQM(JT).NE.2) THEN
0653             ISID=4
0654             IF(KCQM(JT).EQ.-1) ISID=5
0655             IDAU=N-1
0656             IF(KCQ1(JT).EQ.0.OR.KCQ2(JT).EQ.2) IDAU=N
0657             K(ID,ISID)=K(ID,ISID)+IDAU
0658             K(IDAU,ISID)=MSTU(5)*ID
0659 C...Then redirect colour flow if mother octet.
0660           ELSEIF(KCQ1(JT).EQ.0.OR.KCQ2(JT).EQ.0) THEN
0661             IDAU=N-1
0662             IF(KCQ1(JT).EQ.0) IDAU=N
0663             K(ID,4)=K(ID,4)+IDAU
0664             K(ID,5)=K(ID,5)+IDAU
0665             K(IDAU,4)=MSTU(5)*ID
0666             K(IDAU,5)=MSTU(5)*ID
0667           ELSE
0668             ISID=4
0669             IF(KCQ1(JT).EQ.-1) ISID=5
0670             IF(KCQ1(JT).EQ.2) ISID=INT(4.5D0+PYR(0))
0671             K(ID,ISID)=K(ID,ISID)+(N-1)
0672             K(ID,9-ISID)=K(ID,9-ISID)+N
0673             K(N-1,ISID)=MSTU(5)*ID
0674             K(N,9-ISID)=MSTU(5)*ID
0675           ENDIF
0676
0677 C...Insert junction
0678           IF(IABS(KCQ1(JT)+KCQ2(JT)-KCQM(JT)).EQ.3) THEN
0679             N=N+1
0680 C...~q* mother: type 3 junction. ~q mother: type 4.
0681             ITJUNC(JT)=(7+KCQM(JT))/2
0682 C...Specify junction KF and set colour flow from junction
0683             K(N,1)=42
0684             K(N,2)=88
0685             K(N,3)=ID
0686 C...Junction type encoded together with mother:
0687             K(N,4)=ID+ITJUNC(JT)*MSTU(5)
0688             K(N,5)=N-1+MSTU(5)*(N-2)
0689 C...Zero P and V for junction (V filled later)
0690             DO 310 J=1,5
0691               P(N,J)=0D0
0692               V(N,J)=0D0
0693   310       CONTINUE
0694 C...Set colour flow from mother to junction
0695             K(ID,8-ITJUNC(JT))= N + MSTU(5)*(K(ID,8-ITJUNC(JT))/MSTU(5))
0696 C...Set colour flow from daughters to junction
0697             DO 320 II=N-2,N-1
0698               K(II,4) = 0
0699               K(II,5) = 0
0700 C...(Anti-)colour mother is junction.
0701               K(II,1+ITJUNC(JT)) = MSTU(5)*(N)
0702   320       CONTINUE
0703           ENDIF
0704         ENDIF
0705
0706 C...End loop over resonances for daughter flavour and mass selection.
0707         MSTU(10)=MSTU10
0708   330   IF(MWID(KCA).NE.0.AND.(KFL1(JT).EQ.0.OR.KFL3(JT).NE.0))
0709      &  NINH=NINH+1
0710         IF(IRES.GT.0.AND.MWID(KCA).NE.0.AND.MDCY(KCA,1).NE.0.AND.
0711      &  KFL1(JT).EQ.0) THEN
0712           WRITE(CODE,'(I9)') K(ID,2)
0713           WRITE(MASS,'(F9.3)') P(ID,5)
0714           CALL PYERRM(3,'(PYRESD:) Failed to decay particle'//
0715      &    CODE//' with mass'//MASS)
0716           MINT(51)=1
0717           GOTO 720
0718         ENDIF
0719   340 CONTINUE
0720
0721 C...Check for allowed combinations. Skip if no decays.
0722       IF(JTMAX.EQ.1) THEN
0723         IF(KDCY(1).EQ.0) GOTO 710
0724       ELSEIF(JTMAX.EQ.2) THEN
0725         IF(KDCY(1).EQ.0.AND.KDCY(2).EQ.0) GOTO 710
0726         IF(KEQL(1).EQ.4.AND.KEQL(2).EQ.4) GOTO 180
0727         IF(KEQL(1).EQ.5.AND.KEQL(2).EQ.5) GOTO 180
0728       ELSEIF(JTMAX.EQ.3) THEN
0729         IF(KDCY(1).EQ.0.AND.KDCY(2).EQ.0.AND.KDCY(3).EQ.0) GOTO 710
0730         IF(KEQL(1).EQ.4.AND.KEQL(2).EQ.4) GOTO 180
0731         IF(KEQL(1).EQ.4.AND.KEQL(3).EQ.4) GOTO 180
0732         IF(KEQL(2).EQ.4.AND.KEQL(3).EQ.4) GOTO 180
0733         IF(KEQL(1).EQ.5.AND.KEQL(2).EQ.5) GOTO 180
0734         IF(KEQL(1).EQ.5.AND.KEQL(3).EQ.5) GOTO 180
0735         IF(KEQL(2).EQ.5.AND.KEQL(3).EQ.5) GOTO 180
0736       ENDIF
0737
0738 C...Special case: matrix element option for Z0 decay to quarks.
0739       IF(MSTP(48).EQ.1.AND.ISUB.EQ.1.AND.JTMAX.EQ.1.AND.
0740      &IABS(MINT(11)).EQ.11.AND.IABS(KFL1(1)).LE.5) THEN
0741
0742 C...Check consistency of MSTJ options set.
0743         IF(MSTJ(109).EQ.2.AND.MSTJ(110).NE.1) THEN
0744           CALL PYERRM(6,
0745      &    '(PYRESD:) MSTJ(109) value requires MSTJ(110) = 1')
0746           MSTJ(110)=1
0747         ENDIF
0748         IF(MSTJ(109).EQ.2.AND.MSTJ(111).NE.0) THEN
0749           CALL PYERRM(6,
0750      &    '(PYRESD:) MSTJ(109) value requires MSTJ(111) = 0')
0751
0752           MSTJ(111)=0
0753         ENDIF
0754
0755 C...Select alpha_strong behaviour.
0756         MST111=MSTU(111)
0757         PAR112=PARU(112)
0758         MSTU(111)=MSTJ(108)
0759         IF(MSTJ(108).EQ.2.AND.(MSTJ(101).EQ.0.OR.MSTJ(101).EQ.1))
0760      &  MSTU(111)=1
0761         PARU(112)=PARJ(121)
0762         IF(MSTU(111).EQ.2) PARU(112)=PARJ(122)
0763
0764 C...Find axial fraction in total cross section for scalar gluon model.
0765         PARJ(171)=0D0
0766         IF((IABS(MSTJ(101)).EQ.1.AND.MSTJ(109).EQ.1).OR.
0767      &  (MSTJ(101).EQ.5.AND.MSTJ(49).EQ.1)) THEN
0768           POLL=1D0-PARJ(131)*PARJ(132)
0769           SFF=1D0/(16D0*XW*XW1)
0770           SFW=P(ID,5)**4/((P(ID,5)**2-PARJ(123)**2)**2+
0771      &    (PARJ(123)*PARJ(124))**2)
0772           SFI=SFW*(1D0-(PARJ(123)/P(ID,5))**2)
0773           VE=4D0*XW-1D0
0774           HF1I=SFI*SFF*(VE*POLL+PARJ(132)-PARJ(131))
0775           HF1W=SFW*SFF**2*((VE**2+1D0)*POLL+2D0*VE*
0776      &    (PARJ(132)-PARJ(131)))
0777           KFLC=IABS(KFL1(1))
0778           PMQ=PYMASS(KFLC)
0779           QF=KCHG(KFLC,1)/3D0
0780           VQ=1D0
0781           IF(MOD(MSTJ(103),2).EQ.1) VQ=SQRT(MAX(0D0,
0782      &    1D0-(2D0*PMQ/P(ID,5))**2))
0783           VF=SIGN(1D0,QF)-4D0*QF*XW
0784           RFV=0.5D0*VQ*(3D0-VQ**2)*(QF**2*POLL-2D0*QF*VF*HF1I+
0785      &    VF**2*HF1W)+VQ**3*HF1W
0786           IF(RFV.GT.0D0) PARJ(171)=MIN(1D0,VQ**3*HF1W/RFV)
0787         ENDIF
0788
0789 C...Choice of jet configuration.
0790         CALL PYXJET(P(ID,5),NJET,CUT)
0791         KFLC=IABS(KFL1(1))
0792         KFLN=21
0793         IF(NJET.EQ.4) THEN
0794           CALL PYX4JT(NJET,CUT,KFLC,P(ID,5),KFLN,X1,X2,X4,X12,X14)
0795         ELSEIF(NJET.EQ.3) THEN
0796           CALL PYX3JT(NJET,CUT,KFLC,P(ID,5),X1,X3)
0797         ELSE
0798           MSTJ(120)=1
0799         ENDIF
0800
0801 C...Fill jet configuration; return if incorrect kinematics.
0802         NC=N-2
0803         IF(NJET.EQ.2.AND.MSTJ(101).NE.5) THEN
0804           CALL PY2ENT(NC+1,KFLC,-KFLC,P(ID,5))
0805         ELSEIF(NJET.EQ.2) THEN
0806           CALL PY2ENT(-(NC+1),KFLC,-KFLC,P(ID,5))
0807         ELSEIF(NJET.EQ.3) THEN
0808           CALL PY3ENT(NC+1,KFLC,21,-KFLC,P(ID,5),X1,X3)
0809         ELSEIF(KFLN.EQ.21) THEN
0810           CALL PY4ENT(NC+1,KFLC,KFLN,KFLN,-KFLC,P(ID,5),X1,X2,X4,
0811      &    X12,X14)
0812         ELSE
0813           CALL PY4ENT(NC+1,KFLC,-KFLN,KFLN,-KFLC,P(ID,5),X1,X2,X4,
0814      &    X12,X14)
0815         ENDIF
0816         IF(MSTU(24).NE.0) THEN
0817           MINT(51)=1
0818           MSTU(111)=MST111
0819           PARU(112)=PAR112
0820           GOTO 720
0821         ENDIF
0822
0823 C...Angular orientation according to matrix element.
0824         IF(MSTJ(106).EQ.1) THEN
0825           CALL PYXDIF(NC,NJET,KFLC,P(ID,5),CHIZ,THEZ,PHIZ)
0826           IF(MINT(11).LT.0) THEZ=PARU(1)-THEZ
0827           CTHE(1)=COS(THEZ)
0828           CALL PYROBO(NC+1,N,0D0,CHIZ,0D0,0D0,0D0)
0829           CALL PYROBO(NC+1,N,THEZ,PHIZ,0D0,0D0,0D0)
0830         ENDIF
0831
0832 C...Boost partons to Z0 rest frame.
0833         CALL PYROBO(NC+1,N,0D0,0D0,P(ID,1)/P(ID,4),
0834      &  P(ID,2)/P(ID,4),P(ID,3)/P(ID,4))
0835
0836 C...Mark decayed resonance and add documentation lines,
0837         K(ID,1)=K(ID,1)+10
0838         IDOC=MINT(83)+MINT(4)
0839         DO 360 I=NC+1,N
0840           I1=MINT(83)+MINT(4)+1
0841           K(I,3)=I1
0842           IF(MSTP(128).GE.1) K(I,3)=ID
0843           IF(MSTP(128).LE.1.AND.MINT(4).LT.MSTP(126)) THEN
0844             MINT(4)=MINT(4)+1
0845             K(I1,1)=21
0846             K(I1,2)=K(I,2)
0847             K(I1,3)=IREF(IP,4)
0848             DO 350 J=1,5
0849               P(I1,J)=P(I,J)
0850   350       CONTINUE
0851           ENDIF
0852   360   CONTINUE
0853
0854 C...Generate parton shower.
0855         IF(MSTJ(101).EQ.5.AND.MINT(35).LE.1) THEN
0856           CALL PYSHOW(N-1,N,P(ID,5))
0857         ELSEIF(MSTJ(101).EQ.5.AND.MINT(35).GE.2) THEN
0858           NPART=2
0859           IPART(1)=N-1
0860           IPART(2)=N
0861           PTPART(1)=0.5D0*P(ID,5)
0862           PTPART(2)=PTPART(1)
0863           NCT=NCT+1
0864           IF(K(N-1,2).GT.0) THEN
0865             MCT(N-1,1)=NCT
0866             MCT(N,2)=NCT
0867           ELSE
0868             MCT(N-1,2)=NCT
0869             MCT(N,1)=NCT
0870           ENDIF
0871           CALL PYPTFS(2,0.5D0*P(ID,5),0D0,PTGEN)
0872         ENDIF
0873
0874 C... End special case for Z0: skip ahead.
0875         MSTU(111)=MST111
0876         PARU(112)=PAR112
0877         GOTO 700
0878       ENDIF
0879
0880 C...Order incoming partons and outgoing resonances.
0881       IF(JTMAX.EQ.2.AND.ISUB.NE.0.AND.MSTP(47).GE.1.AND.
0882      &NINH.EQ.0) THEN
0883         ILIN(1)=MINT(84)+1
0884         IF(K(MINT(84)+1,2).GT.0) ILIN(1)=MINT(84)+2
0885         IF(K(ILIN(1),2).EQ.21.OR.K(ILIN(1),2).EQ.22)
0886      &  ILIN(1)=2*MINT(84)+3-ILIN(1)
0887         ILIN(2)=2*MINT(84)+3-ILIN(1)
0888         IMIN=1
0889         IF(IREF(IP,7).EQ.25.OR.IREF(IP,7).EQ.35.OR.IREF(IP,7)
0890      &  .EQ.36) IMIN=3
0891         IMAX=2
0892         IORD=1
0893         IF(K(IREF(IP,1),2).EQ.23) IORD=2
0894         IF(K(IREF(IP,1),2).EQ.24.AND.K(IREF(IP,2),2).EQ.-24) IORD=2
0895         IAKIPD=IABS(K(IREF(IP,IORD),2))
0896         IF(IAKIPD.EQ.25.OR.IAKIPD.EQ.35.OR.IAKIPD.EQ.36) IORD=3-IORD
0897         IF(KDCY(IORD).EQ.0) IORD=3-IORD
0898
0899 C...Order decay products of resonances.
0900         DO 370 JT=IORD,3-IORD,3-2*IORD
0901           IF(KDCY(JT).EQ.0) THEN
0902             ILIN(IMAX+1)=NSD(JT)
0903             IMAX=IMAX+1
0904           ELSEIF(K(NSD(JT)+1,2).GT.0) THEN
0905             ILIN(IMAX+1)=N+2*JT-1
0906             ILIN(IMAX+2)=N+2*JT
0907             IMAX=IMAX+2
0908             K(N+2*JT-1,2)=K(NSD(JT)+1,2)
0909             K(N+2*JT,2)=K(NSD(JT)+2,2)
0910           ELSE
0911             ILIN(IMAX+1)=N+2*JT
0912
0913             ILIN(IMAX+2)=N+2*JT-1
0914             IMAX=IMAX+2
0915             K(N+2*JT-1,2)=K(NSD(JT)+1,2)
0916             K(N+2*JT,2)=K(NSD(JT)+2,2)
0917           ENDIF
0918   370   CONTINUE
0919
0920 C...Find charge, isospin, left- and righthanded couplings.
0921         DO 390 I=IMIN,IMAX
0922           DO 380 J=1,4
0923             COUP(I,J)=0D0
0924   380     CONTINUE
0925           KFA=IABS(K(ILIN(I),2))
0926           IF(KFA.EQ.0.OR.KFA.GT.20) GOTO 390
0927           COUP(I,1)=KCHG(KFA,1)/3D0
0928           COUP(I,2)=(-1)**MOD(KFA,2)
0929           COUP(I,4)=-2D0*COUP(I,1)*XWV
0930           COUP(I,3)=COUP(I,2)+COUP(I,4)
0931   390   CONTINUE
0932
0933 C...Full propagator dependence and flavour correlations for 2 gamma*/Z.
0934         IF(ISUB.EQ.22) THEN
0935           DO 420 I=3,5,2
0936             I1=IORD
0937             IF(I.EQ.5) I1=3-IORD
0938             DO 410 J1=1,2
0939               DO 400 J2=1,2
0940                 CORL(I/2,J1,J2)=COUP(1,1)**2*HGZ(I1,1)*COUP(I,1)**2/
0941      &          16D0+COUP(1,1)*COUP(1,J1+2)*HGZ(I1,2)*COUP(I,1)*
0942      &          COUP(I,J2+2)/4D0+COUP(1,J1+2)**2*HGZ(I1,3)*
0943      &          COUP(I,J2+2)**2
0944   400         CONTINUE
0945   410       CONTINUE
0946   420     CONTINUE
0947           COWT12=(CORL(1,1,1)+CORL(1,1,2))*(CORL(2,1,1)+CORL(2,1,2))+
0948      &    (CORL(1,2,1)+CORL(1,2,2))*(CORL(2,2,1)+CORL(2,2,2))
0949           COMX12=(CORL(1,1,1)+CORL(1,1,2)+CORL(1,2,1)+CORL(1,2,2))*
0950      &    (CORL(2,1,1)+CORL(2,1,2)+CORL(2,2,1)+CORL(2,2,2))
0951
0952           IF(COWT12.LT.PYR(0)*COMX12) GOTO 180
0953         ENDIF
0954       ENDIF
0955
0956 C...Select angular orientation type - Z'/W' only.
0957       MZPWP=0
0958       IF(ISUB.EQ.141) THEN
0959         IF(PYR(0).LT.PARU(130)) MZPWP=1
0960         IF(IP.EQ.2) THEN
0961           IF(IABS(K(IREF(2,1),2)).EQ.37) MZPWP=2
0962           IAKIR=IABS(K(IREF(2,2),2))
0963           IF(IAKIR.EQ.25.OR.IAKIR.EQ.35.OR.IAKIR.EQ.36) MZPWP=2
0964           IF(IAKIR.LE.20) MZPWP=2
0965         ENDIF
0966         IF(IP.GE.3) MZPWP=2
0967       ELSEIF(ISUB.EQ.142) THEN
0968         IF(PYR(0).LT.PARU(136)) MZPWP=1
0969         IF(IP.EQ.2) THEN
0970           IAKIR=IABS(K(IREF(2,2),2))
0971           IF(IAKIR.EQ.25.OR.IAKIR.EQ.35.OR.IAKIR.EQ.36) MZPWP=2
0972           IF(IAKIR.LE.20) MZPWP=2
0973         ENDIF
0974         IF(IP.GE.3) MZPWP=2
0975       ENDIF
0976
0977 C...Select random angles (begin of weighting procedure).
0978   430 DO 440 JT=1,JTMAX
0979         IF(KDCY(JT).EQ.0) GOTO 440
0980         IF(JTMAX.EQ.1.AND.ISUB.NE.0.AND.IHDEC.EQ.0) THEN
0981           CTHE(JT)=VINT(13)+(VINT(33)-VINT(13)+VINT(34)-VINT(14))*PYR(0)
0982           IF(CTHE(JT).GT.VINT(33)) CTHE(JT)=CTHE(JT)+VINT(14)-VINT(33)
0983           PHI(JT)=VINT(24)
0984         ELSE
0985           CTHE(JT)=2D0*PYR(0)-1D0
0986           PHI(JT)=PARU(2)*PYR(0)
0987         ENDIF
0988   440 CONTINUE
0989
0990       IF(JTMAX.EQ.2.AND.MSTP(47).GE.1.AND.NINH.EQ.0) THEN
0991 C...Construct massless four-vectors.
0992         DO 460 I=N+1,N+4
0993           K(I,1)=1
0994           DO 450 J=1,5
0995             P(I,J)=0D0
0996             V(I,J)=0D0
0997   450     CONTINUE
0998   460   CONTINUE
0999         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