PRO eae_plot_resids, eae, pselect, LPR = lpr
; Plot the data in the eae file in various ways
;  pselect chooses the following:
;     0 - Compare the BND rates
;     1 - HRMA areas
;     2 - MEG 1 area
;     3 - HEG 1 area
;     4 - MEG 0 area
;     5 - HEG 0 area
;     6 - LEG 1 area
;     7 - LEG 0 area
;    12 - MEG 1 efficiency
;    13 - HEG 1 efficiency
;    14 - MEG 0 efficiency
;    15 - HEG 0 efficiency
;    16 - LEG 1 efficiency
;    17 - LEG 0 efficiency
;
; 6/29/98 dd Modified to use new HRMA-XRCF ea values
; 4/07/99 dd Modified to use the CALDB files
; 6/26/99 dd Rename and include residual plot

; Put the HRMA, grating predictions in common so only need to
; read once...
COMMON eae_preds, meg1,heg1,leg1,meg0,heg0,leg0, $
		meg1e,heg1e,leg1e,meg0e,heg0e,leg0e, $
		hrma, hmeg, hheg, es

; CALDBaxafcal :
; HETG Cal products are in !DDCIPSDIR
; LETG Cal products are in !DDLETGCAL+'/cip/'
; HRMA Cal products are in !DDHRMACAL+'/cip/'
; HXDS Cal products are in !DDHXDSDIR

; LEG model factor: the eae_effic_anal.pro when it applies
; the fraction corrections, explicitly corrects the LEG
; efficiency to the efficiency INCLUDING the support
; structure diffraction.  The LETG CIP 'N0002 is for this
; same efficiency.  So no additional correction is needed:
leg_fact = 1.0

