PRO feat_fraction, eae, NO_GRATING=no_grating, $
	NO_PLOT=no_plot, PS_PLOT=ps_plot, $
	FRACTION_OUT=fraction_out, AP_EVENTS = ap_events, $
	DYFAC = dyfac, DZFAC = dzfac
;
; KEYWORDS:
;    DYFAC - in mm: amount to offset the aperture by (from nominal MARX
;            simulated value) for the as-run condition
;    DZFAC - in mm: ditto
;
;
; The "feature fraction" is given as:
;
;                      all counts due to the feature
;  feature_fraction = -------------------------------
;                      observed counts from analysis
;
; so that multiplying the observed counts by the "feature_fraction"
; gives the feature counts.
;
; The "observed counts" depends on location, aperture, and analysis
; method.  Currently the location is based on .loc files,
; the aperture is circular from CMDB (eae file), and ROI analysis
; is used (mimicing the current eae analysis.)
;
; The fraction is calculated for the nominal aperture and location
; as well as for +/- aperture size variation and +/- [Y|Z]_xrcf offset
; errors.  FRACTION_OUT is an array of these seven values plus the
; number of feature counts and the Observed-feature fraction
; at nominal location:
;  fracNom fracDplus fracDminus fracYplus fracYminus 
;       fracYplus fracYminus Nfeature ObsFeatFrac

; The fractions are based on simulations contained in:
marx_dir = !DDPHASE1SIM

; Input parameters:
; Generate the parameters from the eae structure:
;        feat_fraction, eae(nnn)
; OR pass them explicitly with a structure like:
;        feat_fraction, {trw_id:'D-LXH-dF-16.001',iteration:'2',
;                       energy:1.254,fp_aperture:10000.,order:5,grating:'MEG',
;                       fp_detector:'SSD',focal_analysis:'ROI,1100-1400'}
;
trw_id = eae.trw_id
if KEYWORD_SET(NO_GRATING) then trw_id = trw_id+'NG'
it_str = 'i'+eae.iteration
diameter = eae.fp_aperture/1000.0
energy = eae.energy
order = eae.order
grating = eae.grating
detector = eae.fp_detector
anal_str = eae.focal_analysis  ; e.g., 'ROI,800-1400'

; Internal parameters:
; SSD energy per channel
ssd_const = 0.0056  ; keV/channel
; the spectral feature is E +/- de
de = 0.01 * energy

; Offset parameters
if n_elements(DYFAC) EQ 0 then DYFAC = 0.0
if n_elements(DZFAC) EQ 0 then DZFAC = 0.0

; Variation parameters
; Diameter error:
dia_error = 0.025 ; mm
; Y-pos error:
y_error = 0.050 ; mm
; Z-pos error:
z_error = 0.050 ; mm

; Derived parameters
; If the focal analysis is ROI, then set the limits
; based on the ROI values and the detector type,
; otherwise use a guess:
roi_low = 0.5*energy  ; these are in keV and apply to the 
roi_high = 1.5*energy ; simulated data
; Check for ROI
roiloc = STRPOS(anal_str,'ROI,')
if roiloc GE 0 then begin
  ; Get the remaining string in the field
  roi_str = STRMID(anal_str,roiloc+4,80)
  ; Get the ROI limits if they are there
  ; e.g., in format ROI,12-120
  pieces = STR_SEP(roi_str,'-')
  if n_elements(pieces) GE 2 then begin
    chan_low = 0
    READS, pieces(0), chan_low
    chan_high = 0
    READS, pieces(1), chan_high
    ; Convert these detector-specific limits to keV
    CASE detector of
    'SSD': begin
        roi_low = ssd_const * chan_low
        roi_high = ssd_const * chan_high
      end
    'FPC2': begin
        chan_ave = (chan_low + chan_high)/2.
        if energy GT 4.0 then begin
          ; includes escape peak...
          fpc_const = (energy-1.5)/chan_ave
        end else begin
          fpc_const = energy/chan_ave
        end
        roi_low = fpc_const * chan_low
        roi_high = fpc_const * chan_high
      end
    else: begin
        print, ' * Can not convert chans to keV for detector = '+detector
        RETURN
      end
    endcase
  end  ; of "ROI,12-120" format
  ; Check for "ROI,All" as used for HSI
  if STRPOS(roi_str,'All') GE 0 then begin
    roi_low = 0.0
    roi_high = 30.0
  end
