; $Id: dcs_normalization_counts.pro,v 1.13 2005/09/14 15:54:42 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_normalization_counts,dcs_base,files_selected,nfiles_selected,nruns,$
	normalization_values,fbnormalization_values,davePtr=davePtr,bitflips,maxbitflips,$
	skipmessage,normalizetoBM1,dontsavebminfo,nhedet=nhedet,fbdecision=fbdecision,$
	message_str,skipdefaultchoices=skipdefaultchoices,logentry=logentry
;************************************************************************************************
; Runs are read into "davePtr" and the beam monitor peaks are analyzed.
; Input parameters include files_selected, nruns
; Output parameters include normalization_values and fbnormalization_values
;
compile_opt strictarr
;
omithisto=1
printinfo=0
filetype="binary"
;
;	Initialize a structure and replicate it for each beam monitor, each run.
; -- allcounts is the sum of counts in all time channels.
; -- peakcounts, positions, peakwidths and backgrounds are from fit of Gaussian
;			plus constant background.
; -- postheory is theoretical position of peak.
; -- posmaxmum is experimental position of peak.
;
bmvalues={allcounts:0.0,peakcounts:0.0,positions:0.0,peakwidths:0.0,$
	backgrounds:0.0,postheory:0.0,posmaxmum:0.0}
bmvalues=replicate(bmvalues,nruns,4)
normalization_values=fltarr(nruns)
fbnormalization_values=fltarr(nruns)
other_values=fltarr(nruns,4)
duration=fltarr(nruns)
nterms=4
params=fltarr(nterms)
;
; Loop over all runs.
for irun=0,nruns-1 do begin
	file=files_selected[irun]
;
	dcs_readdatafile,file,omithisto,printinfo,filetype,davePtr=davePtr
;
; Store run duration in the array "duration".
	duration[irun]=(*(*(*davePtr).dataStrPtr).specificPtr).duration
;
; Code that is executed for the first run only STARTS HERE.
	if (irun eq 0) then begin
;
; 		Define quantities needed to obtain beam monitor peak positions and intensities.
;			Master speed, master period, period at sample.
				mspeed=(*(*(*davePtr).dataStrPtr).specificPtr).ch_ms
				mperiod=60000000.0/mspeed
				speriod=mperiod*(*(*(*davePtr).dataStrPtr).specificPtr).ch_srdenom
;
;			Determine "ibm" and "other":
;			"ibm" contains the locations of the beam monitor channels in "histohigh"
;			"other" contains the locations of the EXT, CH5, CH7 and AND channels.
				lines=strarr(5)
				bm_dis=[0.0,(*(*(*davePtr).dataStrPtr).specificPtr).fc_dis]
				lines[0]="Beam monitor DAQ board numbers"
				lines[1]="Assumed distances from sample (mm)"
				ibm=intarr(4)
				other=intarr(4)
;
; 		Following loop is over the 17 Lemo inputs in crate 1 (the lower crate).
				for k=0,16 do begin
					string123=strmid((*(*(*davePtr).dataStrPtr).specificPtr).highsource[k],0,3)
					if (string123 eq "FC0") then begin
						ibm[0]=k
						lines[0]=lines[0]+": BM0,"+strcompress(k)
					endif
					if (string123 eq "FC1") then begin
						ibm[1]=k
						lines[0]=lines[0]+": BM1,"+strcompress(k)
						lines[1]=lines[1]+": BM1,"+strcompress(bm_dis[1])
					endif
					if (string123 eq "FC2") then begin
						ibm[2]=k
						lines[0]=lines[0]+": BM2,"+strcompress(k)
						lines[1]=lines[1]+": BM2,"+strcompress(bm_dis[2])
					endif
					if (string123 eq "FC3") then begin
						ibm[3]=k
						lines[0]=lines[0]+": BM3,"+strcompress(k)
						lines[1]=lines[1]+": BM3,"+strcompress(bm_dis[3])
					endif
					if (string123 eq "EXT") then other[0]=k
					if (string123 eq "CH5") then other[1]=k
					if (string123 eq "CH7") then other[2]=k
					if (string123 eq "AND") then other[3]=k
				endfor
;
;			(The variables "ibm" and "other" have been defined.)
;
;			Incident wavelength, incident inverse velocity, resolution mode
				lambda=(*(*(*davePtr).dataStrPtr).specificPtr).ch_wl; incident wavelength
				tau0=lambda/!dcs_hom; incident inverse velocity
				resmode=(*(*(*davePtr).dataStrPtr).specificPtr).ch_res; resolution mode