if TOTAL(pselect EQ 0) GT 0 then begin

  if KEYWORD_SET(LPR) then begin
    ; For .ps output do 2 per page
    !p.multi = [0,2,1]
    pre_print_landscape
  end else begin
    ; For screen display cram them in
    !p.multi = [0,2,3]
  end

  ; BND rate comparison plot
  hn_valid = where(eae.fpc_hn_counts GT 1.)
  eaev = eae(hn_valid)
  n_rate = eaev.fpc_hn_counts/eaev.fpc_hn_time
  s_rate = eaev.fpc_hs_counts/eaev.fpc_hs_time
  t_rate = eaev.fpc_ht_counts/eaev.fpc_ht_time
  b_rate = eaev.fpc_hb_counts/eaev.fpc_hb_time
  plot_oo, n_rate, s_rate, PSYM=4, $
	TITLE = 'Comparing FPC_Hx Rates to FPC_HN Rate', $
	XTITLE = 'FPC_HN Rate (c/s)', $
	YTITLE = 'Other FPC_H Rates (c/s)', $
	XRANGE=[1.,1.E4], YRANGE=[1.,1.E4], CHARSIZE=1.4
  oplot, n_rate, t_rate, PSYM=1
  oplot, n_rate, b_rate, PSYM=2
  xvals = 10.0^(FLOAT(indgen(17))/4.0)
  ; Equal flux line:
  oplot, xvals, xvals
  xyouts, 140., 40., 'FPC_HN = ALL line'
  xyouts, 3., 600., 'FPC_HN = 36 line'
  ; HN at 36 mm line:
  oplot, xvals, xvals/( eae_bnd_hrma_area('FPC_HN','36') / $
	eae_bnd_hrma_area('FPC_HN','ALL') )

  plot_oo, eaev.energy, s_rate/n_rate, PSYM=4, $
	TITLE = 'Comparing FPC_Hx rates to FPC_HN rate', $
	XTITLE='Energy (keV)', $
	YRANGE = [0.1,10.], $
	YTITLE = 'FPC_Hx rate / FPC_HH rate', CHARSIZE=1.4
  oplot, eaev.energy, t_rate/n_rate, PSYM=1
  oplot, eaev.energy, b_rate/n_rate, PSYM=2
  oplot, [0.1,10.],[1.,1.]
  oplot, [0.1,10.], [1.,1.]/( eae_bnd_hrma_area('FPC_HN','36') / $
	eae_bnd_hrma_area('FPC_HN','ALL') )
  xyouts, 1.2, 0.5, 'FPC_HN = ALL line'
  xyouts, 0.12, 4., 'FPC_HN = 36 line'

  ; Plot the measured fluxes
  hn_flux = fltarr(n_elements(eae))
  f5_flux = fltarr(n_elements(eae))
  s5_flux = fltarr(n_elements(eae))

  ; Calculate them first
  for ir=0,n_elements(eae)-1 do begin
    ; FPC_HN
    area = eae_bnd_hrma_area('FPC_HN',eae(ir).fpc_hn_aperture)
    if area LT 0.0 OR eae(ir).fpc_hn_time LT 0.1 then begin
      hn_flux(ir) = -1.0  ; not valid      
    end else begin
      hn_flux(ir) = eae(ir).fpc_hn_counts/( $
		eae(ir).fpc_hn_time * area)
    end
    ; FPC_5
    area = eae_bnd_hrma_area('FPC_5',eae(ir).fpc_5_aperture)
    if area LT 0.0 OR eae(ir).fpc_5_time LT 0.1 then begin
      f5_flux(ir) = -1.0  ; not valid      
    end else begin
      f5_flux(ir) = eae(ir).fpc_5_counts/( $
		eae(ir).fpc_5_time * area)
    end

    ; SSD_5
    area = eae_bnd_hrma_area('SSD_5',eae(ir).ssd_5_aperture)
    if area LT 0.0 OR eae(ir).ssd_5_time LT 0.1 then begin
      s5_flux(ir) = -1.0  ; not valid      
    end else begin
      s5_flux(ir) = eae(ir).ssd_5_counts/( $
		eae(ir).ssd_5_time * area)
    end
  end

  plot_oo, hn_flux, f5_flux, PSYM=4, $
	TITLE = 'Comparing FPC_5 Flux to FPC_HN Flux', $
	XTITLE = 'FPC_HN Flux ( c / s cm^2 )', $
	YTITLE = 'FPC_5 flux ( c / s cm^2)', CHARSIZE=1.4
  xvals = 10.0^(FLOAT(indgen(8)-3))
  oplot, xvals, xvals

  plot_oo, eae.energy, f5_flux/hn_flux, PSYM=4, $
	TITLE = 'Comparing FPC_5 Flux to FPC_HN', $
	XTITLE='Energy (keV)', $
	YRANGE=[1.E-3, 1.E2], $
	YTITLE = 'FPC_5 flux / FPC_HN flux ( counts / counts )', CHARSIZE=1.4
  oplot, xvals, 1.+ 0.*xvals

  ; Overplot the FPC_5/FPC_HN ratio expected
  ; From CALDB:
  fpc_hn_qe = rdb_read(!DDHXDSDIR+'/fpc_hn.fitQe_N19980902.rdb')
  fpc5_qe = rdb_read(!DDHXDSDIR+'/fpc_5.fitQe_N19980902.rdb')
  oplot, fpc5_qe.energy, fpc5_qe.qe/fpc_hn_qe.qe
  xyouts, 0.35, 0.2, 'FPC_5 / FPC_HN'

  ; plot SSD only above 0.8 keV data sets
  ssd_range = where(eae.energy GE 0.8)

  plot_oo, f5_flux(ssd_range), s5_flux(ssd_range), PSYM=4, $
	TITLE = 'Comparing SSD_5 Flux to FPC_HN Flux', $
	XTITLE = 'FPC_HN Flux ( c / s cm^2 )', $
	YTITLE = 'SSD_5 flux ( c / s cm^2)', CHARSIZE=1.4
  xvals = 10.0^(FLOAT(indgen(8)-3))
  oplot, xvals, xvals

  plot_oo, eae(ssd_range).energy, s5_flux(ssd_range)/hn_flux(ssd_range), PSYM=4, $
	TITLE = 'Comparing SSD_5 Flux to FPC_HN (above 0.8 keV)', $
	XTITLE='Energy (keV)', XRANGE=[0.7,10.0], XSTYLE=1, $
	YRANGE=[1.E-1, 1.E1], $
	YTITLE = 'SSD_5 flux / FPC_HN flux ( counts / counts )', CHARSIZE=1.4
  oplot, xvals, 1.+ 0.*xvals

  ; Overplot the SSD/FPC ratio expected

  ; Get the SSD 5 QE
  ssd_qe = rdb_read(!DDHXDSDIR+'/ssd_5_qe_N19990208.rdb')
  ; convert energy to keV:
  ssd_qe.energy = ssd_qe.energy/1000.0
  ; convert to common ssd grid
  fpc_hn_qe_on_ssd_grid = INTERPOL(fpc_hn_qe.qe, fpc_hn_qe.energy, ssd_qe.energy)
  oplot, ssd_qe.energy, ssd_qe.qe/fpc_hn_qe_on_ssd_grid
  xyouts, 2.5, 0.4, 'SSD_QE(E) / FPC_QE(E)'
  
  if KEYWORD_SET(LPR) then begin
  ; If .ps output close it...
    device, /close
    set_plot, 'X'
  end
  ; in any case
  !p.multi = 0

