PRO hsi_source_spectra, trw_id, LPR = lpr, NO_ANAL=no_anal

@cmdb_common
@labx_common
COMMON internal, trw_run_ids

; Use TRW ID to derrive all other info...
if n_elements(trw_run_ids) LT 2 then begin
  trw_run_ids = rdb_read('/spectra/d6/Runids/trw_mod.rdb')
end

; Setup where the data come from
data_prefix = '/spectra/d6/'

; Assume CMDB is available
if cmdb_using_cal GT 0 then begin
@cmdb_cal_integers
end else begin
@cmdb_trw_integers
end

; Get XRCF parameters
xrcf_params, ANGLES=angles, PERIODS=periods, HC=hc, $ 
	ROWLAND=rowland, NAMES=names

; Spectral extraction width
spec_width = 2.0  ; mm

; Where to keep the resulting files
save_em_dir = '/spectra/d6/Source_spec/'

print, '..............................'
print, '  TRW ID  : '+trw_id

; Parameters needed for analysis
this_meas = where(trw_run_ids.trwid EQ trw_id AND $
	STRPOS(trw_run_ids.irig_start,'NULL') LT 0, nfound)
if nfound GE 1 then begin
  run_id = STRCOMPRESS(STRING( $
	LONG(trw_run_ids(this_meas(0)).runid) ),/REMOVE_ALL)
  irig_str = trw_run_ids(this_meas(0)).irig_start
  pieces = STR_SEP(irig_str,' ')
  date_pieces = STR_SEP(pieces(0),'/')
  if n_elements(date_pieces) EQ 3 then begin
    d_m = date_pieces(0)
    if STRLEN(d_m) EQ 1 then d_m = '0'+date_pieces(0)
    d_d = date_pieces(1)
    if STRLEN(d_d) EQ 1 then d_d = '0'+date_pieces(1)
    d_y = STRMID(date_pieces(2),2,2)
    date_str = d_y + d_m + d_d
  end else begin
    print, 'For '+trwid+ $
		" didn't get 3 pieces from irig date: "+irig_str
    date_str = 'NULL'
  end
  date_code = date_str
end else begin
  print, 'TRW ID '+trw_id+' not found in trw_mod.rdb'
  RETURN
end

print, 'date code : '+date_code
print, '   run id : '+run_id

; Use CMDB to find out the grating type, source, and energy
this_meas = where(cmdb_fields(*,c_id) EQ trw_id, nfound)
grating = -1
if nfound EQ 1 then begin
  if cmdb_fields(this_meas(0),c_grating) EQ 'LETG' then grating = 0
  if cmdb_fields(this_meas(0),c_grating) EQ 'HETG' then begin
    shutter_str = cmdb_fields(this_meas(0),c_shutters)
    if STRPOS(shutter_str,'MEG') GE 0 then grating = 1
    if STRPOS(shutter_str,'HEG') GE 0 then grating = 2
    if STRPOS(shutter_str,'ALL,ALL') GE 0 then grating = -1
  end
  ; get a nice source string
  source_str = cmdb_source_str(this_meas(0))
  print, '   source : '+source_str
  ; and the CMDB line energy
  src_energy = cmdb_data(this_meas(0)).energy
  print, '   energy : '+STRING(src_energy,FORMAT='(F8.4)')+' keV'
end else begin
  print, trw_id+' found no or multiple times in cmdb :'+cmdb_file
  RETURN
end

; Find the HSI location in the focal plane
if STRPOS(trw_id, '-3D-') GT 0 OR STRPOS(trw_id, '-dF-') GT 0 then begin
  ; Read the loc file to know where the HSI images are
  xrcf_loc_file_read, trw_id, Xloc, Yloc, Zloc
end else begin
  ; Use the CMDB location
  ir = this_meas(0)
  if STRPOS(cmdb_fields(ir,c_measurement_config),'HXDA') $
	GE 0 then begin
    ; Using HXDA
    Y_str = cmdb_fields(ir,c_HXDA_Y)
    Z_str = cmdb_fields(ir,c_HXDA_Z)
    X_str = cmdb_fields(ir,c_HXDA_defocus)
  end else begin
    ; Using FAM etc.
    Y_str = cmdb_fields(ir,c_five_axis_y)
    Z_str = cmdb_fields(ir,c_five_axis_Z)
    X_str = cmdb_fields(ir,c_five_axis_X)
  end
  temp = 0.0
  READS, X_str, temp
  Xloc = [temp]
  READS, Y_str, temp
  Yloc = [temp]
  READS, Z_str, temp
  Zloc = [temp]
