 
C*********************************************************************
 
C...PYPREP
C...Rearranges partons along strings. Allows small systems
C...to collapse into one or two particles and checks flavours.
 
      SUBROUTINE PYPREP(IP)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      IMPLICIT INTEGER(I-N)
      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...Local arrays.
      DIMENSION DPS(5),DPC(5),UE(3)
 
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 320
      NS=N
  140 NSIN=N-NS
      PDM=1D0+PARJ(32)
      IC=0
      DO 190 I=MAX(1,IP),NS
        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 320
 
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.
      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 320
        KQ1=KCHG(KC1,2)*ISIGN(1,K(IC1,2))
        KQ2=KCHG(KC2,2)*ISIGN(1,K(IC2,2))
        IF(KQ1+KQ2.NE.0) GOTO 320
C.. Start with qq, if there is one. Only allow for rank 1 popcorn meson
  200   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 200
      ELSE
        IF(IABS(K(IC2,2)).NE.21) GOTO 320
C.. No room for popcorn mesons in closed string -> 2 hadrons.
        MSTU(125)=0
  210   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 210
      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 320
      IF(P(N+2,5)+P(N+3,5)+PARJ(64).GE.PECM) GOTO 260
 
C...Perform two-particle decay of jet system, if possible.
      IF(PECM.GE.0.02D0*DPC(4)) 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 220 J=1,3
          P(N+2,J)=PA*UE(J)
          P(N+3,J)=-PA*UE(J)
  220   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,0D0,0D0,DPC(1)/DPC(4),DPC(2)/DPC(4),
     &  DPC(3)/DPC(4))
      ELSE
        NP=0
        DO 230 I=IC1,IC2
          IF(K(I,1).EQ.1.OR.K(I,1).EQ.2) NP=NP+1
  230   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)) GOTO 260
        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 240 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)
  240   CONTINUE
      ENDIF
      DO 250 J=1,4
        V(N+1,J)=V(IC1,J)
        V(N+2,J)=V(IC1,J)
        V(N+3,J)=V(IC2,J)
  250 CONTINUE
      V(N+1,5)=0D0
      V(N+2,5)=0D0
      V(N+3,5)=0D0
      N=N+3
      GOTO 300
 
C...Else form one particle from the flavours available, if possible.
  260 K(N+1,5)=N+2
      IF(IABS(K(IC1,2)).GT.100.AND.IABS(K(IC2,2)).GT.100) THEN
        GOTO 320
      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 260
      P(N+2,5)=PYMASS(K(N+2,2))
 
C...Find parton/particle which combines to largest extra mass.
      IR=0
      HA=0D0
      HSM=0D0
      DO 280 MCOMB=1,3
        IF(IR.NE.0) GOTO 280
        DO 270 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 270
          IF(MCOMB.EQ.1) KCI=PYCOMP(K(I,2))
          IF(MCOMB.EQ.1.AND.KCI.EQ.0) GOTO 270
          IF(MCOMB.EQ.1.AND.KCHG(KCI,2).EQ.0.AND.I.LE.NS) GOTO 270
          IF(MCOMB.EQ.2.AND.IABS(K(I,2)).GT.10.AND.IABS(K(I,2)).LE.100)
     &    GOTO 270
          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
  270   CONTINUE
  280 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 290 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)
  290   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.
  300 DO 310 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
  310 CONTINUE
      IF(N.LT.MSTU(4)-MSTU(32)-5) GOTO 140
 
C...Check flavours and invariant masses in parton systems.
  320 NP=0
      KFN=0
      KQS=0
      DO 330 J=1,5
        DPS(J)=0D0
  330 CONTINUE
      DO 360 I=MAX(1,IP),N
        IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 360
        KC=PYCOMP(K(I,2))
        IF(KC.EQ.0) GOTO 360
        KQ=KCHG(KC,2)*ISIGN(1,K(I,2))
        IF(KQ.EQ.0) GOTO 360
        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 340 J=1,4
          DPS(J)=DPS(J)+P(I,J)
  340   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) CALL PYERRM(3,
     &    '(PYPREP:) too small mass in jet system')
          NP=0
          KFN=0
          KQS=0
          DO 350 J=1,5
            DPS(J)=0D0
  350     CONTINUE
        ENDIF
  360 CONTINUE
 
      RETURN
      END
