SUBROUTINE H_TOF(ABORT,errmsg)
*--------------------------------------------------------
*-
*-   Purpose and Methods : Analyze scintillator information for each track 
*-
*-      Required Input BANKS     HMS_RAW_SCIN
*-                               HMS_DECODED_SCIN
*-                               HMS_FOCAL_PLANE
*-
*-      Output BANKS             HMS_TRACK_TESTS
*-
*-   Output: ABORT           - success or failure
*-         : err             - reason for failure, if any
*- 
* author: John Arrington
* created: 2/22/94
*
* h_tof finds the time of flight for a particle from
* the hodoscope TDC information.  It corrects for pulse
* height walk, time lag from the hit to the pmt signal,
* and time offsets for each signal.  It requires the
* hodoscope ADC and TDC information, the track, and
* the correction parameters.
*
* $Log: h_tof.f,v $
* Revision 1.19.6.2.2.7  2009/03/31 19:33:00  cdaq
* *** empty log message ***
*
* Revision 1.19.6.2.2.6  2008/12/11 17:11:43  cdaq
* For gfortran compiler the arguements generic functions
* must be the same type.
*
* Revision 1.19.6.2.2.5  2008/11/19 12:46:50  cdaq
* Add line to reutrn if nparam < 1
*
* Revision 1.19.6.2.2.4  2008/11/17 16:03:47  cdaq
* Fixed minor bug in printing of ADCHIST
*
* Revision 1.19.6.2.2.3  2008/11/17 16:00:22  cdaq
* Major revision to do tof calibration internatlly
*
* Revision 1.19.6.2.2.2  2008/10/28 21:03:18  cdaq
* Changed default betap to 1
*
* Revision 1.19.6.2.2.1  2008/10/27 16:34:54  cdaq
* changes for F1 TDCs
*
* Revision 1.19.6.2  2007/10/26 16:44:55  pcarter
* made the arguments to max() match data types -- GCC 4 is picky about that
*
* Revision 1.19.6.1  2007/10/24 16:37:16  cdaq
* *** empty log message ***
*
* Revision 1.19.4.3  2007/05/16 19:50:02  cdaq
* P. Bosted fixed bug in new code to dump TOF data
*
* Revision 1.19.4.2  2007/05/10 21:15:10  cdaq
* changes for writing dump file for Peter's tof code
*
* Revision 1.19.4.1  2007/05/02 21:19:30  jones
* Add new code needed for  adjusting scintillator timing using P Bosted's method.
*
* Revision 1.19  2005/03/15 21:08:08  jones
* Add code to filter the scintillator tdc hits and group them by time. ( P. Bosted)
*
* Revision 1.18  1999/06/10 16:52:12  csa
* (JRA) Cosmetic changes
*
* Revision 1.17  1997/03/19 18:43:45  saw
* (JRA) Don't neglect negative side of hodoscopes
*
* Revision 1.16  1996/09/04 13:36:00  saw
* (JRA) Include actual beta in calculation of focal plane time.
*
* Revision 1.15  1996/01/16 22:00:15  cdaq
* (JRA)
*
* Revision 1.14  1995/05/22 19:39:29  cdaq
* (SAW) Split gen_data_data_structures into gen, hms, sos, and coin parts"
*
* Revision 1.13  1995/02/10  18:59:41  cdaq
* (JRA) Add track index to hgood_plane_time, hgood_scin_time, hgood_tdc_pos,
* and  hgood_tdc_neg
*
* Revision 1.12  1995/02/02  16:35:25  cdaq
* (JRA) Zero out some variables at start, minph variables now per pmt,
*       hscin_adc_pos/neg change to floats.
*
* Revision 1.11  1995/01/31  21:49:32  cdaq
* (JRA) Added count of pmt's firing and cosmetic changes.
*
* Revision 1.10  1995/01/30  22:09:24  cdaq
* (JRA) Cosmetic changes.  Remove commented out code to dump time of
*       flight fitting data.
*
* Revision 1.9  1995/01/27  19:26:13  cdaq
* (JRA) Add calculation of time for each plane.  Add commented out
*       code to dump time of flight fitting data.
*
* Revision 1.8  1995/01/18  16:26:48  cdaq
* (SAW) Catch negative ADC values in argument of square root
*
* Revision 1.7  1994/09/13  21:25:35  cdaq
* (JRA) save extra diagnostic variables, require 2 hits/counter, add dedx
*
* Revision 1.6  1994/08/02  20:11:47  cdaq
* (JRA) Some hacks
*
* Revision 1.5  1994/07/21  13:29:45  cdaq
* (JRA) Correct sign on a time correction
*
* Revision 1.4  1994/07/08  19:43:53  cdaq
* (JRA) Keep list of wether hits are on track or not
*
* Revision 1.3  1994/05/13  02:36:30  cdaq
* (DFG) remove h_prt_track_tests call
*
* Revision 1.2  1994/04/13  16:28:53  cdaq
* (DFG) Add check for zero track
*
* Revision 1.1  1994/02/21  16:06:29  cdaq
* Initial revision
*
*--------------------------------------------------------
      IMPLICIT NONE
*
      character*50 here
      parameter (here= 'H_TOF')
*
      logical ABORT
      character*(*) errmsg
