; $Id: dm_load_dcsdata.pro,v 1.45 2008/02/14 21:00:48 ymqiu Exp $
;#######################################################################
;
; NAME:
;  dm_load_dcsdata
;
; PURPOSE:
;  This program loads dcs file(s) for dcs_mslice and dcs_cryst_align program.
;
; CATEGORY:
;  dcs_mslice
;
; AUTHOR:
;  Yiming Qiu
;  NIST Center for Neutron Research
;  100 Bureau Drive, Gaithersburg, MD 20899-6102
;  United States
;  yiming.qiu@nist.gov
;  February, 2008
;
; 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.
;
;#######################################################################

; Calls:
;   dcs_read_binaryOctavefile   :   the program which actually reads in a dcs data
;   dm_center_kid               :   widget centering
;   dm_define_pointer           :   define a dave pointer
;   dm_ispower2                 :   determine if  a number is power of two
;   dm_locate_datastrptr        :   locate datastrPtr from the dave pointer
;   dm_to_string                :   string conversion
; Parameters:
;   open_path:  file path,string
;   open_file:  file name(s),string or string array
; Keywords:
;  input:
;   avgsum:     0:average 1:sum, effective when diffuse keyword is set,default is 1
;   diffetyp:   0:single energy  1:multiple energy, only needed when diffuse keyword is set
;   diffuse:    if set, the output is intended for dcs_mslice diffuse scattering type
;   eadjust:    if set, change !dcs_dsd to 3995.5
;   eint_yn:    only applicable for diffuse mode, if set, I=sum(I[i,j]*dE[i,j]) over the specified time channels
;   e_range:    energy range, effective when t_chan=1&2
;   monfit:     if set then monitor is fitted by a gaussian, and the monitor count is the integrated intensity, otherwise it's just the sum
;               default is set
;   montype:    0-fc0 1-bm1 2-duration, default is 1
;   parent:     message widget parent
;   raw:        if set, data will not be normalized, the output is intended for dcs_cryst_align program
;   t_chan:     0:all 1:elastic 2:specify erange
;   title:      message title
;  output:  all data types are float
;   qty:        [ntchan,ndet] or [ndet,nfiles]-raw or [ntchan,ndet,nfiles]-diffuse(same as raw if diffetyp=0)
;   dqty:       [ntchan,ndet] or [ndet,nfiles]-raw or [ntchan,ndet,nfiles]-diffuse
;   ang:        [11,nfiles], [*,0] contains [rotation,height,x-tran,y-tran,x-tilt,y-tilt,file index,T_sample,T_control,T_setpoint,monitor_count]
;               available when raw or diffuse keywords are set
;   bgrate:     [ndet] overall background rate in cts/hr/det, available when shutter closed files are present
;   dbgrate:    [ndet] error of bgrate
;   comment:    comment line of first file
;   duration:   duration of first file, available when diffuse or raw keywords are set
;   ehist:      [ntchan+1]
;   einc:       Ei
;   emean:      [ntchan]
;   e_range:    if initial values are not set, will save the e_range info
;   error:      1 if any error occurs, the results should be discarded
;               possible errors: 1.different file type 2.shutter closed for all nonbackground data
;                                3.shutter open for all dark count files 4. non existing filename
;   ewid:       diffuse:ewid[tchans]; powder,single:[ntchan], the corresponding energy width for the time channels
;   info:       {ch_wl:ch_wl,ch_ms:ch_ms,ch_srdenom:ch_srdenom,ch_res:ch_res,tsdmin:ts0,temperature:setpoint} is used to check
;               file sameness
;   kfactor:    [ntchan]-(ki/kf)  or a scalor
;   monrate:    monitor rate for normalizing dark count rate in projection calculation
;   tmean:      [ntchan]    time channel in us
;   weight:     sum of integrated intensity of the beam monitor
;               or
;               [nfile] array of integrated intensity of the beam monitor if diffuse keyword is set
pro dm_load_dcsdata,open_path,open_file,qty=qty,dqty=dqty,ang=ang,avgsum=avgsum,$
    bgrate=bgrate,dbgrate=dbgrate,comment=comment,diffetyp=diffetyp,diffuse=diffuse,$
    duration=duration,eadjust=eadjust,ehist=ehist,einc=einc,eint_yn=eint_yn,emean=emean,$
    e_range=e_range,error=error,ewid=ewid,info=info,kfactor=kfactor,monfit=monfit,monrate=monrate,$
    montype=montype,parent=parent,raw=raw,t_chan=t_chan,title=title,tmean=tmean,weight=weight,debug=debug
    darkcount   = arg_present(bgrate)
    defsysv,'!dcs_dsd',exists = exists
    if exists eq 1 then old_dcs_dsd = !dcs_dsd else old_dcs_dsd = 4010.0
    new_dcs_dsd = 3995.5
    if (keyword_set(diffuse) or keyword_set(raw)) then begin
       if n_elements(avgsum)   eq 0 then avgsum = 1         ;default is sum
       if n_elements(t_chan)   eq 0 then t_chan = 1         ;default is elastic peak
       if n_elements(diffetyp) eq 0 then diffetyp = 0       ;default is single energy
    endif else begin
       if n_elements(t_chan)   eq 0 then t_chan = 0         ;default is all time channels
       if n_elements(diffetyp) eq 0 then diffetyp = 1       ;default is multiple energy
    endelse
    if n_elements(title)   eq 0 then title='Loading data ...'
    if n_elements(monfit)  eq 0 then monfit=1   ;default is to fit monitor with a gaussian, and use the integrated intensity as normalization
    if n_elements(montype) eq 0 then montype=1  ;default is BM1
    n_dkfile    = 0
    error       = 0b                            ;initialization - no error
    sch_tch     = (t_chan gt 0)                 ;flag for locating time channels
    nofirstline = 0 & omithisto=0               ;dcs_read_binaryOctavefile keyword
    id_mon      = ([2,3,2])[montype]            ;index of monitor in histohigh[0:ntchan,*]

    n_files     = n_elements(open_file)
    mesg        = widget_base(title=title,/col,/base_align_center,group_leader=parent,tlb_frame_attr=9,map=0)
    ;get the correct path separator
    if float(!version.release) ge 5.5 then pathsep=call_function("path_sep") $
    else begin
       defsysv,'!delimiter',exists=exists
       if not exists then begin
          case !version.os_family of
            'Windows':  defsysv,'!delimiter','\'
            'unix':     defsysv,'!delimiter','/'
            'MacOS':    defsysv,'!delimiter',':'
            else:       defsysv,'!delimiter','/'
          endcase
       endif
       pathsep = !delimiter
    endelse

    ;unexpected errors
    catch, myerror
    if myerror ne 0 then begin
       error = 1b
       catch,/cancel
       ok = dialog_message(!error_state.msg,/error)
       if widget_info(mesg,/valid_id) then widget_control,mesg,/destroy
       if ptr_valid(davePtr) then heap_free,davePtr & tmp=check_math(mask=32)
       defsysv,'!dcs_dsd',old_dcs_dsd
       return
    end

    for i=0L,n_files-1L do begin
        file = open_path+pathsep+open_file[i]
        if not ptr_valid(davePtr) then davePtr=dm_define_pointer()
        if keyword_set(eadjust) then defsysv,'!dcs_dsd',new_dcs_dsd
        if i eq 0 then begin
           label = widget_label(mesg,value='Loading '+file,/dynamic_resize,/align_center)
           dm_center_kid,mesg,parent
           widget_control,mesg,/realize,map=1
        endif else $
           widget_control,label,set_value='loading '+file,tlb_set_title='Loading...  '+$
              strtrim(string(i+1),2)+'/'+strtrim(string(n_files),2)
        ;read in data file
        dcs_read_binaryOctavefile,file,nofirstline,omithisto,var_read,unknown,readerror,davePtr=davePtr
        if strlen(readerror) ne 0 then begin
           error = 1b                           ;error occurs
           if widget_info(mesg,/valid_id) then widget_control,mesg,/destroy      ;clean up
           heap_free,davePtr & tmp=check_math(mask=32)
           ok = dialog_message([file,readerror],/error,dialog_parent=parent)
           defsysv,'!dcs_dsd',old_dcs_dsd
           return
        end
        this          = dm_locate_datastrptr(davePtr)               ;get the data structure pointer
        ts0           = float((*(*this).specificPtr).tsdmin)        ;tsdmin
        ch_wl         = float((*(*this).specificPtr).ch_wl)         ;wavelength
        ch_ms         = float((*(*this).specificPtr).ch_ms)         ;chopper master speed
        ch_srdenom    = (*(*this).specificPtr).ch_srdenom           ;chopper speed ratio denomenator
        ch_res        = (*(*this).specificPtr).ch_res               ;resolution mode
        setpoint      = float((*(*this).specificPtr).temp_setpoint) ;temperature set point
        motor_pos     = (*(*this).specificPtr).motor_pos            ;motor position
        this_duration = (*(*this).specificPtr).duration             ;duration of run in secs
        shut_stat     = (*(*this).specificPtr).shutter_stat         ;shutter status 0:closed 1:open -1:not defined
        temp_sample   = (*(*this).specificPtr).temp_sample          ;sample temperature
        temp_control  = (*(*this).specificPtr).temp_control
        postemp       = where(temp_sample gt 0,count)
        if count ne 0 then begin
           temp_sample  = temp_sample[postemp]
           temp_control = temp_control[postemp]
        endif
        temp_sample   = mean(temp_sample)
        temp_control  = mean(temp_control)
        ;check sameness
        if n_elements(info) eq 0 then begin                      ;save file information
           info = {ch_wl:ch_wl,ch_ms:ch_ms,ch_srdenom:ch_srdenom,ch_res:ch_res,tsdmin:ts0,temperature:setpoint}
        endif else if not keyword_set(raw) then begin            ;check file sameness except when raw keyword is set
           if (info.ch_wl ne ch_wl) or (info.ch_ms ne ch_ms) or (info.ch_srdenom ne ch_srdenom) or $
           (info.tsdmin ne ts0) or (info.ch_res ne ch_res) then begin
              ok = dialog_message('Selected data files are of different types.',/error,dialog_parent=parent)
              error = 1b                        ;error occurs
              widget_control,mesg,/destroy      ;clean up
              heap_free,davePtr & tmp=check_math(mask=32)
              defsysv,'!dcs_dsd',old_dcs_dsd
              return                            ;exit
           endif
        endif
        if i eq 0 then begin
           nhedet   = (*(*this).specificPtr).nhedet               ;number of H3-3 detectors
           comment  = strtrim((*(*this).specificPtr).comments,2)  ;comment
           ntchan   = (*(*this).specificPtr).ntchan               ;number of time channels
           speriod  = 6.0d7/ch_ms*ch_srdenom                      ;time between pulses at sample in micro second
           ; Define the number of good time channels for the He3 detectors, ngtchan0,
           ; and the number of good time channels for the high end inputs , ngtchan1,
           twid     = (*(*this).specificPtr).tchanlook*0.05
           times0   = [0,reform(total(twid[0,*],/cumulative))]
           times1   = [0,reform(total(twid[1,*],/cumulative))]
           tstart   = times0[0:ntchan-1]
           tstop    = times0[1:ntchan]
           chan     = where(tstart le speriod-0.05 and tstop gt speriod-0.05) + 1
           ngtchan0 = chan[0]
           if ngtchan0 eq 0 then begin
              ngtchan0 = 1000L
              warning  = 1b
           endif
           tstart   = times1[0:ntchan-1]
           tstop    = times1[1:ntchan]
           chan     = where(tstart le speriod-0.05 and tstop gt speriod-0.05) + 1
           ngtchan1 = chan[0]
           if ngtchan1 eq 0 then begin
              ngtchan1 = 1000L
              warning  = 1b
           endif
           if keyword_set(warning) then begin
              ok = dialog_message('The time channel look-up table in this file might not be correct.',dialog_parent=parent,title='Warning: '+open_file[i])
              widget_control,label,/show,/input_focus
           endif
           tstart = tstart[0:ngtchan1-1]
           tstop  = tstop[0:ngtchan1-1]
           tmean  = (tstart+tstop)*0.5
           times  = float((*(*(*this).specificPtr).time_propsPtr)[0:ngtchan0])
           if not keyword_set(raw) then begin     ;no use of twid for dcs_cryst_align data
              twid  = times[1:ngtchan0]-times[0:ngtchan0-1]
              w_max = max(twid,min=w_min)
              if w_max eq w_min then $
                 twid = twid[0] $
              else $
                 twid = twid#(fltarr(nhedet)+1.0)
           endif
           einc       = float(!dcs_hsq2mn/(ch_wl)^2)   ;incident energy
           ehist      = (*(*(*this).specificPtr).energy_propsPtr)[0:ngtchan0]
           emean      = float(0.5*(ehist[0:ngtchan0-1]+ehist[1:ngtchan0]))
           ewid       = (ehist[1:ngtchan0]-ehist[0:ngtchan0-1])
           ef         = einc-emean                ;final energy
           ki_kf      = sqrt(einc/ef)             ;ratio of ki/kf
           ntchan     = ngtchan0
           nitem      = ngtchan0
           if keyword_set(avgsum) then nitem = 1.0  ;1 for sum
           if darkcount then begin
              fb_weight  = 0.0                  ;save the duration in hrs, for shutter-closed files
           endif
           if keyword_set(raw) or keyword_set(diffuse) then begin   ;initialization
              duration = this_duration          ;duration of first file
           endif else begin
              if n_elements(qty) eq 0 then begin
                 qty    = 0.0
                 dqty   = 0.0
                 weight = 0.0
              endif
           endelse
        endif
        tmp_qty = (*(*this).commonStr.histPtr).qty
        tmp_err = (*(*this).commonStr.histPtr).err
        ptr_free,(*this).commonStr.histPtr
        ;correct bit flip errors
        tdh_dsum = [long64(total(tmp_qty,1,/double)),long64(total((*(*this).specificPtr).histohigh,1,/double))]
        tdh_tsum = long64(total(tmp_qty,2,/double))+long64(total((*(*this).specificPtr).histohigh,2,/double))
        bad_t = where((tdh_tsum-(*(*this).specificPtr).timsum) ne 0,tcount)    ;disagreement in time sums over detector
        bad_d = where((tdh_dsum-(*(*this).specificPtr).detsum) ne 0,dcount)    ;disagreement in detector sums over time
        if tcount lt 50 and dcount lt 50 then begin
           ;the following code can only correct isolated bitflip errors
           for ii=0,tcount-1 do begin
               widget_control,label,set_value='correcting bitflip errors...',tlb_set_title=open_file[i]
               for jj=0,dcount-1 do begin
                   tmp1 = tdh_tsum[bad_t[ii]]-((*(*this).specificPtr).timsum)[bad_t[ii]]
                   tmp2 = tdh_dsum[bad_d[jj]]-((*(*this).specificPtr).detsum)[bad_d[jj]]
                   if (tmp1 eq tmp2) and (bad_t[ii] lt ngtchan0) and (bad_d[jj] lt nhedet) then begin
                      if dm_ispower2(abs(tmp1)) eq 0 then begin
                         ok = dialog_message(['Bitflip error at det='+dm_to_string(bad_d[jj])+', tchan='+dm_to_string(bad_t[ii])+$
                           " may involve several bits.",'This bitflip error has been corrected.',$
                           '','Detector number and time channel number start from 0.'],dialog_parent=parent,title='Warning: '+open_file[i])
                         widget_control,label,/show,/input_focus
                      endif
                      tmp_qty[bad_t[ii],bad_d[jj]] = tmp_qty[bad_t[ii],bad_d[jj]]-tmp1
                      tmp_err[bad_t[ii],bad_d[jj]] = sqrt(0>(tmp_qty[bad_t[ii],bad_d[jj]]))
                      tdh_dsum = [long64(total(tmp_qty,1,/double)),long64(total((*(*this).specificPtr).histohigh,1,/double))]
                      tdh_tsum = long64(total(tmp_qty,2,/double))+long64(total((*(*this).specificPtr).histohigh,2,/double))
                   endif else begin
                      if n_elements(uncorr) eq 0 then uncorr = [bad_d[jj],bad_t[ii]] $
                      else uncorr = [[uncorr],[bad_d[jj],bad_t[ii]]]
                   endelse
               endfor
           endfor
           if n_elements(uncorr) ne 0 then begin
              tmp     = n_elements(uncorr[0,*])
              be_info = 'The following bitflip error'+(tmp gt 1?'s are':' is')+' not corrected:'
              for kk=0,(30<(tmp-1)) do $
                  be_info = [be_info,string(uncorr[*,kk],format='("    det num = ",i4,"   tchan num = ",i4)')]
              if kk lt tmp then $
                 be_info = [be_info,'            ...','            ...','            ...'] $
              else  $
                 be_info[kk] = be_info[tmp]+' .'
              be_info = [be_info,'','Detector number and time channel number start from 0.',$
                        'Detector number greater than '+dm_to_string(nhedet-1)+' or time channel number',$
                        'greater than '+dm_to_string(ngtchan0-1)+' is for histohigh data.']
              ok = dialog_message(be_info,dialog_parent=parent,title='Warning: '+open_file[i])
              widget_control,label,/show,/input_focus
              tmp = temporary(uncorr) & tmp = -1
           endif
           tdh_dsum = -1 & tdh_tsum=-1
        endif else $
           ok = dialog_message('There are too many bit flip errors. Skip bit flip error correction.',dialog_parent=parent)
        ;check if shutter is closed
        mon_counts = (*(*this).specificPtr).histohigh[0:ngtchan1-1,id_mon]
        mon_sum = total(mon_counts)
        if shut_stat eq -1 then begin
           shutter_close = total(tmp_qty[0:ngtchan0-1,*])/mon_sum
           shutter_close = (shutter_close gt 200)
        endif else $
           shutter_close = 1-shut_stat
        if montype eq 2 then begin
           mon_sum = this_duration
        endif
        if i eq 0 then first_mon_sum = mon_sum
        if ch_wl eq 0 then shutter_close=1  ;regard lambda=0 as shutter closed measurement
        if sch_tch then begin           ;need to locate time channels
           if n_elements(e_range) eq 2 then $
              if total(finite(e_range)) eq 2 then $
                 if t_chan eq 1 then t_chan=2   ;if e_range is given, regard it as specify
           sch_tch = 0                  ;no need to locate tchan again once this is done
           all_tch = 1                  ;show all time channels flag
           case t_chan of
              1: begin                      ;elastic peak
                 if shutter_close eq 1 then begin
                    sch_tch = 1             ;need to do it again
                    break
                 endif
                 tmp_y  = tmp_qty[0:ngtchan0-1,*]
                 tmp_y  = total(tmp_y,2,/double)
                 tmp_ny = n_elements(tmp_y)
                 tmp_x  = findgen(tmp_ny)
                 tmp_ymax = max(tmp_y,imax)
                 ;shift the peak approximately to the center channel
                 tmp_y  = shift(tmp_y,tmp_ny/2.-imax)
                 estimates = [tmp_ymax,tmp_ny/2.,20.0,0.0]
                 result = gaussfit(tmp_x,tmp_y,params,estimates=estimates,nterms=n_elements(estimates),chisq=chisq)
                 if keyword_set(debug) then begin
                    print,params
                    window,/free
                    plot,tmp_x,tmp_y,psym=6,xr=tmp_ny/2.+[-params[2],params[2]]*3,title='chisq='+strtrim(string(chisq),2)
                    oplot,tmp_x,result
                    wait,2 & wdelete,!D.WINDOW
                 endif
                 tmp_x = -1 & tmp_y = -1
                 params[2] =   abs(params[2])
                 if (abs(params[1]-tmp_ny/2.) gt 100) or (params[2] le 1) then begin;somehow the fitting is not right use all time channels
                    ok = dialog_message(["Can't locate the elastic peak from the first file.","All time channels are used instead."],$
                        dialog_parent=parent)
                    widget_control,label,/show,/input_focus
                 endif else begin
                    all_tch = 0     ;fitting is good
                    tchans  = params[1]-(tmp_ny/2.-imax)-fix(2.5*params[2])+indgen(fix(5.0*params[2])+1)
                    index   = where(tchans lt 0,count)
                    if count ne 0 then tchans[index]=tchans[index]+ngtchan0
                    index   = where(tchans ge ngtchan0,count)
                    if count ne 0 then tchans[index]=tchans[index]-ngtchan0
                    emin    = min(emean[tchans],max=emax)
                    e_range = [emin,emax]
                 endelse
                 end
              2: begin             ;specify energy range
                 if n_elements(e_range) eq 2 then begin
                    if total(finite(e_range)) eq 2 then begin
                       tchans = where((emean ge e_range[0]) and (emean le e_range[1]),count)
                       if count ne 0 then begin
                          all_tch = 0
                          break
                       endif
                    endif
                    ok = dialog_message(['The energy range specified is invalid.','All time channels are used instead.'],$
                         dialog_parent=parent)
                    widget_control,label,/show,/input_focus
                 endif
                 end
              else:
           endcase
           if all_tch eq 1 then begin  ;all time channels
              tchans  = 0 & tmp = temporary(tchans) ;get rid of tchans
              emin    = min(emean,max=emax)
              e_range = [emin,emax]
           endif
           ntchan = n_elements(tchans)
           if (sch_tch eq 0) then begin
              if diffetyp eq 1 then begin
                 if ntchan ne 0 then begin
                    if n_elements(qty) gt 1 then begin
                       qty  = qty[tchans,*,*]
                       dqty = dqty[tchans,*,*]
                    endif
                 endif
              endif
              if ntchan ne 0 then begin
                 if n_elements(twid) ne 1 then twid = twid[tchans,*]
                 ewid  = ewid[tchans]
                 emean = emean[tchans]
                 ki_kf = ki_kf[tchans]
              endif
           endif
           if ntchan eq 0 then ntchan = ngtchan0
           nitem = (avgsum eq 0)?(1.0>ntchan):1.0      ;1 for sum
        endif
        if ntchan ne ngtchan0 then begin
           tmp_qty = tmp_qty[tchans,*]
           tmp_err = tmp_err[tchans,*]
        endif else begin
           tmp_qty = tmp_qty[0:ngtchan0-1,*]
           tmp_err = tmp_err[0:ngtchan0-1,*]
        endelse
        if keyword_set(raw) then begin  ;for dcs_cryst_align, no correction or normalization
           factor  = float(first_mon_sum/mon_sum)/nitem
           if i eq 0 then begin
              qty  = total(tmp_qty,1)*factor
              dqty = sqrt(total(tmp_err^2,1))*factor
           endif else begin
              qty  = [[qty],[total(tmp_qty,1)*factor]]
              dqty = [[dqty],[sqrt(total(tmp_err^2,1))*factor]]
           endelse
           if n_elements(ang) eq 0 then $
              ang = [motor_pos[1:6],i,temp_sample,temp_control,setpoint,mon_sum] $
           else $
              ang = [[ang],[[motor_pos[1:6],i,temp_sample,temp_control,setpoint,mon_sum]]]
        endif else begin                ;for dcs_mslice
           if shutter_close then begin  ;shutter closed
              n_dkfile = n_dkfile+1     ;record the number of shutter closed files
              if darkcount then begin   ;add to darkcount rate if background keyword is set
                 fb_weight = fb_weight+float(this_duration/3600.0)   ;in hours
                 tmp_qty   = total(tmp_qty,1)
                 if n_elements(bgrate) eq 0 then bgrate=temporary(tmp_qty) $
                 else bgrate=temporary(bgrate)+temporary(tmp_qty)
              endif else begin          ;discard the data if not measuring dark count file
                 if ch_wl eq 0 then $
                    tmp = 'Wavelength is 0 for this file. ' $
                 else $
                    tmp = 'Shutter is closed for this file. '
                 widget_control,label,set_value=temporary(tmp)+'This file is discarded.',tlb_set_title=open_file[i]
                 wait,2
              endelse
           endif else begin                 ;shutter open
              if darkcount then begin
                 widget_control,label,set_value='Shutter is open. This file is discarded.',tlb_set_title=open_file[i]
                 wait,2
                 continue
              endif
              if keyword_set(monfit) and (montype eq 1) then begin  ;fit to a gaussian only for bm1
                 mon_max = max(mon_counts,imax) & imax=imax[0]
                 tbm_exp = (tstart[imax]+tstop[imax])/2.0
                 ;shift the peak time to roughly the center of pulse
                 tmp     = (tmean+(speriod/2.-tbm_exp)+speriod) mod speriod
                 params  = [mon_max,speriod/2.,10.0,0.0]
                 if keyword_set(debug) then print,'estimates=',params
                 result  = gaussfit(tmp,mon_counts,params,estimates=params,nterms=n_elements(params))
                 if keyword_set(debug) then begin
                    print,'parameter=',params
                    window,/free
                    plot,tmean,mon_counts,psym=6,xr=tbm_exp+[-params[2],params[2]]*3.
                    oplot,tmean,result
                    wait,2 & wdelete,!D.WINDOW
                 endif
                 normalization_values = float(abs(params[0]*params[2])) ;integrated intensity ~ monrate* npulse
              endif else normalization_values = mon_sum
              ;monrate = float(normalization_values*speriod/this_duration)*1e-6
              if n_elements(monrate) eq 0 then monrate = [normalization_values,speriod,this_duration] $
              else monrate = [[monrate],[[normalization_values,speriod,this_duration]]]
              ;divide (*histPtr).qty and (*histPtr).err by the time channel width
              if keyword_set(diffuse) then begin
                 if n_elements(weight) eq 0 then $
                    weight = normalization_values $
                 else $
                    weight = [weight,normalization_values]
                 if n_elements(ang) eq 0 then $
                    ang = [motor_pos[1:6],i,temp_sample,temp_control,setpoint,mon_sum] $
                 else $
                    ang = [[ang],[[motor_pos[1:6],i,temp_sample,temp_control,setpoint,mon_sum]]]
                 tmp_qty   = tmp_qty/twid    ;~I * npulse
                 tmp_err   = tmp_err/twid
                 if diffetyp eq 0 then begin        ;single energy
                    if keyword_set(eint_yn) then begin
                       for j=0,n_elements(tmp_qty[0,*])-1 do begin
                           tmp_qty[*,j] = tmp_qty[*,j]*ewid
                           tmp_err[*,j] = tmp_err[*,j]*ewid
                       endfor
                       factor = 1e6/normalization_values
                    endif else begin
                       factor = 1e6/nitem/normalization_values
                       nmonrate = n_elements(monrate[0,*])
                       monrate[0,nmonrate-1] = monrate[0,nmonrate-1]*nitem
                    endelse
                    if i eq 0 then begin
                       qty  = total(tmp_qty,1)*factor
                       dqty = sqrt(total(tmp_err^2,1))*factor
                    endif else begin
                       qty  = [[qty],[total(tmp_qty,1)*factor]]
                       dqty = [[dqty],[sqrt(total(tmp_err^2,1))*factor]]
                    endelse
                 endif else begin                   ;multiple energy
                    factor = 1e6/normalization_values
                    if i eq 0 then begin
                       qty  = tmp_qty*factor
                       dqty = tmp_err*factor
                    endif else begin
                       qty  = [[[qty]],[[tmp_qty*factor]]]
                       dqty = [[[dqty]],[[tmp_err*factor]]]
                    endelse
                 endelse
              endif else begin
                 if n_elements(ang) eq 0 then $
                    ang = [motor_pos[1:6],i,temp_sample,temp_control,setpoint,mon_sum]
                 weight = temporary(weight)+normalization_values
                 qty  = temporary(qty)+tmp_qty/twid
                 dqty = sqrt(temporary(dqty)^2+(tmp_err/twid)^2)
              endelse
           endelse
        endelse
        tmp_qty = 0 & tmp_err = 0
    endfor

    defsysv,'!dcs_dsd',old_dcs_dsd
    twid = 0    ;clear memory
    tmp = check_math(mask=32)    ;Clear any floating-point underflows that has occurred
    widget_control,mesg,/destroy
    heap_free,davePtr

    if keyword_set(raw) then return

    ;monitor rate
    if n_elements(monrate) ne 0 then begin
       monrate = float(total(monrate[0,*])*monrate[1,0]/total(monrate[2,*]))*1e-6
    endif

    if darkcount then begin
       if n_elements(bgrate) ne 0 then begin
          dbgrate = sqrt(bgrate)/fb_weight
          bgrate  = bgrate/fb_weight
       endif else begin
          error = 1
          ok    = dialog_message('All files are with shutter open.',/error,dialog_parent=parent)
       endelse
    endif else begin
       if n_dkfile ne n_files then begin
          kfactor = ki_kf
          if keyword_set(diffuse) then begin
             if diffetyp eq 0 then begin
                if t_chan eq 2 then begin
                   emean   = mean(emean)
                   kfactor = mean(kfactor)
                endif else begin
                   emean   = 0.0
                   kfactor = 1.0
                endelse
             endif
          endif else begin
             factor  = float(1d6/weight)
             weight  = float(weight)
             qty     = temporary(qty)*factor
             dqty    = temporary(dqty)*factor
          endelse
       endif else begin
          error = 1b
          ok    = dialog_message('All files are with shutter closed.',/error,dialog_parent=parent)
       endelse
    endelse
