first analysis added

gevp tools/analysis.jl

+using juobs, ADerrors, DelimitedFiles, PyPlot, LaTeXStrings, Revise, LsqFit, DataFrames, CSV
+include("/Users/ale/juobs/gevp tools/tools.jl")
+include("/Users/ale/juobs/constants/juobs_const.jl")
+##
+const path_plot = "/Users/ale/Google Drive/phd/analysis/gevp analysis/plot"
+const path_result = "/Users/ale/Google Drive/phd/analysis/gevp analysis"
+const path_data = "/Users/ale/Desktop/data"
+##
+const ensembles = ["H400", "N200", "N203", "N300", "J303"]
+const deg = [true, false, false, true, false]
+
+const L = [32, 48, 48, 48, 64]
+const beta = [3.46 , 3.55, 3.55, 3.70, 3.70]
+const R = ["H400r001", ["N200r000", "N200r001"], ["N203r000", "N203r001"], "N300r002", "J303r003"]
+
+
+
+mu_pp = Vector{Vector{Vector{Float64}}}(undef, length(ensembles))
+mu_aa = Vector{Vector{Vector{Float64}}}(undef, length(ensembles))
+
+m_lh = Vector{Vector{uwreal}}(undef, length(ensembles))
+
+
+#### choose the ensembles and read data 
+iens =5
+pp = read_dat(ensembles[iens], "G5","G5")
+a1a1 = read_dat(ensembles[iens], "G1G5","G1G5")    
+a2a2 = read_dat(ensembles[iens], "G2G5","G2G5")    
+a3a3 = read_dat(ensembles[iens], "G3G5","G3G5")   
+pa1 = read_dat(ensembles[iens], "G5","G1G5")    
+pa2 = read_dat(ensembles[iens], "G5","G2G5")    
+pa3 = read_dat(ensembles[iens], "G5","G3G5")    
+    
+pp_obs = corr_obs.(pp, L=L[iens])
+a1a1_obs = corr_obs.(a1a1, L=L[iens])
+a2a2_obs = corr_obs.(a2a2, L=L[iens])
+a3a3_obs = corr_obs.(a3a3, L=L[iens])
+pa1_obs = corr_obs.(pa1, L=L[iens])
+pa2_obs = corr_obs.(pa2, L=L[iens])
+pa3_obs = corr_obs.(pa3, L=L[iens])
+
+mu_pp[iens] = getfield.(pp_obs, :mu)
+##############################
+## GEVP ANALYSIS #############
+##############################
+
+
+###################
+## HH SECTOR
+####################
+
+#Initialize result container 
+
+
+##initialize partial  results 
+mul, mus, muh = get_mu(mu_pp[iens], deg[iens])
+
+m_lh = Vector{uwreal}(undef, length(muh))
+m_sh = Vector{uwreal}(undef, length(muh))
+m_lh_star = Vector{uwreal}(undef, length(muh))
+m_sh_star = Vector{uwreal}(undef, length(muh))
+m_hh = Vector{uwreal}(undef, length(muh))
+m_hh_star = Vector{uwreal}(undef, length(muh))
+f_lh = Vector{uwreal}(undef, length(muh))
+f_sh = Vector{uwreal}(undef, length(muh))
+f_lh_star = Vector{uwreal}(undef, length(muh))
+f_sh_star = Vector{uwreal}(undef, length(muh))
+f_hh = Vector{uwreal}(undef, length(muh))
+f_hh_star = Vector{uwreal}(undef, length(muh))
+##
+hh11 = get_hh(mu_pp[iens], [pp_obs], deg[iens])
+hh12 = get_hh(mu_pp[iens], [pa1_obs, pa2_obs, pa3_obs], deg[iens])
+hh22 = get_hh(mu_pp[iens], [a1a1_obs, a2a2_obs, a3a3_obs], deg[iens])
+gevp_hhmat = comp_mat(hh11, hh12, hh22)
+@time hh_en, hh_evec = comp_energy(hh11,hh12,hh22, true)
+
+
+
+##  HH fitting procedure 
+mu_hh = 3
+mu_hs = 3
+mu_hl = 3
+@. model(x,p)=p[1] + p[2] * exp(-(p[3]-p[1])*x)
+ydata1 = hh_en[mu_hh].en_val[1][15:end-15] 
+ydata2 = hh_en[mu_hh].en_val[2][13:end-10] 
+xdata1 = collect(1:length(ydata1))
+xdata2 = collect(1:length(ydata2))
+wpm = Dict{Int64, Vector{Float64}}()
+wpm[303] = [-1.0, 2.0,-1.0, -1.0]
+##
+hhpar1 = fit_routine(model, ydata1, wpm=wpm)
+##
+hhpar2 = fit_routine(model, ydata2, wpm=wpm)
+##
+fig = figure()
+plt.errorbar(xdata1, value.(ydata1), yerr = err.(ydata1), ms=2, fmt="^", mfc="none", mec="cornflowerblue",ecolor="cornflowerblue", capsize=2, label=L"$am_{\eta_c}$") #am_{D_s}
+plt.plot(xdata1, model(xdata1, value.(hhpar1)), lw=0.5, c="midnightblue" )
+plt.errorbar(xdata2, value.(ydata2), yerr = err.(ydata2), ms=2, fmt="^", mfc="none", mec="tomato", ecolor="tomato" ,capsize=2, label=L"$am_{J/\psi}$") #am_{D_s^*}
+plt.plot(xdata2, model(xdata2, value.(hhpar2)), lw=0.5, c="darkred" )
+plt.xlim(1,50)
+#plt.ylim(0.70,0.80)
+plt.xlabel(L"$t$")
+plt.ylabel(L"$E(t, t_0=5)$")
+plt.title(string(ensembles[iens]))#,"  ", L"$\mu_c = $", fen[mu_sec].mu) )
+plt.axvline(x=5, lw = 1, c ="green", ls="dashed", label=L"$t_0$")
+plt.legend()
+name = string(ensembles[iens], "_hh_gevp.pdf")
+#plt.savefig(joinpath(path_plot, name), dpi=500)
+display(fig)
+
+## temporary container for hh masses
+hhps_mass = hhpar1[1]
+hhvec_mass = hhpar2[1]
+m_hh[mu_hh] = hhpar1[1]
+m_hh_star[mu_hh] = hhpar2[1]
+
+## pseudo decay hh
+
+ps_pneff = pseudo_mat_elem(hh_evec[mu_hh], hh11[1][mu_hh].obs, hhps_mass)
+uwerr.(ps_pneff)
+@. dec_mod(x,p) = p[1] + p[2] *x
+ydat = ps_pneff[30:80]
+ps_tofit = fit_routine(dec_mod, ydat, 2, wpm=wpm)
+figu = figure()
+plt.errorbar(collect(1:length(ydat)), value.(ydat), yerr = err.(ydat), fmt="d")
+#plt.ylim(0.26,0.27)
+plt.plot(collect(1:length(ydat)), (dec_mod(collect(1:length(ydat)),value.(ps_tofit)  )))
+display(figu)
+hhps_dec = extract_dec_const(ps_tofit[1], hhps_mass, hh_en[mu_hh].mu)
+uwerr(hhps_dec)
+hhps_dec
+f_hh[mu_hh] = hhps_dec
+
+## vec decay hh
+ak_corr = (hh22[1][mu_hh].obs .+ hh22[2][mu_hh].obs .+ hh22[3][mu_hh].obs)
+aux_vec = vec_mat_elem(hh_evec[mu_hh], ak_corr, gevp_hhmat[mu_hh].mat, hhvec_mass)
+uwerr.(aux_vec)
+ydat = aux_vec[30:80]
+vec_tofit = fit_routine(dec_mod, ydat, 2)
+figu = figure()
+plt.errorbar(collect(1:length(ydat)), value.(ydat), yerr = err.(ydat), fmt="d")
+#plt.ylim(0.26,0.27)
+plt.plot(collect(1:length(ydat)), (dec_mod(collect(1:length(ydat)),value.(vec_tofit)  )))
+display(figu)
+hhvec_dec = extract_dec_const(vec_tofit[1], hhvec_mass, hh_en[mu_hh].mu)
+uwerr(hhvec_dec)
+hhvec_dec
+f_hh_star[mu_hh] = hhvec_dec
+
+#################
+## HS SECTOR ####
+#################
+hs11 = get_sh(mu_pp[iens], [pp_obs], deg[iens])
+hs12 = get_sh(mu_pp[iens], [pa1_obs, pa2_obs, pa3_obs], deg[iens])
+hs22 = get_sh(mu_pp[iens], [a1a1_obs, a2a2_obs, a3a3_obs], deg[iens])
+gevp_hsmat = comp_mat(hs11, hs12, hs22)
+hs_en, hs_evec = comp_energy(hs11, hs12, hs22, true)
+
+## HS fitting procedure 
+
+@. model(x,p)=p[1] + p[2] * exp(-(p[3]-p[1])*x)
+ydata1 = hs_en[mu_hs].en_val[1][11:end-10] 
+ydata2 = hs_en[mu_hs].en_val[2][11:end-10] 
+xdata1 = collect(1:length(ydata1))
+xdata2 = collect(1:length(ydata2))
+wpm = Dict{Int64, Vector{Float64}}()
+wpm[303] = [-1.0, 2.0,-1.0, -1.0]
+println(hs_en[mu_hs].mu)
+##
+hspar1 = fit_routine(model, ydata1, wpm=wpm)
+##
+hspar2 = fit_routine(model, ydata2, wpm=wpm)
+##
+fig = figure()
+plt.errorbar(xdata1, value.(ydata1), yerr = err.(ydata1), ms=2, fmt="^", mfc="none", mec="cornflowerblue",ecolor="cornflowerblue", capsize=2, label=L"$am_{\eta_c}$") #am_{D_s}
+plt.plot(xdata1, model(xdata1, value.(hspar1)), lw=0.5, c="midnightblue" )
+plt.errorbar(xdata2, value.(ydata2), yerr = err.(ydata2), ms=2, fmt="^", mfc="none", mec="tomato", ecolor="tomato" ,capsize=2, label=L"$am_{J/\psi}$") #am_{D_s^*}
+plt.plot(xdata2, model(xdata2, value.(hspar2)), lw=0.5, c="darkred" )
+plt.xlim(1,40)
+#plt.ylim(0.40,0.7)
+plt.xlabel(L"$t$")
+plt.ylabel(L"$E(t, t_0=5)$")
+plt.title(string(ensembles[iens]))#,"  ", L"$\mu_c = $", fen[mu_sec].mu) )
+plt.axvline(x=5, lw = 1, c ="green", ls="dashed", label=L"$t_0$")
+plt.legend()
+name = string(ensembles[iens], "_hl_gevp.pdf")
+#plt.savefig(joinpath(path_plot, name), dpi=500)
+display(fig)
+
+##temporary container for hs masses
+hsps_mass = hspar1[1]
+hsvec_mass = hspar2[1]
+m_sh[mu_hs] = hspar1[1]
+m_sh_star[mu_hs] = hspar2[1]
+
+## pseudo decay hs
+ps_pneff = pseudo_mat_elem(hs_evec[mu_hs], hs11[1][mu_hs].obs, hsps_mass)
+uwerr.(ps_pneff)
+@. dec_mod(x,p) = p[1] + p[2] *x
+ydat = ps_pneff[30:60]
+ps_tofit = fit_routine(dec_mod, ydat, 2)
+figu = figure()
+plt.errorbar(collect(1:length(ydat)), value.(ydat), yerr = err.(ydat), fmt="d")
+#plt.ylim(0.26,0.27)
+plt.plot(collect(1:length(ydat)), (dec_mod(collect(1:length(ydat)),value.(ps_tofit)  )))
+display(figu)
+hsps_dec = extract_dec_const(ps_tofit[1], hsps_mass, hs_en[mu_hs].mu)
+uwerr(hsps_dec)
+hsps_dec
+f_sh[mu_hs] = hsps_dec
+
+## vec decay hs
+ak_corr = (hs22[1][mu_hs].obs .+ hs22[2][mu_hs].obs .+ hs22[3][mu_hs].obs)
+aux_vec = vec_mat_elem(hs_evec[mu_hs], ak_corr, gevp_hsmat[mu_hs].mat, hsvec_mass)
+uwerr.(aux_vec)
+ydat = aux_vec[25:50]
+vec_tofit = fit_routine(dec_mod, ydat, 2)
+figu = figure()
+plt.errorbar(collect(1:length(ydat)), value.(ydat), yerr = err.(ydat), fmt="d")
+#plt.ylim(0.26,0.27)
+plt.plot(collect(1:length(ydat)), (dec_mod(collect(1:length(ydat)),value.(vec_tofit)  )))
+display(figu)
+hsvec_dec = extract_dec_const(vec_tofit[1], hsvec_mass, hs_en[mu_hs].mu)
+uwerr(hsvec_dec)
+hsvec_dec
+f_sh_star[mu_hs] = hsvec_dec
+
+
+#################
+## HL SECTOR ####
+#################
+hl11 = get_lh(mu_pp[iens], [pp_obs], deg[iens])
+hl12 = get_lh(mu_pp[iens], [pa1_obs, pa2_obs, pa3_obs], deg[iens])
+hl22 = get_lh(mu_pp[iens], [a1a1_obs, a2a2_obs, a3a3_obs], deg[iens])
+gevp_hlmat = comp_mat(hl11, hl12, hl22)
+hl_en, hl_evec = comp_energy(hl11, hl12, hl22, true)
+
+## HL fitting procedure 
+
+@. model(x,p)=p[1] + p[2] * exp(-(p[3]-p[1])*x)
+ydata1 = hl_en[mu_hl].en_val[1][10:end-10] 
+ydata2 = hl_en[mu_hl].en_val[2][10:end-10] 
+xdata1 = collect(1:length(ydata1))
+xdata2 = collect(1:length(ydata2))
+wpm = Dict{Int64, Vector{Float64}}()
+wpm[303] = [-1.0, 2.0,-1.0, -1.0]
+println(hl_en[mu_hl].mu)
+##
+hlpar1 = fit_routine(model, ydata1, wpm=wpm)
+##
+hlpar2 = fit_routine(model, ydata2, wpm=wpm)
+##
+fig = figure()
+plt.errorbar(xdata1, value.(ydata1), yerr = err.(ydata1), ms=2, fmt="^", mfc="none", mec="cornflowerblue",ecolor="cornflowerblue", capsize=2, label=L"$am_{\eta_c}$") #am_{D_s}
+plt.plot(xdata1, model(xdata1, value.(hlpar1)), lw=0.5, c="midnightblue" )
+plt.errorbar(xdata2, value.(ydata2), yerr = err.(ydata2), ms=2, fmt="^", mfc="none", mec="tomato", ecolor="tomato" ,capsize=2, label=L"$am_{J/\psi}$") #am_{D_s^*}
+plt.plot(xdata2, model(xdata2, value.(hlpar2)), lw=0.5, c="darkred" )
+plt.xlim(1,40)
+plt.ylim(0.40,0.7)
+plt.xlabel(L"$t$")
+plt.ylabel(L"$E(t, t_0=5)$")
+plt.title(string(ensembles[iens]))#,"  ", L"$\mu_c = $", fen[mu_sec].mu) )
+plt.axvline(x=5, lw = 1, c ="green", ls="dashed", label=L"$t_0$")
+plt.legend()
+name = string(ensembles[iens], "_hl_gevp.pdf")
+#plt.savefig(joinpath(path_plot, name), dpi=500)
+display(fig)
+
+##temporary container for hs masses
+hlps_mass = hlpar1[1]
+hlvec_mass = hlpar2[1]
+m_lh[mu_hl] = hlpar1[1]
+m_lh_star[mu_hl] = hlpar2[1]
+
+## pseudo decay hl
+
+ps_pneff = pseudo_mat_elem(hl_evec[mu_hl], hl11[1][mu_hl].obs, hlps_mass)
+uwerr.(ps_pneff)
+@. dec_mod(x,p) = p[1] + p[2] *x
+ydat = ps_pneff[20:50]
+ps_tofit = fit_routine(dec_mod, ydat, 2)
+figu = figure()
+plt.errorbar(collect(1:length(ydat)), value.(ydat), yerr = err.(ydat), fmt="d")
+#plt.ylim(0.26,0.27)
+plt.plot(collect(1:length(ydat)), (dec_mod(collect(1:length(ydat)),value.(ps_tofit)  )))
+display(figu)
+hlps_dec = extract_dec_const(ps_tofit[1], hlps_mass, hl_en[mu_hl].mu)
+uwerr(hlps_dec)
+hlps_dec
+##store hl dec
+f_lh[mu_hl] = hlps_dec
+
+## vec decay hl
+ak_corr = (hl22[1][mu_hl].obs .+ hl22[2][mu_hl].obs .+ hl22[3][mu_hl].obs)
+aux_vec = vec_mat_elem(hl_evec[mu_hl], ak_corr, gevp_hlmat[mu_hl].mat, hlvec_mass)
+uwerr.(aux_vec)
+ydat = aux_vec[20:50]
+vec_tofit = fit_routine(dec_mod, ydat, 2)
+figu = figure()
+plt.errorbar(collect(1:length(ydat)), value.(ydat), yerr = err.(ydat), fmt="d")
+#plt.ylim(0.26,0.27)
+plt.plot(collect(1:length(ydat)), (dec_mod(collect(1:length(ydat)),value.(vec_tofit)  )))
+display(figu)
+hlvec_dec = extract_dec_const(vec_tofit[1], hlvec_mass, hl_en[mu_hl].mu)
+uwerr(hlvec_dec)
+hlvec_dec
+##store hl_star dec
+f_lh_star[mu_hl] = hlvec_dec
+
+
+
+
+#############
+## MATCHING
+#############
+mul, mus, muh = get_mu(mu_pp[iens], deg[iens])
+##
+target = a(beta[iens]) * (2*M[1] + 6*M[2] + M[3] + 3*M[4]) / (12*hc)
+if !deg[iens] 
+    muh_target = match_muc(muh, m_lh, m_sh, m_lh_star, m_sh_star, target)
+else
+    muh_target = match_muc(muh, m_lh, m_lh_star, target)
+end
+uwerr(muh_target)
+muh_target
+
+##interpolate m_lh, m_lh_star, m_sh, m_sh_star, m_hh, m_hh_star
+#m_lh
+par, chi2exp = lin_fit(muh, m_lh)
+m_lh_match = y_lin_fit(par, muh_target)
+uwerr(m_lh_match)
+
+#m_lh_star
+par, chi2exp = lin_fit(muh, m_lh_star)
+m_lh_star_match = y_lin_fit(par, muh_target)
+uwerr(m_lh_star_match)
+
+#m_sh
+par, chi2exp = lin_fit(muh, m_sh)
+m_sh_match = y_lin_fit(par, muh_target)
+uwerr(m_sh_match)
+
+#m_sh_star
+par, chi2exp = lin_fit(muh, m_sh_star)
+m_sh_star_match = y_lin_fit(par, muh_target)
+uwerr(m_sh_star_match)
+
+#m_hh
+par, chi2exp = lin_fit(muh, m_hh)
+m_hh_match = y_lin_fit(par, muh_target)
+uwerr(m_hh_match)
+
+#m_hh_star
+par, chi2exp = lin_fit(muh, m_hh_star)
+m_hh_star_match = y_lin_fit(par, muh_target)
+uwerr(m_hh_star_match)
+
+#f_lh
+par, chi2exp = lin_fit(muh, f_lh)
+f_lh_match = y_lin_fit(par, muh_target)
+uwerr(f_lh_match)
+
+#f_lh_star
+par, chi2exp = lin_fit(muh, f_lh_star)
+f_lh_star_match = y_lin_fit(par, muh_target)
+uwerr(f_lh_star_match)
+
+#f_sh
+par, chi2exp = lin_fit(muh, f_sh)
+f_sh_match = y_lin_fit(par, muh_target)
+uwerr(f_sh_match)
+
+#f_sh_star
+par, chi2exp = lin_fit(muh, f_sh_star)
+f_sh_star_match = y_lin_fit(par, muh_target)
+uwerr(f_sh_star_match)
+
+#f_hh
+par, chi2exp = lin_fit(muh, f_hh)
+f_hh_match = y_lin_fit(par, muh_target)
+uwerr(f_hh_match)
+
+#f_hh_star
+par, chi2exp = lin_fit(muh, f_hh_star)
+f_hh_star_match = y_lin_fit(par, muh_target)
+uwerr(f_hh_star_match)
+
+###################
+##  PLOTS #########
+##################
+
+## m_lh, m_lh_star
+figure()
+title(ensembles[iens])
+xlabel(L"$a\mu_c$")
+ylabel(L"$aM$")
+errorbar(muh, value.(m_lh), err.(m_lh), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
+errorbar(value(muh_target), value(m_lh_match), err(m_lh_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
+
+errorbar(muh, value.(m_lh_star), err.(m_lh_star), fmt="*",ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
+errorbar(value(muh_target), value(m_lh_star_match), err(m_lh_star_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="navy", lw=1)
+legend([L"$m_D$", L"$\langle m_D \rangle$", L"$m_{D^*}$", L"$ \langle m_{D^*}\rangle $"])  
+display(gcf())
+t = string(ensembles[iens], "_mD.pdf")
+savefig(joinpath(path_plot,t))
+close()
+
+
+## m_sh, m_sh_star
+figure()
+title(ensembles[iens])
+xlabel(L"$a\mu_c$")
+ylabel(L"$aM$")
+errorbar(muh, value.(m_sh), err.(m_sh), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
+errorbar(value(muh_target), value(m_sh_match), err(m_sh_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
+
+errorbar(muh, value.(m_sh_star), err.(m_sh_star), fmt="*",ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
+errorbar(value(muh_target), value(m_sh_star_match), err(m_sh_star_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="navy", lw=1)
+legend([L"$m_{D_s}$", L"$\langle m_{D_s} \rangle$", L"$m_{D^*_s}$", L"$ \langle m_{D^*_s}\rangle $"])  
+display(gcf())
+t = string(ensembles[iens], "_mDs.pdf")
+savefig(joinpath(path_plot,t))
+close()
+
+## m_hh, m_hh_star
+figure()
+title(ensembles[iens])
+xlabel(L"$a\mu_c$")
+ylabel(L"$aM$")
+errorbar(muh, value.(m_hh), err.(m_hh), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
+errorbar(value(muh_target), value(m_hh_match), err(m_hh_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
+
+errorbar(muh, value.(m_hh_star), err.(m_hh_star), fmt="*",ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
+errorbar(value(muh_target), value(m_hh_star_match), err(m_hh_star_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="navy", lw=1)
+legend([L"$m_{\eta_c}$", L"$\langle m_{\eta_c} \rangle$", L"$m_{J/\psi}$", L"$ \langle m_{J/\psi}\rangle $"])  
+display(gcf())
+t = string(ensembles[iens], "_m_eta.pdf")
+savefig(joinpath(path_plot,t))
+close()
+
+## f_lh, f_sh
+figure()
+title(ensembles[iens])
+xlabel(L"$a\mu_c$")
+ylabel(L"$af$")
+errorbar(muh, value.(f_lh), err.(f_lh), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
+errorbar(value(muh_target), value(f_lh_match), err(f_lh_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
+
+errorbar(muh, value.(f_sh), err.(f_sh), fmt="*",ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
+errorbar(value(muh_target), value(f_sh_match), err(f_sh_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="navy", lw=1)
+legend([L"$f_{D}$", L"$\langle f_{D} \rangle$", L"$f_{D_s}$", L"$ \langle f_{D_s}\rangle $"])  
+display(gcf())
+t = string(ensembles[iens], "_fD.pdf")
+savefig(joinpath(path_plot,t))
+close()
+
+## f_lh_star, f_sh_star
+figure()
+title(ensembles[iens])
+xlabel(L"$a\mu_c$")
+ylabel(L"$af$")
+errorbar(muh, value.(f_lh_star), err.(f_lh_star), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
+errorbar(value(muh_target), value(f_lh_star_match), err(f_lh_star_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
+
+errorbar(muh, value.(f_sh_star), err.(f_sh_star), fmt="*",ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
+errorbar(value(muh_target), value(f_sh_star_match), err(f_sh_star_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="navy", lw=1)
+legend([L"$f_{D^*}$", L"$\langle f_{D^*} \rangle$", L"$f_{D_s^*}$", L"$ \langle f_{D_s^*}\rangle $"])  
+display(gcf())
+t = string(ensembles[iens], "_fD*.pdf")
+savefig(joinpath(path_plot,t))
+close()
+
+
+## f_hh, f_hh_star
+figure()
+title(ensembles[iens])
+xlabel(L"$a\mu_c$")
+ylabel(L"$af$")
+errorbar(muh, value.(f_hh), err.(f_hh), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
+errorbar(value(muh_target), value(f_hh_match), err(f_hh_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
+
+errorbar(muh, value.(f_hh_star), err.(f_hh_star), fmt="*",ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
+errorbar(value(muh_target), value(f_hh_star_match), err(f_hh_star_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="navy", lw=1)
+legend([L"$f_{\eta_c}$", L"$\langle f_{\eta_c} \rangle$", L"$f_{J/\psi}$", L"$ \langle f_{J/\psi}\rangle $"])  
+display(gcf())
+t = string(ensembles[iens], "_f_eta.pdf")
+savefig(joinpath(path_plot,t))
+close()
+
+###############
+## RESULTS#####
+###############
+
+#Quark mass
+muc = zm_tm(beta[iens]) * muh_target * sqrt(8 * t0(beta[iens]))
+uwerr(muc)
+#Meson masses
+m_D = m_lh_match * sqrt(8 * t0(beta[iens]))
+uwerr(m_D)
+m_Ds = m_sh_match * sqrt(8 * t0(beta[iens]))
+uwerr(m_Ds)
+m_D_star = m_lh_star_match * sqrt(8 * t0(beta[iens]))
+uwerr(m_D_star)
+m_Ds_star = m_sh_star_match * sqrt(8 * t0(beta[iens]))
+uwerr(m_Ds_star)
+m_eta = m_hh_match * sqrt(8 * t0(beta[iens]))
+uwerr(m_eta)
+m_j = m_hh_star_match * sqrt(8 * t0(beta[iens]))
+uwerr(m_j)
+
+#Decay constant 
+f_D = f_lh_match * sqrt(8 * t0(beta[iens]))
+uwerr(f_D)
+f_Ds = f_sh_match * sqrt(8 * t0(beta[iens]))
+uwerr(f_Ds)
+f_D_star = f_lh_star_match * sqrt(8 * t0(beta[iens]))
+uwerr(f_D_star)
+f_Ds_star = f_sh_star_match * sqrt(8 * t0(beta[iens]))
+uwerr(f_Ds_star)
+f_eta = f_hh_match * sqrt(8 * t0(beta[iens]))
+uwerr(f_eta)
+f_j = f_hh_star_match * sqrt(8 * t0(beta[iens]))
+uwerr(f_j)
+
+
+println("         ", ensembles[iens]) 
+println(L"Z^{tm}_M \mu_c \sqrt{8t_0} =  ", muc)
+println( L"M_D \sqrt{8t_0} =  ", m_D)
+println( L"M_Ds \sqrt{8t_0} =  ", m_Ds)
+println( L"M_D* \sqrt{8t_0} =  ", m_D_star)
+println( L"M_Ds* \sqrt{8t_0} =  ", m_Ds_star)
+println( L"M_eta* \sqrt{8t_0} =  ", m_eta)
+println( L"M_jpsi* \sqrt{8t_0} =  ", m_j)
+println( L"f_D \sqrt{8t_0} =  ", f_D)
+println( L"f_Ds \sqrt{8t_0} =  ", f_Ds)
+println( L"f_D* \sqrt{8t_0} =  ", f_D_star)
+println( L"f_Ds* \sqrt{8t_0} =  ", f_Ds_star)
+println( L"f_eta \sqrt{8t_0} =  ", f_eta)
+println( L"f_jpsi \sqrt{8t_0} =  ", f_j)
+
+
+
+
+
+
+## DataFrames with partial mass data 
+
+mass_df = DataFrame(ens=String[], muc=Float64[], mDs= Float64[], mDserr= Float64[], mDst= Float64[], mDsterr= Float64[], meta= Float64[], metaerr= Float64[], mJ= Float64[], mJerr=Float64[] )
+##
+push!(mass_df , ["N203", 0.173375, 0.61040, 35e-5, 0.6567, 10e-4, 0.93127, 14e-5, 0.96969, 25e-5])
+push!(mass_df , ["N203", 0.18225, 0.62854, 34e-5, 0.6732, 10e-4, 0.96287, 14e-5, 1.00036, 24e-5])
+push!(mass_df , ["N203", 0.191625, 0.64604, 41e-5, 0.6895, 10e-4, 0.99403, 13e-5, 1.03073, 22e-5])
+push!(mass_df , ["N200", 0.17285, 0.61362, 39e-5, 0.66079, 65e-5,0.92917 , 14e-5 ,0.96765 ,31e-5])
+push!(mass_df , ["N200", 0.1819, 0.63145, 44e-5, 0.67644, 66e-5, 0.96076, 13e-5 , 0.99831 , 30e-5])
+push!(mass_df , ["N200", 0.190995, 0.64913, 41e-5, 0.69333, 66e-5, 0.99178, 13e-5 ,1.02860  ,29e-5 ])
+
+##
+path_result
+##
+CSV.write(joinpath(path_result,"uncomplete_result.csv"), mass_df)

gevp tools/results.rtf

+{\rtf1\ansi\ansicpg1252\cocoartf2509
+\cocoatextscaling0\cocoaplatform0{\fonttbl\f0\fswiss\fcharset0 Helvetica;}
+{\colortbl;\red255\green255\blue255;}
+{\*\expandedcolortbl;;}
+\paperw11900\paperh16840\margl1440\margr1440\vieww10800\viewh8400\viewkind0
+\pard\tx566\tx1133\tx1700\tx2267\tx2834\tx3401\tx3968\tx4535\tx5102\tx5669\tx6236\tx6803\pardirnatural\partightenfactor0
+
+\f0\fs24 \cf0 $N200\
+$Z^\{tm\}_M \\mu_c \\sqrt\{8t_0\} =  $3.095508084190539 +/- 0.07424875770953374\
+$M_D \\sqrt\{8t_0\} =  $3.958264727738654 +/- 0.05515406625525418\
+$M_Ds \\sqrt\{8t_0\} =  $4.072538504215807 +/- 0.05467443552294891\
+$M_D* \\sqrt\{8t_0\} =  $4.229272703354319 +/- 0.05006705711946871\
+$M_Ds* \\sqrt\{8t_0\} =  $4.361851498402342 +/- 0.04975162489041616\
+$M_eta* \\sqrt\{8t_0\} =  $6.198773176967105 +/- 0.0967669789483958\
+$M_jpsi* \\sqrt\{8t_0\} =  $6.439877922484902 +/- 0.0942314903150228\
+$f_D \\sqrt\{8t_0\} =  $0.47414683350372805 +/- 0.002821170657997782\
+$f_Ds \\sqrt\{8t_0\} =  $0.5204052951135102 +/- 0.002351052892691432\
+$f_D* \\sqrt\{8t_0\} =  $0.18581645952897405 +/- 0.0032030018851700124\
+$f_Ds* \\sqrt\{8t_0\} =  $0.20920280395715848 +/- 0.0024790842850004667\
+$f_eta \\sqrt\{8t_0\} =  $0.9266388913538548 +/- 0.011566795910907063\
+$f_jpsi \\sqrt\{8t_0\} =  $0.42904806660632133 +/- 0.007884932722702343\
+$\
+$N203\
+$Z^\{tm\}_M \\mu_c \\sqrt\{8t_0\} =  $3.07608073397371 +/- 0.07386718423714976\
+$M_D \\sqrt\{8t_0\} =  $3.966409130806326 +/- 0.054959061347465483\
+$M_Ds \\sqrt\{8t_0\} =  $4.030944118480624 +/- 0.05473249525087825\
+$M_D* \\sqrt\{8t_0\} =  $4.254701932623177 +/- 0.04997144163893854\
+$M_Ds* \\sqrt\{8t_0\} =  $4.320267491882463 +/- 0.049563464976048445\
+$M_eta* \\sqrt\{8t_0\} =  $6.175405773876267 +/- 0.09615293679832428\
+$M_jpsi* \\sqrt\{8t_0\} =  $6.418715171251667 +/- 0.09346001316883264\
+$f_D \\sqrt\{8t_0\} =  $0.48810439044765425 +/- 0.0026045102710287034\
+$f_Ds \\sqrt\{8t_0\} =  $0.513802783972585 +/- 0.002279763193762235\
+$f_D* \\sqrt\{8t_0\} =  $0.19370973720440873 +/- 0.0028186032755263607\
+$f_Ds* \\sqrt\{8t_0\} =  $0.20485862913653624 +/- 0.0024571277939203565\
+$f_eta \\sqrt\{8t_0\} =  $0.9262934445841321 +/- 0.011380030137209169\
+$f_jpsi \\sqrt\{8t_0\} =  $0.4287220496155749 +/- 0.007707612363034976\
+$\
+$J303\
+$Z^\{tm\}_M \\mu_c \\sqrt\{8t_0\} =  $3.106805094851814 +/- 0.07428519131808012\
+$M_D \\sqrt\{8t_0\} =  $3.951817203146415 +/- 0.05634686437630071\
+$M_Ds \\sqrt\{8t_0\} =  $4.091122020600407 +/- 0.05584799044824401\
+$M_D* \\sqrt\{8t_0\} =  $4.220571557429657 +/- 0.050571226554394465\
+$M_Ds* \\sqrt\{8t_0\} =  $4.378430678774368 +/- 0.05075096037526913\
+$M_eta* \\sqrt\{8t_0\} =  $6.222653661097548 +/- 0.09980503048093914\
+$M_jpsi* \\sqrt\{8t_0\} =  $6.45871173311589 +/- 0.09688893871726818\
+$f_D \\sqrt\{8t_0\} =  $0.46790893343380957 +/- 0.0035987447382407898\
+$f_Ds \\sqrt\{8t_0\} =  $0.5204463638699105 +/- 0.0025343572218284698\
+$f_D* \\sqrt\{8t_0\} =  $0.17563003352045548 +/- 0.004603348336097794\
+$f_Ds* \\sqrt\{8t_0\} =  $0.20465609887290984 +/- 0.0030029574421112505\
+$f_eta \\sqrt\{8t_0\} =  $0.876239237941407 +/- 0.008803719053134774\
+$f_jpsi \\sqrt\{8t_0\} =  $0.4070411608133558 +/- 0.006650155013434349\
+}
\ No newline at end of file