---
title: "Obtaining a Dose Rate Calibration Curve"
author: "N. Frerebeau, B. Lebrun, S. Kreutzer"
date: "Last modified: `r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    fig_caption: yes
    toc: true
    number_sections: true
header-includes:
   - \usepackage{amsmath}
   - \usepackage{amssymb}
bibliography: bibliography.bib
vignette: >
  %\VignetteIndexEntry{Obtaining a Dose Rate Calibration Curve}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  message = FALSE,
  fig.width = 7,
  fig.height = 3.5,
  comment = "#>"
)
```

# Scope

This vignette walks you through the steps to obtain a calibration curve for the
gamma dose rate prediction. You will need the calibration curve to calculate 
the gamma dose rate from spectra measured in the field where the true gamma dose
is yet to be determined. In other words, you will need this curve to analyse
your field data. 

# Load library and import files

First things first, let's load the **gamma** package.

```{r packages}
library(gamma)
```

Now, we will use already measured spectra shipped with **gamma** for this example.
If you determine your calibration curve, you will need to use
your own data. All data used here were measured with a NaI probe. 

```{r import}
# Import CNF files for calibration
spc_dir <- system.file("extdata/AIX_NaI_1/calibration", package = "gamma")
spc <- read(spc_dir)
spc

# Import a CNF file of background measurement
bkg_dir <- system.file("extdata/AIX_NaI_1/background", package = "gamma")
bkg <- read(bkg_dir)
bkg
```

The object `spc` is a set of spectra measured in settings around Clermont-Ferrand 
[@miallier2009] and Bordeaux [@richter2010] with a known gamma dose rate, 
where `bgk` is a background curve measured in a lead castle. 

# Energy scale calibration

Before we can further work with the spectra, we have to perform an energy calibration, 
this is assigning values in terms of energy to the channel numbers. 

## Reference Spectra

First, we remove the baseline from the set of spectra for easier peak detection. 
, in the following subsections, we perform and apply the energy calibration 
to each spectrum, including the background spectrum. One subsection for each
spectrum with the corresponding R code.

```{r signal}
# Spectrum pre-processing
# Remove baseline for peak detection
bsl <- spc |>
  signal_slice(-1:-40) |>
  signal_stabilize(f = sqrt) |>
  signal_smooth(method = "savitzky", m = 21) |>
  signal_correct()
```

### BRIQUE
```{r calibrate-BRIQUE}
# Peak detection
pks <- peaks_find(bsl[["BRIQUE"]])
# Set energy values
set_energy(pks) <- c(238, NA, NA, NA, 1461, NA, NA, 2615)
# Adjust the energy scale
BRIQUE <- energy_calibrate(spc[["BRIQUE"]], pks)
```

```{r plot-BRIQUE, echo=FALSE}
plot(BRIQUE, pks) +
  ggplot2::theme_bw()
```

### C341
```{r calibrate-C341}
# Spectrum pre-processing and peak detection
pks <- peaks_find(bsl[["C341"]])
# Set energy values
set_energy(pks) <- c(238, NA, NA, NA, 1461, NA, 2615)
# Adjust the energy scale
C341 <- energy_calibrate(spc[["C341"]], pks)
```

```{r plot-C341, echo=FALSE}
plot(C341, pks) +
  ggplot2::theme_bw()
```

### C347
```{r calibrate-C347}
# Spectrum pre-processing and peak detection
pks <- peaks_find(bsl[["C347"]], span = 10)
# Set energy values
set_energy(pks) <- c(238, NA, NA, NA, NA, 1461, NA, 2615)
# Adjust the energy scale
C347 <- energy_calibrate(spc[["C347"]], pks)
```

```{r plot-C347, echo=FALSE}
plot(C347, pks) +
  ggplot2::theme_bw()
```

### GOU
```{r calibrate-GOU}
# Spectrum pre-processing and peak detection
pks <- peaks_find(bsl[["GOU"]])
# Set energy values
set_energy(pks) <- c(238, NA, NA, NA, 1461, NA, NA, 2615)
# Adjust the energy scale
GOU <- energy_calibrate(spc[["GOU"]], pks)
```

```{r plot-GOU, echo=FALSE}
plot(GOU, pks) +
  ggplot2::theme_bw()
