; $Id: hfbsreduction__define.pro,v 1.32 2007/10/12 17:03:03 azuah Exp $
;###############################################################################
;
; NAME:
;  HFBSREDUCTION__DEFINE
;
; PURPOSE:
;  See description below.
;
; CATEGORY:
;  DAVE, HFBS, data reduction
;
; 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.
;
;###############################################################################
;+
; NAME:
;       HFBSREDUCTION__DEFINE
;
; PURPOSE:
;
;       This object widget program performs the data reduction steps for
;   data sets collected on the High Flux Backscattering Spectrometer.
;   This program is meant to run in the DAVE environment.
;
; 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
;
; CATEGORY:
;
;       Objects, widgets, data reduction, DAVE software
;
; CALLING SEQUENCE:
;
;       object = obj_new('hfbsreduction')
;
;
; INPUT PARAMETERS:
;
;       NONE
;
; INPUT KEYWORDS:
;
;   GROUP_LEADER   - Parent widget of HFBSREDUCTION object widget
;   NOTIFYID    - Vector of TOP and ID of calling widget
;   WORKDIR      - working directory, where DAVE and/or CON files will be put
;   DATDIR      - data directory, directory of raw data
;
; REQUIRED PROGRAMS:
;
;        DREBIN.PRO
;    DREBIN_POINTS_WRAPPER.PRO
;    OPAN_REBIN_WIDGET.PRO
;
; COMMON BLOCKS:
;
;       NONE
;
; RESTRICTIONS
;
;       NONE
;
; OBJECT METHODS:
;
; There are no explicitly private object methods in IDL but the
; methods used in this class are divided into PUBLIC and PRIVATE
; to indicate which ones should be used by users who wish to run
; the program from the command line, for instance.  Execution of
; the program's controls from the command line is described in the
; example below.
;
; PUBLIC OBJECT PROCEDURE METHODS:
;
;
;
; PRIVATE OBJECT PROCEDURE METHODS:
;
;
; EXAMPLE
;
;
; MODIFICATION HISTORY:
;
;       Written by Rob Dimeo, October 23, 2002.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsReductionCleanup,tlb
widget_control,tlb,get_uvalue = self
if obj_valid(self) then begin
  s = size(self->getNotifyIds())
  if s[0] ne 0 then begin
    if s[0] eq 1 then count = 0 else count = s[2]-1
    for j = 0,count do begin
      hfbsInfo = {HFBSReductionEvent,$
                            ID:(self->getNotifyIds())[0,j],$
                            Top:(self->getNotifyIds())[1,j],$
                            Handler:0l,$
                            daveFiles:self->getDaveFiles()}
      if widget_info((self->getNotifyIds())[0,j],/valid_id) then begin $
        widget_control,(self->getNotifyIds())[0,j],send_event = hfbsInfo
      endif
    endfor
  endif else begin
    ;obj_destroy,self
  endelse
endif
obj_destroy,self
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function hfbsreduction::getDaveFiles
return,self.daveFiles
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::cleanup
; First let's restore the colors
tvlct,*self.rPtr,*self.gPtr,*self.bPtr
; Now free up the pointers
ptr_free,self.odataPtr,self.oerrorPtr
ptr_free,self.oxvalsPtr,self.oyvalsPtr
ptr_free,self.dataPtr,self.errorPtr
ptr_free,self.xvalsPtr,self.yvalsPtr
ptr_free,self.notifyIdPtr;,self.daveFiles
ptr_free,self.rPtr,self.gPtr,self.bPtr
ptr_free,self.pgroup_leader,self.sigFilePtr
ptr_free,self.vanFilePtr,self.bgFilePtr
ptr_free,self.selDetPtr,self.tempPtr,self.treatmentPtr
ptr_free,self.tbmPtr,self.wbmPtr,self.dopplerPtr,self.headerPtr
ptr_free,self.titlePtr,self.detPtr
; Delete any pixmap windows
wdelete,self.winPix
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::quit,event = event
; Widget, destroy thyself
widget_control,self.tlb,/destroy
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::selSigFiles,event = event
compile_opt idl2
;files = dialog_pickfile(dialog_parent = event.top,/read, $
;                               filter = ['*.hfbs'],$
;                               title = 'Select raw signal data file(s)',$
;                               path = self.datDir, $
;                               /multiple_files)



;040407
;LRK
;MAKE THE FOLLOWING CHANGE TO ALLEVIATE THE CRASHING PROBLEM OBERVED BY Vicky, Cecile and Craig.
if !version.release ge 6.4 then begin
    files = dialog_pickfile(       /read,$
                                   filter = '*.hfbs',$
                                   title = 'Select raw signal data file(s)',$
                                   path = self.datDir, $
                                   /must_exist,$
                                   /multiple_files)
endif else begin

    files = HFBSFileList( startdir=self.datDir,$
                          filter='*.hfbs',$
                          title='Select raw signal data file(s)',$
                          commentcheck=0)
endelse

nfiles = n_elements(files)
if (1.0*total(file_test(files,/regular)) ne nfiles) then begin
  void = dialog_message(dialog_parent = event.top, $
   'At least one of the selected files is invalid')
  return
endif
*self.sigFilePtr = files
nfiles = n_elements(files)
dispFiles = strarr(nfiles)
for i = 0,nfiles-1 do begin
  text = files[i]
  textLen = strlen(text)
  pos = textLen-16
  dispFiles[i] = strmid(text,pos)
endfor
widget_control,self.fileText,set_value = dispFiles
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::selBgFiles,event = event
compile_opt idl2

;files = dialog_pickfile(dialog_parent = event.top,/read, $
;                               filter = ['*.hfbs'],$
;                               title = 'Select background data file(s)',$
;                               path = self.datDir,$
;                               /multiple_files)

;040407
;LRK
;MAKE THE FOLLOWING CHANGE TO ALLEVIATE THE CRASHING PROBLEM OBERVED BY Vicky, Cecile and Craig.
if !version.release ge 6.4 then begin
    files = dialog_pickfile(       /read,$
                                   filter = '*.hfbs',$
                                   title = 'Select background data file(s)',$
                                   path = self.datDir, $
                                   /must_exist,$
                                   /multiple_files)
endif else begin

    files = HFBSFileList( startdir=self.datDir,$
                          filter='*.hfbs',$
                          title='Select background data file(s)',$
                          commentcheck=0)
endelse


nfiles = n_elements(files)
if (1.0*total(file_test(files,/regular)) ne nfiles) then begin
  void = dialog_message(dialog_parent = event.top, $
   'At least one of the selected files is invalid')
  return
endif
*self.bgFilePtr = files
nfiles = n_elements(files)
dispFiles = strarr(nfiles)
for i = 0,nfiles-1 do begin
  text = files[i]
  textLen = strlen(text)
  pos = textLen-16
  dispFiles[i] = strmid(text,pos)
endfor
widget_control,self.bgText,set_value = dispFiles
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::selVanFile,event = event
compile_opt idl2
;files = dialog_pickfile(dialog_parent = event.top,/read, $
;files = dialog_pickfile(       /read, $
;                               filter = '*.hfbs',$
;                               title = 'Select vanadium data file',$
;                               path = self.datDir,$
;                               /must_exist,$
;                               /multiple_files)

;040407
;LRK
;MAKE THE FOLLOWING CHANGE TO ALLEVIATE THE CRASHING PROBLEM OBERVED BY Vicky, Cecile and Craig.
if !version.release ge 6.4 then begin
    files = dialog_pickfile(       /read,$
                                   filter = '*.hfbs',$
                                   title = 'Select vanadium data file',$
                                   path = self.datDir, $
                                   /must_exist,$
                                   /multiple_files)
endif else begin

    files = HFBSFileList( startdir=self.datDir,$
                          filter='*.hfbs',$
                          title='Select vanadium data file',$
                          commentcheck=0)
endelse


nfiles = n_elements(files)
if (1.0*total(file_test(files,/regular)) ne nfiles) then begin
  void = dialog_message(dialog_parent = event.top,'The selected file is invalid')
  return
endif
*self.vanFilePtr = files
nfiles = n_elements(files)
dispFiles = strarr(nfiles)
for i = 0,nfiles-1 do begin
  text = files[i]
  textLen = strlen(text)
  pos = textLen-16
  dispFiles[i] = strmid(text,pos)
endfor
widget_control,self.vanText,set_value = dispFiles
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function hfbsreduction::getnotifyIds
return,*self.notifyIdPtr
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsReductionEvents,event
;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,dialog_parent=event.top)
        catch, /cancel
        return
    endif
endif
;end error handler-------------------------------------------------------------

if dave_set_focus(event) then return
if tag_names(event,/structure_name) eq 'WIDGET_BASE' then begin
  widget_control,event.top,get_uvalue = self
  self->resize,event = event
  return
endif
widget_control,event.id,get_uvalue = cmd
call_method,cmd.method,cmd.object,event = event
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::plotData,event = event
if n_elements(*self.dataPtr) eq 0 then return
widget_control,self.grpSlider,get_value = val
val = fix(val[0])
z = reform((*self.dataPtr)[*,val])
zerr = reform((*self.errorPtr)[*,val])
x = (*self.xvalsPtr)
y = (*self.yvalsPtr)[val]
if self.autoscale eq 1 then begin
  self.xrange = [min(x),max(x)]
  dz = 0.1*(max(z)-min(z))
  self.yrange = [min(z)-max(zerr)-dz, $
                 max(z)+max(zerr)+dz]
endif
plot,x,z,psym = 4,xrange = self.xrange,yrange = self.yrange, $
     xstyle = 1,ystyle = 1,xtitle = self.xlabel, ytitle = self.zlabel, $
     title = (*self.titlePtr)[val]
errplot,x,z-zerr,z+zerr,width = 0.0

return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::doNothing,event = event
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function hfbsreduction::concatenate_detector_list
selDets = *self.selDetPtr
result = where(total(selDets,2) ne 0,newgrps)
return,selDets[result,*]
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::detgrouping,event = event
yout = hfbsDetGroupWidget(group_leader = event.top,self.selDetPtr)
*self.selDetPtr = yout
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::group_detectors,event = event
if n_elements(*self.dataPtr) eq 0 then return

data = *self.dataPtr
error = *self.errorPtr
y = *self.yvalsPtr   & ny = n_elements(y)
x = *self.xvalsPtr   & nx = n_elements(x)

