; $Id: dm_load_spe.pro,v 1.18 2008/03/19 20:46:50 ymqiu Exp $
;#######################################################################
;
; NAME:
;  dm_load_spe
;
; PURPOSE:
;  load spe file returns zdat[nx,ny], err[nx,ny], xdat[nx or nx+1], ydat[ny or ny+1]
;
; 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.
;
;#######################################################################

pro dm_load_spe,unit,zdat=zdat,err=err,xdat=xdat,ydat=ydat
    nx = 0L & ny = 0L & junk=' '
    format1 = '(8(g10.5))'
    format2 = '(2(i5))'

    readf,unit,junk
    if strlen(junk) eq 11 then begin
       format1 = '(((g11.5),7(g10.5)))'
       format2 = '((i6,i5))'
    endif

    tmp = strtrim(junk,2)
    if (strpos(tmp,' '))[0] eq -1 then begin
       reads,junk,nx,ny,format=format2
    endif else begin
       reads,junk,nx,ny
    endelse

    zdat = fltarr(nx,ny)
    err  = fltarr(nx,ny)
    xdat = fltarr(nx+1)
    ydat = fltarr(ny+1)
    readf,unit,junk    ;read ### phi grid
    readf,unit,xdat,format=format1
    readf,unit,junk    ;read ### Energy Grid
    readf,unit,ydat,format=format1

    ;determine whether the data is point type or histogram type
    if (xdat[nx] eq 0) and (nx ge 2) then begin
       if xdat[nx-1]-xdat[0] ge 0 then $
          tmp = xdat[1:nx]-xdat[0:nx-1] $
       else $
          tmp = xdat[0:nx-1]-xdat[1:nx]
       tmp1 = where(tmp ge 0,count)
       if count ne nx then $    ;not monotonic
          xdat = xdat[0:nx-1] $
       else if tmp[nx-1] gt 5.*tmp[nx-2] then $    ;arbitrary criterion
          xdat = xdat[0:nx-1]
    endif
    if (ydat[ny] eq 0) and (ny gt 2) then begin
       if ydat[ny-1]-ydat[0] ge 0 then $
          tmp = ydat[1:ny]-ydat[0:ny-1] $
       else $
          tmp = ydat[0:ny-1]-ydat[1:ny]
       tmp1 = where(tmp ge 0,count)
       if count ne ny then $    ;not monotonic
          ydat = ydat[0:ny-1] $
       else if tmp[ny-1] gt 5.*tmp[ny-2] then $    ;arbitrary criterion
          ydat = ydat[0:ny-1]
    endif
    tmp  = fltarr(ny)
    for i=0L,nx-1L do begin
        readf,unit,junk   ;read ### S(phi,w)
        readf,unit,tmp,format=format1;format
        zdat[i,*] = tmp
        readf,unit,junk   ;read ### Errors
        readf,unit,tmp,format=format1;format
        err[i,*]  = tmp
    endfor
    free_lun,unit
end

