Calibration Curve

N. Frerebeau, B. Lebrun

2024-04-08

library(gamma)

This vignette shows how to build a calibration curve for gamma dose rate prediction.

Import Files

# Import CNF files for calibration
spc_dir <- system.file("extdata/AIX_NaI_1/calibration", package = "gamma")
(spc <- read(spc_dir))
#> A collection of 5 gamma spectra: BRIQUE, C341, C347, GOU, PEP

# Import a CNF file of background measurement
bkg_dir <- system.file("extdata/AIX_NaI_1/background", package = "gamma")
(bkg <- read(bkg_dir))
#> Gamma spectrum:
#> *  name: PB
#> *  date: 2019-03-27 12:06:02
#> *  live_time: 7707.42
#> *  real_time: 7714.93
#> *  channels: 1024
#> *  energy_min: -7
#> *  energy_max: 3124.91

Energy Scale Calibration

Reference Spectra

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

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)

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)

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)

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)

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)

Signal Integration

spc_scaled <- list(BRIQUE, C341, C347, GOU, PEP)
spc_scaled <- methods::as(spc_scaled, "GammaSpectra")

Count Threshold

# Integration range (in keV)
Ni_range <- c(200, 2800)

# Integrate background spectrum
(Ni_bkg <- signal_integrate(bkg_scaled, range = Ni_range, energy = FALSE))
#>      value      error 
#> 1.40046864 0.01906326

# Integrate reference spectra
(Ni_spc <- signal_integrate(spc_scaled, range = Ni_range, background = Ni_bkg, 
                            energy = FALSE, simplify = TRUE))
#>           value     error
#> BRIQUE 209.2953 0.2135497
#> C341    87.7405 0.2249214
#> C347   150.6917 0.2929845
#> GOU    170.0013 0.3105710
#> PEP    264.7712 0.3969934

Energy Threshold

# Integration range (in keV)
NiEi_range <- c(200, 2800)

# Integrate background spectrum
(NiEi_bkg <- signal_integrate(bkg_scaled, range = NiEi_range, energy = TRUE))
#>        value        error 
#> 1108.0565661    0.5362181

# Integrate reference spectra
(NiEi_signal <- signal_integrate(spc_scaled, range = NiEi_range, 
                                 background = NiEi_bkg, energy = TRUE,
                                 simplify = TRUE))
#>            value    error
#> BRIQUE 110923.34 4.933826
#> C341    46657.02 5.215421
#> C347    81717.33 6.843667
#> GOU     89644.92 7.152993
#> PEP    139711.03 9.136525

Calibration Curve

Summary

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

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

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

# Summary
summarise(AIX_NaI)
#> $Ni
#> $Ni$residuals
#> [1]  33.402986   7.968894   6.460507 -18.670747  78.089036
#> 
#> $Ni$coefficients
#>            Estimate Std. Error
#> Intercept 40.013121 41.8531813
#> Slope      9.140417  0.2721024
#> 
#> $Ni$MSWD
#> [1] 0.8938093
#> 
#> $Ni$df
#> [1] 3
#> 
#> $Ni$p_value
#> [1] 0.4433926
#> 
#> 
#> $NiEi
#> $NiEi$residuals
#> [1]  33.906630  12.519006 -21.932843  -8.121213  86.160947
#> 
#> $NiEi$coefficients
#>              Estimate   Std. Error
#> Intercept 27.88537142 4.202826e+01
#> Slope      0.01735135 5.142132e-04
#> 
#> $NiEi$MSWD
#> [1] 0.9375068
#> 
#> $NiEi$df
#> [1] 3
#> 
#> $NiEi$p_value
#> [1] 0.421443

# Plot curve
plot(AIX_NaI, model = "Ni") +
  ggplot2::theme_bw()
plot(AIX_NaI, model = "NiEi") +
  ggplot2::theme_bw()

Check Model

Count Threshold (Ni)

Energy Threshold (Ni.Ei)

Predict New Dose Rates

# Import CNF files for dose rate prediction
test_dir <- system.file("extdata/AIX_NaI_1/test", package = "gamma")
(test <- read(test_dir))
#> A collection of 5 gamma spectra: NAR19-P2-1, NAR19-P3-1, NAR19-P4-1, NAR19-P5-1, NAR19-P6-1

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


# Dose rate prediction
# (assuming that the energy scale of each spectrum was adjusted first)
(rates <- dose_predict(AIX_NaI, test, sigma = 2))
#>        names   dose_Ni error_Ni dose_NiEi error_NiEi
#> 1 NAR19-P2-1  925.2131 1388.917  898.0336   1348.102
#> 2 NAR19-P3-1 1116.4969 1676.068 1089.8538   1636.056
#> 3 NAR19-P4-1  894.0082 1342.074  869.4770   1305.233
#> 4 NAR19-P5-1 1391.7762 2089.312 1369.9492   2056.527
#> 5 NAR19-P6-1 1171.7514 1759.018 1152.1000   1729.499

R session

#> R version 4.3.3 (2024-02-29)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 22.04.4 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#> 
#> locale:
#>  [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=C              
#>  [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   
#>  [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Europe/Paris
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] gamma_1.0.5
#> 
#> loaded via a namespace (and not attached):
#>  [1] vctrs_0.6.5       cli_3.6.2         knitr_1.45        rlang_1.1.3      
#>  [5] xfun_0.43         highr_0.10        rxylib_0.2.12     generics_0.1.3   
#>  [9] jsonlite_1.8.8    labeling_0.4.3    glue_1.7.0        colorspace_2.1-0 
#> [13] htmltools_0.5.8   sass_0.4.9        fansi_1.0.6       scales_1.3.0     
#> [17] rmarkdown_2.26    grid_4.3.3        evaluate_0.23     munsell_0.5.0    
#> [21] jquerylib_0.1.4   tibble_3.2.1      fastmap_1.1.1     yaml_2.3.8       
#> [25] lifecycle_1.0.4   compiler_4.3.3    dplyr_1.1.4       Rcpp_1.0.12      
#> [29] pkgconfig_2.0.3   rstudioapi_0.16.0 farver_2.1.1      IsoplotR_6.1     
#> [33] digest_0.6.35     R6_2.5.1          tidyselect_1.2.1  utf8_1.2.4       
#> [37] pillar_1.9.0      magrittr_2.0.3    bslib_0.7.0       withr_3.0.0      
#> [41] tools_4.3.3       gtable_0.3.4      ggplot2_3.5.0     cachem_1.0.8