; $Id: dm_calc_deteff.pro,v 1.7 2007/03/15 13:57:17 ymqiu Exp $
;#######################################################################
;
; NAME:
;  dm_calc_deteff
;
; PURPOSE:
;   calculate detector efficiency
;
; CATEGORY:
;  dcs_mslice
;
; AUTHOR:
;  Yiming Qiu
;  NIST Center for Neutron Research
;  100 Bureau Drive, Gaithersburg, MD 20899-8562
;  United States
;  yiming.qiu@nist.gov
;  March, 2007
;
; 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.
;
;#######################################################################

; parameter:
;   qty:  [ne,ndet]
; keyword:
;   fit:  if set, fit by an gaussian and use the data in the gaussian fucntion range
function dm_calc_deteff,qty,fit=fit,debug=debug
    if keyword_set(fit) then begin
       ;first fit the det-sum and find out the elastic range
       y    = total(qty,2,/double)
       ny   = n_elements(y)
       if keyword_set(debug) then print,'n=',ny
       x    = findgen(ny)
       ymax = max(y,imax) &imax=imax[0]
       ymin = min(y)
       ;shift the peak approcimately to the center channel
       y    = shift(y,ny/2-imax)
       fwhm = total(y-ymin)/(ymax-ymin)
       estimates = [ymax,ny/2,fwhm/sqrt(8.*alog(2)),ymin]
       estimates[2] = estimates[2]/2.
       if keyword_set(debug) then print,'estimates=',estimates
       result = gaussfit(x,y,params,estimates=estimates,nterms=n_elements(estimates))
       if keyword_set(debug) then print,'parameter=',params
       tmp  = shift(qty,ny/2-imax,0)
       params[2] = abs(params[2])
       istart = (0>(ny/2-fix(2.5*params[2])))
       iend = ((ny-1)<(ny/2+fix(2.5*params[2])))
       tmp  = tmp[istart:iend,*]
       eff  = total(tmp,1,/double)
       if keyword_set(debug) then begin
         window,0
         plot,y,psym=6,xr=[istart,iend]
         oplot,result
       endif
    endif else begin
       eff  = total(qty,1,/double)
    endelse
    return,float(eff/mean(eff,/double))
end