# Overview

gamma is intended to process in-situ gamma-ray spectrometry measurements for luminescence dating. This package allows to import, inspect and (automatically) correct the energy scale of the spectrum. It provides methods for estimating the gamma dose rate by the use of a calibration curve. This package only supports Canberra CNF and TKA files.

Typical workflow:

1. Import spectrum data with read(),
2. Inspect spectrum with plot(),
3. Calibrate the energy scale with calibrate(),
4. Estimate the gamma dose rate with predict().

The estimation of the gamma dose rate requires a calibration curve. This package contains several built-in curves, but as these are instrument-specific you may need to build your own (see help(fit)). These built-in curves are in use in several labs and can be used to reproduce published results. Check out the following vignettes for an overview of the fitting process.

# IRAMAT-CRP2A LaBr (BDX1)
utils::vignette("CRP2A#1", package = "gamma")

# CEREGE NaI (AIX1)
utils::vignette("CEREGE#1", package = "gamma")

Note that gamma uses the International System of Units: energy values are assumed to be in keV and dose rates in µGy/y.

# Import file

Both Canberra CNF and TKA files can be imported. Several channels can be skipped during import to retain only part of the spectrum.

If skip is TRUE, an attempt is made to define the number of channels to skip at the beginning of the spectrum. This skips all channels before the highest count maximum. This is intended to deal with the artefact produced by the rapid growth of random background noise towards low energies.

# Automatically skip the first chanels
# Import a CNF file
cnf_test <- system.file("extdata/test_CNF.cnf", package = "gamma")
(spc_cnf <- read(cnf_test, skip = TRUE))
#> Gamma spectrum:
#>   Reference: test_CNF
#>   Instrument: InSpector 1000
#>   Date: 2019-02-07 11:48:18
#>   Number of chanels: 989
#>   Energy range (keV): 97.94 - 3124.91
#>   Dose rate: not known

# Import a TKA file
tka_test <- system.file("extdata/test_TKA.tka", package = "gamma")
(spc_tka <- read(tka_test, skip = TRUE))
#> Gamma spectrum:
#>   Reference: test_TKA
#>   Instrument: unknown
#>   Date: 2019-06-06 20:07:46
#>   Number of chanels: 989
#>   Energy range (keV): not calibrated
#>   Dose rate: not known

# Import all files in a given directory
files_test <- system.file("extdata", package = "gamma")
(spc <- read(files_test, skip = TRUE))
#> A collection of 2 gamma spectra: test_CNF, test_TKA

# Inspect spectrum

# Plot TKA spectrum
plot(spc_cnf) +
ggplot2::theme_bw()
#> Registered S3 methods overwritten by 'ggplot2':
#>   method         from
#>   [.quosures     rlang
#>   c.quosures     rlang
#>   print.quosures rlang

# Plot spectra
plot(spc, facet = TRUE) +
ggplot2::theme_bw()

# Calibrate the energy scale

The energy calibration of a spectrum is the most tricky part. To do this, you must specify the position of at least two observed peaks and the corresponding energy value (in keV).

The package allows you to provide the channel-energy pairs you want to use. However, the spectrum can be noisy so it is difficult to properly determine the peak channel. In this case, a better approach may be to fit each peak with a Gaussian curve and then use the parameters estimated in this way.

Regardless of the approach you choose, it is strongly recommended to check the result before proceeding.

## Peak parameters estimation

The automatic estimation of peak parameters is performed in two steps that can influence the final result. First, the baseline of the spectrum is estimated and removed. Then, a non-linear model is fitted for each peak. The model parameters can finally be used for the spectrum energy scale calibration.

These two steps (baseline and peak fitting) are performed in one go. The following sections will help you to fine-tune the individual parameters.

### Estimate and remove baseline

The baseline estimation is performed using the SNIP algorithm (Ryan et al. 1988; Morháč et al. 1997; Morháč and Matoušek 2008). You can apply the LLS operator to your data, use a decreasing clipping window or change the number of iterations (see references).

