SubRoutine CluHad(Nrjet, IErr)
      Implicit None
C
C ************************************************************************
C
C PURPOSE:       To subdivide the hit calorimeter cells into jets/clusters.
C                CDF-like algorithm.
C CALLED FROM:   CccOfl
C COMMUNICATION: Common blocks CluCom, CDFCom and CluPar
C AUTHOR:        Paul de Jong, NIKHEF-H, 08-01-92
C ************************************************************************
C --DATE--:--NAME--:--MODIFICATIONS---------------------------------------
C
C ************************************************************************
C
#include "zrconst.inc"
#include "clucom.inc"
#include "clupar.inc"
#include "cdfcom.inc"
#include "zrunit.inc"
#include "cluout.inc"
C
      Integer K,L,M,Ncel,Iteration
      Real Deta,Deta1,Deta2,Dphi,Dphi1,Dphi2,Dist1,Dist2,SumE,SumEt,
     ,     SumEta,SumPhi,Phimax,Phimin,TPhi,Eoverlap
      Integer IErr,Nrjet,I,J,I1,I2,ITRY1,ITRY2,IREC,IDEL,IORI,IEMP
      REAL PS(5),PSS,RINIT,R2T,R2M,R2J,R2,R2MAX,R2MIN,PEMAX,
     ,     PSJT,TSAV,R2ACC,PMAX,EVISIB
      Integer NSAV,NP,NPRE,NREM,NJET,IJET,ISPL,IMIN1,IMIN2,NLOOP
      Integer NEDIT,IMAX,IMIN,INEW,ITRY
      Real Type,PiMass,Charge,Life,Ub(10)
      Integer Nwb, Name(5), Idum, Icjet, Icjet1, Icjet2
      Logical First, Ready
      Data First /.TRUE./
C
      Call ModIn ('CluHad',Idum)
      IErr=0
C
C Loop over the seed cells and arrange them in preclusters.
C
      Ijet=0
      Do 100 I=1,Nseed
         If (KC(I,2).NE.2) GoTo 100
         Ijet=Ijet+1
         KC(I,1)=Ijet
         KC(I,2)=0
         Do 110 J=I+1,Nseed
            If (KC(J,2).NE.2) GoTo 110
            Deta=Abs(PC(I,1)-PC(J,1))
            Dphi=Abs(PC(I,2)-PC(J,2))
            If (Dphi.GT.Pi) Dphi=Twopi-Dphi
            If (Sqrt(Deta**2+Dphi**2).LE.RadPC) Then
               KC(J,1)=KC(I,1)
               KC(J,2)=0
            Endif
  110    Continue
  100 Continue
C
C Check for safety margin
C
      If ((Npart+2*Ijet).GT.Nsblock) Then
         Write(Lp,1001)
 1001    Format(' CluHad: Overflow in common block!')
         Call MoSetC('Overflow in common block!',1.,1)
         IErr=5
         GoTo 999
      Endif