;
;			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)
;
; 		Define time channel widths for high end data (in us).
				tchanwid=(*(*(*davePtr).dataStrPtr).specificPtr).tchanlook[1,*]*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.
				ngtchan1=(*(*(*davePtr).dataStrPtr).specificPtr).ngtchan1
;
;			Truncate tstart and tstop.
				tstart=tstart[0:ngtchan1-1]
				tstop=tstop[0:ngtchan1-1]
				lines[2]=" Number of time channels = "+strcompress(ngtchan1)
				lines[3]=" Time between pulses at sample = "+string(speriod)+" us."
				lines[4]=" Time channel"+strcompress(ngtchan1)+" ends at "+string(tstop[ngtchan1-1])+" us."
;
; 		Add information to log.
				for k=0,4 do logentry=logentry+string(10B)+lines[k]
;
; 		Centers of time channels (for use in gaussfit)
				tmean=0.5*(tstart+tstop)
;
; 		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: begin
						if (!debug) then print,"Invalid STARTCHOICE"
						nruns=0
						return
					end
				endcase
;
;			If desired, print information message.
				if (not skipmessage) then begin
					blabel=['OK']
					def_button=0
					dcs_multichoice_modal,lines,decision,blabel,def_button,$
						group_leader=dcs_base,alignment=0
				endif
;
; 		If not normalizing to BM1, ask how to normalize among data sets.
				if (not normalizetoBM1) then begin
					dcs_multiplechoice,dcs_base,decision,"Normalization",message_str,80,$
						skipdefaultchoices=skipdefaultchoices,nbuttonrow=[3,4,2],blabel=blabel
				endif else begin
;					if (not lmgr(/runtime)) then begin
						mess=["This new dialog may surprise you.",$
						"It is included to satisfy the special needs of one of our users.",$
						"","Unless you are that person the likely response is ''No'' (or hit Enter).",$
						"","We will endeavor to handle this particular user's need in a more",$
						"transparent fashion some time in the near future. Meanwhile",$
						"our apologies for any confusion that this may be causing.",$
						"","Do you wish to normalize to duration (rather than to the beam monitor)?"]
						res=dialog_message(/question,mess,/default_no)
						if (res eq "No") then decision=2 else decision=7
;					endif else begin
;						decision=2
;					endelse
					blabel=['BM0','BM1 total','BM1 peak sum',$
						'EXT','CH5','CH7','AND','Duration','Main menu']
				endelse
;
				if (decision eq 8) then begin
					nruns=0
					return
				endif
;
; 		If there is background (dark count) data, and if the normalization of sample
; 		runs is other than to "EXT" or "Duration", ask how to normalize the background data.
;
; BEGIN
;	The following lines have been commented out (SIMPLIFICATION PROJECT)
;				case decision of
;					3: fbdecision=0
;					6: fbdecision=1
;					else: begin
;						if (nfiles_selected[4] eq 0) then begin
;							fbdecision=-1
;						endif else begin
;							dcs_multiplechoice,dcs_base,fbdecision,"FBNormalization",message_str,80,$
;								skipdefaultchoices=skipdefaultchoices,blabel=fbblabel
;						endelse
;					end
;				endcase
;	... and the following line has been inserted (SIMPLIFICATION PROJECT)
				fbdecision=0
				fbblabel=['EXT','Duration','Main menu']
; END
;
				if (fbdecision eq 2) then begin
					nruns=0
					return
				endif
;
			endif; Code that is executed for the first run only STOPS HERE.
;
;
	lines=["Treatment of beam monitors...",$
		"Reading file "+strcompress(irun+1)+":"+file]
	dcs_message_display,lines,base,group_leader=dcs_base,alignment=0

; Treat each beam monitor in turn.
	for j=0,3 do begin
;
; 	Compute integrated counts over all time channels.
		allcounts=$
			ulong(total((*(*(*davePtr).dataStrPtr).specificPtr).histohigh[0:ngtchan1-1,ibm[j]],/double))
		bmvalues[irun,j].allcounts=allcounts
;
; 	Examine beam monitor peaks, excluding BM0 and monitors with very few total counts.
		if (j gt 0 and allcounts gt 100) then begin
			counts=(*(*(*davePtr).dataStrPtr).specificPtr).histohigh[0:ngtchan1-1,ibm[j]]
