subroutine hec_beam_fit(mode)
***********************************************************************
***    Reconstruct beam trajectory parameters and put them in common/beam/
***    At most 5 MWPC can be used
***    if mode=1, reconstruct beam particle trajectory
***    if mode=2, print test/tune histograms ID=2051-2062
***
***    Created  20-FEB-1995  A.MINAENKO
***    Modified 25-JUN-1997  A.Minaenko
***
***********************************************************************
      implicit none

      include 'hec_par.inc'              ! shared parameters
      INCLUDE 'hec_beam.inc'
      INCLUDE 'hec_beam_sys_cluset.inc'
      INCLUDE 'hec_beam_sys_clusters.inc'
      INCLUDE 'hec_beam_sys_config.inc'
      INCLUDE 'hec_datacard.inc'
      include 'hec_evth.inc'
      INCLUDE 'hec_trig.inc'
c
C===    Local variables
      CHARACTER*80 CHTIT
      CHARACTER*1 CHN(5)
      SAVE CHN
C===
      LOGICAL FIRST
      SAVE FIRST
C
      real resmin
      integer mode
      
      INTEGER IHZ
      SAVE IHZ
      INTEGER I,I1,I2,I3,I4,I5,NCX,NCY,NCMAX,NCMAX_S
C
      REAL RESMAXMAX    	!  acceptable biggest residual 
      SAVE RESMAXMAX            !  in 0.5*(cluster size) units
      REAL ZX2(5),ZY2(5)
      SAVE ZX2,ZY2
      REAL SLOP,XZERO,RESMEAN,RMAX
      REAL SLOP_S,XZERO_S,RMAX_S
      REAL IVCLU_S(5),CCLU_S(5),WCLU_S(5)
c
c data statement initializations
c
      data first/.true./
      data ihz/5050/
      data resmaxmax/1.5/
      data chn/'1','2','3','4','5'/
C
C ----------------------------------------------------------------------
C
c
c check if package is requested
c
      IF (z_beam .gt. 100000.) RETURN 
c
      IF (mode.NE.1.AND.mode.NE.2) RETURN	! wrong mode
C===    Initialize arrays
      IF (FIRST) THEN
       FIRST=.FALSE.
       DO I=1,5
          ZX2(I)=ZXPLANE(I)**2
          ZY2(I)=ZYPLANE(I)**2
       ENDDO
c
c        Book histograms
c
       IF (hpr_beam.ne.0) THEN
        CHTIT='MEAN RESIDUALS X-PLANES AFTER FIRST LINEFIT CALL$'
        CALL HBOOK1(IHZ+1,CHTIT,100,0.,10.,0.)
        CHTIT='MEAN RESIDUALS Y-PLANES AFTER FIRST LINEFIT CALL$'
        CALL HBOOK1(IHZ+2,CHTIT,100,0.,10.,0.)
        CHTIT='MAX RESIDUALS X-PLANES AFTER FIRST LINEFIT CALL$'
        CALL HBOOK1(IHZ+3,CHTIT,100,0.,10.,0.)
        CHTIT='MAX RESIDUALS Y-PLANES AFTER FIRST LINEFIT CALL$'
        CALL HBOOK1(IHZ+4,CHTIT,100,0.,10.,0.)
        CHTIT='MEAN RESIDUALS X-PLANES AFTER LAST LINEFIT CALL$'
        CALL HBOOK1(IHZ+6,CHTIT,100,0.,10.,0.)
        CHTIT='MEAN RESIDUALS Y-PLANES AFTER LAST LINEFIT CALL$'
        CALL HBOOK1(IHZ+7,CHTIT,100,0.,10.,0.)
        CHTIT='MAX RESIDUALS X-PLANES AFTER LAST LINEFIT CALL$'
        CALL HBOOK1(IHZ+8,CHTIT,100,0.,10.,0.)
        CHTIT='MAX RESIDUALS Y-PLANES AFTER LAST LINEFIT CALL$'
        CALL HBOOK1(IHZ+9,CHTIT,100,0.,10.,0.)
        CHTIT='REJECTED X-PLANES LIST$'
        CALL HBOOK1(IHZ+11,CHTIT,5,1.,6.,0.)
        CHTIT='REJECTED Y-PLANES LIST$'
        CALL HBOOK1(IHZ+12,CHTIT,5,1.,6.,0.)
       ENDIF
c
      ENDIF
C
C===    Print histograms
C
      IF (mode.EQ.2) THEN
        IF (hpr_beam.ne.0) THEN
          DO I=1,4
            CALL HPRINT(IHZ+I)
            CALL HPRINT(IHZ+5+I)
            CALL HDELET(IHZ+I)
            CALL HDELET(IHZ+5+I)
          ENDDO
          CALL HPRINT(IHZ+11)
          CALL HPRINT(IHZ+12)
          CALL HDELET(IHZ+11)
          CALL HDELET(IHZ+12)
        ENDIF
        RETURN
      ENDIF
