C*********************************************************************
 
C...PYPREP
C...Rearranges partons along strings. Allows small systems
C...to collapse into one or two particles and checks flavours.
 
C... Taken from pythia 6.121 and modified to implement:
C... * Multiple attempts to form two hadrons from a cluster
C... * Anisotropic decay of a cluster
C... * New cluster collapse scheme
C... See list of new parameters below
 
      SUBROUTINE PYPREP(IP)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      INTEGER PYK,PYCHGE,PYCOMP
C...Commonblocks.
      COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
      COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
      COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
      COMMON/PYDAT3/MDCY(500,3),MDME(4000,2),BRAT(4000),KFDP(4000,5)
      SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYDAT3/
C...Parameters.
      COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
C...Local arrays.
      DIMENSION DPS(5),DPC(5),UE(3)
C... New (EN).
      DIMENSION PG(5),E1(3),E2(3),E3(3),E4(3),ED(3),ECL(3)
      DIMENSION PS(5)
 
 
C... ************ New parameters ************
C... MSTJ(16): Use new version?, integer (D=0)
C...     (= 0 : old, = 1 : new).
C... MSTJ(17): Number of tries in decay, integer (D=1)
C... PARJ(27): Parameter in non-isotropic decay of a cluster,
C...     double (D=0D0). Best value if using non-isotropic
C...     decay is PARJ(27)=10D0.
C...     The default of PARJ(27)=0D0 is isotropic decay.
 
C... Set the new parameters (EN)
C... You have two options:
C    * Set the parameters in your main program.
C    * Set the parameters here. They will override anything you
C      set in your main program. We suggest using the new parameters
C
C      MSTJ(16)=1
C      MSTJ(17)=2
C      PARJ(27)=0.D0
C
C      together with the following (old) ones:
C
C      PMAS(1,1) = 0.33D0
C      PMAS(2,1) = 0.33D0
C      PMAS(3,1) = 0.50D0
C      PMAS(4,1) = 1.50D0
C      MSTP(92) = 3
C      PARP(91) = 1.D0
C      PARP(93) = 5.D0
C
C      (See hep-ph/9809226 or LU TP 98-24 for the details)
 
C... Do sanity check on parameters (EN)
      IF (PARJ(27).LT.0.D0) THEN
        write(*,*) 'Warning: PARJ(27) negative'
        WRITE(*,*) 'Now set to default (0.D0)'
        PARJ(27)=0.D0
      END IF
      IF (MSTJ(17).LE.0) THEN
        WRITE(*,*) 'Warning: MSTJ(17)<1'
        WRITE(*,*) 'Now set to default (1)'
        MSTJ(17)=1
      END IF
 
C...Rearrange parton shower product listing along strings: begin loop.
      I1=N
      DO 130 MQGST=1,2
        DO 120 I=MAX(1,IP),N
          IF(K(I,1).NE.3) GOTO 120
          KC=PYCOMP(K(I,2))
          IF(KC.EQ.0) GOTO 120
          KQ=KCHG(KC,2)
          IF(KQ.EQ.0.OR.(MQGST.EQ.1.AND.KQ.EQ.2)) GOTO 120
 
C...Pick up loose string end.
          KCS=4
          IF(KQ*ISIGN(1,K(I,2)).LT.0) KCS=5
          IA=I
          NSTP=0
  100     NSTP=NSTP+1
          IF(NSTP.GT.4*N) THEN
            CALL PYERRM(14,'(PYPREP:) caught in infinite loop')
            RETURN
          ENDIF
 
C...Copy undecayed parton.
          IF(K(IA,1).EQ.3) THEN
            IF(I1.GE.MSTU(4)-MSTU(32)-5) THEN
              CALL PYERRM(11,'(PYPREP:) no more memory left in PYJETS')
              RETURN
            ENDIF
            I1=I1+1
            K(I1,1)=2
            IF(NSTP.GE.2.AND.KCHG(PYCOMP(K(IA,2)),2).NE.2) K(I1,1)=1
            K(I1,2)=K(IA,2)
            K(I1,3)=IA
            K(I1,4)=0
            K(I1,5)=0
            DO 110 J=1,5
              P(I1,J)=P(IA,J)
              V(I1,J)=V(IA,J)
  110       CONTINUE
            K(IA,1)=K(IA,1)+10
            IF(K(I1,1).EQ.1) GOTO 120
          ENDIF
 
