; $Id: hfbs_readhscn.pro,v 1.10 2007/06/19 19:38:06 azuah Exp $
;###############################################################################
;
; NAME:
;  HFBS_READHSCN
;
; PURPOSE:
;  Reads in raw HFBS fixed-window scan files.
;
; CATEGORY:
;  DAVE, HFBS, utility
;
; AUTHOR:
;   Robert M. Dimeo, Ph.D.
;   NIST Center for Neutron Research
;   100 Bureau Drive
;   Gaithersburg, MD 20899
;   Phone: (301) 975-8135
;   E-mail: robert.dimeo@nist.gov
;   http://www.ncnr.nist.gov/staff/dimeo
;
; LICENSE:
;  The software in this file is written by an employee of
;  National Institute of Standards and Technology
;  as part of the DAVE software project.
;
;  The DAVE software package is not subject to copyright protection
;  and is in the public domain. It should be considered as an
;  experimental neutron scattering data reduction, visualization, and
;  analysis system. As such, the authors assume no responsibility
;  whatsoever for its use, and make no guarantees, expressed or
;  implied, about its quality, reliability, or any other
;  characteristic. The use of certain trade names or commercial
;  products does not imply any endorsement of a particular product,
;  nor does it imply that the named product is necessarily the best
;  product for the stated purpose. We would appreciate acknowledgment
;  if the DAVE software is used of if the code in this file is
;  included in another product.
;
;###############################################################################
pro hfbs_readhscn,file,time,ctemp,stemp,data,error,comments,$
    TBM,FC,header,xout,titles,$
    WBM = wbm,$
    norm2Mon = norm2Mon,$
    norm2One = norm2One,$
    norm2Time = norm2Time,$
    liveTime = lTime,$
    DETRATES = detrates


;begin error handler------------------------------------------------------------
; RTA - handler to catch misc. errors and provide a graceful exit so
;       that main app does not crash. Remember to switch off debug flag when debugging!
if (n_elements(!debug) && (!debug eq 0)) then begin
    catch, catchError
    if (catchError ne 0) then begin
        ;;print, 'Error handled!'
        eTitle = 'Error encountered'
        eMsg = 'An error or unusual condition was encountered!'
        eMsg = [eMsg,'Please, report the following to the DAVE team:']
        eMsg = [eMsg,!error_state.msg]
;    eMsg = [eMsg,'Function/module name: calledfunc()',!error_state.msg]
        void = dialog_message(/error,eMsg,title=eTitle,/center);,dialog_parent=event.top)
        catch, /cancel
        return
    endif
endif
;end error handler-------------------------------------------------------------

; Modified this routine to read in files in the old data format (created prior
; to 19990901): RMD,11/2/2000
;

if n_elements(norm2Mon) eq 0 then norm2Mon = 0
if n_elements(norm2One) eq 0 then norm2One = 0
if n_elements(norm2Time) eq 0 then norm2Time = 0

res = strpos(file,'.hscn')
;print,'Filename: ',file
;print,'File part: ',strmid(file,res-11,1)
if strmid(file,res-11,1) ne '1' and strmid(file,res-11,1) ne '2' then begin
;  print,'Not a usual scan file'
  present = 1
endif else begin

  compnum = double(strmid(file,res-11,8))
;  print,'Compnum: ',compnum
  if compnum ge 19990901d then begin
    if compnum ge 20051012d then begin
        ;101305
        ;ADDING THIS OPTION FOLLOWING THE ADDITION OF A NEW TEMP PROBE.
        ;HANDLE THIS IN if STATEMENTS BELOW.
        present = 2
    endif else begin
        present = 1
    endelse

  endif else begin
    if compnum lt 19990901d and compnum gt 19990601d then begin
      present = -1
    endif
    if compnum lt 19990601d then begin
      present = 0
    endif
  endelse
endelse
;print,'Present: ',present

OPENR, lun, file, /GET_LUN

header=strarr(5)
comments = ''
dumstr = ''
dumstr1 = ''
dumstrmode = ''
ch = 0
intensity = 0.0

;Begin reading in information

READF, lun, dumstr     ; read comments
header[0] = dumstr
comments = strmid(dumstr,14,50)
;state.comments[buffer]=strmid(dumstr,21)
if present eq 1 or present eq 2 then begin
  READF, lun, dumstrmode   ; read users
  header[1] = dumstrmode
endif else begin
  if present eq -1 then begin
    header[1] = ''
  endif
  if present eq 0 then begin
    header[1] = ''
  endif