C
C Calculate the centres of the preclusters from the contributing cells
C
      Do 120 I=1,Ijet
         SumE=0.
         SumEt=0.
         Sumeta=0.
         Sumphi=0.
         Ncel=0
         Phimax=0.
         Phimin=Twopi
         Do 130 J=1,Nseed
            If (KC(J,1).EQ.I) Then
               SumE=SumE+PC(J,4)
               SumEt=SumEt+PC(J,3)
               If (Imode.EQ.1) Then
                  Sumeta=Sumeta+PC(J,1)*PC(J,3)
                  Sumphi=Sumphi+PC(J,2)*PC(J,3)
               Else if (Imode.EQ.2) Then
                  Sumeta=Sumeta+PC(J,1)*PC(J,4)
                  Sumphi=Sumphi+PC(J,2)*PC(J,4)
               Endif
               If (PC(J,2).GT.Phimax) Phimax=PC(J,2)
               If (PC(J,2).LT.Phimin) Phimin=PC(J,2)
               Ncel=Ncel+1
            Endif
  130    Continue
         If (SumEt.LE.0.) Then
            Write(Lp,1003)
 1003       Format(' CluHad: Et = 0.?')
            SumEt=1.
            SumE=1.
         Endif
         If (Imode.EQ.1) Then
            PC(Npart+I,1)=Sumeta/SumEt
            PC(Npart+I,2)=Sumphi/SumEt
         Else if (Imode.EQ.2) Then
            PC(Npart+I,1)=Sumeta/SumE
            PC(Npart+I,2)=Sumphi/SumE
         Endif
         If ((Phimax-Phimin).GT.(3*Pi/2.)) Then
            Sumphi=0.
            Do 135 J=1,Nseed
               If (KC(J,1).EQ.I) Then
                  If (PC(J,2).GT.Pi) Then
                     TPhi=PC(J,2)-Twopi
                  Else
                     TPhi=PC(J,2)
                  Endif
                  If (Imode.EQ.1) Then
                     Sumphi=Sumphi+TPhi*PC(J,3)
                  Else if (Imode.EQ.2) Then
                     Sumphi=Sumphi+TPhi*PC(J,4)
                  Endif
               Endif
  135       Continue
            If (Imode.EQ.1) Then
               PC(Npart+I,2)=Sumphi/SumEt
            Else if (Imode.EQ.2) Then
               PC(Npart+I,2)=Sumphi/SumE
            Endif
            If (PC(Npart+I,2).LT.0.) PC(Npart+I,2)=
     +            PC(Npart+I,2)+Twopi
         Endif
         PC(Npart+I,3)=SumEt
         PC(Npart+I,4)=SumE
         PC(Npart+I,5)=REAL(I)
         KC(Npart+I,1)=Ncel
         KC(Npart+I,2)=3
  120 Continue
C
C We now have the preclusters with their centres. For each cluster we now
C assign cells to it if they are within the required distance.
C
      Iteration=0
    5 Ready=.TRUE.
      Iteration=Iteration+1
C
C .. reset the counter
C
      Do 10 I=1,Nsblock
         PC(I,5)=0.
   10 Continue
      Do 140 I=1,Ijet
         Do 150 J=1,Npart
            Deta=Abs(PC(Npart+I,1)-PC(J,1))
            Dphi=Abs(PC(Npart+I,2)-PC(J,2))
            If (Dphi.GT.Pi) Dphi=Twopi-Dphi
            If (Sqrt(Deta**2+Dphi**2).LE.Radius) Then
               PC(J,5)=PC(J,5)+1.
               If (PC(J,5).GE.9.) Then
                  PC(J,5)=8.
                  Write(Lp,1004)
 1004             Format(' CluHad: Overflow in counter!')
                  Call MoSetC('Overflow in counter!',1.,0)
                  IErr=6
               Endif
               If (PC(J,5).GE.6.) Then
                  If (INT(VC(J,INT(PC(J,5)-5.))).NE.I) Ready=.FALSE.
                  VC(J,INT(PC(J,5)-5.))=REAL(I)
               Else
                  If (KC(J,INT(PC(J,5))).NE.I) Ready=.FALSE.
                  KC(J,INT(PC(J,5)))=I
               Endif
            Endif
  150    Continue
  140 Continue
C
C .. also check the counter for convergence
C
      Do 155 J=1,Npart
         If (INT(PC(J,5)).NE.INT(VC(J,4))) Ready=.FALSE.
         VC(J,4)=PC(J,5)
  155 Continue
