PRO p1d_anal, p1d, run_id
; Plot, etc. the Mg-K 1D slit scan data...
;

; Read in the data: do this on the command line and
; pass p1d into this routine.
; dd's ROI data:
;;p1d = rdb_read('/nfs/spectra/d6/p1d/p1d_wflux.rdb')
; Pete Ratzlaff's XSPEC fits
;;mgk = rdb_read('/nfs/spectra/d6/p1d/Mg-K.rdb')
; find out which it is to allow either to be used...
; Pete uses RUNID and I use RUN_ID
is_petes = TOTAL(STRPOS(tag_names(p1d),'RUNID') GE 0)

; If run_id is provided then use just that one, otherwise
; loop through all of them...
if n_elements(run_id) GT 0 then begin
  ; Selected run_id's
  run_id = run_id
  run_id_str = STRING(run_id,FORMAT='(I6)')
end else begin
  ; Loop through a list of all of them
  alltids = p1d.trw_id
  uids = UNIQ(alltids, SORT(alltids))
  for iu=0,n_elements(uids)-1 do begin
    if is_petes GT 0.5 then begin
      print, p1d(uids(iu)).trw_id, ' ',p1d(uids(iu)).runid
    end else begin
      print, p1d(uids(iu)).trw_id, ' ',p1d(uids(iu)).run_id
    end
  end
  RETURN
end


; Select the data
if is_petes GT 0.5 then begin
  idata = where(p1d.runid EQ run_id)
end else begin
  idata = where(p1d.run_id EQ run_id)
end


; Get the normalized data, errors, and scan location
; The data is the FP count rate - assume source flux is 
; constant and avoid dividing by noisy BNDs.
if is_petes GT 0.5 then begin
  data_norm = p1d(idata).linecounts_fpc_x2
  norm = TOTAL(1.0*data_norm)
  data_norm = data_norm/norm
  data_perror = p1d(idata).linecounts_err_fpc_x2/norm
  data_merror = p1d(idata).linecounts_err_fpc_x2/norm
  ap_str = p1d(idata(0)).ap_sz
  ; fix the "?"s
  if STRPOS(ap_str,'?') GE 0 then begin  
    ap_str = '10x200h'
    step_size = 0.0
    reads, ap_str, step_size
    ; X,Y,Z not defined either...
    data_loc = step_size*FLOAT(indgen(n_elements(idata)))
    data_loc = data_loc - TOTAL(FLOAT(data_loc))/n_elements(data_loc)
    xsource = ', assumed'
  end else begin
    step_size = 0.0
    reads, ap_str, step_size
    if step_size GT 50. then step_size = 40.0
    if STRPOS(ap_str, 'h') GT 0 then begin
      ; horizontal slit
      data_loc = p1d(idata).Z
      data_loc = data_loc - TOTAL(FLOAT(data_loc))/n_elements(data_loc)
    end else begin
      data_loc = p1d(idata).Y
      data_loc = data_loc - TOTAL(FLOAT(data_loc))/n_elements(data_loc)
    end
    xsource = ', encoder'
    ; By pass the above and use the step size...
    data_loc = step_size*FLOAT(indgen(n_elements(idata)))
    data_loc = data_loc - TOTAL(FLOAT(data_loc))/n_elements(data_loc)
    xsource = ', assumed'
  end
  ; Plot info
  if p1d(idata(0)).grating EQ 'LETG' then grat='LEG' else begin
    if p1d(idata(0)).shells EQ '13' then grat='MEG' else grat='HEG'
  end
  ord_str = string(p1d(idata(0)).order,FORMAT='(I1)')
  if ord_str NE '0' then ord_str = '+'+ord_str
  ptitle = 'Mg-K PSF/1D Slit Scan:  '+grat+' '+ord_str+',  '+ $
	run_id_str +' (XSPEC)'
  pxrange = ((n_elements(idata)-1)/2.)*step_size*[-1.,1.]