```

### PEP
```{r calibrate-PEP}
# Spectrum pre-processing and peak detection
pks <- peaks_find(bsl[["PEP"]])
# Set energy values
set_energy(pks) <- c(238, NA, NA, NA, 1461, NA, NA, 2615)
# Adjust the energy scale
PEP <- energy_calibrate(spc[["PEP"]], pks)
```

```{r plot-PEP, echo=FALSE}
plot(PEP, pks) +
  ggplot2::theme_bw()
```

### Background Spectrum
```{r calibrate-bkg}
# Pb212, K40, Tl208
lines <- data.frame(
  channel = c(86, 496, 870),
  energy = c(238, 1461, 2615)
) 
bkg_scaled <- energy_calibrate(bkg, lines = lines)
```

```{r plot-bkg, echo=FALSE}
plot(bkg_scaled, xaxis = "energy", yaxis = "rate") +
  ggplot2::geom_vline(xintercept = c(238, 1461, 2615), linetype = 3) +
  ggplot2::theme_bw()
```

## Signal integration

Two methods can be used to calculate a gamma dose rate from the measured counts: the 'window' 
method and the 'threshold' method. The 'window' method compares counts under nuclide-specific gamma-ray peaks to  
reference spectra where the nuclide concentration is known. In this case, only 
a small part of the spectrum is used. The 'threshold' method uses the circumstance
that above a certain threshold, the number of counts for a given dose rate changes
only slightly with the emitters [e.g., U, Th, K -- cf. @guerin2011; @mercier2007; @lvborg1974]. 

The exact position of this threshold depends on the used probe and needs to be determined 
either by modelling or measurements of different sites with known but different radionuclides
concentrations. 

For this example (and in the package **gamma** in general), we will focus solely on the 'threshold' method, a crucial tool 
used to derive a calibration curve and derive gamma dose rates. This method's practical application in deriving a calibration curve 
is a key aspect of our exploration.  Our first step is to combine the energy-calibrated spectra into a unified set, 
laying the foundation for the subsequent calculations. 

```{r calibrate-spc}
spc_scaled <- list(BRIQUE, C341, C347, GOU, PEP)
spc_scaled <- methods::as(spc_scaled, "GammaSpectra")
spc_scaled
```

Now we can integrate the spectra using either the cumulative count spectra or
the energy corresponding to a certain energy [cf. @guerin2011 for details]. 

The following two subsections show how to derive the integrated spectrometer
response with each technique. In both cases, however, we assume that the
**threshold itself is already known**! The examples are for educational reasons
only and the here used function `signal_integrate()` is used in the back by the
function `dose_fit()` we will use later to derive the calibration curve. 

### Using the count threshold
```{r integrate-Ni}
# Integration range (in keV)
Ni_range <- c(200, 2800)

# Integrate background spectrum
Ni_bkg <- signal_integrate(
  object = bkg_scaled, 
  range = Ni_range, 
  energy = FALSE)

# Integrate reference spectra
Ni_spc <- signal_integrate(
  object = spc_scaled, 
  range = Ni_range, 
  background = Ni_bkg, 
  energy = FALSE, 
  simplify = TRUE)
```

### Using the energy threshold

```{r integrate-NiEi}
# Integration range (in keV)
NiEi_range <- c(200, 2800)

# Integrate background spectrum
NiEi_bkg <- signal_integrate(
  object = bkg_scaled, 
  range = NiEi_range, 
  energy = TRUE)

# Integrate reference spectra
NiEi_signal <- signal_integrate(
  object = spc_scaled, 
  range = NiEi_range, 
  background = NiEi_bkg, 
  energy = TRUE,
  simplify = TRUE)