C
C .. reevaluate the jet centres based on the new list
C
      Do 160 I=1,Ijet
         SumE=0.
         SumEt=0.
         Sumeta=0.
         Sumphi=0.
         Ncel=0
         Phimax=0.
         Phimin=Twopi
         Do 170 J=1,Npart
            Do 180 K=1,INT(PC(J,5))
               If (K.GE.6) Then
                  Icjet = INT(VC(J,K-5))
               Else
                  Icjet = KC(J,K)
               Endif
               If (Icjet.EQ.I) Then
                  SumE=SumE+PC(J,4)
                  SumEt=SumEt+PC(J,3)
                  If (Imode.EQ.1) Then
                     Sumeta=Sumeta+PC(J,1)*PC(J,3)
                     Sumphi=Sumphi+PC(J,2)*PC(J,3)
                  Else if (Imode.EQ.2) Then
                     Sumeta=Sumeta+PC(J,1)*PC(J,4)
                     Sumphi=Sumphi+PC(J,2)*PC(J,4)
                  Endif
                  If (PC(J,2).GT.Phimax) Phimax=PC(J,2)
                  If (PC(J,2).LT.Phimin) Phimin=PC(J,2)
                  Ncel=Ncel+1
               Endif
  180       Continue
  170    Continue
         If (SumEt.LE.0.) Then
            Write(Lp,1003)
            SumEt=1.
            SumE=1.
         Endif
         If (Imode.EQ.1) Then
            PC(Npart+I,1)=Sumeta/SumEt
            PC(Npart+I,2)=Sumphi/SumEt
         Else if (Imode.EQ.2) Then
            PC(Npart+I,1)=Sumeta/SumE
            PC(Npart+I,2)=Sumphi/SumE
         Endif
         If ((Phimax-Phimin).GT.(3*Pi/2.)) Then
            Sumphi=0.
            Do 175 J=1,Npart
               Do 185 K=1,INT(PC(J,5))
                  If (K.GE.6) Then
                     Icjet = INT(VC(J,K-5))
                  Else
                     Icjet = KC(J,K)
                  Endif
                  If (Icjet.EQ.I) Then
                     If (PC(J,2).GT.Pi) Then
                        TPhi=PC(J,2)-Twopi
                     Else
                        TPhi=PC(J,2)
                     Endif
                     If (Imode.EQ.1) Then
                        Sumphi=Sumphi+TPhi*PC(J,3)
                     Else if (Imode.EQ.2) Then
                        Sumphi=Sumphi+TPhi*PC(J,4)
                     Endif
                  Endif
  185          Continue
  175       Continue
            If (Imode.EQ.1) Then
               PC(Npart+I,2)=Sumphi/SumEt
            Else if (Imode.EQ.2) Then
               PC(Npart+I,2)=Sumphi/SumE
            Endif
            If (PC(Npart+I,2).LT.0.) PC(Npart+I,2)=
     +            PC(Npart+I,2)+Twopi
         Endif
         PC(Npart+I,3)=SumEt
         PC(Npart+I,4)=SumE
         PC(Npart+I,5)=REAL(I)
         KC(Npart+I,1)=Ncel
         KC(Npart+I,2)=3
  160 Continue
C
C Repeat the iteration if list has changed, but cut off at MAXITER iterations
C
      If ((.NOT.Ready).AND.Iteration.LT.MaxIter) GoTo 5
C
C Resolve overlaps. If the amount of overlapping energy is more than some
C fraction (typically 75%) of the energy of the smallest cluster, the
C clusters are merged, otherwise the cells are assigned to the closest cluster
C
C...Reorder jets according to Et
C
      DO 460 I=Npart+1,Npart+IJET
      DO 460 J=1,5
         VC(I,J)=PC(I,J)
         VC(I+IJET,J)=REAL(KC(I,J))
  460 CONTINUE
      DO 490 INEW=Npart+1,Npart+IJET
         PEMAX=0.
         DO 470 ITRY=Npart+1,Npart+IJET
            IF(VC(ITRY,3).LE.PEMAX) GOTO 470
            IMAX=ITRY
            PEMAX=VC(ITRY,3)
  470    CONTINUE
         DO 480 J=1,5
            KC(INEW,J)=INT(VC(IMAX+IJET,J))
            PC(INEW,J)=VC(IMAX,J)
  480    CONTINUE
         VC(IMAX,3)=-1.
  490 CONTINUE
