; $Id: dtas_reduce__define.pro,v 1.31 2008/05/23 19:35:05 azuah Exp $
;###############################################################################
;
; NAME:
;  DTAS_REDUCE__DEFINE
;
; PURPOSE:
;  Class for TAS reduction.
;
; CATEGORY:
;  DAVE, TAS, 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.
;
;###############################################################################
; +
;
; -
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro dtas_reduce::cleanup
nobj = self.intermediate->count()
if nobj gt 0 then begin
   oall = self.intermediate->get(/all)
   obj_destroy,oall
endif
nsignal = self.signal_container->count()
nempty = self.empty_container->count()
nreduced = self.reduced_container->count()
nint = self.intermediate->count()
nzip = self.zip_container->count()
ncombine = self.combine_container->count()
if ncombine gt 0 then begin
   oall = self.combine_container->get(/all)
   obj_destroy,oall
endif
if nint gt 0 then begin
   oall = self.intermediate->get(/all)
   obj_destroy,oall
endif
if nzip gt 0 then begin
   oall = self.zip_container->get(/all)
   obj_destroy,oall
endif
if nsignal gt 0 then begin
   oall = self.signal_container->get(/all)
   obj_destroy,oall
endif
if nempty gt 0 then begin
   oall = self.empty_container->get(/all)
   obj_destroy,oall
endif
if nreduced gt 0 then begin
   oall = self.reduced_container->get(/all)
   obj_destroy,oall
endif
n_here = self->count()
if n_here gt 0 then begin
   oall = self->get(/all)
   obj_destroy,oall
endif
obj_destroy,[self.intermediate,self.signal_container,self.empty_container, $
   self.reduced_container,self.zip_container,self.combine_container]
self->IDL_CONTAINER::cleanup
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function dtas_reduce::destroy_element,o
if not obj_valid(o) then return,0
ret = o->get_property(signal_type = signal_type)
case signal_type of
'SIGNAL': $
   begin
      self.signal_container->remove,o
      obj_destroy,o
   end
'INTERMEDIATE': $
   begin
      self.intermediate->remove,o
      obj_destroy,o
   end
'EMPTY': $
   begin
      self.empty_container->remove,o
      obj_destroy,o
   end
'REDUCED': $
   begin
      self.reduced_container->remove,o
      obj_destroy,o
   end
'COMBINED': $
   begin
      self.combine_container->remove,o
      obj_destroy,o
   end
else:
endcase
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function dtas_reduce::update_plot_structure,oadd,DISPLAY = display,_Extra = extra
if n_params() eq 0 then return,0
if not obj_valid(oadd) then return,0
; Updates the plottable data in the object
; Get out the plot structure
ret = oadd->get_property(plot_str = plot_str,data = data,var_names = var_names)
xvar_name = plot_str.xtitle
x_index = where(var_names eq xvar_name,count_x)
if count_x gt 0 then begin
   *plot_str.x = *data.(x_index[0])
endif
nvars = n_elements(var_names)
*plot_str.y = *data.(nvars-2)
*plot_str.dy = *data.(nvars-1)
ret = oadd->set_property(plot_str = plot_str)
if n_elements(display) ne 0 then $
   ret = oadd->display(psym = 4,_Extra = extra)
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function dtas_reduce::sum_poisson,msg = msg, signal = signal, empty = empty, o = o
if keyword_set(empty) then ocont = self.empty_container
if keyword_set(signal) then ocont = self.signal_container
msg = ''
n_obj = ocont->count()
if n_obj lt 1 then begin
   msg = 'At least two files are required in the summation'
   return,0
endif
if n_obj eq 1 then begin
   msg = 'Only one data file present'
   oall = ocont->get(/all)
   o = oall[0]->clone()
   oadd = oall[0]->clone()
   return,oadd
endif
; Check if the independent variables match up
oall = ocont->get(/all)
for i = 0,n_obj-1 do begin
   ret = oall[i]->get_property(plot_str = plot_str)
   if i eq 0 then xtitle = plot_str.xtitle else $
      xtitle = [xtitle,plot_str.xtitle]
endfor
if n_elements(uniq(xtitle)) ne 1 then begin
   msg = 'The independent variables do not match up!'
   return,0
endif
ret = oall[0]->get_property(data = data,total_mon = m0,  $
      plot_str = plot_str, var_names = var_names)
if m0 eq 0 then begin
   msg = 'No monitor counts have been defined for this file'
   return,0
endif
baseline_mon = m0
nvars = n_elements(var_names)
var_ptr = ptrarr(nvars,/allocate_heap)

for i = 0,n_elements(oall)-1 do begin
   ret = oall[i]->get_property(data = data,total_mon = m0,  $
         plot_str = plot_str, var_names = var_names,filename = filename)
   if i eq 0 then files_to_sum = filename else $
      files_to_sum = [files_to_sum,filename]
   if i eq 0 then begin
      for j = 0,nvars-1 do begin
         *var_ptr[j] = *data.(j)
      endfor

   endif else begin
      sf = double(baseline_mon)/double(m0)
      for j = 0,nvars-1 do begin
         if (j eq (nvars-1)) or (j eq (nvars-2)) then begin
            ; Must scale the intensity and error bars by the monitor
            ; scale factor
            *var_ptr[j] = [*var_ptr[j],sf*(*data.(j))]
         endif else begin
            *var_ptr[j] = [*var_ptr[j],*data.(j)]
         endelse
      endfor

   endelse
endfor

; Create a data structure consistent with the one in DTAS_DATA class
oadd = oall[0]->clone()
ret = oadd->get_property(data = data,time = time)
for i = 0,nvars-1 do begin
   *data.(i) = *var_ptr[i]
endfor
; Get the independent variable out
xvar_name = plot_str.xtitle
x_index = where(strupcase(var_names) eq strupcase(xvar_name),count_index)
if count_index eq 0 then return,0

ret = oadd->set_property(data = data)
ptr_free,var_ptr

y = reform(*data.(nvars-2))
dy = reform(*data.(nvars-1))

x = reform(*data.(x_index[0]))
; Convert data structure to a data array
darray = oadd->convert_data_to_array()
merge_str = merge_arrays(x,y,dy,uniq_index = uniq_index)
nx = n_elements(merge_str.x)
new_array = dblarr(nx,nvars)
; Fill in what we know first
new_array[*,x_index[0]] = merge_str.x
new_array[*,nvars-2] = merge_str.y
new_array[*,nvars-1] = merge_str.dy

; Now must merge the rest of the data array!
new_array[*,0] = dindgen(nx)
for i = 1,nvars-3 do begin
   if i ne x_index[0] then begin
      new_array[*,i] = darray[uniq_index,i]
   endif
endfor

; Now convert the updated data array to a structure and put it back
; in to the summed object
ret = oadd->convert_array_to_data(new_array)
ret = oadd->get_property(data = new_data)


; Add the summed data to the TAS data container and the reduced container
self->add,oadd
self.reduced_container->add,oadd
ret = self->update_plot_structure(oadd)
nfiles = n_elements(files_to_sum)
for i = 0,nfiles-1 do begin
   new_treatment = 'File added (Poisson): '+files_to_sum[i]
   ret = oadd->update_treatment(new_treatment)
endfor
ret = oadd->update_display_name('_add')
if arg_present(o) then o = oadd->clone()
return,oadd
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function dtas_reduce::check_same_tags,oarray
; This function checks to see if the tag names in the data structure are
; the same
nobj = n_elements(oarray)
if nobj eq 1 then return,1B
ret = oarray[0]->get_property(var_names = names)
n = n_elements(names)

