;pro plot_erf_results, ps = ps
;restore, 'erf_results.sav'
; loads X2, energy,width,pixels,date

; Make a histogram of the mean energy and deviation

; Calculate some statistics:
;print, 'Statistics for the Reconstructed Mean Energy:'
loadct, 39

;tsize =3.5
tsize =1.0

;PRINT, 'Mean: ', result[0] & PRINT, '1-Sigma: ', std & $
;   PRINT, 'FWHM: ', fwhm & $
;   PRINT, 'Skewness: ', result[2] & PRINT, 'Kurtosis: ', result[3] 

;print, noutlier


result = MOMENT(allx2)
std = sqrt(result[1])
fwhm = 2.3582 * std
diff = where(allx2 gt 1.5, noutlier)


bmax=max(allx2)
bmin=0.0
nbin=100.

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

stp2=findgen(256)/256*(bmax-bmin)+bmin
bar2=fltarr(n_elements(stp2),2)
bar2[*,0]=stp2
bar2[*,1]=stp2
x2b=bytscl(allx2,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 = 'X2 of ERF Fit', $
      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, xr1[0]+(xr1[1]-xr1[0])*0.05, yr1[1]*0.7, charsize=csize, $
;	'Outliers (>3sigma): '+string(noutlier, format = '(i2)')
xyouts, 0.1,0.9,charsize=csize,/normal,$
	'Statistics for the ERF Fits:'



if not keyword_set(set_ps_plot) then hak


;;;;;;;;;;;;;;;;;;;; Now do the mean energy ;;;;;;;;;;;;;;;;;;;;

; Identify the pixels that failed this test:
badOnes = where(allx2 eq 0, nbad)
if nbad gt 0 then begin
   energy[badOnes] = average(energy[where(allx2 ne 0)])
endif

result = MOMENT(energy)
std = sqrt(result[1])

diff = where(abs(energy - result[0]) / std gt 3, noutlier)
newenergy = energy
if noutlier gt 0 then begin
   goodones = where(abs(energy - result[0]) / std le 3)
   goodmean = average(energy[goodones])
   for i =0, noutlier -1 do newenergy[diff[i]] = goodmean
   result = moment(newenergy)
   std = sqrt(result[1])
endif


fwhm = 2.3582 * std
;bmax=result[0] + 2*std
;bmin=result[0] - 2*std
bmax = max(newenergy)
bmin = min(newenergy)
nbin=100.

eh = histogram(energy, min = bmin, max = bmax, binsize = (bmax-bmin)/nbin)
eind = findgen(n_elements(eh)) * (bmax-bmin)/nbin + bmin

stp2=findgen(256)/256*(bmax-bmin)+bmin
bar2=fltarr(n_elements(stp2),2)
bar2[*,0]=stp2
bar2[*,1]=stp2
energyb=bytscl(newenergy,min=bmin,max=bmax)
energyb[diff] = max(energyb)

plotimage,energyb,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=result[0] + 1*[-sqrt(result[1]), sqrt(result[1])]
xr1=[bmin,bmax]
yr1=[0, max(eh)*1.5]
plot, eind, eh, psym = 10, position = [0.65,0.18,0.95,0.8], $
      xtitle = 'Reconstructed mean energy (keV)', $
      ytitle = 'Number of pixels', $
      charsize=csize1,/noerase,$
      xrange = xr1, yrange=yr1, /xst, /yst

;oplot, [2.25, 2.25], [0, 100], color = 5
;oplot, [result[0], result[0]], [0, 100]

;xyouts, result[0] + 1 * std, 25, 'Actual (red): 2.25 keV', charsize = 1
xyouts, xr1[0]+(xr1[1]-xr1[0])*0.05,yr1[1]*0.9, charsize=csize, $
	'Mean: '+string(result[0], format = '(f4.2)')+' keV'
xyouts, xr1[0]+(xr1[1]-xr1[0])*0.05,yr1[1]*0.8, charsize=csize, $
	'1-Sigma: '+string(std, format = '(f5.3)')+' keV'
xyouts, xr1[0]+(xr1[1]-xr1[0])*0.05,yr1[1]*0.7, charsize=csize, $
	'Outliers: '+string(noutlier, format = '(i2)')

xyouts, 0.1,0.9,charsize=csize,/normal,$
       'Statistics for the Reconstructed Mean Energy:'

if not keyword_set(set_ps_plot) then hak


;;;;;;;;;;;;;; Now do the 1-Sigma width of the ERF Functions ;;;;;;;;;;;

badOnes = where(allx2 eq 0, nbad)
if nbad gt 0 then begin
   width[badOnes] = max(width[where(allx2 ne 0)])
endif

;print, 'Statistics for the Reconstructed 1-Sigma Width:'
result = MOMENT(width)
std = sqrt(result[1])
diff = where(abs(width - result[0]) / std gt 3, noutlier)

; Find real noisy pixels:
outliers = where(width gt 0.4, noutlier)
newWidth = width
if noutlier gt 0 then begin
   goodones = where(abs(width - result[0]) / std le 3)
   goodmean = average(width[goodones])
   for i = 0, noutlier -1 do newWidth[outliers[i]] = goodmean
   result = moment(newwidth)
   std = sqrt(result[1])
endif


fwhm = 2.3582 * std

bmin=min(newwidth)
bmax=max(newwidth) 
nbin=100.


wh = histogram(width, min = bmin, max = bmax, binsize = (bmax-bmin)/nbin)
wInd = findgen(n_elements(wh)) * (bmax-bmin)/nbin + bmin

stp2=findgen(256)/256*(bmax-bmin)+bmin
bar2=fltarr(n_elements(stp2),2)
bar2[*,0]=stp2
bar2[*,1]=stp2
widthb=bytscl(newwidth,min=bmin,max=bmax)
widthb[diff] = max(widthb)

plotimage,widthb,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 = result[0] + 1 * [-sqrt(result[1]), sqrt(result[1])]
xr1 = [bmin, bmax]
yr1 = [0, max(wh)*1.5]
plot, wind, wh, psym = 10, position = [0.65,0.18,0.95,0.8], $
      xtitle = '1-Sigma Width (keV)', $
      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)')+' keV'
