;pro fit_offset

;restore, 'offset_data.sav'
; loads val, which is [frame, col, row] with 
; frame = 0-4, row and col are both 0 - 31
; Stored data are the digitized Vout.

; Convert val to Voltage
range = 10.d
nsteps = 4095.d
res = range / nsteps
volts = dc_offset * res


; volts[0, *, *] is Vout with Vin set to 0 V
; volts[1, *, *] is Vout with Vin set to 2 V
; volts[2, *, *] is Preamp Output Volt (should be uniform and ~ 0.1 V)
; volts[3, *, *] is leakage (should be uniform and ~ 1.6 V)
; volts[4, *, *] is DC level at shaper (?)

; leakage is 150e6 ohms.
; I = V / R

; First thing is to generate a Vout -> Vin
; map, assuming linearity:
; Vin = m * Vout + offset
; 2 = m[*, *] * (volts[1, *, *] - volts[0, *, *])
; m[*, *] = 2 / (volts[1, *, *] - volts[0, *, *])
; offset = - m[*, *] * volts[0, *, *]

slope = 2.d / (reform(volts[1, *, *]) - reform(volts[0, *, *]))
offset = -slope * reform(volts[0, *, *])

;slope = 2.d / (reform(pixeldc[1, *, *]) - reform(pixeldc[0, *, *]))
;offset = -slope * reform(pixeldc[0, *, *])

;print, slope
;stop


for i = 2, 4 do volts[i, *, *] = volts[i, *, *] * slope + offset

zerolevel = reform(volts[0, *, *])
raillevel = reform(volts[1, *, *])
preAmp = reform(volts[2, *, *])
leakage = reform(volts[3, *, *]) / 150e6
dclevel = reform(volts[4, *, *])

; Now plotting for the reports








;;;;;;;;;;;; Zero-level image ;;;;;;;;;;;;;;;;;;;

result = moment(zerolevel)
mean = result[0]
std = sqrt(result[1])
bmax=mean+3*std 
bmin= (mean-3*std)
nbin=100.

x2h = histogram(zerolevel, min = bmin, max = bmax, binsize = (bmax-bmin)/nbin)
x2ind = findgen(n_elements(x2h)) * (bmax-bmin)/nbin + bmin

stp2=findgen(256)/256*(bmax-bmin)+bmin
bar2=fltarr(n_elements(stp2),2)
bar2[*,0]=stp2
bar2[*,1]=stp2

x2b=bytscl(zerolevel,min=bmin,max=bmax)

plotimage,x2b,yr=[32,0],/xst,/yst,$
	pixel_aspect_ratio=1.0,position=[0.14,0.12,0.48,0.75],$
	charsize=csize1,xtit="Column",ytit="Row"

plotimage,bytscl(bar2),position=[0.14,0.77,0.48,0.80],$
        yticks=1,ytickname=[' ',' '],charsize=csize1, /noerase,$
        xthick=tsize, ythick=tsize, charthick=tsize,$
        imgxrange=[bmin,bmax]

xr1 = [bmin, bmax]
yr1 = [0, max(x2h)*1.5]
plot, x2ind, x2h, psym = 10, position = [0.65,0.18,0.95,0.8], $
      xtitle = 'Vout', $
      ytitle = 'Number of pixels', $
      charsize=csize1, /noerase, $
      xrange = xr1, yrange=yr1, /xst, /yst

xyouts, xr1[0]+(xr1[1]-xr1[0])*0.05, yr1[1]*0.9, charsize=csize, $
	'Mean: '+string(result[0], format = '(f4.2)')+' V'
xyouts, xr1[0]+(xr1[1]-xr1[0])*0.05, yr1[1]*0.8, charsize=csize, $
	'1-Sigma: '+string(std, format = '(f4.2)')+' V'

xyouts, 0.1,0.9,charsize=csize,/normal,$
	'Readout Voltage for 0 Volts On-Pixel:'