for i = 1,nobj-1 do begin
   ret = oarray[i]->get_property(var_names = names)
   n1 = n_elements(names)
   if n1 ne n then return,0B
   n = n1
endfor

return,1B
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function dtas_reduce::remove_mismatch,oarray
; This function removes a single column from the data that cannot be found in
; the other data sets.
nobj = n_elements(oarray)
if nobj eq 1 then return,1B
; First go through each object and get out the variable names
index = -1
for i = 0,nobj-1 do begin
   ret = oarray[i]->get_property(var_names = names)
   if i eq 0 then all_names = names else $
      all_names = [all_names,names]
endfor

; What are the unique names?
uniq_names = all_names[uniq(all_names,sort(all_names))]

; Now find out which objects lack one of the elements in the unique
; name list
for i = 0,nobj-1 do begin
   ret = oarray[i]->get_property(var_names = names,data = this_data)
   help,*this_data.(0)
   nunique_names = n_elements(uniq_names)
   for j = 0,nunique_names-1 do begin
      this_name = uniq_names[j]
      match_index = where(names eq this_name,count_match)
      if count_match eq 0 then bad_var_name = this_name
   endfor
endfor
; Now go through each object and eliminate the data member whose
; tag matches bad_var_name
for i = 0,nobj-1 do begin
   ret = oarray[i]->get_property(data = data,var_names = var_names)
   bad_index = where(var_names eq bad_var_name,count_match)
   if count_match gt 0 then begin
      ; Create a new data structure without the offensive data member
      nvars = n_elements(var_names)
      var_ptr = ptrarr(nvars,/allocate_heap)

      s = create_struct(var_names[0],data.(0))
      for j = 1,nvars-1 do begin
         if j ne bad_index then begin
            s = create_struct(s,var_names[j],data.(j))
         endif
      endfor
      ret = oarray[i]->set_property(data = s)
   endif
endfor
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function dtas_reduce::sum_straight,msg = msg, signal = signal, empty = empty, $
   o = o
; This function performs a straight addition of raw counts and monitor counts
; for overlapping points and a simple merge for non-overlapping points (but
; the monitor values are updated).
if keyword_set(empty) then ocont = self.empty_container
if keyword_set(signal) then ocont = self.signal_container
msg = ''
n_obj = ocont->count()
if n_obj lt 1 then begin
   msg = 'At least two files are required in the summation'
   return,0
endif
if n_obj eq 1 then begin
   msg = 'Only one data file present'
   oall = ocont->get(/all)
   oadd = oall[0]->clone()
   o = oadd->clone()
   self->add,oadd
   return,oadd
endif
; Check if the independent variables match up
oall = ocont->get(/all)

; The following statement checks to see if there are identical
; tag names in the independent variables.
ret = self->check_same_tags(oall)
if ~ret then ret = self->remove_mismatch(oall)

for i = 0,n_obj-1 do begin
   ret = oall[i]->get_property(plot_str = plot_str)
   if i eq 0 then xtitle = plot_str.xtitle else $
      xtitle = [xtitle,plot_str.xtitle]
endfor
if n_elements(uniq(xtitle)) ne 1 then begin
   msg = 'The independent variables do not match up!'
   return,0
endif
ret = oall[0]->get_property(data = data,total_mon = m0,  $
      plot_str = plot_str, var_names = var_names)

if m0 eq 0 then begin
   msg = 'No monitor counts have been defined for this file'
   return,0
endif
baseline_mon = m0
nvars = n_elements(var_names)
var_ptr = ptrarr(nvars,/allocate_heap)

for i = 0,n_elements(oall)-1 do begin
   ret = oall[i]->get_property(data = data,total_mon = m0,  $
         plot_str = plot_str, var_names = var_names,        $
         treatment = treatment, filename = filename)
  ; for j = 0,n_elements(var_names)-1 do help,*data.(j)

;   if i gt 0 then begin
;      if i eq 1 then treatment_array = treatment else $
;         treatment_array = [treatment_array,treatment]
;   endif
   treatment_array = (i eq 0)? treatment : [treatment_array,treatment]

   if i eq 0 then begin
      files_to_sum = filename
      tags = tag_names(data)
      for j = 0,nvars-1 do begin
         *var_ptr[j] = *data.(j)
      endfor

   endif else begin
      files_to_sum = [files_to_sum,filename]
      sf = double(baseline_mon)/double(m0)
      tags = tag_names(data)
      mon_index = where(strupcase(tags) eq 'MON',count_mon)
      if (count_mon eq 0) then $ ; then try looking for MONITOR tag
        mon_index = where(strupcase(tags) eq 'MONITOR',count_mon)

      for j = 0,nvars-1 do begin
         if (j eq (nvars-1)) or (j eq (nvars-2)) or $
            (j eq mon_index) then begin
            *var_ptr[j] = [*var_ptr[j],(*data.(j))]
         endif else begin
            *var_ptr[j] = [*var_ptr[j],*data.(j)]
         endelse
      endfor

   endelse
endfor
; Create a data structure consistent with the one in DTAS_DATA class
oadd = oall[0]->clone()
ret = oadd->update_treatment(treatment_array)
ret = oadd->get_property(data = data,time = time)

for i = 0,nvars-1 do begin
   *data.(i) = *var_ptr[i]
endfor
; Get the independent variable out
xvar_name = plot_str.xtitle
x_index = where(strupcase(var_names) eq strupcase(xvar_name),count_index)
if count_index eq 0 then return,0

ret = oadd->set_property(data = data)
ptr_free,var_ptr
tags = strupcase(tag_names(data))
mon_index = where(tags eq 'MON',count_mon)
if (count_mon eq 0) then $      ; then try looking for MONITOR tag
  mon_index = where(tags eq 'MONITOR',count_mon)
  
time_index = where(tags eq 'MIN',count_time)
if (count_time eq 0) then $ ; then look for TIME tag 
  time_index = where(tags eq 'TIME',count_time)

mon = reform(*data.(mon_index)) & dmon = sqrt(mon)
time = reform(*data.(time_index)) & dtime = sqrt(time)
y = reform(*data.(nvars-2))
dy = reform(*data.(nvars-1))

x = reform(*data.(x_index[0]))
; Convert data structure to a data array
darray = oadd->convert_data_to_array()
; Now be very careful here.  We need the contents of x in the original order
; in which they were stored in the data pointer.  The function called
; MERGE_ARRAYS will change the values.  Must make a clean copy here and use
; it in the function call.
xclean = x
xin = xclean
merge_str = merge_arrays(xin,y,dy,uniq_index = uniq_index,/straight)
xin = xclean
mon_merge = merge_arrays(xin,mon,dmon,uniq_index = mon_uniq,/straight)
xin = xclean
time_merge = merge_arrays(xin,time,dtime,uniq_index = time_uniq,/straight)

nx = n_elements(merge_str.x)
new_array = dblarr(nx,nvars)

; Now must insert the rest of the data array!
new_array[*,0] = dindgen(nx)

;for i = 2,nvars-3 do begin
for i = 1,nvars-3 do begin      ; RTA - don't understand why starting from 2, change to 1
;   if ((i ne x_index[0])) then begin
      new_array[*,i] = darray[uniq_index,i]
;   endif
endfor

; Fill in what we know *after* - RTA
new_array[*,x_index[0]] = merge_str.x
new_array[*,nvars-2] = merge_str.y
new_array[*,nvars-1] = merge_str.dy
new_array[*,mon_index[0]] = mon_merge.y
new_array[*,time_index[0]] = time_merge.y