```

## Get dose rate calibration curve

The dose rate calibration curve models detector specific counts against a know gamma-dose so
that such gamma dose can be derived from unknown spectra with the the same detector.

### Model calibration curve

To model the calibration curve we will need have measured sites with known radionuclide concentrations.
`'gamma'` ships a table with values from know sites in the dataset `clermont`. If you, as we have, measured
with your probe in that locations you do not have to manually add values from @miallier2009 or @richter2010.

```{r}
# Get reference dose rates
data("clermont")
doses <- clermont[, c("gamma_dose", "gamma_error")]
```

```{r, echo = FALSE}
knitr::kable(clermont)
```

In the next step we create a list with additional information about the used equipment
to be passed onto the calibration R object. You have to modify those values for your own 
detector.

```{r}
# Metadata
info <- list(
  laboratory = "CEREGE",
  instrument = "InSpector 1000",
  detector = "NaI",
  authors = "CEREGE Luminescence Team"
)
```

In the last step we construct the calibration curve and inspect the results. 

```{r}
# Build the calibration curve
AIX_NaI <- dose_fit(
  object = spc_scaled, 
  background = bkg_scaled, 
  doses = doses,
  range_Ni = Ni_range, 
  range_NiEi = NiEi_range,
  details = info
)
AIX_NaI 
```

If this information seems too opaque , we can further inspect the output numerically 
(`summary()`) and graphically using the package `plot()` methods.

```{r}
# show summary
summarise(AIX_NaI)
```

```{r calibration, fig.width=3.5, fig.show='hold'}
# plot calibration curves
plot(AIX_NaI, energy = FALSE) +
  ggplot2::theme_bw()
plot(AIX_NaI, energy = TRUE) +
  ggplot2::theme_bw()
```

Well, this looks all good and should make a good calibration curve. 
To use this calibration curve in the future, you can save the object using the 
base R function `save()`. 

```{r, eval=FALSE}
save(AIX_NaI, file = "<you_path>/<date>_NaI_DoseRate_Calibration.rda")
```

```{r save, eval=FALSE, echo=FALSE}
# DANGER ZONE
# AIX_NaI_1 <- AIX_NaI
# usethis::use_data(AIX_NaI_1, internal = FALSE, overwrite = FALSE)
```

### Check the model

Although the calibration curve appears reasonable, we can use a few additional 
R plots to further inspect the performed calibration. 

#### Count threshold (`Ni`)

```{r check-Ni, echo=FALSE, fig.width=3.5, fig.show='hold'}
Ni_residuals  <- get_residuals(AIX_NaI[["Ni"]])

# Residuals vs fitted values
ggplot2::ggplot(Ni_residuals, ggplot2::aes(x = fitted, y = residuals)) +
  ggplot2::geom_hline(yintercept = 0, linetype = 3) +
  ggplot2::geom_segment(ggplot2::aes(x = fitted, xend = fitted,
                                     y = 0, yend = residuals)) + 
  ggplot2::geom_point() +
  ggplot2::theme_bw() +
  ggplot2::labs(title = "Residuals vs fitted values",
                x = "Fitted values", y = "Residuals")

# Std. residuals vs fitted values
ggplot2::ggplot(Ni_residuals, ggplot2::aes(x = fitted, y = standardized)) +
  ggplot2::geom_hline(yintercept = 0, linetype = 3) +
  ggplot2::geom_hline(yintercept = c(-2, 2), linetype = 2) +
  ggplot2::geom_segment(ggplot2::aes(x = fitted, xend = fitted,
                                     y = 0, yend = standardized)) + 
  ggplot2::geom_point() +
  ggplot2::theme_bw() +
  ggplot2::labs(title = "Std. residuals vs fitted values",
                x = "Fitted values", y = "Standardized residuals")

# Normal QQ plot of standardized residuals
ggplot2::ggplot(Ni_residuals, ggplot2::aes(sample = standardized)) +
  ggplot2::geom_abline(slope = 1, intercept = 0) +
  ggplot2::geom_qq_line(colour = "red") +
  ggplot2::geom_qq() +
  ggplot2::theme_bw() +
  ggplot2::labs(title = "Normal QQ plot",
                x = "Theoretical quantiles",
                y = "Standardize residuals")

