PRO hsi_snoop_anal, COMBINE_ALL=combine_all, QED_WRITE=qed_write
;
; Having created rdb files from the HSI exposures:
;  - combine all the individual rdb files into one, include
;    numerical values of the runids and cumulative values
;    of the counts. (This is only done once by setting
;    keyword COMBINE_ALL.)
;
;  - Make useful plots from the combined data, etc.
;

; where to get the rdb output files
keep_em_dir = !DDHETGCAL+'/cmp/hsi_snoop'

print, ''
print, ' HSI Snooper Analysis'
print, ''

if KEYWORD_SET(COMBINE_ALL) then begin
; - - - - - - - - - - -
; Loop through the date codes
for iy = 96,97 do begin
  mstart = 12
  mstop = 12
  if iy EQ 97 then begin
    mstart = 1
    mstop = 2
  end
  for im = mstart, mstop do begin
    for id=1, 31 do begin
      month_str = STRCOMPRESS(im,/REMOVE_ALL)
      if STRLEN(month_str) EQ 1 then month_str = '0'+month_str
      day_str = STRCOMPRESS(id,/REMOVE_ALL)
      if STRLEN(day_str) EQ 1 then day_str = '0'+day_str
      date_code = STRCOMPRESS(iy,/REMOVE_ALL) + month_str + day_str

      rdb_file = findfile(keep_em_dir+ $
		'/hsi_snoop_'+date_code+'.rdb', COUNT=n_found)
      if n_found EQ 1 then begin

        print, ' Getting ['+date_code+']'
        data = rdb_read(rdb_file(0))

        ; Create numerical run id values
        nids = FLTARR(n_elements(data))
        temp=0.0
        for it=0,n_elements(data)-1 do begin
          READS, data(it).run_id, temp
          nids(it) = temp
        end

        ; Assemble all the data into arrays
        if n_elements(all_nids) EQ 0 then begin
          all_date = data.mst_date
          all_run_id = data.run_id
          all_iter = data.iteration
          all_total = data.c_total
          all_tile = data.c_tile
          all_depress = data.c_depress
          all_nids = nids
        end else begin
          all_date = [all_date, data.mst_date]
          all_run_id = [all_run_id, data.run_id]
          all_iter = [all_iter, data.iteration]
          all_total = [all_total, data.c_total]
          all_tile = [all_tile, data.c_tile]
          all_depress = [all_depress, data.c_depress]
          all_nids = [all_nids,nids]
        end

      end ; found

    end ; id
  end ; im
end ; iy

; create the cumualtive counts
cum_total = all_total
cum_tile = all_tile
cum_depress = all_depress
for id=1,n_elements(all_depress)-1 do begin
  cum_depress(id) = cum_depress(id-1) + cum_depress(id)
  cum_total(id) = cum_total(id-1) + cum_total(id)
  cum_tile(id) = cum_tile(id-1) + cum_tile(id)
end

; write it all out to one rdb file; include a column for
; the QE depression correction value for each file as well.
out = REPLICATE({mst_date:'', run_id:'', iteration:'', $
	c_total:0.0, c_tile:0.0, c_depress:0.0, $
	nid:0.0, $
	cum_total:0.0, cum_tile:0.0, cum_depress:0.0, $
	qed_val:1.0}, n_elements(all_nids))

out.mst_date = all_date
out.run_id = all_run_id
out.iteration = all_iter
out.c_total = all_total
out.c_tile = all_tile
out.c_depress = all_depress
out.nid = all_nids
out.cum_depress = cum_depress
out.cum_tile = cum_tile
out.cum_total = cum_total

rdb_write, keep_em_dir+'/hsi_snoop_all.rdb', out

; - - - - - - - - - - -
end ; of COMBINE_ALL


; Read in the combined data file and make some plots...

sn = rdb_read(keep_em_dir+'/hsi_snoop_all.rdb')

; Show cumulative counts for the three different regions...

pre_print_portrait
!p.multi=[0,1,3]

plot, sn.nid - 1.E5, sn.cum_total, $
	CHARSIZE=1.6, $
	TITLE = 'TOTAL HSI Counts Detected vs Run ID', $
	XRANGE=[6400.,12000.], XSTYLE=1, $
        XTITLE = 'Runid - 100000', $
	YTITLE='Cummulative Counts'
plot, sn.nid - 1.E5, sn.cum_tile, $
	CHARSIZE=1.6, $
	TITLE = 'Central Tile Counts Detected vs Run ID', $
	XRANGE=[6400.,12000.], XSTYLE=1, $
        XTITLE = 'Runid - 100000', $
	YTITLE='Cummulative Counts'