; Now convert the updated data array to a structure and put it back
; in to the summed object
ret = oadd->convert_array_to_data(new_array)
ret = oadd->get_property(data = new_data)


; Add the summed data to the TAS data container and the reduced container
self->add,oadd
self.reduced_container->add,oadd
ret = self->update_plot_structure(oadd)
nfiles = n_elements(files_to_sum)
if keyword_set(empty) then str = '(empty)'
if keyword_set(signal) then str = '(signal)'
for i = 0,nfiles-1 do begin
   new_treatment = str+'File added (straight): '+files_to_sum[i]
   ret = oadd->update_treatment(new_treatment)
endfor
ret = oadd->update_display_name('_add')
if arg_present(o) then o = oadd->clone()
return,oadd
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function dtas_reduce::subtract_data,oadd,osubtract,msg = msg,o = o
; RTA - Make major changes. The basis of the subtracted results was
;       osubtraction, change this to oadd! It was assummed that the
;       monitor is always at column 2 in the data structure which is
;       not the case.

if n_params() ne 2 then begin
   msg = 'Two data objects must be passed in'
   return,0
endif
msg = ''

; The following statement checks to see if there are identical
; tag names in the independent variables.
ret = self->check_same_tags([oadd,osubtract])
if ~ret then ret = self->remove_mismatch([oadd,osubtract])

; First see if the independent variables match
ret = oadd->get_property(data = add_data,total_mon = m0,  $
   plot_str = plot_str, var_names = var_names)
nvars = n_elements(var_names)

if m0 eq 0 then begin
   msg = 'No monitor counts have been defined for this file'
   return,0
endif

; Get the INTENSITY and ERROR values out
y_baseline = reform(*add_data.(nvars-2))
dy_baseline = reform(*add_data.(nvars-1))
; Get the independent variable out
xvar_name = plot_str.xtitle
x_index = where(strupcase(var_names) eq strupcase(xvar_name),count_index)
if count_index eq 0 then return,0
x_index = x_index[0]

; ************************************************************
; Sort the "add" data out by x_index (added by RMD 05/18/04)
xsort = sort(*add_data.(x_index))
for i = 0,n_tags(add_data)-1 do $
   (*add_data.(i))[*] = (*add_data.(i))[xsort]
y_baseline = y_baseline[xsort]
dy_baseline = dy_baseline[xsort]
; ************************************************************

x_baseline = reform(*add_data.(x_index))

; Make a copy of the first TAS data object for use in the difference.

;-------------START CHANGE - RTA
; ** replace following lines
;mon_baseline = reform(*add_data.(1))
;oresult = osubtract->clone()
;ret = oadd->get_property(treatment = add_treatment,descriptr_str = descriptr_str)
;ret = oresult->set_private(treatment = add_treatment)
;ret = oresult->set_property(descriptr_str = descriptr_str)
;ret = osubtract->get_property(plot_str = plot_str,total_mon = mi, $
;      filename = filename,treatment = treatment,var_names = var_names, $
;      data = sub_data)
;ret = oresult->update_treatment(treatment)
;ret = oresult->get_property(data = result_data)
; * with
mon_baseline = oadd->get_monitor()
oresult = oadd->clone()
ret = osubtract->get_property(total_mon = mi,filename=sub_filename,data = sub_data)
ret = oresult->get_property(data = result_data)
;--------------END CHANGE - RTA


; ************************************************************
; Sort the data on the indexed variable (added by RMD 05/18/04)
x1sort = sort(*result_data.(x_index))
x2sort = sort(*sub_data.(x_index))
for i = 0,n_tags(result_data)-1 do begin
   (*result_data.(i))[*] = (*result_data.(i))[x1sort]
   (*sub_data.(i))[*] = (*sub_data.(i))[x2sort]
endfor
; ************************************************************

if mi eq 0 then begin
      obj_destroy,oresult
      msg = 'No monitor counts have been defined for this file'
      return,0
endif


xi = reform(*sub_data.(x_index))
; Now get the minimal range based on the intersection of x_baseline and xi
xb_min = min(x_baseline,max = xb_max)
xi_min = min(xi,max = xi_max)
xmin = max([xb_min,xi_min])
xmax = min([xb_max,xi_max])
where_xn_baseline = where((x_baseline ge xmin) and (x_baseline le xmax),count)
if count lt 2 then begin
   obj_destroy,oresult
   msg = 'There is insufficient overlap of the files used in the subtraction'
   return,0
endif
x_baseline = x_baseline[where_xn_baseline]

; *************************BUG FIX*************************************
; The following lines were added 2-10-04 by RMD.  These lines fix a bug
; in the original subtraction algorithm where data ranges were not
; adequately accounted for in the subtraction.
y_baseline = y_baseline[where_xn_baseline]
dy_baseline = dy_baseline[where_xn_baseline]
; ***********************END BUG FIX***********************************

xb_min = min(x_baseline,max = xb_max)

if (total(xi-x_baseline) ne 0.0) or $
   (n_elements(xi) ne n_elements(x_baseline)) then begin
; Interpolate the values contained in osubtract to those of the independent
; variable in oadd:
   *result_data.(x_index) = x_baseline
   ntags = n_tags(sub_data)
   for i = 0,ntags-1 do begin
      if i ne x_index[0] then begin
         v = double(reform(*sub_data.(i)))
         x = xi
         u = x_baseline
         ok = where((x ge xb_min) and (x le xb_max))
         v = v[ok]
         x = x[ok]
         case self.interp of
         'BINNING':  $
            begin
               dv = sqrt(v)
               drebin,x,v,dv,u,vi,dvi,/points,/to_points,err = err,emsg = emsg
               if i eq (ntags-2) then dvi_out = dvi
            end
         'LINEAR':   $
            begin
               vi = interpol(v,x,u)
            end
         'SPLINE':   $
            begin
               vi = interpol(v,x,u,/spline)
            end
         endcase
         if ( (i eq (ntags-1)) and (self.interp eq 'BINNING') ) then begin
            *result_data.(i) = dvi_out
         endif else begin
            *result_data.(i) = vi
         endelse
      endif
   endfor
endif

;-------------START CHANGE - RTA
; ** replace following lines
;moni = reform(*result_data.(1))
;; Update the monitor (RMD 03-17-04)
;mon_new = moni*weight
;ret = oresult->set_monitor(mon_new)
;yi = reform(*result_data.(nvars-2))
;dyi = reform(*result_data.(nvars-1))
; * with
moni = osubtract->get_monitor()
yi = reform(*sub_data.(nvars-2))
dyi = reform(*sub_data.(nvars-1))
;--------------END CHANGE - RTA


mon_baseline = mon_baseline[where_xn_baseline]
weight = double(mon_baseline)/double(moni)
y_baseline = y_baseline - weight*yi
dy_baseline = sqrt(dy_baseline^2+(weight*dyi)^2)


;new_treatment = strarr(2)
new_treatment = 'Empty subtracted from signal'
new_treatment = [new_treatment,'Empty file: '+sub_filename]
if self.interp eq 'BINNING' then $
   new_treatment = [new_treatment,self.interp+' used in subtraction'] else $
      new_treatment = [new_treatment,self.interp+' interpolation used in subtraction']

ret = oresult->update_treatment(new_treatment)