end


; Report parameters:
print, ''
print, ' . . . . . . . . . . . . . . . . . . . . '
print, '    TRW_ID : ',trw_id+'/'+it_str
print, '  (Grating : '+grating+' )'
print, '    (Order : '+STRING(FIX(order),FORMAT=I4)+' )'
print, '  Detector : '+detector
print, '  Diameter : ',diameter
if DYFAC NE 0.0 then print, '    dY_fac : ', DYFAC
if DYFAC NE 0.0 then print, '    dZ_fac : ', DZFAC
feature_str = STRING(energy,FORMAT='(F8.4)')+'  +/-'+ $
		STRING(de,FORMAT='(F7.4)')
print, '   Feature : ', feature_str
print, '  ROI(keV) : ',roi_low, roi_high
print, ' . . . . . . . . . . . . . . . . . . . . '
print, ' dDiameter : ',dia_error
print, '    dYxrcf : ',y_error

; get the simulated dataset:
; it will be in the structure ml1 with tags: detx, dety, energy
restore, marx_dir + '/' + trw_id + '/'+ it_str +'/ml1.idlsav'

; Counts in feature.  Does not depend on the aperture.
; use E+/-de and spatially +/- 0.33*order_spacing
xrcf_params, HC=hc, PERIODS=periods, ROWLANDS=rowlands
CASE grating OF
  'LEG': begin
    order_spacing = rowlands(0)*(hc/energy)/periods(0)
  end
  'MEG': begin
    order_spacing = rowlands(1)*(hc/energy)/periods(1)
  end
  'HEG': begin
    order_spacing = rowlands(2)*(hc/energy)/periods(2)
  end
  'NONE': begin
    order_spacing = 6.0/0.33  ; mm everything in simulation
  end
  ELSE: begin
    print, '* feat_fraction: unknown grating: '+grating
    RETURN
  end
ENDCASE

; ALL FEATURE COUNTS
all_feature_counts = where(ABS(ml1.energy - energy) LT de AND $
	ABS(ml1.detx) LT 0.33*order_spacing, nfeature)

print, ' N_feature : ', nfeature
print, ' . . . . . . . . . . . . . . . . . . . . '
print, ''

; Calculate the observed counts and fractions for several
; cases...
;  fracNom fracDplus fracDminus fracYplus fracYminus fracZplus fracZminus
yvals = 0. + y_error*[0.,0.,0.,1.,-1.,0., 0.]
zvals = 0. + z_error*[0.,0.,0.,0., 0.,1.,-1.]
diams = diameter + dia_error*[0.,1.,-1.,0.,0.,0.,0.]

