PRO hsi_source_spectra, trw_id, LPR = lpr, NO_ANAL=no_anal, $
	COMPARE=compare, PS_COMPARE = ps_compare, HEAE_ANAL=heae_anal
;+
; Analyze HSI images to extract source spectra
;   /LPR creates output plot files
;   /NO_ANAL does no processing: useful to see the various parameters
;            of the measurement printed out.
;   /COMPARE converts counts/keV to photons/cm^2 keV and compares
;            result with labxray source model.
;-
; 12/15/97 dd Added /COMPARE: conversion of counts/keV spectra to
;             photons/cm^2 keV using HRMA-Grating-HSI CIPs
; 4/8/99   dd Update the HRMA, Grating, and HSI QE values from CALDBcal.
; 4/13/99  dd Do a QE depression correction when reading in the HSI data,
;             burried below are the lines:
;        qed_value = hsi_qe_depress(run_id_numerical)
;        hsi_data = xrcf_hsi_read(file_name,/FAC,EXPOSURE=hsi_time, $
;		RAW_HSI= raw_hsi, $
;		QE_DEPRESS=qed_value, QE_INFO_STR=qe_info_str)
;             And also make "tile plots" for central HSI tile, using:
;        hsi_tile_plots, raw_hsi, plot_title, QE_INFO_STR = qe_info_str
; 4/23/99 dd Add the HEAE_ANAL option to add counts info to the
;            heae_whcounts.rdb file.
;  - The /HEAE_ANAL keyword to hsi_source_spectra will
;    have it read heae_whcounts.rdb, modify the appropriate line,
;    by adding focal_counts, focal_err_plus, focal_err_minus,
;    focal_time, focal_analysis=ROI,ALL , etc.
;    and then write it back out.
; 5/18,31/99 dd Use new HRMA 'N0004 QEs and HSI QE ('990530).
; - - - - - - - - - - - - - 
@cmdb_common
@labx_common
; Make a common for big, long-to-read-in things
COMMON internal, trw_run_ids, cip_hrma, cip_heg, cip_meg, cip_leg, cip_hsi

; - - - - - - - - - - - - - - - - - - 
; Preparation work...

; Use TRW ID to derrive all other info...
; lookup table to map TRW-ID to runid and date:
if n_elements(trw_run_ids) LT 2 then begin
  trw_run_ids = rdb_read(!DDIDL+'/xrcf/eae_trw_mod.rdb')
end

; Get cal products from CALDB cal areas
; Load them if not already
if (KEYWORD_SET(COMPARE) AND n_elements(cip_hrma) LT 10) then begin
  print, 'Loading CIPs...'
  print, '  HSI .'
  cip_hsi = rdb_read(!DDHXDSDIR+'/hsi_qe_N19990530.rdb')
  ; Use the HRMA N0003 XRCF value to remove artificial 2 keV jump...
  ; Use the HRMA N0004 XRCF values (artificial 0.94 keV jump)
  print, '  HRMA at XRCF ...'
  cip_hrma = rdb_read(!DDHRMACAL+'/cip/hrmaD1996-11-01XeffareaN0004.rdb')
  print, '  MEG .....'
  cip_meg = rdb_read(!DDHETGCAL+'/cip/hetgmegD1996-11-01greffpr005N0002.rdb')
  print, '  HEG .....'
  cip_heg = rdb_read(!DDHETGCAL+'/cip/hetghegD1996-11-01greffpr005N0002.rdb')
  print, '  LETG ..........'
  cip_leg = rdb_read(!DDLETGCAL+'/cip/letgD1996-11-01greffpr005N0002.rdb')
endif

; Setup where the data come from
data_prefix = !DDHXDSDATA+'/phase1/data/'

; 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, $ 
	ROWLANDS=rowlands, NAMES=names

; Spectral extraction width
spec_width = 2.0  ; mm

