; $Id: fcsdata__define.pro,v 1.4 2005/02/08 17:20:27 azuah Exp $
;+
; NAME:
;       FCSDATA__DEFINE
;
; PURPOSE:
;
;       Object class for data collected on the Fermi Chopper
;		Spectrometer (FCS) at the NIST Center for Neutron Research.
;
; AUTHORS:
;       Craig Brown
;		NIST Center for Neutron Research
;       100 Bureau Drive
;		Gaithersburg, MD 20899
;       Phone: (301) 975-5134
;       E-mail: craig.brown@nist.gov
;                and
;       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
;
; CALLING SEQUENCE:
;
;       object = obj_new('FCSDATA',DATADIRECTORY = datadirectory)
;
;
;###############################################################################
; 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 fcsdata::cleanup
ptr_free,self.distPtr
ptr_free,self.GPqPtr
ptr_free,self.monPtr
ptr_free,self.qPtr
ptr_free,self.GPdetPtr
ptr_free,self.vanefficiencyPtr
ptr_free,self.detPtr
ptr_free,self.dettypePtr
ptr_free,self.incdetPtr
ptr_free,self.grpdetPtr
ptr_free,self.timePtr
ptr_free,self.treatmentPtr
ptr_free,self.dataPtr
ptr_free,self.errorPtr
ptr_free,self.rawdataPtr
ptr_free,self.rawerrorPtr
ptr_free,self.energyPtr
ptr_free,self.qmatPtr
ptr_free,self.n_phi_t_ptr
ptr_free,self.ne_phi_t_ptr
ptr_free,self.n_phi_e_ptr
ptr_free,self.ne_phi_e_ptr
ptr_free,self.s_phi_e_ptr
ptr_free,self.se_phi_e_ptr
ptr_free,self.s_q_e_ptr
ptr_free,self.se_q_e_ptr
ptr_free,self.g_e_ptr
ptr_free,self.ge_e_ptr
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::quit
widget_control,self.tlb,/destroy
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function fcsdata::loaddave,filename = filename,group_leader = group_leader
;restore a dave file.
if n_elements(group_leader) eq 0 then group_leader = 0L

catch,errno
if (errno ne 0) then begin
  catch,/cancel
  void = dialog_message(/error,["Problem with your DAVE _raw file?",!err_string])
  return,-1
endif

if n_elements(filename) eq 0 then begin
  filename = DIALOG_PICKFILE(title = 'Select a DAVE "raw" file to reduce',$
                             /read,filter = '*.dave',$
                             path = self.datDir, $
                             dialog_parent = group_leader)
endif
if ((file_test(filename) eq 0) or (filename eq '')) then return,-1
	;strip file list down to the files:
        ;f2=StrSplit(filename,'/',/Extract)
        ;filename=f2[N_elements(f2)-1]
  restore,filename
  self.filename=filename
  xx=(*(*(*davePtr).dataStrPtr).commonStr.histPtr).qty
  ee=(*(*(*davePtr).dataStrPtr).commonStr.histPtr).err
  *self.monPtr = *(*(*(*davePtr).datastrptr).specificptr).monptr
  *self.rawerrorPtr = [[sqrt(*self.monPtr)],[ee]]
  *self.rawdataPtr = [[(*self.monPtr)],[xx]]
  *self.n_phi_t_ptr = xx[*,1:(size(xx))[2]-1]
  *self.ne_phi_t_ptr = ee[*,1:(size(xx))[2]-1]
  *self.n_phi_e_ptr = *self.n_phi_t_ptr
  *self.ne_phi_e_ptr = *self.ne_phi_t_ptr
  *self.timePtr = (*(*(*davePtr).dataStrPtr).commonStr.histPtr).x
  dets=(*(*(*davePtr).dataStrPtr).commonStr.histPtr).y
  *self.detPtr = dets
  self.wavelength = (*(*(*davePtr).dataStrPtr).specificPtr).wavelength
  ;*self.detPtr = *(*(*(*davePtr).dataStrPtr).specificPtr).detectoranglesptr
  self.temperature=(*(*davePtr).descriptr).qty
  *self.dettypePtr=*(*(*(*davePtr).datastrptr).specificptr).dettypeptr
  *self.distPtr=*(*(*(*davePtr).datastrptr).specificptr).distptr
  if ptr_valid(((*(*davePtr).datastrptr).commonstr).treatmentptr) then begin
     (*self.treatmentptr)=*((*(*davePtr).datastrptr).commonstr).treatmentptr
  endif
  *self.dataPtr = [[(*self.monPtr)],[xx]]
  *self.errorPtr = [[sqrt(*self.monPtr)],[ee]]
  self.ndets=(size(*self.dataPtr))[2]
  self.nchan=(size(*self.dataPtr))[1]
  self.run= (*(*(*davePtr).datastrptr).specificptr).run
  self.time=(*(*(*davePtr).datastrptr).specificptr).time
  self.user= (*(*(*davePtr).datastrptr).specificptr).user
  self.comments= (*(*(*davePtr).datastrptr).specificptr).comments
  self.chanwidth= (*(*(*davePtr).datastrptr).specificptr).chanwidth
  self.chop_pulses= (*(*(*davePtr).datastrptr).specificptr).chop_pulses
  self.chop_freq= (*(*(*davePtr).datastrptr).specificptr).chop_freq
  self.monitor= (*(*(*davePtr).datastrptr).specificptr).monitor
  self.date= (*(*(*davePtr).datastrptr).specificptr).date
self.dataLoaded = 1
   ptr_free,(*(*davePtr).datastrptr).commonstr.histptr
   ptr_free,(*(*davePtr).datastrptr).commonstr.commonstrptr
   ptr_free,(*(*davePtr).datastrptr).commonstr.treatmentptr
   ptr_free,(*(*(*davePtr).datastrptr).specificptr).dettypeptr
   ptr_free,(*(*(*davePtr).datastrptr).specificptr).detectoranglesptr
   ptr_free,(*(*(*davePtr).datastrptr).specificptr).distptr
   ptr_free,(*(*(*davePtr).datastrptr).specificptr).monptr
   ptr_free,(*(*davePtr).datastrptr).specificptr
   ptr_free,(*davePtr).datastrptr
   ptr_free,(*davePtr).descriptr
   ptr_free,davePtr
  *self.gpdetptr=*self.detptr
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function fcsdata::loaddata,filename = filename,group_leader = group_leader
if n_elements(group_leader) eq 0 then group_leader = 0L

catch,errno
if (errno ne 0) then begin
  catch,/cancel
  void = dialog_message(/error,["Problem with your ascii .TOF file?",!err_string])
  return,-1
endif

if n_elements(filename) eq 0 then begin
  filename = DIALOG_PICKFILE(title = 'Select an FCS data file to reduce',$
                             /read,filter = '*.tof',$
                             path = self.datDir, $
                             dialog_parent = group_leader)
endif
if ((file_test(filename) eq 0) or (filename eq '')) then return,-1
self.filename = filename
; Read in the header
openr, lun, filename, /get_lun
dumstr = ""
readf,lun,dumstr
colonLocation = strpos(dumstr,':')
self.user = strmid(dumstr,colonLocation+1)
readf,lun,dumstr
colonLocation = strpos(dumstr,':')
self.run = strmid(dumstr,colonLocation+1)

readf,lun,dumstr
colonLocation = strpos(dumstr,':')
self.date = strmid(dumstr,colonLocation+1)

readf,lun,dumstr
colonLocation = strpos(dumstr,':')
self.comments = strmid(dumstr,colonLocation+1)

readf,lun,dumstr
colonLocation = strpos(dumstr,':')
self.time = double(strmid(dumstr,colonLocation+1))

readf,lun,dumstr
colonLocation = strpos(dumstr,':')
self.monitor = double(strmid(dumstr,colonLocation+1))

readf,lun,dumstr
colonLocation = strpos(dumstr,':')
self.chop_pulses = double(strmid(dumstr,colonLocation+1))

readf,lun,dumstr
colonLocation = strpos(dumstr,':')
self.wavelength = double(strmid(dumstr,colonLocation+1))

readf,lun,dumstr
colonLocation = strpos(dumstr,':')
self.chanwidth = double(strmid(dumstr,colonLocation+1))

readf,lun,dumstr
colonLocation = strpos(dumstr,':')
self.chop_freq = double(strmid(dumstr,colonLocation+1))

readf,lun,dumstr
colonLocation = strpos(dumstr,':')
self.temperature = double(strmid(dumstr,colonLocation+1))

readf,lun,dumstr
readf,lun,dumstr

; Read in the detectors
ndet = 107
angle = dblarr(ndet)
group = intarr(ndet)
distance = dblarr(ndet)
type = strarr(ndet)
detector = intarr(ndet)

det = 0.
ang = 0.0
dist = 0.0
typ = ''
grp = 0.

; Parse the string
for i = 0,ndet-1 do begin
  dumstr = ''
  readf,lun,dumstr
  thisLine = strtrim(dumstr,2)
  reads,thisLine,det,ang,dist,dumstr
  newstr = strtrim(dumstr,2)
  spPos = strpos(newstr,' ')
  typ = strmid(newstr,0,1)
  grp = fix(strmid(newstr,spPos+1))
  detector[i] = det
  angle[i] = ang
  distance[i] = dist
  group[i] = grp
  type[i] = typ
endfor



; Ok, now determine the values of each of the detector groups.  Note that there
; are some detectors that are in the same "GROUP".  These should be averaged
; together to yield an angle array smaller than NDET.

phi = angle[0]
d = distance[0]
dettyp=strarr(group[ndet-1])

for i=1,ndet-1 do begin
  ok = where(group eq group[i],count)
  phi = [phi,1.0*total(angle[ok])/count]
  d = [d,1.0*total(distance[ok])/count]
  dettyp[group[i]-2]=type[i]					; keep track of the type of detector
endfor
*self.dettypeptr=dettyp

phi = phi[uniq(phi)]
ndets = n_elements(phi)
*self.detPtr = reverse(phi[0:ndets-2])
self.ndets = ndets - 1
*self.distPtr = d

nchan = 1024

self.nchan = nchan
self.ndets = ndets
grp=1+indgen(ndets)
temp=fltarr(nchan,1)
tof_data=fltarr(nchan,ndets)

; Read raw counts
readf,lun, dumstr
temp1=0
dumnum=0
for i=0,ndets-1 do begin
    readf,lun, dumnum
	readf, lun, temp		;Read in the blocks of data into a temp file
	tof_data[*,i] = temp	;Transfer the temp data into the large array