*
      INCLUDE 'hms_data_structures.cmn'
      INCLUDE 'gen_constants.par'
      INCLUDE 'gen_units.par'
      include 'hms_scin_parms.cmn'
      include 'hms_scin_tof.cmn'
      include 'hms_tracking.cmn'
      integer*4 hit, trk
      integer*4 plane,ind
      integer*4 hntof_pairs
      real*4 adc_ph                     !pulse height (channels)
      real*4 xhit_coord,yhit_coord
      real*4 time
      real*4 p,betap         !momentum and velocity from momentum, assuming desired mass
      real*4 path,zcor,num_fp_time
      real*4 sum_fp_time,sum_plane_time(hnum_scin_planes)
      integer*4 num_plane_time(hnum_scin_planes)
      integer timehist(200),i,j,jmax,maxhit,nfound
      real*4 time_pos(1000),time_neg(1000),tmin,time_tolerance
      logical keep_pos(1000),keep_neg(1000),first/.true./
      integer nsv, idetsv(100)
      real*8 tr0sv(100),psv(100),zcsv(100)
      real*8 tc1sv(100),adcsv(100)
      save
*
*--------------------------------------------------------
*
      ABORT= .FALSE.
      errmsg = ' '

      nsv = 0

c      if(hdumptof.eq.1) write(37,
c     > '(''ntrk,tothits''2i4,f8.3)') 
c     >  hntracks_fp,hscin_tot_hits,hscin_tdc_to_time
      if(hntracks_fp.le.0 .or. hscin_tot_hits.le.0) then
        do trk = 1 , hntracks_fp
          hnum_scin_hit(trk) = 0
          hnum_pmt_hit(trk) = 0
          hbeta(trk) = 0
          hbeta_chisq(trk) = -3
          htime_at_fp(trk) = 0
        enddo
        goto 666
      endif
 
**MAIN LOOP:  Loop over all tracks and get corrected time, tof, beta...
      do trk = 1 , hntracks_fp

** Initialize counter,flags...
        hntof = 0
        hntof_pairs = 0
        sum_fp_time = 0.
        num_fp_time = 0.
        hnum_scin_hit(trk) = 0
        hnum_pmt_hit(trk) = 0
        p = hp_tar(trk)
c if p=0, assume electrons
        if(abs(p).lt.0.1) then
          betap = 1.0
        else
          betap = p/sqrt(p*p+hpartmass*hpartmass)
c put in check for reasonable
          betap = min(1., max(0.3, betap))
        endif
        do plane = 1 , hnum_scin_planes
          hgood_plane_time(trk,plane) = .false.
          sum_plane_time(plane) = 0.
          num_plane_time(plane) = 0
        enddo

! Calculate all corrected hit times and histogram
! This uses a copy of code below. Results are save in time_pos,neg
! including the z-pos. correction assuming nominal value of betap
! Code is currently hard-wired to look for a peak in the
! range of 0 to 100 nsec, with a group of times that all
! agree withing a time_tolerance of time_tolerance nsec. The normal
! peak position appears to be around 35 nsec.
! NOTE: if want to find farticles with beta different than
!       reference particle, need to make sure this is big enough
!       to accomodate difference in TOF for other particles
! Default value in case user hasnt definedd something reasonable
        time_tolerance=20.0
        if(htof_tolerance.gt.0.5.and.htof_tolerance.lt.10000.) then
          time_tolerance=htof_tolerance
        endif
! Use wide window if dumping events for fitting
cc        if(hdumptof.eq.1) time_tolerance=25.0
        if(first) then
           first=.false.
           write(*,'(1x,''Using '',f8.2,'' nsec window for'',
     >     '' hms tof and fp calculations'')') time_tolerance
           if(hdumptof.eq.1) 
     >       write(*,'(/1x,''TOF calibration being done:'',
     >         ''  see output in HTOFCAL directory'')')
        endif
        nfound = 0
        do j=1,200
          timehist(j)=0
        enddo
        do hit = 1 , hscin_tot_hits
          i=min(1000,hit)
          time_pos(i)=-99.
          time_neg(i)=-99.
          keep_pos(i)=.false.
          keep_neg(i)=.false.
          plane = hscin_plane_num(hit)
          xhit_coord = hx_fp(trk) + hxp_fp(trk)*hscin_zpos(hit)
          yhit_coord = hy_fp(trk) + hyp_fp(trk)*hscin_zpos(hit)
          if (plane.eq.1 .or. plane.eq.3) then !x plane
            hscin_trans_coord(hit) = xhit_coord
            hscin_long_coord(hit) = yhit_coord
          else if (plane.eq.2 .or. plane.eq.4) then !y plane
            hscin_trans_coord(hit) = yhit_coord
            hscin_long_coord(hit) = xhit_coord
          else                          !bad plane #.
            abort = .true.
            write(errmsg,*) 'hscin_plane_num(',hit,') = ',plane
            call g_prepend(here,errmsg)
            return
          endif
          if (abs(hscin_center_coord(hit)-hscin_trans_coord(hit))
     &         .lt.(hscin_width(hit)/2.+hscin_slop(hit))) then
            if(hscin_tdc_pos(hit) .ge. hscin_tdc_min .and.  
     &          hscin_tdc_pos(hit) .le. hscin_tdc_max) then
              adc_ph = hscin_adc_pos(hit)
              path = hscin_pos_coord(hit) - hscin_long_coord(hit)
              time = hscin_tdc_pos(hit) * hscin_tdc_to_time
              time = time - (hscin_zpos(hit)/(29.979*betap) *
     &               sqrt(1. + hxp_fp(trk)**2 + hyp_fp(trk)**2))
              time_pos(i) = time - hscin_pos_invadc_offset(hit) -
     >            path / hscin_pos_invadc_linear(hit) -
     >            hscin_pos_invadc_adc(hit)/sqrt(max(20.,adc_ph))
              nfound = nfound + 1
              do j=1,200
                tmin = 0.5*float(j)                
                if(time_pos(i) .gt. tmin .and.
     >             time_pos(i) .lt. tmin + time_tolerance) 
     >            timehist(j) = timehist(j) + 1
              enddo
            endif
            if (hscin_tdc_neg(hit).ge.hscin_tdc_min .and. !good tdc
     1           hscin_tdc_neg(hit).le.hscin_tdc_max) then
              adc_ph = hscin_adc_neg(hit)
              path = hscin_long_coord(hit) - hscin_neg_coord(hit)
              time = hscin_tdc_neg(hit) * hscin_tdc_to_time
              time = time - (hscin_zpos(hit)/(29.979*betap) *
     &               sqrt(1. + hxp_fp(trk)**2 + hyp_fp(trk)**2))
              time_neg(i) = time - hscin_neg_invadc_offset(hit) -
     >            path / hscin_neg_invadc_linear(hit) -
     >            hscin_neg_invadc_adc(hit)/sqrt(max(20.,adc_ph))
              nfound = nfound + 1
              do j=1,200
                tmin = 0.5*float(j)                
                if(time_neg(i) .gt. tmin .and.
     >             time_neg(i) .lt. tmin + time_tolerance) 
     >            timehist(j) = timehist(j)+1
             enddo
            endif
          endif
        enddo