C...Go to next parton in colour space.
          IB=IA
          IF(MOD(K(IB,KCS)/MSTU(5)**2,2).EQ.0.AND.MOD(K(IB,KCS),MSTU(5))
     &    .NE.0) THEN
            IA=MOD(K(IB,KCS),MSTU(5))
            K(IB,KCS)=K(IB,KCS)+MSTU(5)**2
            MREV=0
          ELSE
            IF(K(IB,KCS).GE.2*MSTU(5)**2.OR.MOD(K(IB,KCS)/MSTU(5),
     &      MSTU(5)).EQ.0) KCS=9-KCS
            IA=MOD(K(IB,KCS)/MSTU(5),MSTU(5))
            K(IB,KCS)=K(IB,KCS)+2*MSTU(5)**2
            MREV=1
          ENDIF
          IF(IA.LE.0.OR.IA.GT.N) THEN
            CALL PYERRM(12,'(PYPREP:) colour rearrangement failed')
            RETURN
          ENDIF
          IF(MOD(K(IA,4)/MSTU(5),MSTU(5)).EQ.IB.OR.MOD(K(IA,5)/MSTU(5),
     &    MSTU(5)).EQ.IB) THEN
            IF(MREV.EQ.1) KCS=9-KCS
            IF(MOD(K(IA,KCS)/MSTU(5),MSTU(5)).NE.IB) KCS=9-KCS
            K(IA,KCS)=K(IA,KCS)+2*MSTU(5)**2
          ELSE
            IF(MREV.EQ.0) KCS=9-KCS
            IF(MOD(K(IA,KCS),MSTU(5)).NE.IB) KCS=9-KCS
            K(IA,KCS)=K(IA,KCS)+MSTU(5)**2
          ENDIF
          IF(IA.NE.I) GOTO 100
          K(I1,1)=1
  120   CONTINUE
  130 CONTINUE
      N=I1
      IF(MSTJ(14).LT.0) RETURN
 
C...Find lowest-mass colour singlet jet system, OK if above threshold.
      IF(MSTJ(14).EQ.0) GOTO 500
      NS=N
  140 NSIN=N-NS
      PDM=1D0+PARJ(32)
      IC=0
 
      DO 190 I=MAX(1,IP),N
        IF(K(I,1).NE.1.AND.K(I,1).NE.2) THEN
        ELSEIF(K(I,1).EQ.2.AND.IC.EQ.0) THEN
          NSIN=NSIN+1
          IC=I
          DO 150 J=1,4
            DPS(J)=P(I,J)
  150     CONTINUE
          MSTJ(93)=1
          DPS(5)=PYMASS(K(I,2))
        ELSEIF(K(I,1).EQ.2) THEN
          DO 160 J=1,4
            DPS(J)=DPS(J)+P(I,J)
  160     CONTINUE
        ELSEIF(IC.NE.0.AND.KCHG(PYCOMP(K(I,2)),2).NE.0) THEN
          DO 170 J=1,4
            DPS(J)=DPS(J)+P(I,J)
  170     CONTINUE
          MSTJ(93)=1
          DPS(5)=DPS(5)+PYMASS(K(I,2))
          PD=SQRT(MAX(0D0,DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2))-
     &    DPS(5)
          IF(PD.LT.PDM) THEN
            PDM=PD
            DO 180 J=1,5
              DPC(J)=DPS(J)
  180       CONTINUE
            IC1=IC
            IC2=I
          ENDIF
          IC=0
        ELSE
          NSIN=NSIN+1
        ENDIF
  190 CONTINUE
 
      IF(PDM.GE.PARJ(32)) GOTO 500
 
C...Fill small-mass system as cluster.
      NSAV=N
      PECM=SQRT(MAX(0D0,DPC(4)**2-DPC(1)**2-DPC(2)**2-DPC(3)**2))
      K(N+1,1)=11
      K(N+1,2)=91
      K(N+1,3)=IC1
      K(N+1,4)=N+2
      K(N+1,5)=N+3
      P(N+1,1)=DPC(1)
      P(N+1,2)=DPC(2)
      P(N+1,3)=DPC(3)
      P(N+1,4)=DPC(4)
      P(N+1,5)=PECM
 
