From bin-files to De-distributions

Prerequisites

## load the library
library("Luminescence")

## set the working directory
setwd("/home/micha/Documents/Arbeit/Projects/R Luminescence/Tutorial from bin-file to De-plots")

Import and inspect bin-files

## import the bin-file
SAR.data <- readBIN2R("bin/TA110817N13-14_OSL_4_ed.BIN")
## 
## [readBIN2R.R]
##   >> bin/TA110817N13-14_OSL_4_ed.BIN
## 
   |=================================================================| 100%
##   >>  724 records have been read successfully!
## 
## 

## show a short summary of the data set
SAR.data
## 
## Risoe.BINfileData Object
##  Version:              03
##  Object Date:          060120
##  User:                 krb
##  System ID:            30
##  Overall Records:      724
##  Records Type:         TL=360; OSL=340; IRSL=24;
##  Position Range:       1 : 24
##  Run Range:            1 : 46
##  Set Range:            1 : 2

## show another summary, perhaps similar to what the Analyst shows for ten
## measurements parameters ID (1), SEL (2), LTYPE (7), POSITION (17), RUN
## (18), DTYPE (23) and IRR_TIME (24) are shown
SAR.data@METADATA[1:10, c(1, 2, 7, 17, 18, 23, 24)]
##    ID  SEL LTYPE POSITION RUN       DTYPE IRR_TIME
## 1   1 TRUE    TL        1   1     Natural        0
## 2   2 TRUE   OSL        1   2     Natural        0
## 3   3 TRUE    TL        1   4     Natural        0
## 4   4 TRUE   OSL        1   5 Bleach+dose       80
## 5   5 TRUE    TL        1   7 Bleach+dose        0
## 6   6 TRUE   OSL        1   8 Bleach+dose     1000
## 7   7 TRUE    TL        1  10 Bleach+dose        0
## 8   8 TRUE   OSL        1  11 Bleach+dose       80
## 9   9 TRUE    TL        1  13 Bleach+dose        0
## 10 10 TRUE   OSL        1  14 Bleach+dose     1800
## for a complete list of parameters and respective column numbers use
## this code
rbind(1:5, colnames(SAR.data@METADATA)[1:5], 6:10, colnames(SAR.data@METADATA)[6:10], 
    11:15, colnames(SAR.data@METADATA)[11:15], 16:20, colnames(SAR.data@METADATA)[16:20], 
    21:25, colnames(SAR.data@METADATA)[21:25], 26:30, colnames(SAR.data@METADATA)[26:30], 
    31:35, colnames(SAR.data@METADATA)[31:35], 36:40, colnames(SAR.data@METADATA)[36:40], 
    41:45, colnames(SAR.data@METADATA)[41:45])
##       [,1]          [,2]       [,3]          [,4]       [,5]      
##  [1,] "1"           "2"        "3"           "4"        "5"       
##  [2,] "ID"          "SEL"      "VERSION"     "LENGTH"   "PREVIOUS"
##  [3,] "6"           "7"        "8"           "9"        "10"      
##  [4,] "NPOINTS"     "LTYPE"    "LOW"         "HIGH"     "RATE"    
##  [5,] "11"          "12"       "13"          "14"       "15"      
##  [6,] "TEMPERATURE" "XCOORD"   "YCOORD"      "TOLDELAY" "TOLON"   
##  [7,] "16"          "17"       "18"          "19"       "20"      
##  [8,] "TOLOFF"      "POSITION" "RUN"         "TIME"     "DATE"    
##  [9,] "21"          "22"       "23"          "24"       "25"      
## [10,] "SEQUENCE"    "USER"     "DTYPE"       "IRR_TIME" "IRR_TYPE"
## [11,] "26"          "27"       "28"          "29"       "30"      
## [12,] "IRR_UNIT"    "BL_TIME"  "BL_UNIT"     "AN_TEMP"  "AN_TIME" 
## [13,] "31"          "32"       "33"          "34"       "35"      
## [14,] "NORM1"       "NORM2"    "NORM3"       "BG"       "SHIFT"   
## [15,] "36"          "37"       "38"          "39"       "40"      
## [16,] "SAMPLE"      "COMMENT"  "LIGHTSOURCE" "SET"      "TAG"     
## [17,] "41"          "42"       "43"          "44"       "45"      
## [18,] "GRAIN"       "LPOWER"   "SYSTEMID"    NA         NA

