PRO eae_effic_anal, eae, FRACTION=fraction, ffs
;
; This procedure forms the ratio of the measured
; effective area to the reference effective area
; and puts the result in the efficiency.
; 
; modified 2/23/98 dd Include the "feature_fractions.rdb" corrections
; if available.
;

if KEYWORD_SET(FRACTION) then begin
  ; Correct for the line/feature fractions if available...
  ; read in the "feature_fractions" if not passed
  if n_elements(ffs) LE 0 then begin
    ffs = rdb_read(!DDHETGCAL+'/cmp/eae/feature_fractions.rdb')
  end
  if n_elements(eae) NE n_elements(ffs) then begin
    print, ' * eae_effic_anal.pro: eae and ffs not same length'
    RETURN
  end
  ; Calculate the feature fraction errors
  dff_err = feat_frac_err(ffs)
end

for ir=0, n_elements(eae)-1 do begin
  if eae(ir).reference_value GT 0.0 then begin
    ; Calculate the simple efficiency ratio...
    eae(ir).effic = eae(ir).effective/eae(ir).reference_value
    ; Simple ratio, the errors are the same as for the effective
    ; area because we're ratioing grating/no-grating measurements
    ; and not correcting for the "line fractions"
    eae(ir).effic_method = 'simple ratio (eae_effic_anal.pro)'
    eae(ir).effic_err_plus = eae(ir).effective_err_plus/ $
				eae(ir).reference_value
    eae(ir).effic_err_minus = eae(ir).effective_err_minus/ $
				eae(ir).reference_value
    eae(ir).effic_syst_plus = eae(ir).effective_syst_plus/ $
				eae(ir).reference_value
    eae(ir).effic_syst_minus = eae(ir).effective_syst_minus/ $
				eae(ir).reference_value
    if KEYWORD_SET(FRACTION) then begin
      ; Correct for the line/feature fractions if available...
      if ffs(ir).fracNom GT 0.0 then begin
        ; OK, have this feature_fraction, what about the reference
        ; value feature fraction?        
        ; 
        ; Using logic similar to the eae_update_ref_vals.pro we'll
        ; calculate the "effective feature fraction" and its error:
        ;    ref_effective_ff
        ;    ref_effective_dff
        ;

        ; if there is an = then it is fixed value
        if STRPOS(eae(ir).reference_id,'=') GE 0 then begin
          ; If the reference value is fixed
          value = 0.0
          pieces = STR_SEP(eae(ir).reference_id, '=')
          READS, pieces(1), value
          eae(ir).reference_value = value
          ; nominal feature fractions and error...
          ref_effective_ff = 1.0
          ref_effective_dff = 0.0
          if eae(ir).trw_id EQ 'E-LXF-3D-22.044' then begin
            ; The reference for the LEG Ti-L is the expected
            ; HRMA area.  A no-grating simulation 
            ;  > xrcf_sim, 'E-LXF-3D-22.044', /NO_GRAT, /RUN, ITER=0, NUM=200000
            ;  > feat_fraction, eae(294),/NO_GRAT
            ; gives the following feature fraction:
            ref_effective_ff = 0.333
            ref_effective_dff = 0.25  ; Large error to flag this point
          end
          if eae(ir).trw_id EQ 'E-HXF-3D-22.047' then begin
            ; The reference here is the HRMA area based on
            ; assuming a known LEG 1st order efficiency
            ; the fraction corrections are for that case:
            ref_effective_ff = 0.5*(0.928+0.933)
            ref_effective_dff = 0.5*(0.0332+0.0339)
          end
        end else begin
          ; strip off any factor
          pieces = STR_SEP(eae(ir).reference_id, '*')
          ref_ids = pieces(0)
          ; second piece is the factor if it is there
          ; the factor will not effect the effective feature fraction
          ; but keep it all here for completeness
          factor = 1.0
          if n_elements(pieces) EQ 2 then begin
            READS, pieces(1), factor
          end
          ; See if it is a compound
          pieces = STR_SEP(ref_ids,'+')
          ; value is the summ of the areas
          ; and fvalue is the sum of (feat_frac times area)
          value = 0.0
          fvalue = 0.0
          ; for error, the dffvalue is sum of:
          ;  square of dff_err times fvalue
          dffvalue = 0.0
          f_ok = 1  ; set to zero if fraction is missing for a piece
          for ip=0,n_elements(pieces)-1 do begin
            match = where(eae.trw_id EQ pieces(ip), nmatch)
            if nmatch EQ 1 then begin
              value = value + factor*eae(match(0)).effective
              fvalue = fvalue + factor*ffs(match(0)).fracNom * $
			eae(match(0)).effective
              dffvalue = dffvalue + ( dff_err(match(0)) * $
			(factor*ffs(match(0)).fracNom * $
			eae(match(0)).effective) )^2
              if ffs(match(0)).fracNom LT 0. then begin
                f_ok = 0
                print, ' * eae_effic_anal: f_ok=0 for index ',match(0)
                RETURN
              end
            end else begin
              ; found zero or multiple times...
              print, 'eae_update_ref_vals: zero/multiple of:'+pieces(ip)
              if nmatch GT 1 then print, match
              RETURN
            end
          end ; loop over multiple ids
          ref_effective_ff = fvalue/value
          ref_effective_dff = SQRT(dffvalue)/fvalue
        end ; of non-fixed value reference

        ; Now we have ref_effective_ff and ref_effective_dff
        ; so we can correct the simple efficiency ratio value...
        corr_denom = ref_effective_ff
        dff_denom = ref_effective_dff
        corr_numer = ffs(ir).fracNom
        dff_numer = dff_err(ir)
        eae(ir).effic_method = 'fraction corrected, ' + $
		STRCOMPRESS(STRING(corr_numer, FORMAT='(F8.3)')+'/'+ $
		STRING(corr_denom,FORMAT='(F8.3)') ,/REMOVE_ALL)+ $
		' (eae_effic_anal.pro)'
        eae(ir).effic = eae(ir).effic * corr_numer/corr_denom
        ; The statistical error is still given by the
        ; error in the with-grating effective area
        eae(ir).effic_err_plus = eae(ir).effic_err_plus * $
				corr_numer/corr_denom
        eae(ir).effic_err_minus = eae(ir).effic_err_minus * $
				corr_numer/corr_denom
        ; The systematic error is now the combination of the
        ; errors in the fractions and does not depend on the
        ; previous error estimates.
        dff_total = SQRT(dff_numer^2 + dff_denom^2)
        ; The dff is a fractional error so the actual +/- errors are given
        ; by:
        eae(ir).effic_syst_plus = eae(ir).effic*(1.+dff_total) - $
							eae(ir).effic
        eae(ir).effic_syst_minus = ABS(eae(ir).effic/(1.+dff_total) - $
							eae(ir).effic)
      end ; of check for valid fraction for this one
    end ; FRACTION?
  end ; of check if ref value available
end ; of loop over eae

RETURN
END