*result_data.(nvars-2) = y_baseline
*result_data.(nvars-1) = dy_baseline
; Update the DATA for the object
ret = oresult->set_property(data = result_data)
; Add the summed data to the TAS data container and the reduced container
self->add,oresult
self.reduced_container->add,oresult
; Update the plot structure
ret = self->update_plot_structure(oresult)
;ret = oresult->update_display_name('_subtracted')


; Plot the data
if arg_present(o) then o = oresult->clone()
;-------------START CHANGE - RTA
; Copy the name of the "ADD" file to this one:
;ret = oadd->get_property(display_name = display_name, $
;   filename = filename)
;ret = oresult->set_property(display_name = display_name, $
;   filename = filename)
;ret = oresult->get_property(descriptr_str = d)
;--------------END CHANGE - RTA
ret = oresult->update_display_name('_subtracted')
return,oresult
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function dtas_reduce::set_property, INTERP = interp,  $
                                    XVAR_NAME = xvar_name
if n_elements(interp) ne 0 then self.interp = interp
if n_elements(xvar_name) ne 0 then begin
   nobj = self->count()
   if nobj ne 0 then begin
      oall = self->get(/all)
      for i = 0,nobj-1 do begin
         ret = oall[i]->set_property(xvar_name = xvar_name)
      endfor
   endif
endif
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function dtas_reduce::get_property, INTERMEDIATE = intermediate,  $
                                    REDUCED = reduced,            $
                                    SIGNAL = signal,              $
                                    INTERP = interp,              $
                                    COMBINED = combined,          $
                                    EMPTY = empty
combined = self.combine_container
interp = self.interp
reduced = self.reduced_container
signal = self.signal_container
empty = self.empty_container
intermediate = self.intermediate
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function dtas_reduce::clear_zip
; This function clears the objects from the zip container and destroys them.
nzip = self.zip_container->count()
if nzip gt 0 then begin
   oall = self.zip_container->get(/all)
   obj_destroy,oall
   obj_destroy,self.zip_container
   self.zip_container = obj_new('IDL_CONTAINER')
endif
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function dtas_reduce::add_zip,oarray
; This function adds the data objects contained in OARRAY to the zip container
; object, sorting them by their descriptive variable first.  Clone the objects
; before adding!
n = n_elements(oarray)
for i = 0,n-1 do begin
   if obj_valid(oarray[i]) then begin
      clone = oarray[i]->clone()
      self.zip_container->add,clone
   endif
endfor
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function dtas_reduce::check_indep_var_limits
; Check the limits of the independent variable.
;n = self.zip_container->count()
;if n eq 0 then return,0
;oall = self.zip_container->get(/all)
;for i = 0,n-1 do begin
;   ret = oall[i]->get_property(plot_str = plot_str)
;   x = *plot_str.x
;   if i eq 0 then begin
;
;   endif else begin
;   endelse
;endfor
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function dtas_reduce::zip_data, MSG = msg,DISPLAY = display
; This function "ZIPS" data together based on a descriptive variable such as
; temperature or magnetic field for instance.  The data used in the zipping is
; selected from any of the SIGNAL, EMPTY, or REDUCED containers.
msg = ''
nzip = self.zip_container->count()
if nzip eq 0 then return,0
; Search for the descriptive variable for each and determine if they are the
; same.
oall = self.zip_container->get(/all)

dflag = 1
for i = 0,nzip-1 do begin
   ret = oall[i]->get_property(descriptr_str = descriptr)
   if n_elements(descriptr) eq 0 then dflag = 0
endfor

for i = 0,nzip-1 do begin
   ret = oall[i]->get_property(descriptr_str = descriptr_str,plot_str = plot_str)
   ; If there is no descriptr present then fill it out
   if dflag eq 0 then begin
      ret = oall[i]->set_index_descriptr(float(i))
      ret = oall[i]->get_property(descriptr_str = descriptr_str)
   endif
   ;ret = oall[i]->get_property(descriptr_str = descriptr_str)
   if i eq 0 then begin
      name = descriptr_str.name
      qty = descriptr_str.qty
      xtitles = plot_str.xtitle
   endif else begin
      name = [name,descriptr_str.name]
      qty = [qty,descriptr_str.qty]
      xtitles = [xtitles,plot_str.xtitle]
   endelse
endfor
if n_elements(uniq(name)) gt 1 then begin
   msg = 'Combining variables do not match!'
   return,0
endif
;  See if the first independent variable, in particular the name of x,
;  is the same for all data sets
if n_elements(uniq(xtitles)) gt 1 then begin
   msg = 'Independent variable types do not match!'
   return,0
endif
dsort = sort(qty)
; Perform a check on the independent variables
ret = self->check_indep_var_limits()
if ret eq 0 then return,0
; Reorganize the zip container based on the sorting
for i = 0,nzip-1 do self.zip_container->remove,oall[i]
obj_destroy,self.zip_container
self.zip_container = obj_new('IDL_CONTAINER')
for i = 0,nzip-1 do self.zip_container->add,oall[dsort[i]]
oall = self.zip_container->get(/all)
nzip = self.zip_container->count()
; Now put the final data set together
for i = 0,nzip-1 do begin
   ret = oall[i]->get_property(descriptr_str = descriptr_str, $
      plot_str = plot_str,treatment = this_treatment)
   if i eq 0 then treatment = this_treatment else $
      treatment = [treatment,this_treatment]
   ny = nzip
   if i eq 0 then begin
      nx = n_elements(*plot_str.x)
      x = *plot_str.x
      z = dblarr(nx,ny)
      dz = dblarr(nx,ny)
      y = dblarr(ny)
   endif
   xnew = *plot_str.x
   y[i] = descriptr_str.qty
   if total(xnew-x) ne 0 then begin ; interpolate
      znew = *plot_str.y
      dznew = *plot_str.dy
      nx = n_elements(x)
      nxnew = n_elements(xnew)
      if (x[0] lt xnew[0] or x[nx-1] gt xnew[nxnew-1]) then begin
          ;; avoid any extrapolation by running interpol() only for
          ;; values of x that lie within xnew 
          inIndex = where(x ge xnew[0] and x le xnew[nxnew-1], inCnt $
                        ,complement=outIndex,ncomplement=outCnt)
          if (inCnt gt 0) then begin
              xValid = x[inIndex]
              zinterp = interpol(znew,xnew,xValid)
              dzinterp = interpol(dznew,xnew,xValid)
              z[inIndex,i] = zinterp
              dz[inIndex,i] = dzinterp
          endif
          if (outCnt gt 0) then begin
              z[outIndex,i] = 0.0
              dz[outIndex,i] = 0.0
          endif
      endif
   endif else begin
      z[*,i] = *plot_str.y
      dz[*,i] = *plot_str.dy
      x = *plot_str.x
   endelse
endfor
treatment = [treatment,'The files listed above were combined']
; The return value will be a structure
px = ptr_new(y,/no_copy)
py = ptr_new(x,/no_copy)
pz = ptr_new(transpose(z),/no_copy)
bad_index = where(dz eq 0.0,count_bad)
if count_bad gt 0 then dz[bad_index] = 1.0
pdz = ptr_new(transpose(dz),/no_copy)
xtitle = descriptr_str.legend
ytitle = plot_str.xtitle

str = {  px:px,               $
         py:py,               $
         pz:pz,               $
         pdz:pdz,             $
         xtitle:xtitle,       $
         ytitle:ytitle        }

o = obj_new('dtas_zip')
ret = o->set_property(  xval = *str.px,                  $
                        yval = *str.py,                  $
                        z = *str.pz,                     $
                        dz = *str.pdz,                   $
                        xtitle = xtitle,                 $
                        display_name = 'COMBINED_DATA',  $
                        treatment = treatment,           $
                        ytitle = ytitle                  )
