HIJING/hipyset1_35/pyresd_hijing.F

   1 * $Id$
   2
   3 C*********************************************************************
   4
   5       SUBROUTINE PYRESD_HIJING
   6
   7 C...Allows resonances to decay (including parton showers for hadronic
   8 C...channels).
   9       IMPLICIT DOUBLE PRECISION(D)
  10 #include "lujets_hijing.inc"
  11 #include "ludat1_hijing.inc"
  12 #include "ludat2_hijing.inc"
  13 #include "ludat3_hijing.inc"
  14 #include "pysubs_hijing.inc"
  15 #include "pypars_hijing.inc"
  16 #include "pyint1_hijing.inc"
  17 #include "pyint2_hijing.inc"
  18 #include "pyint4_hijing.inc"
  19       DIMENSION IREF(10,6),KDCY(2),KFL1(2),KFL2(2),NSD(2),ILIN(6),
  20      &COUP(6,4),PK(6,4),PKK(6,6),CTHE(2),PHI(2),WDTP(0:40),
  21      &WDTE(0:40,0:5)
  22       COMPLEX FGK,HA(6,6),HC(6,6)
  23
  24 C...The F, Xi and Xj functions of Gunion and Kunszt
  25 C...(Phys. Rev. D33, 665, plus errata from the authors).
  26       FGK(I1,I2,I3,I4,I5,I6)=4.*HA(I1,I3)*HC(I2,I6)*(HA(I1,I5)*
  27      &HC(I1,I4)+HA(I3,I5)*HC(I3,I4))
  28       DIGK(DT,DU)=-4.*D34*D56+DT*(3.*DT+4.*DU)+DT**2*(DT*DU/(D34*D56)-
  29      &2.*(1./D34+1./D56)*(DT+DU)+2.*(D34/D56+D56/D34))
  30       DJGK(DT,DU)=8.*(D34+D56)**2-8.*(D34+D56)*(DT+DU)-6.*DT*DU-
  31      &2.*DT*DU*(DT*DU/(D34*D56)-2.*(1./D34+1./D56)*(DT+DU)+
  32      &2.*(D34/D56+D56/D34))
  33
  34 C...Define initial two objects, initialize loop.
  35       ISUB=MINT(1)
  36       SH=VINT(44)
  37       IREF(1,5)=0
  38       IREF(1,6)=0
  39       IF(ISET(ISUB).EQ.1.OR.ISET(ISUB).EQ.3) THEN
  40         IREF(1,1)=MINT(84)+2+ISET(ISUB)
  41         IREF(1,2)=0
  42         IREF(1,3)=MINT(83)+6+ISET(ISUB)
  43         IREF(1,4)=0
  44       ELSEIF(ISET(ISUB).EQ.2.OR.ISET(ISUB).EQ.4) THEN
  45         IREF(1,1)=MINT(84)+1+ISET(ISUB)
  46         IREF(1,2)=MINT(84)+2+ISET(ISUB)
  47         IREF(1,3)=MINT(83)+5+ISET(ISUB)
  48         IREF(1,4)=MINT(83)+6+ISET(ISUB)
  49       ENDIF
  50       NP=1
  51       IP=0
  52   100 IP=IP+1
  53       NINH=0
  54
  55 C...Loop over one/two resonances; reset decay rates.
  56       JTMAX=2
  57       IF(IP.EQ.1.AND.(ISET(ISUB).EQ.1.OR.ISET(ISUB).EQ.3)) JTMAX=1
  58       DO 140 JT=1,JTMAX
  59       KDCY(JT)=0
  60       KFL1(JT)=0
  61       KFL2(JT)=0
  62       NSD(JT)=IREF(IP,JT)
  63       ID=IREF(IP,JT)
  64       IF(ID.EQ.0) GOTO 140
  65       KFA=IABS(K(ID,2))
  66       IF(KFA.LT.23.OR.KFA.GT.40) GOTO 140
  67       IF(MDCY(KFA,1).NE.0) THEN
  68         IF(ISUB.EQ.1.OR.ISUB.EQ.141) MINT(61)=1
  69         CALL PYWIDT_HIJING(KFA,P(ID,5),WDTP,WDTE)
  70         IF(KCHG(KFA,3).EQ.0) THEN
  71           IPM=2
  72         ELSE
  73           IPM=(5+ISIGN(1,K(ID,2)))/2
  74         ENDIF
  75         IF(JTMAX.EQ.1.OR.IABS(K(IREF(IP,1),2)).NE.IABS(K(IREF(IP,2),2)))
  76      &  THEN
  77           I12=4
  78         ELSE
  79           IF(JT.EQ.1) I12=INT(4.5+RLU_HIJING(0))
  80           I12=9-I12
  81         ENDIF
  82         RKFL=(WDTE(0,1)+WDTE(0,IPM)+WDTE(0,I12))*RLU_HIJING(0)
  83         DO 120 I=1,MDCY(KFA,3)
  84         IDC=I+MDCY(KFA,2)-1
  85         KFL1(JT)=KFDP(IDC,1)*ISIGN(1,K(ID,2))
  86         KFL2(JT)=KFDP(IDC,2)*ISIGN(1,K(ID,2))
  87         RKFL=RKFL-(WDTE(I,1)+WDTE(I,IPM)+WDTE(I,I12))
  88         IF(RKFL.LE.0.) GOTO 130
  89   120   CONTINUE
  90   130   CONTINUE
  91       ENDIF
  92
  93 C...Summarize result on decay channel chosen.
  94       IF((KFA.EQ.23.OR.KFA.EQ.24).AND.KFL1(JT).EQ.0) NINH=NINH+1
  95       IF(KFL1(JT).EQ.0) GOTO 140
  96       KDCY(JT)=2
  97       IF(IABS(KFL1(JT)).LE.10.OR.KFL1(JT).EQ.21) KDCY(JT)=1
  98       IF((IABS(KFL1(JT)).GE.23.AND.IABS(KFL1(JT)).LE.25).OR.
  99      &(IABS(KFL1(JT)).EQ.37)) KDCY(JT)=3
 100       NSD(JT)=N
 101
 102 C...Fill decay products, prepared for parton showers for quarks.
 103       IF(KDCY(JT).EQ.1) THEN
 104         CALL LU2ENT_HIJING(-(N+1),KFL1(JT),KFL2(JT),P(ID,5))
 105       ELSE
 106         CALL LU2ENT_HIJING(N+1,KFL1(JT),KFL2(JT),P(ID,5))
 107       ENDIF
 108       IF(JTMAX.EQ.1) THEN
 109          CTHE(JT)=VINT(13)+(VINT(33)-VINT(13)+VINT(34)-VINT(14))
 110      $        *RLU_HIJING(0)
 111         IF(CTHE(JT).GT.VINT(33)) CTHE(JT)=CTHE(JT)+VINT(14)-VINT(33)
 112         PHI(JT)=VINT(24)
 113       ELSE
 114         CTHE(JT)=2.*RLU_HIJING(0)-1.
 115         PHI(JT)=PARU(2)*RLU_HIJING(0)
 116       ENDIF
 117   140 CONTINUE
 118       IF(MINT(3).EQ.1.AND.IP.EQ.1) THEN
 119         MINT(25)=KFL1(1)
 120         MINT(26)=KFL2(1)
 121       ENDIF
 122       IF(JTMAX.EQ.1.AND.KDCY(1).EQ.0) GOTO 530
 123       IF(JTMAX.EQ.2.AND.KDCY(1).EQ.0.AND.KDCY(2).EQ.0) GOTO 530
 124       IF(MSTP(45).LE.0.OR.IREF(IP,2).EQ.0.OR.NINH.GE.1) GOTO 500
 125       IF(K(IREF(1,1),2).EQ.25.AND.IP.EQ.1) GOTO 500
 126       IF(K(IREF(1,1),2).EQ.25.AND.KDCY(1)*KDCY(2).EQ.0) GOTO 500
 127
 128 C...Order incoming partons and outgoing resonances.
 129       ILIN(1)=MINT(84)+1
 130       IF(K(MINT(84)+1,2).GT.0) ILIN(1)=MINT(84)+2
 131       IF(K(ILIN(1),2).EQ.21) ILIN(1)=2*MINT(84)+3-ILIN(1)
 132       ILIN(2)=2*MINT(84)+3-ILIN(1)
 133       IMIN=1
 134       IF(IREF(IP,5).EQ.25) IMIN=3
 135       IMAX=2
 136       IORD=1
 137       IF(K(IREF(IP,1),2).EQ.23) IORD=2
 138       IF(K(IREF(IP,1),2).EQ.24.AND.K(IREF(IP,2),2).EQ.-24) IORD=2
 139       IF(IABS(K(IREF(IP,IORD),2)).EQ.25) IORD=3-IORD
 140       IF(KDCY(IORD).EQ.0) IORD=3-IORD
 141
 142 C...Order decay products of resonances.
 143       DO 390 JT=IORD,3-IORD,3-2*IORD
 144       IF(KDCY(JT).EQ.0) THEN
 145         ILIN(IMAX+1)=NSD(JT)
 146         IMAX=IMAX+1
 147       ELSEIF(K(NSD(JT)+1,2).GT.0) THEN
 148         ILIN(IMAX+1)=N+2*JT-1
 149         ILIN(IMAX+2)=N+2*JT
 150         IMAX=IMAX+2
 151         K(N+2*JT-1,2)=K(NSD(JT)+1,2)
 152         K(N+2*JT,2)=K(NSD(JT)+2,2)
 153       ELSE
 154         ILIN(IMAX+1)=N+2*JT
 155         ILIN(IMAX+2)=N+2*JT-1
 156         IMAX=IMAX+2
 157         K(N+2*JT-1,2)=K(NSD(JT)+1,2)
 158         K(N+2*JT,2)=K(NSD(JT)+2,2)
 159       ENDIF
 160   390 CONTINUE
 161
 162 C...Find charge, isospin, left- and righthanded couplings.
 163       XW=PARU(102)
 164       DO 410 I=IMIN,IMAX
 165       DO 400 J=1,4
 166   400 COUP(I,J)=0.
 167       KFA=IABS(K(ILIN(I),2))
 168       IF(KFA.GT.20) GOTO 410
 169       COUP(I,1)=LUCHGE_HIJING(KFA)/3.
 170       COUP(I,2)=(-1)**MOD(KFA,2)
 171       COUP(I,4)=-2.*COUP(I,1)*XW
 172       COUP(I,3)=COUP(I,2)+COUP(I,4)
 173   410 CONTINUE
 174       SQMZ=PMAS(23,1)**2
 175       GZMZ=PMAS(23,1)*PMAS(23,2)
 176       SQMW=PMAS(24,1)**2
 177       GZMW=PMAS(24,1)*PMAS(24,2)
 178       SQMZP=PMAS(32,1)**2
 179       GZMZP=PMAS(32,1)*PMAS(32,2)
 180
 181 C...Select random angles; construct massless four-vectors.
 182   420 DO 430 I=N+1,N+4
 183       K(I,1)=1
 184       DO 430 J=1,5
 185   430 P(I,J)=0.
 186       DO 440 JT=1,JTMAX
 187       IF(KDCY(JT).EQ.0) GOTO 440
 188       ID=IREF(IP,JT)
 189       P(N+2*JT-1,3)=0.5*P(ID,5)
 190       P(N+2*JT-1,4)=0.5*P(ID,5)
 191       P(N+2*JT,3)=-0.5*P(ID,5)
 192       P(N+2*JT,4)=0.5*P(ID,5)
 193       CTHE(JT)=2.*RLU_HIJING(0)-1.
 194       PHI(JT)=PARU(2)*RLU_HIJING(0)
 195       CALL LUDBRB_HIJING(N+2*JT-1,N+2*JT,ACOS(CTHE(JT)),PHI(JT),
 196      &DBLE(P(ID,1)/P(ID,4)),DBLE(P(ID,2)/P(ID,4)),DBLE(P(ID,3)/P(ID,4)))
 197   440 CONTINUE
 198
 199 C...Store incoming and outgoing momenta, with random rotation to
 200 C...avoid accidental zeroes in HA expressions.
 201       DO 450 I=1,IMAX
 202       K(N+4+I,1)=1
 203       P(N+4+I,4)=SQRT(P(ILIN(I),1)**2+P(ILIN(I),2)**2+P(ILIN(I),3)**2+
 204      &P(ILIN(I),5)**2)
 205       P(N+4+I,5)=P(ILIN(I),5)
 206       DO 450 J=1,3
 207   450 P(N+4+I,J)=P(ILIN(I),J)
 208       THERR=ACOS(2.*RLU_HIJING(0)-1.)
 209       PHIRR=PARU(2)*RLU_HIJING(0)
 210       CALL LUDBRB_HIJING(N+5,N+4+IMAX,THERR,PHIRR,0D0,0D0,0D0)
 211       DO 460 I=1,IMAX
 212       DO 460 J=1,4
 213   460 PK(I,J)=P(N+4+I,J)
 214
 215 C...Calculate internal products.
 216       IF(ISUB.EQ.22.OR.ISUB.EQ.23.OR.ISUB.EQ.25) THEN
 217         DO 470 I1=IMIN,IMAX-1
 218         DO 470 I2=I1+1,IMAX
 219         HA(I1,I2)=SQRT((PK(I1,4)-PK(I1,3))*(PK(I2,4)+PK(I2,3))/
 220      &  (1E-20+PK(I1,1)**2+PK(I1,2)**2))*CMPLX(PK(I1,1),PK(I1,2))-
 221      &  SQRT((PK(I1,4)+PK(I1,3))*(PK(I2,4)-PK(I2,3))/
 222      &  (1E-20+PK(I2,1)**2+PK(I2,2)**2))*CMPLX(PK(I2,1),PK(I2,2))
 223         HC(I1,I2)=CONJG(HA(I1,I2))
 224         IF(I1.LE.2) HA(I1,I2)=CMPLX(0.,1.)*HA(I1,I2)
 225         IF(I1.LE.2) HC(I1,I2)=CMPLX(0.,1.)*HC(I1,I2)
 226         HA(I2,I1)=-HA(I1,I2)
 227   470   HC(I2,I1)=-HC(I1,I2)
 228       ENDIF
 229       DO 480 I=1,2
 230       DO 480 J=1,4
 231   480 PK(I,J)=-PK(I,J)
 232       DO 490 I1=IMIN,IMAX-1
 233       DO 490 I2=I1+1,IMAX
 234       PKK(I1,I2)=2.*(PK(I1,4)*PK(I2,4)-PK(I1,1)*PK(I2,1)-
 235      &PK(I1,2)*PK(I2,2)-PK(I1,3)*PK(I2,3))
 236   490 PKK(I2,I1)=PKK(I1,I2)
 237
 238       IF(IREF(IP,5).EQ.25) THEN
 239 C...Angular weight for H0 -> Z0 + Z0 or W+ + W- -> 4 quarks/leptons
 240         WT=16.*PKK(3,5)*PKK(4,6)
 241         IF(IP.EQ.1) WTMAX=SH**2
 242         IF(IP.GE.2) WTMAX=P(IREF(IP,6),5)**4
 243
 244       ELSEIF(ISUB.EQ.1) THEN
 245         IF(KFA.NE.37) THEN
 246 C...Angular weight for gamma*/Z0 -> 2 quarks/leptons
 247           EI=KCHG(IABS(MINT(15)),1)/3.
 248           AI=SIGN(1.,EI+0.1)
 249           VI=AI-4.*EI*XW
 250           EF=KCHG(KFA,1)/3.
 251           AF=SIGN(1.,EF+0.1)
 252           VF=AF-4.*EF*XW
 253           GG=1.
 254           GZ=1./(8.*XW*(1.-XW))*SH*(SH-SQMZ)/((SH-SQMZ)**2+GZMZ**2)
 255           ZZ=1./(16.*XW*(1.-XW))**2*SH**2/((SH-SQMZ)**2+GZMZ**2)
 256           IF(MSTP(43).EQ.1) THEN
 257 C...Only gamma* production included
 258             GZ=0.
 259             ZZ=0.
 260           ELSEIF(MSTP(43).EQ.2) THEN
 261 C...Only Z0 production included
 262             GG=0.
 263             GZ=0.
 264           ENDIF
 265           ASYM=2.*(EI*AI*GZ*EF*AF+4.*VI*AI*ZZ*VF*AF)/(EI**2*GG*EF**2+
 266      &    EI*VI*GZ*EF*VF+(VI**2+AI**2)*ZZ*(VF**2+AF**2))
 267           WT=1.+ASYM*CTHE(JT)+CTHE(JT)**2
 268           WTMAX=2.+ABS(ASYM)
 269         ELSE
 270 C...Angular weight for gamma*/Z0 -> H+ + H-
 271           WT=1.-CTHE(JT)**2
 272           WTMAX=1.
 273         ENDIF
 274
 275       ELSEIF(ISUB.EQ.2) THEN
 276 C...Angular weight for W+/- -> 2 quarks/leptons
 277         WT=(1.+CTHE(JT))**2
 278         WTMAX=4.
 279
 280       ELSEIF(ISUB.EQ.15.OR.ISUB.EQ.19) THEN
 281 C...Angular weight for f + fb -> gluon/gamma + Z0 ->
 282 C...-> gluon/gamma + 2 quarks/leptons
 283         WT=((COUP(1,3)*COUP(3,3))**2+(COUP(1,4)*COUP(3,4))**2)*
 284      &  (PKK(1,3)**2+PKK(2,4)**2)+((COUP(1,3)*COUP(3,4))**2+
 285      &  (COUP(1,4)*COUP(3,3))**2)*(PKK(1,4)**2+PKK(2,3)**2)
 286         WTMAX=(COUP(1,3)**2+COUP(1,4)**2)*(COUP(3,3)**2+COUP(3,4)**2)*
 287      &  ((PKK(1,3)+PKK(1,4))**2+(PKK(2,3)+PKK(2,4))**2)
 288
 289       ELSEIF(ISUB.EQ.16.OR.ISUB.EQ.20) THEN
 290 C...Angular weight for f + fb' -> gluon/gamma + W+/- ->
 291 C...-> gluon/gamma + 2 quarks/leptons
 292         WT=PKK(1,3)**2+PKK(2,4)**2
 293         WTMAX=(PKK(1,3)+PKK(1,4))**2+(PKK(2,3)+PKK(2,4))**2
 294
 295       ELSEIF(ISUB.EQ.22) THEN
 296 C...Angular weight for f + fb -> Z0 + Z0 -> 4 quarks/leptons
 297         S34=P(IREF(IP,IORD),5)**2
 298         S56=P(IREF(IP,3-IORD),5)**2
 299         TI=PKK(1,3)+PKK(1,4)+S34
 300         UI=PKK(1,5)+PKK(1,6)+S56
 301         WT=COUP(1,3)**4*((COUP(3,3)*COUP(5,3)*ABS(FGK(1,2,3,4,5,6)/
 302      &  TI+FGK(1,2,5,6,3,4)/UI))**2+(COUP(3,4)*COUP(5,3)*ABS(
 303      &  FGK(1,2,4,3,5,6)/TI+FGK(1,2,5,6,4,3)/UI))**2+(COUP(3,3)*
 304      &  COUP(5,4)*ABS(FGK(1,2,3,4,6,5)/TI+FGK(1,2,6,5,3,4)/UI))**2+
 305      &  (COUP(3,4)*COUP(5,4)*ABS(FGK(1,2,4,3,6,5)/TI+FGK(1,2,6,5,4,3)/
 306      &  UI))**2)+COUP(1,4)**4*((COUP(3,3)*COUP(5,3)*ABS(
 307      &  FGK(2,1,5,6,3,4)/TI+FGK(2,1,3,4,5,6)/UI))**2+(COUP(3,4)*
 308      &  COUP(5,3)*ABS(FGK(2,1,6,5,3,4)/TI+FGK(2,1,3,4,6,5)/UI))**2+
 309      &  (COUP(3,3)*COUP(5,4)*ABS(FGK(2,1,5,6,4,3)/TI+FGK(2,1,4,3,5,6)/
 310      &  UI))**2+(COUP(3,4)*COUP(5,4)*ABS(FGK(2,1,6,5,4,3)/TI+
 311      &  FGK(2,1,4,3,6,5)/UI))**2)
 312         WTMAX=4.*S34*S56*(COUP(1,3)**4+COUP(1,4)**4)*(COUP(3,3)**2+
 313      &  COUP(3,4)**2)*(COUP(5,3)**2+COUP(5,4)**2)*4.*(TI/UI+UI/TI+
 314      &  2.*SH*(S34+S56)/(TI*UI)-S34*S56*(1./TI**2+1./UI**2))
 315
 316       ELSEIF(ISUB.EQ.23) THEN
 317 C...Angular weight for f + fb' -> Z0 + W +/- -> 4 quarks/leptons
 318         D34=P(IREF(IP,IORD),5)**2
 319         D56=P(IREF(IP,3-IORD),5)**2
 320         DT=PKK(1,3)+PKK(1,4)+D34
 321         DU=PKK(1,5)+PKK(1,6)+D56
 322         CAWZ=COUP(2,3)/SNGL(DT)-2.*(1.-XW)*COUP(1,2)/(SH-SQMW)
 323         CBWZ=COUP(1,3)/SNGL(DU)+2.*(1.-XW)*COUP(1,2)/(SH-SQMW)
 324         WT=COUP(5,3)**2*ABS(CAWZ*FGK(1,2,3,4,5,6)+CBWZ*
 325      &  FGK(1,2,5,6,3,4))**2+COUP(5,4)**2*ABS(CAWZ*
 326      &  FGK(1,2,3,4,6,5)+CBWZ*FGK(1,2,6,5,3,4))**2
 327         WTMAX=4.*D34*D56*(COUP(5,3)**2+COUP(5,4)**2)*(CAWZ**2*
 328      &  DIGK(DT,DU)+CBWZ**2*DIGK(DU,DT)+CAWZ*CBWZ*DJGK(DT,DU))
 329
 330       ELSEIF(ISUB.EQ.24) THEN
 331 C...Angular weight for f + fb -> Z0 + H0 -> 2 quarks/leptons + H0
 332         WT=((COUP(1,3)*COUP(3,3))**2+(COUP(1,4)*COUP(3,4))**2)*
 333      &  PKK(1,3)*PKK(2,4)+((COUP(1,3)*COUP(3,4))**2+(COUP(1,4)*
 334      &  COUP(3,3))**2)*PKK(1,4)*PKK(2,3)
 335         WTMAX=(COUP(1,3)**2+COUP(1,4)**2)*(COUP(3,3)**2+COUP(3,4)**2)*
 336      &  (PKK(1,3)+PKK(1,4))*(PKK(2,3)+PKK(2,4))
 337
 338       ELSEIF(ISUB.EQ.25) THEN
 339 C...Angular weight for f + fb -> W+ + W- -> 4 quarks/leptons
 340         D34=P(IREF(IP,IORD),5)**2
 341         D56=P(IREF(IP,3-IORD),5)**2
 342         DT=PKK(1,3)+PKK(1,4)+D34
 343         DU=PKK(1,5)+PKK(1,6)+D56
 344         CDWW=(COUP(1,3)*SQMZ/(SH-SQMZ)+COUP(1,2))/SH
 345         CAWW=CDWW+0.5*(COUP(1,2)+1.)/SNGL(DT)
 346         CBWW=CDWW+0.5*(COUP(1,2)-1.)/SNGL(DU)
 347         CCWW=COUP(1,4)*SQMZ/(SH-SQMZ)/SH
 348         WT=ABS(CAWW*FGK(1,2,3,4,5,6)-CBWW*FGK(1,2,5,6,3,4))**2+
 349      &  CCWW**2*ABS(FGK(2,1,5,6,3,4)-FGK(2,1,3,4,5,6))**2
 350         WTMAX=4.*D34*D56*(CAWW**2*DIGK(DT,DU)+CBWW**2*DIGK(DU,DT)-CAWW*
 351      &  CBWW*DJGK(DT,DU)+CCWW**2*(DIGK(DT,DU)+DIGK(DU,DT)-DJGK(DT,DU)))
 352
 353       ELSEIF(ISUB.EQ.26) THEN
 354 C...Angular weight for f + fb' -> W+/- + H0 -> 2 quarks/leptons + H0
 355         WT=PKK(1,3)*PKK(2,4)
 356         WTMAX=(PKK(1,3)+PKK(1,4))*(PKK(2,3)+PKK(2,4))
 357
 358       ELSEIF(ISUB.EQ.30) THEN
 359 C...Angular weight for f + g -> f + Z0 -> f + 2 quarks/leptons
 360         IF(K(ILIN(1),2).GT.0) WT=((COUP(1,3)*COUP(3,3))**2+
 361      &  (COUP(1,4)*COUP(3,4))**2)*(PKK(1,4)**2+PKK(3,5)**2)+
 362      &  ((COUP(1,3)*COUP(3,4))**2+(COUP(1,4)*COUP(3,3))**2)*
 363      &  (PKK(1,3)**2+PKK(4,5)**2)
 364         IF(K(ILIN(1),2).LT.0) WT=((COUP(1,3)*COUP(3,3))**2+
 365      &  (COUP(1,4)*COUP(3,4))**2)*(PKK(1,3)**2+PKK(4,5)**2)+
 366      &  ((COUP(1,3)*COUP(3,4))**2+(COUP(1,4)*COUP(3,3))**2)*
 367      &  (PKK(1,4)**2+PKK(3,5)**2)
 368         WTMAX=(COUP(1,3)**2+COUP(1,4)**2)*(COUP(3,3)**2+COUP(3,4)**2)*
 369      &  ((PKK(1,3)+PKK(1,4))**2+(PKK(3,5)+PKK(4,5))**2)
 370
 371       ELSEIF(ISUB.EQ.31) THEN
 372 C...Angular weight for f + g -> f' + W+/- -> f' + 2 quarks/leptons
 373         IF(K(ILIN(1),2).GT.0) WT=PKK(1,4)**2+PKK(3,5)**2
 374         IF(K(ILIN(1),2).LT.0) WT=PKK(1,3)**2+PKK(4,5)**2
 375         WTMAX=(PKK(1,3)+PKK(1,4))**2+(PKK(3,5)+PKK(4,5))**2
 376
 377       ELSEIF(ISUB.EQ.141) THEN
 378 C...Angular weight for gamma*/Z0/Z'0 -> 2 quarks/leptons
 379         EI=KCHG(IABS(MINT(15)),1)/3.
 380         AI=SIGN(1.,EI+0.1)
 381         VI=AI-4.*EI*XW
 382         API=SIGN(1.,EI+0.1)
 383         VPI=API-4.*EI*XW
 384         EF=KCHG(KFA,1)/3.
 385         AF=SIGN(1.,EF+0.1)
 386         VF=AF-4.*EF*XW
 387         APF=SIGN(1.,EF+0.1)
 388         VPF=APF-4.*EF*XW
 389         GG=1.
 390         GZ=1./(8.*XW*(1.-XW))*SH*(SH-SQMZ)/((SH-SQMZ)**2+GZMZ**2)
 391         GZP=1./(8.*XW*(1.-XW))*SH*(SH-SQMZP)/((SH-SQMZP)**2+GZMZP**2)
 392         ZZ=1./(16.*XW*(1.-XW))**2*SH**2/((SH-SQMZ)**2+GZMZ**2)
 393         ZZP=2./(16.*XW*(1.-XW))**2*
 394      &  SH**2*((SH-SQMZ)*(SH-SQMZP)+GZMZ*GZMZP)/
 395      &  (((SH-SQMZ)**2+GZMZ**2)*((SH-SQMZP)**2+GZMZP**2))
 396         ZPZP=1./(16.*XW*(1.-XW))**2*SH**2/((SH-SQMZP)**2+GZMZP**2)
 397         IF(MSTP(44).EQ.1) THEN
 398 C...Only gamma* production included
 399           GZ=0.
 400           GZP=0.
 401           ZZ=0.
 402           ZZP=0.
 403           ZPZP=0.
 404         ELSEIF(MSTP(44).EQ.2) THEN
 405 C...Only Z0 production included
 406           GG=0.
 407           GZ=0.
 408           GZP=0.
 409           ZZP=0.
 410           ZPZP=0.
 411         ELSEIF(MSTP(44).EQ.3) THEN
 412 C...Only Z'0 production included
 413           GG=0.
 414           GZ=0.
 415           GZP=0.
 416           ZZ=0.
 417           ZZP=0.
 418         ELSEIF(MSTP(44).EQ.4) THEN
 419 C...Only gamma*/Z0 production included
 420           GZP=0.
 421           ZZP=0.
 422           ZPZP=0.
 423         ELSEIF(MSTP(44).EQ.5) THEN
 424 C...Only gamma*/Z'0 production included
 425           GZ=0.
 426           ZZ=0.
 427           ZZP=0.
 428         ELSEIF(MSTP(44).EQ.6) THEN
 429 C...Only Z0/Z'0 production included
 430           GG=0.
 431           GZ=0.
 432           GZP=0.
 433         ENDIF
 434         ASYM=2.*(EI*AI*GZ*EF*AF+EI*API*GZP*EF*APF+4.*VI*AI*ZZ*VF*AF+
 435      &  (VI*API+VPI*AI)*ZZP*(VF*APF+VPF*AF)+4.*VPI*API*ZPZP*VPF*APF)/
 436      &  (EI**2*GG*EF**2+EI*VI*GZ*EF*VF+EI*VPI*GZP*EF*VPF+
 437      &  (VI**2+AI**2)*ZZ*(VF**2+AF**2)+(VI*VPI+AI*API)*ZZP*
 438      &  (VF*VPF+AF*APF)+(VPI**2+API**2)*ZPZP*(VPF**2+APF**2))
 439         WT=1.+ASYM*CTHE(JT)+CTHE(JT)**2
 440         WTMAX=2.+ABS(ASYM)
 441
 442       ELSE
 443         WT=1.
 444         WTMAX=1.
 445       ENDIF
 446 C...Obtain correct angular distribution by rejection techniques.
 447       IF(WT.LT.RLU_HIJING(0)*WTMAX) GOTO 420
 448
 449 C...Construct massive four-vectors using angles chosen. Mark decayed
 450 C...resonances, add documentation lines. Shower evolution.
 451   500 DO 520 JT=1,JTMAX
 452       IF(KDCY(JT).EQ.0) GOTO 520
 453       ID=IREF(IP,JT)
 454       CALL LUDBRB_HIJING(NSD(JT)+1,NSD(JT)+2,ACOS(CTHE(JT)),PHI(JT),
 455      &DBLE(P(ID,1)/P(ID,4)),DBLE(P(ID,2)/P(ID,4)),DBLE(P(ID,3)/P(ID,4)))
 456       K(ID,1)=K(ID,1)+10
 457       K(ID,4)=NSD(JT)+1
 458       K(ID,5)=NSD(JT)+2
 459       IDOC=MINT(83)+MINT(4)
 460       DO 510 I=NSD(JT)+1,NSD(JT)+2
 461       MINT(4)=MINT(4)+1
 462       I1=MINT(83)+MINT(4)
 463       K(I,3)=I1
 464       K(I1,1)=21
 465       K(I1,2)=K(I,2)
 466       K(I1,3)=IREF(IP,JT+2)
 467       DO 510 J=1,5
 468   510 P(I1,J)=P(I,J)
 469       IF(JTMAX.EQ.1) THEN
 470         MINT(7)=MINT(83)+6+2*ISET(ISUB)
 471         MINT(8)=MINT(83)+7+2*ISET(ISUB)
 472       ENDIF
 473       IF(MSTP(71).GE.1.AND.KDCY(JT).EQ.1) CALL LUSHOW_HIJING(NSD(JT)+1,
 474      &NSD(JT)+2,P(ID,5))
 475
 476 C...Check if new resonances were produced, loop back if needed.
 477       IF(KDCY(JT).NE.3) GOTO 520
 478       NP=NP+1
 479       IREF(NP,1)=NSD(JT)+1
 480       IREF(NP,2)=NSD(JT)+2
 481       IREF(NP,3)=IDOC+1
 482       IREF(NP,4)=IDOC+2
 483       IREF(NP,5)=K(IREF(IP,JT),2)
 484       IREF(NP,6)=IREF(IP,JT)
 485   520 CONTINUE
 486   530 IF(IP.LT.NP) GOTO 100
 487
 488       RETURN
 489       END