end

;the following procedure is for loading empty can file for dcs_mslice
pro dcs_mslice::dm_load_dcs_emptycan,open_path,open_file,parent=parent,error=error,comment=comment,$
    qty=qty,dqty=dqty,ang=ang,weight=weight
    tmp = self.e_range
    if self.samp_typ ge 2 then begin
       if self.e_bin[self.samp_typ] ne 2 then self.e_range=[!values.f_nan,!values.f_nan]
       dm_load_dcsdata,open_path,open_file,ang=ang,qty=qty,dqty=dqty,weight=weight,$
            comment=comment,ewid=ewid,t_chan=self.e_bin[self.samp_typ],e_range=self.e_range,$
            eint_yn=self.eint_yn,avgsum=self.bin_avgsum,parent=parent,error=error,$
            /diffuse,monfit=self.mon_sumfit,title='loading empty can file...',$
            eadjust=self.eadjust,diffetyp=(self.samp_typ eq 3),emean=emean,montype=self.dcs_mon
       ang = reform(ang[0,*])   ;only need a2 angle
    endif else begin
       ang = 0
       dm_load_dcsdata,open_path,open_file,qty=qty,dqty=dqty,parent=parent,$
            title='loading empty can file...',error=error,monfit=self.mon_sumfit,$
            comment=comment,eadjust=self.eadjust,montype=self.dcs_mon
    endelse
    self.e_range = tmp