Analyse SAR-data

## some convenient settings
signal.integral <- 1:5
background.integral <- 200:250
position <- 1
info.measurement <- "Arbitrary sample 1"

## analyse the dataset for position 1
SAR.results <- Analyse_SAR.OSLdata(input.data = SAR.data, signal.integral = signal.integral, 
    background.integral = background.integral, position = position, info.measurement = info.measurement, 
    output.plot = TRUE)

plot of chunk unnamed-chunk-4

## 
## [Analyse_OSLCurves.R] >> Figure for position 1 produced.

## show some of the output, just for illustration
str(SAR.results)
## List of 3
##  $ LnLxTnTx         :List of 1
##   ..$ :'data.frame': 7 obs. of  15 variables:
##   .. ..$ Name          : chr [1:7] "Natural" "R1" "R2" "R3" ...
##   .. ..$ Dose          : num [1:7] 0 1000 1800 2200 3000 0 1800
##   .. ..$ Repeated      : logi [1:7] FALSE FALSE FALSE FALSE FALSE FALSE ...
##   .. ..$ LnLx          : num [1:7] 8014 7592 11796 14128 18165 ...
##   .. ..$ LnLx.BG       : num [1:7] 659 886 1126 1410 2020 ...
##   .. ..$ TnTx          : num [1:7] 1862 2006 2239 2393 2891 ...
##   .. ..$ TnTx.BG       : num [1:7] 614 748 873 1056 1341 ...
##   .. ..$ Net_LnLx      : num [1:7] 7355 6706 10670 12718 16145 ...
##   .. ..$ Net_LnLx.Error: num [1:7] 143 106 206 200 248 ...
##   .. ..$ Net_TnTx      : num [1:7] 1248 1258 1366 1337 1550 ...
##   .. ..$ Net_TnTx.Error: num [1:7] 56 74.5 58.4 63.5 87.2 ...
##   .. ..$ LxTx          : num [1:7] 4.3 3.78 5.27 5.9 6.28 ...
##   .. ..$ LxTx.Error    : num [1:7] 0.288 0.327 0.366 0.476 0.607 ...
##   .. ..$ LnLx.curveID  : num [1:7] 2 6 10 14 18 22 26
##   .. ..$ TnTx.curveID  : num [1:7] 4 8 12 16 20 24 28
##  $ RejectionCriteria: num [1, 1:3] 1.13 0.17 NA
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:3] "R2/R6" "Recuperation" "IRSL_BOSL"
##  $ SARParameters    :'data.frame':   2 obs. of  3 variables:
##   ..$ readTemp: int [1:2] 0 0
##   ..$ cutheat : num [1:2] 200 160
##   ..$ systemID: int [1:2] 30 30

Get De-values from growth-curves

## extract and show LxTx data, data frame conversion necessary for plot functions
data.LxTx <- as.data.frame(cbind(SAR.results[[1]][[1]][c(2,    # Dose
                                                         12,   # LxTx
                                                         13,   # LxTx.Error
                                                         6)])) # TnTx
data.LxTx
##   Dose  LxTx LxTx.Error TnTx
## 1    0 4.304    0.28838 1862
## 2 1000 3.785    0.32684 2006
## 3 1800 5.268    0.36604 2239
## 4 2200 5.904    0.47588 2393
## 5 3000 6.283    0.60718 2891
## 6    0 0.732    0.07193 2045
## 7 1800 4.654    0.46571 2829

