; File: zero_order_effect
; 
; Execute this IDL comand file with:  IDL> @zero_order_effect
; Note that the xrcf_hsi_read and other data gathering lines
; may be commented out (";;") after the first execution...
;
; Use HSI Al-K data to search for differences between the zero
; order image and the no-grating image.
;

; Point to where the Phase 1 data are:
phase_1_data_dir = '/nfs/atum/d3/dd/dd_d6/'

; histogram bin size
bin_size = 10.0 ; um
; size of region to use in analysis: +/- reg_size in Y and Z
; 9000 um is full HSI area, 200 um excludes LETG fine support orders
reg_size = 320. ; um  


; HRMA data set
; D-IXH-PI-3.004,5,6,7 : images from shells 1,3,4,6
; Date code: 970104, runids: 107546,7,8,9

hrma1 = xrcf_hsi_read(phase_1_data_dir+'970104/hsi107546i0.fits',/FAC,/PBLU)
hrma3 = xrcf_hsi_read(phase_1_data_dir+'970104/hsi107547i0.fits',/FAC,/PBLU)
hrma4 = xrcf_hsi_read(phase_1_data_dir+'970104/hsi107548i0.fits',/FAC,/PBLU)
hrma6 = xrcf_hsi_read(phase_1_data_dir+'970104/hsi107549i0.fits',/FAC,/PBLU)

; ------ this analysis for Marty for HSI FWHM...
; Assemble full HRMA image for Marty...
!p.multi=[0,1,2]
hrma_full_y = [hrma1.y,hrma3.y,hrma4.y,hrma6.y]
hrma_full_z = [hrma1.z,hrma3.z,hrma4.z,hrma6.z]
grating_name = 'Full HRMA'
hrma_reg_size = 150. ; um 
hrma_bin_size = 2.0 ; um
zero_y = hrma_full_y
zero_z = hrma_full_z
sel = where(ABS(zero_y) LT hrma_reg_size AND ABS(zero_z) LT hrma_reg_size)
zero_y = zero_y(sel)
zero_z = zero_z(sel)
plot, zero_y, zero_z, PSYM=3, $
	TITLE='Full HRMA Al-K image from shells 1,3,4,6', $
	XTITLE='Facility Y (um)', XRANGE=hrma_reg_size*[-1.,1.], XSTY=1, $
	YTITLE='Facility Z (um)', YRANGE=hrma_reg_size*[-1.,1.], YSTY=1
lin_hist, zero_y, hrma_bin_size, hist_axis, hist_counts
plot, hist_axis, hist_counts, PSYM=10, $
	XRANGE=hrma_reg_size*[-1.,1.], XSTYLE=1, $
	XTITLE = 'Facility Y axis (um)', $
	YTITLE = 'counts/bin, bin size = 2 um', $
	TITLE = grating_name+' Zero-order Y Projection'
; Now make an rdb file of the Y projection at 10 microns ...
lin_hist, zero_y, 10.0, bins_out, counts_out
struct_defn = {bin:0.0, count:0.0, err:0.0}
hist_out = REPLICATE(struct_defn, n_elements(bins_out))
hist_out.bin = bins_out
hist_out.count = FLOAT(counts_out)
hist_out.err = SQRT(FLOAT(counts_out))
file_out = 'hrma_y_profile.rdb'
hdr = ['# created by xrcf/zero_order_effect.pro, '+SYSTIME(), $
	'# input data is from D-IXH-PI-3.004 through 007', $
	'# date/files = 970104/hsi107546i0.fits through 107549', $
	'# bin size is 10 microns', $
	'#']
rdb_write, file_out, hist_out, HEADER=hdr
; ----------------------------------------------

;;hrma_meg_y = [hrma1.y,hrma3.y]
;;hrma_meg_z = [hrma1.z,hrma3.z]
;;hrma_heg_y = [hrma4.y,hrma6.y]
;;hrma_heg_z = [hrma4.z,hrma6.z]
;;hrma_leg_y = [hrma_meg_y,hrma_heg_y]
;;hrma_leg_z = [hrma_meg_z,hrma_heg_z]


