; $Id: dm_choose_cutdata.pro,v 1.20 2008/03/19 20:46:50 ymqiu Exp $
;#######################################################################
;
; NAME:
;  dm_choose_cutdata
;
; PURPOSE:
;  choose cut data
;
; 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.
;
;#######################################################################

; parameters:
;   None
; keywords:
;   xdat,ydat,yerr: data to be returned
;   (xy)tit:        label of the data
pro dcs_mslice::dm_choose_cutdata,xdat=xdat,ydat=ydat,yerr=yerr,xtit=xtit,ytit=ytit,debug=debug
    xp = [self.view_u1typ, self.view_u2typ]  ;powder x-axis
    x  = [xp[self.cut_x[0]],3+self.cut_x[1],3+self.cut_x[2],3+self.cut_x[3]]
    if keyword_set(debug) then begin
       print,'^^choose cut data...'
       current=systime(/sec)
    endif

    emax = max((*self.projStrPtr).en,min=emin)
    dE   = emax-emin    ;default energy integration width

    ;ewid
    if self.eint_yn and (self.samp_typ ne 2) and (self.cut_y[self.samp_typ] eq 0) then $
       self->dm_choose_trajdata,-2,data=ewid
    ok_ewid = (n_elements(ewid) ne 0)


    ;xdat
    self->dm_choose_trajdata,x[self.samp_typ],data=xdat,label=xtit
    if (xtit eq 'E (meV)') or (xtit eq 'Energy') then begin
        if self.samp_typ eq 3 then $
           uniq_E = xdat[*,0,0] $
        else $
           uniq_E = xdat[*,0]
    endif

    ndat = n_elements(xdat)

    if keyword_set(debug) then print,'choose x done',systime(/sec)-current,' sec'
    viewlbs = [self.view_u1lab[(self.samp_typ ge 1)],self.view_u2lab[(self.samp_typ ge 1)],self.view_u3lab]
    ;thickness
    case self.samp_typ of
       0:begin     ;powder
         if total(finite(self.cut_t1ran)) ne 0 then begin
            self->dm_choose_trajdata,xp[(self.cut_x[0]+1) mod 2],index=index1,range=self.cut_t1ran,empty=empty
            if empty then begin
               index1 = -1 & tmp = temporary(index1)
               emptyitem = viewlbs[(self.cut_x[0]+1) mod 2]
            endif else if n_elements(index1) eq ndat then begin
               index1 = -1 & tmp = temporary(index1)
            endif
         endif
         if ok_ewid then begin  ;figure out dE
            if xp[self.cut_x[0]] eq 0 then dE=(finite(self.cut_step)?(self.cut_step):dE) $
            else if xp[(self.cut_x[0]+1) mod 2] eq 0 then begin
                 if total(finite(self.cut_t1ran)) eq 2 then $
                    if self.cut_t1ran[1] gt self.cut_t1ran[0] then $
                       dE = self.cut_t1ran[1]-self.cut_t1ran[0]
            endif
         endif
         end
       1:begin     ;single crystal
         t1_id = 3+((self.cut_x[1]+1) mod 3)
         t2_id = 3+((self.cut_x[1]+2) mod 3)
         if total(finite(self.cut_t1ran)) ne 0 then begin
            self->dm_choose_trajdata,t1_id,index=index1,range=self.cut_t1ran,empty=empty
            if empty then begin
               index1 = -1 & tmp = temporary(index1)
               emptyitem = viewlbs[(self.cut_x[1]+1) mod 3]
            endif else if n_elements(index1) eq ndat then begin
               index1 = -1 & tmp = temporary(index1)
            endif
         endif
         if total(finite(self.cut_t2ran)) ne 0 then begin
            self->dm_choose_trajdata,t2_id,index=index2,range=self.cut_t2ran,empty=empty
            if empty then begin
               index2 = -1 & tmp = temporary(index2)
               if n_elements(emptyitem) eq 0 then emptyitem = viewlbs[(self.cut_x[1]+2) mod 3] else $
               emptyitem = emptyitem+', '+viewlbs[(self.cut_x[1]+2) mod 3]
            endif else if n_elements(index2) eq ndat then begin
               index2 = -1 & tmp = temporary(index2)
            endif
         endif
         if self.disp_fold then begin
            ;fold xaxis
            x_id = 3+((self.disp_perp+1) mod 3)
            if x_id eq t1_id then begin
               if total(finite(self.cut_t1ran)) ne 0 then begin
                  self->dm_choose_trajdata,t1_id,index=tmp,range=[-self.cut_t1ran[1],-self.cut_t1ran[0]],empty=empty
                  if empty eq 0 then begin
                     if n_elements(index1) eq 0 then index1 = tmp $
                     else index1 = [index1,tmp]
                  endif
               endif
            endif else if x_id eq t2_id then begin
               if total(finite(self.cut_t2ran)) ne 0 then begin
                  self->dm_choose_trajdata,t2_id,index=tmp,range=[-self.cut_t2ran[1],-self.cut_t2ran[0]],empty=empty
                  if empty eq 0 then begin
                     if n_elements(index2) eq 0 then index2 = tmp $
                     else index2 = [index2,tmp]
                  endif
               endif
            endif
         endif
         if ok_ewid then begin ;figure out dE
            ind_e = where([self.view_u1[3],self.view_u2[3],self.view_u3[3]] eq 1,count)
            if count eq 1 then begin
               if self.cut_x[1] eq ind_e[0] then dE=(finite(self.cut_step)?(self.cut_step):dE) $
               else if ((self.cut_x[1]+1) mod 3) eq ind_e[0] then begin
                    if total(finite(self.cut_t1ran)) eq 2 then $
                       if self.cut_t1ran[1] gt self.cut_t1ran[0] then $
                          dE=self.cut_t1ran[1]-self.cut_t1ran[0]
               endif else if ((self.cut_x[1]+2) mod 3) eq ind_e[0] then begin
                    if total(finite(self.cut_t2ran)) eq 2 then $
                       if self.cut_t2ran[1] gt self.cut_t2ran[0] then $
                          dE=self.cut_t2ran[1]-self.cut_t2ran[0]
               endif
            endif
         endif
         end
       2:begin     ;diffuse scattering single E
         if total(finite(self.cut_t1ran)) ne 0 then begin
            self->dm_choose_trajdata,3+((self.cut_x[2]+1) mod 2),index=index1,range=self.cut_t1ran,empty=empty
            if empty then begin
               index1 = -1 & tmp = temporary(index1)
               emptyitem = viewlbs[(self.cut_x[2]+1) mod 2]
            endif else if n_elements(index1) eq ndat then begin
               index1 = -1 & tmp = temporary(index1)
            endif
         endif
         end
       3:begin
         t1_id  = 3+((self.cut_x[3]+1) mod 3)
         t2_id  = 3+((self.cut_x[3]+2) mod 3)
         if total(finite(self.cut_t1ran)) ne 0 then begin
            self->dm_choose_trajdata,t1_id,index=index1,range=self.cut_t1ran,empty=empty
            if empty then begin
               index1 = -1 & tmp = temporary(index1)
               emptyitem = viewlbs[(self.cut_x[3]+1) mod 3]
            endif else if n_elements(index1) eq ndat then begin
               index1 = -1 & tmp = temporary(index1)
            endif
         endif
         if total(finite(self.cut_t2ran)) ne 0 then begin
            self->dm_choose_trajdata,t2_id,index=index2,range=self.cut_t2ran,empty=empty
            if empty then begin
               index2 = -1 & tmp = temporary(index2)
               if n_elements(emptyitem) eq 0 then emptyitem = viewlbs[(self.cut_x[3]+2) mod 3] else $
               emptyitem = emptyitem+', '+viewlbs[(self.cut_x[3]+2) mod 3]
            endif else if n_elements(index2) eq ndat then begin
               index2 = -1 & tmp = temporary(index2)
            endif
         endif
         if self.disp_fold then begin
            ;fold xaxis
            x_id = 3+((self.disp_perp+1) mod 3)
            if x_id eq t1_id then begin
               if total(finite(self.cut_t1ran)) ne 0 then begin
                  self->dm_choose_trajdata,t1_id,index=tmp,range=[-self.cut_t1ran[1],-self.cut_t1ran[0]],empty=empty
                  if (empty eq 0) and (n_elements(index1) ne 0) then begin
                     if index1[0] eq -1 then index1 = tmp $
                     else index1 = [index1,tmp]
                  endif
               endif
            endif else if x_id eq t2_id then begin
               if total(finite(self.cut_t2ran)) ne 0 then begin
                  self->dm_choose_trajdata,t2_id,index=tmp,range=[-self.cut_t2ran[1],-self.cut_t2ran[0]],empty=empty
                  if (empty eq 0) and (n_elements(index2) ne 0) then begin
                     if index2[0] eq -1 then index2 = tmp $
                     else index2 = [index2,tmp]
                  endif
               endif
            endif
         endif
         if ok_ewid then begin ;figure out dE
            if self.cut_x[3] eq 2 then dE=(finite(self.cut_step)?(self.cut_step):dE) $
            else if ((self.cut_x[3]+1) mod 3) eq 2 then begin
               if total(finite(self.cut_t1ran)) eq 2 then $
                  if self.cut_t1ran[1] gt self.cut_t1ran[0] then $
                     dE = self.cut_t1ran[1]-self.cut_t1ran[0]
            endif else if ((self.cut_x[3]+2) mod 3) eq 2 then begin
               if total(finite(self.cut_t2ran)) eq 2 then $
                  if self.cut_t2ran[1] gt self.cut_t2ran[0] then $
                     dE = self.cut_t2ran[1]-self.cut_t2ran[0]
            endif
         endif
         end
       else:
    endcase

    ;check extra hkl range
    n_hklran = total(finite(self.disp_hklran))
    if (self.samp_typ gt 0) and (self.hklran_yn eq 1) and (n_hklran ne 0) then begin
       self->dm_choose_trajdata,self.disp_hkl,index=index3,range=self.disp_hklran,empty=empty,debug=debug
       if empty then begin
          index3 = -1 & tmp = temporary(index3)
          if n_elements(emptyitem) eq 0 then emptyitem = (['H','K','L'])[self.disp_hkl] $
          else emptyitem = emptyitem+', '+(['H','K','L'])[self.disp_hkl]
       endif else if n_elements(index3) eq ndat then begin
          index3 = -1 & tmp = temporary(index3)
       endif
    endif

    if n_elements(emptyitem) ne 0 then $
    ok = dialog_message(['The specified '+emptyitem+' range covers no data.', $
               'The entire possible range is used instead.'],dialog_parent=self.tlb)

    if n_elements(index1) ne 0 then index = temporary(index1)
    if n_elements(index2) ne 0 then begin
       if n_elements(index) ne 0 then index = dm_common(index,index2,/no_copy) $
       else index = temporary(index2)
    endif
    if n_elements(index3) ne 0 then begin
       if n_elements(index) ne 0 then index = dm_common(index,index3,/no_copy) $
       else index = temporary(index3)
    endif

    ;ydat
    self->dm_choose_trajdata, self.cut_y[self.samp_typ]-1,data=ydat,error=yerr,label=ytit
    ok_err = (n_elements(yerr) ne 0)

    if n_elements(index) ne 0 then begin
       if index[0] ne -1 then begin
          xdat = xdat[index]
          ydat = ydat[index]
          if ok_err  then yerr = yerr[index]
          if ok_ewid then ewid = ewid[index]
       endif
       index = -1
    endif

    ;xrange
    if total(finite(self.cut_ran)) eq 2 then begin
       if self.cut_ran[1] gt self.cut_ran[0] then begin
          index = where((xdat le self.cut_ran[1]) and (xdat ge self.cut_ran[0]),count)
          if count gt 0 then begin
             xdat = xdat[index]
             ydat = ydat[index]
             if ok_err  then yerr = yerr[index]
             if ok_ewid then ewid = ewid[index]
             index = -1
          endif
       endif
    endif else begin
       if finite(self.cut_ran[0]) then begin
          index = where(xdat ge self.cut_ran[0],count)
          if count gt 0 then begin
             xdat = xdat[index]
             ydat = ydat[index]
             if ok_err  then yerr = yerr[index]
             if ok_ewid then ewid = ewid[index]
             index = -1
          endif
       endif
       if finite(self.cut_ran[1]) then begin
          index = where(xdat le self.cut_ran[1],count)
          if count gt 0 then begin
             xdat = xdat[index]
             ydat = ydat[index]
             if ok_err  then yerr = yerr[index]
             if ok_ewid then ewid = ewid[index]
             index = -1
          endif
       endif
    endelse

    if keyword_set(debug) then help,xdat,ydat,ewid

    rebinned = 0b
    ;rebinning according to the steps
    if finite(self.cut_step) then begin
       if (self.cut_step ne 0) then begin
          dm_step_bin,self.cut_step,xdat,ydat=ydat,yerr=yerr,debug=debug,avgsum=self.bin_avgsum,$
             conststep=self.bin_conststep,ewid=ewid
          rebinned = 1b
       endif
    endif
    if rebinned eq 0 then begin
       ;rebin if x-axis is E
       if ((xtit eq 'E (meV)') or (xtit eq 'Energy')) and (finite(self.cut_step) eq 0) then $
          dm_step_bin,1e-5,xdat,ydat=ydat,yerr=yerr,uniq_val=uniq_E,debug=debug,avgsum=self.bin_avgsum,$
             ewid=ewid
    endif

    if ok_ewid then begin
       ydat = ydat*dE
       if ok_err then yerr = yerr*dE
    endif

    if keyword_set(debug) then print,'^^choose cut data finished in ',systime(/sec)-current,' sec.'
end