include	"vpthit.h"

#  vpt_remove -- Remove particle events for high speed photometer files
#
#  Description:
#  ------------
#  
#  Date		Author			Description
#  ----		------			-----------
#  29-Oct-1990  J.-C. Hsu		rewrite in SPP
#------------------------------------------------------------------------------

procedure vpt_remove (indata, inmask, scheme, npts, type, 
			samptm, sigmas, dthresh, athresh, lowflag, hitval, 
			mean, err, remove)

						## input
real	indata[npts], inmask[npts]
char	scheme[SZ_SCHEME]
int	npts
int	type
real	samptm
real	sigmas, athresh, dthresh, hitval
bool	lowflag
						## output
real	mean, err
int	remove
						## local
double	sum, sumsq, ave, stndev
real	thresh, low_thresh
int	ngood, i, j, k, xremove

int	strsearch()
#==============================================================================
begin
	# set the threshold level
	if (strsearch (scheme, "stat") > 0)
	    thresh = INDEF
	else {
	    if (type == DIGITAL)
		thresh = dthresh
	    else
		thresh = athresh
	}

	# reset counters
	ngood = 0
	remove = 0
	xremove = 0

	# reset sums
	sum = 0.d0
	sumsq = 0.d0

	# accumulate sum and sum square of the input data and remove particle 
	# events
	do i = 1, npts {
	    if (nint(inmask[i]) == OKVAL) {
		if (indata[i] > thresh) {
		    remove = remove + 1
		    inmask[i] = hitval
		} else {
		    ngood = ngood + 1
		    sum = double (indata[i]) + sum
		    sumsq = double (indata[i])**2 + sumsq
		    inmask[i] = OKVAL
		}
	    }
	}

	# iterate the above process for statistical method until no more 
	# points are removed
	if (strsearch (scheme, "stat") > 0) {
	    do j = 1, npts {
		if (ngood == 0)
		    call error (1, "all points are removed")
		else {
		    ave = sum / double(ngood)
		    stndev = sqrt(abs(sumsq / double(ngood) - ave**2))
		}

		# if the standard deviation is small (<1) set it to 1
		if (stndev < 1.d0) 
		    stndev = 1.d0
		thresh = ave + sigmas * stndev
		if (lowflag)
		    low_thresh = ave - sigmas * stndev
		else
		    low_thresh = -INDEF
		xremove = remove
		sum = 0.d0
		sumsq = 0.d0
		ngood = 0

		# accumulate sum and sum square of the input data and remove 
		# particle events
		do k = 1, npts {
	    	    if (nint(inmask[k]) == OKVAL) {
		        if (indata[k] > thresh || indata[k] < low_thresh) {
			    remove = remove + 1
			    inmask[k] = hitval
		    	} else {
		    	    ngood = ngood + 1
		    	    sum = double (indata[k]) + sum
		    	    sumsq = double (indata[k])**2 + sumsq
			    inmask[k] = OKVAL
		        }
	    	    }
		}
		if (remove == xremove)
		    break
	    }
	}
	# calculate the final mean rate/DN and its mean error
	if (ngood <= 0) 
	    call error (1, "all points are removed, threshold set too low")

	# if only one point left, use square root of the count as error for
	# digital data and half DN for analog data
	ave = sum / double(ngood)
	mean = ave
	if (ngood == 1) {
	    if (type == DIGITAL)
		err = sqrt(ave)
	    else
		err = 0.5
	} else
	    err = sqrt(abs((sumsq/double(ngood) - ave**2)) / double(ngood-1))
	
	# convert counts to count rate for digital data
	if (type == DIGITAL) {
	    mean = ave / samptm
	    err = err / samptm
	}
end