; $Id: dcs_docorrections.pro,v 1.17 2005/09/14 15:56:43 jrdcopley 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_qflagPtr,ngroups,elasticqvals,davePtr=davePtr,qflagPtr=qflagPtr
;************************************************************************************************
; This procedure is called if the grouping is in Q.
; First we construct the array Q(t,d), whose pointer is inelasticqPtr.
;	Then we construct a 2-d grouping array whose pointer is tdgPtr
;		('two dimensional group pointer').
;	Then we construct an array of pointers, qflagPtr, to 2-d arrays of 1s and 0s.
; All of these arrays are dimensioned ngtchan0*nhedet.
;
compile_opt strictarr
;
device,decomposed=0
loadct,13
;
progressBar = Obj_New($
	"SHOWPROGRESS",$
;	/cancelbutton,$;	(How does this work?)
	title="Work in progress. Please wait.",$
	xsize=300,$
	ysize=10,$
	message="Please be patient...",$
	color=184)
progressBar->Start
;
nhedet=(*(*(*davePtr).dataStrPtr).specificPtr).nhedet
dcs_getangles,angles,nhedet=nhedet
cos2th=cos((!pi/180.)*angles)
unitvec_d=intarr(nhedet)+1
;
lambda=(*(*(*davePtr).dataStrPtr).specificPtr).ch_wl
e0=float(!dcs_hsq2mn/lambda^2)
;
ngtchan0=(*(*(*davePtr).dataStrPtr).specificPtr).ngtchan0
efinal=e0-float((*(*(*(*davePtr).dataStrPtr).specificPtr).energy_propsptr)[0:ngtchan0-1])
unitvec_t=intarr(ngtchan0)+1
;
conversion=float(4.*!pi^2/!dcs_hsq2mn)
qsqptr=ptr_new(conversion*(e0+efinal#unitvec_d-2.*sqrt(e0*efinal)#cos2th),/no_copy)
sign_of_2th=(angles gt 0)*2-1
if (elasticqvals[0] lt 0.0) then $
	inelasticqPtr=ptr_new(sqrt(*qsqptr)*(unitvec_t#sign_of_2th),/no_copy) else $
	inelasticqPtr=ptr_new(sqrt(*qsqptr),/no_copy)
ptr_free,qsqptr
;
elasticqvals=float(elasticqvals)
tdgPtr=ptr_new(bytarr(ngtchan0,nhedet))
*tdgPtr=*inelasticqPtr ge elasticqvals[0] and *inelasticqPtr lt elasticqvals[1]
;
progressBar->Update,10
;
for k=1,ngroups-1 do begin
	*tdgPtr=*tdgPtr+(k+1)*$
		(*inelasticqPtr ge elasticqvals[k] and *inelasticqPtr lt elasticqvals[k+1])
	progressBar->Update,10+k/float(ngroups-1)*80
endfor
;
ptr_free,inelasticqPtr
;
qflagPtr=ptrarr(ngroups,/allocate_heap,/nozero)
;
progressBar->Update,100
;
for k=0,ngroups-1 do begin
	*qflagPtr[k]=(*tdgPtr eq k)
endfor
;
ptr_free,tdgPtr
;
progressBar->Destroy
Obj_Destroy, progressBar
;
device,decomposed=1
;
end


;************************************************************************************************
pro dcs_redefine_props,angle_props,Q_props,ngroups,davePtr=davePtr
;************************************************************************************************
; This procedure is called if the grouping is in Q.
; It redefines angle_props and Q_props.
; Note that (*Q_propsPtr)[*,0] will become the y component in histPtr.
;
compile_opt strictarr
;
;
lambda=(*(*(*davePtr).dataStrPtr).specificPtr).ch_wl
e0=float(!dcs_hsq2mn/lambda^2)
conversion=float(4.*!pi^2/!dcs_hsq2mn)
k0=sqrt(conversion*e0)
;
; The variable qrat is a 2-column array containing absolute (Q/k_0) limits for each group.
qrat=transpose(abs(Q_props[*,1:2])/k0)
;
; The variables min2th and max2th contain minimum and maximum values of abs(2theta).
nhedet=(*(*(*davePtr).dataStrPtr).specificPtr).nhedet
dcs_getangles,angles,nhedet=nhedet
min2th=min(abs(angles))
max2th=max(abs(angles))
cmin2th=cos(min2th*!dtor)
cmax2th=cos(max2th*!dtor)
argmin=cmin2th^2+qrat^2-1
argmax=cmax2th^2+qrat^2-1
;
; Expand Q_props to 7 rows.
; Rows 1-3 contain the mean, minimum and maximum Q for the group.
; Rows 4-5 contain the minimum and maximum energy transfer for the minimum Q in the group.
; Rows 6-7 contain the minimum and maximum energy transfer for the maximum Q in the group.
Q_props=[[Q_props],[fltarr(ngroups,4)]]
; Compute the quantities that go in rows 4-7 of Q_props.
jk=-1
for j=0,ngroups-1 do begin
	for k=0,1 do begin
		jk=jk+1
		elim=fltarr(2)
		nr=-1
		if (argmin[jk] gt 0.0) then begin
			roomin=sqrt(argmin[jk])
			if (qrat[jk] gt 1.0) then begin
				nr=nr+1
				elim[nr]=e0*(cmin2th+roomin)^2
			endif
			if (qrat[jk] le 1.0 and cmin2th gt 0.0) then begin
				nr=nr+1
				elim[nr]=e0*(cmin2th-roomin)^2
				nr=nr+1
				elim[nr]=e0*(cmin2th+roomin)^2
			endif
		endif
		if (argmax[jk] gt 0.0) then begin
			roomax=sqrt(argmax[jk])
			if (qrat[jk] gt 1.0) then begin
				nr=nr+1
				elim[nr]=e0*(cmax2th+roomax)^2
			endif
			if (qrat[jk] le 1.0 and cmax2th gt 0.0) then begin
				nr=nr+1
				elim[nr]=e0*(cmax2th-roomax)^2
				nr=nr+1
				elim[nr]=e0*(cmax2th+roomax)^2
			endif
		endif
		elim=e0-elim
		elim=elim[sort(elim)]
		Q_props[j,3+k*2]=elim[0]
		Q_props[j,3+k*2+1]=elim[1]
	endfor
endfor
end


;************************************************************************************************
pro dcs_make_eta,det_eff_qty,det_eff_err,grouping,ngroups,doQgrouping,$
	qflagPtr,ngtchan0,nhedet,fssc,fcsc,fcc,$
	sumeffssc_qty,sumeffssc_err,sumeffcsc_qty,sumeffcsc_err,sumeffcc_qty,sumeffcc_err
;************************************************************************************************
; This procedure constructs the quantity eta-hat(t,g) (now in 3 flavors)
; Quantities returned are dimensioned ngtchan0*ngroups.
;
compile_opt strictarr
;
unitvec_t=intarr(ngtchan0)+1
;
case ngroups of
;	One group.
	1: begin
		sumeffssc_qty = unitvec_t *      total( det_eff_qty*fssc)
		sumeffssc_err = unitvec_t * sqrt(total((det_eff_err*fssc)^2 ))
		sumeffcsc_qty = unitvec_t *      total( det_eff_qty*fcsc)
		sumeffcsc_err = unitvec_t * sqrt(total((det_eff_err*fcsc)^2 ))
		sumeffcc_qty  = unitvec_t *      total( det_eff_qty*fcc )
		sumeffcc_err  = unitvec_t * sqrt(total((det_eff_err*fcc )^2 ))
	end
;
; No grouping (hence nhedet groups).
	nhedet: begin
		sumeffssc_qty = unitvec_t # (det_eff_qty*fssc)
		sumeffssc_err = unitvec_t # (det_eff_err*fssc)
		sumeffcsc_qty = unitvec_t # (det_eff_qty*fcsc)
		sumeffcsc_err = unitvec_t # (det_eff_err*fcsc)
		sumeffcc_qty  = unitvec_t # (det_eff_qty*fcc )
		sumeffcc_err  = unitvec_t # (det_eff_err*fcc )
	end
;
;	Otherwise...
	else: begin
;
; Grouping in 2theta (rather than in Q):
		if (not doQgrouping) then begin
			sumeffssc_qty = fltarr(ngtchan0,ngroups)
			sumeffssc_err = fltarr(ngtchan0,ngroups)
			sumeffcsc_qty = fltarr(ngtchan0,ngroups)
			sumeffcsc_err = fltarr(ngtchan0,ngroups)
			sumeffcc_qty  = fltarr(ngtchan0,ngroups)
			sumeffcc_err  = fltarr(ngtchan0,ngroups)
			for k=0,ngroups-1 do begin
				list=where(grouping eq k+1)
				sumeffssc_qty[*,k] = unitvec_t #      total( det_eff_qty[list]*fssc[list])
				sumeffssc_err[*,k] = unitvec_t # sqrt(total((det_eff_err[list]*fssc[list])^2))
				sumeffcsc_qty[*,k] = unitvec_t #      total( det_eff_qty[list]*fcsc[list])
				sumeffcsc_err[*,k] = unitvec_t # sqrt(total((det_eff_err[list]*fcsc[list])^2))
				sumeffcc_qty[*,k]  = unitvec_t #      total( det_eff_qty[list]*fcc[list])
				sumeffcc_err[*,k]  = unitvec_t # sqrt(total((det_eff_err[list]*fcc[list])^2))
			endfor
		endif else begin

; Grouping in Q (rather than in angle):
			sumeffssc_qty = fltarr(ngtchan0,ngroups)
			sumeffssc_err = fltarr(ngtchan0,ngroups)
			sumeffcsc_qty = fltarr(ngtchan0,ngroups)
			sumeffcsc_err = fltarr(ngtchan0,ngroups)
			sumeffcc_qty  = fltarr(ngtchan0,ngroups)
			sumeffcc_err  = fltarr(ngtchan0,ngroups)
			for k=0,ngroups-1 do begin
				sumeffssc_qty[*,k] =      total((unitvec_t#(det_eff_qty*fssc)  ) * (*qflagPtr[k]),2)
				sumeffssc_err[*,k] = sqrt(total((unitvec_t#(det_eff_err*fssc)^2) * (*qflagPtr[k]),2) )
				sumeffcsc_qty[*,k] =      total((unitvec_t#(det_eff_qty*fcsc)  ) * (*qflagPtr[k]),2)
				sumeffcsc_err[*,k] = sqrt(total((unitvec_t#(det_eff_err*fcsc)^2) * (*qflagPtr[k]),2) )
				sumeffcc_qty[*,k]  =      total((unitvec_t#(det_eff_qty*fcc)   ) * (*qflagPtr[k]),2)
				sumeffcc_err[*,k]  = sqrt(total((unitvec_t#(det_eff_err*fcc)^2 ) * (*qflagPtr[k]),2) )
			endfor
		endelse
	end
endcase
;
end


;************************************************************************************************
pro dcs_make_ibhat,fbdata_qty,fbdata_err,fbnormalization,grouping,ngroups,doQgrouping,$
	speriod,qflagPtr,ngtchan0,nhedet,sumfb_qty,sumfb_err
;************************************************************************************************
; This procedure constructs the quantity I_B-hat(t,g)/tau.
; Quantities returned are dimensioned ngtchan0*ngroups if tnorm ne 0.
; If tnorm=0 (i.e. T_B^tot = 0), scalar zeros are returned.
;
compile_opt strictarr
;
tnorm=max(fbnormalization)
;
if (tnorm eq 0.0) then begin
	sumfb_qty=0.0
	sumfb_err=0.0
	return
endif
;
unitvec_t=intarr(ngtchan0)+1
;
case ngroups of
;	One group.
	1: begin
		sumfb_qty = unitvec_t*total(fbdata_qty)/speriod
		sumfb_err = unitvec_t*sqrt(total(fbdata_err^2))/speriod
	end
;
; No grouping (hence nhedet groups).
	nhedet: begin
		sumfb_qty = (unitvec_t#fbdata_qty)/speriod
		sumfb_err = (unitvec_t#fbdata_err)/speriod
	end
;
;	Otherwise...
	else: begin
;
; Grouping in 2theta (rather than in Q):
		if (not doQgrouping) then begin
			sumfb_qty=fltarr(ngtchan0,ngroups)
			sumfb_err=fltarr(ngtchan0,ngroups)
			for k=0,ngroups-1 do begin
				list=where(grouping eq k+1)
				sumfb_qty[*,k] = (unitvec_t#(total(fbdata_qty[list])))/speriod
				sumfb_err[*,k] = (unitvec_t#(sqrt(total(fbdata_err[list]^2))))/speriod
			endfor
		endif else begin

; Grouping in Q (rather than in angle):
			sumfb_qty=fltarr(ngtchan0,ngroups)
			sumfb_err=fltarr(ngtchan0,ngroups)
			for k=0,ngroups-1 do begin
				sumfb_qty[*,k] = total((unitvec_t#fbdata_qty)*(*qflagPtr[k]),2)/speriod
				sumfb_err[*,k] = sqrt(total((unitvec_t#fbdata_err^2)*(*qflagPtr[k]),2))/speriod
			endfor
		endelse
	end
endcase
;
end


;************************************************************************************************
pro dcs_docorrections,dcs_base,sumfiles,n_sumfiles,det_eff_qty,det_eff_err,grouping,ngroups,$
	angle_props,Q_props,grpfiles,temperatures,logentry=logentry,elasticqvals=elasticqvals,$
	fbdata_qty=fbdata_qty,fbdata_err=fbdata_err,fbnormalization=fbnormalization,doQgrouping,$
	message_str,skipdefaultchoices=skipdefaultchoices

;************************************************************************************************
; Performs fast background subtraction, detector grouping, correction for detector efficiencies,
; data normalization, and corrections using empty can and cadmium runs.
;
; Input variables:
; sumfiles and n_sumfiles contain information about the .sum files to be used in the procedure.
;	det_eff_qty and det_eff_err are 913-long vectors with obvious contents.
;	grouping and ngroups define, along with doQgrouping, how the data are to be grouped.
;	angle_props and Q_props contain angle and Q information about the groups.
; grpfiles are the names of the grouped data files.
; temperatures
;	elasticqvals
; fbdata_qty and fbdata_err contain background count rates and errors for each detector.
; fbnormalization is a set of 1s and 0s depending whether or not the detector is masked.

compile_opt strictarr
;
renorm_factor=1.0e6
;
;res=ptr_valid(count=nptr)
;print,'Entering docorrections, ',nptr,' pointers in use.'
;
; This is the number of files to be treated (S+C+A)
nfiles=fix(n_sumfiles[0]+n_sumfiles[1]+n_sumfiles[3])
;
; Generate the names of the files to be treated.
files=sumfiles[0:n_sumfiles[0]-1]
;
; Extract names of container and absorber files.
if (n_sumfiles[1] gt 0) then begin
	canfile=sumfiles[n_sumfiles[0]]
	files=[files,canfile]
endif else begin
	canfile=""
endelse
;
if (n_sumfiles[3] gt 0) then begin
	absfile=sumfiles[nfiles-1]
	files=[files,absfile]
endif else begin
	absfile=""
endelse

logentry=logentry+string(10B)+$
	"Background subtraction, detector grouping, detector efficiency correction, "+$
	"data normalization and corrections using empty can and cadmium runs have been performed."
logentry=logentry+string(10B)+$
	"Data renormalized. Multiplicative factor = "+strcompress(string(renorm_factor))+"."
;
lines=strarr(2*nfiles+4)
lines[0:2]=[$
	"Performing background subtraction, detector grouping,",$
	"beam monitor normalization, detector efficiency corrections,",$
	"and corrections using empty can and cadmium runs"]
bases=intarr(2*nfiles)
;
grpfiles=strarr(n_sumfiles[0])
temperatures=fltarr(n_sumfiles[0])
;
; Read the first file in order to get at quantities such as ch_wl, nhedet etc.
filename=sumfiles[0]
restore,filename=filename
;
nhedet=(*(*(*davePtr).dataStrPtr).specificPtr).nhedet
ngtchan0=(*(*(*davePtr).dataStrPtr).specificPtr).ngtchan0
;
; Calculate time between pulses at sample.
speriod=6.0e7/(*(*(*davePtr).dataStrPtr).specificPtr).ch_ms*$
	(*(*(*davePtr).dataStrPtr).specificPtr).ch_srdenom
;
if (doQgrouping) then begin
	dcs_redefine_props,angle_props,Q_props,ngroups,davePtr=davePtr
	dcs_define_qflagPtr,ngroups,elasticqvals,davePtr=davePtr,qflagPtr=qflagPtr
endif
;
; The normalization sum is the same for all unmasked detectors and 0 otherwise.
; These two lines reconstruct the detector masking vector M(d).
normalization_sum=(*(*(*davePtr).dataStrPtr).specificPtr).normalization_sum
masking=normalization_sum/max(normalization_sum)
;
; This code has been added to handle self-shielding corrections.
dcs_multiplechoice,dcs_base,ssf_decision,"Sample_selfshielding",$
	message_str,60,skipdefaultchoices=skipdefaultchoices
;
; ssf_decision=0 means angle-dependent self-shielding factors are to be included.
; ssf_decision=1 means there is a single angle-independent self-shielding factor.
; ssf_decision=2 means there is no self-shielding correction.
;
;	The following lines define nhedet-element arrays fssc, fcsc, and fcc.
fssc=fltarr(nhedet)+1.0
fcsc=fltarr(nhedet)+1.0
fcc=fltarr(nhedet)+1.0
;
if (ssf_decision eq 0) then begin
	filename=dialog_pickfile(path=!dcs_outp_directory,$
		title="Define the self-shielding file.")
	cancel=filename eq ""
	if (cancel) then begin
		ssf_decision=1
	endif else begin
		openr,unit,filename,/get_lun
		line=''
		readf,unit,line
		if (line eq 'Slab geometry') then begin
			nlines=0
			readf,unit,nlines
			for i=0,nlines-1 do readf,unit,line
			values=fltarr(7)
			for i=0,nhedet-1 do begin
				if (eof(unit)) then begin
					res=dialog_message("End-of-file encountered."+$
						"Angle-dependent self-shielding corrections will not be applied.")
					ssf_decision=1
					break
				endif
				readf,unit,values
				fssc[i]=values[1]
				fcsc[i]=values[2]
				fcc[i]=values[3]
			endfor
			free_lun,unit
		endif else begin
			res=dialog_message("Angle-dependent self-shielding "+$
				"corrections can only be performed for slab samples")
			ssf_decision=1
		endelse
	endelse
endif
;
if (ssf_decision eq 1) then begin
	title="Enter mean self-shielding factor"
	text="Mean self-shielding factor"
	value='1.0000'
	dave_get_values,group_leader=dcs_base,title,text,value,/checkall
	value=value[0]
	fcsc=value; (fssc and fcc are unity)
endif
;
; (if ssf_decision = 2, fssc, fcsc and fcc are all unity)
;
; Construct eta-hat(t,g) (now in 3 flavors).
dcs_make_eta,det_eff_qty,det_eff_err,grouping,ngroups,doQgrouping,$
	qflagPtr,ngtchan0,nhedet,fssc,fcsc,fcc,$
	sumeffssc_qty,sumeffssc_err,sumeffcsc_qty,sumeffcsc_err,sumeffcc_qty,sumeffcc_err
;
; Construct I_B-hat(t,g)/tau (scalar zeros if T_B^tot=0).
dcs_make_ibhat,fbdata_qty,fbdata_err,fbnormalization,grouping,ngroups,doQgrouping,$
	speriod,qflagPtr,ngtchan0,nhedet,sumfb_qty,sumfb_err
;
;;f_a=0.0
absfile_qty=0.0
absfile_err=0.0
;;fa_tnorm=0.0
;
;;f_c=0.0
canfile_qty=0.0
canfile_err=0.0
;;fc_tnorm=0.0
tcnorm=0.0
;
; Loop over files, decrementing the file index.
for i=nfiles-1,0,-1 do begin
	filename=files[i]
	icomp=nfiles-1-i
	lines[icomp*2+4]="Reading input file "+filename
	dcs_message_display,lines,base,group_leader=dcs_base,alignment=0
	bases[icomp*2]=base
;
	free_all_ptrs,davePtr
	restore,filename=filename
;
;	Put angle_props and Q_props into the DAVE structure.
	free_Ptr,(*(*(*davePtr).dataStrPtr).specificPtr).angle_propsPtr
	free_Ptr,(*(*(*davePtr).dataStrPtr).specificPtr).Q_propsPtr
	(*(*(*davePtr).dataStrPtr).specificPtr).angle_propsPtr=ptr_new(angle_props)
	(*(*(*davePtr).dataStrPtr).specificPtr).Q_propsPtr=ptr_new(Q_props)
;
; The quantities norm and tnorm are N_X^tot and T_X^tot/N_X^tot respectively.
;	We neglect their errors.
	norm=max((*(*(*davePtr).dataStrPtr).specificPtr).normalization_sum)
	tnorm=max((*(*(*davePtr).dataStrPtr).specificPtr).fbnormalization_sum)/norm
;
;	qtyptr points to I_X(t,d), errptr points to its error, xptr points to x.
; Dereference these quantities, free the pointer, redefine the pointer to
; reflect arrays dimensioned ngroups rather than nhedet, populate the structure.
;
	qtyPtr=ptr_new((*(*(*davePtr).dataStrPtr).commonStr.histPtr).qty)
	errPtr=ptr_new((*(*(*davePtr).dataStrPtr).commonStr.histPtr).err)
	xPtr=ptr_new((*(*(*davePtr).dataStrPtr).commonStr.histPtr).x)
;
	free_Ptr,(*(*davePtr).dataStrPtr).commonStr.histPtr
;
	(*(*davePtr).dataStrPtr).commonStr.histPtr=$
		ptr_new({qty:fltarr(ngtchan0,ngroups),$
		err:fltarr(ngtchan0,ngroups),$
		x:*xPtr,$
		y:(*(*(*(*davePtr).dataStrPtr).specificPtr).Q_propsPtr)[*,0]})
;
; In what follows, qty is I_X-hat(t,g)/N_X^tot and err is the
;	corresponding error.
;
	case ngroups of
;	One group.
		1: begin
			(*(*(*davePtr).dataStrPtr).commonStr.histPtr).qty = total((*qtyPtr),2)/norm
			(*(*(*davePtr).dataStrPtr).commonStr.histPtr).err = $
				sqrt(total((*errPtr)^2,2))/norm; this line fixed 5/8/02
		end
;
; No grouping (hence nhedet groups).
		nhedet: begin
			(*(*(*davePtr).dataStrPtr).commonStr.histPtr).qty = *qtyPtr/norm
			(*(*(*davePtr).dataStrPtr).commonStr.histPtr).err = *errPtr/norm
		end
;
;	Otherwise...
		else: begin
;
; Grouping in 2theta (rather than in Q):
			if (not doQgrouping) then begin
				for k=0,ngroups-1 do begin
					list=where(grouping eq k+1)
					(*(*(*davePtr).dataStrPtr).commonStr.histPtr).qty[*,k] = $
						total((*qtyPtr)[*,list],2)/norm
					(*(*(*davePtr).dataStrPtr).commonStr.histPtr).err[*,k] = $
						sqrt(total( (*errPtr)[*,list]^2,2 ))/norm
				endfor
			endif else begin
;
; Grouping in Q (rather than in angle):
				device,decomposed=0
				loadct,13
;
				progressBar = Obj_New($
					"SHOWPROGRESS",$
;					/cancelbutton,$;	(How does this work?)
					title="Work in progress. Please wait.",$
					xsize=300,$
					ysize=10,$
					message="Please be patient...",$
					color=184)
				progressBar->Start
;
				for k=0,ngroups-1 do begin
					(*(*(*davePtr).dataStrPtr).commonStr.histPtr).qty[*,k] = $
						total((*qtyPtr)*(*qflagPtr[k]),2)/norm
					(*(*(*davePtr).dataStrPtr).commonStr.histPtr).err[*,k] = $
						sqrt(total( (*errPtr)^2*(*qflagPtr[k]),2 ))/norm
					progressBar->Update,10+(k+1)/float(ngroups)*90
				endfor
;
				progressBar->Destroy
				Obj_Destroy, progressBar
;
device,decomposed=1
			endelse
		end
	endcase
;
	if (filename eq canfile) then begin
;;		title="Enter mean sample attenuation factor"
;;		text="Mean sample attenuation factor"
;;		value='1.0000'
;;		dave_get_values,group_leader=dcs_base,title,text,value,/checkall
;;		f_c=float(value[0])
;;		canfile_qty=f_c*(*(*(*davePtr).dataStrPtr).commonStr.histPtr).qty
;;		canfile_err=f_c*(*(*(*davePtr).dataStrPtr).commonStr.histPtr).err
;;		fc_tnorm=f_c*tnorm
		canfile_qty=(*(*(*davePtr).dataStrPtr).commonStr.histPtr).qty
		canfile_err=(*(*(*davePtr).dataStrPtr).commonStr.histPtr).err
;;		fc_tnorm=f_c*tnorm
		tcnorm=tnorm
		lines[icomp*2+5]="Processing container file "+filename
		dcs_message_display,lines,base,group_leader=dcs_base,alignment=0
		bases[icomp*2+1]=base
	endif
;
	if (filename eq absfile) then begin
;;		f_a=1.0-f_c
;;		absfile_qty=f_a*(*(*(*davePtr).dataStrPtr).commonStr.histPtr).qty
;;		absfile_err=f_a*(*(*(*davePtr).dataStrPtr).commonStr.histPtr).err
		absfile_qty=(*(*(*davePtr).dataStrPtr).commonStr.histPtr).qty
		absfile_err=(*(*(*davePtr).dataStrPtr).commonStr.histPtr).err
;;		fa_tnorm=f_a*tnorm
		lines[icomp*2+5]="Processing absorber file "+filename
		dcs_message_display,lines,base,group_leader=dcs_base,alignment=0
		bases[icomp*2+1]=base
	endif
;
; The absorber file is no longer included in the correction procedure.
	if (filename ne absfile and filename ne canfile) then begin
		qty=(*(*(*davePtr).dataStrPtr).commonStr.histPtr).qty
		err=(*(*(*davePtr).dataStrPtr).commonStr.histPtr).err
;
;		scorr_qty=(   (qty - f_a*absfile_qty - f_c*canfile_qty) - sumfb_qty*(tnorm - fa_tnorm - fc_tnorm) )  /sumeff_qty
;		scorr_err=sqrt(  (err^2 + (f_a*absfile_err)^2 + (f_c*canfile_err)^2) +	sumfb_err^2*(tnorm-fa_tnorm-fc_tnorm)^2  ) /sumeff_qty
;
		ssfratio=sumeffcsc_qty/sumeffcc_qty
		scorr_qty=  	  ( qty   -  canfile_qty*ssfratio    -  sumfb_qty   * (tnorm - tcnorm*ssfratio)   ) /sumeffssc_qty
		scorr_err=	sqrt( err^2 + (canfile_err*ssfratio)^2 +	sumfb_err^2 * (tnorm - tcnorm*ssfratio)^2 ) /sumeffssc_qty
;
; Test for divide by 0.0; replace +/-NaN with 0.0
		indices=where(finite(scorr_qty) eq 0,count)
		if (count gt 0) then begin
			scorr_qty[indices]=0.0
			scorr_err[indices]=0.0
		endif
		(*(*(*davePtr).dataStrPtr).commonStr.histPtr).qty=scorr_qty*renorm_factor
		(*(*(*davePtr).dataStrPtr).commonStr.histPtr).err=scorr_err*renorm_factor
;
		strput,filename,"cor.dave",strlen(filename)-8
		grpfiles[i]=filename
		lines[icomp*2+5]="Writing sample output file "+filename
		dcs_message_display,lines,base,group_leader=dcs_base,alignment=0
		bases[icomp*2+1]=base
;
;	widget_control,/hourglass
;
		(*(*davePtr).dataStrPtr).commonStr.ylabel = 'mean Q'
		(*(*davePtr).dataStrPtr).commonStr.yunits = 'wavevector:A-1'
		(*(*davePtr).dataStrPtr).commonStr.ytype = 'points'
		(*(*davePtr).dataStrPtr).commonStr.histlabel = 'Grouped Intensity'
;
		(*(*(*davePtr).dataStrPtr).commonStr.treatmentPtr)=$
			strsplit(logentry,string(10B),/extract)
;
		save,davePtr,filename=filename
;
		temperatures[i]=(*(*(*davePtr).dataStrPtr).specificPtr).temp_setpoint > 0
;
	endif
;
	ptr_free,qtyPtr
	ptr_free,errPtr
	ptr_free,xPtr
endfor
;
; NB qflagPtr (if it exists) is a pointer array.
if (total(ptr_valid(qflagPtr)) gt 0) then ptr_free,qflagPtr
;
for icomp=2*nfiles-1,0,-1 do widget_control,bases[icomp],/destroy
;
free_all_ptrs,davePtr
;
;res=ptr_valid(count=nptr)
;print,'Exiting docorrections, ',nptr,' pointers in use.'
;
end