C
      IF (.NOT.trig_physics) RETURN	! not a physical event
      CALL hec_BEAM_GETCLUST(1)		! take clusters from beam region
      NCX=0
      NCY=0
      xbeam_ok=.false.
      ybeam_ok=.false.
      DO I=1,5
         IVCLU(I)=0
         IF (NCLUX(I).GT.0) NCX=NCX+1
         IF (NCLUY(I).GT.0) NCY=NCY+1
      ENDDO
C
C===    X-planes fit
C
      IF (NCX.GE.3) THEN
       RESMIN=1.E+38
       CALL UCOPY(ZXPLANE,ZCLU,5)
       CALL UCOPY(ZX2,ZCLU2,5)
C===    Track match loop
       DO I1=1, MAX0(NCLUX(1),1)
        IF (NCLUX(1).GT.0) THEN
         IVCLU(1)=I1
         CCLU(1)=CLUX(1,I1,1)
         WCLU(1)=CLUX(2,I1,1)
        ENDIF
        DO I2=1, MAX0(NCLUX(2),1)
         IF (NCLUX(2).GT.0) THEN
          IVCLU(2)=I2
          CCLU(2)=CLUX(1,I2,2)
          WCLU(2)=CLUX(2,I2,2)
         ENDIF
         DO I3=1, MAX0(NCLUX(3),1)
          IF (NCLUX(3).GT.0) THEN
           IVCLU(3)=I3
           CCLU(3)=CLUX(1,I3,3)
           WCLU(3)=CLUX(2,I3,3)
          ENDIF
          DO I4=1, MAX0(NCLUX(4),1)
           IF (NCLUX(4).GT.0) THEN
            IVCLU(4)=I4
            CCLU(4)=CLUX(1,I4,4)
            WCLU(4)=CLUX(2,I4,4)
           ENDIF
           DO I5=1, MAX0(NCLUX(5),1)
            IF (NCLUX(5).GT.0) THEN
             IVCLU(5)=I5
             CCLU(5)=CLUX(1,I5,5)
             WCLU(5)=CLUX(2,I5,5)
            ENDIF
            CALL hec_BEAM_LINEFIT(SLOP,XZERO,RESMEAN,RMAX,NCMAX)
            IF (RESMEAN.LT.RESMIN) THEN
             RESMIN=RESMEAN
             SLOP_S=SLOP
             XZERO_S=XZERO
             RMAX_S=RMAX
             NCMAX_S=NCMAX
             DO I=1,5
              IVCLU_S(I)=IVCLU(I)
              CCLU_S(I)=CCLU(I)
              WCLU_S(I)=WCLU(I)
             ENDDO
            ENDIF
           ENDDO
          ENDDO
         ENDDO				
        ENDDO
       ENDDO
C===     Take the best combination of clusters
       RESMEAN=RESMIN
       SLOP=SLOP_S
       XZERO=XZERO_S
       RMAX=RMAX_S
       NCMAX=NCMAX_S
       DO I=1,5
          IVCLU(I)=IVCLU_S(I)
          CCLU(I)=CCLU_S(I)
          WCLU(I)=WCLU_S(I)
       ENDDO
C===
       IF (hpr_beam.ne.0) THEN
        CALL HF1(IHZ+1,RESMEAN,1.)
        CALL HF1(IHZ+3,RMAX,1.)
       ENDIF
C===    Check the fit quality and reject plane, if necessary
       IF (RMAX.LE.RESMAXMAX) THEN
        xbeam_ok=.TRUE.
       ELSE
        NCX=NCX-1
        CALL HF1(IHZ+11,FLOAT(NCMAX),1.)
        IF (NCX.GE.3) THEN
         IVCLU(NCMAX)=0
         CALL hec_BEAM_LINEFIT(SLOP,XZERO,RESMEAN,RMAX,NCMAX)
         IF (RMAX.LE.RESMAXMAX) THEN
          xbeam_ok=.TRUE.
         ELSE
          NCX=NCX-1
          CALL HF1(IHZ+11,FLOAT(NCMAX),1.)
          IF (NCX.GE.3) THEN
           IVCLU(NCMAX)=0
           CALL hec_BEAM_LINEFIT(SLOP,XZERO,RESMEAN,RMAX,NCMAX)
           IF (RMAX.LE.RESMAXMAX) THEN
            xbeam_ok=.TRUE.
           ELSE
            CALL HF1(IHZ+11,FLOAT(NCMAX),1.)
           ENDIF
          ENDIF
         ENDIF				! ELSE RMAX.LE.RESMAXMAX
        ENDIF				! NCX.GE.3
       ENDIF				! ELSE RMAX.LE.RESMAXMAX
C
       xbeam = SLOP*zbeam+XZERO
       xslop = SLOP
