PRO heg_level1_spec
;
; execute at IDL command line with @heg_level1_spec

;--------------------------------------
;--------------------------------------
; Use SAOSAC ray file 8.04778_000_000.fits as
; flight prediction and calculate expected
; line profile...
; MARX run output in /spectra/d7/heg_level1_sim/saosac
xye_image, x, y, e, MARX='/spectra/d7/heg_level1_sim/saosac'
r = sqrt(x^2 + y^2)

; Measure fractions in 1,3,10 arc sec diameter
print, ''
print, ' HRMA Encircled fraction, SAOSAC 8.04 keV'
print, '     Diameter    Fraction'
dias = [1.,3.,10.]
for id=0,2 do begin
  inside = where(r LT 50.E-3*dias(id)/2.0)
  print, dias(id), FLOAT(n_elements(inside))/FLOAT(n_elements(r))
end
print, ''


;window,0
pre_print_landscape

!p.multi=[0,3,1]
print, ' FWHM vs Other blur'
blurs = [0.,10.6, 14.5]  ; um
btitles = ['SAOSAC/Flight, 8.04778 keV', $
	'Optimistic, FWHM ~ 44 um',$
	'Conservative, FWHM ~ 52 um']
for ibs=0,2 do begin
  ;add gaussian blur to x values
  blur = fltarr(n_elements(x))
  for ib=0,n_elements(x)-1 do blur(ib) = blurs(ibs)*0.001*g_rand(seed)
  xb = x+blur
  bin_width = 0.001
  lin_hist, xb(where(ABS(xb) LT 0.1)),bin_width,bin_values_out, bin_counts_out
  plot, bin_width*bin_values_out*1000.0, bin_counts_out,PSYM=10, $
	YRANGE=[0.,400.], $
	title = btitles(ibs), CHARSIZE=1.8, XTITLE='Y location (um)',$
	ytitle = 'counts'
  ; Measure FWHM
  maxofdist = MAX(bin_counts_out)
  maxofdist = bin_counts_out(where(bin_counts_out GE 0.8*maxofdist))
  maxofdist = TOTAL(maxofdist)/FLOAT(n_elements(maxofdist))
  oplot, [-0.05, 0.05], maxofdist*[1.,1.], LINESTYLE=2
  oplot, [-0.05, 0.05], maxofdist*[1.,1.]/2., LINESTYLE=2
  above = where(bin_counts_out GE maxofdist/2)
  print, STRING(blurs(ibs)),' um  -->   ', $
	STRING(1000.*bin_width*(MAX(above)-MIN(above)),FORMAT='(F7.3)'), ' um'
end
device,/close
set_plot, 'X'
!p.multi=0



;--------------------------------------
;--------------------------------------
; Use ebrt routines to see current error budget inputs
; eb_reset, 0
; eb_ps(eb_EkeV) = 8.0
; eb_update
; eb_dump
;
; Find out that currently sigma_H at 8 keV is 32.6683 um
; and other sigmas are 10.6 (14.5) optimistic (conservative).



;--------------------------------------
;--------------------------------------
; Try to analyze XRCF MC data set for as-measured E/dE at Cu-K
; 
;;print, ''
;;print, ' Analyze XRCF test H-HAS-MC-3.001 for HEG E/dE at Cu-K (8.04 keV)'
;;print, ''

;--------------------------------------
; Get the PSU, dd, processed data:
; (takes a while, comment it out after doing it once...)
;;restore, '/spectra/d6/ACIS_anal/H-HAS-MC-3.001_l1.idlsav'

; Plot an "image" of selected events:
;;sel = where( (l1.grade EQ 0 OR l1.grade EQ 3 or l1.grade EQ 4) AND $
;;	l1.energy GT 6.0 AND l1.energy LT 10.0 AND $
;;	ABS(l1.detx) GT 5.0 AND ABS(l1.detx) LT 9.0)

;;print, ' Number of selected events = ',n_elements(sel)
;;print, ''

; See where they are in time...
;;plot, l1(sel).detx, PSYM=3


; Select a stable time period:

; with grade 0 only:
;;l1use = l1(sel(6000:10000))

; including grades 3 and 4:
;;l1use = l1(sel(12000:20000))
; or 
;;l1use = l1(sel(21000:29000))

; and plot the final selected ones
;;plot, l1use.detx, PSYM=3

; Make histograms of the first orders
;;lin_hist, l1use.detx, 0.002, bin_values_out, bin_counts_out
;;plot, bin_values_out, bin_counts_out, PSYM=10,XRANGE=[6.1,6.4]
;;plot, bin_values_out, bin_counts_out, PSYM=10,XRANGE=[-7.4,-7.]


;------------------------
; Need to figure out when NS and TB shutters are...
; Get the PSU and the l1 data...
;;psu_view, 'H-HAS-MC-3.001', psu=psu, l1=l1


; EXPOSURE between 1500 2000 and 2350 2800

;;sel1 = where( l1.energy GT 6.0 AND l1.energy LT 10.0 AND $
;;		(l1.grade EQ 0 OR l1.grade EQ 3 OR l1.grade EQ 4) AND $
	;;	(psu.EXPOSURE GT 1500 AND psu.EXPOSURE LT 2000) )

;;sel2 = where( l1.energy GT 6.0 AND l1.energy LT 10.0 AND $
;;		(l1.grade EQ 0 OR l1.grade EQ 3 OR l1.grade EQ 4) AND $
;;		(psu.EXPOSURE GT 2350 AND psu.EXPOSURE LT 2800) )

;;print, n_elements(sel1)
;;print, n_elements(sel2)

;;pre_print_landscape
;;plot, l1(sel1).detx, PSYM=3, YRANGE=[-10.,10.]
;;plot, l1(sel2).detx, PSYM=3, YRANGE=[-10.,10.]

;;plot, l1(sel1(2000:14000)).detx, PSYM=3, YRANGE=[-10.,10.]
;;xvals = l1(sel1(2000:14000)).detx
;;lin_hist, xvals(where(ABS(xvals) LT 10.)), 0.002, bin_values_out, bin_counts_out
;;;pre_print_landscape
;;plot, bin_values_out, bin_counts_out, PSYM=10,XRANGE=[6.1,6.4], $
	TITLE='EXPOSURE 1500 2000'
;;plot, bin_values_out, bin_counts_out, PSYM=10,XRANGE=[-7.3,-7.], $
	TITLE='EXPOSURE 1500 2000'
;;;post_print


;;plot, l1(sel2(2500:12000)).detx, PSYM=3, YRANGE=[-10.,10.]
;;xvals = l1(sel2(2500:12000)).detx
;;lin_hist, xvals(where(ABS(xvals) LT 10.)), 0.002, $
;;	bin_values_out, bin_counts_out

;pre_print_landscape
;;plot, bin_values_out, bin_counts_out, PSYM=10,XRANGE=[6.3,6.6], $
;;	TITLE='EXPOSURE 2350 2800'
;;plot, bin_values_out, bin_counts_out, PSYM=10,XRANGE=[-7.2,-6.9], $
;;	TITLE='EXPOSURE 2350 2800'
;post_print

;;device,/close
;;set_plot, 'X'

RETURN
END