;pro fit_erf


;restore, 'LOG090209.sav'
;restore, outfile
; loads date, pixels, and counts

;npix = n_elements(counts[0, *])
npix=32*32

!p.multi = [0, 8, 8]
;!y.range = [-0.3, 0.9]
!y.style = 1
!p.charsize=0.8
!p.symsize=0.5

; Just pix off the first one, for now.
;X2 = fltarr(npix)
allfit = fltarr(2, 32, 32)
allx2 = fltarr(32, 32)
ind = findgen(16)
;for i = 0, npix - 1 do begin
for i=0,15 do begin
   for y0=0,7 do begin
       for x0=0,7 do begin
	   x1=x0+((i mod 4)*8)
	   y1=y0+(i/4)*8

           thiscounts = reverse(reform(counts2[x1,y1,*]))
           statErr = sqrt(thiscounts++)
           tooSmall = where(statErr eq 0, nsmall)
           if nsmall gt 0 then statErr[tooSmall]++
           thiscounts /= 200.
           statErr /= 200.
 
          plot, ind, thiscounts,psym=3,$
  	     title=strcompress('X '+string(x1)+', Y '+string(y1)),$
	     yrange=[0,1],xtitl='Thresholds',ytit='Normalized counts'
;           oploterr, ind, thiscounts, statErr, psym = 3

           start = [8.d, 5.d]
           test = myerf(ind, start)
           
           result = MPFITFUN('MYERF', ind, thiscounts, $
               statErr, start, yfit = yfit, bestnorm = best, dof = dof, /quiet)
           oplot,ind,yfit
           ;X2[i] = best / dof
           allfit[*,x1,y1] = result
           allx2[x1, y1] = best / dof
           ;print,i,yfit
           ;print,x1,y1
     endfor
   endfor
endfor
;if not keyword_set(set_ps_plot) then hak

; allfit[0, *] is the mean value
; allfit[1, *] is the deviation (square root of the variance)

; Conversion factors (via Rick Cook):
; The test scans the threshold setting from 20 through 50 in steps of 2.
; The test pulse amplitude was fixed at nominally 2.25 keV.
; The rough scale factor for threshold setting is 2.25 keV/35 steps.

; So the mean value conversion goes: (2.25 / 35) * (mean + 20)
; while the variance conversion goes (2.25 / 35) * var
; to put both in keV.

; Now conversions:
; setting N
; V = (N / 2^16) * 5
; Q = V * C, C = 20e-15 farads
; E = Q / Q_1kev, Q_1keV = 200 * 1.7e-19
; E = V * C / Q_1keV
; Channel-to-energy conversion:
; E = channel * (2^-16 * 5 * 220e-15) / (200. * 1.6e-19)

dE =  2d^(-16) * 5 * 20e-15 / (200. * 1.6e-19)
energy = (25. + 5.*reform(allfit[0, *, *]) ) *dE
width = 5d*allfit[1, *, *] * dE

;save, X2, energy,width,pixels,date, file = 'erf_results.sav'

end