; Where to keep the resulting files
; Not really something to fill up CALDB area with...
save_em_dir = !DDHSISPECTRA


; - - - - - - - - - - - - - - - - - - 
; Get and Report Parameters needed for analysis

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

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)
  run_id_numerical = trw_run_ids(this_meas(0)).runid
  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 xrcf/eae_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(44) the grating is MEG(HEG) and it is offset correctly
; .loc file was calculated at zero angle because of "ALL,ALL" shutters.
; So the Yloc is really the dispersion hypooneus...
if trw_id EQ 'D-HXH-3D-11.043' then begin
  grating = 1
  Zloc = Yloc * SIN(angles(grating)*!DTOR)
  Yloc = Yloc * COS(angles(grating)*!DTOR)
end
if trw_id EQ 'D-HXH-3D-11.044' then begin
  grating = 2
  Zloc = Yloc * SIN(angles(grating)*!DTOR)
  Yloc = Yloc * COS(angles(grating)*!DTOR)
end
if trw_id EQ 'D-HXH-3D-11.037' then begin
  grating = 1
  Zloc = Yloc * SIN(angles(grating)*!DTOR)
  Yloc = Yloc * COS(angles(grating)*!DTOR)
end
if trw_id EQ 'D-HXH-3D-11.038' then begin
  grating = 2
  Zloc = Yloc * SIN(angles(grating)*!DTOR)
  Yloc = Yloc * COS(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
if trw_id EQ 'D-LXH-EA-13.011' then begin
  ; The Y location is slightly off resulting in a split
  ; spectral peak... adjust for aligned 3rd order:
  Yloc = 0.05  ; mm
  ; (the neg first order is on HSI gap it appears)
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 = rowlands(grating)*(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
; Special order selection
if trw_id EQ 'E-HXH-3D-10.008a' then order = 1
if trw_id EQ 'D-LXH-3D-9.006' then order = 5

; 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'

; Do QE depression correction when reading in the HSI data:
qed_value = hsi_qe_depress(run_id_numerical)
hsi_data = xrcf_hsi_read(file_name,/FAC,EXPOSURE=hsi_time, $
		RAW_HSI= raw_hsi, $
		QE_DEPRESS=qed_value, QE_INFO_STR=qe_info_str)

; EXPOSURE keyword is the actual exposure time in seconds from
; the FITS header
print, 'hsi_time (sec) = ', hsi_time

; Create "global" coordinates w.r.t. zero-order
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/rowlands(grating)))


; Plot to .gif
set_plot, 'X'
window, 0, xsi=500, ysi=700

dd_load_ct, /OTHER


; Make the hsi "tile" plots here
; use the raw HSI data

plot_title = trw_id + $
	' ('+date_code+'/hsi'+run_id+'i'+it_id+'.fits) '
hsi_tile_plots, raw_hsi, plot_title, QE_INFO_STR = qe_info_str

if KEYWORD_SET(LPR) then begin
  ; Finish the .gif output
  wshow,0,iconic=0
  image = tvrd()
  tvlct, red, green, blue, /GET
  gif_file = save_em_dir+'/'+trw_id+'_tile.gif'
  write_gif, gif_file, image, red, green, blue
end else begin
  wait, 5.
end

; Plot to .gif
set_plot, 'X'
window, 0, xsi=500, ysi=700

dd_load_ct, /OTHER

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

; plot the events
plot, glob_Y, glob_Z, PSYM=3, BACK=0, COLOR=1, $
	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, COLOR=10
oplot, ys, zs - 1000.0*spec_width, COLOR=10
oplot, ys, zs + 1000.0*spec_width/2.0, COLOR=10
oplot, ys, zs + 1000.0*spec_width, COLOR=10

