Obtaining a Dose Rate Calibration Curve

In gamma: Dose Rate Estimation from in-Situ Gamma-Ray Spectrometry Measurements

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.

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.

# 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.

# 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

# 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)

plot(BRIQUE, pks) +
  ggplot2::theme_bw()

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)

plot(C341, pks) +
  ggplot2::theme_bw()

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)

plot(C347, pks) +
  ggplot2::theme_bw()

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)

plot(GOU, pks) +
  ggplot2::theme_bw()

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)

plot(PEP, pks) +
  ggplot2::theme_bw()

Background Spectrum

# Pb212, K40, Tl208
lines <- data.frame(
  channel = c(86, 496, 870),
  energy = c(238, 1461, 2615)
) 
bkg_scaled <- energy_calibrate(bkg, lines = lines)

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.

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

# 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

# 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.

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

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.

# 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.

# 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.

# show summary
summarise(AIX_NaI)

# 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().

save(AIX_NaI, file = "<you_path>/<date>_NaI_DoseRate_Calibration.rda")

# 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)

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)

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

# 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

# Inspect spectra
plot(test, yaxis = "rate") +
  ggplot2::theme_bw()

Energy calibrate the spectra

# 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

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.

sessionInfo()

References {-}

gamma documentation built on Sept. 24, 2024, 1:07 a.m.