C...Form two particles from flavours of lowest-mass system, if feasible.
C... We could try this several times if it doesn't work the first (EN)
      NTRY = 0
  200 NTRY = NTRY + 1
      K(N+2,1)=1
      K(N+3,1)=1
      IF(MSTU(16).NE.2) THEN
        K(N+2,3)=N+1
        K(N+3,3)=N+1
      ELSE
        K(N+2,3)=IC1
        K(N+3,3)=IC2
      ENDIF
      K(N+2,4)=0
      K(N+3,4)=0
      K(N+2,5)=0
      K(N+3,5)=0
      IF(IABS(K(IC1,2)).NE.21) THEN
        KC1=PYCOMP(K(IC1,2))
        KC2=PYCOMP(K(IC2,2))
        IF(KC1.EQ.0.OR.KC2.EQ.0) GOTO 500
        KQ1=KCHG(KC1,2)*ISIGN(1,K(IC1,2))
        KQ2=KCHG(KC2,2)*ISIGN(1,K(IC2,2))
        IF(KQ1+KQ2.NE.0) GOTO 500
C...Start with qq, if there is one. Only allow for rank 1 popcorn meson
  210   K1=K(IC1,2)
        IF(IABS(K(IC2,2)).GT.10) K1=K(IC2,2)
        MSTU(125)=0
        CALL PYDCYK(K1,0,KFLN,K(N+2,2))
        CALL PYDCYK(K(IC1,2)+K(IC2,2)-K1,-KFLN,KFLDMP,K(N+3,2))
        IF(K(N+2,2).EQ.0.OR.K(N+3,2).EQ.0) GOTO 210
      ELSE
        IF(IABS(K(IC2,2)).NE.21) GOTO 500
C...No room for popcorn mesons in closed string -> 2 hadrons.
        MSTU(125)=0
  220   CALL PYDCYK(1+INT((2D0+PARJ(2))*PYR(0)),0,KFLN,KFDMP)
        CALL PYDCYK(KFLN,0,KFLM,K(N+2,2))
        CALL PYDCYK(-KFLN,-KFLM,KFLDMP,K(N+3,2))
        IF(K(N+2,2).EQ.0.OR.K(N+3,2).EQ.0) GOTO 220
      ENDIF
      P(N+2,5)=PYMASS(K(N+2,2))
      P(N+3,5)=PYMASS(K(N+3,2))
      IF(P(N+2,5)+P(N+3,5)+PARJ(64).GE.PECM.AND.NSIN.EQ.1) GOTO 500
      IF(P(N+2,5)+P(N+3,5)+PARJ(64).GE.PECM) THEN
C... Try again? (EN)
        IF(NTRY.LT.MSTJ(17)) THEN
          GOTO 200
        ELSE
          GOTO 280
        END IF
      END IF
 
C...Perform two-particle decay of jet system, if possible.
C... Non-isotropic decay presently only for charm. (EN)
      ICHARM=0
      IF (PARJ(27).GT.0.D0) THEN
        IF (ABS(K(IC1,2)).EQ.4) THEN
          ICHARM=IC1
        ELSE IF (ABS(K(IC2,2)).EQ.4) THEN
          ICHARM=IC2
        END IF
      END IF
 
      IF(PECM.GE.0.02D0*DPC(4)) THEN
C... Isotropic decay. (EN)
        IF (ICHARM.EQ.0) THEN
          PA=SQRT((PECM**2-(P(N+2,5)+P(N+3,5))**2)*(PECM**2-
     &    (P(N+2,5)-P(N+3,5))**2))/(2D0*PECM)
          UE(3)=2D0*PYR(0)-1D0
          PHI=PARU(2)*PYR(0)
          UE(1)=SQRT(1D0-UE(3)**2)*COS(PHI)
          UE(2)=SQRT(1D0-UE(3)**2)*SIN(PHI)
          DO 230 J=1,3
            P(N+2,J)=PA*UE(J)
            P(N+3,J)=-PA*UE(J)
  230     CONTINUE
          P(N+2,4)=SQRT(PA**2+P(N+2,5)**2)
          P(N+3,4)=SQRT(PA**2+P(N+3,5)**2)
          MSTU(33)=1
          CALL PYROBO(N+2,N+3,0.D0,0.D0,DPC(1)/DPC(4),DPC(2)/DPC(4),
     &    DPC(3)/DPC(4))
        ELSE
