HIJING/hipyset1_35/lushow_hijing.F

   1 * $Id$
   2
   3 C*********************************************************************
   4
   5       SUBROUTINE LUSHOW_HIJING(IP1,IP2,QMAX)
   6
   7 C...Purpose: to generate timelike parton showers from given partons.
   8       IMPLICIT DOUBLE PRECISION(D)
   9 #include "lujets_hijing.inc"
  10 #include "ludat1_hijing.inc"
  11 #include "ludat2_hijing.inc"
  12       DIMENSION PMTH(5,40),PS(5),PMA(4),PMSD(4),IEP(4),IPA(4),
  13      &KFLA(4),KFLD(4),KFL(4),ITRY(4),ISI(4),ISL(4),DP(4),DPT(5,4)
  14
  15 C...Initialization of cutoff masses etc.
  16       IF(MSTJ(41).LE.0.OR.(MSTJ(41).EQ.1.AND.QMAX.LE.PARJ(82)).OR.
  17      &QMAX.LE.MIN(PARJ(82),PARJ(83)).OR.MSTJ(41).GE.3) RETURN
  18       PMTH(1,21)=ULMASS_HIJING(21)
  19       PMTH(2,21)=SQRT(PMTH(1,21)**2+0.25*PARJ(82)**2)
  20       PMTH(3,21)=2.*PMTH(2,21)
  21       PMTH(4,21)=PMTH(3,21)
  22       PMTH(5,21)=PMTH(3,21)
  23       PMTH(1,22)=ULMASS_HIJING(22)
  24       PMTH(2,22)=SQRT(PMTH(1,22)**2+0.25*PARJ(83)**2)
  25       PMTH(3,22)=2.*PMTH(2,22)
  26       PMTH(4,22)=PMTH(3,22)
  27       PMTH(5,22)=PMTH(3,22)
  28       PMQTH1=PARJ(82)
  29       IF(MSTJ(41).EQ.2) PMQTH1=MIN(PARJ(82),PARJ(83))
  30       PMQTH2=PMTH(2,21)
  31       IF(MSTJ(41).EQ.2) PMQTH2=MIN(PMTH(2,21),PMTH(2,22))
  32       DO 100 IF=1,8
  33       PMTH(1,IF)=ULMASS_HIJING(IF)
  34       PMTH(2,IF)=SQRT(PMTH(1,IF)**2+0.25*PMQTH1**2)
  35       PMTH(3,IF)=PMTH(2,IF)+PMQTH2
  36       PMTH(4,IF)=SQRT(PMTH(1,IF)**2+0.25*PARJ(82)**2)+PMTH(2,21)
  37   100 PMTH(5,IF)=SQRT(PMTH(1,IF)**2+0.25*PARJ(83)**2)+PMTH(2,22)
  38       PT2MIN=MAX(0.5*PARJ(82),1.1*PARJ(81))**2
  39       ALAMS=PARJ(81)**2
  40       ALFM=LOG(PT2MIN/ALAMS)
  41
  42 C...Store positions of shower initiating partons.
  43       M3JC=0
  44       IF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.IP2.EQ.0) THEN
  45         NPA=1
  46         IPA(1)=IP1
  47       ELSEIF(MIN(IP1,IP2).GT.0.AND.MAX(IP1,IP2).LE.MIN(N,MSTU(4)-
  48      &MSTU(32))) THEN
  49         NPA=2
  50         IPA(1)=IP1
  51         IPA(2)=IP2
  52       ELSEIF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.IP2.LT.0.
  53      &AND.IP2.GE.-3) THEN
  54         NPA=IABS(IP2)
  55         DO 110 I=1,NPA
  56   110   IPA(I)=IP1+I-1
  57       ELSE
  58         CALL LUERRM_HIJING(12,
  59      &  '(LUSHOW_HIJING:) failed to reconstruct showering system')
  60         IF(MSTU(21).GE.1) RETURN
  61       ENDIF
  62
  63 C...Check on phase space available for emission.
  64       IREJ=0
  65       DO 120 J=1,5
  66   120 PS(J)=0.
  67       PM=0.
  68       DO 130 I=1,NPA
  69       KFLA(I)=IABS(K(IPA(I),2))
  70       PMA(I)=P(IPA(I),5)
  71       IF(KFLA(I).NE.0.AND.(KFLA(I).LE.8.OR.KFLA(I).EQ.21))
  72      &PMA(I)=PMTH(3,KFLA(I))
  73       PM=PM+PMA(I)
  74       IF(KFLA(I).EQ.0.OR.(KFLA(I).GT.8.AND.KFLA(I).NE.21).OR.
  75      &PMA(I).GT.QMAX) IREJ=IREJ+1
  76       DO 130 J=1,4
  77   130 PS(J)=PS(J)+P(IPA(I),J)
  78       IF(IREJ.EQ.NPA) RETURN
  79       PS(5)=SQRT(MAX(0.,PS(4)**2-PS(1)**2-PS(2)**2-PS(3)**2))
  80       IF(NPA.EQ.1) PS(5)=PS(4)
  81       IF(PS(5).LE.PM+PMQTH1) RETURN
  82       IF(NPA.EQ.2.AND.MSTJ(47).GE.1) THEN
  83         IF(KFLA(1).GE.1.AND.KFLA(1).LE.8.AND.KFLA(2).GE.1.AND.
  84      &  KFLA(2).LE.8) M3JC=1
  85         IF(MSTJ(47).GE.2) M3JC=1
  86       ENDIF
  87
  88 C...Define imagined single initiator of shower for parton system.
  89       NS=N
  90       IF(N.GT.MSTU(4)-MSTU(32)-5) THEN
  91          CALL LUERRM_HIJING(11
  92      $        ,'(LUSHOW_HIJING:) no more memory left in LUJETS_HIJING')
  93         IF(MSTU(21).GE.1) RETURN
  94       ENDIF
  95       IF(NPA.GE.2) THEN
  96         K(N+1,1)=11
  97         K(N+1,2)=21
  98         K(N+1,3)=0
  99         K(N+1,4)=0
 100         K(N+1,5)=0
 101         P(N+1,1)=0.
 102         P(N+1,2)=0.
 103         P(N+1,3)=0.
 104         P(N+1,4)=PS(5)
 105         P(N+1,5)=PS(5)
 106         V(N+1,5)=PS(5)**2
 107         N=N+1
 108       ENDIF
 109
 110 C...Loop over partons that may branch.
 111       NEP=NPA
 112       IM=NS
 113       IF(NPA.EQ.1) IM=NS-1
 114   140 IM=IM+1
 115       IF(N.GT.NS) THEN
 116         IF(IM.GT.N) GOTO 380
 117         KFLM=IABS(K(IM,2))
 118         IF(KFLM.EQ.0.OR.(KFLM.GT.8.AND.KFLM.NE.21)) GOTO 140
 119         IF(P(IM,5).LT.PMTH(2,KFLM)) GOTO 140
 120         IGM=K(IM,3)
 121       ELSE
 122         IGM=-1
 123       ENDIF
 124       IF(N+NEP.GT.MSTU(4)-MSTU(32)-5) THEN
 125          CALL LUERRM_HIJING(11
 126      $        ,'(LUSHOW_HIJING:) no more memory left in LUJETS_HIJING')
 127         IF(MSTU(21).GE.1) RETURN
 128       ENDIF
 129
 130 C...Position of aunt (sister to branching parton).
 131 C...Origin and flavour of daughters.
 132       IAU=0
 133       IF(IGM.GT.0) THEN
 134         IF(K(IM-1,3).EQ.IGM) IAU=IM-1
 135         IF(N.GE.IM+1.AND.K(IM+1,3).EQ.IGM) IAU=IM+1
 136       ENDIF
 137       IF(IGM.GE.0) THEN
 138         K(IM,4)=N+1
 139         DO 150 I=1,NEP
 140   150   K(N+I,3)=IM
 141       ELSE
 142         K(N+1,3)=IPA(1)
 143       ENDIF
 144       IF(IGM.LE.0) THEN
 145         DO 160 I=1,NEP
 146   160   K(N+I,2)=K(IPA(I),2)
 147       ELSEIF(KFLM.NE.21) THEN
 148         K(N+1,2)=K(IM,2)
 149         K(N+2,2)=K(IM,5)
 150       ELSEIF(K(IM,5).EQ.21) THEN
 151         K(N+1,2)=21
 152         K(N+2,2)=21
 153       ELSE
 154         K(N+1,2)=K(IM,5)
 155         K(N+2,2)=-K(IM,5)
 156       ENDIF
 157
 158 C...Reset flags on daughers and tries made.
 159       DO 170 IP=1,NEP
 160       K(N+IP,1)=3
 161       K(N+IP,4)=0
 162       K(N+IP,5)=0
 163       KFLD(IP)=IABS(K(N+IP,2))
 164       ITRY(IP)=0
 165       ISL(IP)=0
 166       ISI(IP)=0
 167   170 IF(KFLD(IP).GT.0.AND.(KFLD(IP).LE.8.OR.KFLD(IP).EQ.21)) ISI(IP)=1
 168       ISLM=0
 169
 170 C...Maximum virtuality of daughters.
 171       IF(IGM.LE.0) THEN
 172         DO 180 I=1,NPA
 173         IF(NPA.GE.3) P(N+I,4)=(PS(4)*P(IPA(I),4)-PS(1)*P(IPA(I),1)-
 174      &  PS(2)*P(IPA(I),2)-PS(3)*P(IPA(I),3))/PS(5)
 175         P(N+I,5)=MIN(QMAX,PS(5))
 176         IF(NPA.GE.3) P(N+I,5)=MIN(P(N+I,5),P(N+I,4))
 177   180   IF(ISI(I).EQ.0) P(N+I,5)=P(IPA(I),5)
 178       ELSE
 179         IF(MSTJ(43).LE.2) PEM=V(IM,2)
 180         IF(MSTJ(43).GE.3) PEM=P(IM,4)
 181         P(N+1,5)=MIN(P(IM,5),V(IM,1)*PEM)
 182         P(N+2,5)=MIN(P(IM,5),(1.-V(IM,1))*PEM)
 183         IF(K(N+2,2).EQ.22) P(N+2,5)=PMTH(1,22)
 184       ENDIF
 185       DO 190 I=1,NEP
 186       PMSD(I)=P(N+I,5)
 187       IF(ISI(I).EQ.1) THEN
 188         IF(P(N+I,5).LE.PMTH(3,KFLD(I))) P(N+I,5)=PMTH(1,KFLD(I))
 189       ENDIF
 190   190 V(N+I,5)=P(N+I,5)**2
 191
 192 C...Choose one of the daughters for evolution.
 193   200 INUM=0
 194       IF(NEP.EQ.1) INUM=1
 195       DO 210 I=1,NEP
 196   210 IF(INUM.EQ.0.AND.ISL(I).EQ.1) INUM=I
 197       DO 220 I=1,NEP
 198       IF(INUM.EQ.0.AND.ITRY(I).EQ.0.AND.ISI(I).EQ.1) THEN
 199         IF(P(N+I,5).GE.PMTH(2,KFLD(I))) INUM=I
 200       ENDIF
 201   220 CONTINUE
 202       IF(INUM.EQ.0) THEN
 203         RMAX=0.
 204         DO 230 I=1,NEP
 205         IF(ISI(I).EQ.1.AND.PMSD(I).GE.PMQTH2) THEN
 206           RPM=P(N+I,5)/PMSD(I)
 207           IF(RPM.GT.RMAX.AND.P(N+I,5).GE.PMTH(2,KFLD(I))) THEN
 208             RMAX=RPM
 209             INUM=I
 210           ENDIF
 211         ENDIF
 212   230   CONTINUE
 213       ENDIF
 214
 215 C...Store information on choice of evolving daughter.
 216       INUM=MAX(1,INUM)
 217       IEP(1)=N+INUM
 218       DO 240 I=2,NEP
 219       IEP(I)=IEP(I-1)+1
 220   240 IF(IEP(I).GT.N+NEP) IEP(I)=N+1
 221       DO 250 I=1,NEP
 222   250 KFL(I)=IABS(K(IEP(I),2))
 223       ITRY(INUM)=ITRY(INUM)+1
 224       IF(ITRY(INUM).GT.200) THEN
 225          CALL LUERRM_HIJING(14
 226      $        ,'(LUSHOW_HIJING:) caught in infinite loop')
 227         IF(MSTU(21).GE.1) RETURN
 228       ENDIF
 229       Z=0.5
 230       IF(KFL(1).EQ.0.OR.(KFL(1).GT.8.AND.KFL(1).NE.21)) GOTO 300
 231       IF(P(IEP(1),5).LT.PMTH(2,KFL(1))) GOTO 300
 232
 233 C...Calculate allowed z range.
 234       IF(NEP.EQ.1) THEN
 235         PMED=PS(4)
 236       ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN
 237         PMED=P(IM,5)
 238       ELSE
 239         IF(INUM.EQ.1) PMED=V(IM,1)*PEM
 240         IF(INUM.EQ.2) PMED=(1.-V(IM,1))*PEM
 241       ENDIF
 242       IF(MOD(MSTJ(43),2).EQ.1) THEN
 243         ZC=PMTH(2,21)/PMED
 244         ZCE=PMTH(2,22)/PMED
 245       ELSE
 246         ZC=0.5*(1.-SQRT(MAX(0.,1.-(2.*PMTH(2,21)/PMED)**2)))
 247         IF(ZC.LT.1E-4) ZC=(PMTH(2,21)/PMED)**2
 248         ZCE=0.5*(1.-SQRT(MAX(0.,1.-(2.*PMTH(2,22)/PMED)**2)))
 249         IF(ZCE.LT.1E-4) ZCE=(PMTH(2,22)/PMED)**2
 250       ENDIF
 251       ZC=MIN(ZC,0.491)
 252       ZCE=MIN(ZCE,0.491)
 253       IF((MSTJ(41).EQ.1.AND.ZC.GT.0.49).OR.(MSTJ(41).EQ.2.AND.
 254      &MIN(ZC,ZCE).GT.0.49)) THEN
 255         P(IEP(1),5)=PMTH(1,KFL(1))
 256         V(IEP(1),5)=P(IEP(1),5)**2
 257         GOTO 300
 258       ENDIF
 259
 260 C...Integral of Altarelli-Parisi z kernel for QCD.
 261       IF(MSTJ(49).EQ.0.AND.KFL(1).EQ.21) THEN
 262         FBR=6.*LOG((1.-ZC)/ZC)+MSTJ(45)*(0.5-ZC)
 263       ELSEIF(MSTJ(49).EQ.0) THEN
 264         FBR=(8./3.)*LOG((1.-ZC)/ZC)
 265
 266 C...Integral of Altarelli-Parisi z kernel for scalar gluon.
 267       ELSEIF(MSTJ(49).EQ.1.AND.KFL(1).EQ.21) THEN
 268         FBR=(PARJ(87)+MSTJ(45)*PARJ(88))*(1.-2.*ZC)
 269       ELSEIF(MSTJ(49).EQ.1) THEN
 270         FBR=(1.-2.*ZC)/3.
 271         IF(IGM.EQ.0.AND.M3JC.EQ.1) FBR=4.*FBR
 272
 273 C...Integral of Altarelli-Parisi z kernel for Abelian vector gluon.
 274       ELSEIF(KFL(1).EQ.21) THEN
 275         FBR=6.*MSTJ(45)*(0.5-ZC)
 276       ELSE
 277         FBR=2.*LOG((1.-ZC)/ZC)
 278       ENDIF
 279
 280 C...Integral of Altarelli-Parisi kernel for photon emission.
 281       IF(MSTJ(41).EQ.2.AND.KFL(1).GE.1.AND.KFL(1).LE.8)
 282      &FBRE=(KCHG(KFL(1),1)/3.)**2*2.*LOG((1.-ZCE)/ZCE)
 283
 284 C...Inner veto algorithm starts. Find maximum mass for evolution.
 285   260 PMS=V(IEP(1),5)
 286       IF(IGM.GE.0) THEN
 287         PM2=0.
 288         DO 270 I=2,NEP
 289         PM=P(IEP(I),5)
 290         IF(KFL(I).GT.0.AND.(KFL(I).LE.8.OR.KFL(I).EQ.21)) PM=
 291      &  PMTH(2,KFL(I))
 292   270   PM2=PM2+PM
 293         PMS=MIN(PMS,(P(IM,5)-PM2)**2)
 294       ENDIF
 295
 296 C...Select mass for daughter in QCD evolution.
 297       B0=27./6.
 298       DO 280 IF=4,MSTJ(45)
 299   280 IF(PMS.GT.4.*PMTH(2,IF)**2) B0=(33.-2.*IF)/6.
 300       IF(MSTJ(44).LE.0) THEN
 301          PMSQCD=PMS*EXP(MAX(-100.,LOG(RLU_HIJING(0))*PARU(2)/(PARU(111)
 302      $        *FBR)))
 303       ELSEIF(MSTJ(44).EQ.1) THEN
 304         PMSQCD=4.*ALAMS*(0.25*PMS/ALAMS)**(RLU_HIJING(0)**(B0/FBR))
 305       ELSE
 306         PMSQCD=PMS*RLU_HIJING(0)**(ALFM*B0/FBR)
 307       ENDIF
 308       IF(ZC.GT.0.49.OR.PMSQCD.LE.PMTH(4,KFL(1))**2) PMSQCD=
 309      &PMTH(2,KFL(1))**2
 310       V(IEP(1),5)=PMSQCD
 311       MCE=1
 312
 313 C...Select mass for daughter in QED evolution.
 314       IF(MSTJ(41).EQ.2.AND.KFL(1).GE.1.AND.KFL(1).LE.8) THEN
 315          PMSQED=PMS*EXP(MAX(-100.,LOG(RLU_HIJING(0))*PARU(2)/(PARU(101)
 316      $        *FBRE)))
 317         IF(ZCE.GT.0.49.OR.PMSQED.LE.PMTH(5,KFL(1))**2) PMSQED=
 318      &  PMTH(2,KFL(1))**2
 319         IF(PMSQED.GT.PMSQCD) THEN
 320           V(IEP(1),5)=PMSQED
 321           MCE=2
 322         ENDIF
 323       ENDIF
 324
 325 C...Check whether daughter mass below cutoff.
 326       P(IEP(1),5)=SQRT(V(IEP(1),5))
 327       IF(P(IEP(1),5).LE.PMTH(3,KFL(1))) THEN
 328         P(IEP(1),5)=PMTH(1,KFL(1))
 329         V(IEP(1),5)=P(IEP(1),5)**2
 330         GOTO 300
 331       ENDIF
 332
 333 C...Select z value of branching: q -> qgamma.
 334       IF(MCE.EQ.2) THEN
 335         Z=1.-(1.-ZCE)*(ZCE/(1.-ZCE))**RLU_HIJING(0)
 336         IF(1.+Z**2.LT.2.*RLU_HIJING(0)) GOTO 260
 337         K(IEP(1),5)=22
 338
 339 C...Select z value of branching: q -> qg, g -> gg, g -> qqbar.
 340       ELSEIF(MSTJ(49).NE.1.AND.KFL(1).NE.21) THEN
 341         Z=1.-(1.-ZC)*(ZC/(1.-ZC))**RLU_HIJING(0)
 342         IF(1.+Z**2.LT.2.*RLU_HIJING(0)) GOTO 260
 343         K(IEP(1),5)=21
 344       ELSEIF(MSTJ(49).EQ.0.AND.MSTJ(45)*(0.5-ZC).LT.RLU_HIJING(0)*FBR)
 345      $       THEN
 346         Z=(1.-ZC)*(ZC/(1.-ZC))**RLU_HIJING(0)
 347         IF(RLU_HIJING(0).GT.0.5) Z=1.-Z
 348         IF((1.-Z*(1.-Z))**2.LT.RLU_HIJING(0)) GOTO 260
 349         K(IEP(1),5)=21
 350       ELSEIF(MSTJ(49).NE.1) THEN
 351         Z=ZC+(1.-2.*ZC)*RLU_HIJING(0)
 352         IF(Z**2+(1.-Z)**2.LT.RLU_HIJING(0)) GOTO 260
 353         KFLB=1+INT(MSTJ(45)*RLU_HIJING(0))
 354         PMQ=4.*PMTH(2,KFLB)**2/V(IEP(1),5)
 355         IF(PMQ.GE.1.) GOTO 260
 356         PMQ0=4.*PMTH(2,21)**2/V(IEP(1),5)
 357         IF(MOD(MSTJ(43),2).EQ.0.AND.(1.+0.5*PMQ)*SQRT(1.-PMQ).LT.
 358      &  RLU_HIJING(0)*(1.+0.5*PMQ0)*SQRT(1.-PMQ0)) GOTO 260
 359         K(IEP(1),5)=KFLB
 360
 361 C...Ditto for scalar gluon model.
 362       ELSEIF(KFL(1).NE.21) THEN
 363         Z=1.-SQRT(ZC**2+RLU_HIJING(0)*(1.-2.*ZC))
 364         K(IEP(1),5)=21
 365       ELSEIF(RLU_HIJING(0)*(PARJ(87)+MSTJ(45)*PARJ(88)).LE.PARJ(87))
 366      $       THEN
 367         Z=ZC+(1.-2.*ZC)*RLU_HIJING(0)
 368         K(IEP(1),5)=21
 369       ELSE
 370         Z=ZC+(1.-2.*ZC)*RLU_HIJING(0)
 371         KFLB=1+INT(MSTJ(45)*RLU_HIJING(0))
 372         PMQ=4.*PMTH(2,KFLB)**2/V(IEP(1),5)
 373         IF(PMQ.GE.1.) GOTO 260
 374         K(IEP(1),5)=KFLB
 375       ENDIF
 376       IF(MCE.EQ.1.AND.MSTJ(44).GE.2) THEN
 377         IF(Z*(1.-Z)*V(IEP(1),5).LT.PT2MIN) GOTO 260
 378         IF(ALFM/LOG(V(IEP(1),5)*Z*(1.-Z)/ALAMS).LT.RLU_HIJING(0)) GOTO
 379      $       260
 380       ENDIF
 381
 382 C...Check if z consistent with chosen m.
 383       IF(KFL(1).EQ.21) THEN
 384         KFLGD1=IABS(K(IEP(1),5))
 385         KFLGD2=KFLGD1
 386       ELSE
 387         KFLGD1=KFL(1)
 388         KFLGD2=IABS(K(IEP(1),5))
 389       ENDIF
 390       IF(NEP.EQ.1) THEN
 391         PED=PS(4)
 392       ELSEIF(NEP.GE.3) THEN
 393         PED=P(IEP(1),4)
 394       ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN
 395         PED=0.5*(V(IM,5)+V(IEP(1),5)-PM2**2)/P(IM,5)
 396       ELSE
 397         IF(IEP(1).EQ.N+1) PED=V(IM,1)*PEM
 398         IF(IEP(1).EQ.N+2) PED=(1.-V(IM,1))*PEM
 399       ENDIF
 400       IF(MOD(MSTJ(43),2).EQ.1) THEN
 401         PMQTH3=0.5*PARJ(82)
 402         IF(KFLGD2.EQ.22) PMQTH3=0.5*PARJ(83)
 403         PMQ1=(PMTH(1,KFLGD1)**2+PMQTH3**2)/V(IEP(1),5)
 404         PMQ2=(PMTH(1,KFLGD2)**2+PMQTH3**2)/V(IEP(1),5)
 405         ZD=SQRT(MAX(0.,(1.-V(IEP(1),5)/PED**2)*((1.-PMQ1-PMQ2)**2-
 406      &  4.*PMQ1*PMQ2)))
 407         ZH=1.+PMQ1-PMQ2
 408       ELSE
 409         ZD=SQRT(MAX(0.,1.-V(IEP(1),5)/PED**2))
 410         ZH=1.
 411       ENDIF
 412       ZL=0.5*(ZH-ZD)
 413       ZU=0.5*(ZH+ZD)
 414       IF(Z.LT.ZL.OR.Z.GT.ZU) GOTO 260
 415       IF(KFL(1).EQ.21) V(IEP(1),3)=LOG(ZU*(1.-ZL)/MAX(1E-20,ZL*
 416      &(1.-ZU)))
 417       IF(KFL(1).NE.21) V(IEP(1),3)=LOG((1.-ZL)/MAX(1E-10,1.-ZU))
 418
 419 C...Three-jet matrix element correction.
 420       IF(IGM.EQ.0.AND.M3JC.EQ.1) THEN
 421         X1=Z*(1.+V(IEP(1),5)/V(NS+1,5))
 422         X2=1.-V(IEP(1),5)/V(NS+1,5)
 423         X3=(1.-X1)+(1.-X2)
 424         IF(MCE.EQ.2) THEN
 425           KI1=K(IPA(INUM),2)
 426           KI2=K(IPA(3-INUM),2)
 427           QF1=KCHG(IABS(KI1),1)*ISIGN(1,KI1)/3.
 428           QF2=KCHG(IABS(KI2),1)*ISIGN(1,KI2)/3.
 429           WSHOW=QF1**2*(1.-X1)/X3*(1.+(X1/(2.-X2))**2)+
 430      &    QF2**2*(1.-X2)/X3*(1.+(X2/(2.-X1))**2)
 431           WME=(QF1*(1.-X1)/X3-QF2*(1.-X2)/X3)**2*(X1**2+X2**2)
 432         ELSEIF(MSTJ(49).NE.1) THEN
 433           WSHOW=1.+(1.-X1)/X3*(X1/(2.-X2))**2+
 434      &    (1.-X2)/X3*(X2/(2.-X1))**2
 435           WME=X1**2+X2**2
 436         ELSE
 437           WSHOW=4.*X3*((1.-X1)/(2.-X2)**2+(1.-X2)/(2.-X1)**2)
 438           WME=X3**2
 439         ENDIF
 440         IF(WME.LT.RLU_HIJING(0)*WSHOW) GOTO 260
 441
 442 C...Impose angular ordering by rejection of nonordered emission.
 443       ELSEIF(MCE.EQ.1.AND.IGM.GT.0.AND.MSTJ(42).GE.2) THEN
 444         MAOM=1
 445         ZM=V(IM,1)
 446         IF(IEP(1).EQ.N+2) ZM=1.-V(IM,1)
 447         THE2ID=Z*(1.-Z)*(ZM*P(IM,4))**2/V(IEP(1),5)
 448         IAOM=IM
 449   290   IF(K(IAOM,5).EQ.22) THEN
 450           IAOM=K(IAOM,3)
 451           IF(K(IAOM,3).LE.NS) MAOM=0
 452           IF(MAOM.EQ.1) GOTO 290
 453         ENDIF
 454         IF(MAOM.EQ.1) THEN
 455           THE2IM=V(IAOM,1)*(1.-V(IAOM,1))*P(IAOM,4)**2/V(IAOM,5)
 456           IF(THE2ID.LT.THE2IM) GOTO 260
 457         ENDIF
 458       ENDIF
 459
 460 C...Impose user-defined maximum angle at first branching.
 461       IF(MSTJ(48).EQ.1) THEN
 462         IF(NEP.EQ.1.AND.IM.EQ.NS) THEN
 463           THE2ID=Z*(1.-Z)*PS(4)**2/V(IEP(1),5)
 464           IF(THE2ID.LT.1./PARJ(85)**2) GOTO 260
 465         ELSEIF(NEP.EQ.2.AND.IEP(1).EQ.NS+2) THEN
 466           THE2ID=Z*(1.-Z)*(0.5*P(IM,4))**2/V(IEP(1),5)
 467           IF(THE2ID.LT.1./PARJ(85)**2) GOTO 260
 468         ELSEIF(NEP.EQ.2.AND.IEP(1).EQ.NS+3) THEN
 469           THE2ID=Z*(1.-Z)*(0.5*P(IM,4))**2/V(IEP(1),5)
 470           IF(THE2ID.LT.1./PARJ(86)**2) GOTO 260
 471         ENDIF
 472       ENDIF
 473
 474 C...End of inner veto algorithm. Check if only one leg evolved so far.
 475   300 V(IEP(1),1)=Z
 476       ISL(1)=0
 477       ISL(2)=0
 478       IF(NEP.EQ.1) GOTO 330
 479       IF(NEP.EQ.2.AND.P(IEP(1),5)+P(IEP(2),5).GE.P(IM,5)) GOTO 200
 480       DO 310 I=1,NEP
 481       IF(ITRY(I).EQ.0.AND.KFLD(I).GT.0.AND.(KFLD(I).LE.8.OR.KFLD(I).EQ.
 482      &21)) THEN
 483         IF(P(N+I,5).GE.PMTH(2,KFLD(I))) GOTO 200
 484       ENDIF
 485   310 CONTINUE
 486
 487 C...Check if chosen multiplet m1,m2,z1,z2 is physical.
 488       IF(NEP.EQ.3) THEN
 489         PA1S=(P(N+1,4)+P(N+1,5))*(P(N+1,4)-P(N+1,5))
 490         PA2S=(P(N+2,4)+P(N+2,5))*(P(N+2,4)-P(N+2,5))
 491         PA3S=(P(N+3,4)+P(N+3,5))*(P(N+3,4)-P(N+3,5))
 492         PTS=0.25*(2.*PA1S*PA2S+2.*PA1S*PA3S+2.*PA2S*PA3S-
 493      &  PA1S**2-PA2S**2-PA3S**2)/PA1S
 494         IF(PTS.LE.0.) GOTO 200
 495       ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2.OR.MOD(MSTJ(43),2).EQ.0) THEN
 496         DO 320 I1=N+1,N+2
 497         KFLDA=IABS(K(I1,2))
 498         IF(KFLDA.EQ.0.OR.(KFLDA.GT.8.AND.KFLDA.NE.21)) GOTO 320
 499         IF(P(I1,5).LT.PMTH(2,KFLDA)) GOTO 320
 500         IF(KFLDA.EQ.21) THEN
 501           KFLGD1=IABS(K(I1,5))
 502           KFLGD2=KFLGD1
 503         ELSE
 504           KFLGD1=KFLDA
 505           KFLGD2=IABS(K(I1,5))
 506         ENDIF
 507         I2=2*N+3-I1
 508         IF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN
 509           PED=0.5*(V(IM,5)+V(I1,5)-V(I2,5))/P(IM,5)
 510         ELSE
 511           IF(I1.EQ.N+1) ZM=V(IM,1)
 512           IF(I1.EQ.N+2) ZM=1.-V(IM,1)
 513           PML=SQRT((V(IM,5)-V(N+1,5)-V(N+2,5))**2-
 514      &    4.*V(N+1,5)*V(N+2,5))
 515           PED=PEM*(0.5*(V(IM,5)-PML+V(I1,5)-V(I2,5))+PML*ZM)/V(IM,5)
 516         ENDIF
 517         IF(MOD(MSTJ(43),2).EQ.1) THEN
 518           PMQTH3=0.5*PARJ(82)
 519           IF(KFLGD2.EQ.22) PMQTH3=0.5*PARJ(83)
 520           PMQ1=(PMTH(1,KFLGD1)**2+PMQTH3**2)/V(I1,5)
 521           PMQ2=(PMTH(1,KFLGD2)**2+PMQTH3**2)/V(I1,5)
 522           ZD=SQRT(MAX(0.,(1.-V(I1,5)/PED**2)*((1.-PMQ1-PMQ2)**2-
 523      &    4.*PMQ1*PMQ2)))
 524           ZH=1.+PMQ1-PMQ2
 525         ELSE
 526           ZD=SQRT(MAX(0.,1.-V(I1,5)/PED**2))
 527           ZH=1.
 528         ENDIF
 529         ZL=0.5*(ZH-ZD)
 530         ZU=0.5*(ZH+ZD)
 531         IF(I1.EQ.N+1.AND.(V(I1,1).LT.ZL.OR.V(I1,1).GT.ZU)) ISL(1)=1
 532         IF(I1.EQ.N+2.AND.(V(I1,1).LT.ZL.OR.V(I1,1).GT.ZU)) ISL(2)=1
 533         IF(KFLDA.EQ.21) V(I1,4)=LOG(ZU*(1.-ZL)/MAX(1E-20,ZL*(1.-ZU)))
 534         IF(KFLDA.NE.21) V(I1,4)=LOG((1.-ZL)/MAX(1E-10,1.-ZU))
 535   320   CONTINUE
 536         IF(ISL(1).EQ.1.AND.ISL(2).EQ.1.AND.ISLM.NE.0) THEN
 537           ISL(3-ISLM)=0
 538           ISLM=3-ISLM
 539         ELSEIF(ISL(1).EQ.1.AND.ISL(2).EQ.1) THEN
 540           ZDR1=MAX(0D0,V(N+1,3)/MAX(1D-6,V(N+1,4))-1D0)
 541           ZDR2=MAX(0D0,V(N+2,3)/MAX(1D-6,V(N+2,4))-1D0)
 542           IF(ZDR2.GT.RLU_HIJING(0)*(ZDR1+ZDR2)) ISL(1)=0
 543           IF(ISL(1).EQ.1) ISL(2)=0
 544           IF(ISL(1).EQ.0) ISLM=1
 545           IF(ISL(2).EQ.0) ISLM=2
 546         ENDIF
 547         IF(ISL(1).EQ.1.OR.ISL(2).EQ.1) GOTO 200
 548       ENDIF
 549       IF(IGM.GT.0.AND.MOD(MSTJ(43),2).EQ.1.AND.(P(N+1,5).GE.
 550      &PMTH(2,KFLD(1)).OR.P(N+2,5).GE.PMTH(2,KFLD(2)))) THEN
 551         PMQ1=V(N+1,5)/V(IM,5)
 552         PMQ2=V(N+2,5)/V(IM,5)
 553         ZD=SQRT(MAX(0.,(1.-V(IM,5)/PEM**2)*((1.-PMQ1-PMQ2)**2-
 554      &  4.*PMQ1*PMQ2)))
 555         ZH=1.+PMQ1-PMQ2
 556         ZL=0.5*(ZH-ZD)
 557         ZU=0.5*(ZH+ZD)
 558         IF(V(IM,1).LT.ZL.OR.V(IM,1).GT.ZU) GOTO 200
 559       ENDIF
 560
 561 C...Accepted branch. Construct four-momentum for initial partons.
 562   330 MAZIP=0
 563       MAZIC=0
 564       IF(NEP.EQ.1) THEN
 565         P(N+1,1)=0.
 566         P(N+1,2)=0.
 567         P(N+1,3)=SQRT(MAX(0.,(P(IPA(1),4)+P(N+1,5))*(P(IPA(1),4)-
 568      &  P(N+1,5))))
 569         P(N+1,4)=P(IPA(1),4)
 570         V(N+1,2)=P(N+1,4)
 571       ELSEIF(IGM.EQ.0.AND.NEP.EQ.2) THEN
 572         PED1=0.5*(V(IM,5)+V(N+1,5)-V(N+2,5))/P(IM,5)
 573         P(N+1,1)=0.
 574         P(N+1,2)=0.
 575         P(N+1,3)=SQRT(MAX(0.,(PED1+P(N+1,5))*(PED1-P(N+1,5))))
 576         P(N+1,4)=PED1
 577         P(N+2,1)=0.
 578         P(N+2,2)=0.
 579         P(N+2,3)=-P(N+1,3)
 580         P(N+2,4)=P(IM,5)-PED1
 581         V(N+1,2)=P(N+1,4)
 582         V(N+2,2)=P(N+2,4)
 583       ELSEIF(NEP.EQ.3) THEN
 584         P(N+1,1)=0.
 585         P(N+1,2)=0.
 586         P(N+1,3)=SQRT(MAX(0.,PA1S))
 587         P(N+2,1)=SQRT(PTS)
 588         P(N+2,2)=0.
 589         P(N+2,3)=0.5*(PA3S-PA2S-PA1S)/P(N+1,3)
 590         P(N+3,1)=-P(N+2,1)
 591         P(N+3,2)=0.
 592         P(N+3,3)=-(P(N+1,3)+P(N+2,3))
 593         V(N+1,2)=P(N+1,4)
 594         V(N+2,2)=P(N+2,4)
 595         V(N+3,2)=P(N+3,4)
 596
 597 C...Construct transverse momentum for ordinary branching in shower.
 598       ELSE
 599         ZM=V(IM,1)
 600         PZM=SQRT(MAX(0.,(PEM+P(IM,5))*(PEM-P(IM,5))))
 601         PMLS=(V(IM,5)-V(N+1,5)-V(N+2,5))**2-4.*V(N+1,5)*V(N+2,5)
 602         IF(PZM.LE.0.) THEN
 603           PTS=0.
 604         ELSEIF(MOD(MSTJ(43),2).EQ.1) THEN
 605           PTS=(PEM**2*(ZM*(1.-ZM)*V(IM,5)-(1.-ZM)*V(N+1,5)-
 606      &    ZM*V(N+2,5))-0.25*PMLS)/PZM**2
 607         ELSE
 608           PTS=PMLS*(ZM*(1.-ZM)*PEM**2/V(IM,5)-0.25)/PZM**2
 609         ENDIF
 610         PT=SQRT(MAX(0.,PTS))
 611
 612 C...Find coefficient of azimuthal asymmetry due to gluon polarization.
 613         HAZIP=0.
 614         IF(MSTJ(49).NE.1.AND.MOD(MSTJ(46),2).EQ.1.AND.K(IM,2).EQ.21.
 615      &  AND.IAU.NE.0) THEN
 616           IF(K(IGM,3).NE.0) MAZIP=1
 617           ZAU=V(IGM,1)
 618           IF(IAU.EQ.IM+1) ZAU=1.-V(IGM,1)
 619           IF(MAZIP.EQ.0) ZAU=0.
 620           IF(K(IGM,2).NE.21) THEN
 621             HAZIP=2.*ZAU/(1.+ZAU**2)
 622           ELSE
 623             HAZIP=(ZAU/(1.-ZAU*(1.-ZAU)))**2
 624           ENDIF
 625           IF(K(N+1,2).NE.21) THEN
 626             HAZIP=HAZIP*(-2.*ZM*(1.-ZM))/(1.-2.*ZM*(1.-ZM))
 627           ELSE
 628             HAZIP=HAZIP*(ZM*(1.-ZM)/(1.-ZM*(1.-ZM)))**2
 629           ENDIF
 630         ENDIF
 631
 632 C...Find coefficient of azimuthal asymmetry due to soft gluon
 633 C...interference.
 634         HAZIC=0.
 635         IF(MSTJ(46).GE.2.AND.(K(N+1,2).EQ.21.OR.K(N+2,2).EQ.21).
 636      &  AND.IAU.NE.0) THEN
 637           IF(K(IGM,3).NE.0) MAZIC=N+1
 638           IF(K(IGM,3).NE.0.AND.K(N+1,2).NE.21) MAZIC=N+2
 639           IF(K(IGM,3).NE.0.AND.K(N+1,2).EQ.21.AND.K(N+2,2).EQ.21.AND.
 640      &    ZM.GT.0.5) MAZIC=N+2
 641           IF(K(IAU,2).EQ.22) MAZIC=0
 642           ZS=ZM
 643           IF(MAZIC.EQ.N+2) ZS=1.-ZM
 644           ZGM=V(IGM,1)
 645           IF(IAU.EQ.IM-1) ZGM=1.-V(IGM,1)
 646           IF(MAZIC.EQ.0) ZGM=1.
 647           HAZIC=(P(IM,5)/P(IGM,5))*SQRT((1.-ZS)*(1.-ZGM)/(ZS*ZGM))
 648           HAZIC=MIN(0.95,HAZIC)
 649         ENDIF
 650       ENDIF
 651
 652 C...Construct kinematics for ordinary branching in shower.
 653   340 IF(NEP.EQ.2.AND.IGM.GT.0) THEN
 654         IF(MOD(MSTJ(43),2).EQ.1) THEN
 655           P(N+1,4)=PEM*V(IM,1)
 656         ELSE
 657           P(N+1,4)=PEM*(0.5*(V(IM,5)-SQRT(PMLS)+V(N+1,5)-V(N+2,5))+
 658      &    SQRT(PMLS)*ZM)/V(IM,5)
 659         ENDIF
 660         PHI=PARU(2)*RLU_HIJING(0)
 661         P(N+1,1)=PT*COS(PHI)
 662         P(N+1,2)=PT*SIN(PHI)
 663         IF(PZM.GT.0.) THEN
 664           P(N+1,3)=0.5*(V(N+2,5)-V(N+1,5)-V(IM,5)+2.*PEM*P(N+1,4))/PZM
 665         ELSE
 666           P(N+1,3)=0.
 667         ENDIF
 668         P(N+2,1)=-P(N+1,1)
 669         P(N+2,2)=-P(N+1,2)
 670         P(N+2,3)=PZM-P(N+1,3)
 671         P(N+2,4)=PEM-P(N+1,4)
 672         IF(MSTJ(43).LE.2) THEN
 673           V(N+1,2)=(PEM*P(N+1,4)-PZM*P(N+1,3))/P(IM,5)
 674           V(N+2,2)=(PEM*P(N+2,4)-PZM*P(N+2,3))/P(IM,5)
 675         ENDIF
 676       ENDIF
 677
 678 C...Rotate and boost daughters.
 679       IF(IGM.GT.0) THEN
 680         IF(MSTJ(43).LE.2) THEN
 681           BEX=P(IGM,1)/P(IGM,4)
 682           BEY=P(IGM,2)/P(IGM,4)
 683           BEZ=P(IGM,3)/P(IGM,4)
 684           GA=P(IGM,4)/P(IGM,5)
 685           GABEP=GA*(GA*(BEX*P(IM,1)+BEY*P(IM,2)+BEZ*P(IM,3))/(1.+GA)-
 686      &    P(IM,4))
 687         ELSE
 688           BEX=0.
 689           BEY=0.
 690           BEZ=0.
 691           GA=1.
 692           GABEP=0.
 693         ENDIF
 694         THE=ULANGL_HIJING(P(IM,3)+GABEP*BEZ,SQRT((P(IM,1)+GABEP*BEX)**2+
 695      &  (P(IM,2)+GABEP*BEY)**2))
 696         PHI=ULANGL_HIJING(P(IM,1)+GABEP*BEX,P(IM,2)+GABEP*BEY)
 697         DO 350 I=N+1,N+2
 698         DP(1)=COS(THE)*COS(PHI)*P(I,1)-SIN(PHI)*P(I,2)+
 699      &  SIN(THE)*COS(PHI)*P(I,3)
 700         DP(2)=COS(THE)*SIN(PHI)*P(I,1)+COS(PHI)*P(I,2)+
 701      &  SIN(THE)*SIN(PHI)*P(I,3)
 702         DP(3)=-SIN(THE)*P(I,1)+COS(THE)*P(I,3)
 703         DP(4)=P(I,4)
 704         DBP=BEX*DP(1)+BEY*DP(2)+BEZ*DP(3)
 705         DGABP=GA*(GA*DBP/(1D0+GA)+DP(4))
 706         P(I,1)=DP(1)+DGABP*BEX
 707         P(I,2)=DP(2)+DGABP*BEY
 708         P(I,3)=DP(3)+DGABP*BEZ
 709   350   P(I,4)=GA*(DP(4)+DBP)
 710       ENDIF
 711
 712 C...Weight with azimuthal distribution, if required.
 713       IF(MAZIP.NE.0.OR.MAZIC.NE.0) THEN
 714         DO 360 J=1,3
 715         DPT(1,J)=P(IM,J)
 716         DPT(2,J)=P(IAU,J)
 717   360   DPT(3,J)=P(N+1,J)
 718         DPMA=DPT(1,1)*DPT(2,1)+DPT(1,2)*DPT(2,2)+DPT(1,3)*DPT(2,3)
 719         DPMD=DPT(1,1)*DPT(3,1)+DPT(1,2)*DPT(3,2)+DPT(1,3)*DPT(3,3)
 720         DPMM=DPT(1,1)**2+DPT(1,2)**2+DPT(1,3)**2
 721         DO 370 J=1,3
 722         DPT(4,J)=DPT(2,J)-DPMA*DPT(1,J)/DPMM
 723   370   DPT(5,J)=DPT(3,J)-DPMD*DPT(1,J)/DPMM
 724         DPT(4,4)=SQRT(DPT(4,1)**2+DPT(4,2)**2+DPT(4,3)**2)
 725         DPT(5,4)=SQRT(DPT(5,1)**2+DPT(5,2)**2+DPT(5,3)**2)
 726         IF(MIN(DPT(4,4),DPT(5,4)).GT.0.1*PARJ(82)) THEN
 727           CAD=(DPT(4,1)*DPT(5,1)+DPT(4,2)*DPT(5,2)+
 728      &    DPT(4,3)*DPT(5,3))/(DPT(4,4)*DPT(5,4))
 729           IF(MAZIP.NE.0) THEN
 730             IF(1.+HAZIP*(2.*CAD**2-1.).LT.RLU_HIJING(0)*(1.+ABS(HAZIP)))
 731      &      GOTO 340
 732           ENDIF
 733           IF(MAZIC.NE.0) THEN
 734             IF(MAZIC.EQ.N+2) CAD=-CAD
 735             IF((1.-HAZIC)*(1.-HAZIC*CAD)/(1.+HAZIC**2-2.*HAZIC*CAD).
 736      &      LT.RLU_HIJING(0)) GOTO 340
 737           ENDIF
 738         ENDIF
 739       ENDIF
 740
 741 C...Continue loop over partons that may branch, until none left.
 742       IF(IGM.GE.0) K(IM,1)=14
 743       N=N+NEP
 744       NEP=2
 745       IF(N.GT.MSTU(4)-MSTU(32)-5) THEN
 746          CALL LUERRM_HIJING(11
 747      $        ,'(LUSHOW_HIJING:) no more memory left in LUJETS_HIJING')
 748         IF(MSTU(21).GE.1) N=NS
 749         IF(MSTU(21).GE.1) RETURN
 750       ENDIF
 751       GOTO 140
 752
 753 C...Set information on imagined shower initiator.
 754   380 IF(NPA.GE.2) THEN
 755         K(NS+1,1)=11
 756         K(NS+1,2)=94
 757         K(NS+1,3)=IP1
 758         IF(IP2.GT.0.AND.IP2.LT.IP1) K(NS+1,3)=IP2
 759         K(NS+1,4)=NS+2
 760         K(NS+1,5)=NS+1+NPA
 761         IIM=1
 762       ELSE
 763         IIM=0
 764       ENDIF
 765
 766 C...Reconstruct string drawing information.
 767       DO 390 I=NS+1+IIM,N
 768       IF(K(I,1).LE.10.AND.K(I,2).EQ.22) THEN
 769         K(I,1)=1
 770       ELSEIF(K(I,1).LE.10) THEN
 771         K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))
 772         K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))
 773       ELSEIF(K(MOD(K(I,4),MSTU(5))+1,2).NE.22) THEN
 774         ID1=MOD(K(I,4),MSTU(5))
 775         IF(K(I,2).GE.1.AND.K(I,2).LE.8) ID1=MOD(K(I,4),MSTU(5))+1
 776         ID2=2*MOD(K(I,4),MSTU(5))+1-ID1
 777         K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))+ID1
 778         K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))+ID2
 779         K(ID1,4)=K(ID1,4)+MSTU(5)*I
 780         K(ID1,5)=K(ID1,5)+MSTU(5)*ID2
 781         K(ID2,4)=K(ID2,4)+MSTU(5)*ID1
 782         K(ID2,5)=K(ID2,5)+MSTU(5)*I
 783       ELSE
 784         ID1=MOD(K(I,4),MSTU(5))
 785         ID2=ID1+1
 786         K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))+ID1
 787         K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))+ID1
 788         K(ID1,4)=K(ID1,4)+MSTU(5)*I
 789         K(ID1,5)=K(ID1,5)+MSTU(5)*I
 790         K(ID2,4)=0
 791         K(ID2,5)=0
 792       ENDIF
 793   390 CONTINUE
 794
 795 C...Transformation from CM frame.
 796       IF(NPA.GE.2) THEN
 797         BEX=PS(1)/PS(4)
 798         BEY=PS(2)/PS(4)
 799         BEZ=PS(3)/PS(4)
 800         GA=PS(4)/PS(5)
 801         GABEP=GA*(GA*(BEX*P(IPA(1),1)+BEY*P(IPA(1),2)+BEZ*P(IPA(1),3))
 802      &  /(1.+GA)-P(IPA(1),4))
 803       ELSE
 804         BEX=0.
 805         BEY=0.
 806         BEZ=0.
 807         GABEP=0.
 808       ENDIF
 809       THE=ULANGL_HIJING(P(IPA(1),3)+GABEP*BEZ,SQRT((P(IPA(1),1)
 810      &+GABEP*BEX)**2+(P(IPA(1),2)+GABEP*BEY)**2))
 811       PHI=ULANGL_HIJING(P(IPA(1),1)+GABEP*BEX,P(IPA(1),2)+GABEP*BEY)
 812       IF(NPA.EQ.3) THEN
 813          CHI=ULANGL_HIJING(COS(THE)*COS(PHI)*(P(IPA(2),1)+GABEP*BEX)
 814      $        +COS(THE)*SIN(PHI)*(P(IPA(2),2)+GABEP*BEY)-SIN(THE)
 815      $        *(P(IPA(2),3)+GABEP*BEZ),-SIN(PHI)*(P(IPA(2),1)+GABEP*BEX)
 816      $        +COS(PHI)*(P(IPA(2),2)+GABEP*BEY))
 817         CALL LUDBRB_HIJING(NS+1,N,0.,CHI,0D0,0D0,0D0)
 818       ENDIF
 819       DBEX=DBLE(BEX)
 820       DBEY=DBLE(BEY)
 821       DBEZ=DBLE(BEZ)
 822       CALL LUDBRB_HIJING(NS+1,N,THE,PHI,DBEX,DBEY,DBEZ)
 823
 824 C...Decay vertex of shower.
 825       DO 400 I=NS+1,N
 826       DO 400 J=1,5
 827   400 V(I,J)=V(IP1,J)
 828
 829 C...Delete trivial shower, else connect initiators.
 830       IF(N.EQ.NS+NPA+IIM) THEN
 831         N=NS
 832       ELSE
 833         DO 410 IP=1,NPA
 834         K(IPA(IP),1)=14
 835         K(IPA(IP),4)=K(IPA(IP),4)+NS+IIM+IP
 836         K(IPA(IP),5)=K(IPA(IP),5)+NS+IIM+IP
 837         K(NS+IIM+IP,3)=IPA(IP)
 838         IF(IIM.EQ.1.AND.MSTU(16).NE.2) K(NS+IIM+IP,3)=NS+1
 839         K(NS+IIM+IP,4)=MSTU(5)*IPA(IP)+K(NS+IIM+IP,4)
 840   410   K(NS+IIM+IP,5)=MSTU(5)*IPA(IP)+K(NS+IIM+IP,5)
 841       ENDIF
 842
 843       RETURN
 844       END