heap_free,str
self.combine_container->add,o
return,o
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function dtas_reduce::zip_non_orth_data, MSG = msg,DISPLAY = display
; This function "ZIPS" data together based on a descriptive variable such as
; temperature or magnetic field for instance.  The data used in the zipping is
; selected from any of the SIGNAL, EMPTY, or REDUCED containers.
msg = ''
nzip = self.zip_container->count()
if nzip eq 0 then return,0
; Search for the descriptive variable for each and determine if they are the
; same.
oall = self.zip_container->get(/all)
first_var = strarr(nzip)
second_var = strarr(nzip)
for i = 0,nzip-1 do begin
   ret = oall[i]->get_property(descriptr_str = descriptr,   $
                               var_names = var_names,data = data,plot_str = plot_str $
                              ,second_var_name=sndVarName)
   ;col_pos = strpos(descriptr.name,':')
   ;index = where(var_names eq strmid(descriptr.name,0,col_pos))
   if i eq 0 then begin
      first_var_name = plot_str.xtitle
      second_var_name = descriptr.name ;var_names[index[0]]
   endif else begin
      first_var_name = [first_var_name,plot_str.xtitle]
      second_var_name = [second_var_name,descriptr.name];[second_var_name,var_names[index[0]]]
   endelse
   if n_elements(uniq(first_var_name)) ne 1  then begin
      msg = 'First independent variable names do not match up'
      return,0
  endif
;   if n_elements(uniq(second_var_name)) ne 1  then begin
;      msg = 'Combining variable names do not match up'
;      return,0
;  endif

endfor

dflag = 1
for i = 0,nzip-1 do begin
   ret = oall[i]->get_property(descriptr_str = descriptr,   $
                               var_names = var_names,data = data,plot_str = plot_str $
                              ,second_var_name=sndVarName)
   if n_elements(descriptr) eq 0 then dflag = 0
   ;col_pos = strpos(descriptr.name,':')
   ;index = where(var_names eq strmid(descriptr.name,0,col_pos))
   first_var_name = plot_str.xtitle
   second_var_name = descriptr.name ;var_names[index[0]]
   ; Get the second independent variable
   if i eq 0 then begin
      y = descriptr.qty ;*data.(index)
      x = *plot_str.x
      z = *plot_str.y
      dz = *plot_str.dy
   endif else begin
      y = [y,descriptr.qty]    ;[y,*data.(index)]
      x = [x,*plot_str.x]
      z = [z,*plot_str.y]
      dz = [dz,*plot_str.dy]
   endelse
endfor

; Now determine the grid
xmin = min(x,max = xmax)
ymin = min(y,max = ymax)
nygrid = n_elements(y)          ;fix(1.0*n_elements(y)/(1.0*nzip))
nxgrid = fix(1.0*n_elements(x)/(1.0*nzip))

dx = (xmax-xmin)/(nxgrid-1.0)
xgrid = xmin+dx*findgen(nxgrid)
dy = (ymax-ymin)/(nygrid-1.0)
ygrid = ymin+dy*findgen(nygrid)
limits = [xmin,ymin,xmax,ymax]
; Now go through and do the interpolation to a regular grid
triangulate, x, y, triangles, boundary
grid = trigrid(x, y, z, triangles,XOUT = xgrid,       $
   YOUT = ygrid,missing = 0.0);,max = max(z))
dgrid = trigrid(x, y, dz, triangles, XOUT = xgrid,    $
   YOUT = ygrid,missing = 0.0)
;grid = trigrid(x, y, z, triangles, limits,XGRID = xgrid, $
;   YGRID = ygrid,missing = 0.0,nx = 10);,max = max(z))
;dgrid = trigrid(x, y, dz, triangles, limits,XGRID = xgrid, $
;   YGRID = ygrid,missing = 0.0,ny = 10)
zout = grid
dzout = dgrid
xout = xgrid
yout = ygrid
;tvlct,r,g,b,/get
;loadct,1,/silent
;;tvimage,bytscl(zout)
;plotimage, bytscl((zout)), imgxrange=[min(xout),max(xout)], imgyrange=[min(yout),max(yout)]
;tvlct,r,g,b
;return,0

;for i = 0,nzip-1 do begin
;   ret = oall[i]->get_property(descriptr_str = descriptr_str,plot_str = plot_str)
;   ; If there is no descriptr present then fill it out
;   if dflag eq 0 then begin
;      ret = oall[i]->set_index_descriptr(float(i))
;      ret = oall[i]->get_property(descriptr_str = descriptr_str)
;   endif
;   ;ret = oall[i]->get_property(descriptr_str = descriptr_str)
;   if i eq 0 then begin
;      name = descriptr_str.name
;      qty = descriptr_str.qty
;      xtitles = plot_str.xtitle
;   endif else begin
;      name = [name,descriptr_str.name]
;      qty = [qty,descriptr_str.qty]
;      xtitles = [xtitles,plot_str.xtitle]
;   endelse
;endfor
;if n_elements(uniq(name)) gt 1 then begin
;   msg = 'Combining variables do not match!'
;   return,0
;endif
;  See if the first independent variable, in particular the name of x,
;  is the same for all data sets
;if n_elements(uniq(xtitles)) gt 1 then begin
;   msg = 'Independent variable types do not match!'
;   return,0
;endif
;dsort = sort(qty)
;; Perform a check on the independent variables
;ret = self->check_indep_var_limits()
;if ret eq 0 then return,0
;; Reorganize the zip container based on the sorting
;for i = 0,nzip-1 do self.zip_container->remove,oall[i]
;obj_destroy,self.zip_container
;self.zip_container = obj_new('IDL_CONTAINER')
;for i = 0,nzip-1 do self.zip_container->add,oall[dsort[i]]
;oall = self.zip_container->get(/all)
;nzip = self.zip_container->count()
;; Now put the final data set together
;for i = 0,nzip-1 do begin
;   ret = oall[i]->get_property(descriptr_str = descriptr_str, $
;      plot_str = plot_str,treatment = this_treatment)
;   if i eq 0 then treatment = this_treatment else $
;      treatment = [treatment,this_treatment]
;   ny = nzip
;   if i eq 0 then begin
;      nx = n_elements(*plot_str.x)
;      x = *plot_str.x
;      z = dblarr(nx,ny)
;      dz = dblarr(nx,ny)
;      y = dblarr(ny)
;   endif
;   xnew = *plot_str.x
;   y[i] = descriptr_str.qty
;   if total(xnew-x) ne 0 then begin ; interpolate
;      znew = *plot_str.y
;      dznew = *plot_str.dy
;      zinterp = interpol(znew,xnew,x)
;      dzinterp = interpol(dznew,xnew,x)
;      z[*,i] = zinterp
;      dz[*,i] = dzinterp
;   endif else begin
;      z[*,i] = *plot_str.y
;      dz[*,i] = *plot_str.dy
;      x = *plot_str.x
;   endelse
;endfor
;treatment = [treatment,'The files listed above were combined']
; The return value will be a structure
px = ptr_new(xout,/no_copy)
py = ptr_new(yout,/no_copy)
pz = ptr_new(transpose(zout),/no_copy)
bad_index = where(dzout eq 0.0,count_bad)
if count_bad gt 0 then dzout[bad_index] = 1.0
pdz = ptr_new(transpose(dzout),/no_copy)
xtitle = first_var_name;descriptr_str.legend
ytitle = second_var_name;plot_str.xtitle