C
      Do 190 I=1,Ijet
         Do 200 J=Ijet,I+1,-1
            Eoverlap=0.
            Do 210 K=1,Npart
               If (PC(K,5).LE.1.) GoTo 210
               Do 220 L=1,INT(PC(K,5))
               Do 220 M=1,INT(PC(K,5))
                  If (L.GE.6) Then
                     Icjet1 = INT(VC(K,L-5))
                  Else
                     Icjet1 = KC(K,L)
                  Endif
                  If (M.GE.6) Then
                     Icjet2 = INT(VC(K,M-5))
                  Else
                     Icjet2 = KC(K,M)
                  Endif
                  If ((Icjet1.EQ.INT(PC(Npart+I,5))).AND.
     .                (Icjet2.EQ.INT(PC(Npart+J,5)))) Then
                     Eoverlap=Eoverlap+PC(K,4)
                     GoTo 210
                  Endif
  220          Continue
  210       Continue
            If (Eoverlap.LE.0.) GoTo 200
            If (Eoverlap.GE.Frac*MIN(PC(Npart+I,4),
     ,            PC(Npart+J,4))) Then
C
C .. merge !
C
               Do 230 K=1,Npart
                  Do 240 L=1,INT(PC(K,5))
                     If (L.GE.6) Then
                        If (INT(VC(K,L-5)).EQ.INT(PC(Npart+J,5))) Then
                           VC(K,L-5)=PC(Npart+I,5)
                        Endif
                     Else
                        If (KC(K,L).EQ.INT(PC(Npart+J,5))) Then
                           KC(K,L)=INT(PC(Npart+I,5))
                        Endif
                     Endif
  240             Continue
  230          Continue
            Else
C
C .. assign each cell to the closest jet
C
               Do 250 K=1,Npart
                  Deta1=Abs(PC(Npart+I,1)-PC(K,1))
                  Dphi1=Abs(PC(Npart+I,2)-PC(K,2))
                  If (Dphi1.GT.Pi) Dphi1=Twopi-Dphi1
                  Dist1=Sqrt(Deta1**2+Dphi1**2)
                  Deta2=Abs(PC(Npart+J,1)-PC(K,1))
                  Dphi2=Abs(PC(Npart+J,2)-PC(K,2))
                  If (Dphi2.GT.Pi) Dphi2=Twopi-Dphi2
                  Dist2=Sqrt(Deta2**2+Dphi2**2)
                  If (Dist1.LE.Dist2) Then
                     Do 260 L=1,INT(PC(K,5))
                        If (L.GE.6) Then
                         If (INT(VC(K,L-5)).EQ.INT(PC(Npart+J,5))) Then
                            VC(K,L-5)=PC(Npart+I,5)
                         Endif
                        Else
                           If (KC(K,L).EQ.INT(PC(Npart+J,5))) Then
                              KC(K,L)=INT(PC(Npart+I,5))
                           Endif
                        Endif
  260                Continue
                  Else
                     Do 270 L=1,INT(PC(K,5))
                        If (L.GE.6) Then
                         If (INT(VC(K,L-5)).EQ.INT(PC(Npart+I,5))) Then
                            VC(K,L-5)=PC(Npart+J,5)
                         Endif
                        Else
                           If (KC(K,L).EQ.INT(PC(Npart+I,5))) Then
                              KC(K,L)=INT(PC(Npart+J,5))
                           Endif
                        Endif
  270                Continue
                  Endif
  250          Continue
            Endif
  200    Continue
  190 Continue
