; $Id: dcs_define_det_eff.pro,v 1.7 2005/02/14 17:50:14 azuah Exp $
; Written by J.R.D. Copley.
;************************************************************************************************
;###############################################################################
;
; 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_define_det_eff,dcs_base,det_eff_qty,det_eff_err,method,nhedet=nhedet,$
	vansumfile=vansumfile,logentry=logentry
;************************************************************************************************
;
compile_opt strictarr
;
;
restore,filename=vansumfile
;
; Much of what follows was lifted from "dcs_normalization_counts".
;
; Define quantities needed to obtain detector peak positions and intensities.
; Master speed, master period, period at sample.
mspeed=(*(*(*davePtr).dataStrPtr).specificPtr).ch_ms
mperiod=60000000.0/mspeed; master period
speriod=mperiod*(*(*(*davePtr).dataStrPtr).specificPtr).ch_srdenom
;
;	Detector distance, incident wavelength, incident inverse velocity, resolution mode
det_dis=(*(*(*davePtr).dataStrPtr).specificPtr).det_dis
lambda=(*(*(*davePtr).dataStrPtr).specificPtr).ch_wl
tau0=lambda/!dcs_hom
resmode=(*(*(*davePtr).dataStrPtr).specificPtr).ch_res
;
;	Chopper distances and slot positions.
distances=(*(*(*davePtr).dataStrPtr).specificPtr).ch_dis
angles=(*(*(*davePtr).dataStrPtr).specificPtr).ch_slots
;
; 3-character STARTCHOICE
startchoice=strmid((*(*(*davePtr).dataStrPtr).specificPtr).startchoice,0,3)
;
; Time channel widths for detectors (in us).
tchanwid=(*(*(*davePtr).dataStrPtr).specificPtr).tchanlook[0,*]*0.05d0
;
;	Start and stop times for each channel (in us). Channel 0 starts at 0.
tstop=reform(total(tchanwid,/cumulative))
tstart=shift(tstop,1)
tstart[0]=0
;
;	Number of time channels.
nch=(*(*(*davePtr).dataStrPtr).specificPtr).ngtchan0
;
;	Truncate tstart and tstop.
tstart=tstart[0:nch-1]
tstop=tstop[0:nch-1]
;
; Centers of time channels (for use in gaussfit)
tmean=0.5*(tstart+tstop); (for use in gaussfit)
;
; Pointer to input data.
indataPtr=ptr_new((*(*(*davePtr).dataStrPtr).commonStr.histPtr).qty[0:nch-1,*])
;
; Method 0 uses the total counts for each detector.
;
if (method eq 0) then begin
	det_eff_qty=total(*indataPtr,1)
	help,det_eff_qty
	norm=mean(det_eff_qty)
	det_eff_qty=det_eff_qty/norm
	ptr_free,indataPtr
	logentry=logentry+string(10B)+$
		"Detector efficiencies computed using the total counts for each detector."
	return
endif
;
; Methods 1 and 2 use elastic peak sums.
;
; Compute ts0, which is the time that the DAQ start pulse occurs,
; 	relative to the time that neutrons cross the sample axis.
;
case startchoice of
	"AND": ts0=(*(*(*davePtr).dataStrPtr).specificPtr).tsdmin
	"CH1": ts0=speriod+tau0*distances[0]-(30-angles[0,resmode-1])*(speriod/360.0)
	"CH2": ts0=speriod+tau0*distances[1]-(30-angles[1,resmode-1])*(speriod/360.0)
	"CH3": ts0=speriod+tau0*distances[2]-(30-angles[2,resmode-1])*(speriod/360.0)
	"CH4": ts0=speriod+tau0*distances[3]-(30-angles[3,resmode-1])*(speriod/360.0)
	"CH5": ts0=speriod+tau0*distances[4]-(30-angles[4,resmode-1])*(speriod/360.0)
	"CH6": ts0=speriod+tau0*distances[5]-(30-angles[5,resmode-1])*(speriod/360.0)
	"CH7": ts0=speriod+tau0*distances[6]-(30-angles[6,resmode-1])*(speriod/360.0)
	else: stop,"Invalid STARTCHOICE"
endcase
;
;
;	Compute elastic time of flight from sample to detector.
t_sd=tau0*det_dis
;
; Hence elastic time of arrival at detector, relative to DAQ start.
;td_the=dave_modulus(t_sd-ts0,speriod) (OLD VERSION)
td_the=(t_sd-ts0) mod speriod + speriod * ((t_sd-ts0) lt 0.0)
;
; Compute corresponding time channel.
tch_the=where(tstart le td_the and td_the lt tstop)
tch_the=tch_the[0]
if (!debug) then print,'Theoretical elastic time channel = ',tch_the
;
; Renumber channels so that peak count is roughly in channel nch/2.
; The number of channels shifted is tch_shift.
;tch_shift=dave_modulus(round(-tch_the+nch/2),nch) (OLD VERSION)
tch_shift=round(-tch_the+nch/2) mod nch + nch * (round(-tch_the+nch/2) lt 0.0)
if (!debug) then print,'tch_shift = ',tch_shift
;
; Rotate data so it peaks in channel nch/2.
*indataPtr=shift(*indataPtr,tch_shift,0)
;
; Redefine mean times for time channels.
tmean_rot=shift(tmean,tch_shift)
addition=fltarr(nch)+speriod
if (tch_shift ne 0) then addition[0:tch_shift-1]=0
tmean_rot=tmean_rot+addition
;
; Compute time shift corresponding to tch_shift.