end

;the following procedure is for loading dark count files for dcs_mslice
pro dcs_mslice::dm_load_dcs_darkcount,open_path,open_file,parent=parent,error=error,comment=comment,$
    bgrate=bgrate,dbgrate=dbgrate
    dm_load_dcsdata,open_path,open_file,bgrate=bgrate,dbgrate=dbgrate,error=error,comment=comment,$
    parent=parent,title='loading dark count file...'
end

;the following procedure is for loading detector dependent background file for dcs_mslice
pro dcs_mslice::dm_load_dcs_detbackground,open_path,open_file,parent=parent,error=error,comment=comment,$
    qty=qty,dqty=dqty
    dm_load_dcsdata,open_path,open_file,qty=qty,dqty=dqty,t_chan=2,$
        e_range=self.bgeran,avgsum=0,title='loading background file...',comment=comment,$
        parent=parent,error=error,monfit=self.mon_sumfit,$
        eadjust=self.eadjust,montype=self.dcs_mon
    ntch = n_elements(qty[*,0])
    qty  = total(qty,1)
    dqty = total(dqty^2,1)
    qty  = qty/ntch
    dqty = sqrt(dqty)/ntch
end

;the following procedure is for loading time channel background file for dcs_mslice
pro dcs_mslice::dm_load_dcs_tchanbackground,open_path,open_file,parent=parent,error=error,comment=comment,$
    qty=qty,dqty=dqty,emean=emean
    dm_load_dcsdata,open_path,open_file,qty=qty,dqty=dqty,t_chan=2,$
        e_range=self.dcor_eqran[0:1],avgsum=0,title='loading background file...',comment=comment,$
        parent=parent,error=error,diffetyp=1,monfit=self.mon_sumfit,$
        eadjust=self.eadjust,montype=self.dcs_mon,emean=emean
