library(ggplot2)library(Rose)count=1countfree=1df<-data.frame("momentum"=c(0),"E2"=c(0),"err"=c(0),"corr"=c(0),"delta"=c(0),"deltaerr"=c(0), "k"=c(0),"kerr"=c(0),"L"=c(0), "kcotd"=c(0),"kcotderr"=c(0), "T"=c(0))dfree<-data.frame("momentum"=c(0),"E2"=c(0),"err"=c(0),"Efreelatt"=c(0), "p"=c(0))TLs<-list( c(32,32) )gg<-ggplot()iTL=1for (TL in TLs){T<-TL[1]L<-TL[2]dir<-paste("../out")file=sprintf("%s/G2t_T%d_L%d_msq0-4.900000_msq1-4.650000_l02.500000_l12.500000_mu5.000000_g0.000000_rep0_output",dir,T,L)cat("\n####", file,'\n\n')mt<-read_df(file)all_obs<-get_all_corr(mt)string=sprintf("\\b%s\\b","E1_0")# need to put the delimiters on the word to grep l<-grep(string,all_obs[,"corr"])label<-paste(gsub('\\\\b','',string),"L=",L)n<-all_obs[l,"n"]fit<-get_fit_n(mt,n)m0<-fit[1,1]momenta<-data.frame(c(0,0,0,0,0,0),c(1,0,0,0,0,0),c(1,1,0,0,0,0),c(1,1,1,0,0,0),c(1,0,0,-1,0,0))for (p in momenta){ norm=p[1]^2+p[2]^2+p[3]^2 norm=norm*(2*pi/L)^2 norm1=p[4]^2+p[5]^2+p[6]^2 norm1=norm1*(2*pi/L)^2 norm_tot=(p[1]+p[4])^2+(p[2]+p[5])^2+(p[3]+p[6])^2 norm_tot=norm_tot*(2*pi/L)^2 E2=sqrt(m0^2+norm1)+sqrt(m0^2+norm ) E2cm=E2#sqrt(E2^2-norm_tot) s=sprintf("(%d,%d,%d)(%d,%d,%d)",p[1],p[2],p[3],p[4],p[5],p[6]) hatp2=4*( sin(pi*p[1]/L)^2+sin(pi*p[2]/L)^2+sin(pi*p[3]/L)^2 ) Efl<-acosh(cosh(m0) +0.5*hatp2 ) hatp2=4*(sin(pi*p[4]/L)^2+sin(pi*p[5]/L)^2+sin(pi*p[6]/L)^2 ) Efl<-Efl +acosh(cosh(m0) +0.5*hatp2 ) Efl=Efl#sqrt(Efl^2-norm_tot) dfree[countfree,]<-list(s, E2cm , 0, Efl, norm_tot*4 ) countfree=countfree+1}mom=1for ( s inc("E2_0")){ string=sprintf("\\b%s\\b",s)# need to put the delimiters on the word to grep l<-grep(string,all_obs[,"corr"]) label<-paste(gsub('\\\\b','',string),"L=",L) n<-all_obs[l,"n"] fit<-get_fit_n(mt,n) df[count,]<-list(dfree[1,1],fit[1,1],fit[1,2] ,label,fit[3,5],fit[3,6], fit[3,7],fit[3,8], L, fit[3,9],fit[3,10],T) count=count+1}for ( s inc("E2_0_p1","E2_0_p11","E2_0_p111")){ string=sprintf("\\b%s\\b",s)# need to put the delimiters on the word to grep l<-grep(string,all_obs[,"corr"]) label<-paste(gsub('\\\\b','',string),"L=",L) n<-all_obs[l,"n"] fit<-get_fit_n(mt,n) df[count,]<-list(dfree[mom+1,1],fit[1,1],fit[1,2] ,label,fit[2,5],fit[2,6], fit[2,7],fit[2,8], L, fit[2,9],fit[2,10],T) count=count+1 mom=mom+1}for ( s inc("E2_0_A1")){ string=sprintf("\\b%s\\b",s)# need to put the delimiters on the word to grep l<-grep(string,all_obs[,"corr"]) label<-paste(gsub('\\\\b','',string),"L=",L) n<-all_obs[l,"n"] fit<-get_fit_n(mt,n) df[count,]<-list(dfree[5,1],fit[1,1],fit[1,2] ,label,fit[2,5],fit[2,6], fit[2,7],fit[2,8],L, fit[2,9],fit[2,10],T) count=count+1}iTL=iTL+1}lab1<-"$E_2-E_2^{free,latt}$"lab<-"$E_2-E_2^{free,cont}$"n<-"a"gg<-gg+scale_color_manual(labels=c(lab,lab1),values=c("red","blue"))gg<-gg+scale_shape_manual(labels=c(lab,lab1),values=c(1,2))gg<-gg+scale_fill_manual(labels=c(lab,lab1),values=c("red","blue")) gg<-gg+geom_point(data=df, mapping=aes(x=df$momentum, y=df$E2-dfree$E2,color=n, shape=n) ) gg<-gg+geom_errorbar(data=df, mapping=aes(x=df$momentum, ymin=df$E2-df$err-dfree$E2,ymax=df$E2+df$err-dfree$E2,color=n, shape=n ) , width =0.2 )###########plot first gg<-gg+geom_abline(slope=0, itercept=0, size =1. )gg <-gg+ggplot2::theme_bw()+theme(axis.text.x =element_text(angle =45, vjust =1, hjust=1))+theme(legend.title =element_blank(),legend.position =c(0.73,0.81))gg<-gg+ggplot2::ylab("$\\Delta E_2$")+ggplot2::xlab("$p_1,p_2$")gg<-gg +coord_cartesian(ylim=c(-0.001,0.006))fig<-ggsave_pdf = "Delta_E"texfile=paste0(save_pdf,".tex")tikzDevice::tikz(texfile,standAlone=TRUE, width =5/1.618033, height =5/1.618033)plot(fig)dev.off()tools::texi2dvi(paste0(save_pdf,".tex"),pdf=TRUE)# ############## E_latt n1<-"b"gg<-gg+geom_point(data=df, mapping=aes(x=df$momentum, y=df$E2-dfree$Efreelatt, color=n1, shape=n1) , width =0.2 ) gg<-gg+geom_errorbar(data=df, mapping=aes(x=df$momentum, ymin=df$E2-df$err-dfree$Efreelatt,ymax=df$E2+df$err-dfree$Efreelatt, color=n1,shape=n1) , width =0.2 )fig<-ggsave_pdf = "Delta_E_latt"texfile=paste0(save_pdf,".tex")tikzDevice::tikz(texfile,standAlone=TRUE, width =5/1.618033, height =5/1.618033)plot(fig)dev.off()tools::texi2dvi(paste0(save_pdf,".tex"),pdf=TRUE)