PRO hist_play, hist_file, $
	XCENTER=xcenter, REG_SIZE=reg_size, $
	PHA_RANGE = pha_range, $
	NORM=norm, $
	CUMMULATIVE=cummulative, $
        FLIP_X = flip_x, $
	CENTER_PEAK=center_peak, $ 
	OPLOT_HIST=oplot_hist, $
	GAUSSFIT=gaussfit, $
	XRANGE=xrange, YRANGE=yrange, $
	LOG=log, $
	CHAR_SCALE=char_scale
;
; KEYWORDS:
;  XCENTER, REG_SIZE - defines relevant region for
;                      normalization, fitting, plotting...
;  NORM - normalize to the total area for plotting, fitting
;  CUMMULATIVE - plots fraction-of-y-less-than-x vs x
;  FLIP_X - change the sign of the X-axis values read in
;  CENTER_PEAK - if present the bins are offset to have the peak at 0.
;  OPLOT_HIST - if set the histogram is plotted with diamonds
;          using the OPLOT command (assumes hist_play has
;          been previously run).
;          Note: use of OPLOT disables the following options:
;  GAUSSFIT - if present gaussfit is performed on the data...
;  XRANGE, YRANGE - passed on to plot routine if present
;  LOG  - if present selects LOG counts axis
;  CHAR_SCALE - if present it scales the default char size 
;
; This procedure reads in a 1D histogram file in rdb format
; and does stuff to it...
;
; The rdb file has columns: bin, count, err
; rdb parameters that may be included are: title, bin_axis, bin_unit
; (all strings when returned by rdb_param)

print, ''
if n_elements(hist_file) EQ 0 then hist_file = 'temp_hist.rdb'
print, ' file: '+hist_file

; Get the histogram data
hist = rdb_read(hist_file)
print, ' TOTAL(counts): '+STRING(TOTAL(hist.count))
;
; and check for optional parameters, setting defaults
hist_title = rdb_param(hist_file, 'title')
if STRLEN(hist_title) EQ 0 then hist_title = 'Histogram from '+hist_file
bin_axis = rdb_param(hist_file, 'bin_axis')
if STRLEN(bin_axis) EQ 0 then bin_axis = 'Bins'
bin_unit = rdb_param(hist_file, 'bin_unit')
if STRLEN(bin_unit) EQ 0 then bin_unit = ' '

count_axis = 'Counts per bin'
count_unit = ' '

if n_elements(xcenter) EQ 0 then xcenter=0.0
hist.bin = hist.bin - xcenter
; reg_size applies about zero hist.bin
if n_elements(reg_size) EQ 0 then reg_size= MAX(ABS(hist.bin))
; select only the bins in the region
sel = where(ABS(hist.bin) LE reg_size)
hist = hist(sel)

; Normalize the bin values?
if KEYWORD_SET(NORM) then begin
  norm_value = TOTAL(hist.count)
  hist.count = hist.count/norm_value
  hist.err = hist.err/norm_value
  count_unit = 'normalized'
end

; Flip the X-axis?
if KEYWORD_SET(FLIP_X) then begin
  hist.bin = -1.0*hist.bin
end

; Shift the bin values to have the peak be at zero?
if KEYWORD_SET(CENTER_PEAK) then begin
  peak_bins = where(hist.count GE 0.6*MAX(hist.count))
  peak_loc = TOTAL(hist(peak_bins).bin*hist(peak_bins).count)/ $
		TOTAL(hist(peak_bins).count)
  hist.bin = hist.bin - peak_loc
end

; Character size modified?
if n_elements(char_scale) LE 0 then char_scale = 1.0