;this procedure reads a spe file for dcs_mslice program
;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 sest, default is self.bin_avgsum
;   background: if set, the loaded files are regarded as background files and will be summed
;   diffetyp:   only needed when diffuse keyword is set 0:single energy  1:multiple energy
;   diffuse:    if set, diffuse scattering mode
;   parent:     message widget parent,  default is self.tlb
;   title:      message title
;   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 - diffuse keyword set, default is self.e_range
;   t_chan:     0:all 1:elastic 2:specify erange  - efffective when diffuse keyword set, default is self.e_bin[self.samp_typ]
;   phxfile:    flag for reading phxfile, if set the detector pointer will be created no matter what
;   eief:       if present and given a value, the procedure will not prompt for Ei or Ef value from user
;   masklimit:  a number indicating the threshhold below which the detector is masked, default:-1e20
; output:
;   error:      error flag 1 or 0
;   qty:        [ne,ndet] or [ne,ndet,nfiles]-diffuse keywords set [ne,ndet] if nfiles=1
;   dqty:       [ne,ndet] or [ne,ndet,nfiles]-diffuse keywords set [ne,ndet] if nfiles=1
;   ang:        [nfiles] - diffuse keyword set, not in spe file, need to be typed in by the user
;   eief:       Ei or Ef(depends on the instrument geometry), not in the spe file, need to be typed in by the user
;   emean:      [ne]
;   ewid:       diffuse:ewid[tchans] ; powder,single:[ntchan], the corresponding energy width for the time channels
;   info:       {ne:ne, ndet:ndet}  used for checking file compatibility
;   weight:     weight or duration info , not in the file, need to be typed in by the user
;   kfactor:    [ntchan]-(ki/kf)  or a scalor
pro dcs_mslice::dm_load_spe,open_path,open_file,ang=ang,avgsum=avgsum,background=background,$
    diffetyp=diffetyp,diffuse=diffuse,eief=eief,e_range=e_range,eint_yn=eint_yn,emean=emean,$
    error=error,ewid=ewid,info=info,kfactor=kfactor,parent=parent,phxfile=phxfile,qty=qty,dqty=dqty,$
    t_chan=t_chan,title=title,masklimit=masklimit

    if n_elements(masklimit)   eq 0 then masklimit = -1e20
    if keyword_set(diffuse) then begin
       if n_elements(avgsum)   eq 0 then avgsum   = self.bin_avgsum
       if n_elements(t_chan)   eq 0 then t_chan   = self.e_bin[self.samp_typ]
       if n_elements(e_range)  ne 2 then e_range  = self.e_range
       if n_elements(diffetyp) eq 0 then diffetyp = (self.samp_typ eq 3)
    endif
    if n_elements(parent) eq 0 then parent = self.tlb
    if n_elements(title)  eq 0 then title  = 'Loading data ...'
    error   = 0b  ;initialize no error
    n_files = n_elements(open_file)

    mesg = widget_base(title=title,/col,/base_align_center,group_leader=parent)

    ;catch and clear possible io errors
    catch, myerror
    if myerror ne 0 then begin
       error = 1b
       catch,/cancel
       if widget_info(mesg,/valid_id) then widget_control,mesg,/destroy
       ok = dialog_message(!error_state.msg,/error)
       if n_elements(unit) ne 0 then free_lun,unit
       return
    end

    ;ask for ei
    if arg_present(eief) then begin
       prompt = (['Ei (meV):','Ef (meV):'])[0>(self.instrgeom)<1]
       if n_elements(eief) ne 1 then eief = !values.f_nan
       while(finite(eief) eq 0) do $
          eief = dm_dialog_input(prompt,/float,dialog_parent=parent,/e2lambda)
    endif else eief = 1e3  ;arbitrary number

    ;ask for weight info
    myweight = 1.0
    if arg_present(weight) or (n_files ne 1) then begin
       myweight = !values.f_nan
       while(total(finite(myweight)) ne n_files) do $
            myweight = dm_dialog_listinput(open_file,title='Please input weight info:',default=1,$
                 label='weight',dialog_parent=parent,reset=1,info=['The weight is for summing files.',$
                 "Press OK if you don't know the information."])
    endif

    ;ask for rotation angles
    if keyword_set(diffuse) and (not keyword_set(background)) then begin
       ang = !values.f_nan
       while(total(finite(ang)) ne n_files) do $
          ang = dm_dialog_listinput(open_file,title='Please input rotation angles:',default=0,$
                label='angle',dialog_parent=parent,reset=0)
    endif

    for i=0L,n_files-1L do begin
        file = open_path+self.pathsep+open_file[i]
        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=strtrim(title,2)+'  '+$
                strtrim(string(i+1),2)+'/'+strtrim(string(n_files),2)

        openr,unit,file,/get_lun,error=fileerror
        if fileerror ne 0 then begin
           error = 1b
           widget_control,mesg,/destroy
           ok = dialog_message("Can't open "+file,/error,dialog_parent=parent)
           return
        endif

        ;read the spe file
        dm_load_spe,unit,zdat=zdat,err=dzdat,xdat=xdat,ydat=ydat    ;xdat:tth, ydat:en

        if n_elements(mskcount) eq 0 then $
           mask = where(zdat[*,0] lt masklimit,mskcount,complement=nomask)
        if mskcount gt 0 then begin
           zdat  = zdat[nomask,*]
           dzdat = dzdat[nomask,*]
           xdat  = xdat[nomask]
        endif

        ndet = n_elements(xdat)
        nen  = n_elements(ydat)
        zdim = size(zdat,/dimension)

        if n_elements(info) eq 0 then begin
           info={ndet:ndet,nen:nen,nx:zdim[0],ny:zdim[1]}
        endif else begin
           if (info.ndet ne ndet) or (info.nen ne nen) or (info.nx ne zdim[0]) or (info.ny ne zdim[1]) 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
              return
           endif
        endelse

        if i eq 0 then begin
           emean = temporary(ydat)
           tth   = temporary(xdat)
           dtth  = fltarr(zdim[0])
           psi   = fltarr(zdim[0])
           dpsi  = fltarr(zdim[0])
           ewid  = emean[1:nen-1]-emean[0:nen-2]
           if n_elements(ewid) ne zdim[1] then ewid = [ewid,2*ewid[nen-2]-ewid[nen-3]]
           if ndet eq zdim[0]+1 then tth   = (tth[0:ndet-2]+tth[1:ndet-1])/2.
           if nen  eq zdim[1]+1 then emean = (emean[0:nen-2]+emean[1:nen-1])/2.

           if keyword_set(diffuse) then begin
              weight  = myweight
              all_tch = 1
              kfactor = 1.0
              ;search the proper time channels to be summed or averaged
              case t_chan of
                   0:          ;all time channels
                   1:  begin   ;elastic peak
                       tmp_y    = total(zdat,1)
                       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 tmp_ny/5.) or (params[2] le 1) then begin
                          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]+tmp_ny
                          index   = where(tchans ge tmp_ny,count)
                          if count ne 0 then tchans[index] = tchans[index]-tmp_ny
                          emin    = min(emean[tchans],max=emax)
                          e_range = [emin,emax]
                       endelse
                       end
                   2:  begin   ;specified energy range
                       if n_elements(e_range) eq 2 then $
                          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
                                if arg_present(kfactor) then begin
                                   if self.instrgeom eq 0 then begin   ;direct geometry
                                      kfactor = sqrt(eief/(eief-mean(e_range))) ;ki/kf_avg
                                   endif else begin   ;inverse geometry
                                      kfactor = sqrt((eief+mean(e_range))/eief) ;ki_avg/kf
                                   endelse
                                endif
                                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
                       end
                   else:
                endcase
                if all_tch eq 1 then begin
                   tchans = indgen(zdim[1])
                   emin   = min(emean,max=emax)
                   e_range = [emin,emax]
                endif
                ewid = ewid[tchans]
                nitem = (avgsum eq 0)?(1.0>(n_elements(ewid))):1.0
                if keyword_set(eint_yn) then factor=1.0 else factor = 1.0/nitem
                if (diffetyp eq 1) and arg_present(kfactor) then begin
                   if self.instrgeom eq 0 then begin   ;direct geometry
                      kfactor = (sqrt(eief/(eief-emean)))[tchans]
                   endif else begin   ;inverse geometry
                      kfactor = (sqrt((eief+emean)/eief))[tchans]
                   endelse
                endif
            endif else begin
                qty    = 0.0
                dqty   = 0.0
                weight = 0.0
            endelse
        endif

        ;add or append data
        if keyword_set(diffuse) then begin
           zdat  = zdat[*,tchans]
           dzdat = dzdat[*,tchans]
           if diffetyp eq 0 then begin
              if keyword_set(eint_yn) then begin
                 for j=0L,n_elements(zdat[*,0])-1 do begin
                     zdat[j,*] = zdat[j,*]*ewid
                     zdat[j,*] = zdat[j,*]*ewid
                 endfor
              endif
              if i eq 0 then begin
                 qty  = total(zdat,2)*factor
                 dqty = sqrt(total(dzdat^2,2))*factor
              endif else begin
                 qty  = [[qty],[total(zdat,2)*factor]]
                 dqty = [[dqty],[sqrt(total(dzdat^2,2))*factor]]
              endelse
           endif else begin
              if i eq 0 then begin
                 qty  = transpose(zdat)
                 dqty = transpose(dzdat)
              endif else begin
                 qty  = [[[qty]],[[transpose(zdat)]]]
                 dqty = [[[dqty]],[[transpose(dzdat)]]]
              endelse
           endelse
        endif else begin
           qty    = temporary(qty)+temporary(zdat*myweight[i])
           dqty   = temporary(dqty)+temporary(dzdat*myweight[i])
           weight = temporary(weight)+myweight[i]
        endelse
    endfor

    if not keyword_set(diffuse) then begin
       if arg_present(kfactor) then begin
          if self.instrgeom eq 0 then begin   ;direct geometry
             kfactor = sqrt(eief/(eief-emean))
          endif else begin   ;inverse geometry
             kfactor = sqrt((eief+emean)/eief)
          endelse
       endif
       if weight eq 0 then weight = 1.0
       qty  = qty/weight
       dqty = dqty/weight
       qty  = transpose(qty)   ;->[nen,ndet]
       dqty = transpose(dqty)  ;->[nen,ndet]
    endif else begin
       if diffetyp eq 1 then emean=emean[tchans] $
       else begin
          if t_chan eq 2 then emean=mean(emean[tchans]) $
          else emean = 0.0
       endelse
       if keyword_set(background) then begin    ;sum the files
          ind  = (size(qty))[0]
          totw = total(weight)
          if totw le 0 then totw = 1.0
          weight = weight/totw
          for i=0,n_elements(weight)-1 do begin
              qty[*,i]  = qty[*,i]*weight[i]
              dqty[*,i] = (dqty[*,i]*weight[i])^2
          endfor
          qty    = total(qty,ind)
          dqty   = sqrt(total(dqty,ind))
          weight = totw
       endif
    endelse

    if mskcount eq 0 then tmp =temporary(nomask)  ;erase nomask

    ;read the phx file
    if keyword_set(phxfile) then begin
       widget_control,label,set_value='Loading detector file...'
       self->dm_load_phx,/phx,tth=tth,psi=psi,dtth=dtth,dpsi=dpsi,nomask=nomask
    endif
    widget_control,mesg,/destroy
end

;the following procedure is for loading detector dependent background file for dcs_mslice
pro dcs_mslice::dm_load_spe_detbackground,open_path,open_file,error=error,qty=qty,dqty=dqty,$
    parent=parent
    self->dm_load_spe,open_path,open_file,qty=qty,dqty=dqty,t_chan=2,e_range=self.bgeran,avgsum=0,$
        title='loading background file...',parent=parent,error=error,/diffuse,/background,diffetyp=0
end

;the following procedure is for loading time channel background file for dcs_mslice
pro dcs_mslice::dm_load_spe_tchanbackground,open_path,open_file,error=error,qty=qty,dqty=dqty,$
    parent=parent,emean=emean
    self->dm_load_spe,open_path,open_file,qty=qty,dqty=dqty,t_chan=2,e_range=self.dcor_eqran[0:1],$
        avgsum=0,title='loading background file...',parent=parent,error=error,/diffuse,/background,$
        diffetyp=1,emean=emean
end