end

;the following procedure is for loading vanadium file for dcs_mslice
pro dcs_mslice::dm_load_dcs_vanadium,open_path,open_file,parent=parent,error=error,$
    comment=comment,qty=qty,ei=ei
    dm_load_dcsdata,open_path,open_file,qty=qty,parent=parent,title='loading vanadium file...',error=error,$
                monfit=self.mon_sumfit,comment=comment,info=info,eadjust=self.eadjust,montype=self.dcs_mon,$
                einc=ei
end


;the following procedure is for loading dcs data when "Load Data" button is pressed in dcs_mslice
pro dcs_mslice::dm_load_dcs,open_path,open_file,parent=parent,error=error,comment=comment,$
    qty=qty,dqty=dqty,weight=weight,ang=ang,ei=ei,info=info,ewid=ewid,emean=emean,monrate=monrate,$
    kfactor=kfactor,title=title,tmean=tmean
    if self.samp_typ ge 2 then begin
       e_range = self.e_range
       dm_load_dcsdata,open_path,open_file,ang=ang,qty=qty,dqty=dqty,weight=weight,info=info,comment=comment,$
           monrate=monrate,ewid=ewid,kfactor=kfactor,t_chan=self.e_bin[self.samp_typ],e_range=e_range,eint_yn=self.eint_yn,$
           avgsum=self.bin_avgsum,parent=parent,error=error,/diffuse,monfit=self.mon_sumfit,title=title,$
           eadjust=self.eadjust,diffetyp=(self.samp_typ eq 3),emean=emean,montype=self.dcs_mon,einc=ei
       if n_elements(e_range) ne 0 then self.e_range=e_range
       ang = reform(ang[0,*])   ;only need a2 angle
    endif else begin
       dm_load_dcsdata,open_path,open_file,qty=qty,dqty=dqty,weight=weight,emean=emean,info=info,ewid=ewid,$
           monrate=monrate,kfactor=kfactor,parent=parent,comment=comment,error=error,monfit=self.mon_sumfit,$
           title=title,eadjust=self.eadjust,montype=self.dcs_mon,ang=ang,tmean=tmean,einc=ei
       ang = ang[0]             ;only need a2 angle
    endelse
    if self.samp_typ ne 0 then tmean=0
end