Lvalues <- c(32,24,20)
ulim <- c(0.8, 1.0, 1.2)
load("Energies.Rdata")
for(L in Lvalues) {
i <- which(Lvalues == L)
tikzfiles <- tikz.init(paste("energylevels", L, sep=""),width=5,height=5)
ii <- which(Epipi$L == L)
mirrep <- max(Epipi[ii,]$irrep)
mlevel <- max(Epipi[ii,]$level)
mframe <- max(Epipi[ii,]$frame)
N <- length(Epipi[Epipi$L==L & Epipi$level==1,]$E)
xv <- (Epipi[ii,]$irrep-1)*mframe + Epipi[ii,]$frame
xv[xv==9] <- 6
E <- Epipi[Epipi$L==L,]$E
dE <- Epipi[Epipi$L==L,]$dE
Mpi <- Epi[Epipi$L==L,]$E[1]
p1 <- array(c(0,0,0, 0,0,1, 0,1,1, 1,1,1, 0,0,1, 0,1,1), dim=c(3,6))
p2 <- array(c(0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,1,1), dim=c(3,6))
p12 <- array(c(0,0,1, 1,0,0, 0,1,0, 1,0,0, 0,1,1), dim=c(3,5))
p22 <- array(c(0,0,1, -1,0,1, 0,0,1, 0,1,1, 0,1,1), dim=c(3,5))
plot(NA, type="n", ylab="$aE$", xlab="", ylim=c(0,ulim[i]), xlim=c(0.5,N+0.5), xaxt="n")
rect(xleft=xv-0.2, xright=xv, ytop=E+dE, ybottom=E-dE, col="gray")
legend("topright", legend=c("$E_{\\pi\\pi}$"), bty="n", fill=c("gray"))
labels <- c("A1, $\\vec d=[0,0,0]$", "A1, $\\vec d=[0,0,1]$", "A1, $\\vec d=[0,1,1]$", "A1, $\\vec d=[1,1,1]$", "E, $\\vec d=[0,0,0]$", "T2, $\\vec d=[0,0,0]$")
axis(1, at=c(1:N), labels=rep("", times=N))
text(x=c(1:N), y=-0.09, labels=labels, srt = 90, adj = 1, xpd = TRUE,cex=0.8,srt=60)
legend("topleft", legend=c(paste("A40 L =", L)), bty="n", cex=1.2)
disp <- c()
disp2 <- c()
for(i in c(1:6)) {
disp[i] <- acosh( cosh(Mpi) + 2*sum(sin(pi*p1[,i]/L)^2)) + acosh( cosh(Mpi) + 2*sum(sin(pi*p2[,i]/L)^2))
}
for(i in c(1:5)) {
disp2[i] <- acosh( cosh(Mpi) + 2*sum(sin(pi*p12[,i]/L)^2)) + acosh( cosh(Mpi) + 2*sum(sin(pi*p22[,i]/L)^2))
}
arrows(y0=disp, y1=disp, x0=c(1:6)-0.2, x1=c(1:6)+0.2, length=0, col="blue")
arrows(y0=disp2, y1=disp2, x0=c(1:5)-0.2, x1=c(1:5)+0.2, length=0, col="blue")
legend("bottomright", legend=c("$E_{\\pi}(\\vec p_1)+E_{\\pi}(\\vec p_2)$"), bty="n", lty=c(1), col=c("blue"))
tikz.finalize(tikzfiles=tikzfiles,clean=TRUE)
}
cut <- 3.
load(paste("fit-detEq-result-cut", cut, ".Rdata", sep=""))
qsq <- numeric()
## determine the rough q^2 values for applying a cut
for(i in c(1:length(Epipi$E))) {
dvec=c(0,0,0)
if(Epipi$frame[i] == 2) dvec=c(0,0,1)
if(Epipi$frame[i] == 3) dvec=c(0,1,1)
if(Epipi$frame[i] == 4) dvec=c(1,1,1)
qtsq <- compute.qtildesq(Epipi$E[i], dvec=dvec, L=Epipi$L[i], mpi=Epi$E[i])
qsq[i] <- qtsq$q^2
}
Epipi$qsq <- qsq
Epi$qsq <- qsq
k <- which(Epi$qsq/Epi$E^2 <= cut)
Epipicut <- Epipi[k,]
Epicut <- Epi[k,]
Epipicut <- Epipicut[iiEdata,]
Epicut <- Epicut[iiEdata,]
for(L in Lvalues) {
i <- which(Lvalues == L)
tikzfiles <- tikz.init(paste("energylevels-pfit", L, "-cut", cut, sep=""),width=5,height=6)
op <- options(cex=1.3)
ii <- which(Epipi$L == L)
mirrep <- max(Epipi[ii,]$irrep)
mlevel <- max(Epipi[ii,]$level)
mframe <- max(Epipi[ii,]$frame)
N <- length(Epipi[Epipi$L==L & Epipi$level==1,]$E)
xv <- (Epipi[ii,]$irrep-1)*mframe + Epipi[ii,]$frame
xv[xv==9] <- 6
E <- Epipi[Epipi$L==L,]$E
dE <- Epipi[Epipi$L==L,]$dE
Mpi <- Epi[Epipi$L==L,]$E[1]
a <- 0.2
plot(NA, type="n", ylab="$aE$", xlab="", ylim=c(a,ulim[i]), xlim=c(0.5,N+0.5), xaxt="n")
rect(xleft=xv-0.2, xright=xv, ytop=E+dE, ybottom=E-dE, col="gray")
legend("topright", legend=c("data $E_{\\pi\\pi}$", "fitted $E_{\\pi\\pi}$"), bty="n", fill=c("gray", "red"), border=c("black", "red"))
legend("bottomleft", legend=c(paste("$q^2/M_\\pi^2\\leq$", cut, sep="")), bty="n", cex=1.3)
labels <- c("A1, $\\vec d=[0,0,0]$", "A1, $\\vec d=[0,0,1]$", "A1, $\\vec d=[0,1,1]$", "A1, $\\vec d=[1,1,1]$", "E, $\\vec d=[0,0,0]$", "T2, $\\vec d=[0,0,0]$")
axis(1, at=c(1:N), labels=rep("", times=N))
if(L==32) a <- a-0.03/(ulim[i]-a)
if(L==24) a <- a-0.04/(ulim[i]-a)
if(L==20) a <- a-0.07/(ulim[i]-a)
text(x=c(1:N), y=a, labels=labels, srt = 90, adj = 1, xpd = TRUE,cex=0.8,srt=60)
legend("topleft", legend=c(paste("A40 ensemble, L =", L)), bty="n", cex=1.2)
p1 <- array(c(0,0,0, 0,0,1, 0,1,1, 1,1,1, 0,0,1, 0,1,1), dim=c(3,6))
p2 <- array(c(0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,1,1), dim=c(3,6))
p12 <- array(c(0,0,1, 1,0,0, 0,1,0, 1,0,0, 0,1,1), dim=c(3,5))
p22 <- array(c(0,0,1, -1,0,1, 0,0,1, 0,1,1, 0,1,1), dim=c(3,5))
disp <- c()
disp2 <- c()
for(i in c(1:6)) {
disp[i] <- acosh( cosh(Mpi) + 2*sum(sin(pi*p1[,i]/L)^2)) + acosh( cosh(Mpi) + 2*sum(sin(pi*p2[,i]/L)^2))
}
for(i in c(1:5)) {
disp2[i] <- acosh( cosh(Mpi) + 2*sum(sin(pi*p12[,i]/L)^2)) + acosh( cosh(Mpi) + 2*sum(sin(pi*p22[,i]/L)^2))
}
arrows(y0=disp, y1=disp, x0=c(1:6)-0.2, x1=c(1:6)+0.2, length=0, col="blue")
arrows(y0=disp2, y1=disp2, x0=c(1:5)-0.2, x1=c(1:5)+0.2, length=0, col="blue")
legend("bottomright", legend=c("$E_{\\pi}(\\vec p_1)+E_{\\pi}(\\vec p_2)$"), bty="n", lty=c(1), col=c("blue"))
## add the fitted values to the plot
ii <- which(Epipicut$L == L)
mirrep <- max(Epipicut[ii,]$irrep)
mlevel <- max(Epipicut[ii,]$level)
mframe <- max(Epipicut[ii,]$frame)
N <- length(Epipicut[Epipicut$L==L & Epipicut$level==1,]$E)
xv <- (Epipicut[ii,]$irrep-1)*mframe + Epipicut[ii,]$frame
xv[xv==9] <- 6
E <- bootevalues[1,ii]
dE <- apply(bootevalues[c(1:10),ii], 2, sd)
rect(xleft=xv, xright=xv+0.2, ytop=E+dE, ybottom=E-dE, col="red", border="red")
tikz.finalize(tikzfiles=tikzfiles,clean=TRUE, crop=FALSE)
options(op)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.