! Find bin with most hits
        jmax=0
        maxhit=0
        do j=1,200
          if(timehist(j) .gt. maxhit) then
            jmax = j
            maxhit = timehist(j)
          endif
        enddo
c$$$        if(hdumptof.eq.1) then
c$$$          write(37,'(''trk='',2i3,8f8.3)') trk,jmax,
c$$$     >      hx_fp(trk),hxp_fp(trk),hy_fp(trk),hyp_fp(trk),
c$$$     >      hp_tar(trk)
c$$$        endif
        if(jmax.gt.0) then
          tmin = 0.5*float(jmax) 
          do hit = 1 , hscin_tot_hits
            i=min(1000,hit)
            if(time_pos(i) .gt. tmin .and.
     >         time_pos(i) .lt. tmin + time_tolerance) then
              keep_pos(i) = .true.
           endif
            if(time_neg(i) .gt. tmin .and.
     >         time_neg(i) .lt. tmin + time_tolerance) then
              keep_neg(i) = .true.
            endif
          enddo
        endif

! Resume regular tof code, now using time filer from above
        do hit = 1 , hscin_tot_hits
          hgood_scin_time(trk,hit) = .false.
          hgood_tdc_pos(trk,hit) = .false.
          hgood_tdc_neg(trk,hit) = .false.
          hscin_time(hit) = 0
          hscin_sigma(hit) = 100.
        enddo

        do hit = 1 , hscin_tot_hits
          plane = hscin_plane_num(hit)

**    Find hit position
          xhit_coord = hx_fp(trk) + hxp_fp(trk)*hscin_zpos(hit)
          yhit_coord = hy_fp(trk) + hyp_fp(trk)*hscin_zpos(hit)
          if (plane.eq.1 .or. plane.eq.3) then !x plane
            hscin_trans_coord(hit) = xhit_coord
            hscin_long_coord(hit) = yhit_coord
          else if (plane.eq.2 .or. plane.eq.4) then !y plane
            hscin_trans_coord(hit) = yhit_coord
            hscin_long_coord(hit) = xhit_coord
          else                          !bad plane #.
            abort = .true.
            write(errmsg,*) 'hscin_plane_num(',hit,') = ',plane
            call g_prepend(here,errmsg)
            return
          endif

**    Check if scin is on track
          if (abs(hscin_center_coord(hit)-hscin_trans_coord(hit))
     &         .gt.(hscin_width(hit)/2.+hscin_slop(hit))) then

            hscin_on_track(trk,hit) = .false.
          else
            hscin_on_track(trk,hit) = .true.
**    Check for good TDC
            if (hscin_tdc_pos(hit) .ge. hscin_tdc_min .and.  
     &          hscin_tdc_pos(hit) .le. hscin_tdc_max .and.
     >          keep_pos(hit)) then

**    Calculate time for each tube with a good tdc. 'pos' side first.
              hgood_tdc_pos(trk,hit) = .true.
              hntof = hntof + 1
              adc_ph = hscin_adc_pos(hit)
              path = hscin_pos_coord(hit) - hscin_long_coord(hit)

*     Convert TDC value to time, do pulse height correction, correction for
*     propogation of light thru scintillator, and offset.
              time = hscin_tdc_pos(hit) * hscin_tdc_to_time
              hscin_pos_time(hit)=time - hscin_pos_invadc_offset(hit) -
     >            path / hscin_pos_invadc_linear(hit) -
     >            hscin_pos_invadc_adc(hit)/sqrt(max(20.,adc_ph))
              zcor =  (hscin_zpos(hit)/(29.979*betap) * sqrt(1.+
     >               hxp_fp(trk)*hxp_fp(trk)+hyp_fp(trk)*hyp_fp(trk)))
              if(hdumptof.eq.1) then
c$$$                write(37,'(1x,''1'',2i3,5f10.3)') 
c$$$     >             hscin_plane_num(hit),
c$$$     >             hscin_counter_num(hit),
c$$$     >             hscin_tdc_pos(hit) * hscin_tdc_to_time,
c$$$     >             path,zcor,
c$$$     >             hscin_pos_time(hit)-zcor,adc_ph
                nsv = min(100, nsv + 1)
                idetsv(nsv) = 20 * (hscin_plane_num(hit)-1) +
     >            hscin_counter_num(hit)
                tr0sv(nsv) = hscin_tdc_pos(hit) * 
     >            hscin_tdc_to_time
                psv(nsv) = path
                zcsv(nsv) = zcor
                tc1sv(nsv) = hscin_pos_time(hit)-zcor
                adcsv(nsv) = adc_ph
              endif
            endif