; OK, jump in here and do the HEAE_ANAL stuff if requested...
if KEYWORD_SET(HEAE_ANAL) then begin
  ; read in the heae_whcounts.rdb structure
  heae_dir = !DDHETGCAL+'/cmp/eae/'
  heae = rdb_read(heae_dir+'heae_whcounts.rdb', /SILENT)
  ; find this trw_id in the heae
  this_one = where(heae.trw_id EQ trw_id, n_found)
  if n_found EQ 1 then begin
    this_one = this_one(0)
    ; load in mst_date, run_id, and iteration
    heae(this_one).mst_date = date_code
    heae(this_one).run_id = run_id
    heae(this_one).iteration = STRCOMPRESS(STRING(iter),/REMOVE_ALL)
    ; load in focal_time, focal_analysis=ROI,All
    heae(this_one).focal_time = hsi_time
    heae(this_one).focal_analysis = "ROI,All"
    ; load in focal_counts, focal_err_plus, focal_err_minus
    ; Find the number of events in circle of diameter heae.fp_aperture microns.
    ; about hsiY,Z = 0,0:
    ap_radius = heae(this_one).fp_aperture/2.0
    in_circ = where ( SQRT(hsi_data.Y^2 + hsi_data.Z^2) LT $
		ap_radius , n_found)
    heae(this_one).focal_counts = FLOAT(n_found)
    heae(this_one).focal_err_plus = SQRT(heae(this_one).focal_counts)
    heae(this_one).focal_err_minus = heae(this_one).focal_err_plus
    rdb_write, heae_dir+'heae_whcounts.rdb', heae, /SILENT
    ; overplot the HEAE extraction region circle
    facet_circle, 0.0, 360.0, circX, circY
    oplot, 1000.*Yloc(iter)+circX*ap_radius, $
		1000.*Zloc(iter)+circY*ap_radius, COLOR=100
    ; and label the plot
    xyouts, 1000.*Yloc(iter), 1000.*Zloc(iter) + 1.15*5000.0, $
	STRCOMPRESS(n_found)+' counts in ' + $
	STRING(hsi_time,FORMAT='(F6.1)')+' s, ' + $
	STRCOMPRESS(FIX(heae(this_one).fp_aperture)) + ' micron diameter', $
	ALIGNMENT=0.5, CHARSIZE=1.2, color=100
    ; all done, write out the heae_whcounts.rdb file
  end else begin
    print, ''
    print, ' * Did not find '+trw_id+' in heae_wcounts.rdb file ...'
    print, ''
  end
end

; 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*2.*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*2.*spec_width/2.0) )

; plot the selected events
plot, glob_Y(spec_photons), glob_Z(spec_photons), PSYM=3, BACK=0, COLOR=1, $
	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
oplot, ys, zs - 1000.0*spec_width/2.0, COLOR=10
oplot, ys, zs + 1000.0*spec_width/2.0, COLOR=10

; plot the background events
plot, glob_Y(bkg_photons), glob_Z(bkg_photons), PSYM=3, BACK=0, COLOR=1,$
	XRANGE = 1000.0*Yloc(iter)+[-9000.,9000.], $
	TITLE = 'Events in Background-extraction Regions', $
	XTITLE = 'Facility Y (um)', YTITLE = 'Facility Z (um)', $
	CHARSIZE = char_size
oplot, ys, zs - 1000.0*spec_width/2.0, COLOR=10
oplot, ys, zs - 1000.0*spec_width, COLOR=10
oplot, ys, zs + 1000.0*spec_width/2.0, COLOR=10
oplot, ys, zs + 1000.0*spec_width, COLOR=10

; 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/rowlands(grating)))
lam_min = ( ((periods(grating)/ABS(FLOAT(order)))* $
		(disp_min/rowlands(grating))) > 1.24 )
dlam = (lam_max - lam_min)/500.0

lam_spec_hist = histogram(wavelengths(spec_photons), $
		BINSIZE=dlam, MIN=lam_min, MAX=lam_max)
lam_bkg_hist = histogram(wavelengths(bkg_photons), $
		BINSIZE=dlam, MIN=lam_min, MAX=lam_max)