; Create output array,
; include the number of feature counts, and the
; observed_feature fraction at nominal
; parameters.
fraction_out = FLTARR(n_elements(diams)+2)
; loop over the cases
for ic=0, n_elements(diams)-1 do begin
  ; how far are the events from the center?
  radius = SQRT( (ml1.detx-(DYFAC+yvals(ic)))^2 + $
			(ml1.dety-(DZFAC+zvals(ic)))^2 )

  ; OBSERVED COUNTS
  observed_counts = where(ml1.energy GT roi_low AND $
				ml1.energy LT roi_high AND  $
				radius LT diams(ic)/2.0, nobs)

  ; OBSERVED FEATURE COUNTS are all counts
  ; in the feature that are observed
  obsfeat_counts = where(  (ABS(ml1.energy - energy) LT de AND $
	ABS(ml1.detx) LT 0.33*order_spacing)  AND  $
			(ml1.energy GT roi_low AND $
				ml1.energy LT roi_high AND $
				radius LT diams(ic)/2.0) $
			, nobsfeat)
  print, ' feat/obsfeat = ', FLOAT(nfeature)/FLOAT(nobsfeat)
  print, ' obsfeat/obs = ', FLOAT(nobsfeat)/FLOAT(nobs)
  print, ''
  fraction_out(ic) = FLOAT(nfeature)/FLOAT(nobs)
  if ic EQ 0 then begin
    ; Add the other info
    fraction_out(n_elements(fraction_out)-1) = FLOAT(nobsfeat)/FLOAT(nobs)
    fraction_out(n_elements(fraction_out)-2) = FLOAT(nfeature)
    ; Return the Events in the nominal aperture
    ap_events = ml1(where(radius LT diams(ic)/2.0))
  end
end

print, fraction_out
print, ''


; Make a plot to show things... if desired...
if NOT(KEYWORD_SET(NO_PLOT)) then begin
; - - - - - - - - - - - - - - - - - - -
if KEYWORD_SET(PS_PLOT) then begin
; Plot to .ps

; color table:
; 0 - black
; 1 - white
; 2 - dark blue for background
; 3 to 254 go from red to white-blue
; 255 - white
dd_load_ct, /OTHER

!p.multi = [0,1,4]

; Colors for the events based on energy
ecolor_low = 0.2
energy_colors = 3+FIX(252.*((ml1.energy - ecolor_low)/(10.0-ecolor_low)))
energy_colors = ((energy_colors > 3) < 254 )

!p.multi = [0,1,2]
char_size = 1.2

; Most useful plot is in E vs X space:
;
plot, ml1.detx, ml1.energy, PSYM=3, /NODATA, $
	TITLE=trw_id+'/'+it_str+',  Feature is '+feature_str, $
	XTITLE='Facility Y (mm w.r.t. detector center)', $
	XRANGE=[-6.,6.], XSTYLE=1,$
	YTITLE='Photon Energy (keV)',  $
	CHARSIZE=char_size, BACK=1, COLOR=0
plots, ml1.detx, ml1.energy, PSYM=3, COLOR=energy_colors

; Plot solid lines to show the FEATURE:
; The feature is region around energy and within range of the
; order:
oplot, order_spacing*0.33*[-1.,1.], (energy+de)*[1.,1.]
oplot, order_spacing*0.33*[-1.,1.], (energy-de)*[1.,1.]

; OBSERVED events depend on spatial-spectral analysis method...
; For the nominally placed aperture with ROI:
; show ROI limits:
oplot, [-10.,10.], roi_low*[1.,1.], LINESTYLE=2
oplot, [-10.,10.], roi_high*[1.,1.], LINESTYLE=2
; show Aperture limits (in X):
oplot,  DYFAC+0.5*diameter*[1.,1.], [0.,10.], LINESTYLE=2
oplot, DYFAC-0.5*diameter*[1.,1.], [0.,10.], LINESTYLE=2

; Calculate and show the nominal fraction to go with the plot
; ALL FEATURE COUNTS
all_feature_counts = where(ABS(ml1.energy - energy) LT de AND $
	ABS(ml1.detx) LT 0.33*order_spacing, nfeature)
print, ''
print, '   Regarding the plot:'
print, 'All feature counts = ', nfeature

; ALL LINE-ORDER REGION COUNTS
all_counts = where(ABS(ml1.detx) LT 0.33*order_spacing, nallcounts)
print, 'All counts in 0.33*order_spacing = ', nallcounts

radius = SQRT((ml1.detx-DYFAC)^2+(ml1.dety-DZFAC)^2)
; OBSERVED COUNTS
observed_counts = where(ml1.energy GT roi_low AND $
				ml1.energy LT roi_high AND $
				radius LT diameter/2.0, nobs)