C
C Recalculate the centres from the new lists, and clean up.
C
      Nrjet=0
      Do 280 I=1,Ijet
         SumE=0.
         SumEt=0.
         Sumeta=0.
         Sumphi=0.
         Ncel=0
         Phimax=0.
         Phimin=Twopi
         Do 290 J=1,Npart
            If (INT(PC(J,5)).LE.0) GoTo 290
            If (KC(J,1).EQ.INT(PC(Npart+I,5))) Then
               SumE=SumE+PC(J,4)
               SumEt=SumEt+PC(J,3)
               If (Imode.EQ.1) Then
                  Sumeta=Sumeta+PC(J,1)*PC(J,3)
                  Sumphi=Sumphi+PC(J,2)*PC(J,3)
               Else if (Imode.EQ.2) Then
                  Sumeta=Sumeta+PC(J,1)*PC(J,4)
                  Sumphi=Sumphi+PC(J,2)*PC(J,4)
               Endif
               If (PC(J,2).GT.Phimax) Phimax=PC(J,2)
               If (PC(J,2).LT.Phimin) Phimin=PC(J,2)
               Ncel=Ncel+1
            Endif
  290    Continue
         If (Ncel.LE.0) Then
            Do 300 J=1,5
               PC(Npart+I,J)=0.
               KC(Npart+I,J)=0
  300       Continue
         Else
            Nrjet=Nrjet+1
            If (SumEt.LE.0.) Then
               Write(Lp,1003)
               SumEt=1.
               SumE=1.
            Endif
            If (Imode.EQ.1) Then
               PC(Npart+I,1)=Sumeta/SumEt
               PC(Npart+I,2)=Sumphi/SumEt
            Else if (Imode.EQ.2) Then
               PC(Npart+I,1)=Sumeta/SumE
               PC(Npart+I,2)=Sumphi/SumE
            Endif
            If ((Phimax-Phimin).GT.(3*Pi/2.)) Then
               Sumphi=0.
               Do 305 J=1,Npart
                  If (KC(J,1).EQ.INT(PC(Npart+I,5))) Then
                     If (PC(J,2).GT.Pi) Then
                        TPhi=PC(J,2)-Twopi
                     Else
                        TPhi=PC(J,2)
                     Endif
                     If (Imode.EQ.1) Then
                        Sumphi=Sumphi+TPhi*PC(J,3)
                     Else if (Imode.EQ.2) Then
                        Sumphi=Sumphi+TPhi*PC(J,4)
                     Endif
                  Endif
  305          Continue
               If (Imode.EQ.1) Then
                  PC(Npart+I,2)=Sumphi/SumEt
               Else if (Imode.EQ.2) Then
                  PC(Npart+I,2)=Sumphi/SumE
               Endif
               If (PC(Npart+I,2).LT.0.) PC(Npart+I,2)=
     +               PC(Npart+I,2)+Twopi
            Endif
            PC(Npart+I,3)=SumEt
            PC(Npart+I,4)=SumE
            KC(Npart+I,1)=Ncel
            KC(Npart+I,2)=3
         Endif
  280 Continue
C
C...Finally reorder jets according to Et
C
      DO 310 I=Npart+1,Npart+IJET
      DO 310 J=1,5
         VC(I,J)=PC(I,J)
         VC(I+IJET,J)=REAL(KC(I,J))
  310 CONTINUE
      DO 320 INEW=Npart+1,Npart+IJET
         PEMAX=-0.001
         DO 330 ITRY=Npart+1,Npart+IJET
            IF(VC(ITRY,3).LE.PEMAX) GOTO 330
            IMAX=ITRY
            PEMAX=VC(ITRY,3)
  330    CONTINUE
         DO 340 J=1,5
            KC(INEW,J)=INT(VC(IMAX+IJET,J))
            PC(INEW,J)=VC(IMAX,J)
  340    CONTINUE
         VC(IMAX,3)=-1.
  320 CONTINUE

      Npart=Npart+Nrjet
C
  999 Call ModOut ('CluHad')
      RETURN
      END