end

if TOTAL(pselect GE 1) GT 0 then begin
  ; These plots include the predictions
  ; so read in the predicted effic.s etc. if not
  ; already in common...
  if n_elements(meg1) LE 0 then begin
    ; get hrma areas from CALDBaxafcal
      hrma_cip = rdb_read(!DDHRMACAL+'/cip/'+ $
	'hrmaD1996-11-01XeffareaN0004.rdb')
      es = hrma_cip.energy
      hrma = hrma_cip.ea_1 + hrma_cip.ea_3 + hrma_cip.ea_4 + $
		hrma_cip.ea_6
      hmeg = hrma_cip.ea_1 + hrma_cip.ea_3
      hheg = hrma_cip.ea_4 + hrma_cip.ea_6

    ; get grating predictions
      meg_effic = rdb_read(!DDCIPSDIR+ $
		'/hetgmegD1996-11-01greffN0003.rdb')
      meg_effic_err = rdb_read(!DDCIPSDIR+ $
		'/hetgmegD1996-11-01grerrN0003.rdb')
      heg_effic = rdb_read(!DDCIPSDIR+ $
		'/hetghegD1996-11-01greffN0003.rdb')
      heg_effic_err = rdb_read(!DDCIPSDIR+ $
		'/hetghegD1996-11-01grerrN0003.rdb')
      leg_effic = rdb_read(!DDLETGCAL+ $
		'/cip/letgD1996-11-01greffpr001N0003.rdb')
      leg_effic_err = rdb_read(!DDLETGCAL+ $
		'/cip/letgD1996-11-01grerrpr001N0003.rdb')

      meg1 = INTERPOL(meg_effic.op1, meg_effic.energy, es)
      heg1 = INTERPOL(heg_effic.op1, heg_effic.energy, es)
      leg1 = INTERPOL(leg_effic.op1, leg_effic.energy, es)
      meg0 = INTERPOL(meg_effic.oz, meg_effic.energy, es)
      heg0 = INTERPOL(heg_effic.oz, heg_effic.energy, es)
      leg0 = INTERPOL(leg_effic.oz, leg_effic.energy, es)
      meg1e = INTERPOL(meg_effic_err.op1, meg_effic.energy, es)
      heg1e = INTERPOL(heg_effic_err.op1, heg_effic.energy, es)
      leg1e = INTERPOL(leg_effic_err.op1, leg_effic.energy, es)
      meg0e = INTERPOL(meg_effic_err.oz, meg_effic.energy, es)
      heg0e = INTERPOL(heg_effic_err.oz, heg_effic.energy, es)
      leg0e = INTERPOL(leg_effic_err.oz, leg_effic.energy, es)

  end
end


if TOTAL(pselect EQ 1) GT 0 then begin
  ; HRMA no-grating areas

  ; Predictions
  plot_oo, es, hrma, XRANGE=[0.1,10.], $
	TITLE = 'HRMA-only Effective Areas', $
	YTITLE = 'Effective Area (cm^2)', $
	XTITLE = 'Energy (keV)', LINESTYLE=1, CHARSIZE=1.4
  oplot, es, hmeg
  oplot, es, hheg, LINESTYLE = 2

  ; Plot no-grating areas (MEG)
  no_grat = where(eae.order EQ 999 AND $
		STRPOS(eae.shutters,'MEG') GE 0, nwhere)
print, 'MEG no-grats:',nwhere
  if nwhere GT 0 then begin
    oplot, [eae(no_grat).energy], [eae(no_grat).effective], PSYM=4
    plot_errors, 0.98*[eae(no_grat).energy], [eae(no_grat).energy], $
	1.02*[eae(no_grat).energy], $
	[eae(no_grat).effective - eae(no_grat).effective_err_minus], $
	[eae(no_grat).effective + eae(no_grat).effective_err_plus]
    plot_errors, [eae(no_grat).energy], [eae(no_grat).energy], $
	[eae(no_grat).energy], $
	[eae(no_grat).effective - eae(no_grat).effective_syst_minus], $
	[eae(no_grat).effective + eae(no_grat).effective_syst_plus]
  end

  no_grat = where(eae.order EQ 999 AND $
		STRPOS(eae.shutters,'HEG') GE 0, nwhere)
