; Written by  Craig Brown
; adapted from J.R.D. Copley
;************************************************************************************************
;###############################################################################
;
; 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 or if the code in this file is
;  included in another product.
;
;###############################################################################
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro ddcs_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()
ndarkcount = self.darkcount_container->count()
nvanadium = self.vanadium_container->count()
nabsorber = self.absorber_container->count()
nreduced = self.reduced_container->count()
nint = self.intermediate->count()
nzip = self.zip_container->count()
ncorrected = self.corrected_container->count()
if ncorrected gt 0 then begin
   oall = self.corrected_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 nvanadium gt 0 then begin
   oall = self.vanadium_container->get(/all)
   obj_destroy,oall
endif
if nabsorber gt 0 then begin
   oall = self.absorber_container->get(/all)
   obj_destroy,oall
endif
if ndarkcount gt 0 then begin
   oall = self.darkcount_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.corrected_container,$
   self.vanadium_container,self.darkcount_container,self.absorber_container]
ptr_free,self.signal_normalization_values_ptr,$
   self.vanadium_normalization_values_ptr,$
   self.empty_normalization_values_ptr,$
   self.absorber_normalization_values_ptr,$
   self.darkcount_normalization_values_ptr,$
   self.signal_fbnormalization_values_ptr,$
   self.vanadium_fbnormalization_values_ptr,$
   self.empty_fbnormalization_values_ptr,$
   self.absorber_fbnormalization_values_ptr,$
   self.darkcount_fbnormalization_values_ptr,$
   self.grouping_ptr,       $
   self.angle_props_ptr,    $
   self.elasticqvals_ptr,   $
   self.Q_props_ptr,        $
   self.qflagPtr

self->IDL_CONTAINER::cleanup
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ddcs_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
'DARKCOUNT': $
   begin
      self.darkcount_container->remove,o
      obj_destroy,o
   end
'ABSORBER': $
   begin
      self.absorber_container->remove,o
      obj_destroy,o
   end
'VANADIUM': $
   begin
      self.vanadium_container->remove,o
      obj_destroy,o
   end
'REDUCED': $
   begin
      self.reduced_container->remove,o
      obj_destroy,o
   end
'CORRECTED': $
   begin
      self.corrected_container->remove,o
      obj_destroy,o
   end
else:
endcase
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ddcs_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 ddcs_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 ddcs_reduce::get_property, INTERMEDIATE = intermediate,  $
                                    REDUCED = reduced,            $
                                    SIGNAL = signal,              $
                                    VANADIUM = vanadium,          $
                                    DARKCOUNT = darkcount,        $
                                    ABSORBER = absorber,          $
                                    INTERP = interp,              $
                                    CORRECTED = CORRECTED,          $
                                    EMPTY = empty
corrected = self.corrected_container
interp = self.interp
reduced = self.reduced_container
signal = self.signal_container
empty = self.empty_container
darkcount = self.darkcount_container
vanadium = self.vanadium_container
intermediate = self.intermediate
absorber = self.absorber_container
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ddcs_reduce::clear_int
; This function clears the objects from the int container and destroys them.
nint = self.intermediate->count()
if nint gt 0 then begin
   oall = self.intermediate->get(/all)
   obj_destroy,oall
   obj_destroy,self.intermediate
   self.intermediate = obj_new('IDL_CONTAINER')
endif
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ddcs_reduce::clear_com
; This function clears the objects from the int container and destroys them.
ncom = self.corrected_container->count()
if ncom gt 0 then begin
   oall = self.corrected_container->get(/all)
   obj_destroy,oall
   obj_destroy,self.corrected_container
   self.corrected_container = obj_new('IDL_CONTAINER')
endif
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ddcs_reduce::Scatter, _EXTRA = extra,$
                     type=type,$
                     temperatures=temperatures,$
                     dodetailed=dodetailed,$
                     dodeteff=dodeteff,$
                     dozeroenergy=dozeroenergy,$
                     qty=qty,err=err,$
                     qelmean=qelmean,tthmean=tthmean,$
                     xaxis=xaxis,xunits=xunits,xlabel=xlabel
sample_index=-99
o_corrected = self.corrected_container->get(/all)
ncorrected = self.corrected_container->count()
for i = 0 , ncorrected-1 do begin
    ret=o_corrected[i]->get_property(display_name=display_name)
    if display_name eq 'CORRECTED_SIGNAL' then sample_index=i
endfor
if sample_index lt 0 then  return,0

logentry=""
; Values of kT in meV.
if type eq 'GDOS' then begin
   ;make the temperature some value
   if (size(temperatures))[2] eq 0 then temperature =300.0
endif
if (size(temperatures))[2] ne 0 then begin
if temperatures le 0.0 then temperatures=300.
   if temperatures gt 0.0 then begin
      ktemperatures=temperatures*!dcs_boltzmann/!dcs_meVcon
   endif
endif

if (dodeteff) then begin
    message="Performing detector efficiency corrections."
endif else begin
    message="SKIPPING detector efficiency corrections."
endelse
logentry=logentry+string(10B)+message

ret = o_corrected[sample_index]->get_property(n0ngtchan=n0ngtchan,time_props=time_props,ch_wl=ch_wl,$
                                angle_props=angle_props,Q_props=Q_props, qty=qty, err=err)

;ch_wl=ch_wl,n0ngtchan=n0ngtchan,ch_ms=ch_ms,ch_srdenom=ch_srdenom,$
 ;                           normalization_sum=normalization_sum, nhedet=nhedet)
ngtchan0=n0ngtchan

; Uncorrected times of flight and time channel boundaries
times_uncorr=time_props[0:ngtchan0]
tstart_uncorr=times_uncorr[0:ngtchan0-1]
tstop_uncorr=times_uncorr[1:ngtchan0]
;
e0=!dcs_hsq2mn/ch_wl^2
;
tsd0=ch_wl/!dcs_hom*!dcs_dsd

; If requested, determine zero energy channel time correction.
    if (dozeroenergy) then begin
       peakfrac=0.5
; Check to see if dimension is 1 or 2 XXX will die if 1D Q vector
       if (size(qty))[0] eq  1 then summed_int=qty else summed_int=total(qty,2)
       skip=0
       maxsummed_int=max(summed_int)
       maxch=where(maxsummed_int eq summed_int,count)
       if (count ne 1) then begin
           print,"More than one channel has peak intensity"
           tcorr=0.0
           skip=1
       endif
       channels=where(summed_int ge peakfrac*maxsummed_int,nchannels)
       if (nchannels le 3) then begin
           print,"Too few channels in range determined by peakfrac."
           tcorr=0.0
           skip=1
       endif
       if (skip ne 1) then begin
          lchannel=channels[0]
          rchannel=channels[nchannels-1]
          tmean=(tstart_uncorr+tstop_uncorr)*0.5
          twid=tstop_uncorr-tstart_uncorr
          moment0=total(summed_int[lchannel:rchannel]*twid[lchannel:rchannel])
          moment1=total(summed_int[lchannel:rchannel]*twid[lchannel:rchannel]*(tmean[lchannel:rchannel]-tmean[maxch[0]]))
          tcentroid=moment1/moment0+tmean[maxch[0]]
          tcorr=tsd0-tcentroid
       endif
    endif else begin
       tcorr=0.0
    endelse

; Corrected times and energies.
    times=times_uncorr+tcorr
    energies=e0-!dcs_hsq2mn*!dcs_dsd^2/!dcs_hom^2/times^2
;
; Boundaries of corrected time and energy channels.
    tstart=times[0:ngtchan0-1]
    tstop=times[1:ngtchan0]
    etstart=energies[0:ngtchan0-1]
    etstop=energies[1:ngtchan0]
;
; Mean times of flight from sample to detector
    tmean=0.5*(tstart+tstop)
;
; Mean energy transfers
    etmean=0.5*(etstart+etstop)

; Mean scattering angles and elastic Q values
    tthmean=angle_props[*,0]
    qelmean=Q_props[*,0]
    ntth=n_elements(tthmean)

; Values of 2*cos(2theta) and sin(2theta)
    twocostth=2.0*cos(tthmean*!dtor)
    sintth=sin(tthmean*!dtor)
;
;
    ef=e0-etmean
    root=sqrt(e0*ef)
;
; Ratios (k_i/k_f)
    kikfratio=root/ef
;
; Inverse of detector efficiency as function of Ef.
    if (dodeteff) then begin
       det_thick=1.0; assumed detector thickness in cm.
       sigma_acoeff=0.4466; (from Jeremy in cm-1)
       sigma_s=0.0008; (from Jeremy in cm-1)
       wl0=ch_wl
       deteff_inv=1.0/(1.0-exp(-det_thick*(sigma_s+sigma_acoeff*wl0*kikfratio)))
    endif else begin
       deteff_inv=fltarr(ngtchan0)+1
    endelse

    unitvectth=fltarr(ntth)+1

