      SUBROUTINE PYCLUS(NJET)
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)
      SAVE /PYJETS/,/PYDAT1/,/PYDAT2/
C...Local arrays and saved variables.
      DIMENSION PS(5)
      SAVE NSAV,NP,PS,PSS,RINIT,NPRE,NREM
C...Functions: distance measure in pT, (pseudo)mass or Durham pT.
      R2T(I1,I2)=(P(I1,5)*P(I2,5)-P(I1,1)*P(I2,1)-P(I1,2)*P(I2,2)-
     &P(I1,3)*P(I2,3))*2D0*P(I1,5)*P(I2,5)/(0.0001D0+P(I1,5)+P(I2,5))**2
      R2M(I1,I2)=2D0*P(I1,4)*P(I2,4)*(1D0-(P(I1,1)*P(I2,1)+P(I1,2)*
     &P(I2,2)+P(I1,3)*P(I2,3))/(P(I1,5)*P(I2,5)))
      R2D(I1,I2)=2D0*MIN(P(I1,4),P(I2,4))**2*(1D0-(P(I1,1)*P(I2,1)+
     &P(I1,2)*P(I2,2)+P(I1,3)*P(I2,3))/(P(I1,5)*P(I2,5)))
C...If first time, reset. If reentering, skip preliminaries.
      IF(MSTU(48).LE.0) THEN
        NP=0
        DO 100 J=1,5
          PS(J)=0D0
  100   CONTINUE
        PSS=0D0
        PIMASS=PMAS(PYCOMP(211),1)
      ELSE
        NJET=NSAV
        IF(MSTU(43).GE.2) N=N-NJET
        DO 110 I=N+1,N+NJET
          P(I,5)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
  110   CONTINUE
        IF(MSTU(46).LE.3.OR.MSTU(46).EQ.5) THEN
          R2ACC=PARU(44)**2
        ELSE
          R2ACC=PARU(45)*PS(5)**2
        ENDIF
        NLOOP=0
        GOTO 300
      ENDIF
C...Find which particles are to be considered in cluster search.
      DO 140 I=1,N
        IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 140
        IF(MSTU(41).GE.2) THEN
          KC=PYCOMP(K(I,2))
          IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
     &    KC.EQ.18) GOTO 140
          IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.PYCHGE(K(I,2)).EQ.0)
     &    GOTO 140
        ENDIF
        IF(N+2*NP.GE.MSTU(4)-MSTU(32)-5) THEN
          CALL PYERRM(11,'(PYCLUS:) no more memory left in PYJETS')
          NJET=-1
          RETURN
        ENDIF