end
; Some kludges
; For D-HXH-3D-11.043 the grating is MEG and it is offset correctly
; .loc file was calculated "wrong" because of "ALL,ALL" shutters.
if trw_id EQ 'D-HXH-3D-11.043' then begin
  grating = 1
  Zloc = Yloc * TAN(angles(grating)*!DTOR)
end
if trw_id EQ 'D-HXH-EA-13.012' then begin
  ; Both MEG and HEG first orders are visible
  grating = 2
end

; Use location to help uncover grating being used
; if ALL,ALL
if grating LT 0 then begin
  if cmdb_fields(this_meas(0),c_grating) EQ 'HETG' then begin
    if ABS(Yloc(0)) GT 0.5 then begin
      if TOTAL(Zloc/Yloc GT 0.05) GT 0 then begin
        ; MEG
        grating = 1
      end
      if TOTAL(Zloc/Yloc LT -0.05) GT 0 then begin
        ; HEG
        grating = 2
      end
    end
  end  
end

; print, location(s)
for iloc = 0, n_elements(Xloc)-1 do begin
print, '    Y,Z,X : '+STRING(Yloc(iloc)) + $
	STRING(Zloc(iloc))+STRING(Xloc(iloc))
end

if grating LT 0 then begin
  print, 'TRW ID '+trw_id+' not well-defined grating type'
  RETURN
end else begin
  print, '  grating : '+names(grating)
  ; Calculate the order based on the location, energy, and 
  ; grating
  disp_first = rowland*(hc/src_energy)/periods(grating)
  sign = FIX(1.001*Yloc(0)/(ABS(Yloc(0)) > 1.0) )  ; gives -1,0,+1
  order = sign * FIX((ABS(Yloc(0))+disp_first/2.0 )/disp_first)
  ; if the order is 0 then use 1 (probably a high-energy
  ; line so image was zero-centered...)
  if order EQ 0 then order = 1  
  print, '    order : '+ STRING(order,FORMAT='(I3)')
end

; The iteration is generally 0 but in case it is not...
iter = 0
; Special cases...
if trw_id EQ 'D-LXH-dF-16.001' then iter = 2
if trw_id EQ 'D-HXH-dF-16.002' then iter = 2
if trw_id EQ 'D-HXH-dF-16.003' then iter = 2
print, '     iter : '+ STRING(iter,FORMAT='(I3)')

; Done with preliminaries...
if KEYWORD_SET(NO_ANAL) then RETURN

; Do only one iteration
it_id = STRCOMPRESS(STRING(iter),/REMOVE_ALL)
; Read the data file
file_name = data_prefix+date_code+'/hsi'+run_id+'i'+it_id+'.fits'
hsi_data = xrcf_hsi_read(file_name,/FAC)

; Single iteration
glob_Y = hsi_data.Y + 1000.*Yloc(iter)
glob_Z = hsi_data.Z + 1000.*Zloc(iter)

; Create wavelengths of each photon
disp_dist = (ABS(glob_Y)/1000.0)/COS(!DTOR*angles(grating))
wavelengths = ((periods(grating)/ABS(FLOAT(order)))*(disp_dist/rowland))


if KEYWORD_SET(LPR) then begin
  pre_print_portrait
end

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

; plot the events
plot, glob_Y, glob_Z, PSYM=3, $
	XRANGE = 1000.0*Yloc(iter)+[-9000.,9000.], $
	YRANGE = 1000.0*Zloc(iter)+[-9000.,9000.], $
	TITLE = trw_id+': '+source_str+',  '+names(grating)+ $
	'['+STRCOMPRESS(STRING(order),/REMOVE_ALL)+'] ('+ $
		date_code+'/hsi'+run_id+'i'+ $
		STRCOMPRESS(STRING(iter),/REMOVE_ALL)+'.fits)', $
	XTITLE = 'Facility Y (um)', YTITLE = 'Facility Z (um)', $
	CHARSIZE = char_size

