PRO eae_sim, eae, SIM=sim, ANAL=anal, ORDER_CHECK=order_check, SELECT=select
;
; Simulate and analyze  the eae measurements from Phase 1
; Use these simulations to calculate "feature fractions"
; to correct the eae measurements...
;
; 4/21/99 dd Add fracZplus and fracZminus

; Parameters in code:
marx_sim_dir = !DDPHASE1SIM

; Read in the eae file if not supplied
if n_elements(eae) LE 0 then begin
  eae = rdb_read(!DDHETGCAL+'/cmp/eae/eae_weffics.rdb')
end

; Select those lines to do
; - - - - - - - - - - - -
; Initial simulation of all NONE, 0, and +/-1 order measurements
select1 = where( eae.source NE 'HIREF-W' AND $
	(eae.grating EQ 'NONE' OR ABS(eae.order) EQ 1 OR eae.order EQ 0)  )

; Note: in earlier version the -1 orders were not simulated because
; of symmetry reasons and to save time.  However because of the HXDS
; stage angle the orders may not be symmetric...  This time (April 99)
; simulate them all...
; - - -
; Many of the "3D" measurements have +1 and -1 so we can
; use feat_frac_trans to copy the +1 fractions into -1 (generally
; the simulations would be symmetric...)
; The "EE" measurements were done one at a time so better to 
; simulate these for -1 orders (rather than copy them by trying
; to match up with +1 orders...)
;;select2 = where( eae.source NE 'HIREF-W' AND $
;;	STRPOS(eae.trw_id,'-EE-') GE 0 AND eae.order EQ -1 , nselect2 )
; and include some 3D's that have only -1 order:
;;other_ids = ['D-HXS-3D-29.005','D-HXS-3D-29.003','D-HXF-3D-22.050']
;;for nid=0,n_elements(other_ids)-1 do begin
;;  add_this = where(eae.trw_id EQ other_ids(nid) AND eae.order EQ -1,nwhere)
;;  if nwhere EQ 1 then select2 = [select2,add_this]
;;end
;;if nselect2 GT 0 then begin
;;  select = [select1,select2]
;;end else begin
;;  select = select1
;;end

; - - - - - - - - - - - -
; Nominal simulations:
select = select1
want_valid_events = 40000L
; - - - - - - - - - - - -

; - - - - - - - - - - - -
; Test the s/w with self-contained Al-K set: 10-22
; Ultra10 running time: 9:50 - 10:25  or 2.7 minutes each.
; For 320 would take: 14.4 hours ...
;;select = select1(10:22)
;;want_valid_events = 40000L
; - - - - - - - - - - - -

; - - - - - - - - - - - -
; For some measurements we'll want more counts to reduce the statistical
; error due to the simulation itself...
; These were simulated to get 160K counts
; to get stat errors lower...
;;select = [67,69,73,75,91,154,157,158,267,269,349,351,354,356]
;;want_valid_events =  160000L
; - - - - - - - - - - - -

; Show the selected list
eae_list, eae(select)

; - - - - Simulation - - - - - 
; 6/22/99 dd Disable simulation to prevent overwriting
;            existing simulations (add "AND (0 EQ 1)" to below:)
if KEYWORD_SET(SIM) AND (0 EQ 1) then begin
  ; Do them...
  for is=0,n_elements(select)-1 do begin
    ie = select(is)
    ; simulate it once with few rays...
    it_num =0
    READS, eae(ie).iteration, it_num
    xrcf_sim, eae(ie).trw_id, /RUN, /POST, iter=it_num, $
	NUMRAYS=100000, N_POST_EVENTS = nout, MARX_SIM_DIR=marx_sim_dir
    ; nout is the number of events saved - want it to be want_valid_events
    if nout LT (3*want_valid_events/4) OR  $
		nout GT (5*want_valid_events/4) then begin
      print, ''
      print, ' * Repeat sim *  Nout was ',nout
      print, ''
      xrcf_sim, eae(ie).trw_id, /RUN, /POST, iter=it_num, $
	NUMRAYS=LONG(1.E5*want_valid_events/nout), MARX_SIM_DIR=marx_sim_dir
    end
    ; Now remove all the MARX output vectors to save disk space...
    it_dir = marx_sim_dir+'/'+eae(ie).trw_id + '/i' + eae(ie).iteration
    SPAWN, 'rm '+it_dir+'/marx_out/*'
  end
end ; simulation

; - - - - - Analysis - - - - -
if KEYWORD_SET(ANAL) then begin
  ; The result of analysis of the simulations is an rdb table
  ; giving the "feature_fractions" for different measurements
  ; for a method of analysis (e.g., ROI).  This feature_fraction
  ; is calculated by the routine feat_fraction.
  ; 
  ; This table is then used to correct the grating efficiency
  ; as calculated by eae_effic_anal.
  ;

  ; format for rdb table...
  ;  eaeIndex fracNom fracDplus fracDminus fracYplus fracYminus
  ;    fracZplus fracZminus Nfeature ObsFeatFrac
  ; 
  frac_struct = {eaeIndex:-1., fracNom:-1., $
		fracDplus:-1., fracDminus:-1., $
		fracYplus:-1., fracYminus:-1., $
		fracZplus:-1., fracZminus:-1., $
		Nfeature:-1., ObsFeatFrac:-1.}
  ffs = REPLICATE(frac_struct, n_elements(eae))
  ; fill the index values
  ffs.eaeIndex = INDGEN(n_elements(ffs))
  ; for the simulated ones, 
  ; carry out feat_fraction analysis
  ; and fill in the feature fraction information...
  for is=0,n_elements(select)-1 do begin
    ie = select(is)
    feat_fraction, eae(ie), FRACTION_OUT = fracs
    for itag = 0, n_elements(fracs)-1 do begin
      ffs(ie).(itag+1) = fracs(itag)
    end
  end

  ; Transfer the +1 order values to -1's ...
  ; Not any more!
  ;;feat_frac_trans, ffs, eae

  ; Now save the file
  rdb_write, 'feature_fractions.rdb', ffs

end ; analysis


; - - - - Order Check - - - - 
if KEYWORD_SET(ORDER_CHECK) then begin
  ; Check that the eae file order and the simulated
  ; detector location (in MARX/trw_id/iN/det_name_xyz.txt) agree.
  for is=0,n_elements(select)-1 do begin
    ie = select(is)
    it_dir = marx_sim_dir+'/'+eae(ie).trw_id + '/i' + eae(ie).iteration+'/'
    ; read the det_name_xyz.txt file
    OPENR, det_unit, it_dir+'det_name_xyz.txt', /GET_LUN
    sim_det_name = ''
    readf, det_unit, sim_det_name
    simx=0. & simy = 0. & simz = 0.
    readf, det_unit, simx, simy, simz
    CLOSE, det_unit
    FREE_LUN, det_unit

    grating = eae(ie).grating
    energy = eae(ie).energy

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 = 1000.0  ; mm everything in simulation
  end
  ELSE: begin
    print, '* eae_sim/ORDER_CHECK: unknown grating: '+grating
    RETURN
  end
ENDCASE
    ord_compare = -1.*eae(ie).order
    if ord_compare LT -900. then ord_compare = 0.
    flag = ''
    if NOT(approx_equal(ord_compare,simy/order_spacing,0.01)) then flag = '***'
    print, ie, ':  i'+eae(ie).iteration, eae(ie).order, simy/order_spacing, $
	flag
  end 
end ; order check

RETURN
END