detGroups = self->concatenate_detector_list()
detSize = size(detGroups)
newgrps = detSize[1]
oldGroups = detSize[2]

dout = fltarr(nx,newgrps+1)
errout = fltarr(nx,newgrps+1)
; Stuff the monitor information into the new arrays
dout[*,0] = data[*,0]
errout[*,0] = error[*,0]

; Calculate the new group values based on the grouping selected
newY = fltarr(newgrps+1)
newY[0] = 0.0  ; beam monitor
for i = 0,newgrps-1 do begin
  newY[i+1] = total(y[1:ny-1]*reform(detGroups[i,*]))/(total(detGroups[i,*]))
endfor

; Apply grouping operation on the data
for i = 0,newgrps-1 do begin
  valid = where(detGroups[i,*] eq 1,nsum)
  for j = 0,nsum-1 do begin
    dout[*,i+1] = dout[*,i+1] + data[*,valid[j]+1]
    errout[*,i+1] = sqrt((errout[*,i+1])^2 + (error[*,valid[j]+1])^2)
  endfor
endfor

*self.dataPtr = dout
*self.errorPtr = errout
*self.yvalsPtr = newY

self.ndet = newgrps
widget_control,self.grpSlider,set_slider_max = self.ndet

; Create the new title list

unit = !dave_invAngstromSym
detectorLabels = strarr(self.ndet+1)
thisFormat = '(f10.2)'
detectorLabels[0] = 'Monitor'
for i = 1,self.ndet do begin
  detectorLabels[i] = 'Q='+strtrim(string(newY[i], $
                       format = thisFormat),2)+' '+unit
endfor

*self.titlePtr = detectorLabels

treat = '16 detectors grouped into '+ $
        strtrim(string(self.ndet),2)+' groups'
*self.treatmentPtr = [*self.treatmentPtr,treat]

return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::oldebinning,event = event
if n_elements(*self.dataPtr) eq 0 then return
*self.dataPtr = *self.odataPtr
*self.xvalsPtr = *self.oxvalsPtr
*self.errorPtr = *self.oerrorPtr
treat = 'Old energy binning restored'
;*self.treatmentPtr = [*self.treatmentPtr,treat]

self->display,event = event
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::ebinning,event = event
if n_elements(*self.dataPtr) eq 0 then return
result = opan_rebin_widget(*self.xvalsPtr,group_leader = event.top)
data = *self.dataPtr
error = *self.errorPtr
datSize = size(data)
nx = datSize[1]
if result.cancel then return
if result.nbins gt nx then begin
  strout = 'Number of desired bins cannot exceed the number of input bins'
  void = dialog_message(dialog_parent = event.top,strout)
  return
endif

x = *self.xvalsPtr
oxrange = '('+strtrim(string(min(x)),2)+','+strtrim(string(max(x)),2)+')'
xrange = '('+strtrim(string(result.xlo),2)+','+strtrim(string(result.xhi),2)+')'
treat = ['Data has been rebinned in energy.', $
        'Old bins:'+strtrim(string(nx),2), $
        'Old x-range: '+oxrange, $
        'New bins:'+strtrim(string(result.nbins),2), $
        'New x-range: ',xrange]
;*self.treatmentPtr = [*self.treatmentPtr,treat]

; Ok, we've done the error checking...now do the rebinning
xlo = result.xlo & xhi = result.xhi
nbins = result.nbins
dx = (xhi-xlo)/(nbins-1.0)
x_out = xlo+dx*dindgen(nbins)

ngrps = datSize[2]
old = *self.dataPtr
errold = *self.errorPtr
err = dblarr(nbins,ngrps)
new = dblarr(nbins,ngrps)

for i = 0,ngrps-1 do begin
      z_in = old[*,i]
    dz_in = errold[*,i]
    x_in = *self.xvalsPtr
    drebin_points_wrapper,x_in,z_in,dz_in,x_out,z_out,dz_out
   new[*,i] = z_out[*]
   err[*,i] = dz_out[*]
endfor

*self.dataPtr = new
*self.errorPtr = err
*self.xvalsPtr = x_out
self->display,event = event

return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::normalize_to_vanadium,event = event
if n_elements(*self.dataPtr) eq 0 then return
if n_elements(*self.vanFilePtr) eq 0 then return
; First pull out the current data values
sigData = *self.dataPtr
sigError = *self.errorPtr
sigX = *self.xvalsPtr
sigY = *self.yvalsPtr
sigTBM = *self.tbmPtr
sigWBM = *self.wbmPtr
sigTemp = *self.tempPtr

vanFile = *self.vanFilePtr

; Wed Oct 11 15:39:54 EDT 2006
; Goran.Gasparovic@nist.gov
; Enable multiple vanadium normalizaton files.
; Compare with pro hfbsreduction::sum_runs,event = event,err = err
; Ref Begin gasparovic-061011.01
begin
	gg_nfiles=n_elements(vanFile)
	self->read_raw_hfbs,([vanFile])[0],event = event

	; Do this below only if there is more than one vanadium file.
	if n_elements(vanFile) gt 1 then begin
	; Copy first file data into temporary variables
		gg_data = *self.dataptr
		gg_temperature = *self.tempptr
		gg_doparray = *self.dopplerptr
		gg_wbm = *self.wbmptr
		gg_tbm = *self.tbmptr

		; Loop over remaining files (first one is done already).
		for gg_i=1,gg_nfiles-1 do begin
			; Read next file
			self->read_raw_hfbs,vanFile[gg_i],event = event

			; Add data, append the rest, to temporary variables
			gg_data = gg_data + (*self.dataptr)
			gg_temperature = [gg_temperature, *self.tempptr]
			gg_doparray = [gg_doparray, *self.dopplerptr]
			gg_wbm = [gg_wbm, *self.wbmptr]
			gg_tbm = [gg_tbm, *self.tbmptr]
		end

		; Good.  Now copy back from temporary variables into self pointers.
		*self.dataptr = gg_data
		*self.tempptr = gg_temperature
		*self.dopplerptr = gg_doparray
		*self.wbmptr = gg_wbm
		*self.tbmptr = gg_tbm

		; Erase gg_data for memory conservation.  The rest is small; can stay.
		gg_junk = ((temporary(gg_data))[0])

		;Excellent.  All done.
	end
end
; End gasparovic-061011.01

self->convert_to_energy,event = event
self->normalize_to_pv,event = event
self->normalize_to_monitor,event = event,nomessage = 1

xfitlo = -5.0 & xfithi = 5.0
vanX = *self.xvalsPtr
vanData = *self.dataPtr
region = where((vanX ge xfitLo) and (vanX le xfitHi))

for i = 0,self.ndet-1 do begin
  grp = i + 1
  x = vanX[region]
  y = reform(vanData[region,grp])
  ; Get a good estimate of the starting parameters for the fits
  ymax = max(y)
  sig = 1.0/2.354
  po = [ymax*sqrt(2.0*!pi*sig^2),0.0,1.0,0.0]
  ; Fit it...
  result = gaussfit(x,y,p,estimates = po,nterms = 4)
  a = p[0]*sqrt(2.0*!pi*p[2]^2)
  if i eq 0 then area = a else area = [area,a]
endfor

stats = moment(area[3:self.ndet-1])
area = temporary(area)/stats[0]
ux = 1+bytarr(n_elements(sigX))
dwf = ux#(exp(0.066*sigY^2))  ; Debye-Waller factor for vanadium
a = ux#([1.0,area])
dwf = 1.0
sigData = sigData*dwf/a
sigError = sigError*dwf/a
*self.dataPtr = sigData
*self.errorPtr = sigError
*self.xvalsPtr = sigX
*self.yvalsPtr = sigY
*self.tbmPtr = sigTBM
*self.wbmPtr = sigWBM
*self.tempPtr = sigTemp
treat = 'Normalized to vanadium file: '+strtrim(string(vanFile),2)
*self.treatmentPtr = [*self.treatmentPtr,treat]

return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function normFitFunction,parm,_Extra = extra
i = extra.index
x = *extra.xPtr
data = (*extra.dataPtr)[*,i]
mon = (*extra.dataPtr)[*,0]
b = parm[0];extra.b

; Pull data down by b
data = data - b
; Normalize to monitor
data = data/mon
; Determine the counts on the left and right sides
; of +/-5 ueV
eValue = 5.0
range_plus = where(x le -eValue)
range_minus = where(x ge eValue)
n_plus = total(data[range_plus])
n_minus = total(data[range_minus])
f = abs(n_plus-n_minus)
return,f
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::normalize_to_monitor,event = event,nomessage = nomessage, $
   barray = barray
if n_elements(*self.dataPtr) eq 0 then return
if n_elements(nomessage) eq 0 then nomessage = 0

ndet = self.ndet
dataPtr = self.dataPtr
xPtr = self.xvalsPtr
x = *xPtr
data = (*self.dataPtr)[*,1:*]
mon = (*self.dataPtr)[*,0]

if n_elements(barray) eq 0 then begin
   bArray = dblarr(self.ndet)
   wt = float(x GE 5) - float(x LE (-5))
   bArray = dblarr(ndet)
   for i = 0,ndet-1 do begin
       bArray[i] = total(wt*data[*,i]/mon)/total(wt/mon)
   ;    bArray[i] = total(wt*data[*,i]/mon^2)/total(wt/mon^2)
   endfor
endif

data = *self.dataPtr
error = *self.errorPtr
x = *self.xvalsPtr
datSize = size(data)
nchan = datSize[1]
monitor = reform(data[*,0])
monerr = reform(error[*,0])

for i = 1,self.ndet do begin
  data[*,i] = data[*,i]-bArray[i-1]
endfor

; Replace the loop with array operations for monitor normalization
uvec = 1+bytarr(self.ndet)
u = data[0:nchan-1,1:self.ndet]
v = monitor#uvec
sigu = error[0:nchan-1,1:self.ndet]
sigv = monerr#uvec

z = u/v
; Need to catch the error if either u or v are
; zeros.
uzero = where(u eq 0.0,count_uz)
vzero = where(v eq 0.0,count_vz)
if (count_uz gt 0) or (count_vz gt 0) then begin
   sigz = 1.0+fltarr(nchan,self.ndet)
endif else begin
   sigz = z*sqrt((sigu/u)^2+(sigv/v)^2)
