subroutine hec_beam_survey(mode)
************************************************************************
***   if mode=1 store data for alignment, reading current event
***   if mode=2 calculate alignment constants
***
***   created  17-feb-1995  a.minaenko
***   updated  02-mar-1995  a.kiryunin
***   modified 24-jun-1997  a.minaenko
***
************************************************************************
      implicit none
c
      include 'hec_par.inc'
      include 'hec_beam_sys_clusters.inc'
      include 'hec_beam_sys_config.inc'
      include 'hec_beam_sys_shifts.inc'
      include 'hec_datacard.inc'
      include 'hec_runh.inc'
c
      integer mode
c===    local variables
      character*80 chtit
      character*1 chn(5)
      save chn
c
      integer numev, ihz
      save numev, ihz
      integer i, ipla, idum, nbin, nbmax, lvmax, nblow, nbup, iclu
c
      real y_tabl_ref
      save y_tabl_ref
      real dx_refrun, dy_refrun, dzx, dzy
      save dx_refrun, dy_refrun, dzx, dzy
      real hi, hstati, hsum
      real vect(200),x1,x2,bwidth,dum,x_slop,xzero,y_slop,yzero
      real nevx(5),nevy(5),cont,out,x,y,dx,dy,xpred,ypred
c
      data chn/'1','2','3','4','5'/
      data numev/0/
      data ihz/5000/
      data y_tabl_ref/0.0/
c
      if (mode.lt.1.or.mode.gt.2) return	! incorrect mode
      if (mode.eq.1) then
       if (numev.eq.0) then			! initialization
         if (nrfl.eq.5) y_tabl_ref=y_tabl
         dzx=zxplane(nrfl)-zxplane(nrff)
         dzy=zyplane(nrfl)-zyplane(nrff)
c===    save shifts of the last ref.mwpc from reference run
         dx_refrun=shiftx(nrfl)
         dy_refrun=shifty(nrfl)
         do i=1,5
            chtit='(xpred-xreal) x-plane chamber #'//chn(i)
            call hbook1(ihz+i,chtit,120,-6.,6.,0.)
            chtit='(ypred-yreal) y-plane chamber #'//chn(i)
            call hbook1(ihz+5+i,chtit,120,-6.,6.,0.)
         enddo
         if (hpr_beam.ne.0) then
          do i=21,25
            chtit='beam profile: x-plane chamber #'//chn(i-20)
            call hbook1(ihz+i,chtit,120,-6.,6.,0.)
            chtit='beam profile: y-plane chamber #'//chn(i-20)
            call hbook1(ihz+5+i,chtit,120,-6.,6.,0.)
          enddo
         endif
       endif					! numev=0
c
c===      first 40 events are used to find precise histogram boundaries
c===      redefine alignment histograms with precize boundaries
       if (numev.eq.40) then
        do i=1,5
           call hgive(ihz+i,chtit,nbin,x1,x2,idum,dum,dum,
     &                idum,idum)
           call hunpak(ihz+i,vect,' ',0)
           nbmax=lvmax(vect,nbin)
           bwidth=(x2-x1)/nbin
           x1=x1+bwidth*(nbmax-0.5)
           if (i.eq.nrff.or.i.eq.nrfl) then
              x2=x1+2.0
              x1=x1-2.0
              chtit=' beam profile x-plane chamber #'//chn(i)
           else
              x2=x1+2.5*step(i)
              x1=x1-2.5*step(i)
              chtit='(xpred-xreal) x-plane chamber #'//chn(i)
           endif
           call hdelet(ihz+i)
           call hbook1(ihz+i,chtit,50,x1,x2,0.)
           call hidopt(ihz+i,'stat')
c===
           call hgive(ihz+5+i,chtit,nbin,x1,x2,idum,dum,dum,
     &                idum,idum)
           call hunpak(ihz+5+i,vect,' ',0)
           nbmax=lvmax(vect,nbin)
           bwidth=(x2-x1)/nbin
           x1=x1+bwidth*(nbmax-0.5)
           if (i.eq.nrff.or.i.eq.nrfl) then
              x2=x1+2.0
              x1=x1-2.0
              chtit=' beam profile y-plane chamber #'//chn(i)
           else
              x2=x1+2.5*step(i)
              x1=x1-2.5*step(i)
              chtit='(ypred-yreal) y-plane chamber #'//chn(i)
           endif
           call hdelet(ihz+5+i)
           call hbook1(ihz+5+i,chtit,50,x1,x2,0.)
           call hidopt(ihz+5+i,'stat')
        enddo
       endif					! numev=40
c===    begin alignment histogram filling
       numev=numev+1
       call hec_beam_getclust(0)		! take all clusters
