;example program that reads SphinX FITS file, 
;calculates spectrum from raw event lists
;Then the program uses sphinx_select.pro to filter out good solar events 
;and overplots in green cleaned spectrum

P_init = !P
!P.background=255
!P.color=0
!P.charsize=1.5

;---- specify the FITS file name to read ---- 
;---- this file should be in working directory 
f = 'SPHINX_090704_044307_095331_evn_D1_L1.fts'         
;---- read SphinX FITS file ----
; read events extension
events_raw = mrdfits(f, 1, hdr_events_raw)              
; read exposureure extension
exposure_raw = mrdfits(f, 2, hdr_exposure_raw)          

;---- calculate raw lightcurve from event list -----
sphinx_spectrum, events_raw, exposure_raw, spectrum_raw ;, time , /rate        

;---- determine exposures for raw lightcurve -----
exp_raw = (exposure_raw.tstop - exposure_raw.tstart)*(exposure_raw.fracexp)

;---- divide by exposures to obtain rate spectrum [c/s] ----
spectrum_rate_raw = spectrum_raw*0D
for i=0L, n_elements(exp_raw) - 1 do begin
  spectrum_rate_raw[*, i] = spectrum_raw[*, i]/exp_raw[i]  
endfor

;---- accumulate spectra over entire file time range -----
spectrum_rawi = total(spectrum_raw, 2)
spectrum_rate_rawi = total(spectrum_rate_raw, 2)/n_elements(exp_raw)

;---- plot, spectrum vs SphinX channels ---- 
window, 0, xs=600, ys=400
erase

plot, spectrum_rate_rawi, /yl, yr=[0.001, 1.5*max(spectrum_rate_rawi)], $
/xs, /ys, ps=10, title = f, xtitle = '#channel', ytitle = 'counts/s', col=0  
;oplot, spectrum_rawi, ps=10

;---- Filter event list ----
;---- with this call sphinx_select should filter out only events of solar origin 
sphinx_select, events_raw, hdr_events_raw, $
events_filtered, exposure_filtered, min_GTI_event_num = 5 ;, light_curve = lcx  ;, /show_event_selection_plots

;---- calculate spectrum from filtered event list -----
sphinx_spectrum, events_filtered, exposure_filtered, spectrum_filtered

;---- determine exposureures for filtered lightcurve -----
exp_filtered = (exposure_filtered.tstop - exposure_filtered.tstart)*exposure_filtered.fracexp

;---- divide by exposures to obtain filtered rate spectrum [c/s] ---- 
spectrum_rate_filtered  = spectrum_filtered *0D
for i=0L, n_elements(exp_filtered) - 1 do begin
  spectrum_rate_filtered[*, i] = spectrum_filtered[*, i]/exp_filtered[i]  
endfor

;---- accumulate filtered spectra over entire file time range -----
spectrum_filteredi = total(spectrum_filtered , 2)
spectrum_rate_filteredi = total(spectrum_rate_filtered, 2)/n_elements(exp_filtered)

;---- oplot filtered rate spectrum plot ----
wh = where(spectrum_rate_filteredi gt 0, count)
if count ne 0 then begin 
  tvlct, 0, 255, 0, 128
  oplot, spectrum_rate_filteredi[wh], ps=10, col=128
endif else print, 'No valid events found during filtration' 

stop
;---- obtain DRM channel energy bounds from SphinX calibration FITS file ----
fDRM = 'SPHINX_RSP_256_nom_D1.fts'
DRM_struct = mrdfits(fDRM, 1, hdr_events_raw)
DRM = DRM_struct.matrix
ebounds = mrdfits(fDRM, 2, hdr_events_raw)


;---- calculate energy bin centers ----
energy_bin_cen = (ebounds.E_MAX + ebounds.E_MIN)/2D   

;---- plot SphinX spectrum vs energy ---- 
plot, energy_bin_cen, spectrum_rate_rawi, $
/yl, yr=[0.001, 1.5*max(spectrum_rate_rawi)], $
/xs, /ys, ps=10, title = f, xtitle = 'Energy [keV]', ytitle = 'counts/s', col=0  

;---- oplot filtered rate spectrum plot ----
wh = where(spectrum_rate_filteredi gt 0, count)
if count ne 0 then begin 
  tvlct, 0, 255, 0, 128
  oplot, energy_bin_cen[wh], spectrum_rate_filteredi[wh], ps=10, col=128
endif else print, 'No valid events found during filtration' 

stop
;---- calculate energy bin widths ----
energy_bin_widths = (ebounds.E_MAX - ebounds.E_MIN)   

;---- plot SphinX spectrum vs energy for 1 kev energy bins ---- 
plot, energy_bin_cen, spectrum_rate_rawi/energy_bin_widths , $
/yl, yr=[0.001, 1.5*max(spectrum_rate_rawi/energy_bin_widths)], $
/xs, /ys, ps=10, title = f, xtitle = 'Energy [keV]', ytitle = 'counts/s/keV', col=0  

;---- oplot filtered rate spectrum plot ----
wh = where(spectrum_rate_filteredi gt 0, count)
if count ne 0 then begin 
  tvlct, 0, 255, 0, 128
  oplot, energy_bin_cen[wh], spectrum_rate_filteredi[wh]/energy_bin_widths[wh], ps=10, col=128
endif else print, 'No valid events found during filtration' 


END