; $Id: dcs_bitflip_correct.pro,v 1.3 2005/02/14 17:50:32 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.
;
;###############################################################################
function dcs_bitflip_correct,histodata,histohigh,detsum,timsum,bitflips=bitflips,maxbf=maxbf,$
	inputfilename=inputfilename,onlydetect=onlydetect,printlevel=printlevel
;************************************************************************************************
;
compile_opt strictarr
;
if (n_elements(inputfilename) eq 0) then begin
	filename=""
endif else begin
	filename="File"+inputfilename
endelse
;
if (n_elements(onlydetect) eq 0) then onlydetect=0
;
if (n_elements(printlevel) eq 0) then printlevel=0
;
; Maximum number of bitflips.
dtcountmax=maxbf
;
;	Sum He3 data and high end data over detectors and time channels
datatimsum=long(total(histodata,2,/double))+long(total(histohigh,2,/double))
datadetsum=[long(total(histodata,1,/double)),long(total(histohigh,1,/double))]
;
; Form differences (which would all be 0 if there were no bitflips).
tdh1=datadetsum-long(detsum)
tdh2=datatimsum-long(timsum)
;
;	Determine detector channels and time channels corresponding to bitflips.
bad_t=where(tdh2 ne 0,tcount); disagreement in time sums over detector
bad_d=where(tdh1 ne 0,dcount); disagreement in detector sums over time
;
;	If the number of bitflip detector channels is not the same as
;	the number of bitflip time channels we return, since the bitflips
; cannot be determined unambiguously.
if (dcount ne tcount) then begin
	message=[filename,$
		"Number of bitflip detector channels differs",$
		"from number of bitflip time channels"]
	if (printlevel ge 1) then x=dialog_message(message)
	return,-1
endif
;
dtcount=dcount
;
; If there are no bitflips we return.
if (dtcount eq 0) then begin
	message=[filename,strcompress("No bitflips detected.")]
	if (printlevel ge 3) then x=dialog_message(message)
	return,0
endif
;
;	If there are too many bitflips we return.
if (dtcount gt dtcountmax) then begin
	message=[filename,strcompress("More than"+string(dtcountmax)+" bitflips.")]
	if (printlevel ge 1) then x=dialog_message(message)
	return,-dtcount
endif
;
;	Determine numbers of unique differences and ensure that they are
; both identical to the number of bitflips, since if two or more
; differences are identical the individual bitflips cannot be
;	determined unambiguously.
nd=n_elements(uniq(tdh1,sort(tdh1)))-1
nt=n_elements(uniq(tdh2,sort(tdh2)))-1
if (nd ne dtcount or nd ne dtcount) then begin
	message=[filename,strcompress("Cannot unambiguously correct for bitflips.")]
	if (printlevel ge 1) then x=dialog_message(message)
	return,-dtcount
endif
;
;	Start correcting for bitflips.
message=filename
if (dtcount gt 1) then begin
	message=[message,strcompress(string(dtcount),/remove_all)+" bitflips found.",""]
endif else message=[message,'1 bitflip found.',""]
;
bitflips=lonarr(3,maxbf)
for k=0,dtcount-1 do begin
	badt=bad_t[k]
	difference=tdh2[badt]
	power=alog(difference)/alog(2)
	if (abs(2.0^power-difference) lt 0.1) then power=round(power)
	badd=where(tdh1 eq difference)
	message=[message,strcompress("Bitflip detected at [t,d] = ["+$
		string(badt)+","+string(badd)+"].")]
	message=[message,$
		strcompress("Magnitude ="+string(difference)+" = 2^"+strcompress(power,/remove_all)+".")]
	if (not onlydetect) then begin
		if (badd le 912) then begin
			histodata[badt,badd]=histodata[badt,badd]-difference
		endif else begin
			histohigh[badt,badd-913]=histohigh[badt,badd-913]-difference
		endelse
		message=[message,$
			strcompress("Bitflip has been corrected."),""]
	endif else begin
		message=[message,$
			strcompress("Bitflip has **NOT** been corrected."),""]
	endelse
	bitflips[0,k]=badt
	bitflips[1,k]=badd
	bitflips[2,k]=difference
endfor
if (printlevel ge 2) then x=dialog_message(message)
;
return,dtcount
end