if not keyword_set(set_ps_plot) then hak

;;;;;;;;;;;; 2V Rail-level image ;;;;;;;;;;;;;;;;;;;

result = moment(raillevel)
mean = result[0]
std = sqrt(result[1])
bmax=mean+3*std
bmin=mean-3*std
nbin=100.

x2h = histogram(raillevel, min = bmin, max = bmax, binsize = (bmax-bmin)/nbin)
x2ind = findgen(n_elements(x2h)) * (bmax-bmin)/nbin + bmin

stp2=findgen(256)/256*(bmax-bmin)+bmin
bar2=fltarr(n_elements(stp2),2)
bar2[*,0]=stp2
bar2[*,1]=stp2

x2b=bytscl(raillevel,min=bmin,max=bmax)

plotimage,x2b,yr=[32,0],/xst,/yst,$
	pixel_aspect_ratio=1.0,position=[0.14,0.12,0.48,0.75],$
	charsize=csize1,xtit="Column",ytit="Row"

plotimage,bytscl(bar2),position=[0.14,0.77,0.48,0.80],$
        yticks=1,ytickname=[' ',' '],charsize=csize1, /noerase,$
        xthick=tsize, ythick=tsize, charthick=tsize,$
        imgxrange=[bmin,bmax]

xr1 = [bmin, bmax]
yr1 = [0, max(x2h)*1.5]
plot, x2ind, x2h, psym = 10, position = [0.65,0.18,0.95,0.8], $
      xtitle = 'Vout', $
      ytitle = 'Number of pixels', $
      charsize=csize1, /noerase, $
      xrange = xr1, yrange=yr1, /xst, /yst

xyouts, xr1[0]+(xr1[1]-xr1[0])*0.05, yr1[1]*0.9, charsize=csize, $
	'Mean: '+string(result[0], format = '(f4.2)')+' V'
xyouts, xr1[0]+(xr1[1]-xr1[0])*0.05, yr1[1]*0.8, charsize=csize, $
	'1-Sigma: '+string(std, format = '(f4.2)')+' V'

xyouts, 0.1,0.9,charsize=csize,/normal,$
	'Readout Voltage for 2 Volts On-Pixel:'

if not keyword_set(set_ps_plot) then hak

;;;;;;;;;;;; Slope for conversion to On-Pixel Voltage ;;;;;;;;;;;;;;;;;;;

result = moment(slope)
mean = result[0]
std = sqrt(result[1])
bmax=mean+3*std
bmin=mean-3*std
nbin=100.

x2h = histogram(slope, min = bmin, max = bmax, binsize = (bmax-bmin)/nbin)
x2ind = findgen(n_elements(x2h)) * (bmax-bmin)/nbin + bmin

stp2=findgen(256)/256*(bmax-bmin)+bmin
bar2=fltarr(n_elements(stp2),2)
bar2[*,0]=stp2
bar2[*,1]=stp2

x2b=bytscl(slope,min=bmin,max=bmax)

plotimage,x2b,yr=[32,0],/xst,/yst,$
	pixel_aspect_ratio=1.0,position=[0.14,0.12,0.48,0.75],$
	charsize=csize1,xtit="Column",ytit="Row"

plotimage,bytscl(bar2),position=[0.14,0.77,0.48,0.80],$
        yticks=1,ytickname=[' ',' '],charsize=csize1, /noerase,$
        xthick=tsize, ythick=tsize, charthick=tsize,$
        imgxrange=[bmin,bmax]

xr1 = [bmin, bmax]
yr1 = [0, max(x2h)*1.5]
plot, x2ind, x2h, psym = 10, position = [0.65,0.18,0.95,0.8], $
      xtitle = 'Slope', $
      ytitle = 'Number of pixels', $
      charsize=csize1, /noerase, $
      xrange = xr1, yrange=yr1, /xst, /yst

xyouts, xr1[0]+(xr1[1]-xr1[0])*0.05, yr1[1]*0.9, charsize=csize, $
	'Mean: '+string(result[0], format = '(f4.2)')