# Estimate the baseline of a single file
# Default parameters:
# * The LLS operator is not used,
# * An increasing clipping window is used,
# * 100 iterations are performed.
baseline <- estimateBaseline(spc_tka)

# Plot spectrum + baseline
plot(spc_tka, baseline) +
ggplot2::theme_bw()

# Substract the previously estimated baseline
spc_clean1 <- spc_tka - baseline
# Or, remove the baseline in one go
spc_clean2 <- removeBaseline(spc_tka)

identical(spc_clean1, spc_clean2)
#> [1] TRUE

# Plot the corrected spectrum
plot(spc_clean1) +
ggplot2::theme_bw()

# You can change the default parameters (see help for details)
# help(SNIP)
spc_clean3 <- removeBaseline(spc_tka, LLS = TRUE)

# The default settings don't seem too bad
plot(spc_clean1, spc_clean3) +
ggplot2::theme_bw()

# You can remove the baseline of multiple spectra in one go
# Note that the same parameters will be used for all spectra
spc_clean <- removeBaseline(spc)

### Automatic peak detection

Once you have fixed the parameters for the baseline estimation, you can try to automatically find peaks in the spectrum. A local maximum has to be the highest one in the given window and has to be higher than $$k$$ times the noise to be recognized as peak.

Again, peak finding is performed on the baseline corrected spectrum, so you may want to pass extra arguments for the baseline estimation.

# Detect peaks in a single file
peak_index <- findPeaks(spc_tka)

# Plot spectrum
plot(peak_index) +
ggplot2::theme_bw()

Automatic peak detection
Chanel Energy [keV] Count Count rate [1/s]
peak #1 86 NA 2376 0.702
peak #2 210 NA 840 0.248
peak #3 316 NA 297 0.088
peak #4 496 NA 910 0.269
peak #5 876 NA 110 0.032

The result of this automatic detection can be used as starting positions for peak fitting.

### Peak fitting

If you know the approximate position of some peaks, you can refine it by fitting each peak with a Gaussian. The parameters of these non-linear models are determined using the port algorithm (see help(nls) from the {stats} package).

You only need to provide the starting peak positions (chanel). If problems arise, it may be useful to constrain the parameters (see help(fit)) or improve the baseline estimate (peak fitting is performed on the baseline corrected spectrum).

# Fit peaks at a given chanel
peak_fit <- fitPeaks(spc_tka, peaks = c(86, 496, 876))

## Plot results
plot(peak_fit) +
ggplot2::labs(title = "Fitted peaks from fixed starting positions") +
ggplot2::theme_bw()

Peak positions
Mean (chanel) Std. dev. (chanel) Height (count)
peak #1 85.806 2.686 2447.413
peak #2 494.862 9.744 877.721
peak #3 877.186 14.328 98.544

## Energy scale calibration

### The careful way

If you know the correspondence between several channels and the energy scale, you can use these pairs of values to calibrate the spectrum. A second order polynomial model is fitted on these energy vs chanel values, then the chanel values are used to predict the new energy scale.

# Create a list of chanel-enegy pairs
# Names are optional, but values must always be set in this order:
# channel then energy
calib_lines <- list(
Pb212 = c(chanel = 84, energy = 238),
K40 = c(chanel = 492, energy = 1461),
Tl208 = c(chanel = 865, energy = 2615)
)

# Calibrate the spectrum using these fixed lines
tka_scaled1 <- calibrate(spc_tka, lines = calib_lines)

Alternatively, you can use the results of the peak fitting and pass the expected energy values.

# Calibrate the spectrum using the peak models
# Lines should be in keV
tka_scaled2 <- calibrate(peak_fit, lines = c(238, 1461, 2615))

# Plot the spectrum (note the top scale)
plot(tka_scaled2) +
ggplot2::labs(title = "Calibrated spectra") +
ggplot2::theme_bw()

### The fast way

If you want to calibrate several spectra at once, you can try the following approach.