endelse

READF,lun, dumstrmode   ; read mode
header[2] = dumstrmode
mode = strmid(dumstrmode,14,4)
READF,lun, dumstrmode
header[3] = dumstrmode
;state.start_time[buffer]=strmid(dumstr,21)
READF, lun, dumstr     ; read in the time/pt
header[4] = dumstr
READF, lun, dumstr
titles = dumstr

;present = strpos(dumstr,'CTEMP')
;present = present[0]
;print,present
;present = 0

if present eq -1 then begin     ; old file

  ndet1 = 23
  nbig  = 10000
  x     = dblarr(ndet1,nbig)
  x1    = dblarr(ndet1)
  count = 0
  data = fltarr(ndet1)
  WHILE (NOT EOF(lun)) DO BEGIN
     READF, lun, data,format='(23F0)'
     x[0:ndet1-1,count]=data[0:ndet1-1]
     count = count + 1
     ;print,x[0,count-1]
  ENDWHILE
  FREE_LUN, lun
  nchan = count

  nchan = count
  xout = fltarr(ndet1,nchan)
  xout[0:ndet1-1,0:nchan-1] = x[0:ndet1-1,0:nchan-1]

  ctemp = fltarr(nchan) & ctemp[0:nchan-1] = x[5,0:nchan-1]
  stemp = fltarr(nchan) & stemp[0:nchan-1] = x[5,0:nchan-1]
  time = fltarr(nchan) & time[0:nchan-1] = x[0,0:nchan-1]


  ;scan_livetimePtr = ptr_new(x[1,0:nchan-1])
  wbm = (x[2,0:nchan-1])
  tbm = (x[3,0:nchan-1])
  FC = (x[22,0:nchan-1])
  sumFC = total(FC)


  ndet = 16
  data = dblarr(ndet1,nchan)
  for i = 0,nchan-1 do begin
    data[1:ndet,i]=x[6:ndet1-2,i];/x[1,i]
  endfor
  data[0,0:nchan-1] = x[ndet1-1,0:nchan-1]
  data[17,0:nchan-1] = x[2,0:nchan-1]
  data[18,0:nchan-1] = x[3,0:nchan-1]
  error = dblarr(ndet1,nchan)
  error[0:ndet1-1,0:nchan-1] = sqrt(data[0:ndet1-1,0:nchan-1])
  for i=0,nchan-1 do begin
    data[0:ndet-1,i]=data[0:ndet-1,i]
    error[0:ndet-1,i]=error[0:ndet-1,i]
  endfor
  energy = FLTARR(nchan)

  if mode eq 'time' then energy[0:nchan-1] = x[5,0:nchan-1]
  if mode eq 'temp' then energy[0:nchan-1] = x[5,0:nchan-1]
  xlabel = 'temp'
  BuffLabel = 'temp'
  moncount = 0l
  moncount = TOTAL(x[22,0:nchan-1])
