pythiaeRHIC/pythia/pystrf.f

0001
0002 C*********************************************************************
0003
0004 C...PYSTRF
0005 C...Handles the fragmentation of an arbitrary colour singlet
0006 C...jet system according to the Lund string fragmentation model.
0007
0008       SUBROUTINE PYSTRF(IP)
0009
0010 C...Double precision and integer declarations.
0011       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0012       IMPLICIT INTEGER(I-N)
0013       INTEGER PYK,PYCHGE,PYCOMP
0014 C...Commonblocks.
0015       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0016       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0017       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0018       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/
0019 C...Local arrays. All MOPS variables ends with MO
0020       DIMENSION DPS(5),KFL(3),PMQ(3),PX(3),PY(3),GAM(3),IE(2),PR(2),
0021      &IN(9),DHM(4),DHG(4),DP(5,5),IRANK(2),MJU(4),IJU(6),PJU(5,5),
0022      &TJU(5),KFJH(2),NJS(2),KFJS(2),PJS(4,5),MSTU9T(8),PARU9T(8),
0023      &INMO(9),PM2QMO(2),XTMO(2),EJSTR(2),IJUORI(2),IBARRK(2),
0024      &PBST(3,5),TJUOLD(5)
0025
0026 C...Function: four-product of two vectors.
0027       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)
0028       DFOUR(I,J)=DP(I,4)*DP(J,4)-DP(I,1)*DP(J,1)-DP(I,2)*DP(J,2)-
0029      &DP(I,3)*DP(J,3)
0030
0031 C...Reset counters.
0032       MSTJ(91)=0
0033       NSAV=N
0034       MSTU90=MSTU(90)
0035       NP=0
0036       KQSUM=0
0037       DO 100 J=1,5
0038         DPS(J)=0D0
0039   100 CONTINUE
0040       MJU(1)=0
0041       MJU(2)=0
0042       NTRYFN=0
0043       IJUORI(1)=0
0044       IJUORI(2)=0
0045
0046 C...Identify parton system.
0047       I=IP-1
0048   110 I=I+1
0049       IF(I.GT.MIN(N,MSTU(4)-MSTU(32))) THEN
0050         CALL PYERRM(12,'(PYSTRF:) failed to reconstruct jet system')
0051         IF(MSTU(21).GE.1) RETURN
0052       ENDIF
0053       IF(K(I,1).NE.1.AND.K(I,1).NE.2.AND.K(I,1).NE.41) GOTO 110
0054       KC=PYCOMP(K(I,2))
0055       IF(KC.EQ.0) GOTO 110
0056       KQ=KCHG(KC,2)*ISIGN(1,K(I,2))
0057       IF(KQ.EQ.0.AND.K(I,1).NE.41) GOTO 110
0058       IF(N+5*NP+11.GT.MSTU(4)-MSTU(32)-5) THEN
0059         CALL PYERRM(11,'(PYSTRF:) no more memory left in PYJETS')
0060         IF(MSTU(21).GE.1) RETURN
0061       ENDIF
0062
0063 C...Take copy of partons to be considered. Check flavour sum.
0064       NP=NP+1
0065       DO 120 J=1,5
0066         K(N+NP,J)=K(I,J)
0067         P(N+NP,J)=P(I,J)
0068         IF(J.NE.4) DPS(J)=DPS(J)+P(I,J)
0069   120 CONTINUE
0070       DPS(4)=DPS(4)+SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
0071       K(N+NP,3)=I
0072       IF(KQ.NE.2) KQSUM=KQSUM+KQ
0073       IF(K(I,1).EQ.41) THEN
0074         IF(MOD(KQSUM,2).EQ.0.AND.MJU(1).EQ.0) THEN
0075           MJU(1)=N+NP
0076           IJUORI(1)=I
0077         ELSE
0078           MJU(2)=N+NP
0079           IJUORI(2)=I
0080         ENDIF
0081       ENDIF
0082       IF(K(I,1).EQ.2.OR.K(I,1).EQ.41) GOTO 110
0083       IF(MOD(KQSUM,3).NE.0) THEN
0084         CALL PYERRM(12,'(PYSTRF:) unphysical flavour combination')
0085         IF(MSTU(21).GE.1) RETURN
0086       ENDIF
0087       IF(MJU(1).GT.0.OR.MJU(2).GT.0) MSTU(29)=1
0088
0089 C...Boost copied system to CM frame (for better numerical precision).
0090       IF(ABS(DPS(3)).LT.0.99D0*DPS(4)) THEN
0091         MBST=0
0092         MSTU(33)=1
0093         CALL PYROBO(N+1,N+NP,0D0,0D0,-DPS(1)/DPS(4),-DPS(2)/DPS(4),
0094      &  -DPS(3)/DPS(4))
0095       ELSE
0096         MBST=1
0097         HHBZ=SQRT(MAX(1D-6,DPS(4)+DPS(3))/MAX(1D-6,DPS(4)-DPS(3)))
0098         DO 130 I=N+1,N+NP
0099           HHPMT=P(I,1)**2+P(I,2)**2+P(I,5)**2
0100           IF(P(I,3).GT.0D0) THEN
0101             HHPEZ=MAX(1D-10,(P(I,4)+P(I,3))/HHBZ)
0102             P(I,3)=0.5D0*(HHPEZ-HHPMT/HHPEZ)
0103             P(I,4)=0.5D0*(HHPEZ+HHPMT/HHPEZ)
0104           ELSE
0105             HHPEZ=MAX(1D-10,(P(I,4)-P(I,3))*HHBZ)
0106             P(I,3)=-0.5D0*(HHPEZ-HHPMT/HHPEZ)
0107             P(I,4)=0.5D0*(HHPEZ+HHPMT/HHPEZ)
0108           ENDIF
0109   130   CONTINUE
0110       ENDIF
0111
0112 C...Search for very nearby partons that may be recombined.
0113       NTRYR=0
0114       NTRYWR=0
0115       PARU12=PARU(12)
0116       PARU13=PARU(13)
0117       MJU(3)=MJU(1)
0118       MJU(4)=MJU(2)
0119       NR=NP
0120       NRMIN=2
0121       IF(MJU(1).GT.0) NRMIN=NRMIN+2
0122       IF(MJU(2).GT.0) NRMIN=NRMIN+2
0123   140 IF(NR.GT.NRMIN) THEN
0124         PDRMIN=2D0*PARU12
0125         DO 150 I=N+1,N+NR
0126           IF(I.EQ.N+NR.AND.IABS(K(N+1,2)).NE.21) GOTO 150
0127           I1=I+1
0128           IF(I.EQ.N+NR) I1=N+1
0129           IF(K(I,1).EQ.41.OR.K(I1,1).EQ.41) GOTO 150
0130           IF(MJU(1).NE.0.AND.I1.LT.MJU(1).AND.IABS(K(I1,2)).NE.21)
0131      &    GOTO 150
0132           IF(MJU(2).NE.0.AND.I.GT.MJU(2).AND.IABS(K(I,2)).NE.21)
0133      &    GOTO 150
0134           PAP=SQRT((P(I,1)**2+P(I,2)**2+P(I,3)**2)*(P(I1,1)**2+
0135      &    P(I1,2)**2+P(I1,3)**2))
0136           PVP=P(I,1)*P(I1,1)+P(I,2)*P(I1,2)+P(I,3)*P(I1,3)
0137           PDR=4D0*(PAP-PVP)**2/MAX(1D-6,PARU13**2*PAP+2D0*(PAP-PVP))
0138           IF(PDR.LT.PDRMIN) THEN
0139             IR=I
0140             PDRMIN=PDR
0141           ENDIF
0142   150   CONTINUE
0143
0144 C...Recombine very nearby partons to avoid machine precision problems.
0145         IF(PDRMIN.LT.PARU12.AND.IR.EQ.N+NR) THEN
0146           DO 160 J=1,4
0147             P(N+1,J)=P(N+1,J)+P(N+NR,J)
0148   160     CONTINUE
0149           P(N+1,5)=SQRT(MAX(0D0,P(N+1,4)**2-P(N+1,1)**2-P(N+1,2)**2-
0150      &    P(N+1,3)**2))
0151           NR=NR-1
0152           GOTO 140
0153         ELSEIF(PDRMIN.LT.PARU12) THEN
0154           DO 170 J=1,4
0155             P(IR,J)=P(IR,J)+P(IR+1,J)
0156   170     CONTINUE
0157           P(IR,5)=SQRT(MAX(0D0,P(IR,4)**2-P(IR,1)**2-P(IR,2)**2-
0158      &    P(IR,3)**2))
0159           IF(MJU(2).NE.0.AND.IR.GT.MJU(2)) K(IR,2)=K(IR+1,2)
0160           DO 190 I=IR+1,N+NR-1
0161             K(I,1)=K(I+1,1)
0162             K(I,2)=K(I+1,2)
0163             DO 180 J=1,5
0164               P(I,J)=P(I+1,J)
0165   180       CONTINUE
0166   190     CONTINUE
0167           IF(IR.EQ.N+NR-1) K(IR,2)=K(N+NR,2)
0168           NR=NR-1
0169           IF(MJU(1).GT.IR) MJU(1)=MJU(1)-1
0170           IF(MJU(2).GT.IR) MJU(2)=MJU(2)-1
0171           GOTO 140
0172         ENDIF
0173       ENDIF
0174       NTRYR=NTRYR+1
0175
0176 C...Reset particle counter. Skip ahead if no junctions are present;
0177 C...this is usually the case!
0178       NRS=MAX(5*NR+11,NP)
0179       NTRY=0
0180   200 NTRY=NTRY+1
0181       IF(NTRY.GT.100.AND.NTRYR.LE.8.AND.NR.GT.NRMIN) THEN
0182         PARU12=4D0*PARU12
0183         PARU13=2D0*PARU13
0184         GOTO 140
0185       ELSEIF(NTRY.GT.100.OR.NTRYR.GT.100) THEN
0186         CALL PYERRM(14,'(PYSTRF:) caught in infinite loop')
0187         IF(MSTU(21).GE.1) RETURN
0188       ENDIF
0189       I=N+NRS
0190       MSTU(90)=MSTU90
0191       IF(MJU(1).EQ.0.AND.MJU(2).EQ.0) GOTO 650
0192       IF(MSTJ(12).GE.4) CALL PYERRM(29,'(PYSTRF:) sorry,'//
0193      &     ' junction strings not handled by MSTJ(12)>3 options')
0194       DO 640 JT=1,2
0195         NJS(JT)=0
0196         IF(MJU(JT).EQ.0) GOTO 640
0197         JS=3-2*JT
0198
0199 C++SKANDS
0200 C...Find and sum up momentum on three sides of junction.
0201 C...Begin with previous boost = zero.
0202         IJRFIT=0
0203         DO 210 IX=1,3
0204           TJUOLD(IX)=0D0
0205   210   CONTINUE
0206         TJUOLD(4)=1D0
0207   220   IU=0
0208 C...Beginning and end of string system in event record.
0209         I1BEG=N+1+(JT-1)*(NR-1)
0210         I1END=N+NR+(JT-1)*(1-NR)
0211 C...Look for junction string piece end points
0212         DO 230 I1=I1BEG,I1END,JS
0213           IF(K(I1,2).NE.21.AND.IU.LE.5.AND.IJRFIT.EQ.0) THEN
0214 C...Store junction string piece end points.
0215 C                 1-junction systems        2-junction systems
0216 C           IU :  1     2     3   4     1     2   3     4   5     6
0217 C       IJU(IU):  q-g-g-q-g-g-j-g-q     q-g-g-q-g-j-g-g-j-g-q-g-g-q
0218             IU=IU+1
0219             IJU(IU)=I1
0220           ENDIF
0221 C...Sum over momenta, from junction outwards.
0222   230   CONTINUE
0223         DO 280 IU=1,3
0224           PWT=0D0
0225 C...Initialize junction drag and string piece 4-vectors.
0226           DO 240 J=1,5
0227             PBST(IU,J)=0D0
0228             PJU(IU,J)=0D0
0229   240     CONTINUE
0230 C...First two branches. Inwards out means opposite direction to JS.
0231 C...(JS is 1 for JT=1, -1 for JT=2)
0232           IF (IU.LT.3) THEN
0233             I1A=IJU(IU+1)-JS
0234             I1B=IJU(IU)
0235             IDIR=-JS
0236 C...Last branch (gq or gjgqgq). Direction now reversed.
0237           ELSE
0238             I1A=IJU(IU)+JS
0239             I1B=I1END
0240             IDIR=JS
0241           ENDIF
0242           DO 270 I1=I1A,I1B,IDIR
0243 C...Sum up momentum directions with exponential suppression
0244 C...for use in finding junction rest frame below.
0245             IF (K(I1,2).EQ.88) THEN
0246 C...gjgqgq type system encountered. Use current PWT as start
0247 C...for both strings.
0248               PWTOLD=PWT
0249             ELSE
0250               IF (I1.EQ.IJU(5)+IDIR) PWT=PWTOLD
0251 C...Sum up string piece (boosted) 4-momenta.
0252               DO 250 J=1,4
0253                 PJU(IU,J)=PJU(IU,J)+P(I1,J)
0254   250         CONTINUE
0255 C...Compute "junction drag" vectors from (boosted) 4-momenta (initial
0256 C...boost is zero, see above). Skip parton if suppression factor large.
0257               IF (PWT.GT.10D0) GOTO 270
0258 C...Compute momentum in current frame:
0259               TDP=TJUOLD(1)*P(I1,1)+TJUOLD(2)*P(I1,2)+TJUOLD(3)*P(I1,3)
0260               BFC=TDP/(1D0+TJUOLD(4))+P(I1,4)
0261               DO 260 J=1,3
0262                 PTMP=P(I1,J)+TJUOLD(J)*BFC
0263                 PBST(IU,J)=PBST(IU,J)+PTMP*EXP(-PWT)
0264   260         CONTINUE
0265 C...Boosted energy
0266               PTMP=TJUOLD(4)*P(I1,4)+TDP
0267               PBST(IU,4)=PBST(IU,J)+PTMP*EXP(-PWT)
0268               PWT=PWT+PTMP/PARJ(48)
0269             ENDIF
0270   270     CONTINUE
0271 C...Put |p| rather than m in 5th slot.
0272           PBST(IU,5)=SQRT(PBST(IU,1)**2+PBST(IU,2)**2+PBST(IU,3)**2)
0273           PJU(IU,5)=SQRT(PJU(IU,1)**2+PJU(IU,2)**2+PJU(IU,3)**2)
0274   280   CONTINUE
0275
0276 C...Calculate boost from present frame to next JRF candidate.
0277         IJRFIT=IJRFIT+1
0278         CALL PYJURF(PBST,TJU)
0279
0280 C...After some iterations do not take full step in new direction.
0281         IF(IJRFIT.GT.5) THEN
0282           REDUCE=0.8D0**(IJRFIT-5)
0283           TJU(1)=REDUCE*TJU(1)
0284           TJU(2)=REDUCE*TJU(2)
0285           TJU(3)=REDUCE*TJU(3)
0286           TJU(4)=SQRT(1D0+TJU(1)**2+TJU(2)**2+TJU(3)**2)
0287         ENDIF
0288
0289 C...Combine new boost (TJU) with old boost (TJUOLD)
0290         TMP=TJU(1)*TJUOLD(1)+TJU(2)*TJUOLD(2)+TJU(3)*TJUOLD(3)
0291         DO 290 IX=1,3
0292           TJUOLD(IX)=TJU(IX)+TJUOLD(IX)*(TMP/(1D0+TJUOLD(4))+TJU(4))
0293   290   CONTINUE
0294         TJUOLD(4)=SQRT(1D0+TJUOLD(1)**2+TJUOLD(2)**2+TJUOLD(3)**2)
0295
0296 C...If last boost small, accept JRF, else iterate.
0297 C...Also prevent possibility of infinite loop.
0298         IF (ABS((TJU(4)-1D0)/TJUOLD(4)).GT.0.01D0.AND.
0299      &  IJRFIT.LT.MSTJ(18)) THEN
0300           GOTO 220
0301         ELSEIF (IJRFIT.GE.MSTJ(18)) THEN
0302           CALL PYERRM(1,'(PYSTRF:) failed to converge on JRF')
0303         ENDIF
0304
0305 C...Now store total boost in TJU and change perception.
0306 C...TJUOLD = boost vector from CM of string syst -> JRF. Henceforth,
0307 C...TJU = junction motion vector in string CM, so the sign changes.
0308         DO 300 J=1,3
0309           TJU(J)=-TJUOLD(J)
0310   300   CONTINUE
0311         TJU(4)=SQRT(1D0+TJU(1)**2+TJU(2)**2+TJU(3)**2)
0312
0313 C--SKANDS
0314
0315 C...Calculate string piece energies in junction rest frame.
0316         DO 310 IU=1,3
0317           PJU(IU,5)=TJU(4)*PJU(IU,4)-TJU(1)*PJU(IU,1)-TJU(2)*PJU(IU,2)-
0318      &    TJU(3)*PJU(IU,3)
0319           PBST(IU,5)=TJU(4)*PBST(IU,4)-TJU(1)*PBST(IU,1)-
0320      &    TJU(2)*PBST(IU,2)-TJU(3)*PBST(IU,3)
0321   310   CONTINUE
0322
0323 C...Start preparing for fragmentation of two strings from junction.
0324         ISTA=I
0325         NTRYER=0
0326   320   NTRYER=NTRYER+1
0327         I=ISTA
0328         DO 620 IU=1,2
0329           NS=IABS(IJU(IU+1)-IJU(IU))
0330
0331 C...Junction strings: find longitudinal string directions.
0332           DO 350 IS=1,NS
0333             IS1=IJU(IU)+JS*(IS-1)
0334             IS2=IJU(IU)+JS*IS
0335             DO 330 J=1,5
0336               DP(1,J)=0.5D0*P(IS1,J)
0337               IF(IS.EQ.1) DP(1,J)=P(IS1,J)
0338               DP(2,J)=0.5D0*P(IS2,J)
0339               IF(IS.EQ.NS) DP(2,J)=(-PBST(IU,J)+2D0*PBST(IU,5)*TJU(J))*
0340      &        (PJU(IU,5)/PBST(IU,5))
0341   330       CONTINUE
0342             IF(IS.EQ.NS) DP(2,5)=SQRT(MAX(0D0,PJU(IU,4)**2-
0343      &      PJU(IU,1)**2-PJU(IU,2)**2-PJU(IU,3)**2))
0344             DP(3,5)=DFOUR(1,1)
0345             DP(4,5)=DFOUR(2,2)
0346             DHKC=DFOUR(1,2)
0347             IF(DP(3,5)+2D0*DHKC+DP(4,5).LE.0D0) THEN
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               DP(3,5)=0D0
0351               DP(4,5)=0D0
0352               DHKC=DFOUR(1,2)
0353             ENDIF
0354             DHKS=SQRT(DHKC**2-DP(3,5)*DP(4,5))
0355             DHK1=0.5D0*((DP(4,5)+DHKC)/DHKS-1D0)
0356             DHK2=0.5D0*((DP(3,5)+DHKC)/DHKS-1D0)
0357             IN1=N+NR+4*IS-3
0358             P(IN1,5)=SQRT(DP(3,5)+2D0*DHKC+DP(4,5))
0359             DO 340 J=1,4
0360               P(IN1,J)=(1D0+DHK1)*DP(1,J)-DHK2*DP(2,J)
0361               P(IN1+1,J)=(1D0+DHK2)*DP(2,J)-DHK1*DP(1,J)
0362   340       CONTINUE
0363   350     CONTINUE
0364
0365 C...Junction strings: initialize flavour, momentum and starting pos.
0366           ISAV=I
0367           MSTU91=MSTU(90)
0368   360     NTRY=NTRY+1
0369           IF(NTRY.GT.100.AND.NTRYR.LE.8.AND.NR.GT.NRMIN) THEN
0370             PARU12=4D0*PARU12
0371             PARU13=2D0*PARU13
0372             GOTO 140
0373           ELSEIF(NTRY.GT.100) THEN
0374             CALL PYERRM(14,'(PYSTRF:) caught in infinite loop')
0375             IF(MSTU(21).GE.1) RETURN
0376           ENDIF
0377           I=ISAV
0378           MSTU(90)=MSTU91
0379           IRANKJ=0
0380           IE(1)=K(N+1+(JT/2)*(NP-1),3)
0381           IF (MOD(JT+IU,2).NE.0) THEN
0382             IE(1)=K(IJU(IU),3)
0383             IF (NP-NR.NE.0) THEN
0384 C...If gluons have disappeared. Original IJU must be used.
0385               IT=IP
0386               NE=1
0387   370         IT=IT+1
0388               IF (K(IT,2).NE.21) THEN
0389                 NE=NE+1
0390               ENDIF
0391               IF (NE.EQ.IU+4*(JT-1)) THEN
0392                 IE(1)=IT
0393               ELSEIF (IT.LE.IP+NP) THEN
0394                 GOTO 370
0395               ELSE
0396                 CALL PYERRM(14,'(PYSTRF:) '//
0397      &               'Original IJU could not be reconstructed!')
0398               ENDIF
0399             ENDIF
0400           ENDIF
0401           IN(4)=N+NR+1
0402           IN(5)=IN(4)+1
0403           IN(6)=N+NR+4*NS+1
0404           DO 390 JQ=1,2
0405             DO 380 IN1=N+NR+2+JQ,N+NR+4*NS-2+JQ,4
0406               P(IN1,1)=2-JQ
0407               P(IN1,2)=JQ-1
0408               P(IN1,3)=1D0
0409   380       CONTINUE
0410   390     CONTINUE
0411           KFL(1)=K(IJU(IU),2)
0412           PX(1)=0D0
0413           PY(1)=0D0
0414           GAM(1)=0D0
0415           DO 400 J=1,5
0416             PJU(IU+3,J)=0D0
0417   400     CONTINUE
0418
0419 C...Junction strings: find initial transverse directions.
0420           DO 410 J=1,4
0421             DP(1,J)=P(IN(4),J)
0422             DP(2,J)=P(IN(4)+1,J)
0423             DP(3,J)=0D0
0424             DP(4,J)=0D0
0425   410     CONTINUE
0426           DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)
0427           DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)
0428           DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4)
0429           DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4)
0430           DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4)
0431           IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1D0
0432           IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1D0
0433           IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1D0
0434           IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1D0
0435           DHC12=DFOUR(1,2)
0436           DHCX1=DFOUR(3,1)/DHC12
0437           DHCX2=DFOUR(3,2)/DHC12
0438           DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12)
0439           DHCY1=DFOUR(4,1)/DHC12
0440           DHCY2=DFOUR(4,2)/DHC12
0441           DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12
0442           DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2)
0443           DO 420 J=1,4
0444             DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J))
0445             P(IN(6),J)=DP(3,J)
0446             P(IN(6)+1,J)=DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)-
0447      &      DHCYX*DP(3,J))
0448   420     CONTINUE
0449
0450 C...Junction strings: produce new particle, origin.
0451   430     I=I+1
0452           IF(2*I-NSAV.GE.MSTU(4)-MSTU(32)-5) THEN
0453             CALL PYERRM(11,'(PYSTRF:) no more memory left in PYJETS')
0454             IF(MSTU(21).GE.1) RETURN
0455           ENDIF
0456           IRANKJ=IRANKJ+1
0457           K(I,1)=1
0458           K(I,3)=IE(1)
0459           K(I,4)=0
0460           K(I,5)=0
0461
0462 C...Junction strings: generate flavour, hadron, pT, z and Gamma.
0463   440     CALL PYKFDI(KFL(1),0,KFL(3),K(I,2))
0464           IF(K(I,2).EQ.0) GOTO 360
0465           IF(IRANKJ.EQ.1.AND.IABS(KFL(1)).LE.10.AND.
0466      &    IABS(KFL(3)).GT.10) THEN
0467             IF(PYR(0).GT.PARJ(19)) GOTO 440
0468           ENDIF
0469           P(I,5)=PYMASS(K(I,2))
0470           CALL PYPTDI(KFL(1),PX(3),PY(3))
0471           PR(1)=P(I,5)**2+(PX(1)+PX(3))**2+(PY(1)+PY(3))**2
0472           CALL PYZDIS(KFL(1),KFL(3),PR(1),Z)
0473           IF(IABS(KFL(1)).GE.4.AND.IABS(KFL(1)).LE.8.AND.
0474      &    MSTU(90).LT.8) THEN
0475             MSTU(90)=MSTU(90)+1
0476             MSTU(90+MSTU(90))=I
0477             PARU(90+MSTU(90))=Z
0478           ENDIF
0479           GAM(3)=(1D0-Z)*(GAM(1)+PR(1)/Z)
0480           DO 450 J=1,3
0481             IN(J)=IN(3+J)
0482   450     CONTINUE
0483
0484 C...Junction strings: stepping within 'low' string region.
0485           IF(IN(1)+1.EQ.IN(2).AND.Z*P(IN(1)+2,3)*P(IN(2)+2,3)*
0486      &    P(IN(1),5)**2.GE.PR(1)) THEN
0487             P(IN(1)+2,4)=Z*P(IN(1)+2,3)
0488             P(IN(2)+2,4)=PR(1)/(P(IN(1)+2,4)*P(IN(1),5)**2)
0489             DO 460 J=1,4
0490               P(I,J)=(PX(1)+PX(3))*P(IN(3),J)+(PY(1)+PY(3))*P(IN(3)+1,J)
0491   460       CONTINUE
0492             GOTO 560
0493 C...Has used up energy of junction string, i.e. no more hadrons in it.
0494           ELSEIF(IN(1)+1.EQ.IN(2).AND.IN(1).EQ.N+NR+4*NS-3) THEN
0495             DO 470 J=1,5
0496               P(I,J)=0D0
0497   470       CONTINUE
0498             GOTO 600
0499 C...Stepping from 'low' string region
0500           ELSEIF(IN(1)+1.EQ.IN(2)) THEN
0501             P(IN(2)+2,4)=P(IN(2)+2,3)
0502             P(IN(2)+2,1)=1D0
0503             IN(2)=IN(2)+4
0504             IF(IN(2).GT.N+NR+4*NS) GOTO 360
0505             IF(FOUR(IN(1),IN(2)).LE.1D-2) THEN
0506               P(IN(1)+2,4)=P(IN(1)+2,3)
0507               P(IN(1)+2,1)=0D0
0508               IN(1)=IN(1)+4
0509             ENDIF
0510           ENDIF
0511
0512 C...Junction strings: find new transverse directions.
0513   480     IF(IN(1).GT.N+NR+4*NS.OR.IN(2).GT.N+NR+4*NS.OR.
0514      &    IN(1).GT.IN(2)) GOTO 360
0515           IF(IN(1).NE.IN(4).OR.IN(2).NE.IN(5)) THEN
0516             DO 490 J=1,4
0517               DP(1,J)=P(IN(1),J)
0518               DP(2,J)=P(IN(2),J)
0519               DP(3,J)=0D0
0520               DP(4,J)=0D0
0521   490       CONTINUE
0522             DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)
0523             DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)
0524             DHC12=DFOUR(1,2)
0525             IF(DHC12.LE.1D-2) THEN
0526               P(IN(1)+2,4)=P(IN(1)+2,3)
0527               P(IN(1)+2,1)=0D0
0528               IN(1)=IN(1)+4
0529               GOTO 480
0530             ENDIF
0531             IN(3)=N+NR+4*NS+5
0532             DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4)
0533             DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4)
0534             DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4)
0535             IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1D0
0536             IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1D0
0537             IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1D0
0538             IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1D0
0539             DHCX1=DFOUR(3,1)/DHC12
0540             DHCX2=DFOUR(3,2)/DHC12
0541             DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12)
0542             DHCY1=DFOUR(4,1)/DHC12
0543             DHCY2=DFOUR(4,2)/DHC12
0544             DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12
0545             DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2)
0546             DO 500 J=1,4
0547               DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J))
0548               P(IN(3),J)=DP(3,J)
0549               P(IN(3)+1,J)=DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)-
0550      &        DHCYX*DP(3,J))
0551   500       CONTINUE
0552 C...Express pT with respect to new axes, if sensible.
0553             PXP=-(PX(3)*FOUR(IN(6),IN(3))+PY(3)*FOUR(IN(6)+1,IN(3)))
0554             PYP=-(PX(3)*FOUR(IN(6),IN(3)+1)+PY(3)*FOUR(IN(6)+1,IN(3)+1))
0555             IF(ABS(PXP**2+PYP**2-PX(3)**2-PY(3)**2).LT.0.01D0) THEN
0556               PX(3)=PXP
0557               PY(3)=PYP
0558             ENDIF
0559           ENDIF
0560
0561 C...Junction strings: sum up known four-momentum, coefficients for m2.
0562           DO 530 J=1,4
0563             DHG(J)=0D0
0564             P(I,J)=PX(1)*P(IN(6),J)+PY(1)*P(IN(6)+1,J)+PX(3)*P(IN(3),J)+
0565      &      PY(3)*P(IN(3)+1,J)
0566             DO 510 IN1=IN(4),IN(1)-4,4
0567               P(I,J)=P(I,J)+P(IN1+2,3)*P(IN1,J)
0568   510       CONTINUE
0569             DO 520 IN2=IN(5),IN(2)-4,4
0570               P(I,J)=P(I,J)+P(IN2+2,3)*P(IN2,J)
0571   520       CONTINUE
0572   530     CONTINUE
0573           DHM(1)=FOUR(I,I)
0574           DHM(2)=2D0*FOUR(I,IN(1))
0575           DHM(3)=2D0*FOUR(I,IN(2))
0576           DHM(4)=2D0*FOUR(IN(1),IN(2))
0577
0578 C...Junction strings: find coefficients for Gamma expression.
0579           DO 550 IN2=IN(1)+1,IN(2),4
0580             DO 540 IN1=IN(1),IN2-1,4
0581               DHC=2D0*FOUR(IN1,IN2)
0582               DHG(1)=DHG(1)+P(IN1+2,1)*P(IN2+2,1)*DHC
0583               IF(IN1.EQ.IN(1)) DHG(2)=DHG(2)-P(IN2+2,1)*DHC
0584               IF(IN2.EQ.IN(2)) DHG(3)=DHG(3)+P(IN1+2,1)*DHC
0585               IF(IN1.EQ.IN(1).AND.IN2.EQ.IN(2)) DHG(4)=DHG(4)-DHC
0586   540       CONTINUE
0587   550     CONTINUE
0588
0589 C...Junction strings: solve (m2, Gamma) equation system for energies.
0590           DHS1=DHM(3)*DHG(4)-DHM(4)*DHG(3)
0591           IF(ABS(DHS1).LT.1D-4) GOTO 360
0592           DHS2=DHM(4)*(GAM(3)-DHG(1))-DHM(2)*DHG(3)-DHG(4)*
0593      &    (P(I,5)**2-DHM(1))+DHG(2)*DHM(3)
0594           DHS3=DHM(2)*(GAM(3)-DHG(1))-DHG(2)*(P(I,5)**2-DHM(1))
0595           P(IN(2)+2,4)=0.5D0*(SQRT(MAX(0D0,DHS2**2-4D0*DHS1*DHS3))/
0596      &    ABS(DHS1)-DHS2/DHS1)
0597           IF(DHM(2)+DHM(4)*P(IN(2)+2,4).LE.0D0) GOTO 360
0598           P(IN(1)+2,4)=(P(I,5)**2-DHM(1)-DHM(3)*P(IN(2)+2,4))/
0599      &    (DHM(2)+DHM(4)*P(IN(2)+2,4))
0600
0601 C...Junction strings: step to new region if necessary.
0602           IF(P(IN(2)+2,4).GT.P(IN(2)+2,3)) THEN
0603             P(IN(2)+2,4)=P(IN(2)+2,3)
0604             P(IN(2)+2,1)=1D0
0605             IN(2)=IN(2)+4
0606             IF(IN(2).GT.N+NR+4*NS) GOTO 360
0607             IF(FOUR(IN(1),IN(2)).LE.1D-2) THEN
0608               P(IN(1)+2,4)=P(IN(1)+2,3)
0609               P(IN(1)+2,1)=0D0
0610               IN(1)=IN(1)+4
0611             ENDIF
0612             GOTO 480
0613           ELSEIF(P(IN(1)+2,4).GT.P(IN(1)+2,3)) THEN
0614             P(IN(1)+2,4)=P(IN(1)+2,3)
0615             P(IN(1)+2,1)=0D0
0616             IN(1)=IN(1)+4
0617             GOTO 480
0618           ENDIF
0619
0620 C...Junction strings: particle four-momentum, remainder, loop back.
0621   560     DO 570 J=1,4
0622             P(I,J)=P(I,J)+P(IN(1)+2,4)*P(IN(1),J)+
0623      &      P(IN(2)+2,4)*P(IN(2),J)
0624             PJU(IU+3,J)=PJU(IU+3,J)+P(I,J)
0625   570     CONTINUE
0626           IF(P(I,4).LT.P(I,5)) GOTO 360
0627           PJU(IU+3,5)=TJU(4)*PJU(IU+3,4)-TJU(1)*PJU(IU+3,1)-
0628      &    TJU(2)*PJU(IU+3,2)-TJU(3)*PJU(IU+3,3)
0629           IF(PJU(IU+3,5).LT.PJU(IU,5)) THEN
0630             KFL(1)=-KFL(3)
0631             PX(1)=-PX(3)
0632             PY(1)=-PY(3)
0633             GAM(1)=GAM(3)
0634             IF(IN(3).NE.IN(6)) THEN
0635               DO 580 J=1,4
0636                 P(IN(6),J)=P(IN(3),J)
0637                 P(IN(6)+1,J)=P(IN(3)+1,J)
0638   580         CONTINUE
0639             ENDIF
0640             DO 590 JQ=1,2
0641               IN(3+JQ)=IN(JQ)
0642               P(IN(JQ)+2,3)=P(IN(JQ)+2,3)-P(IN(JQ)+2,4)
0643               P(IN(JQ)+2,1)=P(IN(JQ)+2,1)-(3-2*JQ)*P(IN(JQ)+2,4)
0644   590       CONTINUE
0645             GOTO 430
0646           ENDIF
0647
0648 C...Junction strings: save quantities left after each string.
0649           IF(IABS(KFL(1)).GT.10) GOTO 360
0650   600     I=I-1
0651           KFJH(IU)=KFL(1)
0652           DO 610 J=1,4
0653             PJU(IU+3,J)=PJU(IU+3,J)-P(I+1,J)
0654   610     CONTINUE
0655
0656 C...Junction strings: loopback if much unused energy in both strings.
0657           PJU(IU+3,5)=TJU(4)*PJU(IU+3,4)-TJU(1)*PJU(IU+3,1)-
0658      &    TJU(2)*PJU(IU+3,2)-TJU(3)*PJU(IU+3,3)
0659           EJSTR(IU)=PJU(IU,5)-PJU(IU+3,5)
0660   620   CONTINUE
0661         IF((MIN(EJSTR(1),EJSTR(2)).GT.PARJ(49).OR.
0662      &  EJSTR(1).GT.PARJ(49)+PYR(0)*PARJ(50).OR.
0663      &  EJSTR(2).GT.PARJ(49)+PYR(0)*PARJ(50))
0664      &  .AND.NTRYER.LT.10) GOTO 320
0665
0666 C...Junction strings: put together to new effective string endpoint.
0667         NJS(JT)=I-ISTA
0668         KFLS=2*INT(PYR(0)+3D0*PARJ(4)/(1D0+3D0*PARJ(4)))+1
0669         IF(KFJH(1).EQ.KFJH(2)) KFLS=3
0670         KFJS(JT)=ISIGN(1000*MAX(IABS(KFJH(1)),IABS(KFJH(2)))+
0671      &  100*MIN(IABS(KFJH(1)),IABS(KFJH(2)))+KFLS,KFJH(1))
0672         DO 630 J=1,4
0673           PJS(JT,J)=PJU(1,J)+PJU(2,J)+P(MJU(JT),J)
0674           PJS(JT+2,J)=PJU(4,J)+PJU(5,J)
0675   630   CONTINUE
0676         PJS(JT,5)=SQRT(MAX(0D0,PJS(JT,4)**2-PJS(JT,1)**2-PJS(JT,2)**2-
0677      &  PJS(JT,3)**2))
0678         PJS(JT+2,5)=0D0
0679   640 CONTINUE
0680
0681 C...Open versus closed strings. Choose breakup region for latter.
0682   650 IF(MJU(1).NE.0.AND.MJU(2).NE.0) THEN
0683         NS=MJU(2)-MJU(1)
0684         NB=MJU(1)-N
0685       ELSEIF(MJU(1).NE.0) THEN
0686         NS=N+NR-MJU(1)
0687         NB=MJU(1)-N
0688       ELSEIF(MJU(2).NE.0) THEN
0689         NS=MJU(2)-N
0690         NB=1
0691       ELSEIF(IABS(K(N+1,2)).NE.21) THEN
0692         NS=NR-1
0693         NB=1
0694       ELSE
0695         NS=NR+1
0696         W2SUM=0D0
0697         DO 660 IS=1,NR
0698           P(N+NR+IS,1)=0.5D0*FOUR(N+IS,N+IS+1-NR*(IS/NR))
0699           W2SUM=W2SUM+P(N+NR+IS,1)
0700   660   CONTINUE
0701         W2RAN=PYR(0)*W2SUM
0702         NB=0
0703   670   NB=NB+1
0704         W2SUM=W2SUM-P(N+NR+NB,1)
0705         IF(W2SUM.GT.W2RAN.AND.NB.LT.NR) GOTO 670
0706       ENDIF
0707
0708 C...Find longitudinal string directions (i.e. lightlike four-vectors).
0709       DO 700 IS=1,NS
0710         IS1=N+IS+NB-1-NR*((IS+NB-2)/NR)
0711         IS2=N+IS+NB-NR*((IS+NB-1)/NR)
0712         DO 680 J=1,5
0713           DP(1,J)=P(IS1,J)
0714           IF(IABS(K(IS1,2)).EQ.21) DP(1,J)=0.5D0*DP(1,J)
0715           IF(IS1.EQ.MJU(1)) DP(1,J)=PJS(1,J)-PJS(3,J)
0716           DP(2,J)=P(IS2,J)
0717           IF(IABS(K(IS2,2)).EQ.21) DP(2,J)=0.5D0*DP(2,J)
0718           IF(IS2.EQ.MJU(2)) DP(2,J)=PJS(2,J)-PJS(4,J)
0719   680   CONTINUE
0720         IF(IS1.EQ.MJU(1)) DP(1,5)=SQRT(MAX(0D0,DP(1,4)**2-DP(1,1)**2-
0721      &  DP(1,2)**2-DP(1,3)**2))
0722         IF(IS2.EQ.MJU(2)) DP(2,5)=SQRT(MAX(0D0,DP(2,4)**2-DP(2,1)**2-
0723      &  DP(2,2)**2-DP(2,3)**2))
0724         DP(3,5)=DFOUR(1,1)
0725         DP(4,5)=DFOUR(2,2)
0726         DHKC=DFOUR(1,2)
0727         IF(DP(3,5)+2D0*DHKC+DP(4,5).LE.0D0) GOTO 200
0728         DHKS=SQRT(DHKC**2-DP(3,5)*DP(4,5))
0729         DHK1=0.5D0*((DP(4,5)+DHKC)/DHKS-1D0)
0730         DHK2=0.5D0*((DP(3,5)+DHKC)/DHKS-1D0)
0731         IN1=N+NR+4*IS-3
0732         P(IN1,5)=SQRT(DP(3,5)+2D0*DHKC+DP(4,5))
0733         DO 690 J=1,4
0734           P(IN1,J)=(1D0+DHK1)*DP(1,J)-DHK2*DP(2,J)
0735           P(IN1+1,J)=(1D0+DHK2)*DP(2,J)-DHK1*DP(1,J)
0736   690   CONTINUE
0737   700 CONTINUE
0738
0739 C...Begin initialization: sum up energy, set starting position.
0740       ISAV=I
0741       MSTU91=MSTU(90)
0742   710 NTRY=NTRY+1
0743       IF(NTRY.GT.100.AND.NTRYR.LE.8.AND.NR.GT.NRMIN) THEN
0744         PARU12=4D0*PARU12
0745         PARU13=2D0*PARU13
0746         GOTO 140
0747       ELSEIF(NTRY.GT.100) THEN
0748         CALL PYERRM(14,'(PYSTRF:) caught in infinite loop')
0749         IF(MSTU(21).GE.1) RETURN
0750       ENDIF
0751       I=ISAV
0752       MSTU(90)=MSTU91
0753       DO 730 J=1,4
0754         P(N+NRS,J)=0D0
0755         DO 720 IS=1,NR
0756           P(N+NRS,J)=P(N+NRS,J)+P(N+IS,J)
0757   720   CONTINUE
0758   730 CONTINUE
0759       DO 750 JT=1,2
0760         IRANK(JT)=0
0761         IF(MJU(JT).NE.0) IRANK(JT)=NJS(JT)
0762         IF(NS.GT.NR) IRANK(JT)=1
0763         IBARRK(JT)=0
0764         IE(JT)=K(N+1+(JT/2)*(NP-1),3)
0765         IN(3*JT+1)=N+NR+1+4*(JT/2)*(NS-1)
0766         IN(3*JT+2)=IN(3*JT+1)+1
0767         IN(3*JT+3)=N+NR+4*NS+2*JT-1
0768         DO 740 IN1=N+NR+2+JT,N+NR+4*NS-2+JT,4
0769           P(IN1,1)=2-JT
0770           P(IN1,2)=JT-1
0771           P(IN1,3)=1D0
0772   740   CONTINUE
0773   750 CONTINUE
0774
0775 C.. MOPS variables and switches
0776       NRVMO=0
0777       XBMO=1D0
0778       MSTU(121)=0
0779       MSTU(122)=0
0780
0781 C...Initialize flavour and pT variables for open string.
0782       IF(NS.LT.NR) THEN
0783         PX(1)=0D0
0784         PY(1)=0D0
0785         IF(NS.EQ.1.AND.MJU(1)+MJU(2).EQ.0) CALL PYPTDI(0,PX(1),PY(1))
0786         PX(2)=-PX(1)
0787         PY(2)=-PY(1)
0788         DO 760 JT=1,2
0789           KFL(JT)=K(IE(JT),2)
0790           IF(MJU(JT).NE.0) KFL(JT)=KFJS(JT)
0791           IF(MJU(JT).NE.0.AND.IABS(KFL(JT)).GT.1000) IBARRK(JT)=1
0792           MSTJ(93)=1
0793           PMQ(JT)=PYMASS(KFL(JT))
0794           GAM(JT)=0D0
0795   760   CONTINUE
0796
0797 C...Closed string: random initial breakup flavour, pT and vertex.
0798       ELSE
0799         KFL(3)=INT(1D0+(2D0+PARJ(2))*PYR(0))*(-1)**INT(PYR(0)+0.5D0)
0800         IBMO=0
0801   770   CALL PYKFDI(KFL(3),0,KFL(1),KDUMP)
0802 C.. Closed string: first vertex diq attempt => enforced second
0803 C.. vertex diq
0804         IF(IABS(KFL(1)).GT.10)THEN
0805            IBMO=1
0806            MSTU(121)=0
0807            GOTO 770
0808         ENDIF
0809         IF(IBMO.EQ.1) MSTU(121)=-1
0810         KFL(2)=-KFL(1)
0811         CALL PYPTDI(KFL(1),PX(1),PY(1))
0812         PX(2)=-PX(1)
0813         PY(2)=-PY(1)
0814         PR3=MIN(25D0,0.1D0*P(N+NR+1,5)**2)
0815   780   CALL PYZDIS(KFL(1),KFL(2),PR3,Z)
0816         ZR=PR3/(Z*P(N+NR+1,5)**2)
0817         IF(ZR.GE.1D0) GOTO 780
0818         DO 790 JT=1,2
0819           MSTJ(93)=1
0820           PMQ(JT)=PYMASS(KFL(JT))
0821           GAM(JT)=PR3*(1D0-Z)/Z
0822           IN1=N+NR+3+4*(JT/2)*(NS-1)
0823           P(IN1,JT)=1D0-Z
0824           P(IN1,3-JT)=JT-1
0825           P(IN1,3)=(2-JT)*(1D0-Z)+(JT-1)*Z
0826           P(IN1+1,JT)=ZR
0827           P(IN1+1,3-JT)=2-JT
0828           P(IN1+1,3)=(2-JT)*(1D0-ZR)+(JT-1)*ZR
0829   790   CONTINUE
0830       ENDIF
0831 C.. MOPS variables
0832       DO 800 JT=1,2
0833          XTMO(JT)=1D0
0834          PM2QMO(JT)=PMQ(JT)**2
0835          IF(IABS(KFL(JT)).GT.10) PM2QMO(JT)=0D0
0836   800 CONTINUE
0837
0838 C...Find initial transverse directions (i.e. spacelike four-vectors).
0839       DO 840 JT=1,2
0840         IF(JT.EQ.1.OR.NS.EQ.NR-1.OR.MJU(1)+MJU(2).NE.0) THEN
0841           IN1=IN(3*JT+1)
0842           IN3=IN(3*JT+3)
0843           DO 810 J=1,4
0844             DP(1,J)=P(IN1,J)
0845             DP(2,J)=P(IN1+1,J)
0846             DP(3,J)=0D0
0847             DP(4,J)=0D0
0848   810     CONTINUE
0849           DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)
0850           DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)
0851           DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4)
0852           DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4)
0853           DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4)
0854           IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1D0
0855           IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1D0
0856           IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1D0
0857           IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1D0
0858           DHC12=DFOUR(1,2)
0859           DHCX1=DFOUR(3,1)/DHC12
0860           DHCX2=DFOUR(3,2)/DHC12
0861           DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12)
0862           DHCY1=DFOUR(4,1)/DHC12
0863           DHCY2=DFOUR(4,2)/DHC12
0864           DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12
0865           DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2)
0866           DO 820 J=1,4
0867             DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J))
0868             P(IN3,J)=DP(3,J)
0869             P(IN3+1,J)=DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)-
0870      &      DHCYX*DP(3,J))
0871   820     CONTINUE
0872         ELSE
0873           DO 830 J=1,4
0874             P(IN3+2,J)=P(IN3,J)
0875             P(IN3+3,J)=P(IN3+1,J)
0876   830     CONTINUE
0877         ENDIF
0878   840 CONTINUE
0879
0880 C...Remove energy used up in junction string fragmentation.
0881       IF(MJU(1)+MJU(2).GT.0) THEN
0882         DO 860 JT=1,2
0883           IF(NJS(JT).EQ.0) GOTO 860
0884           DO 850 J=1,4
0885             P(N+NRS,J)=P(N+NRS,J)-PJS(JT+2,J)
0886   850     CONTINUE
0887   860   CONTINUE
0888         PARJST=PARJ(33)
0889         IF(MSTJ(11).EQ.2) PARJST=PARJ(34)
0890         WMIN=PARJST+PMQ(1)+PMQ(2)
0891         WREM2=FOUR(N+NRS,N+NRS)
0892         IF(P(N+NRS,4).LT.0D0.OR.WREM2.LT.WMIN**2) THEN
0893           NTRYWR=NTRYWR+1
0894           IF(MOD(NTRYWR,20).NE.0) NTRYR=NTRYR-1
0895           GOTO 140
0896         ENDIF
0897       ENDIF
0898
0899 C...Produce new particle: side, origin.
0900   870 I=I+1
0901       IF(2*I-NSAV.GE.MSTU(4)-MSTU(32)-5) THEN
0902         CALL PYERRM(11,'(PYSTRF:) no more memory left in PYJETS')
0903         IF(MSTU(21).GE.1) RETURN
0904       ENDIF
0905 C.. New side priority for popcorn systems
0906       IF(MSTU(121).LE.0)THEN
0907          JT=1.5D0+PYR(0)
0908          IF(IABS(KFL(3-JT)).GT.10) JT=3-JT
0909          IF(IABS(KFL(3-JT)).GE.4.AND.IABS(KFL(3-JT)).LE.8) JT=3-JT
0910       ENDIF
0911       JR=3-JT
0912       JS=3-2*JT
0913       IRANK(JT)=IRANK(JT)+1
0914       K(I,1)=1
0915       K(I,4)=0
0916       K(I,5)=0
0917
0918 C...Generate flavour, hadron and pT.
0919   880 K(I,3)=IE(JT)
0920       CALL PYKFDI(KFL(JT),0,KFL(3),K(I,2))
0921       IF(K(I,2).EQ.0) GOTO 710
0922       MU90MO=MSTU(90)
0923       IF(MSTU(121).EQ.-1) GOTO 910
0924       IF(IRANK(JT).EQ.1.AND.IABS(KFL(JT)).LE.10.AND.
0925      &IABS(KFL(3)).GT.10) THEN
0926         IF(PYR(0).GT.PARJ(19)) GOTO 880
0927       ENDIF
0928       IF(IBARRK(JT).EQ.1.AND.MOD(IABS(K(I,2)),10000).GT.1000)
0929      &K(I,3)=IJUORI(JT)
0930       P(I,5)=PYMASS(K(I,2))
0931       CALL PYPTDI(KFL(JT),PX(3),PY(3))
0932       PR(JT)=P(I,5)**2+(PX(JT)+PX(3))**2+(PY(JT)+PY(3))**2
0933
0934 C...Final hadrons for small invariant mass.
0935       MSTJ(93)=1
0936       PMQ(3)=PYMASS(KFL(3))
0937       PARJST=PARJ(33)
0938       IF(MSTJ(11).EQ.2) PARJST=PARJ(34)
0939       WMIN=PARJST+PMQ(1)+PMQ(2)+PARJ(36)*PMQ(3)
0940       IF(IABS(KFL(JT)).GT.10.AND.IABS(KFL(3)).GT.10) WMIN=
0941      &WMIN-0.5D0*PARJ(36)*PMQ(3)
0942       WREM2=FOUR(N+NRS,N+NRS)
0943       IF(WREM2.LT.0.10D0) GOTO 710
0944       IF(WREM2.LT.MAX(WMIN*(1D0+(2D0*PYR(0)-1D0)*PARJ(37)),
0945      &PARJ(32)+PMQ(1)+PMQ(2))**2) GOTO 1080
0946
0947 C...Choose z, which gives Gamma. Shift z for heavy flavours.
0948       CALL PYZDIS(KFL(JT),KFL(3),PR(JT),Z)
0949       IF(IABS(KFL(JT)).GE.4.AND.IABS(KFL(JT)).LE.8.AND.
0950      &MSTU(90).LT.8) THEN
0951         MSTU(90)=MSTU(90)+1
0952         MSTU(90+MSTU(90))=I
0953         PARU(90+MSTU(90))=Z
0954       ENDIF
0955       KFL1A=IABS(KFL(1))
0956       KFL2A=IABS(KFL(2))
0957       IF(MAX(MOD(KFL1A,10),MOD(KFL1A/1000,10),MOD(KFL2A,10),
0958      &MOD(KFL2A/1000,10)).GE.4) THEN
0959         PR(JR)=(PMQ(JR)+PMQ(3))**2+(PX(JR)-PX(3))**2+(PY(JR)-PY(3))**2
0960         PW12=SQRT(MAX(0D0,(WREM2-PR(1)-PR(2))**2-4D0*PR(1)*PR(2)))
0961         Z=(WREM2+PR(JT)-PR(JR)+PW12*(2D0*Z-1D0))/(2D0*WREM2)
0962         PR(JR)=(PMQ(JR)+PARJST)**2+(PX(JR)-PX(3))**2+(PY(JR)-PY(3))**2
0963         IF((1D0-Z)*(WREM2-PR(JT)/Z).LT.PR(JR)) GOTO 1080
0964       ENDIF
0965       GAM(3)=(1D0-Z)*(GAM(JT)+PR(JT)/Z)
0966
0967 C.. MOPS baryon model modification
0968       XTMO3=(1D0-Z)*XTMO(JT)
0969       IF(IABS(KFL(3)).LE.10) NRVMO=0
0970       IF(IABS(KFL(3)).GT.10.AND.MSTJ(12).GE.4) THEN
0971          GTSTMO=1D0
0972          PTSTMO=1D0
0973          RTSTMO=PYR(0)
0974          IF(IABS(KFL(JT)).LE.10)THEN
0975             XBMO=MIN(XTMO3,1D0-(2D-10))
0976             GBMO=GAM(3)
0977             PMMO=0D0
0978             PGMO=GBMO+LOG(1D0-XBMO)*PM2QMO(JT)
0979             GTSTMO=1D0-PARF(192)**PGMO
0980          ELSE
0981             IF(IRANK(JT).EQ.1) THEN
0982                GBMO=GAM(JT)
0983                PMMO=0D0
0984                XBMO=1D0
0985             ENDIF
0986             IF(XBMO.LT.1D0-(1D-10))THEN
0987                PGNMO=GBMO*XTMO3/XBMO+PM2QMO(JT)*LOG(1D0-XTMO3)
0988                GTSTMO=(1D0-PARF(192)**PGNMO)/(1D0-PARF(192)**PGMO)
0989                PGMO=PGNMO
0990             ENDIF
0991             IF(MSTJ(12).GE.5)THEN
0992                PMNMO=SQRT((XBMO-XTMO3)*(GAM(3)/XTMO3-GBMO/XBMO))
0993                PMMO=PMMO+PMAS(PYCOMP(K(I,2)),1)-PMAS(PYCOMP(K(I,2)),3)
0994                PTSTMO=EXP((PMMO-PMNMO)*PARF(193))
0995                PMMO=PMNMO
0996             ENDIF
0997          ENDIF
0998
0999 C.. MOPS Accepting popcorn system hadron.
1000          IF(PTSTMO*GTSTMO.GT.RTSTMO) THEN
1001             IF(IRANK(JT).EQ.1.OR.IABS(KFL(JT)).LE.10) THEN
1002                NRVMO=I-N-NR
1003                IF(I+NRVMO.GT.MSTU(4)-MSTU(32)-5) THEN
1004                   CALL PYERRM(11,
1005      &                 '(PYSTRF:) no more memory left in PYJETS')
1006                   IF(MSTU(21).GE.1) RETURN
1007                ENDIF
1008                IMO=I
1009                KFLMO=KFL(JT)
1010                PMQMO=PMQ(JT)
1011                PXMO=PX(JT)
1012                PYMO=PY(JT)
1013                GAMMO=GAM(JT)
1014                IRMO=IRANK(JT)
1015                XMO=XTMO(JT)
1016                DO 900 J=1,9
1017                   IF(J.LE.5) THEN
1018                      DO 890 LINE=1,I-N-NR
1019                         P(MSTU(4)-MSTU(32)-LINE,J)=P(N+NR+LINE,J)
1020                         K(MSTU(4)-MSTU(32)-LINE,J)=K(N+NR+LINE,J)
1021   890                CONTINUE
1022                   ENDIF
1023                   INMO(J)=IN(J)
1024   900          CONTINUE
1025             ENDIF
1026          ELSE
1027 C..Reject popcorn system, flag=-1 if enforcing new one
1028             MSTU(121)=-1
1029             IF(PTSTMO.GT.RTSTMO) MSTU(121)=-2
1030          ENDIF
1031       ENDIF
1032
1033
1034 C..Lift restoring string outside MOPS block
1035   910 IF(MSTU(121).LT.0) THEN
1036          IF(MSTU(121).EQ.-2) MSTU(121)=0
1037          MSTU(90)=MU90MO
1038          NRVMO=0
1039          IF(IRANK(JT).EQ.1.OR.IABS(KFL(JT)).LE.10) GOTO 880
1040          I=IMO
1041          KFL(JT)=KFLMO
1042          PMQ(JT)=PMQMO
1043          PX(JT)=PXMO
1044          PY(JT)=PYMO
1045          GAM(JT)=GAMMO
1046          IRANK(JT)=IRMO
1047          XTMO(JT)=XMO
1048          DO 930 J=1,9
1049             IF(J.LE.5) THEN
1050                DO 920 LINE=1,I-N-NR
1051                   P(N+NR+LINE,J)=P(MSTU(4)-MSTU(32)-LINE,J)
1052                   K(N+NR+LINE,J)=K(MSTU(4)-MSTU(32)-LINE,J)
1053   920          CONTINUE
1054             ENDIF
1055             IN(J)=INMO(J)
1056   930    CONTINUE
1057          GOTO 880
1058       ENDIF
1059       XTMO(JT)=XTMO3
1060 C.. MOPS end of modification
1061
1062       DO 940 J=1,3
1063         IN(J)=IN(3*JT+J)
1064   940 CONTINUE
1065
1066 C...Stepping within or from 'low' string region easy.
1067       IF(IN(1)+1.EQ.IN(2).AND.Z*P(IN(1)+2,3)*P(IN(2)+2,3)*
1068      &P(IN(1),5)**2.GE.PR(JT)) THEN
1069         P(IN(JT)+2,4)=Z*P(IN(JT)+2,3)
1070         P(IN(JR)+2,4)=PR(JT)/(P(IN(JT)+2,4)*P(IN(1),5)**2)
1071         DO 950 J=1,4
1072           P(I,J)=(PX(JT)+PX(3))*P(IN(3),J)+(PY(JT)+PY(3))*P(IN(3)+1,J)
1073   950   CONTINUE
1074         GOTO 1040
1075       ELSEIF(IN(1)+1.EQ.IN(2)) THEN
1076         P(IN(JR)+2,4)=P(IN(JR)+2,3)
1077         P(IN(JR)+2,JT)=1D0
1078         IN(JR)=IN(JR)+4*JS
1079         IF(JS*IN(JR).GT.JS*IN(4*JR)) GOTO 710
1080         IF(FOUR(IN(1),IN(2)).LE.1D-2) THEN
1081           P(IN(JT)+2,4)=P(IN(JT)+2,3)
1082           P(IN(JT)+2,JT)=0D0
1083           IN(JT)=IN(JT)+4*JS
1084         ENDIF
1085       ENDIF
1086
1087 C...Find new transverse directions (i.e. spacelike string vectors).
1088   960 IF(JS*IN(1).GT.JS*IN(3*JR+1).OR.JS*IN(2).GT.JS*IN(3*JR+2).OR.
1089      &IN(1).GT.IN(2)) GOTO 710
1090       IF(IN(1).NE.IN(3*JT+1).OR.IN(2).NE.IN(3*JT+2)) THEN
1091         DO 970 J=1,4
1092           DP(1,J)=P(IN(1),J)
1093           DP(2,J)=P(IN(2),J)
1094           DP(3,J)=0D0
1095           DP(4,J)=0D0
1096   970   CONTINUE
1097         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)
1098         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)
1099         DHC12=DFOUR(1,2)
1100         IF(DHC12.LE.1D-2) THEN
1101           P(IN(JT)+2,4)=P(IN(JT)+2,3)
1102           P(IN(JT)+2,JT)=0D0
1103           IN(JT)=IN(JT)+4*JS
1104           GOTO 960
1105         ENDIF
1106         IN(3)=N+NR+4*NS+5
1107         DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4)
1108         DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4)
1109         DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4)
1110         IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1D0
1111         IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1D0
1112         IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1D0
1113         IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1D0
1114         DHCX1=DFOUR(3,1)/DHC12
1115         DHCX2=DFOUR(3,2)/DHC12
1116         DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12)
1117         DHCY1=DFOUR(4,1)/DHC12
1118         DHCY2=DFOUR(4,2)/DHC12
1119         DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12
1120         DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2)
1121         DO 980 J=1,4
1122           DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J))
1123           P(IN(3),J)=DP(3,J)
1124           P(IN(3)+1,J)=DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)-
1125      &    DHCYX*DP(3,J))
1126   980   CONTINUE
1127 C...Express pT with respect to new axes, if sensible.
1128         PXP=-(PX(3)*FOUR(IN(3*JT+3),IN(3))+PY(3)*
1129      &  FOUR(IN(3*JT+3)+1,IN(3)))
1130         PYP=-(PX(3)*FOUR(IN(3*JT+3),IN(3)+1)+PY(3)*
1131      &  FOUR(IN(3*JT+3)+1,IN(3)+1))
1132         IF(ABS(PXP**2+PYP**2-PX(3)**2-PY(3)**2).LT.0.01D0) THEN
1133           PX(3)=PXP
1134           PY(3)=PYP
1135         ENDIF
1136       ENDIF
1137
1138 C...Sum up known four-momentum. Gives coefficients for m2 expression.
1139       DO 1010 J=1,4
1140         DHG(J)=0D0
1141         P(I,J)=PX(JT)*P(IN(3*JT+3),J)+PY(JT)*P(IN(3*JT+3)+1,J)+
1142      &  PX(3)*P(IN(3),J)+PY(3)*P(IN(3)+1,J)
1143         DO 990 IN1=IN(3*JT+1),IN(1)-4*JS,4*JS
1144           P(I,J)=P(I,J)+P(IN1+2,3)*P(IN1,J)
1145   990   CONTINUE
1146         DO 1000 IN2=IN(3*JT+2),IN(2)-4*JS,4*JS
1147           P(I,J)=P(I,J)+P(IN2+2,3)*P(IN2,J)
1148  1000   CONTINUE
1149  1010 CONTINUE
1150       DHM(1)=FOUR(I,I)
1151       DHM(2)=2D0*FOUR(I,IN(1))
1152       DHM(3)=2D0*FOUR(I,IN(2))
1153       DHM(4)=2D0*FOUR(IN(1),IN(2))
1154
1155 C...Find coefficients for Gamma expression.
1156       DO 1030 IN2=IN(1)+1,IN(2),4
1157         DO 1020 IN1=IN(1),IN2-1,4
1158           DHC=2D0*FOUR(IN1,IN2)
1159           DHG(1)=DHG(1)+P(IN1+2,JT)*P(IN2+2,JT)*DHC
1160           IF(IN1.EQ.IN(1)) DHG(2)=DHG(2)-JS*P(IN2+2,JT)*DHC
1161           IF(IN2.EQ.IN(2)) DHG(3)=DHG(3)+JS*P(IN1+2,JT)*DHC
1162           IF(IN1.EQ.IN(1).AND.IN2.EQ.IN(2)) DHG(4)=DHG(4)-DHC
1163  1020   CONTINUE
1164  1030 CONTINUE
1165
1166 C...Solve (m2, Gamma) equation system for energies taken.
1167       DHS1=DHM(JR+1)*DHG(4)-DHM(4)*DHG(JR+1)
1168       IF(ABS(DHS1).LT.1D-4) GOTO 710
1169       DHS2=DHM(4)*(GAM(3)-DHG(1))-DHM(JT+1)*DHG(JR+1)-DHG(4)*
1170      &(P(I,5)**2-DHM(1))+DHG(JT+1)*DHM(JR+1)
1171       DHS3=DHM(JT+1)*(GAM(3)-DHG(1))-DHG(JT+1)*(P(I,5)**2-DHM(1))
1172       P(IN(JR)+2,4)=0.5D0*(SQRT(MAX(0D0,DHS2**2-4D0*DHS1*DHS3))/
1173      &ABS(DHS1)-DHS2/DHS1)
1174       IF(DHM(JT+1)+DHM(4)*P(IN(JR)+2,4).LE.0D0) GOTO 710
1175       P(IN(JT)+2,4)=(P(I,5)**2-DHM(1)-DHM(JR+1)*P(IN(JR)+2,4))/
1176      &(DHM(JT+1)+DHM(4)*P(IN(JR)+2,4))
1177
1178 C...Step to new region if necessary.
1179       IF(P(IN(JR)+2,4).GT.P(IN(JR)+2,3)) THEN
1180         P(IN(JR)+2,4)=P(IN(JR)+2,3)
1181         P(IN(JR)+2,JT)=1D0
1182         IN(JR)=IN(JR)+4*JS
1183         IF(JS*IN(JR).GT.JS*IN(4*JR)) GOTO 710
1184         IF(FOUR(IN(1),IN(2)).LE.1D-2) THEN
1185           P(IN(JT)+2,4)=P(IN(JT)+2,3)
1186           P(IN(JT)+2,JT)=0D0
1187           IN(JT)=IN(JT)+4*JS
1188         ENDIF
1189         GOTO 960
1190       ELSEIF(P(IN(JT)+2,4).GT.P(IN(JT)+2,3)) THEN
1191         P(IN(JT)+2,4)=P(IN(JT)+2,3)
1192         P(IN(JT)+2,JT)=0D0
1193         IN(JT)=IN(JT)+4*JS
1194         GOTO 960
1195       ENDIF
1196
1197 C...Four-momentum of particle. Remaining quantities. Loop back.
1198  1040 DO 1050 J=1,4
1199         P(I,J)=P(I,J)+P(IN(1)+2,4)*P(IN(1),J)+P(IN(2)+2,4)*P(IN(2),J)
1200         P(N+NRS,J)=P(N+NRS,J)-P(I,J)
1201  1050 CONTINUE
1202       IF(P(IN(1)+2,4).GT.1D0+PARU(14).OR.P(IN(1)+2,4).LT.-PARU(14).OR.
1203      &P(IN(2)+2,4).GT.1D0+PARU(14).OR.P(IN(2)+2,4).LT.-PARU(14))
1204      &GOTO 200
1205       IF(P(I,4).LT.P(I,5)) GOTO 710
1206       KFL(JT)=-KFL(3)
1207       PMQ(JT)=PMQ(3)
1208       PX(JT)=-PX(3)
1209       PY(JT)=-PY(3)
1210       GAM(JT)=GAM(3)
1211       IF(IN(3).NE.IN(3*JT+3)) THEN
1212         DO 1060 J=1,4
1213           P(IN(3*JT+3),J)=P(IN(3),J)
1214           P(IN(3*JT+3)+1,J)=P(IN(3)+1,J)
1215  1060   CONTINUE
1216       ENDIF
1217       DO 1070 JQ=1,2
1218         IN(3*JT+JQ)=IN(JQ)
1219         P(IN(JQ)+2,3)=P(IN(JQ)+2,3)-P(IN(JQ)+2,4)
1220         P(IN(JQ)+2,JT)=P(IN(JQ)+2,JT)-JS*(3-2*JQ)*P(IN(JQ)+2,4)
1221  1070 CONTINUE
1222       IF(IBARRK(JT).EQ.1.AND.MOD(IABS(K(I,2)),10000).GT.1000)
1223      &IBARRK(JT)=0
1224       GOTO 870
1225
1226 C...Final hadron: side, flavour, hadron, mass.
1227  1080 I=I+1
1228       K(I,1)=1
1229       K(I,3)=IE(JR)
1230       K(I,4)=0
1231       K(I,5)=0
1232       CALL PYKFDI(KFL(JR),-KFL(3),KFLDMP,K(I,2))
1233       IF(K(I,2).EQ.0) GOTO 710
1234       IF(IBARRK(JT).EQ.1.AND.MOD(IABS(K(I-1,2)),10000).GT.1000)
1235      &IBARRK(JT)=0
1236       IF(IBARRK(JT).EQ.1.AND.MOD(IABS(K(I,2)),10000).GT.1000)
1237      &K(I,3)=IJUORI(JT)
1238       IF(IBARRK(JR).EQ.1.AND.MOD(IABS(K(I,2)),10000).GT.1000)
1239      &K(I,3)=IJUORI(JR)
1240       P(I,5)=PYMASS(K(I,2))
1241       PR(JR)=P(I,5)**2+(PX(JR)-PX(3))**2+(PY(JR)-PY(3))**2
1242
1243 C...Final two hadrons: find common setup of four-vectors.
1244       JQ=1
1245       IF(P(IN(4)+2,3)*P(IN(5)+2,3)*FOUR(IN(4),IN(5)).LT.
1246      &P(IN(7)+2,3)*P(IN(8)+2,3)*FOUR(IN(7),IN(8))) JQ=2
1247       DHC12=FOUR(IN(3*JQ+1),IN(3*JQ+2))
1248       DHR1=FOUR(N+NRS,IN(3*JQ+2))/DHC12
1249       DHR2=FOUR(N+NRS,IN(3*JQ+1))/DHC12
1250       IF(IN(4).NE.IN(7).OR.IN(5).NE.IN(8)) THEN
1251         PX(3-JQ)=-FOUR(N+NRS,IN(3*JQ+3))-PX(JQ)
1252         PY(3-JQ)=-FOUR(N+NRS,IN(3*JQ+3)+1)-PY(JQ)
1253         PR(3-JQ)=P(I+(JT+JQ-3)**2-1,5)**2+(PX(3-JQ)+(2*JQ-3)*JS*
1254      &  PX(3))**2+(PY(3-JQ)+(2*JQ-3)*JS*PY(3))**2
1255       ENDIF
1256
1257 C...Solve kinematics for final two hadrons, if possible.
1258       WREM2=2D0*DHR1*DHR2*DHC12
1259       FD=(SQRT(PR(1))+SQRT(PR(2)))/SQRT(WREM2)
1260       IF(MJU(1)+MJU(2).NE.0.AND.I.EQ.ISAV+2.AND.FD.GE.1D0) GOTO 200
1261       IF(FD.GE.1D0) GOTO 710
1262       FA=WREM2+PR(JT)-PR(JR)
1263       FB=SQRT(MAX(0D0,FA**2-4D0*WREM2*PR(JT)))
1264       PREVCF=PARJ(42)
1265       IF(MSTJ(11).EQ.2) PREVCF=PARJ(39)
1266       PREV=1D0/(1D0+EXP(MIN(50D0,PREVCF*FB*PARJ(40))))
1267       FB=SIGN(FB,JS*(PYR(0)-PREV))
1268       KFL1A=IABS(KFL(1))
1269       KFL2A=IABS(KFL(2))
1270       IF(MAX(MOD(KFL1A,10),MOD(KFL1A/1000,10),MOD(KFL2A,10),
1271      &MOD(KFL2A/1000,10)).GE.6) FB=SIGN(SQRT(MAX(0D0,FA**2-
1272      &4D0*WREM2*PR(JT))),DBLE(JS))
1273       DO 1090 J=1,4
1274         P(I-1,J)=(PX(JT)+PX(3))*P(IN(3*JQ+3),J)+(PY(JT)+PY(3))*
1275      &  P(IN(3*JQ+3)+1,J)+0.5D0*(DHR1*(FA+FB)*P(IN(3*JQ+1),J)+
1276      &  DHR2*(FA-FB)*P(IN(3*JQ+2),J))/WREM2
1277         P(I,J)=P(N+NRS,J)-P(I-1,J)
1278  1090 CONTINUE
1279       IF(P(I-1,4).LT.P(I-1,5).OR.P(I,4).LT.P(I,5)) GOTO 710
1280       DM2F1=P(I-1,4)**2-P(I-1,1)**2-P(I-1,2)**2-P(I-1,3)**2-P(I-1,5)**2
1281       DM2F2=P(I,4)**2-P(I,1)**2-P(I,2)**2-P(I,3)**2-P(I,5)**2
1282       IF(DM2F1.GT.1D-10*P(I-1,4)**2.OR.DM2F2.GT.1D-10*P(I,4)**2) THEN
1283         NTRYFN=NTRYFN+1
1284         IF(NTRYFN.LT.100) GOTO 140
1285         CALL PYERRM(13,'(PYSTRF:) bad energies for final two hadrons')
1286       ENDIF
1287
1288 C...Mark jets as fragmented and give daughter pointers.
1289       N=I-NRS+1
1290       DO 1100 I=NSAV+1,NSAV+NP
1291         IM=K(I,3)
1292         K(IM,1)=K(IM,1)+10
1293         IF(MSTU(16).NE.2) THEN
1294           K(IM,4)=NSAV+1
1295           K(IM,5)=NSAV+1
1296         ELSE
1297           K(IM,4)=NSAV+2
1298           K(IM,5)=N
1299         ENDIF
1300  1100 CONTINUE
1301
1302 C...Document string system. Move up particles.
1303       NSAV=NSAV+1
1304       K(NSAV,1)=11
1305       K(NSAV,2)=92
1306       K(NSAV,3)=IP
1307       K(NSAV,4)=NSAV+1
1308       K(NSAV,5)=N
1309       DO 1110 J=1,4
1310         P(NSAV,J)=DPS(J)
1311         V(NSAV,J)=V(IP,J)
1312  1110 CONTINUE
1313       P(NSAV,5)=SQRT(MAX(0D0,DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2))
1314       V(NSAV,5)=0D0
1315       DO 1130 I=NSAV+1,N
1316         DO 1120 J=1,5
1317           K(I,J)=K(I+NRS-1,J)
1318           P(I,J)=P(I+NRS-1,J)
1319           V(I,J)=0D0
1320  1120   CONTINUE
1321  1130 CONTINUE
1322       MSTU91=MSTU(90)
1323       DO 1140 IZ=MSTU90+1,MSTU91
1324         MSTU9T(IZ)=MSTU(90+IZ)-NRS+1-NSAV+N
1325         PARU9T(IZ)=PARU(90+IZ)
1326  1140 CONTINUE
1327       MSTU(90)=MSTU90
1328
1329 C...Order particles in rank along the chain. Update mother pointer.
1330       DO 1160 I=NSAV+1,N
1331         DO 1150 J=1,5
1332           K(I-NSAV+N,J)=K(I,J)
1333           P(I-NSAV+N,J)=P(I,J)
1334  1150   CONTINUE
1335  1160 CONTINUE
1336       I1=NSAV
1337       DO 1190 I=N+1,2*N-NSAV
1338         IF(K(I,3).NE.IE(1).AND.K(I,3).NE.IJUORI(1)) GOTO 1190
1339         I1=I1+1
1340         DO 1170 J=1,5
1341           K(I1,J)=K(I,J)
1342           P(I1,J)=P(I,J)
1343  1170   CONTINUE
1344         IF(MSTU(16).NE.2) K(I1,3)=NSAV
1345         DO 1180 IZ=MSTU90+1,MSTU91
1346           IF(MSTU9T(IZ).EQ.I) THEN
1347             MSTU(90)=MSTU(90)+1
1348             MSTU(90+MSTU(90))=I1
1349             PARU(90+MSTU(90))=PARU9T(IZ)
1350           ENDIF
1351  1180   CONTINUE
1352  1190 CONTINUE
1353       DO 1220 I=2*N-NSAV,N+1,-1
1354         IF(K(I,3).EQ.IE(1).OR.K(I,3).EQ.IJUORI(1)) GOTO 1220
1355         I1=I1+1
1356         DO 1200 J=1,5
1357           K(I1,J)=K(I,J)
1358           P(I1,J)=P(I,J)
1359  1200   CONTINUE
1360         IF(MSTU(16).NE.2) K(I1,3)=NSAV
1361         DO 1210 IZ=MSTU90+1,MSTU91
1362           IF(MSTU9T(IZ).EQ.I) THEN
1363             MSTU(90)=MSTU(90)+1
1364             MSTU(90+MSTU(90))=I1
1365             PARU(90+MSTU(90))=PARU9T(IZ)
1366           ENDIF
1367  1210   CONTINUE
1368  1220 CONTINUE
1369
1370 C...Boost back particle system. Set production vertices.
1371       IF(MBST.EQ.0) THEN
1372         MSTU(33)=1
1373         CALL PYROBO(NSAV+1,N,0D0,0D0,DPS(1)/DPS(4),DPS(2)/DPS(4),
1374      &  DPS(3)/DPS(4))
1375       ELSE
1376         DO 1230 I=NSAV+1,N
1377           HHPMT=P(I,1)**2+P(I,2)**2+P(I,5)**2
1378           IF(P(I,3).GT.0D0) THEN
1379             HHPEZ=(P(I,4)+P(I,3))*HHBZ
1380             P(I,3)=0.5D0*(HHPEZ-HHPMT/HHPEZ)
1381             P(I,4)=0.5D0*(HHPEZ+HHPMT/HHPEZ)
1382           ELSE
1383             HHPEZ=(P(I,4)-P(I,3))/HHBZ
1384             P(I,3)=-0.5D0*(HHPEZ-HHPMT/HHPEZ)
1385             P(I,4)=0.5D0*(HHPEZ+HHPMT/HHPEZ)
1386           ENDIF
1387  1230   CONTINUE
1388       ENDIF
1389       DO 1250 I=NSAV+1,N
1390         DO 1240 J=1,4
1391           V(I,J)=V(IP,J)
1392  1240   CONTINUE
1393  1250 CONTINUE
1394
1395       RETURN
1396       END