PYTHIA/jetset/lustrf.F

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