PHOS/shaker/lustrf.f

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