fit_OSLcurve {OSLdecomposition}

Multi-exponential CW-OSL curve fitting

Description

Fitting function for multi-exponentially decaying CW-OSL measurements, based on the algorithm described by Bluszcz & Adamiec (2006).

Usage

fit_OSLcurve(
  curve,
  K.max = 5,
  F.threshold = 150,
  stimulation.intensity = 30,
  stimulation.wavelength = 470,
  verbose = TRUE,
  output.complex = FALSE,
  parallel.computing = FALSE
)

Arguments

curve

RLum.Data.Curve or data.frame or matrix (required): CW-OSL record or average CW-OSL curve created by sum_OSLcurves. If no column ⁠$time⁠ exists, the first column is defined as measurement time (x-axis). Time intervals must be constant. If no column ⁠$signal⁠ exists, the second column is defined as signal values (y-axis). Further columns will be ignored

K.max

numeric (with default): Maximum number of components K. The computing time increases exponentially with the component number. K < 7 is recommended

F.threshold

numeric (with default): Fitting stop criterion. If the F-value is lower than this threshold, the fitting procedure stops and the K - 1 fit is returned

stimulation.intensity

numeric (with default): Intensity of optical stimulation in mW / cm². Used to calculate photoionisation cross sections.

stimulation.wavelength

numeric (with default): Wavelength of optical stimulation in nm. Used to calculate photoionisation cross sections. If a wavelength between 465 and 480 nm is chosen, the cross sections are set into relation with literature values to name the signal components automatically.

verbose

logical (with default): Enables console text output.

output.complex

logical (with default): If TRUE, the function returns a list of objects, see section Value for further information. If FALSE, the function returns a data.frame with the CW-OSL model parameters of the fitting chosen by the F-test. Setting the parameter to FALSE is not recommended when fitting a global average curve created by sum_OSLcurves as over-fitting is likely in such cases.

parallel.computing

logical (with default): Enables the use of multiple CPU cores. This increases the execution speed significantly but may need administrator rights and/or a firewall exception. See DEoptim::DEoptim.control for further information.

Details

The function assumes multiple exponentially decaying signal components with first-order kinetics:

I(t) = n_1 \lambda_1 exp(-\lambda_1 t) + n_2 \lambda_2 exp(-\lambda_2 t) + ... + n_K \lambda_K exp(-\lambda_K t)

with I(t) the CW-OSL signal, n the signal component intensity, \lambda the signal component decay constant and K the number of signal components. For actual fitting, the integrated version of this formula is used, see Mittelstraß et al. (2021) for details.

The fitting algorithm is an implementation of the hybrid evolutionary-linear algorithm (HELA) by Bluszcz & Adamiec (2006). See there or Mittelstraß et al. (in preparation) for details. The differential evolution part of HELA is performed by DEoptim::DEoptim. The linear regression part of HELA is performed by decompose_OSLcurve. The parameter refinement by Levenberg-Marquardt fitting is performed by minpack.lm::nlsLM.

F-test

Bluszcz & Adamiec (2006) suggest the use of an F-test to determine the correct number of signal components. This function compares the residual square sum (RSS_K) value of each fitting with the RSS_K-1 value of the previous fitting and calculates an Improvement-in-fitting-quality criterion:

F_K = {(RSS_{K-1} - RSS_K)/2} / {RSS_K(N - 2K)}

Here, N is the number data points (channels) of the measurement and K is the number of OSL components in the fitting model. If F_K falls below the threshold value (F.threshold), the fitting model with K components is apparently not significantly better than the K - 1 model of the previous fitting cycle. Thus, the K - 1 model will be recommended as fitting solution.

Photoionisation cross sections

While the function is suited for the analysis of a wide variety of multi-exponential decay problems, it is targeted to CW-OSL measurements of quartz under SAR protocol conditions (470 nm stimulation at 125 °C). To compare the calculated OSL components with OSL components reported in published literature, photoionisation cross sections are calculated using the stimulation.wavelength \lambda_{stim} and stimulation.intensity \Phi_{stim}:

\sigma_k=\lambda_k {hc / \Phi_{stim}\lambda_{stim}}

Here \sigma_k is the photoionisation cross section of component k in cm^2, \lambda_k the CW-OSL decay constant in s^-1, h the Planck constant and c the speed of light.

If a stimulation.intensity between 460 nm and 485 nm is defined, the components are named automatically in accordance to the cross-sections published by Durcan and Duller (2011), Jain et al. (2003) and Singarayer and Bailey (2003). For the Ultrafast and the Slow4 component, no consistent literature values could be found, so their range is tentatively assigned:

Component	Lower limit (cm^2)	Upper limit (cm^2)
Ultrafast	1e-16	1e-15
Fast	1.9e-17	3.1e-17
Medium	3e-18	9e-18
Slow1	1e-18	1.85e-18
Slow2	1.1e-19	4e-19
Slow3	1e-20	4.67e-20
Slow4	1e-21	1e-20

Value

If output.complex = FALSE, a data.frame is returned. It contains the signal decay rates and signal intensities of the best fit. The best fit was either chosen by the F-test or the last successful fitting iteration.

If output.complex = TRUE, a list of objects is returned:

Element	Type	Description
decay.rates	numeric	vector of the best suiting decay rates
K.selected	numeric	number of components of the best fit
F.test	data.frame	table containing the F-test parameter and the decay rates of each fitting model
F.test.print	data.frame	the same table as above, but formated for pretty console and report output
info.text	character	collected messages from the algorithms
component.tables	list	result data.frames for all tested models
curve	data.frame	fitted time-signal-curve
components	data.frame	best fit; same object as output.complex = FALSE returns
fit.results	list	list of nls objects for all tested models
plot.data	data.frame	factorized results for overview plotting with plot_PhotoCrosssections
parameters	list	function arguments and the needed computing time

Last update

2022-07-27, DM: Moved residual sum of squares (RSS) calculation during DE-optimization cycle to decompose_OSLcurve() to improve computing time by factor 3 to 4

Author(s)

Dirk Mittelstraß, [e-mail]

Please cite the package the following way:

Mittelstraß, D., Schmidt, C., Beyer, J., Heitmann, J. and Straessner, A.: R package OSLdecomposition: Automated identification and separation of quartz CW-OSL signal components, in preparation.

References

Bluszcz, A., Adamiec, G., 2006. Application of differential evolution to fitting OSL decay curves. Radiation Measurements 41, 886–891.

Durcan, J.A., Duller, G.A.T., 2011. The fast ratio: A rapid measure for testing the dominance of the fast component in the initial OSL signal from quartz. Radiation Measurements 46, 1065–1072.

Jain, M., Murray, A.S., Bøtter-Jensen, L., 2003. Characterisation of blue-light stimulated luminescence components in different quartz samples: implications for dose measurement. Radiation Measurements 37, 441–449.

Mittelstraß, D., 2019. Decomposition of weak optically stimulated luminescence signals and its application in retrospective dosimetry at quartz (Master thesis). TU Dresden, Dresden.

Singarayer, J.S., Bailey, R.M., 2003. Further investigations of the quartz optically stimulated luminescence components using linear modulation. Radiation Measurements, Proceedings of the 10th international Conference on Luminescence and Electron-Spin Resonance Dating (LED 2002) 37, 451–458.

See Also

RLum.OSL_decomposition, sum_OSLcurves, decompose_OSLcurve, plot_OSLcurve, plot_PhotoCrosssections, minpack.lm::nlsLM, DEoptim::DEoptim

Examples

# Create a simple curve with just one component
curve <- data.frame(
  X = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12),
  Y = c(377, 244, 163, 93, 59, 28, 17, 13, 10, 8, 9, 5))
# Perform fitting
components <- fit_OSLcurve(curve, F.threshold = 3)

# Display results
plot_OSLcurve(curve, components)