endfor
free_lun,lun

*self.monPtr = tof_data[*,0]
*self.dataPtr = tof_data
whereZero = where(tof_data eq 0.0,countZero)
error = sqrt(tof_data)
if countZero gt 0 then error[whereZero] = 1.0
*self.errorPtr = error

*self.timePtr = 0.5*self.chanwidth+self.chanwidth*findgen(nchan)
ok = where(*self.dataPtr gt 0.0,count)


;Having determined detector type,
; now need to normalise to detector height *self.dettypeptr

halfheight=where(*self.dettypeptr eq 'A')
newdata=(*self.dataptr)[*,halfheight]*2
newerr=(*self.errorptr)[*,halfheight]*2
(*self.dataPtr)[*,(n_elements(*self.dettypeptr)-n_elements(halfheight))-1:ndets-2]=newdata
(*self.errorPtr)[*,(n_elements(*self.dettypeptr)-n_elements(halfheight))-1:ndets-2]=newerr


*self.n_phi_t_ptr = (*self.dataPtr)[*,1:ndets-1]
*self.ne_phi_t_ptr = sqrt((*self.errorPtr)[*,1:ndets-1])
*self.n_phi_t_ptr = reverse(*self.n_phi_t_ptr,2)
*self.ne_phi_t_ptr = reverse(*self.ne_phi_t_ptr,2)

*self.n_phi_e_ptr = *self.n_phi_t_ptr
*self.ne_phi_e_ptr = *self.ne_phi_t_ptr

*self.rawdataptr=[[(*self.dataPtr)[*,0]],[(*self.dataPtr)[*,1:ndets-1]]]
*self.rawerrorptr=[[sqrt((*self.dataPtr)[*,0])],[(*self.errorPtr)[*,1:ndets-1]]]

self.dataLoaded = 1
*self.gpdetptr=*self.detptr

return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::removeDarkCounts,lo_chan = lo_chan,hi_chan = hi_chan
data = *self.dataPtr
dsize = size(data)
nchan = dsize[1]
ndets = dsize[2]

;hard wire in a better background choice
;just use detectors 40-62
;scan for lowest data that is not 1?
;somewhat problematic really.... rely on visual inspection????

for i = 0,ndets-1 do begin
  bg = (moment(data[lo_chan:hi_chan,i]))[0]
  data[*,i] = data[*,i] - bg
endfor
*self.dataPtr = data > 0.0
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::calcQMatrix
; Compute the Q-matrix for use later
lambda = self.wavelength
y = *self.detPtr	&	ny = n_elements(y)
x = *self.energyPtr	&	nx = n_elements(x)
uy = 1+bytarr(ny)
ux = 1+bytarr(nx)
xmat = x#uy
phimat = ux#y
eo = 81.81/lambda^2
alpha = 2.0717
qmat = sqrt((2.0*eo-xmat-2.0*sqrt(eo*(eo-xmat))*cos(phimat*!dtor))/alpha)
*self.qmatPtr = qmat
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function fcsdata::timeToEnergy,	elastic_channel = elastic_channel, $
								lo_chan = lo_chan, $
								hi_chan = hi_chan, $
								sigma = sigma

;make sure start with the original data every time!

catch,errno
if (errno ne 0) then begin
  catch,/cancel
  void = dialog_message(/error,["Something is wrong in TimeToEnergy conversion",!err_string])
  return,-1
endif


ndets=(size(*self.rawdataPtr))[2]
*self.n_phi_t_ptr =(*self.rawdataPtr)[*,1:(ndets-1)]
*self.ne_phi_t_ptr =(*self.rawerrorPtr)[*,1:(ndets-1)]
*self.n_phi_t_ptr = reverse(*self.n_phi_t_ptr,2)
*self.ne_phi_t_ptr = reverse(*self.ne_phi_t_ptr,2)
*self.n_phi_e_ptr =*self.n_phi_t_ptr
*self.ne_phi_e_ptr =*self.ne_phi_t_ptr
;if n_elements(*self.gpdetptr) gt 0 then *self.detptr=*self.gpdetptr
*self.detptr=*self.gpdetptr

;put the reordered data back into self.data
*self.dataPtr = [[*self.monPtr],[*self.n_phi_t_ptr]]
*self.errorPtr = [[sqrt(*self.monPtr)],[*self.ne_phi_t_ptr]]

if n_elements(sigma) eq 0 then sigma = 1.0
; Remove dark counts
self->removeDarkCounts,lo_chan = lo_chan,hi_chan = hi_chan

data = *self.dataPtr
error = *self.errorPtr
time = *self.timePtr
lambda = self.wavelength

; Determine time at which peak occurs in the monitor spectrum
lsd = 1.e-3*(moment(*self.distPtr))[0]	; average sample-detector distance in m
dms=0.483				; monitor-sample distance in m
h=0.457					; detector height in m
vo = 3.958e3/lambda		; neutron speed in m/s
hbarc = 1973e3			; meV.angstroms
mnc2 = 939.6e9			; meV

eo = 4.0*!dpi^2*(hbarc^2)/(2.0*mnc2*lambda^2)	; meV

; Next let's determine the approximate location of the
; elastic peak in time in each of the detectors.  What we'll
; do is search for the channel with the most counts in it and
; then assign that channel number as the elastic channel.
;dsize = size(data)
;nchan = dsize[1] & ndets = 44; dsize[2] ; CB MADE THIS TO 44 DETECTORS
;telastic = fltarr(ndets-1)
;datmax = fltarr(ndets)
;bg = fltarr(ndets)
;imax = intarr(ndets)
;ratio = fltarr(ndets-1)
;num = 200
;nchan = self.nchan
;for i = 1,ndets-1 do begin
;  bg[i-1] = (moment(data[nchan-num:nchan-1,i]))[0]
;  datmax[i] = max(data[0:600,i],i_max)
;  ratio[i-1] = bg[i-1]/datmax[i]
;  imax[i] = i_max
;endfor
;
;; Where is the ratio less than 0.1 ?
;ok = where(ratio le 0.1,count)
;
; if (count le 1) then elastic_channel=400
;dmax = max(datmax,i_max)
;;imax = imax[i_max]					; elastic channel identified!
;imax_index = where(imax gt 0.0)
;if n_elements(elastic_channel) eq 0 then begin
;  imax = fix((moment(imax[imax_index]))[0])
;endif else begin
;  imax = elastic_channel
;endelse

imax=self->findelastic()
if imax eq -1 then imax=400

telastic = time[imax]				; elastic time in micro-seconds
dt_elastic = lsd/vo					; elastic time difference in seconds
tsd = 1.0e-6*time - (1.0e-6*telastic-dt_elastic)	; time in seconds
c = 3.0e8	; m/s
vn = lsd/tsd

