2.3 Equivalent dose calculation

2.3.1 Minimalist approach

In most cases the analysed data are from quartz SAR measurements, which requires the analyse_SAR.CWOSL() to be used for equivalent dose calculation. Even though analyse_SAR.CWOSL() offers many different options (i.e., arguments) to modify the fitting procedure, there are only few arguments that are actually required to run this function.

The object argument takes an RLum.Analysis object as input (or a list of them). Other than the input data themselves we only need to define the signal integration intervals.

de <- analyse_SAR.CWOSL(object = rlum[[3]], 
                        signal.integral.min = 1, 
                        signal.integral.max = 3, 
                        background.integral.min = 200, 
                        background.integral.max = 250, 
                        plot.single = 6)
## [plot_GrowthCurve()] Fit: EXP (interpolation) | De = 495.45 | D01 = 447.38

Note that we also used plot.single = 6 in the function call, which restricts the graphical output to the dose response curve. By default, the function would also plot the TL (if available), Tx and Lx curves as well as a series of plots for the rejection criteria.

2.3.2 Applying further criteria

When you have a look at the documentation of ?analyse_SAR.CWOSL you may notice that there is no obvious argument for setting the applied function to be fitted (single saturating exponential, linear, etc.). This is because this and related options are part of the plot_GrowthCurve() function, the actual workhorse of analyse_SAR.CWOSL(). Arguments provided by plot_GrowthCurve() can also be used for analyse_SAR.CWOSL(). This works because of the magical ... argument, which always indicates that a function accepts an indeterminate number of further arguments that are internally passed on to other functions. From ?plot_GrowthCurve we can learn that we must use the fit.method argument to change the fitting function. When we now use fit.method = "EXP+LIN" as input in our analyse_SAR.CWOSL() call, the argument will be automatically passed to the internally called plot_GrowthCurve() function.

In the following code block we want to further customise our fitting procedure and also adjust the rejection criteria. From the documentation we can learn that the argument rejection.criteria takes a named list with numerical values. You can use the following keywords as names for the list elements:

  • recycling.ratio
  • recuperation.rate
  • testdose.error
  • palaeodose.error
  • exceed.max.regpoint

Before adjusting the threshold values we first have a look at the performance of our previous fit.

print(get_RLum(de, "rejection.criteria"))
##                  Criteria       Value Threshold Status
## 1 Recycling ratio (R5/R1)   1.0337000       0.1     OK
## 2     Recuperation rate 1   0.0165000       0.1     OK
## 3          Testdose error   0.0119357       0.1     OK
## 4        Palaeodose error   0.0290200       0.1     OK
## 5    De > max. dose point 495.4500000     850.0     OK
##                                  UID
## 1 2018-02-16-03:38.0.465975960250944
## 2 2018-02-16-03:38.0.465975960250944
## 3 2018-02-16-03:38.0.465975960250944
## 4 2018-02-16-03:38.0.465975960250944
## 5 2018-02-16-03:38.0.465975960250944

The data.frame tells us that the aliquot passed all rejection criteria (“OK” status). For the purpose of demonstrating we will apply much stricter threshold values so that our aliquot fails at least one of the criteria.

rc <- list(
  recycling.ratio = 2,
  recuperation.rate = 5,
  testdose.error = 5,
  palaeodose.error = 5,
  exceed.max.regpoint = TRUE)

de <- analyse_SAR.CWOSL(object = rlum[[3]], 
                        signal.integral.min = 1, 
                        signal.integral.max = 3, 
                        background.integral.min = 200, 
                        background.integral.max = 250, 
                        plot = FALSE, 
                        verbose = FALSE,
                        rejection.criteria = rc,
                        # arguments inherited from plot_GrowthCurve()
                        fit.method = "EXP+LIN",
                        mode = "interpolation",
                        fit.force_through_origin = FALSE,
                        fit.includingRepeatedRegPoints = TRUE,
                        fit.weights = TRUE,
                        NumberIterations.MC = 100)

print(get_RLum(de, "rejection.criteria"))
##                  Criteria       Value Threshold Status
## 1 Recycling ratio (R5/R1)   1.0337000      0.02 FAILED
## 2     Recuperation rate 1   0.0165000      0.05     OK
## 3          Testdose error   0.0119357      0.05     OK
## 4        Palaeodose error   0.0287400      0.05     OK
## 5    De > max. dose point 532.1000000    850.00     OK
##                                  UID
## 1 2018-02-16-03:38.0.712970336666331
## 2 2018-02-16-03:38.0.712970336666331
## 3 2018-02-16-03:38.0.712970336666331
## 4 2018-02-16-03:38.0.712970336666331
## 5 2018-02-16-03:38.0.712970336666331