plot, sn.nid - 1.E5, sn.cum_depress, $
	CHARSIZE=1.6, $
	TITLE = 'Depressed Region Counts Detected vs Run ID', $
	XRANGE=[6400.,12000.], XSTYLE=1, $
        XTITLE = 'Runid - 100000', $
	YTITLE='Cummulative Counts'

device, /close
SPAWN, 'cp idl.ps '+keep_em_dir+'/hsi_snoop_cumul.ps'

set_plot, 'X'


; Get the measured QE values ("table" keywords)
dummy = hsi_qe_depress(sn(0).nid, $
	TABLE_IDS=table_ids, TABLE_QEDS=table_qeds)

; Find the cum counts for each of the tabulated
; QED measurements (skipping first and last)
table_counts = table_ids  ; define the array
table_counts(0) = 0.0
last = n_elements(table_counts)-1
table_counts(last) = max(sn.cum_depress)
for it = 1, last-1 do begin
  these_counts = where(sn.nid EQ table_ids(it))
  table_counts(it) = sn(these_counts(0)).cum_depress
end

; Fit the data to a quadratic
fit_qeds = table_qeds(2:5)
fit_counts = table_counts(2:5)

coefs = poly_fit(fit_counts, fit_qeds, 2)

calc_counts = findgen(100)*40000.0
calc_qed = coefs(0) + coefs(1)*calc_counts + $
	coefs(2)*calc_counts*calc_counts

; Evaluate the QE depression for each file using the counts to qed
; fit; use the cumulative counts minus 1/2 of this measurements
; counts to get average QED for the measurement...
ave_counts = sn.cum_depress  - sn.c_depress/2.0
sn.qed_val = coefs(0) + coefs(1)*ave_counts + $
	coefs(2)*ave_counts*ave_counts

; Make these plots
pre_print_portrait
!p.multi = [0,1,4]

; Cumulative counts plot
plot, sn.nid - 1.E5, sn.cum_depress, $
	CHARSIZE=1.6, $
	TITLE = 'HSI Depressed Region Counts Detected vs Run ID', $
	XRANGE=[6400.,12000.], XSTYLE=1, $
        XTITLE = 'Runid - 100000', $
	YTITLE='Cummulative Counts'

; "Measured" QE depression values vs Run ID
plot, table_ids - 1.E5, table_qeds, PSYM=2, $
	CHARSIZE=1.6, $
	TITLE = 'Measured QE Depression value vs Run ID', $
	XRANGE=[6400.,12000.], XSTYLE=1, $
        XTITLE = 'Runid - 100000', $
	YTITLE='HSI QE Depression', $
	YRANGE=[0.65, 1.05], YSTYLE=1

; "Measured" QE depression values vs Cum Counts
plot, table_counts, table_qeds, PSYM=2, $
	CHARSIZE=1.6, $
	TITLE = 'Measured QE Depression value vs Cummulative Depressed counts', $
	YTITLE='HSI QE Depression', $
	XTITLE = 'Cummulative Counts', $
	YRANGE=[0.65, 1.05], YSTYLE=1
oplot, calc_counts, calc_qed, LINESTYLE=2

; "Predicted" QE depression values vs Run ID
plot, sn.nid - 1.E5, sn.qed_val, PSYM=3, $
	CHARSIZE=1.6, $
	TITLE = 'Predicted QE Depression value vs Run ID', $
	XRANGE=[6400.,12000.], XSTYLE=1, $
        XTITLE = 'Runid - 100000', $
	YTITLE='HSI QE Depression', $
	YRANGE=[0.65, 1.05], YSTYLE=1
oplot, table_ids - 1.E5, table_qeds, PSYM=2

device, /close
SPAWN, 'cp idl.ps '+keep_em_dir+'/hsi_snoop_model.ps'

set_plot, 'X'

; Write out the file with the QED values
if KEYWORD_SET(QED_WRITE) then begin
  rdb_write, keep_em_dir+'/hsi_snoop_all.rdb', sn
end

; "Prune" these values for a smaller set needed to get
; qed value accurate to 1% from runid lookup table:
xs_needed = interp_pruned_x(sn.nid, [[sn.qed_val]], [0.01])

; Print these values so that they can be cut and pasted into
; the hsi_qe_depress routine lookup table:
print, ''
for il=0,n_elements(xs_needed)-1 do print, sn(xs_needed(il)).nid, ', $'
print, ''
for il=0,n_elements(xs_needed)-1 do print, sn(xs_needed(il)).qed_val, ', $'

RETURN
END