print, 'HEG no-grats:',nwhere
  if nwhere GT 0 then begin
    oplot, [eae(no_grat).energy], [eae(no_grat).effective], PSYM=7
    plot_errors, 0.98*[eae(no_grat).energy], [eae(no_grat).energy], $
	1.02*[eae(no_grat).energy], $
	[eae(no_grat).effective - eae(no_grat).effective_err_minus], $
	[eae(no_grat).effective + eae(no_grat).effective_err_plus]
    plot_errors, [eae(no_grat).energy], [eae(no_grat).energy], $
	[eae(no_grat).energy], $
	[eae(no_grat).effective - eae(no_grat).effective_syst_minus], $
	[eae(no_grat).effective + eae(no_grat).effective_syst_plus]
  end

  no_grat = where(eae.order EQ 999 AND $
		STRPOS(eae.shutters,'ALL,ALL') GE 0,nwhere)
print, 'ALL no-grats:',nwhere
  if nwhere GT 0 then begin
    oplot, [eae(no_grat).energy], [eae(no_grat).effective], PSYM=5
    plot_errors, 0.98*[eae(no_grat).energy], [eae(no_grat).energy], $
	1.02*[eae(no_grat).energy], $
	[eae(no_grat).effective - eae(no_grat).effective_err_minus], $
	[eae(no_grat).effective + eae(no_grat).effective_err_plus]
    plot_errors, [eae(no_grat).energy], [eae(no_grat).energy], $
	[eae(no_grat).energy], $
	[eae(no_grat).effective - eae(no_grat).effective_syst_minus], $
	[eae(no_grat).effective + eae(no_grat).effective_syst_plus]
  end

end

if TOTAL(pselect EQ 2) GT 0 then begin
  ; HRMA grating areas

  ; Predictions
  plot_oo, es, meg1*hmeg, XRANGE=[0.1,10.], $
	TITLE = 'MEG First-Order Effective Areas', $
	YTITLE = 'Effective Area (cm^2)', YRANGE=[1.0,100.], $
	XTITLE = 'Energy (keV)', CHARSIZE=1.4
  oplot, es, (meg1+meg1e)*hmeg, LINESTYLE=1
  oplot, es, (meg1-meg1e)*hmeg, LINESTYLE=1
  ; MEG data
  w_grat = where(ABS(eae.order) EQ 1 AND $
		STRPOS(eae.shutters,'MEG') GE 0, nwhere)
print, 'MEG 1st order meas.s:',nwhere
  if nwhere GT 0 then begin
    oplot, [eae(w_grat).energy], [eae(w_grat).effective], PSYM=4
    plot_errors, 0.951*[eae(w_grat).energy], [eae(w_grat).energy], $
	1.049*[eae(w_grat).energy], $
	[eae(w_grat).effective - eae(w_grat).effective_err_minus], $
	[eae(w_grat).effective + eae(w_grat).effective_err_plus]
    plot_errors, [eae(w_grat).energy], [eae(w_grat).energy], $
	[eae(w_grat).energy], $
	[eae(w_grat).effective - eae(w_grat).effective_syst_minus], $
	[eae(w_grat).effective + eae(w_grat).effective_syst_plus]
  end

end

if TOTAL(pselect EQ 3) GT 0 then begin
  ; HRMA grating areas

  ; Plot HRMA-grating areas (FPC only)

  ; Predictions
  plot_oo, es, heg1*hheg, XRANGE=[0.1,10.], $
	TITLE = 'HEG First-Order Effective Area', $
	YTITLE = 'Effective Area (cm^2)', YRANGE=[1.0,100.], $
	XTITLE = 'Energy (keV)', CHARSIZE=1.4
  oplot, es, (heg1+heg1e)*hheg, LINESTYLE=1
  oplot, es, (heg1-heg1e)*hheg, LINESTYLE=1

  ; HEG data
  w_grat = where(ABS(eae.order) EQ 1 AND $
		STRPOS(eae.shutters,'HEG') GE 0, nwhere)