str = {  px:px,               $
         py:py,               $
         pz:pz,               $
         pdz:pdz,             $
         xtitle:xtitle,       $
         ytitle:ytitle        }

o = obj_new('dtas_zip')
ret = o->set_property(  xval = *str.px,                  $
                        yval = *str.py,                  $
                        z = *str.pz,                     $
                        dz = *str.pdz,                   $
                        xtitle = xtitle,                 $
                        display_name = 'COMBINED_DATA',  $
                        treatment = treatment,           $
                        ytitle = ytitle                  )
heap_free,str
self.combine_container->add,o
return,o
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function dtas_reduce::reduce, SUM_SIGNAL = sum_signal,                        $
                              SUM_BG = sum_bg,                                $
                              STRAIGHT_ADD = straight_add,                    $
                              DETAILED_BALANCE = detailed_balance,            $
                              RESOLUTION_CORRECTION = resolution_correction,  $
                              MONITOR_CORRECTION = monitor_correction,        $
                              FAST_BACKGROUND = fast_background,              $
                              FBG_SET = fbg_set,                              $
                              INTERMEDIATE = intermediate,                    $
                              SUBTRACTION = subtraction,                      $
                              NORMALIZE_TO_MONITOR = normalize_to_monitor,    $
                              NORMALIZE_TO_time = normalize_to_time,    $
                              INTERP = interp,                                $
                              SF_SET = sf_set,                                $
                              SCALE_FACTOR = scale_factor,                    $
                              _EXTRA = extra,                                 $
                              MSG = msg
msg = ''
; The reduction steps are performed in the following order:
; (1).   If the SUM keyword is set then two possible sums are performed: (a). the
;        data from the SIGNAL container are summed and (b). the data from the
;        EMPTY container are summed.  If SUM is not set then all of the data objects
;        in the SIGNAL and EMPTY containers are treated identically.  If the
;        STRAIGHT_ADD keyword is set then both the contents of the SIGNAL container
;        and the EMPTY container are summed as raw counts (with no scaling to
;        monitor).  If STRAIGHT_ADD is not set then the data files are summed
;        in accordance with POISSON statistics.  The default is for straight
;        addition since we are assuming that the input data files are in terms of
;        raw counts.
; (2).   If the SUBTRACTION keyword is set then the sum of the EMPTY container is
;        subtracted from the sum of the SIGNAL container.  If the independent axes
;        do not line up then either LINEAR or (cubic) SPLINE interpolation are used.
;        The INTERP keyword can be set to "LINEAR" or "SPLINE" at the user's discretion
;        but LINEAR is the default.
; (3).   If FAST_BACKGROUND keyword is set equal to a number (assumed in counts/minute)
;        then this subtraction is applied for each value of the dependent variable.
;        Many of the subsequent reduction operations require the FAST_BACKGROUND
;        subtraction.  Such operations are noted explicitly below.
; (4).   If RESOLUTION_CORRECTION is set (Boolean) then the resolution correction
;        is applied.  Requires FAST_BACKGROUND subtraction.
; (5).   If DETAILED_BALANCE is set (Boolean) then the symmetrization of the
;        scattering profile is performed (if the temperature is provided in the
;        raw data file).  Requires FAST_BACKGROUND subtraction.
; (6).   If MONITOR_CORRECTION_FACTOR is applied then the data is corrected for
;        variations in the beam monitor which is provided in an external file
;        provided in this software.  Requires FAST_BACKGROUND subtraction.
; (7).   If NORMALIZE_TO_MONITOR is set then the scattering profile is normalized
;        to the beam monitor.
; (8).   If SCALE_FACTOR is set equal to a number then the dependent variable and
;        its uncertainty get scaled up by the scale factor provided.
;
signal_cont = self.signal_container & nsignal = signal_cont->count()
empty_cont = self.empty_container & nempty = empty_cont->count()

; Determine if there are any ICE files.  If so then provide a warning
; message and indicate that ICE files cannot be reduced yet.
;if nsignal gt 0 then begin
;   osignal = signal_cont->get(/all)
;   for i = 0,nsignal-1 do begin
;      ret = osignal[i]->get_property(filename = filename)
;      if is_ice(filename) then begin
;         msg = 'ICE files cannot be reduced yet.'
;         return,0B
;      endif
;   endfor
;endif
;if nempty gt 0 then begin
;   oempty = empty_cont->get(/all)
;   for i = 0,nempty-1 do begin
;      ret = oempty[i]->get_property(filename = filename)
;      if is_ice(filename) then begin
;         msg = 'ICE files cannot be reduced yet.'
;         return,0B
;      endif
;   endfor
;endif


if (nempty+nsignal) eq 0 then return,0

if n_elements(interp) ne 0 then self.interp = interp

; Sum the signal files?
if keyword_set(sum_signal) then begin
   ; If present then sum the SIGNAL files
   if nsignal gt 0 then begin
      if keyword_set(straight_add) then begin
         add_object = self->sum_straight(MSG = msg, o = oadd,/signal)
      endif else begin
         add_object = self->sum_poisson(MSG = msg, o = oadd,/signal)
      endelse
      if not obj_valid(add_object) then begin
         msg = msg
         return,0
      endif
      ret = self->update_plot_structure(oadd)
      ret = self->update_plot_structure(add_object)
      if intermediate then self.intermediate->add,oadd else obj_destroy,oadd
      red_object = add_object
   endif
endif

; Sum the background files?
if keyword_set(sum_bg) then begin
   if nempty gt 0 then begin
      if keyword_set(straight_add) then begin
         sub_object = self->sum_straight(MSG = msg, o = osub,/empty)
      endif else begin
         sub_object = self->sum_poisson(MSG = msg, o = osub,/empty)
      endelse
      if not obj_valid(sub_object) then begin
         msg = msg
         return,0
      endif
      ret = self->update_plot_structure(osub)
      ret = self->update_plot_structure(sub_object)
      if intermediate then self.intermediate->add,osub else obj_destroy,osub
      ;if n_elements(red_object) gt 0 then red_object = [red_object,sub_object] else $
      red_object = sub_object
   endif
endif

if obj_valid(add_object) then begin
   if ~obj_valid(sub_object) and nempty gt 0 then begin
      red_object = objarr(nempty+1)
      red_object[0] = add_object
      red_object[1:nempty] = (empty_cont->get(/all))->clone()
   endif
endif
if obj_valid(sub_object) then begin
   if ~obj_valid(add_object) and nsignal gt 0 then begin
      red_object = objarr(nsignal+1)
      oall = signal_cont->get(/all)
      for i = 0,n_elements(oall)-1 do red_object[i] = (oall[i])->clone()
      ;red_object[0:nsignal-1] = (signal_cont->get(/all))->clone()
      red_object[n_elements(oall)] = sub_object
   endif
endif

; If we have not summed any of the files then do this...
;if not (obj_valid(red_object))[0] then begin
if ((~obj_valid(add_object)) and (~obj_valid(sub_object))) then begin
   red_object = objarr(nsignal+nempty)
   if nsignal gt 0 then begin
      sig_all = signal_cont->get(/all)
      for i = 0,nsignal-1 do begin
         red_object[i] = sig_all[i]->clone()
      endfor
   endif
   if nempty gt 0 then begin
      empty_all = empty_cont->get(/all)
      for i = 0,nempty-1 do begin
         red_object[nsignal+i] = empty_all[i]->clone()
      endfor
   endif
   self->add,red_object
