---
title: "RLumModel - Getting started with RLumModel"
author: "Johannes Friedrich, University of Bayreuth"
date: '`r Sys.Date()`'
output: 
  rmarkdown::html_vignette:
    fig_caption: yes
    fig_height: 5
    fig_width: 5
    number_sections: yes
    toc: yes
bibliography: RLumModel.bib
csl: elsevier-harvard.csl
vignette: |
  %\VignetteIndexEntry{RLumModel - Getting started with RLumModel} 
  %\VignetteEncoding{UTF-8}{inputenc}
  %\VignetteEngine{knitr::rmarkdown} 
---

```{r global_options, include=FALSE}
knitr::opts_chunk$set(fig.pos = 'H', fig.align = 'center')
```

```{r, echo=FALSE, message = FALSE}
library(RLumModel)
library(knitr)
library(kableExtra)
```


# Introduction

This vignette shows a few examples for the **R**-package 'RLumModel'. The main function `model_LuminescenceSignals()` and their arguments will be explained.
All calculations were done with 'RLumModel' 0.2.2 and 'Luminescence' 0.7.5.


# Object structure of RLumModel

The output from the main function `model_LuminescenceSignals()` is of class `RLum.Analysis` [@Kreutzer_2012] and contains data of class `RLum.Data.Curve` in the slot 'records'. The advantage of this infrastructure is that the package 'Luminescence' offers a lot of methods to visualize and manipulate data. 

All simualted data are stored in the slot 'records': TL/OSL/RF curves as well as the concentrations of every energy level from every step.

The following code loads a data set provided by the 'RLumModel' package and shows how to seperate TL/OSL/RF data from concentrations and how to visualize them.

```{r, eval = FALSE}
data("ExampleData.ModelOutput", package = "RLumModel")

##show class
class(model.output)

##show structure
Luminescence::structure_RLum(model.output)

##seperate TL-curve from TL-concentrations
TL_curve <- Luminescence::get_RLum(model.output, recordType = "TL$")
TL_conc <- Luminescence::get_RLum(model.output, recordType = "(TL)", drop = FALSE)

##also possible: TL_curve <- get_RLum(model.output, record.id = 1)

##plot results
Luminescence::plot_RLum(TL_curve)
Luminescence::plot_RLum(TL_conc)

```

Some notes to the code example above:

* in 'TL_curve <- ...' appears "TL$". This is necessary to match the pattern "TL" without any sign after "TL", e.g. a bracket. The brackets are used (by default) for the concentrations.
* in 'TL_conc <- ...' the pattern "(TL)" will match all concentrations with "(TL)", see structure.
* `drop = FALSE` was used to keep the `RLum.Analysis` class for 'TL_conc'.
* To see a single plot of every energy-level, use the option `plot.single = TRUE` in `plot_RLum()`. For more details see the manual of 'Luminescence'.

```{r, eval = FALSE, fig.align = "center"}
##plot every energy-level by an extra plot
Luminescence::plot_RLum(TL_conc, plot.single = TRUE)
```

It is also possible to choose a `RLum.Data.Curve` by their 'record.id', which can be seen with:

```{r, eval = FALSE}
##see structure of model.output
Luminescence::structure_RLum(model.output)
```




# Selecting a quartz luminescence model

The first argument required for the function `model_LuminescenceSignals()` is the name of a 
quartz luminescence model to be used, respectively the used parameter set in this quartz luminescence model. 
All currently implemented quartz luminescence models were described in @Friedrich_2016.
The command to select a set of parameters from a specific model in `RLumModel` is a character string
with the name of the author and the year, e.g.

```{r}
model <- "Bailey2001"
```

The available models are "Bailey2001", "Bailey2002", "Bailey2004", "Pagonis2007", "Pagonis2008","Friedrich2017" and "customized" (@Bailey_2001, @Bailey_2002, @Bailey_2004, @Pagonis_2007a, @Pagonis_2008b, @Friedrich_2017). For customized or own parameter sets, see vignette [RLumModel - Using own parameter sets](RLumModel_-_Using_own_parameter_sets.html)

The corresponding parameter set will be loaded automatically with the function call. 

# Creating a sequence