c===     x-plane histogram filling
       if (nclux(nrff).eq.1.and.nclux(nrfl).eq.1) then
        x_slop=(clux(1,1,nrfl)-clux(1,1,nrff))/dzx
        xzero=clux(1,1,nrff)-x_slop*zxplane(nrff)
        call hf1(ihz+nrff,clux(1,1,nrff),1.)
        call hf1(ihz+nrfl,clux(1,1,nrfl),1.)
        do ipla=1,5
        if (ipla.ne.nrff.and.ipla.ne.nrfl) then
         if (nclux(ipla).ge.1) then
          xpred=xzero+x_slop*zxplane(ipla)
          do iclu=1,nclux(ipla)
          if (clux(2,iclu,ipla).le.(0.6*step(ipla))) then
           dx=clux(1,iclu,ipla)-xpred
           call hf1(ihz+ipla,dx,1.)
          endif
          if (hpr_beam.ne.0) call hf1(ihz+20+ipla,clux(1,1,ipla),1.)
          enddo
         endif
        endif				! ipla.ne.nrff.and.ipla.ne.nrfl
        enddo
       endif				! nclux(nrff).eq.1.and.nclux(nrfl).le.1
c===     y-plane histogram filling
       if (ncluy(nrff).eq.1.and.ncluy(nrfl).eq.1) then
        y_slop=(cluy(1,1,nrfl)-cluy(1,1,nrff)-y_tabl_ref)/dzy
        yzero=cluy(1,1,nrff)+y_tabl-y_slop*zyplane(nrff)
        call hf1(ihz+5+nrff,cluy(1,1,nrff),1.)
        call hf1(ihz+5+nrfl,cluy(1,1,nrfl),1.)
        do ipla=1,5
        if (ipla.ne.nrff.and.ipla.ne.nrfl) then
         if (ncluy(ipla).ge.1) then
          do iclu=1,ncluy(ipla)
          if (cluy(2,iclu,ipla).le.(0.6*step(ipla))) then
           ypred=yzero+y_slop*zyplane(ipla)
           dy=cluy(1,iclu,ipla)-ypred
           if (ipla.le.4) dy=dy+y_tabl
           call hf1(ihz+5+ipla,dy,1.)
           if (hpr_beam.ne.0) call hf1(ihz+25+ipla,cluy(1,1,ipla),1.)
          endif
          enddo
         endif
        endif				! ilpa.ne.nrff.and.ipla.ne.nrfl
        enddo
       endif				! ncluy(nrff).eq.1.and.ncluy(nrfl).le.1
      else				! mode.eq.1
c
c===    calculate alignment constants
c
c===    get beam position in reference chambers
       nevx(nrff) = hsum(ihz+nrff)
       if (nevx(nrff).ne.0) then
         shiftx(nrff)=hstati(ihz+nrff,1,' ',0)
       else
         shiftx(nrff)=0.0
       endif
c
       nevy(nrff) = hsum(ihz+nrff+5)
       if (nevy(nrff).ne.0) then
         shifty(nrff)=hstati(ihz+5+nrff,1,' ',0)
       else
         shifty(nrff)=0.0
       endif
c
       nevx(nrfl) = hsum(ihz+nrfl)
       if (nevx(nrfl).ne.0) then
         shiftx(nrfl)=hstati(ihz+nrfl,1,' ',0)
       else
         shiftx(nrfl)=0.0
       endif
c
       nevy(nrfl) = hsum(ihz+nrfl+5)
       if (nevy(nrfl).ne.0) then
         shifty(nrfl)=hstati(ihz+5+nrfl,1,' ',0)
       else
         shifty(nrfl)=0.0
       endif
c
       beamx=shiftx(nrff)
       beamy=shifty(nrff)
c
       do ipla=1,5
       if (ipla.ne.nrff.and.ipla.ne.nrfl) then
        cont=hsum(ihz+ipla)
        nevx(ipla)=cont
        if (nevx(ipla).ge.1) then