endelse
data[0:nchan-1,1:self.ndet] = z
error[0:nchan-1,1:self.ndet] = sigz
*self.dataPtr = data
*self.errorPtr = error

*self.odataPtr = *self.dataPtr
*self.oerrorPtr = *self.errorPtr

if nomessage eq 0 then begin
  treat = 'Normalized to incident beam monitor'
  *self.treatmentPtr = [*self.treatmentPtr,treat]
endif

return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::old_normalize_to_monitor,event = event,nomessage = nomessage
if n_elements(*self.dataPtr) eq 0 then return
if n_elements(nomessage) eq 0 then nomessage = 0
data = *self.dataPtr
error = *self.errorPtr
x = *self.xvalsPtr
datSize = size(data)
nchan = datSize[1]
monitor = reform(data[*,0])
monerr = reform(error[*,0])

;; Fit a straight line to the first few points and subtract this off from
;; the entire data set
fewPts = 10
for i = 1,self.ndet do begin
  xline = fltarr(fewPts)
  yline = fltarr(fewPts)
  eline = fltarr(fewPts)
  xline[0:fewPts-1] = x[0:fewPts-1]
  yline[0:fewPts-1] = data[0:fewPts-1,i]
  eline[0:fewPts-1] = error[0:fewPts-1,i]
  result = linfit(xline,yline,measure_errors = eline)
  a = result[0] & b = result[1]
  ymin = b*xline[0]+a
  data[*,i] = data[*,i]-ymin
endfor

; Replace the loop with array operations for monitor normalization
uvec = 1+bytarr(self.ndet)
u = data[0:nchan-1,1:self.ndet]
v = monitor#uvec
sigu = error[0:nchan-1,1:self.ndet]
sigv = monerr#uvec

z = u/v
sigz = z*sqrt((sigu/u)^2+(sigv/v)^2)
data[0:nchan-1,1:self.ndet] = z
error[0:nchan-1,1:self.ndet] = sigz
*self.dataPtr = data
*self.errorPtr = error

*self.odataPtr = *self.dataPtr
*self.oerrorPtr = *self.errorPtr

if nomessage eq 0 then begin
  treat = 'Normalized to incident beam monitor'
  *self.treatmentPtr = [*self.treatmentPtr,treat]
endif

return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::read_raw_hfbs,filename,event = event,err = err
err = 0
if n_elements(filename) eq 0 then return
self.curFile = filename
nheader = 9
header = strarr(nheader)
dum = ''
openr,lun,filename,/get_lun
; Read in the header
readf,lun,header

; What is the mode of data acquisition?
modeStatement = header[2]
pos = strpos(modeStatement,':')
mode = strtrim(strmid(modeStatement,pos+1),2)
self.mode = strupcase(mode)
if strupcase(mode) ne 'VELOCITY' then begin
  self.mode = 'TIME'
endif

*self.headerPtr = header

readf,lun,dum
hash_pos = strpos(dum,'#')
if hash_pos[0] ne -1 then begin
   ; Read in the doppler frequencies
   count = 0L
   hash_pos[0] = -1
   while hash_pos[0] eq -1 do begin
      readf,lun,dum
      hash_pos = strpos(dum,'#')
      if hash_pos[0] eq -1 then begin
         if count eq 0L then dop_freq = double(dum) else $
            dop_freq = [dop_freq,double(dum)]
         count = count+1L
      endif
   endwhile
endif

dop_freq = dop_freq[0:count-1]
*self.dopplerPtr = dop_freq


hash_pos = strpos(dum,'#')
if hash_pos[0] ne -1 then begin
   ; Read in the white beam monitor counts
   count = 0L
   hash_pos[0] = -1
   while hash_pos[0] eq -1 do begin
      readf,lun,dum
      hash_pos = strpos(dum,'#')
      if hash_pos[0] eq -1 then begin
         if count eq 0L then wbm = double(dum) else $
            wbm = [wbm,double(dum)]
         count = count+1L
      endif
   endwhile
endif
wbm = wbm[0:count-1]
*self.wbmPtr = wbm

hash_pos = strpos(dum,'#')
if hash_pos[0] ne -1 then begin
   ; Read in the transmitted beam monitor counts
   count = 0L
   hash_pos[0] = -1
   while hash_pos[0] eq -1 do begin
      readf,lun,dum
      hash_pos = strpos(dum,'#')
      if hash_pos[0] eq -1 then begin
         if count eq 0L then tbm = double(dum) else $
            tbm = [tbm,double(dum)]
         count = count+1L
      endif
   endwhile
endif
tbm = tbm[0:count-1]
*self.tbmPtr = tbm

compare = dum eq '# Detector 0'
if not compare then begin
   ; There must be a temperature field here so we will read it in
   hash_pos = strpos(dum,'#')
   if hash_pos[0] ne -1 then begin
      ; Read in the transmitted beam monitor counts
      count = 0L
      hash_pos[0] = -1
      while hash_pos[0] eq -1 do begin
         readf,lun,dum
         hash_pos = strpos(dum,'#')
         if hash_pos[0] eq -1 then begin
            if count eq 0L then temp = double(dum) else $
               temp = [temp,double(dum)]
            count = count+1L
         endif
      endwhile
   endif
   temp = temp[0:count-1]
endif else begin
   ; No temperature in file so just fill up the temperature field
   ; with dummy values.
   temp = 300.+dblarr(count)
endelse
*self.tempPtr = temp

; Now read in the detectors 4096 channels at a time.
nchan = 4096
ngroup = self.ngroup
; Now read the data in big chunks, 4096 channels at a time
firstDat = {data:fltarr(2,nchan)}
readf,lun,firstDat
datStr = {instring:'',data:fltarr(2,nchan)}
str_arr = replicate(datStr,ngroup-1)
readf,lun,str_arr
data = fltarr(nchan,ngroup)
data[*,0] = reform(firstDat.data[1,*])
data[*,1:19] = reform(str_arr.data[1,0:nchan-1,*])
free_lun,lun,/force

error = sqrt(data)
wherezero = where(error eq 0,thisnum)
if thisnum gt 0 then begin
   error[wherezero] = 1.0
endif

*self.dataPtr = data
*self.errorPtr = error
*self.xvalsPtr = indgen(nchan)

*self.odataPtr = *self.dataPtr
*self.oerrorPtr = *self.errorPtr
*self.oxvalsPtr = *self.xvalsPtr
*self.oyvalsPtr = *self.yvalsPtr

self.xlabel = 'Channel'
self.ylabel = 'Group'
self.zlabel = 'Intensity'
self.ndet = 16

; Calculate the q-values for future use here
detAngles = [14.46,20.98,27.08,32.31,36.0,43.75,51.5,59.25,67.0,$
             74.75,82.5,90.25,98.0,105.75,113.5,121.5]
lamo = 6.27
q = (4.0*!pi/lamo)*sin(0.5*!dtor*detAngles)
*self.yvalsPtr = 0.0*indgen(self.ngroup)
(*self.yvalsPtr)[1:self.ndet] = q[*]

unit = !dave_invAngstromSym
detectorLabels = strarr(self.ndet)
thisFormat = '(f10.2)'
for i = 0,15 do begin
  detectorLabels[i] = 'Q='+strtrim(string(q[i],format = thisFormat),2)+' '+unit
endfor
titleArray = ['Monitor','Detector 1', 'Detector 2', 'Detector 3', $
              'Detector 4', 'Detector 5', 'Detector 6','Detector 7', $
              'Detector 8', 'Detector 9','Detector 10', 'Detector 11', $
              'Detector 12', 'Detector 13','Detector 14', 'Detector 15', $
              'Detector 16', 'Arbitrary','Monitor P(v)','Detector P(v)']
titleArray[1:16] = detectorLabels[0:15]
*self.titlePtr = titleArray




end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::old_read_raw_hfbs,filename,event = event,err = err
; This method reads a raw HFBS data file.  Note that this only works with
; files that were created after the data acquisition system converted over
; to 'CAM BINNING'.
err = 0
if n_elements(filename) eq 0 then return
self.curFile = filename
nbig = 5000
openr,lun,filename,/get_lun
; Read header information
header = strarr(10)
dummy = ''
for i = 0,9 do begin
  readf,lun,dummy
  header[i] = dummy
endfor

; What is the mode of data acquisition?
modeStatement = header[2]
pos = strpos(modeStatement,':')
mode = strtrim(strmid(modeStatement,pos+1),2)
self.mode = strupcase(mode)
if strupcase(mode) ne 'VELOCITY' then begin
  self.mode = 'TIME'
;  free_lun,lun
;  strout = 'Cannot read data file'
;  void = dialog_message(dialog_parent = event.top,strout)
;  err = 1
;  return
endif

*self.headerPtr = header
; Read in Doppler frequency
hpos = -1
dopArray = strarr(nbig)
count = 0
while hpos eq (-1) do begin
  readf,lun,dummy
  dopArray[count] = dummy
  count = count+1
  hpos = strpos(dummy,'#')
endwhile
dopArray = float(dopArray[0:count-2])
*self.dopplerPtr = dopArray

hpos = -1
wbmArray = strarr(nbig)
count = 0
while hpos eq (-1) do begin
  readf,lun,dummy
  wbmArray[count] = dummy
  count = count+1
  hpos = strpos(dummy,'#')
endwhile
wbmArray = float(wbmArray[0:count-2])
*self.wbmPtr = wbmArray

hpos = -1
tbmArray = strarr(nbig)
count = 0
while hpos eq (-1) do begin
  readf,lun,dummy
  tbmArray[count] = dummy
  count = count+1
  hpos = strpos(dummy,'#')
endwhile
tbmArray = float(tbmArray[0:count-2])
*self.tbmPtr = tbmArray

hpos = -1
tempArray = strarr(nbig)
count = 0
while hpos eq (-1) do begin
  readf,lun,dummy
  tempArray[count] = dummy
  count = count+1
  hpos = strpos(dummy,'#')
endwhile
tempArray = float(tempArray[0:count-2])
*self.tempPtr = tempArray
nchan = 4096

; Now read the data in big chunks, 4096 channels at a time
firstDat = {data:fltarr(2,nchan)}
readf,lun,firstDat
datStr = {instring:'',data:fltarr(2,nchan)}
str_arr = replicate(datStr,self.ngroup-1)
readf,lun,str_arr
data = fltarr(nchan,self.ngroup)
data[*,0] = reform(firstDat.data[1,*])
data[*,1:19] = reform(str_arr.data[1,0:nchan-1,*])
free_lun,lun