print, 'Observed counts  = ', nobs
; OBSERVED FEATURE COUNTS are all counts
; in the feature that are observed
obsfeat_counts = where(  (ABS(ml1.energy - energy) LT de AND $
	ABS(ml1.detx) LT 0.33*order_spacing)  AND  $
			(ml1.energy GT roi_low AND $
				ml1.energy LT roi_high AND $
				radius LT diameter/2.0) $
			, nobsfeat)
print, 'Observed feature counts  = ', nobsfeat

nom_fraction = FLOAT(nfeature)/FLOAT(nobs)
frac_str = ' Feature-fraction = '+ $
		STRING(nom_fraction,FORMAT='(F7.4)')
xyouts, -5.5, 0.5*(roi_low+energy-de), frac_str, CHARSIZE=1.2
xyouts, (0.33*order_spacing-0.75 < 3.0), 1.05*(energy+de), $
		STRING(nfeature,FORMAT='(I6)'), CHARSIZE=1.2
xyouts, DYFAC + 0.5*diameter - 0.1, 1.05*roi_low, STRING(nobs,FORMAT='(I6)'), $
		ALIGNMENT=1.0, CHARSIZE=1.2

; Now two little plots...
!p.multi = [2,2,2]

; Energy plot: detailed feature histogram... to show how the
; feature relates to the surrounding spectrum.  Does not show
; how the feature and aperture are related.
featehist = HISTOGRAM(ml1(all_feature_counts).energy - energy,$
	BINSIZE = de/10.0,MIN=-20.*de,MAX=20.*de)
allehist = HISTOGRAM(ml1(all_counts).energy - energy,$
	BINSIZE = de/10.0,MIN=-20.*de,MAX=20.*de)
histes = (findgen(n_elements(featehist))-0.5*FLOAT(n_elements(featehist)))*de/10.0
hist_color = 3+FIX(252.*((energy - ecolor_low)/(10.0-ecolor_low)))
pha_colors = 3+FIX(252.*((energy+histes - ecolor_low)/(10.0-ecolor_low)))
plot_io, histes, featehist, PSYM=10, YRANGE=[1.,MAX(featehist)], $
	XRANGE=[MIN(histes),MAX(histes)], XSTYLE=1, $
	CHARSIZE=char_size, TITLE='Feature in Spectrum', $
	XTITLE = 'dEnergy (keV)', YTITLE='Number of Events'
plots, histes, featehist, PSYM=1, COLOR = hist_color
plots, histes, allehist, PSYM=3, COLOR=pha_colors

; and a plot of the X,Y location of the events...
xyxr = (3.0*diameter*0.5 < 6.0)
xyyr = (3.0*diameter*0.5 < 6.0)
plot, ml1.detx, ml1.dety, PSYM=3, COLOR=0, /CLIP, $
	XRANGE= xyxr*[-1.,1.], XSTYLE=1, $
	YRANGE= xyyr*[-1.,1.], YSTYLE=1, /NODATA, $
	CHARSIZE = char_size, TITLE='Y-Z Plot of Events', $
	XTITLE='Facility Y (mm w.r.t. ...)', $
	YTITLE='Facility Z (mm w.r.t. detector center)'
; need to clip it ourselves?!?:
in_plot = where(ABS(ml1.detx) LT 3.0*diameter*0.5 AND $
		ABS(ml1.dety) LT 3.0*diameter*0.5)
plots, ml1(in_plot).detx, ml1(in_plot).dety, PSYM=3, $
	COLOR=energy_colors(in_plot)
