asic='H26'
;---------------------------------------------------
loadct,39
datadir='~/nustar/asic_screening/data/protohybrid/'+asic+'PRE/'
reportdir='~/nustar/asic_screening/reports/compare_pre-post/'+asic+'/'
; Make the directory, if necessary
file_mkdir, reportdir

sfx_ndisc='DISC_'+asic

csize =1.0
csize1 = 1.0

;-----------------   NDISC Test  ------------------
   
filename=strcompress(datadir+sfx_ndisc+'PRE.TXT',/remove)
outps1 =strcompress(reportdir+asic+'_thresh_comp.ps',/remove)
outpdf1 =strcompress(reportdir+asic+'_thresh_comp.pdf',/remove)

.run setps

print,filename
.run parse_ndisc
.run fit_erf_comp
.run plot_erf_comp

prethresh = thresh


datadir='~/nustar/asic_screening/data/protohybrid/'+asic+'POST/'
filename=strcompress(datadir+sfx_ndisc+'POST.TXT',/remove)
print,filename
.run parse_ndisc
.run fit_erf_comp
.run plot_erf_comp



   diff = 100. * (thresh - prethresh) / (0.5 * (thresh + prethresh))

   bmin = min(diff)
   bmax = max(diff)
   nbin = 100
   
   dh = histogram(diff, min = bmin, max = bmax, binsize = (bmax - bmin) / nbin)
   dInd = findgen(n_elements(dh))* (bmax - bmin) / nbin + bmin

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

   
   plotimage,diffB,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 = 'Percent difference in Threshold'


   xr1 = [bmin, bmax]
   yr1 = [0, max(dh)*1.5]
   result = moment(diff)
   plot, dind, dh, psym = 10, position = [0.65,0.18,0.95,0.8], $
   xtitle = 'Percent difference in Threshold', $
   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)')


.run endps