;;PRO hsi_B_tile_anal
;
; USAGE: to command line type:
; IDL> isel = 0  ; 0-7
; IDL> @hsi_pha_changes
;
; Look at the Y-axis projection of smooth spectra to see the
; QE dip!  isel = 1, 6, 7, 10
;
; To test for QE changes look for PHA distribution changes
; in the region which post-XRCF shows QE depression:
; Y in range 2030 to 2060, Z in range 2033, 2063
;


if n_elements(isel) EQ 0 then isel = 7

names = ['313: Al-K, Alx4 filter', '536: C-K, dispersed', $
		'748: Penning, Bex2', '586: Mg-K, dispersed', $
		'587: Mg-K, dispersed', '694: Si-K, dispersed', $
		'891: Be-K, dispersed', '897: B-K, dispersed', $
		'927: C-K, filtered', '987: C-K filtered', $
		'1208: HIREFS 0.206 keV, dispersed', '1279: Al-K, dispersed']
date_code_strs = ['970104', '970113', $
			'970120', '970116', $
			'970116', '970119', $
			'970125', '970126', $
			'970128', '970129', $
			'970207', '970208']
runid_strs = ['107545', '108577', $
			'109635', '109002', $
			'109003', '109313', $
			'110192', '110212', $
			'110299', '110391', $
			'111402', '111492']
it_strs = 'i'+['0', '0', $
		'0', '2', $
		'2', '0', $
		'0', '0', $
		'0','0', $
		'0', '0']


name = names(isel)
date_code_str = date_code_strs(isel)
runid_str = runid_strs(isel)
it_str = it_strs(isel)

runid_numerical = 0.0
READS, runid_str, runid_numerical

; Get the QE Depression value and read in the data file
qed = hsi_qe_depress(runid_numerical)
;;qed = 1.0  ; this line "turns off" the correction
hsi_ignore = xrcf_hsi_read(!DDHXDSDATA+'/phase1/data/'+date_code_str+ $
	'/hsi'+runid_str+it_str+'.fits', RAW_HSI=hsi, $
	QE_DEPRESS=qed, QE_INFO_STR=qe_info_str)

; make the HSI "tile" plots

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

plot_title = name + '  ' + date_code_str + '/' + $
		runid_str + it_str
hsi_tile_plots, hsi, plot_title, z_center = 2045.5 , z_width = 30.0, $
	YCHANS=ychans, YCOUNTS=ycounts, QE_INFO_STR = qe_info_str

; Finish the .gif output
wshow,0,iconic=0
image = tvrd()
tvlct, red, green, blue, /GET
write_gif, 'hsi'+runid_str+'_tile.gif', image, red, green, blue



; Fit a line to the B data
!p.multi = [0,1,3]
plot, ychans, ycounts, PSYM=10
ysel = where ( (ychans GT 1980 and ychans LT 2010) OR $
		(ychans GT 2090 AND ychans LT 2120) )
yallsel = where(ychans GT 1980 AND ychans LT 2120)
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 ...
@df_common
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


; and plot the Gaussian...
plot, ychans, ycounts/line_value, PSYM=10, $
	XRANGE=[min(ychans(ysel)), max(ychans(ysel))], XSTYLE=1, $
	YRANGE=[0.6,1.2]
oplot, [min(ychans(ysel)), max(ychans(ysel))], [1.0,1.0], $
	LINESTYLE=2

yamp = a(0)
ycent = a(1)
ysig = a(2)
gauss = EXP(-1.0*(ychans - ycent)^2/(2.*ysig^2))
oplot, ychans, 1.0 + yamp * gauss

; Finish the .gif output
gif_file = 'hsi'+runid_str+'_tile_anal.gif'
if qed LE 0.99 then gif_file = 'hsi'+runid_str+'_QED_tile_anal.gif'
wshow,0,iconic=0
image = tvrd()
tvlct, red, green, blue, /GET
write_gif, gif_file, image, red, green, blue