; and the diameter of the aperture...
facet_circle, 0., 360., xcirc, ycirc
oplot, DYFAC+diameter*0.5*xcirc, DZFAC+diameter*0.5*ycirc, COLOR=0
oplot, [DYFAC], [DZFAC], PSYM=1, COLOR=0
xyouts, -0.8*xyxr, 0.8*xyyr, 'EE correction = ' + $
	STRING(FLOAT(nfeature)/FLOAT(nobsfeat),FORMAT='(F5.3)'), $
	CHARSIZE = 1.1

print, ''

!p.multi=0

end else begin
; GIF PLOTS

; Plot to .gif
set_plot, 'X'
gif_window = 0
window, gif_window, xsi=550, ysi=700

; color table:
; 0 - black
; 1 - white
; 2 - dark blue for background
; 3 to 254 go from red to white-blue
; 255 - white
dd_load_ct, /OTHER

!p.multi = [0,1,4]

; Colors for the events based on energy
ecolor_low = 0.2
energy_colors = 3+FIX(252.*((ml1.energy - ecolor_low)/(10.0-ecolor_low)))
energy_colors = ((energy_colors > 3) < 254 )

!p.multi = [0,1,2]
char_size = 1.2

; Most useful plot is in E vs X space:
;
plot, ml1.detx, ml1.energy, PSYM=3, /NODATA, $
	TITLE=trw_id+'/'+it_str+',  Feature is '+feature_str, $
	XTITLE='Facility Y (mm w.r.t. detector center)', $
	XRANGE=[-6.,6.], XSTYLE=1,$
	YTITLE='Photon Energy (keV)',  $
	CHARSIZE=char_size, BACK=2, COLOR=1
plots, ml1.detx, ml1.energy, PSYM=3, COLOR=energy_colors


; Plot solid lines to show the FEATURE:
; The feature is region around energy and within range of the
; order:
oplot, order_spacing*0.33*[-1.,1.], (energy+de)*[1.,1.]
oplot, order_spacing*0.33*[-1.,1.], (energy-de)*[1.,1.]

; OBSERVED events depend on spatial-spectral analysis method...
; For the nominally placed aperture with ROI:
; show ROI limits:
oplot, [-10.,10.], roi_low*[1.,1.], LINESTYLE=2
oplot, [-10.,10.], roi_high*[1.,1.], LINESTYLE=2
; show Aperture limits (in X):
oplot,  DYFAC+0.5*diameter*[1.,1.], [0.,10.], LINESTYLE=2
oplot, DYFAC-0.5*diameter*[1.,1.], [0.,10.], LINESTYLE=2

; Calculate and show the nominal fraction to go with the plot
; ALL FEATURE COUNTS
all_feature_counts = where(ABS(ml1.energy - energy) LT de AND $
	ABS(ml1.detx) LT 0.33*order_spacing, nfeature)
print, ''
print, '   Regarding the plot:'
print, 'All feature counts = ', nfeature

; ALL LINE-ORDER REGION COUNTS
all_counts = where(ABS(ml1.detx) LT 0.33*order_spacing, nallcounts)
print, 'All counts in 0.33*order_spacing = ', nallcounts

radius = SQRT((ml1.detx-DYFAC)^2+(ml1.dety-DZFAC)^2)
; OBSERVED COUNTS
observed_counts = where(ml1.energy GT roi_low AND $
				ml1.energy LT roi_high AND $
				radius LT diameter/2.0, nobs)
print, 'Observed counts  = ', nobs
; OBSERVED FEATURE COUNTS are all counts
; in the feature that are observed
obsfeat_counts = where(  (ABS(ml1.energy - energy) LT de AND $
	ABS(ml1.detx) LT 0.33*order_spacing)  AND  $
			(ml1.energy GT roi_low AND $
				ml1.energy LT roi_high AND $
				radius LT diameter/2.0) $
			, nobsfeat)
print, 'Observed feature counts  = ', nobsfeat

nom_fraction = FLOAT(nfeature)/FLOAT(nobs)
frac_str = ' Feature-fraction = '+ $
		STRING(nom_fraction,FORMAT='(F7.4)')
