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.

utils::vignette("CRP2A#1", package = "gamma")

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.

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.

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.

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") +

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.

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

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.

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


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.