xyouts, xr1[0]+(xr1[1]-xr1[0])*0.05, yr1[1]*0.8, charsize=csize, $
	'1-Sigma: '+string(std, format = '(f4.2)')

xyouts, 0.1,0.9,charsize=csize,/normal,$
	'Readout to On-pixel Scale Factor:'

if not keyword_set(set_ps_plot) then hak

;;;;;;;;;;;; Pre-amp image ;;;;;;;;;;;;;;;;;;;

result = moment(preamp)
mean = result[0]
std = sqrt(result[1])
bmax=mean+3*std
bmin=0
nbin=100.

outliers = where( preamp lt 20e-3, noutliers)

x2h = histogram(preamp, min = bmin, max = bmax, binsize = (bmax-bmin)/nbin)
x2ind = findgen(n_elements(x2h)) * (bmax-bmin)/nbin + bmin

stp2=findgen(256)/256*(bmax-bmin)+bmin
bar2=fltarr(n_elements(stp2),2)
bar2[*,0]=stp2
bar2[*,1]=stp2

x2b= 255B - bytscl(preamp,min=bmin,max=bmax)
plotimage,x2b,yr=[32,0],/xst,/yst,$
	pixel_aspect_ratio=1.0,position=[0.14,0.12,0.48,0.75],$
	charsize=csize1,xtit="Column",ytit="Row"

plotimage, 255B - bytscl(bar2),position=[0.14,0.77,0.48,0.80],$
        yticks=1,ytickname=[' ',' '],charsize=csize1, /noerase,$
        xthick=tsize, ythick=tsize, charthick=tsize,$
        imgxrange=[bmin,bmax]

xr1 = [0, bmax]
yr1 = [0, max(x2h)*1.5]
plot, x2ind, x2h, psym = 10, position = [0.65,0.18,0.95,0.8], $
      xtitle = 'Vout', $
      ytitle = 'Number of pixels', $
      charsize=csize1, /noerase, $
      xrange = xr1, yrange=yr1, /xst, /yst

xyouts, xr1[0]+(xr1[1]-xr1[0])*0.05, yr1[1]*0.9, charsize=csize, $
	'Mean: '+string(result[0], format = '(f4.2)')+' V'
xyouts, xr1[0]+(xr1[1]-xr1[0])*0.05, yr1[1]*0.8, charsize=csize, $
	'1-Sigma: '+string(std, format = '(f4.2)')+' V'
xyouts, xr1[0]+(xr1[1]-xr1[0])*0.05, yr1[1]*0.7, charsize=csize, $
	'Low-level (< 20 mV): '+string(noutliers, format = '(i2)')

xyouts, 0.1,0.9,charsize=csize,/normal,$
	'Pre-Amp Output:'

if not keyword_set(set_ps_plot) then hak

;;;;;;;;;;;; Leakage image ;;;;;;;;;;;;;;;;;;;

leakage *= 1e9
result = moment(leakage)
mean = result[0]
std = sqrt(result[1])
bmax=mean+3*std
bmin=mean-3*std
nbin=100.

outliers = where( (leakage - result[0]) / std gt 3, noutliers)

x2h = histogram(leakage, min = bmin, max = bmax, binsize = (bmax-bmin)/nbin)
x2ind = findgen(n_elements(x2h)) * (bmax-bmin)/nbin + bmin

stp2=findgen(256)/256*(bmax-bmin)+bmin
bar2=fltarr(n_elements(stp2),2)
bar2[*,0]=stp2
bar2[*,1]=stp2

x2b=bytscl(leakage,min=bmin,max=bmax)

plotimage,x2b,yr=[32,0],/xst,/yst,$
	pixel_aspect_ratio=1.0,position=[0.14,0.12,0.48,0.75],$
	charsize=csize1,xtit="Column",ytit="Row"