endif else begin  ; new files which record both temperatures
  if present eq 2 then begin ;NEW FILES WITH 3 TEMPERATURES
      ndet1 = 25
      nbig  = 10000
      x     = dblarr(ndet1,nbig)
      x1    = dblarr(ndet1)
      count = 0
      data = fltarr(ndet1)
      WHILE (NOT EOF(lun)) DO BEGIN
         READF, lun, data,format='(25F0)'
         x[0:ndet1-1,count]=data[0:ndet1-1]
         count = count + 1
         ;print,x[0,count-1]
      ENDWHILE
      FREE_LUN, lun
      nchan = count

      nchan = count
      xout = fltarr(ndet1,nchan)
      xout[0:ndet1-1,0:nchan-1] = x[0:ndet1-1,0:nchan-1]

      ;THE NEW TEMPERATURE DATA IS IN COLUMN 7
      ;SHIFT ALL LATER COLUMNS ONE COLUMN TO THE RIGHT.
      auxtemp = fltarr(nchan) & auxtemp[0:nchan-1] = x[7,0:nchan-1]

      ctemp = fltarr(nchan) & ctemp[0:nchan-1] = x[5,0:nchan-1]
      stemp = fltarr(nchan) & stemp[0:nchan-1] = x[6,0:nchan-1]
      time = fltarr(nchan) & time[0:nchan-1] = x[0,0:nchan-1]


      ;scan_livetimePtr = ptr_new(x[1,0:nchan-1])
      ltime = reform(x[1,0:nchan-1])
      wbm = (x[2,0:nchan-1])
      tbm = (x[3,0:nchan-1])


      ;SHIFT fc COLUMN
      fc = (x[24,0:nchan-1])
      sumFC = total(fc)

      ndet = 16
      data = dblarr(ndet1,nchan)
      for i = 0,nchan-1 do begin
        ;SHIFT DATA COLUMN
        data[1:ndet,i]=x[8:ndet1-2,i];/x[1,i]
      endfor
      data[0,0:nchan-1] = x[ndet1-1,0:nchan-1]

      ;DON'T NEED TO CHANGE THIS
      data[17,0:nchan-1] = x[2,0:nchan-1]
      data[18,0:nchan-1] = x[3,0:nchan-1]
      error = dblarr(ndet1,nchan)
      error[0:ndet1-1,0:nchan-1] = sqrt(data[0:ndet1-1,0:nchan-1])
      for i=0,nchan-1 do begin
        data[0:ndet-1,i]=data[0:ndet-1,i]
        error[0:ndet-1,i]=error[0:ndet-1,i]
      endfor
    ;;;;;;;;;;;;
    ;  normalize to fission chamber
      ;for i = 0,nchan-1 do begin
      ;  data[1:ndet,i] = data[1:ndet,i]/x[23,i]
      ;  error[1:ndet,i] = error[1:ndet,i]/x[23,i]
      ;endfor
    ;;;;;;;;;;;;
      energy = FLTARR(nchan)

      if mode eq 'time' then energy[0:nchan-1] = x[5,0:nchan-1]
      if mode eq 'temp' then energy[0:nchan-1] = x[5,0:nchan-1]
      xlabel = 'temp'
      BuffLabel = 'temp'
      moncount = 0L

      ;SHIFT THIS BY 1 COLUMN
      moncount = TOTAL(x[24,0:nchan-1])
  endif;present eq 2

  if present eq 1 then begin
  ndet1 = 24
  nbig  = 10000
  x     = dblarr(ndet1,nbig)
  x1    = dblarr(ndet1)
  count = 0
  data = fltarr(ndet1)
  WHILE (NOT EOF(lun)) DO BEGIN
     READF, lun, data,format='(24F0)'
     x[0:ndet1-1,count]=data[0:ndet1-1]
     count = count + 1
     ;print,x[0,count-1]
  ENDWHILE
  FREE_LUN, lun
  nchan = count

  nchan = count
  xout = fltarr(ndet1,nchan)
  xout[0:ndet1-1,0:nchan-1] = x[0:ndet1-1,0:nchan-1]

  ctemp = fltarr(nchan) & ctemp[0:nchan-1] = x[5,0:nchan-1]
  stemp = fltarr(nchan) & stemp[0:nchan-1] = x[6,0:nchan-1]
  time = fltarr(nchan) & time[0:nchan-1] = x[0,0:nchan-1]


  ;scan_livetimePtr = ptr_new(x[1,0:nchan-1])
  ltime = reform(x[1,0:nchan-1])
  wbm = (x[2,0:nchan-1])
  tbm = (x[3,0:nchan-1])
  fc = (x[23,0:nchan-1])
  sumFC = total(fc)

  ndet = 16
  data = dblarr(ndet1,nchan)
  for i = 0,nchan-1 do begin
    data[1:ndet,i]=x[7:ndet1-2,i];/x[1,i]
  endfor
  data[0,0:nchan-1] = x[ndet1-1,0:nchan-1]
  data[17,0:nchan-1] = x[2,0:nchan-1]
  data[18,0:nchan-1] = x[3,0:nchan-1]
  error = dblarr(ndet1,nchan)
  error[0:ndet1-1,0:nchan-1] = sqrt(data[0:ndet1-1,0:nchan-1])
  for i=0,nchan-1 do begin
    data[0:ndet-1,i]=data[0:ndet-1,i]
    error[0:ndet-1,i]=error[0:ndet-1,i]
  endfor
;;;;;;;;;;;;
;  normalize to fission chamber
  ;for i = 0,nchan-1 do begin
  ;  data[1:ndet,i] = data[1:ndet,i]/x[23,i]
  ;  error[1:ndet,i] = error[1:ndet,i]/x[23,i]
  ;endfor