print, 'HEG 1st order meas.s:',nwhere
  if nwhere GT 0 then begin
    oplot, [eae(w_grat).energy], [eae(w_grat).effective], PSYM=4
    plot_errors, 0.951*[eae(w_grat).energy], [eae(w_grat).energy], $
	1.049*[eae(w_grat).energy], $
	[eae(w_grat).effective - eae(w_grat).effective_err_minus], $
	[eae(w_grat).effective + eae(w_grat).effective_err_plus]
    plot_errors, [eae(w_grat).energy], [eae(w_grat).energy], $
	[eae(w_grat).energy], $
	[eae(w_grat).effective - eae(w_grat).effective_syst_minus], $
	[eae(w_grat).effective + eae(w_grat).effective_syst_plus]
  end
end

if TOTAL(pselect EQ 4) GT 0 then begin
  ; HRMA grating areas

  ; Plot HRMA-grating areas (FPC only)

  ; Predictions
  plot_oo, es, meg0*hmeg, XRANGE=[0.1,10.], $
	TITLE = 'MEG Zero-Order Effective Area', $
	YTITLE = 'Effective Area (cm^2)', YRANGE=[4.,400.], YSTYLE=1,$
	XTITLE = 'Energy (keV)', CHARSIZE=1.4
  oplot, es, (meg0+meg0e)*hmeg, LINESTYLE=1
  oplot, es, (meg0-meg0e)*hmeg, LINESTYLE=1

  ; MEG data
  w_grat = where(eae.order EQ 0 AND $
		STRPOS(eae.shutters,'MEG') GE 0, nwhere)
print, 'MEG 0th order meas.s:',nwhere
  if nwhere GT 0 then begin
    oplot, [eae(w_grat).energy], [eae(w_grat).effective], PSYM=4
    plot_errors, 0.951*[eae(w_grat).energy], [eae(w_grat).energy], $
	1.049*[eae(w_grat).energy], $
	[eae(w_grat).effective - eae(w_grat).effective_err_minus], $
	[eae(w_grat).effective + eae(w_grat).effective_err_plus]
    plot_errors, [eae(w_grat).energy], [eae(w_grat).energy], $
	[eae(w_grat).energy], $
	[eae(w_grat).effective - eae(w_grat).effective_syst_minus], $
	[eae(w_grat).effective + eae(w_grat).effective_syst_plus]
  end

end

if TOTAL(pselect EQ 5) GT 0 then begin
  ; HRMA grating areas

  ; Predictions
  plot_oo, es, heg0*hheg, XRANGE=[0.1,10.], $
	TITLE = 'HEG Zero-Order Effective Area', $
	YTITLE = 'Effective Area (cm^2)', YRANGE=[2.,200.], YSTYLE=1, $
	XTITLE = 'Energy (keV)', CHARSIZE=1.4
  oplot, es, (heg0+heg0e)*hheg, LINESTYLE=1
  oplot, es, (heg0-heg0e)*hheg, LINESTYLE=1

  ; HEG data
  w_grat = where(eae.order EQ 0 AND $
		STRPOS(eae.shutters,'HEG') GE 0, nwhere)
print, 'HEG 0th order meas.s:',nwhere
  if nwhere GT 0 then begin
    oplot, [eae(w_grat).energy], [eae(w_grat).effective], PSYM=4
    plot_errors, 0.951*[eae(w_grat).energy], [eae(w_grat).energy], $
	1.049*[eae(w_grat).energy], $
	[eae(w_grat).effective - eae(w_grat).effective_err_minus], $
	[eae(w_grat).effective + eae(w_grat).effective_err_plus]
    plot_errors, [eae(w_grat).energy], [eae(w_grat).energy], $
	[eae(w_grat).energy], $
	[eae(w_grat).effective - eae(w_grat).effective_syst_minus], $
	[eae(w_grat).effective + eae(w_grat).effective_syst_plus]
  end
end


if TOTAL(pselect EQ 6) GT 0 then begin
  ; HRMA grating areas

  ; Plot HRMA-grating areas (FPC only)

  ; Predictions
  plot_oo, es, leg_fact*leg1*hrma, XRANGE=[0.1,10.], $
	TITLE = 'LEG First-Order Effective Areas', $
	YTITLE = 'Effective Area (cm^2)', $
	YRANGE=[2.,200.], YSTYLE=1, $
	XTITLE = 'Energy (keV)', CHARSIZE=1.4
  oplot, es, (leg_fact*leg1+leg_fact*leg1e)*hrma, LINESTYLE=1
  oplot, es, (leg_fact*leg1-leg_fact*leg1e)*hrma, LINESTYLE=1

  ; LEG data
  w_grat = where(ABS(eae.order) EQ 1 AND $
		STRPOS(eae.grating,'LE') GE 0, nwhere)
