; Copyright (c) 1998-2010 A.P. Hitchcock  All rights reserved
;+
;NAME:
;	READ_XRF_SPECTRA
;
;LAST CHANGED: ----------------------------------- 24-Jul-10 (aph)
;
;PURPOSE:
;	This function reads in XGLabs X-ray fluorescence spectra from
; the  *.dta or *.hdf files written by XGLabs software
;
;CATEGORY:
;	STAND ALONE: read in utility
;
;CALLING SEQUENCE:
;	Result = READ_XRF([file=file, NOF=nof, group=group,COUNTS=COUNTS,DETECTOR=DETECTOR, ALL=ALL, $
;           HDF5=HDF5, VERBOSE=VERBOSE, MAXE=MAXE,UPPERE=UPPERE, SILENT=SILENT, PARAM=PARAM, _extra=e])
;
;CALLED FROM AXIS:
;	Read->Spectra->Twinmic
;
;INPUTS: none. All input parameters are passed as keywords.
;
;KEYWORDS:
; ALL   - read as an ALL file with multiple detector signals
; COUNTS  if set, return only 1d array with counts with zeros filled
; DETECTOR number of detector to read (0-7)
;	FILE 	filename to read
;	GROUP	group leader (Axis_ID if called from aXis2000 and axis)
;	HDF5  if set, read in from hdf file
;	MAXE	maximum energy (use as estimator of photon E max)
; NOF   no filter (default filter = '*.axb')
;	OLD		flag indicating if pre-Dec09 or post Dec-09 software / detectors
; SILENT  if set, do not display feedback on progress
; upperE  if set, only data below upperE is retained
;	VERBOSE print log information about processing
;	EXTRA 	passed on parameters
;
;OUTPUTS:
;	RESULT = AXIS 1d  structure
;
;COMMON BLOCKS: none
;
;MODIFICATION HISTORY:
; (11-Dec-09 aph) first version
; (04-mar-10 aph) allow read-in of full 8092 spectrum
; (24-jul-10 aph) adapt for new multi-column files; hdf5 files

Function read_xrf_spectra, file=file, NOF=nof, group=group,detector=detector, all=all, silent=silent, $
       hdf5=hdf5, maxE=maxE, upperE=upperE, counts=counts, verbose=verbose, old=old, _extra=e
on_error,2

IF NOT keyword_set(file) then begin
  if keyword_Set(hdf5) then fltr = '*.hdf' else $
    if keyword_set(all) then fltr='*_All_*.dta' else fltr='*.dta'
	if keyword_set(nof) then fltr='*'
	file = pickfile2(/READ, FILTER=fltr, /LPATH, DEFPATH=defpath)
ENDIF
s = ''
if file EQ '' THEN RETURN, s  ; bail-out if no filename
check = findfile(file)		 ; bail-out if non-existent file
if strlen(check(0)) EQ 0 then begin
	axis_log, 'Cannot find ' + file
	return, s
endif
t = ax_name(file)					; file format: yymmdd_hhmmss_(type).ssv)
fileshort = t(1)

