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