;
; 		If necessary perform bitflip correction	to data.
			nbitflips=0
			for k=0,maxbitflips-1 do begin
				if (bitflips[1,k,irun] eq nhedet+ibm[j]) then begin
					nbitflips=nbitflips+1
					if (nbitflips gt 2) then begin
						result=dialog_message("Too many bitflips in dcs_normalization_counts.",/error)
						stop
					endif
					bitflipchannel=bitflips[0,k,irun-1]
					bitflipamount=bitflips[2,k,irun-1]
					result=dialog_message(["A bitflip has been corrected in dcs_normalization_counts.",$
					"Beam monitor"+strcompress(j)+", time channel"+strcompress(bitflips[0,k,irun-1])+$
					", correction"+strcompress(bitflips[2,k,irun-1])])
					allcounts=allcounts-bitflipamount
					counts[bitflipchannel]=counts[bitflipchannel]-bitflipamount
				endif
			endfor
;
;			Compute time of flight from sample to beam monitor.
			t_sbm=tau0*bm_dis[j]
;
; 		Hence time of arrival at beam monitor, relative to DAQ start.
;			tbm_the=dave_modulud(t_sbm-ts0,speriod) (OLD VERSION)
			tbm_the=(t_sbm-ts0) mod speriod + speriod * ((t_sbm-ts0) lt 0.0)
;
; 		Compute corresponding time channel.
			tch_the=where(tstart le tbm_the and tbm_the lt tstop)
;
;			Determine channel(s) with maximum number of counts, corresponding times,
;				the mean of these times, and the corresponding time channel.
			maxcounts=max(counts)
			tch_exp=where(counts eq maxcounts)
			tbm_exp=(tstart[tch_exp]+tstop[tch_exp])/2.0
			tbm_exp=mean(tbm_exp)
			tch_exp=where(tstart le tbm_exp and tbm_exp lt tstop)
;
; 		Renumber channels so that peak count is roughly in channel ngtchan1/2.
; 		The number of channels shifted is tch_shift.
;			tch_shift=dave_modulus(round(-tch_exp[0]+ngtchan1/2),ngtchan1) (OLD VERSION)
			tch_shift=round(-tch_exp[0]+ngtchan1/2) mod ngtchan1 + ngtchan1 * (round(-tch_exp[0]+ngtchan1/2) lt 0.0)
			counts_rot=shift(counts,tch_shift)
			tmean_rot=shift(tmean,tch_shift)
			addition=fltarr(ngtchan1)+speriod
			if (tch_shift ne 0) then addition[0:tch_shift-1]=0
			tmean_rot=tmean_rot+addition
;
; 		Check whether there have been any floating-point underflows.
			if (check_math(mask=32)) then if (!debug) then print, "NB: floating-point underflows"
;
; 		Compute time shift corresponding to tch_shift.
;			tshift=tstop[dave_modulus(ngtchan1-1-tch_shift,ngtchan1)] (OLD VERSION)
			tshift=tstop[(ngtchan1-1-tch_shift) mod ngtchan1 + ngtchan1 * ((ngtchan1-1-tch_shift) lt 0.0)]
;
; 		Create estimates of Gaussian fit parameters.
			estimates=[maxcounts,tshift+0.5*speriod,20.0,0.0]
;
; 		Fit a Gaussian plus constant background to the rotated BM data.
			result=gaussfit(tmean_rot,counts_rot,params,estimates=estimates,nterms=nterms)
;
; 		Ensure that peak position lies within period at sample.
;			params[1]=dave_modulus(params[1],speriod) (OLD VERSION)
			params[1]=(params[1] mod speriod) + speriod * (params[1] lt 0.0)
;
; 		Populate the bmvalues array.
			bmvalues[irun,j].peakcounts =params[0]
			bmvalues[irun,j].positions  =params[1]
			bmvalues[irun,j].peakwidths =params[2]
			bmvalues[irun,j].backgrounds=params[3]
			bmvalues[irun,j].postheory  =tbm_the
			bmvalues[irun,j].posmaxmum  =tbm_exp
;
; 		Clear any floating-point underflows that occurred in "gaussfit",
; 		or in calculating the theoretical counts.
			dummy=check_math(mask=32)
;
; 		If detailed information is to be saved, compute the fitted counts and save to file with
; 		fit parameters, etc.
			if (not dontsavebminfo) then begin
				counts_calc=params[0]*exp(-0.5*((tmean-params[1])/params[2])^2)+params[3]

;				Generate file name.
				filename="_"+strcompress(irun,/remove_all)+"_0"+strcompress(j,/remove_all)+".bmd"
				filename=!dcs_outp_directory+!dcs_rootname+filename
				save,filename=filename,allcounts,params,tbm_the,tbm_exp,counts,counts_calc
			endif
