subroutine h_tof_fit(abort,errmsg,trk)

*-------------------------------------------------------------------
* author: John Arrington
* created: 2/22/94
*
* h_tof_fit fits the velocity of the paritcle from the corrected
* times generated by h_tof.
*
* modifications:
* $Log: h_tof_fit.f,v $
* Revision 1.10  1996/09/04 13:36:24  saw
* (JRA) Include actual beta in calculation of focal plane time.
*
* Revision 1.9  1995/05/22 19:39:29  cdaq
* (SAW) Split gen_data_data_structures into gen, hms, sos, and coin parts"
*
* Revision 1.8  1995/02/23  13:29:15  cdaq
* (SAW) Cosmetic Changes
*
* Revision 1.7  1995/02/10  18:49:57  cdaq
* (JRA) Add track index to hgood_scin_time
*
* Revision 1.6  1994/09/13  21:26:53  cdaq
* (JRA) fix chisq calculation
*
* Revision 1.5  1994/07/13  15:05:08  cdaq
* (SAW) Add abs around tmpdenom that I left out last update
*
* Revision 1.4  1994/07/11  18:34:35  cdaq
* (JRA) Increase comparison of tmpdenom from 1e-15 to 1e-10
*
* Revision 1.3  1994/07/08  19:42:31  cdaq
* (JRA) Change fit from velocity to beta.  Bad fits give beta=0
*
* Revision 1.2  1994/06/14  04:53:41  cdaq
* (DFG) Protect against divide by 0 in beta calc
*
* Revision 1.1  1994/04/13  16:29:15  cdaq
* Initial revision
*
*-------------------------------------------------------------------

      implicit none

      include 'hms_data_structures.cmn'
      include 'hms_scin_parms.cmn'
      include 'hms_scin_tof.cmn'

      logical abort
      character*1024 errmsg
      character*20 here
      parameter (here = 'h_tof_fit')
      
      real*4 sumw, sumt, sumz, sumzz, sumtz
      real*4 scin_weight
      real*4 tmp, t0 ,tmpdenom
      real*4 pathnorm
      integer*4 hit, trk
      save
      
      sumw = 0.
      sumt = 0.
      sumz = 0.
      sumzz = 0.
      sumtz = 0.

      do hit = 1 , hscin_tot_hits
        if (hgood_scin_time(trk,hit)) then
          scin_weight = 1./hscin_sigma(hit)**2
          sumw = sumw + scin_weight
          sumt = sumt + scin_weight * hscin_time(hit)
          sumz = sumz + scin_weight * hscin_zpos(hit)
          sumzz = sumzz + scin_weight * hscin_zpos(hit)**2
          sumtz = sumtz + scin_weight * hscin_zpos(hit) *
     1         hscin_time(hit)
        endif
      enddo

* The formula for beta (and t0) come from taking chi-squared (as
* defined below), and differentiating  with respect to each
* of the fit paramters (beta and t0 for fit to z=beta*(t-t0)).
* Setting both of these derivatives to zero gives the minumum
* chisquared (since they are quadratic in beta and t0), and
* gives a solution for beta in terms of sums of z, t, and w.

      tmp = sumw*sumzz - sumz*sumz
      t0 = (sumt*sumzz - sumz*sumtz) / tmp
      tmpdenom = sumw*sumtz - sumz*sumt
      if(abs(tmpdenom) .gt. 1.e-10) then
        hbeta(trk) = tmp / tmpdenom     !velocity in cm/ns.
        hbeta_chisq(trk) = 0.
        do hit = 1 , hscin_tot_hits
          if (hgood_scin_time(trk,hit)) then
            hbeta_chisq(trk) = hbeta_chisq(trk) + 
     1           (hscin_zpos(hit)/hbeta(trk) -
     1           (hscin_time(hit) - t0))**2 / hscin_sigma(hit)**2
          endif
        enddo

        pathnorm = sqrt(1 + hxp_fp(trk)**2 + hyp_fp(trk)**2)
        hbeta(trk) = hbeta(trk) * pathnorm !take angle into account
        hbeta(trk) = hbeta(trk) / 29.979 !velocity/c
      else
        hbeta(trk) = 0.                 ! set unphysical beta
        hbeta_chisq(trk) = -2
      endif                             ! end if on denomimator = 0.

      return
      end