**    Repeat for pmts on 'negative' side
            if (hscin_tdc_neg(hit).ge.hscin_tdc_min .and. !good tdc
     1           hscin_tdc_neg(hit).le.hscin_tdc_max.and.
     >           keep_neg(hit)) then

              hgood_tdc_neg(trk,hit) = .true.
              hntof = hntof + 1
              adc_ph = hscin_adc_neg(hit)
              path = hscin_long_coord(hit) - hscin_neg_coord(hit)
              time = hscin_tdc_neg(hit) * hscin_tdc_to_time
              hscin_neg_time(hit)=time - hscin_neg_invadc_offset(hit) -
     >            path / hscin_neg_invadc_linear(hit) -
     >            hscin_neg_invadc_adc(hit)/sqrt(max(20.,adc_ph))
              zcor =  (hscin_zpos(hit)/(29.979*betap) * sqrt(1.+
     >               hxp_fp(trk)*hxp_fp(trk)+hyp_fp(trk)*hyp_fp(trk)))
c$$$              if(hdumptof.eq.1.and.hntracks_fp.eq.1.and.
c$$$     $             timehist(max(1,jmax)).gt.6) then
               if(hdumptof.eq.1) then
c$$$                write(37,'(1x,''2'',2i3,5f10.3)') 
c$$$     >             hscin_plane_num(hit),
c$$$     >             hscin_counter_num(hit),
c$$$     >             hscin_tdc_neg(hit) * hscin_tdc_to_time,
c$$$     >             path,zcor,
c$$$     >             hscin_neg_time(hit)-zcor,adc_ph
                nsv = min(100, nsv + 1)
                idetsv(nsv) = 20 * (hscin_plane_num(hit)-1) +
     >            hscin_counter_num(hit) + 100
                tr0sv(nsv) = hscin_tdc_neg(hit) * 
     >            hscin_tdc_to_time
                psv(nsv) = path
                zcsv(nsv) = zcor
                tc1sv(nsv) = hscin_neg_time(hit)-zcor
                adcsv(nsv) = adc_ph
              endif
            endif

**    Calculate ave time for scintillator and error.
            if (hgood_tdc_pos(trk,hit)) then
              if (hgood_tdc_neg(trk,hit)) then
                hscin_time(hit) = (hscin_neg_time(hit) + hscin_pos_time(hit))/2.
                hscin_sigma(hit) = max(0.1,sqrt(hscin_neg_sigma(hit)**2 + 
     1               hscin_pos_sigma(hit)**2)/2.)
                hgood_scin_time(trk,hit) = .true.
                hntof_pairs = hntof_pairs + 1
              else
                hscin_time(hit) = hscin_pos_time(hit)
                hscin_sigma(hit) = max(0.1,hscin_pos_sigma(hit))
                hgood_scin_time(trk,hit) = .true.
*                hgood_scin_time(trk,hit) = .false.
              endif
            else                        ! if hgood_tdc_neg = .false.
              if (hgood_tdc_neg(trk,hit)) then
                hscin_time(hit) = hscin_neg_time(hit)
                hscin_sigma(hit) = max(0.1,hscin_neg_sigma(hit))
                hgood_scin_time(trk,hit) = .true.
*                hgood_scin_time(trk,hit) = .false.
              endif
            endif
c     Get time at focal plane
            if (hgood_scin_time(trk,hit)) then
              hscin_time_fp(hit) = hscin_time(hit)
     &             - (hscin_zpos(hit)/(29.979*betap) *
     &             sqrt(1.+hxp_fp(trk)*hxp_fp(trk)+hyp_fp(trk)*hyp_fp(trk)))
              sum_fp_time = sum_fp_time + hscin_time_fp(hit) /
     >          hscin_sigma(hit)**2
              num_fp_time = num_fp_time + 1./hscin_sigma(hit)**2
              sum_plane_time(plane)=sum_plane_time(plane)
     &             +hscin_time_fp(hit)
              num_plane_time(plane)=num_plane_time(plane)+1
              hnum_scin_hit(trk) = hnum_scin_hit(trk) + 1
              hscin_hit(trk,hnum_scin_hit(trk)) = hit
              hscin_fptime(trk,hnum_scin_hit(trk)) = hscin_time_fp(hit)

              if (hgood_tdc_pos(trk,hit) .and. hgood_tdc_neg(trk,hit)) then
                hnum_pmt_hit(trk) = hnum_pmt_hit(trk) + 2
              else
                hnum_pmt_hit(trk) = hnum_pmt_hit(trk) + 1
              endif
              if (hgood_tdc_pos(trk,hit)) then
                if (hgood_tdc_neg(trk,hit)) then
                  hdedx(trk,hnum_scin_hit(trk)) = sqrt(max(0.,
     &                 hscin_adc_pos(hit)*hscin_adc_neg(hit)))
                else
                  hdedx(trk,hnum_scin_hit(trk))=max(0.,hscin_adc_pos(hit))
                endif
              else
                if (hgood_tdc_neg(trk,hit)) then
                  hdedx(trk,hnum_scin_hit(trk))=max(0.,hscin_adc_neg(hit))
                else
                  hdedx(trk,hnum_scin_hit(trk)) = 0.
                endif
              endif
            endif
            
          endif                 !end of 'if scintillator was on the track'
          