C... Anisotropic decay for charm (EN)
C... Boost charm quark to restframe of the cluster (EN)
          CALL PYROBO(ICHARM,ICHARM,0D0,0D0,-DPC(1)/DPC(4),
     &    -DPC(2)/DPC(4),-DPC(3)/DPC(4))
C... Calculate angles (EN)
          THE=PYANGL(P(ICHARM,3),SQRT(P(ICHARM,1)**2+P(ICHARM,2)**2))
          PHI=PYANGL(P(ICHARM,1),P(ICHARM,2)**2)
C... Perform anisotropic cluster decay in restframe of cluster,
C... along charm momentum. (EN)
          AA=PARJ(27)*(P(N+1,5)-P(N+2,5)-P(N+3,5))
          PA=SQRT((PECM**2-(P(N+2,5)+P(N+3,5))**2)*(PECM**2-
     &    (P(N+2,5)-P(N+3,5))**2))/(2D0*PECM)
C... Generate the distribution dN/dx \prop exp(AA*(x-1)),
C... where x=cos(THETA)
          UE(3)=1+(1/AA)*log(PYR(0)*(1-exp(-2*AA))+exp(-2*AA))
          PHI2=PARU(2)*PYR(0)
          UE(1)=SQRT(1D0-UE(3)**2)*COS(PHI2)
          UE(2)=SQRT(1D0-UE(3)**2)*SIN(PHI2)
          DO 240 J=1,3
            P(N+2,J)=PA*UE(J)
            P(N+3,J)=-PA*UE(J)
  240     CONTINUE
          P(N+2,4)=SQRT(PA**2+P(N+2,5)**2)
          P(N+3,4)=SQRT(PA**2+P(N+3,5)**2)
          MSTU(33)=1
C... Boost everything back to original system.
          CALL PYROBO(N+2,N+3,THE,PHI,DPC(1)/DPC(4),DPC(2)/DPC(4),
     &    DPC(3)/DPC(4))
          CALL PYROBO(ICHARM,ICHARM,0D0,0D0,DPC(1)/DPC(4),DPC(2)/DPC(4),
     &    DPC(3)/DPC(4))
        END IF
      ELSE
        NP=0
        DO 250 I=IC1,IC2
          IF(K(I,1).EQ.1.OR.K(I,1).EQ.2) NP=NP+1
  250   CONTINUE
        HA=P(IC1,4)*P(IC2,4)-P(IC1,1)*P(IC2,1)-P(IC1,2)*P(IC2,2)-
     &  P(IC1,3)*P(IC2,3)
        IF(NP.GE.3.OR.HA.LE.1.25D0*P(IC1,5)*P(IC2,5)) THEN
C... Try again? (EN)
          IF(NTRY.LT.MSTJ(17)) THEN
            GOTO 200
          ELSE
            GOTO 280
          END IF
        END IF
        HD1=0.5D0*(P(N+2,5)**2-P(IC1,5)**2)
        HD2=0.5D0*(P(N+3,5)**2-P(IC2,5)**2)
        HR=SQRT(MAX(0D0,((HA-HD1-HD2)**2-(P(N+2,5)*P(N+3,5))**2)/
     &  (HA**2-(P(IC1,5)*P(IC2,5))**2)))-1D0
        HC=P(IC1,5)**2+2D0*HA+P(IC2,5)**2
        HK1=((P(IC2,5)**2+HA)*HR+HD1-HD2)/HC
        HK2=((P(IC1,5)**2+HA)*HR+HD2-HD1)/HC
        DO 260 J=1,4
          P(N+2,J)=(1D0+HK1)*P(IC1,J)-HK2*P(IC2,J)
          P(N+3,J)=(1D0+HK2)*P(IC2,J)-HK1*P(IC1,J)
  260   CONTINUE
      ENDIF
      DO 270 J=1,4
        V(N+1,J)=V(IC1,J)
        V(N+2,J)=V(IC1,J)
        V(N+3,J)=V(IC2,J)
  270 CONTINUE
      V(N+1,5)=0D0
      V(N+2,5)=0D0
      V(N+3,5)=0D0
      N=N+3
      GOTO 480
 
 
