Nothing
#' @title Al2O3:C Passive Dosimeter Measurement Analysis
#'
#' @description The function provides the analysis routines for measurements on a
#' FI lexsyg SMART reader using Al2O3:C chips according to Kreutzer et al., 2018
#'
#' @details
#'
#' **Working with a travel dosimeter**
#'
#' The function allows to define particular aliquots as travel dosimeters. For example:
#' `travel_dosimeter = c(1,3,5)` sets aliquots 1, 3 and 5 as travel dosimeters. These dose values
#' of this dosimeters are combined and automatically subtracted from the obtained dose values
#' of the other dosimeters.
#'
#' **Calculate TL dose**
#'
#' The argument `calculate_TL_dose` provides the possibility to experimentally calculate a TL-dose,
#' i.e. an apparent dose value derived from the TL curve ratio. However, it should be noted that
#' this value is only a fall back in case something went wrong during the measurement of the optical
#' stimulation. The TL derived dose value is corrected for cross-talk and for the irradiation time,
#' but not considered if a travel dosimeter is defined.
#'
#' Calculating the palaeodose is possible without **any TL** curve in the sequence!
#'
#' **Test parameters**
#'
#' `TL_peak_shift` [numeric] (default: `15`):
#'
#' Checks whether the TL peak shift is bigger > 15 K, indicating a problem with the
#' thermal contact of the chip.
#'
#' `stimulation_power` [numeric] (default: `0.05`):
#'
#' So far available, information on the delivered optical stimulation are compared. Compared are
#' the information from the first curves with all others. If the ratio differs more from
#' unity than the defined by the threshold, a warning is returned.
#'
#' @param object [RLum.Analysis-class] (**required**): measurement input
#'
#' @param signal_integral [numeric] (*optional*): signal integral, used for the signal
#' and the background. Example: `c(1:10)` for the first 10 channels.
#' If nothing is provided the full range is used
#'
#' @param dose_points [numeric] (*with default*):
#' vector with dose points, if dose points are repeated, only the general
#' pattern needs to be provided. Default values follow the suggestions
#' made by Kreutzer et al., 2018
#'
#' @param recordType [character] (*with default*): input curve selection, which is passed to
#' function [get_RLum]. To deactivate the automatic selection set the argument to `NULL`
#'
#' @param irradiation_time_correction [numeric] or [RLum.Results-class] (*optional*):
#' information on the used irradiation time correction obtained by another experiments.
#' I a `numeric` is provided it has to be of length two: mean, standard error
#'
#' @param calculate_TL_dose [logical] (*with default*): Enables/disables experimental dose estimation
#' based on the TL curves. Taken is the ratio of the peak sums of each curves +/- 5 channels.
#'
#' @param cross_talk_correction [numeric] or [RLum.Results-class] (*optional*):
#' information on the used irradiation time correction obtained by another experiments.
#' If a `numeric` vector is provided it has to be of length three:
#' mean, 2.5 % quantile, 97.5 % quantile.
#'
#' @param travel_dosimeter [numeric] (*optional*): specify the position of the travel dosimeter
#' (so far measured a the same time). The dose of travel dosimeter will be subtracted from all
#' other values.
#'
#' @param test_parameters [list] (*with default*):
#' set test parameters. Supported parameters are: `TL_peak_shift` All input: [numeric]
#' values, `NA` and `NULL` (s. Details)
#'
#' @param verbose [logical] (*with default*):
#' enable/disable verbose mode
#'
#' @param plot [logical] (*with default*): enable/disable plot output, if `object` is of type [list],
#' a [numeric] vector can be provided to limit the plot output to certain aliquots
#'
#' @param ... further arguments that can be passed to the plot output, supported are `norm`, `main`, `mtext`,
#' `title` (for self-call mode to specify, e.g., sample names)
#'
#' @return Function returns results numerically and graphically:
#'
#' -----------------------------------\cr
#' `[ NUMERICAL OUTPUT ]`\cr
#' -----------------------------------\cr
#'
#' **`RLum.Results`**-object
#'
#' **slot:** **`@data`**
#'
#' \tabular{lll}{
#' **Element** \tab **Type** \tab **Description**\cr
#' `$data` \tab `data.frame` \tab the estimated equivalent dose \cr
#' `$data_table` \tab `data.frame` \tab full dose and signal table \cr
#' `test_parameters` \tab `data.frame` \tab results with test parameters \cr
#' `data_TDcorrected` \tab `data.frame` \tab travel dosimeter corrected results (only if TD was provided)\cr
#' }
#'
#' *Note: If correction the irradiation time and the cross-talk correction method is used, the De
#' values in the table `data` table are already corrected, i.e. if you want to get an uncorrected value,
#' you can use the column `CT_CORRECTION` remove the correction*
#'
#'**slot:** **`@info`**
#'
#' The original function call
#'
#' ------------------------\cr
#' `[ PLOT OUTPUT ]`\cr
#' ------------------------\cr
#'
#' - OSL and TL curves, combined on two plots.
#'
#'
#' @section Function version: 0.2.6
#'
#' @author Sebastian Kreutzer, Institute of Geography, Heidelberg University (Germany)
#'
#' @seealso [analyse_Al2O3C_ITC]
#'
#' @references
#'
#' Kreutzer, S., Martin, L., Guérin, G., Tribolo, C., Selva, P., Mercier, N., 2018. Environmental Dose Rate
#' Determination Using a Passive Dosimeter: Techniques and Workflow for alpha-Al2O3:C Chips.
#' Geochronometria 45, 56-67.
#'
#' @keywords datagen
#'
#' @examples
#' ##load data
#' data(ExampleData.Al2O3C, envir = environment())
#'
#' ##run analysis
#' analyse_Al2O3C_Measurement(data_CrossTalk)
#'
#' @md
#' @export
analyse_Al2O3C_Measurement <- function(
object,
signal_integral = NULL,
dose_points = c(0,4),
recordType = c("OSL (UVVIS)", "TL (UVVIS)"),
calculate_TL_dose = FALSE,
irradiation_time_correction = NULL,
cross_talk_correction = NULL,
travel_dosimeter = NULL,
test_parameters = NULL,
verbose = TRUE,
plot = TRUE,
...
){
# Self call -----------------------------------------------------------------------------------
if(is(object, "list")){
if(!all(unlist(lapply(object, function(x){is(x, "RLum.Analysis")})))){
stop("[analyse_Al2O3C_Measurement()] The elements in 'object' are not all of type 'RLum.Analsyis'", call. = FALSE)
}
##expand input arguments
if(!is.null(signal_integral)){
signal_integral <- rep(list(signal_integral), length = length(object))
}
##dose points
if(is(dose_points, "list")){
dose.points <- rep(dose_points, length = length(object))
}else{
dose_points <- rep(list(dose_points), length = length(object))
}
##irradiation time correction
if(is(irradiation_time_correction, "list")){
irradiation_time_correction <- rep(irradiation_time_correction, length = length(object))
}else{
irradiation_time_correction <- rep(list(irradiation_time_correction), length = length(object))
}
##cross talk correction
if(is( cross_talk_correction, "list")){
cross_talk_correction <- rep( cross_talk_correction, length = length(object))
}else{
cross_talk_correction <- rep(list( cross_talk_correction), length = length(object))
}
##test_parameters
if(is(test_parameters[[1]], "list")){
test_parameters <- rep(test_parameters, length = length(object))
}else{
test_parameters <- rep(list(test_parameters), length = length(object))
}
##verbose
##plot
if(is(plot, "logical")){
plot <- rep(x = plot, length(object))
}else{
plot <- 1:length(object)%in%plot
}
##run analyis
results <- lapply(1:length(object), function(x) {
temp <- analyse_Al2O3C_Measurement(
object = object[[x]],
signal_integral = signal_integral[[x]],
dose_points = dose_points[[x]],
irradiation_time_correction = irradiation_time_correction[[x]],
cross_talk_correction = cross_talk_correction[[x]],
test_parameters = test_parameters[[x]],
calculate_TL_dose = calculate_TL_dose,
verbose = verbose,
plot = plot[x],
...
)
##adjusting the terminal output, to avoid confusions
if(verbose)
cat(" ... (#",x, " | ALQ POS: ", temp$data$POSITION,")\n", sep = "")
##add running number to the plot, but only of we had a plot here...
if(plot[x]){
title(main = paste0(list(...)$title[x], " ","#", x), adj = 1, line = 3)
}
return(temp)
})
##merge results
results <- merge_RLum(results)
##correct sys.call, otherwise it gets a little bit strange
##why this is not implemented in the merge_RLum() method ... because here it would be wrong!
results@info[names(results@info) == "call"] <- NULL
results@info$call <- sys.call()
##travel dosimeter
##check for travel dosimeter and subtract the values so far this is meaningful at all
if(!is.null(travel_dosimeter)){
##check data type
if(!is(travel_dosimeter, "numeric"))
stop("[analyse_Al2O3C_Measurement()] Input for `travel_dosimeter` is not numeric!",
call. = FALSE)
##check whether everything is subtracted from everything ... you never know, users do weird stuff
if(length(travel_dosimeter) == nrow(results$data))
try(stop("[analyse_Al2O3C_Measurement()] You specified every position as travel dosimeter, nothing corrected!",
call. = FALSE))
##check if the position is valid
if(!any(travel_dosimeter%in%results$data$POSITION))
try(stop("[analyse_Al2O3C_Measurement()] Invalid position in 'travel_dosimeter', nothing corrected!",
call. = FALSE))
##correct for the travel dosimeter calculating the weighted mean and the sd (as new error)
##if only one value is given just take it
if(length(travel_dosimeter) == 1 && nrow(results$data[travel_dosimeter==results$data$POSITION,c(1,2)]) == 1){
correction <- as.numeric(results$data[travel_dosimeter==results$data$POSITION,c(1,2)])
}else{
temp.correction <- results$data[results$data$POSITION%in%travel_dosimeter,c(1,2)]
correction <- c(
stats::weighted.mean(
x = temp.correction[[1]],
w = if(all(temp.correction[[2]]==0)){rep(1, length(temp.correction[[2]]))} else {temp.correction[[2]]}),
sd(temp.correction[,1]))
rm(temp.correction)
}
##subtract all the values, in a new data frame, we do not touch the original data
data_TDcorrected <- data.frame(
DE = results@data$data[!results$data$POSITION%in%travel_dosimeter,1] - correction[1],
DE_ERROR = sqrt(results@data$data[!results$data$POSITION%in%travel_dosimeter,2]^2 + correction[2]^2),
POSITION = results@data$data[!results$data$POSITION%in%travel_dosimeter, "POSITION"]
)
##however, we set information on the travel dosimeter in the corresponding column
results@data$data$TRAVEL_DOSIMETER <- results$data$POSITION%in%travel_dosimeter
##attach the new element to the results output
results@data <- c(results@data, list(data_TDcorrected = data_TDcorrected))
##return message
if(verbose)
cat("\n ...+ travel dosimeter correction applied.\n ...+ results stored in object $data_TDcorrected.\n\n")
} ##end travel dosimeter
##return results
return(results)
}
# Integrity check ---------------------------------------------------------------------------
##TODO ... do more, push harder
##Add sufficient unit tests
# Preparation ---------------------------------------------------------------------------------
##select curves based on the recordType selection; if not NULL
if(!is.null(recordType)){
object_raw <- object
object <- get_RLum(object, recordType = recordType, drop = FALSE)
}
##set signal integral
if(is.null(signal_integral)){
signal_integral <- c(1:nrow(object[[1]][]))
}else{
##check whether the input is valid, otherwise make it valid
if(min(signal_integral) < 1 | max(signal_integral) > nrow(object[[1]][])){
signal_integral <- c(1:nrow(object[[1]][]))
warning(
paste0(
"[analyse_Al2O3C_Measurement()] Input for 'signal_integral' corrected to 1:", nrow(object[[1]][])
),
call. = FALSE
)
}
}
## Set Irradiation Time Correction ---------------
if (!is.null(irradiation_time_correction)) {
if (is(irradiation_time_correction, "RLum.Results")) {
if (irradiation_time_correction@originator == "analyse_Al2O3C_ITC") {
irradiation_time_correction <- get_RLum(irradiation_time_correction)
##consider the case for more than one observation ...
if(nrow(irradiation_time_correction)>1){
irradiation_time_correction <- c(mean(irradiation_time_correction[[1]]), sd(irradiation_time_correction[[1]]))
}else{
irradiation_time_correction <- c(irradiation_time_correction[[1]], irradiation_time_correction[[2]])
}
} else{
stop(
"[analyse_Al2O3C_Measurement()] The object provided for the argument 'irradiation_time_correction' was created by an unsupported function!",
call. = FALSE
)
}
}
}
## Set Cross Talk Correction ---------------
##check wehther the information on the position was stored in the input
##object
if(!is.null(get_RLum(object = object[[1]], info.object = "position"))){
POSITION <- get_RLum(object = object[[1]], info.object = "position")
}else{
message("[analyse_Al2O3_Measurement()] Aliquot position number was not found. No cross talk correction was applied!")
cross_talk_correction <- c(0,0,0)
POSITION <- NA
}
if(is.null(cross_talk_correction)){
cross_talk_correction <- c(0,0,0)
}else{
##check whether the input is of type RLum.Results and check orignator
if (is(cross_talk_correction, "RLum.Results") &&
cross_talk_correction@originator == "analyse_Al2O3C_CrossTalk") {
##grep cross talk correction and calculate values for
##this particular carousel position
cross_talk_correction <-
as.numeric(predict(cross_talk_correction$fit,
newdata = data.frame(x = POSITION),
interval = "confidence"))
}else{
stop(
"[analyse_Al2O3C_Measurement()] The object provided for the argument 'cross_talk_correction' was created by an unsupported function or has a wrong originator!",
call. = FALSE
)
}
}
# Calculation ---------------------------------------------------------------------------------
##we have two dose points, and one background curve, we do know only the 2nd dose
##set test parameters
test_parameters.default <- list(
TL_peak_shift = 15,
stimulation_power = 0.05
)
##modify default values by given input
if(!is.null(test_parameters)){
test_parameters <- modifyList(test_parameters.default, test_parameters)
##remove NULL elements from list
test_parameters <- test_parameters[!sapply(test_parameters, is.null)]
}else{
test_parameters <- test_parameters.default
}
##calculate integrated light values
NATURAL <- sum(get_RLum(object, recordType = "OSL")[[1]]@data[signal_integral, 2])
REGENERATED <- sum(get_RLum(object, recordType = "OSL")[[2]]@data[signal_integral, 2])
BACKGROUND <- sum(get_RLum(object, recordType = "OSL")[[3]]@data[signal_integral, 2])
##do the same for the TL
if (calculate_TL_dose[1] && any(grepl("TL", names(object)))){
NATURAL_TL <- try(sum(
object@records[[2]]@data[
(which.max(object@records[[2]]@data[,2])-5):(which.max(object@records[[2]]@data[,2])+5),2]), silent = TRUE)
REGENERATED_TL <- try(sum(
object@records[[4]]@data[
(which.max(object@records[[4]]@data[,2])-5):(which.max(object@records[[4]]@data[,2])+5),2]), silent = TRUE)
##catch errors if the integration fails
if(inherits(NATURAL_TL, "try-error")){
NATURAL_TL <- NA
warning("[analyse_Al2O3_Measurement()] Natural TL signal out of bounds, NA returned!", call. = FALSE, immediate. = TRUE)
}
if(inherits(REGENERATED_TL, "try-error")){
REGENERATED_TL <- NA
warning("[analyse_Al2O3_Measurement()] Regenerated TL signal out of bounds, NA returned!", call. = FALSE, immediate. = TRUE)
}
}else{
NATURAL_TL <- NA
REGENERATED_TL <- NA
}
##combine into data.frame
temp_df <- data.frame(
POSITION = POSITION,
DOSE = if(!is.null(irradiation_time_correction)){
dose_points + irradiation_time_correction[1]
}else{
dose_points
},
DOSE_ERROR = if(!is.null(irradiation_time_correction)){
dose_points * irradiation_time_correction[2]/irradiation_time_correction[1]
}else{
0
},
STEP = c("NATURAL", "REGENERATED"),
INTEGRAL = c(NATURAL, REGENERATED),
BACKGROUND = c(BACKGROUND, BACKGROUND),
NET_INTEGRAL = c(NATURAL - BACKGROUND, REGENERATED - BACKGROUND),
NATURAL_TL = NATURAL_TL,
REGENERATED_TL = REGENERATED_TL,
row.names = NULL
)
##0 dose points should not be biased by the correction ..
##Note: it does not mean that 0 s beneath the source has a dose of 0, however, in the certain
##case aliquot was never moved under the source
id_zero <- which(dose_points == 0)
temp_df$DOSE[id_zero] <- 0
temp_df$DOSE_ERROR[id_zero] <- 0
##calculate DE by using the irradiation time correction AND the cross talk correction
##(1) sample dose point values with irradiation time corrections (random)
if(!is.null(irradiation_time_correction)){
DOSE_MC <- rnorm(1000, mean = temp_df$DOSE[2], sd = temp_df$DOSE_ERROR[2])
}else{
DOSE_MC <- temp_df$DOSE[2]
}
##(2) random sampling from cross-irradiation
CT <- runif(1000, min = cross_talk_correction[2], max = cross_talk_correction[3])
##(3) signal ratio
INTEGRAL_RATIO <- temp_df$NET_INTEGRAL[1]/temp_df$NET_INTEGRAL[2]
##(4) calculate DE
temp_DE <- (DOSE_MC * INTEGRAL_RATIO)
##(5) substract cross-talk value from DE
temp_DE <- temp_DE - CT
##(5.1) calculate TL based DE
##calculate a dose based on TL
##Note: we use irradiation time correction and CT correction based on GSL measurements
if(calculate_TL_dose){
TL_Ratio <- NATURAL_TL/REGENERATED_TL
temp_TL_DE <- (DOSE_MC * TL_Ratio) - CT
TL_DE <- mean(temp_TL_DE)
TL_DE.ERROR <- sd(temp_TL_DE)
}else{
TL_DE <- NA
TL_DE.ERROR <- NA
}
##(6) create final data.frame
data <- data.frame(
DE = mean(temp_DE),
DE_ERROR = sd(temp_DE),
POSITION,
INTEGRAL_RATIO,
TRAVEL_DOSIMETER = NA,
CT_CORRECTION = cross_talk_correction[1],
CT_CORRECTION_Q2.5 = cross_talk_correction[2],
CT_CORRECTION_Q97.5 = cross_talk_correction[3],
TL_DE = TL_DE,
TL_DE.ERROR = TL_DE.ERROR,
row.names = NULL
)
##calculate test parameters
##TL_peak_shift
##check TL peak positions, if it differers more than the threshold, return a message
##can be done better, but should be enough here.
if (any("TL_peak_shift"%in%names(test_parameters)) && any(grepl("TL", names(object)))){
##calculate value
TP_TL_peak_shift.value <- abs((object[[2]][which.max(object[[2]][,2]),1] -
object[[4]][which.max(object[[4]][,2]),1]))
##compare
TP_TL_peak_shift.status <- TP_TL_peak_shift.value > test_parameters$TL_peak_shift
##return warning
if(TP_TL_peak_shift.status)
warning("TL peak shift detected for aliquot position ",POSITION, "! Check curves!", call. = FALSE)
##set data.frame
TP_TL_peak_shift <- data.frame(
CRITERIA = "TL_peak_shift",
THRESHOLD = test_parameters$TL_peak_shift,
VALUE = TP_TL_peak_shift.value,
STATUS = TP_TL_peak_shift.status,
stringsAsFactors = FALSE)
}else{
TP_TL_peak_shift <- data.frame(stringsAsFactors = FALSE)
}
##stimulation_power
if(any("stimulation_power"%in%names(test_parameters))){
##get curves ids holding the information on the stimulation power
temp_curves_OSL <- get_RLum(object_raw, recordType = "OSL", curveType = "measured")
temp_curves_OSL <- lapply(temp_curves_OSL, function(o){
if("stimulator"%in%names(o@info)){
if(grepl(o@info$stimulator, pattern = "LED", fixed = TRUE)){
return(o)
}
}
return(NULL)
})
##remove NULL
temp_curves_OSL <- temp_curves_OSL[!sapply(temp_curves_OSL, is.null)]
##check whether something is left
if(length(temp_curves_OSL) < 2){
TP_stimulation_power.value <- NA
TP_stimulation_power.status <- FALSE
}else{
##calculate sum of the power
TP_stimulation_power.value <- vapply(temp_curves_OSL, function(x){
sum(x@data[,2])
}, numeric(1))
##estimate a theoretical value based on the first value ... it does not
##matter which value is correct or not
TP_stimulation_power.value <- abs(1 -
sum(TP_stimulation_power.value)/(TP_stimulation_power.value[1] * length(TP_stimulation_power.value)))
TP_stimulation_power.status <- TP_stimulation_power.value > test_parameters$stimulation_power
if(TP_stimulation_power.status)
warning("Stimulation power was not stable for ALQ ",POSITION, "! Results are likely to be wrong!", call. = FALSE)
}
##remove object
rm(temp_curves_OSL)
##set data.frame
TP_stimulation_power <- data.frame(
CRITERIA = "stimulation_power",
THRESHOLD = test_parameters$stimulation_power,
VALUE = TP_stimulation_power.value,
STATUS = TP_stimulation_power.status,
stringsAsFactors = FALSE)
}else{
TP_stimulation_power <- data.frame(stringsAsFactors = FALSE)
}
##compile all test parameter df
df_test_parameters <- rbind(
TP_TL_peak_shift,
TP_stimulation_power)
# Terminal output -----------------------------------------------------------------------------
if(verbose){
cat(" [analyse_Al2O3_Measurement()] #",POSITION, " ", "DE: ",
round(data$DE, 2), " \u00B1 ", round(data$DE_ERROR,2), "\n", sep = "")
}
# Plotting ------------------------------------------------------------------------------------
##enable or disable plot ... we cannot put the condition higher, because we here
##calculate something we are going to need later
if (plot) {
##get plot settings
par.default <- par()$mfrow
on.exit(par(mfrow = par.default))
##settings
plot_settings <- list(
main = c(paste("ALQ POS:", POSITION, "| OSL"), paste("ALQ POS:", POSITION, "| TL")),
norm = TRUE,
mtext = ""
)
##modify on request
plot_settings <- modifyList(x = plot_settings, val = list(...),)
##plot curves
if(any(grepl("TL", names(object))))
par(mfrow = c(1,2))
plot_RLum(
object,
plot.single = TRUE,
combine = TRUE,
mtext = list(paste0("DE: ", round(data$DE,2), " \u00b1 ", round(data$DE_ERROR,2)), ""),
xlab = list("Simulation [s]", "Temperature [\u00B0C]"),
legend.text = list(list("#1 NAT", "#3 REG", "#5 BG"), list("#2 NAT", "#4 REG")),
legend.pos = list("topright", "topleft"),
main = as.list(plot_settings$main),
norm = plot_settings$norm
)
}
# Output --------------------------------------------------------------------------------------
UID <- create_UID()
output <- set_RLum(
class = "RLum.Results",
data = list(
data = cbind(data, UID),
data_table = cbind(temp_df, UID),
test_parameters = cbind(df_test_parameters, UID)
),
info = list(
call = sys.call()
)
)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.