**    See if there are any good time measurements in the plane.
          if (hgood_scin_time(trk,hit)) then
            hgood_plane_time(trk,plane) = .true. !still in loop over hits.
          endif
          
        enddo                           !end of loop over hit scintillators

**    Fit beta if there are enough time measurements (one upper, one lower)
        if ((hgood_plane_time(trk,1) .or. hgood_plane_time(trk,2)) .and.
     1       (hgood_plane_time(trk,3) .or. hgood_plane_time(trk,4))) then
          call h_tof_fit(abort,errmsg,trk) !fit velocity of particle
          if (abort) then
            call g_prepend(here,errmsg)
            return
          endif
        else             !cannot fit beta from given time measurements
          hbeta(trk) = 0.
          hbeta_chisq(trk) = -1.
        endif
        if (num_fp_time .ne. 0) then
          htime_at_fp(trk) = sum_fp_time / num_fp_time
        endif
        
        do ind=1,4
          if (num_plane_time(ind) .ne. 0) then
            h_fptime(ind)=sum_plane_time(ind)/float(num_plane_time(ind))
          else
            h_fptime(ind)=1000.*ind
          endif
        enddo

        h_fptimedif(1)=h_fptime(1)-h_fptime(2)
        h_fptimedif(2)=h_fptime(1)-h_fptime(3)
        h_fptimedif(3)=h_fptime(1)-h_fptime(4)
        h_fptimedif(4)=h_fptime(2)-h_fptime(3)
        h_fptimedif(5)=h_fptime(2)-h_fptime(4)
        h_fptimedif(6)=h_fptime(3)-h_fptime(4)   
