1 cdaq  1.2       subroutine h_tof_fit(abort,errmsg,trk)

  2 cdaq  1.1 
  3           *-------------------------------------------------------------------
  4           * author: John Arrington
  5           * created: 2/22/94
  6           *
  7           * h_tof_fit fits the velocity of the paritcle from the corrected
  8           * times generated by h_tof.
  9           *
 10           * modifications:

 11 cdaq  1.2 * $Log: h_tof_fit.f,v $

 12 saw   1.10 * Revision 1.9  1995/05/22 19:39:29  cdaq
 13            * (SAW) Split gen_data_data_structures into gen, hms, sos, and coin parts"
 14            *

 15 cdaq  1.9  * Revision 1.8  1995/02/23  13:29:15  cdaq
 16            * (SAW) Cosmetic Changes
 17            *

 18 cdaq  1.8  * Revision 1.7  1995/02/10  18:49:57  cdaq
 19            * (JRA) Add track index to hgood_scin_time
 20            *

 21 cdaq  1.7  * Revision 1.6  1994/09/13  21:26:53  cdaq
 22            * (JRA) fix chisq calculation
 23            *

 24 cdaq  1.6  * Revision 1.5  1994/07/13  15:05:08  cdaq
 25            * (SAW) Add abs around tmpdenom that I left out last update
 26            *

 27 cdaq  1.5  * Revision 1.4  1994/07/11  18:34:35  cdaq
 28            * (JRA) Increase comparison of tmpdenom from 1e-15 to 1e-10
 29            *

 30 cdaq  1.4  * Revision 1.3  1994/07/08  19:42:31  cdaq
 31            * (JRA) Change fit from velocity to beta.  Bad fits give beta=0
 32            *

 33 cdaq  1.3  * Revision 1.2  1994/06/14  04:53:41  cdaq
 34            * (DFG) Protect against divide by 0 in beta calc
 35            *

 36 cdaq  1.2  * Revision 1.1  1994/04/13  16:29:15  cdaq
 37            * Initial revision
 38            *

 39 cdaq  1.1  *-------------------------------------------------------------------
 40            

 41 cdaq  1.2        implicit none

 42 cdaq  1.1  

 43 cdaq  1.9        include 'hms_data_structures.cmn'

 44 cdaq  1.2        include 'hms_scin_parms.cmn'
 45                  include 'hms_scin_tof.cmn'
 46            
 47                  logical abort
 48                  character*1024 errmsg
 49                  character*20 here
 50                  parameter (here = 'h_tof_fit')
 51                  
 52                  real*4 sumw, sumt, sumz, sumzz, sumtz
 53                  real*4 scin_weight
 54                  real*4 tmp, t0 ,tmpdenom
 55                  real*4 pathnorm
 56                  integer*4 hit, trk
 57                  save
 58                  
 59                  sumw = 0.
 60                  sumt = 0.
 61                  sumz = 0.
 62                  sumzz = 0.
 63                  sumtz = 0.

 64 cdaq  1.1  

 65 cdaq  1.2        do hit = 1 , hscin_tot_hits

 66 cdaq  1.8          if (hgood_scin_time(trk,hit)) then
 67                      scin_weight = 1./hscin_sigma(hit)**2
 68                      sumw = sumw + scin_weight
 69                      sumt = sumt + scin_weight * hscin_time(hit)
 70                      sumz = sumz + scin_weight * hscin_zpos(hit)
 71                      sumzz = sumzz + scin_weight * hscin_zpos(hit)**2
 72                      sumtz = sumtz + scin_weight * hscin_zpos(hit) *
 73                 1         hscin_time(hit)
 74                    endif

 75 cdaq  1.2        enddo

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

 84 cdaq  1.2        tmp = sumw*sumzz - sumz*sumz
 85                  t0 = (sumt*sumzz - sumz*sumtz) / tmp
 86                  tmpdenom = sumw*sumtz - sumz*sumt

 87 cdaq  1.5        if(abs(tmpdenom) .gt. 1.e-10) then

 88 cdaq  1.8          hbeta(trk) = tmp / tmpdenom     !velocity in cm/ns.
 89                    hbeta_chisq(trk) = 0.
 90                    do hit = 1 , hscin_tot_hits
 91                      if (hgood_scin_time(trk,hit)) then
 92                        hbeta_chisq(trk) = hbeta_chisq(trk) + 
 93                 1           (hscin_zpos(hit)/hbeta(trk) -
 94                 1           (hscin_time(hit) - t0))**2 / hscin_sigma(hit)**2
 95                      endif
 96                    enddo
 97            

 98 saw   1.10         pathnorm = sqrt(1 + hxp_fp(trk)**2 + hyp_fp(trk)**2)

 99 cdaq  1.8          hbeta(trk) = hbeta(trk) * pathnorm !take angle into account

100 saw   1.10         hbeta(trk) = hbeta(trk) / 29.979 !velocity/c

101 cdaq  1.2        else

102 cdaq  1.8          hbeta(trk) = 0.                 ! set unphysical beta
103                    hbeta_chisq(trk) = -2

104 cdaq  1.2        endif                             ! end if on denomimator = 0.

105 cdaq  1.1  

106 cdaq  1.2        return
107                  end