;tshift=tstop[dave_modulus(nch-1-tch_shift,nch)] (OLD VERSION)
tshift=tstop[(nch-1-tch_shift) mod nch + nch * ((nch-1-tch_shift) lt 0)]
;
; Check whether there have been any floating-point underflows.
if (check_math(mask=32)) then if (!debug) then print, "NB: floating-point underflows"
;
nterms=5
params=fltarr(nterms)
;

; Method 1 uses elastic peak sums based on an overall fit of a Gaussian plus
; constant background to the rotated counts summed over all detectors.
;
; mult determines how far out we go (how many multiples of sigma), either side of
; the peak position; the background is then estimated by giong out a further
; mult*sigma in either direction.
mult=3
if (method eq 1) then begin
	sumdata_rot=total(*indataPtr,2)
	maxcounts=sumdata_rot[round(nch/2)]
	estimates=[maxcounts,tmean_rot[round(nch/2)],20.0,0.0,0.0]
	result=gaussfit(tmean_rot,sumdata_rot,params,estimates=estimates,nterms=nterms)
	lotime=params[1]-mult*params[2]
	loch=where(tmean_rot-10.0 lt lotime and tmean_rot+10.0 ge lotime)
	loch=round(mean(loch))
	hitime=params[1]+mult*params[2]
	hich=where(tmean_rot-10.0 lt hitime and tmean_rot+10.0 ge hitime)
	hich=round(mean(hich))
	nchans=hich-loch+1
	intint=total((*indataPtr)[loch:hich,*],1)
	lo_int=total((*indataPtr)[loch-nchans+1:loch,*],1)
	hi_int=total((*indataPtr)[hich:hich+nchans-1,*],1)
	det_eff_qty=intint-0.5*(lo_int+hi_int)
	norm=mean(det_eff_qty)
	det_eff_qty=det_eff_qty/norm
	logentry=logentry+string(10B)+$
		"Detector efficiencies computed using elastic peak sums."
endif
;
; Method 2 uses elastic peak sums based on individual fits of a Gaussian plus
; constant background to the rotated counts for each of the detectors.
if (method eq 2) then begin
	lines=["Detector"+strcompress(1)+" of"+strcompress(nhedet)]
	dcs_message_display,lines,base,group_leader=dcs_base,alignment=0
	for j=0,nhedet-1 do begin
		maxcounts=(*indataPtr)[round(nch/2),j]
		estimates=[maxcounts,tmean_rot[round(nch/2)],20.0,0.0,0.0]
		result=gaussfit(tmean_rot,(*indataPtr)[*,j],params,estimates=estimates,nterms=nterms)
		lotime=params[1]-mult*params[2]
		loch=where(tmean_rot-10.0 lt lotime and tmean_rot+10.0 ge lotime)
		loch=round(mean(loch))
		hitime=params[1]+mult*params[2]
		hich=where(tmean_rot-10.0 lt hitime and tmean_rot+10.0 ge hitime)
		hich=round(mean(hich))
		nchans=hich-loch+1
		intint=total((*indataPtr)[loch:hich,j])
		lo_int=total((*indataPtr)[loch-nchans+1:loch,j])
		hi_int=total((*indataPtr)[hich:hich+nchans-1,j])
		det_eff_qty[j]=intint-0.5*(lo_int+hi_int)
		if ((j+1) mod 25 eq 0) then begin
			widget_control,base,/destroy
			lines=["Processing detector"+strcompress(j+1)+" of"+strcompress(nhedet),$
				"Please be patient..."]
			dcs_message_display,lines,base,group_leader=dcs_base,alignment=0
		endif
	endfor
	widget_control,base,/destroy
	norm=mean(det_eff_qty)
	det_eff_qty=det_eff_qty/norm
	logentry=logentry+string(10B)+$
		"Detector efficiencies computed using individual elastic peak fits."
endif
Q=*(*(*(*davePtr).dataStrPtr).specificPtr).Q_propsPtr
dwi=exp(Q^2*0.066)
det_eff_qty=det_eff_qty*dwi
;dcs_plot,Q,det_eff_qty
;
ptr_free,indataPtr
;
; Clear any floating-point underflows that may have occurred in "gaussfit",
dummy=check_math(mask=32)
;
end