; Create 2-d array of inelastic Q values
; Next line creates matrices with dimensions [ngtchan0,ntth].
    hsqqsq=dcs_operate((e0+ef)#unitvectth,root#twocostth,"-")
    conversion=2*!pi/sqrt(!dcs_hsq2mn)
    qinel=sqrt(hsqqsq)*conversion
;
; Create 1-d array of (maximum Q^4 - minimum Q^4) values
    tthmin=10.0
    tthmax=140.0
    q4maxmin=(e0+ef-2.0*root*cos(tthmax*!dtor))^2 - (e0+ef-2.0*root*cos(tthmin*!dtor))^2

; Values of exp(-E/2kT)
if dodetailed  then begin
     if ((size(temperatures))[2] ne 0) then begin
       dbal=exp(-etmean/(2.0*ktemperatures))
     endif else dbal=1.0
endif else begin
    dbal=1.0
endelse
;
; Lines containing matrix multiplication (#) create matrices with
; dimensions [ngtchan0,ntth].

;     "GDOS": factor=((deteff_inv*kikfratio^4*(1-dbal^2)*etmean)#unitvectth)/hsqqsq
; NOTE VDOS IS DIFFERENT!

case type of
      "DODT": $
         begin
             factor=deteff_inv#unitvectth
             qty=factor * qty
             err=factor * err
             xaxis=times
             xunits="time:usec"
             xlabel="Time-of-flight"

         end
      "DODE":  $
         begin
             factor=(deteff_inv*kikfratio^3)#unitvectth
             qty=factor * qty
             err=factor * err
             xaxis=energies
             xunits="energy:meV"
             xlabel="Energy transfer"

         end
      "SQOM":  $
         begin
             factor=(deteff_inv*kikfratio^4)#unitvectth
             qty=factor * qty
             err=factor * err
             xaxis=energies
             xunits="energy:meV"
             xlabel="Energy transfer"
         end
      "SSQW":  $
         begin
             factor=(deteff_inv*kikfratio^4*dbal)#unitvectth
             qty=factor * qty
             err=factor * err
             xaxis=energies
             xunits="energy:meV"
             xlabel="Energy transfer"
         end
      "KIQW":  $
         begin
             factor=(deteff_inv*kikfratio^4*(1-dbal^2))#unitvectth
             qty=factor * qty
             err=factor * err
             xaxis=energies
             xunits="energy:meV"
             xlabel="Energy transfer"
         end
      "GDOS":  $
         begin
             factor=((deteff_inv*kikfratio^4*(1-dbal^2)*etmean)#unitvectth)/hsqqsq
             qty=factor * qty
             err=factor * err
             xaxis=energies
             xunits="energy:meV"
             xlabel="Energy transfer"
             if (size(qty))[0] eq  1 then qty=qty else  qty=total(qty,2)
             if (size(err))[0] eq  1 then err=err else err=sqrt(total(err^2,2))
         end
      "VDOS":  $
         begin
             factor=(deteff_inv*kikfratio^3*(1-dbal^2)*etmean/q4maxmin)#sintth
             qty=factor * qty
             err=factor * err
             xaxis=energies
             xunits="energy:meV"
             xlabel="Energy transfer"
             qty=total(qty,2)
             err=sqrt(total(err^2,2))
         end
      "SGDOS":  $
         begin
             factor=((deteff_inv*kikfratio^4*(1-dbal^2)*etmean)#unitvectth)/hsqqsq
             qty=factor * qty
             err=factor * err
             xaxis=energies
             xunits="energy:meV"
             xlabel="Energy transfer"
             qty=total(qty,2)
             err=sqrt(total(err^2,2))
         end
      else:
    endcase


    if (dozeroenergy) then begin
       message="Zero energy channel correction,"+strcompress(string(tcorr))+" us."
    endif else begin
       message="NO zero energy channel correction."
    endelse
    logentry=logentry+string(10B)+message

    bigfloat=exp(80.0)
    indices=where(1-finite(qty),countq)
    if (countq gt 0) then begin
       qty[indices]=bigfloat
    endif
    indices=where(1-finite(err),counte)
    if (counte gt 0) then begin
       err[indices]=bigfloat
    endif

return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ddcs_reduce::reduce, _EXTRA = extra,                                 $
                              mask = mask,                                    $
                              msg = msg,                                      $
                              dark=dark,                                      $
                              vanarr = vanarr,                                $
                              q_group = q_group,                              $
                              no_det_grouping=no_det_grouping,                $
                              attenuation = attenuation,                      $
                              norm=norm,                                      $
                              log=log,                                        $
                              slider=slider
self.log=log
msg = ''
; The COPLEY reduction steps are performed in the following order:
; (1).  load all data files  (DONE IN THE GUI)
;       -move them between groups, do not let the same file be repeated in a group.
;       -this has only been implemented for signal files
; (2).  perform sanity checks, unlike John, you dont have a choice about the level! (PART OF THE DATA OBJECT)
; (3).  check their "sameness", "dcs_sameness" or skip (I ALWAYS DO HERE)
; (4).  Call a procedure that typically returns the beam monitor counts and durations for each run.
;       These quantities are used for normalization purposes. No errors are calculated for these
;       quantities.  dcs_normalization_counts-
;       normalization_values{nruns} and fbnormalization_values{nruns} John creates.
; (5)   Have the user define a detector mask. dcs_choose_mask. mask{nhedet,nruns} defined
;       IMPLEMENTED IN THE PLOT WINDOW. CHOOSING MASK ON PLOT?
; (6)   Sum "identical" files.  Also sum the two sets of normalization values.Mask(s) are also applied.
;       The wavelength is the same for all runs other than the background file(s).
;       dcs_sumidenticals
; (7)   At this stage we generally wish to save the summed file names and wavelength
;       information for possible use in a later reduction session. We need a name for
;       the file where this information is to be saved. We either use the default file
;       name or request a file name (with the option to not save).
;           fullfilename=!dcs_outp_directory+!dcs_rootname+".sfi"
;           save,filename=fullfilename,lambda,sumfiles,n_sumfiles
; (8)   Construct normalized background, if there is one.
;       fbdata_qty=fltarr(nhedet)    fbdata_err=fltarr(nhedet)   fbnormalization=fltarr(nhedet)
;       dcs_bgd_subtract
; (9)   Have the user define detector efficiencies.
;       dcs_choose_det_eff  "det_eff_qty{nhedet}" and "det_eff_err{nhedet}"
; (10)  Have the user define detector groupings.
;       dcs_choose_grouping   "grouping{nhedet}","angle_props{ngroups,3}" and
;       "Q_props{ngroups,3}" have been defined (rows 1, 2, 3 are mean, min and max).
;       Also doQgrouping is defined.]
; (11)  Perform corrections, i.e. "fast background" subtraction, detector grouping, detector
;       efficiency correction, data normalization, corrections for container and absorber runs.
;       dcs_docorrections  "samfiles{n}" has been defined: n is the number of corrected sample files created
;     and samfiles contains the names of the corrected sample files created.
; (12)  Convert data to d2sigma/domegadt or some other quantityand save
;       Convert data to d2sigma/domegadt or some other quantity


; dark is a [float,float] darkcount rate/error  /detector/second

signal_cont = self.signal_container & nsignal = signal_cont->count()
empty_cont = self.empty_container & nempty = empty_cont->count()
vanadium_cont = self.vanadium_container & nvanadium = vanadium_cont->count()
darkcount_cont = self.darkcount_container & ndarkcount = darkcount_cont->count()
absorber_cont = self.absorber_container & nabsorber = absorber_cont->count()

if (nsignal) eq 0 then return,0

*self.attenuation_ptr=attenuation

;check that all files in each folder are similar tags:
ret=self->check_contents()

;check that all pertinent data is consistent
ret=self->check_sameness(lines=lines)
;need to do something besides just dump the program if they are not the same.
;possible put in an override flag?
if (n_elements(lines) ge 1) and (lines[0] ne '') then begin
   msg='Mismatch of  '+lines
   return,0
endif

;get normalization values
ret=self->dcs_normalization(norm=norm)

;sum identicals
ret=self->dcs_sumindeticals(mask=mask)

; Construct the time-averaged epithermal background for later subtraction.
ret=self->dcs_bgd_subtract(dark=dark)


;Detector efficiencies
ret = self->dcs_det_eff(vanarr=vanarr)

; Q or 2th groupings from the q_group structure
;q_group={use_det_grouping:use_det_grouping,NEGPOS:NEGPOS,LOMIUP:LOMIUP,LOWUPP:LOWUPP,DONOTC:DONOTC,$
;        EQUAL2:EQUAL2, EQUALQ:EQUALQ, TTSTEP:TTSTEP, QSTEP:QSTEP}
;
ret = self -> dcs_choose_q(q_group=q_group,no_det_grouping=no_det_grouping,slider=slider)


; Do corrections
ret = self -> dcs_corrections()
;this is where the corrected data object is created.

if ptr_valid(self.corrected_container) eq 0 then begin
    msg=[msg+'Done!']
    return,msg
endif else begin
    return,0
endelse

return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ddcs_reduce::dcs_corrections

renorm_factor=1.0e6
logentry=''
; This is the number of files to be treated (S+C+A)
o_int = self.intermediate->get(/all)
nint = self.intermediate->count()
ngroups = self.ngroups

sample_index=-1
abs_index=-1
empty_index=-1
nfiles=0
for i = 0 , nint-1 do begin
    ret=o_int[i]->get_property(display_name=display_name)
    if display_name eq 'SUMMED_SIGNAL' then begin
        sample_index=i
        nfiles=nfiles+1
    endif
    if display_name eq 'SUMMED_ABSORBER' then begin
         abs_index=i
         nfiles=nfiles+1
    endif
    if display_name eq 'SUMMED_EMPTY' then begin
          empty_index=i
          nfiles=nfiles+1
    endif
endfor
;

logentry=logentry+string(10B)+$
    "Background subtraction, detector grouping, detector efficiency correction, "+$
    "data normalization and corrections using empty can and cadmium runs have been performed."
logentry=logentry+string(10B)+$
    "Data renormalized. Multiplicative factor = "+strcompress(string(renorm_factor))+"."

grpfiles=strarr(1)       ; JOHN HAS AN ARRAY HERE FOR MULTIPLE FILE CONDITIONS
temperatures=fltarr(1)   ;    "                        "                     "

; Read the first file in order to get at quantities such as ch_wl, nhedet etc.
ret = o_int[sample_index]->get_property(ch_wl=ch_wl,n0ngtchan=n0ngtchan,ch_ms=ch_ms,ch_srdenom=ch_srdenom,$
                            normalization_sum=normalization_sum, nhedet=nhedet)
ngtchan0=n0ngtchan

; Calculate time between pulses at sample.
speriod=6.0e7/ch_ms*ch_srdenom

if (self.doQgrouping) then begin
     ret = self ->dcs_redefine_props(ch_wl=ch_wl,nhedet=nhedet)
     ret = self ->dcs_define_qflagPtr(ch_wl=ch_wl,nhedet=nhedet,ngtchan0=ngtchan0)
endif

;
; The normalization sum is the same for all unmasked detectors and 0 otherwise.
; These two lines reconstruct the detector masking vector M(d).
masking=normalization_sum/max(normalization_sum)
;
; Construct eta-hat(t,g).
ret = self -> dcs_make_eta(ngtchan0=ngtchan0,nhedet=nhedet,$
            sumeff_qty=sumeff_qty,sumeff_err=sumeff_err)
;will this return the values needed sumeff_qty,sumeff_err ? XXX


; Construct I_B-hat(t,g)/tau (scalar zeros if T_B^tot=0).
ret = self -> dcs_make_ibhat(ngtchan0=ngtchan0,nhedet=nhedet,speriod=speriod,$
            sumfb_qty=sumfb_qty,sumfb_err=sumfb_err)

;
f_a=0.0
absfile_qty=0.0
absfile_err=0.0
fa_tnorm=0.0
;
f_c=0.0
canfile_qty=0.0
canfile_err=0.0
fc_tnorm=0.0
;
; Loop over files, E- A- S-
order=[empty_index,abs_index,sample_index]

for i=0,n_elements(order)-1 do begin
    if order[i] le -1 then CONTINUE
    ret = o_int[order[i]]->get_property(display_name=display_name,normalization_sum=normalization_sum,$
                fbnormalization_sum=fbnormalization_sum,xval=xval,qty=qty,err=err)

    icomp=nfiles-1-i
;
; LETS DUPLICATE THE SUMMED_SAMPLE FILE FOR NOW AND SEE.
;
;;;     ret = o_new->set_property(angle_props=*self.angle_props_ptr,Q_props=*self.Q_props_ptr)
;
; The quantities norm and tnorm are N_X^tot and T_X^tot/N_X^tot respectively.
;   We neglect their errors.
    norm=max(normalization_sum)
    tnorm=max(fbnormalization_sum)/norm
;
;   qtyptr points to I_X(t,d), errptr points to its error, xptr points to x.
; Dereference these quantities, free the pointer, redefine the pointer to
; reflect arrays dimensioned ngroups rather than nhedet, populate the structure.
    qtyPtr=ptr_new(qty)
    errPtr=ptr_new(err)
    xPtr=ptr_new(xval)

    qty=fltarr(ngtchan0,ngroups)
    err=fltarr(ngtchan0,ngroups)
    x=*xPtr
    y=(*self.Q_props_Ptr)[*,0]

;
; In what follows, qty is I_X-hat(t,g)/N_X^tot and err is the
;   corresponding error.
;
    case ngroups of
;   One group.
       1: begin
         qty = total((*qtyPtr),2)/norm
         err = sqrt(total((*errPtr)^2,2))/norm; this line fixed 5/8/02
       end
;
; No grouping (hence nhedet groups).
       nhedet: begin
         qty = *qtyPtr/norm
         err = *errPtr/norm
       end
;
;   Otherwise...
       else: begin
;
; Grouping in 2theta (rather than in Q):
         if (not self.doQgrouping) then begin
          for k=0,ngroups-1 do begin
              list=where(*self.grouping_ptr eq k+1)
              qty[*,k] = total((*qtyPtr)[*,list],2)/norm
              err[*,k] = sqrt(total( (*errPtr)[*,list]^2,2 ))/norm
          endfor
         endif else begin
;
; Grouping in Q (rather than in angle):
;          device,decomposed=0
;          loadct,13
;
          progressBar = Obj_New($
              "SHOWPROGRESS",$
              title="Work in progress. Please wait.",$
              xsize=300,$
              ysize=10,$
              message="Please be patient...",$
              color=184)
          progressBar->Start
;
          for k=0,ngroups-1 do begin
              qty[*,k] = total((*qtyPtr)*((*self.qflagPtr)[k,*,*]),2)/norm
              err[*,k] = sqrt(total( (*errPtr)^2*((*self.qflagPtr)[k,*,*]),2 ))/norm
              progressBar->Update,10+(k+1)/float(ngroups)*90
          endfor
;
          progressBar->Destroy
          Obj_Destroy, progressBar
;
;device,decomposed=1
         endelse
       end
    endcase
;
    if (display_name eq 'SUMMED_EMPTY') then begin
          f_c=*self.attenuation_ptr
          canfile_qty=f_c*qty
          canfile_err=f_c*err
          fc_tnorm=f_c*tnorm
    endif
;
    if (display_name eq 'SUMMED_ABSORBER') then begin
          f_c=*self.attenuation_ptr
          f_a=1.0-f_c
          absfile_qty=f_a*qty
          absfile_err=f_a*err
          fa_tnorm=f_a*tnorm
   endif
;
    if (display_name ne 'SUMMED_ABSORBER' and display_name ne 'SUMMED_EMPTY') then begin

;it seems that the attenuation factor is used twice.....WHAT is the formula?

; LETS DUPLICATE THE SUMMED_SAMPLE FILE FOR NOW AND SEE.
;create an object in which to copy the relevent data.
    newtype='CORRECTED_SIGNAL'
    ;create a new object to contain everything:
     o_new=o_int[order[i]]->clone()
     ret = o_new->set_property(signal_type = newtype,display_name=newtype)
;   Put angle_props and Q_props into the DAVE structure.
     ret = o_new->set_property(angle_props=*self.angle_props_ptr,Q_props=*self.Q_props_ptr)

;--------------------------------------------------------------------------------
          scorr_qty=((qty - absfile_qty - canfile_qty) - $
            sumfb_qty*(tnorm - fa_tnorm - fc_tnorm))/sumeff_qty
          scorr_err=sqrt(((err^2 + (f_a*absfile_err)^2 + (f_c*canfile_err)^2) + $
             sumfb_err^2*(tnorm-fa_tnorm-fc_tnorm)^2))/sumeff_qty
;---------------------------------------------------------------------------------
;this is Johns original. It does f_c^2
;--------------------------------------------------------------------------------
;          scorr_qty=((qty - f_a*absfile_qty - f_c*canfile_qty) - $
;            sumfb_qty*(tnorm - fa_tnorm - fc_tnorm))/sumeff_qty
;          scorr_err=sqrt(((err^2 + (f_a*absfile_err)^2 + (f_c*canfile_err)^2) + $
;             sumfb_err^2*(tnorm-fa_tnorm-fc_tnorm)^2))/sumeff_qty
;---------------------------------------------------------------------------------
; Test for divide by 0.0; replace +/-NaN with 0.0
       indices=where(finite(scorr_qty) eq 0,count)
       if (count gt 0) then begin
         scorr_qty[indices]=0.0
         scorr_err[indices]=0.0
       endif
       qty=scorr_qty*renorm_factor
       err=scorr_err*renorm_factor
;
       ret = o_new->set_property(qty=qty,err=err)
;
       ylabel = 'mean Q'
       yunits = 'wavevector:A-1'
       ytype = 'points'
       histlabel = 'Grouped Intensity'
       ret = o_new->get_property(treatment=treatment)
       lines=strsplit(logentry,string(10B),/extract)
       treatment=[treatment,lines]
;
       ret = o_new->set_property(ylabel=ylabel,yunits=yunits,ytype=ytype,histlabel=histlabel,treatment=treatment)
;
    endif
;
    ptr_free,qtyPtr
    ptr_free,errPtr
    ptr_free,xPtr
endfor
;
; NB qflagPtr (if it exists) is a pointer to an array (ngroups,1000,913).
if (total(ptr_valid(qflagPtr)) gt 0) then ptr_free,qflagPtr
;
self.log=self.log+string(10B)+logentry
;this is where the log file is added to the treatment pointer

ret = o_new->get_property(treatment=treatment)

lines=strsplit(self.log,string(10B),/extract)
treatment=[[treatment],lines]

ret = o_new->set_property(treatment=treatment)

;MOVE THIS TO THE OVERALL STRUCTURE
     if obj_valid(o_new) then begin
        print,'WE HAVE A CORRECTED SIGNAL OBJECT'
        self.corrected_container->add,o_new
     endif else begin
        obj_destroy,o_new
        return,0
     endelse

return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ddcs_reduce::dcs_make_ibhat,ngtchan0=ngtchan0,nhedet=nhedet,speriod=speriod,$
                                    sumfb_qty=sumfb_qty,sumfb_err=sumfb_err

; This procedure constructs the quantity I_B-hat(t,g)/tau.
; Quantities returned are dimensioned ngtchan0*ngroups if tnorm ne 0.
; If tnorm=0 (i.e. T_B^tot = 0), scalar zeros are returned.
;
;get other parameters

compile_opt strictarr
;
grouping=*self.grouping_ptr
ngroups=max(*self.grouping_ptr)
doQgrouping=self.doQgrouping
fbdata_qty=self.fbdata_qty
fbdata_err=self.fbdata_err
fbnormalization=self.fbnormalization

tnorm=max(fbnormalization)
;
if (tnorm eq 0.0) then begin
    sumfb_qty=0.0
    sumfb_err=0.0
    return,1
endif
;
unitvec_t=intarr(ngtchan0)+1
;
case ngroups of
;   One group.
    1: begin
       sumfb_qty = unitvec_t*total(fbdata_qty)/speriod
       sumfb_err = unitvec_t*sqrt(total(fbdata_err^2))/speriod
    end
;
; No grouping (hence nhedet groups).
    nhedet: begin
       sumfb_qty = (unitvec_t#fbdata_qty)/speriod
       sumfb_err = (unitvec_t#fbdata_err)/speriod
    end
;
;   Otherwise...
    else: begin
;
; Grouping in 2theta (rather than in Q):
       if (not doQgrouping) then begin
         sumfb_qty=fltarr(ngtchan0,ngroups)
         sumfb_err=fltarr(ngtchan0,ngroups)
         for k=0,ngroups-1 do begin
          list=where(grouping eq k+1)
          sumfb_qty[*,k] = (unitvec_t#(total(fbdata_qty[list])))/speriod
          sumfb_err[*,k] = (unitvec_t#(sqrt(total(fbdata_err[list]^2))))/speriod
         endfor
       endif else begin

; Grouping in Q (rather than in angle):
         sumfb_qty=fltarr(ngtchan0,ngroups)
         sumfb_err=fltarr(ngtchan0,ngroups)
         for k=0,ngroups-1 do begin
          sumfb_qty[*,k] = total((unitvec_t#fbdata_qty)*((*self.qflagPtr)[k,*,*]),2)/speriod
          sumfb_err[*,k] = sqrt(total((unitvec_t#fbdata_err^2)*((*self.qflagPtr)[k,*,*]),2))/speriod
         endfor
       endelse
    end
endcase
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ddcs_reduce::dcs_make_eta,ngtchan0=ngtchan0,nhedet=nhedet,$
                                    sumeff_qty=sumeff_qty,sumeff_err=sumeff_err
;Quantities returned are dimensioned ngtchan0*ngroups.
;get other parameters
;det_eff_qty, det_eff_qty, qflagPtr


unitvec_t=intarr(ngtchan0)+1
grouping=*self.grouping_ptr
ngroups=max(*self.grouping_ptr)
det_eff_qty=self.det_eff_qty
det_eff_err=self.det_eff_err

doQgrouping=self.doQgrouping

case ngroups of
;   One group.
    1: begin
       sumeff_qty = unitvec_t*total(det_eff_qty)
       sumeff_err = unitvec_t*sqrt(total(det_eff_err^2 ))
    end
;
; No grouping (hence nhedet groups).
    nhedet: begin
       sumeff_qty = unitvec_t#det_eff_qty
       sumeff_err = unitvec_t#det_eff_err
    end
;
;   Otherwise...
    else: begin
;
; Grouping in 2theta (rather than in Q):
       if (not doQgrouping) then begin
         sumeff_qty=fltarr(ngtchan0,ngroups)
         sumeff_err=fltarr(ngtchan0,ngroups)
         for k=0,ngroups-1 do begin
          list=where(grouping eq k+1)
          sumeff_qty[*,k] = unitvec_t#total(det_eff_qty[list])
          sumeff_err[*,k] = unitvec_t#sqrt(total(det_eff_err[list]^2))
         endfor
       endif else begin

; Grouping in Q (rather than in angle):
         sumeff_qty=fltarr(ngtchan0,ngroups)
         sumeff_err=fltarr(ngtchan0,ngroups)
         for k=0,ngroups-1 do begin
          sumeff_qty[*,k] = total((unitvec_t#det_eff_qty)*((*self.qflagPtr)[k,*,*]),2)
          sumeff_err[*,k] = sqrt(total((unitvec_t#det_eff_err^2)*((*self.qflagPtr)[k,*,*]),2 ))
         endfor
       endelse
    end
endcase

return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ddcs_reduce::dcs_redefine_props,ch_wl=ch_wl,nhedet=nhedet
; This procedure is called if the grouping is in Q.
; It redefines angle_props and Q_props.
; Note that (*Q_propsPtr)[*,0] will become the y component in histPtr.
;
compile_opt strictarr
;
;
lambda=ch_wl
e0=float(!dcs_hsq2mn/lambda^2)
conversion=float(4.*!pi^2/!dcs_hsq2mn)
k0=sqrt(conversion*e0)
;
Q_props=*self.Q_props_ptr
angles=self.angles
ngroups=(self.ngroups)



; The variable qrat is a 2-column array containing absolute (Q/k_0) limits for each group.
qrat=transpose(abs(Q_props[*,1:2])/k0)
;
; The variables min2th and max2th contain minimum and maximum values of abs(2theta).
;nhedet=(*(*(*davePtr).dataStrPtr).specificPtr).nhedet
;dcs_getangles,angles,nhedet=nhedet
;get angles from elsewhere...

min2th=min(abs(angles))
max2th=max(abs(angles))
cmin2th=cos(min2th*!dtor)
cmax2th=cos(max2th*!dtor)
argmin=cmin2th^2+qrat^2-1
argmax=cmax2th^2+qrat^2-1
;
; Expand Q_props to 7 rows.
; Rows 1-3 contain the mean, minimum and maximum Q for the group.
; Rows 4-5 contain the minimum and maximum energy transfer for the minimum Q in the group.
; Rows 6-7 contain the minimum and maximum energy transfer for the maximum Q in the group.
Q_props=[[Q_props],[fltarr(ngroups,4)]]
; Compute the quantities that go in rows 4-7 of Q_props.
jk=-1
for j=0,ngroups-1 do begin
    for k=0,1 do begin
       jk=jk+1
       elim=fltarr(2)
       nr=-1
       if (argmin[jk] gt 0.0) then begin
         roomin=sqrt(argmin[jk])
         if (qrat[jk] gt 1.0) then begin
          nr=nr+1
          elim[nr]=e0*(cmin2th+roomin)^2
         endif
         if (qrat[jk] le 1.0 and cmin2th gt 0.0) then begin
          nr=nr+1
          elim[nr]=e0*(cmin2th-roomin)^2
          nr=nr+1
          elim[nr]=e0*(cmin2th+roomin)^2
         endif
       endif
       if (argmax[jk] gt 0.0) then begin
         roomax=sqrt(argmax[jk])
         if (qrat[jk] gt 1.0) then begin
          nr=nr+1
          elim[nr]=e0*(cmax2th+roomax)^2
         endif
         if (qrat[jk] le 1.0 and cmax2th gt 0.0) then begin
          nr=nr+1
          elim[nr]=e0*(cmax2th-roomax)^2
          nr=nr+1
          elim[nr]=e0*(cmax2th+roomax)^2
         endif
       endif
       elim=e0-elim
       elim=elim[sort(elim)]
       Q_props[j,3+k*2]=elim[0]
       Q_props[j,3+k*2+1]=elim[1]
    endfor
endfor

*self.Q_props_ptr=Q_props


return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ddcs_reduce::dcs_define_qflagPtr,nhedet=nhedet,ch_wl=ch_wl,ngtchan0=ngtchan0
;ngroups,elasticqvals,davePtr=davePtr,
;what we want qflagPtr=qflagPtr

;get the old colors:


;device,decomposed=0
;loadct,13
;
progressBar = Obj_New($
    "SHOWPROGRESS",$
;   /cancelbutton,$;  (How does this work?)
    title="Work in progress. Please wait.",$
    xsize=300,$
    ysize=10,$
    message="Please be patient...",$
    color=184)
progressBar->Start
;
;nhedet=(*(*(*davePtr).dataStrPtr).specificPtr).nhedet
;dcs_getangles,angles,nhedet=nhedet
angles=self.angles
ngroups=max(*self.grouping_ptr)

cos2th=cos((!pi/180.)*angles)
unitvec_d=intarr(nhedet)+1
;
lambda=ch_wl
e0=float(!dcs_hsq2mn/lambda^2)
;
;ngtchan0=(*(*(*davePtr).dataStrPtr).specificPtr).ngtchan0
;efinal=e0-float((*(*(*(*davePtr).dataStrPtr).specificPtr).energy_propsptr)[0:ngtchan0-1])
;need to fetch this guy...
o_int = self.intermediate->get(/all)
nint = self.intermediate->count()
for i = 0 , nint-1 do begin
    ret=o_int[i]->get_property(display_name=display_name)
    if display_name eq 'SUMMED_SIGNAL' then index=i
endfor
ret = o_int[index]->get_property(energy_props=energy_props)

efinal=e0-float(energy_props[0:ngtchan0-1])
unitvec_t=intarr(ngtchan0)+1
;
elasticqvals=*self.elasticqvals_ptr
;
conversion=float(4.*!pi^2/!dcs_hsq2mn)
qsqptr=ptr_new(conversion*(e0+efinal#unitvec_d-2.*sqrt(e0*efinal)#cos2th),/no_copy)
sign_of_2th=(angles gt 0)*2-1
if (elasticqvals[0] lt 0.0) then $
    inelasticqPtr=ptr_new(sqrt(*qsqptr)*(unitvec_t#sign_of_2th),/no_copy) else $
    inelasticqPtr=ptr_new(sqrt(*qsqptr),/no_copy)
ptr_free,qsqptr
;
elasticqvals=float(elasticqvals)
tdgPtr=ptr_new(bytarr(ngtchan0,nhedet))
*tdgPtr=*inelasticqPtr ge elasticqvals[0] and *inelasticqPtr lt elasticqvals[1]
;
progressBar->Update,10
;
for k=1,ngroups-1 do begin
    *tdgPtr=*tdgPtr+(k+1)*$
       (*inelasticqPtr ge elasticqvals[k] and *inelasticqPtr lt elasticqvals[k+1])
    progressBar->Update,10+k/float(ngroups-1)*80
endfor
;
ptr_free,inelasticqPtr
;
qflagPtrs=ptrarr(ngroups,/allocate_heap,/nozero)
;
test=fltarr(ngroups,(size(*tdgPtr))[1],(size(*tdgPtr))[2])

progressBar->Update,100
;
for k=0,ngroups-1 do begin
    *qflagPtrs[k]=(*tdgPtr eq k)
    test[k,*,*]=(*tdgPtr eq k)
endfor
;
ptr_free,tdgPtr
;
progressBar->Destroy
Obj_Destroy, progressBar
;
;device,decomposed=1
;
;self.qflagPtr=*qflagPtrs
*self.qflagPtr=test
ptr_free,qflagPtrs

return,1
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ddcs_reduce::dcs_choose_q,q_group=q_group,no_det_grouping=no_det_grouping,$
                        slider=slider

;load the detector info
inputfile=!DAVE_AUXILIARY_DIR+'dcs_DetPos2.dat'
doQgrouping=self.doQgrouping  ; default -1

content=self.content
if content[0] EQ '' then begin  ; i.e. we need only read this once.
   openr,unit,inputfile,/get_lun
   readf,unit,content
   free_lun,unit
   self.content=content
endif

o_int = self.intermediate->get(/all)
nint = self.intermediate->count()
for i = 0 , nint-1 do begin
    ret=o_int[i]->get_property(display_name=display_name)
    if display_name eq 'SUMMED_SIGNAL' then index=i
endfor

ret = o_int[index]->get_property(ch_wl=ch_wl,nhedet=nhedet)

grouping=indgen(nhedet)+1
;
tthmin=-29.99
tthmax=139.99
qmin=4.0*!pi*sin(tthmin*0.5/!radeg)/ch_wl
qmax=4.0*!pi*sin(tthmax*0.5/!radeg)/ch_wl
;

bank=[" 1"," 0","-1"]
contents=strarr(6,nhedet+1)
contents[0:4,0]=["Detector","Bank","Angle","''Q''","Group"]
contents[0,1:nhedet]=[string(grouping)]
;
angles=fltarr(nhedet)
Q_values=fltarr(nhedet)
;
for i=1,nhedet do begin
    pieces=strsplit(content[i]," ",/extract)
    contents[1,i]=bank[fix(pieces[6])+1]
    contents[2,i]=pieces[8]
    angles[i-1]=float(pieces[8])
    Q_values[i-1]=4.0*!pi*sin(angles[i-1]*0.5/!radeg)/ch_wl
    contents[3,i]=string(Q_values[i-1],format="(f7.3)")
    contents[4,i]=contents[0,i]
endfor
self.angles=angles
self.Q_values=Q_values

;convert my values into things similar to Johns.

; np=0 to NOT "Combine negative and positive angle detectors"
; np=1 to     "Combine negative and positive angle detectors"
; uml=1 to "Combine lower, middle and upper banks"
; uml=2 to "Combine lower and upper banks only")
; uml=3 to "Do not combine detector banks"
; qthchce=1 to "Combine into groups with equal widths in 2theta" (calc_vals contains 2theta step)
; qthchce=2 to "Combine into groups with equal widths in ''Q''" (calc_vals contains Q step)
; qthchce=3 to "Specify 2theta boundaries of detector groups" (calc_vals contains 2theta values)
; qthchce=4 to "Specify ''Q'' boundaries of detector groups" (calc_vals contains Q values)
; qthchce=1 2theta groups
; qthchce=2 Q grouping
;   calc_vals contains quantities described above
;   contents contains the contents of "detpos2.dat"
; nhedet=913
; lambda=incident wavelength
; tthmin,tthmax are basically -30 and +140
; qmin,qmax are corresponding Q limits
;

;noq=q_group.use_det_grouping
np=q_group.NEGPOS
if q_group.LOMIUP eq 1 then uml=1
if q_group.LOWUPP eq 1 then uml=2
if q_group.DONOTC eq 1 then uml=3
if q_group.EQUAL2 eq 1 then begin
    qthchce=1
    calc_vals=q_group.TTSTEP
endif else begin
    if q_group.EQUALQ eq 1 then begin
       qthchce=2
       calc_vals=q_group.QSTEP
    endif
endelse
;
lambda=ch_wl;
;
;;;;ADDED CB

if  (no_det_grouping eq 1) then begin
;     Treat case that "DO NOT GROUP" button was checked, hence no grouping.
;   grouping=transpose(grouping); turn grouping into a row vector.

   grouping=indgen(nhedet)+1
   ngroups=max(grouping)
   angle_props=fltarr(ngroups,3)
   Q_props=fltarr(ngroups,3)

   ;additionally, check the attenuation value to see if we need to interpolate.
   if n_elements(*self.attenuation_ptr) gt 1 then begin
      if n_elements(*self.attenuation_ptr) ne nhedet then begin
         ;need to interpolate
         ; for now send a message and make value =1.0
         msg='Interpolation of SSFs not implemented, Value set to 1.0'
         void = dialog_message(dialog_parent = self.group_leader,msg,  $
         title = 'Problems with DCS Data Reduction',/information)
         *self.attenuation_ptr=1.0
      endif
         ;we can then do direct subtractions in the corrections function
         ;it seems that there needs to be an additional inout saying if we are using
         ;plate or cylinder geometry as they do not have the same function.
   endif

   if (np) then begin
     angles=abs(angles)
     Q_values=abs(Q_values)
   endif

   for ig=0,ngroups-1 do begin
       list=where(grouping eq ig+1)
       angle_props[ig,0]=mean(angles[list])
       Q_props[ig,0]=mean(Q_values[list])
       angle_props[ig,1]=min(angles[list])
       Q_props[ig,1]=min(Q_values[list])
       angle_props[ig,2]=max(angles[list])
       Q_props[ig,2]=max(Q_values[list])
   endfor
   elasticqvals=q_props[*,0]
   if (ngroups eq 1 or ngroups eq nhedet or doQgrouping eq 1) then doQgrouping=0 else begin
      doQgrouping = 1
   endelse
endif else begin
;     Treat case that "DO NOT GROUP" button was unchecked, hence  grouping.
    grouping=fix(contents[0,1:nhedet])
    bank=fix(contents[1,1:nhedet])
    angle=float(contents[2,1:nhedet])

    case qthchce of
    1: begin
       tthstep=float(calc_vals[0])
       if (tthstep le 0.0) then begin
         res=dialog_message("Invalid 2theta step size:"+calc_vals[0],/error)
         grouping=0
         return,0
       endif
       if (np eq 0) then nneg=fix(-tthmin/tthstep+1.0) else nneg=0
       npos=fix(tthmax/tthstep+1.0)
       nvals=nneg+1+npos
       vals=fltarr(nvals)
       if (np eq 0) then vals[0]=tthmin else vals[0]=0.0
       if (nvals gt 2) then for i=1,nvals-2 do vals[i]=tthstep*(i-nneg)
       vals[nvals-1]=tthmax
    end
    2: begin
       qstep=float(calc_vals[0])
       if (qstep le 0.0) then begin
         res=dialog_message("Invalid Q step size:"+calc_vals[0],/error)
         grouping=0
         return,0
       endif
         if (np eq 0) then nneg=fix(-qmin/qstep+1.0) else nneg=0
       npos=fix(qmax/qstep+1.0)
       nvals=nneg+1+npos
       vals=fltarr(nvals)
       if (np eq 0) then vals[0]=tthmin else vals[0]=0.0
       if (nvals gt 2) then for i=1,nvals-2 do vals[i]=2.*asin(qstep*(i-nneg)*lambda/(4.0*!pi))*!radeg
       vals[nvals-1]=tthmax
    end
    else:
   endcase
;
; Hence elasticqvals, which contains the limits on Q(elastic).
   elasticqvals=4.0*!pi*sin(vals*!pi/360.0)/lambda

  if (np eq 0) then begin
    for i=0,nneg-1 do begin
       vmin=vals[nneg-i-1]
       vmax=vals[nneg-i]
       indices=where((angle ge vmin) and (angle lt vmax))
       if (indices[0] gt -1) then grouping[indices]=1+i*uml
    endfor
    for i=nneg,nvals-2 do begin
       vmin=vals[i]
       vmax=vals[i+1]
       indices=where((angle ge vmin) and (angle lt vmax))
       if (indices[0] gt -1) then grouping[indices]=1+i*uml
    endfor
   endif else begin
    for i=0,nvals-2 do begin
       vmin=vals[i]
       vmax=vals[i+1]
       indices=where((abs(angle) ge vmin) and (abs(angle) lt vmax))
       if (indices[0] gt -1) then grouping[indices]=1+i*uml
     endfor
   endelse

   if (uml eq 2) then grouping=grouping+abs(bank)
;
   if (uml eq 3) then grouping=grouping+1+bank
;
   grouping=transpose(grouping); turn grouping into a row vector.
;
; The final grouping array should include at least one instance of each integer from 1 to ngroups.
   dcs_general_collapse,grouping,cgrouping
   if (n_elements(cgrouping) eq 1) then message,'Fatal error in dcs_calculate_grouping.'
   grouping=cgrouping+1
   ngroups=max(grouping)
;we have groups, np, elasticqvals, ngroups,

   angle_props=fltarr(ngroups,3)
   Q_props=fltarr(ngroups,3)
;
   if (np) then begin
       angles=abs(angles)
       Q_values=abs(Q_values)
   endif

   for ig=0,ngroups-1 do begin
       list=where(grouping eq ig+1)
       angle_props[ig,0]=mean(angles[list])
       Q_props[ig,0]=mean(Q_values[list])
       angle_props[ig,1]=min(angles[list])
       Q_props[ig,1]=min(Q_values[list])
       angle_props[ig,2]=max(angles[list])
       Q_props[ig,2]=max(Q_values[list])
   endfor
;

   if (ngroups eq 1 or ngroups eq nhedet or qthchce eq 1) then doQgrouping=0 else begin
        doQgrouping = 0
   endelse

endelse

   ; q_group.EQUAL2 = 1
   ; qthchce=1
   ; q_group.EQUALQ=0

self.ngroups=ngroups
*self.grouping_ptr=grouping
*self.angle_props_ptr=angle_props
*self.Q_props_ptr=Q_props
self.doQgrouping=doQgrouping
*self.elasticqvals_ptr=elasticqvals

log='Number of detector groups = '+string(ngroups)
log=log+string(10B)+'Angle Mean,        Min,           Max'
for i=0,n_elements(angle_props)/3-1 do begin
   log=log+string(10B)+string(angle_props[i,0])+string(angle_props[i,1])+string(angle_props[i,2])
endfor
log=log+string(10B)+'   Q  Mean,       Min,              Max'
for i=0,n_elements(Q_props)/3-1 do begin
   log=log+string(10B)+string(Q_props[i,0])+string(Q_props[i,1])+string(Q_props[i,2])
endfor
log=log+string(10B)+'Elastic Q values'
for i=0,n_elements(elasticqvals)-1 do begin
   log=log+string(10B)+string(elasticqvals[i])
endfor

self.log=self.log+string(10B)+log

;just for the GUI slider
slider=ngroups
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ddcs_reduce::dcs_det_eff,vanarr=vanarr

logentry=""
;get some kind of preference for vanadium.
;set it for now

if n_elements(vanarr) eq 913 then begin ;use this array
    logentry=logentry+string(10B)+"Detector efficiencies are from a file."
    self.det_eff_qty=vanarr
    self.det_eff_err=vanarr*0
    return,1
endif else begin
    ; we will check to see if there is a file, if not use unity values
    vanchoice=1
endelse

;check to see if there is a summed_vanadium file
if (vanchoice eq 1) then begin
      o_int = self.intermediate->get(/all)
      nint = self.intermediate->count()
    for i = 0 , nint-1 do begin
        ret=o_int[i]->get_property(display_name=display_name)
        if display_name eq 'SUMMED_VANADIUM' then index=i
    endfor
    test=size(index,/type)
    if (test eq 0) then begin
       logentry=logentry+string(10B)+$
       "There is no vanadium file so efficiencies cannot be computed."
       vanchoice=0
    endif
endif

; The variable "method" controls how the detector efficiencies are calculated.
; method=0 means the total counts are computed for each detector. (NOT AVAILABLE)
; method=1 means the sum of the elastic peaks is fitted to a Gaussian.
; method=2 means elastic peaks are individually fitted to Gaussians. (NOT AVAILABLE)
method=1

case vanchoice of
    0: begin
       logentry=logentry+string(10B)+"Detector efficiencies have been set to unity."
       osignal = self.signal_container->get(/all)
       ret = osignal[0]->get_property(nhedet=nhedet)
       self.det_eff_qty=fltarr(nhedet)+1
       self.det_eff_err=fltarr(nhedet)
    end
    1: begin
       ret=self->dcs_define_det_eff(method=method)
       logentry=logentry+string(10B)+"Detector efficiencies have been computed."
    end
    else:
endcase

self.log=self.log+string(10B)+logentry
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ddcs_reduce::dcs_define_det_eff,method=method
o_int = self.intermediate->get(/all)
nint = self.intermediate->count()
for i = 0 , nint-1 do begin
    ret=o_int[i]->get_property(display_name=display_name)
    if display_name eq 'SUMMED_VANADIUM' then index=i
endfor
test=size(index,/type)
if (test eq 0) then return,0
logentry=''
; Define quantities needed to obtain detector peak positions and intensities.
; Master speed, master period, period at sample.

ret = o_int[index]->get_property(ch_ms=ch_ms,ch_srdenom=ch_srdenom,$
                    ch_res=ch_res,ch_wl=ch_wl,det_dis=det_dis,$
                    ch_dis=ch_dis,ch_slots=ch_slots,startchoice=startchoice,$
                    tchanlook=tchanlook,n0ngtchan=n0ngtchan,qty=qty,$
                    tsdmin=tsdmin,ngtchan1=ngtchan1,Q_props=Q_props)

mspeed=ch_ms
mperiod=60000000.0/mspeed; master period
speriod=mperiod*ch_srdenom

;   Detector distance, incident wavelength, incident inverse velocity, resolution mode
lambda=ch_wl; incident wavelength
tau0=lambda/!dcs_hom; incident inverse velocity
resmode=ch_res; resolution mode

;   Chopper distances and slot positions.
distances=ch_dis
angles=ch_slots
;
;    3-character STARTCHOICE
startchoice=strmid(startchoice,0,3)
;
;    Define time channel widths for high end data (in us).
tchanwid=tchanlook[1,*]*0.05d0
;
;      Start and stop times for each channel (in us). Channel 0 starts at 0.
tstop=reform(total(tchanwid,/cumulative))
tstart=shift(tstop,1)
tstart[0]=0

;   Number of time channels.
nch=n0ngtchan
;
;   Truncate tstart and tstop.
tstart=tstart[0:nch-1]
tstop=tstop[0:nch-1]
;
; Centers of time channels (for use in gaussfit)
tmean=0.5*(tstart+tstop); (for use in gaussfit)

; Pointer to input data.
indataPtr=ptr_new(double(qty[0:nch-1,*]))

; Method 0 uses the total counts for each detector.
;
if (method eq 0) then begin
    det_eff_qty=total(*indataPtr,1)
    help,det_eff_qty
    norm=mean(det_eff_qty)
    self.det_eff_qty=det_eff_qty/norm
    ptr_free,indataPtr
    return,1
endif

; Methods 1 and 2 use elastic peak sums.
;
; Compute ts0, which is the time that the DAQ start pulse occurs,
;   relative to the time that neutrons cross the sample axis.
;
case startchoice of
    "AND": ts0=tsdmin
    "CH1": ts0=speriod+tau0*distances[0]-(30-angles[0,resmode-1])*(speriod/360.0)
    "CH2": ts0=speriod+tau0*distances[1]-(30-angles[1,resmode-1])*(speriod/360.0)
    "CH3": ts0=speriod+tau0*distances[2]-(30-angles[2,resmode-1])*(speriod/360.0)
    "CH4": ts0=speriod+tau0*distances[3]-(30-angles[3,resmode-1])*(speriod/360.0)
    "CH5": ts0=speriod+tau0*distances[4]-(30-angles[4,resmode-1])*(speriod/360.0)
    "CH6": ts0=speriod+tau0*distances[5]-(30-angles[5,resmode-1])*(speriod/360.0)
    "CH7": ts0=speriod+tau0*distances[6]-(30-angles[6,resmode-1])*(speriod/360.0)
    else: return,0
endcase
;
;
;   Compute elastic time of flight from sample to detector.
t_sd=tau0*det_dis

;
; Hence elastic time of arrival at detector, relative to DAQ start.
;td_the=dave_modulus(t_sd-ts0,speriod) (OLD VERSION)
td_the=(t_sd-ts0) mod speriod + speriod * ((t_sd-ts0) lt 0.0)
;
; Compute corresponding time channel.
tch_the=where(tstart le td_the and td_the lt tstop)
tch_the=tch_the[0]
;if (!debug) then print,'Theoretical elastic time channel = ',tch_the
;
; Renumber channels so that peak count is roughly in channel nch/2.
; The number of channels shifted is tch_shift.
;tch_shift=dave_modulus(round(-tch_the+nch/2),nch) (OLD VERSION)
tch_shift=round(-tch_the+ngtchan1/2) mod ngtchan1 + ngtchan1 * (round(-tch_the+ngtchan1/2) lt 0.0)
;if (!debug) then print,'tch_shift = ',tch_shift
;
; Rotate data so it peaks in channel nch/2.
*indataPtr=shift(*indataPtr,tch_shift,0)
;
; Redefine mean times for time channels.
tmean_rot=shift(tmean,tch_shift)
addition=fltarr(ngtchan1)+speriod
if (tch_shift ne 0) then addition[0:tch_shift-1]=0
tmean_rot=tmean_rot+addition
;
; Compute time shift corresponding to tch_shift.

;tshift=tstop[dave_modulus(ngtchan1-1-tch_shift,ngtchan1)] (OLD VERSION)
tshift=tstop[(ngtchan1-1-tch_shift) mod ngtchan1 + ngtchan1 * ((ngtchan1-1-tch_shift) lt 0)]
;
nterms=5
params=fltarr(nterms)
;
; Method 1 uses elastic peak sums based on an overall fit of a Gaussian plus
; constant background to the rotated counts summed over all detectors.
;
; mult determines how far out we go (how many multiples of sigma), either side of
; the peak position; the background is then estimated by giong out a further
; mult*sigma in either direction.
mult=3
if (method eq 1) then begin
    sumdata_rot=total(*indataPtr,2)
    maxcounts=sumdata_rot[round(ngtchan1/2)]
    estimates=[maxcounts,tmean_rot[round(ngtchan1/2)],20.0,0.0,0.0]
    result=gaussfit(tmean_rot,sumdata_rot,params,estimates=estimates,nterms=nterms)
    lotime=params[1]-mult*params[2]
    loch=where(tmean_rot-10.0 lt lotime and tmean_rot+10.0 ge lotime)
    loch=round(mean(loch))
    hitime=params[1]+mult*params[2]
    hich=where(tmean_rot-10.0 lt hitime and tmean_rot+10.0 ge hitime)
    hich=round(mean(hich))
    nchans=hich-loch+1
    intint=total((*indataPtr)[loch:hich,*],1)
    lo_int=total((*indataPtr)[loch-nchans+1:loch,*],1)
    hi_int=total((*indataPtr)[hich:hich+nchans-1,*],1)
    det_eff_qty=intint-0.5*(lo_int+hi_int)
    norm=mean(det_eff_qty)
    det_eff_qty=det_eff_qty/norm
    logentry=logentry+string(10B)+$
     "Detector efficiencies computed using elastic peak sums."
endif
;
; Method 2 uses elastic peak sums based on individual fits of a Gaussian plus
; constant background to the rotated counts for each of the detectors.
if (method eq 2) then begin
    lines=["Detector"+strcompress(1)+" of"+strcompress(nhedet)]
    ;dcs_message_display,lines,base,group_leader=dcs_base,alignment=0
    for j=0,nhedet-1 do begin
       maxcounts=(*indataPtr)[round(ngtchan1/2),j]
       estimates=[maxcounts,tmean_rot[round(ngtchan1/2)],20.0,0.0,0.0]
       result=gaussfit(tmean_rot,(*indataPtr)[*,j],params,estimates=estimates,nterms=nterms)
       lotime=params[1]-mult*params[2]
       loch=where(tmean_rot-10.0 lt lotime and tmean_rot+10.0 ge lotime)
       loch=round(mean(loch))
       hitime=params[1]+mult*params[2]
       hich=where(tmean_rot-10.0 lt hitime and tmean_rot+10.0 ge hitime)
       hich=round(mean(hich))
       nchans=hich-loch+1
       intint=total((*indataPtr)[loch:hich,j])
       lo_int=total((*indataPtr)[loch-nchans+1:loch,j])
       hi_int=total((*indataPtr)[hich:hich+nchans-1,j])
       det_eff_qty[j]=intint-0.5*(lo_int+hi_int)
    endfor
    norm=mean(det_eff_qty)
    det_eff_qty=det_eff_qty/norm
    logentry=logentry+string(10B)+$
      "Detector efficiencies computed using individual elastic peak fits."
endif
dwi=exp(Q_props^2*0.066)
self.det_eff_qty=det_eff_qty*dwi


ptr_free,indataPtr
; Clear any floating-point underflows that may have occurred in "gaussfit",
dummy=check_math(mask=32)
self.log=self.log+logentry
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ddcs_reduce::dcs_bgd_subtract,dark=dark
; fbdata_qty --- an nhedet-element array containing the total background count rate for each
;   detector per external pulser count; the external pulser is assumed to be running at
;   !ext_pgen Hz (see code at the end of "dcs_normalization_counts"), in which case
;   !ext_pgen*fbdata_qty is the background count rate in counts/sec.
; fbdata_err --- the corresponding error.
; fbnormalization --- an nhedet-element array containing 1s and 0s depending whether or not
;   the detector is masked.
;
;'dark' is a user supplied array that is cnts/det/sec and err
; add it to any array calculated here.
; if the array has any value it should be used.
;
logentry=''

fbdecision='' ;start with some value
o_signal = self.signal_container->get(/all)
nsignal = self.signal_container->count()
o_darkcount = self.darkcount_container->get(/all)
ndarkcount = self.darkcount_container->count()
ret = o_signal[0]->get_property(ntchan=ntchan,nhedet=nhedet)
fbdata_qty=fltarr(nhedet)
fbdata_err=fltarr(nhedet)
fbnormalization=fltarr(nhedet)


o_int = self.intermediate->get(/all)
nint = self.intermediate->count()
for i = 0 , nint-1 do begin
   ret=o_int[i]->get_property(display_name=display_name)
   if display_name eq 'SUMMED_DARKCOUNT' then index=i
endfor
test=size(index,/type)

if (test eq 0) and  (dark[0] eq 0) then fbdecision='DO NOTHING'
if (test eq 0) and  (dark[0] ne 0) then fbdecision='CONSTANT'


case fbdecision of
    'DO NOTHING': begin
        self.log = self.log +string(10B)+'No fast background defined'
        return,0
    end
    'CONSTANT':   begin
       ;should get some preferences if possible
       fbdata_qty=(double(fltarr(nhedet)+1))*float(dark[0])/3600.0/!ext_pgen  ;make DOUBLE CB
       fbdata_err=(double(fltarr(nhedet)+1))*float(dark[1])/3600.0/!ext_pgen  ; should this be !ext_pgen?
       fbnormalization=(fltarr(nhedet)+1)*1.0
       logentry=logentry+$
         "A mean background count rate has been supplied: ("+strcompress(string(dark[0]))+$
                         "+/-"+strcompress(string(dark[1]))+"/detector/h."
    end
    else: begin
       ;should not have to check this again, but...
       o_int = self.intermediate->get(/all)
       nint = self.intermediate->count()
    for i = 0 , nint-1 do begin
        ret=o_int[i]->get_property(display_name=display_name)
        if display_name eq 'SUMMED_DARKCOUNT' then index=i
    endfor
    test=size(index,/type)
    if test eq 0 then return,0

; Fast background normalization counts for fast background run

    ret = o_int[index]->get_property(tchanlook=tchanlook,$
                    nhedet=nhedet,qty=qty)
;need some values from the actual data files
    ret= o_darkcount[0]->get_property(ch_ms=ch_ms,$
                    ch_srdenom=ch_srdenom,n0ngtchan=n0ngtchan)

;    Determine fast background per pulse of the external pulser, I_B-bar(d).
    fbnormalization=*self.darkcount_fbnormalization_values_ptr
    speriod=(60000000.0/ch_ms)*ch_srdenom
    twid=transpose(tchanlook[0,0:n0ngtchan-1]*0.05)
    unitvec_nhedet=intarr(nhedet)+1
    twidPtr=ptr_new(twid#unitvec_nhedet,/no_copy)
    fbdata_qty=total(double(qty)*(*twidPtr),1)/max(fbnormalization); get some more accuracy CB
    fbdata_err=sqrt(fbdata_qty/max(fbnormalization))
    ; add the array to it too.
    if dark[0] ne 0 then fbdata_qty=fbdata_qty+(double(fltarr(nhedet)+1))*float(dark[0])/3600.0/!ext_pgen  ;
    if dark[1] ne 0 then fbdata_err=fbdata_err+(double(fltarr(nhedet)+1))*float(dark[1])/3600.0/!ext_pgen  ; should this be !ext_pgen?
    fbnormalization=fbnormalization/max(fbnormalization)
    ptr_free,twidptr
;;;;
    end
endcase

logentry=logentry+string(10B)+"Dark count backgrounds have been calculated."
if dark[0] ne 0 then begin
   logentry=logentry+"A mean background count rate has been supplied: ("$
         +strcompress(string(dark[0]))+"+/-"+strcompress(string(dark[1]))+"/detector/h."
endif
self.log=self.log + logentry + string(10B) + "Average counts/detector/hour = "+$
     strcompress(string(total(fbdata_qty)/913*3600*!ext_pgen))+" +/- "+strcompress(string(total(fbdata_err)/913*3600*!ext_pgen))
;;;;

;we should make these values global now.
self.fbnormalization=fbnormalization
self.fbdata_err=fbdata_err
self.fbdata_qty=fbdata_qty
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ddcs_reduce::dcs_sumindeticals,mask=mask
o_signal = self.signal_container->get(/all)
o_empty = self.empty_container->get(/all)
o_vanadium = self.vanadium_container->get(/all)
o_absorber = self.absorber_container->get(/all)
o_darkcount = self.darkcount_container->get(/all)

nsignal = self.signal_container->count()
nempty = self.empty_container->count()
nvanadium = self.vanadium_container->count()
nabsorber = self.absorber_container->count()
ndarkcount = self.darkcount_container->count()

;sample files first;
if nsignal ge 1 then  begin    ; this should be redundent at this point.
     ret = self->clear_int()       ; remove previous intermediates if existing
     o_new_s = self->dcs_do_sumindeticals(mask=mask,type='SIGNAL')
     ; add it to the data leaves....
     if obj_valid(o_new_s) then begin
        print,'WE HAVE A SUMMED SIGNAL OBJECT'
        self.intermediate->add,o_new_s
     endif else begin
        obj_destroy,o_new_s
        return,0
     endelse
endif

;other files next;
if nempty ge 1 then  begin
     o_new_e = self->dcs_do_sumindeticals(mask=mask,type='EMPTY')
     ; add it to the data leaves....
     if obj_valid(o_new_e) then begin
        print,'WE HAVE A SUMMED EMPTY OBJECT'
        self.intermediate->add,o_new_e
     endif else begin
        obj_destroy,o_new_e
     endelse
endif
if nabsorber ge 1 then  begin
     o_new_a = self->dcs_do_sumindeticals(mask=mask,type='ABSORBER')
     ; add it to the data leaves....
     if obj_valid(o_new_a) then begin
        print,'WE HAVE A SUMMED ABSORBER OBJECT'
        self.intermediate->add,o_new_a
     endif else begin
        obj_destroy,o_new_a
     endelse
endif
if nvanadium ge 1 then  begin
     o_new_v = self->dcs_do_sumindeticals(mask=mask,type='VANADIUM')
     ; add it to the data leaves....
     if obj_valid(o_new_v) then begin
        print,'WE HAVE A SUMMED VANADIUM OBJECT'
        self.intermediate->add,o_new_v
     endif else begin
        obj_destroy,o_new_v
     endelse
endif
if ndarkcount ge 1 then  begin
     o_new_d = self->dcs_do_sumindeticals(mask=mask,type='DARKCOUNT')
     ; add it to the data leaves....
     if obj_valid(o_new_d) then begin
        print,'WE HAVE A SUMMED DARKCOUNT OBJECT'
        self.intermediate->add,o_new_d
     endif else begin
        obj_destroy,o_new_d
     endelse
endif


return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ddcs_reduce::dcs_do_sumindeticals,mask=mask,type=type,msg = msg

case type of
'SIGNAL' : $
    begin
           o = self.signal_container->get(/all)
           numo = self.signal_container->count()
           normalization_values=*self.signal_normalization_values_ptr
           fbnormalization_values=*self.signal_fbnormalization_values_ptr
    end
'EMPTY' : $
    begin
           o = self.empty_container->get(/all)
           numo = self.empty_container->count()
           normalization_values=*self.empty_normalization_values_ptr
           fbnormalization_values=*self.empty_fbnormalization_values_ptr
    end
'VANADIUM' : $
    begin
           o = self.vanadium_container->get(/all)
           numo = self.vanadium_container->count()
           normalization_values=*self.vanadium_normalization_values_ptr
           fbnormalization_values=*self.vanadium_fbnormalization_values_ptr
    end
'ABSORBER' : $
    begin
           o = self.absorber_container->get(/all)
           numo = self.absorber_container->count()
           normalization_values=*self.absorber_normalization_values_ptr
           fbnormalization_values=*self.absorber_fbnormalization_values_ptr
    end
'DARKCOUNT' : $
    begin
           o = self.darkcount_container->get(/all)
           numo = self.darkcount_container->count()
           normalization_values=*self.darkcount_normalization_values_ptr
           fbnormalization_values=*self.darkcount_fbnormalization_values_ptr
    end
  else: return,0
endcase

logentry=''
lines=''

ret = o[0]->get_property(ntchan=ntchan,nhedet=nhedet)
unitvec_ntchan=bytarr(ntchan)+1
unitvec_nhedet=bytarr(nhedet)+1

lines="Summing identical files, applying masks, correcting bitflips,"
lines=[[lines],"dividing each channel by its width in time."]


; Quantities from raw data files (other than histodata) that are to be summed.
squantity=["badbrd","baddet","duration","nframes","ncycles",$
    "detsum","grandsum","histohigh","repeats","resets",$
    "timsum","totals"]

nsq=n_elements(squantity)

for j=0,numo-1 do begin
   ret = o[j]->get_property(maxbitflips=maxbitflips,bitflips=bitflips,$
          qty=qty,histohigh=histohigh,duration=duration,ncycles=ncycles,nframes=nframes,$
          badbrd=badbrd,baddet=baddet,detsum=detsum,grandsum=grandsum,repeats=repeats,$
          resets=resets,timsum=timsum,totals=totals,temp_setpoint=temp_setpoint,$
          command=command, comments=comments, sample_desc=sample_desc,start_date=start_date,$
          stop_date=stop_date, user=user, cduration_csum=duration_csum, he3_csum=he3_csum, $
          high_csum=high_csum, temp_control=temp_control, temp_sample=temp_sample,$
          treatment=treatment)
   ;get values

; If necessary perform bitflip correction   to data.
    for k=0,maxbitflips-1 do begin
       if (bitflips[1,k] lt 1000) then begin
         bftim=bitflips[0,k]
         bfdet=bitflips[1,k]
         bfamt=bitflips[2,k]
         lines=[[lines],"A bitflip has been corrected in dcs_sumidenticals in "+type]
         ;if (!debug) then print,j,bftim,bfdet,bfamt
         if (bfdet lt nhedet) then begin
          qty[bftim,bfdet] = qty[bftim,bfdet]-bfamt
         endif else begin
          histohigh[bftim,bfdet-nhedet] = histohigh[bftim,bfdet-nhedet]-bfamt
         endelse
       endif
    endfor
;


; Next line creates matrix with dimensions [ntchan,nhedet].
    maskj=unitvec_ntchan#mask

    if (j eq 0) then begin

;create an object in which to copy the relevent data.
newtype='SUMMED_'+type
;create a new object to contain everything:
o_new=o[j]->clone()
ret = o_new->set_property(signal_type = newtype)
;
; These quantities are to be summed.
       pdata_qty_0=ptr_new(qty*maskj)
       normalization_sum_0=(normalization_values[j]*mask)
       fbnormalization_sum_0=(fbnormalization_values[j]*mask)

;
; These quantities, from the first file, will be placed in the output sum file.
       temp_setpoint_0=temp_setpoint
       command_0=command
       comments_0=comments
       sample_desc_0=sample_desc
       start_date_0=start_date
       stop_date_0=stop_date
       user_0=user

; These quantities are to be concatenated.
       duration_csum_0=duration_csum
       he3_csum_0=he3_csum
       high_csum_0=high_csum
       temp_control_0=temp_control
       temp_sample_0=temp_sample

; These quantities are to be summed.
       for k=0,nsq-1 do begin
         com=squantity[k]+"sum="+squantity[k]
         res=execute(com)
       endfor;

;   add treatments together
    tlog=treatment
    endif

    if (j gt 0) then begin
;
; Quantities to be summed.
       *pdata_qty_0=*pdata_qty_0+qty*maskj
       normalization_sum_0=normalization_sum_0+normalization_values[j]*mask
       fbnormalization_sum_0=fbnormalization_sum_0+fbnormalization_values[j]*mask
       for k=0,nsq-1 do begin
         com=squantity[k]+"sum="+squantity[k]+$
          "sum+"+squantity[k]
         res=execute(com)
       endfor
;
; Quantities to be concatenated.
;
       void=where(duration_csum_0 ge 0,nt0)
       void=where(duration_csum ge 0,nt)
       if (nt gt 0) then begin
         duration_csum_0[nt0:nt0+nt-1]=duration_csum[0:nt-1]
         he3_csum_0[nt0:nt0+nt-1]=he3_csum[0:nt-1]
         high_csum_0[*,nt0:nt0+nt-1]=high_csum[*,0:nt-1]
       endif
;
       void=where(temp_control_0 gt 0,ntc0)
       void=where(temp_control gt 0,ntc)
       if (ntc gt 0) then $
         temp_control_0[ntc0:ntc0+ntc-1]=temp_control[0:ntc-1]
;
       void=where(temp_sample_0 gt 0,nts0)
       void=where(temp_sample gt 0,nts)
       if (nts gt 0) then $
         temp_sample_0[nts0:nts0+nts-1]=temp_sample[0:nts-1]
;
      tlog=[[tlog],treatment]

     endif
endfor

;newtype='SUMMED_'+type
;;create a new object to contain everything:
;o_new = obj_new('ddcs_data',nhedet=nhedet, ntchan =ntchan, $
;                   signal_type = newtype)

ret = o_new->set_property(qty=*pdata_qty_0, err=sqrt(*pdata_qty_0), badbrd=badbrdsum,       $
             baddet=baddetsum, duration=durationsum, nframes=nframessum, ncycles=ncyclessum,   $
             detsum=detsumsum, grandsum=grandsumsum,histohigh=histohighsum, repeats=repeatssum,$
             resets=resetssum, timsum=timsumsum, totals=totalssum, datamax=max(*pdata_qty_0),    $
             highmax=max(histohighsum),temp_setpoint=temp_setpoint_0, command=command_0,$
             comments=comments_0, sample_desc=sample_desc_0, start_date=start_date_0,   $
             stop_date=stop_date_0, user=user_0,cduration_csum=duration_csum_0,          $
             he3_csum=he3_csum_0,high_csum=high_csum_0,temp_control=temp_control_0,     $
             temp_sample=temp_sample_0, normalization_sum=normalization_sum_0 , $
             fbnormalization_sum=fbnormalization_sum_0,display_name=newtype,$
             treatment=tlog)


; Place (most of) qty, (most of) err and x and y
;   into new pointers; qty and err are first divided by the appropriate time
; channel width.
; The x array is redefined as times of flight from sample to detector.
ret = o_new->get_property(n0ngtchan=n0ngtchan,time_props=time_props)
times=(time_props)[0:n0ngtchan]
qty=*pdata_qty_0
err=sqrt(*pdata_qty_0)
twid=(times[1:n0ngtchan]-times[0:n0ngtchan-1])#unitvec_nhedet
qty=qty[0:n0ngtchan-1,*]/twid
err=err[0:n0ngtchan-1,*]/twid
x=times
;y=y still


ptr_free,pdata_qty_0
;
xlabel = 'Time-of-flight from sample to detector'
xunits = 'time:us'
xtype = 'histogram'
histlabel = 'Summed intensity'
;
ret = o_new->set_property(qty=qty,err=err,xval=x,xlabel=xlabel,ylabel=ylabel,xtype=xtype,histlabel=histlabel,xunits=xunits)

newline=''
for i=0,n_elements(lines)-1 do begin
   newline=newline+string(10B)+lines[i]
endfor
self.log=self.log+string(10B)+string(newline)
self.log=self.log+string(10B)+'Nominally similar files have been summed'

if  obj_valid(o_new) then begin
   return,o_new
endif

return,0
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ddcs_reduce::dcs_normalization,norm=norm
;

o_signal = self.signal_container->get(/all)
o_empty = self.empty_container->get(/all)
o_vanadium = self.vanadium_container->get(/all)
o_absorber = self.absorber_container->get(/all)
o_darkcount = self.darkcount_container->get(/all)

nsignal = self.signal_container->count()
nempty = self.empty_container->count()
nvanadium = self.vanadium_container->count()
nabsorber = self.absorber_container->count()
ndarkcount = self.darkcount_container->count()

nruns=(nsignal+nempty+nvanadium+nabsorber+ndarkcount)
logentry=''
;
;   Initialize a structure and replicate it for each beam monitor, each run.
; -- allcounts is the sum of counts in all time channels.
; -- peakcounts, positions, peakwidths and backgrounds are from fit of Gaussian
;      plus constant background.
; -- postheory is theoretical position of peak.
; -- posmaxmum is experimental position of peak.
;
decision=norm
bmvalues={allcounts:0.0,peakcounts:0.0,positions:0.0,peakwidths:0.0,$
    backgrounds:0.0,postheory:0.0,posmaxmum:0.0}
signal_bmvalues=replicate(bmvalues,nsignal,4)
signal_normalization_values=fltarr(nsignal)
signal_fbnormalization_values=fltarr(nsignal)
signal_other_values=fltarr(nsignal,4)
signal_duration=fltarr(nsignal)

if nempty ne 0 then begin
 empty_bmvalues=replicate(bmvalues,nempty,4)
 empty_normalization_values=fltarr(nempty)
 empty_fbnormalization_values=fltarr(nempty)
 empty_other_values=fltarr(nempty,4)
 empty_duration=fltarr(nempty)
endif

if nvanadium ne 0 then begin
 vanadium_bmvalues=replicate(bmvalues,nvanadium,4)
 vanadium_normalization_values=fltarr(nvanadium)
 vanadium_fbnormalization_values=fltarr(nvanadium)
 vanadium_other_values=fltarr(nvanadium,4)
 vanadium_duration=fltarr(nvanadium)
endif

if nabsorber ne 0 then begin
 absorber_bmvalues=replicate(bmvalues,nabsorber,4)
 absorber_normalization_values=fltarr(nabsorber)
 absorber_fbnormalization_values=fltarr(nabsorber)
 absorber_other_values=fltarr(nabsorber,4)
 absorber_duration=fltarr(nabsorber)
endif

if ndarkcount ne 0 then begin
 darkcount_bmvalues=replicate(bmvalues,ndarkcount,4)
 darkcount_normalization_values=fltarr(ndarkcount)
 darkcount_fbnormalization_values=fltarr(ndarkcount)
 darkcount_other_values=fltarr(ndarkcount,4)
 darkcount_duration=fltarr(ndarkcount)
endif

nterms=4
params=fltarr(nterms)


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

for irun=0,nsignal-1 do begin
   ret = o_signal[irun]->get_property(ch_res=ch_res,ch_wl=ch_wl,ch_dis=ch_dis,$
         ch_slots=ch_slots,startchoice=startchoice,tchanlook=tchanlook,tsdmin=tsdmin,$
         detsum=detsum,nframes=nframes,ch_ms=ch_ms,ch_srdenom=ch_srdenom,histohigh=histohigh,$
          fc_dis=fc_dis,highsource=highsource,ngtchan1=ngtchan1,nhedet=nhedet,maxbitflips=maxbitflips,$
          bitflips=bitflips)

   ret = o_signal[irun]->get_property(duration=duration)

   signal_duration[irun]=duration

; Code that is executed for the first run only STARTS HERE.
    if (irun eq 0) then begin
;
;    Define quantities needed to obtain beam monitor peak positions and intensities.
;      Master speed, master period, period at sample.
;          ret = o_signal[irun]->get_property(ch_ms=ch_ms,ch_srdenom=ch_srdenom,$
;          fc_dis=fc_dis,highsource=highsource,ngtchan1=ngtchan1)
          mspeed=ch_ms
          mperiod=60000000.0/mspeed
          speriod=mperiod*ch_srdenom
;
;      Determine "ibm" and "other":
;      "ibm" contains the locations of the beam monitor channels in "histohigh"
;      "other" contains the locations of the EXT, CH5, CH7 and AND channels.
          lines=strarr(5)
          bm_dis=[0.0,fc_dis]
          lines[0]="Beam monitor DAQ board numbers"
          lines[1]="Assumed distances from sample (mm)"
          ibm=intarr(4)
          other=intarr(4)

;    Following loop is over the 17 Lemo inputs in crate 1 (the lower crate).
          for k=0,16 do begin
              string123=strmid(highsource[k],0,3)
              if (string123 eq "FC0") then begin
                 ibm[0]=k
                 lines[0]=lines[0]+": BM0,"+strcompress(k)
              endif
              if (string123 eq "FC1") then begin
                 ibm[1]=k
                 lines[0]=lines[0]+": BM1,"+strcompress(k)
                 lines[1]=lines[1]+": BM1,"+strcompress(bm_dis[1])
              endif
              if (string123 eq "FC2") then begin
                 ibm[2]=k
                 lines[0]=lines[0]+": BM2,"+strcompress(k)
                 lines[1]=lines[1]+": BM2,"+strcompress(bm_dis[2])
              endif
              if (string123 eq "FC3") then begin
                 ibm[3]=k
                 lines[0]=lines[0]+": BM3,"+strcompress(k)
                 lines[1]=lines[1]+": BM3,"+strcompress(bm_dis[3])
              endif
              if (string123 eq "EXT") then other[0]=k
              if (string123 eq "CH5") then other[1]=k
              if (string123 eq "CH7") then other[2]=k
              if (string123 eq "AND") then other[3]=k
          endfor
;
;      (The variables "ibm" and "other" have been defined.)
;
;      Incident wavelength, incident inverse velocity, resolution mode
          lambda=ch_wl; incident wavelength
          tau0=lambda/!dcs_hom; incident inverse velocity
          resmode=ch_res; resolution mode
;
;      Chopper distances and slot positions.
          distances=ch_dis
          angles=ch_slots
;
;    3-character STARTCHOICE
          startchoice=strmid(startchoice,0,3)
;
;    Define time channel widths for high end data (in us).
          tchanwid=tchanlook[1,*]*0.05d0
;
;      Start and stop times for each channel (in us). Channel 0 starts at 0.
          tstop=reform(total(tchanwid,/cumulative))
          tstart=shift(tstop,1)
          tstart[0]=0
;
;      Number of time channels.
          ngtchan1=ngtchan1
;
;      Truncate tstart and tstop.
          tstart=tstart[0:ngtchan1-1]
          tstop=tstop[0:ngtchan1-1]
          lines[2]=" Number of time channels = "+strcompress(ngtchan1)
          lines[3]=" Time between pulses at sample = "+string(speriod)+" us."
          lines[4]=" Time channel"+strcompress(ngtchan1)+" ends at "+string(tstop[ngtchan1-1])+" us."

;
;    Add information to log.
;          for k=0,4 do logentry=logentry+string(10B)+lines[k]
;
;    Centers of time channels (for use in gaussfit)
          tmean=0.5*(tstart+tstop)
;
;    Compute ts0, which is the time that the DAQ start pulse occurs,
;        relative to the time that neutrons cross the sample axis.
;
          case startchoice of
              "AND": ts0=tsdmin
              "CH1": ts0=speriod+tau0*distances[0]-(30-angles[0,resmode-1])*(speriod/360.0)
              "CH2": ts0=speriod+tau0*distances[1]-(30-angles[1,resmode-1])*(speriod/360.0)
              "CH3": ts0=speriod+tau0*distances[2]-(30-angles[2,resmode-1])*(speriod/360.0)
              "CH4": ts0=speriod+tau0*distances[3]-(30-angles[3,resmode-1])*(speriod/360.0)
              "CH5": ts0=speriod+tau0*distances[4]-(30-angles[4,resmode-1])*(speriod/360.0)
              "CH6": ts0=speriod+tau0*distances[5]-(30-angles[5,resmode-1])*(speriod/360.0)
              "CH7": ts0=speriod+tau0*distances[6]-(30-angles[6,resmode-1])*(speriod/360.0)
              else: begin
               ;  if (!debug) then print,"Invalid STARTCHOICE"
               ;  nruns=0
                 return,0
              end
          endcase
;    If there is background (dark count) data, and if the normalization of sample
;    runs is other than to "EXT" or "Duration", ask how to normalize the background data.
;          decision=2    ; 7 would be time
          fbdecision=0
          blabel=['BM0','BM1 total','BM1 peak sum',$
                 'EXT','CH5','CH7','AND','Duration']
            fbblabel=['EXT','Duration']

   endif; Code that is executed for the first run only STOPS HERE.

; Treat each beam monitor in turn.;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    for j=0,3 do begin
;
;   Compute integrated counts over all time channels.
       allcounts=$
         ulong(total(histohigh[0:ngtchan1-1,ibm[j]],/double))
       signal_bmvalues[irun,j].allcounts=allcounts
;
;   Examine beam monitor peaks, excluding BM0 and monitors with very few total counts.
       if (j gt 0 and allcounts gt 100) then begin
         counts=histohigh[0:ngtchan1-1,ibm[j]]
;
;    If necessary perform bitflip correction to data.
         nbitflips=0
         for k=0,maxbitflips-1 do begin
          if (bitflips[1,k] eq nhedet+ibm[j]) then begin
              nbitflips=nbitflips+1
              if (nbitflips gt 2) then begin
                 result=dialog_message("Too many bitflips in dcs_normalization.",/error)
                 stop
              endif
              bitflipchannel=bitflips[0,k]
              bitflipamount=bitflips[2,k]
              result=dialog_message(["A bitflip has been corrected in dcs_normalization_counts.",$
              "Beam monitor"+strcompress(j)+", time channel"+strcompress(bitflips[0,k])+$
              ", correction"+strcompress(bitflips[2,k])])
              allcounts=allcounts-bitflipamount
              counts[bitflipchannel]=counts[bitflipchannel]-bitflipamount
          endif
         endfor
;
;      Compute time of flight from sample to beam monitor.
         t_sbm=tau0*bm_dis[j]
;
;    Hence time of arrival at beam monitor, relative to DAQ start.
;      tbm_the=dave_modulud(t_sbm-ts0,speriod) (OLD VERSION)
         tbm_the=(t_sbm-ts0) mod speriod + speriod * ((t_sbm-ts0) lt 0.0)
;
;    Compute corresponding time channel.
         tch_the=where(tstart le tbm_the and tbm_the lt tstop)
;
;      Determine channel(s) with maximum number of counts, corresponding times,
;          the mean of these times, and the corresponding time channel.
         maxcounts=max(counts)
         tch_exp=where(counts eq maxcounts)
         tbm_exp=(tstart[tch_exp]+tstop[tch_exp])/2.0
         tbm_exp=mean(tbm_exp)
         tch_exp=where(tstart le tbm_exp and tbm_exp lt tstop)
;
;    Renumber channels so that peak count is roughly in channel ngtchan1/2.
;    The number of channels shifted is tch_shift.
;      tch_shift=dave_modulus(round(-tch_exp[0]+ngtchan1/2),ngtchan1) (OLD VERSION)
         tch_shift=round(-tch_exp[0]+ngtchan1/2) mod ngtchan1 + ngtchan1 * (round(-tch_exp[0]+ngtchan1/2) lt 0.0)
         counts_rot=shift(counts,tch_shift)
         tmean_rot=shift(tmean,tch_shift)
         addition=fltarr(ngtchan1)+speriod
         if (tch_shift ne 0) then addition[0:tch_shift-1]=0
         tmean_rot=tmean_rot+addition
;
;    Check whether there have been any floating-point underflows.
;;;;         if (check_math(mask=32)) then if (!debug) then print, "NB: floating-point underflows"
;
;    Compute time shift corresponding to tch_shift.
;      tshift=tstop[dave_modulus(ngtchan1-1-tch_shift,ngtchan1)] (OLD VERSION)
         tshift=tstop[(ngtchan1-1-tch_shift) mod ngtchan1 + ngtchan1 * ((ngtchan1-1-tch_shift) lt 0.0)]
;
;    Create estimates of Gaussian fit parameters.
         estimates=[maxcounts,tshift+0.5*speriod,20.0,0.0]
;
;    Fit a Gaussian plus constant background to the rotated BM data.
         result=gaussfit(tmean_rot,counts_rot,params,estimates=estimates,nterms=nterms)
;
;    Ensure that peak position lies within period at sample.
;      params[1]=dave_modulus(params[1],speriod) (OLD VERSION)
         params[1]=(params[1] mod speriod) + speriod * (params[1] lt 0.0)
;
;    Populate the bmvalues array.
         signal_bmvalues[irun,j].peakcounts =params[0]
         signal_bmvalues[irun,j].positions  =params[1]
         signal_bmvalues[irun,j].peakwidths =params[2]
         signal_bmvalues[irun,j].backgrounds=params[3]
         signal_bmvalues[irun,j].postheory  =tbm_the
         signal_bmvalues[irun,j].posmaxmum  =tbm_exp
;
;    Clear any floating-point underflows that occurred in "gaussfit",
;    or in calculating the theoretical counts.
         dummy=check_math(mask=32)
;
;    If detailed information is to be saved, compute the fitted counts and save to file with
;    fit parameters, etc.
;;;         if (not dontsavebminfo) then begin
;;;          counts_calc=params[0]*exp(-0.5*((tmean-params[1])/params[2])^2)+params[3]
;;;
;;;;          Generate file name.
;;;          filename="_"+strcompress(irun,/remove_all)+"_0"+strcompress(j,/remove_all)+".bmd"
;;;          filename=!dcs_outp_directory+!dcs_rootname+filename
;;;          save,filename=filename,allcounts,params,tbm_the,tbm_exp,counts,counts_calc
;;;         endif


       endif; end of condition that it is not BM0 and there are at least a few counts.
    endfor; end of loop over beam monitors;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;

;   The array "other_values" contains counts for "EXT", "CH5", CH7" and "AND".
    signal_other_values[irun,0]=detsum[913+other[0]]
;
    if (startchoice eq "CH5") then $
       signal_other_values[irun,1]=nframes else $
       signal_other_values[irun,1]=detsum[913+other[1]]
;
    if (startchoice eq "CH7") then $
       signal_other_values[irun,2]=nframes else $
       signal_other_values[irun,2]=detsum[913+other[2]]
;
    if (startchoice eq "AND") then $
       signal_other_values[irun,3]=nframes else $
       signal_other_values[irun,3]=detsum[913+other[3]]
;

endfor  ; end loop over sample runs

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;-------------------------------------;

if nvanadium ge 1 then begin ; just check that there is something
for vrun=0,nvanadium-1 do begin
   ret = o_vanadium[vrun]->get_property(duration=duration,histohigh=histohigh,highsource=highsource,$
                            nframes=nframes,detsum=detsum,bitflips=bitflips)
   vanadium_duration[vrun]=duration

; Treat each beam monitor in turn.;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    for j=0,3 do begin
;
;   Compute integrated counts over all time channels.
       allcounts=$
         ulong(total(histohigh[0:ngtchan1-1,ibm[j]],/double))
       vanadium_bmvalues[vrun,j].allcounts=allcounts
;
;   Examine beam monitor peaks, excluding BM0 and monitors with very few total counts.
       if (j gt 0 and allcounts gt 100) then begin
         counts=histohigh[0:ngtchan1-1,ibm[j]]
;
;    If necessary perform bitflip correction to data.
         nbitflips=0
         for k=0,maxbitflips-1 do begin
          if (bitflips[1,k] eq nhedet+ibm[j]) then begin
              nbitflips=nbitflips+1
              if (nbitflips gt 2) then begin
                 result=dialog_message("Too many bitflips in dcs_normalization.",/error)
                 stop
              endif
              bitflipchannel=bitflips[0,k]
              bitflipamount=bitflips[2,k]
              result=dialog_message(["A bitflip has been corrected in dcs_normalization_counts.",$
              "Beam monitor"+strcompress(j)+", time channel"+strcompress(bitflips[0,k])+$
              ", correction"+strcompress(bitflips[2,k])])
              allcounts=allcounts-bitflipamount
              counts[bitflipchannel]=counts[bitflipchannel]-bitflipamount
          endif
         endfor
;
;      Compute time of flight from sample to beam monitor.
         t_sbm=tau0*bm_dis[j]
;
;    Hence time of arrival at beam monitor, relative to DAQ start.
;      tbm_the=dave_modulud(t_sbm-ts0,speriod) (OLD VERSION)
         tbm_the=(t_sbm-ts0) mod speriod + speriod * ((t_sbm-ts0) lt 0.0)
;
;    Compute corresponding time channel.
         tch_the=where(tstart le tbm_the and tbm_the lt tstop)
;
;      Determine channel(s) with maximum number of counts, corresponding times,
;          the mean of these times, and the corresponding time channel.
         maxcounts=max(counts)
         tch_exp=where(counts eq maxcounts)
         tbm_exp=(tstart[tch_exp]+tstop[tch_exp])/2.0
         tbm_exp=mean(tbm_exp)
         tch_exp=where(tstart le tbm_exp and tbm_exp lt tstop)
;
;    Renumber channels so that peak count is roughly in channel ngtchan1/2.
;    The number of channels shifted is tch_shift.
;      tch_shift=dave_modulus(round(-tch_exp[0]+ngtchan1/2),ngtchan1) (OLD VERSION)
         tch_shift=round(-tch_exp[0]+ngtchan1/2) mod ngtchan1 + ngtchan1 * (round(-tch_exp[0]+ngtchan1/2) lt 0.0)
         counts_rot=shift(counts,tch_shift)
         tmean_rot=shift(tmean,tch_shift)
         addition=fltarr(ngtchan1)+speriod
         if (tch_shift ne 0) then addition[0:tch_shift-1]=0
         tmean_rot=tmean_rot+addition
;
;    Check whether there have been any floating-point underflows.
;;;;         if (check_math(mask=32)) then if (!debug) then print, "NB: floating-point underflows"
;
;    Compute time shift corresponding to tch_shift.
;      tshift=tstop[dave_modulus(ngtchan1-1-tch_shift,ngtchan1)] (OLD VERSION)
         tshift=tstop[(ngtchan1-1-tch_shift) mod ngtchan1 + ngtchan1 * ((ngtchan1-1-tch_shift) lt 0.0)]
;
;    Create estimates of Gaussian fit parameters.
         estimates=[maxcounts,tshift+0.5*speriod,20.0,0.0]
;
;    Fit a Gaussian plus constant background to the rotated BM data.
         result=gaussfit(tmean_rot,counts_rot,params,estimates=estimates,nterms=nterms)
;
;    Ensure that peak position lies within period at sample.
;      params[1]=dave_modulus(params[1],speriod) (OLD VERSION)
         params[1]=(params[1] mod speriod) + speriod * (params[1] lt 0.0)
;
;    Populate the bmvalues array.
         vanadium_bmvalues[vrun,j].peakcounts =params[0]
         vanadium_bmvalues[vrun,j].positions  =params[1]
         vanadium_bmvalues[vrun,j].peakwidths =params[2]
         vanadium_bmvalues[vrun,j].backgrounds=params[3]
         vanadium_bmvalues[vrun,j].postheory  =tbm_the
         vanadium_bmvalues[vrun,j].posmaxmum  =tbm_exp
;
;    Clear any floating-point underflows that occurred in "gaussfit",
;    or in calculating the theoretical counts.
         dummy=check_math(mask=32)
;
;    If detailed information is to be saved, compute the fitted counts and save to file with
;    fit parameters, etc.
;;;         if (not dontsavebminfo) then begin
;;;          counts_calc=params[0]*exp(-0.5*((tmean-params[1])/params[2])^2)+params[3]
;;;
;;;;          Generate file name.
;;;          filename="_"+strcompress(irun,/remove_all)+"_0"+strcompress(j,/remove_all)+".bmd"
;;;          filename=!dcs_outp_directory+!dcs_rootname+filename
;;;          save,filename=filename,allcounts,params,tbm_the,tbm_exp,counts,counts_calc
;;;         endif


       endif; end of condition that it is not BM0 and there are at least a few counts.
    endfor; end of loop over beam monitors;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;   The array "other_values" contains counts for "EXT", "CH5", CH7" and "AND".
    vanadium_other_values[vrun,0]=detsum[913+other[0]]
    if (startchoice eq "CH5") then $
       vanadium_other_values[vrun,1]=nframes else $
       vanadium_other_values[vrun,1]=detsum[913+other[1]]
    if (startchoice eq "CH7") then $
       vanadium_other_values[vrun,2]=nframes else $
       vanadium_other_values[vrun,2]=detsum[913+other[2]]
    if (startchoice eq "AND") then $
       vanadium_other_values[vrun,3]=nframes else $
       vanadium_other_values[vrun,3]=detsum[913+other[3]]
endfor  ; end loop over vanadium runs
endif ; done checking

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


if nabsorber ge 1 then begin ; just check that there is something
for irun=0,nabsorber-1 do begin
   ret = o_absorber[irun]->get_property(duration=duration,histohigh=histohigh,highsource=highsource,$
                            nframes=nframes,detsum=detsum)
   absorber_duration[irun]=duration

; Treat each beam monitor in turn.;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    for j=0,3 do begin
;
;   Compute integrated counts over all time channels.
       allcounts=$
         ulong(total(histohigh[0:ngtchan1-1,ibm[j]],/double))
       absorber_bmvalues[irun,j].allcounts=allcounts
;
;   Examine beam monitor peaks, excluding BM0 and monitors with very few total counts.
       if (j gt 0 and allcounts gt 100) then begin
         counts=histohigh[0:ngtchan1-1,ibm[j]]
;
;    If necessary perform bitflip correction to data.
         nbitflips=0
         for k=0,maxbitflips-1 do begin
          if (bitflips[1,k] eq nhedet+ibm[j]) then begin
              nbitflips=nbitflips+1
              if (nbitflips gt 2) then begin
                 result=dialog_message("Too many bitflips in dcs_normalization.",/error)
                 stop
              endif
              bitflipchannel=bitflips[0,k]
              bitflipamount=bitflips[2,k]
              result=dialog_message(["A bitflip has been corrected in dcs_normalization_counts.",$
              "Beam monitor"+strcompress(j)+", time channel"+strcompress(bitflips[0,k])+$
              ", correction"+strcompress(bitflips[2,k])])
              allcounts=allcounts-bitflipamount
              counts[bitflipchannel]=counts[bitflipchannel]-bitflipamount
          endif
         endfor
;
;      Compute time of flight from sample to beam monitor.
         t_sbm=tau0*bm_dis[j]
;
;    Hence time of arrival at beam monitor, relative to DAQ start.
;      tbm_the=dave_modulud(t_sbm-ts0,speriod) (OLD VERSION)
         tbm_the=(t_sbm-ts0) mod speriod + speriod * ((t_sbm-ts0) lt 0.0)
;
;    Compute corresponding time channel.
         tch_the=where(tstart le tbm_the and tbm_the lt tstop)
;
;      Determine channel(s) with maximum number of counts, corresponding times,
;          the mean of these times, and the corresponding time channel.
         maxcounts=max(counts)
         tch_exp=where(counts eq maxcounts)
         tbm_exp=(tstart[tch_exp]+tstop[tch_exp])/2.0
         tbm_exp=mean(tbm_exp)
         tch_exp=where(tstart le tbm_exp and tbm_exp lt tstop)
;
;    Renumber channels so that peak count is roughly in channel ngtchan1/2.
;    The number of channels shifted is tch_shift.
;      tch_shift=dave_modulus(round(-tch_exp[0]+ngtchan1/2),ngtchan1) (OLD VERSION)
         tch_shift=round(-tch_exp[0]+ngtchan1/2) mod ngtchan1 + ngtchan1 * (round(-tch_exp[0]+ngtchan1/2) lt 0.0)
         counts_rot=shift(counts,tch_shift)
         tmean_rot=shift(tmean,tch_shift)
         addition=fltarr(ngtchan1)+speriod
         if (tch_shift ne 0) then addition[0:tch_shift-1]=0
         tmean_rot=tmean_rot+addition
;
;    Check whether there have been any floating-point underflows.
;;;;         if (check_math(mask=32)) then if (!debug) then print, "NB: floating-point underflows"
;
;    Compute time shift corresponding to tch_shift.
;      tshift=tstop[dave_modulus(ngtchan1-1-tch_shift,ngtchan1)] (OLD VERSION)
         tshift=tstop[(ngtchan1-1-tch_shift) mod ngtchan1 + ngtchan1 * ((ngtchan1-1-tch_shift) lt 0.0)]
;
;    Create estimates of Gaussian fit parameters.
         estimates=[maxcounts,tshift+0.5*speriod,20.0,0.0]
;
;    Fit a Gaussian plus constant background to the rotated BM data.
         result=gaussfit(tmean_rot,counts_rot,params,estimates=estimates,nterms=nterms)
;
;    Ensure that peak position lies within period at sample.
;      params[1]=dave_modulus(params[1],speriod) (OLD VERSION)
         params[1]=(params[1] mod speriod) + speriod * (params[1] lt 0.0)
;
;    Populate the bmvalues array.
         absorber_bmvalues[irun,j].peakcounts =params[0]
         absorber_bmvalues[irun,j].positions  =params[1]
         absorber_bmvalues[irun,j].peakwidths =params[2]
         absorber_bmvalues[irun,j].backgrounds=params[3]
         absorber_bmvalues[irun,j].postheory  =tbm_the
         absorber_bmvalues[irun,j].posmaxmum  =tbm_exp
;
;    Clear any floating-point underflows that occurred in "gaussfit",
;    or in calculating the theoretical counts.
         dummy=check_math(mask=32)
;
;    If detailed information is to be saved, compute the fitted counts and save to file with
;    fit parameters, etc.
;;;         if (not dontsavebminfo) then begin
;;;          counts_calc=params[0]*exp(-0.5*((tmean-params[1])/params[2])^2)+params[3]
;;;
;;;;          Generate file name.
;;;          filename="_"+strcompress(irun,/remove_all)+"_0"+strcompress(j,/remove_all)+".bmd"
;;;          filename=!dcs_outp_directory+!dcs_rootname+filename
;;;          save,filename=filename,allcounts,params,tbm_the,tbm_exp,counts,counts_calc
;;;         endif


       endif; end of condition that it is not BM0 and there are at least a few counts.
    endfor; end of loop over beam monitors;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;   The array "other_values" contains counts for "EXT", "CH5", CH7" and "AND".
    absorber_other_values[irun,0]=detsum[913+other[0]]
    if (startchoice eq "CH5") then $
       absorber_other_values[irun,1]=nframes else $
       absorber_other_values[irun,1]=detsum[913+other[1]]
    if (startchoice eq "CH7") then $
       absorber_other_values[irun,2]=nframes else $
       absorber_other_values[irun,2]=detsum[913+other[2]]
    if (startchoice eq "AND") then $
       absorber_other_values[irun,3]=nframes else $
       absorber_other_values[irun,3]=detsum[913+other[3]]
endfor  ; end loop over absorber runs
endif ; done checking

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


if nempty ge 1 then begin ; just check that there is something
for irun=0,nempty-1 do begin
   ret = o_empty[irun]->get_property(duration=duration,histohigh=histohigh,highsource=highsource,$
                            nframes=nframes,detsum=detsum)
   empty_duration[irun]=duration

; Treat each beam monitor in turn.;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    for j=0,3 do begin
;
;   Compute integrated counts over all time channels.
       allcounts=$
         ulong(total(histohigh[0:ngtchan1-1,ibm[j]],/double))
       empty_bmvalues[irun,j].allcounts=allcounts
;
;   Examine beam monitor peaks, excluding BM0 and monitors with very few total counts.
       if (j gt 0 and allcounts gt 100) then begin
         counts=histohigh[0:ngtchan1-1,ibm[j]]
;
;    If necessary perform bitflip correction to data.
         nbitflips=0
         for k=0,maxbitflips-1 do begin
          if (bitflips[1,k] eq nhedet+ibm[j]) then begin
              nbitflips=nbitflips+1
              if (nbitflips gt 2) then begin
                 result=dialog_message("Too many bitflips in dcs_normalization.",/error)
                 stop
              endif
              bitflipchannel=bitflips[0,k]
              bitflipamount=bitflips[2,k]
              result=dialog_message(["A bitflip has been corrected in dcs_normalization_counts.",$
              "Beam monitor"+strcompress(j)+", time channel"+strcompress(bitflips[0,k])+$
              ", correction"+strcompress(bitflips[2,k])])
              allcounts=allcounts-bitflipamount
              counts[bitflipchannel]=counts[bitflipchannel]-bitflipamount
          endif
         endfor
;
;      Compute time of flight from sample to beam monitor.
         t_sbm=tau0*bm_dis[j]
;
;    Hence time of arrival at beam monitor, relative to DAQ start.
;      tbm_the=dave_modulud(t_sbm-ts0,speriod) (OLD VERSION)
         tbm_the=(t_sbm-ts0) mod speriod + speriod * ((t_sbm-ts0) lt 0.0)
;
;    Compute corresponding time channel.
         tch_the=where(tstart le tbm_the and tbm_the lt tstop)
;
;      Determine channel(s) with maximum number of counts, corresponding times,
;          the mean of these times, and the corresponding time channel.
         maxcounts=max(counts)
         tch_exp=where(counts eq maxcounts)
         tbm_exp=(tstart[tch_exp]+tstop[tch_exp])/2.0
         tbm_exp=mean(tbm_exp)
         tch_exp=where(tstart le tbm_exp and tbm_exp lt tstop)
;
;    Renumber channels so that peak count is roughly in channel ngtchan1/2.
;    The number of channels shifted is tch_shift.
;      tch_shift=dave_modulus(round(-tch_exp[0]+ngtchan1/2),ngtchan1) (OLD VERSION)
         tch_shift=round(-tch_exp[0]+ngtchan1/2) mod ngtchan1 + ngtchan1 * (round(-tch_exp[0]+ngtchan1/2) lt 0.0)
         counts_rot=shift(counts,tch_shift)
         tmean_rot=shift(tmean,tch_shift)
         addition=fltarr(ngtchan1)+speriod
         if (tch_shift ne 0) then addition[0:tch_shift-1]=0
         tmean_rot=tmean_rot+addition
;
;    Check whether there have been any floating-point underflows.
;;;;         if (check_math(mask=32)) then if (!debug) then print, "NB: floating-point underflows"
;
;    Compute time shift corresponding to tch_shift.
;      tshift=tstop[dave_modulus(ngtchan1-1-tch_shift,ngtchan1)] (OLD VERSION)
         tshift=tstop[(ngtchan1-1-tch_shift) mod ngtchan1 + ngtchan1 * ((ngtchan1-1-tch_shift) lt 0.0)]
;
;    Create estimates of Gaussian fit parameters.
         estimates=[maxcounts,tshift+0.5*speriod,20.0,0.0]
;
;    Fit a Gaussian plus constant background to the rotated BM data.
         result=gaussfit(tmean_rot,counts_rot,params,estimates=estimates,nterms=nterms)
;
;    Ensure that peak position lies within period at sample.
;      params[1]=dave_modulus(params[1],speriod) (OLD VERSION)
         params[1]=(params[1] mod speriod) + speriod * (params[1] lt 0.0)
;
;    Populate the bmvalues array.
         empty_bmvalues[irun,j].peakcounts =params[0]
         empty_bmvalues[irun,j].positions  =params[1]
         empty_bmvalues[irun,j].peakwidths =params[2]
         empty_bmvalues[irun,j].backgrounds=params[3]
         empty_bmvalues[irun,j].postheory  =tbm_the
         empty_bmvalues[irun,j].posmaxmum  =tbm_exp
;
;    Clear any floating-point underflows that occurred in "gaussfit",
;    or in calculating the theoretical counts.
         dummy=check_math(mask=32)
;
;    If detailed information is to be saved, compute the fitted counts and save to file with
;    fit parameters, etc.
;;;         if (not dontsavebminfo) then begin
;;;          counts_calc=params[0]*exp(-0.5*((tmean-params[1])/params[2])^2)+params[3]
;;;
;;;;          Generate file name.
;;;          filename="_"+strcompress(irun,/remove_all)+"_0"+strcompress(j,/remove_all)+".bmd"
;;;          filename=!dcs_outp_directory+!dcs_rootname+filename
;;;          save,filename=filename,allcounts,params,tbm_the,tbm_exp,counts,counts_calc
;;;         endif


       endif; end of condition that it is not BM0 and there are at least a few counts.
    endfor; end of loop over beam monitors;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;   The array "other_values" contains counts for "EXT", "CH5", CH7" and "AND".
    empty_other_values[irun,0]=detsum[913+other[0]]
    if (startchoice eq "CH5") then $
       empty_other_values[irun,1]=nframes else $
       empty_other_values[irun,1]=detsum[913+other[1]]
    if (startchoice eq "CH7") then $
       empty_other_values[irun,2]=nframes else $
       empty_other_values[irun,2]=detsum[913+other[2]]
    if (startchoice eq "AND") then $
       empty_other_values[irun,3]=nframes else $
       empty_other_values[irun,3]=detsum[913+other[3]]
endfor  ; end loop over empty container runs
endif ; done checking


if ndarkcount ge 1 then begin ; just check that there is something
for irun=0,ndarkcount-1 do begin
   ret = o_darkcount[irun]->get_property(duration=duration,histohigh=histohigh,highsource=highsource,$
                            nframes=nframes,detsum=detsum)
   darkcount_duration[irun]=duration

; Treat each beam monitor in turn.;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    for j=0,3 do begin
;
;   Compute integrated counts over all time channels.
       allcounts=$
         ulong(total(histohigh[0:ngtchan1-1,ibm[j]],/double))
       darkcount_bmvalues[irun,j].allcounts=allcounts
;
;   Examine beam monitor peaks, excluding BM0 and monitors with very few total counts.
       if (j gt 0 and allcounts gt 100) then begin
         counts=histohigh[0:ngtchan1-1,ibm[j]]
;
;    If necessary perform bitflip correction to data.
         nbitflips=0
         for k=0,maxbitflips-1 do begin
          if (bitflips[1,k] eq nhedet+ibm[j]) then begin
              nbitflips=nbitflips+1
              if (nbitflips gt 2) then begin
                 result=dialog_message("Too many bitflips in dcs_normalization.",/error)
                 stop
              endif
              bitflipchannel=bitflips[0,k]
              bitflipamount=bitflips[2,k]
              result=dialog_message(["A bitflip has been corrected in dcs_normalization_counts.",$
              "Beam monitor"+strcompress(j)+", time channel"+strcompress(bitflips[0,k])+$
              ", correction"+strcompress(bitflips[2,k])])
              allcounts=allcounts-bitflipamount
              counts[bitflipchannel]=counts[bitflipchannel]-bitflipamount
          endif
         endfor
;
;      Compute time of flight from sample to beam monitor.
         t_sbm=tau0*bm_dis[j]
;
;    Hence time of arrival at beam monitor, relative to DAQ start.
;      tbm_the=dave_modulud(t_sbm-ts0,speriod) (OLD VERSION)
         tbm_the=(t_sbm-ts0) mod speriod + speriod * ((t_sbm-ts0) lt 0.0)
;
;    Compute corresponding time channel.
         tch_the=where(tstart le tbm_the and tbm_the lt tstop)
;
;      Determine channel(s) with maximum number of counts, corresponding times,
;          the mean of these times, and the corresponding time channel.
         maxcounts=max(counts)
         tch_exp=where(counts eq maxcounts)
         tbm_exp=(tstart[tch_exp]+tstop[tch_exp])/2.0
         tbm_exp=mean(tbm_exp)
         tch_exp=where(tstart le tbm_exp and tbm_exp lt tstop)
;
;    Renumber channels so that peak count is roughly in channel ngtchan1/2.
;    The number of channels shifted is tch_shift.
;      tch_shift=dave_modulus(round(-tch_exp[0]+ngtchan1/2),ngtchan1) (OLD VERSION)
         tch_shift=round(-tch_exp[0]+ngtchan1/2) mod ngtchan1 + ngtchan1 * (round(-tch_exp[0]+ngtchan1/2) lt 0.0)
         counts_rot=shift(counts,tch_shift)
         tmean_rot=shift(tmean,tch_shift)
         addition=fltarr(ngtchan1)+speriod
         if (tch_shift ne 0) then addition[0:tch_shift-1]=0
         tmean_rot=tmean_rot+addition
;
;    Check whether there have been any floating-point underflows.
;;;;         if (check_math(mask=32)) then if (!debug) then print, "NB: floating-point underflows"
;
;    Compute time shift corresponding to tch_shift.
;      tshift=tstop[dave_modulus(ngtchan1-1-tch_shift,ngtchan1)] (OLD VERSION)
         tshift=tstop[(ngtchan1-1-tch_shift) mod ngtchan1 + ngtchan1 * ((ngtchan1-1-tch_shift) lt 0.0)]
;
;    Create estimates of Gaussian fit parameters.
         estimates=[maxcounts,tshift+0.5*speriod,20.0,0.0]
;
;    Fit a Gaussian plus constant background to the rotated BM data.
         result=gaussfit(tmean_rot,counts_rot,params,estimates=estimates,nterms=nterms)
;
;    Ensure that peak position lies within period at sample.
;      params[1]=dave_modulus(params[1],speriod) (OLD VERSION)
         params[1]=(params[1] mod speriod) + speriod * (params[1] lt 0.0)
;
;    Populate the bmvalues array.
         darkcount_bmvalues[irun,j].peakcounts =params[0]
         darkcount_bmvalues[irun,j].positions  =params[1]
         darkcount_bmvalues[irun,j].peakwidths =params[2]
         darkcount_bmvalues[irun,j].backgrounds=params[3]
         darkcount_bmvalues[irun,j].postheory  =tbm_the
         darkcount_bmvalues[irun,j].posmaxmum  =tbm_exp
;
;    Clear any floating-point underflows that occurred in "gaussfit",
;    or in calculating the theoretical counts.
         dummy=check_math(mask=32)
;
;    If detailed information is to be saved, compute the fitted counts and save to file with
;    fit parameters, etc.
;;;         if (not dontsavebminfo) then begin
;;;          counts_calc=params[0]*exp(-0.5*((tmean-params[1])/params[2])^2)+params[3]
;;;
;;;;          Generate file name.
;;;          filename="_"+strcompress(irun,/remove_all)+"_0"+strcompress(j,/remove_all)+".bmd"
;;;          filename=!dcs_outp_directory+!dcs_rootname+filename
;;;          save,filename=filename,allcounts,params,tbm_the,tbm_exp,counts,counts_calc
;;;         endif


       endif; end of condition that it is not BM0 and there are at least a few counts.
    endfor; end of loop over beam monitors;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;   The array "other_values" contains counts for "EXT", "CH5", CH7" and "AND".
    darkcount_other_values[irun,0]=detsum[913+other[0]]
    if (startchoice eq "CH5") then $
       darkcount_other_values[irun,1]=nframes else $
       darkcount_other_values[irun,1]=detsum[913+other[1]]
    if (startchoice eq "CH7") then $
       darkcount_other_values[irun,2]=nframes else $
       darkcount_other_values[irun,2]=detsum[913+other[2]]
    if (startchoice eq "AND") then $
       darkcount_other_values[irun,3]=nframes else $
       darkcount_other_values[irun,3]=detsum[913+other[3]]
endfor  ; end loop over empty container runs
endif ; done checking
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;--------------------------------------;
mincts=25; Peak intensity criterion for deciding whether or not there is a peak.
;
; Add information to log.
logentry=logentry+string(10B)+"The beam monitor peaks have been analyzed."
;if (not dontsavebminfo) then logentry=logentry+string(10B)+$
;     "Beam monitor information has been saved in a series of Run_BM_ files."
;

;------ Signal files
dummy=where(signal_bmvalues.(1) gt mincts,count)
if (count eq 0) then $
    logentry=logentry+string(10B)+"No peaks have peak intensity >"+strcompress(mincts)+"."
;
if (count gt 0) then begin
    logentry=logentry+string(10B)+"The following lines contain the results of Gaussian fits: "+$
       string(10B)+"Results are only printed for peaks with peak intensity >"+strcompress(mincts)+"."+$
       string(10B)+"Filename, BM_ID, peak intensity, position (us), "+$
       "width (us), background, and ''integrated intensity''."



    for ibm=0,3 do begin
       for irun=0,nsignal-1 do begin
       ret = o_signal[irun]->get_property(display_name=display_name)
         if (signal_bmvalues[irun,ibm].(1) gt 10) then begin
          pvals=string(signal_bmvalues[irun,ibm].(1),signal_bmvalues[irun,ibm].(2),$
              signal_bmvalues[irun,ibm].(3),signal_bmvalues[irun,ibm].(4),$
              signal_bmvalues[irun,ibm].(1)*signal_bmvalues[irun,ibm].(3)*sqrt(2.*!pi),$
              format='(f12.0,f9.1,f7.1,f7.1,f15.1)')
              logentry=logentry+string(10B)+$
                 "File "+display_name+": BM"+strcompress(ibm,/remove_all)+pvals
         endif
       endfor
    endfor
endif
;

;make a 1d array
arr=fltarr(nhedet)+1

case decision of ; (decision 2 is BM1 "integrated intensities")
    0: *self.signal_normalization_values_ptr=arr*signal_bmvalues[*,0].allcounts; BM0 total counts
    1: *self.signal_normalization_values_ptr=arr*signal_bmvalues[*,1].allcounts; BM1 total counts
    2: *self.signal_normalization_values_ptr=arr*signal_bmvalues[*,1].peakcounts*signal_bmvalues[*,1].peakwidths
    3: *self.signal_normalization_values_ptr=arr*signal_other_values[*,0]; EXT counts
    4: *self.signal_normalization_values_ptr=arr*signal_other_values[*,1]; CH5 counts
    5: *self.signal_normalization_values_ptr=arr*signal_other_values[*,2]; CH7 counts
    6: *self.signal_normalization_values_ptr=arr*signal_other_values[*,3]; AND counts
    7: *self.signal_normalization_values_ptr=arr*signal_duration; run duration
    else:
endcase
;
;------vanadium files
if nvanadium ge 1 then begin
 dummy=where(vanadium_bmvalues.(1) gt mincts,count)
 if (count eq 0) then $
    logentry=logentry+string(10B)+"No peaks have peak intensity >"+strcompress(mincts)+"."
;
 if (count gt 0) then begin
    logentry=logentry+string(10B)+"The following lines contain the results of Gaussian fits: "+$
       string(10B)+"Results are only printed for peaks with peak intensity >"+strcompress(mincts)+"."+$
       string(10B)+"Filename, BM_ID, peak intensity, position (us), "+$
       "width (us), background, and ''integrated intensity''."



    for ibm=0,3 do begin
       for irun=0,nvanadium-1 do begin
       ret = o_vanadium[irun]->get_property(display_name=display_name)
         if (vanadium_bmvalues[irun,ibm].(1) gt 10) then begin
          pvals=string(vanadium_bmvalues[irun,ibm].(1),vanadium_bmvalues[irun,ibm].(2),$
              vanadium_bmvalues[irun,ibm].(3),vanadium_bmvalues[irun,ibm].(4),$
              vanadium_bmvalues[irun,ibm].(1)*vanadium_bmvalues[irun,ibm].(3)*sqrt(2.*!pi),$
              format='(f12.0,f9.1,f7.1,f7.1,f15.1)')
              logentry=logentry+string(10B)+$
                 "File "+display_name+": BM"+strcompress(ibm,/remove_all)+pvals
         endif
       endfor
    endfor
 endif
;
 case decision of ; (decision 2 is BM1 "integrated intensities")
    0: *self.vanadium_normalization_values_ptr=arr*vanadium_bmvalues[*,0].allcounts; BM0 total counts
    1: *self.vanadium_normalization_values_ptr=arr*vanadium_bmvalues[*,1].allcounts; BM1 total counts
    2: *self.vanadium_normalization_values_ptr=arr*vanadium_bmvalues[*,1].peakcounts*vanadium_bmvalues[*,1].peakwidths
    3: *self.vanadium_normalization_values_ptr=arr*vanadium_other_values[*,0]; EXT counts
    4: *self.vanadium_normalization_values_ptr=arr*vanadium_other_values[*,1]; CH5 counts
    5: *self.vanadium_normalization_values_ptr=arr*vanadium_other_values[*,2]; CH7 counts
    6: *self.vanadium_normalization_values_ptr=arr*vanadium_other_values[*,3]; AND counts
    7: *self.vanadium_normalization_values_ptr=arr*vanadium_duration; run duration
    else:
 endcase
endif
;
;------empty files
if nempty ge 1 then begin
 dummy=where(empty_bmvalues.(1) gt mincts,count)
 if (count eq 0) then $
    logentry=logentry+string(10B)+"No peaks have peak intensity >"+strcompress(mincts)+"."
;
 if (count gt 0) then begin
    logentry=logentry+string(10B)+"The following lines contain the results of Gaussian fits: "+$
       string(10B)+"Results are only printed for peaks with peak intensity >"+strcompress(mincts)+"."+$
       string(10B)+"Filename, BM_ID, peak intensity, position (us), "+$
       "width (us), background, and ''integrated intensity''."

    for ibm=0,3 do begin
       for irun=0,nempty-1 do begin
       ret = o_empty[irun]->get_property(display_name=display_name)
         if (empty_bmvalues[irun,ibm].(1) gt 10) then begin
          pvals=string(empty_bmvalues[irun,ibm].(1),empty_bmvalues[irun,ibm].(2),$
              empty_bmvalues[irun,ibm].(3),empty_bmvalues[irun,ibm].(4),$
              empty_bmvalues[irun,ibm].(1)*empty_bmvalues[irun,ibm].(3)*sqrt(2.*!pi),$
              format='(f12.0,f9.1,f7.1,f7.1,f15.1)')
              logentry=logentry+string(10B)+$
                 "File "+display_name+": BM"+strcompress(ibm,/remove_all)+pvals
         endif
       endfor
    endfor
 endif
;
 case decision of ; (decision 2 is BM1 "integrated intensities")
    0: *self.empty_normalization_values_ptr=arr*empty_bmvalues[*,0].allcounts; BM0 total counts
    1: *self.empty_normalization_values_ptr=arr*empty_bmvalues[*,1].allcounts; BM1 total counts
    2: *self.empty_normalization_values_ptr=arr*empty_bmvalues[*,1].peakcounts*empty_bmvalues[*,1].peakwidths
    3: *self.empty_normalization_values_ptr=arr*empty_other_values[*,0]; EXT counts
    4: *self.empty_normalization_values_ptr=arr*empty_other_values[*,1]; CH5 counts
    5: *self.empty_normalization_values_ptr=arr*empty_other_values[*,2]; CH7 counts
    6: *self.empty_normalization_values_ptr=arr*empty_other_values[*,3]; AND counts
    7: *self.empty_normalization_values_ptr=arr*empty_duration; run duration
    else:
 endcase
endif

;------absorber files
if nabsorber ge 1 then begin
 dummy=where(absorber_bmvalues.(1) gt mincts,count)
 if (count eq 0) then $
    logentry=logentry+string(10B)+"No peaks have peak intensity >"+strcompress(mincts)+"."
;
 if (count gt 0) then begin
    logentry=logentry+string(10B)+"The following lines contain the results of Gaussian fits: "+$
       string(10B)+"Results are only printed for peaks with peak intensity >"+strcompress(mincts)+"."+$
       string(10B)+"Filename, BM_ID, peak intensity, position (us), "+$
       "width (us), background, and ''integrated intensity''."

    for ibm=0,3 do begin
       for irun=0,nabsorber-1 do begin
       ret = o_absorber[irun]->get_property(display_name=display_name)
         if (absorber_bmvalues[irun,ibm].(1) gt 10) then begin
          pvals=string(absorber_bmvalues[irun,ibm].(1),absorber_bmvalues[irun,ibm].(2),$
              absorber_bmvalues[irun,ibm].(3),absorber_bmvalues[irun,ibm].(4),$
              absorber_bmvalues[irun,ibm].(1)*absorber_bmvalues[irun,ibm].(3)*sqrt(2.*!pi),$
              format='(f12.0,f9.1,f7.1,f7.1,f15.1)')
              logentry=logentry+string(10B)+$
                 "File "+display_name+": BM"+strcompress(ibm,/remove_all)+pvals
         endif
       endfor
    endfor
 endif

 case decision of ; (decision 2 is BM1 "integrated intensities")
    0: *self.absorber_normalization_values_ptr=arr*absorber_bmvalues[*,0].allcounts; BM0 total counts
    1: *self.absorber_normalization_values_ptr=arr*absorber_bmvalues[*,1].allcounts; BM1 total counts
    2: *self.absorber_normalization_values_ptr=arr*absorber_bmvalues[*,1].peakcounts*absorber_bmvalues[*,1].peakwidths
    3: *self.absorber_normalization_values_ptr=arr*absorber_other_values[*,0]; EXT counts
    4: *self.absorber_normalization_values_ptr=arr*absorber_other_values[*,1]; CH5 counts
    5: *self.absorber_normalization_values_ptr=arr*absorber_other_values[*,2]; CH7 counts
    6: *self.absorber_normalization_values_ptr=arr*absorber_other_values[*,3]; AND counts
    7: *self.absorber_normalization_values_ptr=arr*absorber_duration; run duration
    else:
 endcase
endif

;------darkcount files
if ndarkcount ge 1 then begin
 dummy=where(darkcount_bmvalues.(1) gt mincts,count)
 if (count eq 0) then $
    logentry=logentry+string(10B)+"No peaks have peak intensity >"+strcompress(mincts)+"."
;
 if (count gt 0) then begin
    logentry=logentry+string(10B)+"The following lines contain the results of Gaussian fits: "+$
       string(10B)+"Results are only printed for peaks with peak intensity >"+strcompress(mincts)+"."+$
       string(10B)+"Filename, BM_ID, peak intensity, position (us), "+$
       "width (us), background, and ''integrated intensity''."

    for ibm=0,3 do begin
       for irun=0,ndarkcount-1 do begin
       ret = o_darkcount[irun]->get_property(display_name=display_name)
         if (darkcount_bmvalues[irun,ibm].(1) gt 10) then begin
          pvals=string(darkcount_bmvalues[irun,ibm].(1),darkcount_bmvalues[irun,ibm].(2),$
              darkcount_bmvalues[irun,ibm].(3),darkcount_bmvalues[irun,ibm].(4),$
              darkcount_bmvalues[irun,ibm].(1)*darkcount_bmvalues[irun,ibm].(3)*sqrt(2.*!pi),$
              format='(f12.0,f9.1,f7.1,f7.1,f15.1)')
              logentry=logentry+string(10B)+$
                 "File "+display_name+": BM"+strcompress(ibm,/remove_all)+pvals
         endif
       endfor
    endfor
 endif
;
;
 case decision of ; (decision 2 is BM1 "integrated intensities")
    0: *self.darkcount_normalization_values_ptr=arr*darkcount_bmvalues[*,0].allcounts; BM0 total counts
    1: *self.darkcount_normalization_values_ptr=arr*darkcount_bmvalues[*,1].allcounts; BM1 total counts
    2: *self.darkcount_normalization_values_ptr=arr*darkcount_bmvalues[*,1].peakcounts*darkcount_bmvalues[*,1].peakwidths
    3: *self.darkcount_normalization_values_ptr=arr*darkcount_other_values[*,0]; EXT counts
    4: *self.darkcount_normalization_values_ptr=arr*darkcount_other_values[*,1]; CH5 counts
    5: *self.darkcount_normalization_values_ptr=arr*darkcount_other_values[*,2]; CH7 counts
    6: *self.darkcount_normalization_values_ptr=arr*darkcount_other_values[*,3]; AND counts
    7: *self.darkcount_normalization_values_ptr=arr*darkcount_duration; run duration
    else:
 endcase
endif


case fbdecision of
    0: begin
       *self.signal_fbnormalization_values_ptr=arr#signal_other_values[*,0]; EXT counts
       if nempty gt 0 then  *self.empty_fbnormalization_values_ptr=arr#empty_other_values[*,0]; EXT counts
       if nvanadium gt 0 then *self.vanadium_fbnormalization_values_ptr=arr#vanadium_other_values[*,0]; EXT counts
       if nabsorber gt 0 then *self.absorber_fbnormalization_values_ptr=arr#absorber_other_values[*,0]; EXT counts
       if ndarkcount  gt 0 then *self.darkcount_fbnormalization_values_ptr=arr#darkcount_other_values[*,0]; EXT counts
       end
    1: begin
       *self.signal_fbnormalization_values_ptr=arr*signal_duration[*]; run duration
       if nempty gt 0 then  *self.empty_fbnormalization_values_ptr=arr*empty_duration[*]; EXT counts
       if nvanadium gt 0 then *self.vanadium_fbnormalization_values_ptr=arr*vanadium_other_duration[*]; EXT counts
       if nabsorber gt 0 then *self.absorber_fbnormalization_values_ptr=arr*absorber_other_duration[*]; EXT counts
       if ndarkcount  gt 0 then *self.darkcount_fbnormalization_values_ptr=arr*darkcount_duration[*]; EXT counts

       end
    else:
endcase
;
logentry=logentry+string(10B)+$
    "Normalization other than to background will be performed using "+$
    blabel[decision]+" counts."
;
if (ndarkcount ne 0) then begin
    logentry=logentry+string(10B)+"Normalization to background will be performed using "+$
       fbblabel[fbdecision]+" counts."
       sample_ratio_ext=*self.signal_fbnormalization_values_ptr/signal_duration
       ;should not need these.... just check the sample values is enough
;       if nvanadium gt 0 then vanadium_ratio_ext=vanadium_fbnormalization_values/vanadium_duration
;       if nempty gt 0 then empty_ratio_ext=empty_fbnormalization_values/empty_duration
;       if nabsorber gt 0 then absorber_ratio_ext=absorber_fbnormalization_values/absorber_duration
;       if ndarkcount  gt 0 then darkcount_ratio_ext=darkcount_fbnormalization_values/darkcount_duration

;    if (!debug) then print,ratio_ext
    extpgmin=!ext_pgen*0.98
    extpgmax=!ext_pgen*1.02
    for irun=0,nsignal-1 do begin
       if (sample_ratio_ext[irun] lt extpgmin or sample_ratio_ext[irun] gt extpgmax) then begin
         res=dialog_message("It appears that the external"+strcompress(!ext_pgen)+$
          " Hz pulser was not in use.  The ''fast background'' run MUST be removed from the tree")
;;;; need to remove fast background from the calling tree!;;;;;;;;;;;XXX
         ;nfiles_selected[4]=0
         ;files_selected=files_selected[0:nruns-2]
         return,0
       endif
    endfor
endif

self.log=self.log+string(10B)+logentry

return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ddcs_reduce::check_sameness,lines=lines
; These are the quantities that are required to be identical.
; Quantities removed from this list: "motor_pos","coll_osc"...
;
quantity=["ch_ms","ch_res","ch_srdenom","ch_srmode","ch_wl",$
    "coll_amp","coll_mean","det_dis","tsdmin",$
    "ch_bid","ch_delay","ch_dis","ch_input","ch_phase","ch_slots",$
    "fc_dis","tchanlook","startchoice","highsource"]
ntests=n_elements(quantity)

; number of files other than background files.
o_signal = self.signal_container->get(/all)
o_empty = self.empty_container->get(/all)
o_vanadium = self.vanadium_container->get(/all)
o_absorber = self.absorber_container->get(/all)
nsignal = self.signal_container->count()
nempty = self.empty_container->count()
nvanadium = self.vanadium_container->count()
nabsorber = self.absorber_container->count()

lines=strarr(nsignal+nempty+nvanadium+nabsorber)
;
nline=0

notsame=fltarr(ntests,nsignal)
for i = 0,nsignal-1 do begin
   ret = o_signal[i]->get_property(display_name=display_name,$
         ch_srdenom=ch_srdenom,ch_srmode=ch_srmode,ch_input=ch_input)
   if (ch_srdenom eq 2) then begin
      ch_srmode = 2
      ch_input[3] = 2.000
   endif
    if (i eq 0) then begin ; All files are compared with file 1.
       for k=0,ntests-1 do begin
         com1="ret=o_signal[i]->get_property("+quantity[k]+"="+quantity[k]+")"
         res=execute(com1)
         com=quantity[k]+"_0="+quantity[k]
         res=execute(com)
       endfor
    endif else begin
       for k=0,ntests-1 do begin
         com1="ret=o_signal[i]->get_property("+quantity[k]+"="+quantity[k]+")"
         res=execute(com1)
         com="tst=where("+quantity[k]+" ne "+quantity[k]+"_0,count)"
         res=execute(com)
         com="notsame["+strcompress(k,/remove_all)+",[i]]=count"
         res=execute(com)
       endfor
       if (total(notsame[*,i],1) gt 0) then begin
         ;nline=nline+2
         lines[nline]=display_name+" o_signal["
         lines[nline]=lines[nline]+strcompress(string(i),/remove_all)+"]  /"
         for k=0,ntests-1 do begin
          if (notsame[k,i] gt 0) then lines[nline]=lines[nline]+quantity[k]+"/
         endfor
         nline=nline+1
       endif
    endelse
endfor

if nempty ge 1 then begin
   notsame=fltarr(ntests,nempty)
   for i = 0,nempty-1 do begin
      ret = o_empty[i]->get_property(ch_srdenom=ch_srdenom,ch_srmode=ch_srmode,ch_input=ch_input)
      if (ch_srdenom eq 2) then begin
         ch_srmode = 2
         ch_input[3] = 2.000
      endif
      for k=0,ntests-1 do begin
         com1="ret=o_empty[i]->get_property("+quantity[k]+"="+quantity[k]+")"
         res=execute(com1)

    ; stop those problems with samples differing from empty
       if quantity[k] eq "ch_phase" then begin
         chtest=0
         for ii=0,n_elements(ch_phase)-1 do begin
           if ch_phase[ii] ne ch_phase_0[ii] then chtest=1
         endfor
         if chtest eq 1 then begin
            stringcc="Your EMPTY 'ch_phase' does not match the sample. Do you want to ignore these errors?"
            ctest= dialog_message(stringcc, /question,title='Question')
            if ctest eq "Yes" then  ch_phase = ch_phase_0
         endif
       endif


         com="tst=where("+quantity[k]+" ne "+quantity[k]+"_0,count)"
         res=execute(com)
         com="notsame["+strcompress(k,/remove_all)+",[i]]=count"
         res=execute(com)
      endfor
      if (total(notsame[*,i],1) gt 0) then begin
         ;nline=nline+2
         lines[nline]=display_name+" o_empty["
         lines[nline]=lines[nline]+strcompress(string(i),/remove_all)+"]  /"
         for k=0,ntests-1 do begin
          if (notsame[k,i] gt 0) then lines[nline]=lines[nline]+quantity[k]+"/
         endfor
         nline=nline+1
      endif
   endfor
endif

if nvanadium ge 1 then begin
   notsame=fltarr(ntests,nvanadium)
   for i = 0,nvanadium-1 do begin
      ret = o_vanadium[i]->get_property(ch_srdenom=ch_srdenom,ch_srmode=ch_srmode,ch_input=ch_input)
      if (ch_srdenom eq 2) then begin
         ch_srmode = 2
         ch_input[3] = 2.000
      endif


      for k=0,ntests-1 do begin
         com1="ret=o_vanadium[i]->get_property("+quantity[k]+"="+quantity[k]+")"
         res=execute(com1)

    ; stop those problems with samples differing from vanadium
       if quantity[k] eq "coll_amp" then begin
         if coll_amp ne coll_amp_0 then begin
            stringcc="Your vanadium 'coll_amp' does not match the sample. Do you want to ignore this error?"
            ctest= dialog_message(stringcc, /question,title='Question')
            if ctest eq "Yes" then  coll_amp = coll_amp_0
         endif
       endif

       if quantity[k] eq "highsource" then begin
         htst=where(highsource ne highsource_0,counth)
         if counth ne 0 then begin
          stringhh="Your vanadium 'highsource' header does not match the sample. Do you want to ignore this error?"
            htest= dialog_message(stringhh,/question,title='Question')
            if htest eq "Yes" then  highsource = highsource_0
         endif
       endif

         com="tst=where("+quantity[k]+" ne "+quantity[k]+"_0,count)"
         res=execute(com)
         com="notsame["+strcompress(k,/remove_all)+",[i]]=count"
         res=execute(com)
      endfor
      if (total(notsame[*,i],1) gt 0) then begin
         ;nline=nline+2
         lines[nline]=display_name+" o_vanadium["
         lines[nline]=lines[nline]+strcompress(string(i),/remove_all)+"]  /"
         for k=0,ntests-1 do begin
          if (notsame[k,i] gt 0) then lines[nline]=lines[nline]+quantity[k]+"/
         endfor
         nline=nline+1
      endif
   endfor
endif

if nabsorber ge 1 then begin
   notsame=fltarr(ntests,nabsorber)
   for i = 0,nabsorber-1 do begin
      ret = o_absorber[i]->get_property(ch_srdenom=ch_srdenom,ch_srmode=ch_srmode,ch_input=ch_input)
      if (ch_srdenom eq 2) then begin
         ch_srmode = 2
         ch_input[3] = 2.000
      endif
      for k=0,ntests-1 do begin
         com1="ret=o_absorber[i]->get_property("+quantity[k]+"="+quantity[k]+")"
         res=execute(com1)
         com="tst=where("+quantity[k]+" ne "+quantity[k]+"_0,count)"
         res=execute(com)
         com="notsame["+strcompress(k,/remove_all)+",[i]]=count"
         res=execute(com)
      endfor
      if (total(notsame[*,i],1) gt 0) then begin
         ;nline=nline+2
         lines[nline]=display_name+" o_absorber["
         lines[nline]=lines[nline]+strcompress(string(i),/remove_all)+"]  /"
         for k=0,ntests-1 do begin
          if (notsame[k,i] gt 0) then lines[nline]=lines[nline]+quantity[k]+"/
         endfor
         nline=nline+1
      endif
   endfor
endif

;now what to do?
;return lines and process the errors?

test=where(lines ne '',count)
if count ge 1 then begin
   newlines=strarr(n_elements(test))
   for i=0,count-1 do begin
      newlines[i]=lines[test[i]]
   endfor
   lines=newlines
endif else begin
   lines=''
endelse

end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ddcs_reduce::add,oarray
for i = 0,n_elements(oarray)-1 do begin
   if obj_class(oarray[i]) eq 'DDCS_DATA' and obj_valid(oarray[i]) 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]
   'VANADIUM':   self.vanadium_container->add,oarray[i]
   'DARKCOUNT':   self.darkcount_container->add,oarray[i]
   'ABSORBER':   self.absorber_container->add,oarray[i]
   else:
   endcase
endfor
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ddcs_reduce::check_contents
; This function goes through the contents of the SIGNAL container and the EMPTY
; container etc.and determines if each of the contents' signal_types matches the
; container type (EMPTY, VANADIUM, DARKCOUNT or SIGNAL).
oall = self->get(/all)
o_signal = self.signal_container->get(/all)
o_empty = self.empty_container->get(/all)
o_vanadium = self.vanadium_container->get(/all)
o_darkcount = self.darkcount_container->get(/all)
o_absorber = self.absorber_container->get(/all)
nsignal = self.signal_container->count()
nempty = self.empty_container->count()
nvanadium = self.vanadium_container->count()
ndarkcount = self.darkcount_container->count()
nabsorber = self.absorber_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
   if signal_type eq 'VANADIUM' then begin
      self.signal_container->remove,o_signal[i]
      self.vanadium_container->add,o_signal[i]
   endif
   if signal_type eq 'DARKCOUNT' then begin
      self.signal_container->remove,o_signal[i]
      self.darkcount_container->add,o_signal[i]
   endif
   if signal_type eq 'ABSORBER' then begin
      self.signal_container->remove,o_signal[i]
      self.absorber_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
   if signal_type eq 'VANADIUM' then begin
      self.empty_container->remove,o_empty[i]
      self.vanadium_container->add,o_empty[i]
   endif
   if signal_type eq 'DARKCOUNT' then begin
      self.empty_container->remove,o_empty[i]
      self.darkcount_container->add,o_empty[i]
   endif
   if signal_type eq 'ABSORBER' then begin
      self.empty_container->remove,o_empty[i]
      self.absorber_container->add,o_empty[i]
   endif
endfor
for i = 0,nvanadium-1 do begin
   ret = o_vanadium[i]->get_property(signal_type = signal_type)
   if signal_type eq 'SIGNAL' then begin
      self.vanadium_container->remove,o_vanadium[i]
      self.signal_container->add,o_vanadium[i]
   endif
   if signal_type eq 'EMPTY' then begin
      self.vanadium_container->remove,o_vanadium[i]
      self.empty_container->add,o_vanadium[i]
   endif
   if signal_type eq 'DARKCOUNT' then begin
      self.vanadium_container->remove,o_vanadium[i]
      self.darkcount_container->add,o_vanadium[i]
   endif
   if signal_type eq 'ABSORBER' then begin
      self.vanadium_container->remove,o_vanadium[i]
      self.absorber_container->add,o_vanadium[i]
   endif

endfor
for i = 0,nabsorber-1 do begin
   ret = o_absorber[i]->get_property(signal_type = signal_type)
   if signal_type eq 'SIGNAL' then begin
      self.absorber_container->remove,o_absorber[i]
      self.signal_container->add,o_absorber[i]
   endif
   if signal_type eq 'EMPTY' then begin
      self.absorber_container->remove,o_absorber[i]
      self.empty_container->add,o_absorber[i]
   endif
   if signal_type eq 'DARKCOUNT' then begin
      self.absorber_container->remove,o_absorber[i]
      self.darkcount_container->add,o_absorber[i]
   endif
   if signal_type eq 'VANADIUM' then begin
      self.absorber_container->remove,o_absorber[i]
      self.vanadium_container->add,o_absorber[i]
   endif

endfor

for i = 0,ndarkcount-1 do begin
   ret = o_darkcount[i]->get_property(signal_type = signal_type)
   if signal_type eq 'SIGNAL' then begin
      self.darkcount_container->remove,o_darkcount[i]
      self.signal_container->add,o_darkcount[i]
   endif
   if signal_type eq 'VANADIUM' then begin
      self.darkcount_container->remove,o_darkcount[i]
      self.vanadium_container->add,o_darkcount[i]
   endif
   if signal_type eq 'EMPTY' then begin
      self.darkcount_container->remove,o_darkcount[i]
      self.empty_container->add,o_darkcount[i]
   endif
   if signal_type eq 'ABSORBER' then begin
      self.darkcount_container->remove,o_darkcount[i]
      self.absorber_container->add,o_darkcount[i]
   endif

endfor

return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ddcs_reduce::init,obj_array,group_leader=group_leader
if not self->idl_container::init() then return,0
if keyword_set(group_leader) then begin
   self.group_leader=group_leader
endif else begin
   self.group_leader=0L
endelse
self.signal_container = obj_new('IDL_CONTAINER')
self.empty_container = obj_new('IDL_CONTAINER')
self.vanadium_container = obj_new('IDL_CONTAINER')
self.darkcount_container = obj_new('IDL_CONTAINER')
self.absorber_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]
      'DARKCOUNT':self.darkcount_container->add,obj_array[i]
      'VANADIUM': self.vanadium_container->add,obj_array[i]
      'ABSORBER': self.absorber_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.corrected_container = obj_new('IDL_CONTAINER')

;stores beam monitor counts and durations for each run
self.signal_normalization_values_ptr=ptr_new(/allocate_heap)
self.vanadium_normalization_values_ptr=ptr_new(/allocate_heap)
self.empty_normalization_values_ptr=ptr_new(/allocate_heap)
self.absorber_normalization_values_ptr=ptr_new(/allocate_heap)
self.darkcount_normalization_values_ptr=ptr_new(/allocate_heap)
self.signal_fbnormalization_values_ptr=ptr_new(/allocate_heap)
self.vanadium_fbnormalization_values_ptr=ptr_new(/allocate_heap)
self.empty_fbnormalization_values_ptr=ptr_new(/allocate_heap)
self.absorber_fbnormalization_values_ptr=ptr_new(/allocate_heap)
self.darkcount_fbnormalization_values_ptr=ptr_new(/allocate_heap)

nhedet=913
self.det_eff_qty=fltarr(nhedet)+1
self.det_eff_err=fltarr(nhedet)
self.fbnormalization=fltarr(nhedet)
self.fbdata_err=fltarr(nhedet)
self.fbdata_qty=fltarr(nhedet)+1

;;;
self.ngroups=-1
self.grouping_ptr=ptr_new(/allocate_heap)
self.angle_props_ptr=ptr_new(/allocate_heap)
self.Q_props_ptr=ptr_new(/allocate_heap)
self.elasticqvals_ptr=ptr_new(/allocate_heap)
self.qflagPtr=ptr_new(/allocate_heap)
self.doQgrouping=-1
self.attenuation_ptr=ptr_new(/allocate_heap)
self.angles=fltarr(nhedet)
self.q_values=fltarr(nhedet)
self.content=strarr(nhedet+1)
;;;
self.log=''
self.interp = 'LINEAR'  ; (default)
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro ddcs_reduce__define
define = {  ddcs_reduce,                  $
            group_leader:0L,              $
            interp:'',                    $
            intermediate:obj_new(),       $
            signal_container:obj_new(),   $
            empty_container:obj_new(),    $
            absorber_container:obj_new(),    $
            vanadium_container:obj_new(), $
            darkcount_container:obj_new(), $
            reduced_container:obj_new(),  $
            zip_container:obj_new(),      $
            corrected_container:obj_new(),  $
            signal_normalization_values_ptr:ptr_new(),  $
            vanadium_normalization_values_ptr:ptr_new(),  $
            empty_normalization_values_ptr:ptr_new(),  $
            absorber_normalization_values_ptr:ptr_new(),  $
            darkcount_normalization_values_ptr:ptr_new(),  $
            signal_fbnormalization_values_ptr:ptr_new(),  $
            vanadium_fbnormalization_values_ptr:ptr_new(),  $
            empty_fbnormalization_values_ptr:ptr_new(),  $
            absorber_fbnormalization_values_ptr:ptr_new(),  $
            darkcount_fbnormalization_values_ptr:ptr_new(),  $
            det_eff_qty:fltarr(913),      $
            det_eff_err:fltarr(913),      $
            fbnormalization:fltarr(913),  $
            fbdata_err:fltarr(913),       $
            fbdata_qty:fltarr(913),       $
            content:strarr(914),          $
            angles:fltarr(913),           $
            q_values:fltarr(913),         $
            ngroups:0,                    $
            grouping_ptr:ptr_new(),       $
            angle_props_ptr:ptr_new(),    $
            elasticqvals_ptr:ptr_new(),   $
            Q_props_ptr:ptr_new(),        $
            qflagPtr:ptr_new(),           $
            doQgrouping: -1       ,       $
            attenuation_ptr:ptr_new(),    $
            log:'',                       $
inherits IDL_CONTAINER  }
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;