; MEG zero-order
;;megT = xrcf_hsi_read(phase_1_data_dir+'970103/hsi107301i0.fits',/FAC,/PBLU)
;;megN = xrcf_hsi_read(phase_1_data_dir+'970103/hsi107302i0.fits',/FAC,/PBLU)
;;megB = xrcf_hsi_read(phase_1_data_dir+'970103/hsi107303i0.fits',/FAC,/PBLU)
;;megS = xrcf_hsi_read(phase_1_data_dir+'970103/hsi107304i0.fits',/FAC,/PBLU)
;;z_meg_y = [megT.y,megN.y,megB.y,megS.y]
;;z_meg_z = [megT.z,megN.z,megB.z,megS.z]

; HEG zero-order
;;hegT = xrcf_hsi_read(phase_1_data_dir+'970103/hsi107310i0.fits',/FAC,/PBLU)
;;hegN = xrcf_hsi_read(phase_1_data_dir+'970103/hsi107311i0.fits',/FAC,/PBLU)
;;hegB = xrcf_hsi_read(phase_1_data_dir+'970103/hsi107312i0.fits',/FAC,/PBLU)
;;hegS = xrcf_hsi_read(phase_1_data_dir+'970103/hsi107313i0.fits',/FAC,/PBLU)
;;z_heg_y = [hegT.y,hegN.y,hegB.y,hegS.y]
;;z_heg_z = [hegT.z,hegN.z,hegB.z,hegS.z]

; LEG zero-order
;;legT = xrcf_hsi_read(phase_1_data_dir+'970103/hsi107296i0.fits',/FAC,/PBLU)
;;legN = xrcf_hsi_read(phase_1_data_dir+'970103/hsi107297i0.fits',/FAC,/PBLU)
;;legB = xrcf_hsi_read(phase_1_data_dir+'970103/hsi107299i0.fits',/FAC,/PBLU)
;;legS = xrcf_hsi_read(phase_1_data_dir+'970103/hsi107300i0.fits',/FAC,/PBLU)
;;z_leg_y = [legT.y,legN.y,legB.y,legS.y]
;;z_leg_z = [legT.z,legN.z,legB.z,legS.z]

print, ''
print, ' HRMA-MEG events: ',n_elements(hrma_meg_y)
print, ' MEG-zerord events: ',n_elements(z_meg_y)
print, ' HRMA-HEG events: ',n_elements(hrma_heg_y)
print, ' HEG-zerord events: ',n_elements(z_heg_y)
;;print, ' HRMA-LEG events: ',n_elements(hrma_leg_y)
;;print, ' LEG-zerord events: ',n_elements(z_leg_y)
print, ''

; Set up the desired data set...
grating_name = 'MEG'  ; HEG,MEG,LEG
hrma_shells = '13'  ; 46,13,1346
hrma_y = hrma_meg_y
hrma_z = hrma_meg_z
zero_y = z_meg_y
zero_z = z_meg_z

; Use only events in the region...
sel = where(ABS(zero_y) LT reg_size AND ABS(zero_z) LT reg_size)
zero_y = zero_y(sel)
zero_z = zero_z(sel)

; Set output to postscript
pre_print_widerect
; Number of plots
;  [0,2,2] if spot plots included, [0,2,2] if not
!p.multi = [0,2,1]

; Can show a spot diagram...  or not! (comment out with ;;)
;;plot, hrma_y, hrma_z, psym=3, $
;;	XRANGE=reg_size*[-1.,1.],YRANGE=reg_size*[1.,-1.], XSTYLE=1,YSTYLE=1

ytitlestr = 'Counts (fraction per'+STRCOMPRESS(STRING(FIX(bin_size))) + $
	' um bin)'
; HRMA-only projections:
lin_hist, hrma_y, bin_size, hist_axis_hy, hist_counts_hy
;;plot_io, hist_axis_hy, hist_counts_hy/TOTAL(1.0*hist_counts_hy), PSYM=10, $
;;	XRANGE=reg_size*[-1.,1.], XSTYLE=1, YRANGE=[0.0003,0.3], YSTYLE=1, $
;;	XTITLE = 'Facility Y axis (um)', $
;;	YTITLE = ytitlestr, $
;;	TITLE = 'HRMA-'+hrma_shells+' Y Projection'
lin_hist, hrma_z, bin_size, hist_axis_hz, hist_counts_hz
;;plot_io, hist_axis_hz, hist_counts_hz/TOTAL(1.0*hist_counts_hz), PSYM=10, $
;;	XRANGE=reg_size*[-1.,1.], XSTYLE=1, YRANGE=[0.0003,0.3], YSTYLE=1, $
;;	XTITLE = 'Facility Z axis (um)', $
;;	YTITLE = ytitlestr, $
;;	TITLE = 'HRMA-'+hrma_shells+' Z Projection'