;;;;;;;;;;;;
  energy = FLTARR(nchan)

  if mode eq 'time' then energy[0:nchan-1] = x[5,0:nchan-1]
  if mode eq 'temp' then energy[0:nchan-1] = x[5,0:nchan-1]
  xlabel = 'temp'
  BuffLabel = 'temp'
  moncount = 0l
  moncount = TOTAL(x[23,0:nchan-1])
  endif;present eq 1
  if present eq 0 then begin  ; can read in Ivkov's files (really old!!!)

      ndet1 = 23
      nbig  = 2500
      x     = dblarr(ndet1,nbig)
      x1    = dblarr(ndet1)
      count = 0
      data = fltarr(ndet1)
      WHILE (NOT EOF(lun)) DO BEGIN
         READF, lun, data,format='(23F0)'
         x[0:ndet1-1,count]=data[0:ndet1-1]
         count = count + 1
         ;print,x[0,count-1]
      ENDWHILE
      FREE_LUN, lun
      nchan = count
      xout = fltarr(ndet1,nchan)
      xout[0:ndet1-1,0:nchan-1] = x[0:ndet1-1,0:nchan-1]

      ctemp = fltarr(nchan) & ctemp[0:nchan-1] = x[0,0:nchan-1]
      stemp = fltarr(nchan) & stemp[0:nchan-1] = x[0,0:nchan-1]
      time = fltarr(nchan) & time[0:nchan-1] = x[0,0:nchan-1]

      ;scan_livetimePtr = ptr_new(x[0,0:nchan-1])
      wbm = (x[1,0:nchan-1])
      TBM = (x[1,0:nchan-1])
      FC = (x[1,0:nchan-1])
      sumFC = total(FC)

      ndet = 20
      data = dblarr(ndet1,nchan)
      for i = 0,nchan-1 do begin
        data[1:16,i]=x[2:17,i];/x[1,i]
      endfor
      data[0,0:nchan-1] = x[1,0:nchan-1]
      data[17,0:nchan-1] = x[2,0:nchan-1]
      data[18,0:nchan-1] = x[3,0:nchan-1]
      error = dblarr(ndet1,nchan)
      error[0:ndet1-1,0:nchan-1] = sqrt(data[0:ndet1-1,0:nchan-1])
      for i=0,nchan-1 do begin
        data[0:ndet-1,i]=data[0:ndet-1,i]
        error[0:ndet-1,i]=error[0:ndet-1,i]
      endfor
      energy = FLTARR(nchan)

      if mode eq 'time' then energy[0:nchan-1] = x[0,0:nchan-1]
      if mode eq 'temp' then energy[0:nchan-1] = x[0,0:nchan-1]
      xlabel = 'temp'
      BuffLabel = 'temp'
      moncount = 0l
      moncount = TOTAL(x[1,0:nchan-1])
  endif
endelse

liveTime = (nchan gt 1)? time[n_elements(time)-1] : ltime

; calculate the count rates in each detector (cps)
detRates = (nchan gt 1)? $
           (total(data[1:16,*],2))/float(liveTime) : $
           data[1:16,*]/float(liveTime[0])

if n_elements(norm2Mon) eq 0 then norm2Mon = 1
if norm2Time eq 1 then begin
  ndet = 17
  u = 1+bytarr(ndet-1)
  data[1:ndet-1,*] = data[1:ndet-1,0:nchan-1]/(u#ltime)
  error[1:ndet-1,*] = error[1:ndet-1,0:nchan-1]/(u#ltime)
  ;for i = 1,16 do begin
  ;  data[i,0:nchan-1] = data[i,0:nchan-1]/ltime[0:nchan-1]
  ;  error[i,0:nchan-1] = error[i,0:nchan-1]/ltime[0:nchan-1]
  ;endfor
endif
if norm2Mon eq 1 then begin
  ndet = 17
  u = 1+bytarr(ndet-1)
  data[1:ndet-1,*] = data[1:ndet-1,*]/(u#reform(data[0,*]))
  bad_data = where(~finite(data),count_bad)
  if count_bad gt 0 then data[bad_data] = 0.0
  error[1:ndet-1,*] = error[1:ndet-1,*]/(u#reform(data[0,*]))
  if count_bad gt 0 then error[bad_data] = 0.0
;  for i = 1,16 do begin
;    data[i,0:nchan-1] = data[i,0:nchan-1]/data[0,0:nchan-1]
;    error[i,0:nchan-1] = error[i,0:nchan-1]/data[0,0:nchan-1]
;  endfor
endif
if norm2One eq 1 then begin
  for i = 1,16 do begin
    dataMax = max(data[i,0:nchan-1])
    data[i,0:nchan-1] = data[i,0:nchan-1]/dataMax
    error[i,0:nchan-1] = error[i,0:nchan-1]/dataMax
  endfor
endif

return
end