*     
*     Dump tof common blocks if (hdebugprinttoftracks is set

        if(hdebugprinttoftracks.ne.0 ) call h_prt_tof(trk)


        if(hntracks_fp.gt.1000) then
          if(trk.eq.1) write(*,'(/1x,''hms tol='',f8.2)') time_tolerance
          write(*,'(5i3,4L2,7f7.2)') trk,nfound,jmax,timehist(max(1,jmax)),
     >      hnum_pmt_hit(trk),
     >      hgood_plane_time(trk,1),hgood_plane_time(trk,3),
     >      hgood_plane_time(trk,2),hgood_plane_time(trk,4),
     >      htime_at_fp(trk),hbeta(trk),hbeta_chisq(trk),
     >      hdelta_tar(trk),hy_tar(trk),hxp_tar(trk),hyp_tar(trk)
        endif
        if(hntracks_fp.eq.1.and.
     >    hdumptof.eq.1.and.
     >     timehist(max(1,jmax)).gt.6) then
c        if(hdumptof.eq.1) then
           write(37,'(1x,''0'',2i3,5f10.3)') trk,hntracks_fp,
     >      p, betap
          call h_tofcal_fill(nsv,idetsv,tr0sv,psv,zcsv,
     >       tc1sv,adcsv)
        endif
      enddo                             !end of loop over tracks

 666  continue

      RETURN
      END



! Fit TOF for Hall C HMS with the form for each PMT:
! tcorr = time - offset - path * velocity - adccor / sqrt(ADC)
! where offset, velocity, and adccor are parameters
! September 20085 P. Bosted
! Modified to run automatically during replay: no longer
! any need to dump large text files
! To activate, set hdumptof = 1 in hdebug.param in
! the PARAM directory (can also do from command line)
! The output parameters will be in tof/hodoxxxx.param, where
! xxxxx is the run number
! Normal values of invadc_offset are between -50 and 50,
! normal values of invadc_velocity are between 12 and 17 (50 is
! default if not enough data for the fit), and normal values of
! shodo_pos_invadc_adc are 20 to 50. Normal values of the sigmas
! are 0.3 to 0.8 nsec. 

      subroutine h_tofcal_init
! initialize common block variables at begin run
      implicit none
      integer i,j
! common block variables
      integer thist(200,10),adchist(200,18),phist(200,18)
      integer nhit(200),ip1(200),ip2(200),ip3(200)
      integer ipdet(600),iptyp(600),nparam
      real*8 ax(1000,1000),bx(1000),avtime,ntime,avsig(200)
      common/htofcal/ ax,bx,thist,adchist,phist,nhit,
     >   ip1,ip2,ip3,ipdet,iptyp,avtime,ntime,avsig,nparam

      nparam=0
      do i=1,200
        nhit(i)=0
        ip1(i)=0
        ip2(i)=0
        ip3(i)=0
        avsig(i)=0.
        do j=1,10
          thist(i,j)=0
        enddo
        do j=1,18
          adchist(i,j)=0
          phist(i,j)=0
        enddo
      enddo
      do i=1,600
        ipdet(i)=0
        iptyp(i)=0
      enddo

! Initialize the fitting arrays
      do i=1,1000
        bx(i)=0.
        do j=1,1000
          ax(i,j)=0.
        enddo
      enddo

      avtime = 0.
      ntime = 0.

      return
      end

      subroutine h_tofcal_fill(n,idet,tr0,p,zc,tc1,adc)
! Fill in the arrays for HMS tof cal
! Inputs are:
! n  number of PMTs n
! idet detector code (from 1 to 200)
! tr0 TDCtime
! p path length
! zc time correction due to z 
! tc1 corrected time using current variables
! adc ADC

      implicit none

! local and input variables
      integer i,j,k,n,idt,idet(100)
      integer k1,k2,k3,k4,k5,k6
      real*8 tr0(100),p(100),zc(100),tc1(100),p2(100),adc(100),tr(100)
      real*8 dt,avval
      logical first_time/.true./
! common block variables
      integer thist(200,10),adchist(200,18),phist(200,18)
      integer nhit(200),ip1(200),ip2(200),ip3(200)
      integer ipdet(600),iptyp(600),nparam
      real*8 ax(1000,1000),bx(1000),avtime,ntime,avsig(200)
      common/htofcal/ ax,bx,thist,adchist,phist,nhit,
     >   ip1,ip2,ip3,ipdet,iptyp,avtime,ntime,avsig,nparam

      if(first_time) then
        first_time = .false.
        call h_tofcal_init
      endif

! need at least 6 PMTs for fitting
      if(n.le.5) return

c      write(6,'(/i3,10f6.1)') n,(tc1(i),i=1,min(10,n))
c      write(6,'(i3,10i6)') n,(idet(i),i=1,min(10,n))
! Loop over all PMTs
      avval = 0. 
      do i=1,n
! Check for valid detector code
        if(idet(i).le.0.or.idet(i).gt.200) then
          write(6,'(''ERROR, in h_tofcal, idet='',2i8)') i,idet(i)
          return
        endif

! Fill in ADC histograms
        k = min(18., max(1., (adc(i)/20.)+1))
        adchist(idet(i), k) = adchist(idet(i), k) + 1 

! correct raw times for zpos using betap
        tr(i) = tr0(i)- zc(i)

! Put 1./sqrt(ADC) in p2 variable
        p2(i) = 1./sqrt(max(20., adc(i)))

! Histogram path length variable
        k = min(18, max(1, int(p(i)/7.)+1))
        phist(idet(i), k) = phist(idet(i), k) + 1 
! average time
        avval = avval + tc1(i)
      enddo
      avval = avval / float(n)

! Loop over PMTS again
! Get average h_start_time and sigmas for each PMT
! THESE SHOULD BE DONE on a 2nd iteration of the
! TOF calibration for a given run, setting
! h_tof_tolerance to something small like 3 nsec
      do j=1,n
        nhit(idet(j))=nhit(idet(j))+1
        avsig(idet(j)) = avsig(idet(j)) + 
     >    (tc1(j) - avval)**2

! If first time detector used, assign corresponding parameters
! Note that detector 4 had has a fixed time offset (ip1) of zero
! since all times are relative. 
        if(nhit(idet(j)).eq.1) then
          if(idet(j).eq.4) then
            ip1(idet(j))=0
          else
! fixed time offsets
            nparam=nparam+1
            ip1(idet(j))=nparam
            ipdet(nparam)=idet(j)
            iptyp(nparam)=1
          endif

! linear term in path
! Changed 11/08 to make same for both pos. and neg. ends! 
          if(idet(j).le.100) then
           nparam=nparam+1
           ip2(idet(j))=nparam
           ip2(idet(j)+100)=nparam
           ipdet(nparam)=idet(j)
           iptyp(nparam)=2
          endif

! 1/sqrt(adc) terms (or could be path length**2 if desired)
          nparam=nparam+1
          ip3(idet(j))=nparam
          ipdet(nparam)=idet(j)
          iptyp(nparam)=3
          k=idet(j)
c          write(6,'(''h_tofcal_fill'',i3,4i5)') 
c     >      k,nhit(k),ip1(k),ip2(k),ip3(k)
        endif
        avtime = avtime + tc1(j)
        ntime = ntime + 1.
       enddo ! loop over n

! now loop over all pairs in fill in the matrix
! Also, histogram time differences using current corrections
       do j=1,n-1
          do k=j+1,n
            if(ip2(idet(j)).gt.0.and.ip2(idet(k)).gt.0) then
              dt = tc1(j)-tc1(k)
              idt = min(10,max(1,int((dt+5.))+1))
              thist(idet(j),idt) = thist(idet(j),idt) + 1
              dt = tc1(k)-tc1(j)
              idt = min(10,max(1,int((dt+5.))+1))
              thist(idet(k),idt) = thist(idet(k),idt) + 1
              k1 = idet(j)
              k2 = idet(k)
              k1=ip1(idet(j))
              k2=ip1(idet(k))
              k3=ip2(idet(j))
              k4=ip2(idet(k))
              k5=ip3(idet(j))
              k6=ip3(idet(k))
              if(k1.gt.0) then
                bx(k1) = bx(k1) - (tr(j)-tr(k))
                ax(k1,k1) = ax(k1,k1) + 1.
                ax(k1,k3) = ax(k1,k3) + p(j)
                ax(k1,k4) = ax(k1,k4) - p(k)
                ax(k1,k5) = ax(k1,k5) + p2(j)
                ax(k1,k6) = ax(k1,k6) - p2(k)
                ax(k3,k1) = ax(k3,k1) + p(j)
                ax(k4,k1) = ax(k4,k1) - p(k)
                ax(k5,k1) = ax(k5,k1) + p2(j)
                ax(k6,k1) = ax(k6,k1) - p2(k)
              endif
              if(k1.gt.0.and.k2.gt.0) then
                ax(k1,k2) = ax(k1,k2) - 1.
                ax(k2,k1) = ax(k2,k1) - 1.
              endif
              if(k2.gt.0) then
                bx(k2) = bx(k2) + (tr(j)-tr(k))
                ax(k2,k2) = ax(k2,k2) + 1.
                ax(k2,k3) = ax(k2,k3) - p(j)
                ax(k2,k4) = ax(k2,k4) + p(k)
                ax(k2,k5) = ax(k2,k5) - p2(j)
                ax(k2,k6) = ax(k2,k6) + p2(k)
                ax(k3,k2) = ax(k3,k2) - p(j)    
                ax(k4,k2) = ax(k4,k2) + p(k)    
                ax(k5,k2) = ax(k5,k2) - p2(j)    
                ax(k6,k2) = ax(k6,k2) + p2(k)    
              endif
              bx(k3) = bx(k3) - (tr(j)-tr(k)) * p(j)
              bx(k4) = bx(k4) + (tr(j)-tr(k)) * p(k)
              bx(k5) = bx(k5) - (tr(j)-tr(k)) * p2(j)
              bx(k6) = bx(k6) + (tr(j)-tr(k)) * p2(k)
              ax(k3,k3) = ax(k3,k3) + p(j)*p(j)
              ax(k3,k4) = ax(k3,k4) - p(k)*p(j)
              ax(k3,k5) = ax(k3,k5) + p2(j)*p(j)
              ax(k3,k6) = ax(k3,k6) - p2(k)*p(j)
              ax(k4,k3) = ax(k4,k3) - p(j)*p(k)
              ax(k4,k4) = ax(k4,k4) + p(k)*p(k)
              ax(k4,k5) = ax(k4,k5) - p2(j)*p(k)
              ax(k4,k6) = ax(k4,k6) + p2(k)*p(k)
              ax(k5,k3) = ax(k5,k3) + p(j)*p2(j)
              ax(k5,k4) = ax(k5,k4) - p(k)*p2(j)
              ax(k5,k5) = ax(k5,k5) + p2(j)*p2(j)
              ax(k5,k6) = ax(k5,k6) - p2(k)*p2(j)
              ax(k6,k3) = ax(k6,k3) - p(j)*p2(k)
              ax(k6,k4) = ax(k6,k4) + p(k)*p2(k)
              ax(k6,k5) = ax(k6,k5) - p2(j)*p2(k)
              ax(k6,k6) = ax(k6,k6) + p2(k)*p2(k)
           endif
         enddo
      enddo


      return
      end


      subroutine h_tofcal_endrun(runno)
! Solve simultaneous linear equations for best values of
! tof parameters, and write to file
! local and input variables
      implicit none

      INCLUDE 'hms_data_structures.cmn'
      include 'hms_scin_parms.cmn'
      include 'hms_scin_tof.cmn'
      integer i,j,runno,iwork(1000),ifail
      real*8 toff(200),vel(200),quad(200)
      character*80 fname
! common block variables
      integer thist(200,10),adchist(200,18),phist(200,18)
      integer nhit(200),ip1(200),ip2(200),ip3(200)
      integer ipdet(600),iptyp(600),nparam
      real*8 ax(1000,1000),bx(1000),avtime,ntime,avsig(200)
      common/htofcal/ ax,bx,thist,adchist,phist,nhit,
     >   ip1,ip2,ip3,ipdet,iptyp,avtime,ntime,avsig,nparam


! find the solutions
      if (nparam .lt. 1) return
      call deqn (nparam,ax,1000,iwork,ifail,1,bx)

! association of parameters with detectors
      do i=1,200
        toff(i)=0
        vel(i)=0.
        quad(i)=0.
      enddo

      do i=1,nparam
        if(iptyp(i).eq.1) toff(ipdet(i))=bx(i)
        if(iptyp(i).eq.2)  vel(ipdet(i))=bx(i)
        if(iptyp(i).eq.2)  vel(ipdet(i)+100)=bx(i)
        if(iptyp(i).eq.3) quad(ipdet(i))=bx(i)
      enddo


! write solutions
      write(fname,'(''HTOFCAL/htofcal'',i5.5,''.param'')') 
     >  runno

      open(unit=10,file=fname)

      write(10,'(''; This parameter determines how close'',
     >  '' in time the ''/
     >  ''; corrected scint. have to be to each other. '',
     >  '' For initial calibrations, use 50.  ''/
     >  ''; For final calibration,use 3. For regular '',
     >  '' running, use about 10.''/
     >  ''; Used in h_trans_scin.f and h_tof.f''/
     >  ''   htof_tolerance = '',f6.1)') 
     >  htof_tolerance

      write(10,'(/''; This is default (average) value'',
     >  '' of start time for ''/
     >  ''; drift chambers. It is used in h_trans_scin.f''/
     >  ''   hstart_time_center = '',f6.1)') avtime/
     >  max(1.,ntime)

! copied from previous hhodo.param
      write(10,'(/''; This is 1/2 width of winow'',
     >  '' on hstart_time_center''/
     >  ''; it is used in h_trans_scin.f''/
     >  ''   hstart_time_slop = '',f6.1)') 
     >  hstart_time_slop

      write(10,'(/''; Minimum and Maximum raw TDC'',
     >  '' that will be used ''/
     >  ''; Check raw TDC spectra to make sure ok''/
     >  ''   hscin_tdc_min = '',i6/ 
     >  ''   hscin_tdc_max = '',i6)') 
     >  int(hscin_tdc_min),int(hscin_tdc_max)

      write(10,'(/''; TDC time in nsec per channel''/
     >  ''   hscin_tdc_to_time = '',f8.4)') 
     >  hscin_tdc_to_time

      write(10,'(/''; Position tolerance in cm'',
     >  '' used in efficiency calcultions ''/
     >  ''; used in hms_scin_eff.f''/
     >  ''   hhodo_slop = '',
     >  f5.0,'','',f5.0,'','',f5.0,'','',f5.0)') 
     >  (hhodo_slop(i),i=1,4)

      write(10,'(/''hhodo_pos_invadc_offset ='',3(f8.2,'',''),
     >  f8.2)') (-1.0*toff(i),i=1,80,20)
      do j=2,16
       write(10,'(1x,''                        '',3(f8.2,'',''),
     >  f8.2)')(-1.0*toff(i),i=j,79+j,20)
      enddo

      write(10,'(/''hhodo_neg_invadc_offset ='',3(f8.2,'',''),
     >  f8.2)')(-1.0*toff(i),i=101,180,20)
      do j=2,16
       write(10,'(1x,''                        '',3(f8.2,'',''),
     >  f8.2)')(-1.0*toff(i),i=100+j,179+j,20)
      enddo

      write(10,'(/''hhodo_pos_invadc_linear ='',3(f8.2,'',''),
     >  f8.2)')( -1./min(-0.02,vel(i)),i=1,80,20)
      do j=2,16
       write(10,'(1x,''                        '',3(f8.2,'',''),
     >  f8.2)')(-1./min(-0.02,vel(i)),i=j,79+j,20)
      enddo

      write(10,'(/''hhodo_neg_invadc_linear ='',3(f8.2,'',''),
     >  f8.2)')( -1./min(-0.02,vel(i)),i=101,180,20)
      do j=2,16
       write(10,'(1x,''                        '',3(f8.2,'',''),
     >  f8.2)')(-1./min(-0.02,vel(i)),i=100+j,179+j,20)
      enddo

      write(10,'(/''hhodo_pos_invadc_adc='',3(f9.2,'',''),
     >  f9.2)')(-1.*quad(i),i=1,80,20)
      do j=2,16
       write(10,'(1x,''                   '',3(f9.2,'',''),
     >  f9.2)')(-1.*quad(i),i=j,79+j,20)
      enddo

      write(10,'(/''hhodo_neg_invadc_adc='',3(f9.2,'',''),
     >  f9.2)')(-1.0*quad(i),i=101,180,20)
      do j=2,16
       write(10,'(1x,''                   '',3(f9.2,'',''),
     >  f9.2)')(-1.*quad(i),i=100+j,179+j,20)
      enddo

      do i=1,200
        avsig(i) = avsig(i) / max(1,nhit(i))
      enddo

      write(10,'(/''hhodo_pos_sigma ='',3(f8.2,'',''),
     >  f8.2)') (max(0.1,min(100.,avsig(i))),i=1,80,20)
      do j=2,16
       write(10,'(1x,''               '',3(f8.2,'',''),
     >  f8.2)')(max(0.1,min(100.,avsig(i))),i=j,79+j,20)
      enddo

      write(10,'(/''hhodo_neg_sigma ='',3(f8.2,'',''),
     >  f8.2)') (max(0.3,min(100.,avsig(i))),i=101,180,20)
      do j=2,16
       write(10,'(1x,''               '',3(f8.2,'',''),
     >  f8.2)')(max(0.3,min(100.,avsig(i))),i=100+j,179+j,20)
      enddo

      write(10,'(/''hhodo_pos_ped_limit = 1000,'',
     >  ''1000,1000,1000,1000,1000,1000,1000''/
     >  22x,''1000,1000,1000,1000,1000,1000,1000,1000''/   
     >  22x,''1000,1000,1000,1000,1000,1000,1000,1000''/   
     >  22x,''1000,1000,1000,1000,1000,1000,1000,1000''/   
     >  22x,''1000,1000,1000,1000,1000,1000,1000,1000''/   
     >  22x,''1000,1000,1000,1000,1000,1000,1000,1000''/   
     >  22x,''1000,1000,1000,1000,1000,1000,1000,1000''/   
     >  22x,''1000,1000,1000,1000,1000,1000,1000,1000'')')

      write(10,'(/''hhodo_neg_ped_limit = 1000,'',
     >  ''1000,1000,1000,1000,1000,1000,1000''/
     >  22x,''1000,1000,1000,1000,1000,1000,1000,1000''/   
     >  22x,''1000,1000,1000,1000,1000,1000,1000,1000''/   
     >  22x,''1000,1000,1000,1000,1000,1000,1000,1000''/   
     >  22x,''1000,1000,1000,1000,1000,1000,1000,1000''/   
     >  22x,''1000,1000,1000,1000,1000,1000,1000,1000''/   
     >  22x,''1000,1000,1000,1000,1000,1000,1000,1000''/   
     >  22x,''1000,1000,1000,1000,1000,1000,1000,1000'')')

      close(unit=10)

! Diagnositc information 
      write(fname,'(''HTOFCAL/htofcal'',i5.5,''.diag'')') 
     >  runno
      open(unit=10,file=fname)

      write(10,'(1x,''ifail='',i10,
     >  '' (desired value is 0 if fit worked)'')') ifail

      do i=1,100
       if(nhit(i).gt.0 .or.nhit(100+i).gt.0.) then
        write(10,'(i3,2i6,6f7.1)') i,nhit(i),nhit(100+i),
     >    -toff(i),-toff(100+i),-1./vel(i),-1./vel(100+i),
     >    -quad(i),-quad(100+i) 
       endif
      enddo
      write(10,'(''ADCHIST'')')
      do i=1,200
        if(nhit(i).gt.0) write(10,'(i4,16i4)') 
     >    i,(adchist(i,j)/10,j=1,16)
      enddo
      write(10,'(''THIST'')')
      do i=1,200
        if(nhit(i).gt.0) write(10,'(i4,10i5)') 
     >    i,(thist(i,j)/100,j=1,10)
      enddo
      write(10,'('' PHIST'')')
      do i=1,200
        if(nhit(i).gt.0) write(10,'(i4,16i4)') 
     >    i,(phist(i,j)/20,j=1,16)
      enddo

      close(unit=10)

      return
      end