print, 'LEG 1st order meas.s:',nwhere
  if nwhere GT 0 then begin
    oplot, [eae(w_grat).energy], [eae(w_grat).effective], PSYM=4
    plot_errors, 0.951*[eae(w_grat).energy], [eae(w_grat).energy], $
	1.049*[eae(w_grat).energy], $
	[eae(w_grat).effective - eae(w_grat).effective_err_minus], $
	[eae(w_grat).effective + eae(w_grat).effective_err_plus]
    plot_errors, [eae(w_grat).energy], [eae(w_grat).energy], $
	[eae(w_grat).energy], $
	[eae(w_grat).effective - eae(w_grat).effective_syst_minus], $
	[eae(w_grat).effective + eae(w_grat).effective_syst_plus]
  end
end

if TOTAL(pselect EQ 7) GT 0 then begin
  ; HRMA grating areas

  ; Plot HRMA-grating areas (FPC only)

  ; Predictions
  plot_oo, es, leg_fact*leg0*hrma, XRANGE=[0.1,10.], $
	TITLE = 'LEG Zero-Order Effective Areas', $
	YTITLE = 'Effective Area (cm^2)', $
	YRANGE=[4.,400.], YSTYLE=1, $
	XTITLE = 'Energy (keV)', CHARSIZE=1.4
  oplot, es, (leg_fact*leg0+leg_fact*leg0e)*hrma, LINESTYLE=1
  oplot, es, (leg_fact*leg0-leg_fact*leg0e)*hrma, LINESTYLE=1

  ; LEG data
  w_grat = where(eae.order EQ 0 AND $
		STRPOS(eae.grating,'LE') GE 0, nwhere)
print, 'LEG 0th order meas.s:',nwhere
  if nwhere GT 0 then begin
    oplot, [eae(w_grat).energy], [eae(w_grat).effective], PSYM=4
    plot_errors, 0.951*[eae(w_grat).energy], [eae(w_grat).energy], $
	1.049*[eae(w_grat).energy], $
	[eae(w_grat).effective - eae(w_grat).effective_err_minus], $
	[eae(w_grat).effective + eae(w_grat).effective_err_plus]
    plot_errors, [eae(w_grat).energy], [eae(w_grat).energy], $
	[eae(w_grat).energy], $
	[eae(w_grat).effective - eae(w_grat).effective_syst_minus], $
	[eae(w_grat).effective + eae(w_grat).effective_syst_plus]
  end
end


if TOTAL(pselect EQ 12) GT 0 then begin
  ; Efficiency

  ; Select the data
  w_grat = where(ABS(eae.order) EQ 1 AND $
		STRPOS(eae.shutters,'MEG') GE 0 AND $
		eae.effic GT 0.0, nwhere)
  print, 'MEG 1st order meas.s:',nwhere

  ; Set parameters for plotting
  this_title = 'MEG First-Order Efficiency'
  x_title = 'Energy (keV)'
  y_title = 'Efficiency'
  x_range = [0.1, 10.]
  y_range = [0.003, 0.3]
  char_size = 0.8
  resid_frac = 0.3
  ;
  model_es = es
  model = meg1
  model_errp = meg1e
  model_errm = meg1e
  ;
  data_es = [eae(w_grat).energy]
  data = [eae(w_grat).effic]
  data_errp = [eae(w_grat).effic_err_minus]
  data_errm = [eae(w_grat).effic_err_plus]
  data_systp = [eae(w_grat).effic_syst_minus]
  data_systm = [eae(w_grat).effic_syst_plus]
  resid_range = [-60.0, 60.0]

end

if TOTAL(pselect EQ 13) GT 0 then begin
  ; Efficiency

  ; Select the data
  ; HEG data
  w_grat = where(ABS(eae.order) EQ 1 AND $
		STRPOS(eae.shutters,'HEG') GE 0 AND $
		eae.effic GT 0.0, nwhere)
  print, 'HEG 1st order meas.s:',nwhere

  ; Set parameters for plotting
  this_title = 'HEG First-Order Efficiency'
  x_title = 'Energy (keV)'
  y_title = 'Efficiency'
  x_range = [0.1, 10.]
  y_range = [0.003, 0.3]
  char_size = 0.8
  resid_frac = 0.3
  ;
  model_es = es
  model = heg1
  model_errp = heg1e
  model_errm = heg1e
  ;
  data_es = [eae(w_grat).energy]
  data = [eae(w_grat).effic]
  data_errp = [eae(w_grat).effic_err_minus]
  data_errm = [eae(w_grat).effic_err_plus]
  data_systp = [eae(w_grat).effic_syst_minus]
  data_systm = [eae(w_grat).effic_syst_plus]
  resid_range = [-60.0, 60.0]