xyouts, -5.5, 0.5*(roi_low+energy-de), frac_str, CHARSIZE=1.2
xyouts, (0.33*order_spacing-0.75 < 3.0), 1.05*(energy+de), $
		STRING(nfeature,FORMAT='(I6)'), CHARSIZE=1.2
xyouts, DYFAC + 0.5*diameter - 0.1, 1.05*roi_low, STRING(nobs,FORMAT='(I6)'), $
		ALIGNMENT=1.0, CHARSIZE=1.2

; Now two little plots...
!p.multi = [2,2,2]

; Energy plot: detailed feature histogram... to show how the
; feature relates to the surrounding spectrum.  Does not show
; how the feature and aperture are related.
featehist = HISTOGRAM(ml1(all_feature_counts).energy - energy,$
	BINSIZE = de/10.0,MIN=-20.*de,MAX=20.*de)
allehist = HISTOGRAM(ml1(all_counts).energy - energy,$
	BINSIZE = de/10.0,MIN=-20.*de,MAX=20.*de)
histes = (findgen(n_elements(featehist))-0.5*FLOAT(n_elements(featehist)))*de/10.0
hist_color = 3+FIX(252.*((energy - ecolor_low)/(10.0-ecolor_low)))
pha_colors = 3+FIX(252.*((energy+histes - ecolor_low)/(10.0-ecolor_low)))
plot_io, histes, featehist, PSYM=10, YRANGE=[1.,MAX(featehist)], $
	XRANGE=[MIN(histes),MAX(histes)], XSTYLE=1, $
	CHARSIZE=char_size, TITLE='Feature in Spectrum', $
	XTITLE = 'dEnergy (keV)', YTITLE='Number of Events'
plots, histes, featehist, PSYM=1, COLOR = hist_color
plots, histes, allehist, PSYM=3, COLOR=pha_colors

; and a plot of the X,Y location of the events...
xyxr = (3.0*diameter*0.5 < 6.0)
xyyr = (3.0*diameter*0.5 < 6.0)
plot, ml1.detx, ml1.dety, PSYM=3, COLOR=1, /CLIP, $
	XRANGE= xyxr*[-1.,1.], XSTYLE=1, $
	YRANGE= xyyr*[-1.,1.], YSTYLE=1, /NODATA, $
	CHARSIZE = char_size, TITLE='Y-Z Plot of Events', $
	XTITLE='Facility Y (mm w.r.t. ...)', $
	YTITLE='Facility Z (mm w.r.t. detector center)'
; need to clip it ourselves?!?:
in_plot = where(ABS(ml1.detx) LT 3.0*diameter*0.5 AND $
		ABS(ml1.dety) LT 3.0*diameter*0.5)
plots, ml1(in_plot).detx, ml1(in_plot).dety, PSYM=3, $
	COLOR=energy_colors(in_plot)
; and the diameter of the aperture...
facet_circle, 0., 360., xcirc, ycirc
oplot, DYFAC+diameter*0.5*xcirc, DZFAC+diameter*0.5*ycirc, COLOR=1
oplot, [DYFAC], [DZFAC], PSYM=1, COLOR=1
xyouts, -0.8*xyxr, 0.8*xyyr, 'EE correction = ' + $
	STRING(FLOAT(nfeature)/FLOAT(nobsfeat),FORMAT='(F5.3)'), $
	CHARSIZE = 1.1
; Finish the .gif output
wshow, gif_window,iconic=0
image = tvrd()
tvlct, red, green, blue, /GET
; write it into the directory the ml1 file came from:
gif_file = marx_dir + '/' + trw_id + '/'+ it_str +'/'+ $
		trw_id + '_'+ it_str +'_ff.gif'
write_gif, gif_file, image, red, green, blue

print, ''

!p.multi=0
end ; of .gif plot
; - - - no_plot - - -
end

RETURN
END