HIJING/hipyset1_35/lustrf_hijing.F

   1 * $Id$
   2
   3 C*********************************************************************
   4
   5       SUBROUTINE LUSTRF_HIJING(IP)
   6 C...Purpose: to handle the fragmentation of an arbitrary colour singlet
   7 C...jet system according to the Lund string fragmentation model.
   8       IMPLICIT DOUBLE PRECISION(D)
   9 #include "lujets_hijing.inc"
  10 #include "ludat1_hijing.inc"
  11 #include "ludat2_hijing.inc"
  12       DIMENSION DPS(5),KFL(3),PMQ(3),PX(3),PY(3),GAM(3),IE(2),PR(2),
  13      &IN(9),DHM(4),DHG(4),DP(5,5),IRANK(2),MJU(4),IJU(3),PJU(5,5),
  14      &TJU(5),KFJH(2),NJS(2),KFJS(2),PJS(4,5)
  15
  16 C...Function: four-product of two vectors.
  17       FOUR(I,J)=P(I,4)*P(J,4)-P(I,1)*P(J,1)-P(I,2)*P(J,2)-P(I,3)*P(J,3)
  18       DFOUR(I,J)=DP(I,4)*DP(J,4)-DP(I,1)*DP(J,1)-DP(I,2)*DP(J,2)-
  19      &DP(I,3)*DP(J,3)
  20
  21 C...Reset counters. Identify parton system.
  22       MSTJ(91)=0
  23       NSAV=N
  24       NP=0
  25       KQSUM=0
  26       DO 100 J=1,5
  27   100 DPS(J)=0.
  28       MJU(1)=0
  29       MJU(2)=0
  30       I=IP-1
  31   110 I=I+1
  32       IF(I.GT.MIN(N,MSTU(4)-MSTU(32))) THEN
  33          CALL LUERRM_HIJING(12
  34      $        ,'(LUSTRF_HIJING:) failed to reconstruct jet system')
  35         IF(MSTU(21).GE.1) RETURN
  36       ENDIF
  37       IF(K(I,1).NE.1.AND.K(I,1).NE.2.AND.K(I,1).NE.41) GOTO 110
  38       KC=LUCOMP_HIJING(K(I,2))
  39       IF(KC.EQ.0) GOTO 110
  40       KQ=KCHG(KC,2)*ISIGN(1,K(I,2))
  41       IF(KQ.EQ.0) GOTO 110
  42       IF(N+5*NP+11.GT.MSTU(4)-MSTU(32)-5) THEN
  43          CALL LUERRM_HIJING(11
  44      $        ,'(LUSTRF_HIJING:) no more memory left in LUJETS_HIJING')
  45         IF(MSTU(21).GE.1) RETURN
  46       ENDIF
  47
  48 C...Take copy of partons to be considered. Check flavour sum.
  49       NP=NP+1
  50       DO 120 J=1,5
  51       K(N+NP,J)=K(I,J)
  52       P(N+NP,J)=P(I,J)
  53   120 DPS(J)=DPS(J)+P(I,J)
  54       K(N+NP,3)=I
  55       IF(P(N+NP,4)**2.LT.P(N+NP,1)**2+P(N+NP,2)**2+P(N+NP,3)**2) THEN
  56         P(N+NP,4)=SQRT(P(N+NP,1)**2+P(N+NP,2)**2+P(N+NP,3)**2+
  57      &  P(N+NP,5)**2)
  58         DPS(4)=DPS(4)+MAX(0.,P(N+NP,4)-P(I,4))
  59       ENDIF
  60       IF(KQ.NE.2) KQSUM=KQSUM+KQ
  61       IF(K(I,1).EQ.41) THEN
  62         KQSUM=KQSUM+2*KQ
  63         IF(KQSUM.EQ.KQ) MJU(1)=N+NP
  64         IF(KQSUM.NE.KQ) MJU(2)=N+NP
  65       ENDIF
  66       IF(K(I,1).EQ.2.OR.K(I,1).EQ.41) GOTO 110
  67       IF(KQSUM.NE.0) THEN
  68          CALL LUERRM_HIJING(12
  69      $        ,'(LUSTRF_HIJING:) unphysical flavour combination')
  70         IF(MSTU(21).GE.1) RETURN
  71       ENDIF
  72
  73 C...Boost copied system to CM frame (for better numerical precision).
  74       CALL LUDBRB_HIJING(N+1,N+NP,0.,0.,-DPS(1)/DPS(4),-DPS(2)/DPS(4),
  75      &-DPS(3)/DPS(4))
  76
  77 C...Search for very nearby partons that may be recombined.
  78       NTRYR=0
  79       PARU12=PARU(12)
  80       PARU13=PARU(13)
  81       MJU(3)=MJU(1)
  82       MJU(4)=MJU(2)
  83       NR=NP
  84   130 IF(NR.GE.3) THEN
  85         PDRMIN=2.*PARU12
  86         DO 140 I=N+1,N+NR
  87         IF(I.EQ.N+NR.AND.IABS(K(N+1,2)).NE.21) GOTO 140
  88         I1=I+1
  89         IF(I.EQ.N+NR) I1=N+1
  90         IF(K(I,1).EQ.41.OR.K(I1,1).EQ.41) GOTO 140
  91         IF(MJU(1).NE.0.AND.I1.LT.MJU(1).AND.IABS(K(I1,2)).NE.21)
  92      &  GOTO 140
  93         IF(MJU(2).NE.0.AND.I.GT.MJU(2).AND.IABS(K(I,2)).NE.21) GOTO 140
  94         PAP=SQRT((P(I,1)**2+P(I,2)**2+P(I,3)**2)*(P(I1,1)**2+
  95      &  P(I1,2)**2+P(I1,3)**2))
  96         PVP=P(I,1)*P(I1,1)+P(I,2)*P(I1,2)+P(I,3)*P(I1,3)
  97         PDR=4.*(PAP-PVP)**2/(PARU13**2*PAP+2.*(PAP-PVP))
  98         IF(PDR.LT.PDRMIN) THEN
  99           IR=I
 100           PDRMIN=PDR
 101         ENDIF
 102   140   CONTINUE
 103
 104 C...Recombine very nearby partons to avoid machine precision problems.
 105         IF(PDRMIN.LT.PARU12.AND.IR.EQ.N+NR) THEN
 106           DO 150 J=1,4
 107   150     P(N+1,J)=P(N+1,J)+P(N+NR,J)
 108           P(N+1,5)=SQRT(MAX(0.,P(N+1,4)**2-P(N+1,1)**2-P(N+1,2)**2-
 109      &    P(N+1,3)**2))
 110           NR=NR-1
 111           GOTO 130
 112         ELSEIF(PDRMIN.LT.PARU12) THEN
 113           DO 160 J=1,4
 114   160     P(IR,J)=P(IR,J)+P(IR+1,J)
 115           P(IR,5)=SQRT(MAX(0.,P(IR,4)**2-P(IR,1)**2-P(IR,2)**2-
 116      &    P(IR,3)**2))
 117           DO 170 I=IR+1,N+NR-1
 118           K(I,2)=K(I+1,2)
 119           DO 170 J=1,5
 120   170     P(I,J)=P(I+1,J)
 121           IF(IR.EQ.N+NR-1) K(IR,2)=K(N+NR,2)
 122           NR=NR-1
 123           IF(MJU(1).GT.IR) MJU(1)=MJU(1)-1
 124           IF(MJU(2).GT.IR) MJU(2)=MJU(2)-1
 125           GOTO 130
 126         ENDIF
 127       ENDIF
 128       NTRYR=NTRYR+1
 129
 130 C...Reset particle counter. Skip ahead if no junctions are present;
 131 C...this is usually the case!
 132       NRS=MAX(5*NR+11,NP)
 133       NTRY=0
 134   180 NTRY=NTRY+1
 135       IF(NTRY.GT.100.AND.NTRYR.LE.4) THEN
 136         PARU12=4.*PARU12
 137         PARU13=2.*PARU13
 138         GOTO 130
 139       ELSEIF(NTRY.GT.100) THEN
 140          CALL LUERRM_HIJING(14
 141      $        ,'(LUSTRF_HIJING:) caught in infinite loop')
 142         IF(MSTU(21).GE.1) RETURN
 143       ENDIF
 144       I=N+NRS
 145       IF(MJU(1).EQ.0.AND.MJU(2).EQ.0) GOTO 500
 146       DO 490 JT=1,2
 147       NJS(JT)=0
 148       IF(MJU(JT).EQ.0) GOTO 490
 149       JS=3-2*JT
 150
 151 C...Find and sum up momentum on three sides of junction. Check flavours.
 152       DO 190 IU=1,3
 153       IJU(IU)=0
 154       DO 190 J=1,5
 155   190 PJU(IU,J)=0.
 156       IU=0
 157       DO 200 I1=N+1+(JT-1)*(NR-1),N+NR+(JT-1)*(1-NR),JS
 158       IF(K(I1,2).NE.21.AND.IU.LE.2) THEN
 159         IU=IU+1
 160         IJU(IU)=I1
 161       ENDIF
 162       DO 200 J=1,4
 163   200 PJU(IU,J)=PJU(IU,J)+P(I1,J)
 164       DO 210 IU=1,3
 165   210 PJU(IU,5)=SQRT(PJU(IU,1)**2+PJU(IU,2)**2+PJU(IU,3)**2)
 166       IF(K(IJU(3),2)/100.NE.10*K(IJU(1),2)+K(IJU(2),2).AND.
 167      &K(IJU(3),2)/100.NE.10*K(IJU(2),2)+K(IJU(1),2)) THEN
 168          CALL LUERRM_HIJING(12
 169      $        ,'(LUSTRF_HIJING:) unphysical flavour combination')
 170         IF(MSTU(21).GE.1) RETURN
 171       ENDIF
 172
 173 C...Calculate (approximate) boost to rest frame of junction.
 174       T12=(PJU(1,1)*PJU(2,1)+PJU(1,2)*PJU(2,2)+PJU(1,3)*PJU(2,3))/
 175      &(PJU(1,5)*PJU(2,5))
 176       T13=(PJU(1,1)*PJU(3,1)+PJU(1,2)*PJU(3,2)+PJU(1,3)*PJU(3,3))/
 177      &(PJU(1,5)*PJU(3,5))
 178       T23=(PJU(2,1)*PJU(3,1)+PJU(2,2)*PJU(3,2)+PJU(2,3)*PJU(3,3))/
 179      &(PJU(2,5)*PJU(3,5))
 180       T11=SQRT((2./3.)*(1.-T12)*(1.-T13)/(1.-T23))
 181       T22=SQRT((2./3.)*(1.-T12)*(1.-T23)/(1.-T13))
 182       TSQ=SQRT((2.*T11*T22+T12-1.)*(1.+T12))
 183       T1F=(TSQ-T22*(1.+T12))/(1.-T12**2)
 184       T2F=(TSQ-T11*(1.+T12))/(1.-T12**2)
 185       DO 220 J=1,3
 186   220 TJU(J)=-(T1F*PJU(1,J)/PJU(1,5)+T2F*PJU(2,J)/PJU(2,5))
 187       TJU(4)=SQRT(1.+TJU(1)**2+TJU(2)**2+TJU(3)**2)
 188       DO 230 IU=1,3
 189   230 PJU(IU,5)=TJU(4)*PJU(IU,4)-TJU(1)*PJU(IU,1)-TJU(2)*PJU(IU,2)-
 190      &TJU(3)*PJU(IU,3)
 191
 192 C...Put junction at rest if motion could give inconsistencies.
 193       IF(PJU(1,5)+PJU(2,5).GT.PJU(1,4)+PJU(2,4)) THEN
 194         DO 240 J=1,3
 195   240   TJU(J)=0.
 196         TJU(4)=1.
 197         PJU(1,5)=PJU(1,4)
 198         PJU(2,5)=PJU(2,4)
 199         PJU(3,5)=PJU(3,4)
 200       ENDIF
 201
 202 C...Start preparing for fragmentation of two strings from junction.
 203       ISTA=I
 204       DO 470 IU=1,2
 205       NS=IJU(IU+1)-IJU(IU)
 206
 207 C...Junction strings: find longitudinal string directions.
 208       DO 260 IS=1,NS
 209       IS1=IJU(IU)+IS-1
 210       IS2=IJU(IU)+IS
 211       DO 250 J=1,5
 212       DP(1,J)=0.5*P(IS1,J)
 213       IF(IS.EQ.1) DP(1,J)=P(IS1,J)
 214       DP(2,J)=0.5*P(IS2,J)
 215   250 IF(IS.EQ.NS) DP(2,J)=-PJU(IU,J)
 216       IF(IS.EQ.NS) DP(2,4)=SQRT(PJU(IU,1)**2+PJU(IU,2)**2+PJU(IU,3)**2)
 217       IF(IS.EQ.NS) DP(2,5)=0.
 218       DP(3,5)=DFOUR(1,1)
 219       DP(4,5)=DFOUR(2,2)
 220       DHKC=DFOUR(1,2)
 221       IF(DP(3,5)+2.*DHKC+DP(4,5).LE.0.) THEN
 222         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)
 223         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)
 224         DP(3,5)=0D0
 225         DP(4,5)=0D0
 226         DHKC=DFOUR(1,2)
 227       ENDIF
 228       DHKS=SQRT(DHKC**2-DP(3,5)*DP(4,5))
 229       DHK1=0.5*((DP(4,5)+DHKC)/DHKS-1.)
 230       DHK2=0.5*((DP(3,5)+DHKC)/DHKS-1.)
 231       IN1=N+NR+4*IS-3
 232       P(IN1,5)=SQRT(DP(3,5)+2.*DHKC+DP(4,5))
 233       DO 260 J=1,4
 234       P(IN1,J)=(1.+DHK1)*DP(1,J)-DHK2*DP(2,J)
 235   260 P(IN1+1,J)=(1.+DHK2)*DP(2,J)-DHK1*DP(1,J)
 236
 237 C...Junction strings: initialize flavour, momentum and starting pos.
 238       ISAV=I
 239   270 NTRY=NTRY+1
 240       IF(NTRY.GT.100.AND.NTRYR.LE.4) THEN
 241         PARU12=4.*PARU12
 242         PARU13=2.*PARU13
 243         GOTO 130
 244       ELSEIF(NTRY.GT.100) THEN
 245          CALL LUERRM_HIJING(14
 246      $        ,'(LUSTRF_HIJING:) caught in infinite loop')
 247         IF(MSTU(21).GE.1) RETURN
 248       ENDIF
 249       I=ISAV
 250       IRANKJ=0
 251       IE(1)=K(N+1+(JT/2)*(NP-1),3)
 252       IN(4)=N+NR+1
 253       IN(5)=IN(4)+1
 254       IN(6)=N+NR+4*NS+1
 255       DO 280 JQ=1,2
 256       DO 280 IN1=N+NR+2+JQ,N+NR+4*NS-2+JQ,4
 257       P(IN1,1)=2-JQ
 258       P(IN1,2)=JQ-1
 259   280 P(IN1,3)=1.
 260       KFL(1)=K(IJU(IU),2)
 261       PX(1)=0.
 262       PY(1)=0.
 263       GAM(1)=0.
 264       DO 290 J=1,5
 265   290 PJU(IU+3,J)=0.
 266
 267 C...Junction strings: find initial transverse directions.
 268       DO 300 J=1,4
 269       DP(1,J)=P(IN(4),J)
 270       DP(2,J)=P(IN(4)+1,J)
 271       DP(3,J)=0.
 272   300 DP(4,J)=0.
 273       DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)
 274       DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)
 275       DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4)
 276       DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4)
 277       DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4)
 278       IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1.
 279       IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1.
 280       IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1.
 281       IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1.
 282       DHC12=DFOUR(1,2)
 283       DHCX1=DFOUR(3,1)/DHC12
 284       DHCX2=DFOUR(3,2)/DHC12
 285       DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12)
 286       DHCY1=DFOUR(4,1)/DHC12
 287       DHCY2=DFOUR(4,2)/DHC12
 288       DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12
 289       DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2)
 290       DO 310 J=1,4
 291       DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J))
 292       P(IN(6),J)=DP(3,J)
 293   310 P(IN(6)+1,J)=DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)-
 294      &DHCYX*DP(3,J))
 295
 296 C...Junction strings: produce new particle, origin.
 297   320 I=I+1
 298       IF(2*I-NSAV.GE.MSTU(4)-MSTU(32)-5) THEN
 299          CALL LUERRM_HIJING(11
 300      $        ,'(LUSTRF_HIJING:) no more memory left in LUJETS_HIJING')
 301         IF(MSTU(21).GE.1) RETURN
 302       ENDIF
 303       IRANKJ=IRANKJ+1
 304       K(I,1)=1
 305       K(I,3)=IE(1)
 306       K(I,4)=0
 307       K(I,5)=0
 308
 309 C...Junction strings: generate flavour, hadron, pT, z and Gamma.
 310   330 CALL LUKFDI_HIJING(KFL(1),0,KFL(3),K(I,2))
 311       IF(K(I,2).EQ.0) GOTO 270
 312       IF(MSTJ(12).GE.3.AND.IRANKJ.EQ.1.AND.IABS(KFL(1)).LE.10.AND.
 313      &IABS(KFL(3)).GT.10) THEN
 314         IF(RLU_HIJING(0).GT.PARJ(19)) GOTO 330
 315       ENDIF
 316       P(I,5)=ULMASS_HIJING(K(I,2))
 317       CALL LUPTDI_HIJING(KFL(1),PX(3),PY(3))
 318       PR(1)=P(I,5)**2+(PX(1)+PX(3))**2+(PY(1)+PY(3))**2
 319       CALL LUZDIS_HIJING(KFL(1),KFL(3),PR(1),Z)
 320       GAM(3)=(1.-Z)*(GAM(1)+PR(1)/Z)
 321       DO 340 J=1,3
 322   340 IN(J)=IN(3+J)
 323
 324 C...Junction strings: stepping within or from 'low' string region easy.
 325       IF(IN(1)+1.EQ.IN(2).AND.Z*P(IN(1)+2,3)*P(IN(2)+2,3)*
 326      &P(IN(1),5)**2.GE.PR(1)) THEN
 327         P(IN(1)+2,4)=Z*P(IN(1)+2,3)
 328         P(IN(2)+2,4)=PR(1)/(P(IN(1)+2,4)*P(IN(1),5)**2)
 329         DO 350 J=1,4
 330   350   P(I,J)=(PX(1)+PX(3))*P(IN(3),J)+(PY(1)+PY(3))*P(IN(3)+1,J)
 331         GOTO 420
 332       ELSEIF(IN(1)+1.EQ.IN(2)) THEN
 333         P(IN(2)+2,4)=P(IN(2)+2,3)
 334         P(IN(2)+2,1)=1.
 335         IN(2)=IN(2)+4
 336         IF(IN(2).GT.N+NR+4*NS) GOTO 270
 337         IF(FOUR(IN(1),IN(2)).LE.1E-2) THEN
 338           P(IN(1)+2,4)=P(IN(1)+2,3)
 339           P(IN(1)+2,1)=0.
 340           IN(1)=IN(1)+4
 341         ENDIF
 342       ENDIF
 343
 344 C...Junction strings: find new transverse directions.
 345   360 IF(IN(1).GT.N+NR+4*NS.OR.IN(2).GT.N+NR+4*NS.OR.
 346      &IN(1).GT.IN(2)) GOTO 270
 347       IF(IN(1).NE.IN(4).OR.IN(2).NE.IN(5)) THEN
 348         DO 370 J=1,4
 349         DP(1,J)=P(IN(1),J)
 350         DP(2,J)=P(IN(2),J)
 351         DP(3,J)=0.
 352   370   DP(4,J)=0.
 353         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)
 354         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)
 355         DHC12=DFOUR(1,2)
 356         IF(DHC12.LE.1E-2) THEN
 357           P(IN(1)+2,4)=P(IN(1)+2,3)
 358           P(IN(1)+2,1)=0.
 359           IN(1)=IN(1)+4
 360           GOTO 360
 361         ENDIF
 362         IN(3)=N+NR+4*NS+5
 363         DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4)
 364         DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4)
 365         DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4)
 366         IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1.
 367         IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1.
 368         IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1.
 369         IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1.
 370         DHCX1=DFOUR(3,1)/DHC12
 371         DHCX2=DFOUR(3,2)/DHC12
 372         DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12)
 373         DHCY1=DFOUR(4,1)/DHC12
 374         DHCY2=DFOUR(4,2)/DHC12
 375         DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12
 376         DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2)
 377         DO 380 J=1,4
 378         DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J))
 379         P(IN(3),J)=DP(3,J)
 380   380   P(IN(3)+1,J)=DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)-
 381      &  DHCYX*DP(3,J))
 382 C...Express pT with respect to new axes, if sensible.
 383         PXP=-(PX(3)*FOUR(IN(6),IN(3))+PY(3)*FOUR(IN(6)+1,IN(3)))
 384         PYP=-(PX(3)*FOUR(IN(6),IN(3)+1)+PY(3)*FOUR(IN(6)+1,IN(3)+1))
 385         IF(ABS(PXP**2+PYP**2-PX(3)**2-PY(3)**2).LT.0.01) THEN
 386           PX(3)=PXP
 387           PY(3)=PYP
 388         ENDIF
 389       ENDIF
 390
 391 C...Junction strings: sum up known four-momentum, coefficients for m2.
 392       DO 400 J=1,4
 393       DHG(J)=0.
 394       P(I,J)=PX(1)*P(IN(6),J)+PY(1)*P(IN(6)+1,J)+PX(3)*P(IN(3),J)+
 395      &PY(3)*P(IN(3)+1,J)
 396       DO 390 IN1=IN(4),IN(1)-4,4
 397   390 P(I,J)=P(I,J)+P(IN1+2,3)*P(IN1,J)
 398       DO 400 IN2=IN(5),IN(2)-4,4
 399   400 P(I,J)=P(I,J)+P(IN2+2,3)*P(IN2,J)
 400       DHM(1)=FOUR(I,I)
 401       DHM(2)=2.*FOUR(I,IN(1))
 402       DHM(3)=2.*FOUR(I,IN(2))
 403       DHM(4)=2.*FOUR(IN(1),IN(2))
 404
 405 C...Junction strings: find coefficients for Gamma expression.
 406       DO 410 IN2=IN(1)+1,IN(2),4
 407       DO 410 IN1=IN(1),IN2-1,4
 408       DHC=2.*FOUR(IN1,IN2)
 409       DHG(1)=DHG(1)+P(IN1+2,1)*P(IN2+2,1)*DHC
 410       IF(IN1.EQ.IN(1)) DHG(2)=DHG(2)-P(IN2+2,1)*DHC
 411       IF(IN2.EQ.IN(2)) DHG(3)=DHG(3)+P(IN1+2,1)*DHC
 412   410 IF(IN1.EQ.IN(1).AND.IN2.EQ.IN(2)) DHG(4)=DHG(4)-DHC
 413
 414 C...Junction strings: solve (m2, Gamma) equation system for energies.
 415       DHS1=DHM(3)*DHG(4)-DHM(4)*DHG(3)
 416       IF(ABS(DHS1).LT.1E-4) GOTO 270
 417       DHS2=DHM(4)*(GAM(3)-DHG(1))-DHM(2)*DHG(3)-DHG(4)*
 418      &(P(I,5)**2-DHM(1))+DHG(2)*DHM(3)
 419       DHS3=DHM(2)*(GAM(3)-DHG(1))-DHG(2)*(P(I,5)**2-DHM(1))
 420       P(IN(2)+2,4)=0.5*(SQRT(MAX(0D0,DHS2**2-4.*DHS1*DHS3))/ABS(DHS1)-
 421      &DHS2/DHS1)
 422       IF(DHM(2)+DHM(4)*P(IN(2)+2,4).LE.0.) GOTO 270
 423       P(IN(1)+2,4)=(P(I,5)**2-DHM(1)-DHM(3)*P(IN(2)+2,4))/
 424      &(DHM(2)+DHM(4)*P(IN(2)+2,4))
 425
 426 C...Junction strings: step to new region if necessary.
 427       IF(P(IN(2)+2,4).GT.P(IN(2)+2,3)) THEN
 428         P(IN(2)+2,4)=P(IN(2)+2,3)
 429         P(IN(2)+2,1)=1.
 430         IN(2)=IN(2)+4
 431         IF(IN(2).GT.N+NR+4*NS) GOTO 270
 432         IF(FOUR(IN(1),IN(2)).LE.1E-2) THEN
 433           P(IN(1)+2,4)=P(IN(1)+2,3)
 434           P(IN(1)+2,1)=0.
 435           IN(1)=IN(1)+4
 436         ENDIF
 437         GOTO 360
 438       ELSEIF(P(IN(1)+2,4).GT.P(IN(1)+2,3)) THEN
 439         P(IN(1)+2,4)=P(IN(1)+2,3)
 440         P(IN(1)+2,1)=0.
 441         IN(1)=IN(1)+JS
 442         GOTO 710
 443       ENDIF
 444
 445 C...Junction strings: particle four-momentum, remainder, loop back.
 446   420 DO 430 J=1,4
 447       P(I,J)=P(I,J)+P(IN(1)+2,4)*P(IN(1),J)+P(IN(2)+2,4)*P(IN(2),J)
 448   430 PJU(IU+3,J)=PJU(IU+3,J)+P(I,J)
 449       IF(P(I,4).LE.0.) GOTO 270
 450       PJU(IU+3,5)=TJU(4)*PJU(IU+3,4)-TJU(1)*PJU(IU+3,1)-
 451      &TJU(2)*PJU(IU+3,2)-TJU(3)*PJU(IU+3,3)
 452       IF(PJU(IU+3,5).LT.PJU(IU,5)) THEN
 453         KFL(1)=-KFL(3)
 454         PX(1)=-PX(3)
 455         PY(1)=-PY(3)
 456         GAM(1)=GAM(3)
 457         IF(IN(3).NE.IN(6)) THEN
 458           DO 440 J=1,4
 459           P(IN(6),J)=P(IN(3),J)
 460   440     P(IN(6)+1,J)=P(IN(3)+1,J)
 461         ENDIF
 462         DO 450 JQ=1,2
 463         IN(3+JQ)=IN(JQ)
 464         P(IN(JQ)+2,3)=P(IN(JQ)+2,3)-P(IN(JQ)+2,4)
 465   450   P(IN(JQ)+2,1)=P(IN(JQ)+2,1)-(3-2*JQ)*P(IN(JQ)+2,4)
 466         GOTO 320
 467       ENDIF
 468
 469 C...Junction strings: save quantities left after each string.
 470       IF(IABS(KFL(1)).GT.10) GOTO 270
 471       I=I-1
 472       KFJH(IU)=KFL(1)
 473       DO 460 J=1,4
 474   460 PJU(IU+3,J)=PJU(IU+3,J)-P(I+1,J)
 475   470 CONTINUE
 476
 477 C...Junction strings: put together to new effective string endpoint.
 478       NJS(JT)=I-ISTA
 479       KFJS(JT)=K(K(MJU(JT+2),3),2)
 480       KFLS=2*INT(RLU_HIJING(0)+3.*PARJ(4)/(1.+3.*PARJ(4)))+1
 481       IF(KFJH(1).EQ.KFJH(2)) KFLS=3
 482       IF(ISTA.NE.I) KFJS(JT)=ISIGN(1000*MAX(IABS(KFJH(1)),
 483      &IABS(KFJH(2)))+100*MIN(IABS(KFJH(1)),IABS(KFJH(2)))+
 484      &KFLS,KFJH(1))
 485       DO 480 J=1,4
 486       PJS(JT,J)=PJU(1,J)+PJU(2,J)+P(MJU(JT),J)
 487   480 PJS(JT+2,J)=PJU(4,J)+PJU(5,J)
 488       PJS(JT,5)=SQRT(MAX(0.,PJS(JT,4)**2-PJS(JT,1)**2-PJS(JT,2)**2-
 489      &PJS(JT,3)**2))
 490   490 CONTINUE
 491
 492 C...Open versus closed strings. Choose breakup region for latter.
 493   500 IF(MJU(1).NE.0.AND.MJU(2).NE.0) THEN
 494         NS=MJU(2)-MJU(1)
 495         NB=MJU(1)-N
 496       ELSEIF(MJU(1).NE.0) THEN
 497         NS=N+NR-MJU(1)
 498         NB=MJU(1)-N
 499       ELSEIF(MJU(2).NE.0) THEN
 500         NS=MJU(2)-N
 501         NB=1
 502       ELSEIF(IABS(K(N+1,2)).NE.21) THEN
 503         NS=NR-1
 504         NB=1
 505       ELSE
 506         NS=NR+1
 507         W2SUM=0.
 508         DO 510 IS=1,NR
 509         P(N+NR+IS,1)=0.5*FOUR(N+IS,N+IS+1-NR*(IS/NR))
 510   510   W2SUM=W2SUM+P(N+NR+IS,1)
 511         W2RAN=RLU_HIJING(0)*W2SUM
 512         NB=0
 513   520   NB=NB+1
 514         W2SUM=W2SUM-P(N+NR+NB,1)
 515         IF(W2SUM.GT.W2RAN.AND.NB.LT.NR) GOTO 520
 516       ENDIF
 517
 518 C...Find longitudinal string directions (i.e. lightlike four-vectors).
 519       DO 540 IS=1,NS
 520       IS1=N+IS+NB-1-NR*((IS+NB-2)/NR)
 521       IS2=N+IS+NB-NR*((IS+NB-1)/NR)
 522       DO 530 J=1,5
 523       DP(1,J)=P(IS1,J)
 524       IF(IABS(K(IS1,2)).EQ.21) DP(1,J)=0.5*DP(1,J)
 525       IF(IS1.EQ.MJU(1)) DP(1,J)=PJS(1,J)-PJS(3,J)
 526       DP(2,J)=P(IS2,J)
 527       IF(IABS(K(IS2,2)).EQ.21) DP(2,J)=0.5*DP(2,J)
 528   530 IF(IS2.EQ.MJU(2)) DP(2,J)=PJS(2,J)-PJS(4,J)
 529       DP(3,5)=DFOUR(1,1)
 530       DP(4,5)=DFOUR(2,2)
 531       DHKC=DFOUR(1,2)
 532       IF(DP(3,5)+2.*DHKC+DP(4,5).LE.0.) THEN
 533         DP(3,5)=DP(1,5)**2
 534         DP(4,5)=DP(2,5)**2
 535         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2+DP(1,5)**2)
 536         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2+DP(2,5)**2)
 537         DHKC=DFOUR(1,2)
 538       ENDIF
 539       DHKS=SQRT(DHKC**2-DP(3,5)*DP(4,5))
 540       DHK1=0.5*((DP(4,5)+DHKC)/DHKS-1.)
 541       DHK2=0.5*((DP(3,5)+DHKC)/DHKS-1.)
 542       IN1=N+NR+4*IS-3
 543       P(IN1,5)=SQRT(DP(3,5)+2.*DHKC+DP(4,5))
 544       DO 540 J=1,4
 545       P(IN1,J)=(1.+DHK1)*DP(1,J)-DHK2*DP(2,J)
 546   540 P(IN1+1,J)=(1.+DHK2)*DP(2,J)-DHK1*DP(1,J)
 547
 548 C...Begin initialization: sum up energy, set starting position.
 549       ISAV=I
 550   550 NTRY=NTRY+1
 551       IF(NTRY.GT.100.AND.NTRYR.LE.4) THEN
 552         PARU12=4.*PARU12
 553         PARU13=2.*PARU13
 554         GOTO 130
 555       ELSEIF(NTRY.GT.100) THEN
 556          CALL LUERRM_HIJING(14
 557      $        ,'(LUSTRF_HIJING:) caught in infinite loop')
 558         IF(MSTU(21).GE.1) RETURN
 559       ENDIF
 560       I=ISAV
 561       DO 560 J=1,4
 562       P(N+NRS,J)=0.
 563       DO 560 IS=1,NR
 564   560 P(N+NRS,J)=P(N+NRS,J)+P(N+IS,J)
 565       DO 570 JT=1,2
 566       IRANK(JT)=0
 567       IF(MJU(JT).NE.0) IRANK(JT)=NJS(JT)
 568       IF(NS.GT.NR) IRANK(JT)=1
 569       IE(JT)=K(N+1+(JT/2)*(NP-1),3)
 570       IN(3*JT+1)=N+NR+1+4*(JT/2)*(NS-1)
 571       IN(3*JT+2)=IN(3*JT+1)+1
 572       IN(3*JT+3)=N+NR+4*NS+2*JT-1
 573       DO 570 IN1=N+NR+2+JT,N+NR+4*NS-2+JT,4
 574       P(IN1,1)=2-JT
 575       P(IN1,2)=JT-1
 576   570 P(IN1,3)=1.
 577
 578 C...Initialize flavour and pT variables for open string.
 579       IF(NS.LT.NR) THEN
 580         PX(1)=0.
 581         PY(1)=0.
 582         IF(NS.EQ.1.AND.MJU(1)+MJU(2).EQ.0) CALL LUPTDI_HIJING(0,PX(1)
 583      $       ,PY(1))
 584         PX(2)=-PX(1)
 585         PY(2)=-PY(1)
 586         DO 580 JT=1,2
 587         KFL(JT)=K(IE(JT),2)
 588         IF(MJU(JT).NE.0) KFL(JT)=KFJS(JT)
 589         MSTJ(93)=1
 590         PMQ(JT)=ULMASS_HIJING(KFL(JT))
 591   580   GAM(JT)=0.
 592
 593 C...Closed string: random initial breakup flavour, pT and vertex.
 594       ELSE
 595          KFL(3)=INT(1.+(2.+PARJ(2))*RLU_HIJING(0))*(-1)
 596      $        **INT(RLU_HIJING(0)+0.5)
 597         CALL LUKFDI_HIJING(KFL(3),0,KFL(1),KDUMP)
 598         KFL(2)=-KFL(1)
 599         IF(IABS(KFL(1)).GT.10.AND.RLU_HIJING(0).GT.0.5) THEN
 600           KFL(2)=-(KFL(1)+ISIGN(10000,KFL(1)))
 601         ELSEIF(IABS(KFL(1)).GT.10) THEN
 602           KFL(1)=-(KFL(2)+ISIGN(10000,KFL(2)))
 603         ENDIF
 604         CALL LUPTDI_HIJING(KFL(1),PX(1),PY(1))
 605         PX(2)=-PX(1)
 606         PY(2)=-PY(1)
 607         PR3=MIN(25.,0.1*P(N+NR+1,5)**2)
 608   590   CALL LUZDIS_HIJING(KFL(1),KFL(2),PR3,Z)
 609         ZR=PR3/(Z*P(N+NR+1,5)**2)
 610         IF(ZR.GE.1.) GOTO 590
 611         DO 600 JT=1,2
 612         MSTJ(93)=1
 613         PMQ(JT)=ULMASS_HIJING(KFL(JT))
 614         GAM(JT)=PR3*(1.-Z)/Z
 615         IN1=N+NR+3+4*(JT/2)*(NS-1)
 616         P(IN1,JT)=1.-Z
 617         P(IN1,3-JT)=JT-1
 618         P(IN1,3)=(2-JT)*(1.-Z)+(JT-1)*Z
 619         P(IN1+1,JT)=ZR
 620         P(IN1+1,3-JT)=2-JT
 621   600   P(IN1+1,3)=(2-JT)*(1.-ZR)+(JT-1)*ZR
 622       ENDIF
 623
 624 C...Find initial transverse directions (i.e. spacelike four-vectors).
 625       DO 640 JT=1,2
 626       IF(JT.EQ.1.OR.NS.EQ.NR-1) THEN
 627         IN1=IN(3*JT+1)
 628         IN3=IN(3*JT+3)
 629         DO 610 J=1,4
 630         DP(1,J)=P(IN1,J)
 631         DP(2,J)=P(IN1+1,J)
 632         DP(3,J)=0.
 633   610   DP(4,J)=0.
 634         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)
 635         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)
 636         DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4)
 637         DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4)
 638         DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4)
 639         IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1.
 640         IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1.
 641         IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1.
 642         IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1.
 643         DHC12=DFOUR(1,2)
 644         DHCX1=DFOUR(3,1)/DHC12
 645         DHCX2=DFOUR(3,2)/DHC12
 646         DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12)
 647         DHCY1=DFOUR(4,1)/DHC12
 648         DHCY2=DFOUR(4,2)/DHC12
 649         DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12
 650         DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2)
 651         DO 620 J=1,4
 652         DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J))
 653         P(IN3,J)=DP(3,J)
 654   620   P(IN3+1,J)=DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)-
 655      &  DHCYX*DP(3,J))
 656       ELSE
 657         DO 630 J=1,4
 658         P(IN3+2,J)=P(IN3,J)
 659   630   P(IN3+3,J)=P(IN3+1,J)
 660       ENDIF
 661   640 CONTINUE
 662
 663 C...Remove energy used up in junction string fragmentation.
 664       IF(MJU(1)+MJU(2).GT.0) THEN
 665         DO 660 JT=1,2
 666         IF(NJS(JT).EQ.0) GOTO 660
 667         DO 650 J=1,4
 668   650   P(N+NRS,J)=P(N+NRS,J)-PJS(JT+2,J)
 669   660   CONTINUE
 670       ENDIF
 671
 672 C...Produce new particle: side, origin.
 673   670 I=I+1
 674       IF(2*I-NSAV.GE.MSTU(4)-MSTU(32)-5) THEN
 675          CALL LUERRM_HIJING(11
 676      $        ,'(LUSTRF_HIJING:) no more memory left in LUJETS_HIJING')
 677         IF(MSTU(21).GE.1) RETURN
 678       ENDIF
 679       JT=1.5+RLU_HIJING(0)
 680       IF(IABS(KFL(3-JT)).GT.10) JT=3-JT
 681       JR=3-JT
 682       JS=3-2*JT
 683       IRANK(JT)=IRANK(JT)+1
 684       K(I,1)=1
 685       K(I,3)=IE(JT)
 686       K(I,4)=0
 687       K(I,5)=0
 688
 689 C...Generate flavour, hadron and pT.
 690   680 CALL LUKFDI_HIJING(KFL(JT),0,KFL(3),K(I,2))
 691       IF(K(I,2).EQ.0) GOTO 550
 692       IF(MSTJ(12).GE.3.AND.IRANK(JT).EQ.1.AND.IABS(KFL(JT)).LE.10.AND.
 693      &IABS(KFL(3)).GT.10) THEN
 694         IF(RLU_HIJING(0).GT.PARJ(19)) GOTO 680
 695       ENDIF
 696       P(I,5)=ULMASS_HIJING(K(I,2))
 697       CALL LUPTDI_HIJING(KFL(JT),PX(3),PY(3))
 698       PR(JT)=P(I,5)**2+(PX(JT)+PX(3))**2+(PY(JT)+PY(3))**2
 699
 700 C...Final hadrons for small invariant mass.
 701       MSTJ(93)=1
 702       PMQ(3)=ULMASS_HIJING(KFL(3))
 703       WMIN=PARJ(32+MSTJ(11))+PMQ(1)+PMQ(2)+PARJ(36)*PMQ(3)
 704       IF(IABS(KFL(JT)).GT.10.AND.IABS(KFL(3)).GT.10) WMIN=
 705      &WMIN-0.5*PARJ(36)*PMQ(3)
 706       WREM2=FOUR(N+NRS,N+NRS)
 707       IF(WREM2.LT.0.10) GOTO 550
 708       IF(WREM2.LT.MAX(WMIN*(1.+(2.*RLU_HIJING(0)-1.)*PARJ(37)),
 709      &PARJ(32)+PMQ(1)+PMQ(2))**2) GOTO 810
 710
 711 C...Choose z, which gives Gamma. Shift z for heavy flavours.
 712       CALL LUZDIS_HIJING(KFL(JT),KFL(3),PR(JT),Z)
 713       KFL1A=IABS(KFL(1))
 714       KFL2A=IABS(KFL(2))
 715       IF(MAX(MOD(KFL1A,10),MOD(KFL1A/1000,10),MOD(KFL2A,10),
 716      &MOD(KFL2A/1000,10)).GE.4) THEN
 717         PR(JR)=(PMQ(JR)+PMQ(3))**2+(PX(JR)-PX(3))**2+(PY(JR)-PY(3))**2
 718         PW12=SQRT(MAX(0.,(WREM2-PR(1)-PR(2))**2-4.*PR(1)*PR(2)))
 719         Z=(WREM2+PR(JT)-PR(JR)+PW12*(2.*Z-1.))/(2.*WREM2)
 720         PR(JR)=(PMQ(JR)+PARJ(32+MSTJ(11)))**2+(PX(JR)-PX(3))**2+
 721      &  (PY(JR)-PY(3))**2
 722         IF((1.-Z)*(WREM2-PR(JT)/Z).LT.PR(JR)) GOTO 810
 723       ENDIF
 724       GAM(3)=(1.-Z)*(GAM(JT)+PR(JT)/Z)
 725       DO 690 J=1,3
 726   690 IN(J)=IN(3*JT+J)
 727
 728 C...Stepping within or from 'low' string region easy.
 729       IF(IN(1)+1.EQ.IN(2).AND.Z*P(IN(1)+2,3)*P(IN(2)+2,3)*
 730      &P(IN(1),5)**2.GE.PR(JT)) THEN
 731         P(IN(JT)+2,4)=Z*P(IN(JT)+2,3)
 732         P(IN(JR)+2,4)=PR(JT)/(P(IN(JT)+2,4)*P(IN(1),5)**2)
 733         DO 700 J=1,4
 734   700   P(I,J)=(PX(JT)+PX(3))*P(IN(3),J)+(PY(JT)+PY(3))*P(IN(3)+1,J)
 735         GOTO 770
 736       ELSEIF(IN(1)+1.EQ.IN(2)) THEN
 737         P(IN(JR)+2,4)=P(IN(JR)+2,3)
 738         P(IN(JR)+2,JT)=1.
 739         IN(JR)=IN(JR)+4*JS
 740         IF(JS*IN(JR).GT.JS*IN(4*JR)) GOTO 550
 741         IF(FOUR(IN(1),IN(2)).LE.1E-2) THEN
 742           P(IN(JT)+2,4)=P(IN(JT)+2,3)
 743           P(IN(JT)+2,JT)=0.
 744           IN(JT)=IN(JT)+4*JS
 745         ENDIF
 746       ENDIF
 747
 748 C...Find new transverse directions (i.e. spacelike string vectors).
 749   710 IF(JS*IN(1).GT.JS*IN(3*JR+1).OR.JS*IN(2).GT.JS*IN(3*JR+2).OR.
 750      &IN(1).GT.IN(2)) GOTO 550
 751       IF(IN(1).NE.IN(3*JT+1).OR.IN(2).NE.IN(3*JT+2)) THEN
 752         DO 720 J=1,4
 753         DP(1,J)=P(IN(1),J)
 754         DP(2,J)=P(IN(2),J)
 755         DP(3,J)=0.
 756   720   DP(4,J)=0.
 757         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)
 758         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)
 759         DHC12=DFOUR(1,2)
 760         IF(DHC12.LE.1E-2) THEN
 761           P(IN(JT)+2,4)=P(IN(JT)+2,3)
 762           P(IN(JT)+2,JT)=0.
 763           IN(JT)=IN(JT)+4*JS
 764           GOTO 710
 765         ENDIF
 766         IN(3)=N+NR+4*NS+5
 767         DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4)
 768         DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4)
 769         DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4)
 770         IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1.
 771         IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1.
 772         IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1.
 773         IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1.
 774         DHCX1=DFOUR(3,1)/DHC12
 775         DHCX2=DFOUR(3,2)/DHC12
 776         DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12)
 777         DHCY1=DFOUR(4,1)/DHC12
 778         DHCY2=DFOUR(4,2)/DHC12
 779         DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12
 780         DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2)
 781         DO 730 J=1,4
 782         DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J))
 783         P(IN(3),J)=DP(3,J)
 784   730   P(IN(3)+1,J)=DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)-
 785      &  DHCYX*DP(3,J))
 786 C...Express pT with respect to new axes, if sensible.
 787         PXP=-(PX(3)*FOUR(IN(3*JT+3),IN(3))+PY(3)*
 788      &  FOUR(IN(3*JT+3)+1,IN(3)))
 789         PYP=-(PX(3)*FOUR(IN(3*JT+3),IN(3)+1)+PY(3)*
 790      &  FOUR(IN(3*JT+3)+1,IN(3)+1))
 791         IF(ABS(PXP**2+PYP**2-PX(3)**2-PY(3)**2).LT.0.01) THEN
 792           PX(3)=PXP
 793           PY(3)=PYP
 794         ENDIF
 795       ENDIF
 796
 797 C...Sum up known four-momentum. Gives coefficients for m2 expression.
 798       DO 750 J=1,4
 799       DHG(J)=0.
 800       P(I,J)=PX(JT)*P(IN(3*JT+3),J)+PY(JT)*P(IN(3*JT+3)+1,J)+
 801      &PX(3)*P(IN(3),J)+PY(3)*P(IN(3)+1,J)
 802       DO 740 IN1=IN(3*JT+1),IN(1)-4*JS,4*JS
 803   740 P(I,J)=P(I,J)+P(IN1+2,3)*P(IN1,J)
 804       DO 750 IN2=IN(3*JT+2),IN(2)-4*JS,4*JS
 805   750 P(I,J)=P(I,J)+P(IN2+2,3)*P(IN2,J)
 806       DHM(1)=FOUR(I,I)
 807       DHM(2)=2.*FOUR(I,IN(1))
 808       DHM(3)=2.*FOUR(I,IN(2))
 809       DHM(4)=2.*FOUR(IN(1),IN(2))
 810
 811 C...Find coefficients for Gamma expression.
 812       DO 760 IN2=IN(1)+1,IN(2),4
 813       DO 760 IN1=IN(1),IN2-1,4
 814       DHC=2.*FOUR(IN1,IN2)
 815       DHG(1)=DHG(1)+P(IN1+2,JT)*P(IN2+2,JT)*DHC
 816       IF(IN1.EQ.IN(1)) DHG(2)=DHG(2)-JS*P(IN2+2,JT)*DHC
 817       IF(IN2.EQ.IN(2)) DHG(3)=DHG(3)+JS*P(IN1+2,JT)*DHC
 818   760 IF(IN1.EQ.IN(1).AND.IN2.EQ.IN(2)) DHG(4)=DHG(4)-DHC
 819
 820 C...Solve (m2, Gamma) equation system for energies taken.
 821       DHS1=DHM(JR+1)*DHG(4)-DHM(4)*DHG(JR+1)
 822       IF(ABS(DHS1).LT.1E-4) GOTO 550
 823       DHS2=DHM(4)*(GAM(3)-DHG(1))-DHM(JT+1)*DHG(JR+1)-DHG(4)*
 824      &(P(I,5)**2-DHM(1))+DHG(JT+1)*DHM(JR+1)
 825       DHS3=DHM(JT+1)*(GAM(3)-DHG(1))-DHG(JT+1)*(P(I,5)**2-DHM(1))
 826       P(IN(JR)+2,4)=0.5*(SQRT(MAX(0D0,DHS2**2-4.*DHS1*DHS3))/ABS(DHS1)-
 827      &DHS2/DHS1)
 828       IF(DHM(JT+1)+DHM(4)*P(IN(JR)+2,4).LE.0.) GOTO 550
 829       P(IN(JT)+2,4)=(P(I,5)**2-DHM(1)-DHM(JR+1)*P(IN(JR)+2,4))/
 830      &(DHM(JT+1)+DHM(4)*P(IN(JR)+2,4))
 831
 832 C...Step to new region if necessary.
 833       IF(P(IN(JR)+2,4).GT.P(IN(JR)+2,3)) THEN
 834         P(IN(JR)+2,4)=P(IN(JR)+2,3)
 835         P(IN(JR)+2,JT)=1.
 836         IN(JR)=IN(JR)+4*JS
 837         IF(JS*IN(JR).GT.JS*IN(4*JR)) GOTO 550
 838         IF(FOUR(IN(1),IN(2)).LE.1E-2) THEN
 839           P(IN(JT)+2,4)=P(IN(JT)+2,3)
 840           P(IN(JT)+2,JT)=0.
 841           IN(JT)=IN(JT)+4*JS
 842         ENDIF
 843         GOTO 710
 844       ELSEIF(P(IN(JT)+2,4).GT.P(IN(JT)+2,3)) THEN
 845         P(IN(JT)+2,4)=P(IN(JT)+2,3)
 846         P(IN(JT)+2,JT)=0.
 847         IN(JT)=IN(JT)+4*JS
 848         GOTO 710
 849       ENDIF
 850
 851 C...Four-momentum of particle. Remaining quantities. Loop back.
 852   770 DO 780 J=1,4
 853       P(I,J)=P(I,J)+P(IN(1)+2,4)*P(IN(1),J)+P(IN(2)+2,4)*P(IN(2),J)
 854   780 P(N+NRS,J)=P(N+NRS,J)-P(I,J)
 855       IF(P(I,4).LE.0.) GOTO 550
 856       KFL(JT)=-KFL(3)
 857       PMQ(JT)=PMQ(3)
 858       PX(JT)=-PX(3)
 859       PY(JT)=-PY(3)
 860       GAM(JT)=GAM(3)
 861       IF(IN(3).NE.IN(3*JT+3)) THEN
 862         DO 790 J=1,4
 863         P(IN(3*JT+3),J)=P(IN(3),J)
 864   790   P(IN(3*JT+3)+1,J)=P(IN(3)+1,J)
 865       ENDIF
 866       DO 800 JQ=1,2
 867       IN(3*JT+JQ)=IN(JQ)
 868       P(IN(JQ)+2,3)=P(IN(JQ)+2,3)-P(IN(JQ)+2,4)
 869   800 P(IN(JQ)+2,JT)=P(IN(JQ)+2,JT)-JS*(3-2*JQ)*P(IN(JQ)+2,4)
 870       GOTO 670
 871
 872 C...Final hadron: side, flavour, hadron, mass.
 873   810 I=I+1
 874       K(I,1)=1
 875       K(I,3)=IE(JR)
 876       K(I,4)=0
 877       K(I,5)=0
 878       CALL LUKFDI_HIJING(KFL(JR),-KFL(3),KFLDMP,K(I,2))
 879       IF(K(I,2).EQ.0) GOTO 550
 880       P(I,5)=ULMASS_HIJING(K(I,2))
 881       PR(JR)=P(I,5)**2+(PX(JR)-PX(3))**2+(PY(JR)-PY(3))**2
 882
 883 C...Final two hadrons: find common setup of four-vectors.
 884       JQ=1
 885       IF(P(IN(4)+2,3)*P(IN(5)+2,3)*FOUR(IN(4),IN(5)).LT.P(IN(7),3)*
 886      &P(IN(8),3)*FOUR(IN(7),IN(8))) JQ=2
 887       DHC12=FOUR(IN(3*JQ+1),IN(3*JQ+2))
 888       DHR1=FOUR(N+NRS,IN(3*JQ+2))/DHC12
 889       DHR2=FOUR(N+NRS,IN(3*JQ+1))/DHC12
 890       IF(IN(4).NE.IN(7).OR.IN(5).NE.IN(8)) THEN
 891         PX(3-JQ)=-FOUR(N+NRS,IN(3*JQ+3))-PX(JQ)
 892         PY(3-JQ)=-FOUR(N+NRS,IN(3*JQ+3)+1)-PY(JQ)
 893         PR(3-JQ)=P(I+(JT+JQ-3)**2-1,5)**2+(PX(3-JQ)+(2*JQ-3)*JS*
 894      &  PX(3))**2+(PY(3-JQ)+(2*JQ-3)*JS*PY(3))**2
 895       ENDIF
 896
 897 C...Solve kinematics for final two hadrons, if possible.
 898       WREM2=WREM2+(PX(1)+PX(2))**2+(PY(1)+PY(2))**2
 899       FD=(SQRT(PR(1))+SQRT(PR(2)))/SQRT(WREM2)
 900       IF(MJU(1)+MJU(2).NE.0.AND.I.EQ.ISAV+2.AND.FD.GE.1.) GOTO 180
 901       IF(FD.GE.1.) GOTO 550
 902       FA=WREM2+PR(JT)-PR(JR)
 903       IF(MSTJ(11).EQ.2) PREV=0.5*FD**PARJ(37+MSTJ(11))
 904       IF(MSTJ(11).NE.2) PREV=0.5*EXP(MAX(-100.,LOG(FD)*
 905      &PARJ(37+MSTJ(11))*(PR(1)+PR(2))**2))
 906       FB=SIGN(SQRT(MAX(0.,FA**2-4.*WREM2*PR(JT))),JS*(RLU_HIJING(0)-PREV
 907      $     ))
 908       KFL1A=IABS(KFL(1))
 909       KFL2A=IABS(KFL(2))
 910       IF(MAX(MOD(KFL1A,10),MOD(KFL1A/1000,10),MOD(KFL2A,10),
 911      &MOD(KFL2A/1000,10)).GE.6) FB=SIGN(SQRT(MAX(0.,FA**2-
 912      &4.*WREM2*PR(JT))),FLOAT(JS))
 913       DO 820 J=1,4
 914       P(I-1,J)=(PX(JT)+PX(3))*P(IN(3*JQ+3),J)+(PY(JT)+PY(3))*
 915      &P(IN(3*JQ+3)+1,J)+0.5*(DHR1*(FA+FB)*P(IN(3*JQ+1),J)+
 916      &DHR2*(FA-FB)*P(IN(3*JQ+2),J))/WREM2
 917   820 P(I,J)=P(N+NRS,J)-P(I-1,J)
 918
 919 C...Mark jets as fragmented and give daughter pointers.
 920       N=I-NRS+1
 921       DO 830 I=NSAV+1,NSAV+NP
 922       IM=K(I,3)
 923       K(IM,1)=K(IM,1)+10
 924       IF(MSTU(16).NE.2) THEN
 925         K(IM,4)=NSAV+1
 926         K(IM,5)=NSAV+1
 927       ELSE
 928         K(IM,4)=NSAV+2
 929         K(IM,5)=N
 930       ENDIF
 931   830 CONTINUE
 932
 933 C...Document string system. Move up particles.
 934       NSAV=NSAV+1
 935       K(NSAV,1)=11
 936       K(NSAV,2)=92
 937       K(NSAV,3)=IP
 938       K(NSAV,4)=NSAV+1
 939       K(NSAV,5)=N
 940       DO 840 J=1,4
 941       P(NSAV,J)=DPS(J)
 942   840 V(NSAV,J)=V(IP,J)
 943       P(NSAV,5)=SQRT(MAX(0D0,DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2))
 944       V(NSAV,5)=0.
 945       DO 850 I=NSAV+1,N
 946       DO 850 J=1,5
 947       K(I,J)=K(I+NRS-1,J)
 948       P(I,J)=P(I+NRS-1,J)
 949   850 V(I,J)=0.
 950
 951 C...Order particles in rank along the chain. Update mother pointer.
 952       DO 860 I=NSAV+1,N
 953       DO 860 J=1,5
 954       K(I-NSAV+N,J)=K(I,J)
 955   860 P(I-NSAV+N,J)=P(I,J)
 956       I1=NSAV
 957       DO 880 I=N+1,2*N-NSAV
 958       IF(K(I,3).NE.IE(1)) GOTO 880
 959       I1=I1+1
 960       DO 870 J=1,5
 961       K(I1,J)=K(I,J)
 962   870 P(I1,J)=P(I,J)
 963       IF(MSTU(16).NE.2) K(I1,3)=NSAV
 964   880 CONTINUE
 965       DO 900 I=2*N-NSAV,N+1,-1
 966       IF(K(I,3).EQ.IE(1)) GOTO 900
 967       I1=I1+1
 968       DO 890 J=1,5
 969       K(I1,J)=K(I,J)
 970   890 P(I1,J)=P(I,J)
 971       IF(MSTU(16).NE.2) K(I1,3)=NSAV
 972   900 CONTINUE
 973
 974 C...Boost back particle system. Set production vertices.
 975       CALL LUDBRB_HIJING(NSAV+1,N,0.,0.,DPS(1)/DPS(4),DPS(2)/DPS(4),
 976      &DPS(3)/DPS(4))
 977       DO 910 I=NSAV+1,N
 978       DO 910 J=1,4
 979   910 V(I,J)=V(IP,J)
 980
 981       RETURN
 982       END