hijing/src/lustrf.f

0001
0002 C*********************************************************************
0003
0004       SUBROUTINE LUSTRF(IP)
0005 C...Purpose: to handle the fragmentation of an arbitrary colour singlet
0006 C...jet system according to the Lund string fragmentation model.
0007       IMPLICIT DOUBLE PRECISION(D)
0008       COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
0009       SAVE /LUJETS/
0010       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0011       SAVE /LUDAT1/
0012       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
0013       SAVE /LUDAT2/
0014       DIMENSION DPS(5),KFL(3),PMQ(3),PX(3),PY(3),GAM(3),IE(2),PR(2),
0015      &IN(9),DHM(4),DHG(4),DP(5,5),IRANK(2),MJU(4),IJU(3),PJU(5,5),
0016      &TJU(5),KFJH(2),NJS(2),KFJS(2),PJS(4,5)
0017
0018 C...Function: four-product of two vectors.
0019       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)
0020       DFOUR(I,J)=DP(I,4)*DP(J,4)-DP(I,1)*DP(J,1)-DP(I,2)*DP(J,2)-
0021      &DP(I,3)*DP(J,3)
0022
0023 C...Reset counters. Identify parton system.
0024       MSTJ(91)=0
0025       NSAV=N
0026       NP=0
0027       KQSUM=0
0028       DO 100 J=1,5
0029   100 DPS(J)=0.
0030       MJU(1)=0
0031       MJU(2)=0
0032       I=IP-1
0033   110 I=I+1
0034       IF(I.GT.MIN(N,MSTU(4)-MSTU(32))) THEN
0035         CALL LUERRM(12,'(LUSTRF:) failed to reconstruct jet system')
0036         IF(MSTU(21).GE.1) RETURN
0037       ENDIF
0038       IF(K(I,1).NE.1.AND.K(I,1).NE.2.AND.K(I,1).NE.41) GOTO 110
0039       KC=LUCOMP(K(I,2))
0040       IF(KC.EQ.0) GOTO 110
0041       KQ=KCHG(KC,2)*ISIGN(1,K(I,2))
0042       IF(KQ.EQ.0) GOTO 110
0043       IF(N+5*NP+11.GT.MSTU(4)-MSTU(32)-5) THEN
0044         CALL LUERRM(11,'(LUSTRF:) no more memory left in LUJETS')
0045         IF(MSTU(21).GE.1) RETURN
0046       ENDIF
0047
0048 C...Take copy of partons to be considered. Check flavour sum.
0049       NP=NP+1
0050       DO 120 J=1,5
0051       K(N+NP,J)=K(I,J)
0052       P(N+NP,J)=P(I,J)
0053   120 DPS(J)=DPS(J)+P(I,J)
0054       K(N+NP,3)=I
0055       IF(P(N+NP,4)**2.LT.P(N+NP,1)**2+P(N+NP,2)**2+P(N+NP,3)**2) THEN
0056         P(N+NP,4)=SQRT(P(N+NP,1)**2+P(N+NP,2)**2+P(N+NP,3)**2+
0057      &  P(N+NP,5)**2)
0058         DPS(4)=DPS(4)+MAX(0.,P(N+NP,4)-P(I,4))
0059       ENDIF
0060       IF(KQ.NE.2) KQSUM=KQSUM+KQ
0061       IF(K(I,1).EQ.41) THEN
0062         KQSUM=KQSUM+2*KQ
0063         IF(KQSUM.EQ.KQ) MJU(1)=N+NP
0064         IF(KQSUM.NE.KQ) MJU(2)=N+NP
0065       ENDIF
0066       IF(K(I,1).EQ.2.OR.K(I,1).EQ.41) GOTO 110
0067       IF(KQSUM.NE.0) THEN
0068         CALL LUERRM(12,'(LUSTRF:) unphysical flavour combination')
0069         IF(MSTU(21).GE.1) RETURN
0070       ENDIF
0071
0072 C...Boost copied system to CM frame (for better numerical precision).
0073       CALL LUDBRB(N+1,N+NP,0.,0.,-DPS(1)/DPS(4),-DPS(2)/DPS(4),
0074      &-DPS(3)/DPS(4))
0075
0076 C...Search for very nearby partons that may be recombined.
0077       NTRYR=0
0078       PARU12=PARU(12)
0079       PARU13=PARU(13)
0080       MJU(3)=MJU(1)
0081       MJU(4)=MJU(2)
0082       NR=NP
0083   130 IF(NR.GE.3) THEN
0084         PDRMIN=2.*PARU12
0085         DO 140 I=N+1,N+NR
0086         IF(I.EQ.N+NR.AND.IABS(K(N+1,2)).NE.21) GOTO 140
0087         I1=I+1
0088         IF(I.EQ.N+NR) I1=N+1
0089         IF(K(I,1).EQ.41.OR.K(I1,1).EQ.41) GOTO 140
0090         IF(MJU(1).NE.0.AND.I1.LT.MJU(1).AND.IABS(K(I1,2)).NE.21)
0091      &  GOTO 140
0092         IF(MJU(2).NE.0.AND.I.GT.MJU(2).AND.IABS(K(I,2)).NE.21) GOTO 140
0093         PAP=SQRT((P(I,1)**2+P(I,2)**2+P(I,3)**2)*(P(I1,1)**2+
0094      &  P(I1,2)**2+P(I1,3)**2))
0095         PVP=P(I,1)*P(I1,1)+P(I,2)*P(I1,2)+P(I,3)*P(I1,3)
0096         PDR=4.*(PAP-PVP)**2/(PARU13**2*PAP+2.*(PAP-PVP))
0097         IF(PDR.LT.PDRMIN) THEN
0098           IR=I
0099           PDRMIN=PDR
0100         ENDIF
0101   140   CONTINUE
0102
0103 C...Recombine very nearby partons to avoid machine precision problems.
0104         IF(PDRMIN.LT.PARU12.AND.IR.EQ.N+NR) THEN
0105           DO 150 J=1,4
0106   150     P(N+1,J)=P(N+1,J)+P(N+NR,J)
0107           P(N+1,5)=SQRT(MAX(0.,P(N+1,4)**2-P(N+1,1)**2-P(N+1,2)**2-
0108      &    P(N+1,3)**2))
0109           NR=NR-1
0110           GOTO 130
0111         ELSEIF(PDRMIN.LT.PARU12) THEN
0112           DO 160 J=1,4
0113   160     P(IR,J)=P(IR,J)+P(IR+1,J)
0114           P(IR,5)=SQRT(MAX(0.,P(IR,4)**2-P(IR,1)**2-P(IR,2)**2-
0115      &    P(IR,3)**2))
0116           DO 170 I=IR+1,N+NR-1
0117           K(I,2)=K(I+1,2)
0118           DO 170 J=1,5
0119   170     P(I,J)=P(I+1,J)
0120           IF(IR.EQ.N+NR-1) K(IR,2)=K(N+NR,2)
0121           NR=NR-1
0122           IF(MJU(1).GT.IR) MJU(1)=MJU(1)-1
0123           IF(MJU(2).GT.IR) MJU(2)=MJU(2)-1
0124           GOTO 130
0125         ENDIF
0126       ENDIF
0127       NTRYR=NTRYR+1
0128
0129 C...Reset particle counter. Skip ahead if no junctions are present;
0130 C...this is usually the case!
0131       NRS=MAX(5*NR+11,NP)
0132       NTRY=0
0133   180 NTRY=NTRY+1
0134       IF(NTRY.GT.100.AND.NTRYR.LE.4) THEN
0135         PARU12=4.*PARU12
0136         PARU13=2.*PARU13
0137         GOTO 130
0138       ELSEIF(NTRY.GT.100) THEN
0139         CALL LUERRM(14,'(LUSTRF:) caught in infinite loop')
0140         IF(MSTU(21).GE.1) RETURN
0141       ENDIF
0142       I=N+NRS
0143       IF(MJU(1).EQ.0.AND.MJU(2).EQ.0) GOTO 500
0144       DO 490 JT=1,2
0145       NJS(JT)=0
0146       IF(MJU(JT).EQ.0) GOTO 490
0147       JS=3-2*JT
0148
0149 C...Find and sum up momentum on three sides of junction. Check flavours.
0150       DO 190 IU=1,3
0151       IJU(IU)=0
0152       DO 190 J=1,5
0153   190 PJU(IU,J)=0.
0154       IU=0
0155       DO 200 I1=N+1+(JT-1)*(NR-1),N+NR+(JT-1)*(1-NR),JS
0156       IF(K(I1,2).NE.21.AND.IU.LE.2) THEN
0157         IU=IU+1
0158         IJU(IU)=I1
0159       ENDIF
0160       DO 200 J=1,4
0161   200 PJU(IU,J)=PJU(IU,J)+P(I1,J)
0162       DO 210 IU=1,3
0163   210 PJU(IU,5)=SQRT(PJU(IU,1)**2+PJU(IU,2)**2+PJU(IU,3)**2)
0164       IF(K(IJU(3),2)/100.NE.10*K(IJU(1),2)+K(IJU(2),2).AND.
0165      &K(IJU(3),2)/100.NE.10*K(IJU(2),2)+K(IJU(1),2)) THEN
0166         CALL LUERRM(12,'(LUSTRF:) unphysical flavour combination')
0167         IF(MSTU(21).GE.1) RETURN
0168       ENDIF
0169
0170 C...Calculate (approximate) boost to rest frame of junction.
0171       T12=(PJU(1,1)*PJU(2,1)+PJU(1,2)*PJU(2,2)+PJU(1,3)*PJU(2,3))/
0172      &(PJU(1,5)*PJU(2,5))
0173       T13=(PJU(1,1)*PJU(3,1)+PJU(1,2)*PJU(3,2)+PJU(1,3)*PJU(3,3))/
0174      &(PJU(1,5)*PJU(3,5))
0175       T23=(PJU(2,1)*PJU(3,1)+PJU(2,2)*PJU(3,2)+PJU(2,3)*PJU(3,3))/
0176      &(PJU(2,5)*PJU(3,5))
0177       T11=SQRT((2./3.)*(1.-T12)*(1.-T13)/(1.-T23))
0178       T22=SQRT((2./3.)*(1.-T12)*(1.-T23)/(1.-T13))
0179       TSQ=SQRT((2.*T11*T22+T12-1.)*(1.+T12))
0180       T1F=(TSQ-T22*(1.+T12))/(1.-T12**2)
0181       T2F=(TSQ-T11*(1.+T12))/(1.-T12**2)
0182       DO 220 J=1,3
0183   220 TJU(J)=-(T1F*PJU(1,J)/PJU(1,5)+T2F*PJU(2,J)/PJU(2,5))
0184       TJU(4)=SQRT(1.+TJU(1)**2+TJU(2)**2+TJU(3)**2)
0185       DO 230 IU=1,3
0186   230 PJU(IU,5)=TJU(4)*PJU(IU,4)-TJU(1)*PJU(IU,1)-TJU(2)*PJU(IU,2)-
0187      &TJU(3)*PJU(IU,3)
0188
0189 C...Put junction at rest if motion could give inconsistencies.
0190       IF(PJU(1,5)+PJU(2,5).GT.PJU(1,4)+PJU(2,4)) THEN
0191         DO 240 J=1,3
0192   240   TJU(J)=0.
0193         TJU(4)=1.
0194         PJU(1,5)=PJU(1,4)
0195         PJU(2,5)=PJU(2,4)
0196         PJU(3,5)=PJU(3,4)
0197       ENDIF
0198
0199 C...Start preparing for fragmentation of two strings from junction.
0200       ISTA=I
0201       DO 470 IU=1,2
0202       NS=IJU(IU+1)-IJU(IU)
0203
0204 C...Junction strings: find longitudinal string directions.
0205       DO 260 IS=1,NS
0206       IS1=IJU(IU)+IS-1
0207       IS2=IJU(IU)+IS
0208       DO 250 J=1,5
0209       DP(1,J)=0.5*P(IS1,J)
0210       IF(IS.EQ.1) DP(1,J)=P(IS1,J)
0211       DP(2,J)=0.5*P(IS2,J)
0212   250 IF(IS.EQ.NS) DP(2,J)=-PJU(IU,J)
0213       IF(IS.EQ.NS) DP(2,4)=SQRT(PJU(IU,1)**2+PJU(IU,2)**2+PJU(IU,3)**2)
0214       IF(IS.EQ.NS) DP(2,5)=0.
0215       DP(3,5)=DFOUR(1,1)
0216       DP(4,5)=DFOUR(2,2)
0217       DHKC=DFOUR(1,2)
0218       IF(DP(3,5)+2.*DHKC+DP(4,5).LE.0.) THEN
0219         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)
0220         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)
0221         DP(3,5)=0D0
0222         DP(4,5)=0D0
0223         DHKC=DFOUR(1,2)
0224       ENDIF
0225       DHKS=SQRT(DHKC**2-DP(3,5)*DP(4,5))
0226       DHK1=0.5*((DP(4,5)+DHKC)/DHKS-1.)
0227       DHK2=0.5*((DP(3,5)+DHKC)/DHKS-1.)
0228       IN1=N+NR+4*IS-3
0229       P(IN1,5)=SQRT(DP(3,5)+2.*DHKC+DP(4,5))
0230       DO 260 J=1,4
0231       P(IN1,J)=(1.+DHK1)*DP(1,J)-DHK2*DP(2,J)
0232   260 P(IN1+1,J)=(1.+DHK2)*DP(2,J)-DHK1*DP(1,J)
0233
0234 C...Junction strings: initialize flavour, momentum and starting pos.
0235       ISAV=I
0236   270 NTRY=NTRY+1
0237       IF(NTRY.GT.100.AND.NTRYR.LE.4) THEN
0238         PARU12=4.*PARU12
0239         PARU13=2.*PARU13
0240         GOTO 130
0241       ELSEIF(NTRY.GT.100) THEN
0242         CALL LUERRM(14,'(LUSTRF:) caught in infinite loop')
0243         IF(MSTU(21).GE.1) RETURN
0244       ENDIF
0245       I=ISAV
0246       IRANKJ=0
0247       IE(1)=K(N+1+(JT/2)*(NP-1),3)
0248       IN(4)=N+NR+1
0249       IN(5)=IN(4)+1
0250       IN(6)=N+NR+4*NS+1
0251       DO 280 JQ=1,2
0252       DO 280 IN1=N+NR+2+JQ,N+NR+4*NS-2+JQ,4
0253       P(IN1,1)=2-JQ
0254       P(IN1,2)=JQ-1
0255   280 P(IN1,3)=1.
0256       KFL(1)=K(IJU(IU),2)
0257       PX(1)=0.
0258       PY(1)=0.
0259       GAM(1)=0.
0260       DO 290 J=1,5
0261   290 PJU(IU+3,J)=0.
0262
0263 C...Junction strings: find initial transverse directions.
0264       DO 300 J=1,4
0265       DP(1,J)=P(IN(4),J)
0266       DP(2,J)=P(IN(4)+1,J)
0267       DP(3,J)=0.
0268   300 DP(4,J)=0.
0269       DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)
0270       DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)
0271       DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4)
0272       DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4)
0273       DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4)
0274       IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1.
0275       IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1.
0276       IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1.
0277       IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1.
0278       DHC12=DFOUR(1,2)
0279       DHCX1=DFOUR(3,1)/DHC12
0280       DHCX2=DFOUR(3,2)/DHC12
0281       DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12)
0282       DHCY1=DFOUR(4,1)/DHC12
0283       DHCY2=DFOUR(4,2)/DHC12
0284       DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12
0285       DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2)
0286       DO 310 J=1,4
0287       DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J))
0288       P(IN(6),J)=DP(3,J)
0289   310 P(IN(6)+1,J)=DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)-
0290      &DHCYX*DP(3,J))
0291
0292 C...Junction strings: produce new particle, origin.
0293   320 I=I+1
0294       IF(2*I-NSAV.GE.MSTU(4)-MSTU(32)-5) THEN
0295         CALL LUERRM(11,'(LUSTRF:) no more memory left in LUJETS')
0296         IF(MSTU(21).GE.1) RETURN
0297       ENDIF
0298       IRANKJ=IRANKJ+1
0299       K(I,1)=1
0300       K(I,3)=IE(1)
0301       K(I,4)=0
0302       K(I,5)=0
0303
0304 C...Junction strings: generate flavour, hadron, pT, z and Gamma.
0305   330 CALL LUKFDI(KFL(1),0,KFL(3),K(I,2))
0306       IF(K(I,2).EQ.0) GOTO 270
0307       IF(MSTJ(12).GE.3.AND.IRANKJ.EQ.1.AND.IABS(KFL(1)).LE.10.AND.
0308      &IABS(KFL(3)).GT.10) THEN
0309         IF(RLU(0).GT.PARJ(19)) GOTO 330
0310       ENDIF
0311       P(I,5)=ULMASS(K(I,2))
0312       CALL LUPTDI(KFL(1),PX(3),PY(3))
0313       PR(1)=P(I,5)**2+(PX(1)+PX(3))**2+(PY(1)+PY(3))**2
0314       CALL LUZDIS(KFL(1),KFL(3),PR(1),Z)
0315       GAM(3)=(1.-Z)*(GAM(1)+PR(1)/Z)
0316       DO 340 J=1,3
0317   340 IN(J)=IN(3+J)
0318
0319 C...Junction strings: stepping within or from 'low' string region easy.
0320       IF(IN(1)+1.EQ.IN(2).AND.Z*P(IN(1)+2,3)*P(IN(2)+2,3)*
0321      &P(IN(1),5)**2.GE.PR(1)) THEN
0322         P(IN(1)+2,4)=Z*P(IN(1)+2,3)
0323         P(IN(2)+2,4)=PR(1)/(P(IN(1)+2,4)*P(IN(1),5)**2)
0324         DO 350 J=1,4
0325   350   P(I,J)=(PX(1)+PX(3))*P(IN(3),J)+(PY(1)+PY(3))*P(IN(3)+1,J)
0326         GOTO 420
0327       ELSEIF(IN(1)+1.EQ.IN(2)) THEN
0328         P(IN(2)+2,4)=P(IN(2)+2,3)
0329         P(IN(2)+2,1)=1.
0330         IN(2)=IN(2)+4
0331         IF(IN(2).GT.N+NR+4*NS) GOTO 270
0332         IF(FOUR(IN(1),IN(2)).LE.1E-2) THEN
0333           P(IN(1)+2,4)=P(IN(1)+2,3)
0334           P(IN(1)+2,1)=0.
0335           IN(1)=IN(1)+4
0336         ENDIF
0337       ENDIF
0338
0339 C...Junction strings: find new transverse directions.
0340   360 IF(IN(1).GT.N+NR+4*NS.OR.IN(2).GT.N+NR+4*NS.OR.
0341      &IN(1).GT.IN(2)) GOTO 270
0342       IF(IN(1).NE.IN(4).OR.IN(2).NE.IN(5)) THEN
0343         DO 370 J=1,4
0344         DP(1,J)=P(IN(1),J)
0345         DP(2,J)=P(IN(2),J)
0346         DP(3,J)=0.
0347   370   DP(4,J)=0.
0348         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)
0349         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)
0350         DHC12=DFOUR(1,2)
0351         IF(DHC12.LE.1E-2) THEN
0352           P(IN(1)+2,4)=P(IN(1)+2,3)
0353           P(IN(1)+2,1)=0.
0354           IN(1)=IN(1)+4
0355           GOTO 360
0356         ENDIF
0357         IN(3)=N+NR+4*NS+5
0358         DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4)
0359         DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4)
0360         DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4)
0361         IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1.
0362         IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1.
0363         IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1.
0364         IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1.
0365         DHCX1=DFOUR(3,1)/DHC12
0366         DHCX2=DFOUR(3,2)/DHC12
0367         DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12)
0368         DHCY1=DFOUR(4,1)/DHC12
0369         DHCY2=DFOUR(4,2)/DHC12
0370         DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12
0371         DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2)
0372         DO 380 J=1,4
0373         DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J))
0374         P(IN(3),J)=DP(3,J)
0375   380   P(IN(3)+1,J)=DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)-
0376      &  DHCYX*DP(3,J))
0377 C...Express pT with respect to new axes, if sensible.
0378         PXP=-(PX(3)*FOUR(IN(6),IN(3))+PY(3)*FOUR(IN(6)+1,IN(3)))
0379         PYP=-(PX(3)*FOUR(IN(6),IN(3)+1)+PY(3)*FOUR(IN(6)+1,IN(3)+1))
0380         IF(ABS(PXP**2+PYP**2-PX(3)**2-PY(3)**2).LT.0.01) THEN
0381           PX(3)=PXP
0382           PY(3)=PYP
0383         ENDIF
0384       ENDIF
0385
0386 C...Junction strings: sum up known four-momentum, coefficients for m2.
0387       DO 400 J=1,4
0388       DHG(J)=0.
0389       P(I,J)=PX(1)*P(IN(6),J)+PY(1)*P(IN(6)+1,J)+PX(3)*P(IN(3),J)+
0390      &PY(3)*P(IN(3)+1,J)
0391       DO 390 IN1=IN(4),IN(1)-4,4
0392   390 P(I,J)=P(I,J)+P(IN1+2,3)*P(IN1,J)
0393       DO 400 IN2=IN(5),IN(2)-4,4
0394   400 P(I,J)=P(I,J)+P(IN2+2,3)*P(IN2,J)
0395       DHM(1)=FOUR(I,I)
0396       DHM(2)=2.*FOUR(I,IN(1))
0397       DHM(3)=2.*FOUR(I,IN(2))
0398       DHM(4)=2.*FOUR(IN(1),IN(2))
0399
0400 C...Junction strings: find coefficients for Gamma expression.
0401       DO 410 IN2=IN(1)+1,IN(2),4
0402       DO 410 IN1=IN(1),IN2-1,4
0403       DHC=2.*FOUR(IN1,IN2)
0404       DHG(1)=DHG(1)+P(IN1+2,1)*P(IN2+2,1)*DHC
0405       IF(IN1.EQ.IN(1)) DHG(2)=DHG(2)-P(IN2+2,1)*DHC
0406       IF(IN2.EQ.IN(2)) DHG(3)=DHG(3)+P(IN1+2,1)*DHC
0407   410 IF(IN1.EQ.IN(1).AND.IN2.EQ.IN(2)) DHG(4)=DHG(4)-DHC
0408
0409 C...Junction strings: solve (m2, Gamma) equation system for energies.
0410       DHS1=DHM(3)*DHG(4)-DHM(4)*DHG(3)
0411       IF(ABS(DHS1).LT.1E-4) GOTO 270
0412       DHS2=DHM(4)*(GAM(3)-DHG(1))-DHM(2)*DHG(3)-DHG(4)*
0413      &(P(I,5)**2-DHM(1))+DHG(2)*DHM(3)
0414       DHS3=DHM(2)*(GAM(3)-DHG(1))-DHG(2)*(P(I,5)**2-DHM(1))
0415       P(IN(2)+2,4)=0.5*(SQRT(MAX(0D0,DHS2**2-4.*DHS1*DHS3))/ABS(DHS1)-
0416      &DHS2/DHS1)
0417       IF(DHM(2)+DHM(4)*P(IN(2)+2,4).LE.0.) GOTO 270
0418       P(IN(1)+2,4)=(P(I,5)**2-DHM(1)-DHM(3)*P(IN(2)+2,4))/
0419      &(DHM(2)+DHM(4)*P(IN(2)+2,4))
0420
0421 C...Junction strings: step to new region if necessary.
0422       IF(P(IN(2)+2,4).GT.P(IN(2)+2,3)) THEN
0423         P(IN(2)+2,4)=P(IN(2)+2,3)
0424         P(IN(2)+2,1)=1.
0425         IN(2)=IN(2)+4
0426         IF(IN(2).GT.N+NR+4*NS) GOTO 270
0427         IF(FOUR(IN(1),IN(2)).LE.1E-2) THEN
0428           P(IN(1)+2,4)=P(IN(1)+2,3)
0429           P(IN(1)+2,1)=0.
0430           IN(1)=IN(1)+4
0431         ENDIF
0432         GOTO 360
0433       ELSEIF(P(IN(1)+2,4).GT.P(IN(1)+2,3)) THEN
0434         P(IN(1)+2,4)=P(IN(1)+2,3)
0435         P(IN(1)+2,1)=0.
0436         IN(1)=IN(1)+JS
0437         GOTO 710
0438       ENDIF
0439
0440 C...Junction strings: particle four-momentum, remainder, loop back.
0441   420 DO 430 J=1,4
0442       P(I,J)=P(I,J)+P(IN(1)+2,4)*P(IN(1),J)+P(IN(2)+2,4)*P(IN(2),J)
0443   430 PJU(IU+3,J)=PJU(IU+3,J)+P(I,J)
0444       IF(P(I,4).LE.0.) GOTO 270
0445       PJU(IU+3,5)=TJU(4)*PJU(IU+3,4)-TJU(1)*PJU(IU+3,1)-
0446      &TJU(2)*PJU(IU+3,2)-TJU(3)*PJU(IU+3,3)
0447       IF(PJU(IU+3,5).LT.PJU(IU,5)) THEN
0448         KFL(1)=-KFL(3)
0449         PX(1)=-PX(3)
0450         PY(1)=-PY(3)
0451         GAM(1)=GAM(3)
0452         IF(IN(3).NE.IN(6)) THEN
0453           DO 440 J=1,4
0454           P(IN(6),J)=P(IN(3),J)
0455   440     P(IN(6)+1,J)=P(IN(3)+1,J)
0456         ENDIF
0457         DO 450 JQ=1,2
0458         IN(3+JQ)=IN(JQ)
0459         P(IN(JQ)+2,3)=P(IN(JQ)+2,3)-P(IN(JQ)+2,4)
0460   450   P(IN(JQ)+2,1)=P(IN(JQ)+2,1)-(3-2*JQ)*P(IN(JQ)+2,4)
0461         GOTO 320
0462       ENDIF
0463
0464 C...Junction strings: save quantities left after each string.
0465       IF(IABS(KFL(1)).GT.10) GOTO 270
0466       I=I-1
0467       KFJH(IU)=KFL(1)
0468       DO 460 J=1,4
0469   460 PJU(IU+3,J)=PJU(IU+3,J)-P(I+1,J)
0470   470 CONTINUE
0471
0472 C...Junction strings: put together to new effective string endpoint.
0473       NJS(JT)=I-ISTA
0474       KFJS(JT)=K(K(MJU(JT+2),3),2)
0475       KFLS=2*INT(RLU(0)+3.*PARJ(4)/(1.+3.*PARJ(4)))+1
0476       IF(KFJH(1).EQ.KFJH(2)) KFLS=3
0477       IF(ISTA.NE.I) KFJS(JT)=ISIGN(1000*MAX(IABS(KFJH(1)),
0478      &IABS(KFJH(2)))+100*MIN(IABS(KFJH(1)),IABS(KFJH(2)))+
0479      &KFLS,KFJH(1))
0480       DO 480 J=1,4
0481       PJS(JT,J)=PJU(1,J)+PJU(2,J)+P(MJU(JT),J)
0482   480 PJS(JT+2,J)=PJU(4,J)+PJU(5,J)
0483       PJS(JT,5)=SQRT(MAX(0.,PJS(JT,4)**2-PJS(JT,1)**2-PJS(JT,2)**2-
0484      &PJS(JT,3)**2))
0485   490 CONTINUE
0486
0487 C...Open versus closed strings. Choose breakup region for latter.
0488   500 IF(MJU(1).NE.0.AND.MJU(2).NE.0) THEN
0489         NS=MJU(2)-MJU(1)
0490         NB=MJU(1)-N
0491       ELSEIF(MJU(1).NE.0) THEN
0492         NS=N+NR-MJU(1)
0493         NB=MJU(1)-N
0494       ELSEIF(MJU(2).NE.0) THEN
0495         NS=MJU(2)-N
0496         NB=1
0497       ELSEIF(IABS(K(N+1,2)).NE.21) THEN
0498         NS=NR-1
0499         NB=1
0500       ELSE
0501         NS=NR+1
0502         W2SUM=0.
0503         DO 510 IS=1,NR
0504         P(N+NR+IS,1)=0.5*FOUR(N+IS,N+IS+1-NR*(IS/NR))
0505   510   W2SUM=W2SUM+P(N+NR+IS,1)
0506         W2RAN=RLU(0)*W2SUM
0507         NB=0
0508   520   NB=NB+1
0509         W2SUM=W2SUM-P(N+NR+NB,1)
0510         IF(W2SUM.GT.W2RAN.AND.NB.LT.NR) GOTO 520
0511       ENDIF
0512
0513 C...Find longitudinal string directions (i.e. lightlike four-vectors).
0514       DO 540 IS=1,NS
0515       IS1=N+IS+NB-1-NR*((IS+NB-2)/NR)
0516       IS2=N+IS+NB-NR*((IS+NB-1)/NR)
0517       DO 530 J=1,5
0518       DP(1,J)=P(IS1,J)
0519       IF(IABS(K(IS1,2)).EQ.21) DP(1,J)=0.5*DP(1,J)
0520       IF(IS1.EQ.MJU(1)) DP(1,J)=PJS(1,J)-PJS(3,J)
0521       DP(2,J)=P(IS2,J)
0522       IF(IABS(K(IS2,2)).EQ.21) DP(2,J)=0.5*DP(2,J)
0523   530 IF(IS2.EQ.MJU(2)) DP(2,J)=PJS(2,J)-PJS(4,J)
0524       DP(3,5)=DFOUR(1,1)
0525       DP(4,5)=DFOUR(2,2)
0526       DHKC=DFOUR(1,2)
0527       IF(DP(3,5)+2.*DHKC+DP(4,5).LE.0.) THEN
0528         DP(3,5)=DP(1,5)**2
0529         DP(4,5)=DP(2,5)**2
0530         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2+DP(1,5)**2)
0531         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2+DP(2,5)**2)
0532         DHKC=DFOUR(1,2)
0533       ENDIF
0534       DHKS=SQRT(DHKC**2-DP(3,5)*DP(4,5))
0535       DHK1=0.5*((DP(4,5)+DHKC)/DHKS-1.)
0536       DHK2=0.5*((DP(3,5)+DHKC)/DHKS-1.)
0537       IN1=N+NR+4*IS-3
0538       P(IN1,5)=SQRT(DP(3,5)+2.*DHKC+DP(4,5))
0539       DO 540 J=1,4
0540       P(IN1,J)=(1.+DHK1)*DP(1,J)-DHK2*DP(2,J)
0541   540 P(IN1+1,J)=(1.+DHK2)*DP(2,J)-DHK1*DP(1,J)
0542
0543 C...Begin initialization: sum up energy, set starting position.
0544       ISAV=I
0545   550 NTRY=NTRY+1
0546       IF(NTRY.GT.100.AND.NTRYR.LE.4) THEN
0547         PARU12=4.*PARU12
0548         PARU13=2.*PARU13
0549         GOTO 130
0550       ELSEIF(NTRY.GT.100) THEN
0551         CALL LUERRM(14,'(LUSTRF:) caught in infinite loop')
0552         IF(MSTU(21).GE.1) RETURN
0553       ENDIF
0554       I=ISAV
0555       DO 560 J=1,4
0556       P(N+NRS,J)=0.
0557       DO 560 IS=1,NR
0558   560 P(N+NRS,J)=P(N+NRS,J)+P(N+IS,J)
0559       DO 570 JT=1,2
0560       IRANK(JT)=0
0561       IF(MJU(JT).NE.0) IRANK(JT)=NJS(JT)
0562       IF(NS.GT.NR) IRANK(JT)=1
0563       IE(JT)=K(N+1+(JT/2)*(NP-1),3)
0564       IN(3*JT+1)=N+NR+1+4*(JT/2)*(NS-1)
0565       IN(3*JT+2)=IN(3*JT+1)+1
0566       IN(3*JT+3)=N+NR+4*NS+2*JT-1
0567       DO 570 IN1=N+NR+2+JT,N+NR+4*NS-2+JT,4
0568       P(IN1,1)=2-JT
0569       P(IN1,2)=JT-1
0570   570 P(IN1,3)=1.
0571
0572 C...Initialize flavour and pT variables for open string.
0573       IF(NS.LT.NR) THEN
0574         PX(1)=0.
0575         PY(1)=0.
0576         IF(NS.EQ.1.AND.MJU(1)+MJU(2).EQ.0) CALL LUPTDI(0,PX(1),PY(1))
0577         PX(2)=-PX(1)
0578         PY(2)=-PY(1)
0579         DO 580 JT=1,2
0580         KFL(JT)=K(IE(JT),2)
0581         IF(MJU(JT).NE.0) KFL(JT)=KFJS(JT)
0582         MSTJ(93)=1
0583         PMQ(JT)=ULMASS(KFL(JT))
0584   580   GAM(JT)=0.
0585
0586 C...Closed string: random initial breakup flavour, pT and vertex.
0587       ELSE
0588         KFL(3)=INT(1.+(2.+PARJ(2))*RLU(0))*(-1)**INT(RLU(0)+0.5)
0589         CALL LUKFDI(KFL(3),0,KFL(1),KDUMP)
0590         KFL(2)=-KFL(1)
0591         IF(IABS(KFL(1)).GT.10.AND.RLU(0).GT.0.5) THEN
0592           KFL(2)=-(KFL(1)+ISIGN(10000,KFL(1)))
0593         ELSEIF(IABS(KFL(1)).GT.10) THEN
0594           KFL(1)=-(KFL(2)+ISIGN(10000,KFL(2)))
0595         ENDIF
0596         CALL LUPTDI(KFL(1),PX(1),PY(1))
0597         PX(2)=-PX(1)
0598         PY(2)=-PY(1)
0599         PR3=MIN(25.,0.1*P(N+NR+1,5)**2)
0600   590   CALL LUZDIS(KFL(1),KFL(2),PR3,Z)
0601         ZR=PR3/(Z*P(N+NR+1,5)**2)
0602         IF(ZR.GE.1.) GOTO 590
0603         DO 600 JT=1,2
0604         MSTJ(93)=1
0605         PMQ(JT)=ULMASS(KFL(JT))
0606         GAM(JT)=PR3*(1.-Z)/Z
0607         IN1=N+NR+3+4*(JT/2)*(NS-1)
0608         P(IN1,JT)=1.-Z
0609         P(IN1,3-JT)=JT-1
0610         P(IN1,3)=(2-JT)*(1.-Z)+(JT-1)*Z
0611         P(IN1+1,JT)=ZR
0612         P(IN1+1,3-JT)=2-JT
0613   600   P(IN1+1,3)=(2-JT)*(1.-ZR)+(JT-1)*ZR
0614       ENDIF
0615
0616 C...Find initial transverse directions (i.e. spacelike four-vectors).
0617       DO 640 JT=1,2
0618       IF(JT.EQ.1.OR.NS.EQ.NR-1) THEN
0619         IN1=IN(3*JT+1)
0620         IN3=IN(3*JT+3)
0621         DO 610 J=1,4
0622         DP(1,J)=P(IN1,J)
0623         DP(2,J)=P(IN1+1,J)
0624         DP(3,J)=0.
0625   610   DP(4,J)=0.
0626         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)
0627         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)
0628         DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4)
0629         DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4)
0630         DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4)
0631         IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1.
0632         IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1.
0633         IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1.
0634         IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1.
0635         DHC12=DFOUR(1,2)
0636         DHCX1=DFOUR(3,1)/DHC12
0637         DHCX2=DFOUR(3,2)/DHC12
0638         DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12)
0639         DHCY1=DFOUR(4,1)/DHC12
0640         DHCY2=DFOUR(4,2)/DHC12
0641         DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12
0642         DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2)
0643         DO 620 J=1,4
0644         DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J))
0645         P(IN3,J)=DP(3,J)
0646   620   P(IN3+1,J)=DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)-
0647      &  DHCYX*DP(3,J))
0648       ELSE
0649         DO 630 J=1,4
0650         P(IN3+2,J)=P(IN3,J)
0651   630   P(IN3+3,J)=P(IN3+1,J)
0652       ENDIF
0653   640 CONTINUE
0654
0655 C...Remove energy used up in junction string fragmentation.
0656       IF(MJU(1)+MJU(2).GT.0) THEN
0657         DO 660 JT=1,2
0658         IF(NJS(JT).EQ.0) GOTO 660
0659         DO 650 J=1,4
0660   650   P(N+NRS,J)=P(N+NRS,J)-PJS(JT+2,J)
0661   660   CONTINUE
0662       ENDIF
0663
0664 C...Produce new particle: side, origin.
0665   670 I=I+1
0666       IF(2*I-NSAV.GE.MSTU(4)-MSTU(32)-5) THEN
0667         CALL LUERRM(11,'(LUSTRF:) no more memory left in LUJETS')
0668         IF(MSTU(21).GE.1) RETURN
0669       ENDIF
0670       JT=1.5+RLU(0)
0671       IF(IABS(KFL(3-JT)).GT.10) JT=3-JT
0672       JR=3-JT
0673       JS=3-2*JT
0674       IRANK(JT)=IRANK(JT)+1
0675       K(I,1)=1
0676       K(I,3)=IE(JT)
0677       K(I,4)=0
0678       K(I,5)=0
0679
0680 C...Generate flavour, hadron and pT.
0681   680 CALL LUKFDI(KFL(JT),0,KFL(3),K(I,2))
0682       IF(K(I,2).EQ.0) GOTO 550
0683       IF(MSTJ(12).GE.3.AND.IRANK(JT).EQ.1.AND.IABS(KFL(JT)).LE.10.AND.
0684      &IABS(KFL(3)).GT.10) THEN
0685         IF(RLU(0).GT.PARJ(19)) GOTO 680
0686       ENDIF
0687       P(I,5)=ULMASS(K(I,2))
0688       CALL LUPTDI(KFL(JT),PX(3),PY(3))
0689       PR(JT)=P(I,5)**2+(PX(JT)+PX(3))**2+(PY(JT)+PY(3))**2
0690
0691 C...Final hadrons for small invariant mass.
0692       MSTJ(93)=1
0693       PMQ(3)=ULMASS(KFL(3))
0694       WMIN=PARJ(32+MSTJ(11))+PMQ(1)+PMQ(2)+PARJ(36)*PMQ(3)
0695       IF(IABS(KFL(JT)).GT.10.AND.IABS(KFL(3)).GT.10) WMIN=
0696      &WMIN-0.5*PARJ(36)*PMQ(3)
0697       WREM2=FOUR(N+NRS,N+NRS)
0698       IF(WREM2.LT.0.10) GOTO 550
0699       IF(WREM2.LT.MAX(WMIN*(1.+(2.*RLU(0)-1.)*PARJ(37)),
0700      &PARJ(32)+PMQ(1)+PMQ(2))**2) GOTO 810
0701
0702 C...Choose z, which gives Gamma. Shift z for heavy flavours.
0703       CALL LUZDIS(KFL(JT),KFL(3),PR(JT),Z)
0704       KFL1A=IABS(KFL(1))
0705       KFL2A=IABS(KFL(2))
0706       IF(MAX(MOD(KFL1A,10),MOD(KFL1A/1000,10),MOD(KFL2A,10),
0707      &MOD(KFL2A/1000,10)).GE.4) THEN
0708         PR(JR)=(PMQ(JR)+PMQ(3))**2+(PX(JR)-PX(3))**2+(PY(JR)-PY(3))**2
0709         PW12=SQRT(MAX(0.,(WREM2-PR(1)-PR(2))**2-4.*PR(1)*PR(2)))
0710         Z=(WREM2+PR(JT)-PR(JR)+PW12*(2.*Z-1.))/(2.*WREM2)
0711         PR(JR)=(PMQ(JR)+PARJ(32+MSTJ(11)))**2+(PX(JR)-PX(3))**2+
0712      &  (PY(JR)-PY(3))**2
0713         IF((1.-Z)*(WREM2-PR(JT)/Z).LT.PR(JR)) GOTO 810
0714       ENDIF
0715       GAM(3)=(1.-Z)*(GAM(JT)+PR(JT)/Z)
0716       DO 690 J=1,3
0717   690 IN(J)=IN(3*JT+J)
0718
0719 C...Stepping within or from 'low' string region easy.
0720       IF(IN(1)+1.EQ.IN(2).AND.Z*P(IN(1)+2,3)*P(IN(2)+2,3)*
0721      &P(IN(1),5)**2.GE.PR(JT)) THEN
0722         P(IN(JT)+2,4)=Z*P(IN(JT)+2,3)
0723         P(IN(JR)+2,4)=PR(JT)/(P(IN(JT)+2,4)*P(IN(1),5)**2)
0724         DO 700 J=1,4
0725   700   P(I,J)=(PX(JT)+PX(3))*P(IN(3),J)+(PY(JT)+PY(3))*P(IN(3)+1,J)
0726         GOTO 770
0727       ELSEIF(IN(1)+1.EQ.IN(2)) THEN
0728         P(IN(JR)+2,4)=P(IN(JR)+2,3)
0729         P(IN(JR)+2,JT)=1.
0730         IN(JR)=IN(JR)+4*JS
0731         IF(JS*IN(JR).GT.JS*IN(4*JR)) GOTO 550
0732         IF(FOUR(IN(1),IN(2)).LE.1E-2) THEN
0733           P(IN(JT)+2,4)=P(IN(JT)+2,3)
0734           P(IN(JT)+2,JT)=0.
0735           IN(JT)=IN(JT)+4*JS
0736         ENDIF
0737       ENDIF
0738
0739 C...Find new transverse directions (i.e. spacelike string vectors).
0740   710 IF(JS*IN(1).GT.JS*IN(3*JR+1).OR.JS*IN(2).GT.JS*IN(3*JR+2).OR.
0741      &IN(1).GT.IN(2)) GOTO 550
0742       IF(IN(1).NE.IN(3*JT+1).OR.IN(2).NE.IN(3*JT+2)) THEN
0743         DO 720 J=1,4
0744         DP(1,J)=P(IN(1),J)
0745         DP(2,J)=P(IN(2),J)
0746         DP(3,J)=0.
0747   720   DP(4,J)=0.
0748         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)
0749         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)
0750         DHC12=DFOUR(1,2)
0751         IF(DHC12.LE.1E-2) THEN
0752           P(IN(JT)+2,4)=P(IN(JT)+2,3)
0753           P(IN(JT)+2,JT)=0.
0754           IN(JT)=IN(JT)+4*JS
0755           GOTO 710
0756         ENDIF
0757         IN(3)=N+NR+4*NS+5
0758         DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4)
0759         DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4)
0760         DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4)
0761         IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1.
0762         IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1.
0763         IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1.
0764         IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1.
0765         DHCX1=DFOUR(3,1)/DHC12
0766         DHCX2=DFOUR(3,2)/DHC12
0767         DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12)
0768         DHCY1=DFOUR(4,1)/DHC12
0769         DHCY2=DFOUR(4,2)/DHC12
0770         DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12
0771         DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2)
0772         DO 730 J=1,4
0773         DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J))
0774         P(IN(3),J)=DP(3,J)
0775   730   P(IN(3)+1,J)=DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)-
0776      &  DHCYX*DP(3,J))
0777 C...Express pT with respect to new axes, if sensible.
0778         PXP=-(PX(3)*FOUR(IN(3*JT+3),IN(3))+PY(3)*
0779      &  FOUR(IN(3*JT+3)+1,IN(3)))
0780         PYP=-(PX(3)*FOUR(IN(3*JT+3),IN(3)+1)+PY(3)*
0781      &  FOUR(IN(3*JT+3)+1,IN(3)+1))
0782         IF(ABS(PXP**2+PYP**2-PX(3)**2-PY(3)**2).LT.0.01) THEN
0783           PX(3)=PXP
0784           PY(3)=PYP
0785         ENDIF
0786       ENDIF
0787
0788 C...Sum up known four-momentum. Gives coefficients for m2 expression.
0789       DO 750 J=1,4
0790       DHG(J)=0.
0791       P(I,J)=PX(JT)*P(IN(3*JT+3),J)+PY(JT)*P(IN(3*JT+3)+1,J)+
0792      &PX(3)*P(IN(3),J)+PY(3)*P(IN(3)+1,J)
0793       DO 740 IN1=IN(3*JT+1),IN(1)-4*JS,4*JS
0794   740 P(I,J)=P(I,J)+P(IN1+2,3)*P(IN1,J)
0795       DO 750 IN2=IN(3*JT+2),IN(2)-4*JS,4*JS
0796   750 P(I,J)=P(I,J)+P(IN2+2,3)*P(IN2,J)
0797       DHM(1)=FOUR(I,I)
0798       DHM(2)=2.*FOUR(I,IN(1))
0799       DHM(3)=2.*FOUR(I,IN(2))
0800       DHM(4)=2.*FOUR(IN(1),IN(2))
0801
0802 C...Find coefficients for Gamma expression.
0803       DO 760 IN2=IN(1)+1,IN(2),4
0804       DO 760 IN1=IN(1),IN2-1,4
0805       DHC=2.*FOUR(IN1,IN2)
0806       DHG(1)=DHG(1)+P(IN1+2,JT)*P(IN2+2,JT)*DHC
0807       IF(IN1.EQ.IN(1)) DHG(2)=DHG(2)-JS*P(IN2+2,JT)*DHC
0808       IF(IN2.EQ.IN(2)) DHG(3)=DHG(3)+JS*P(IN1+2,JT)*DHC
0809   760 IF(IN1.EQ.IN(1).AND.IN2.EQ.IN(2)) DHG(4)=DHG(4)-DHC
0810
0811 C...Solve (m2, Gamma) equation system for energies taken.
0812       DHS1=DHM(JR+1)*DHG(4)-DHM(4)*DHG(JR+1)
0813       IF(ABS(DHS1).LT.1E-4) GOTO 550
0814       DHS2=DHM(4)*(GAM(3)-DHG(1))-DHM(JT+1)*DHG(JR+1)-DHG(4)*
0815      &(P(I,5)**2-DHM(1))+DHG(JT+1)*DHM(JR+1)
0816       DHS3=DHM(JT+1)*(GAM(3)-DHG(1))-DHG(JT+1)*(P(I,5)**2-DHM(1))
0817       P(IN(JR)+2,4)=0.5*(SQRT(MAX(0D0,DHS2**2-4.*DHS1*DHS3))/ABS(DHS1)-
0818      &DHS2/DHS1)
0819       IF(DHM(JT+1)+DHM(4)*P(IN(JR)+2,4).LE.0.) GOTO 550
0820       P(IN(JT)+2,4)=(P(I,5)**2-DHM(1)-DHM(JR+1)*P(IN(JR)+2,4))/
0821      &(DHM(JT+1)+DHM(4)*P(IN(JR)+2,4))
0822
0823 C...Step to new region if necessary.
0824       IF(P(IN(JR)+2,4).GT.P(IN(JR)+2,3)) THEN
0825         P(IN(JR)+2,4)=P(IN(JR)+2,3)
0826         P(IN(JR)+2,JT)=1.
0827         IN(JR)=IN(JR)+4*JS
0828         IF(JS*IN(JR).GT.JS*IN(4*JR)) GOTO 550
0829         IF(FOUR(IN(1),IN(2)).LE.1E-2) THEN
0830           P(IN(JT)+2,4)=P(IN(JT)+2,3)
0831           P(IN(JT)+2,JT)=0.
0832           IN(JT)=IN(JT)+4*JS
0833         ENDIF
0834         GOTO 710
0835       ELSEIF(P(IN(JT)+2,4).GT.P(IN(JT)+2,3)) THEN
0836         P(IN(JT)+2,4)=P(IN(JT)+2,3)
0837         P(IN(JT)+2,JT)=0.
0838         IN(JT)=IN(JT)+4*JS
0839         GOTO 710
0840       ENDIF
0841
0842 C...Four-momentum of particle. Remaining quantities. Loop back.
0843   770 DO 780 J=1,4
0844       P(I,J)=P(I,J)+P(IN(1)+2,4)*P(IN(1),J)+P(IN(2)+2,4)*P(IN(2),J)
0845   780 P(N+NRS,J)=P(N+NRS,J)-P(I,J)
0846       IF(P(I,4).LE.0.) GOTO 550
0847       KFL(JT)=-KFL(3)
0848       PMQ(JT)=PMQ(3)
0849       PX(JT)=-PX(3)
0850       PY(JT)=-PY(3)
0851       GAM(JT)=GAM(3)
0852       IF(IN(3).NE.IN(3*JT+3)) THEN
0853         DO 790 J=1,4
0854         P(IN(3*JT+3),J)=P(IN(3),J)
0855   790   P(IN(3*JT+3)+1,J)=P(IN(3)+1,J)
0856       ENDIF
0857       DO 800 JQ=1,2
0858       IN(3*JT+JQ)=IN(JQ)
0859       P(IN(JQ)+2,3)=P(IN(JQ)+2,3)-P(IN(JQ)+2,4)
0860   800 P(IN(JQ)+2,JT)=P(IN(JQ)+2,JT)-JS*(3-2*JQ)*P(IN(JQ)+2,4)
0861       GOTO 670
0862
0863 C...Final hadron: side, flavour, hadron, mass.
0864   810 I=I+1
0865       K(I,1)=1
0866       K(I,3)=IE(JR)
0867       K(I,4)=0
0868       K(I,5)=0
0869       CALL LUKFDI(KFL(JR),-KFL(3),KFLDMP,K(I,2))
0870       IF(K(I,2).EQ.0) GOTO 550
0871       P(I,5)=ULMASS(K(I,2))
0872       PR(JR)=P(I,5)**2+(PX(JR)-PX(3))**2+(PY(JR)-PY(3))**2
0873
0874 C...Final two hadrons: find common setup of four-vectors.
0875       JQ=1
0876       IF(P(IN(4)+2,3)*P(IN(5)+2,3)*FOUR(IN(4),IN(5)).LT.P(IN(7),3)*
0877      &P(IN(8),3)*FOUR(IN(7),IN(8))) JQ=2
0878       DHC12=FOUR(IN(3*JQ+1),IN(3*JQ+2))
0879       DHR1=FOUR(N+NRS,IN(3*JQ+2))/DHC12
0880       DHR2=FOUR(N+NRS,IN(3*JQ+1))/DHC12
0881       IF(IN(4).NE.IN(7).OR.IN(5).NE.IN(8)) THEN
0882         PX(3-JQ)=-FOUR(N+NRS,IN(3*JQ+3))-PX(JQ)
0883         PY(3-JQ)=-FOUR(N+NRS,IN(3*JQ+3)+1)-PY(JQ)
0884         PR(3-JQ)=P(I+(JT+JQ-3)**2-1,5)**2+(PX(3-JQ)+(2*JQ-3)*JS*
0885      &  PX(3))**2+(PY(3-JQ)+(2*JQ-3)*JS*PY(3))**2
0886       ENDIF
0887
0888 C...Solve kinematics for final two hadrons, if possible.
0889       WREM2=WREM2+(PX(1)+PX(2))**2+(PY(1)+PY(2))**2
0890       FD=(SQRT(PR(1))+SQRT(PR(2)))/SQRT(WREM2)
0891       IF(MJU(1)+MJU(2).NE.0.AND.I.EQ.ISAV+2.AND.FD.GE.1.) GOTO 180
0892       IF(FD.GE.1.) GOTO 550
0893       FA=WREM2+PR(JT)-PR(JR)
0894       IF(MSTJ(11).EQ.2) PREV=0.5*FD**PARJ(37+MSTJ(11))
0895       IF(MSTJ(11).NE.2) PREV=0.5*EXP(MAX(-100.,LOG(FD)*
0896      &PARJ(37+MSTJ(11))*(PR(1)+PR(2))**2))
0897       FB=SIGN(SQRT(MAX(0.,FA**2-4.*WREM2*PR(JT))),JS*(RLU(0)-PREV))
0898       KFL1A=IABS(KFL(1))
0899       KFL2A=IABS(KFL(2))
0900       IF(MAX(MOD(KFL1A,10),MOD(KFL1A/1000,10),MOD(KFL2A,10),
0901      &MOD(KFL2A/1000,10)).GE.6) FB=SIGN(SQRT(MAX(0.,FA**2-
0902      &4.*WREM2*PR(JT))),FLOAT(JS))
0903       DO 820 J=1,4
0904       P(I-1,J)=(PX(JT)+PX(3))*P(IN(3*JQ+3),J)+(PY(JT)+PY(3))*
0905      &P(IN(3*JQ+3)+1,J)+0.5*(DHR1*(FA+FB)*P(IN(3*JQ+1),J)+
0906      &DHR2*(FA-FB)*P(IN(3*JQ+2),J))/WREM2
0907   820 P(I,J)=P(N+NRS,J)-P(I-1,J)
0908
0909 C...Mark jets as fragmented and give daughter pointers.
0910       N=I-NRS+1
0911       DO 830 I=NSAV+1,NSAV+NP
0912       IM=K(I,3)
0913       K(IM,1)=K(IM,1)+10
0914       IF(MSTU(16).NE.2) THEN
0915         K(IM,4)=NSAV+1
0916         K(IM,5)=NSAV+1
0917       ELSE
0918         K(IM,4)=NSAV+2
0919         K(IM,5)=N
0920       ENDIF
0921   830 CONTINUE
0922
0923 C...Document string system. Move up particles.
0924       NSAV=NSAV+1
0925       K(NSAV,1)=11
0926       K(NSAV,2)=92
0927       K(NSAV,3)=IP
0928       K(NSAV,4)=NSAV+1
0929       K(NSAV,5)=N
0930       DO 840 J=1,4
0931       P(NSAV,J)=DPS(J)
0932   840 V(NSAV,J)=V(IP,J)
0933       P(NSAV,5)=SQRT(MAX(0D0,DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2))
0934       V(NSAV,5)=0.
0935       DO 850 I=NSAV+1,N
0936       DO 850 J=1,5
0937       K(I,J)=K(I+NRS-1,J)
0938       P(I,J)=P(I+NRS-1,J)
0939   850 V(I,J)=0.
0940
0941 C...Order particles in rank along the chain. Update mother pointer.
0942       DO 860 I=NSAV+1,N
0943       DO 860 J=1,5
0944       K(I-NSAV+N,J)=K(I,J)
0945   860 P(I-NSAV+N,J)=P(I,J)
0946       I1=NSAV
0947       DO 880 I=N+1,2*N-NSAV
0948       IF(K(I,3).NE.IE(1)) GOTO 880
0949       I1=I1+1
0950       DO 870 J=1,5
0951       K(I1,J)=K(I,J)
0952   870 P(I1,J)=P(I,J)
0953       IF(MSTU(16).NE.2) K(I1,3)=NSAV
0954   880 CONTINUE
0955       DO 900 I=2*N-NSAV,N+1,-1
0956       IF(K(I,3).EQ.IE(1)) GOTO 900
0957       I1=I1+1
0958       DO 890 J=1,5
0959       K(I1,J)=K(I,J)
0960   890 P(I1,J)=P(I,J)
0961       IF(MSTU(16).NE.2) K(I1,3)=NSAV
0962   900 CONTINUE
0963
0964 C...Boost back particle system. Set production vertices.
0965       CALL LUDBRB(NSAV+1,N,0.,0.,DPS(1)/DPS(4),DPS(2)/DPS(4),
0966      &DPS(3)/DPS(4))
0967       DO 910 I=NSAV+1,N
0968       DO 910 J=1,4
0969   910 V(I,J)=V(IP,J)
0970
0971       RETURN
0972       END