; over plot extraction region for spectrum and background
ys = 1000.0*Yloc(iter)+[-9000.,9000.]
zs = ys*TAN(angles(grating)*!DTOR)
oplot, ys, zs - 1000.0*spec_width/2.0
oplot, ys, zs - 1000.0*spec_width
oplot, ys, zs + 1000.0*spec_width/2.0
oplot, ys, zs + 1000.0*spec_width

; select spectrum and background photons
expected_z = glob_Y * TAN(angles(grating)*!DTOR)
spec_photons = where(glob_Z GT expected_z - 1000.0*spec_width/2.0 AND $
		glob_Z LT expected_z + 1000.0*spec_width/2.0)
bkg_photons = where( (glob_Z GT expected_z - 1000.0*3.*spec_width/2.0 AND $
		glob_Z LT expected_z - 1000.0*spec_width/2.0) OR $
		(glob_Z GT expected_z + 1000.0*spec_width/2.0 AND $
		glob_Z LT expected_z + 1000.0*3.*spec_width/2.0) )

; plot the selected events
plot, glob_Y(spec_photons), glob_Z(spec_photons), PSYM=3, $
	XRANGE = 1000.0*Yloc(iter)+[-9000.,9000.], $
	TITLE = 'Events in Spectral-extraction Region,' + $
	STRING(spec_width,FORMAT='(F5.2)')+' mm wide', $
	XTITLE = 'Facility Y (um)', YTITLE = 'Facility Z (um)', $
	CHARSIZE = char_size

; plot the background events
plot, glob_Y(bkg_photons), glob_Z(bkg_photons), PSYM=3, $
	XRANGE = 1000.0*Yloc(iter)+[-9000.,9000.], $
	TITLE = 'Events in Background-extraction Regions', $
	XTITLE = 'Facility Y (um)', YTITLE = 'Facility Z (um)', $
	CHARSIZE = char_size

; Set the wavelength scale...
disp_ave = ABS(Yloc(iter))
disp_max = disp_ave + 9.0
disp_min = (disp_ave - 9.0) > 1.0
lam_max = ((periods(grating)/ABS(FLOAT(order)))*(disp_max/rowland))
lam_min = ( ((periods(grating)/ABS(FLOAT(order)))*(disp_min/rowland)) > 1.24 )
dlam = (lam_max - lam_min)/500.0

lam_spec_hist = histogram(wavelengths(spec_photons), $
		BINSIZE=dlam, MIN=0., MAX=lam_max)
lam_bkg_hist = histogram(wavelengths(bkg_photons), $
		BINSIZE=dlam, MIN=0., MAX=lam_max)
lams = lam_max*(INDGEN(n_elements(lam_spec_hist))+0.5)/ $
			FLOAT(n_elements(lam_spec_hist))

if KEYWORD_SET(LPR) then begin
  device,/close
  set_plot, 'X'
  SPAWN, 'mv idl.ps '+trw_id+'_image.ps'  
  SPAWN, 'ps2gif '+trw_id+'_image.ps'  
  SPAWN, 'cp '+trw_id+'_image.gif '+save_em_dir+trw_id+'_image.gif'
  SPAWN, 'rm '+trw_id+'_image.*'  
  pre_print_portrait
end

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

plot_io, lams, lam_spec_hist, PSYM=10, $
	TITLE = trw_id+': '+source_str+',  '+names(grating)+ $
		'['+STRCOMPRESS(STRING(order),/REMOVE_ALL)+']', $
		XTITLE = 'Wavelength (A)', YTITLE = 'counts/bin '+ $
		'(bin size = '+STRING(dlam,FORMAT='(F7.4)')+' A)', $
		YRANGE=[1.,MAX(lam_spec_hist(where(lams GE 0.5)))], $
		XRANGE = [lam_min,lam_max], $
		CHARSIZE = char_size
oplot, lams, lam_bkg_hist, PSYM=1
; Add indication of CMDB assumed energy
oplot, hc/src_energy*[1.,1.], [1.,1.E6], LINESTYLE=2

; Make a plot vs Energy
spec_es = hc/lams