error = sqrt(data)
wherezero = where(error eq 0,thisnum)
if thisnum gt 0 then begin
   error[wherezero] = 1.0
endif

*self.dataPtr = data
*self.errorPtr = error
*self.xvalsPtr = indgen(nchan)

*self.odataPtr = *self.dataPtr
*self.oerrorPtr = *self.errorPtr
*self.oxvalsPtr = *self.xvalsPtr
*self.oyvalsPtr = *self.yvalsPtr

self.xlabel = 'Channel'
self.ylabel = 'Group'
self.zlabel = 'Intensity'
self.ndet = 16

; Calculate the q-values for future use here
detAngles = [14.46,20.98,27.08,32.31,36.0,43.75,51.5,59.25,67.0,$
             74.75,82.5,90.25,98.0,105.75,113.5,121.5]
lamo = 6.27
q = (4.0*!pi/lamo)*sin(0.5*!dtor*detAngles)
*self.yvalsPtr = 0.0*indgen(self.ngroup)
(*self.yvalsPtr)[1:self.ndet] = q[*]

unit = !dave_invAngstromSym
detectorLabels = strarr(self.ndet)
thisFormat = '(f10.2)'
for i = 0,15 do begin
  detectorLabels[i] = 'Q='+strtrim(string(q[i],format = thisFormat),2)+' '+unit
endfor
titleArray = ['Monitor','Detector 1', 'Detector 2', 'Detector 3', $
              'Detector 4', 'Detector 5', 'Detector 6','Detector 7', $
              'Detector 8', 'Detector 9','Detector 10', 'Detector 11', $
              'Detector 12', 'Detector 13','Detector 14', 'Detector 15', $
              'Detector 16', 'Arbitrary','Monitor P(v)','Detector P(v)']
titleArray[1:16] = detectorLabels[0:15]
*self.titlePtr = titleArray

return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::sum_runs,event = event,err = err
err = 0
if n_elements(*self.sigFilePtr) eq 0 then begin
  err = -1
  return
endif

files = *self.sigFilePtr
nfiles = n_elements(files)
treat = 'Files used in sum: '
for i = 0,nfiles-1 do begin
  treat = [treat,files[i]]
  self->read_raw_hfbs,files[i],event = event,err = err
  if err eq 1 then return
  if i eq 0 then begin
    data = *self.dataPtr
    temperature = *self.tempPtr
    dopArray = *self.dopplerPtr
    wbm = *self.wbmPtr
    tbm = *self.tbmPtr
  endif else begin
    data = data+*self.dataPtr
    temperature = [temperature,*self.tempPtr]
    dopArray = [dopArray,*self.dopplerPtr]
    wbm = [wbm,*self.wbmPtr]
    tbm = [tbm,*self.tbmPtr]
  endelse
endfor

*self.treatmentPtr = [*self.treatmentPtr,treat]
error = sqrt(data)
where_zero = where(data eq 0,count_zero)
if count_zero gt 0 then error[where_zero] = 1.0
*self.errorPtr = error
*self.dataPtr = data
*self.wbmPtr = wbm
*self.tbmPtr = tbm
*self.tempPtr = temperature
*self.dopplerPtr = dopArray

return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::subtract_background,event = event
;
; In order to subtract the background we must do a few things first.....
;
; 1.  Pull out the signal data and store it in local variables.  We are assuming
;    that the data has been converted to energy and normalized to P(v).
; 2.  Read in and sum up all of the background files.
; 3.  Convert to energy and normalize to P(v).
; 4.  Rebin the background onto the same energy grid as the signal file.
; 5.  Perform the subtraction scaling the background up to the signal file.
;
if n_elements(*self.dataPtr) eq 0 then return
if n_elements(*self.bgFilePtr) eq 0 then return
nbgFiles = n_elements(*self.bgFilePtr)
bgFiles = *self.bgFilePtr

; Store the signal data locally while determining the background signal
sigData = *self.dataPtr
sigError = *self.errorPtr
sigx = *self.xvalsPtr
sigy = *self.yvalsPtr
sigDoppler = *self.dopplerPtr
sigTemp = *self.tempPtr
sigWBM = *self.wbmPtr
sigTBM = *self.tbmPtr

treat = 'Background file(s) used in subtraction:'
; Sum the background files if there are more than one
for i = 0,nbgfiles-1 do begin
  treat = [treat,bgfiles[i]]
  self->read_raw_hfbs,bgfiles[i],event = event
  if i eq 0 then begin
    bgdata = *self.dataPtr
    bgtemperature = *self.tempPtr
    bgdopArray = *self.dopplerPtr
    bgwbm = *self.wbmPtr
    bgtbm = *self.tbmPtr
  endif else begin
    bgdata = bgdata+*self.dataPtr
    bgtemperature = [bgtemperature,*self.tempPtr]
    bgdopArray = [bgdopArray,*self.dopplerPtr]
    bgwbm = [bgwbm,*self.wbmPtr]
    bgtbm = [bgtbm,*self.tbmPtr]
  endelse
endfor

*self.treatmentPtr = [*self.treatmentPtr,treat]
*self.dataPtr = bgData
bgError = sqrt(bgData)
*self.errorPtr = bgError
*self.dopplerPtr = bgDopArray
*self.tempPtr = bgtemperature
*self.wbmPtr = bgwbm
*self.tbmPtr = bgtbm
; Convert the background files to energy and normalize to P(v)
self->convert_to_energy,event = event
self->normalize_to_pv,event = event

bgx = *self.xvalsPtr
bgData = *self.dataPtr
bgError = *self.errorPtr

; Now rebin the background to the same energy axis as the signal file
datSize = size(bgData)
nbins = n_elements(sigx)
ngrps = datSize[2]

drebin,bgx,bgData,bgError,sigx,data_out,ddata_out,err=err,emsg=emsg,$
   /points,/to_points

new = data_out
err = ddata_out

; Do the background subtraction
bgData = new
bgError = err
bgx = sigx

sigMon = total(sigWBM)
bgMon = total(bgwbm)
sf = (1.0*sigMon/bgMon)
newData = sigData
newError = sigError

; Add the monitors
bg_monitor = new[*,0]
sig_monitor = sigData[*,0]
newData[*,0] = bg_monitor + sig_monitor
newError[*,0] = sqrt(bgError[*,0]^2+sigError[*,0]^2)

; Cast the scale factor into an appropriate form
datSize = size(sigData)
nchan = datSize[1]
ndet = datSize[2]
sfm = (1+bytarr(nchan))#(self.sf)
newData[*,1:self.ndet] = 1.0*sigData[*,1:self.ndet] - $
                           sf*((bgData[*,1:self.ndet])*sfm)
newError[*,1:self.ndet] = sqrt(sigError[*,1:self.ndet]^2+ $
                                 (sf*(bgError[*,1:self.ndet])*sfm)^2)

;print,sigData[fix(nchan/2.0),10]
;print,bgData[fix(nchan/2.0),10]
;print,newData[fix(nchan/2.0),10]

*self.xvalsPtr = sigx
*self.dataPtr = newData
*self.errorPtr = newError
*self.dopplerPtr = [sigDoppler,bgDopArray]
*self.tempPtr = [sigTemp];[sigTemp,bgtemperature]
*self.wbmPtr = [sigWBM,bgwbm]
*self.tbmPtr = [sigTBM,bgtbm]

return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::convertTimeToEnergy,event = event

if n_elements(*self.dataPtr) eq 0 then return
datSize = size(*self.dataPtr)
nchan = datSize[1]
ndet = datSize[2]
data = *self.dataPtr
error = *self.errorPtr

res = strpos(self.curFile,'.hfbs')
compnum = double(strmid(self.curFile,res-11,8))
if compnum ge 20000601d or compnum le 19990824d then begin
  camType = 1  ; triangle cam
endif else begin
  camType = 0  ; sine cam
endelse
timebin = float(strmid((*self.headerPtr)[8],21,2))
DopFreq = (moment(*self.dopplerPtr))[0]
channels = findgen(nchan)
guesschannels = fix(1.0/(1.0e-6*DopFreq*timebin))
realchannels = guesschannels
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; begin the time to energy conversion
lamo = 6.27
eo = 81.81/lamo^2
vo = 3956.0/lamo
dms = 2.25
dsad = 4.15
dsafc = 4.24    ; distance to fission chamber
nchannels = fix(dsafc/(vo*TimeBin*1.0e-6))
if realchannels lt nchan then begin
          data[0:realchannels-1,0]=shift(data[0:realchannels-1,0],nchannels)
          error[0:realchannels-1,0]=shift(error[0:realchannels-1,0],nchannels)
endif
camData = fltarr(2,nchan)
dataNew = fltarr(nchan,16)
energy = fltarr(nchan)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Perform a circular shift of the data up to realchannels
; in order to center the data in the effective time window
ytest = fltarr(nchan)
for i=5,13 do begin
          ytest[0:nchan-1]=ytest[0:nchan-1]+data[0:nchan-1,i]
endfor
numtest = 1300
y1max = max(ytest[0:numtest-1]) & y2max = max(ytest[numtest:4095])
y1ind = where(ytest[0:numtest-1] eq y1max)
y2ind = where(ytest[numtest:4095] eq y2max)
y2ind = y2ind + numtest

y1ind = y1ind[0] & y2ind = y2ind[0]

nhalf = fix((y2ind-y1ind)/2.0)
realhalf = fix(0.5*realchannels)
halfshift = (nhalf-realhalf)

if DopFreq gt 8.0 and compnum ge 19990824d then begin
        for i = 0,16 do begin
          data[0:realchannels-1,i]=shift(data[0:realchannels-1,i],halfshift)
          error[0:realchannels-1,i]=shift(error[0:realchannels-1,i],halfshift)
        endfor
endif
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
ytest = fltarr(nchan)
for i=5,13 do begin
          ytest[0:nchan-1]=ytest[0:nchan-1]+data[0:nchan-1,i]
endfor
; Added the following line for numtest on 11/1/2000
;       numtest = 1300
numtest = 1000
y1max = max(ytest[0:numtest-1]) & y2max = max(ytest[numtest:4095])
y1ind = where(ytest[0:numtest-1] eq y1max)
y2ind = where(ytest[numtest:4095] eq y2max)
y2ind = y2ind + numtest

