;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


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

energy = 100. * (postenergy - preenergy) / (0.5 *(postenergy + preenergy))
result = MOMENT(energy)
std = sqrt(result[1])

diff = where(abs(energy - result[0]) / std gt 3, noutlier)
if noutlier gt 0 then begin
   goodones = where(abs(energy - result[0]) / std le 3)
   newenergy = energy[goodones]
   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(energy,min=bmin,max=bmax)

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 = 'Percent Difference in Reconstructed mean energy', $
      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 = '(f8.2)')+' %'
xyouts, xr1[0]+(xr1[1]-xr1[0])*0.05,yr1[1]*0.8, charsize=csize, $
	'1-Sigma: '+string(std, format = '(f5.3)')+' %'
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 ;;;;;;;;;;;


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

if noutlier gt 0 then begin
   goodones = where(abs(width - result[0]) / std le 3)
   newwidth = width[goodones]
   result = moment(newwidth)
   std = sqrt(result[1])
endif


fwhm = 2.3582 * std

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

; Redefine outliers as those with width > 0.4 keV

outliers = where(width gt 0.4, noutlier)

;print, noutlier

;PRINT, 'Mean: ', result[0] & PRINT, '1-Sigma: ', sqrt(result[1]) & $
;   PRINT, 'FWHM: ', 2.3582 * sqrt(result[1]) & $
;   PRINT, 'Skewness: ', result[2] & PRINT, 'Kurtosis: ', result[3] 


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(width,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, /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 = 'Percent Difference in 1-Sigma Width', $
      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 = '(f8.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,$
	'Statistics for the Reconstructed 1-Sigma Width:'

end