# Cook's distance
# ggplot2::ggplot(Ni_residuals, ggplot2::aes(x = name, y = cook)) +
#   ggplot2::geom_hline(yintercept = 0, linetype = 3) +
#   ggplot2::geom_hline(yintercept = 1, linetype = 2) +
#   ggplot2::geom_segment(ggplot2::aes(x = name, xend = name,
#                                      y = 0, yend = cook)) +
#   ggplot2::geom_point() +
#   ggplot2::theme_bw() +
#   ggplot2::labs(title = "Cook's distance",
#                 x = "Observation", y = "D")
```

#### Energy threshold (`Ni.Ei`)

```{r check-NiEi, echo=FALSE, fig.width=3.5, fig.show='hold'}
NiEi_residuals  <- get_residuals(AIX_NaI[["NiEi"]])

# Residuals vs fitted values
ggplot2::ggplot(NiEi_residuals, ggplot2::aes(x = fitted, y = residuals)) +
  ggplot2::geom_hline(yintercept = 0, linetype = 3) +
  ggplot2::geom_segment(ggplot2::aes(x = fitted, xend = fitted,
                                     y = 0, yend = residuals)) + 
  ggplot2::geom_point() +
  ggplot2::theme_bw() +
  ggplot2::labs(title = "Residuals vs fitted values",
                x = "Fitted values", y = "Residuals")

# Std. residuals vs fitted values
ggplot2::ggplot(NiEi_residuals, ggplot2::aes(x = fitted, y = standardized)) +
  ggplot2::geom_hline(yintercept = 0, linetype = 3) +
  ggplot2::geom_hline(yintercept = c(-2, 2), linetype = 2) +
  ggplot2::geom_segment(ggplot2::aes(x = fitted, xend = fitted,
                                     y = 0, yend = standardized)) + 
  ggplot2::geom_point() +
  ggplot2::theme_bw() +
  ggplot2::labs(title = "Std. residuals vs fitted values",
                x = "Fitted values", y = "Standardized residuals")

# Normal QQ plot of standardized residuals
ggplot2::ggplot(NiEi_residuals, ggplot2::aes(sample = standardized)) +
  ggplot2::geom_abline(slope = 1, intercept = 0) +
  ggplot2::geom_qq_line(colour = "red") +
  ggplot2::geom_qq() +
  ggplot2::theme_bw() +
  ggplot2::labs(title = "Normal QQ plot",
                x = "Theoretical quantiles",
                y = "Standardize residuals")

# Cook's distance
# ggplot2::ggplot(NiEi_residuals, ggplot2::aes(x = name, y = cook)) +
#   ggplot2::geom_hline(yintercept = 0, linetype = 3) +
#   ggplot2::geom_hline(yintercept = 1, linetype = 2) +
#   ggplot2::geom_segment(ggplot2::aes(x = name, xend = name,
#                                      y = 0, yend = cook)) +
#   ggplot2::geom_point() +
#   ggplot2::theme_bw() +
#   ggplot2::labs(title = "Cook's distance",
#                 x = "Observation", y = "D")
```

## Application example: opredict new dose rates

After calibrating our detector, it is time to derive gamma dose rates from 
sites with unknown radionuclide concentrations. As mentioned above, we will use 
data we ship with `'gamma'`. Our routine analysis consists of 
only three steps: 

### Import measured data

```{r}
# Import CNF files for dose rate prediction
test_dir <- system.file("extdata/AIX_NaI_1/test", package = "gamma")
test <- read(test_dir)
```

### Inspect the spectra 

```{r predict}
# Inspect spectra
plot(test, yaxis = "rate") +
  ggplot2::theme_bw()
```

### Energy calibrate the spectra

```{r}
# Pb212, K40, Tl208
pks <- data.frame(
  channel = c(86, 490, 870),
  energy = c(238, 1461, 2615)
) |> as("PeakPosition")

## energy calibrate
test <- energy_calibrate(test, pks)

## check the calibration for one curve
plot(test[[1]], pks) +
  ggplot2::theme_bw()

## show all energy calibrated spectra 
# Inspect spectra
plot(test, xaxis = "energy", yaxis = "rate") +
  ggplot2::theme_bw()
```

### Predict the dose rates

```{r}
rates <- dose_predict(AIX_NaI, test, sigma = 1.96)
rates
```

*Note: We assume here that the energy scale of each spectrum was adjusted first.*

# Outro: our R session {-}

Finally, for transparency, the R session setting used for rendering this vignette. 
```{r session-info, echo=FALSE}
sessionInfo()
```

# References {-}