c
c===     x-planes alignment     
c
         call hgive(ihz+ipla,chtit,nbin,x1,x2,idum,dum,dum,
     &              idum,idum)
         out=hi(ihz+ipla,0)+hi(ihz+ipla,nbin+1)
         if (out.gt.cont) then
          write(*,1020)ipla,nint(cont),nint(out)
 1020   format(//2x,'!!!!!!!!! x-plane #',i1,': Warning - too many '
     &  ,'entries out of survey histogram; cont=',i5,'  out=',i6//)
          print *,' ihz,ipla,nbin=',ihz,ipla,nbin,' cont=',cont,
     &            ' out=',out
         endif
         call hunpak(ihz+ipla,vect,' ',0)
         x=hstati(ihz+ipla,1,' ',0)
         call hxi(ihz+ipla,x,nblow)
         nbup=nblow+10
         nblow=nblow-10
         if (nblow.ge.1.and.nbup.le.50) then
          bwidth=(x2-x1)/nbin
          x=x1+(nblow-0.5)*bwidth
          cont=0.
          do i=nblow,nbup
          x=x+bwidth
          cont=cont+vect(i)
          shiftx(ipla)=shiftx(ipla)+x*vect(i)
          enddo
          shiftx(ipla)=shiftx(ipla)/cont
         else
          shiftx(ipla)=1001.
          print *,' *** survey failed for x-plane #',ipla
          print *,' <x>=',x,' nblow=',nblow,' nbup=',nbup
         endif
        endif				! nevx(ipla).ge.1
c
c===     y-planes alignment
c
        cont=hsum(ihz+5+ipla)
        nevy(ipla)=cont
        if (nevy(ipla).ge.1) then
         call hgive(ihz+5+ipla,chtit,nbin,x1,x2,idum,dum,dum,
     &              idum,idum)
         out=hi(ihz+5+ipla,0)+hi(ihz+5+ipla,nbin+1)
         if (out.gt.cont) then
          write(*,1021)ipla,ifix(cont),ifix(out)
 1021   format(//2x,'!!!!!!!!! y-plane #',i1,': worning - too many '
     &  ,'entries out of survey histogram; cont=',i5,'  out=',i6//)
          print *,' ihz, ipla=',ihz,ipla,' cont=',cont,' out=',out
         endif
         call hunpak(ihz+5+ipla,vect,' ',0)
         y=hstati(ihz+5+ipla,1,' ',0)
         call hxi(ihz+5+ipla,y,nblow)
         nbup=nblow+10
         nblow=nblow-10
         if (nblow.ge.1.and.nbup.le.50) then
          bwidth=(x2-x1)/nbin
          x=x1+(nblow-0.5)*bwidth
          cont=0.
          do i=nblow,nbup
          x=x+bwidth
          cont=cont+vect(i)
          shifty(ipla)=shifty(ipla)+x*vect(i)
          enddo
          shifty(ipla)=shifty(ipla)/cont
         else
          shifty(ipla)=1001.
          print *,' *** survey failed for y-plane #',ipla
          print *,' <y>=',x,' nblow=',nblow,' nbup=',nbup
         endif
        endif				! nevy(ipla).ge.1
       endif				! ipla.ne.nrff.and.ipla.ne.nrfl
       enddo
c
c      write (*,1011) shiftx
 1011 format(/5x,'***** alignment parameters (cm) before rotation',
     ,'*****'/,5x,'x-planes: ',5(f7.4,2x))
c      write (*,1002) shifty
c
c===    take into account reference chambers x,y positions
c
       dx=shiftx(nrfl)-shiftx(nrff)
       dy=shifty(nrfl)-shifty(nrff)
       shiftx(nrfl)=dx_refrun
       shifty(nrfl)=dy_refrun
       shiftx(nrff)=0.
       shifty(nrff)=0.
c      slopx=shiftx(nrfl)/dzx
c      xzero=-slopx*zxplane(nrff)
c      slopy=shifty(nrfl)/dzy
c      yzero=-slopy*zyplane(nrff)
c      do ipla=1,5
c      if (ipla.eq.nrff.or.ipla.eq.nrfl) go to 106
c      if (nevx(ipla).gt.0) then
c        shiftx(ipla)=shiftx(ipla)-slopx*zxplane(ipla)-xzero
c      else
c        shiftx(ipla)=0.0
c      endif
c      if (nevy(ipla).gt.0) then
c        shifty(ipla)=shifty(ipla)-slopy*zyplane(ipla)-yzero
c      else
c        shifty(ipla)=0.0
c      endif
c  106 enddo
c
c===    print resulting alignment parameters
c
       write (*,1001) shiftx
 1001  format(/5x,'***** alignment parameters (cm) *****'/,
     , 5x,'x-planes: ',5(f8.3,2x))
       write (*,1002) shifty
 1002  format(5x,'y-planes: ',5(f8.3,2x))
       write (*,1003) numev
 1003  format(/5x,'***** ',i5,' events were used for alignment')
       write (*,1210) nevx
 1210  format(/5x,'***** number of events for alignment *****'/
     & 5x,'x-planes: ',5(f8.0,2x))
       write (*,1220) nevy
 1220  format(5x,'y-planes: ',5(f8.0,2x))
       write (*,1012) dx,dy
 1012  format(/5x,'last reference mwpc displacements: dx=',f6.3,
     , '  dy=',f6.3//)
c
c===     write the computed constants to the file align.const
       if ( run_mode .eq. 1) then
         open (unit=io_beam(2), file='align.const', status='new')
       endif
       write (io_beam(2), *) run_no,
     ,     shiftx,shifty,beamx,beamy
c
       if (hpr_beam .ne. 0 .and. mk_his .gt. 0) then
        do ipla=1,5
        call hprint(ihz+ipla)
        call hprint(ihz+5+ipla)
        enddo
        do ipla=21,25
        call hprint(ihz+ipla)
        call hprint(ihz+5+ipla)
        enddo
       endif
       if (mk_his .eq. 0 .or. hpr_beam .eq. 0) then 
         do ipla=1,5
           call hdelet(ihz+ipla)
           call hdelet(ihz+5+ipla)
           call hdelet(ihz+20+ipla)
           call hdelet(ihz+25+ipla)
         enddo
       endif
      endif                     ! else mode.eq.1
c
      return
      end