#' @title Calculate `Lx/Tx` ratio for CW-OSL curves
#'
#' @description
#' Calculate `Lx/Tx` ratios from a given set of CW-OSL curves assuming late light
#' background subtraction.
#'
#' @details
#' The integrity of the chosen values for the signal and background integral is
#' checked by the function; the signal integral limits have to be lower than
#' the background integral limits. If a [vector] is given as input instead
#' of a [data.frame], an artificial [data.frame] is produced. The
#' error calculation is done according to Galbraith (2002).
#'
#' **Please note:** In cases where the calculation results in `NaN` values (for
#' example due to zero-signal, and therefore a division of 0 by 0), these `NaN`
#' values are replaced by 0.
#'
#' **`sigmab`**
#'
#' The default value of `sigmab` is calculated assuming the background is
#' constant and **would not** applicable when the background varies as,
#' e.g., as observed for the early light subtraction method.
#'
#' **sig0**
#'
#' This argument allows to add an extra component of error to the final `Lx/Tx`
#' error value. The input will be treated as factor that is multiplied with
#' the already calculated `LxTx` and the result is add up by:
#'
#' \deqn{se(LxTx) = \sqrt(se(LxTx)^2 + (LxTx * sig0)^2)}
#'
#'
#' **background.count.distribution**
#'
#' This argument allows selecting the distribution assumption that is used for
#' the error calculation. According to Galbraith (2002, 2014) the background
#' counts may be overdispersed (i.e. do not follow a Poisson distribution,
#' which is assumed for the photomultiplier counts). In that case (might be the
#' normal case) it has to be accounted for the overdispersion by estimating
#' \eqn{\sigma^2} (i.e. the overdispersion value). Therefore the relative
#' standard error is calculated as:
#'
#' - `poisson`
#' \deqn{rse(\mu_{S}) \approx \sqrt(Y_{0} + Y_{1}/k^2)/Y_{0} - Y_{1}/k}
#'
#' - `non-poisson`
#' \deqn{rse(\mu_{S}) \approx \sqrt(Y_{0} + Y_{1}/k^2 + \sigma^2(1+1/k))/Y_{0} - Y_{1}/k}
#'
#' **Please note** that when using the early background subtraction method in
#' combination with the 'non-poisson' distribution argument, the corresponding `Lx/Tx` error
#' may considerably increase due to a high `sigmab` value.
#' Please check whether this is valid for your data set and if necessary
#' consider to provide an own `sigmab` value using the corresponding argument `sigmab`.
#'
#' @param Lx.data [RLum.Data.Curve-class] or [data.frame] (**required**):
#' requires a CW-OSL shine down curve (x = time, y = counts)
#'
#' @param Tx.data [RLum.Data.Curve-class] or [data.frame] (*optional*):
#' requires a CW-OSL shine down curve (x = time, y = counts). If no
#' input is given the `Tx.data` will be treated as `NA` and no `Lx/Tx` ratio
#' is calculated.
#'
#' @param signal.integral [numeric] (**required**): vector with the limits for the signal integral.
#' Can be set to `NA` than now integrals are considered and all other integrals are set to `NA` as well.
#'
#' @param signal.integral.Tx [numeric] (*optional*):
#' vector with the limits for the signal integral for the `Tx`-curve. If nothing is provided the
#' value from `signal.integral` is used.
#'
#' @param background.integral [numeric] (**required**):
#' vector with the bounds for the background integral.
#' Can be set to `NA` than now integrals are considered and all other integrals are set to `NA` as well.
#'
#' @param background.integral.Tx [numeric] (*optional*):
#' vector with the limits for the background integral for the `Tx` curve.
#' If nothing is provided the value from `background.integral` is used.
#'
#' @param background.count.distribution [character] (*with default*):
#' sets the count distribution assumed for the error calculation.
#' Possible arguments `poisson` or `non-poisson`. See details for further information
#'
#' @param use_previousBG [logical] (*with default*):
#' If set to `TRUE` the background of the `Lx`-signal is subtracted also
#' from the `Tx`-signal. Please note that in this case separate
#' signal integral limits for the `Tx`-signal are not allowed and will be reset.
#'
#' @param sigmab [numeric] (*optional*):
#' option to set a manual value for the overdispersion (for `LnTx` and `TnTx`),
#' used for the `Lx/Tx` error calculation. The value should be provided as
#' absolute squared count values, e.g. `sigmab = c(300,300)`.
#' **Note:** If only one value is provided this value is taken for both (`LnTx` and `TnTx`) signals.
#'
#' @param sig0 [numeric] (*with default*):
#' allow adding an extra component of error to the final `Lx/Tx` error value
#' (e.g., instrumental error, see details).
#'
#' @param digits [integer] (*with default*):
#' round numbers to the specified digits.
#' If digits is set to `NULL` nothing is rounded.
#'
#' @return
#' Returns an S4 object of type [RLum.Results-class].
#'
#' Slot `data` contains a [list] with the following structure:
#'
#' **@data**
#' ```
#' $LxTx.table (data.frame)
#' .. $ LnLx
#' .. $ LnLx.BG
#' .. $ TnTx
#' .. $ TnTx.BG
#' .. $ Net_LnLx
#' .. $ Net_LnLx.Error
#' .. $ Net_TnTx
#' .. $ Net_TnTx.Error
#' .. $ LxTx
#' .. $ LxTx.Error
#' $ calc.parameters (list)
#' .. $ sigmab.LnTx
#' .. $ sigmab.TnTx
#' .. $ k
#' ```
#'
#' **@info**
#' ```
#' $ call (original function call)
#' ```
#'
#' @note
#' The results of this function have been cross-checked with the Analyst
#' (version 3.24b). Access to the results object via [get_RLum].
#'
#' **Caution:** If you are using early light subtraction (EBG), please either provide your
#' own `sigmab` value or use `background.count.distribution = "poisson"`.
#'
#' @section Function version: 0.8.0
#'
#' @author
#' Sebastian Kreutzer, Institute of Geography, Heidelberg University (Germany)
#'
#' @seealso [RLum.Data.Curve-class], [Analyse_SAR.OSLdata], [plot_GrowthCurve],
#' [analyse_SAR.CWOSL]
#'
#' @references Duller, G., 2018. Analyst v4.57 - User Manual.
#' `https://users.aber.ac.uk/ggd`\cr
#'
#' Galbraith, R.F., 2002. A note on the variance of a background-corrected OSL
#' count. Ancient TL, 20 (2), 49-51.
#'
#' Galbraith, R.F., 2014. A further note on the variance of a
#' background-corrected OSL count. Ancient TL, 31 (2), 1-3.
#'
#' @keywords datagen
#'
#' @examples
#'
#' ##load data
#' data(ExampleData.LxTxOSLData, envir = environment())
#'
#' ##calculate Lx/Tx ratio
#' results <- calc_OSLLxTxRatio(
#' Lx.data = Lx.data,
#' Tx.data = Tx.data,
#' signal.integral = c(1:2),
#' background.integral = c(85:100))
#'
#' ##get results object
#' get_RLum(results)
#'
#' @md
#' @export
calc_OSLLxTxRatio <- function(
Lx.data,
Tx.data = NULL,
signal.integral,
signal.integral.Tx = NULL,
background.integral,
background.integral.Tx = NULL,
background.count.distribution = "non-poisson",
use_previousBG = FALSE,
sigmab = NULL,
sig0 = 0,
digits = NULL
){
# Test input data ---------------------------------------------------------
##(1) - integrity checks
if(!is.null(Tx.data)){
##(a) - check data type
if(is(Lx.data)[1]!=is(Tx.data)[1]){
stop("[calc_OSLLxTxRatio()] Data type of Lx and Tx data differs!",
call. = FALSE)
}
##(b) - test if data.type is valid in general
if(inherits(Lx.data, "RLum.Data.Curve")){
Lx.data <- as(Lx.data, "data.frame")
Tx.data <- as(Tx.data, "data.frame")
}else{
##go further
if((is(Lx.data)[1] != "data.frame" &
is(Lx.data)[1] != "numeric") &
is(Lx.data)[1] != "matrix"){
stop("[calc_OSLLxTxRatio()] Data type error! Required types are data.frame or numeric vector.", call. = FALSE)
}
}
##(c) - convert vector to data.frame if necessary
if(is(Lx.data)[1] != "data.frame" &
is(Lx.data)[1] != "matrix"){
Lx.data <- data.frame(x = 1:length(Lx.data), y = Lx.data)
Tx.data <- data.frame(x = 1:length(Tx.data), y = Tx.data)
}
##(d) - check if Lx and Tx curves have the same channel length
if(length(Lx.data[,2]) != length(Tx.data[,2]))
stop("[calc_OSLLxTxRatio()] Channel numbers of Lx and Tx data differ!",
call. = FALSE)
}else{
Tx.data <- data.frame(x = NA,y = NA)
##support RLum.objects
if(inherits(Lx.data, "RLum.Data.Curve")){
Lx.data <- as(Lx.data, "data.frame")
}
##check for matrix
if(is(Lx.data)[1] == "matrix"){
Lx.data <- as.data.frame(Lx.data)
}
##no it should be a data.frame, if not, try to produce one
if(is(Lx.data)[1]!="data.frame") {
Lx.data <- data.frame(x = 1:length(Lx.data),y = Lx.data)
}
}#endif::missing Tx.data
# Alternate mode ----------------------------------------------------------
if(any(is.na(c(signal.integral, background.integral)))){
signal.integral <- background.integral <- NA
LnLx <- sum(Lx.data[,2])
TnTx <- sum(Tx.data[,2])
LnLxTnTx <- data.frame(
LnLx = LnLx,
LnLx.BG = 0,
TnTx = TnTx,
TnTx.BG = 0,
Net_LnLx = LnLx,
Net_LnLx.Error = 0,
Net_TnTx = TnTx,
Net_TnTx.Error = 0,
LxTx = LnLx/TnTx,
LxTx.Error = 0)
return(set_RLum(
class = "RLum.Results",
data = list(
LxTx.table = LnLxTnTx,
calc.parameters = NULL,
info = list(call = sys.call())
)))
}
# Continue checks ---------------------------------------------------------
##(e) - check if signal integral is valid
if(min(signal.integral) < 1 | max(signal.integral>length(Lx.data[,2]))){
stop("[calc_OSLLxTxRatio()] signal.integral is not valid!", call. = FALSE)}
##(f) - check if background integral is valid
if(min(background.integral)<1 | max(background.integral>length(Lx.data[,2]))){
stop(paste("[calc_OSLLxTxRatio()] background.integral is not valid! Max: ",length(Lx.data[,2]),sep=""), call. = FALSE)}
##(g) - check if signal and background integral overlapping
if(min(background.integral)<=max(signal.integral)){
stop("[calc_OSLLxTxRatio()] Overlapping of 'signal.integral' and 'background.integral' is not permitted!", call. = FALSE)}
##(h) - similar procedure for the Tx limits
if(all(c(!is.null(signal.integral.Tx),!is.null(background.integral.Tx)))){
if(use_previousBG){
warning("[calc_OSLLxTxRatio()] For option use_previousBG = TRUE independent Lx and Tx integral limits are not allowed. Integral limits of Lx used for Tx.", call. = FALSE)
signal.integral.Tx <- signal.integral
background.integral.Tx <- background.integral
}
if(min(signal.integral.Tx) < 1 | max(signal.integral.Tx>length(Tx.data[,2]))){
stop("[calc_OSLLxTxRatio()] signal.integral.Tx is not valid!", call. = FALSE)}
if(min(background.integral.Tx)<1 | max(background.integral.Tx>length(Tx.data[,2]))){
stop(paste("[calc_OSLLxTxRatio()] background.integral.Tx is not valid! Max: ",length(Tx.data[,2]),sep=""),
call. = FALSE)}
if(min(background.integral.Tx)<=max(signal.integral.Tx)){
stop("[calc_OSLLxTxRatio()] Overlapping of 'signal.integral.Tx' and 'background.integral.Tx' is not permitted!",
call. = FALSE)}
}else if(!all(c(is.null(signal.integral.Tx),is.null(background.integral.Tx)))){
stop("[calc_OSLLxTxRatio()] You have to provide both: signal.integral.Tx and background.integral.Tx!",
call. = FALSE)
}else{
signal.integral.Tx <- signal.integral
background.integral.Tx <- background.integral
}
##check sigmab
if (!is.null(sigmab)) {
if (!is(sigmab, "numeric")) {
stop("[calc_OSLLxTxRatio()] 'sigmab' has to be of type numeric.", call. = FALSE)
}
if (length(sigmab) > 2) {
stop("[calc_OSLLxTxRatio()] Maximum allowed vector length for 'sigmab' is 2.", call. = FALSE)
}
}
##--------------------------------------------------------------------------##
##(2) - read data and produce background subtracted values
## calculate k value - express the background as mutiple value from the number
## of signal integral channels, however, it can be < 1 also
n <- length(signal.integral)
m <- length(background.integral)
k <- m/n
n.Tx <- length(signal.integral.Tx)
##use previous BG and account for the option to set different integral limits
if(use_previousBG){
m.Tx <- m
}else{
m.Tx <- length(background.integral.Tx)
}
k.Tx <- m.Tx/n.Tx
##LnLx (comments are corresponding variables to Galbraith, 2002)
Lx.curve <- Lx.data[,2]
Lx.signal <- sum(Lx.curve[signal.integral]) #Y.0
Lx.background <- sum(Lx.curve[background.integral]) #Y.1
Lx.background <- Lx.background*1/k #mu.B
LnLx <- Lx.signal - Lx.background
##TnTx
Tx.curve <- ifelse(is.na(Tx.data[,1])==FALSE, Tx.data[,2], NA)
Tx.signal <- sum(Tx.curve[signal.integral.Tx])
##use previous BG
if(use_previousBG){
Tx.background <- Lx.background
}else{
Tx.background <- sum(Tx.curve[background.integral.Tx])*1/k.Tx
}
TnTx <- (Tx.signal-Tx.background)
##--------------------------------------------------------------------------##
##(3)
## calculate Lx/Tx Errors according Galbraith (2002) and the personal
## communication of Galbraith (2014) via e-mail
## Nomenclature as stated in the articles
##(a)
## set Y.0 (sum OSL signal including the background) and
## Y.1 (total counts over m later channels)
Y.0 <- Lx.signal
Y.0_TnTx <- Tx.signal
Y.1 <- sum(Lx.curve[background.integral])
Y.1_TnTx <- sum(Tx.curve[background.integral.Tx])
##(b) estimate overdispersion (here called sigmab), see equation (4) in
## Galbraith (2002), Galbraith (2014)
## If else condition for the case that k < 2
if(round(k,digits = 1) >= 2 & ((min(background.integral) + length(signal.integral)*(2+1)) <= length(Lx.curve))){
##(b)(1)(1)
## note that m = n*k = multiple of background.integral from signal.integral
Y.i <- vapply(0:round(k,digits=0), function(i){
sum(Lx.curve[
(min(background.integral)+length(signal.integral)*i):
(min(background.integral)+length(signal.integral)+length(signal.integral)*i)])
}, FUN.VALUE = vector(mode = "numeric", length = 1L))
Y.i <- na.exclude(Y.i)
sigmab.LnLx <- abs(var(Y.i) - mean(Y.i)) ##sigmab is denoted as sigma^2 = s.Y^2-Y.mean
##therefore here absolute values are given
}else{
## provide warning if m is < 25, as suggested by Rex Galbraith
## low number of degree of freedom
if (m < 25) {
warning(
"[calc_OSLLxTxRatio()] Number of background channels for Lx < 25; error estimation might not be reliable!",
call. = FALSE)
}
sigmab.LnLx <- abs((var(Lx.curve[background.integral]) -
mean(Lx.curve[background.integral])) * n)
}
if (round(k.Tx, digits = 1) >= 2 &
((
min(background.integral.Tx) + length(signal.integral.Tx) * (2 + 1)
) <= length(Tx.curve))) {
##(b)(1)(1)
## note that m.Tx = n.Tx*k.Tx = multiple of background.integral.Tx from signal.integral.Tx
## also for the TnTx signal
Y.i_TnTx <- vapply(0:round(k.Tx, digits = 0), function(i) {
sum(Tx.curve[(min(background.integral.Tx) + length(signal.integral.Tx) *
i):(
min(background.integral.Tx) + length(signal.integral.Tx) + length(signal.integral.Tx) *
i
)])
}, FUN.VALUE = vector(mode = "numeric", length = 1L))
Y.i_TnTx <- na.exclude(Y.i_TnTx)
sigmab.TnTx <- abs(var(Y.i_TnTx) - mean(Y.i_TnTx))
} else{
## provide warning if m is < 25, as suggested by Rex Galbraith
## low number of degree of freedom
if (m.Tx < 25 && use_previousBG == FALSE) {
warning("[calc_OSLLxTxRatio()] Number of background channels for Tx < 25; error estimation might not be reliable!",
call. = FALSE)
}
sigmab.TnTx <- abs((var(Tx.curve[background.integral.Tx]) -
mean(Tx.curve[background.integral.Tx])) * n.Tx)
}
##account for a manually set sigmab value
if (!is.null(sigmab)) {
if (length(sigmab) == 2) {
sigmab.LnLx <- sigmab[1]
sigmab.TnTx <- sigmab[2]
}else{
sigmab.LnLx <- sigmab[1]
sigmab.TnTx <- sigmab[1]
}
}
##(c)
## Calculate relative error of the background subtracted signal
## according to Galbraith (2002), equation (6) with changes
## from Galbraith (2014), equation 6
## Discussion with Rex Galbraith via e-mail (2014-02-27):
## Equation 6 is appropriate to be implemented as standard
if(background.count.distribution == "poisson"){
##(c.1) estimate relative standard error for assuming a poisson distribution
LnLx.relError <- sqrt((Y.0 + Y.1/k^2))/(Y.0-Y.1/k) ## rse(mu.s)
TnTx.relError <- sqrt((Y.0_TnTx + Y.1_TnTx/k^2))/(Y.0_TnTx-Y.1_TnTx/k)
}else{
##(c.2) estimate relative standard error for a non-poisson distribution
if(background.count.distribution != "non-poisson"){
warning("Unknown method for background.count.distribution. A non-poisson distribution is assumed!", call. = FALSE)}
LnLx.relError <- sqrt(Y.0 + Y.1/k^2 + sigmab.LnLx*(1+1/k))/
(Y.0 - Y.1/k)
TnTx.relError <- sqrt(Y.0_TnTx + Y.1_TnTx/k^2 + sigmab.TnTx*(1+1/k))/
(Y.0_TnTx - Y.1_TnTx/k)
}
##(d)
##calculate absolute standard error
LnLx.Error <- abs(LnLx*LnLx.relError)
TnTx.Error <- abs(TnTx*TnTx.relError)
##we do not want to have NaN values, as they are mathematically correct, but make
##no sense and would result in aliquots that become rejected later
if(is.nan(LnLx.Error)) LnLx.Error <- 0
if(is.nan(TnTx.Error)) TnTx.Error <- 0
##combine results
LnLxTnTx <- cbind(
Lx.signal,
Lx.background,
Tx.signal,
Tx.background,
LnLx,
LnLx.Error,
TnTx,
TnTx.Error
)
##--------------------------------------------------------------------------##
##(4) Calculate LxTx error according Galbraith (2014)
#transform results in a data.frame
LnLxTnTx <- as.data.frame((LnLxTnTx))
#add col names
colnames(LnLxTnTx)<-c("LnLx", "LnLx.BG",
"TnTx", "TnTx.BG",
"Net_LnLx", "Net_LnLx.Error",
"Net_TnTx", "Net_TnTx.Error")
##calculate Ln/Tx
LxTx <- LnLxTnTx$Net_LnLx/LnLxTnTx$Net_TnTx
##set NaN
if(is.nan(LxTx)) LxTx <- 0
##calculate Ln/Tx error
LxTx.relError <- sqrt(LnLx.relError^2 + TnTx.relError^2)
LxTx.Error <- abs(LxTx * LxTx.relError)
##set NaN
if(is.nan(LxTx.Error)) LxTx.Error <- 0
##add an extra component of error
LxTx.Error <- sqrt(LxTx.Error^2 + (sig0 * LxTx)^2)
##return combined values
temp <- cbind(LnLxTnTx,LxTx,LxTx.Error)
##apply digits if wanted
if(!is.null(digits)){
temp[1,] <- round(temp[1,], digits = digits)
}
calc.parameters <- list(
sigmab.LnLx = sigmab.LnLx,
sigmab.TnTx = sigmab.TnTx,
k = k)
##set results object
return(set_RLum(
class = "RLum.Results",
data = list(
LxTx.table = temp,
calc.parameters = calc.parameters),
info = list(call = sys.call())
))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.