C...Else form one particle from the flavours available, if possible.
  280 K(N+1,5)=N+2
      IF(IABS(K(IC1,2)).GT.100.AND.IABS(K(IC2,2)).GT.100) THEN
        GOTO 500
      ELSEIF(IABS(K(IC1,2)).NE.21) THEN
        CALL PYKFDI(K(IC1,2),K(IC2,2),KFLDMP,K(N+2,2))
      ELSE
        KFLN=1+INT((2D0+PARJ(2))*PYR(0))
        CALL PYKFDI(KFLN,-KFLN,KFLDMP,K(N+2,2))
      ENDIF
      IF(K(N+2,2).EQ.0) GOTO 280
      P(N+2,5)=PYMASS(K(N+2,2))
 
C... Use old algorithm for E/p conservation? (EN)
      IF (MSTJ(16).EQ.0) GOTO 440
 
C... Find the stringpiece closest to the cluster. (EN)
C    IS1 = Row of current string beginning
C    IS2 = Row of current string end
C    NPI(ECS) = Number of partons in the current (closest) string
C    IPI(ECE) = Current string piece (closest to cluster)
C    row N+1: Cluster
C    row N+2: Hadron
C    MINF   :flag
C    DGLOMI : Global minimum
      DGLOMI=1.0D30
      MINF=0
      IS1=0
      IS2=0
      NPIECS=0
      IPIECE=0
C... Loop over the undecayed partons not in present cluster (EN)
      DO 330 I=N,MAX(2,IP+1),-1
C... not a parton in cluster ?
        IF (I.LT.IC1.OR.I.GT.IC2) THEN
C... Beginning of a string
          IF (K(I,1).EQ.1.AND.K(I-1,1).EQ.2) THEN
            IS1TMP=I
            NPI=1
            IPI=I
C... Continuing string
          ELSE IF (K(I,1).EQ.2.AND.K(I-1,1).EQ.2) THEN
            NPI=NPI+1
            IPI=I
C... End of a string
          ELSE IF (K(I,1).EQ.2.AND.K(I-1,1).NE.2) THEN
            NPI=NPI+1
            IF (MINF.eq.1) THEN
              NPIECS=NPI
              MINF=0
            END IF
            GOTO 330
          ELSE
            GOTO 330
          END IF
 
C... Now, is this piece the closest so far? Well let's find out...
C... Define normalized vectors e3 and e4 (e1 and e2 are auxiliary)
          DO 290 J=1,3
            E1(J)=P(IPI,J)/P(IPI,4)
            E2(J)=P(IPI-1,J)/P(IPI-1,4)
            ECL(J)=P(N+1,J)/P(N+1,4)
            E3(J)=E2(J)-E1(J)
            E4(J)=ECL(J)-E1(J)
  290     CONTINUE
 
C... Calculate alfa that minimizes D=(e4-alfa*e3)**2
          ALFA=(E3(1)*E4(1)+E3(2)*E4(2)+E3(3)*E4(3))/
     &         (E3(1)**2+E3(2)**2+E3(3)**2)
C... Calculate minimal distance to cluster
C... Three possibilities: dD/da=0 (0<ALFA<1) ...
          IF (ALFA.LT.1D0.AND.ALFA.GT.0D0) THEN
            DO 300 J=1,3
              ED(J) = E4(J)-ALFA*E3(J)
  300       CONTINUE
            DDMIN=ED(1)**2+ED(2)**2+ED(3)**2
          ELSE
            DDMIN=1.D30
          END IF
C... or ALFA=0 and 1
          DO 320 JALFA=0,1
            DO 310 J=1,3
              ED(J) = E4(J)-JALFA*E3(J)
  310       CONTINUE
            DD=ED(1)**2+ED(2)**2+ED(3)**2
            IF(DD.LT.DDMIN) DDMIN=DD
  320     CONTINUE