Due to the much stricter threshold values the aliquot failed the recycling ratio criteria. While it is certainly nice to know why the aliquot failed the rejection criteria, for the subsequent analysis (calculating a mean dose) it is often enough to know that it failed in the first place. This is why the summary table that also contains the equivalent dose also provides the field RC.Status.

print(get_RLum(de))
##      De De.Error D01 D01.ERROR D02 D02.ERROR  De.MC     Fit RC.Status
## 1 532.1    15.29  NA        NA  NA        NA 532.93 EXP+LIN    FAILED
##   signal.range background.range signal.range.Tx background.range.Tx
## 1        1 : 3        200 : 250         NA : NA             NA : NA
##                                  UID
## 1 2018-02-16-03:38.0.712970336666331

As per the previous output table RC.Status contains the character string "FAILED", which we can later use to filter the data before applying any further analysis.

2.3.3 Iterating over all aliquots/grains

Usually we do not want to calculate the equivalent dose of an individual aliquot, but for all of the aliquots of a bin file. analyse_SAR.CWOSL() allows us to also provide a list of RLum.Analysis objects, so we do not need to bother with a for-loop or the *apply() function family.

de <- analyse_SAR.CWOSL(object = rlum, 
                        signal.integral.min = 1, 
                        signal.integral.max = 3, 
                        background.integral.min = 200, 
                        background.integral.max = 250, 
                        plot = FALSE, 
                        verbose = FALSE)

de_df <- get_RLum(de)

print(de_df)
##        De De.Error    D01 D01.ERROR D02 D02.ERROR  De.MC Fit RC.Status
## 1  582.87    55.36 302.71 14.161031  NA        NA 585.60 EXP        OK
## 2  411.10     5.06 710.58 16.355284  NA        NA 411.19 EXP        OK
## 3  495.45    14.47 447.38 14.352758  NA        NA 496.01 EXP        OK
## 4  407.79     7.68 523.85 14.783003  NA        NA 407.39 EXP        OK
## 5  440.33     7.15 619.88 21.670377  NA        NA 440.02 EXP        OK
## 6  569.63     5.41 646.92 10.472932  NA        NA 570.53 EXP        OK
## 7  425.10    16.63 341.88  9.799217  NA        NA 424.20 EXP        OK
## 8  337.63    18.49 536.04 54.248322  NA        NA 338.24 EXP        OK
## 9  709.25    19.55 530.91 15.232269  NA        NA 712.73 EXP        OK
## 10 475.90     8.61 502.61 13.271291  NA        NA 477.30 EXP        OK
##    signal.range background.range signal.range.Tx background.range.Tx
## 1         1 : 3        200 : 250         NA : NA             NA : NA
## 2         1 : 3        200 : 250         NA : NA             NA : NA
## 3         1 : 3        200 : 250         NA : NA             NA : NA
## 4         1 : 3        200 : 250         NA : NA             NA : NA
## 5         1 : 3        200 : 250         NA : NA             NA : NA
## 6         1 : 3        200 : 250         NA : NA             NA : NA
## 7         1 : 3        200 : 250         NA : NA             NA : NA
## 8         1 : 3        200 : 250         NA : NA             NA : NA
## 9         1 : 3        200 : 250         NA : NA             NA : NA
## 10        1 : 3        200 : 250         NA : NA             NA : NA
##                                    UID
## 1   2018-02-16-03:38.0.452163838082924
## 2   2018-02-16-03:38.0.753062088275328
## 3  2018-02-16-03:38.0.0190066744107753
## 4  2018-02-16-03:38.0.0946386174764484
## 5    2018-02-16-03:39.0.20289219962433
## 6   2018-02-16-03:39.0.245586855802685
## 7   2018-02-16-03:39.0.722360961837694
## 8   2018-02-16-03:39.0.642113087465987
## 9   2018-02-16-03:39.0.891926837153733
## 10  2018-02-16-03:39.0.668478885432705

Even though we provided a list of 10 RLum.Analysis object (one for each aliquot) we are only returned a single RLum.Results object. The data table now contains 10 rows, however, with the results for each of the aliquots.