end

if TOTAL(pselect EQ 14) GT 0 then begin
  ; Efficiency

  ; Select the data
  ; MEG data
  w_grat = where(eae.order EQ 0 AND $
		STRPOS(eae.shutters,'MEG') GE 0 AND $
		eae.effic GT 0.0, nwhere)
  print, 'MEG 0th order meas.s:',nwhere

  ; Set parameters for plotting
  this_title = 'MEG Zero-Order Efficiency'
  x_title = 'Energy (keV)'
  y_title = 'Efficiency'
  x_range = [0.1, 10.]
  y_range = [0.01, 1.0]
  char_size = 0.8
  resid_frac = 0.3
  ;
  model_es = es
  model = meg0
  model_errp = meg0e
  model_errm = meg0e
  ;
  data_es = [eae(w_grat).energy]
  data = [eae(w_grat).effic]
  data_errp = [eae(w_grat).effic_err_minus]
  data_errm = [eae(w_grat).effic_err_plus]
  data_systp = [eae(w_grat).effic_syst_minus]
  data_systm = [eae(w_grat).effic_syst_plus]
  resid_range = [-60.0, 60.0]

end

if TOTAL(pselect EQ 15) GT 0 then begin
  ; Efficiency

  ; Select the data
  ; HEG data
  w_grat = where(eae.order EQ 0 AND $
		STRPOS(eae.shutters,'HEG') GE 0 AND $
		eae.effic GT 0.0, nwhere)
  print, 'HEG 0th order meas.s:',nwhere

  ; Set parameters for plotting
  this_title = 'HEG Zero-Order Efficiency'
  x_title = 'Energy (keV)'
  y_title = 'Efficiency'
  x_range = [0.1, 10.]
  y_range = [0.01, 1.0]
  char_size = 0.8
  resid_frac = 0.3
  ;
  model_es = es
  model = heg0
  model_errp = heg0e
  model_errm = heg0e
  ;
  data_es = [eae(w_grat).energy]
  data = [eae(w_grat).effic]
  data_errp = [eae(w_grat).effic_err_minus]
  data_errm = [eae(w_grat).effic_err_plus]
  data_systp = [eae(w_grat).effic_syst_minus]
  data_systm = [eae(w_grat).effic_syst_plus]
  resid_range = [-60.0, 60.0]

end


if TOTAL(pselect EQ 16) GT 0 then begin
  ; Efficiency

  ; Select the data
  ; LEG data
  w_grat = where(ABS(eae.order) EQ 1 AND $
		STRPOS(eae.grating,'LE') GE 0 AND $
		eae.effic GT 0.0, nwhere)
  print, 'LEG 1st order meas.s:',nwhere

  ; Set parameters for plotting
  this_title = 'LEG First-Order Efficiency'
  x_title = 'Energy (keV)'
  y_title = 'Efficiency'
  x_range = [0.1, 10.]
  y_range = [0.003, 0.3]
  char_size = 0.8
  resid_frac = 0.3
  ;
  model_es = es
  model = leg1
  model_errp = leg1e
  model_errm = leg1e
  ;
  data_es = [eae(w_grat).energy]
  data = [eae(w_grat).effic]
  data_errp = [eae(w_grat).effic_err_minus]
  data_errm = [eae(w_grat).effic_err_plus]
  data_systp = [eae(w_grat).effic_syst_minus]
  data_systm = [eae(w_grat).effic_syst_plus]
  resid_range = [-60.0, 60.0]

end

if TOTAL(pselect EQ 17) GT 0 then begin
  ; Efficiency

  ; Select the data
  ; LEG data
  w_grat = where(eae.order EQ 0 AND $
		STRPOS(eae.grating,'LE') GE 0 AND $
		eae.effic GT 0.0, nwhere)

  print, 'LEG 0th order meas.s:',nwhere

  ; Set parameters for plotting
  this_title = 'LEG Zero-Order Efficiency'
  x_title = 'Energy (keV)'
  y_title = 'Efficiency'
  x_range = [0.1, 10.]
  y_range = [0.01, 1.0]
  char_size = 0.8
  resid_frac = 0.3
  ;
  model_es = es
  model = leg0
  model_errp = leg0e
  model_errm = leg0e
  ;
  data_es = [eae(w_grat).energy]
  data = [eae(w_grat).effic]
  data_errp = [eae(w_grat).effic_err_minus]
  data_errm = [eae(w_grat).effic_err_plus]
  data_systp = [eae(w_grat).effic_syst_minus]
  data_systm = [eae(w_grat).effic_syst_plus]
  resid_range = [-60.0, 60.0]