plotimage,bytscl(bar2),position=[0.14,0.77,0.48,0.80],$
        yticks=1,ytickname=[' ',' '],charsize=csize1, /noerase,$
        xthick=tsize, ythick=tsize, charthick=tsize,$
        imgxrange=[bmin,bmax], title = 'nA'

xr1 = [bmin, bmax]
yr1 = [0, max(x2h)*1.5]
plot, x2ind, x2h, psym = 10, position = [0.65,0.18,0.95,0.8], $
      xtitle = 'nA', $
      ytitle = 'Number of pixels', $
      charsize=csize1, /noerase, $
      xrange = xr1, yrange=yr1, /xst, /yst

xyouts, xr1[0]+(xr1[1]-xr1[0])*0.05, yr1[1]*0.9, charsize=csize, $
	'Mean: '+string(result[0], format = '(f6.2)')+' nA'
xyouts, xr1[0]+(xr1[1]-xr1[0])*0.05, yr1[1]*0.8, charsize=csize, $
	'1-Sigma: '+string(std, format = '(f6.2)')+' nA'
xyouts, xr1[0]+(xr1[1]-xr1[0])*0.05, yr1[1]*0.7, charsize=csize, $
	'Outliers (>3sigma): '+string(noutliers, format = '(i2)')

xyouts, 0.1,0.9,charsize=csize,/normal,$
	'Leakage Baseline for Normal Mode:'

if not keyword_set(set_ps_plot) then hak


;;;;;;;;;;;; DC-Level image ;;;;;;;;;;;;;;;;;;;

result = moment(dclevel)
mean = result[0]
std = sqrt(result[1])
bmax=mean+3*std
bmin=mean-3*std
nbin=100.

outliers = where( (dclevel - result[0]) / std gt 3, noutliers)

x2h = histogram(dclevel, min = bmin, max = bmax, binsize = (bmax-bmin)/nbin)
x2ind = findgen(n_elements(x2h)) * (bmax-bmin)/nbin + bmin

stp2=findgen(256)/256*(bmax-bmin)+bmin
bar2=fltarr(n_elements(stp2),2)
bar2[*,0]=stp2
bar2[*,1]=stp2

x2b=bytscl(dclevel,min=bmin,max=bmax)

plotimage,x2b,yr=[32,0],/xst,/yst,$
	pixel_aspect_ratio=1.0,position=[0.14,0.12,0.48,0.75],$
	charsize=csize1,xtit="Column",ytit="Row"

plotimage,bytscl(bar2),position=[0.14,0.77,0.48,0.80],$
        yticks=1,ytickname=[' ',' '],charsize=csize1, /noerase,$
        xthick=tsize, ythick=tsize, charthick=tsize,$
        imgxrange=[bmin,bmax]

xr1 = [bmin, bmax]
yr1 = [0, max(x2h)*1.5]
plot, x2ind, x2h, psym = 10, position = [0.65,0.18,0.95,0.8], $
      xtitle = 'Vout', $
      ytitle = 'Number of pixels', $
      charsize=csize1, /noerase, $
      xrange = xr1, yrange=yr1, /xst, /yst

xyouts, xr1[0]+(xr1[1]-xr1[0])*0.05, yr1[1]*0.9, charsize=csize, $
	'Mean: '+string(result[0], format = '(f4.2)')+' V'
xyouts, xr1[0]+(xr1[1]-xr1[0])*0.05, yr1[1]*0.8, charsize=csize, $
	'1-Sigma: '+string(std, format = '(f4.2)')+' V'
xyouts, xr1[0]+(xr1[1]-xr1[0])*0.05, yr1[1]*0.7, charsize=csize, $
	'Outliers (>3sigma): '+string(noutliers, format = '(i2)')

xyouts, 0.1,0.9,charsize=csize,/normal,$
	'DC Level Data:'

;save, slope, offset, zerolevel, raillevel, leakage, leakage, dclevel, file = 'offset_volts.sav'

end