EPR Spectra Data Analysis

Calibration Curve Co-60

# Integral simples dos espectros
# ref.  http://stackoverflow.com/questions/30969034/how-can-i-do-vector-integral-in-r #Jellen Vermeir
discreteIntegrationTrapeziumRule <- function(v,lower=1,upper=length(v),stepsize=1)
{
  if(upper > length(v))
    upper=length(v)
  if(lower < 1)
    lower=1
  integrand <- v[lower:upper]
  l <- length(integrand)
  stepsize*(0.5*integrand[1]+sum(integrand[2:(l-1)])+0.5*v[l])
}
############

1. Espectro Alanina-Gs (d141.DSC, d142.DSC, d143.DSC, d144.DSC, d145.DSC, d146.DSC, d147.DSC, d148.DSC, d149.DSC, d150.DSC, d151.DSC, d152.DSC, d153.DSC, d154.DSC, d155.DSC, d156.DSC, d157.DSC, d158.DSC, d159.DSC, d160.DSC, d161.DSC, d162.DSC, d163.DSC, d164.DSC, d165.DSC, d166.DSC, d167.DSC, d168.DSC, d169 50kGy.DSC, d169.DSC, d170 50kGy.DSC, d170.DSC, d171 50kGy.DSC, d171.DSC, d172 50kGy.DSC, d172.DSC, d173 100kGy.DSC, d174 100kGy.DSC, d175 100kGy.DSC, d176 100kGy.DSC)

• Espectro de EPR normalizado utilizando a equação da Bruker.

library(baseline)
#carregando a lista de arquivos DSC do dir
fileDSC<-list.files(pattern="*.DSC$")
fileDTA<-list.files(pattern="*.DTA$")
nn<-length(fileDSC)
length(fileDTA)

[1] 40

fileDSC

 [1] "d141.DSC"        "d142.DSC"        "d143.DSC"        "d144.DSC"        "d145.DSC"       
 [6] "d146.DSC"        "d147.DSC"        "d148.DSC"        "d149.DSC"        "d150.DSC"       
[11] "d151.DSC"        "d152.DSC"        "d153.DSC"        "d154.DSC"        "d155.DSC"       
[16] "d156.DSC"        "d157.DSC"        "d158.DSC"        "d159.DSC"        "d160.DSC"       
[21] "d161.DSC"        "d162.DSC"        "d163.DSC"        "d164.DSC"        "d165.DSC"       
[26] "d166.DSC"        "d167.DSC"        "d168.DSC"        "d169 50kGy.DSC"  "d169.DSC"       
[31] "d170 50kGy.DSC"  "d170.DSC"        "d171 50kGy.DSC"  "d171.DSC"        "d172 50kGy.DSC" 
[36] "d172.DSC"        "d173 100kGy.DSC" "d174 100kGy.DSC" "d175 100kGy.DSC" "d176 100kGy.DSC"

df<-data.frame(Arq=character(nn),Pos=character(nn),Dose=double(nn),Hpap=double(nn),In=double(nn),Irb=double(nn),IIRB=double(nn),stringsAsFactors = FALSE)
nome<-vector(mode="numeric", length=nn)
hpap<-vector(mode = "double",length = nn)
i1<-vector(mode="double",length = nn)
irb<-vector(mode="double",length = nn)
iirb<-vector(mode="double",length = nn)
pos<-vector(mode="character",length=nn)
for (i in 1:nn)
  {
  dsc <- readLines(fileDSC[i],encoding="UTF-8")
  #
  nome<-(strsplit(fileDSC[i],"\\."))[[1]]
  df$Arq[i]<-nome[1]
  #
  fileDTA<-paste0(nome[1],".DTA")
  pontos<-grep("XPTS",dsc,value=T);campos<-strsplit(pontos,'\t')[[1]];XPTS<-as.numeric(campos[2])
  pontos<-grep("XMIN",dsc,value=T);campos<-strsplit(pontos,'\t')[[1]];XMIN<-as.numeric(campos[2])
  pontos<-grep("XWID",dsc,value=T);campos<-strsplit(pontos,'\t')[[1]];XWID<-as.numeric(campos[2])
  pontos<-grep("MWFQ",dsc,value=T);campos<-strsplit(pontos,'    ')[[1]];MF<-as.numeric(campos[2])/1e9#GHz
  pontos<-grep("RCAG",dsc,value=T);campos<-strsplit(pontos,'    ')[[1]];RG<-as.numeric(campos[2])
  pontos<-grep("SPTP",dsc,value=T);campos<-strsplit(pontos,'    ')[[1]];CT<-as.numeric(campos[2])*1000#ms
  pontos<-grep("AVGS",dsc,value=T);campos<-strsplit(pontos,'    ')[[1]];SCAN<-as.numeric(campos[2])
  info<-grep("DATE", dsc,value=T);campos<-strsplit(info,'    ')[[1]];data1<-(campos[2])
  passo<-round(XWID/XPTS,1)
  b<-seq(1:XPTS)
  for(j in 1:XPTS){
    b[j]<-XMIN+(passo*(j-1))
  }
  b<-data.frame(b)
  g<-714.55*MF/b
  s<-readBin(fileDTA, double(), n = XPTS, size = 8, endian = "big");s<-data.frame(s);
  #Data normalization - Bruker Manual cap.5.6
  sn<-s/(CT*SCAN*20*10^(RG/20))
  #
  dados<-cbind(b,g,s,sn);colnames(dados)[2] <- "g";colnames(dados)[4] <- "sn"
  #
  df$Hpap[i]<-max(dados[,4])-min(dados[,4])
  df$In[i] <- discreteIntegrationTrapeziumRule(dados$sn,stepsize = 0.1)
  #
  #
  # Para a correção da linha de base foi utilizada a biblioteca "baseline"
  sn1t<-t(dados$sn)#
  #b1 <- baseline(sn1t[1,,drop=FALSE],method='irls')
  #b1 <- baseline(sn1t[1,, drop=FALSE], method='lowpass')#nao funciona
  #b1 <- baseline(sn1t[1,, drop=FALSE], method='rfbaseline',span=NULL, NoXP=1000)
  b1 <- baseline(sn1t[1,, drop=FALSE], hwm=300, method='medianWindow')
  #b1 <- baseline(sn1t[1,, drop=FALSE], lambda=10, method='als')
  #b1 <- baseline(sn1t[1,, drop=FALSE], method='modpolyfit', deg=6)
  #b1 <- baseline(sn1t[1,, drop=FALSE], wm=20,ws=20,method='rollingBall')#não serve
  #
  cor1<-t(getCorrected(b1))
  soma1<-data.frame(cumsum(cor1))
  df$Irb[i] <- discreteIntegrationTrapeziumRule(soma1$cumsum.cor1.,stepsize = 0.1)  
  #
  sn2<-data.frame(cumsum(dados$sn))
  sn2t<-t(sn2$cumsum.dados.sn.)
  #b2<-baseline(sn2t)
  #b2 <- baseline(sn2t[1,, drop=FALSE], method='lowpass')#nao funciona
  #b2 <- baseline(sn2t[1,, drop=FALSE], method='rfbaseline',span=NULL, NoXP=1000)
  b2 <- baseline(sn2t[1,, drop=FALSE], hwm=300, method='medianWindow')
  #b2 <- baseline(sn2t[1,, drop=FALSE], lambda=10, method='als')
  #b2 <- baseline(sn2t[1,, drop=FALSE], method='modpolyfit', deg=6)
  #b2<-b2mpol
  #
  cor2<-t(getCorrected(b2))
  df$IIRB[i] <- discreteIntegrationTrapeziumRule(cor2,stepsize = 0.1)  
}
df$Dose[1:4]<-1
df$Dose[5:8]<-10
df$Dose[9:12]<-50
df$Dose[13:16]<-100
df$Dose[17:20]<-500
df$Dose[21:24]<-1000
df$Dose[25:28]<-10000
df$Dose[29]<-50000
df$Dose[31]<-50000
df$Dose[33]<-50000
df$Dose[35]<-50000
df$Dose[37:40]<-100000
# retirado do set de irradiacoes no ACC
df$Dose[30]<-0
df$Dose[32]<-0
df$Dose[34]<-0
df$Dose[36]<-0
#pos
pos<-c("A","B","C","D","A","B","C","D","A","B","C","D","A","B","C","D","A","B","C","D","A","B","C","D","A","B","C","D","A","A","B","B","C","C","D","D","A","B","C","D")
df$Pos<-pos
plot(df$Hpap,df$In,col="blue")

