FUNCTION feat_frac_err, ffs, eae, PLOT=plot, STAT_REL = stat_rel
; Calculate and return the error for the feature-fraction
;
; 4/21/99 dd Added in Z plus/minus error and slight change
;            to equation for aperture errors, "dia_dff".

; Parameter:
; how accurate are the correction factors?
frac_syst_error = 0.10  ; 10% of deviation from 1.0

; Calculate the feature_fraction error as a combination
; of:

; SIMULATION STATISTICS
;  - stat error based on number of counts in simulation
stat_dff = SQRT(1.+(ffs.fracNom > 0.))/SQRT(ffs.nfeature > 1.)

; SYSTEMATICS of CORRECTIONS
; - correction due to fraction observed (contamination):
syst_obs_dff = frac_syst_error*(EXP(ABS(ALOG((ffs.ObsFeatFrac > 0.1)))) - 1.0)

; - correction due to fraction not observed (encircled-energy correction):
syst_ee_dff = frac_syst_error*(EXP(ABS(ALOG( $
	(ffs.fracNom/ffs.ObsFeatFrac > 0.1) ))) - 1.0)

; total systematic
syst_dff = SQRT( (syst_obs_dff)^2 + (syst_ee_dff)^2 )

; APERTURE 
;  - correction for diameter/location variation
;; old equation
;;dia_dff =  SQRT( ABS(ffs.fracDminus-ffs.fracDplus)^2 + $
;;		ABS(ffs.fracYminus-ffs.fracYplus)^2 + $
;;		ABS(ffs.fracZminus-ffs.fracZplus)^2 )/ffs.fracNom
; better estimate
; for backward compatibility don't require the fracZ's:
ffs_tags = tag_names(ffs)
if TOTAL(ffs_tags EQ 'FRACZPLUS') GE 1 then begin
  dia_dff =  SQRT( ( (ffs.fracDminus-ffs.fracNom)^2 > $
		(ffs.fracDplus-ffs.fracNom)^2 ) + $
		 ( (ffs.fracYminus-ffs.fracNom)^2 > $
		(ffs.fracYplus-ffs.fracNom)^2 ) + $
		 ( (ffs.fracZminus-ffs.fracNom)^2 > $
		(ffs.fracZplus-ffs.fracNom)^2 ) )/ffs.fracNom
end else begin
  print, ''
  print, ' Old ffs file - no Z information...'
  dia_dff =  SQRT( ( (ffs.fracDminus-ffs.fracNom)^2 > $
		(ffs.fracDplus-ffs.fracNom)^2 ) + $
		 ( (ffs.fracYminus-ffs.fracNom)^2 > $
		(ffs.fracYplus-ffs.fracNom)^2 ) )/ffs.fracNom
end

; Total error (as df/f):
dff_err = SQRT(stat_dff^2 + syst_dff^2 + dia_dff^2)

; Plot the error and contributions...
if KEYWORD_SET(PLOT) then begin
  !p.multi = [0,1,2]
  valid = where(ffs.fracNom GT 0.)
  plot_oo, eae(valid).energy, stat_dff(valid), PSYM=1, $
	YRANGE=[0.0001,10.0], XRANGE=[0.1,10.], $
	XTITLE='Test Energy (keV)', $
	YTITLE='df/f Error', $
	TITLE='dF/F Error Contributions: +=Stat, trangle=Syst, square=Dia,Y,Z'
  oplot, eae(valid).energy, syst_dff(valid), PSYM=5
  oplot, eae(valid).energy, dia_dff(valid), PSYM=6

  plot_oo, eae(valid).energy, dff_err(valid), PSYM=4, $
	YRANGE=[0.001,10.0], XRANGE=[0.1,10.0], $
	XTITLE='Test Energy (keV)', $
	YTITLE='df/f Error', $
	TITLE='Total dF/F Error'
  !p.multi=0

  ; Flag where stat error is relevant - total error is over 2%
  ; and stat error is > 0.7 times this.
  stat_rel = where(stat_dff(valid)/dff_err(valid) GE 0.7 AND $
		dff_err(valid) GT 0.02, nval)
  if nval GE 1 then begin
    print, ''
    print, ' Found some tests where the statistical (simulation) error'
    print, ' is significant: '
    print, ''
    stat_rel = valid(stat_rel)
    ; List these measurements:
    eae_list, eae(stat_rel)
    ; and the stat and total errors;  Counts Factor gives the
    ; factor increase in simulation counts needed to have
    ; stat error equal the other errors...
    print, ' index       TRW_ID         Stat Error   Other Error  Total Error',$
		'  Counts Factor'
    for il=0,n_elements(stat_rel)-1 do begin
      print, STRING(stat_rel(il),FORMAT='(I6)'), $
	'  ', eae(stat_rel(il)).trw_id, $
		stat_dff(stat_rel(il)), $
		SQRT(dff_err(stat_rel(il))^2 - $
		stat_dff(stat_rel(il))^2 ), $
		dff_err(stat_rel(il)), $
		(stat_dff(stat_rel(il))/SQRT(dff_err(stat_rel(il))^2 - $
		stat_dff(stat_rel(il))^2 ) )^2
    end
  end
end

RETURN, dff_err
END