## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
## Title:   Simulation Equivalent Dose Difference - Kreutzer & Mittelstraß, 2025
## Author:  Sebastian Kreutzer, Institute of Geography, Heidelberg University (DE)
## Contact: sebastian.kreutzer@uni-heidelberg.de
## Init:    Mon Aug 11 08:52:51 2025
## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Load libraries ----------------------------------------------------------
library(RLumModel)
library(Luminescence)

# Set parameters ----------------------------------------------------------
## set parameters for iteration
given_dose <- seq(1,1000,25)
## the signal resolution returned by RLumModel is cts/0.1 s
signal_reduction <- seq(0.0001, 0.0011, 0.0001)

## what are the scenarios we want to test
## we want to plot D0 against signal strength and then then see the
## difference
# Simulate signals --------------------------------------------------------
df <- lapply(given_dose, \(d) {
  cli::cli_progress_step(paste0("Runing dose: ", d))
  dose_points <- c(0,50,100,200,300,400,500,600,700,800,900,1000,0,50)
  sequence <- list(
    RegDose = dose_points,
    TestDose = 25,
    PH = 220,
    CH = 220,
    OSL_temp = 125,
    Irr_2recover = d,
    OSL_duration = 70)
  
  # own_parameters <- RLumModel::.set_pars("Bailey2001")
  # own_parameters$N[[3]] <- 1e+7
  # own_parameters$N[[9]] <- 1e+11
  
  ##model
  set.seed(1234)
  model.output <- model_LuminescenceSignals(
    sequence = sequence,
    simulate_sample_history = TRUE,
    model = "Bailey2001",
    #own_parameters = own_parameters,
    lab.dose_rate = 1,
    verbose = FALSE,
    plot = FALSE)
  
  lapply(signal_reduction, \(s) {
    ## lower count rates
    model.output@records <- lapply(model.output@records, \(x) {
      x@data[,2] <- x@data[,2] * s
      x
    })
    
    ## get max signal
    signal_max <- max(vapply(get_RLum(model.output, recordType = "^OSL$"), \(x) max(x@data[,2]), numeric(1)))
    
    ## get only curves of interest
    #object <- model.output
    set.seed(1234)
    results_normal <- get_RLum( model.output, recordType = c("^OSL$", "^TL$"), drop = FALSE) |>
      analyse_SAR.CWOSL(
        signal.integral.min = 1,
        signal.integral.max = 15,
        background.integral.min = 601,
        background.integral.max = 701,
        fit.method = "EXP",
        plot_onePage = TRUE,
        xlab = "Dose [a.u.]",
        verbose = FALSE,
        n.MC = 10,
        plot = FALSE,
        dose.points = dose_points
      )
    
    ## results corrected
    set.seed(1234)
    results_corrected <- get_RLum( model.output, recordType = c("^OSL$", "^TL$"), drop = FALSE) |>
      correct_PMTLinearity(PMT_pulse_pair_resolution = 18) |>
      analyse_SAR.CWOSL(
        signal.integral.min = 1,
        signal.integral.max = 15,
        background.integral.min = 601,
        background.integral.max = 701,
        fit.method = "EXP",
        xlab = "Dose [a.u.]",
        plot_onePage = TRUE,
        plot = FALSE,
        n.MC = 10,
        verbose = FALSE,
        dose.points = dose_points)
    
    ## combine in data.frame
    df <- data.frame(
      DOSE = d,
      SIGNAL = s,
      SIGNAL_MAX = signal_max,
      DE_NORMAL = results_normal$data$De,
      DE_NORMAL_X = results_normal$data$De.Error,
      DE_NORMAL_D0 = results_normal$data$D01,
      DE_CORRECTED = results_corrected$data$De,
      DE_CORRECTED_X = results_corrected$data$De.Error,
      DE_CORRECTED_D0 = results_corrected$data$D01
    )
  }) |> data.table::rbindlist()
}) |> data.table::rbindlist()


# Export ------------------------------------------------------------------
#write.csv(x = df, file = "data/Summary.csv")