PRO hsi_flats, SAVE_CENTRAL=save_central, $
	TILE_PLOTS=tile_plots, ZY_FLIP=zy_flip
;
; Examine the post-XRCF HSI flat-field data sets...
;
; /SAVE_CENTRAL Creates idl save files of the hsi structure
;               consisting of only the central (X,Y in 2000-2100) tile
;               to reduce processing burden of 200 Mb files!
;
; /TILE_PLOTS Creates various .gif output plots...
;
; /ZY_FLIP Switches the HSI-Y and HSI-Z values so that hsi_tile_plots
;          will create Z projections instead of Y projections...
;
; 5/1,3/99 dd
;
@df_common

;------- The data ... ----
; Data files are from SAO : /data/ils/evans_i/hsiflat/events
; Stored locally at MIT : /nfs/spectra/d2/hsiflat
flat_loc = '/nfs/spectra/d2/hsiflat'

; Array of the file names:
flat_files = ['hsi300318.fits', 'hsi300340.fits', 'hsi300348.fits', $
		'hsi300357.fits', 'hsi300363.fits', 'hsi300368.fits', $
		'hsi300395.fits', 'hsi300400.fits', 'hsi300407.fits' ]

flat_targs = ['Fe-La', 'Fe-Ka', 'Mg-Ka', $
		'Al-Ka', 'Sn-La', 'C-Ka', $
		'Cu-La', 'Ti-Ka', 'Mo-La' ]

flat_energies = [ 0.704, 6.403, 1.254, $
		1.486, 3.444, 0.277, $
		0.928, 4.510, 2.293 ]

; List the files...
for iff=0,n_elements(flat_files)-1 do begin
  print, flat_files(iff), ' : ', flat_targs(iff) ,' ,  E = ',flat_energies(iff)
end


if KEYWORD_SET(SAVE_CENTRAL) then begin
for iff=0,n_elements(flat_files)-1 do begin
  ;------- Go through the files..
  ; Select a file
  flat_file = flat_files(iff)

  ; read in the data... (38 M events !?!)
  flat = mrdfits(flat_loc+'/'+flat_file, 1)

  ; Take a peek at the file
  sel = lindgen(100000)
  pre_print_sqr
  plot, flat(sel).Y, flat(sel).Z, PSYM=3, $
	XRANGE=[500.,3500.], XSTYLE=1, $
	YRANGE=[500.,3500.], YSTYLE=1, $
	XTITLE='HSI Y', YTITLE='HSI Z', $
	TITLE='HSI Flatfield Data,  '+flat_file+',  '+flat_targs(iff)+', '+$
		STRCOMPRESS(STRING(flat_energies(iff),FORMAT='(F8.3)'))+' keV'
  device, /close
  SPAWN, 'cp idl.ps '+flat_loc+'/hsi_flat_'+flat_targs(iff)+'.ps'

  ; Select the central tile events:
  ; Define the central tile:
  ytmin = 1920.
  ytmax = 2180.
  ztmin = 1920.
  ztmax = 2180.
  ;
  sel = where(flat.Y GT ytmin)
  flat = flat(sel)
  sel = where(flat.Y LT ytmax)
  flat = flat(sel)
  sel = where(flat.Z GT ztmin)
  flat = flat(sel)
  sel = where(flat.Z LT ztmax)
  flat = flat(sel)

  ; Now save these sub_events to an IDL save file
  save, flat, FILENAME = flat_loc+'/hsi_flat_'+flat_targs(iff)+'.idlsav'
end
end