end else begin
  data_norm = p1d(idata).focal_line
  norm = TOTAL(1.0*data_norm)
  data_norm = data_norm/norm
  data_perror = p1d(idata).focal_line_plus/norm
  data_merror = p1d(idata).focal_line_minus/norm
  step_size = p1d(idata(0)).fp_aperture
  if step_size GT 50. then step_size = 40.0
  data_loc = step_size*FLOAT(indgen(n_elements(idata)))
  data_loc = data_loc - TOTAL(FLOAT(data_loc))/n_elements(data_loc)
  xsource = ', assumed'
  ; Plot info
  grat = p1d(idata(0)).grating
  ord_str = string(p1d(idata(0)).order,FORMAT='(I1)')
  if ord_str NE '0' then ord_str = '+'+ord_str
  ptitle = 'Mg-K PSF/1D Slit Scan:  '+grat+' '+ord_str+',  '+ $
	run_id_str+' (ROI)'
  pxrange = ((n_elements(idata)-1)/2.)*step_size*[-1.,1.]
end
print, ''
print, p1d(idata(0)).trw_id, ' :', $
	n_elements(data_norm), ' points,', step_size, ' um'
print, 'data_loc: '
print, data_loc
print, 'data_norm: '
print, data_norm
print, ''

; START - - - - - - hist_write_template.pro - - - - -
; Include and modify the following code to output data
; to an rdb histogram file
;
; column names
hist_struct = {bin: 0.0, count:0.0, err:0.0}
hist = REPLICATE(hist_struct, n_elements(data_loc))
;
; fill the strucuture from desired arrays...
hist.bin = data_loc/1000.0 ; mm
hist.count = data_norm*norm
hist.err = data_perror*norm
;
; create header information
; each line is an element in the string array
; each line should start with a "#"
; optional parameters have the format:
; "# param_name[:][TAB]value_string"
;
; customize for application...
your_header=['# created by p1d_anal.pro, '+SYSTIME(),$
	'# AXAF XRCF Mg-K PSF/1D scan data', $
	'# trw_id:	'+p1d(idata(0)).trw_id, $
	'# runid:	'+run_id_str, $
	'# grating:	'+grat, $
	'# order:	'+ord_str, $
	'# ap_str:	'+ap_str ]
;
; parameters that hist routines might use:
hist_header=['# for histogram routines use:', $
	'# title:	'+ptitle, $
	'# bin_axis:	Position', $
	'# bin_unit:	mm' ]
header = [your_header, hist_header]
;
; and write it out to desired file
hist_file = '/nfs/spectra/d6/p1d/p1d'+run_id_str+'_hist.rdb' 
rdb_write, hist_file, hist, HEAD = header
; END - - - - - - hist_write_template.pro - - - - -

; Fit it
fit1 = GAUSSFIT(data_loc, data_norm, fit1_params)
print, 'Data fit params:'
print, fit1_params

gaus1 = fit1_params(0)*EXP(-1.*  $
	( (data_loc-fit1_params(1))/fit1_params(2) )^2 $
	 /2.0)


; Plot

plot_io, data_loc, data_norm, PSYM=10, YRANGE=[1.E-4,1.], $
	TITLE = ptitle, $
	YTITLE = 'Count rate (relative)', $
	XTITLE = 'Position (um'+xsource+')', XRANGE = pxrange, XSTYLE=1, $
	CHARSIZE=1.1
plot_errors, data_loc, data_loc, $
	data_loc, $
	data_norm - data_merror, $
	data_norm + data_perror
oplot, data_loc, fit1, LINESTYLE=1
oplot, data_loc, gaus1, LINESTYLE=2
oplot, data_loc, fit1-gaus1, LINESTYLE=2
xyouts, 0.70*pxrange(0), 0.5, $
	'sigma core ='+STRING(fit1_params(2),FORMAT='(F6.2)')+' um'


RETURN
END