subroutine hec_cub_evt
c
c for each event, performs a simple cubic fit to find the 
c max adc and time of multiple sample data
c
      implicit none
c
      include'hec_par.inc'      !shared parameters
      include'hec_adc.inc'      !adc values
      include'hec_cub.inc'      !max adc and time for simple cubic fit
      include'hec_datacard.inc' !datacard values
      include'hec_geo.inc'      !geometry correspondence tables
      include'hec_ped.inc'      !pedestals
c
c set inverse of time powers matrix
c
      real tinv(4,4)
c
      integer ic
      integer i_t
      integer j_adc_max, imxtim
      real adcmax, tmax
      real tmax_prime
      integer i_adctmp
      real tzero
      real adc(4), coef(4)
      integer i, j
      real sq
c
      data tinv(1,1)/1./
      data tinv(1,2)/0./
      data tinv(1,3)/0./
      data tinv(1,4)/0./
      data tinv(2,1)/-1.83333333/
      data tinv(2,2)/3./
      data tinv(2,3)/-1.5/
      data tinv(2,4)/0.33333333/
      data tinv(3,1)/1./
      data tinv(3,2)/-2.5/
      data tinv(3,3)/2./
      data tinv(3,4)/-0.5/      
      data tinv(4,1)/-0.16666666/
      data tinv(4,2)/0.5/
      data tinv(4,3)/-0.5/
      data tinv(4,4)/0.16666666/
c
c reset the array for max adc and time for simple cubic fit
c
      call vzero(adc_cub_ic, 2*i_adc_dim)
c
c if 1996 (which is silly) just output the adc and set time
c to 1 (first and only time slice)
c
      if (irunpd .le. 3) then
        do ic = 1, i_adc_dim
          adc_cub_ic(ic, 1) = real(i_adc_ic(ic))
        enddo
        return
      endif
c
c loop over adcs
c
      do ic = 1, i_adc_used
c     
c find time slice with maximum adc, only considering the selected
c time slice window
c
        j_adc_max = 0.
        imxtim = 1
        do i_t = t1_cub, t2_cub
          if (i_adc_ic_t(ic, i_t) .gt. j_adc_max) then
            j_adc_max = i_adc_ic_t(ic, i_t)
            imxtim = i_t
          endif
        enddo 
c     
c check boundary conditions
c     
        if (imxtim .le. 2 .or. 
     +  imxtim .ge. (i_tim_used - 1)) then
          adcmax = real(i_adc_ic_t(ic, imxtim))
          tmax = real(imxtim)
c     
c if not near boundary look for points surrounding max
c
        else
c
c first find tzero
c
          i_adctmp = max(i_adc_ic_t(ic, imxtim - 2),
     +    i_adc_ic_t(ic, imxtim + 2))
          if (i_adctmp .eq. i_adc_ic_t(ic, imxtim - 2)) then
            tzero = real(imxtim - 2)
          else
            tzero = real(imxtim - 1)
          endif
c
c now get adc for four points used in fit window
c            
          adc(1) = real(i_adc_ic_t(ic, tzero))
          adc(2) = real(i_adc_ic_t(ic, tzero + 1))
          adc(3) = real(i_adc_ic_t(ic, tzero + 2))
          adc(4) = real(i_adc_ic_t(ic, tzero + 3))
c
c now do fit
c
          call vzero(coef, 4)
          do i = 1, 4
            do j = 1, 4
              coef(i) = coef(i) + tinv(i, j)*adc(j)
            enddo
          enddo
c
c use coefficients and cubic polynomial to get tmax and adcmax
c for bad fit, set the fit results to the time slice with 
c maximum adc count in fit window
c
          sq = coef(3)**2 - 3.*coef(2)*coef(4)
          if (sq .lt. 0. .or. coef(4) .eq. 0.) then
            adcmax = real(i_adc_ic_t(ic, imxtim))
            tmax = real(imxtim)
          else
            tmax_prime = (-1.*coef(3) - sqrt(sq))/(3.*coef(4))
c
c make sure the time of maximum is in the fit window
c
            if (tmax_prime .ge. 0. .and. tmax_prime .le. 3) then
              tmax = tmax_prime + tzero
              adcmax = coef(1) + tmax_prime*(coef(2) + 
     +                 tmax_prime*(coef(3) + tmax_prime*coef(4))) 
            else
              adcmax = real(i_adc_ic_t(ic, imxtim))
              tmax = real(imxtim)
            endif
          endif
        endif
c
c fill common for this channel
c
        if (sc_cub .lt. 0. .or. 
     +      adcmax .ge. (adc_ped(ic) + sc_cub*adc_rms(ic))) then
          adc_cub_ic(ic, 1) = adcmax
          adc_cub_ic(ic, 2) = tmax
        else
          adc_cub_ic(ic, 1) = real(i_adc_ic_t(ic, tm_cub))
          adc_cub_ic(ic, 2) = real(tm_cub)
        endif
      enddo
c
      end