endif

; Perform the subtraction?
if keyword_set(subtraction) then begin
   if nsignal eq 1 and nempty eq 1 then begin
      add_object = (signal_cont->get(position = 0))->clone()
      sub_object = (empty_cont->get(position = 0))->clone()
      self->add,add_object
      self->add,sub_object
      if obj_valid(add_object) and obj_valid(sub_object) then begin
            subtraction_object = self->subtract_data(add_object,sub_object,   $
            MSG = msg, o = o_sub_object)
            ret = o_sub_object->get_property(plot_str = plot_str,data = data)
            if obj_valid(subtraction_object) then begin
               ret = self->update_plot_structure(o_sub_object)
               if intermediate then self.intermediate->add,o_sub_object else obj_destroy,o_sub_object
               red_object = subtraction_object
               ret = self->update_plot_structure(subtraction_object)
            endif
      endif
   endif
   if (nsignal gt 1) and (nempty eq 1) then begin
   ; Consider the case where we add the signal files first.
      sub_object = (empty_cont->get(position = 0))->clone()
      self->add,sub_object
      if obj_valid(add_object) then begin
         if obj_valid(add_object) and obj_valid(sub_object) then begin
            subtraction_object = self->subtract_data(add_object,sub_object,   $
            MSG = msg, o = o_sub_object)
            ret = o_sub_object->get_property(plot_str = plot_str,data = data)
            if obj_valid(subtraction_object) then begin
               ret = self->update_plot_structure(o_sub_object)
               if intermediate then self.intermediate->add,o_sub_object else obj_destroy,o_sub_object
               red_object = subtraction_object
               ret = self->update_plot_structure(subtraction_object)
            endif
         endif
      endif else begin
   ; Now consider the case where the one background file is subtracted
   ; from each individual signal file.
         n = nsignal + nempty
         new_red_object = objarr(n-1)
         for i = 0,n-2 do begin
            add_object = (signal_cont->get(position = i))->clone()
            self->add,add_object
            if obj_valid(add_object) and obj_valid(sub_object) then begin
               subtraction_object = self->subtract_data(add_object,sub_object,   $
                  MSG = msg, o = o_sub_object)


               ret = o_sub_object->get_property(plot_str = plot_str,data = data)
               if obj_valid(subtraction_object) then begin
                  ret = self->update_plot_structure(o_sub_object)
                  if intermediate then self.intermediate->add,o_sub_object else obj_destroy,o_sub_object
                  new_red_object[i] = subtraction_object
                  ret = self->update_plot_structure(subtraction_object)
               endif
            endif
         endfor
         red_object = new_red_object
      endelse
   endif

;;   if (nsignal gt 1) and (nempty gt 1) then begin
   if (nsignal ge 1) and (nempty gt 1) then begin
       ;; RTA - add option to perform subtraction of equivalent bkgd
       ;;       from signal files when multiple signal/bkgd are loaded
       if (nsignal eq nempty) and ~keyword_set(sum_bg) then begin
           ;; if same nos of signal and bkgd files then do direct
           ;; subtraction of equivalent files in display order
           err_msg = ['Subtracting corresponding background from signal datasets', $
                     'in the order in which they are listed!']
           void = dialog_message(err_msg,/information)


           n = nsignal + 1
           new_red_object = objarr(n-1)
           for i = 0,n-2 do begin
               add_object = (signal_cont->get(position = i))->clone()
               self->add,add_object
               sub_object = empty_cont->get(position = i)
               if obj_valid(add_object) and obj_valid(sub_object) then begin
                   subtraction_object = self->subtract_data(add_object,sub_object,   $
                                                            MSG = msg, o = o_sub_object)
                   ret = o_sub_object->get_property(plot_str = plot_str,data = data)
                   if obj_valid(subtraction_object) then begin
                       ret = self->update_plot_structure(o_sub_object)
                       if intermediate then self.intermediate->add,o_sub_object else obj_destroy,o_sub_object
                       new_red_object[i] = subtraction_object
                       ret = self->update_plot_structure(subtraction_object)
                   endif
               endif
           endfor
           red_object = new_red_object
       endif else begin
           
           
           if ~obj_valid(sub_object) then begin
               sub_object = (empty_cont->get(position = 0))->clone()
               self->add,sub_object
               err_msg = 'Using the first file in the background folder only in'+$
                         'the subtraction'
               void = dialog_message(err_msg,/information)
           endif
           ;; Consider the case where we add the signal files first.
           if obj_valid(add_object) then begin
               subtraction_object = self->subtract_data(add_object,sub_object,   $
                                                        MSG = msg, o = o_sub_object)
               ret = o_sub_object->get_property(plot_str = plot_str,data = data)
               if obj_valid(subtraction_object) then begin
                   ret = self->update_plot_structure(o_sub_object)
                   if intermediate then self.intermediate->add,o_sub_object else obj_destroy,o_sub_object
                   red_object = subtraction_object
                   ret = self->update_plot_structure(subtraction_object)
               endif
           endif else begin
               ;; Now consider the case where the sum of the background files is subtracted
               ;; from each individual signal file.
               n = nsignal + 1
               new_red_object = objarr(n-1)
               for i = 0,n-2 do begin
                   add_object = (signal_cont->get(position = i))->clone()
                   self->add,add_object
                   if obj_valid(add_object) and obj_valid(sub_object) then begin
                       subtraction_object = self->subtract_data(add_object,sub_object,   $
                                                                MSG = msg, o = o_sub_object)
                       ret = o_sub_object->get_property(plot_str = plot_str,data = data)
                       if obj_valid(subtraction_object) then begin
                           ret = self->update_plot_structure(o_sub_object)
                           if intermediate then self.intermediate->add,o_sub_object else obj_destroy,o_sub_object
                           new_red_object[i] = subtraction_object
                           ret = self->update_plot_structure(subtraction_object)
                       endif
                   endif
               endfor
               red_object = new_red_object
           endelse
       endelse
   endif
   
endif else begin

   if obj_valid(add_object) and obj_valid(sub_object) then begin
      red_object = objarr(2)
      red_object[0] = add_object & red_object[1] = sub_object
   endif

endelse


n_total = n_elements(red_object)
if keyword_set(fbg_set) then begin
   for i = 0,n_total-1 do begin
   ; For all subsequent corrections the fast background correction
   ; is essential.
      ret = red_object[i]->fast_background_subtraction(fast_background,o = o1)
      if ((not intermediate) and (obj_valid(o1))) then obj_destroy,o1
      ret = self->update_plot_structure(red_object[i])
      if obj_valid(o1) then begin
         if intermediate then self.intermediate->add,o1
         ret = self->update_plot_structure(red_object[i])
      endif
      if keyword_set(resolution_correction) then begin
         ret = red_object[i]->resolution_correction(o = o2)
         if (not intermediate) and obj_valid(o2) then obj_destroy,o2
      endif
      ret = self->update_plot_structure(red_object[i])
      if obj_valid(o2) then begin
         if intermediate then self.intermediate->add,o2
         ret = self->update_plot_structure(o2)
      endif
      if keyword_set(detailed_balance) then begin
         ret = red_object[i]->detailed_balance(o = o3)
         if (not intermediate) and obj_valid(o3) then obj_destroy,o3
      endif
         ret = self->update_plot_structure(red_object[i])
      if obj_valid(o3) then begin
         if intermediate then self.intermediate->add,o3
         ret = self->update_plot_structure(o3)
      endif
      if keyword_set(monitor_correction) then begin
         ret = red_object[i]->monitor_correction_factor(o = o4)
         if ((not intermediate) and (obj_valid(o4))) then obj_destroy,o4
      endif

      ret = self->update_plot_structure(red_object[i])
      if obj_valid(o4) then begin
         if intermediate then self.intermediate->add,o4
         ret = self->update_plot_structure(o4)
      endif
      if keyword_set(normalize_to_monitor) then begin
         ret = red_object[i]->normalize_to_monitor(o = o5)
         if (not intermediate) and obj_valid(o5) then obj_destroy,o5
     endif else begin
         if (keyword_set(normalize_to_time)) then begin
             ret = red_object[i]->normalize_to_time(o = o5)
             if (not intermediate) and obj_valid(o5) then obj_destroy,o5
         endif
     endelse

     ret = self->update_plot_structure(red_object[i])
     if obj_valid(o5) then begin
         if intermediate then self.intermediate->add,o5
         ret = self->update_plot_structure(o5)
     endif
   endfor