arbchannels = 4096
nfold = fix((y2ind-y1ind)/2)+y1ind+arbchannels/2

nfold = nfold[0]

channels = findgen(nchan+arbchannels)
energy = findgen(nchan)
newdata = fltarr(nchan+arbchannels,ndet)
jfinal = abs(nfold-fix(nchan))
longdata = fltarr(nchan+arbchannels,ndet)
longdata[arbchannels/2:nchan+arbchannels/2-1,0:ndet-1]=data[0:nchan-1,0:ndet-1]

for j=0,nfold do begin
          if nfold+j lt nchan+arbchannels then begin
            newdata[j+arbchannels/2,0:ndet-1]=longdata[nfold+j,0:ndet-1]+$
                                              longdata[nfold-j,0:ndet-1]
          endif else begin
            newdata[j+arbchannels/2,0:ndet-1]=0.0
          endelse
endfor

nfold2 = arbchannels/2+fix(realchannels/2)
numchannels = nfold2-nfold
ynew = fltarr(nchan+arbchannels,ndet)
ynew = newdata
y = fltarr(nchan+arbchannels,ndet)
y = newdata
for j=1,numchannels-1 do begin
       ynew[nfold+j,0:ndet-1]=y[nfold+j,0:ndet-1]+y[nfold2+numchannels-j,0:ndet-1]
endfor

ynew[nfold2:nchan+arbchannels-1,0:ndet-1]=0.0*findgen(nchan+arbchannels-nfold2,ndet)

indlo = nfold-y1ind-(y2ind-y1ind)/2
indhi = nfold2
indlo = indlo[0]
indhi = indhi[0]
numind = indhi-indlo+1

ydatnew = fltarr(numind,ndet)
ydatnew[0:numind-1,0:ndet-1]=ynew[indlo:indhi,0:ndet-1]

realchannels = numind
artchannels = round(1.0/(2.0e-6*TimeBin*DopFreq))
ndiff = realchannels-artchannels+1

realchannels = realchannels+ndiff
camdata = fltarr(2,2*(realchannels))

camdata = sine_cam_funQ(2*(realchannels),DopFreq)

      ;;;Determine energy from cam profile
      ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
y1new = fltarr(realchannels,ndet)
y1new[0:realchannels-1-ndiff,0:ndet-1] = ydatnew[0:realchannels-1-ndiff,0:ndet-1]

velocity = fltarr(realchannels)

if CamType eq 0 then begin               ; sine cam
         scalefactor = 1.0
         camdata = sine_cam_funQ(2*(realchannels),DopFreq)
         velocity[0:realchannels-1] = scalefactor*camData[1,0:realchannels-1]
endif else begin                         ; triangle cam
         camData = triangle_cam_funQ(2*realchannels,DopFreq)
         velocity[0:realchannels-1] = camData[1,0:realchannels-1]
endelse
energy = fltarr(realchannels)

k1 = 5.22707e-3 & k2 = 1261.64

energy[0:realchannels-1] = k1*(k2*velocity[0:realchannels-1]+$
                                  (velocity[0:realchannels-1])^2)

if realchannels mod 2 eq 0 then begin
         finalchannels = realchannels/2
         newchannels = realchannels-finalchannels+1
endif else begin
         finalchannels = (realchannels+1)/2
         newchannels = realchannels-finalchannels
endelse

newchannels = realchannels-ndiff
newenergy = fltarr(newchannels)
newenergy[0:newchannels-1] = energy[0:newchannels-1]
newdata = fltarr(newchannels,ndet)
newdata[0:newchannels-1,0:ndet-1] = y1new[0:newchannels-1,0:ndet-1]

nchan = newchannels
energy = fltarr(nchan)
energy[0:nchan-1] = newenergy[0:nchan-1]
data = fltarr(nchan,ndet)
error = fltarr(nchan,ndet)


        ;  shift the low angle detectors the appropriate amount so that
        ;  their elastic peaks line up with the high angle detectors

for i = 1,16 do begin
          ymax = max(newdata[0:nchan-1,i]) & ymax = ymax[0]
          yind = where(newdata[0:nchan-1,i] eq ymax) & yind = yind[0]
          emin = min(abs(energy[0:nchan-1])) & emin = emin[0]
          eind = where(abs(energy[0:nchan-1]) eq abs(emin)) & eind = eind[0]
          newdata[0:nchan-1,i] = shift(newdata[0:nchan-1,i],(eind-yind))
endfor
data  = newdata
error = sqrt(data)
totdetcounts = dblarr(ndet)
for i = 0,ndet-1 do begin
          totdetcounts[i] = total(double(data[0:nchan-1,i]))
endfor
;thisDetectorRates = dblarr(ndet)
;thisDetectorRates = totdetcounts/timeMinutes

*self.xvalsPtr = energy;[erange]
*self.dataPtr = data;[erange,*]
*self.errorPtr = error;[erange,*]
self.xlabel = '!6E (!4l!6eV)'

*self.oxvalsPtr = *self.xvalsPtr
*self.odataPtr = *self.dataPtr
*self.oerrorPtr = *self.errorPtr


return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::convert_to_energy,event = event
if n_elements(*self.dataPtr) eq 0 then return
datSize = size(*self.dataPtr)
nchan = datSize[1]
data = *self.dataPtr
error = *self.errorPtr

case self.mode of
'VELOCITY': begin

        energy = dblarr(nchan)
        ;de = 100.0/(nchan-1.0)
        ;energy = -50.0+de*dindgen(nchan)
        ;energy = 4.0*temporary(energy)
      de = (200.0+200.0)/(nchan - 1.0)
      energy = -200.0+de*dindgen(nchan)
        erange = where(energy le 50. and energy ge -50.,nchan)
        *self.xvalsPtr = energy[erange]
        *self.dataPtr = data[erange,*]
        *self.errorPtr = error[erange,*]
        self.xlabel = '!6E (!4l!6eV)'

        *self.oxvalsPtr = *self.xvalsPtr
        *self.odataPtr = *self.dataPtr
        *self.oerrorPtr = *self.errorPtr
      end
'TIME':   begin
            self->convertTimeToEnergy,event = event
      end
else:
endcase
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::normalize_to_pv,event = event
if n_elements(*self.dataPtr) eq 0 then return
datSize = size(*self.dataPtr)
nchan = datSize[1]
data = *self.dataPtr
error = *self.errorPtr
x = *self.xvalsPtr

fcounts = where(data[*,19] ne 0.0)
uvec = 1+bytarr(self.ndet)
pe = data[fcounts,19]/(total(data[fcounts,19]*deriv(x[fcounts])))