# The 'GammaSpectra' class inherits from 'list'
# You can apply a function over it

# First, fit peaks at a given chanel
# The same starting peak positions are used for all spectra
multi_fit <- lapply(X = spc, FUN = fitPeaks, peaks = c(86, 496, 876))

# Then, calibrate the spectra using the peak models
multi_scaled <- lapply(X = multi_fit, FUN = calibrate,
lines = c(238, 1461, 2615))

# Finally, coerce the result to a 'GammaSpectra' object
multi_spc <- as(multi_scaled, "GammaSpectra")

# Plot the spectra
plot(multi_spc, xaxis = "energy") +
ggplot2::labs(title = "Calibrated spectra") +
ggplot2::theme_bw()

# Estimate the gamma dose rate

To estimate the gamma dose rate, you can either use one of the calibration curves distributed with this package or build your own.

# Load one of the built-in curves
data(BDX1) # IRAMAT-CRP2A (Bordeaux)
data(AIX1) # CEREGE (Aix-en-Provence)

The construction of a calibration curve requires a set of reference spectra for which the gamma dose rate is known and a background noise measurement. First, each reference spectrum is integrated over a given interval, then normalized to active time and corrected for background noise. The dose rate is finally modelled by the integrated signal value used as a linear predictor.

# Import the reference spectra
calib_dir <- system.file("extdata/crp2a/calibration", package = "gamma")
calib_spc <- read(calib_dir, skip = TRUE)

# Import the spectra for background noise measurement
noise_dir <- system.file("extdata/crp2a/background", package = "gamma")
noise_spc <- read(noise_dir, skip = TRUE)

# Set the known dose rates
setDoseRate(calib_spc) <- list(
BRIQUE = c(1986, 36),
C341 = c(850, 21),
C347 = c(1424, 24),
GOU = c(1575, 17),
LMP = c(642, 18),
MAZ = c(1141, 12),
PEP = c(2538, 112)
)

# Build the calibration curve
calib_curve <- fit(
calib_spc,
noise_spc,
range = c(165, 2800), # Set the integration range (in keV)
intercept = TRUE
)

# The linear model can be accessed for further investigations
# summary(calib_curve[["model"]])

# Plot the calibration curve
plot(calib_curve) +
ggplot2::theme_bw()

New dose rate values can be predicted.

# Import spectra
test_dir <- system.file("extdata/crp2a/test", package = "gamma")
test_spc <- read(test_dir, skip = TRUE)

# Estimate the gamma dose rate
dose_rate <- predict(calib_curve, test_spc, simplify = TRUE)
Gamma dose rate estimates
Dose Error
BR2011-08 313.632 11.137
BR2011-15A 195.302 6.939
BR2011-15B 193.735 6.881
BR2011-16 141.665 5.033
BR2011-20 1358.085 48.209
BR2011-23 393.568 13.972
BR2011-24 349.229 12.400
BR2011-26 338.581 12.021
BR2011-27 272.910 9.692
BR2011-29 235.469 8.361

# References

Morháč, Miroslav, Ján Kliman, Vladislav Matoušek, Martin Veselský, and Ivan Turzo. 1997. “Background Elimination Methods for Multidimensional Coincidence $$\gamma$$-Ray Spectra.” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 401 (1): 113–32. https://doi.org/10.1016/S0168-9002(97)01023-1.

Morháč, Miroslav, and Vladislav Matoušek. 2008. “Peak Clipping Algorithms for Background Estimation in Spectroscopic Data.” Applied Spectroscopy 62 (1): 91–106. https://doi.org/10.1366/000370208783412762.

Ryan, C.G., E. Clayton, W.L. Griffin, S.H. Sie, and D.R. Cousens. 1988. “SNIP, a Statistics-Sensitive Background Treatment for the Quantitative Analysis of PIXE Spectra in Geoscience Applications.” Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 34 (3): 396–402. https://doi.org/10.1016/0168-583X(88)90063-8.