lams = lam_min + (lam_max-lam_min)*(INDGEN(n_elements(lam_spec_hist))+0.5)/ $
			FLOAT(n_elements(lam_spec_hist))
print, 'Lambdas: MIN, MAX: ',MIN(lams), MAX(lams)

if KEYWORD_SET(LPR) then begin
  ; Finish the .gif output
  wshow,0,iconic=0
  image = tvrd()
  tvlct, red, green, blue, /GET
  gif_file = save_em_dir+'/'+trw_id+'_image.gif'
  write_gif, gif_file, image, red, green, blue
end else begin
  wait, 5.
end

; Plot to .gif
set_plot, 'X'
window, 0, xsi=500, ysi=700
dd_load_ct, /OTHER

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

plot_io, lams, lam_spec_hist, PSYM=10, BACK=0, COLOR=1, $
	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, COLOR=100

; 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, BACK=0, COLOR=1, $
	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, COLOR=100

if KEYWORD_SET(LPR) then begin
  ; Finish the .gif output
  wshow,0,iconic=0
  image = tvrd()
  tvlct, red, green, blue, /GET
  gif_file = save_em_dir+'/'+trw_id+'_spec.gif'
  write_gif, gif_file, image, red, green, blue
end

if KEYWORD_SET(COMPARE) OR KEYWORD_SET(PS_COMPARE) then begin
  ; OK, now convert to photons/s cm^2 keV and compare with
  ; xss_sim model...
  ; HRMA area: cip_hrma - energy, shell1, shell3, shell4, shell6
  ; HSI QE: cip_hsi - energy qe
  ; grating effic: cip_heg, cip_meg, cip_leg - energy, om1, oz, op1, etc.
  ;
  ; Expected "area" (cm^2 counts/photon)
  ; Use grating energy grid
  ; Calculate it on the spec_es energy grid
    CASE grating OF
    0: begin  ; LEG
        ; y = INTERPOL(Y X x)
        hrma_area = INTERPOL(cip_hrma.ea_hrma, cip_hrma.energy, spec_es)
        grating_effic = INTERPOL(cip_leg.(((n_tags(cip_leg))/2)+order), $
				cip_leg.energy, spec_es)
      end
    1: begin  ; MEG
        hrma_area = INTERPOL((cip_hrma.ea_1 + cip_hrma.ea_3), $
				 cip_hrma.energy, spec_es)
        grating_effic = INTERPOL(cip_meg.(((n_tags(cip_meg))/2)+order), $
				cip_meg.energy, spec_es)
      end
    2: begin  ; HEG
        hrma_area = INTERPOL((cip_hrma.ea_4 + cip_hrma.ea_6), $
				 cip_hrma.energy, spec_es)
        grating_effic = INTERPOL(cip_heg.(((n_tags(cip_heg))/2)+order), $
				cip_heg.energy, spec_es)
      end
  endCASE
  hsi_qe = INTERPOL(cip_hsi.qe, cip_hsi.energy, spec_es)
  spec_sys_area = hrma_area * grating_effic * hsi_qe
  phot_spec = (spec_per_kev / (spec_sys_area * hsi_time)) > 0.0

  ; Simulate the spectrum
  xss_sim, this_meas(0) + 1, /NO_BND, /NO_PLOT
  ; This creates trw_id.spec and trw_id.spec.ps files in the current
  ; directory.
  if KEYWORD_SET(LPR) then begin
    SPAWN, 'cp '+trw_id+'.spec '+save_em_dir+'/'+trw_id+'.spec'
  end
  SPAWN, 'rm '+trw_id+'.spec*'

  if KEYWORD_SET(PS_COMPARE) then begin
    ; Plot to .ps, e.g. for SPIE paper!
    ; Parameters for desired aspect ratio
    SET_PLOT, 'PS'
    device, /portrait, font_size = 12, XSIZE=17.5, YSIZE=7.0, $
	YOFFSET=8.4, XOFFSET=2.0
    device,/COLOR
    dd_load_ct, /OTHER
    back_color = 1  ; white
    plot_color = 0  ; black
    data_color = 0
    data_style = 2
    cmdb_color = 0
    model_color = 0
    ; Fit two plots on top of each other...
    !p.multi = [0,1,2]
    char_size = 0.7
  end else begin
    ; Plot to .gif
    set_plot, 'X'
    window, 0, xsi=500, ysi=700

    dd_load_ct, /OTHER
    back_color = 0  ; black
    plot_color = 1  ; white
    data_color = 10 ; red
    data_style = 0
    cmdb_color = 100 ; yellow
    model_color = 250 ; green
    !p.multi = [0,1,2]
    char_size = 1.1
  end

  ; Limit the range for the BOTH plots based on
  ; i) having data (lam_spec_hist > 2)
  min_e_meas = MIN(spec_es(where(lam_spec_hist GT 2)))
  max_e_meas = MAX(spec_es(where(lam_spec_hist GT 2)))
  ; ii) not being too far from the X-ray line
  min_e_meas = min_e_meas > 0.6*src_energy ; 0.5*src_energy
  max_e_meas = max_e_meas < 1.7*src_energy ; 2.*src_energy

  plot_io, spec_es, phot_spec, PSYM=10, BACK=back_color, $
	COLOR=plot_color, $
	TITLE = trw_id+': '+source_str+',  '+names(grating)+ $
		'['+STRCOMPRESS(STRING(order),/REMOVE_ALL)+']', $ 
	XTITLE = 'Energy (keV) '+ $
		'[bin size = '+STRING(dlam,FORMAT='(F7.4)')+' A]', $
	YTITLE = 'photons / cm^2 s keV ', $
	YRANGE=[MIN(phot_spec(where(spec_per_kev GT 0.))), $
			MAX(phot_spec(where(spec_es LT 10.0)))], $
	XRANGE = [min_e_meas, max_e_meas], XSTYLE=1, $
	CHARSIZE = char_size, /NODATA
  oplot, spec_es, phot_spec, PSYM=10, COLOR=data_color
  ; Add indication of CMDB assumed energy
  oplot, src_energy*[1.,1.], [1.E-3,1.E9], LINESTYLE=2, COLOR=cmdb_color
  ; and model spectrum
  oplot, lx_es, lx_tubespec, COLOR=model_color  ; photons/s cm^2 keV
  ; indicate the limits
  oplot, min_e_meas*[1.,1.], [0.001,1.E6],LINESTYLE=1, COLOR=cmdb_color
  oplot, max_e_meas*[1.,1.], [0.001,1.E6],LINESTYLE=1, COLOR=cmdb_color

  ; - - - - - -
  ; Make and Plot normalized cumulative plots


  plot, spec_es, phot_spec, BACK=back_color, COLOR=plot_color, $
		TITLE = 'Normalized cummulative plots ('+ $
		date_code+'/hsi'+run_id+'i'+ $
		STRCOMPRESS(STRING(iter),/REMOVE_ALL)+'.fits)', $
		XTITLE = 'Energy (keV)', YTITLE = 'Cummulative Fraction', $
		XRANGE = [min_e_meas, max_e_meas], XSTYLE=1, $
		YRANGE=[0.,1.0], $
		CHARSIZE = char_size, /NODATA

  ; Add indication of CMDB assumed energy
  oplot, src_energy*[1.,1.], [1.E-3,1.E9], LINESTYLE=2, COLOR=cmdb_color
  ; indicate the limits
  oplot, min_e_meas*[1.,1.], [0.001,1.E9],LINESTYLE=1, COLOR=cmdb_color
  oplot, max_e_meas*[1.,1.], [0.001,1.E9],LINESTYLE=1, COLOR=cmdb_color

  meas_indices = where(spec_es GT min_e_meas AND $
			spec_es LT max_e_meas, nwhere)
  if nwhere GT 2 then begin
    cum_meas = phot_spec(meas_indices)
    cum_meas_es = spec_es(meas_indices)
    cum_meas(0) = cum_meas(0)*ABS(cum_meas_es(1)-cum_meas_es(0))
    for ib=1, n_elements(cum_meas)-1 do begin
      cum_meas(ib) = cum_meas(ib-1) + $
	cum_meas(ib)*ABS(cum_meas_es(ib)-cum_meas_es(ib-1))    
    end
    cum_meas = 1.0 - cum_meas/MAX(cum_meas)
    oplot, cum_meas_es, cum_meas, PSYM=10, COLOR=data_color, $
	LINESTYLE=data_style
  end

  model_indices = where(lx_es GT min_e_meas AND $
			lx_es LT max_e_meas, nwhere)
  if nwhere GT 2 then begin
    cum_model = lx_tubespec(model_indices)
    cum_es = lx_es(model_indices)
    cum_model(0) = cum_model(0)*ABS(cum_es(1)-cum_es(0))
    for ib=1, n_elements(cum_model)-1 do begin
      cum_model(ib) = cum_model(ib-1) + $
	cum_model(ib)*ABS(cum_es(ib)-cum_es(ib-1))    
    end
    cum_model = cum_model/MAX(cum_model)
    oplot, cum_es, cum_model, PSYM=10, COLOR=model_color
  end

  if KEYWORD_SET(PS_COMPARE) then begin
    ; finish the .ps plot
    device, /close
    set_plot, 'X'
    SPAWN, 'cp idl.ps '+save_em_dir+'/'+trw_id+'_compare.ps'
  end
  if KEYWORD_SET(LPR) then begin
    ; Finish the .gif output
    wshow,0,iconic=0
    image = tvrd()
    tvlct, red, green, blue, /GET
    gif_file = save_em_dir+'/'+trw_id+'_compare.gif'
    write_gif, gif_file, image, red, green, blue
  end