endif

if keyword_set(sf_set) ne 0 then begin
;   i = n_total - 1
   for i = 0,n_total-1 do begin
      ret = red_object[i]->calibration(scale_factor,o = o6)
      ret = self->update_plot_structure(red_object[i])
      if (not intermediate) and obj_valid(o6) then obj_destroy,o6
      if obj_valid(o6) then begin
         if intermediate then self.intermediate->add,o6
            ret = self->update_plot_structure(o6)
      endif
   endfor
endif

if intermediate then nobj = self.intermediate->count()
if intermediate then oall = self.intermediate->get(/all)
ret = red_object[n_total-1]->get_property(treatment = treatment)
reduced = self.reduced_container
reduced_objects = reduced->get(/all)
for i = 0,n_elements(red_object)-1 do begin
   if obj_valid(reduced_objects[0]) then begin
      where_same = where(reduced_objects eq red_object[i],count_same)
      if (count_same eq 0) then $
         self.reduced_container->add,red_object[i]
   endif else begin
      for i = 0,n_total-1 do self.reduced_container->add,red_object[i]
   endelse
endfor

;intermediate = self.intermediate
if n_elements(red_object) ne 0 then begin
   for i = 0,n_elements(red_object)-1 do begin
      if obj_valid(red_object[i]) then begin
         ret = self->update_plot_structure(red_object[i])
         ret = red_object[i]->get_property(data = data)
         ncols = n_tags(data)
         slen = max(strlen(tag_names(data)))
         nx = n_elements(*data.(0))
         contents = strarr(nx+1)
         fstr = strtrim(string(slen),2)
         str = strtrim(string(tag_names(data),format = fs),2)
         fs = '('+strtrim(string(ncols),2)+'A'+fstr+')';'('+strtrim(string(ncols),2)+'A22)'
         str = strtrim(string(tag_names(data),format = fs),2)
         contents[0] = str
         f = '('+strtrim(string(ncols),2)+'e'+fstr+'.5)';'('+strtrim(string(ncols),2)+'e22.5)'
         for j = 0,nx-1 do begin
            for k = 0,ncols-1 do begin
               if k eq 0 then this_array = (*data.(k))[j] else $
                  this_array = [this_array,(*data.(k))[j]]
            endfor
            str = strtrim(string(this_array,format = f),2)
            contents[j+1] = str
         endfor
         ret = red_object[i]->set_property(signal_type = 'REDUCED',  $
            contents = contents)
      endif
   endfor
endif
if intermediate then begin
   nint = self.intermediate->count()
   oall = self.intermediate->get(/all)
   for i = 0,nint-1 do begin
      if obj_valid(oall[i]) then begin
         ret = self->update_plot_structure(oall[i])
         ret = oall[i]->get_property(data = data)
         ncols = n_tags(data)
         nx = n_elements(*data.(0))
         contents = strarr(nx+1)
         fstr = strtrim(string(slen),2)
         str = strtrim(string(tag_names(data),format = fs),2)
         fs = '('+strtrim(string(ncols),2)+'A'+fstr+')';'('+strtrim(string(ncols),2)+'A22)'
         str = strtrim(string(tag_names(data),format = fs),2)
         contents[0] = str
         f = '('+strtrim(string(ncols),2)+'e'+fstr+'.5)';'('+strtrim(string(ncols),2)+'e22.5)'
         for j = 0,nx-1 do begin
            for k = 0,ncols-1 do begin
               if k eq 0 then this_array = (*data.(k))[j] else $
                  this_array = [this_array,(*data.(k))[j]]
            endfor
            str = strtrim(string(this_array,format = f),2)
            contents[j+1] = str
         endfor
         ret = oall[i]->set_property(signal_type = 'INTERMEDIATE',  $
            contents = contents)
      endif
   endfor
endif
if n_elements(red_object) eq 0 then red_object = 0
return,red_object
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function dtas_reduce::add,oarray
for i = 0,n_elements(oarray)-1 do begin
   if obj_class(oarray[i]) eq 'DTAS_DATA' then $
      self->add,oarray[i]
   ret = oarray[i]->get_property(signal_type = signal_type)
   case signal_type of
   'SIGNAL':  self.signal_container->add,oarray[i]
   'EMPTY':   self.empty_container->add,oarray[i]
   else:
   endcase
endfor
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function dtas_reduce::check_contents
; This function goes through the contents of the SIGNAL container and the EMPTY
; container and determines if each of the contents' signal_types matches the
; container type (EMPTY or SIGNAL).
oall = self->get(/all)
o_signal = self.signal_container->get(/all)
o_empty = self.empty_container->get(/all)
nsignal = self.signal_container->count()
nempty = self.empty_container->count()
for i = 0,nsignal-1 do begin
   ret = o_signal[i]->get_property(signal_type = signal_type)
   if signal_type eq 'EMPTY' then begin
      self.signal_container->remove,o_signal[i]
      self.empty_container->add,o_signal[i]
   endif
endfor
for i = 0,nempty-1 do begin
   ret = o_empty[i]->get_property(signal_type = signal_type)
   if signal_type eq 'SIGNAL' then begin
      self.empty_container->remove,o_empty[i]
      self.signal_container->add,o_empty[i]
   endif
endfor
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function dtas_reduce::init,obj_array
if not self->idl_container::init() then return,0
self.signal_container = obj_new('IDL_CONTAINER')
self.empty_container = obj_new('IDL_CONTAINER')
if n_elements(obj_array) ne 0 then begin
; Put the objects into the appropriate container
   nobj = n_elements(obj_array)
   for i = 0,nobj-1 do begin
      ret = obj_array[i]->get_property(signal_type = signal_type)
      case signal_type of
      'SIGNAL':   self.signal_container->add,obj_array[i]
      'EMPTY':    self.empty_container->add,obj_array[i]
      else:
      endcase
   endfor
   self->add,obj_array
endif
self.intermediate = obj_new('IDL_CONTAINER')
self.reduced_container = obj_new('IDL_CONTAINER')
self.zip_container = obj_new('IDL_CONTAINER')
self.combine_container = obj_new('IDL_CONTAINER')
self.interp = 'LINEAR'  ; (default)
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro dtas_reduce__define
define = {  dtas_reduce,                  $
            interp:'',                    $
            intermediate:obj_new(),       $
            signal_container:obj_new(),   $
            empty_container:obj_new(),    $
            reduced_container:obj_new(),  $
            zip_container:obj_new(),      $
            combine_container:obj_new(),  $
            inherits IDL_CONTAINER  }
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