; Fit it?
;
; Using GAUSSFIT
if KEYWORD_SET(GAUSSFIT) then begin
  ; for some reason I have to run GAUSSFIT once to get it so I do not
  ; the errors like:
  ; % All array subscripts must be same size. Var = GAUSSFIT
  ; so put this line in... arg
  tmp1 = GAUSSFIT(indgen(10), (indgen(10)-4.5)^2 )
  ; "just a guess ma'am"
  fit1_params = [MAX(hist.count), 0.0, 10.0, 0., 0., 0.]
  fit1 = GAUSSFIT(hist.bin, hist.count, fit1_params)
  print, ''
  print, 'GAUSSFIT fit params:'
  print, fit1_params
  ; evaluate just the gaussian
  gaus1 = fit1_params(0)*EXP(-1.*  $
	( (hist.bin-fit1_params(1))/fit1_params(2) )^2 $
	 /2.0)
  ; for plot...
  if bin_unit EQ 'mm' then begin
    sigma_str = 'sigma core ='+ $
	STRING(1000.0*fit1_params(2),FORMAT='(F6.1)')+' micron'
  end else begin
    if fit1_params(2) LT 100.0 then begin
      sigma_str = 'sigma core ='+ $
	STRING(fit1_params(2),FORMAT='(F6.3)')+' '+bin_unit
    end else begin
      sigma_str = 'sigma core ='+ $
	STRING(fit1_params(2),FORMAT='(F6.1)')+' '+bin_unit
    end
  end
end

; Cummulative plot?
if KEYWORD_SET(CUMMULATIVE) then begin
  ; convert hist.count to cummulative...
    nels = n_elements(hist.count)
    for ic=1,nels-1 do hist(ic).count = hist(ic).count + hist(ic-1).count
  if KEYWORD_SET(GAUSSFIT) then begin
    ; convert fit1 and gaus1 to cummulative
    nels = n_elements(fit1)
    for ic=1,nels-1 do fit1(ic) = fit1(ic) + fit1(ic-1)
    for ic=1,nels-1 do gaus1(ic) = gaus1(ic) + gaus1(ic-1)
  end
end

; Make a plot 
hist_xrange = [MIN(hist.bin), MAX(hist.bin)]
hist_xsty = 0
if KEYWORD_SET(XRANGE) then begin
  hist_xrange = xrange
  hist_xsty = 1
end
hist_yrange = [MIN(hist.count), MAX(hist.count)]
hist_ysty = 0
if KEYWORD_SET(YRANGE) then begin
  hist_yrange = yrange
  hist_ysty = 1
end

if KEYWORD_SET(LOG) AND NOT(KEYWORD_SET(OPLOT_HIST)) then begin
  plot_io, hist.bin, hist.count, PSYM=10, $
	TITLE = hist_title, $
	XTITLE = bin_axis+' ('+bin_unit+')', $
	XRANGE = hist_xrange, XSTYLE=hist_xsty, $
	YTITLE = count_axis+' ('+count_unit+')', $
	YRANGE = hist_yrange, YSTYLE=hist_ysty, $
	CHARSIZE = 1.1 * char_scale
  plot_errors, hist.bin, hist.bin, hist.bin, $
	hist.count - hist.err, $
	hist.count + hist.err
  if KEYWORD_SET(GAUSSFIT) then begin
    oplot, hist.bin, fit1, LINESTYLE=1
    oplot, hist.bin, gaus1, LINESTYLE=2
    oplot, hist.bin, fit1-gaus1, LINESTYLE=2
    xyouts, 0.15*(hist_xrange(1)-hist_xrange(0))+hist_xrange(0), $
	0.6*hist_yrange(1), $
	sigma_str
  end
end else begin
  if KEYWORD_SET(OPLOT_HIST) then begin
    oplot, hist.bin, hist.count, PSYM=-4,LINESTYLE=1
  end else begin
    plot, hist.bin, hist.count, PSYM=10, $
	TITLE = hist_title, $
	XTITLE = bin_axis+' ('+bin_unit+')', $
	XRANGE = hist_xrange, XSTYLE=hist_xsty, $
	YTITLE = count_axis+' ('+count_unit+')', $
	YRANGE = hist_yrange, YSTYLE=hist_ysty, $
	CHARSIZE = 1.1 * char_scale
    plot_errors, hist.bin, hist.bin, hist.bin, $
	hist.count - hist.err, $
	hist.count + hist.err
    if KEYWORD_SET(GAUSSFIT) then begin
      oplot, hist.bin, fit1, LINESTYLE=1
      oplot, hist.bin, gaus1, LINESTYLE=2
      oplot, hist.bin, fit1-gaus1, LINESTYLE=2
      xyouts, 0.05*(hist_xrange(1)-hist_xrange(0))+hist_xrange(0), $
	0.9*(hist_yrange(1)-hist_yrange(0))+hist_yrange(0), $
	sigma_str
    end
  end
end

print, ''

RETURN
END