R/baseline.shirley.R

In baseline: Baseline Correction of Spectra

#' @title Shirley Background Estimation
#' 
#' @description Shirley Background correction for X-ray Photoelectron Spectroscopy.
#' 
#' @details The shape of the spectrum background or baseline is affected by inelastic
#' energy loss processes, secondary electrons and nearby peaks. A reasonable
#' approximation is essential for a qualitative and quantitative analysis of
#' XPS data especially if several components interfere in one spectrum. The
#' choice of an adequate background model is determined by the physical and
#' chemical conditions of the measurements and the significance of the
#' background to the information to be obtained. The subtraction of the
#' baseline before entering the fit iterations or the calculation of the peak
#' area can be an acceptable approximation for simple analytical problems. In
#' order to obtain chemical and physical parameters in detail, however, it is
#' absolutely necessary to include the background function in the iterative
#' peak fit procedure.  The primary function F(E) results from the
#' experimentally obtained function M(E) and the background function U(E) as
#' \deqn{F(E) = M(E)-U(E)}
#' 
#' The kinetic energy E of the spectra can be described as
#' 
#' \deqn{E = SE + SW * (i-1)}
#' 
#' SE means the start energy in eV, SW is the step width in eV and i the
#' channel number. i can assume values between 1 and N with N as the number of
#' data points.
#' 
#' In case of baseline calculation before initiating the fit procedure, the
#' background is set to the averaged experimental function M(E) in a sector
#' around the chosen start and end channels. With \eqn{i{_1}} as left channel
#' (\eqn{E{_1}}: low energy side) and \eqn{i{_2}} as right channel
#' (\eqn{E{_2}}: high energy side) the simulation of the baseline is obtained
#' as
#' 
#' \deqn{U(E_{1})=M(E_{1})}
#' 
#' and
#' 
#' \deqn{U(E_{2})=M(E_{2})}
#' 
#' If ZAP is the number of points used for averaging (can be set in the
#' preferences), the intensity of the averaged measuring function at the low
#' energy side is calculated by
#' 
#' \deqn{M(i_{1})=\frac{\sum_{i=0}^{ZAP-1}M(i_{1}+i)}{ZAP}}
#' 
#' and at the high energy side by
#' \deqn{M(i_{2})=\frac{\sum_{i=0}^{ZAP-1}M(i_{2}+i)}{ZAP}}
#' 
#' In many cases the Shirley model turned out to be a successful approximation
#' for the inelastic background of core level peaks of buried species, which
#' suffered significantly from inelastic losses of the emitted photoelectrons.
#' The calculation of the baseline is an iterative procedure. The number of
#' iteration cycles should be chosen high enough so that the shape of the
#' obtained background function does not change anymore. The analytical
#' expression for the Shirley background is
#' 
#' \deqn{U(E)= \int_{E}^{\infty}F(E')dE'+c}
#' 
#' The algorithm of Proctor and Sherwood ([1] A. Proctor, P.M.A. Sherwood,
#' Anal. Chem. 54 (1982) 13) is based on the assumption that for every point of
#' the spectrum the background intensity generated by a photoelectron line is
#' proportional to the number of all photoelectrons with higher kinetic energy.
#' The intensity of the background U(i) in channel i is given by
#' 
#' \deqn{U(i)=\frac{(a-b)Q(i)}{P(i)+Q(i)}+b}
#' 
#' where a and b are the measured intensities in channel \eqn{i{_1}} and
#' \eqn{i{_2}}, respectively, and P(i) and Q(i) represent the effective peak
#' areas to lower and higher kinetic energies relative to the channel under
#' consideration. An iterative procedure is necessary because P, Q, and U(i)
#' are unknown. In first approximation U(i) = b is used.
#' 
#' The function \code{baseline.shirley} implements the \code{shirley baseline}.
#' It is an iterative algorithm. The iteration stops when the deviation between
#' two consequent iteration is lower than \code{err} or when the max number of
#' iterations \code{maxit} is reached.
#' 
#' @aliases baseline.shirley shirley
#' @param spectra matrix with only 1 y-coordinates by rows (i.e.: y =
#' spectra[1,])
#' @param t Optional vector of spectrum abscissa
#' @param limits list with the y coordinates between calculation of background.
#' Ususally these are the extreme point of the data range.
#' @param maxit max number of iteration
#' @param err Tolerance of difference between iterations
#' @return The baseline function return an object of class \code{baseline}.
#' @seealso \link{baseline}
#' @references A. Proctor, P.M.A. Sherwood, Anal. Chem. 54 (1982) 13.
#' @examples
#' 
#' data("O1s")
#' Data <- O1s
#' 
#' ## The same example with C1s data 
#' # data("C1s")
#' # Data <- C1s
#' 
#' Y <- Data[2,, drop = FALSE]
#' X <- Data[1,]
#' 
#' corrected <- baseline(Y, method = "shirley", t = X)
#' plot(corrected, rev.x = TRUE, labels = X)
#' 
#' \dontrun{
#' # Dependent on external software
#' baselineGUI(Y, labels=X, method="shirley")
#' }
#' @export
baseline.shirley <- function(spectra, t = NULL, limits = NULL, maxit = 50, err = 1e-6) {
# INPUT:
# spectra - matrix with y in first row (y = spectra[1,])
# t       - vector with x coordinates (OPTIONAL)
# limits  - spectrum limits
# maxit   - maximum number of iterations
# err     - cut-off error
#
# OUTPUT:
# baseline  - proposed baseline
# corrected - baseline corrected spectra

   	object <- spectra[1,]
	npt <- length(object)
    if ( is.null(t) ) t <- 1:npt
    
	if ( is.null(limits) ) {
		limits <- list( x = c(t[1], t[npt]), y = c(spectra[1,1], spectra[1,npt]) )
		}
	
	lowLim 	<- min(limits$y)
	BGND 	<- rep.int(lowLim, npt)
	SumBGND <- sum(BGND)
	SumRTF 	<- sum(object)
	RangeY 	<- diff(range(limits$y))
	
	if ( diff(limits$y) > 0 ) {
		nloop <- 0
		repeat {
			nloop <- nloop + 1
			for ( idx in rev(seq_len(npt))) {
				BGND[npt-idx+1] <- ( (RangeY/(SumRTF-sum(BGND)))*(sum(object[npt:idx])-sum(BGND[npt:idx])) ) + lowLim	
			}
			if ( ( abs( (sum(BGND)-SumBGND)/SumBGND ) < err ) || nloop > maxit ){ break }				
			SumBGND <- abs(sum(BGND))
		}

	}
	else {	
		nloop <- 0
		repeat {
			nloop <- nloop + 1
			for ( idx in seq_len(npt)) {
				BGND[idx] <- ( (RangeY/(SumRTF-sum(BGND)))*(sum(object[idx:npt])-sum(BGND[idx:npt])) ) + lowLim	
			}
			if ( ( abs( (sum(BGND)-SumBGND)/SumBGND ) < err ) || nloop > maxit ){ break }				
			SumBGND <- abs(sum(BGND))
		}
	}
		
	baseline <- matrix(data = BGND, nrow = 1)
	
	list(baseline = baseline, corrected = spectra - baseline)
}