if KEYWORD_SET(TILE_PLOTS) then begin
for iff=0,n_elements(flat_files)-1 do begin

  ; Get the HSI event structure for the central tile, "flat":
  restore, FILENAME = flat_loc+'/hsi_flat_'+flat_targs(iff)+'.idlsav'

  ; Switch the Y and Z axes if ZY_FLIP is set...
  flip_str = ''
  ; hsi_tile_plots current default Z_CENTER, Z_WIDTH is 2047.5, 30.0
  z_center = 2045.5
  z_width = 15.0
  if KEYWORD_SET(ZY_FLIP) then begin
    temp = flat.Y
    flat.Y = flat.Z
    flat.Z = temp
    flip_str = '_ZY'
    z_center = 2046.7   ; really "y_center"
  end

  ; Process using hsi_tile_plots just as for XRCF data sets

  ; Plot to .gif
  set_plot, 'X'
  window, 0, xsi=500, ysi=700
  dd_load_ct, /OTHER

  hsi_tile_plots, flat, 'HSI Flat Field, ' + $
	flat_targs(iff)+flip_str+' : Central Tile Plots', $
	QE='No QE depression correction', $
	YCHANS = ychans, YCOUNTS = ycounts, $
	Z_CENTER = z_center, Z_WIDTH = z_width

  ; Finish the .gif output
  wshow,0,iconic=0
  image = tvrd()
  tvlct, red, green, blue, /GET
  gif_file = flat_loc+'/hsi_flat_'+flat_targs(iff)+flip_str+'_tile.gif'
  write_gif, gif_file, image, red, green, blue


  ; Fit the depression bump... (as in hsi_B_tile_anal.pro)
  !p.multi = [0,1,3]
  plot, ychans, ycounts, PSYM=10
  ysel = where ( (ychans GT 1940 and ychans LT 1990) OR $
		(ychans GT 2100 AND ychans LT 2150) )
  yallsel = where(ychans GT 1940 AND ychans LT 2150)
  oplot, ychans(ysel), ycounts(ysel), PSYM=2

  fit_result = linfit(ychans(ysel), ycounts(ysel))
  oplot, [min(ychans(ysel)), max(ychans(ysel))], $
	fit_result(0)+fit_result(1)*[min(ychans(ysel)), max(ychans(ysel))], $
	LINESTYLE=2

  line_value = fit_result(0)+fit_result(1)*ychans


  ; Determine the Gaussian parameters with df_fit ...
  ngauss = 1
  x = ychans(yallsel)
  y = ycounts(yallsel) / line_value(yallsel)
  ; rms error for each bin
  rms = SQRT(ycounts(yallsel)) / line_value(yallsel)  ; counting statistics

  ; initial guess for parameters
  ; a is three params per gaussian plus the continuum level
  ; the gaussian params are: height, center, sigma
  ; a = [gauss_amplitude, gauss_center, gauss_sigma, constant]

  a = [-0.2, 2048.0, 15.0, 1.0]

  ; constant continuum
  df_continuum = 1.0 + 0.0*x  ; tricky way to get all 1.'s

  df_title = 'Fitting the "HSI QE depression"'
  ; Fit it
  df_fit,ngauss,x,y,a,rms, sig,yfit

  ; Show the fit result 
  df_showfit,x,y,a,rms, sig,yfit


  ; create and plot the Gaussian...

  yamp = a(0)
  ycent = a(1)
  ysig = a(2)
  gauss = EXP(-1.0*(ychans - ycent)^2/(2.*ysig^2))
  plttitle = flat_targs(iff)+flip_str + $
	': depress = '+STRING(yamp,FORMAT='(F8.3)')+ $
	', center = '+STRING(ycent,FORMAT='(F8.1)')+ $
	', fwhm = '+ STRING(2.35*ysig,FORMAT='(F7.2)')
  plot, ychans, ycounts/line_value, PSYM=10, $
	XRANGE=[min(ychans(ysel)), max(ychans(ysel))], XSTYLE=1, $
	YRANGE=[0.6,1.2], CHARSIZE=1.8, $
	TITLE = plttitle
  oplot, [min(ychans(ysel)), max(ychans(ysel))], [1.0,1.0], $
	LINESTYLE=2
  oplot, ychans, 1.0 + yamp * gauss

  ; Finish the .gif output
  gif_file = flat_loc+'/hsi_flat_'+flat_targs(iff)+flip_str+'_tile_anal.gif'
  wshow,0,iconic=0
  image = tvrd()
  tvlct, red, green, blue, /GET
  write_gif, gif_file, image, red, green, blue


  ; Finally, make a PHA color-coded image...
  ;
  ; Use events in a key PHA range (where the differences show up)
  pha_max = 100.0
  pha_min = 60.0
  sel = where(flat.PHA LE pha_max AND flat.PHA GT pha_min)
  xye_image, (flat(sel).Y-2048.)/10.0, (flat(sel).Z-2048.)/10.0, $
	flat(sel).PHA, $
	ONE=0.3, GIF_FILE=flat_loc+'/hsi_flat_'+flat_targs(iff)+flip_str+ $
	'_phacolor.gif'

  window, 0

end
end


RETURN
END