C... Is this the smallest so far?
          IF(DDMIN.LT.DGLOMI) THEN
            DGLOMI=DDMIN
            IPIECE = IPI
            MINF = 1
            IS1 = IS1TMP
          END IF
        END IF
  330 CONTINUE
      IS2=IS1-NPIECS+1
 
C... Check if there are any strings to connect to the new gluon. (EN)
      IF (IPIECE.EQ.0) GOTO 440
 
C... Check the sign of the gluon 'mass'.Mg=Mclus-Mhad (EN)
      PG(5)=P(N+1,5)-P(N+2,5)
      IF (PG(5).LT.0.D0) THEN
 
C... Have to absorb a 'gluon' instead.
        I1=IPIECE
        I2=IPIECE-1
 
C... Calculate alfa and beta.
        ALFA=P(N+1,4)*P(I2,4)
     &       -P(N+1,3)*P(I2,3)-P(N+1,2)*P(I2,2)-P(N+1,1)*P(I2,1)
        ALFA=ALFA/(P(I1,4)*P(I2,4)
     &       -P(I1,3)*P(I2,3)-P(I1,2)*P(I2,2)-P(I1,1)*P(I2,1))
        BETA=P(N+1,4)*P(I1,4)
     &       -P(N+1,3)*P(I1,3)-P(N+1,2)*P(I1,2)-P(N+1,1)*P(I1,1)
        BETA=BETA/(P(I1,4)*P(I2,4)
     &       -P(I1,3)*P(I2,3)-P(I1,2)*P(I2,2)-P(I1,1)*P(I2,1))
 
        IF (ALFA.LT.0.D0.OR.BETA.LT.0.D0) THEN
C... Let the string ends emit gluon instead. (EN)
          I1=IS1
          I2=IS2
 
C... Calculate alfa and beta. (EN)
          ALFA=P(N+1,4)*P(I2,4)
     &       -P(N+1,3)*P(I2,3)-P(N+1,2)*P(I2,2)-P(N+1,1)*P(I2,1)
          ALFA=ALFA/(P(I1,4)*P(I2,4)
     &       -P(I1,3)*P(I2,3)-P(I1,2)*P(I2,2)-P(I1,1)*P(I2,1))
          BETA=P(N+1,4)*P(I1,4)
     &       -P(N+1,3)*P(I1,3)-P(N+1,2)*P(I1,2)-P(N+1,1)*P(I1,1)
          BETA=BETA/(P(I1,4)*P(I2,4)
     &       -P(I1,3)*P(I2,3)-P(I1,2)*P(I2,2)-P(I1,1)*P(I2,1))
 
C... Will probably never happen. (EN)
          IF(ALFA.LT.0.D0.OR.BETA.LT.0.D0) GOTO 440
        END IF
 
C... Calc. Pg, a part of which will be added to Phad later. (EN)
        DO 340 I=1,4
          PG(I)=ALFA*P(I1,I)+BETA*P(I2,I)
  340   CONTINUE
        PG(5)=SQRT(MAX(0.D0,PG(4)**2-PG(3)**2-PG(2)**2-PG(1)**2))
 
C... Will probably never happen. (EN)
        IF (PG(5).eq.0.D0) GOTO 440
 
C... Calc. the two deltas. 2nd order equation (EN)
        PCLPG=P(N+1,4)*PG(4)
     &       -P(N+1,3)*PG(3)-P(N+1,2)*PG(2)-P(N+1,1)*PG(1)
        PCLPG=PCLPG/(PG(5)**2)
        ROOT=SQRT(PCLPG**2+(P(N+2,5)**2-P(N+1,5)**2)/(PG(5)**2))
 
C... Use the best (smallest) solution.
        DELTA=-PCLPG+ROOT
 
C... Construct the collapsed hadron: Phad=Pclus + DELTA*Pg.
        DO 350 I=1,4
          P(N+2,I)=P(N+1,I)+DELTA*PG(I)
  350   CONTINUE
 
C... Now, this is the gluon that will be added to the string
C... in order to conserve energy/momentum.
        DO 360 I=1,5
          PG(I) = -1.D0*DELTA*PG(I)
  360   CONTINUE
 
