
### MBT added 12-2-08

userError<-function (msg) 
{
    stop(msg, call. = FALSE)
}
### MBT end

findBsEnv <- function(specDat, idsWindow){

t1 <- 1:length(specDat);

trs<-seq(from=1, to=length(t1), by=idsWindow);
k<-1; sigr<-0*trs;

ks <- (idsWindow-1)/2; ## assuming ODD idsWindow
sigr[1]<-sum(specDat[1:ks]); sigr[length(sigr)]<-sum(specDat[(trs[length(trs)]-ks):trs[length(trs)]]);

for (jk in seq(from=trs[2], to=trs[length(trs)-1], by=idsWindow)){
	k <- k+1; sigr[k] <- sum(specDat[(jk-ks):(jk+ks)]);
}

x<- diff(sigr);
xl<-x[1:(length(x)-1)]; xr<-x[2:length(x)];
nearbot <- (xl<0 & xr>=0) | (xl==0 & xr >0);
xx <- 1:length(x);
nb <- c(1, xx[nearbot]+1); ## tentative minimum
smax <- max(sigr[nb]); imax <- which(sigr[nb]==smax);
ts <- nb[imax];
#print(trs[ts]) ## debug
nxt <- imax+1;

while ((nxt <=length(nb)) & (sigr[nb[nxt]]>0.01*sigr[nb[imax]]) & (sigr[nb[nxt]]>1)) {
	if (sigr[nb[nxt]]<sigr[ts[length(ts)]]) {
		ts<-c(ts, nb[nxt]);
	}
	nxt<-nxt+1;
}
tSampl <- trs[ts]; sigSampl <- sigr[ts]/idsWindow;

envBS <- matrix(tSampl, nrow=length(tSampl), ncol=2);
envBS[,1]<-tSampl; envBS[,2]<-sigSampl;

return(as.matrix(envBS));
}

ExpBkgr_call <- function(tofList, bsSamples, offsetSamples, modelChoice, dataStart)
{

## Usage: Correct Background
## 
## Input Parameters:
## tofList – list of TOF spectra with meta-data
## bsSamples – TOF indices of data to use for model baseline fit 
##             (index is in original, non-truncated record), e.g.: 
##             c([t1:t2],[t3:t4],[t5:t6])
## When automated fit is preferred, provide two values in this field:
## c(Window, corPercent), where 
## idsWindow –- down-sampling step for finding baseline envelope, should be ODD
##           default for Bruker @ 500 MHz sampling rate is 199
## corPercent -– percent correction of baseline envelope for STD of noise
##               default for PCA’07 settings is 0.05 
## offsetSamples –- TOF interval for the estimate of constant offset, e.g.:
##                 c(120000:122000) 
## modelChoice --  choice of analytical baseline model:
##                 1 – linear, 2 – exponential 3 -- Gaussian
##

## use with new meta-data
# dataStart<-as.integer(tofListMetaData[1, "timeOfFlightData.start" ])

if (length(bsSamples)==2) {
	idsWindow <- bsSamples[1];
      corPercent <- bsSamples[2];
} else {
idsWindow <- 199; ## default for Bruker, 500 MHz sampling
corPercent <- 0.05; 
if ( bsSamples[1] <= dataStart )
{
	userError( "'Baseline Start Time' must be greater than 'dataStart' in dbSelector." );
}
}

# Get the input data info from tofList object
spectraName <- names(tofList)
spectraCount <- length(tofList)
#print("Spectra, modelChoice, coef(z)[2], coef(z)[3], coef(z)[1], cOffset")
for (j in 1:spectraCount) {
    spectrum <- tofList[[j]]
    timeLength <- length(spectrum)
	if ( timeLength <= offsetSamples[length(offsetSamples)]  )
	{
		userError( "Time range for offset calculation must be less than the record length." );
	}
        ## estimate and subtract baseline & offset

	t0 <- dataStart:(timeLength+dataStart-1);
	cOffset<-mean(spectrum[offsetSamples]);
	bsTrendModel <- 0*spectrum;

	#cSTD<-sd(spectrum[(timeLength-2000):(timeLength-100)]);
	if (length(bsSamples)<4) {
		#print(idsWindow) ## debug
		envSignl<-findBsEnv(spectrum-cOffset, idsWindow);
		tSamples<-envSignl[,1]; sigSamples<-envSignl[,2];
		#print(tSamples) ## debug
	} else {
		tSamples<- bsSamples ## locations of signal-free regions for baseline model
		sigSamples <- spectrum[tSamples-dataStart]-cOffset;
	}

	# negLength <- length(which(sigSamples <=0));
	if (modelChoice == 1) {
		z <- line(tSamples, sigSamples);
		bsTrendModel <- coef(z)[1]+t0*coef(z)[2];
	} else if (modelChoice == 2) {
		sigSamples[sigSamples < 1] <- 1;
		z<-line(tSamples, log(sigSamples));
		bsTrendModel<- exp(coef(z)[1])*exp(t0*coef(z)[2]);
	} else { ## Gaussian model
		sigSamples[sigSamples < 1] <- 1;
		a1 <- matrix(tSamples, nrow=length(tSamples), ncol=2);
		a1[,1] <- tSamples^2; a1[,2] <- tSamples;
		z<-lsfit(a1,log(sigSamples), intercept=TRUE);
		bsTrendModel<-exp(coef(z)[2]*t0*t0+coef(z)[3]*t0+coef(z)[1]);
	}

#print(c(spectraName[j], modelChoice, coef(z)[2], coef(z)[3], coef(z)[1], cOffset))
	bsTrendModel<-(1-corPercent)*bsTrendModel;

	spectrum<-spectrum-cOffset-bsTrendModel; ## subtract baseline

    ## update the tofList record
    	tofList[spectraName[j]] <- list(spectrum)
}

return(tofList);

}