plot(df$Hpap,df$Irb,col="blue",log="xy")

3 y values <= 0 omitted from logarithmic plot

points(df$Hpap,df$IIRB,col="red")

points(df$Hpap,df$In,col="green")

pairs(data.matrix(df[3:6]))

plot(df$Hpap,df$Irb,col="blue")

plot(df$Hpap,df$IIRB,col="green")

plot(df$Hpap,df$In,col="red")

cores<-c(A="red",B="blue",C="green",D="yellow")
plot(df$Hpap,df$Irb,col="blue",log="xy")

3 y values <= 0 omitted from logarithmic plot

plot(df$Hpap,df$IIRB,col="green",log="xy")

plot(df$Hpap,df$In,col="red",log="xy")

39 y values <= 0 omitted from logarithmic plot

cores

       A        B        C        D 
   "red"   "blue"  "green" "yellow" 

x1<-jitter(df$Irb,2,2)/10000
x2<-jitter(df$In,1,2)
plot(x1,x2,type="n")
text(x1,x2,labels = df$Pos)

# ordenando pela dose
df2<-df[order(df$Dose),]
library(data.table)
as.data.table(df2)[,mean(Hpap), by=.(Dose)]

as.data.table(df2)[,sd(Hpap), by=.(Dose)]

df3<-data.frame(as.data.table(df2)[,mean(Hpap), by=.(Dose)])
df3<-cbind(df3,as.data.table(df2)[,sd(Hpap), by=.(Dose)][,2])
colnames(df3)[1]<-"Dose"
colnames(df3)[2]<-"Hpapm"
colnames(df3)[3]<-"sd"
df3["sdr"]<-df3$sd/df3$Hpapm*100
library(ggplot2)
boxplot(df3$sdr~df3$Dose,
        xlab = "Dose",
        ylab = "Desvio Padrão Relativo (%)")

# library(wordcloud)
# library (calibrate)
# plot(log(df$Dose),log(df$Hpap),col="blue")
# plot(df$Dose,df$Hpap,col="red")
# plot(df$Dose,df$Hpap,col="green",log="xy")
# 
# j1<-df$Dose#jitter(df$Dose)
# j2<-df$Hpap# jitter(df$Hpap)
# plot(j1,j2,col="red",cex=4,xlim=c(0,2),ylim=c(0,0.5))
# 
# plot(jitter(df$Dose),jitter(df$Hpap),col="red",cex=4,xlim=c(0,100),ylim=c(0,15))
# text(df$Dose,df$Hpap,df$Arq,xlim=c(0,2),ylim=c(0,0.5))
#plot(NA, xlim=c(0,10), ylim=c(1, 10^4), xlab="x", ylab="y", log="y", yaxt="n")
#at.y <- outer(1:9, 10^(0:4))
#lab.y <- ifelse(log10(at.y) %% 1 == 0, at.y, NA)
#axis(2, at=at.y, labels=lab.y, las=1)

Orlando (nov/2016)