end

if ( (TOTAL(pselect EQ 12) GT 0) OR $
	(TOTAL(pselect EQ 13) GT 0) OR $
	(TOTAL(pselect EQ 14) GT 0) OR $
	(TOTAL(pselect EQ 15) GT 0) OR $
	(TOTAL(pselect EQ 16) GT 0) OR $
	(TOTAL(pselect EQ 17) GT 0)  ) then begin
; ------------------------------------------
  ; Two plots
  !p.multi = [0,1,2]

  ; Set plot region for main plot
  !p.region = [0.,resid_frac, 1.0, 1.0]

  ; Plot the Model
  if n_elements(model) GT 0 then begin
    plot_oo, model_es, model, XRANGE=x_range, XSTYLE=1, $
	TITLE = this_title, $
	YTITLE = y_title, YRANGE=y_range, YSTYLE=1, $
	XTITLE = x_title, CHARSIZE=char_size
    if n_elements(model_errp) GT 0 then begin
      if n_elements(model_errm) EQ 0 then model_errm = model_errp
      oplot, model_es, model+model_errm, LINESTYLE=1
      oplot, model_es, model-model_errm, LINESTYLE=1
    end
  end else begin
    plot_oo, [0.,1.], [0.,1.], /NODATA, XRANGE=x_range, XSTYLE=1, $
	TITLE = this_title, $
	YTITLE = y_title, YRANGE=y_range, YSTYLE=1, $
	XTITLE = x_title, CHARSIZE=char_size
  end

  ; Plot the Data w/errors
  if n_elements(data) GT 0 then begin
    oplot, data_es, data, PSYM=4
    if n_elements(data_errp) GT 0 then begin
      if n_elements(data_errm) EQ 0 then data_errm = data_errp
      bar_frac = EXP(ALOG(x_range(1)/x_range(0))/80.0)
      plot_errors, data_es/bar_frac, data_es, $
	data_es*bar_frac, $
	[data - data_errm], $
	[data + data_errp]
    end
    if n_elements(data_systp) GT 0 then begin
      if n_elements(data_systm) EQ 0 then data_systm = data_systp
      plot_errors, data_es, data_es, $
	data_es, $
	[data - data_systm], $
	[data + data_systp]
    end
  end

  ; Make the residuals plot
  !p.region = [0.0,0.02, 1.0, resid_frac]
  plot_oi, x_range, [0.0,0.0], XRANGE=x_range, XSTYLE=1, $
	YTITLE = 'Residuals (%)', YRANGE=resid_range, YSTYLE=1, $
	XTITLE = x_title, $
	LINESTYLE=1, $
	CHARSIZE=char_size

  ; Put the data points on it...
  if n_elements(data) GT 0 then begin
    ; Evaluate the model at the data
    if n_elements(model) GT 0 then begin
      model_at_data = INTERPOL_SORT(model, model_es, data_es)
    end else begin
      ; set the model_at_data to the data to plot the errors on linear
      ; scale in resid plot...
      model_at_data = data
    end
    data_model_error = 100.0*(data - model_at_data)/model_at_data
    oplot, data_es, data_model_error, PSYM=4
    if n_elements(data_errp) GT 0 then begin
      plot_errors, data_es/bar_frac, data_es, $
	data_es*bar_frac, $
	[100.0*((data-data_errm) - model_at_data)/model_at_data], $
	[100.0*((data+data_errp) - model_at_data)/model_at_data]
    end
    if n_elements(data_systp) GT 0 then begin
      plot_errors, data_es, data_es, $
	data_es, $
	[100.0*((data-data_systm) - model_at_data)/model_at_data], $
	[100.0*((data+data_systp) - model_at_data)/model_at_data]
    end
  end

  ; and add the model error lines
  if n_elements(model) GT 0 then begin
    oplot, model_es, 100.0*model_errp/model, LINESTYLE=1
    oplot, model_es, -100.0*model_errm/model, LINESTYLE=1
  end

  plot_creator_label, 'ddidl/xrcf/eae_plot_resids.pro'

  !p.region = 0
  !p.multi = 0
; ------------------------------------------
end

RETURN
END