chargeBkgr_call <- function(tofList, offsetSamples, amp, decay) 
{

##----------------------------------------------------
## Calculate baseline from charge accumulation
##
## INPUT:
## tofList - list of the observed mass spectra (raw TOF)
## amp - amplitude fraction of accumulated signal
## decay - decay constant for the accumulated signal -- exp(-1/tau)
## offsetSamples –- TOF interval for the estimate of constant offset, e.g.:
##                 c(120000:122000) 
##
## OUTPUT:
## tofList - baseline subtracted lsi of mass spectra
##
## NOTES:
## --1-- Created by Prof.Cooke (wecook@wm.edu) 10/10/04
## --2-- No self-optimization code is included yet, but defaults work
## Slow with long spectra (> 50K points). Faster version with C-function 
## is available
#
# USAGE:
# tofList_noBS <- chargeBkgr_call(tofList,offsetSamples, amp, decay);
# Dependency: none
#

spectraName <- names(tofList)
spectraCount <- length(tofList)

for (j in 1:spectraCount) {
    spectrum <- tofList[[j]]
    lensp <- length(spectrum)
    basecon <- mean(spectrum[offsetSamples]);
    ytemp <- rep(0,lensp);
#amp=0.0005; decay = 0.9988; basecon = 2114;
#amp=0.0007; decay = 0.9989; basecon = 2140;
    for (i in 2:lensp) {
	ytemp[i] <- (1-amp)*decay*(ytemp[i-1])+amp*(spectrum[i]-basecon);
    }
#print(length(y)) #debug
#print(length(ytemp))
    ybs <- ytemp + basecon;
    spectrum <- spectrum-ybs; ## subtract baseline
    tofList[spectraName[j]] <- list(spectrum); ## update the tofList record
}

return(tofList);
}

offsetBkgr_call <- function(tofList, offsetSamples) 
{

##----------------------------------------------------
## Calculate baseline from ADC offset
##
## INPUT:
## tofList - list of the observed mass spectra (raw TOF)
## offsetSamples –- TOF interval for the estimate of constant offset, e.g.:
##                 c(120000:122000) 
##
## OUTPUT:
## tofList - baseline subtracted lsi of mass spectra
##
## NOTES:
## --1-- Created by Prof.Cooke (wecook@wm.edu) 10/10/04
## --2-- No self-optimization code is included yet, but defaults work
## Slow with long spectra (> 50K points). Faster version with C-function 
## is available
#
# USAGE:
# tofList_noBS <- chargeBkgr_call(tofList,offsetSamples);
# Dependency: none
#

spectraName <- names(tofList)
spectraCount <- length(tofList)

for (j in 1:spectraCount) {
    spectrum <- tofList[[j]]
    lensp <- length(spectrum)
    basecon <- mean(spectrum[offsetSamples]);
    spectrum <- spectrum-basecon; ## subtract baseline
    tofList[spectraName[j]] <- list(spectrum); ## update the tofList record
}

return(tofList);
}


pedRmMAV_call <- function(tofList,pw,nf){

###------------------------------------------------------------------ 
### Function pedRmMAV_call manages pedestal removal 
###  It performs the following:
###	Estimates the baseline using local minima
###	Calls mavSmoothing to smooth baseline estimate
###	Subtracts smooth baseline estimate from signal
###
###  INPUT:
###  tofList
### 	pw - peak half-width (used for window width for MAV ~6*HWFM is a 
###       good start) Assumes that peaks are the same width – valid for 
###       resampled data
###	nf – number/fraction of width in window
###
###  OUTPUT:
### 	tofList with pedestals removed
###
###  DEPENDENCIES:
###	mavSmoothing function (from "LibIDSfilts")
###
###------------------------------------------------------------------
ww<-floor(nf*pw)
#print(c(ww,pw))
numSpec<-length(tofList)
lastPt<-length(tofList[[1]])
if (ww < 2) { ww <- 2 }

for (ispec in 1:numSpec) {

sigt<-tofList[[ispec]]
bs0<-array(0,(lastPt-ww-ww))
 for (i in (ww+1):(lastPt-ww)) {
	bs0[i-ww]<-min(sigt[(i-ww):(i+ww)])
  }

zeroVector<-array(0,ww)
bs0<-c(zeroVector+min(sigt[1:ww]),bs0,zeroVector+min(sigt[(lastPt-ww):lastPt]))

bst <- mavSmoothing(bs0,ww)
sigpr <- sigt-bst

#sigpr[1:(ww)]<-0 # Set first ww points to zero if desired

tofList[[ispec]]<-sigpr
}
return(tofList)
}