C
       IF (hpr_beam.ne.0.AND.xbeam_ok) THEN
          CALL HF1(IHZ+6,RESMEAN,1.)
          CALL HF1(IHZ+8,RMAX,1.)
       ENDIF
      ENDIF				! NCX.GE.3
C
C===    Y-planes fit
C
      IF (NCY.GE.3) THEN
       RESMIN=1.E+38
       CALL UCOPY(ZYPLANE,ZCLU,5)
       CALL UCOPY(ZY2,ZCLU2,5)
       DO I=1,5
          IVCLU(I)=0
       ENDDO
C===    Track match loop
       DO I1=1, MAX0(NCLUY(1),1)
        IF (NCLUY(1).GT.0) THEN
         IVCLU(1)=I1
         CCLU(1)=CLUY(1,I1,1)
         WCLU(1)=CLUY(2,I1,1)
        ENDIF
        DO I2=1, MAX0(NCLUY(2),1)
         IF (NCLUY(2).GT.0) THEN
          IVCLU(2)=I2
          CCLU(2)=CLUY(1,I2,2)
          WCLU(2)=CLUY(2,I2,2)
         ENDIF
         DO I3=1, MAX0(NCLUY(3),1)
          IF (NCLUY(3).GT.0) THEN
           IVCLU(3)=I3
           CCLU(3)=CLUY(1,I3,3)
           WCLU(3)=CLUY(2,I3,3)
          ENDIF
          DO I4=1, MAX0(NCLUY(4),1)
           IF (NCLUY(4).GT.0) THEN
            IVCLU(4)=I4
            CCLU(4)=CLUY(1,I4,4)
            WCLU(4)=CLUY(2,I4,4)
           ENDIF
           DO I5=1, MAX0(NCLUY(5),1)
            IF (NCLUY(5).GT.0) THEN
             IVCLU(5)=I5
             CCLU(5)=CLUY(1,I5,5)
             WCLU(5)=CLUY(2,I5,5)
            ENDIF
            CALL hec_BEAM_LINEFIT(SLOP,XZERO,RESMEAN,RMAX,NCMAX)
            IF (RESMEAN.LT.RESMIN) THEN
             RESMIN=RESMEAN
             SLOP_S=SLOP
             XZERO_S=XZERO
             RMAX_S=RMAX
             NCMAX_S=NCMAX
             DO I=1,5
                IVCLU_S(I)=IVCLU(I)
                CCLU_S(I)=CCLU(I)
                WCLU_S(I)=WCLU(I)
             ENDDO
            ENDIF
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDDO
       RESMEAN=RESMIN
       SLOP=SLOP_S
       XZERO=XZERO_S
       RMAX=RMAX_S
       NCMAX=NCMAX_S
       DO I=1,5
          IVCLU(I)=IVCLU_S(I)
          CCLU(I)=CCLU_S(I)
          WCLU(I)=WCLU_S(I)
       ENDDO
C===
       IF (hpr_beam.ne.0) THEN
          CALL HF1(IHZ+2,RESMEAN,1.)
          CALL HF1(IHZ+4,RMAX,1.)
       ENDIF
C===    Check the fit quality and reject plane, if necessary
       IF (RMAX.LE.RESMAXMAX) THEN
        ybeam_ok=.TRUE.
       ELSE
        NCY=NCY-1
        CALL HF1(IHZ+12,FLOAT(NCMAX),1.)
        IF (NCY.GE.3) THEN
         IVCLU(NCMAX)=0
         CALL hec_BEAM_LINEFIT(SLOP,XZERO,RESMEAN,RMAX,NCMAX)
         IF (RMAX.LE.RESMAXMAX) THEN
          ybeam_ok=.TRUE.
         ELSE
          CALL HF1(IHZ+12,FLOAT(NCMAX),1.)
          NCY=NCY-1
          IF (NCY.GE.3) THEN
           IVCLU(NCMAX)=0
           CALL hec_BEAM_LINEFIT(SLOP,XZERO,RESMEAN,RMAX,NCMAX)
           IF (RMAX.LE.RESMAXMAX) THEN
            ybeam_ok=.TRUE.
           ELSE
            CALL HF1(IHZ+12,FLOAT(NCMAX),1.)
           ENDIF
          ENDIF
         ENDIF				! ELSE RMAX.LE.RESMAXMAX
        ENDIF				! NCY.GE.3
       ENDIF				! ELSE RMAX.LE.RESMAXMAX
C
       ybeam = SLOP*zbeam+XZERO
       yslop = SLOP
C
       IF (hpr_beam.ne.0.AND.ybeam_ok) THEN
          CALL HF1(IHZ+7,RESMEAN,1.)
          CALL HF1(IHZ+9,RMAX,1.)
       ENDIF
      ENDIF				! NCY.GE.3

      RETURN
      END