ef = 0.5*mnc2*(vn/c)^2
de = eo - ef
*self.energyPtr = de
; Limit the energy range
;eok = where(*self.energyPtr le 5.0 and *self.energyPtr ge -50.,countE)
*self.energyPtr = (*self.energyPtr);[eok]
; Finally we have to correct the scattering data for the Jacobian of the
; nonlinear transformation from time to energy.
dsize = size(data)
nchan = n_elements(*self.energyPtr)
ndets = dsize[2]
u = 1+bytarr(ndets-1)
jacobian = abs(((0.5*tsd/ef)#u))
newdata = (data[*,1:ndets-1]*jacobian)
newerror = (error[*,1:ndets-1]*jacobian)




; Now must put in the final correction to go from the partial differential
; cross-section to the scattering function.


; need other scaling factors:
; TXFLUX=RAWMON/beamarea/EFFICIENCY
; XBUFF=4.0*PI/AVOGADRO/RMOLES/SIGMA/DELW/RAWDELTAT/TXFLUX/TSO/RNEUTMASS/RLENGTH/RLENGTH*RJMEV
; HISTO_DATA(IA,J,I)=XBUFF*HISTO_DATA(IB,J,I)
; where:
; Avogadro= 0.602252 	;Avogadro's number (atoms/mole) x 1.E-24
; DELW/2.222E-03 	;18 in. detector solid angle
; RJMEV/0.160217733 	;conversion factor (J/meV) x 1.E+21
; RNEUTMASS/1.6749286 	;neutron mass (kg) x 1.E27
; RLENGTH/2284.8 	;sample-to-detector distance (mm)
; beamarea/25.8/ 	;beam cross-sectional area (cm**2)
; efficiency=3.3e-04*wavelength(1)/4.08
;
;
; mostly just scaling factors, shall ignore here for now. except sigma.
;


ko = 2.0*!dpi/lambda
kf = mnc2*lsd/(hbarc*c*tsd)
term = ((4.0*!dpi*ko/kf)/sigma)#u
newdata = newdata*term
newerror = newerror*term

*self.dataPtr = newdata/self.monitor
*self.errorPtr = newerror/self.monitor

; Use an arbitrary multiplicative factor to get the
; scattering values into a reasonable range.
*self.dataPtr = 1.e6*(*self.dataPtr)
*self.errorPtr = 1.e6*(*self.errorPtr)

; Determine the wavevector transfer
*self.qPtr = 4.0*!pi*sin(0.5*(!dtor)*(*self.detPtr))/self.wavelength

*self.s_phi_e_ptr = *self.dataPtr
*self.se_phi_e_ptr = *self.errorPtr
*self.gpqptr=*self.qptr

self->calcQMatrix

return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::symmetrize
T = self.temperature
e = *self.energyPtr
data = *self.s_phi_e_ptr
error = *self.se_phi_e_ptr

uchan = 1+bytarr(self.nchan)
udet  = 1+bytarr(self.ndets)
symfact = (exp(0.5*11.6*e/T))#udet
*self.s_phi_e_ptr = data*symfact
*self.se_phi_e_ptr = error*symfact
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function fcsdata::findelastic, elastic_channel=elastic_channel
data = *self.dataPtr
; Let's determine the approximate location of the
; elastic peak in time in each of the detectors.  What we'll
; do is search for the channel with the most counts in it and
; then assign that channel number as the elastic channel.


catch,errno
if (errno ne 0) then begin
  catch,/cancel
  void = dialog_message(/error,["I cant seem to find your elastic peak (or other error)"$
						,!err_string])
  return,-1
endif


;ignore the first detectors, always bad, 2:45 usually good, so newdets=44
newdets=44
telastic = fltarr(newdets-1)
time = *self.timePtr
datmax = fltarr(newdets)
bg = fltarr(newdets-1)
imax = intarr(newdets)
ratio = fltarr(newdets-1)
num = 200
nchan = self.nchan
for i = 1,newdets-1 do begin
  bg[i-1] = (moment(data[nchan-num:nchan-1,i]))[0]
  datmax[i] = max(data[0:600,i],i_max)
  ratio[i-1] = bg[i-1]/datmax[i]
  imax[i] = i_max
endfor

; Where is the ratio less than 0.01 ?
ok = where(ratio le 0.01,count)

dmax = max(datmax,i_max)
;imax = imax[i_max]					; elastic channel identified!
imax_index = where(imax gt 0.0)
if n_elements(elastic_channel) eq 0 then begin
  ;imax = fix((moment(imax[imax_index]))[0])
   imax=fix(median(imax[imax_index]))		;median seems to be working OK-better than moment
endif else begin
  imax = elastic_channel
endelse
print,'Elastic Channel calculated to be:',imax

return,imax
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function fcsdata::normalize, factor=factor
;normalize to vanadium elastic peaks intensity for detector norm.

uy=1+bytarr((size(*self.n_phi_t_ptr,/dimension))[0])
fmat=uy#factor
*self.s_phi_e_ptr=*self.s_phi_e_ptr/fmat
*self.se_phi_e_ptr=*self.se_phi_e_ptr/fmat

return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function fcsdata::detectoreff
;    This is based on Treat by T.Udovic et. al.
;         EFFARRAY(0)->EFFARRAY(500) CORRESPONDS TO 0->50.0 MEV BY .1 MEV STEPS
;         EFFARRAY(500)->EFFARRAY(700) CORRESPONDS TO 50->250 MEV BY 1 MEV STEPS
;         ABOVE 250 MEV, EFF IS APPROXIMATED BY AN EQUATION .
;         Detector efficiencies were calculated by numerical integration
;          (101 steps) using program Detefftreat written by T. Udovic. It
;          was found that the effect of sample and detector heights on
;          efficiency was insignificant and could be ignored.
;
;
; CMB took T.Udovics integrated data and fitted that to a 5 parameter exponential
; y= a0 + a1*exp(-a2*x) + a3*exp(-a4*x)
;
; for a very good fit (better than 0.5% difference), the params are

; 0.2801, 0.2739, 0.0697, 0.4489, 0.0103


x=*self.energyptr
y=*self.s_phi_e_ptr
erry=*self.se_phi_e_ptr

; GET ALL THE DETECTORS ON THE SAME SIZE SCALE DEPENDING ON TYPE:

;y(1024,63) x is (1024)

uy=1+bytarr((size(*self.s_phi_e_ptr,/dimension))[1])
xmat=x#uy

*self.s_phi_e_ptr = y / (0.2801 + 0.2739*exp(-0.0697*abs(xmat)) $
			+ 0.4489*exp(-0.0103*abs(xmat)))
*self.se_phi_e_ptr = erry / (0.2801 + 0.2739*exp(-0.0697*abs(xmat)) + $
			0.4489*exp(-0.0103*abs(xmat)))

return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function fcsdata::fitvanelastic, elastic_channel= elastic_channel, width=width, dt = dt

;vaneff only valid if between limits 0.5 to 1.5

catch,errno
if (errno ne 0) then begin
  catch,/cancel
  void = dialog_message(/error,["I cant fit the Vanadium elastic peak (or other error)"$
						,!err_string])
  return,-1
endif


data = *self.n_phi_t_Ptr
time = *self.timePtr
time_channel=elastic_channel*dt
twidth = width * dt
coeffs=fltarr(n_elements(*self.detPtr))			;already removed the monitor
guess=[100,time_channel,twidth,1.0]
for i=0,(n_elements(*self.detPtr)-1) do begin
   yfit=GAUSSFIT(time[elastic_channel-30:elastic_channel+30],$
	data[elastic_channel-30:elastic_channel+30,i],A,estimates=guess,nterms=4)
   coeffs(i)=A[0]
;   window
;   plot,time[elastic_channel-30:elastic_channel+30],data[elastic_channel-30:elastic_channel+30,i]
;   oplot,time[elastic_channel-30:elastic_channel+30],yfit
endfor

median=median(coeffs[where(coeffs gt 0.0)])
vanefficiency=coeffs/median

;make a change so that
;vaneff only valid if between limits
;make those limits 0.5 to 1.5

whereless= where(vanefficiency lt 0.5, countless)
if countless gt 0 then vanefficiency[whereless] = 1.0
wheremore= where(vanefficiency gt 1.5, countmore)
if countmore gt 0 then vanefficiency[wheremore] = 1.0

*self.vanefficiencyPtr=abs(vanefficiency)

return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

pro fcsdata::calcDOS
sqe = *self.s_phi_e_ptr
sqerr = *self.se_phi_e_ptr
e = *self.energyPtr
q = *self.qmatPtr
data = sqe/q^2
error = sqerr/q^2
qsize = size(q)
incdet=*self.incdetptr

nx = qsize[1]
ny = qsize[2]
dos = fltarr(nx)
doserr = fltarr(nx)

delta=intarr(ny)

for i=0,n_elements(incdet)-1 do begin
  delta[incdet[i]]=1
endfor

for i=0,ny-1 do begin
   data[*,i]=data[*,i]*delta[i]
   error[*,i]=error[*,i]*delta[i]
endfor

for i = 0,nx-1 do begin
  dos[i] = int_tabulated(q[i,*],data[i,*])
  dq = reform(q[i,1:ny-1] - q[i,0:ny-2])
  dq = [dq[0],dq]
  doserr[i]  = (sqrt(total(error[i,*]*error[i,*]*dq)))
endfor

dos = (1.0-exp(-e*11.6/self.temperature))*e*dos
doserr = (1.0-exp(-e*11.6/self.temperature))*e*doserr
*self.g_e_ptr = dos
*self.ge_e_ptr = doserr
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::average_dets

;get detectors to include:
include=*self.incdetptr
;get q groups to include:
if n_elements(*self.grpdetptr) gt 0 then begin
   groups=*self.grpdetptr
endif else begin
   groups=intarr(63,63)
   for j=0,62 do begin
     groups[j,j]=1
   endfor
endelse


;at this point, we have n_phi_t, n_phi_e, s_phi_e but not s_q_e
;take the included detectors and the number of group matrix, and
;recreate this data, in the averaged arrays
;
;
;self.GPqPtr..63  double q values
;self.GPdetPtr..63 angles
;GPn_phi_t_ptr..1024,63 float
;GPne_phi_t_ptr..
;GPn_phi_e_ptr..
;GPne_phi_e_ptr..
;GPs_phi_e_ptr..
;GPse_phi_e_ptr..
;
;

GPdet=fltarr(((size(groups))[1]))    ;equal weighting of detector angles:


uy=1+bytarr((size(groups))[1])
xmat=uy#*self.detptr
temp=xmat*groups
for i=0,((size(groups))[1]-1) do begin
 gpdet[i]=total(temp(i,*))/n_elements(where(groups[i,*] ne 0))
endfor


*self.detptr=gpdet

GPq=fltarr(((size(groups))[1]))    ;equal weighting of 	q (as fast as recalculating?)
xmat=uy#*self.gpqptr
temp=xmat*groups
for i=0,((size(groups))[1]-1) do begin
 gpq[i]=total(temp(i,*))/n_elements(where(groups[i,*] ne 0))
endfor


;*self.qptr=gpq


GPne_phi_t=*self.ne_phi_t_ptr^2
GPn_phi_t=*self.n_phi_t_ptr#transpose(groups)
GPne_phi_t=GPne_phi_t#transpose(groups)
GPne_phi_t=sqrt(GPne_phi_t)
for i=0,((size(groups))[1])-1 do begin
   GPn_phi_t[*,i]=GPn_phi_t[*,i]/n_elements(where(groups[i,*] ne 0))
   GPne_phi_t[*,i]=GPne_phi_t[*,i]/n_elements(where(groups[i,*] ne 0))
endfor


GPne_phi_e=*self.ne_phi_e_ptr^2
GPn_phi_e=*self.n_phi_e_ptr#transpose(groups)
GPne_phi_e=GPne_phi_e#transpose(groups)
GPne_phi_e=sqrt(GPne_phi_e)
for i=0,((size(groups))[1])-1 do begin
   GPn_phi_e[*,i]=GPn_phi_e[*,i]/n_elements(where(groups[i,*] ne 0))
   GPne_phi_e[*,i]=GPne_phi_e[*,i]/n_elements(where(groups[i,*] ne 0))
endfor


GPse_phi_e=*self.se_phi_e_ptr^2
GPs_phi_e=*self.s_phi_e_ptr#transpose(groups)
GPse_phi_e=GPse_phi_e#transpose(groups)
GPse_phi_e=sqrt(GPse_phi_e)
for i=0,((size(groups))[1])-1 do begin
   GPs_phi_e[*,i]=GPs_phi_e[*,i]/n_elements(where(groups[i,*] ne 0))
   GPse_phi_e[*,i]=GPse_phi_e[*,i]/n_elements(where(groups[i,*] ne 0))
endfor


*self.n_phi_t_ptr=gpn_phi_t
*self.ne_phi_t_ptr=gpne_phi_t
*self.n_phi_e_ptr=gpn_phi_e
*self.ne_phi_e_ptr=gpne_phi_e
*self.s_phi_e_ptr=gps_phi_e
*self.se_phi_e_ptr=gpse_phi_e

self->calcQMatrix


end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::draw_n_phi_t,_Extra = extra,detector = detector
thisformat = '(f8.2)'
x = *self.timePtr
y = (*self.n_phi_t_ptr)[*,detector]
yerr = (*self.ne_phi_t_ptr)[*,detector]
if extra.autoScale eq 1 then begin
  dy = 0.1*(max(y)-0.0)
  ymin = min(y)-dy
  yrange = [ymin < 0,max(y)+dy]
  xrange = [min(x),max(x)]
  errorflag=extra.errorflag
endif else begin
  xrange = extra.xrange
  yrange = extra.yrange
  errorflag=extra.errorflag
endelse
thisTitle = '!3Scattering angle = '+ $
	strtrim(string((*self.detPtr)[detector],thisFormat),2) + $
	'!E0!N'
plot,x,y,psym = 4,title = thisTitle,xrange = xrange, $
	yrange = yrange,/xsty,/ysty,ytitle = '!3N(!4u!3,t)', xtitle = '!3time (!4l!3s)'
if (errorflag ne 0) then begin
 errplot,x,y-yerr,y+yerr,width = 0.0
endif
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::draw_n_phi_e,_Extra = extra,detector = detector
thisformat = '(f8.2)'
x = *self.energyPtr
y = (*self.n_phi_e_ptr)[*,detector]
yerr = (*self.ne_phi_e_ptr)[*,detector]
if extra.autoScale eq 1 then begin
  dy = 0.1*(max(y)-0.0)
  ymin = min(y)-dy
  yrange = [ymin < 0,max(y)+dy]
  xrange = [min(x),max(x)]
  errorflag=extra.errorflag
endif else begin
  xrange = extra.xrange
  yrange = extra.yrange
  errorflag=extra.errorflag
endelse
thisTitle = '!3Scattering angle = '+ $
	strtrim(string((*self.detPtr)[detector],thisFormat),2) + $
	'!E0!N'
plot,x,y,psym = 4,title = thisTitle,xrange = xrange, $
	yrange = yrange,/xsty,/ysty,ytitle = '!3N(!4u!3,E)', xtitle = 'E (meV)'
if (errorflag ne 0) then begin
 errplot,x,y-yerr,y+yerr,width = 0.0
endif
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::draw_s_phi_e,_Extra = extra,detector = detector
thisformat = '(f8.2)'
x = *self.energyPtr
y = (*self.s_phi_e_ptr)[*,detector]
yerr = (*self.se_phi_e_ptr)[*,detector]
if extra.autoScale eq 1 then begin
  dy = 0.1*(max(y)-0.0)
  ymin = min(y)-dy
  yrange = [ymin < 0,max(y)+dy]
  xrange = [min(x),max(x)]
  errorflag=extra.errorflag
endif else begin
  xrange = extra.xrange
  yrange = extra.yrange
  errorflag=extra.errorflag
endelse

thisTitle = '!3Scattering angle = '+ $
	strtrim(string((*self.detPtr)[detector],thisFormat),2) + $
	'!E0!N'
plot,x,y,psym = 4,title = thisTitle,xtitle = '!3E (meV)', $
	ytitle = '!3S(!4u!3,E)',xrange = xrange,yrange = yrange,/xsty,/ysty
if (errorflag ne 0) then begin
 errplot,x,y-yerr,y+yerr,width = 0.0
endif
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::draw_s_q_e,_Extra = extra,detector = detector
thisformat = '(f8.2)'
x = *self.energyPtr
y = (*self.s_q_e_ptr)[*,detector]
yerr = (*self.se_q_e_ptr)[*,detector]
if extra.autoScale eq 1 then begin
  dy = 0.1*(max(y)-0.0)
  ymin = min(y)-dy
  yrange = [ymin < 0,max(y)+dy]
  xrange = [min(x),max(x)]
  errorflag=extra.errorflag
endif else begin
  xrange = extra.xrange
  yrange = extra.yrange
  errorflag=extra.errorflag
endelse

qVal = (*self.qPtr)[detector]
thisTitle = '!3Q!D!N='+strtrim(string(qVal,thisFormat),2)+ $
	' !3!sA!r!u!9 %!3!n!E-1!N'
plot,x,y,psym = 4,title = thisTitle,xtitle = '!3E (meV)', $
	ytitle = '!3S(Q,E)',xrange = xrange,yrange = yrange,/xsty,/ysty
if (errorflag ne 0) then begin
 errplot,x,y-yerr,y+yerr,width = 0.0
endif
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::draw_dos,_Extra = extra
x = *self.energyPtr
y = *self.g_e_ptr
yerr = *self.ge_e_ptr
if extra.autoScale eq 1 then begin
  dy = 0.1*(max(y)-min(y))
  ymin = min(y)-dy
  yrange = [ymin < 0,max(y)+dy]
  xrange = [min(x),max(x)]
  errorflag=extra.errorflag
endif else begin
  xrange = extra.xrange
  yrange = extra.yrange
  errorflag=extra.errorflag
endelse
thisTitle = '!3Density of States'
plot,x,y,psym = 4,title = thisTitle, xtitle = '!3E (meV)', $
	ytitle = '!3G(E)',xrange = xrange,yrange = yrange,/xsty,/ysty
if (errorflag ne 0) then begin
 errplot,x,y-yerr,y+yerr,width = 0.0
endif
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::draw_elastic,_Extra = extra

data = (*self.s_phi_e_ptr)
energy = (*self.energyptr)

fitlow=-1.0
fithi=1.0

lowvec=where((energy) le fitlow,count)
if count ge 0 then begin
  low=lowvec[count-1]
endif else begin
   low=350
endelse

highvec=where((energy) ge fithi,count)
if count ge 0 then begin
  hi=highvec[count-1]
endif else begin
  hi=600
endelse

intcoeffs=fltarr(n_elements(*self.detPtr))
bgcoeffs=fltarr(n_elements(*self.detPtr))
guess=[1.0,0.001,0.1,1.0]

for i=0,(n_elements(*self.detPtr)-1) do begin
   yfit=GAUSSFIT(energy[low:hi],$
	data[low:hi,i],A,estimates=guess,nterms=4)
   intcoeffs(i)=A[0]
   bgcoeffs(i)=A[3]
endfor

x = *self.detPtr
y1 = intcoeffs
y2 = bgcoeffs


if extra.autoScale eq 1 then begin
  dy = 0.1*(max(y1)-min(y2))
  ymin = min(y2)-dy
  yrange = [ymin < 0,max(y1)+dy]
  xrange = [min(x),max(x)]
  errorflag=extra.errorflag
endif else begin
  xrange = extra.xrange
  yrange = extra.yrange
  errorflag=extra.errorflag
endelse
thisTitle = '!3Integrated Elastic Intensity'
plot,x,y1,psym = -4,title = thisTitle, xtitle = '!3Detector Angle', $
	ytitle = '!3Gaussian Area',xrange = xrange,yrange = yrange,/xsty,/ysty
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::image_s_q_e,_Extra = extra
x = *self.energyPtr
y = *self.qmatPtr
ysize = size(y)
ny = ysize[2]
uy = 1+bytarr(ny)
xmat = uy#x
ymat = transpose(y)
z = transpose(*self.s_phi_e_ptr) < extra.image_max
if extra.autoScale eq 1 then begin
  yrange = [min(x),max(x)]
  xrange = [min(y),max(y)]
endif else begin
  xrange = extra.xrange
  yrange = extra.yrange
endelse
contour,z,ymat,xmat,/fill, $
	nlevels = 30,ytitle = 'E (meV)',xtitle = 'Q (!3!sA!r!u!9 %!3!n!E-1!N)', $
	xrange = xrange,yrange = yrange,/xsty,/ysty
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::image_s_phi_e,_Extra = extra
x = *self.energyPtr
y = *self.detPtr
z = *self.s_phi_e_ptr
if extra.autoScale eq 1 then begin
  yrange = [min(x),max(x)]
  xrange = [min(y),max(y)]
endif else begin
  xrange = extra.xrange
  yrange = extra.yrange
endelse
zdisp = transpose(z) < extra.image_max
contour,zdisp,y,x,/fill, $
	nlevels = 30,ytitle = '!3E (meV)', $
	xtitle = '!3Scattering angle (!E0!N)', $
	xrange = xrange,yrange = yrange,/xsty,/ysty
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::image_n_phi_e,_Extra = extra
x = *self.energyPtr
y = *self.detPtr
z = *self.n_phi_e_ptr
zsize = size(z)
ny = zsize[2]
uy = 1+bytarr(ny)
nx = zsize[1]
ux = 1+bytarr(nx)
if extra.autoScale eq 1 then begin
  yrange = [min(x),max(x)]
  xrange = [min(y),max(y)]
endif else begin
  xrange = extra.xrange
  yrange = extra.yrange
endelse
contour,(transpose(z) < extra.image_max),y,x,/fill, $
	nlevels = 30,ytitle = '!3E (meV)',xtitle = '!3Scattering angle (!E0!N)', $
	xrange = xrange,yrange = yrange, /xsty,/ysty
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::image_n_phi_t,_Extra = extra
x = *self.timePtr
y = *self.detPtr
z = *self.n_phi_t_ptr
if extra.autoScale eq 1 then begin
  yrange = [min(x),max(x)]
  xrange = [min(y),max(y)]
endif else begin
  xrange = extra.xrange
  yrange = extra.yrange
endelse
contour,(transpose(z) < extra.image_max),y,x,/fill, $
	nlevels = 30,ytitle = '!3time channel', $
	xtitle = '!3Scattering angle (!E0!N)', $
	xrange = xrange,yrange = yrange, /xsty,/ysty
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::draw,_Extra = extra,detector = detector,qSelect = qSelect
if n_elements(*self.energyPtr) ne 0 then $
   x = *self.energyPtr else $
   x = *self.timePtr

if n_elements(detector) eq 0 then detector = 17
detector = (detector lt 0) ? 0:detector
detector = (detector gt self.ndets-1) ? (self.ndets-1):detector
if n_elements(qSelect) eq 0 then qSelect = 0

thisformat = '(f8.2)'
if qSelect then begin
	qVal = (*self.qPtr)[detector]
	thisTitle = '!3Q='+strtrim(string(qVal,thisFormat),2)
endif else begin
	qVal = (4.0*!pi*sin(0.5*((*self.detPtr)[detector])*!dtor))/self.wavelength
	thisTitle = '!3Q!DELASTIC!N='+strtrim(string(qVal,thisFormat),2)
endelse
y = (*self.dataPtr)[*,detector]
yerr = (*self.errorPtr)[*,detector]
xsize = 400 & ysize = 400
plot,x,y,psym = 4,_Extra = extra,title = thisTitle
errplot,x,y-yerr,y+yerr,width = 0.0

return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::drawMonitor,_Extra = extra
x = *self.timePtr
y = *self.monPtr
yerr = sqrt(y)
plot,x,y,psym = 4,_Extra = extra,title = "!6Beam Monitor",$
     xtitle = "!6time (!4l!6s)"
errplot,x,y-yerr,y+yerr,width = 0.0
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::animate,_Extra = extra,duration = duration
if n_elements(duration) eq 0 then duration = 0.01
if n_elements(*self.energyPtr) ne 0 then $
   x = *self.energyPtr else $
   x = *self.timePtr

data = *self.dataPtr
error = *self.errorPtr
thisformat = '(f8.2)'
xsize = 400 & ysize = 400
window,0,xsize = xsize,ysize = ysize
winVis = !d.window
window,/free,/pixmap,xsize = xsize,ysize = ysize
winPix = !d.window

for i = 0,self.ndets-1 do begin
  wset,winPix
  y = data[*,i] & yerr = error[*,i]
  thisTitle = 'Detector angle: '+$
               strtrim(string((*self.detPtr)[i],format = thisFormat),2)
  plot,x,y,psym = 4,_Extra = extra,title = thisTitle
  errplot,x,y-yerr,y+yerr,width = 0.0
  wset,winVis
  device,copy = [0,0,!d.x_size,!d.y_size,0,0,winPix]
  wait,duration
endfor
wdelete,winPix
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::drawimage,threshold = threshold, $
                       colortable = colortable, $
                       _Extra = extra
if n_elements(*self.energyPtr) eq 0 then return

if n_elements(colortable) eq 0 then colortable = 6
thisPosition = [0.15,0.15,0.9,0.90]
chsize = 1.0

y = *self.energyPtr
x = *self.detPtr
z = transpose(*self.dataPtr)

; Now interpolate to a uniform grid so scaling of the image
; is correct.
nx = 70
ny = 500
xlo = min(x) & xhi = max(x)
ylo = min(y) & yhi = max(y)
xx = hfbs_makepoints(xlo = xlo,xhi = xhi,npts = nx)
yy = hfbs_makepoints(xlo = ylo,xhi = yhi,npts = ny)
ix = interpol(indgen(n_elements(x)),x,xx)
iy = interpol(indgen(n_elements(y)),y,yy)
zz = interpolate(z,ix,iy,/grid)
z = zz & y = yy & x = xx

if n_elements(threshold) eq 0 then threshold = max(z)
z = temporary(z) > 0
zdisplay = bytscl(temporary(z) < threshold)


loadct,colortable,/silent
xr = [min(x),max(x)]
yr = [min(y),max(y)]
xir = [0,n_elements(x)-1]
yir = [0,n_elements(y)-1]
tvimage, zdisplay, position=thisPosition,/minus_one,/erase,$
         xrange = xir,yrange = yir,xstyle = 1,ystyle = 1
loadct,0,/silent
plot,x,y,/nodata,/noerase,$
     xstyle = 1,ystyle = 1,$
     position = thisPosition, $
     xrange = xr,yrange = yr, $
     xtitle = '!6Angle (!Eo!N)',ytitle = '!6E (meV)',$
     title = self.filename,charsize = chsize

return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::surfaceObject
FSC_SURFACE, transpose(*self.dataPtr),*self.detPtr,*self.energyPtr,$
             xtitle = 'Detector Angle',ytitle = 'E (meV)',/block,$
             /elevation_shading,/shaded
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::setproperty,ndets = ndets, $
			 wavelength = wavelength, $
             detectors = detectors, $
             filename = filename, $
             rawdata = rawdata , $
             rawerror = rawerror, $
             alldata = alldata , $
             error = error, $
		temperature = temperature, $
			dettype = dettype,$
			time = time,$
			dist = dist,$
			mon = mon,$
			nchan = nchan,$
			run = run,$
			rtime = rtime,$
			user = user,$
			comments = comments,$
			chanwidth = chanwidth,$
			chop_pulses = chop_pulses,$
			chop_freq = chop_freq ,$
			totalmonitor = totalmonitor , $
			date = date,$
			treatment=treatment,$
 			loaded=loaded,$
			incdet=incdet,$
			grpdet=grpdet

;if arg_present(ndets) then self.ndets= ndets
if arg_present(incdet) then *self.incdetptr= incdet
if arg_present(grpdet) then *self.grpdetptr= grpdet
if arg_present(loaded) then self.dataloaded=1
if arg_present(wavelength) then self.wavelength=wavelength
if arg_present(date) then self.date=date
if arg_present(dettype) then *self.dettypeptr= dettype
if arg_present(time) then *self.timeptr= time
if arg_present(dist) then *self.distptr= dist
if arg_present(mon) then *self.monptr= mon
if arg_present(nchan) then self.nchan= nchan
if arg_present(run) then self.run= run
if arg_present(rtime) then self.time= rtime
if arg_present(user) then self.user= user
if arg_present(comments) then self.comments= comments
if arg_present(chanwidth) then self.chanwidth= chanwidth
if arg_present(chop_pulses) then self.chop_pulses= chop_pulses
if arg_present(chop_freq) then self.chop_freq= chop_freq
if arg_present(totalmonitor) then self.monitor= totalmonitor
if arg_present(temperature) then self.temperature= temperature
if arg_present(treatment) then  *self.treatmentptr=treatment

if arg_present(ndets) then begin
  if ndets eq 0 then begin
	ndets=n_elements(dettype)
	self.ndets=ndets
  endif
endif
if arg_present(rawdata) then begin
  *self.rawdataPtr=rawdata
endif
if arg_present(rawerror) then begin
  *self.rawerrorPtr=rawerror
endif
if arg_present(alldata) then begin
  *self.dataPtr= alldata
  *self.n_phi_t_ptr = (*self.dataPtr)[*,1:ndets-1]
  *self.n_phi_t_ptr = reverse(*self.n_phi_t_ptr,2)
  *self.n_phi_e_ptr = *self.n_phi_t_ptr
endif
if arg_present(error) then begin
  *self.errorPtr= error
  *self.ne_phi_t_ptr = (*self.errorPtr)[*,1:ndets-1]
  *self.ne_phi_t_ptr = reverse(*self.ne_phi_t_ptr,2)
  *self.ne_phi_e_ptr = *self.ne_phi_t_ptr
endif
if arg_present(detectors) then *self.detptr= detectors
if arg_present(detectors) then *self.gpdetptr= detectors
if arg_present(filename) then self.filename= filename



end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::getproperty,ndets = ndets, $
			 wavelength = wavelength, $
                         detectors = detectors, $
                         filename = filename, $
                         erange = erange, $
                         datarange = datarange, $
                         alldata = alldata , $
                         rawdata = rawdata , $
			 data = data, $
                         rawerror = rawerror, $
                         error = error, $
                         n_phi_t_max = n_phi_t_max, $
                         n_phi_e_max = n_phi_e_max, $
                         s_phi_e_max = s_phi_e_max, $
                         s_q_e_max = s_q_e_max, $
                         g_e_max = g_e_max, $
                         temperature = temperature, $
                         nqbins = nqbins, $
			 vanefficiency=vanefficiency ,$
			colorTable = colorTable,$
			dettype = dettype,$
			energy = energy,$
			time = time,$
			dist = dist,$
			mon = mon,$
			nchan = nchan,$
			run = run,$
			rtime = rtime,$
			user = user,$
			comments = comments,$
			chanwidth = chanwidth,$
			chop_pulses = chop_pulses,$
			chop_freq = chop_freq ,$
			totalmonitor = totalmonitor , $
			date = date,$
			treatment=treatment

if arg_present(date) then date=self.date
if arg_present(colorTable) then colorTable= self.colorTable
if arg_present(dettype) then dettype= *self.dettypePtr
if arg_present(time) then time= *self.timePtr
if arg_present(dist) then dist= *self.distPtr
if arg_present(mon) then mon= *self.monPtr
if arg_present(nchan) then nchan= self.nchan
if arg_present(run) then run= self.run
if arg_present(rtime) then rtime= self.time
if arg_present(user) then user= self.user
if arg_present(comments) then comments= self.comments
if arg_present(chanwidth) then chanwidth= self.chanwidth
if arg_present(chop_pulses) then chop_pulses= self.chop_pulses
if arg_present(chop_freq) then chop_freq= self.chop_freq
if arg_present(totalmonitor) then totalmonitor= self.monitor
if arg_present(treatment) then begin
  if n_elements(*self.treatmentptr) ne 0 then begin
	 treatment= *self.treatmentptr
  endif else begin
	treatment=''
  endelse
endif

if arg_present(nqbins) then nqbins = n_elements(*self.qPtr)
if arg_present(energy) then energy = *self.energyPtr
if arg_present(temperature) then temperature= self.temperature
if arg_present(rawdata) then rawdata= *self.rawdataPtr
if arg_present(rawerror) then rawerror= *self.rawerrorPtr
if arg_present(data) then data= *self.dataPtr
if arg_present(alldata) then alldata= *self.dataPtr
if arg_present(error) then error= *self.errorPtr
if arg_present(ndets) then ndets= self.ndets
if arg_present(detectors) then detectors= *self.detPtr
if arg_present(filename) then filename= self.filename
if arg_present(vanefficiency) then vanefficiency = *self.vanefficiencyPtr

if arg_present(datarange) then datarange= [min(*self.dataPtr),max(*self.dataPtr)]

if arg_present(erange) then $
	erange = [min(*self.energyPtr),max(*self.energyPtr)]
if arg_present(n_phi_t_max) then $
	n_phi_t_max = max(*self.n_phi_t_ptr)
if arg_present(n_phi_e_max) then $
	n_phi_e_max = max(*self.n_phi_e_ptr)
if arg_present(s_phi_e_max) then $
	s_phi_e_max = max(*self.s_phi_e_ptr)
if arg_present(g_e_max) then $
	g_e_max = max(*self.g_e_ptr)
if arg_present(s_q_e_max) then $
	s_q_e_max = max(*self.s_phi_e_ptr)
if arg_present(wavelength) then $
	wavelength = self.wavelength

end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::saveDave,group_leader = group_leader,davePtr = davePtr, $
	n_phi_t = n_phi_t, $
	n_phi_e = n_phi_e, $
	s_phi_e = s_phi_e, $
	s_q_e = s_q_e, $
	g_e = g_e ,$
        raw=raw,$
	file = file
if n_elements(group_leader) eq 0 then begin
   group_leader = 0L
endif else begin
   widget_control,group_leader,get_uvalue = pState
endelse
if	n_phi_t and n_phi_e and $
	s_phi_e and s_q_e and $
	g_e eq 0 then return
true = 0

if n_elements(file) eq 0 or file eq "" then begin
      fileout = DIALOG_PICKFILE(title = 'Enter the name of the DAVE file',$
                             path=(*!dave_defaults).datdir,$
			     /write,filter = '*.dave',$
                             dialog_parent = group_leader)
     if n_elements(fileout) eq 0 or fileout eq " " then return
endif else begin
	fileout=file
endelse

;check to see what 'outfile' is ...
;f2=StrSplit(fileout[0],'/',/Extract)
;file2=f2[N_elements(f2)-1]
file2=fileout[0]
wheredave2=where(strmatch(file2,'*.dave',/fold_case) eq 1,dave2)

if dave2 ge 1 then begin
   file3=StrSplit(file2,'.dave',/Extract,/fold_case,/regex)
   f3=StrSplit(file3,'_',/Extract)
   if f3(n_elements(f3)-1) eq 'raw' or f3(n_elements(f3)-1) eq 'RAW' then begin
   fileout=strmid(file3,0,(strlen(file3)-4))
   endif
endif else begin
  fileout=file2
endelse




timestamp='Time: '+ strmid(systime(0),4,3)+'-'$
		+strmid(systime(0),8,2)+'-'+strmid(systime(0),20,4)+' '$
		+strmid(systime(0),11,5)

hbar_omega = '!6!sh!r!e/!n !7x!6'

if raw ne 0 then begin
  data = (*self.rawdataptr)[*,1:(size(*self.rawdataptr))[2]-1]
  error = (*self.rawerrorptr)[*,1:(size(*self.rawdataptr))[2]-1]
  x = *self.timePtr
  y = *self.gpdetPtr
  xlabel = 't'
  xunits = 'time:us'
  ylabel = '!4u!3'
  yunits = 'angle:degrees'
  histlabel = 'Raw !3I (!4u!3,t)'
  histunits = 'Raw Counts'
  treatment = ['File:'+ self.filename]
  treatment = [[treatment],timestamp ]
  treatment = [[treatment],'Raw data converted to DAVE']
  treatment = [[treatment],'--------------------------']
  if ptr_valid(self.treatmentPtr) then begin
      if n_elements(*self.treatmentPtr) ne 0 then begin
          treatment=[[treatment],*self.treatmentPtr]
      endif
  endif

  thisFile = fileout + '_raw.dave'

  *(*(*davePtr).dataStrPtr).commonStr.histPtr = $
				      		{qty:data,$
				       		 err:error,$
				      		 x:x,$
				      		 y:y}
 (*(*(*davePtr).dataStrPtr).specificPtr).wavelength=self.wavelength
 *(*(*(*davePtr).dataStrPtr).specificPtr).detectorAnglesptr=*self.gpdetPtr
 *(*(*(*davePtr).dataStrPtr).specificPtr).dettypeptr=*self.dettypeptr
 *(*(*(*davePtr).dataStrPtr).specificPtr).distptr=*self.distptr
 *(*(*(*davePtr).dataStrPtr).specificPtr).monptr=*self.monptr
 (*(*(*davePtr).dataStrPtr).specificPtr).run = self.run
 (*(*(*davePtr).dataStrPtr).specificPtr).time = self.time
 (*(*(*davePtr).dataStrPtr).specificPtr).user = self.user
 (*(*(*davePtr).dataStrPtr).specificPtr).comments= self.comments
 (*(*(*davePtr).dataStrPtr).specificPtr).chanwidth= self.chanwidth
 (*(*(*davePtr).dataStrPtr).specificPtr).chop_pulses= self.chop_pulses
 (*(*(*davePtr).dataStrPtr).specificPtr).chop_freq= self.chop_freq
 (*(*(*davePtr).dataStrPtr).specificPtr).monitor = self.monitor
 (*(*(*davePtr).dataStrPtr).specificPtr).date = self.date

  (*(*davePtr).dataStrPtr).commonStr.instrument = 'FCS'
  (*(*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
  *(*(*davePtr).dataStrPtr).commonStr.treatmentPtr = treatment
  (*(*davePtr).descriptr).name='Temperature'
  (*(*davePtr).descriptr).units='K'
  (*(*davePtr).descriptr).legend='Nominal sample temperature'
  (*(*davePtr).descriptr).qty=self.temperature
  (*(*davePtr).descriptr).err=0.0




  save, davePtr, filename = thisFile
endif

if n_phi_t ne 0 then begin
  data = *self.n_phi_t_ptr
  error = *self.ne_phi_t_ptr
  x = *self.timePtr
  y = *self.detPtr
  xlabel = 't'
  xunits = 'time:us'
  ylabel = '!4u!3'
  yunits = 'angle:degrees'
  histlabel = '!3I (!4u!3,t)'
  histunits = ''
  treatment =	['File:'+ self.filename]
  treatment = [[treatment],timestamp]
  if (*(*pState).prefsPtr).VAN_NORM then treatment = [[treatment],'Detectors Normalized to Vanadium Elastic Scattering']
  if (*(*pState).prefsPtr).SYMM_APP then treatment = [[treatment],'Symmetrized function']
  treatment = [[treatment],'--------------------------']
  if ptr_valid(self.treatmentPtr) then begin
      if n_elements(*self.treatmentPtr) ne 0 then begin
          treatment=[[treatment],*self.treatmentPtr]
      endif
  endif


  thisFile = fileout + '_n_phi_t.dave'

  *(*(*davePtr).dataStrPtr).commonStr.histPtr = $
				      		{qty:data,$
				       		 err:error,$
				      		 x:x,$
				      		 y:y}
(*(*(*davePtr).dataStrPtr).specificPtr).wavelength=self.wavelength
 *(*(*(*davePtr).dataStrPtr).specificPtr).detectorAnglesptr=*self.detPtr
 *(*(*(*davePtr).dataStrPtr).specificPtr).dettypeptr=*self.dettypeptr
 *(*(*(*davePtr).dataStrPtr).specificPtr).distptr=*self.distptr
 *(*(*(*davePtr).dataStrPtr).specificPtr).monptr=*self.monptr
 (*(*(*davePtr).dataStrPtr).specificPtr).run = self.run
 (*(*(*davePtr).dataStrPtr).specificPtr).time = self.time
 (*(*(*davePtr).dataStrPtr).specificPtr).user = self.user
 (*(*(*davePtr).dataStrPtr).specificPtr).comments= self.comments
 (*(*(*davePtr).dataStrPtr).specificPtr).chanwidth= self.chanwidth
 (*(*(*davePtr).dataStrPtr).specificPtr).chop_pulses= self.chop_pulses
 (*(*(*davePtr).dataStrPtr).specificPtr).chop_freq= self.chop_freq
 (*(*(*davePtr).dataStrPtr).specificPtr).monitor = self.monitor
 (*(*(*davePtr).dataStrPtr).specificPtr).date = self.date
  (*(*davePtr).dataStrPtr).commonStr.instrument = 'FCS'
  (*(*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
  *(*(*davePtr).dataStrPtr).commonStr.treatmentPtr = treatment
  (*(*davePtr).descriptr).name='Temperature'
  (*(*davePtr).descriptr).units='K'
  (*(*davePtr).descriptr).legend='Nominal sample temperature'
  (*(*davePtr).descriptr).qty=self.temperature
  (*(*davePtr).descriptr).err=0.0

  save, davePtr, filename = thisFile
endif

if n_phi_e ne 0 then begin
  data = *self.n_phi_e_ptr
  error = *self.ne_phi_e_ptr
  x = *self.energyPtr
  y = *self.detPtr
  xlabel = hbar_omega
  xunits = 'energy:meV'
  ylabel = '!4u!3'
  yunits = 'angle:degrees'
  histlabel = '!3I (!4u!3,E)'
  histunits = ''
  treatment =	['File:'+ self.filename, $
  				'converted from time to energy']
  treatment = [[treatment],timestamp]
  if (*(*pState).prefsPtr).VAN_NORM then treatment = [[treatment],'Detectors Normalized to Vanadium Elastic Scattering']
  if (*(*pState).prefsPtr).SYMM_APP then treatment = [[treatment],'Symmetrized function']
  if (*(*pState).prefsPtr).E_APP then treatment = [[treatment],'Energy Binned to ' +  string((*(*pState).prefsPtr).ELO) +',' $
		+  string((*(*pState).prefsPtr).EHI) +','+  string((*(*pState).prefsPtr).NEBINS)+' steps']
  treatment = [[treatment],'--------------------------']
  if ptr_valid(self.treatmentPtr) then begin
      if n_elements(*self.treatmentPtr) ne 0 then begin
          treatment=[[treatment],*self.treatmentPtr]
      endif
  endif


  thisFile = fileout + '_n_phi_e.dave'
  *(*(*davePtr).dataStrPtr).commonStr.histPtr = $
				      		{qty:data,$
				       		 err:error,$
				      		 x:x,$
				      		 y:y}
(*(*(*davePtr).dataStrPtr).specificPtr).wavelength=self.wavelength
 *(*(*(*davePtr).dataStrPtr).specificPtr).detectorAnglesptr=*self.detPtr
 *(*(*(*davePtr).dataStrPtr).specificPtr).dettypeptr=*self.dettypeptr
 *(*(*(*davePtr).dataStrPtr).specificPtr).distptr=*self.distptr
 *(*(*(*davePtr).dataStrPtr).specificPtr).monptr=*self.monptr
 (*(*(*davePtr).dataStrPtr).specificPtr).run = self.run
 (*(*(*davePtr).dataStrPtr).specificPtr).time = self.time
 (*(*(*davePtr).dataStrPtr).specificPtr).user = self.user
 (*(*(*davePtr).dataStrPtr).specificPtr).comments= self.comments
 (*(*(*davePtr).dataStrPtr).specificPtr).chanwidth= self.chanwidth
 (*(*(*davePtr).dataStrPtr).specificPtr).chop_pulses= self.chop_pulses
 (*(*(*davePtr).dataStrPtr).specificPtr).chop_freq= self.chop_freq
 (*(*(*davePtr).dataStrPtr).specificPtr).monitor = self.monitor
 (*(*(*davePtr).dataStrPtr).specificPtr).date = self.date

  (*(*davePtr).dataStrPtr).commonStr.instrument = 'FCS'
  (*(*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
  *(*(*davePtr).dataStrPtr).commonStr.treatmentPtr = treatment

  (*(*davePtr).descriptr).name='Temperature'
  (*(*davePtr).descriptr).units='K'
  (*(*davePtr).descriptr).legend='Nominal sample temperature'
  (*(*davePtr).descriptr).qty=self.temperature
  (*(*davePtr).descriptr).err=0.0
  save, davePtr, filename = thisFile
endif

if s_phi_e ne 0 then begin
  data = *self.s_phi_e_ptr
  error = *self.se_phi_e_ptr
  x = *self.energyPtr
  y = *self.detPtr
  xlabel = hbar_omega
  xunits = 'energy:meV'
  ylabel = '!4u!3'
  yunits = 'angle:degrees'
  histlabel = '!3S (!4u!3,E)'
  histunits = ''
  treatment =	['File:'+ self.filename, $
  				'converted from time to energy', $
  				'converted to scattering function']
  treatment = [[treatment],timestamp]
  if (*(*pState).prefsPtr).VAN_NORM then treatment = [[treatment],'Detectors Normalized to Vanadium ']
  if (*(*pState).prefsPtr).SYMM_APP then treatment = [[treatment],'Symmetrized function']
  if (*(*pState).prefsPtr).E_APP then treatment = [[treatment],'Energy Binned to ' +  string((*(*pState).prefsPtr).ELO) +',' $
		+  string((*(*pState).prefsPtr).EHI) +','+  string((*(*pState).prefsPtr).NEBINS)+' steps']
  treatment = [[treatment],'--------------------------']
  if ptr_valid(self.treatmentPtr) then begin
      if n_elements(*self.treatmentPtr) ne 0 then begin
          treatment=[[treatment],*self.treatmentPtr]
      endif
  endif

  thisFile = fileout + '_s_phi_e.dave'
  *(*(*davePtr).dataStrPtr).commonStr.histPtr = $
				      		{qty:data,$
				       		 err:error,$
				      		 x:x,$
				      		 y:y}
(*(*(*davePtr).dataStrPtr).specificPtr).wavelength=self.wavelength
 *(*(*(*davePtr).dataStrPtr).specificPtr).detectorAnglesptr=*self.detPtr
 *(*(*(*davePtr).dataStrPtr).specificPtr).dettypeptr=*self.dettypeptr
 *(*(*(*davePtr).dataStrPtr).specificPtr).distptr=*self.distptr
 *(*(*(*davePtr).dataStrPtr).specificPtr).monptr=*self.monptr
 (*(*(*davePtr).dataStrPtr).specificPtr).run = self.run
 (*(*(*davePtr).dataStrPtr).specificPtr).time = self.time
 (*(*(*davePtr).dataStrPtr).specificPtr).user = self.user
 (*(*(*davePtr).dataStrPtr).specificPtr).comments= self.comments
 (*(*(*davePtr).dataStrPtr).specificPtr).chanwidth= self.chanwidth
 (*(*(*davePtr).dataStrPtr).specificPtr).chop_pulses= self.chop_pulses
 (*(*(*davePtr).dataStrPtr).specificPtr).chop_freq= self.chop_freq
 (*(*(*davePtr).dataStrPtr).specificPtr).monitor = self.monitor
 (*(*(*davePtr).dataStrPtr).specificPtr).date = self.date

  (*(*davePtr).dataStrPtr).commonStr.instrument = 'FCS'
  (*(*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
  *(*(*davePtr).dataStrPtr).commonStr.treatmentPtr = treatment
  (*(*davePtr).descriptr).name='Temperature'
  (*(*davePtr).descriptr).units='K'
  (*(*davePtr).descriptr).legend='Nominal sample temperature'
  (*(*davePtr).descriptr).qty=self.temperature
  (*(*davePtr).descriptr).err=0.0
  save, davePtr, filename = thisFile
endif

if (s_q_e ne 0) and $
	(n_elements(*self.s_q_e_ptr) ne 0) then begin
  data = *self.s_q_e_ptr
  error = *self.se_q_e_ptr
  x = *self.energyPtr
  y = *self.qPtr
  xlabel = hbar_omega
  xunits = 'energy:meV'
  ylabel = '!3Q'
  yunits = 'wavevector:A-1'
  histlabel = '!3S (Q,E)'
  histunits = ''
    treatment =	['File:'+ self.filename, $
  				'converted from time to energy', $
  				'converted to scattering function', $
  				'rebinned to constant Q bins']
  treatment = [[treatment],timestamp]
  if (*(*pState).prefsPtr).VAN_NORM then treatment = [[treatment],'Detectors Normalized to Vanadium ']
  if (*(*pState).prefsPtr).SYMM_APP then treatment = [[treatment],'Symmetrized function']
  if (*(*pState).prefsPtr).E_APP then begin
                string1='Energy Binned to ' + string((*(*pState).prefsPtr).ELO) +',' $
		+ string((*(*pState).prefsPtr).EHI) +', with '+ string((*(*pState).prefsPtr).NEBINS)$
                + ' steps'
                treatment = [[treatment],string1]
  endif
  if (*(*pState).prefsPtr).Q_APP then treatment = [[treatment],'Q Binned to ' + string((*(*pState).prefsPtr).QLO) +',' $
		+ string((*(*pState).prefsPtr).QHI) +', with '+ string((*(*pState).prefsPtr).NQBINS)+' steps']
  treatment = [[treatment],'--------------------------']
  if ptr_valid(self.treatmentPtr) then begin
      if n_elements(*self.treatmentPtr) ne 0 then begin
          treatment=[[treatment],*self.treatmentPtr]
      endif
  endif

  thisFile = fileout + '_s_q_e.dave'
  *(*(*davePtr).dataStrPtr).commonStr.histPtr = $
				      		{qty:data,$
				       		 err:error,$
				      		 x:x,$
				      		 y:y}
(*(*(*davePtr).dataStrPtr).specificPtr).wavelength=self.wavelength
 *(*(*(*davePtr).dataStrPtr).specificPtr).detectorAnglesptr=*self.detPtr
 *(*(*(*davePtr).dataStrPtr).specificPtr).dettypeptr=*self.dettypeptr
 *(*(*(*davePtr).dataStrPtr).specificPtr).distptr=*self.distptr
 *(*(*(*davePtr).dataStrPtr).specificPtr).monptr=*self.monptr
 (*(*(*davePtr).dataStrPtr).specificPtr).run = self.run
 (*(*(*davePtr).dataStrPtr).specificPtr).time = self.time
 (*(*(*davePtr).dataStrPtr).specificPtr).user = self.user
 (*(*(*davePtr).dataStrPtr).specificPtr).comments= self.comments
 (*(*(*davePtr).dataStrPtr).specificPtr).chanwidth= self.chanwidth
 (*(*(*davePtr).dataStrPtr).specificPtr).chop_pulses= self.chop_pulses
 (*(*(*davePtr).dataStrPtr).specificPtr).chop_freq= self.chop_freq
 (*(*(*davePtr).dataStrPtr).specificPtr).monitor = self.monitor
 (*(*(*davePtr).dataStrPtr).specificPtr).date = self.date

  (*(*davePtr).dataStrPtr).commonStr.instrument = 'FCS'
  (*(*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
  *(*(*davePtr).dataStrPtr).commonStr.treatmentPtr = treatment
  (*(*davePtr).descriptr).name='Temperature'
  (*(*davePtr).descriptr).units='K'
  (*(*davePtr).descriptr).legend='Nominal sample temperature'
  (*(*davePtr).descriptr).qty=self.temperature
  (*(*davePtr).descriptr).err=0.0
  save, davePtr, filename = thisFile
endif

if (g_e ne 0) then begin
  data = *self.g_e_ptr
  error = *self.ge_e_ptr
  x = *self.energyPtr
  y = 0
  xlabel = hbar_omega
  xunits = 'energy:meV'
  ylabel = ''
  yunits = ''
  histlabel = '!3g(E)'
  histunits = ''
  treatment =	['File:'+ self.filename, $
  				'converted from time to energy', $
  				'converted to scattering function', $
  				'converted to density of states']
  treatment = [[treatment],timestamp]
  if (*(*pState).prefsPtr).VAN_NORM then treatment = [[treatment],'Detectors Normalized to Vanadium ']
  if (*(*pState).prefsPtr).SYMM_APP then treatment = [[treatment],'Symmetrized function']
  if (*(*pState).prefsPtr).E_APP then treatment = [[treatment],'Energy Binned to ' + string((*(*pState).prefsPtr).ELO) +',' $
		+ string((*(*pState).prefsPtr).EHI) +', with '+ string((*(*pState).prefsPtr).NEBINS) +' steps']
  if (*(*pState).prefsPtr).Q_APP then treatment = [[treatment],'Q Binned to ' + string((*(*pState).prefsPtr).QLO) +',' $
		+ string((*(*pState).prefsPtr).QHI) +', with '+ string((*(*pState).prefsPtr).NQBINS)+' steps']
  treatment = [[treatment],'--------------------------']
  if ptr_valid(self.treatmentPtr) then begin
      if n_elements(*self.treatmentPtr) ne 0 then begin
          treatment=[[treatment],*self.treatmentPtr]
      endif
  endif

  thisFile = fileout + '_g_e.dave'
  *(*(*davePtr).dataStrPtr).commonStr.histPtr = $
				      		{qty:data,$
				       		 err:error,$
				      		 x:x,$
				      		 y:y}
(*(*(*davePtr).dataStrPtr).specificPtr).wavelength=self.wavelength
 *(*(*(*davePtr).dataStrPtr).specificPtr).detectorAnglesptr=*self.detPtr
 *(*(*(*davePtr).dataStrPtr).specificPtr).dettypeptr=*self.dettypeptr
 *(*(*(*davePtr).dataStrPtr).specificPtr).distptr=*self.distptr
 *(*(*(*davePtr).dataStrPtr).specificPtr).monptr=*self.monptr
 (*(*(*davePtr).dataStrPtr).specificPtr).run = self.run
 (*(*(*davePtr).dataStrPtr).specificPtr).time = self.time
 (*(*(*davePtr).dataStrPtr).specificPtr).user = self.user
 (*(*(*davePtr).dataStrPtr).specificPtr).comments= self.comments
 (*(*(*davePtr).dataStrPtr).specificPtr).chanwidth= self.chanwidth
 (*(*(*davePtr).dataStrPtr).specificPtr).chop_pulses= self.chop_pulses
 (*(*(*davePtr).dataStrPtr).specificPtr).chop_freq= self.chop_freq
 (*(*(*davePtr).dataStrPtr).specificPtr).monitor = self.monitor
 (*(*(*davePtr).dataStrPtr).specificPtr).date = self.date

  (*(*davePtr).dataStrPtr).commonStr.instrument = 'FCS'
  (*(*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
  *(*(*davePtr).dataStrPtr).commonStr.treatmentPtr = treatment
  (*(*davePtr).descriptr).name='Temperature'
  (*(*davePtr).descriptr).units='K'
  (*(*davePtr).descriptr).legend='Nominal sample temperature'
  (*(*davePtr).descriptr).qty=self.temperature
  (*(*davePtr).descriptr).err=0.0
  save, davePtr, filename = thisFile
endif
if true then heap_free,davePtr
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::rebindata,group_leader = group_leader
if n_elements(group_leader) eq 0 then group_leader = 0L
if n_elements(*self.dataPtr) eq 0 then return
ndets = self.ndets
nchan = n_elements(*self.energyPtr)
datSize = size(*self.dataPtr)
val = 0

dat = (*self.dataPtr)[*,val]
daterr = (*self.errorPtr)[*,val]
x = (*self.energyPtr)[*,val]
y = (*self.detPtr)[0,val]

nx = n_elements(x)
result = opan_rebin_widget(x,group_leader = group_leader)
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
; 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

x_in = reform(*self.energyPtr)
zin = old
dzin = errold
drebin,*self.energyPtr,*self.dataPtr,*self.errorPtr,x_out,z_out,dz_out,$
	/points,/to_points,err=err,emsg=emsg

*self.dataPtr = z_out
*self.errorPtr = dz_out
*self.energyPtr = x_out

end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::energy_rebin,elo = elo, ehi = ehi, nbins = nbins
if n_elements(*self.dataPtr) eq 0 then return
ndets = self.ndets
nchan = n_elements(*self.energyPtr)
datSize = size(*self.dataPtr)
val = 0
dat = (*self.dataPtr)[*,val]
daterr = (*self.errorPtr)[*,val]
x = (*self.energyPtr)[*,val]
y = (*self.detPtr)[0,val]

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

ngrps = datSize[2]
old = *self.dataPtr
errold = *self.errorPtr

x_in = reform(*self.energyPtr)
zin = old
dzin = errold
if n_elements(*self.s_phi_e_ptr) gt 0 then begin
  drebin,x_in,*self.s_phi_e_ptr,*self.se_phi_e_ptr,x_out,z_out,dz_out,$
  	  /points,/to_points,err=err,emsg=emsg
  *self.energyPtr = x_out
  *self.s_phi_e_ptr = z_out
  *self.se_phi_e_ptr = dz_out
endif
if n_elements(*self.s_q_e_ptr) gt 0 then begin
  drebin,x_in,*self.s_q_e_ptr,*self.se_q_e_ptr,x_out,z_out,dz_out,$
  	  /points,/to_points,err=err,emsg=emsg
  *self.energyPtr = x_out
  *self.s_q_e_ptr = z_out
  *self.se_q_e_ptr = dz_out
endif
if n_elements(*self.n_phi_e_ptr) gt 0 then begin
  drebin,x_in,*self.n_phi_e_ptr,*self.ne_phi_e_ptr,x_out,z_out,dz_out,$
  	  /points,/to_points,err=err,emsg=emsg
  *self.energyPtr = x_out
  *self.n_phi_e_ptr = z_out
  *self.ne_phi_e_ptr = dz_out
endif
if n_elements(*self.g_e_ptr) gt 0 then begin
  drebin,x_in,*self.g_e_ptr,*self.ge_e_ptr,x_out,z_out,dz_out,$
  	  /points,/to_points,err=err,emsg=emsg
  *self.energyPtr = x_out
  *self.g_e_ptr = z_out
  *self.ge_e_ptr = dz_out
endif

self->calcQMatrix
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata::phi_to_q,nqbins = nqbins,qrange = qrange
; This uses rebinning to perform the conversion from constant scattering
; angle to constant wavevector transfer.
if n_elements(*self.energyPtr) eq 0 then return
lambda = self.wavelength

;include only active detectors!
incdet=*self.incdetptr
ny=n_elements(incdet)   ; how big is the new matrix?
y=fltarr(ny)
for i=0,ny-1 do begin      ; copy over to new array
  y[i]=(*self.detPtr)[incdet[i]]
endfor


;y = *self.detPtr	&	ny = n_elements(y)
x = *self.energyPtr	&	nx = n_elements(x)
uy = 1+bytarr(ny)
ux = 1+bytarr(nx)
xmat = x#uy
phimat = ux#y
eo = 81.81/lambda^2

alpha = 2.0717
qmat = sqrt((2.0*eo-xmat-2.0*sqrt(eo*(eo-xmat))*cos(phimat*!dtor))/alpha)
qmin = min(qmat,max = qmax)
if n_elements(qrange) eq 0 then qrange = [qmin,qmax]
qmin = qrange[0] & qmax = qrange[1]

if n_elements(nqbins) eq 0 then nqbins = fix(0.5*ny)
uq = 1+bytarr(nqbins)
dq = (qmax-qmin)/(nqbins-1.)
qbins = qmin+dq*findgen(nqbins)

dataReb = fltarr(nx,nqbins)
errReb = fltarr(nx,nqbins)

;do the same for the actual values....


zin = *self.s_phi_e_ptr
dzin = *self.se_phi_e_ptr


z=fltarr(n_elements(x),n_elements(y))
dz=fltarr(n_elements(x),n_elements(y))
for i=0,ny-1 do begin      ; copy over to new array
  z[*,i]=zin[*,incdet[i]]
  dz[*,i]=dzin[*,incdet[i]]
endfor



; Now do the rebinning in Q
for i = 0,nx-1 do begin
  z_in = reform(z[i,*])
  dz_in = reform(dz[i,*])
  x_in = reform(qmat[i,*])
  x_out = qbins
  drebin,x_in,z_in,dz_in,x_out,z_out,dz_out,/points,/to_points, $
  		err = err,emsg = emsg
  dataReb[i,*] = z_out
  errReb[i,*] = dz_out
endfor
*self.s_q_e_ptr = dataReb
*self.se_q_e_ptr = errReb
*self.qPtr = qbins
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function fcsdata::dataLoaded
return,self.dataLoaded
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function fcsdata::init,filename = filename,dataDirectory = dataDirectory
if n_elements(filename) eq 0 then filename = " "
if n_elements(dataDirectory) eq 0 then dataDirectory = $
	'c:\dimeo\idlprograms\FCS\data'
self.datDir = dataDirectory
self.filename = filename

self.dataLoaded = 0
self.colorTable = 5
self.dataPtr = ptr_new(/allocate_heap)
self.errorPtr = ptr_new(/allocate_heap)
self.rawdataPtr = ptr_new(/allocate_heap)
self.rawerrorPtr = ptr_new(/allocate_heap)
self.detPtr = ptr_new(/allocate_heap)
self.dettypePtr = ptr_new(/allocate_heap)
self.energyPtr = ptr_new(/allocate_heap)
self.timePtr = ptr_new(/allocate_heap)
self.distPtr = ptr_new(/allocate_heap)
self.qPtr = ptr_new(/allocate_heap)
self.monPtr = ptr_new(/allocate_heap)
self.vanefficiencyPtr = ptr_new(/allocate_heap)
; CONTAINERS FOR AVERAGED DATA
self.GPdetPtr = ptr_new(/allocate_heap)
self.GPqPtr = ptr_new(/allocate_heap)

self.qmatPtr = ptr_new(/allocate_heap)
self.n_phi_t_ptr=ptr_new(/allocate_heap)
self.ne_phi_t_ptr=ptr_new(/allocate_heap)
self.n_phi_e_ptr=ptr_new(/allocate_heap)
self.ne_phi_e_ptr=ptr_new(/allocate_heap)
self.s_phi_e_ptr=ptr_new(/allocate_heap)
self.se_phi_e_ptr=ptr_new(/allocate_heap)
self.s_q_e_ptr=ptr_new(/allocate_heap)
self.se_q_e_ptr=ptr_new(/allocate_heap)
self.g_e_ptr=ptr_new(/allocate_heap)
self.ge_e_ptr=ptr_new(/allocate_heap)
self.treatmentptr=ptr_new(/allocate_heap)
self.incdetptr=ptr_new(/allocate_heap)
self.grpdetptr=ptr_new(/allocate_heap)
self.temperature = 0D
self.nchan = 0L
self.ndets = 0L
self.run = 0L
self.time = 0D
self.user = ""
self.comments = ""
self.chanwidth = 0D
self.chop_pulses = 0D
self.chop_freq = 0D
self.wavelength = 0D
self.monitor = 0D
if self.filename ne " " then begin
  self->loadData,filename = self.filename
endif
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcsdata__define
; Definition module for FCSDATA class
define = {fcsdata,$
	  treatmentptr:ptr_new(), $
          datDir:"", $
          dataLoaded:0,	$
          temperature:0D, $
          run:0L, $
          time:0D, $
          user:"", $
          comments:"", $
          chanwidth:0D, $
          chop_pulses:0D, $
          chop_freq:0D, $
          wavelength:0D, $
          monitor:0D, $
          date:"", $
          ndets:0L, $
          filename:"", $
          nchan:0L, $
          monPtr:ptr_new(), $
          qPtr:ptr_new(), $
          distPtr:ptr_new(), $
          detPtr:ptr_new(), $
          dettypePtr:ptr_new(), $
          incdetptr:ptr_new(),$
          grpdetptr:ptr_new(),$
	  timePtr:ptr_new(), $
          energyPtr:ptr_new(), $
;containers for grouped data
          GPqPtr:ptr_new(), $
          GPdetPtr:ptr_new(), $
;
          qmatPtr:ptr_new(),	$
          n_phi_t_ptr:ptr_new(),	$
          ne_phi_t_ptr:ptr_new(),	$
          n_phi_e_ptr:ptr_new(),	$
          ne_phi_e_ptr:ptr_new(),	$
          s_phi_e_ptr:ptr_new(),	$
          se_phi_e_ptr:ptr_new(),	$
          s_q_e_ptr:ptr_new(),		$
          se_q_e_ptr:ptr_new(),		$
          g_e_ptr:ptr_new(),		$
          ge_e_ptr:ptr_new(),		$
          vanefficiencyPtr:ptr_new(),	$		; only useful for vanadium...
	  colorTable:0,				$
          dataPtr:ptr_new(), $
          errorPtr:ptr_new(), $
          rawdataPtr:ptr_new(), $
          rawerrorPtr:ptr_new() $
          }

return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro fcs_example
filename = "c:\bhd\dave_development\dave\programs\testing\FCS\data\bulk.fcs"
o = obj_new('fcsdata',filename = filename)
;o->loaddata,$
;   filename = "c:\bhd\dave_development\dave\programs\testing\FCS\data\bulk.fcs"
dummy = o->timetoenergy()
window,1,xsize = 500,ysize = 500
o->draw,detector = 45;drawimage
obj_destroy,o
end