data[fcounts,1:16] = data[fcounts,1:16]/(pe#uvec)
error[fcounts,1:16] = error[fcounts,1:16]/(pe#uvec)

; normalize the monitor to its P(v)
fcounts = where(data[*,18] ne 0.0)
data[fcounts,0] = (data[fcounts,0]/data[fcounts,18])* $
                        total(data[fcounts,18]*deriv(x[fcounts]))
error[fcounts,0] = (error[fcounts,0]/error[fcounts,18])*$
                        total(error[fcounts,18]*deriv(x[fcounts]))
*self.dataPtr = data[fcounts,0:16]
*self.errorPtr = error[fcounts,0:16]
*self.xvalsPtr = x[fcounts]
*self.yvalsPtr = (*self.yvalsPtr)[0:16]
widget_control,self.grpSlider,set_slider_max = self.ndet

*self.oyvalsPtr = *self.yvalsPtr
*self.oxvalsPtr = *self.xvalsPtr
*self.odataPtr = *self.dataPtr
*self.oerrorPtr = *self.errorPtr

return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::display,event = event
; This method plots the data (if present) to a pixmap then to the screen.
wset,self.winPix
self->plotData,event = event
wset,self.winVis
device,copy = [0,0,!d.x_size,!d.y_size,0,0,self.winPix]
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::resize,event = event
; This method allows the user to resize the application interface while changing
; the size of the plot window accordingly.

ctrlgeom = widget_info(self.ctrlbase,/geometry)
tlbgeom = widget_info(self.tlb,/geometry)
xsize = event.x
ysize = event.y

; New data window dimensions
newxsize = xsize-ctrlgeom.xsize
newysize = ysize > ctrlgeom.ysize

widget_control,self.win,draw_xsize = newxsize, $
               draw_ysize = newysize
wdelete,self.winPix
window,/free,/pixmap,xsize = newxsize,ysize = newysize
self.winPix = !d.window
self->display,event = event
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::treatmentHandler,event = event
; This method responds to changes in the state of the
; treatment control by sensitizing or desensitizing appropriate
; buttons as necessary.
widget_control,self.treatGroup,get_value = vals
widget_control,self.vanButton,sensitive = vals[1]
widget_control,self.bgButton,sensitive = vals[2]
widget_control,self.grpButton,sensitive = vals[3]
;widget_control,self.binButton,sensitive = vals[4]
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::flashMessage_create,msg,tlb
; This method simply flashes a widget with a message that remains present
; until the flashMessage_destroy method is invoked.
;
; Center it.
geom = widget_info(self.tlb, /geometry)
xpos = geom.xoffset + geom.xsize/2 - 100
ypos = geom.yoffset + geom.ysize/2 - 50

tlb = widget_base(title='Reduction status:',/row,xoffset=xpos,yoffset=ypos, $
      tlb_frame_attr = 3)
void = widget_text(tlb,value = msg,xsize = strlen(msg),/editable)
widget_control,tlb,/realize
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::flashMessage_destroy,tlb
widget_control,tlb,/destroy
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::reduce,event = event


catch, catchError
if (catchError ne 0) then begin
    ;;print, 'Error handled!'
    eTitle = 'HFBS Reduction: Error encountered in hfbsreduction::reduce'
    eMsg = 'An error or unusual condition was encountered!'
    eMsg = [eMsg,'Please, report the following to the DAVE team:']
    eMsg = [eMsg,!error_state.msg]
    void = dialog_message(/error,eMsg,title=eTitle,dialog_parent=event.top)
    self->flashMessage_destroy,tlb
    catch, /cancel
    return
endif



if n_elements(*self.sigFilePtr) eq 0 then return
files = *self.sigFilePtr
if n_elements(*self.treatmentPtr) ne 0 then begin
  ptr_free,self.treatmentPtr
  self.treatmentPtr = ptr_new(/allocate_heap)
endif

; How do we treat the data?
widget_control,self.treatGroup,get_value = treatVal

; Do we treat these individually or sum them?
widget_control,self.runGroup,get_value = thisValue
;*self.treatmentPtr = 'Data reduction/treatment details:'
if thisValue[0] eq 0 then begin   ; treat individually

  nfiles = n_elements(files)

  for i = 0,nfiles-1 do begin
   *self.treatmentPtr = 'Data reduction/treatment details:'
   msg = 'Reading in a raw HFBS file'
   self->flashMessage_create,msg,tlb
   self->read_raw_hfbs,files[i],event = event,err = err ; always
   self->flashMessage_destroy,tlb
   if err eq 1 then return

; The following code was added August 1, 2003 (RMD) to correct a problem
; in which I originally assumed that the HFBS filename would have a standard
; format based on the date.  This allows the user to use any file with a
; .HFBS extension.
;
; First extract the filename
  filename_result = strsplit(files[i],path_sep(),/extract)
  n_comp = n_elements(filename_result)
  filename = filename_result[n_comp-1] ; filename with .HFBS attached
  extPos = strpos(filename,'.hfbs')
  filename = strmid(filename,0,extPos)


;  extPos = strpos(files[i],'.hfbs')
;  filename = strmid(files[i],extPos-11,11)
   fileOut = filepath(filename,root_dir = self.workDir)
   fileOut = fileOut + 'dyn.dave'

   treat = 'Filename: '+strtrim(files[i],2)
   *self.treatmentPtr = [*self.treatmentPtr,treat]

   msg = 'Converting to energy'
   self->flashMessage_create,msg,tlb
   self->convert_to_energy,event = event    ; always
   self->flashMessage_destroy,tlb
   msg = 'Normalizing to P(v)'
   self->flashMessage_create,msg,tlb
   self->normalize_to_pv,event = event      ; always
   self->flashMessage_destroy,tlb

    if n_elements(*self.bgFilePtr) ne 0 and $
      self.mode eq 'VELOCITY' then begin  ; subtract background
     msg = 'Subtracting background file(s)'
     self->flashMessage_create,msg,tlb
     self->subtract_background,event = event
     self->flashMessage_destroy,tlb
    endif
    if treatVal[0] eq 1 then begin
     msg = 'Normalizing to beam monitor'
     self->flashMessage_create,msg,tlb
     self->normalize_to_monitor,event = event     ; normalize to beam monitor
     self->flashMessage_destroy,tlb
    endif
    if treatVal[1] eq 1 and $
      self.mode eq 'VELOCITY' then begin
     msg = 'Normalizing to vanadium run'
     self->flashMessage_create,msg,tlb
     self->normalize_to_vanadium,event = event   ; normalize to vanadium
     self->flashMessage_destroy,tlb
    endif

    msg = 'Grouping detectors'
   self->flashMessage_create,msg,tlb
   self->group_detectors,event = event    ; group detectors
   self->flashMessage_destroy,tlb

    msg = 'Writing to a DAVE file...'
    self->flashMessage_create,msg,tlb
    self->writeDave,event = event,fileout = fileout
    self->flashMessage_destroy,tlb

   if n_elements(*self.treatmentPtr) ne 0 then begin
      ptr_free,self.treatmentPtr
      self.treatmentPtr = ptr_new(/allocate_heap)
   endif
  endfor

endif else begin         ; sum runs together
   *self.treatmentPtr = 'Data reduction/treatment details:'
   msg = 'Reading in files to sum'
   self->flashMessage_create,msg,tlb
   self->sum_runs, event = event, err = err
   self->flashMessage_destroy,tlb
   if err eq 1 then return
   msg = 'Converting to energy'
   self->flashMessage_create,msg,tlb
   self->convert_to_energy,event = event    ; always
   self->flashMessage_destroy,tlb
   msg = 'Normalizing to P(v)'
   self->flashMessage_create,msg,tlb
   self->normalize_to_pv,event = event      ; always
   self->flashMessage_destroy,tlb

; Ok, we've performed all of the essential reduction treatment.  Now
; let's do any further steps that the user wishes.

  if n_elements(*self.bgFilePtr) ne 0 and $
      self.mode eq 'VELOCITY' then begin  ; subtract background
   msg = 'Subtracting background file(s)'
   self->flashMessage_create,msg,tlb
   self->subtract_background,event = event
   self->flashMessage_destroy,tlb
  endif

  if treatVal[0] eq 1 then begin
   msg = 'Normalizing to beam monitor'
   self->flashMessage_create,msg,tlb
   self->normalize_to_monitor,event = event   ; normalize to beam monitor
   self->flashMessage_destroy,tlb
  endif

  if treatVal[1] eq 1 and $
    self.mode eq 'VELOCITY' then begin
   msg = 'Normalizing to vanadium run'
   self->flashMessage_create,msg,tlb
   self->normalize_to_vanadium,event = event    ; normalize to vanadium
   self->flashMessage_destroy,tlb
  endif

  msg = 'Grouping detectors'
  self->flashMessage_create,msg,tlb
  self->group_detectors,event = event   ; group detectors
  self->flashMessage_destroy,tlb
;
; The following code was added August 1, 2003 (RMD) to correct a problem
; in which I originally assumed that the HFBS filename would have a standard
; format based on the date.  This allows the user to use any file with a
; .HFBS extension.
;
; First extract the filename
  filename_result = strsplit(files[0],path_sep(),/extract)
  n_comp = n_elements(filename_result)
  filename = filename_result[n_comp-1] ; filename with .HFBS attached
  extPos = strpos(filename,'.hfbs')
  filename = strmid(filename,0,extPos)
  ;extPos = strpos(files[0],'.hfbs')
  ;filename = strmid(files[0],extPos-11,11)

  fileExample = filename
;  fileOut = dialog_pickfile(dialog_parent = event.top,/write, $
;                            file = filename, $
;                            filter = ['*.dave'],$
;                            title = 'DAVE output filename',$
;                            path = self.workDir)
;
  fileOut = dialog_pickfile(/write, $
                            file = filename, $
                            filter = '*.dave',$
                            title = 'DAVE output filename',$
                            path = self.workDir)


  fileOut = fileOut+'dyn.dave'


  msg = 'Writing to a DAVE file...'
  self->flashMessage_create,msg,tlb
  self->writeDave,event = event,fileout = fileout
  self->flashMessage_destroy,tlb
endelse

self->display,event = event
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::setScaleFactors,event = event
sf = hfbsdetsf(group_leader = event.top,sf = self.sf)
self.sf = sf
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::writeDave_old,event = event,fileout = fileout
; This method populates the DAVE pointer in memory and saves the data into
; a "*.DAVE" file.
treat = 'Converted to DAVE format'
*self.treatmentPtr = [*self.treatmentPtr,treat]

n = n_elements(*self.treatmentPtr)
;for i = 0,n-1 do print,(*self.treatmentPtr)[i]
davePtr = self.davePtr

data = *self.dataPtr
error = *self.errorPtr
x = *self.xvalsPtr
y = *self.yvalsPtr
;print,'Number of groups:',self.ndet
; Ok, start filling up the dave pointer...
; First, fill up the data...
*(*(*davePtr).dataStrPtr).commonStr.histPtr = $
                  {qty:data[*,1:self.ndet],$
                     err:error[*,1:self.ndet],$
                   x:x,$
                   y:y[1:self.ndet]}

datSize = size(data[*,1:self.ndet])
ndet = datSize[2]
nchan = datSize[1]

uP = strpos((*self.headerPtr)[1],':')
user = strmid((*self.headerPtr)[1],uP+1)
cP = strpos((*self.headerPtr)[0],':')
comments = strmid((*self.headerPtr)[0],cP+1)
cycP = strpos((*self.headerPtr)[4],':')
ncycles = strmid((*self.headerPtr)[4],cycP+1)
if n_elements(*self.dopplerPtr) eq 1 then begin
  fdoppler = *self.dopplerPtr
endif else begin
  fdoppler = (moment(*self.dopplerPtr))[0]
endelse
temperature = *self.tempPtr
ltP = strpos((*self.headerPtr)[7],':')
livetime = strmid((*self.headerPtr)[7],ltP+1)

nh = n_elements(*self.headerPtr)
header = (*self.headerPtr)[0:nh-2]
detSize = (size(data))[2]
crates = 60.0*float(total(data[*,1:detSize-1],1))/float(livetime)


specific = {nchan:nchan,$
            ndet:ndet,$
            header:header, $
            monPtr:ptr_new(data[*,0]),$
            monErrPtr:ptr_new(error[*,0]),$
            ncycles:fix(ncycles),$
            phiPtr:ptr_new(*self.detPtr),$
            instrument:'HFBS',$
            wavelength:6.271,$
            camType:'triangular',$
            daqmode:'inelastic',$
            dopplerFrequency:fdoppler,$
            date:systime(),$
            user:user,$
            temp_sample:temperature,$
            temp_setpoint:temperature,$
            temp_control:temperature,$
      livetime:float(livetime), $
            comments:comments,$
            sumfc:total(data[*,0]),$
            sumtbm:total(*self.tbmPtr),$
            sumwbm:total(*self.wbmPtr),$
            wbmRate:60.0*float(total(*self.tbmPtr))/float(livetime),$;
      wbmRateErr:60.0*sqrt(float(total(*self.tbmPtr)))/float(livetime),$;
      cRates:crates}

*(*(*davePtr).dataStrPtr).specificPtr = specific

hbar_omega = '!6!sh!r!e/!n !7x!6'
xlabel = hbar_omega
ylabel = 'Q '
xunits = 'energy:ueV'
yunits = 'wavevector:A-1'
histlabel = 'I (arbitrary units)'
histunits = ''

(*(*davePtr).dataStrPtr).commonStr.instrument = 'HFBS'
(*(*davePtr).dataStrPtr).commonStr.xtype = 'POINTS'
(*(*davePtr).dataStrPtr).commonStr.ytype = 'POINTS'
(*(*davePtr).dataStrPtr).commonStr.xlabel = xlabel
(*(*davePtr).dataStrPtr).commonStr.ylabel = ylabel
(*(*davePtr).dataStrPtr).commonStr.histlabel = histlabel
(*(*davePtr).dataStrPtr).commonStr.xunits = xunits
(*(*davePtr).dataStrPtr).commonStr.yunits = yunits
(*(*davePtr).dataStrPtr).commonStr.histunits = histunits
ltTreat = 'Livetime: '+strtrim(string(livetime),2)+' seconds'
wbmTreat = 'Total white beam monitor counts: '+ $
           strtrim(string(total(*self.wbmPtr)),2)
tbmTreat = 'Total transmitted beam monitor counts: '+ $
           strtrim(string(total(*self.tbmPtr)),2)
fcTreat = 'Total incident beam monitor counts: '+ $
           strtrim(string(total(data[*,0])),2)
*self.treatmentPtr = [*self.treatmentPtr,wbmTreat,tbmTreat,fcTreat]
msg = 'Number of channels: '+strtrim(string(nchan),2)
*self.treatmentPtr = [*self.treatmentPtr,msg]
msg = 'Number of groups: '+strtrim(string(ndet),2)
*self.treatmentPtr = [*self.treatmentPtr,msg]

*(*(*davePtr).dataStrPtr).commonStr.treatmentPtr = *self.treatmentPtr

if n_elements(fileout) eq 0 then return
if (*self.daveFiles)[0] eq '' then begin
  *self.daveFiles = fileOut
endif else begin
  *self.daveFiles = [*self.daveFiles,fileOut]
endelse
save, davePtr, filename = fileout

return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::writeDave,event = event,fileout = fileout
; This method populates the DAVE pointer in memory and saves the data into
; a "*.DAVE" file.
treat = 'Converted to DAVE format'
*self.treatmentPtr = [*self.treatmentPtr,treat]

n = n_elements(*self.treatmentPtr)
;for i = 0,n-1 do print,(*self.treatmentPtr)[i]
davePtr = self.davePtr

data = *self.dataPtr
error = *self.errorPtr
x = *self.xvalsPtr
y = *self.yvalsPtr
;print,'Number of groups:',self.ndet
; Ok, start filling up the dave pointer...
; First, fill up the data...


;*(*(*davePtr).dataStrPtr).commonStr.histPtr = $
;                  {qty:data[*,1:self.ndet],$
;                     err:error[*,1:self.ndet],$
;                   x:x,$
;                   y:y[1:self.ndet]}

datSize = size(data[*,1:self.ndet])
ndet = datSize[2]
nchan = datSize[1]

uP = strpos((*self.headerPtr)[1],':')
user = strmid((*self.headerPtr)[1],uP+1)
cP = strpos((*self.headerPtr)[0],':')
comments = strmid((*self.headerPtr)[0],cP+1)
cycP = strpos((*self.headerPtr)[4],':')
ncycles = strmid((*self.headerPtr)[4],cycP+1)
if n_elements(*self.dopplerPtr) eq 1 then begin
  fdoppler = *self.dopplerPtr
endif else begin
  fdoppler = (moment(*self.dopplerPtr))[0]
endelse
temperature = *self.tempPtr
ltP = strpos((*self.headerPtr)[7],':')
livetime = strmid((*self.headerPtr)[7],ltP+1)

nh = n_elements(*self.headerPtr)
header = (*self.headerPtr)[0:nh-2]
detSize = (size(data))[2]
crates = 60.0*float(total(data[*,1:detSize-1],1))/float(livetime)

specific = {nchan:nchan,$
            ndet:ndet,$
            header:header, $
            monPtr:ptr_new(data[*,0]),$
            monErrPtr:ptr_new(error[*,0]),$
            ncycles:fix(ncycles),$
            phiPtr:ptr_new(*self.detPtr),$
            instrument:'HFBS',$
            wavelength:6.271,$
            camType:'triangular',$
            daqmode:'inelastic',$
            dopplerFrequency:fdoppler,$
            date:systime(),$
            user:user,$
            temp_sample:temperature,$
            temp_setpoint:temperature,$
            temp_control:temperature,$
      livetime:float(livetime), $
            comments:comments,$
            sumfc:total(data[*,0]),$
            sumtbm:total(*self.tbmPtr),$
            sumwbm:total(*self.wbmPtr),$
            wbmRate:60.0*float(total(*self.tbmPtr))/float(livetime),$;
      wbmRateErr:60.0*sqrt(float(total(*self.tbmPtr)))/float(livetime),$;
      cRates:crates}

;*(*(*davePtr).dataStrPtr).specificPtr = specific

hbar_omega = '!6!sh!r!e/!n !7x!6'
xlabel = hbar_omega
ylabel = 'Q '
xunits = 'energy:ueV'
yunits = 'wavevector:A-1'
histlabel = 'I (arbitrary units)'
histunits = ''

;(*(*davePtr).dataStrPtr).commonStr.instrument = 'HFBS'
;(*(*davePtr).dataStrPtr).commonStr.xtype = 'POINTS'
;(*(*davePtr).dataStrPtr).commonStr.ytype = 'POINTS'
;(*(*davePtr).dataStrPtr).commonStr.xlabel = xlabel
;(*(*davePtr).dataStrPtr).commonStr.ylabel = ylabel
;(*(*davePtr).dataStrPtr).commonStr.histlabel = histlabel
;(*(*davePtr).dataStrPtr).commonStr.xunits = xunits
;(*(*davePtr).dataStrPtr).commonStr.yunits = yunits
;(*(*davePtr).dataStrPtr).commonStr.histunits = histunits
ltTreat = 'Livetime: '+strtrim(string(livetime),2)+' seconds'
wbmTreat = 'Total white beam monitor counts: '+ $
           strtrim(string(total(*self.wbmPtr)),2)
tbmTreat = 'Total transmitted beam monitor counts: '+ $
           strtrim(string(total(*self.tbmPtr)),2)
fcTreat = 'Total incident beam monitor counts: '+ $
           strtrim(string(total(data[*,0])),2)
*self.treatmentPtr = [*self.treatmentPtr,wbmTreat,tbmTreat,fcTreat]
msg = 'Number of channels: '+strtrim(string(nchan),2)
*self.treatmentPtr = [*self.treatmentPtr,msg]
msg = 'Number of groups: '+strtrim(string(ndet),2)
*self.treatmentPtr = [*self.treatmentPtr,msg]

;*(*(*davePtr).dataStrPtr).commonStr.treatmentPtr = *self.treatmentPtr
if n_elements(temperature) gt 1 then temperature = (moment(temperature))[0]
retval = create_dave_pointer(self.davePtr,   $
     qty = data[*,1:self.ndet],   $
     err = error[*,1:self.ndet],    $
     xvals = x,            $
     yvals = y[1:self.ndet],      $
     specificstr = specific,      $
     instrument = 'HFBS',      $
     xtype = 'POINTS',        $
     ytype = 'POINTS',        $
     xlabel = xlabel,         $
     ylabel = ylabel,         $
     qtlabel = histlabel,      $
     xunits = xunits,         $
     yunits = yunits,         $
     qtunits = histunits,      $
     dname = 'T',          $
     dunits = 'temperature:K',     $
     dlegend = 'Temperature',    $
     dqty = temperature,         $
     derr = 0.0,            $
     treatment = *self.treatmentptr,   $
     ermsg = ermsg           )



if n_elements(fileout) eq 0 then return
if (*self.daveFiles)[0] eq '' then begin
  *self.daveFiles = fileOut
endif else begin
  *self.daveFiles = [*self.daveFiles,fileOut]
endelse
save, davePtr, filename = fileout

return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::zoomEvents,event = event
if n_elements(*self.dataPtr) eq 0 then return
case event.type of
0: begin    ; button press
     self.mouse = event.press
     if self.mouse eq 4 then begin
       self.autoscale = 1
        self->display,event = event
     endif
     if self.mouse eq 1 then begin
       self.xbox[0] = event.x
       self.ybox[0] = event.y
        self->display,event = event
       empty
       self.autoscale = 0
       widget_control,self.win,/draw_motion_events
     endif
   end
1: begin ; button release
    if self.mouse eq 1 then begin
     xll = self.xbox[0] < self.xbox[1]
     yll = self.ybox[0] < self.ybox[1]
     w = abs(self.xbox[1] - self.xbox[0])
     h = abs(self.ybox[1] - self.ybox[0])
     xur = xll + w
     yur = yll + h
     ll = convert_coord(xll,yll,/device,/to_data)
     ur = convert_coord(xur,yur,/device,/to_data)
     self.xrange = [ll[0],ur[0]]
     self.yrange = [ll[1],ur[1]]
      self->display,event = event
     self.mouse = 0B
     widget_control,self.win,draw_motion_events = 0
    endif
    if self.mouse eq 4 then begin
     self->display,event = event
     self.mouse = 0B
     widget_control,self.win,draw_motion_events = 0
    endif
   end
2: begin ; mouse motion
     if self.mouse eq 1 then begin
        self.xbox[1] = event.x
        self.ybox[1] = event.y
        xc = [self.xbox[0],event.x,event.x,$
              self.xbox[0],$
              self.xbox[0]]
        yc = [self.ybox[0],self.ybox[0],$
              event.y,event.y,$
              self.ybox[0]]
        wset,self.winVis
        device,copy = [0,0,!d.x_size,!d.y_size,0,0,self.winPix]
        plots,xc,yc,/device
        empty
     endif
   end
else:
endcase
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::reduce_help,event = event
;pdf_file = file_which('HFBS_data_reduction.pdf',/include_current_dir)
pdf_file = !DAVE_AUXILIARY_DIR+'HFBS_data_reduction.pdf'
void = launch_help(pdf_file,tlb = event.top)
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction::createwidgets
; Widget definition method
if n_elements(*self.pgroup_leader) eq 0 then begin
  tlb = widget_base(/row, $
        title = 'HFBS Data Reduction',mbar = bar,/tlb_size_events)
endif else begin
  tlb = widget_base(group_leader = *self.pgroup_leader,/row, $
        title = 'HFBS Data Reduction',mbar = bar,/tlb_size_events)
endelse
self.tlb = tlb
self.ctrlBase = widget_base(self.tlb,/row)
lbase = widget_base(self.ctrlbase,/col)
treatTypes = ['Monitor normalization','Vanadium normalization',$
              'Background subtraction','Detector grouping'];, $
;              'Energy binning']
treatGroup = cw_bgroup(lbase,treatTypes,/col,/nonexclusive,$
           set_value = [1,0,0,0],/return_index, /frame,$
           uvalue = {object:self,method:'treatmentHandler'})
self.treatGroup = treatGroup
runTypes = ['Single sample runs','Sum all runs']
runGroup = cw_bgroup(lbase,runTypes,/col,/exclusive,$
           /no_release,set_value = 0,/return_index, /frame,$
           uvalue = {object:self,method:'doNothing'})
self.runGroup = runGroup

void = widget_button(lbase,value = 'Select signal file(s)', $
      uvalue = {object:self,method:'selSigFiles'})
self.vanButton = widget_button(lbase,value = 'Select vanadium file(s)', $
                sensitive = 0, $
               uvalue = {object:self,method:'selVanFile'})
self.bgButton = widget_button(lbase,value = 'Select background file(s)', $
                sensitive = 0, $
               uvalue = {object:self,method:'selBgFiles'})
self.grpButton = widget_button(lbase,value = 'Specify detector grouping', $
                sensitive = 0, $
               uvalue = {object:self,method:'detgrouping'})
void = widget_button(lbase,value = 'Set scale factors', $
               uvalue = {object:self,method:'setScaleFactors'})
void = widget_button(lbase,value = 'Quit', $
       uvalue = {object:self,method:'quit'})

rbase = widget_base(self.ctrlbase,/col)
void = widget_label(rbase,value = 'Signal file(s)')
xtextSize = 20
self.fileText = widget_text(rbase,ysize = 5,xsize = xtextSize,/scroll,/editable)
void = widget_label(rbase,value = 'Vanadium file(s)')
self.vanText = widget_text(rbase,ysize = 3,xsize = xtextSize,/scroll,/editable)
void = widget_label(rbase,value = 'Background file(s)')
self.bgText = widget_text(rbase,xsize = xtextSize,ysize = 5,/scroll,/editable)
void = widget_button(rbase,value = 'reduce data', $
               uvalue = {object:self,method:'reduce'})
void = widget_button(rbase,value = '?help?', $
               uvalue = {object:self,method:'reduce_help'})

plotBase = widget_base(self.tlb,/col)
winxsize = 400 & winysize = 200
self.win = widget_draw(plotBase,xsize = winxsize,ysize = winysize, $
           /button_events,uvalue = {object:self,method:'zoomEvents'})
self.grpSlider = widget_slider(plotBase,minimum = 0,maximum = 19, $
          uvalue = {object:self,method:'display'})

; Get the vertical size right
cgeom = widget_info(self.ctrlBase,/geometry)
winysize = cgeom.ysize > winysize
winxsize = winysize
widget_control,self.win,ysize = winysize,xsize = winxsize

; Center the widget
geom = widget_info(self.tlb,/geometry)
device,get_screen_size = sz
sx = sz[0] & sy = sz[1]
xoff = fix(0.5*(sx-geom.xsize))
yoff = fix(0.5*(sy-geom.ysize))
widget_control,self.tlb,xoffset = xoff,yoffset = yoff

widget_control,self.tlb,/realize

window,/free,/pixmap,xsize = winxsize,ysize = winysize
self.winPix = !d.window
widget_control,self.win,get_value = winVis
self.winVis = winVis

widget_control,tlb,set_uvalue = self
ret = dave_set_focus(self.tlb)
xmanager,'hfbsreduction::createwidgets',self.tlb, $
         event_handler = 'hfbsReductionEvents',$
         cleanup = 'hfbsReductionCleanup', /no_block
return
end;hfbsreduction::createwidgets
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function hfbsreduction::init, group_leader = group_leader, $
                 notifyId = notifyId, $
                 workDir = workDir, $
                 datDir = datDir

tvlct,r,g,b,/get
self.rPtr = ptr_new(r)
self.gPtr = ptr_new(g)
self.bPtr = ptr_new(b)
loadct,0,/silent

self.notifyIdPtr = ptr_new(/allocate_heap)
self.pgroup_leader = ptr_new(/allocate_heap)

if n_elements(notifyId) ne 0 then $
  *self.notifyIdPtr = notifyId

if n_elements(group_leader) ne 0 then begin
  *self.pgroup_leader = group_leader
  ; Get the DAVE pointer information from group leader
  ; Note that the DAVE pointer is not valid if there is no
  ; group leader specified here!
  widget_control,group_leader,get_uvalue = statePtr
  self.davePtr = (*statePtr).davePtr
  ; DO NOT FREE THE STATEPTR HERE!!!! IT BELONGS TO THE
  ; PARENT (CALLING) WIDGET!!!!
endif

cd,current = curDir
if n_elements(workDir) eq 0 then self.workDir = curDir else self.workDir = workDir
if n_elements(datDir) eq 0 then self.datDir = '' else self.datDir = datDir

self->createwidgets
self.autoscale = 1
self.xrange = [0.0,1.0]
self.yrange = [0.0,1.0]
self.xbox = self.xrange
self.ybox = self.yrange
self.mouse = 0B

; Initialize all of the pointers
self.dataPtr = ptr_new(/allocate_heap)
self.errorPtr = ptr_new(/allocate_heap)
self.xvalsPtr = ptr_new(/allocate_heap)
self.yvalsPtr = ptr_new(/allocate_heap)
self.odataPtr = ptr_new(/allocate_heap)
self.oerrorPtr = ptr_new(/allocate_heap)
self.oxvalsPtr = ptr_new(/allocate_heap)
self.oyvalsPtr = ptr_new(/allocate_heap)
self.sigFilePtr = ptr_new(/allocate_heap)
self.bgFilePtr = ptr_new(/allocate_heap)
self.vanFilePtr = ptr_new(/allocate_heap)
self.selDetPtr = ptr_new(/allocate_heap)
self.tbmPtr = ptr_new(/allocate_heap)
self.wbmPtr = ptr_new(/allocate_heap)
self.dopplerPtr = ptr_new(/allocate_heap)
self.tempPtr = ptr_new(/allocate_heap)
self.treatmentPtr = ptr_new(/allocate_heap)
self.headerPtr = ptr_new(/allocate_heap)
self.daveFiles = ptr_new(/allocate_heap)
*self.daveFiles = ''

self.ndet = 16
self.sf = 1+fltarr(self.ndet)
array = intarr(self.ndet,self.ndet)
for i = 0,self.ndet-1 do begin
  array[i,i] = 1
endfor
*self.selDetPtr = array
self.ngroup = 20

self.xlabel = ''
self.ylabel = ''
self.zlabel = ''
self.mode = 'VELOCITY'
self.curFile = ''

; Calculate the q-values for future use here
detAngles = [14.46,20.98,27.08,32.31,36.0,43.75,51.5,59.25,67.0,$
             74.75,82.5,90.25,98.0,105.75,113.5,121.5]
self.detPtr = ptr_new(/allocate_heap)
*self.detPtr = detAngles
lamo = 6.27
q = (4.0*!pi/lamo)*sin(0.5*!dtor*detAngles)
*self.yvalsPtr = 0.0*indgen(self.ngroup)
(*self.yvalsPtr)[1:self.ndet] = q[*]

self.titlePtr = ptr_new(/allocate_heap)
unit = !dave_invAngstromSym
detectorLabels = strarr(16)
thisFormat = '(f10.2)'
for i = 0,15 do begin
  detectorLabels[i] = 'Q='+strtrim(string(q[i],format = thisFormat),2)+' '+unit
endfor
titleArray = ['Monitor','Detector 1', 'Detector 2', 'Detector 3', $
              'Detector 4', 'Detector 5', 'Detector 6','Detector 7', $
              'Detector 8', 'Detector 9','Detector 10', 'Detector 11', $
              'Detector 12', 'Detector 13','Detector 14', 'Detector 15', $
              'Detector 16', 'Arbitrary','Monitor P(v)','Detector P(v)']
titleArray[1:16] = detectorLabels[0:15]

*self.titlePtr = titleArray

return,1
end;hfbsreduction::init
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro hfbsreduction__define
define = {hfbsreduction,      $

          tlb:0L,          $
          ctrlBase:0L,        $
          grpSlider:0L,       $
          runGroup:0L,        $
          treatGroup:0L,      $
          vanButton:0L,       $
          bgButton:0L,        $
          grpButton:0L,       $
          fileText:0L,        $
          vanText:0L,         $
          bgText:0L,       $

          ; Window/display variables
          win:0L,          $
          winVis:0L,       $
          winPix:0L,       $
          xrange:fltarr(2),      $
          yrange:fltarr(2),      $
          xbox:fltarr(2),      $
          ybox:fltarr(2),      $
          autoscale:1,        $
          mouse:0B,         $
          xlabel:'',       $
          ylabel:'',       $
          zlabel:'',       $
          titlePtr:ptr_new(),    $
          daveFiles:ptr_new(),     $

          ; Directory variables
          workDir:'',         $
          datDir:'',       $

          ; Raw data file variables
          tbmPtr:ptr_new(),      $
          wbmPtr:ptr_new(),      $
          dopplerPtr:ptr_new(),   $
          headerPtr:ptr_new(),     $
          detPtr:ptr_new(),      $
          tempPtr:ptr_new(),   $
          ngroup:0,         $
          ndet:0,          $
          sf:fltarr(16),      $
          mode:'',           $
          curFile:'',         $

          ; Original data set
          odataPtr:ptr_new(),    $
          oerrorPtr:ptr_new(),     $
          oxvalsPtr:ptr_new(),     $
          oyvalsPtr:ptr_new(),     $

       ; Working data set
          dataPtr:ptr_new(),   $
          errorPtr:ptr_new(),    $
          xvalsPtr:ptr_new(),    $
          yvalsPtr:ptr_new(),    $

       ; Colors for restoration upon quitting
          rPtr:ptr_new(),      $
          gPtr:ptr_new(),      $
          bPtr:ptr_new(),      $
          ; Detector pointer
          selDetPtr:ptr_new(),     $
          ; Treatment pointer
          treatmentPtr:ptr_new(),   $

          ; File pointers
          sigFilePtr:ptr_new(),   $
          bgFilePtr:ptr_new(),     $
          vanFilePtr:ptr_new(),   $

          ; Dave pointer
          davePtr:ptr_new(),   $

          ; Info about parent
          notifyIdPtr:ptr_new(), $
          pgroup_leader:ptr_new()   $

          }


return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