## build growth curve
growth.curve <- plot_GrowthCurve(data.LxTx, 
                                 output.plot = TRUE)
## [plot_GrowthCurve.R] >> D01 = 1736.48

plot of chunk unnamed-chunk-5


## show fit parameters
growth.curve$Fit
## Nonlinear regression model
##   model:  y ~ fit.functionEXP(a, b, c, x) 
##    data:  data 
##       a       b       c 
##    7.55 1736.48  177.54 
##  weighted residual sum-of-squares: 0.0281
## 
## Algorithm "port", convergence message: relative convergence (4)

## extract and show De-value and delta De
De.data <- cbind(growth.curve$De[1:2])
De.data
##     De De.Error
## 1 1288    136.9

And now this entire work in a loop

## some convenient settings to analyse several positions/aliquots at once
signal.integral <- 1:5
background.integral <- 200:250
position <- 1:20
info.measurement <- "Arbitrary sample 1"

## analyse the dataset for positions 1 through 20
SAR.results <- Analyse_SAR.OSLdata(input.data = SAR.data,
                                   signal.integral = signal.integral,
                                   background.integral = background.integral,
                                   position = position,
                                   info.measurement = info.measurement,
                                   output.plot = FALSE)

## extract LxTx data and create De-values
De.data <- data.frame(De = NA, De.Error = NA)
for(i in 1:max(position)) {
data.LxTx <- as.data.frame(cbind(SAR.results[[1]][[i]][c(2,    # Dose
                                                         12,   # LxTx
                                                         13,   # LxTx.Error
                                                         6)])) # TnTx
growth.curve <- plot_GrowthCurve(data.LxTx, 
                                 output.plot = FALSE)

## extract and show De-value and delta De
De.data[i,] <- as.numeric(growth.curve$De[1:2])
}
## [plot_GrowthCurve.R] >> D01 = 1736.48
## [plot_GrowthCurve.R] >> D01 = 1207.56
## [plot_GrowthCurve.R] >> D01 = 1061.56
## [plot_GrowthCurve.R] >> D01 = 1588.12
## [plot_GrowthCurve.R] >> D01 = 1540.72
## [plot_GrowthCurve.R] >> D01 = 1406.09
## [plot_GrowthCurve.R] >> D01 = 1395.99
## [plot_GrowthCurve.R] >> D01 = 1856.74
## [plot_GrowthCurve.R] >> D01 = 1625.59
## [plot_GrowthCurve.R] >> D01 = 1361.68
## [plot_GrowthCurve.R] >> D01 = 1694.63
## [plot_GrowthCurve.R] >> D01 = 1723.4
## [plot_GrowthCurve.R] >> D01 = 1366.87
## [plot_GrowthCurve.R] >> D01 = 1375.56
## [plot_GrowthCurve.R] >> D01 = 1794.62
## [plot_GrowthCurve.R] >> D01 = 1300.68
## [plot_GrowthCurve.R] >> D01 = 1501.94
## [plot_GrowthCurve.R] >> D01 = 1257.66
## [plot_GrowthCurve.R] >> D01 = 1864.51
## [plot_GrowthCurve.R] >> D01 = 1600.62

Display De-values

## remove potential missing values from De-data set
De.data.na.rm <- De.data[complete.cases(De.data), ]

## plot KDE with some further information
plot_KDE(values = De.data.na.rm, distribution.parameters = c("median", "qr"), 
    summary = c("mean", "median", "sdabs"))

plot of chunk unnamed-chunk-7


## plot radial plot
plot_RadialPlot(sample = De.data.na.rm)

plot of chunk unnamed-chunk-7

Save the data

## save all R-internal data
save(SAR.data, SAR.results, De.data, De.data.na.rm, file = "Analysis.data.RData")

## save external ASCII output
write.table(x = De.data, file = "De-data.txt", row.names = FALSE)

## to read this or other external data use the following function (not
## run) De.data2 <- read.table('De-data.txt', header = TRUE)