;
		endif; end of condition that it is not BM0 and there are at least a few counts.
	endfor; end of loop over beam monitors
;
;	The array "other_values" contains counts for "EXT", "CH5", CH7" and "AND".
	other_values[irun,0]=(*(*(*davePtr).dataStrPtr).specificPtr).detsum[913+other[0]]
;
	if (startchoice eq "CH5") then $
		other_values[irun,1]=(*(*(*davePtr).dataStrPtr).specificPtr).nframes else $
		other_values[irun,1]=(*(*(*davePtr).dataStrPtr).specificPtr).detsum[913+other[1]]
;
	if (startchoice eq "CH7") then $
		other_values[irun,2]=(*(*(*davePtr).dataStrPtr).specificPtr).nframes else $
		other_values[irun,2]=(*(*(*davePtr).dataStrPtr).specificPtr).detsum[913+other[2]]
;
	if (startchoice eq "AND") then $
		other_values[irun,3]=(*(*(*davePtr).dataStrPtr).specificPtr).nframes else $
		other_values[irun,3]=(*(*(*davePtr).dataStrPtr).specificPtr).detsum[913+other[3]]
;
	widget_control,base,/destroy
endfor; end loop over runs
;
mincts=25; Peak intensity criterion for deciding whether or not there is a peak.
;
; Add information to log.
logentry=logentry+string(10B)+"The beam monitor peaks have been analyzed."
if (not dontsavebminfo) then logentry=logentry+string(10B)+$
		"Beam monitor information has been saved in a series of Run_BM_ files."
;
dummy=where(bmvalues.(1) gt mincts,count)
if (count eq 0) then $
	logentry=logentry+string(10B)+"No peaks have peak intensity >"+strcompress(mincts)+"."
;
if (count gt 0) then begin
	logentry=logentry+string(10B)+"The following lines contain the results of Gaussian fits: "+$
		string(10B)+"Results are only printed for peaks with peak intensity >"+strcompress(mincts)+"."+$
		string(10B)+"Filename, BM_ID, peak intensity, position (us), "+$
		"width (us), background, and ''integrated intensity''."
	;
	for ibm=0,3 do begin
		for irun=0,nruns-1 do begin
			if (bmvalues[irun,ibm].(1) gt 10) then begin
				pvals=string(bmvalues[irun,ibm].(1),bmvalues[irun,ibm].(2),$
					bmvalues[irun,ibm].(3),bmvalues[irun,ibm].(4),$
					bmvalues[irun,ibm].(1)*bmvalues[irun,ibm].(3)*sqrt(2.*!pi),$
					format='(f12.0,f9.1,f7.1,f7.1,f15.1)')
					logentry=logentry+string(10B)+$
						"File "+files_selected[irun]+": BM"+strcompress(ibm,/remove_all)+pvals
			endif
		endfor
	endfor
endif
;
case decision of ; (decision 2 is BM1 "integrated intensities")
	0: normalization_values=bmvalues[*,0].allcounts; BM0 total counts
	1: normalization_values=bmvalues[*,1].allcounts; BM1 total counts
	2: normalization_values=bmvalues[*,1].peakcounts*bmvalues[*,1].peakwidths
	3: normalization_values=other_values[*,0]; EXT counts
	4: normalization_values=other_values[*,1]; CH5 counts
	5: normalization_values=other_values[*,2]; CH7 counts
	6: normalization_values=other_values[*,3]; AND counts
	7: normalization_values=duration; run duration
	else:
endcase
;
case fbdecision of
	0: fbnormalization_values=other_values[*,0]; EXT counts
	1: fbnormalization_values=duration; run duration
	else:
endcase
;
logentry=logentry+string(10B)+$
	"Normalization other than to background will be performed using "+$
	blabel[decision]+" counts."
;
if (nfiles_selected[4] ne 0) then begin
	logentry=logentry+string(10B)+"Normalization to background will be performed using "+$
		fbblabel[fbdecision]+" counts."
		ratio_ext=fbnormalization_values/duration
	if (!debug) then print,ratio_ext
	extpgmin=!ext_pgen*0.98
	extpgmax=!ext_pgen*1.02
	for irun=0,nruns-1 do begin
		if (ratio_ext[irun] lt extpgmin or ratio_ext[irun] gt extpgmax) then begin
			res=dialog_message("It appears that the external"+strcompress(!ext_pgen)+$
				" Hz pulser was not in use.  The ''fast background'' run will be ignored.")
			nfiles_selected[4]=0
			files_selected=files_selected[0:nruns-2]
			return
		endif
	endfor
endif
;
end