The second argument in the `model_LuminescenceSignals()` function is the sequence to be simulated. 
There are three different ways of creating a sequence. 

For all sequences, temperature differences between sequence steps are automatically simulated with a heating or 
cooling step in between. Also, after irradiating the sample, it is automatically kept at irradiation temperature 
for further 5 s to allow the system to relax prior to the next step [@Bailey_2001].

## Risø SEQ files

The first one is to use the popular and freely available Risø \textit{Sequence Editor version 
4.36\footnote{http://www.nutech.dtu.dk/english/Products-and-Services/Dosimetry/Radiation-Measurement-Instruments/TL\_OSL\_reader/Software; 2016-04-11}} 
to build a personal sequence and to save it as a SEQ-file (*.seq). 
Files created by the Sequence Editor can be imported directly using the path of the SEQ-file. 
The package comes along with an example SEQ-file in the package folder in 'extdata'. 
Thus, a potential sequence is created with


```{r}
sequence <- system.file(
  "extdata",
  "example_SAR_cycle.SEQ",
  package = "RLumModel")
```

or wherever the SEQ-file is stored. While in the Sequence Editor irradiation is commonly defined in seconds, 
performing the simulation requires a dose transformation to Gray. 
Therefore, the function `model_LuminescenceSignals()` offers a special argument 
called `lab.dose_rate`, representing the dose rate of the irradiation unit in the laboratory. 
By default, this dose rate is 1 Gy/s, but can be modified, e.g.,

```{r}
lab.dose_rate <-  0.105
```


## Keywords

The second way of creating a sequence is by referring to a list with keywords 
and a certain order of code numbers or named values, which are shown
in Table 1.
With these keywords, it is possible to create quickly an R object of type list, 
which can be read by the `model_LuminescenceSignals()`
function.

```{r, echo = FALSE}
keywords <- data.frame(
  ARGUMENTS = c("TL","OSL","ILL","LM_OSL", "RF", "RF_heating", "IRR", "CH", "PH", "PAUSE"),
DESCRIPTION = c("Thermally stimulated luminescence", "Optically stimulated luminescence", "Illumination", "Linear modulated OSL", "Radiofluorescence", "RF during heating/cooling","Irradiation", "Cutheat", "Preheat", "Pause"),
SUB_ARGUMENTS = c("’temp_begin’ [°C], ’temp_end’ [°C], ’heating_rate’ [°C/s]",
"’temp’ [°C], ’duration’ [s], ’optical_power’ [%]",
"’temp’ [°C], ’duration’ [s], ’optical_power’ [%]",
"’temp’ [°C], ’duration’ [s], optional: ’start_power’ [%], ’end_power’ [%]",
"’temp’ [°C], ’dose’ [Gy], ’dose_rate’ [Gy/s]",
"’temp_begin’ [°C], ’temp_end’ [°C], 'heating rate' [°C/s], 'dose_rate' [Gy/s]",
"’temp’ [°C], ’dose’ [Gy], ’dose_rate’ [Gy/s]",
"’temp’ [°C], optional: ’duration’ [s], ’heating_rate’ [°C/s]",
"’temp’ [°C], ’duration’ [s], optional: ’heating_rate’ [°C/s]",
"’temp’ [°C], ’duration [s]’"))

kable(keywords,
      format = "html",
      col.names = c("ARGUMENTS", "DESCRIPTION", "SUB-ARGUMENTS"),
      caption = "Keywords in RLumModel for creating sequences") %>%
    kable_styling(bootstrap_options = c("striped", "hover"), full_width = F)

```

Some examples to this kind of sequence creating:

```{r}
sequence <- list(
 IRR = c(temp = 20, dose = 10, dose_rate = 1),
 TL = c(temp_begin = 20, temp_end = 400 , heating_rate = 5))

```

This sequences describes an irradiation simulation at 20 °C with a dose of
10 Gy and a dose rate of 1 Gy/s, which is followed by a TL simulation from
20 °C to 400 °C with a heating rate of 5 °C/s. Note that it is important
that for each sequence keyword like 'IRR' or 'TL' either the vector has to
be named or the correct order of arguments is used, see 'sub-arguments' in
Table 1. Thus the above mentioned code is equivalent to the following one:

```{r}
sequence <- list(
  IRR = c(20, 10, 1),
  TL = c(20, 400, 5))
```

## Creating a SAR/DRT sequence

However, to create a SAR or dose-recovery-test (DRT) sequence with
the Risø Sequence Editor or with keywords is time-consuming, because it
contains a lot of individual sequence steps (preheat, optical stimulation, irradiation, ...). 
Therefore, a third way was implemented in 'RLumModel' to
create a (SAR) sequence after Murray and Wintle (2000) with the (required)
keywords RegDose, TestDose, PH, CH and OSL temp. In addition to these keywords, 
the user is able to set more detailed parameters for the SAR sequence,
see Table 2:

```{r, echo = FALSE}
SAR_sequence_table <- data.frame(
  ABBREVIATION = c("RegDose", "TestDose", "PH", "CH", "OSL_temp", "OSL_duration", "Irr_temp", "PH_duration", "dose_rate", "optical_power", "Irr_2recover"),
  DESCRIPTION = c("Dose points of the regenerative cycles [Gy]", 
                  "Test dose for the SAR cycles [Gy]",
                  "Temperature of the preheat [°C]",
                  "Temperature of the cutheat [°C]",
                  "Temperature of OSL read out [°C]",
                  "Duration of OSL read out [s]",
                  "Temperature of irradiation [°C]",
                  "Duration of the preheat [s]",
                  "Dose rate of the laboratory irradiation source [Gy/s]",
                  "Percentage of the full illumination power [%]",
                  "Dose to be recovered in a dose-recovery-test [Gy]"),
  EXAMPLE_ARGUMENTS = c("c(0, 80, 140, 260, 320, 0, 80)",
                        "50",
                        "240",
                        "200",
                        "125",
                        "default: 40",
                        "default: 20",
                        "default: 10",
                        "default: 1",
                        "default: 90",
                        "20"))

kable(SAR_sequence_table,
      format = "html",
      col.names = c("ABBREVIATION", "DESCRIPTION", "EXAMPLE ARGUMENTS"),
      caption = "Keywords in RLumModel for creating SAR sequences") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = F)
```

So a possible DRT sequence could be the next code example:

```{r}
sequence <- list(
  RegDose = c(0,10,20,50,90,0,10),
  TestDose = 2,
  PH = 220,
  CH = 220,
  OSL_temp = 125,
  Irr_2recover = 20)
```

This sequence describes a DRT, where a dose of 20 Gy will be recovered with this test. The regenerative doses are definded as 0 (natural), 10 Gy, 20 Gy, 50 Gy, 90 Gy and for recuperation and recycling ratio 0 Gy and 10 Gy, respectively.
The test dose is definded as 2 Gy. Preheat and Cutheat are at 220 °C and all OSL measurements are simualted at 125 °C. There are more options to set, see Table 2.

The RLumModel function `model_LuminescenceSignals()` is able to interprete this (sequence-) list as a DRT sequence.

# Working examples

## Simulate a TL measurement

First of all, a sequence is needed, which produces a TL signal after the sample has received a dose:

```{r}
sequence <- list (
IRR = c (20 , 10 , 1) ,
TL = c (20 , 400 , 5))
```

Here, at a temperature of 20 °C a dose of 10 Gy was applied with a dose rate of 1 Gy/s followed by a TL measurement from 20 °C to 400 °C with a heating rate of 5 °C/s. Running this sequence with the `model_LuminescenceSignals()` function produces a model output:

```{r, fig.cap = "TL curve with parameter set 'Bailey2001' after 10 Gy laboratory dose"}
model.output <- model_LuminescenceSignals(
  model = "Bailey2001",
  sequence = sequence,
  verbose = FALSE)
```

This results in a TL curve like the one published in (@Bailey_2001, Fig. 1), see figure above. In a further step, it is easy to produce known TL phenomena like the shift of the TL peak with varying heating rate. For this purpose, a loop over a TL simulation changes the heating rate in each run.

```{r, fig.cap = "TL signal with different heating rates"}
##set heating rate
heating.rate <- seq(from = 2, to = 10, by = 2)

##model signals
##"verbose = FALSE" for no terminal output
## "TL$" for exact matching TL and not (TL)
model.output <- lapply(heating.rate, function(x){
  sequence <- list(
   IRR = c(20, 10, 1),
   TL = c(20, 400, x))


  TL_data <- model_LuminescenceSignals(
      sequence = sequence,
      model = "Bailey2001",
      plot = FALSE,
      verbose = FALSE)
  
  return(Luminescence::get_RLum(TL_data, recordType = "TL$", drop = FALSE))

})

##merge output
model.output.merged <- merge_RLum(model.output)

##plot results
plot_RLum(
 object = model.output.merged,
 xlab = "Temperature [\u00B0C]",
 ylab = "TL signal [a.u.]",
 main = "TL signal with different heating rates",
 legend.text = paste(heating.rate, "°C/s"),
 combine = TRUE)
```

Some notes to the code above:

* the return of the lapply function is a `RLum.Analysis` object, because of drop = FALSE
* `recordType = TL$` is necessary to match the character TL exact and not the concentrations
* `merge_RLum` is necessary to merge all the single `RLum.Analysis` objects in the list 'model.output' together to one `RLum.Analysis` object
* to see the results with another parameter set, only `model = "..."` has to be changed (see Sec. 2)

## Simulating thermal activation characteristics (TACs)

Another frequently simulated phenomenon is the sensitisation of quartz TL by $\beta$- or $\gamma$-irradiation 
followed by activation at high temperatures [@Zimmerman_1971, @Pagonis_2003, @Pagonis_2008a, @Adamiec_2004], termed as thermal activation characteristics (TACs). 
For a simulation sequence, the reader is referred to @Pagonis_2008a, Tab. 1. To simulate this phenomenon
with the `model_LuminescenceSignals()` function, a loop causing a stepwise increase 
of the activation temperature is needed. The resulting intensity of the 110 °C TL peak can be plotted 
against the activation temperature, which shows TAC for the model parameters 
of "Pagonis2007".

```{r TACs, fig.align="center"}
##set temperature
act.temp <- seq(from = 80, to = 600, by = 20)

##loop over temperature
model.output <- vapply(X = act.temp, FUN = function(x) {
  
  ##set sequence, note: sequence includes sample history
  sequence <- list(
    IRR = c(20, 1, 1e-11),
    IRR = c(20, 10, 1),
    PH = c(x, 1),
    IRR = c(20, 0.1, 1),
    TL = c(20, 150, 5)
    )

  ##run simulation
  temp <- model_LuminescenceSignals(
    sequence = sequence,
    model = "Pagonis2007",
    simulate_sample_history = TRUE,
    plot = FALSE,
    verbose = FALSE
    )
  
  ## "TL$" for exact matching TL and not (TL)
  TL_curve <- Luminescence::get_RLum(temp, recordType = "TL$")

  ##return max value in TL curve
  return(max(get_RLum(TL_curve)[,2]))

}, FUN.VALUE = 1)
```

```{r, echo=FALSE, fig.cap = "TAC with parameter set of 'Pagonis2007'"}
##plot resutls
plot(
  act.temp[-(1:3)],
  model.output[-(1:3)],
  type = "b",
  xlab = "Temperature [\u00B0C]",
  ylab = "TL [a.u.]"
  )
```

## Simulating dependency of the OSL signal on the illumination power density

The function `model_LuminescenceSignals()` is also capable of simulating OSL phenomena.
The next figure shows the dependency of the OSL signal on the power density of illumination for the model "Bailey2004". 

```{r, fig.cap = "OSL measurement with different optical power densities with the parameter set of 'Bailey2004'"}
##set optical power [%]
optical_power <- c(0,20,40,60,80,100)

##loop over power 
model.output <- lapply(optical_power, function(x){

  ##set sequence
  sequence <- list(
    IRR = c(20, 50, 1),
    PH = c(220, 10, 5),
    OSL = c(125, 50, x))

  
  data <-  model_LuminescenceSignals(
    sequence = sequence,
    model = "Bailey2004",
    plot = FALSE,
    verbose = FALSE)
  
  ##"OSL$" for exact matching OSL and not (OSL)
  return(Luminescence::get_RLum(data, recordType = "OSL$", drop = FALSE))
  
})

##merge output
model.output.merged <- Luminescence::merge_RLum(model.output)

##plot results
Luminescence::plot_RLum(
  object = model.output.merged,
  xlab = "Illumination time [s]",
  ylab = "OSL signal [a.u.]",
  legend.text = paste("Optical power ", 20 * optical_power / 100," mW/cm^2"),
  combine = TRUE 
  )
```

## Simulating and analysing SAR measurements

For simulating a DRT, it is necessary to define a sequence with the keyword `Irr_2recover`, 
as  mentioned in Section 4.3. 

It should be mentioned that a simulation of a combinded PHPT and DRT may be very time-consuming, 
because for every preheat temperature a complete SAR cycle has to be run. 
A typical DRT sequence featuring various PH temperatures in 'RLumModel' is listed below. 
Note that in such a DRT simulation a loop over different preheat temperatures has to be written, 
utilising characteristic parameters from the literature. The test dose is set to 10\% and the regeneration 
dose points to 40\%, 70\%, 130\%, 160\%, 0\% and 40\% of the recovery dose.

The data created by 'RLumModel' can be directly passed to the functions 
`Luminescence::analyse_SAR.CWOSL()` and `Luminescence::plot_DRTResults()` 
for routine analyses and plotting.

```{r DRT, fig.cap = "Dose recovery test (DRT) with the parameter set of 'Pagonis2008'"}
##set PH temperatures
PH_temp <- seq(from = 160, to = 300, by = 20)

##set regeneration doses
RegDose = c(0, 80, 140, 260, 320, 0, 80)

##loop over PH temperatures
DRT.output <- lapply(PH_temp, function(x){

  sequence <- list(
       RegDose = RegDose,
       TestDose = 20,
       PH = x,
       CH = x,
       OSL_temp = 125,
       Irr_2recover = 200)

  model.output <- model_LuminescenceSignals(
       sequence = sequence,
       model = "Pagonis2008",
       plot = FALSE,
       verbose = FALSE)

  results <- Luminescence::analyse_SAR.CWOSL(
       object = model.output,
       signal.integral.min = 1,
       signal.integral.max = 7,
       background.integral.min = 301,
       background.integral.max = 401,
       fit.method = "EXP",
       dose.points = RegDose,
       plot = FALSE)
															
  temp <- get_RLum(results)
  out <- data.frame(
    De = temp$De, 
    De.error = temp$De.Error)
    
  return(out)
})

##output as data.frame for plot_DRTResults
DRT.result <- as.data.frame(do.call(rbind, DRT.output))

##plot DRT.results
Luminescence::plot_DRTResults(
     DRT.result, 
     preheat = PH_temp,
     given.dose = 200)
```

In the code above, `plot = FALSE` was chosen, because a single OSL plot is not necessary to 
analyse a SAR sequence. To calculate a D<sub>e</sub> from the 
produced `RLum.Analysis` object 'model.output', the 
function `Luminescence::analyse_SAR.CWOSL()` is suitable.
After specifying a number of evaluation parameters (signal and background integration interval, 
dose points and fit function for the dose response curve) and the analysis process, the reduced data are stored 
in an `RLum.Results` object, here termed 'results' . 
A background integration interval from 301 to 401 translates to the signal from 30 s to 40 s, 
because a channel has the default width of 0.1 s. Accordingly, the signal integral ranges from 0.1 s to 0.7 s. 

```{r SAR, fig.cap = "SAR protocol with the parameter set of 'Pagonis2008'"}
##set RegDose
RegDose = c(0, 80, 140, 260, 320, 0, 80)

##set sequence
sequence <- list(
  RegDose = RegDose,
  TestDose = 20,
  PH = 220,
  CH = 220,
  OSL_temp = 125
  )

##model 
model.output <- model_LuminescenceSignals(
  sequence = sequence,
  model = "Pagonis2008",
  plot = FALSE,
  verbose = FALSE
  )
		
##analyse SAR sequence and plot only the resulting growth curve	
results <-Luminescence::analyse_SAR.CWOSL(
  model.output,
  signal.integral.min = 1,
  signal.integral.max = 7,
  background.integral.min = 301,
  background.integral.max = 401,
  fit.method = "EXP",
  dose.points = RegDose, 
  verbose = FALSE, 
  plot.single = c(6)
  )
```


# References 