; $Id: dm_calc_projection.pro,v 1.45 2008/04/02 17:05:53 ymqiu Exp $
;#######################################################################
;
; NAME:
;  dm_calc_projection
;
; PURPOSE:
;  calculate projection
;  apply masking, rebinning, and detector efficiency adjustment to 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
;  March, 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 dcs_mslice::dm_calc_projection,debug=debug
    self.error = 1  ;initialize the error flag, it will be cleared at the end of the program
    if keyword_set(debug) then begin
       print,'**Calculate projection...'
       current=systime(/sec)
    end

    det_thick = 1.0     ;DCS He3 detector thickness in cm
    det_p     = 6.0     ;DCS He3 detector pressure in bar
    is_spe    = 0b      ;type of psi defintion

    case self.ftype of
        0:  begin   ;dcs
            if not ptr_valid(self.detPosPtr) then self.detPosPtr=dm_load_detpos(/dcs)
            end
        1:  begin   ;spe
            is_spe = 1b      ;spe type detector psi definition
            if not ptr_valid(self.detPosPtr) then self->dm_load_phx,/phx
            end
        2:  begin   ;inx
            if not ptr_valid(self.detPosPtr) then self.detPosPtr=dm_load_detpos(/neat)
            end
        3:  begin   ;macs
            if not ptr_valid(self.detPosPtr) then self.detPosPtr=dm_load_detpos(/macs)
            end
        else:
    endcase

    if not ptr_valid(self.detPosPtr) then begin
       ok = dialog_message(['Detector information is not available.','Projection calculation is aborted.'],$
                /error,dialog_parent=self.tlb)
       return
    endif

    ;check if the common input parameters are all set for single crystal and diffuse scattering type
    if self.samp_typ ne 0 then begin
       if total(finite(self.latt_parm)) ne 6 then begin
          ok = dialog_message('Missing lattice parameters.',/error,dialog_parent=self.tlb)
          return
       endif
       if total([finite(self.orie_u),finite(self.orie_v),finite(self.orie_psi)]) ne 7 then begin
          ok = dialog_message('Missing crystal orientation parameters.',/error,dialog_parent=self.tlb)
          return
       endif
       if self.samp_typ eq 1 then begin
          if total([finite(self.view_u1),finite(self.view_u2),finite(self.view_u3)]) ne 12 then begin
             ok = dialog_message('Missing viewing axes parameters.',/error,dialog_parent=self.tlb)
             return
          endif
       endif else begin
          if total([finite(self.view_u1[0:2]),finite(self.view_u2[0:2])]) ne 6  then begin
             ok = dialog_message('Missing viewing axes parameters.',/error,dialog_parent=self.tlb)
             return
          endif
       endelse
    endif

    if not ptr_valid(self.dataStrPtr) then begin    ;no data loaded yet
       ;ask for Ei or Ef according to instrument geometry
       eief = 5.0
       if (finite(self.eief)) then $
          if self.eief gt 0 then eief=self.eief
       if self.samp_typ eq 2 then begin             ;diffuse scattering single E
          case self.e_bin[self.samp_typ] of
            0:  en = 0.0                                      ;use elastic energy for all time channels
            1:  en = 0.0                                      ;use elastic energy for elastic peak
            2:  en = (self.e_range[0]+self.e_range[1])/2.0    ;use the mean energy for specified energy range
            else: en = 0.0
          endcase
          info=['E='+dm_to_string(en)+' meV']
       endif else if self.samp_typ eq 3 then begin  ;diffuse scattering multiple E
          case self.e_bin[3] of
            0:  info = (['-Ei< E <Ei','-Ef< E <Ef'])[self.instrgeom]
            1:  info = ['E~0 meV']
            2:  info = [dm_to_string(self.e_range[0])+'< E <'+dm_to_string(self.e_range[1])]
            else:
          endcase
       endif else info=(['-Ei< E <Ei','-Ef< E <Ef'])[self.instrgeom]
       if self.ftype eq 0 then begin    ;for DCS convert energy(meV) to wavelength(angstrom)
          eief  = sqrt(81.8042/eief)
          title = "Incident wavelength:"
       endif else $
          title = (["Incident energy:","Final energy:"])[self.instrgeom]
       notok = 1
       while(notok) do begin
          if self.ftype eq 0 then begin
             eief = dm_dialog_input('wavelength ('+string('c5'XB)+'):',/float,default=fix(eief*100.0+0.5,type=3)/100.0,$
                    dialog_parent=self.tlb,info=info,/lambda2e,title=title)
          endif else begin
             eief = dm_dialog_input((['Ei','Ef'])[self.instrgeom]+'(meV):',/float,title=title,$
                    default=fix(eief*1000.0+0.5,type=3)/1000.0,dialog_parent=self.tlb,info=info,/e2lambda)
          endelse
          if (finite(eief)) then $
             if eief gt 0 then notok = 0
       endwhile
       if self.ftype eq 0 then $    ;for DCS convert wavelength to energy
          eief =  81.8042/(eief^2)
       self.eief = eief     ;save ei or ef information
       if self.samp_typ eq 2 then begin             ;diffuse scattering single E
          if (self.instrgeom eq 0) and (en gt eief) then begin
             ok = dialog_message('Incident energy is less than the energy transfer: E>Ei. E=0 is used.',/error)
             en = 0.0
          endif
          if (self.instrgeom eq 1) and (en lt -eief) then begin
             ok = dialog_message('Neutron energy gain is more than the final energy: E<-Ef. E=0 is used.',/error)
             en = 0.0
          endif
       endif else if self.samp_typ eq 3 then begin  ;diffuse scattering multiple E
          case self.e_bin[3] of
            0:  en = (findgen(41)/10.-1)*eief+(2*self.instrgeom-1)*1e-6         ;[-Ei,...,Ei] or [-Ef,...Ef]
            1:  begin
                ;use empirical DCS energy resolution as the energy range
                dE = 0.02008*eief^(1.4735)
                en = (findgen(41)/10.-1)*dE
                end
            2:  begin
                if self.instrgeom eq 0 then begin
                   if min(self.e_range) gt eief then begin
                      ok = dialog_message('Incident energy is smaller than the specified energy range. Elastic peak is used.',/error)
                      ;use empirical DCS energy resolution as the energy range
                      dE = 0.02008*eief^(1.4735)
                      en = (findgen(41)/20.-1)*dE
                   endif else begin
                      eran = [self.e_range[0],min([self.e_range[1],eief])]
                      en = eran[0] + findgen(41)/40*(eran[1]-eran[0])-1e-6
                   endelse
                endif else begin
                   if max(self.e_range) lt -eief then begin
                      ok = dialog_message('Neutron energy gain in the specified energy range is more than the final energy. Elastic peak is used.',/error)
                      ;use empirical DCS energy resolution as the energy range
                      dE = 0.02008*eief^(1.4735)
                      en = (findgen(41)/20.-1)*dE
                   endif else begin
                      eran = [max(self.e_range[0],-eief),self.e_range[1]]
                      en = eran[0] + findgen(41)/40*(eran[1]-eran[0])+1e-6
                   endelse
                endelse
                end
            else:
          endcase
       endif else en = (findgen(81)/40.-1)*eief +(2*self.instrgeom-1)*1e-6         ;[-Ei,...,Ei] or [-Ef,...Ef]
       if self.samp_typ ge 2 then begin ;A2 range for diffuse scattering type
          notok = 1
          if total(finite(self.a2[0:2])) eq 3 then a2=self.a2[0:2] else a2 = [0,50,2]
          while(notok) do begin
             a2  = dm_dialog_input(['Start:','End:','Step:'],title='Psi range:',default=a2,/float,dialog_parent=self.tlb)
             if (total(finite(a2)) eq 3) then $
                if a2[2] ne 0 then notok = 0
          endwhile
          a2[2]  = abs(a2[2])
          min_a  = min(a2[0:1],max=max_a)
          self.a2[0:2] = [min_a,max_a,a2[2]]
          a2     = min_a+findgen(fix((max_a-min_a)/a2[2])+1)*a2[2]
          n_file = n_elements(a2)
       endif
       qty  = 0.
       err  = 0.
       ewid = 0.
       time = 0.
    endif else begin
       ewid = (*self.dataStrPtr).ewid
       eief = (*self.dataStrPtr).eief & self.eief = eief
       qty  = ((*self.dataStrPtr).qty)
       err  = ((*self.dataStrPtr).err)
       while (not finite(eief)) do begin
            eief = dm_dialog_input((['Ei','Ef'])[self.instrgeom]+'(meV):',/float,default=eief,dialog_parent=self.tlb,/e2lambda)
            (*self.dataStrPtr).eief = eief
            self.eief = eief
       endwhile
       en   = (*self.dataStrPtr).e
       time = (*self.dataStrPtr).time
       if self.samp_typ ge 2 then begin                    ;diffuse scattering type
          a2 = (*self.dataStrPtr).rotation                 ;angle 2-rotation
          self.orie_offset = self.orie_psi-self.a2[3]      ;save the offset information
          a2 = a2+self.orie_offset                         ;adjust the offset
          n_file  = n_elements(a2)
          tmpmin  = min(a2,max=tmpmax)
          tmpstep = (n_file gt 1)?((tmpmax-tmpmin)/(n_file-1)):(0)
          self.a2[0:2] = [tmpmin,tmpmax,tmpstep]
       endif else if self.samp_typ eq 1 then begin
          if finite((*self.dataStrPtr).rotation) then $
             self.orie_offset = self.orie_psi-(*self.dataStrPtr).rotation
       endif
       if ptr_valid(self.bgdata) then begin
          bgqty = (*self.bgdata).qty
          bgerr = (*self.bgdata).err
          bgang = (*self.bgdata).rotation
          if self.samp_typ ge 2 then n_bgfile = n_elements(bgang)
       endif
    endelse
    det_tt   = (*self.detPosPtr).two_theta
    det_psi  = (*self.detPosPtr).psi
    det_dtt  = (*self.detPosPtr).dtwo_theta
    det_dpsi = (*self.detPosPtr).dpsi
    nen      = n_elements(en)
    ndet     = n_elements(det_tt)
    widget_control,/hourglass
    if ptr_valid(self.dataStrPtr) then begin
       ;subtract overall background rate
       if ptr_valid(self.bgratePtr) then begin
          bgrate = (*self.bgratePtr)
          nrate  = n_elements(bgrate)
          if nrate gt 1 then begin
             dbgrate = bgrate[*,1]
             bgrate  = bgrate[*,0]
          endif else $
             dbgrate = 1.0      ;arbitrary error bar
          case self.ftype of
            0:  begin   ;dcs
                if nrate eq 1 then begin    ;overall rate
                   if ptr_valid(self.eff) then tmp_eff=(*self.eff) else tmp_eff=fltarr(ndet)+1.
                   tmp_value = bgrate/3600./(self.monrate)*temporary(tmp_eff)
                endif else $
                   tmp_value = bgrate/3600./(self.monrate)
                if self.samp_typ ge 2 then begin      ;diffuse scattering type
                   if self.samp_typ eq 2 then begin   ;single energy
                      if self.eint_yn then $
                         tmp_value = tmp_value*total(ewid) $
                      else $
                         tmp_value = tmp_value
                      for i=0L,n_file-1 do qty[*,i]=qty[*,i]-tmp_value
                      if ptr_valid(self.bgdata) then $
                         for i=0L,n_bgfile-1 do bgqty[*,i]=bgqty[*,i]-tmp_value
                   endif else begin                    ;multiple energy
                      for i=0L,n_file-1L do begin
                          for j=0L,nen-1L do $
                              qty[j,*,i] = qty[j,*,i]-tmp_value
                      endfor
                      if ptr_valid(self.bgdata) then begin
                         for i=0L,n_bgfile-1L do begin
                             for j=0L,nen-1L do $
                                 bgqty[j,*,i] = bgqty[j,*,i]-tmp_value
                         endfor
                      endif
                   endelse
                endif else begin
                   for i=0L,nen-1L do $
                       qty[i,*] = qty[i,*]-tmp_value
                   if ptr_valid(self.bgdata) then $
                      for i=0L,nen-1L do $
                          bgqty[i,*] = bgqty[i,*]-tmp_value
                endelse
                tmp_value = 1     ;free up the memory
                end
            1:  begin   ;spe
                ;no implemented
                end
            2:  begin   ;inx
                ;no implemented
                end
            else:
          endcase
       endif
       ;subtract detector dependent background data
       if ptr_valid(self.bgdetsPtr) then begin
          tmp_value = (*self.bgdetsPtr)
          if self.samp_typ ge 2 then begin
             if self.eint_yn then $
                tmp_value = tmp_value*total(ewid)
             if self.samp_typ eq 2 then begin
                for i=0L,n_file-1L do $
                    qty[*,i] = qty[*,i]-tmp_value
             endif else begin
                tmp_value = (fltarr(nen)+1.0)#tmp_value
                for i=0L,n_file-1L do $
                    qty[*,*,i] = qty[*,*,i]-tmp_value
             endelse
          endif else begin
             qty = qty-(fltarr(nen)+1.0)#tmp_value
          endelse
          tmp_value = 1     ;free up the memory
       endif
       ;subtract time channel background data
       if ptr_valid(self.bgtchnPtr) then begin
          tmp_emin  = min((*self.bgtchnPtr).emean,max=tmp_emax)
          tmp_nen   = n_elements((*self.bgtchnPtr).emean)
          tchan_ind = where(en ge tmp_emin and en le tmp_emax,tchan_count)
          if tchan_count gt tmp_nen then tchan_ind = tchan_ind[0:(tmp_nen-1)]   ;errorproof
          if tchan_count ne 0 then begin
             if self.samp_typ eq 3 then begin
                for i=0L,n_file-1L do $
                    for j=0L,ndet-1L do $
                        qty[tchan_ind,j,i] = qty[tchan_ind,j,i]-(*self.bgtchnPtr).qty
             endif else begin
                for i=0L,ndet-1L do $
                    qty[tchan_ind,i] = qty[tchan_ind,i]-(*self.bgtchnPtr).qty
             endelse
          endif
       endif else tchan_count = 0
       ;absorption correction
       if ptr_valid(self.absp_info) then begin
          info = [(*(self.absp_info))[0],(*(self.absp_info))[2:(n_elements(*(self.absp_info))-1)]]
          tmp  = dm_calc_shielding(info,eief,en,det_tt,interp=(*(self.absp_info))[1],instrgeom=self.instrgeom,$
                 group_leader=self.tlb)
          if self.samp_typ ge 2 then begin    ;diffuse scattering type
             if self.samp_typ eq 2 then begin
                for i=0L,n_file-1L do begin
                    qty[*,i] = qty[*,i]*tmp
                    err[*,i] = err[*,i]*tmp
                endfor
             endif else begin
                for i=0L,n_file-1L do begin
                    qty[*,*,i] = qty[*,*,i]*tmp
                    err[*,*,i] = err[*,*,i]*tmp
                endfor
             endelse
          endif else begin
             qty = temporary(qty)*tmp
             err = temporary(err)*tmp
          endelse
          widget_control,/hourglass
          tmp = 1    ;free up the memory
       endif
       ;subtract empty can file
       if ptr_valid(self.bgdata) then begin
          if finite(self.ssfac) then ssfac=self.ssfac $
          else ssfac = 1.0     ;self-shielding factor
          if self.samp_typ ge 2 then begin    ;diffuse scattering type
             if n_bgfile eq 1 then begin
                if self.samp_typ eq 2 then begin
                   for i=0L,n_file-1L do begin
                       qty[*,i] = qty[*,i]-bgqty[*,0]*ssfac
                       err[*,i] = sqrt(err[*,i]^2+(bgerr[*,0])^2*(ssfac*ssfac))
                   endfor
                endif else begin
                   for i=0L,n_file-1L do begin
                       qty[*,*,i] = qty[*,*,i]-bgqty[*,*,0]*ssfac
                       err[*,*,i] = sqrt(err[*,*,i]^2+(bgerr[*,*,0])^2*(ssfac*ssfac))
                       if tchan_count ne 0 then begin   ;avoid double-subtraction
                          qty[tchan_ind,*,i] = qty[tchan_ind,*,i]+bgqty[tchan_ind,*,0]*ssfac
                          ;error bar is not changed
                       endif
                   endfor
                endelse
             endif else begin
                isort = sort(bgang)
                bgang = bgang[isort]
                if self.samp_typ eq 2 then begin
                   for i=0L,n_file-1L do begin
                       index    = 0>(value_locate(bgang,(*self.dataStrPtr).rotation[i]))
                       qty[*,i] = qty[*,i]-bgqty[*,isort[index]]*ssfac
                       err[*,i] = sqrt(err[*,i]^2+(bgerr[*,isort[index]])^2*(ssfac*ssfac))
                   endfor
                endif else begin
                   for i=0L,n_file-1L do begin
                       index      = 0>(value_locate(bgang,(*self.dataStrPtr).rotation[i]))
                       qty[*,*,i] = qty[*,*,i]-bgqty[*,*,isort[index]]*ssfac
                       err[*,*,i] = sqrt(err[*,*,i]^2+(bgerr[*,*,isort[index]])^2*(ssfac*ssfac))
                       if tchan_count ne 0 then begin   ;avoid double-subtraction
                          qty[tchan_ind,*,i] = qty[tchan_ind,*,i]+bgqty[tchan_ind,*,isort[index]]*ssfac
                          ;error bar is not changed
                       endif
                   endfor
                endelse
             endelse
          endif else begin
             qty = temporary(qty)-(bgqty)*ssfac
             err = sqrt(temporary(err)^2+(bgerr)^2*(ssfac*ssfac))
             if tchan_count ne 0 then begin
                qty[tchan_ind,*] = qty[tchan_ind,*]+(bgqty[tchan_ind,*])*ssfac
                ;error bar is not changed
             endif
          endelse
          bgqty = 0 & bgerr = 0    ;free up the memory
       endif
       ;efficiency adjustment
       if ptr_valid(self.eff) then begin
          if self.samp_typ ge 2 then begin    ;diffuse scattering type
             if self.samp_typ eq 2 then begin
                for i=0L,n_file-1L do begin
                    qty[*,i] = qty[*,i]/(*self.eff)
                    err[*,i] = err[*,i]/(*self.eff)
                endfor
             endif else begin
                for i=0L,n_file-1L do begin
                    for j=0L,nen-1L do begin
                        qty[j,*,i] = qty[j,*,i]/(*self.eff)
                        err[j,*,i] = err[j,*,i]/(*self.eff)
                    endfor
                endfor
             endelse
          endif else begin
             for i=0L,nen-1L do begin
                 qty[i,*] = qty[i,*]/(*self.eff)
                 err[i,*] = err[i,*]/(*self.eff)
             endfor
          endelse
       endif
       ;kfactor adjustment
       if (self.intn_typ eq 1) or (self.intn_typ eq 3) then begin   ;detailed balance factor
          ;ask for a sample temperature
          old = 300. & t_existed = 0b
          if n_tags((*self.dataStrPtr).info) ne 0 then begin
             tnames = tag_names((*self.dataStrPtr).info)
             index  = where(tnames eq 'TEMPERATURE',count)
             if count ne 0 then begin
                old = (*self.dataStrPtr).info.temperature
                t_existed = 1b
             endif
          endif
          notvalid = 1
          while notvalid do begin
             temp = dm_dialog_input('sample temperature:',/float,default=old,dialog_parent=self.tlb,$
                    title='Detailed blance calculation')
             widget_control,/hourglass
             if finite(temp) then $
                if temp gt 0 then begin
                   if t_existed then (*self.dataStrPtr).info.temperature=temp
                   notvalid = 0 ;valid temperature given
                endif
          endwhile
          ;calculate detailed balance factor
          kb   = 8.617342e-2    ;Boltzmann constant in mev.k^-1
          dbal = exp(-en/(2.0*kb*temp))
       endif
       case self.intn_typ of
          0:  tmp = (*self.kfactor)^4               ;S(Q,omega)
          1:  tmp = (*self.kfactor)^4*dbal          ;Symmetrized S(Q,omega)
          2:  tmp = (*self.kfactor)^3               ;d2sigma/dOmegadEf
          3:  tmp = (*self.kfactor)^4*(1-dbal^2)    ;Chi(Q,omega)
          4:  tmp = 1.0                             ;no kfactor adjustment
         else:tmp = 1.0                             ;no kfactor adjustment
       endcase
       if self.eint_yn and (self.samp_typ ne 2) then $
          tmp = temporary(tmp)*ewid
       ;energy dependent detector efficiency adjustment
       if (self.edep_yn eq 1) and (self.instrgeom eq 0) then begin
          sigma_acoeff = 0.4466*det_p/6.0  ; (from Jeremy in cm-1/lambda_f) absorption factor of He3 gas
          sigma_s      = 0.0008*det_p/6.0  ; (from Jeremy in cm-1) sigma_s+
          wl0   = sqrt(81.8042/(*self.dataStrPtr).eief) ;incident wavelength
          sumtd = sigma_acoeff*wl0*(*self.kfactor)+sigma_s
          ;deteff_inv=1.0/(1.0-exp(-det_thick*sumtd))
          deteff_inv = sumtd/(sumtd-sigma_s)/(1.0-exp(-det_thick*sumtd))    ;from Jeremy's Equation 9
          sumtd = 0
          tmp   = temporary(tmp)*temporary(deteff_inv)
       endif
       if self.samp_typ ge 2 then begin
          if self.samp_typ eq 2 then begin  ;tmp is a scalor here, qty is a [ndet,nfile] arrray
             qty = temporary(qty)*tmp
             err = temporary(err)*tmp
          endif else begin                  ;tmp is a [nen] array, qty is a [nen,ndet,nfile] array
             for i=0L,ndet-1L do begin
                 for j=0L,n_file-1L do begin
                     qty[*,i,j] = qty[*,i,j]*tmp
                     err[*,i,j] = err[*,i,j]*tmp
                 endfor
             endfor
          endelse
       endif else begin
          for i=0L,ndet-1L do begin         ;tmp is a [nen] array, qty is a [nen,ndet] array
              qty[*,i] = qty[*,i]*tmp
              err[*,i] = err[*,i]*tmp
          endfor
       endelse
       tmp = 0    ;free up the memory

       ;energy rebin
       if (self.e_bin[self.samp_typ] gt 1) and (self.samp_typ lt 2) then begin
          dm_grid_ebin,self.e_bin[self.samp_typ],qty,err,en,debug=debug,avgsum=(self.bin_avgsum eq 1)
          dm_grid_ebin,self.e_bin[self.samp_typ],ewid
          if (self.samp_typ eq 0) and (self.ftype eq 0) then $
             dm_grid_ebin,self.e_bin[self.samp_typ],time
       endif

       ;overall intensity factor correction
       if self.samp_typ ge 2 then begin
          ;calculate dark angle correction factor and save as overall intensity factor
          if self.dcor_yn and ptr_valid(self.dcor_info) then begin
             tmp_qty = (*self.dcor_info).factor
             tmp_ang = (*self.dcor_info).a2
             isort   = sort(tmp_ang)
             tmp_ang = tmp_ang[isort]
             tmp_qty = tmp_qty[isort]
             tmp_fac = fltarr(n_file)
             for i=0L,n_file-1L do begin
                 index      = (0L)>(value_locate(tmp_ang,(*self.dataStrPtr).rotation[i]))
                 tmp_fac[i] = tmp_qty[index]
             endfor
             (*self.oifac) = temporary(tmp_fac)
          endif
          if self.samp_typ eq 2 then begin
             for i=0L,n_file-1L do begin
                 qty[*,i] = qty[*,i]*((*self.oifac)[i])
                 err[*,i] = err[*,i]*((*self.oifac)[i])
             endfor
          endif else begin
             for i=0L,n_file-1L do begin
                 qty[*,*,i] = qty[*,*,i]*((*self.oifac)[i])
                 err[*,*,i] = err[*,*,i]*((*self.oifac)[i])
             endfor
          endelse
       endif else begin
          qty = qty*(*self.oifac)
          err = err*(*self.oifac)
       endelse
    endif

    ;mask detectors
    if ptr_valid(self.mask) then begin
       tmp = bytarr(ndet)
       tmp[*self.mask] = 1b
       ind = where(tmp eq 0,count) & tmp=0
       if count ne 0 then begin
          if ptr_valid(self.dataStrPtr) then begin
             if self.samp_typ ge 2 then begin
                if self.samp_typ eq 2 then begin
                   qty = qty[ind,*]
                   err = err[ind,*]
                endif else begin
                   qty = qty[*,ind,*]
                   err = err[*,ind,*]
                endelse
             endif else begin
                qty = qty[*,ind]
                err = err[*,ind]
             endelse
          endif
          det_tt   = det_tt[ind]
          det_dtt  = det_dtt[ind]
          det_psi  = det_psi[ind]
          det_dpsi = det_dpsi[ind]
          ndet     = count
       endif
    endif

    if self.samp_typ ge 1 then begin
       ;direct -> reciprocal basis vectors
       abc_r = dm_basis_r(self.latt_parm,error=error)
       if error then return  ;error encountered, no further calcuation
       ;orthogonalize and normalize the crystal orientation axes
       if norm(crossp(self.orie_u,self.orie_v)) eq 0 then begin    ;parallel
          ok = dialog_message('Please check cyrstal orientation.',/error,dialog_parent=self.tlb)
          return          ;no further calculation
       endif
       uvw      = dblarr(3,3)
       uvw[*,0] = abc_r#self.orie_u
       uvw[*,1] = abc_r#self.orie_v
       uvw[*,0] = uvw[*,0]/norm(uvw[*,0],/double)
       uvw[*,1] = uvw[*,1]-uvw[*,0]#(transpose(uvw[*,1])#uvw[*,0])
       uvw[*,1] = uvw[*,1]/norm(uvw[*,1],/double)
       for i=0,2 do uvw[i,2]=uvw[(i+1) mod 3,0]*uvw[(i+2) mod 3,1]-uvw[(i+2) mod 3,0]*uvw[(i+1) mod 3,1]
                 ;cross product, not using crossp() for higher precision
       ;rewrite reciprocal lattice parameters in the (u,v,w) frame
       abc_r = transpose(uvw)#(temporary(abc_r))
    endif

    ptr_free,self.projStrPtr,self.powd_line

    case self.samp_typ of
     0:begin  ;powder
       ;combine equivalent detectors
       if self.negative2th eq 0 then det_tt=abs(temporary(det_tt))
       ind      = sort(det_tt)
       ndat     = n_elements(det_tt)
       if ptr_valid(self.dataStrPtr) then begin
          qty   = qty[*,ind]
          err   = err[*,ind]
       endif
       det_tt   = det_tt[ind]
       det_dtt  = det_dtt[ind]
       det_psi  = det_psi[ind]
       det_dpsi = det_dpsi[ind]
       ind      = uniq(det_tt)
       nun      = n_elements(ind)
       if nun ne ndat then begin
          j = 0L
          for i=0L,nun-1 do begin
              if ptr_valid(self.dataStrPtr) then begin
                 qtmp = 0.0
                 etmp = 0.0
              endif
              itmp = 0L
              while(det_tt[j] eq det_tt[ind[i]]) do begin
                 if det_dtt[i] lt det_dtt[j] then det_dtt[i]=det_dtt[j]
                 if det_dpsi[i] lt det_dpsi[j] then det_dpsi[i]=det_dpsi[j]
                 if ptr_valid(self.dataStrPtr) then begin
                    qtmp = temporary(qtmp)+qty[*,j]
                    etmp = temporary(etmp)+err[*,j]*err[*,j]
                    itmp = itmp+1L
                 endif
                 j = j+1L
                 if j eq ndat then break
              endwhile
              if ptr_valid(self.dataStrPtr) then begin
                 qty[*,i] = qtmp/(itmp>1)
                 err[*,i] = sqrt(etmp)/(itmp>1)
              endif
          endfor
          if ptr_valid(self.dataStrPtr) then begin
             qtmp  = 0 & etmp = 0  ;free up the memory
             qty   = qty[*,0:nun-1]
             err   = err[*,0:nun-1]
          endif
          det_dtt  = det_dtt[0:nun-1]
          det_dpsi = det_dpsi[0:nun-1]
          det_tt   = det_tt[ind]
          det_psi  = det_psi[ind]
       endif
       Q = dm_proj_spec_samp(eief,en,det_tt,det_dtt,det_psi,det_dpsi,type=0,debug=debug,instrgeom=self.instrgeom)   ;Q[nen,ndet]
       self.projStrPtr=ptr_new({Q:temporary(Q),en:temporary(en),two_theta:temporary(det_tt),$
          qty:temporary(qty),err:temporary(err),ewid:temporary(ewid),time:temporary(time)})
       end
     1:begin  ;single crystal
       ;reconstruct the viewing axis
       U1 = [abc_r#self.view_u1[0:2],self.view_u1[3]]
       U2 = [abc_r#self.view_u2[0:2],self.view_u2[3]]
       U3 = [abc_r#self.view_u3[0:2],self.view_u3[3]]
       norm1=norm(U1,/double) & norm2=norm(U2,/double) & norm3=norm(U3,/double)
       if self.view_ortharbi eq 0 then begin    ;check orthogonality
          if max(abs([transpose(U1)#U2/norm1/norm2,transpose(U2)#U3/norm2/norm3,$
          transpose(U3)#U1/norm3/norm1])) gt 1e-8 then begin   ;not orthogonal
             ok = dialog_message('Viewing axes are not orthogonal.',/error,dialog_parent=self.tlb)
             return          ;no further calculation
          endif
       endif else begin
          ;make sure no viewing axes are parallel to each other
          tmp = [[self.view_u1[0:2]],[self.view_u2[0:2]],[self.view_u3[0:2]]]
          for i=0,2 do begin
              if total(abs(tmp[*,i])) eq 0 then begin
                 tmp[*,i]=crossp(tmp[*,(i+1) mod 3],tmp[*,(i+2) mod 3])
                 if norm(tmp[*,i]) eq 0 then begin
                    ok = dialog_message('Parallel viewing axes',/error,dialog_parent=self.tlb)
                    return  ;no further calculation
                 endif
                 break
              endif
          endfor
       endelse
       Q   = dm_proj_spec_samp(eief,en,det_tt,det_dtt,det_psi,det_dpsi,psi=self.orie_psi,type=1,debug=debug,is_spe=is_spe,instrgeom=self.instrgeom)
                                                          ;Q(3,ne,ndet) in sample frame
       V   = dm_proj_view(Q,en,U1,U2,U3,debug=debug)      ;project onto viewing axis
       hkl = dm_q2hkl(Q,abc_r,debug=debug)                ;hkl values
       self.projStrPtr=ptr_new({V:temporary(V),hkl:temporary(hkl),en:temporary(en),$
         two_theta:temporary(det_tt),qty:temporary(qty),err:temporary(err),ewid:temporary(ewid)})
       end
     else:begin  ;diffuse scattering
       ;reconstruct the viewing axis and check orthogonality
       U1 = [abc_r#self.view_u1[0:2]]
       U2 = [abc_r#self.view_u2[0:2]]
       norm1 = norm(U1,/double) & norm2 = norm(U2,/double)
       if self.view_ortharbi eq 0 then begin    ;check orthogonality
          if abs(transpose(U1)#U2/norm1/norm2) gt 1e-8 then begin      ;not orthogonal
             ok = dialog_message('Viewing axes are not orthogonal.',/error,dialog_parent=self.tlb)
             return          ;no further calculation
          endif
       endif else begin
          ;make sure the viewing axes are not parallel to each other
          if norm(crossp(self.view_u1[0:2],self.view_u2[0:2])) eq 0 then begin
             ok = dialog_message('Parallel viewing axes',/error,dialog_parent=self.tlb)
             return  ;no further calculation
          endif
       endelse
       Q   = dm_proj_spec_samp(eief,en,det_tt,det_dtt,det_psi,det_dpsi,psi=a2,type=2,debug=debug,is_spe=is_spe,diffetyp=(self.samp_typ eq 3),instrgeom=self.instrgeom)
                                                                            ;Q[3,ndet,nfile] in sample frame
       V   = dm_proj_view(Q,en,U1,U2,[0.0,0.0,0.0],/diffuse,debug=debug)    ;project onto viewing axis
       tmp = size(V,/dimension)
       case n_elements(tmp) of  ;throw away the third projection
            2:    V = V[0:1,*]
            3:    V = V[0:1,*,*]
            4:    V = V[0:1,*,*,*]
            else: V = V[0:1,*]
       endcase
       hkl = dm_q2hkl(Q,abc_r,debug=debug)                                  ;hkl values
       self.projStrPtr = ptr_new({V:temporary(V),hkl:temporary(hkl),en:temporary(en),$
                two_theta:temporary(det_tt),qty:temporary(qty),err:temporary(err)})
       ;Al, Cu, or Stainless Steel powder line calculation
       if (self.samp_typ eq 2) or (self.disp_perp eq 2) or (not ptr_valid(self.dataStrPtr)) then begin
          case self.powd_typ of
               1: q_value = [sqrt(3),2,sqrt(8),sqrt(11),sqrt(12)]*2*!pi/4.04912  ;Al a=4.04912 FCC
               2: q_value = [sqrt(3),2,sqrt(8),sqrt(11),sqrt(12)]*2*!pi/3.6148   ;Cu a=3.6148  FCC
               3: q_value = [sqrt(3),2,sqrt(8),sqrt(11),sqrt(12)]*2*!pi/3.594    ;Stainless Steel a=3.594 FCC
               else:
          endcase
          tmp = transpose([[cos(findgen(361)*!dtor)],[sin(findgen(361)*!dtor)],[fltarr(361)]])
          for i=0L,n_elements(q_value)-1 do begin
             Vpow = dm_proj_view(q_value[i]*tmp,0,U1,U2,[0.0,0.0,0.0],/diffuse)
             x = reform(Vpow[0,*]) & y = reform(Vpow[1,*])
             if n_elements(lines) ne 0 then lines=[[[lines]],[[[x],[y]]]] $
             else lines=[[x],[y]]
          endfor
          if n_elements(lines) ne 0 then self.powd_line=ptr_new(lines,/no_copy)
       endif
       end
    endcase
    self.error = 0       ;clear the error flag
    ;allow plots
    if not keyword_set(debug) then $
       self->set_display_button $
    else $
       print,'**Projection calculation finished in ',systime(/sec)-current,' sec.'
    if ptr_valid(self.dataStrPtr) and (self.bin_avgsum eq 1) and (self.e_bin[self.samp_typ] gt 1) and (self.samp_typ lt 2)then $
       ok = dialog_message(['The current binning method is SUM.','This message is only a reminder.'],dialog_parent=self.tlb)
end