endIF else begin ; COMPARE
  phot_spec = 0.0 * spec_per_kev
end

;- - - - - - - - - - - - - -
; Output a file!
; RDB format
TAB = STRING(9B)
OPENW, out_unit, save_em_dir+'/'+trw_id+'_spec.rdb',/GET_LUN
lam_fmt = '(F9.4)'
cnt_fmt = '(I6)'
e_fmt = '(F10.6)'
flux_fmt = '(F10.3)'
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, '#exposure:'+TAB+STRCOMPRESS(STRING(hsi_time),/REMOVE_ALL)
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.4)')
printf, out_unit, '#Y_hsi:'+TAB+STRING(Yloc(iter),FORMAT='(F10.4)')
printf, out_unit, '#Z_hsi:'+TAB+STRING(Zloc(iter),FORMAT='(F10.4)')
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, '# photspec   photons / s cm^2 keV based on area, time, QE'
printf, out_unit, '#'
printf, out_unit, 'wavelen'+TAB+'counts'+TAB+'bkgcounts'+TAB+'energy'+ $
			TAB+'photspec'
printf, out_unit, 'N'+TAB+'N'+TAB+'N'+TAB+'N'+TAB+'N'
for io=0,n_elements(lams)-1 do begin
  printf, out_unit, STRING(lams(io),FORMAT=lam_fmt) + TAB  + $
	STRING(FIX(lam_spec_hist(io)),FORMAT=cnt_fmt)+ TAB + $
	STRING(FIX(lam_bkg_hist(io)),FORMAT=cnt_fmt)+ TAB + $
	STRING(spec_es(io),FORMAT=e_fmt)+ TAB + $
	STRING(phot_spec(io),FORMAT=flux_fmt)

end
close, out_unit
free_lun, out_unit

!p.multi = 0

RETURN
END