; Can show a spot diagram...  or not!
;;plot, zero_y, zero_z, psym=3, $
;;	XRANGE=reg_size*[-1.,1.],YRANGE=reg_size*[1.,-1.], XSTYLE=1,YSTYLE=1

lin_hist, zero_y, bin_size, hist_axis, hist_counts
plot_io, hist_axis, hist_counts/TOTAL(1.0*hist_counts), PSYM=10, $
	XRANGE=reg_size*[-1.,1.], XSTYLE=1, YRANGE=[0.0003,0.3], YSTYLE=1, $
	XTITLE = 'Facility Y axis (um)', $
	YTITLE = ytitlestr, $
	TITLE = grating_name+' Zero-order Y Projection'
; overplot the HRMA curve with peak channels lined up
hrma_peaks = (where(hist_counts_hy GT 0.5*MAX(hist_counts_hy)))
hrma_peak = TOTAL(1.0*hist_axis_hy(hrma_peaks)*hist_counts_hy(hrma_peaks)) / $
		TOTAL(1.0*hist_counts_hy(hrma_peaks))
zero_peaks = (where(hist_counts GT 0.5*MAX(hist_counts)))
num=TOTAL(1.0*hist_axis(zero_peaks)*hist_counts(zero_peaks))
denom = TOTAL(1.0*hist_counts(zero_peaks))
zero_peak = num / denom
print, zero_peak, hrma_peak
offset = zero_peak - hrma_peak
oplot, hist_axis_hy+offset, $
	hist_counts_hy/TOTAL(1.0*hist_counts_hy), PSYM=4
oplot, hist_axis_hy+offset, $
	hist_counts_hy/TOTAL(1.0*hist_counts_hy), PSYM=10, LINESTYLE=1
; label HRMA curve
oplot, [0.35*reg_size], [0.1], PSYM=4
xyouts, 0.40*reg_size, 0.093, '= HRMA-only', CHARSIZE=0.7

lin_hist, zero_z, bin_size, hist_axis, hist_counts
plot_io, hist_axis, hist_counts/TOTAL(1.0*hist_counts), PSYM=10, $
	XRANGE=reg_size*[-1.,1.], XSTYLE=1, YRANGE=[0.0003,0.3], YSTYLE=1, $
	XTITLE = 'Facility Z axis (um)', $
	YTITLE = ytitlestr, $
	TITLE = grating_name+' Zero-order Z Projection'
; overplot the HRMA curve with peak channels lined up
hrma_peaks = (where(hist_counts_hz GT 0.5*MAX(hist_counts_hz)))
hrma_peak = TOTAL(1.0*hist_axis_hz(hrma_peaks)*hist_counts_hz(hrma_peaks)) / $
		TOTAL(1.0*hist_counts_hz(hrma_peaks))
zero_peaks = (where(hist_counts GT 0.5*MAX(hist_counts)))
zero_peak = TOTAL(1.0*hist_axis(zero_peaks)*hist_counts(zero_peaks)) / $
		TOTAL(1.0*hist_counts(zero_peaks))
print, zero_peak, hrma_peak
offset = zero_peak - hrma_peak
oplot, hist_axis_hz+offset, $
	hist_counts_hz/TOTAL(1.0*hist_counts_hz), PSYM=4
oplot, hist_axis_hz+offset, $
	hist_counts_hz/TOTAL(1.0*hist_counts_hz), PSYM=10,LINESTYLE=1
; label HRMA curve
oplot, [0.35*reg_size], [0.1], PSYM=4
xyouts, 0.40*reg_size, 0.093, '= HRMA-only', CHARSIZE=0.7

device,/close
SPAWN, 'cp idl.ps zoe_'+grating_name+'.ps'

set_plot, 'X'
!p.multi = 0