C... Calc. new string-momenta as a check (EN)
        DO 380 I=1,4
          PS(I)=PG(I)
          DO 370 J=IS2,IS1
            PS(I)=PS(I)+P(J,I)
  370     CONTINUE
  380   CONTINUE
 
C... Check new string energy (EN)
C    Must not be smaller than the rest-energy (mass) of the string ends
        IF (PS(4)-P(IS1,5)-P(IS2,5).LE.0.D0) THEN
          write(*,*) 'E(string)-Mq-Mqbar < 0', PS(4), P(IS1,5)+P(IS2,5)
          GOTO 440
        END IF
 
C... Check that new string mass is larger than the mass of the
C    quark endpoints.
        PS(5)=SQRT(MAX(0.D0,PS(4)**2-PS(3)**2-PS(2)**2-PS(1)**2))
        IF (PS(5)-P(IS1,5)-P(IS2,5).LE.0.D0) GOTO 440
 
      ELSE
C... dm>0: Emit a 'gluon' (EN)
C... Find "gluon" that is needed to put hadron on the mass shell.
        DO 390 I=1, 4
          P(N+2,I)=(P(N+2,5)/P(N+1,5))*P(N+1,I)
          PG(I)=(PG(5)/P(N+1,5))*P(N+1,I)
  390  CONTINUE
      END IF
 
C... Copy string with new gluon put in (EN)
C... Loop over the string to copy
      ICOUNT=0
      DO 420 I=IS2,IS1
        IF (I.EQ.IPIECE) THEN
C... This is where g is inserted
          DO 400 J=1,5
            P(N+3+ICOUNT,J)=PG(J)
  400     CONTINUE
          K(N+3+ICOUNT,1)=2
          K(N+3+ICOUNT,2)=21
          K(N+3+ICOUNT,3)=N+1
          K(N+3+ICOUNT,4)=0
          K(N+3+ICOUNT,5)=0
          K(N+1,5)=N+3+ICOUNT
          ICOUNT=ICOUNT+1
        END IF
 
C... Insert  the old partons.
        K(N+3+ICOUNT,1)=K(I,1)
        K(N+3+ICOUNT,2)=K(I,2)
C... Parent particle
        K(N+3+ICOUNT,3)=I
        K(N+3+ICOUNT,4)=0
        K(N+3+ICOUNT,5)=0
        DO 410 J=1,5
          P(N+3+ICOUNT,J)=P(I,J)
          V(N+3+ICOUNT,J)=V(I,J)
  410   CONTINUE
 
C... Fix the old partons.
C... Mark as decayed
        K(I,1)=K(I,1)+10
C... First and last daughters are the same.
        K(I,4)=N+3+ICOUNT
        K(I,5)=N+3+ICOUNT
        ICOUNT=ICOUNT+1
  420 CONTINUE
 
C...Mark collapsed system and store daughter pointers. Iterate.
      DO 430 I=IC1,IC2
        IF((K(I,1).EQ.1.OR.K(I,1).EQ.2).AND.KCHG(PYCOMP(K(I,2)),2).NE.0)
     &  THEN
          K(I,1)=K(I,1)+10
          IF(MSTU(16).NE.2) THEN
            K(I,4)=NSAV+1
            K(I,5)=NSAV+1
          ELSE
            K(I,4)=NSAV+2
            K(I,5)=N
         ENDIF
        ENDIF
  430 CONTINUE
 
      N=N+2+ICOUNT
      GOTO 140
 
  440 CONTINUE
 