; ---------------- READ FROM hdf5 FILE
 if keyword_Set(hdf5) then begin
  a = H5_browser(file, /dialog_read)
  d = a._data
  t = size(d)
  e = findgen(t(1))
  xl = 'X-ray energy (eV)'
  yl = 'XGLabs XRF spectrum (HDF5)'
  dl = fileshort   
  s = {t:'1d', x:e, d:d, dn:d, xl:xl, yl:yl, dl:dl}
 
 endif else begin

  ; -------- check if file name contains '_all_' - if so set switch
  if not keyword_set(all) then begin
  	spos = strpos(fileshort, '_all_')		; NB ax_name forces all strings to LOWER CASE !!!
  	if spos GT 1 then all = 1 else all = 0
  endif
  
  ; -------- set maximum number of channels - nb above 5000 (~ 5keV) there is noise (- rollover ?)
  max_ch = 8092
  xrfs = fltarr(3,max_ch)
  xrfs(0,*)=findgen(max_ch)
  
  old = 0     ; flag for pre-Dec09 data formats ; need VERSION numbers !!
  ; ---------------------------------------------------------------------------------------------
  
  if all NE 1 then begin
  ;  ------------ using the word 'all' in file name to differentiate 
  ; ----- single detector (or SUM) from multi-detector files 
  ; 					
  ; ---------------  reading in single detector (may be called _Sum_)
  	a = read_ascii(file, data_start=2)
  	tmp=a.(0)
  	t=size(tmp)
  	if keyword_set(verbose) then print, 'multi-column records are ', t 
  	if t(0) GT 1 then begin
      if t(1) LE 4 then begin
  ; assumes this is a SUM  or 1-detector with 4-columns [det,ch,E,cts] (post Dec-09 version of XGLab)
    		ch_min = tmp(1,0) &  ch_max = tmp(1,0)
    		for i=0, n_elements(tmp(0,*))-1 do begin
    			if tmp(1,i) LT max_ch then begin	; can have non-zero counts above max_ch
    				xrfs(1,tmp(1,i))=tmp(2,i)			; done line by line to accomodate save only non-zero channels
    				xrfs(2,tmp(1,i))=tmp(3,i)
    				if tmp(3,i) NE 0. then ch_max=tmp(1,i)
    			endif
    		endfor
      endif else begin
   ; if more than 4 columns this must be multi-detector, either 4 or 8 channel
      all = 1  
   ; ---------- allow user to select 1 or all channels
   ; ------- if not specified, ask user which detector number to read
      if t(1) GT 4  then num_det_max = t(1)/4 else num_det_max = 4
      if NOT keyword_set(detector) then begin
        num_det_prompt = 'detector (1-' + strtrim(string(fix(num_det_max)),2) +')'
        if not keyword_set(detector) then begin
          if keyword_Set(group) then detector = get_num(prompt=num_det_prompt, val=1, group=group) $
                                else detector = get_num(prompt=num_det_prompt, val=1)
        endif
   ;     if detector LT 1 OR detector GT num_det_max then detector = 2  ;KLUGE ! channel used in Dec-09
      endif
      if detector LT 1 OR detector GT num_det_max then begin
        axis_log, 'Detector number is outside allowed values of 1 - ' + strtrim(string(num_det_max),2)
        return, s 
      endif
      det_num = detector-1
    
  ; multi-column format - each detector (det,ch,E,cts) format (post Dec-09 version of XGLab)
    
      if not keyword_set(silent) then print, 'Reading in a multi-detector data file - post-Dec09 format'
        n_ch  = 1 + 4*det_num
        n_E   = 2 + 4*det_num
        n_cts = 3 + 4*det_num
        ch_min = tmp(n_ch,0) &  ch_max = min([n_elements(tmp(0,*)), max_ch])-1
        for i=0,n_elements(tmp(0,*))-1 do begin
          if tmp(n_ch,i) LT max_ch then begin ; can have non-zero counts above max_ch for save only non-zero
            xrfs(1,tmp(n_ch,i))=tmp(n_E,i)          ; energy
            xrfs(2,tmp(n_ch,i))=tmp(n_cts,i)        ; counts
            if tmp(n_cts,i) GT 0. then ch_max=tmp(n_ch,i)
          endif
        endfor
    endelse
  		
  	endif else begin					; SUM 1-column format (pre Dec-09 version of XGLab)
  		old = 1
  		npts = min([n_elements(tmp), max_ch])
  		xrfs(1,*) = findgen(max_ch)
  		xrfs(2,0:npts-1)=a.field1(0:npts-1)
  		ch_min = 0
  		ch_max = npts-1
  	endelse
  
  	if not keyword_set(silent) then begin
  		axis_log, 'Largest non-zero channel ' + string(fix(ch_max))
  		axis_log, 'at Energy = ' +string(xrfs(1,ch_max), format='(f6.1)')
  	endif
  	IF keyword_set(maxE) then return, xrfs(1,ch_max)
  
  	; generate X-ray energy scale
  	slope = (xrfs(1,ch_max) - xrfs(1,ch_min))/float((ch_max-ch_min))
  	int   =  xrfs(1,ch_min) - slope*ch_min
  	xrfs(1,*)= int + slope*findgen(max_ch)
  
  	IF keyword_set(upperE) then begin
  		tmp  = reform(xrfs(1,*))
  		q=where(tmp LE upperE, count)
  		if count NE 0 then ch_max=count
  	;	print, ' count, ch_max = ', count, ch_max
  	ENDIF
  
  endif else begin	
  ; --------- READ in for All files 
  	a = read_ascii(file, data_start=2)
  
  ; - selecting 1 [or all channels if counts option is selected)
  
  	tagn = strlowcase(tag_names(a))
  	if tagn(0) EQ 'field01' then begin
  	 	ts=size(a.field01)			;	POST Dec09 format
  	 	tmp = a.field01
  	endif else begin
  		old = 1						;	PRE Dec09 format
  		tmp = a.field1
  		ts=size(a.field1)
  	endelse
  
  ; ------- if not specified, ask user which detector number to read
  	if ts(1) GT 4  then num_det_max = ts(1)/4 else num_det_max = 4
  	num_det_prompt = 'detector (1-' + strtrim(string(fix(num_det_max)),2) +')'
  	if not keyword_set(detector) then begin
  		if keyword_Set(group) then detector = get_num(prompt=num_det_prompt, val=1, group=group) $
  		                      else detector = get_num(prompt=num_det_prompt, val=1)
  	endif
  ;	if detector LT 1 OR detector GT num_det_max then detector = 2  ;KLUGE ! channel used in Dec-09
    if detector LT 1 OR detector GT num_det_max then begin
      axis_log, 'Detector number is outside allowed values of 1 - ' + strtrim(string(num_det_max),2)
      return, s 
    endif
  	det_num = detector-1
  
  	if ts(1) GT 4 then begin 			; SUM  16-column (det,ch,E,cts) format (post Dec-09 version of XGLab)
  
  		if not keyword_set(silent) then print, 'Reading in a multi-detector data file - post-Dec09 format'
  		n_ch  = 1 + 4*det_num
  		n_E   = 2 + 4*det_num
  		n_cts = 3 + 4*det_num
  
  		ch_min = tmp(n_ch,0) &  ch_max = min([n_elements(tmp(0,*)), max_ch])-1
  		for i=0,n_elements(tmp(0,*))-1 do begin
  			if tmp(n_ch,i) LT max_ch then begin	; can have non-zero counts above max_ch for save only non-zero
  				xrfs(1,tmp(n_ch,i))=tmp(n_E,i)					; energy
  				xrfs(2,tmp(n_ch,i))=tmp(n_cts,i)				; counts
  				if tmp(n_cts,i) GT 0. then ch_max=tmp(n_ch,i)
  			endif
  		endfor
  	endif else begin					; SUM 1-column format (pre Dec-09 version of XGLab)
  
  		if not keyword_set(silent) then	print, 'Reading in a multi-detector data file - pre-Dec09 format'
  		npts = min([n_elements(tmp(0,*)), max_ch])
  		xrfs(1,*) = findgen(max_ch)
  		xrfs(2,0:npts-1)=tmp(det_num,0:npts-1)
  		ch_min = 0
  		ch_max = npts - 1
  	endelse
  
  	if not keyword_set(silent) then begin
  		if old EQ 1 then print, 'Data from old software (original detectors)' else $
   	print, 'Data from new software (post Dec-09 detectors)'
  		axis_log, 'Largest non-zero channel ' + string(fix(ch_max))
  		axis_log, 'at Energy = ' +string(xrfs(1,ch_max), format='(f6.1)')
  	endif
  	IF keyword_set(maxE) then return, xrfs(1,ch_max)
  
  	; generate X-ray energy scale
  	slope = (xrfs(1,ch_max) - xrfs(1,ch_min))/float((ch_max-ch_min))
  	int   =  xrfs(1,ch_min) - slope*ch_min
  	xrfs(1,*)= int + slope*findgen(max_ch)
  
  	IF keyword_set(upperE) then begin
  		tmp  = reform(xrfs(1,*))
  		q=where(tmp LE upperE, count)
  		if count NE 0 then ch_max=count
  	;	print, ' count, ch_max = ', count, ch_max
  	ENDIF
  
  endelse
  
  
  x = fltarr(ch_max)
  d= fltarr(ch_max)
  x = reform(xrfs(1,0:ch_max-1))
  d = reform(xrfs(2,0:ch_max-1))
  
  
  IF keyword_set(counts) then return, reform(d)
  
  ;  ------------ convert to axis 1d structure
  
  yl = 'XGLabs XRF spectrum'
  if old EQ 1 then begin
  	xl = 'channel #'
  	yl = yl + ' old dets'
  endif else  begin
  	xl = 'X-ray energy (eV)'
  	yl = yl + ' new dets'
  endelse
  if n_elements(det_num) EQ 0 then det_num ='1' else $
         det_num=strtrim(string(fix(det_num+1)),2)
  if all EQ 1 then dl = fileshort + '   Detector# ' + det_num $
    else dl = fileshort
  s = {t:'1d', x:x, d:d, dn:d, xl:xl, yl:yl, dl:dl}

ENDELSE  

  if not keyword_set(silent) then $
       axis_log, 'XRF spectrum read from ' + dl

return, s

END
