hijing/src/pyhiresd.f

0001
0002 C*********************************************************************
0003
0004       SUBROUTINE PYHIRESD
0005
0006 C...Allows resonances to decay (including parton showers for hadronic
0007 C...channels).
0008       IMPLICIT DOUBLE PRECISION(D)
0009       COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
0010       SAVE /LUJETS/
0011       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0012       SAVE /LUDAT1/
0013       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
0014       SAVE /LUDAT2/
0015       COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
0016       SAVE /LUDAT3/
0017       COMMON/PYHISUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
0018       SAVE /PYHISUBS/
0019       COMMON/PYHIPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0020       SAVE /PYHIPARS/
0021       COMMON/PYHIINT1/MINT(400),VINT(400)
0022       SAVE /PYHIINT1/
0023       COMMON/PYHIINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
0024       SAVE /PYHIINT2/
0025       COMMON/PYHIINT4/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3)
0026       SAVE /PYHIINT4/
0027       DIMENSION IREF(10,6),KDCY(2),KFL1(2),KFL2(2),NSD(2),ILIN(6),
0028      &COUP(6,4),PK(6,4),PKK(6,6),CTHE(2),PHI(2),WDTP(0:40),
0029      &WDTE(0:40,0:5)
0030       COMPLEX FGK,HA(6,6),HC(6,6)
0031
0032 C...The F, Xi and Xj functions of Gunion and Kunszt
0033 C...(Phys. Rev. D33, 665, plus errata from the authors).
0034       FGK(I1,I2,I3,I4,I5,I6)=4.*HA(I1,I3)*HC(I2,I6)*(HA(I1,I5)*
0035      &HC(I1,I4)+HA(I3,I5)*HC(I3,I4))
0036       DIGK(DT,DU)=-4.*D34*D56+DT*(3.*DT+4.*DU)+DT**2*(DT*DU/(D34*D56)-
0037      &2.*(1./D34+1./D56)*(DT+DU)+2.*(D34/D56+D56/D34))
0038       DJGK(DT,DU)=8.*(D34+D56)**2-8.*(D34+D56)*(DT+DU)-6.*DT*DU-
0039      &2.*DT*DU*(DT*DU/(D34*D56)-2.*(1./D34+1./D56)*(DT+DU)+
0040      &2.*(D34/D56+D56/D34))
0041
0042 C...Define initial two objects, initialize loop.
0043       ISUB=MINT(1)
0044       SH=VINT(44)
0045       IREF(1,5)=0
0046       IREF(1,6)=0
0047       IF(ISET(ISUB).EQ.1.OR.ISET(ISUB).EQ.3) THEN
0048         IREF(1,1)=MINT(84)+2+ISET(ISUB)
0049         IREF(1,2)=0
0050         IREF(1,3)=MINT(83)+6+ISET(ISUB)
0051         IREF(1,4)=0
0052       ELSEIF(ISET(ISUB).EQ.2.OR.ISET(ISUB).EQ.4) THEN
0053         IREF(1,1)=MINT(84)+1+ISET(ISUB)
0054         IREF(1,2)=MINT(84)+2+ISET(ISUB)
0055         IREF(1,3)=MINT(83)+5+ISET(ISUB)
0056         IREF(1,4)=MINT(83)+6+ISET(ISUB)
0057       ENDIF
0058       NP=1
0059       IP=0
0060   100 IP=IP+1
0061       NINH=0
0062
0063 C...Loop over one/two resonances; reset decay rates.
0064       JTMAX=2
0065       IF(IP.EQ.1.AND.(ISET(ISUB).EQ.1.OR.ISET(ISUB).EQ.3)) JTMAX=1
0066       DO 140 JT=1,JTMAX
0067       KDCY(JT)=0
0068       KFL1(JT)=0
0069       KFL2(JT)=0
0070       NSD(JT)=IREF(IP,JT)
0071       ID=IREF(IP,JT)
0072       IF(ID.EQ.0) GOTO 140
0073       KFA=IABS(K(ID,2))
0074       IF(KFA.LT.23.OR.KFA.GT.40) GOTO 140
0075       IF(MDCY(KFA,1).NE.0) THEN
0076         IF(ISUB.EQ.1.OR.ISUB.EQ.141) MINT(61)=1
0077         CALL PYHIWIDT(KFA,P(ID,5),WDTP,WDTE)
0078         IF(KCHG(KFA,3).EQ.0) THEN
0079           IPM=2
0080         ELSE
0081           IPM=(5+ISIGN(1,K(ID,2)))/2
0082         ENDIF
0083         IF(JTMAX.EQ.1.OR.IABS(K(IREF(IP,1),2)).NE.IABS(K(IREF(IP,2),2)))
0084      &  THEN
0085           I12=4
0086         ELSE
0087           IF(JT.EQ.1) I12=INT(4.5+RLU(0))
0088           I12=9-I12
0089         ENDIF
0090         RKFL=(WDTE(0,1)+WDTE(0,IPM)+WDTE(0,I12))*RLU(0)
0091         DO 120 I=1,MDCY(KFA,3)
0092         IDC=I+MDCY(KFA,2)-1
0093         KFL1(JT)=KFDP(IDC,1)*ISIGN(1,K(ID,2))
0094         KFL2(JT)=KFDP(IDC,2)*ISIGN(1,K(ID,2))
0095         RKFL=RKFL-(WDTE(I,1)+WDTE(I,IPM)+WDTE(I,I12))
0096         IF(RKFL.LE.0.) GOTO 130
0097   120   CONTINUE
0098   130   CONTINUE
0099       ENDIF
0100
0101 C...Summarize result on decay channel chosen.
0102       IF((KFA.EQ.23.OR.KFA.EQ.24).AND.KFL1(JT).EQ.0) NINH=NINH+1
0103       IF(KFL1(JT).EQ.0) GOTO 140
0104       KDCY(JT)=2
0105       IF(IABS(KFL1(JT)).LE.10.OR.KFL1(JT).EQ.21) KDCY(JT)=1
0106       IF((IABS(KFL1(JT)).GE.23.AND.IABS(KFL1(JT)).LE.25).OR.
0107      &(IABS(KFL1(JT)).EQ.37)) KDCY(JT)=3
0108       NSD(JT)=N
0109
0110 C...Fill decay products, prepared for parton showers for quarks.
0111       IF(KDCY(JT).EQ.1) THEN
0112         CALL LU2ENT(-(N+1),KFL1(JT),KFL2(JT),P(ID,5))
0113       ELSE
0114         CALL LU2ENT(N+1,KFL1(JT),KFL2(JT),P(ID,5))
0115       ENDIF
0116       IF(JTMAX.EQ.1) THEN
0117         CTHE(JT)=VINT(13)+(VINT(33)-VINT(13)+VINT(34)-VINT(14))*RLU(0)
0118         IF(CTHE(JT).GT.VINT(33)) CTHE(JT)=CTHE(JT)+VINT(14)-VINT(33)
0119         PHI(JT)=VINT(24)
0120       ELSE
0121         CTHE(JT)=2.*RLU(0)-1.
0122         PHI(JT)=PARU(2)*RLU(0)
0123       ENDIF
0124   140 CONTINUE
0125       IF(MINT(3).EQ.1.AND.IP.EQ.1) THEN
0126         MINT(25)=KFL1(1)
0127         MINT(26)=KFL2(1)
0128       ENDIF
0129       IF(JTMAX.EQ.1.AND.KDCY(1).EQ.0) GOTO 530
0130       IF(JTMAX.EQ.2.AND.KDCY(1).EQ.0.AND.KDCY(2).EQ.0) GOTO 530
0131       IF(MSTP(45).LE.0.OR.IREF(IP,2).EQ.0.OR.NINH.GE.1) GOTO 500
0132       IF(K(IREF(1,1),2).EQ.25.AND.IP.EQ.1) GOTO 500
0133       IF(K(IREF(1,1),2).EQ.25.AND.KDCY(1)*KDCY(2).EQ.0) GOTO 500
0134
0135 C...Order incoming partons and outgoing resonances.
0136       ILIN(1)=MINT(84)+1
0137       IF(K(MINT(84)+1,2).GT.0) ILIN(1)=MINT(84)+2
0138       IF(K(ILIN(1),2).EQ.21) ILIN(1)=2*MINT(84)+3-ILIN(1)
0139       ILIN(2)=2*MINT(84)+3-ILIN(1)
0140       IMIN=1
0141       IF(IREF(IP,5).EQ.25) IMIN=3
0142       IMAX=2
0143       IORD=1
0144       IF(K(IREF(IP,1),2).EQ.23) IORD=2
0145       IF(K(IREF(IP,1),2).EQ.24.AND.K(IREF(IP,2),2).EQ.-24) IORD=2
0146       IF(IABS(K(IREF(IP,IORD),2)).EQ.25) IORD=3-IORD
0147       IF(KDCY(IORD).EQ.0) IORD=3-IORD
0148
0149 C...Order decay products of resonances.
0150       DO 390 JT=IORD,3-IORD,3-2*IORD
0151       IF(KDCY(JT).EQ.0) THEN
0152         ILIN(IMAX+1)=NSD(JT)
0153         IMAX=IMAX+1
0154       ELSEIF(K(NSD(JT)+1,2).GT.0) THEN
0155         ILIN(IMAX+1)=N+2*JT-1
0156         ILIN(IMAX+2)=N+2*JT
0157         IMAX=IMAX+2
0158         K(N+2*JT-1,2)=K(NSD(JT)+1,2)
0159         K(N+2*JT,2)=K(NSD(JT)+2,2)
0160       ELSE
0161         ILIN(IMAX+1)=N+2*JT
0162         ILIN(IMAX+2)=N+2*JT-1
0163         IMAX=IMAX+2
0164         K(N+2*JT-1,2)=K(NSD(JT)+1,2)
0165         K(N+2*JT,2)=K(NSD(JT)+2,2)
0166       ENDIF
0167   390 CONTINUE
0168
0169 C...Find charge, isospin, left- and righthanded couplings.
0170       XW=PARU(102)
0171       DO 410 I=IMIN,IMAX
0172       DO 400 J=1,4
0173   400 COUP(I,J)=0.
0174       KFA=IABS(K(ILIN(I),2))
0175       IF(KFA.GT.20) GOTO 410
0176       COUP(I,1)=LUCHGE(KFA)/3.
0177       COUP(I,2)=(-1)**MOD(KFA,2)
0178       COUP(I,4)=-2.*COUP(I,1)*XW
0179       COUP(I,3)=COUP(I,2)+COUP(I,4)
0180   410 CONTINUE
0181       SQMZ=PMAS(23,1)**2
0182       GZMZ=PMAS(23,1)*PMAS(23,2)
0183       SQMW=PMAS(24,1)**2
0184       GZMW=PMAS(24,1)*PMAS(24,2)
0185       SQMZP=PMAS(32,1)**2
0186       GZMZP=PMAS(32,1)*PMAS(32,2)
0187
0188 C...Select random angles; construct massless four-vectors.
0189   420 DO 430 I=N+1,N+4
0190       K(I,1)=1
0191       DO 430 J=1,5
0192   430 P(I,J)=0.
0193       DO 440 JT=1,JTMAX
0194       IF(KDCY(JT).EQ.0) GOTO 440
0195       ID=IREF(IP,JT)
0196       P(N+2*JT-1,3)=0.5*P(ID,5)
0197       P(N+2*JT-1,4)=0.5*P(ID,5)
0198       P(N+2*JT,3)=-0.5*P(ID,5)
0199       P(N+2*JT,4)=0.5*P(ID,5)
0200       CTHE(JT)=2.*RLU(0)-1.
0201       PHI(JT)=PARU(2)*RLU(0)
0202       CALL LUDBRB(N+2*JT-1,N+2*JT,ACOS(CTHE(JT)),PHI(JT),
0203      &DBLE(P(ID,1)/P(ID,4)),DBLE(P(ID,2)/P(ID,4)),DBLE(P(ID,3)/P(ID,4)))
0204   440 CONTINUE
0205
0206 C...Store incoming and outgoing momenta, with random rotation to
0207 C...avoid accidental zeroes in HA expressions.
0208       DO 450 I=1,IMAX
0209       K(N+4+I,1)=1
0210       P(N+4+I,4)=SQRT(P(ILIN(I),1)**2+P(ILIN(I),2)**2+P(ILIN(I),3)**2+
0211      &P(ILIN(I),5)**2)
0212       P(N+4+I,5)=P(ILIN(I),5)
0213       DO 450 J=1,3
0214   450 P(N+4+I,J)=P(ILIN(I),J)
0215       THERR=ACOS(2.*RLU(0)-1.)
0216       PHIRR=PARU(2)*RLU(0)
0217       CALL LUDBRB(N+5,N+4+IMAX,THERR,PHIRR,0D0,0D0,0D0)
0218       DO 460 I=1,IMAX
0219       DO 460 J=1,4
0220   460 PK(I,J)=P(N+4+I,J)
0221
0222 C...Calculate internal products.
0223       IF(ISUB.EQ.22.OR.ISUB.EQ.23.OR.ISUB.EQ.25) THEN
0224         DO 470 I1=IMIN,IMAX-1
0225         DO 470 I2=I1+1,IMAX
0226         HA(I1,I2)=SQRT((PK(I1,4)-PK(I1,3))*(PK(I2,4)+PK(I2,3))/
0227      &  (1E-20+PK(I1,1)**2+PK(I1,2)**2))*CMPLX(PK(I1,1),PK(I1,2))-
0228      &  SQRT((PK(I1,4)+PK(I1,3))*(PK(I2,4)-PK(I2,3))/
0229      &  (1E-20+PK(I2,1)**2+PK(I2,2)**2))*CMPLX(PK(I2,1),PK(I2,2))
0230         HC(I1,I2)=CONJG(HA(I1,I2))
0231         IF(I1.LE.2) HA(I1,I2)=CMPLX(0.,1.)*HA(I1,I2)
0232         IF(I1.LE.2) HC(I1,I2)=CMPLX(0.,1.)*HC(I1,I2)
0233         HA(I2,I1)=-HA(I1,I2)
0234   470   HC(I2,I1)=-HC(I1,I2)
0235       ENDIF
0236       DO 480 I=1,2
0237       DO 480 J=1,4
0238   480 PK(I,J)=-PK(I,J)
0239       DO 490 I1=IMIN,IMAX-1
0240       DO 490 I2=I1+1,IMAX
0241       PKK(I1,I2)=2.*(PK(I1,4)*PK(I2,4)-PK(I1,1)*PK(I2,1)-
0242      &PK(I1,2)*PK(I2,2)-PK(I1,3)*PK(I2,3))
0243   490 PKK(I2,I1)=PKK(I1,I2)
0244
0245       IF(IREF(IP,5).EQ.25) THEN
0246 C...Angular weight for H0 -> Z0 + Z0 or W+ + W- -> 4 quarks/leptons
0247         WT=16.*PKK(3,5)*PKK(4,6)
0248         IF(IP.EQ.1) WTMAX=SH**2
0249         IF(IP.GE.2) WTMAX=P(IREF(IP,6),5)**4
0250
0251       ELSEIF(ISUB.EQ.1) THEN
0252         IF(KFA.NE.37) THEN
0253 C...Angular weight for gamma*/Z0 -> 2 quarks/leptons
0254           EI=KCHG(IABS(MINT(15)),1)/3.
0255           AI=SIGN(1.,EI+0.1)
0256           VI=AI-4.*EI*XW
0257           EF=KCHG(KFA,1)/3.
0258           AF=SIGN(1.,EF+0.1)
0259           VF=AF-4.*EF*XW
0260           GG=1.
0261           GZ=1./(8.*XW*(1.-XW))*SH*(SH-SQMZ)/((SH-SQMZ)**2+GZMZ**2)
0262           ZZ=1./(16.*XW*(1.-XW))**2*SH**2/((SH-SQMZ)**2+GZMZ**2)
0263           IF(MSTP(43).EQ.1) THEN
0264 C...Only gamma* production included
0265             GZ=0.
0266             ZZ=0.
0267           ELSEIF(MSTP(43).EQ.2) THEN
0268 C...Only Z0 production included
0269             GG=0.
0270             GZ=0.
0271           ENDIF
0272           ASYM=2.*(EI*AI*GZ*EF*AF+4.*VI*AI*ZZ*VF*AF)/(EI**2*GG*EF**2+
0273      &    EI*VI*GZ*EF*VF+(VI**2+AI**2)*ZZ*(VF**2+AF**2))
0274           WT=1.+ASYM*CTHE(JT)+CTHE(JT)**2
0275           WTMAX=2.+ABS(ASYM)
0276         ELSE
0277 C...Angular weight for gamma*/Z0 -> H+ + H-
0278           WT=1.-CTHE(JT)**2
0279           WTMAX=1.
0280         ENDIF
0281
0282       ELSEIF(ISUB.EQ.2) THEN
0283 C...Angular weight for W+/- -> 2 quarks/leptons
0284         WT=(1.+CTHE(JT))**2
0285         WTMAX=4.
0286
0287       ELSEIF(ISUB.EQ.15.OR.ISUB.EQ.19) THEN
0288 C...Angular weight for f + fb -> gluon/gamma + Z0 ->
0289 C...-> gluon/gamma + 2 quarks/leptons
0290         WT=((COUP(1,3)*COUP(3,3))**2+(COUP(1,4)*COUP(3,4))**2)*
0291      &  (PKK(1,3)**2+PKK(2,4)**2)+((COUP(1,3)*COUP(3,4))**2+
0292      &  (COUP(1,4)*COUP(3,3))**2)*(PKK(1,4)**2+PKK(2,3)**2)
0293         WTMAX=(COUP(1,3)**2+COUP(1,4)**2)*(COUP(3,3)**2+COUP(3,4)**2)*
0294      &  ((PKK(1,3)+PKK(1,4))**2+(PKK(2,3)+PKK(2,4))**2)
0295
0296       ELSEIF(ISUB.EQ.16.OR.ISUB.EQ.20) THEN
0297 C...Angular weight for f + fb' -> gluon/gamma + W+/- ->
0298 C...-> gluon/gamma + 2 quarks/leptons
0299         WT=PKK(1,3)**2+PKK(2,4)**2
0300         WTMAX=(PKK(1,3)+PKK(1,4))**2+(PKK(2,3)+PKK(2,4))**2
0301
0302       ELSEIF(ISUB.EQ.22) THEN
0303 C...Angular weight for f + fb -> Z0 + Z0 -> 4 quarks/leptons
0304         S34=P(IREF(IP,IORD),5)**2
0305         S56=P(IREF(IP,3-IORD),5)**2
0306         TI=PKK(1,3)+PKK(1,4)+S34
0307         UI=PKK(1,5)+PKK(1,6)+S56
0308         WT=COUP(1,3)**4*((COUP(3,3)*COUP(5,3)*ABS(FGK(1,2,3,4,5,6)/
0309      &  TI+FGK(1,2,5,6,3,4)/UI))**2+(COUP(3,4)*COUP(5,3)*ABS(
0310      &  FGK(1,2,4,3,5,6)/TI+FGK(1,2,5,6,4,3)/UI))**2+(COUP(3,3)*
0311      &  COUP(5,4)*ABS(FGK(1,2,3,4,6,5)/TI+FGK(1,2,6,5,3,4)/UI))**2+
0312      &  (COUP(3,4)*COUP(5,4)*ABS(FGK(1,2,4,3,6,5)/TI+FGK(1,2,6,5,4,3)/
0313      &  UI))**2)+COUP(1,4)**4*((COUP(3,3)*COUP(5,3)*ABS(
0314      &  FGK(2,1,5,6,3,4)/TI+FGK(2,1,3,4,5,6)/UI))**2+(COUP(3,4)*
0315      &  COUP(5,3)*ABS(FGK(2,1,6,5,3,4)/TI+FGK(2,1,3,4,6,5)/UI))**2+
0316      &  (COUP(3,3)*COUP(5,4)*ABS(FGK(2,1,5,6,4,3)/TI+FGK(2,1,4,3,5,6)/
0317      &  UI))**2+(COUP(3,4)*COUP(5,4)*ABS(FGK(2,1,6,5,4,3)/TI+
0318      &  FGK(2,1,4,3,6,5)/UI))**2)
0319         WTMAX=4.*S34*S56*(COUP(1,3)**4+COUP(1,4)**4)*(COUP(3,3)**2+
0320      &  COUP(3,4)**2)*(COUP(5,3)**2+COUP(5,4)**2)*4.*(TI/UI+UI/TI+
0321      &  2.*SH*(S34+S56)/(TI*UI)-S34*S56*(1./TI**2+1./UI**2))
0322
0323       ELSEIF(ISUB.EQ.23) THEN
0324 C...Angular weight for f + fb' -> Z0 + W +/- -> 4 quarks/leptons
0325         D34=P(IREF(IP,IORD),5)**2
0326         D56=P(IREF(IP,3-IORD),5)**2
0327         DT=PKK(1,3)+PKK(1,4)+D34
0328         DU=PKK(1,5)+PKK(1,6)+D56
0329         CAWZ=COUP(2,3)/SNGL(DT)-2.*(1.-XW)*COUP(1,2)/(SH-SQMW)
0330         CBWZ=COUP(1,3)/SNGL(DU)+2.*(1.-XW)*COUP(1,2)/(SH-SQMW)
0331         WT=COUP(5,3)**2*ABS(CAWZ*FGK(1,2,3,4,5,6)+CBWZ*
0332      &  FGK(1,2,5,6,3,4))**2+COUP(5,4)**2*ABS(CAWZ*
0333      &  FGK(1,2,3,4,6,5)+CBWZ*FGK(1,2,6,5,3,4))**2
0334         WTMAX=4.*D34*D56*(COUP(5,3)**2+COUP(5,4)**2)*(CAWZ**2*
0335      &  DIGK(DT,DU)+CBWZ**2*DIGK(DU,DT)+CAWZ*CBWZ*DJGK(DT,DU))
0336
0337       ELSEIF(ISUB.EQ.24) THEN
0338 C...Angular weight for f + fb -> Z0 + H0 -> 2 quarks/leptons + H0
0339         WT=((COUP(1,3)*COUP(3,3))**2+(COUP(1,4)*COUP(3,4))**2)*
0340      &  PKK(1,3)*PKK(2,4)+((COUP(1,3)*COUP(3,4))**2+(COUP(1,4)*
0341      &  COUP(3,3))**2)*PKK(1,4)*PKK(2,3)
0342         WTMAX=(COUP(1,3)**2+COUP(1,4)**2)*(COUP(3,3)**2+COUP(3,4)**2)*
0343      &  (PKK(1,3)+PKK(1,4))*(PKK(2,3)+PKK(2,4))
0344
0345       ELSEIF(ISUB.EQ.25) THEN
0346 C...Angular weight for f + fb -> W+ + W- -> 4 quarks/leptons
0347         D34=P(IREF(IP,IORD),5)**2
0348         D56=P(IREF(IP,3-IORD),5)**2
0349         DT=PKK(1,3)+PKK(1,4)+D34
0350         DU=PKK(1,5)+PKK(1,6)+D56
0351         CDWW=(COUP(1,3)*SQMZ/(SH-SQMZ)+COUP(1,2))/SH
0352         CAWW=CDWW+0.5*(COUP(1,2)+1.)/SNGL(DT)
0353         CBWW=CDWW+0.5*(COUP(1,2)-1.)/SNGL(DU)
0354         CCWW=COUP(1,4)*SQMZ/(SH-SQMZ)/SH
0355         WT=ABS(CAWW*FGK(1,2,3,4,5,6)-CBWW*FGK(1,2,5,6,3,4))**2+
0356      &  CCWW**2*ABS(FGK(2,1,5,6,3,4)-FGK(2,1,3,4,5,6))**2
0357         WTMAX=4.*D34*D56*(CAWW**2*DIGK(DT,DU)+CBWW**2*DIGK(DU,DT)-CAWW*
0358      &  CBWW*DJGK(DT,DU)+CCWW**2*(DIGK(DT,DU)+DIGK(DU,DT)-DJGK(DT,DU)))
0359
0360       ELSEIF(ISUB.EQ.26) THEN
0361 C...Angular weight for f + fb' -> W+/- + H0 -> 2 quarks/leptons + H0
0362         WT=PKK(1,3)*PKK(2,4)
0363         WTMAX=(PKK(1,3)+PKK(1,4))*(PKK(2,3)+PKK(2,4))
0364
0365       ELSEIF(ISUB.EQ.30) THEN
0366 C...Angular weight for f + g -> f + Z0 -> f + 2 quarks/leptons
0367         IF(K(ILIN(1),2).GT.0) WT=((COUP(1,3)*COUP(3,3))**2+
0368      &  (COUP(1,4)*COUP(3,4))**2)*(PKK(1,4)**2+PKK(3,5)**2)+
0369      &  ((COUP(1,3)*COUP(3,4))**2+(COUP(1,4)*COUP(3,3))**2)*
0370      &  (PKK(1,3)**2+PKK(4,5)**2)
0371         IF(K(ILIN(1),2).LT.0) WT=((COUP(1,3)*COUP(3,3))**2+
0372      &  (COUP(1,4)*COUP(3,4))**2)*(PKK(1,3)**2+PKK(4,5)**2)+
0373      &  ((COUP(1,3)*COUP(3,4))**2+(COUP(1,4)*COUP(3,3))**2)*
0374      &  (PKK(1,4)**2+PKK(3,5)**2)
0375         WTMAX=(COUP(1,3)**2+COUP(1,4)**2)*(COUP(3,3)**2+COUP(3,4)**2)*
0376      &  ((PKK(1,3)+PKK(1,4))**2+(PKK(3,5)+PKK(4,5))**2)
0377
0378       ELSEIF(ISUB.EQ.31) THEN
0379 C...Angular weight for f + g -> f' + W+/- -> f' + 2 quarks/leptons
0380         IF(K(ILIN(1),2).GT.0) WT=PKK(1,4)**2+PKK(3,5)**2
0381         IF(K(ILIN(1),2).LT.0) WT=PKK(1,3)**2+PKK(4,5)**2
0382         WTMAX=(PKK(1,3)+PKK(1,4))**2+(PKK(3,5)+PKK(4,5))**2
0383
0384       ELSEIF(ISUB.EQ.141) THEN
0385 C...Angular weight for gamma*/Z0/Z'0 -> 2 quarks/leptons
0386         EI=KCHG(IABS(MINT(15)),1)/3.
0387         AI=SIGN(1.,EI+0.1)
0388         VI=AI-4.*EI*XW
0389         API=SIGN(1.,EI+0.1)
0390         VPI=API-4.*EI*XW
0391         EF=KCHG(KFA,1)/3.
0392         AF=SIGN(1.,EF+0.1)
0393         VF=AF-4.*EF*XW
0394         APF=SIGN(1.,EF+0.1)
0395         VPF=APF-4.*EF*XW
0396         GG=1.
0397         GZ=1./(8.*XW*(1.-XW))*SH*(SH-SQMZ)/((SH-SQMZ)**2+GZMZ**2)
0398         GZP=1./(8.*XW*(1.-XW))*SH*(SH-SQMZP)/((SH-SQMZP)**2+GZMZP**2)
0399         ZZ=1./(16.*XW*(1.-XW))**2*SH**2/((SH-SQMZ)**2+GZMZ**2)
0400         ZZP=2./(16.*XW*(1.-XW))**2*
0401      &  SH**2*((SH-SQMZ)*(SH-SQMZP)+GZMZ*GZMZP)/
0402      &  (((SH-SQMZ)**2+GZMZ**2)*((SH-SQMZP)**2+GZMZP**2))
0403         ZPZP=1./(16.*XW*(1.-XW))**2*SH**2/((SH-SQMZP)**2+GZMZP**2)
0404         IF(MSTP(44).EQ.1) THEN
0405 C...Only gamma* production included
0406           GZ=0.
0407           GZP=0.
0408           ZZ=0.
0409           ZZP=0.
0410           ZPZP=0.
0411         ELSEIF(MSTP(44).EQ.2) THEN
0412 C...Only Z0 production included
0413           GG=0.
0414           GZ=0.
0415           GZP=0.
0416           ZZP=0.
0417           ZPZP=0.
0418         ELSEIF(MSTP(44).EQ.3) THEN
0419 C...Only Z'0 production included
0420           GG=0.
0421           GZ=0.
0422           GZP=0.
0423           ZZ=0.
0424           ZZP=0.
0425         ELSEIF(MSTP(44).EQ.4) THEN
0426 C...Only gamma*/Z0 production included
0427           GZP=0.
0428           ZZP=0.
0429           ZPZP=0.
0430         ELSEIF(MSTP(44).EQ.5) THEN
0431 C...Only gamma*/Z'0 production included
0432           GZ=0.
0433           ZZ=0.
0434           ZZP=0.
0435         ELSEIF(MSTP(44).EQ.6) THEN
0436 C...Only Z0/Z'0 production included
0437           GG=0.
0438           GZ=0.
0439           GZP=0.
0440         ENDIF
0441         ASYM=2.*(EI*AI*GZ*EF*AF+EI*API*GZP*EF*APF+4.*VI*AI*ZZ*VF*AF+
0442      &  (VI*API+VPI*AI)*ZZP*(VF*APF+VPF*AF)+4.*VPI*API*ZPZP*VPF*APF)/
0443      &  (EI**2*GG*EF**2+EI*VI*GZ*EF*VF+EI*VPI*GZP*EF*VPF+
0444      &  (VI**2+AI**2)*ZZ*(VF**2+AF**2)+(VI*VPI+AI*API)*ZZP*
0445      &  (VF*VPF+AF*APF)+(VPI**2+API**2)*ZPZP*(VPF**2+APF**2))
0446         WT=1.+ASYM*CTHE(JT)+CTHE(JT)**2
0447         WTMAX=2.+ABS(ASYM)
0448
0449       ELSE
0450         WT=1.
0451         WTMAX=1.
0452       ENDIF
0453 C...Obtain correct angular distribution by rejection techniques.
0454       IF(WT.LT.RLU(0)*WTMAX) GOTO 420
0455
0456 C...Construct massive four-vectors using angles chosen. Mark decayed
0457 C...resonances, add documentation lines. Shower evolution.
0458   500 DO 520 JT=1,JTMAX
0459       IF(KDCY(JT).EQ.0) GOTO 520
0460       ID=IREF(IP,JT)
0461       CALL LUDBRB(NSD(JT)+1,NSD(JT)+2,ACOS(CTHE(JT)),PHI(JT),
0462      &DBLE(P(ID,1)/P(ID,4)),DBLE(P(ID,2)/P(ID,4)),DBLE(P(ID,3)/P(ID,4)))
0463       K(ID,1)=K(ID,1)+10
0464       K(ID,4)=NSD(JT)+1
0465       K(ID,5)=NSD(JT)+2
0466       IDOC=MINT(83)+MINT(4)
0467       DO 510 I=NSD(JT)+1,NSD(JT)+2
0468       MINT(4)=MINT(4)+1
0469       I1=MINT(83)+MINT(4)
0470       K(I,3)=I1
0471       K(I1,1)=21
0472       K(I1,2)=K(I,2)
0473       K(I1,3)=IREF(IP,JT+2)
0474       DO 510 J=1,5
0475   510 P(I1,J)=P(I,J)
0476       IF(JTMAX.EQ.1) THEN
0477         MINT(7)=MINT(83)+6+2*ISET(ISUB)
0478         MINT(8)=MINT(83)+7+2*ISET(ISUB)
0479       ENDIF
0480       IF(MSTP(71).GE.1.AND.KDCY(JT).EQ.1) CALL LUSHOW(NSD(JT)+1,
0481      &NSD(JT)+2,P(ID,5))
0482
0483 C...Check if new resonances were produced, loop back if needed.
0484       IF(KDCY(JT).NE.3) GOTO 520
0485       NP=NP+1
0486       IREF(NP,1)=NSD(JT)+1
0487       IREF(NP,2)=NSD(JT)+2
0488       IREF(NP,3)=IDOC+1
0489       IREF(NP,4)=IDOC+2
0490       IREF(NP,5)=K(IREF(IP,JT),2)
0491       IREF(NP,6)=IREF(IP,JT)
0492   520 CONTINUE
0493   530 IF(IP.LT.NP) GOTO 100
0494
0495       RETURN
0496       END