R/analyse_Al2O3C_CrossTalk.R

In Luminescence: Comprehensive Luminescence Dating Data Analysis

#' @title Al2O3:C Reader Cross-Talk Analysis
#'
#' @description The function provides the analysis of cross-talk measurements
#' on a FI lexsyg SMART reader using Al2O3:C chips.
#'
#' @param object [Luminescence::RLum.Analysis-class] or [list] (**required**):
#' measurement input
#'
#' @param signal_integral [numeric] (*optional*):
#' signal integral, used for the signal and the background.
#' If nothing is provided, the full range is used.
#'
#' @param integral_input [character] (*with default*):
#' input type for `signal_integral`, one of `"channel"` (default) or
#' `"measurement"`. If set to `"measurement"`, the best matching channels
#' corresponding to the given time range (in seconds) are selected.
#'
#' @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 [Luminescence::get_RLum]. To deactivate the
#' automatic selection set the argument to `NULL`.
#'
#' @param irradiation_time_correction [numeric] or [Luminescence::RLum.Results-class] (*optional*):
#' information on the used irradiation time correction obtained by another
#' experiment.
#'
#' @param method_control [list] (*optional*):
#' optional parameters to control the calculation.
#' See details for further explanations.
#'
#' @param plot [logical] (*with default*):
#' enable/disable the plot output.
#'
#' @param ... further arguments and graphical parameters to control the plot
#' output. Supported are: `main`, `mtext`, and `pt.cex` (point size).
#'
#' @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 summed apparent dose table \cr
#'   `$data_full` \tab `data.frame` \tab full apparent dose table \cr
#'   `$fit` \tab `lm` \tab the linear model obtained from fitting \cr
#'   `$col.seq` \tab `numeric` \tab the used colour vector \cr
#'  }
#'
#' **slot:** **`@info`**
#'
#' The original function call
#'
#' ------------------------\cr
#' `[ PLOT OUTPUT ]`\cr
#' ------------------------\cr
#'
#' - An overview of the obtained apparent dose values
#'
#' @section Function version: 0.1.5
#'
#' @author Sebastian Kreutzer, F2.1 Geophysical Parametrisation/Regionalisation, LIAG - Institute for Applied Geophysics (Germany)
#'
#' @seealso [Luminescence::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. doi: 10.1515/geochr-2015-0086
#'
#' @keywords datagen
#'
#' @examples
#'
#' ##load data
#' data(ExampleData.Al2O3C, envir = environment())
#'
#' ##run analysis
#' analyse_Al2O3C_CrossTalk(data_CrossTalk)
#'
#' @export
analyse_Al2O3C_CrossTalk <- function(
  object,
  signal_integral = NULL,
  integral_input = c("channel", "measurement"),
  dose_points = c(0,4),
  recordType = "OSL (UVVIS)",
  irradiation_time_correction = NULL,
  method_control = NULL,
  plot = TRUE,
  ...
) {
  .set_function_name("analyse_Al2O3C_CrossTalk")
  on.exit(.unset_function_name(), add = TRUE)

  ## Integrity checks -------------------------------------------------------

  .validate_class(object, c("RLum.Analysis", "list"))
  .validate_not_empty(object, class(object)[1])
  if (is.list(object)) {
    lapply(object, .validate_class, class = "RLum.Analysis",
           name = "All elements of 'object'")
  } else {
    object <- list(object)
  }
  integral_input <- .validate_args(integral_input, c("channel", "measurement"))
  .validate_class(dose_points, c("numeric", "integer"))
  .validate_not_empty(dose_points)
  if (length(dose_points) != 1 && length(dose_points) %% 2 != 0) {
    .throw_error("'dose_points' should have length 1 or divisible by 2")
  }
  .validate_class(recordType, "character", null.ok = TRUE)
  .validate_class(irradiation_time_correction, c("numeric", "RLum.Results"),
                  null.ok = TRUE)
  .validate_class(method_control, "list", null.ok = TRUE)

  ##TODO ... do more, push harder
  ##Accept the entire sequence ... including TL and extract
  ##Add sufficient unit tests

  ## Preparation ------------------------------------------------------------

  ##select curves based on the recordType selection; if not NULL
  if(!is.null(recordType)){
    object <- suppressWarnings(get_RLum(object, recordType = recordType,
                                        drop = FALSE))
  }
  if (is.null(object) || all(lengths(object) == 0)) {
    .throw_error("'object' contains no records with recordType = '", recordType, "'")
  }

  #set method control
  method_control_settings <- list(
    fit.method = "SSE"
  )

  ##modify on request
  if(!is.null(method_control)){
    method_control_settings <- modifyList(x = method_control_settings, val = method_control)
  }

  ## signal integral
  x.range <- object[[1]]@records[[1]][, 1]
  if (integral_input == "measurement") {
    signal_integral <- .convert_to_channels(x.range, signal_integral,
                                            "time", null.ok = TRUE)
  }
  signal_integral <- .validate_integral(signal_integral,
                                        max = length(x.range), null.ok = TRUE)
  if (is.null(signal_integral)) {
    signal_integral <- seq_along(x.range)
  }

  ##check irradiation time correction
    if (inherits(irradiation_time_correction, "RLum.Results")) {
      .validate_originator(irradiation_time_correction, "analyse_Al2O3C_ITC")
      irradiation_time_correction <- get_RLum(irradiation_time_correction)

        ##insert 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]])
        }
    }

  # Calculation ---------------------------------------------------------------------------------
  ##we have two dose points, and one background curve, we do know only the 2nd dose

  ##create signal table list
  signal_table_list <- lapply(object, function(x) {
    ##calculate all the three signals needed
    NATURAL <- sum(x[[1]][signal_integral, 2])
    REGENERATED <- sum(x[[2]][signal_integral, 2])
    BACKGROUND <- sum(x[[3]][signal_integral, 2])

    temp_df <- data.frame(
      POSITION = get_RLum(x[[1]], info.object = "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, REGENERATED) - BACKGROUND,
      row.names = NULL
    )

    ##0 dose points should not be biased by the correction ..
    id_zero <- which(dose_points == 0)
    temp_df$DOSE[id_zero] <- 0
    temp_df$DOSE_ERROR[id_zero] <- 0

    return(temp_df)
  })

  APPARENT_DOSE <- data.table::rbindlist(lapply(signal_table_list, function(x) {
    ##run in MC run
    DOSE <- if (!is.null(irradiation_time_correction)) {
              rnorm(1000, mean = x$DOSE[2], sd = x$DOSE_ERROR[2])
            } else {
              x$DOSE[2]
            }

    ##calculation
    temp <- (DOSE * x$NET_INTEGRAL[1]) / x$NET_INTEGRAL[2]

    data.frame(
      POSITION = x$POSITION[1],
      AD = mean(temp),
      AD_ERROR = sd(temp))
  }))

  ##combine
  data_full <- as.data.frame(cbind(data.table::rbindlist(signal_table_list),
                                   APPARENT_DOSE[rep(1:.N, each = 2), 2:3]))

  # Plotting ------------------------------------------------------------------------------------
  par.default <- .par_defaults()
  on.exit(par(par.default), add = TRUE)

    ## set colours
    col_pal <- grDevices::hcl.colors(100, palette = "RdYlGn", rev = TRUE)

    ##settings
    plot_settings <- list(
      main = "Sample Carousel Crosstalk",
      mtext = "",
      pt.cex = 1
    )

      ##modify on request
      plot_settings <- modifyList(x = plot_settings, list(...))


    ##pre-calculations for graphical parameters
    n.positions <- length(unique(APPARENT_DOSE$POSITION))
    arc.step <- (2 * pi) / n.positions
    step <- 0

  ## calculate mean and standard deviation for similar positions
  AD <- POSITION <- NULL  # silence notes raised by R CMD check
  AD_matrix <- APPARENT_DOSE[, list(AD = mean(AD), AD_ERROR = sd(AD)),
                             by = POSITION]

    ##create colour ramp
    col.seq <- data.frame(
      POSITION = AD_matrix[order(AD_matrix[,2]),1],
      COLOUR = col_pal[seq(1,100, length.out = nrow(AD_matrix))],
      stringsAsFactors = FALSE)

    col.seq <- col.seq[["COLOUR"]][order(col.seq[["POSITION"]])]

    ##calculate model
    fit <- stats::lm(
      formula = y ~ poly(x, 2, raw=TRUE),
      data = data.frame(y = APPARENT_DOSE$AD[order(APPARENT_DOSE$POSITION)], x = sort(APPARENT_DOSE$POSITION)))

    ##enable or disable plot ... we cannot put the condition higher, because we here
    ##calculate something we are going to need later
    if (plot) {

      ##set layout matrix
      graphics::layout(mat = matrix(
        c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 3, 1, 1, 1, 1, 3),
        5,
        5,
        byrow = TRUE
      ))

      ##create empty plot
      par(
        mar = c(1, 1, 1, 1),
        omi = c(1, 1, 1, 1),
        oma = c(0.2, 0.2, 0.2, 0.2),
        cex = 1.1
      )
      shape::emptyplot(c(-1.15, 1.15), main = plot_settings$main, frame.plot = FALSE)

      ##add outher circle
      shape::plotcircle(r = 1.1, col = rgb(0.9, 0.9, 0.9, 1))

      ##add inner octagon
      shape::filledcircle(
        r1 = 0.6,
        mid = c(0, 0),
        lwd = 1,
        lcol = "black",
        col = "white"
      )

      ##add circles
      for (i in 1:n.positions) {
        shape::plotcircle(
          r = 0.05,
          mid = c(cos(step), sin(step)),
          cex = 6,
          pch = 20,
          col = col.seq[i]
        )
        text(x = cos(step) * 0.85,
             y = sin(step) * 0.85,
             labels = i)
        step <- step + arc.step
      }

      ##add center plot with position
      plot(NA, NA,
           xlim = range(AD_matrix[,1]),
           ylim = range(APPARENT_DOSE[,2]),
           frame.plot = FALSE,
           type = "l")

      ## add points
      points(x = APPARENT_DOSE,
             cex = plot_settings$pt.cex,
             pch = 20,
             col = rgb(0,0,0,0.3))

      ## add linear model
      lines(sort(APPARENT_DOSE$POSITION), stats::predict(fit), col = "red")

      ##add colour legend
      shape::emptyplot(c(-1.2, 1.2), frame.plot = FALSE)
      graphics::rect(
        xleft = rep(-0.6, 100),
        ybottom = seq(-1.2,1.1,length.out = 100),
        xright = rep(0, 100),
        ytop = seq(-1.1,1.2,length.out = 100),
        col = col_pal,
        lwd = 0,
        border = FALSE
      )

      ##add scale text
      text(
        x = -0.3,
        y = 1.2,
        label = "[s]",
        pos = 3,
        cex = 1.1
      )
      text(
        x = 0.4,
        y = 1,
        label = round(max(AD_matrix[, 2]),2),
        pos = 3,
        cex = 1.1
      )
      text(
        x = 0.4,
        y = -1.5,
        label = 0,
        pos = 3,
        cex = 1.1
      )
    }

  # Output --------------------------------------------------------------------------------------
  set_RLum(
    class = "RLum.Results",
    data = list(
      data = as.data.frame(AD_matrix),
      data_full = data_full,
      fit = fit,
      col.seq = col.seq
      ),
    info = list(
      call = sys.call()
    )
  )
}