PRO meg_misalign_sim, SIM=sim
; Commands, etc. to simulate MEG measurements...

; Need any commons?
; not yet

; Parameters common to analysis
reg_size = 0.25
bin_size = 0.010


;
; Simulate the HSI MEG+3 defocus set, i2 in particular...
;
if 0 EQ 1 then begin
; - - - - - - 
; MEG +3 order
; 
; Setup directory...
; done once:
; xrcf_sim, 'D-HXH-dF-16.002', MARX_SIM='/spectra/d8/MARX'

; D-HXH-dF-16.002i2 :
; QEIDEAL prevents application of the HSI QE to the result
; (more efficient ray-trace for PSF purposes...)
;;xrcf_sim, 'D-HXH-dF-16.002',/RUN, MARX_SIM='/spectra/d8/MARX', $
;;		N_POST=n_post, /POST, ITER=2, /QEIDEAL

;;num_to_trace = LONG((30000.)*(100000./213.))
;;xrcf_sim, 'D-HXH-dF-16.002',/RUN, MARX_SIM='/spectra/d8/MARX', $
;;		N_POST=n_post, /POST, ITER=2, /QEIDEAL, NUM=num_to_trace

; The result is available in /spectra/d8/MARX/D-HXH-dF-16.002/i2/ml1.idlsav
restore, '/spectra/d8/MARX/D-HXH-dF-16.002/i2/ml1.idlsav'
print, tag_names(ml1)
print, n_elements(ml1)

; This rotation can be useful...  But note that
; the 1D cross-dispersion measurements were done with
; a horizontal slit!
; rotate to grating axes
; MEG is 4.74 degrees
;;rot_angle = -1.0*4.74*!DTOR
;;grat_y = ml1.detx * COS(rot_angle) - ml1.dety * SIN(rot_angle)
;;grat_z = ml1.detx * SIN(rot_angle) + ml1.dety * COS(rot_angle)
; or don't rotate...
grat_y = ml1.detx
grat_z = ml1.dety

; select the events
sel = where( (ABS(grat_y) LE reg_size) AND $
		(ABS(grat_z) LE reg_size) )
grat_y = grat_y(sel)
grat_z = grat_z(sel)

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

plot, grat_y, grat_z, PSYM=3, $
	XRANGE=reg_size*[-1.,1.], XSTYLE=1, $
	YRANGE=reg_size*[-1.,1.], YSTYLE=1, $
	TITLE = 'MARX simulation'

; make the histogram along Z
lin_hist, grat_z, bin_size, hist_axis, hist_counts

plot_io, hist_axis, hist_counts, PSYM=10, $
	XRANGE=reg_size*[-1.,1.], XSTYLE=1, $
	YRANGE=[0.5, 6000.0], YSTYLE=1, $
	TITLE = 'MARX simulation'

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


;
; Simulate the three PSF/1D scans
; (these simulations can be analyzed to produce
;  any of the three types of slit scans...)
;
; MEGs
trw_ids = ['D-HXF-P1-19.016','D-HXF-P1-19.017a','D-HXF-P1-19.018']
postrays = 8.*[30000.,30000.,10000.]
raysper = [3754.,1918.,152.]
;
; Do the HEG too:
trw_ids = ['D-HXF-P1-19.013','D-HXF-P1-19.014','D-HXF-P1-19.015']
postrays = 1.0*[30000.,30000.,10000.]
raysper = [2700.,700.,40.]
; 
; and LEGs for fun...
trw_ids = ['D-LXF-P1-17.007','D-LXF-P1-17.008','D-LXF-P1-17.009']
postrays = 0.2*[30000.,30000.,10000.]
raysper = [6000.,850.,16.]

;
if KEYWORD_SET(SIM) then begin
; - - - - - - 
; 
; Setup directory...
; done once:
; trw_id = 'D-HXF-P1-19.016'
; trw_id = 'D-HXF-P1-19.017a'
; trw_id = 'D-HXF-P1-19.018'
; xrcf_sim, trw_id, MARX_SIM='/spectra/d8/MARX'

; Find out how many rays come out for 100K traced...
; QEIDEAL prevents application of the HSI/FPC QE to the result
; (more efficient ray-trace for PSF purposes...)
; xrcf_sim, trw_id,/RUN, MARX_SIM='/spectra/d8/MARX', $
;		N_POST=n_post, /POST, ITER=0, /QEIDEAL
; print, n_post  ; '16: 3754/100K, '17a: 1918/100K, 18: 152/100K

; Ray trace 'em all
for it = 0, 2 do begin
  trw_id = trw_ids(it)
  raysper100k = raysper(it)
  num_to_trace = LONG((postrays(it))*(100000./raysper100k))
  xrcf_sim, trw_id, /RUN, MARX_SIM='/spectra/d8/MARX', $
  		N_POST=n_post, /POST, ITER=0, /QEIDEAL, NUM=num_to_trace
end
; - - - - - - 
end ; of simulation
;
; and Analyze them...
;
if 1 EQ 1 then begin
; - - - - - - 
; The results are available in /spectra/d8/MARX/'+trw_id+'/i0/ml1.idlsav

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

for it = 0, 2 do begin
  trw_id = trw_ids(it)

  evt_file = '/spectra/d8/MARX/'+trw_id+'/i0/ml1.idlsav'
  print, ''
  print, trw_id, ' : '
  print, ''

  hist_from_events, evt_file, 'temp_hist.rdb', AXIS_ANGLE= 90.0, $
		reg_size=reg_size

  hist_play, 'temp_hist.rdb', /NORM, /LOG, /CENTER, /GAUSS, $
		XRANGE=[-0.25,0.25], YRANGE=[0.0001,0.5]

;  wait, 5

end

!p.multi=0

; - - - - - - 
end

RETURN
END