source("ratio.R")
read.pipidata <- function(path="./", T, Rformat=FALSE, Tformat=FALSE) {
pion.cor <- cf()
pipi.cor <- cf()
if(Rformat) {
load(paste(path, "pion.cor.Rdata", sep="") )
load(paste(path, "pipi.cor.Rdata", sep="") )
}
else if(Tformat) {
ind.vector <- c(2,2)
pipi.cor <- readtextcf("pipi_pipi_A1_corr_TP0_00.dat", T=T, check.t=1, path=path, ind.vector=ind.vector)
pion.cor <- readtextcf("pi_corr_p0.dat", T=T, check.t=1, path=path, ind.vector=ind.vector)
}
else {
files <- Sys.glob( paste(path, "C2_pi+-*", sep="") )
pion.cor <- readbinarycf(files, obs=0, T=T)
files <- Sys.glob( paste(path, "C4_1*", sep="") )
pipi1 <- readbinarycf(files, obs=0, T=T)
files <- Sys.glob( paste(path, "C4_2*", sep="") )
pipi2 <- readbinarycf(files, obs=0, T=T)
files <- Sys.glob( paste(path, "C4_3*", sep="") )
pipi3 <- readbinarycf(files, obs=0, T=T)
pipi3 <- mul.cf(pipi3, a=2.)
rm(files)
save(pion.cor, file=paste(path, "pion.cor.Rdata", sep="") )
save(pipi1, file=paste(path, "pipi1.Rdata", sep="") )
save(pipi2, file=paste(path, "pipi2.Rdata", sep="") )
save(pipi3, file=paste(path, "pipi3.Rdata", sep="") )
pipi.cor <- pipi1 + pipi2 - pipi3
save(pipi.cor, file=paste(path, "pipi.cor.Rdata", sep="") )
}
return(invisible(list(pion.cor=pion.cor, pipi.cor=pipi.cor)))
}
## Luescher Formula to Order 1/L^5
## a0 from deltaE and L
## using
## deltaE = - 4 pi a0 / m/L^3 (1 -2.837297 a0/L + 6.375183 a0^2/L^2)
deltaEvL <- function(a0, deltaE, L, m, debug=FALSE) {
if(debug) cat("1/L^4:", -2.837297*a0/L, " 1/L^5:", 6.375183*a0^2/L^2, "\n")
return(deltaE + 4*pi*a0/(m*L^3)*(1 - 2.837297*a0/L + 6.375183*a0^2/L^2))
##return(deltaE + 4*pi*a0/(m*L^3)*(1 - 2.837297*a0/L + 0*6.375183*a0^2/L^2))
##return(deltaE + 4*pi*a0/(m*L^3))
}
## extrac scattering length Luescher formula using
## the ratio to determine energy shift values
run.pipi.analysis.ratio <- function(data, boot.R=999, boot.l=1, useCov=TRUE,
t1=13, t2=31, tr1, tr2, L, ens="") {
if(missing(t1) || missing(t2)) {
stop("t1 and t2 must be specified, aborting...\n")
}
if(missing(tr1)) {
tr1 <- t1
}
if(missing(tr2)) {
tr2 <- t2
}
if(missing(data)) {
stop("data is missing, aborting...\n")
}
pion.cor <- data$pion.cor
pipi.cor <- data$pipi.cor
Thalf <- pion.cor$Time/2
if(missing(L)) {
L <- Thalf
}
if(t2 > Thalf || tr2 > Thalf) {
stop("t2 or tr2 not in range, aborting...\n")
}
pipi.cor <- bootstrap.cf(pipi.cor, boot.R=boot.R, boot.l=boot.l)
pion.effmass <- pion.cor
if(!inherits(pion.cor, "effectivemassfit")) {
pion.cor <- bootstrap.cf(pion.cor, boot.R=boot.R, boot.l=boot.l)
pion.effmass <- bootstrap.effectivemass(pion.cor, boot.R=boot.R, boot.l=boot.l, type="solve")
pion.effmass <- fit.effectivemass(pion.effmass, t1=t1, t2=t2, useCov=useCov)
}
## this is R(t+1/2)
Rpipi <- compRpipi(c4=pipi.cor$cf0, c2=pion.cor$cf0, Thalf=Thalf)
Rpipi.tsboot <- array(NA, dim=c(boot.R, Thalf))
for(i in c(1:boot.R)) {
Rpipi.tsboot[i,] <- compRpipi(pipi.cor$cf.tsboot$t[i,], pion.cor$cf.tsboot$t[i,], Thalf=Thalf)
}
## error of the ratio
dRpipi <- apply(Rpipi.tsboot, 2, sd)
## chi^2 function includes covariance matrix M
fitfn <- function(par, y, t, m, Thalf, M) {
tp <- t - Thalf
z <- Rfn(par, tp, m)
return((z-y) %*% M %*% (z-y))
}
## the same, but in the format for nls.lm
fitfn2 <- function(par, y, t, m, Thalf, L) {
z <- Rfn(par, t-Thalf, m)
return(L %*% (z-y))
}
## fit range
tt <- c(tr1:tr2)
## ratio is displaced by 1/2
tphys <- tt-0.5
## build the covariance matrix from bootstrap samples, if wanted
M <- diag(1./dRpipi[tt]^2)
if(useCov) M <- invertCovMatrix(Rpipi.tsboot[,tt], boot.samples=TRUE)
par <- c(1.6,0.01)
## fit on the original data
lm.avail <- require(minpack.lm)
if(lm.avail) {
ML <- chol(M)
opt.res <- nls.lm(par, fn=fitfn2, L=ML, Thalf=Thalf, t=tphys, y=Rpipi[tt], m=pion.effmass$opt.res$par[1])
}
else {
opt.res <- optim(par, fn = fitfn, method="BFGS", M=M, Thalf=Thalf, t=tphys, y=Rpipi[tt], m=pion.effmass$opt.res$par[1])
opt.res <- optim(opt.res$par, fn = fitfn, method="BFGS", M=M, Thalf=Thalf, t=tphys, y=Rpipi[tt],
m=pion.effmass$opt.res$par[1], control=list(parscale=1./opt.res$par,ndeps=rep(1.e-6, times=length(par))))
}
data$chisq <- fitfn(opt.res$par, M=M, Thalf=Thalf, t=tphys, y=Rpipi[tt], m=pion.effmass$opt.res$par[1])
a0 <- uniroot(deltaEvL, c(0.,-6.), deltaE=opt.res$par[2], L=L, m=pion.effmass$opt.res$par[1])$root
mpia0 <- pion.effmass$opt.res$par[1]*a0
qtilde <- compute.qtildesq(E=opt.res$par[2]+2*pion.effmass$opt.res$par[1], L=L, dvec=c(0,0,0), mpi=pion.effmass$opt.res$par[1])
qcotdelta <- 2.*Re( LuescherZeta(qtilde$qtsq, l=0, m=0, gamma=qtilde$gammaboost, dvec=c(0,0,0)) )/(qtilde$gammaboost*L*sqrt(pi))
par <- opt.res$par
opt.tsboot <- array(NA, dim=c(boot.R, 7))
for(i in 1:boot.R) {
if(lm.avail) {
opt <- nls.lm(par, fn=fitfn2, L=ML, Thalf=Thalf, t=tphys, y=Rpipi.tsboot[i,tt], m=pion.effmass$massfit.tsboot[i,1])
}
else {
opt <- optim(par, fn = fitfn,
t=tphys, Thalf=Thalf, control=list(ndeps=rep(1.e-6, times=length(par)), parscale=1/par),
method="BFGS", M=M, y = Rpipi.tsboot[i,tt], m=pion.effmass$massfit.tsboot[i,1])
}
opt.tsboot[i, c(1,2)] <- opt$par
opt.tsboot[i, 3] <- fitfn(opt$par, t=tphys, Thalf=Thalf, M=M, y = Rpipi.tsboot[i,tt], m=pion.effmass$massfit.tsboot[i,1])
if(opt$par[2] < 0) opt.tsboot[i, 5] <- NA
else opt.tsboot[i, 5] <- uniroot(deltaEvL, c(0.,-6.), deltaE=opt$par[2], L=L, m=pion.effmass$massfit.tsboot[i,1])$root
opt.tsboot[i, 4] <- pion.effmass$massfit.tsboot[i,1]*opt.tsboot[i, 5]
}
qtildeboot <- compute.qtildesq(E=opt.tsboot[,2]+2*pion.effmass$massfit.tsboot[,1], L=L, dvec=c(0,0,0), mpi=pion.effmass$massfit.tsboot[,1])
qcotdeltaboot <- 2.*Re( LuescherZeta(qtildeboot$qtsq, l=0, m=0, gamma=qtildeboot$gammaboost, dvec=c(0,0,0)) )/(qtildeboot$gammaboost*L*sqrt(pi))
## this is deltaE as a function of a0, m and L
lfn2 <- function(a0, L, m, debug=FALSE) {
if(debug) cat("1/L^4:", -2.837297*a0/L, " 1/L^5:", 6.375183*a0^2/L^2, "\n")
return(- 4*pi*a0/(m*L^3)*(1 - 2.837297*a0/L + 6.375183*a0^2/L^2))
##return(deltaE + 4*pi*a0/(m*L^3)*(1 - 2.837297*a0/L + 0*6.375183*a0^2/L^2))
##return(deltaE + 4*pi*a0/(m*L^3))
}
data$pion.cor <- pion.cor
data$pipi.cor <- pipi.cor
data$Rpipi <- Rpipi
data$dRpipi <- dRpipi
data$opt.res <- opt.res
data$opt.tsboot <- opt.tsboot
data$t1 <- t1
data$t2 <- t2
data$useCov <- useCov
data$tr1 <- tr1
data$tr2 <- tr2
data$boot.R <- boot.R
data$boot.l <- boot.l
data$L <- L
data$T <- pion.cor$Time
data$ens <- ens
data$Ecm <- qtilde$Ecm
data$Ecmboot <- qtildeboot$Ecm
data$qcotdelta <- qcotdelta
data$qcotdeltaboot <- qcotdeltaboot
## data$chisq <- opt.res$value
data$dof <- length(tt)-2
data$Qval <- 1-pchisq(data$chisq, data$dof)
data$pionQval <- pion.effmass$Qval
data$pion.effmass <- pion.effmass
data$mpia0 <- mpia0
data$a0 <- a0
data$mpi <- pion.effmass$opt.res$par[1]
data$lfn <- lfn2
data$tphys <- tphys
attr(data, "class") <- c("pipi", class(data))
return(invisible(data))
}
summary.pipi <- function(pipi) {
summary(pipi$pion.effmass)
cat("\n ** Luescher Analysis **\n\n")
cat("fitrange from", pipi$tr1-0.5, "to", pipi$tr2-0.5, "\n")
cat("correlated fit = ", pipi$useCov, "\n")
chisq <- pipi$opt.res$value
dof <- pipi$dof
cat("chi^2 = ", pipi$chisq, "\n")
cat("dof = ", pipi$dof, "\n")
cat("chi^2/dof = ", pipi$chisq/pipi$dof, "\n")
cat("Qval = ", pipi$Qval, "\n")
cat("Qval of pion fit = ", pipi$pionQval, "\n\n")
# cat("contributions by the orders in L relativ to L cubed\n")
# data$lfn(a0=pipi$mpia0/pipi$mpi, m=pipi$mpi, L=pipi$L, debug=TRUE)
# cat("\n")
cat("deltaE = ", pipi$opt.res$par[2], "+-", sd(pipi$opt.tsboot[,2]), "(", 100*sd(pipi$opt.tsboot[,2])/pipi$opt.res$par[2], "%)", "\n")
cat("a_0 = ", pipi$a0, "+-", sd(pipi$opt.tsboot[, 5], na.rm=TRUE), "(", 100*sd(pipi$opt.tsboot[, 5], na.rm=TRUE)/abs(pipi$a0), "%)\n")
cat("mpi*a_0 = ", pipi$mpia0, "+-", sd(pipi$opt.tsboot[, 4], na.rm=TRUE), "(", 100*sd(pipi$opt.tsboot[, 4], na.rm=TRUE)/abs(pipi$mpia0), "%)\n")
cat("qcotdelta = ", pipi$qcotdelta, "+-", sd(pipi$qcotdeltaboot, na.rm=TRUE), "(", 100*sd(pipi$qcotdeltaboot, na.rm=TRUE)/abs(pipi$qcotdelta), "%)\n")
cat("Ecm = ", pipi$Ecm, "+-", sd(pipi$Ecmboot, na.rm=TRUE), "(", 100*sd(pipi$Ecmboot, na.rm=TRUE)/abs(pipi$Ecm), "%)\n")
}
plot.pipi <- function(pipi, ylim, ...) {
plot(pipi$pion.effmass, ylim=c(0.9*pipi$pion.effmass$opt.res$par[1], 1.1*pipi$pion.effmass$opt.res$par[1]), xlab=c("t/a"), ylab=c("aMeff"), main=paste("pion effectivemass ", pipi$t1, "-", pipi$t2, " p=", pipi$pionQval, sep=""))
if(interactive()) X11()
qqnorm(pipi$opt.tsboot[,1], main=paste("qqnorm amplitude ", pipi$tr1, "-", pipi$tr2, " p=", pipi$Qval, sep=""))
if(interactive()) X11()
qqnorm(pipi$opt.tsboot[,2], main=paste("qqnorm deltaE ", pipi$tr1, "-", pipi$tr2, " p=", pipi$Qval, sep=""))
if(interactive()) X11()
qqplot(pipi$opt.tsboot[,3], rchisq(n=pipi$boot.R, df=pipi$dof), main=paste("qqplot chisq", pipi$tr1, "-", pipi$tr2, " p=", pipi$Qval, sep=""))
Thalf <- pipi$T/2
if(missing(ylim)) {
rg <- range(pipi$Rpipi[c(5:Thalf)])
# ylim <- c(rg[1]-0.1*(rg[2]-rg[1]), rg[2]+0.1*(rg[2]-rg[1]))
ylim <- c(0.96*rg[1], 1.04*rg[2])
}
if(interactive()) X11()
plotwitherror(c(1:(Thalf))-0.5, pipi$Rpipi, pipi$dRpipi, ylim=ylim, xlim=c(5,Thalf), ylab=c("R(t)"), xlab=c("t"), main=paste("Ratio", pipi$tr1, "-", pipi$tr2, " p=", pipi$Qval, sep=""), ...)
x <- seq(pipi$tphys[1], pipi$tphys[length(pipi$tphys)]+1, 0.01)
rat <- Rfn(pipi$opt.res$par, tp=x-Thalf, m=pipi$mpi)
lines(x=x, y=rat, col="red")
write.table(data.frame(x=x, rat=rat), file=paste(pipi$ens, ".", pipi$tr1, "-", pipi$tr2,".ratiofit.dat", sep=""), quote=FALSE, row.names=FALSE, col.names=FALSE)
write.table(data.frame(t=c(1:Thalf)-0.5, pipi$Rpipi, pipi$dRpipi), file=paste(pipi$ens, ".", pipi$tr1, "-", pipi$tr2, ".ratio.dat", sep=""), quote=FALSE, row.names=FALSE, col.names=FALSE)
x <- seq(0, 0.06, 0.001)
dE <- pipi$lfn(pipi$mpia0/pipi$mpi, 1/x, pipi$mpi)
write.table(data.frame(x=x, dE=dE), file=paste(pipi$ens, ".", pipi$tr1, "-", pipi$tr2, ".dEovL.dat", sep=""), quote=FALSE, row.names=FALSE, col.names=FALSE)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.