C...Take copy of these particles, with space left for jets later on.
        NP=NP+1
        K(N+NP,3)=I
        DO 120 J=1,5
          P(N+NP,J)=P(I,J)
  120   CONTINUE
        IF(MSTU(42).EQ.0) P(N+NP,5)=0D0
        IF(MSTU(42).EQ.1.AND.K(I,2).NE.22) P(N+NP,5)=PIMASS
        P(N+NP,4)=SQRT(P(N+NP,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
        P(N+NP,5)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
        DO 130 J=1,4
          PS(J)=PS(J)+P(N+NP,J)
  130   CONTINUE
        PSS=PSS+P(N+NP,5)
  140 CONTINUE
      DO 160 I=N+1,N+NP
        K(I+NP,3)=K(I,3)
        DO 150 J=1,5
          P(I+NP,J)=P(I,J)
  150   CONTINUE
  160 CONTINUE
      PS(5)=SQRT(MAX(0D0,PS(4)**2-PS(1)**2-PS(2)**2-PS(3)**2))
C...Very low multiplicities not considered.
      IF(NP.LT.MSTU(47)) THEN
        CALL PYERRM(8,'(PYCLUS:) too few particles for analysis')
        NJET=-1
        RETURN
      ENDIF
C...Find precluster configuration. If too few jets, make harder cuts.
      NLOOP=0
      IF(MSTU(46).LE.3.OR.MSTU(46).EQ.5) THEN
        R2ACC=PARU(44)**2
      ELSE
        R2ACC=PARU(45)*PS(5)**2
      ENDIF
      RINIT=1.25D0*PARU(43)
      IF(NP.LE.MSTU(47)+2) RINIT=0D0
  170 RINIT=0.8D0*RINIT
      NPRE=0
      NREM=NP
      DO 180 I=N+NP+1,N+2*NP
        K(I,4)=0
  180 CONTINUE
C...Sum up small momentum region. Jet if enough absolute momentum.
      IF(MSTU(46).LE.2) THEN
        DO 190 J=1,4
          P(N+1,J)=0D0
  190   CONTINUE
        DO 210 I=N+NP+1,N+2*NP
          IF(P(I,5).GT.2D0*RINIT) GOTO 210
          NREM=NREM-1
          K(I,4)=1
          DO 200 J=1,4
            P(N+1,J)=P(N+1,J)+P(I,J)
  200     CONTINUE
  210   CONTINUE
        P(N+1,5)=SQRT(P(N+1,1)**2+P(N+1,2)**2+P(N+1,3)**2)
        IF(P(N+1,5).GT.2D0*RINIT) NPRE=1
        IF(RINIT.GE.0.2D0*PARU(43).AND.NPRE+NREM.LT.MSTU(47)) GOTO 170
        IF(NREM.EQ.0) GOTO 170
      ENDIF
C...Find fastest remaining particle.
  220 NPRE=NPRE+1
      PMAX=0D0
      DO 230 I=N+NP+1,N+2*NP
        IF(K(I,4).NE.0.OR.P(I,5).LE.PMAX) GOTO 230
        IMAX=I
        PMAX=P(I,5)
  230 CONTINUE
      DO 240 J=1,5
        P(N+NPRE,J)=P(IMAX,J)
  240 CONTINUE
      NREM=NREM-1
      K(IMAX,4)=NPRE
C...Sum up precluster around it according to pT separation.
      IF(MSTU(46).LE.2) THEN
        DO 260 I=N+NP+1,N+2*NP
          IF(K(I,4).NE.0) GOTO 260
          R2=R2T(I,IMAX)
          IF(R2.GT.RINIT**2) GOTO 260
          NREM=NREM-1
          K(I,4)=NPRE
          DO 250 J=1,4
            P(N+NPRE,J)=P(N+NPRE,J)+P(I,J)
  250     CONTINUE
  260   CONTINUE
        P(N+NPRE,5)=SQRT(P(N+NPRE,1)**2+P(N+NPRE,2)**2+P(N+NPRE,3)**2)
C...Sum up precluster around it according to mass or
C...Durham pT separation.
      ELSE
  270   IMIN=0
        R2MIN=RINIT**2
        DO 280 I=N+NP+1,N+2*NP
          IF(K(I,4).NE.0) GOTO 280
          IF(MSTU(46).LE.4) THEN
            R2=R2M(I,N+NPRE)
          ELSE
            R2=R2D(I,N+NPRE)
          ENDIF
          IF(R2.GE.R2MIN) GOTO 280
          IMIN=I
          R2MIN=R2
  280   CONTINUE
        IF(IMIN.NE.0) THEN
          DO 290 J=1,4
            P(N+NPRE,J)=P(N+NPRE,J)+P(IMIN,J)
  290     CONTINUE
          P(N+NPRE,5)=SQRT(P(N+NPRE,1)**2+P(N+NPRE,2)**2+P(N+NPRE,3)**2)
          NREM=NREM-1
          K(IMIN,4)=NPRE
          GOTO 270
        ENDIF
      ENDIF
C...Check if more preclusters to be found. Start over if too few.
      IF(RINIT.GE.0.2D0*PARU(43).AND.NPRE+NREM.LT.MSTU(47)) GOTO 170
      IF(NREM.GT.0) GOTO 220
      NJET=NPRE
C...Reassign all particles to nearest jet. Sum up new jet momenta.
  300 TSAV=0D0
      PSJT=0D0
  310 IF(MSTU(46).LE.1) THEN
        DO 330 I=N+1,N+NJET
          DO 320 J=1,4
            V(I,J)=0D0
  320     CONTINUE
  330   CONTINUE
        DO 360 I=N+NP+1,N+2*NP
          R2MIN=PSS**2
          DO 340 IJET=N+1,N+NJET
            IF(P(IJET,5).LT.RINIT) GOTO 340
            R2=R2T(I,IJET)
            IF(R2.GE.R2MIN) GOTO 340
            IMIN=IJET
            R2MIN=R2
  340     CONTINUE
          K(I,4)=IMIN-N
          DO 350 J=1,4
            V(IMIN,J)=V(IMIN,J)+P(I,J)
  350     CONTINUE
  360   CONTINUE
        PSJT=0D0
        DO 380 I=N+1,N+NJET
          DO 370 J=1,4
            P(I,J)=V(I,J)
  370     CONTINUE
          P(I,5)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
          PSJT=PSJT+P(I,5)
  380   CONTINUE
      ENDIF
C...Find two closest jets.
      R2MIN=2D0*MAX(R2ACC,PS(5)**2)
      DO 400 ITRY1=N+1,N+NJET-1
        DO 390 ITRY2=ITRY1+1,N+NJET
          IF(MSTU(46).LE.2) THEN
            R2=R2T(ITRY1,ITRY2)
          ELSEIF(MSTU(46).LE.4) THEN
            R2=R2M(ITRY1,ITRY2)
          ELSE
            R2=R2D(ITRY1,ITRY2)
          ENDIF
          IF(R2.GE.R2MIN) GOTO 390
          IMIN1=ITRY1
          IMIN2=ITRY2
          R2MIN=R2
  390   CONTINUE
  400 CONTINUE
C...If allowed, join two closest jets and start over.
      IF(NJET.GT.MSTU(47).AND.R2MIN.LT.R2ACC) THEN
        IREC=MIN(IMIN1,IMIN2)
        IDEL=MAX(IMIN1,IMIN2)
        DO 410 J=1,4
          P(IREC,J)=P(IMIN1,J)+P(IMIN2,J)
  410   CONTINUE
        P(IREC,5)=SQRT(P(IREC,1)**2+P(IREC,2)**2+P(IREC,3)**2)
        DO 430 I=IDEL+1,N+NJET
          DO 420 J=1,5
            P(I-1,J)=P(I,J)
  420     CONTINUE
  430   CONTINUE
        IF(MSTU(46).GE.2) THEN
          DO 440 I=N+NP+1,N+2*NP
            IORI=N+K(I,4)
            IF(IORI.EQ.IDEL) K(I,4)=IREC-N
            IF(IORI.GT.IDEL) K(I,4)=K(I,4)-1
  440     CONTINUE
        ENDIF
        NJET=NJET-1
        GOTO 300
C...Divide up broad jet if empty cluster in list of final ones.
      ELSEIF(NJET.EQ.MSTU(47).AND.MSTU(46).LE.1.AND.NLOOP.LE.2) THEN
        DO 450 I=N+1,N+NJET
          K(I,5)=0
  450   CONTINUE
        DO 460 I=N+NP+1,N+2*NP
          K(N+K(I,4),5)=K(N+K(I,4),5)+1
  460   CONTINUE
        IEMP=0
        DO 470 I=N+1,N+NJET
          IF(K(I,5).EQ.0) IEMP=I
  470   CONTINUE
        IF(IEMP.NE.0) THEN
          NLOOP=NLOOP+1
          ISPL=0
          R2MAX=0D0
          DO 480 I=N+NP+1,N+2*NP
            IF(K(N+K(I,4),5).LE.1.OR.P(I,5).LT.RINIT) GOTO 480
            IJET=N+K(I,4)
            R2=R2T(I,IJET)
            IF(R2.LE.R2MAX) GOTO 480
            ISPL=I
            R2MAX=R2
  480     CONTINUE
          IF(ISPL.NE.0) THEN
            IJET=N+K(ISPL,4)
            DO 490 J=1,4
              P(IEMP,J)=P(ISPL,J)
              P(IJET,J)=P(IJET,J)-P(ISPL,J)
  490       CONTINUE
            P(IEMP,5)=P(ISPL,5)
            P(IJET,5)=SQRT(P(IJET,1)**2+P(IJET,2)**2+P(IJET,3)**2)
            IF(NLOOP.LE.2) GOTO 300
          ENDIF
        ENDIF
      ENDIF
C...If generalized thrust has not yet converged, continue iteration.
      IF(MSTU(46).LE.1.AND.NLOOP.LE.2.AND.PSJT/PSS.GT.TSAV+PARU(48))
     &THEN
        TSAV=PSJT/PSS
        GOTO 310
      ENDIF
C...Reorder jets according to energy.
      DO 510 I=N+1,N+NJET
        DO 500 J=1,5
          V(I,J)=P(I,J)
  500   CONTINUE
  510 CONTINUE
      DO 540 INEW=N+1,N+NJET
        PEMAX=0D0
        DO 520 ITRY=N+1,N+NJET
          IF(V(ITRY,4).LE.PEMAX) GOTO 520
          IMAX=ITRY
          PEMAX=V(ITRY,4)
  520   CONTINUE
        K(INEW,1)=31
        K(INEW,2)=97
        K(INEW,3)=INEW-N
        K(INEW,4)=0
        DO 530 J=1,5
          P(INEW,J)=V(IMAX,J)
  530   CONTINUE
        V(IMAX,4)=-1D0
        K(IMAX,5)=INEW
  540 CONTINUE
C...Clean up particle-jet assignments and jet information.
      DO 550 I=N+NP+1,N+2*NP
        IORI=K(N+K(I,4),5)
        K(I,4)=IORI-N
        IF(K(K(I,3),1).NE.3) K(K(I,3),4)=IORI-N
        K(IORI,4)=K(IORI,4)+1
  550 CONTINUE
      IEMP=0
      PSJT=0D0
      DO 570 I=N+1,N+NJET
        K(I,5)=0
        PSJT=PSJT+P(I,5)
        P(I,5)=SQRT(MAX(P(I,4)**2-P(I,5)**2,0D0))
        DO 560 J=1,5
          V(I,J)=0D0
  560   CONTINUE
        IF(K(I,4).EQ.0) IEMP=I
  570 CONTINUE
C...Select storing option. Output variables. Check for failure.
      MSTU(61)=N+1
      MSTU(62)=NP
      MSTU(63)=NPRE
      PARU(61)=PS(5)
      PARU(62)=PSJT/PSS
      PARU(63)=SQRT(R2MIN)
      IF(NJET.LE.1) PARU(63)=0D0
      IF(IEMP.NE.0) THEN
        CALL PYERRM(8,'(PYCLUS:) failed to reconstruct as requested')
        NJET=-1
        RETURN
      ENDIF
      IF(MSTU(43).LE.1) MSTU(3)=MAX(0,NJET)
      IF(MSTU(43).GE.2) N=N+MAX(0,NJET)
      NSAV=NJET
      RETURN
      END