xyouts, xr1[0]+(xr1[1]-xr1[0])*0.05, yr1[1]*0.8, charsize=csize, $
	'1-Sigma: '+string(std, format = '(f4.2)')+' keV'
xyouts, xr1[0]+(xr1[1]-xr1[0])*0.05, yr1[1]*0.7, charsize=csize, $
	'Noisy pizels (>0.4 keV): '+string(noutlier, format = '(i5)')

xyouts, 0.1,0.9,charsize=csize,/normal,$
	'Statistics for the Reconstructed 1-Sigma Width:'





;; Okay, now going to do the "threshold" measurement:

; mask out outliers:

aveE = average(newenergy)
effE = newenergy - aveE
thresh = 5 * newwidth - effE
bmin = min(thresh)
bmax = max(thresh)
nbin = 100
thresh[diff]  = max(thresh)



wh = histogram(thresh, min = bmin, max = bmax, binsize = (bmax-bmin)/nbin)
wInd = findgen(n_elements(wh)) * (bmax-bmin)/nbin + bmin

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

plotimage,widthb,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 + 2., /noerase,$
        xthick=tsize, ythick=tsize, charthick=tsize,$
        imgxrange=[bmin,bmax], title = 'Threshold Energy (keV)'


cumu = reverse(findgen(n_elements(thresh) + 1))
xr1 = [min(thresh), max(thresh)]
yr1 = [0, max(cumu)]
thresh = thresh[sort(thresh)]
plot, thresh, cumu, position = [0.65,0.58,0.95,0.90], $
      xtitle = 'Threshold Energy (keV)', $
      ytitle = 'Number of noisy pixels', $
      charsize=csize1, /noerase, $
      yrange = yr1,xrange = xr1, /xst, /yst, psym = 4

xyouts, xr1[0]+(xr1[1]-xr1[0])*0.25, yr1[1]*1.05, charsize=csize, $
	'Max threshold: '+string(max(thresh), format = '(f4.2)')+' keV'

xr2 = [1.8, max(thresh)]
yr2 = [0, max(cumu[where(thresh gt xr2[0])])]
plot, thresh, cumu, position = [0.65,0.10,0.95,0.42], $
      xtitle = 'Threshold Energy (keV)', $
      ytitle = 'Number of noisy pixels', $
      charsize=csize1, /noerase, $
      yrange = yr2,xrange = xr2, /xst, /yst, psym = 4

num2 = where(thresh ge 2., noisy)
xyouts, xr2[0]+(xr2[1]-xr2[0])*0.25, yr2[1]*1.05, charsize=csize, $
	'Noisy pixels at 2 keV: '+string(noisy, format = '(f6.2)')

;limit = min(where(cumu lt 10))
;oplot, [thresh[limit], thresh[limit]], yr1

;; xyouts, xr1[0]+(xr1[1]-xr1[0])*0.25, yr1[1]*1.05, charsize=csize, $
;; 	'10 pixel threshold: '+string(thresh[limit], format = '(f4.2)')+' keV'



end