C... Use old algorithm
C...Find parton/particle which combines to largest extra mass.
      IR=0
      HA=0D0
      HSM=0D0
      DO 460 MCOMB=1,3
        IF(IR.NE.0) GOTO 460
        DO 450 I=MAX(1,IP),N
          IF(K(I,1).LE.0.OR.K(I,1).GT.10.OR.(I.GE.IC1.AND.I.LE.IC2
     &    .AND.K(I,1).GE.1.AND.K(I,1).LE.2)) GOTO 450
          IF(MCOMB.EQ.1) KCI=PYCOMP(K(I,2))
          IF(MCOMB.EQ.1.AND.KCI.EQ.0) GOTO 450
          IF(MCOMB.EQ.1.AND.KCHG(KCI,2).EQ.0.AND.I.LE.NS) GOTO 450
          IF(MCOMB.EQ.2.AND.IABS(K(I,2)).GT.10.AND.IABS(K(I,2)).LE.100)
     &    GOTO 450
          HCR=DPC(4)*P(I,4)-DPC(1)*P(I,1)-DPC(2)*P(I,2)-DPC(3)*P(I,3)
          HSR=2D0*HCR+PECM**2-P(N+2,5)**2-2D0*P(N+2,5)*P(I,5)
          IF(HSR.GT.HSM) THEN
            IR=I
            HA=HCR
            HSM=HSR
          ENDIF
  450   CONTINUE
  460 CONTINUE
 
C...Shuffle energy and momentum to put new particle on mass shell.
      IF(IR.NE.0) THEN
        HB=PECM**2+HA
        HC=P(N+2,5)**2+HA
        HD=P(IR,5)**2+HA
        HK2=0.5D0*(HB*SQRT(MAX(0D0,((HB+HC)**2-4D0*(HB+HD)*P(N+2,5)**2)/
     &  (HA**2-(PECM*P(IR,5))**2)))-(HB+HC))/(HB+HD)
        HK1=(0.5D0*(P(N+2,5)**2-PECM**2)+HD*HK2)/HB
        DO 470 J=1,4
          P(N+2,J)=(1D0+HK1)*DPC(J)-HK2*P(IR,J)
          P(IR,J)=(1D0+HK2)*P(IR,J)-HK1*DPC(J)
          V(N+1,J)=V(IC1,J)
          V(N+2,J)=V(IC1,J)
  470   CONTINUE
        V(N+1,5)=0D0
        V(N+2,5)=0D0
        N=N+2
      ELSE
        CALL PYERRM(3,'(PYPREP:) no match for collapsing cluster')
        RETURN
      ENDIF
 
C...Mark collapsed system and store daughter pointers. Iterate.
  480 DO 490 I=IC1,IC2
        IF((K(I,1).EQ.1.OR.K(I,1).EQ.2).AND.KCHG(PYCOMP(K(I,2)),2).NE.0)
     &  THEN
          K(I,1)=K(I,1)+10
          IF(MSTU(16).NE.2) THEN
            K(I,4)=NSAV+1
            K(I,5)=NSAV+1
          ELSE
            K(I,4)=NSAV+2
            K(I,5)=N
          ENDIF
        ENDIF
  490 CONTINUE
 
 
      IF(N.LT.MSTU(4)-MSTU(32)-5) GOTO 140
 
C...Check flavours and invariant masses in parton systems.
  500 NP=0
      KFN=0
      KQS=0
      DO 510 J=1,5
        DPS(J)=0D0
  510 CONTINUE
      DO 540 I=MAX(1,IP),N
        IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 540
        KC=PYCOMP(K(I,2))
        IF(KC.EQ.0) GOTO 540
        KQ=KCHG(KC,2)*ISIGN(1,K(I,2))
        IF(KQ.EQ.0) GOTO 540
        NP=NP+1
        IF(KQ.NE.2) THEN
          KFN=KFN+1
          KQS=KQS+KQ
          MSTJ(93)=1
          DPS(5)=DPS(5)+PYMASS(K(I,2))
        ENDIF
        DO 520 J=1,4
          DPS(J)=DPS(J)+P(I,J)
  520   CONTINUE
        IF(K(I,1).EQ.1) THEN
          IF(NP.NE.1.AND.(KFN.EQ.1.OR.KFN.GE.3.OR.KQS.NE.0)) CALL
     &    PYERRM(2,'(PYPREP:) unphysical flavour combination')
          IF(NP.NE.1.AND.DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2.LT.
     &    (0.9D0*PARJ(32)+DPS(5))**2) THEN
            CALL PYERRM(3,'(PYPREP:) too small mass in jet system')
          END IF
          NP=0
          KFN=0
          KQS=0
          DO 530 J=1,5
            DPS(J)=0D0
  530     CONTINUE
        ENDIF
  540 CONTINUE
 
      RETURN
      END