Fitting function for multi-exponentially decaying CW-OSL measurements, based on the algorithm described by Bluszcz & Adamiec (2006).
fit_OSLcurve(
curve,
K.max = 5,
F.threshold = 150,
stimulation.intensity = 30,
stimulation.wavelength = 470,
verbose = TRUE,
output.complex = FALSE,
parallel.computing = FALSE
)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 |
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
numeric (with default): Maximum number of components K. The computing time increases exponentially with the component number. K < 7 is recommended
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
numeric (with default): Intensity of optical stimulation in mW / cm². Used to calculate photoionisation cross sections.
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.
logical (with default): Enables console text output.
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.
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.
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
Dirk Mittelstraß, dirk.mittelstrass@luminescence.de
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.
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 |
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.
RLum.OSL_decomposition, sum_OSLcurves, decompose_OSLcurve, plot_OSLcurve, plot_PhotoCrosssections, minpack.lm::nlsLM, DEoptim::DEoptim
# 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)
Run the code above in your browser using DataLab