subroutine s_tof_fit(abort,errmsg,trk)

*-------------------------------------------------------------------
* author: John Arrington
* created: 2/22/94
*
* s_tof_fit fits the velocity of the paritcle from the corrected
* times generated by s_tof.
*
* modifications:
* $Log: s_tof_fit.f,v $
* Revision 1.6  1996/09/04 20:36:55  saw
* (JRA) Don't forget the  sqrt in pathnorm
*
* Revision 1.5  1995/05/22 19:45:59  cdaq
* (SAW) Split gen_data_data_structures into gen, hms, sos, and coin parts"
*
* Revision 1.4  1995/02/23  13:28:35  cdaq
* (JRA) Add track index to sgood_scin_time
*
* Revision 1.3  1994/11/23  14:15:49  cdaq
* * (SPB) Recopied from hms file and modified names for SOS
*
* Revision 1.2  1994/06/14  04:36:30  cdaq
* (DFG) Protect against divide by 0 in beta calc
*
* Revision 1.1  1994/04/13  18:45:19  cdaq
* Initial revision
*
*-------------------------------------------------------------------

      implicit none

      include 'sos_data_structures.cmn'
      include 'sos_scin_parms.cmn'
      include 'sos_scin_tof.cmn'

      logical abort
      character*1024 errmsg
      character*20 here
      parameter (here = 's_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 , sscin_tot_hits
        if (sgood_scin_time(trk,hit)) then
          scin_weight = 1./sscin_sigma(hit)**2
          sumw = sumw + scin_weight
          sumt = sumt + scin_weight * sscin_time(hit)
          sumz = sumz + scin_weight * sscin_zpos(hit)
          sumzz = sumzz + scin_weight * sscin_zpos(hit)**2
          sumtz = sumtz + scin_weight * sscin_zpos(hit) *
     1         sscin_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(tmpdenom .gt. 1.e-10) then
        sbeta(trk) = tmp / tmpdenom     !velocity in cm/ns.
        sbeta_chisq(trk) = 0.
        do hit = 1 , sscin_tot_hits
          if (sgood_scin_time(trk,hit)) then
            sbeta_chisq(trk) = sbeta_chisq(trk) + 
     1           (sscin_zpos(hit)/sbeta(trk) -
     1           (sscin_time(hit) - t0))**2 / sscin_sigma(hit)**2
          endif
        enddo

        pathnorm = sqrt(1 + sxp_fp(trk)**2 + syp_fp(trk)**2)
        sbeta(trk) = sbeta(trk) * pathnorm !take angle into account
        sbeta(trk) = sbeta(trk) / 29.979 !velocity/c
      else
        sbeta(trk) = 0.                 !set unphysical beta
        sbeta_chisq(trk) = -2
      endif                             !end if on denomimator = 0.

      return
      end