; Convert the spectrum to units of counts/keV
spec_per_kev = ((lams^2)/hc) * (lam_spec_hist - lam_bkg_hist)/dlam

plot_oo, spec_es, spec_per_kev, PSYM=10, $
	TITLE = '('+ $
		date_code+'/hsi'+run_id+'i'+ $
		STRCOMPRESS(STRING(iter),/REMOVE_ALL)+'.fits)', $
		XTITLE = 'Energy (keV)', YTITLE = 'counts/keV '+ $
		'(bin size = '+STRING(dlam,FORMAT='(F7.4)')+' A)', $
		YRANGE=[( MIN(spec_per_kev(where(spec_per_kev GT 0.))) > 1.), $
			MAX(spec_per_kev(where(spec_es LT 10.0)))], $
		XRANGE = [0.1,10.], $
		CHARSIZE = char_size
; Add indication of CMDB assumed energy
oplot, src_energy*[1.,1.], [1.E-3,1.E9], LINESTYLE=2

; Output a file!
; RDB format
TAB = STRING(9B)
OPENW, out_unit, save_em_dir+trw_id+'_spec.rdb',/GET_LUN
fmt = '(F12.4)'
printf, out_unit, '#filename:'+TAB+trw_id+'_spec.rdb'
printf, out_unit, '#fileCreationDate:'+TAB+SYSTIME()
printf, out_unit, '#swVersion:'+TAB+'hsi_source_spectra.pro'
printf, out_unit, '#trw_id:'+TAB+trw_id
printf, out_unit, '#date_code:'+TAB+date_code
printf, out_unit, '#run_id:'+TAB+run_id
printf, out_unit, '#iteration:'+TAB+STRCOMPRESS(STRING(iter),/REMOVE_ALL)
printf, out_unit, '#source_string:'+TAB+source_str
printf, out_unit, '#energy:'+TAB+STRING(src_energy,FORMAT='(F10.6)')
printf, out_unit, '#X_hsi:'+TAB+STRING(Xloc(iter),FORMAT='(F10.6)')
printf, out_unit, '#Y_hsi:'+TAB+STRING(Yloc(iter),FORMAT='(F10.6)')
printf, out_unit, '#Z_hsi:'+TAB+STRING(Zloc(iter),FORMAT='(F10.6)')
printf, out_unit, '#grating:'+TAB+names(grating)
printf, out_unit, '#order:'+TAB+STRING(order,FORMAT='(I3)')
printf, out_unit, '#spec_width:'+TAB+STRING(spec_width,FORMAT='(F10.6)')
printf, out_unit, '#'
printf, out_unit, '#    Column definitions for humans:'
printf, out_unit, '# wavelen    wavelength of bin in Angstroms'
printf, out_unit, '#            (bins are all the same wavelength width)'
printf, out_unit, '# counts     counts in bin in extraction region'
printf, out_unit, '# bkgcounts  counts in bin in background regions'
printf, out_unit, '# energy     energy of the bin in keV'
printf, out_unit, '#'
printf, out_unit, 'wavelen'+TAB+'counts'+TAB+'bkgcounts'+TAB+'energy'
printf, out_unit, 'N'+TAB+'N'+TAB+'N'+TAB+'N'
for io=FIX(lam_min/dlam -1),n_elements(lams)-1 do begin
  printf, out_unit, STRING(lams(io),FORMAT=fmt) + TAB  + $
	STRING(lam_spec_hist(io),FORMAT=fmt)+ TAB + $
	STRING(lam_bkg_hist(io),FORMAT=fmt)+ TAB + $
	STRING(spec_es(io),FORMAT=fmt)

end
close, out_unit
free_lun, out_unit

if KEYWORD_SET(LPR) then begin
  device,/close
  set_plot, 'X'
  SPAWN, 'mv idl.ps '+trw_id+'_spec.ps'  
  SPAWN, 'ps2gif '+trw_id+'_spec.ps'  
  SPAWN, 'cp '+trw_id+'_spec.gif '+save_em_dir+trw_id+'_spec.gif'
  SPAWN, 'rm '+trw_id+'_spec.ps'  
  SPAWN, 'rm '+trw_id+'_spec.gif'  
end

!p.multi = 0

RETURN
END