initial commit

.vscode/settings.json

+{}
\ No newline at end of file

const.jl

+#========= ENSEMBLE DATABASE ========#
+ens_db = Dict(
+    #"ens_id"=>[L, beta, is_deg?, m_pi]
+    "H400" => [32, 3.46, true, 0.16345],
+    "N200" => [48, 3.55, false, 0.09222],
+    "N203" => [48, 3.55, false, 0.11224],
+    "N300" => [48, 3.70, true, 0.10630],
+    "J303" => [64, 3.70, false, 0.06514]
+)
+ttl_obs_db = Dict(
+    "muc"=> ["mu_c"],
+    "ll" => ["m_pi", "m_rho", "f_pi", "f_rho"],
+    "ls" => ["m_k", "m_k_star", "f_k", "f_k_star"],
+    "lh" => ["m_D", "m_D_star", "f_D", "f_D_star"],
+    "ss" => ["m_etaprime", "m_phi", "f_etaprime", "f_phi"],
+    "sh" => ["m_Ds", "m_Ds_star", "f_Ds", "f_Ds_star"],
+    "hh" => ["m_etac", "m_jpsi", "f_etac", "f_jpsi"]
+)
+ylbl_db = Dict(
+    "muc" => [L"$Z^{tm}_M \mu_c \sqrt{8t_0}$"],
+    "ll" => [L"$m_{π} \sqrt{8t_0}$", L"$m_{\rho} \sqrt{8t_0}$", L"$f_{\pi} \sqrt{8t_0}$", L"$f_{\rho} \sqrt{8t_0}$"],
+    "ls" => [L"$m_K \sqrt{8t_0}$", L"$m_{K^*} \sqrt{8t_0}$", L"$f_K \sqrt{8t_0}$", L"$f_{K^*} \sqrt{8t_0}$"],
+    "lh" => [L"$m_D \sqrt{8t_0}$", L"$m_{D^*} \sqrt{8t_0}$", L"$f_D \sqrt{8t_0}$", L"$f_{D^*} \sqrt{8t_0}$"],
+    "ss" => [L"$m_{\eta'} \sqrt{8t_0}$", L"$m_{\phi} \sqrt{8t_0}$", L"$f_{η'} \sqrt{8t_0}$", L"$f_{\phi} \sqrt{8t_0}$"],
+    "sh" => [L"$m_{D_s} \sqrt{8t_0}$", L"$m_{D_s^*} \sqrt{8t_0}$", L"$f_{D_s} \sqrt{8t_0}$", L"$f_{D_s^*} \sqrt{8t_0}$"],
+    "hh" => [L"$m_etac \sqrt{8t_0}$", L"$m_{J/ψ} \sqrt{8t_0}$", L"$f_{η_c} \sqrt{8t_0}$", L"$f_{J/ψ} \sqrt{8t_0}$"]
+)
+
+#PDG
+const hc = 197.3269804 #MeV fm
+const M_values = [1869.65, 2010.26, 1968.34, 2112.2, 2980.3, 3096.916] #MD, MD*, MDs, MDs*, \eta_c, J/\psi (MeV)
+const M_error = [0.05, 0.05, 0.07, 0.4, 1.2, 0.011]
+#1802.05243 
+const b_values = [3.40, 3.46, 3.55, 3.70, 3.85]
+const b_values2 = [3.40, 3.46, 3.55, 3.70]
+const ZM_data = [1.9684, 1.9935, 2.0253, 2.0630, 2.0814]
+const ZM_error = [35, 27, 33, 38, 45] .* 1e-4
+const ZM_tm_data = [2.6047, 2.6181, 2.6312, 2.6339, 2.6127]
+const ZM_tm_error = [42, 33, 42, 48, 55] .* 1e-4
+#1608.08900
+const t0_data = [2.86, 3.659, 5.164, 8.595]
+const t0_error = [11, 16, 18, 29] .* 1e-3
+const t0_ph_value = [0.413]
+const t0_ph_error = ones(1,1) .* 5e-3
+#Covariance matrices
+const C1 = zeros(5, 5)
+const C2 = zeros(5, 5)
+const C3 = zeros(4, 4)
+const C4 = zeros(6,6)
+for i = 1:6
+    C4[i,i] = M_error[i] ^ 2
+    if i<= 5
+        C1[i, i] = ZM_error[i] ^ 2
+        C2[i, i] = ZM_tm_error[i] ^ 2
+        if i<=4
+            C3[i,i] = t0_error[i] ^ 2
+        end
+    end
+end
+
+const ZM = cobs(ZM_data, C1, "ZM values")
+const ZM_tm = cobs(ZM_tm_data, C2, "ZM_tm values")
+const M = cobs(M_values, C4, "charmed meson masses")
+const t0_ph = cobs(t0_ph_value, t0_ph_error .^ 2, "sqrt(8 t0) (fm)")
+const t0_ = cobs(t0_data, C3, "t0")
+const a_ = t0_ph ./ sqrt.(8 .* t0_)
+
+zm(beta::Float64) = ZM[b_values .== beta][1]
+zm_tm(beta::Float64) = ZM_tm[b_values .== beta][1]
+t0(beta::Float64) = t0_[b_values2 .== beta][1]
+a(beta::Float64) = a_[b_values2 .== beta][1]

gevp.jl

+using Revise, juobs, BDIO, ADerrors, DelimitedFiles, PyPlot, LaTeXStrings, LsqFit
+
+
+#============= SET UP VARIABLES ===========#
+const path_data = "/Users/ale/Desktop/data"
+const path_plat = "/Users/ale/automation/plat.txt"
+const path_results = "/Users/ale/Desktop/results_gevp"
+const ensembles = [ "J303", "N300", "H400"]
+const sector = Dict("ll"=>false, "ls"=>false, "lh"=>true, "ss"=>false, "sh"=>true, "hh"=>true )
+const rwf = false
+const compute_t0 = false
+const mass_shift = false #not implemented yet
+#go to line 443 to change CL+chiral extrapolation fit model
+
+#========= DO NOT CHANGE AFTER THIS LINE ==========#
+include("types.jl")
+include("tools.jl")
+include("const.jl")
+
+#======== GET ENSEMBLE INFORMATION FROM DATABASE ==========#
+ensinfo = Vector{EnsInfo}(undef, length(ensembles))
+for i in 1:length(ensembles)
+    ens = ensembles[i]
+    try
+        ensinfo[i]= EnsInfo(ens, ens_db[ens])
+    catch
+        error("The ensemble id ", ens, " was not found in the const.jl ens_db database.
+        Please check the ensemble id or update the database")
+    end
+end
+
+##
+#======== ANALYSIS ==========#
+gen_mulist = Vector{Vector{Vector{Float64}}}(undef, length(ensembles))
+println("\n Building matrices \n")
+@time begin
+    matinfo = Vector{Vector{MatInfo}}(undef,length(ensembles))
+    matinfo_vec = Vector{Vector{MatInfo}}(undef,length(ensembles))
+    for i in 1:length(ensinfo)
+        @time begin
+            ens = ensinfo[i]
+            println("\n Ensemble: ", ens.id)
+            pp = get_corr(ens,"G5", "G5", rw=rwf)
+            gen_mulist[i] = getfield.(pp, :mu)
+            a1a1 = get_corr(ens, "G1G5", "G1G5", rw=rwf)
+            matinfo[i] = comp_mat(ens, pp, 4)
+            matinfo_vec[i] = comp_mat(ens, a1a1, 4)
+            println("Time:")
+        end
+    end
+    println("Total time:")
+end
+
+##
+wpmm = Dict{Int64, Vector{Float64}}()
+wpmm[303] = [-1.0, 1.5,-1.0, -1.0]
+wpmm[300] = [-1.0, 1.5,-1.0, -1.0]
+println("\n Computing observables \n")
+@time begin
+    ensobs = Vector{EnsObs}(undef, length(ensembles))
+    for i in 1:length(ensinfo)
+        @time begin
+            ens = ensinfo[i]
+            println("\n Ensemble: ", ens.id)
+            en0_ps = gevp_mass(matinfo[i], t0=2, wpm=wpmm)
+            en0_vec = gevp_mass(matinfo_vec[i], t0=2, wpm=wpmm)
+            println( "\nground  ps state: ", en0_ps)
+            println( "\nground  vec state: ", en0_vec)
+            ensobs[i] = EnsObs(ens, gen_mulist[i], en0_ps, en0_vec)
+            if ens.deg
+                ensobs[i].m_ls = ensobs[i].m_ll
+                ensobs[i].m_ls_vec = ensobs[i].m_ll_vec
+                ensobs[i].m_ss = ensobs[i].m_ll
+                ensobs[i].m_ss_vec = ensobs[i].m_ll_vec
+                ensobs[i].m_sh = ensobs[i].m_lh
+                ensobs[i].m_sh_vec = ensobs[i].m_lh_vec
+            end
+            println("Time:")
+        end
+    end
+    println("Total time:")
+end
+## ##################################### SO FAR SO GOOD
+
+println("\n Matching Procedure \n")
+
+@time begin
+    for i in 1:length(ensobs)
+        ens = ensobs[i]
+        mul, mus, muh = get_mu(ens.mu_list, ens.ensinfo.deg)
+        if compute_t0
+            est_t0 = get_t0(ens.ensinfo, rw=rwf)
+            ens.t0 = est_t0
+            est_a = t0_ph[1] / sqrt(8 * ens.t0)
+            ens.a = est_a
+        else
+            ens.t0 = t0(ens.ensinfo.beta)
+            ens.a = a(ens.ensinfo.beta)
+        end
+        target = ens.a * (2/3*M[1] + 1/3* M[3])/ hc
+        if !ens.ensinfo.deg
+            ens.muh_target = match_muc(muh, ens.m_lh, ens.m_sh, target)
+        else
+            ens.muh_target = match_muc(muh, ens.m_lh, ens.m_sh, target)
+        end
+        uwerr(ens.muh_target)
+        println(t0, "\n", a, "\n", ens.muh_target)
+        
+        # interpolate lh
+        if sector["lh"]
+            # m_lh
+            par, chi2exp = lin_fit(muh, ens.m_lh)
+            ens.m_lh_match = y_lin_fit(par, value(ens.muh_target))
+            uwerr(ens.m_lh_match)
+            #m_lh_star
+            par, chi2exp = lin_fit(muh, ens.m_lh_vec)
+            ens.m_lh_vec_match = y_lin_fit(par, value(ens.muh_target))
+            uwerr(ens.m_lh_vec_match)
+            #=
+            # f_lh
+            par, chi2exp = lin_fit(muh, ens.f_lh)
+            ens.f_lh_match = y_lin_fit(par, value(ens.muh_target))
+            uwerr(ens.f_lh_match)
+            #f_lh_star
+            par, chi2exp = lin_fit(muh, ens.f_lh_vec)
+            ens.f_lh_vec_match = y_lin_fit(par, value(ens.muh_target))
+            uwerr(ens.f_lh_vec_match)
+            =#
+        end
+        # interpolate sh
+        if sector["sh"]
+            if ens.ensinfo.deg
+                ens.m_sh_match = ens.m_lh_match
+                ens.m_sh_vec_match = ens.m_lh_vec_match
+                #ens.f_sh_match = ens.f_lh_match
+                #ens.f_sh_vec_match = ens.f_lh_vec_match
+            else
+                # m_sh
+                par, chi2exp = lin_fit(muh, ens.m_sh)     
+                ens.m_sh_match = y_lin_fit(par, value(ens.muh_target))
+                uwerr(ens.m_sh_match)
+                # m_sh_star
+                par, chi2exp = lin_fit(muh, ens.m_sh_vec)     
+                ens.m_sh_vec_match = y_lin_fit(par, value(ens.muh_target))
+                uwerr(ens.m_sh_vec_match)
+                #=
+                # f_sh
+                par, chi2exp = lin_fit(muh, ens.f_sh)     
+                ens.f_sh_match = y_lin_fit(par, value(ens.muh_target))
+                uwerr(ens.f_sh_match)
+                # f_sh_star
+                par, chi2exp = lin_fit(muh, ens.f_sh_vec)     
+                ens.f_sh_vec_match = y_lin_fit(par, value(ens.muh_target))
+                uwerr(ens.f_sh_vec_match)
+                =#
+            end
+        end
+        # interpolate hh
+        if sector["hh"]
+            # m_hh
+            par, chi2exp = lin_fit(muh, ens.m_hh)
+            ens.m_hh_match = y_lin_fit(par, value(ens.muh_target))
+            uwerr(ens.m_hh_match)
+            # m_hh_vec
+            par, chi2exp = lin_fit(muh, ens.m_hh_vec)
+            ens.m_hh_vec_match = y_lin_fit(par, value(ens.muh_target))
+            uwerr(ens.m_hh_vec_match)
+            #=
+            # f_hh
+            par, chi2exp = lin_fit(muh, ens.f_hh)
+            ens.f_hh_match = y_lin_fit(par, value(ens.muh_target))
+            uwerr(ens.f_hh_match)
+            # f_hh_vec
+            par, chi2exp = lin_fit(muh, ens.f_hh_vec)
+            ens.f_hh_vec_match = y_lin_fit(par, value(ens.muh_target))
+            uwerr(ens.f_hh_vec_match)
+            =#
+        end
+    end
+    println(" Total time:")
+end
+    
+    
+#========== MATCHING PLOTS ==========#
+
+println("\n Matchin Plots \n ")
+check_folder = filter(x-> occursin("plots", x), readdir(path_results))
+if !("plots" in check_folder)
+    mkdir(joinpath(path_results,"plots"))
+end
+path_plot = joinpath(path_results,"plots")
+for i in 1:length(ensobs)
+    ens = ensobs[i]
+    mul, mus, muh = get_mu(ens.mu_list, ens.ensinfo.deg)
+    if sector["lh"]
+        # masses
+        figure()
+        title(ens.ensinfo.id)
+        xlabel(L"$a\mu$")
+        ylabel(L"$aM$")
+        #m_lh
+        errorbar(muh, value.(ens.m_lh), yerr=err.(ens.m_lh), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.m_lh_match), yerr=err(ens.m_lh_match), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
+        #m_lh_star
+        errorbar(muh, value.(ens.m_lh_vec), yerr=err.(ens.m_lh_vec), fmt="*", ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.m_lh_vec_match), yerr=err(ens.m_lh_vec_match), 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())
+        name = string(ens.ensinfo.id,"_match_mD_mDvec.pdf")
+        savefig(joinpath(path_plot, name))
+        close()
+        #=
+        # decays
+        figure()
+        title(ens.ensinfo.id)
+        xlabel(L"$a\mu$")
+        ylabel(L"$af$")
+        #m_lh
+        errorbar(muh, value.(ens.f_lh), yerr=err.(ens.f_lh), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.f_lh_match), yerr=err(ens.f_lh_match), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
+        #m_lh_star
+        errorbar(muh, value.(ens.f_lh_vec), yerr=err.(ens.f_lh_vec), fmt="*", ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.f_lh_vec_match), yerr=err(ens.f_lh_vec_match), fmt="*", ms=5, mfc="none", capsize=2, c="navy", lw=1)
+        legend([L"$f_D$", L"$\langle f_D \rangle$", L"$f_{D^*}$", L"$ \langle f_{D^*}\rangle $"])  
+        display(gcf())
+        name = string(ens.ensinfo.id,"_match_fD_fDvec.pdf")
+        savefig(joinpath(path_plot, name))
+        close()
+        =#
+    end
+    if sector["sh"] && !ens.ensinfo.deg
+        # masses
+        figure()
+        title(ens.ensinfo.id)
+        xlabel(L"$a\mu$")
+        ylabel(L"$aM$")
+        #m_sh
+        errorbar(muh, value.(ens.m_sh), yerr=err.(ens.m_sh), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.m_sh_match), yerr=err(ens.m_sh_match), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
+        #m_sh_star
+        errorbar(muh, value.(ens.m_sh_vec), yerr=err.(ens.m_sh_vec), fmt="*", ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.m_sh_vec_match), yerr=err(ens.m_sh_vec_match), 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())
+        name = string(ens.ensinfo.id,"_match_mDs_mDsvec.pdf")
+        savefig(joinpath(path_plot, name))
+        close()
+        #= 
+        # decays
+        figure()
+        title(ens.ensinfo.id)
+        xlabel(L"$a\mu$")
+        ylabel(L"$af$")
+        #f_sh
+        errorbar(muh, value.(ens.f_sh), yerr=err.(ens.f_sh), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.f_sh_match), yerr=err(ens.f_sh_match), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
+        #fsh_star
+        errorbar(muh, value.(ens.f_sh_vec), yerr=err.(ens.f_sh_vec), fmt="*", ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.f_sh_vec_match), yerr=err(ens.f_sh_vec_match), fmt="*", ms=5, mfc="none", capsize=2, c="navy", lw=1)
+        legend([L"$f_{D_s}$", L"$\langle f_{D_s} \rangle$", L"$f_{D^*_s}$", L"$ \langle f_{D^*_s}\rangle $"])  
+        display(gcf())
+        name = string(ens.ensinfo.id,"_match_fDs_fDsvec.pdf")
+        savefig(joinpath(path_plot, name))
+        close()
+        =#
+    end
+    if sector["hh"] 
+        # masses
+        figure()
+        title(ens.ensinfo.id)
+        xlabel(L"$a\mu$")
+        ylabel(L"$aM$")
+        #m_hh
+        errorbar(muh, value.(ens.m_hh), yerr=err.(ens.m_hh), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.m_hh_match), yerr=err(ens.m_hh_match), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
+        #m_hh_star
+        errorbar(muh, value.(ens.m_hh_vec), yerr=err.(ens.m_hh_vec), fmt="*", ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.m_hh_vec_match), yerr=err(ens.m_hh_vec_match), 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())
+        name = string(ens.ensinfo.id,"_match_metac_mJpsi.pdf")
+        savefig(joinpath(path_plot, name))
+        close()
+        
+        #=
+        # decays
+        figure()
+        title(ens.ensinfo.id)
+        xlabel(L"$a\mu$")
+        ylabel(L"$af$")
+        #f_hh
+        errorbar(muh, value.(ens.f_hh), yerr=err.(ens.f_hh), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.f_hh_match), yerr=err(ens.f_hh_match), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
+        #f_hh_star
+        errorbar(muh, value.(ens.f_hh_vec), yerr=err.(ens.f_hh_vec), fmt="*", ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.f_hh_vec_match), yerr=err(ens.f_hh_vec_match), 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())
+        name = string(ens.ensinfo.id,"_match_fetac_fJpsi.pdf")
+        savefig(joinpath(path_plot, name))
+        close()
+        =#
+    end
+end
+
+
+
+## #####################################
+
+
+#======= RESULTS ==========#
+# charm quark mass 
+muc = zm_tm.(getfield.(ensinfo,:beta)) .* getfield.(ensobs,:muh_target) .* sqrt.(8 * getfield.(ensobs,:t0))
+uwerr.(muc)
+
+# ll
+m_pi = Vector{uwreal}(undef, length(ensembles))
+m_rho =Vector{uwreal}(undef, length(ensembles))
+#=
+f_pi = Vector{uwreal}(undef, length(ensembles))
+f_rho =Vector{uwreal}(undef, length(ensembles))
+=#
+if sector["ll"]
+    try 
+        m_pi .= getfield.(ensobs, :m_ll) .* sqrt.(8 * getfield.(ensobs,:t0))
+        m_rho .= getfield.(ensobs, :m_ll_vec) .* sqrt.(8 * getfield.(ensobs,:t0))
+        #f_pi .= getfield.(ensobs, :f_ll) .* sqrt.(8 * getfield.(ensobs,:t0))
+        #f_rho .= getfield.(ensobs, :f_ll_vec) .* sqrt.(8 * getfield.(ensobs,:t0))        
+        uwerr.(m_pi)
+        uwerr.(m_rho)
+        #uwerr.(f_pi)
+        #uwerr.(f_rho)
+    catch
+        error("Check ensobs[i].mu_list to see whether all the ensembles have a light-light sector.")
+    end
+end
+# ls
+m_k = Vector{uwreal}(undef, length(ensembles))
+m_k_star =Vector{uwreal}(undef, length(ensembles))
+#f_k= Vector{uwreal}(undef, length(ensembles))
+#f_k_star =Vector{uwreal}(undef, length(ensembles))
+if sector["ls"] 
+    try 
+        m_k .= getfield.(ensobs, :m_ls) .* sqrt.(8 * getfield.(ensobs,:t0))
+        m_k_star .= getfield.(ensobs, :m_ls_vec) .* sqrt.(8 * getfield.(ensobs,:t0))
+        #f_k .= getfield.(ensobs, :f_ls) .* sqrt.(8 * getfield.(ensobs,:t0))
+        #f_k_star .= getfield.(ensobs, :f_ls_vec) .* sqrt.(8 * getfield.(ensobs,:t0))        
+        uwerr.(m_k)
+        uwerr.(m_k_star)
+        #uwerr.(f_k)
+        #uwerr.(f_k_star)
+    catch
+        error("Check ensobs[i].mu_list to see whether all the ensembles have a light-stange sector.")
+    end
+end
+# lh
+m_D = Vector{uwreal}(undef, length(ensembles))
+m_D_star =Vector{uwreal}(undef, length(ensembles))
+#f_D= Vector{uwreal}(undef, length(ensembles))
+#f_D_star =Vector{uwreal}(undef, length(ensembles))
+if sector["lh"] 
+    try 
+        m_D .= getfield.(ensobs, :m_lh_match) .* sqrt.(8 * getfield.(ensobs,:t0))
+        m_D_star .= getfield.(ensobs, :m_lh_vec_match) .* sqrt.(8 * getfield.(ensobs,:t0))
+        #f_D .= getfield.(ensobs, :f_lh_match) .* sqrt.(8 * getfield.(ensobs,:t0))
+        #f_D_star .= getfield.(ensobs, :f_lh_vec_match) .* sqrt.(8 * getfield.(ensobs,:t0))        
+        uwerr.(m_D)
+        uwerr.(m_D_star)
+        #uwerr.(f_D)
+        #uwerr.(f_D_star)
+catch
+        error("Check ensobs[i].mu_list to see whether all the ensembles have a light-heavy sector.")
+    end
+end
+# ss
+m_etaprime = Vector{uwreal}(undef, length(ensembles))
+m_phi =Vector{uwreal}(undef, length(ensembles))
+#f_etaprime= Vector{uwreal}(undef, length(ensembles))
+#f_phi =Vector{uwreal}(undef, length(ensembles))
+if sector["ss"] 
+    try 
+        m_etaprime .= getfield.(ensobs, :m_ss) .* sqrt.(8 * getfield.(ensobs,:t0))
+        m_phi .= getfield.(ensobs, :m_ss_vec) .* sqrt.(8 * getfield.(ensobs,:t0))
+        #f_etaprime .= getfield.(ensobs, :f_ss) .* sqrt.(8 * getfield.(ensobs,:t0))
+        #f_phi .= getfield.(ensobs, :f_ss_vec) .* sqrt.(8 * getfield.(ensobs,:t0))        
+        uwerr.(m_etaprime)
+        uwerr.(m_phi)
+        #uwerr.(f_etaprime)
+        #uwerr.(f_phi)
+catch
+        error("Check ensobs[i].mu_list to see whether all the ensembles have a strange-strange sector.")
+    end
+end
+# sh 
+m_Ds = Vector{uwreal}(undef, length(ensembles))
+m_Ds_star =Vector{uwreal}(undef, length(ensembles))
+#f_Ds= Vector{uwreal}(undef, length(ensembles))
+#f_Ds_star =Vector{uwreal}(undef, length(ensembles))
+if sector["sh"] 
+    try 
+        m_Ds .= getfield.(ensobs, :m_sh_match) .* sqrt.(8 * getfield.(ensobs,:t0))
+        m_Ds_star .= getfield.(ensobs, :m_sh_vec_match) .* sqrt.(8 * getfield.(ensobs,:t0))
+        #f_Ds .= getfield.(ensobs, :f_sh_match) .* sqrt.(8 * getfield.(ensobs,:t0))
+        #f_Ds_star .= getfield.(ensobs, :f_sh_vec_match) .* sqrt.(8 * getfield.(ensobs,:t0))        
+        uwerr.(m_Ds)
+        uwerr.(m_Ds_star)
+        #uwerr.(f_Ds)
+        #uwerr.(f_Ds_star)
+catch
+        error("Check ensobs[i].mu_list to see whether all the ensembles have a strange-heavy sector.")
+    end
+end
+# hh 
+m_etac = Vector{uwreal}(undef, length(ensembles))
+m_jpsi =Vector{uwreal}(undef, length(ensembles))
+#f_etac= Vector{uwreal}(undef, length(ensembles))
+#f_jpsi =Vector{uwreal}(undef, length(ensembles))
+if sector["hh"] 
+    try 
+        m_etac .= getfield.(ensobs, :m_hh_match) .* sqrt.(8 * getfield.(ensobs,:t0))
+        m_jpsi .= getfield.(ensobs, :m_hh_vec_match) .* sqrt.(8 * getfield.(ensobs,:t0))
+        #f_etac .= getfield.(ensobs, :f_hh_match) .* sqrt.(8 * getfield.(ensobs,:t0))
+        #f_jpsi .= getfield.(ensobs, :f_hh_vec_match) .* sqrt.(8 * getfield.(ensobs,:t0))        
+        uwerr.(m_etac)
+        uwerr.(m_jpsi)
+        #uwerr.(f_etac)
+        #uwerr.(f_jpsi)
+catch
+        error("Check ensobs[i].mu_list to see whether all the ensembles have a heavy-heavy sector.")
+    end
+end
+
+
+#===== CONTINUUM AND CHIRAL EXTRAPOLATION =======#
+
+println("\n Continuum and Chiral Extrapolation \n")
+#=
+obs_db =Dict(
+    "muc"=> [muc],
+    "ll" => [m_pi, m_rho, f_pi, f_rho],
+    "ls" => [m_k, m_k_star, f_k, f_k_star],
+    "lh" => [m_D, m_D_star, f_D, f_D_star],
+    "ss" => [m_etaprime, m_phi, f_etaprime, f_phi],
+    "sh" => [m_Ds, m_Ds_star, f_Ds, f_Ds_star],
+    "hh" => [m_etac, m_jpsi, f_etac, f_jpsi]
+)
+=#
+obs_db =Dict(
+    "muc"=> [muc],
+    "ll" => [m_pi, m_rho],
+    "ls" => [m_k, m_k_star],
+    "lh" => [m_D, m_D_star],
+    "ss" => [m_etaprime, m_phi],
+    "sh" => [m_Ds, m_Ds_star],
+    "hh" => [m_etac, m_jpsi]
+)
+
+
+#preparing observables of interest (data, name, ylabel)
+obs = Vector{Vector{uwreal}}(undef,0)
+ttl_obs = Vector{String}(undef,0)
+ylbl = Vector{String}(undef,0)
+push!(obs, obs_db["muc"][1])
+push!(ttl_obs, ttl_obs_db["muc"][1])
+push!(ylbl, ylbl_db["muc"][1])
+
+label = findall(sector)
+for lab in label
+    for i in 1:length(obs_db[lab])
+        push!(obs, obs_db[lab][i])
+        push!(ttl_obs, ttl_obs_db[lab][i])
+        push!(ylbl, ylbl_db[lab][i])
+    end
+end
+# estimate phi2 from database or form analysis if ll sector is computed 
+phi2 = Vector{uwreal}(undef,length(ensembles))
+sector["ll"] ? phi2 = 8 .* getfield.(ensobs,:t0) .* getfield.(ensobs,:m_ll).^2 : phi2 = 8 .* getfield.(ensobs,:t0) .* getfield.(getfield.(ensobs,:ensinfo),:mpi).^2
+phi2_ph = (t0_ph[1] * 139.57039 / hc)^2
+uwerr(phi2_ph)
+#fit routine 
+x = [1 ./(8 .* getfield.(ensobs,:t0)) phi2] #x[1] = a^2 / (8t0), x[2] = 8t0 mpi^2
+uwerr.(x)
+#modify CL+chiral fit  model
+@. cl_chir_model(x, p) = (p[1] + p[2] * x[:, 1]) + (p[3]) * x[:, 2] #linear fits
+
+obs_t0 = Vector{uwreal}(undef, length(obs)) #sqrt(8t0)obs @ CL & phi2_phys
+xarr = Float64.(range(0.0, stop=0.05, length=100))
+for k = 1:length(obs)
+    par = fit_routine(cl_chir_model, value.(x), obs[k], 4)
+    obs_t0[k] = par[1] + par[3] * phi2_ph
+    uwerr(obs_t0[k])    
+    y = Vector{uwreal}(undef, 100)
+    for j=1:100
+        y[j] = par[1] + par[2] * xarr[j] + par[3] * phi2_ph
+    end
+    uwerr.(y)
+    figure()
+    for b in unique(phi2) #select point with same beta
+        nn = findall(x-> x == b, phi2)
+        lgnd = string(L"$m_\pi = $", Int32(round(sqrt(value(b))*hc/t0_ph_value[1], sigdigits=3)), " MeV")
+        errorbar(value.(x[nn,1]), value.(obs[k][nn]), err.(obs[k][nn]), ms=5, fmt="s", mfc="none", mew=0.8, elinewidth=0.8, capsize=2, label=lgnd)
+    end
+    fill_between(xarr, value.(y) .+ err.(y),value.(y) .- err.(y) , color="tomato", alpha=0.4)
+    lgnd=L"$\mathrm{CL}$"
+    errorbar(0, value(obs_t0[k]), err(obs_t0[k]), fmt="s", zorder=2, ms=4,  mfc="red", c = "red", mew=0.8, elinewidth=1, capsize=2, label=lgnd)
+    axvline(0, ls="--", color="black", zorder=1, lw=0.6, label="")
+    xlabel(L"$a^2/8t_0$")
+    ylabel(ylbl[k])
+    xlim(-0.002,0.04)
+    #ylim(2.8,3.4)
+    legend(ncol=2)
+    display(gcf())
+    t = string("fit_", ttl_obs[k], ".pdf")
+    savefig(joinpath(path_plot, t))
+    #savefig(joinpath("/Users/ale/Desktop/plottest", t))
+    close()
+
+end
+
+obs_ph = Vector{uwreal}(undef, length(obs))
+for k = 1:length(obs)
+    println(ylbl[k], " = ", obs_t0[k])
+    #phys
+    obs_ph[k] = obs_t0[k] * hc / t0_ph[1]
+    uwerr(obs_ph[k])
+    println(ttl_obs[k], "(MeV) = ", obs_ph[k])
+end
+
+#======== STORE RESULTS ========#
+
+open(joinpath(path_results, "results.txt"), "w") do f 
+    print(f, "GEVP analysis results", "\n \n")
+    print(f, "#================================# \n \n")
+    print(f, "Setup:\n")
+    print(f, "Ensembles analysed: ", ensembles, "\n")
+    print(f, "Data storage: ", path_data, "\n")
+    print(f, "Plateaux file: ", path_plat, "\n")
+    print(f, "Results storage: ", path_results, "\n")
+    print(f, "Sector analysed: ", findall(sector), "\n")
+    print(f, "Reweighting factor: ", rwf, "\n")
+    print(f, "Mass shift: ", mass_shift, "\n")
+    if compute_t0
+        print(f, "t0 values: computed from the ms.dat files \n")
+    else
+        print(f, "t0 values: extracted from const.jl \n")
+    end
+    if sector["ll"]
+        print(f, "Phi2 = 8t0m_pi^2: m_pi computed from mesons.dat files \n")
+    else
+        print(f, "Phi2 = 8t0m_pi^2: m_pi extracted from ens_db in const.jl \n")
+    end
+    print(f, " \n#================================# \n \n")
+    print(f, "\n Results at finite lattice spacings \n \n")
+
+    lab = findall(sector)
+    for i in 1:length(ensembles)
+        print(f, "Ensemble ", ensobs[i].ensinfo.id, "\n")
+        for j in 1:length(obs)
+            print(f, ylbl[j], " = ", obs[j][i], "\n")
+        end
+        print(f, "\n")
+    end
+
+    print(f, " \n#================================# \n \n")
+    print(f, "\n Results in the continuum limit \n \n")
+    for k = 1:length(obs)
+        print(f, ylbl[k], " = ", obs_t0[k], "\n")
+        #phys
+        print(f, ttl_obs[k], "(MeV) = ", obs_ph[k], "\n\n")
+    end
+end
+
+
+
+
+## ####################################
+## NOT CONSIDER 
+data = read_mesons(path_dat, "G5", "G5")
+corr_pp = corr_obs.(data, L=64)
+
+##
+hh_pp = corr_pp[6].obs
+##
+pp11 = hh_pp[97:end-7]
+pp12 = hh_pp[98:end-6]
+pp13 = hh_pp[99:end-5]
+pp14 = hh_pp[100:end-4]
+pp22 = hh_pp[99:end-5]
+pp23 = hh_pp[100:end-4]
+pp24 = hh_pp[101:end-3]
+pp33 = hh_pp[101:end-3]
+pp34 = hh_pp[102: end-2]
+pp44 = hh_pp[103:end-1]
+pp_sing = get_matrix([pp11, pp22, pp33, pp44], [pp12, pp13, pp14, pp23, pp24, pp34])
+##
+tau=4
+pp11 = hh_pp[97:end-1-4*tau]
+pp12 = hh_pp[97+tau:end-1-3*tau]
+pp13 = hh_pp[97+2*tau:end-1-2*tau]
+pp14 = hh_pp[100:end-4]
+pp22 = hh_pp[97+2*tau:end-1-2*tau]
+pp23 = hh_pp[97+3*tau:end-1-tau]
+pp24 = hh_pp[101:end-1]
+pp33 = hh_pp[97+4*tau:end-1]
+pp34 = hh_pp[102: end-2]
+pp44 = hh_pp[103:end-1]
+pp_sing = get_matrix([pp11, pp22], [pp12])
+
+##
+
+##test lukas
+evecs  = uweigvecs(pp_sing[3], pp_sing[2])
+proj_pp = Vector{uwreal}(undef, length(pp_sing))
+for i = 1:length(pp_sing)
+    svec = evecs[:,1]
+    proj_pp[i]  = uwdot(svec, uwdot(pp_sing[i], svec))[1]
+    uwerr(proj_pp[i])
+end
+uwerr.(pp11)
+fig = figure()
+    errorbar(collect(1:length(proj_pp)), value.(proj_pp), err.(proj_pp), fmt="*",mec="blue", ecolor="blue")
+    errorbar(collect(1:length(pp11)), value.(pp11), err.(pp11), fmt="+",mec="red", ecolor="red")
+    yscale("log")
+    #ylim(-0.00000000001,0.0000000001)
+display(fig )
+##lukas eff mass
+@.model(x,p) = p[1]*exp(-x*p[2]) + p[3]*exp(-x*p[4])
+ydata =  proj_pp[6:end-1]
+uwerr.(ydata)
+xdata = collect(1:length(ydata))
+par = fit_routine(model, xdata, ydata,4)
+
+##
+mlmass = meff(proj_pp, [120-96, 155-96])
+uwerr(lmass)
+
+##
+evals = getall_eigvals(pp_sing, 1)
+evecs = getall_eigvecs(pp_sing, 1)
+##
+en = energies(evals)
+fig = figure()
+    errorbar( collect(1:length(en[1])),value.(en[1]), yerr=err.(en[1]), fmt="+")
+    ylim(0.71,0.76)
+display(fig)
+## fitting the energy, very good
+@.model(x,p) = p[1] + p[2]*exp(-x*(p[3]))
+en1 = en[1][2:end-4]
+xdata = collect(1:length(en1))
+wpmm = Dict{Int64, Vector{Float64}}()
+wpmm[303] = [-1.0, 1.5,-1.0, -1.0]
+par = fit_routine(model,  xdata, en1, 3, wpm=wpmm)
+
+## try to extract matrix elements
+
+function norm_evec(evec::Matrix{uwreal}, ctnot::Matrix{uwreal})
+    n=size(evec,1)
+    res_evec=Matrix{uwreal}(undef, n, n)
+    for i =1:n
+        aux_norm = (uwdot(evec[:,i], uwdot(ctnot, evec[:,i])).^2).^0.25
+        res_evec[i,:] = 1 ./aux_norm .* evec[:,i]
+    end
+    return res_evec
+end
+
+#1)normalise evecs and extract Ptilda for all times
+
+norm = Array{Matrix{uwreal}}(undef, length(evecs))
+for i=1:length(evecs)
+    norm[i] = norm_evec(evecs[i], pp_sing[4])
+end
+
+#2) construct matrix M(t,t_0)=c(t_0)Ptilda for all times 
+mtt0 = Array{Matrix{uwreal}}(undef, length(evecs))
+for i=1:length(evecs)
+    mtt0[i] = uwdot(pp_sing[4], juobs.transpose(norm[i]))
+end
+
+## extract decays 
+function pseudo_mat_elem(evec::Array{Matrix{uwreal}}, mass::uwreal, mu::Array{Float64}, t0::Int64)
+    Rn = [exp(mass * (t0)/2) * evec[t][1] for t =1:length(evec)]
+    aux = sqrt(2)*sum(mu)  / mass^1.5 
+    return uwdot(aux,Rn)
+end
+
+test_dec = pseudo_mat_elem(mtt0, par[1], [0.14,0.14],4)
+## Standard 
+mass = meff(hh_pp, [123, 155] )
+f = dec_const_pcvc(hh_pp, [123,155], mass, [0.14,0.14],96)
+
+
+## testing generalised
+todiag = uwdot(juobs.invert(pp_sing[4]),pp_sing[10])
+evals, evecs = uweigen(todiag)
+##
+#test = uwdot(juobs.invert(evecs), uwdot(uwdot(juobs.invert(pp_sing[4]),pp_sing[10]), evecs))
+test = uwdot(evecs, uwdot(evals,juobs.invert(evecs)))
+##
+evals2, evecs2 = uweigen(pp_sing[10], pp_sing[4])
+#test2 = uwdot(juobs.invert(evecs2), uwdot(uwdot(juobs.invert(pp_sing[4]),pp_sing[10]), evecs2))
+test2 = uwdot(pp_sing[4],uwdot(evecs2, uwdot(evals2,juobs.invert(evecs2))))
\ No newline at end of file

tools.jl

+@doc raw"""
+    read_data(ens::String, g1::String="G5", g2::String="G5")
+Given the ensemble id and the Dirac structure, this method reads the data in path_data and
+returns a CData object. Replicas are automatically detected and loaded if present.
+"""
+function read_data(ens::String, g1::String="G5", g2::String="G5"; legacy::Bool=false)
+    path = joinpath(path_data, ens)
+    rep = filter(x->occursin("mesons.dat", x), readdir(path, join=true))
+    if length(rep)!=0    
+        length(rep)==1 ? (return read_mesons(rep[1], g1, g2, legacy=legacy)) : (return read_mesons(rep, g1, g2, legacy=legacy))
+    else
+        error("mesons.dat file not found for ensemble ", ens, " in path ", path)
+    end
+end
+function read_rw(ens::String)
+    path = joinpath(path_data, ens)
+    rep = filter(x->occursin("ms1.dat", x), readdir(path, join=true))
+    if length(rep)!=0
+        try
+            length(rep) == 1 ? (return read_ms1(rep[1])) : (return read_ms1.(rep)) 
+        catch
+            length(rep) == 1 && length(rep)!=0 ? (return read_ms1(rep[1], v="1.4")) : (return read_ms1.(rep, v="1.4")) 
+        end
+    else
+        error("ms1.dat file not found for ensemble ", ens, " in path ", path)
+    end
+end
+function read_t0(ens::String)
+    path = joinpath(path_data, ens)
+    rep = filter(x->occursin("ms.dat", x), readdir(path, join=true))
+    if length(rep)!=0
+        length(rep) == 1 ? (return read_ms(rep[1])) : (return read_ms.(rep)) 
+    else
+        error("ms.dat file not found for ensemble ", ens, " in path ", path)
+    end
+end
+function get_corr(ens::EnsInfo, g1::String="G5", g2::String="G5"; rw::Bool=false, legacy::Bool=false)
+    aux = read_data(ens.id, g1, g2, legacy=legacy)
+    !rw ? obs = corr_obs.(aux, L=ens.L) : obs = corr_obs.(aux, L=ens.L, rw=read_rw(ens.id))
+    return obs
+end
+function get_t0(ens::EnsInfo; rw::Bool=false, pl::Bool=false)
+    data = read_t0(ens.id)
+    plat = select_plateau(ens) 
+    !rw ? t0 = comp_t0(data, plat, L=ens.L, pl=pl) : t0 = comp_t0(data, plat, L=ens.L, pl=pl, rw=read_rw(ens.id))
+    return t0
+end
+function get_mu(mu_list::Vector{Vector{Float64}}, deg::Bool)
+    mu_sorted = unique(sort(minimum.(mu_list)))
+
+    mul = mu_sorted[1]
+    deg ? mus = 0.0 : mus = mu_sorted[2]
+    muh = unique(maximum.(mu_list))
+    muh = filter(x-> x > mul && x > mus, muh)
+
+    return mul, mus, muh
+end
+function get_ll(mu_list, obs::Vector{uwreal}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mul in mu && mu[1] == mu[2] #l-l
+            return obs[k]
+        end
+    end
+end
+function get_ls(mu_list, obs::Vector{uwreal}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mul in mu && mus in mu  #l-l
+            return obs[k]
+        end
+    end
+end
+function get_lh(mu_list, obs::Vector{uwreal}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    obs_lh = Vector{uwreal}(undef, 0)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mul in mu && mu[1] != mu[2] && !(mus in mu) #l-h
+            push!(obs_lh, obs[k])
+        end
+    end
+    return obs_lh
+end
+function get_ss(mu_list, obs::Vector{uwreal}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mus in mu && mu[1] == mu[2] #s-s
+            return obs[k]
+        end
+    end
+end
+function get_sh(mu_list, obs::Vector{uwreal}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    obs_sh = Vector{uwreal}(undef, 0)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mus in mu && mu[1] != mu[2] && !(mul in mu)#s-h
+            push!(obs_sh, obs[k])
+        end
+    end
+    return obs_sh
+end
+function get_hh(mu_list, obs::Vector{uwreal}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    obs_hh = Vector{uwreal}(undef,0)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        for i =1:length(muh)
+            if muh[i] in mu && mu[1] == mu[2]#h-h
+                push!(obs_hh, obs[k])
+            end
+        end
+    end
+    return obs_hh
+end
+function comp_mat(ensinfo::EnsInfo, tot_obs::Vector{juobs.Corr}, tau::Int64)
+    #wanted_corr = relevant_corr(ensinfo, tot_obs)
+    mat_info = Vector{MatInfo}(undef, length(tot_obs))
+    for i = 1:length(tot_obs)
+        a11 = tot_obs[i].obs[1:end-1-2*tau]
+        a12 = tot_obs[i].obs[1+tau:end-1-tau]
+        a22 = tot_obs[i].obs[1+2*tau:end-1]
+        aux = juobs.get_matrix([a11, a22], [a12])
+        mat_info[i] = MatInfo(ensinfo, aux, tot_obs[i].mu, tot_obs[i].y0)
+    end
+    return mat_info
+end
+function relevant_corr(ensinfo::EnsInfo, tot_obs::Vector{juobs.Corr})
+    mulist = getfield.(tot_obs, :mu)
+    res_obs = Vector{juobs.Corr}(undef, 0)
+    if sector["ll"]
+        push!(res_obs, get_ll(mulist, tot_obs, ensinfo.deg))
+    end
+    if sector["ls"]
+        push!(res_obs, get_ls(mulist, tot_obs, ensinfo.deg))
+    end
+    if sector["lh"]
+        toadd = get_lh(mulist, tot_obs, ensinfo.deg)
+        [push!(res_obs, toadd[i]) for i=1:length(toadd)]
+    end
+    if sector["ss"]
+        push!(res_obs, get_ss(mulist, tot_obs, ensinfo.deg))
+    end
+    if sector["sh"]
+        toadd = get_sh(mulist, tot_obs, ensinfo.deg)
+        [push!(res_obs, toadd[i]) for i=1:length(toadd)]  
+    end
+    if sector["hh"]
+        toadd = get_hh(mulist, tot_obs, ensinfo.deg)
+        [push!(res_obs, toadd[i]) for i=1:length(toadd)]
+    end
+    return res_obs
+end
+function gevp_mass(mat_obs::Vector{MatInfo}; t0::Int64=2, delta_t_in::Array{Int64}=[2,5], wpm::Union{Nothing, Dict{Int64,Vector{Float64}}}=nothing)
+    y0s = getfield.(mat_obs, :y0)
+    mu_list = getfield.(mat_obs, :mu) 
+    plat_en0 = Vector{uwreal}(undef, length(mat_obs))
+    plat_en1 = Vector{uwreal}(undef, length(mat_obs))
+    for i in 1:length(mat_obs)
+        mat = getfield(mat_obs[i], :mat_list)[y0s[i]+1:end-1] #matrices to diagonalise for a given sector [i]
+        evals = getall_eigvals(mat, t0)
+        en = energies(evals) #ground and excited state energy
+        aux_mass = Vector{uwreal}(undef,0)
+        for t in delta_t_in[1]:delta_t_in[2]
+            try
+                temp = fit_mass(en[1],t_in=t,  wpm=wpm)
+                push!(aux_mass, temp)
+            catch
+                plat = select_plateau(mat_obs[i].ensinfo, mu_list)[end-2] .- y0s[i]
+                temp = plat_av(en[1], plat)
+                uwerr(temp)
+                push!(aux_mass, temp)
+                println("\n Error analysis failed for ensemble ", mat_obs[i].ensinfo.id, " in the mu sector ", mat_obs[i].mu, " with time fit values :", t)
+                println("\n The mass has been computed with a plateau average. Plateaux has been taken from select plateau file.\n")
+            end
+        end
+        plat_en0[i] = param_av(aux_mass)
+    end
+    return plat_en0
+end
+function fit_mass(en_i::Array{uwreal}; t_in::Int64=3, wpm::Union{Nothing, Dict{Int64,Vector{Float64}}}=nothing)
+    @.model(x,p) = p[1] + p[2]*exp(-x*p[3])
+    ydata = en_i[t_in:end-5]
+    xdata = collect(1:length(ydata))
+    par = fit_routine(model, xdata, ydata, 3, wpm=wpm)
+    println("\nPar[1] is ", par[1] )
+    return par[1]
+end
+function param_av(fit_masses::Vector{uwreal})
+    w = 1 ./ err.(fit_masses).^2
+    av = sum(w .* fit_masses) / sum(w)
+    return av  
+end
+function get_meff(obs::Vector{juobs.Corr}, ens::EnsInfo; pl::Bool=false)
+    mu_list  = getfield.(obs, :mu)
+    plat = select_plateau(ens, mu_list)
+    return meff.(obs, plat, pl=pl)
+end
+function get_f(obs::Vector{juobs.Corr}, m::Vector{uwreal}, ens::EnsInfo; pl::Bool=false)
+    mu_list = getfield.(obs, :mu)
+    plat = select_plateau(ens, mu_list)
+    return dec_const_pcvc.(obs, plat, m, pl=pl)
+end
+function select_plateau(ensinfo::EnsInfo)
+    ens = ensinfo.id
+    f = readdlm(path_plat)
+    head  = String.(f[:,1])
+    delim = findall(x->occursin("#", x), head)
+    edelim = findfirst(x->occursin(ens, x), head)
+    cdelim = findfirst(x-> x.== edelim, delim)
+    del = delim[cdelim]+1 : delim[cdelim+1]-1
+    head_ = String.(f[del, 1])
+    plat_ = Int64.(f[del, 2:3])
+    plat = Vector{Int64}(undef, 2)
+    n = findfirst(x-> occursin("t0", x), head_)
+    plat[1] = plat_[n,1]
+    plat[2] = plat_[n,2]
+    return plat
+end
+function select_plateau(ensinf::EnsInfo, mu_list)
+    ens =ensinf.id
+    deg = ensinf.deg
+    mul, mus, muh = get_mu(mu_list, deg)
+    f = readdlm(path_plat)
+    head = String.(f[:, 1])
+
+    delim = findall(x-> occursin("#", x ), head)
+    edelim = findfirst(x-> occursin(ens, x), head)
+    cdelim = findfirst(x-> x.== edelim, delim)
+
+    del = delim[cdelim]+1 : delim[cdelim+1]-1
+    #del = delim
+    plat = Vector{Vector{Int64}}(undef, 0)
+    #plat = Dict()
+    plat_ = Int64.(f[del, 2:3])
+    head_ = String.(f[del, 1])
+    
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mul in mu && mu[1] != mu[2] && !(mus in mu)  #heavy-light
+            n = findfirst(x-> occursin("lh", x), head_)
+            push!(plat,  [plat_[n, 1], plat_[n, 2]]  ) 
+        elseif mul in mu && mu[1] == mu[2] #light-light
+            n = findfirst(x-> occursin("ll", x), head_)
+            push!(plat,  [plat_[n, 1], plat_[n, 2]] )               
+        elseif mul in mu && mus in mu#strange-light
+            n = findfirst(x-> occursin("ls", x), head_)
+            push!(plat,  [plat_[n, 1], plat_[n, 2]] )
+        elseif mus in mu && mu[1] != mu[2] && !(mul in mu)#heavy-strange
+            n = findfirst(x-> occursin("sh", x ), head_)
+            push!(plat, [plat_[n, 1], plat_[n, 2]] )
+        elseif mus in mu && mu[1] == mu[2] && !(mul in mu)#strange-strange
+            n = findfirst(x-> occursin("ss", x), head_)
+            push!(plat, [plat_[n, 1], plat_[n, 2]] )
+        elseif !(mul in mu) && !(mus in mu) #heavy-heavy
+            n = findfirst(x-> occursin("hh", x), head_)
+            push!(plat, [plat_[n, 1], plat_[n, 2]] )
+        end
+    end
+    return plat
+end
+
+function match_muc(muh, m_lh, m_sh, target)
+    M = (2/3 .* m_lh .+ 1/3 .* m_sh) 
+    par, chi2exp = lin_fit(muh, M)
+    muh_target = x_lin_fit(par, target)
+    return muh_target
+end
+
+function get_ll(mu_list, obs::Array{juobs.Corr}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mul in mu && mu[1] == mu[2] #l-l
+            return obs[k]
+        end
+    end
+end
+function get_ls(mu_list, obs::Array{juobs.Corr}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mul in mu && mus in mu  #l-l
+            return obs[k]
+        end
+    end
+end
+function get_lh(mu_list, obs::Array{juobs.Corr}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    obs_lh = Array{juobs.Corr}(undef, 0)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mul in mu && mu[1] != mu[2] && !(mus in mu) #l-h
+            push!(obs_lh, obs[k])
+        end
+    end
+    return obs_lh
+end
+function get_ss(mu_list, obs::Array{juobs.Corr}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mus in mu && mu[1] == mu[2] #s-s
+            return obs[k]
+        end
+    end
+end
+function get_sh(mu_list, obs::Array{juobs.Corr}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    obs_sh = Array{juobs.Corr}(undef, 0)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mus in mu && mu[1] != mu[2] && !(mul in mu)#s-h
+            push!(obs_sh, obs[k])
+        end
+    end
+    return obs_sh
+end
+function get_hh(mu_list, obs::Array{juobs.Corr}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    obs_hh = Array{juobs.Corr}(undef,0)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        for i =1:length(muh)
+            if muh[i] in mu && mu[1] == mu[2]#h-h
+                push!(obs_hh, obs[k])
+            end
+        end
+    end
+    return obs_hh
+end
\ No newline at end of file

types.jl

+mutable struct EnsInfo
+    id::String
+    L::Int64
+    beta::Float64
+    deg::Bool
+    mpi::Float64
+    function EnsInfo(ens_id::String, info::Vector{Float64})
+        id = ens_id
+        L = info[1]
+        beta = info[2]
+        deg = info[3]
+        mpi = info[4]
+        return new(id, L, beta, deg, mpi)
+    end
+end
+mutable struct EnsObs
+    ensinfo::EnsInfo
+    mu_list:: Vector{Vector{Float64}}
+    is_pseudo::Bool
+
+    m_ll::Union{Nothing,uwreal}
+    m_ls::Union{Nothing,uwreal}
+    m_lh::Union{Nothing,Vector{uwreal}}
+    m_ss::Union{Nothing,uwreal}
+    m_sh::Union{Nothing,Vector{uwreal}}
+    m_hh::Union{Nothing,Vector{uwreal}}
+    f_ll::Union{Nothing,uwreal}
+    f_ls::Union{Nothing,uwreal}
+    f_lh::Union{Nothing,Vector{uwreal}}
+    f_ss::Union{Nothing,uwreal}
+    f_sh::Union{Nothing,Vector{uwreal}}
+    f_hh::Union{Nothing,Vector{uwreal}}
+
+    m_ll_vec::Union{Nothing,uwreal}
+    m_ls_vec::Union{Nothing,uwreal}
+    m_lh_vec::Union{Nothing,Vector{uwreal}}
+    m_ss_vec::Union{Nothing,uwreal}
+    m_sh_vec::Union{Nothing,Vector{uwreal}}
+    m_hh_vec::Union{Nothing,Vector{uwreal}}
+    f_ll_vec::Union{Nothing,uwreal}
+    f_ls_vec::Union{Nothing,uwreal}
+    f_lh_vec::Union{Nothing,Vector{uwreal}}
+    f_ss_vec::Union{Nothing,uwreal}
+    f_sh_vec::Union{Nothing,Vector{uwreal}}
+    f_hh_vec::Union{Nothing,Vector{uwreal}}
+
+    m_lh_match::Union{Nothing,uwreal}
+    m_sh_match::Union{Nothing,uwreal}
+    m_hh_match::Union{Nothing,uwreal}
+    f_lh_match::Union{Nothing,uwreal}
+    f_sh_match::Union{Nothing,uwreal}
+    f_hh_match::Union{Nothing,uwreal}
+
+    m_lh_vec_match::Union{Nothing,uwreal}
+    m_sh_vec_match::Union{Nothing,uwreal}
+    m_hh_vec_match::Union{Nothing,uwreal}
+    f_lh_vec_match::Union{Nothing,uwreal}
+    f_sh_vec_match::Union{Nothing,uwreal}
+    f_hh_vec_match::Union{Nothing,uwreal}
+
+    muh_target::Union{Nothing,uwreal}
+    a::Union{Nothing,uwreal}
+    t0::Union{Nothing,uwreal}
+
+    
+    function EnsObs(ens::EnsInfo, mu::Vector{Vector{Float64}}, m_ps::Vector{uwreal}, f_ps::Vector{uwreal}, m_vec::Vector{uwreal}, f_vec::Vector{uwreal})
+        a = new()
+        a.ensinfo = ens
+        a.mu_list  = mu
+        a.is_pseudo = true
+
+        a.m_ll = get_ll(a.mu_list, m_ps,a.ensinfo.deg)
+        a.ensinfo.deg ? a.m_ls = a.m_ll : a.m_ls = get_ls(a.mu_list, m_ps, a.ensinfo.deg)
+        a.ensinfo.deg ? a.m_ss = a.m_ll : a.m_ss = get_ss(a.mu_list, m_ps, a.ensinfo.deg)
+        a.m_lh = get_lh(a.mu_list, m_ps,a.ensinfo.deg)
+        a.ensinfo.deg ? a.m_sh = a.m_lh : a.m_sh = get_sh(a.mu_list, m_ps, a.ensinfo.deg)
+        a.m_hh = get_hh(a.mu_list, m_ps,a.ensinfo.deg)
+        a.f_ll = get_ll(a.mu_list, f_ps, a.ensinfo.deg)
+        a.ensinfo.deg ? a.f_ls = a.f_ll : a.f_ls = get_ls(a.mu_list, f_ps, a.ensinfo.deg)
+        a.ensinfo.deg ? a.f_ss = a.f_ll : a.f_ss = get_ss(a.mu_list, f_ps, a.ensinfo.deg)
+        a.f_lh = get_lh(a.mu_list, f_ps,a.ensinfo.deg)
+        a.ensinfo.deg ? a.f_sh = a.f_lh : a.f_sh = get_sh(a.mu_list, f_ps, a.ensinfo.deg)
+        a.f_hh = get_hh(a.mu_list, f_ps,a.ensinfo.deg)
+
+        a.m_ll_vec = get_ll(a.mu_list, m_vec,a.ensinfo.deg)
+        a.ensinfo.deg ? a.m_ls_vec = a.m_ll_vec : a.m_ls_vec = get_ls(a.mu_list, m_vec, a.ensinfo.deg)
+        a.ensinfo.deg ? a.m_ss_vec = a.m_ll_vec : a.m_ss_vec = get_ss(a.mu_list, m_vec, a.ensinfo.deg)
+        a.m_lh_vec = get_lh(a.mu_list, m_vec, a.ensinfo.deg)
+        a.ensinfo.deg ? a.m_sh_vec = a.m_lh_vec : a.m_sh_vec = get_sh(a.mu_list, m_vec, a.ensinfo.deg)
+        a.m_hh_vec = get_hh(a.mu_list, m_vec, a.ensinfo.deg)
+        a.f_ll_vec = get_ll(a.mu_list, f_vec, a.ensinfo.deg)
+        a.ensinfo.deg ? a.f_ls_vec = a.f_ll_vec : a.f_ls_vec = get_ls(a.mu_list, f_vec, a.ensinfo.deg)
+        a.ensinfo.deg ? a.f_ss_vec = a.f_ll_vec : a.f_ss_vec = get_ss(a.mu_list, f_vec, a.ensinfo.deg)
+        a.f_lh_vec = get_lh(a.mu_list, f_vec,a.ensinfo.deg)
+        a.ensinfo.deg ? a.f_sh_vec = a.f_lh_vec : a.f_sh_vec = get_sh(a.mu_list, f_vec, a.ensinfo.deg)
+        a.f_hh_vec = get_hh(a.mu_list, f_vec, a.ensinfo.deg)
+        
+        a.m_lh_match = nothing
+        a.m_sh_match = nothing
+        a.m_hh_match = nothing
+        a.f_lh_match = nothing
+        a.f_sh_match = nothing
+        a.f_hh_match = nothing
+        
+        a.m_lh_vec_match = nothing
+        a.m_sh_vec_match = nothing
+        a.m_hh_vec_match = nothing
+        a.f_lh_vec_match = nothing
+        a.f_sh_vec_match = nothing
+        a.f_hh_vec_match = nothing
+        
+        a.muh_target = nothing
+        a.a = nothing
+        a.t0 = nothing
+        return a 
+    end
+    function EnsObs(ens::EnsInfo, mu::Vector{Vector{Float64}}, m_ps::Vector{uwreal}, m_vec::Vector{uwreal})
+        a = new()
+        a.ensinfo = ens
+        a.mu_list  = mu
+        a.is_pseudo = true
+
+        a.m_ll = get_ll(a.mu_list, m_ps,a.ensinfo.deg)
+        a.ensinfo.deg ? a.m_ls = a.m_ll : a.m_ls = get_ls(a.mu_list, m_ps, a.ensinfo.deg)
+        a.ensinfo.deg ? a.m_ss = a.m_ll : a.m_ss = get_ss(a.mu_list, m_ps, a.ensinfo.deg)
+        a.m_lh = get_lh(a.mu_list, m_ps,a.ensinfo.deg)
+        a.ensinfo.deg ? a.m_sh = a.m_lh : a.m_sh = get_sh(a.mu_list, m_ps, a.ensinfo.deg)
+        a.m_hh = get_hh(a.mu_list, m_ps,a.ensinfo.deg)
+        
+        a.m_ll_vec = get_ll(a.mu_list, m_vec,a.ensinfo.deg)
+        a.ensinfo.deg ? a.m_ls_vec = a.m_ll_vec : a.m_ls_vec = get_ls(a.mu_list, m_vec, a.ensinfo.deg)
+        a.ensinfo.deg ? a.m_ss_vec = a.m_ll_vec : a.m_ss_vec = get_ss(a.mu_list, m_vec, a.ensinfo.deg)
+        a.m_lh_vec = get_lh(a.mu_list, m_vec, a.ensinfo.deg)
+        a.ensinfo.deg ? a.m_sh_vec = a.m_lh_vec : a.m_sh_vec = get_sh(a.mu_list, m_vec, a.ensinfo.deg)
+        a.m_hh_vec = get_hh(a.mu_list, m_vec, a.ensinfo.deg)
+        
+        a.m_lh_match = nothing
+        a.m_sh_match = nothing
+        a.m_hh_match = nothing
+        
+        a.m_lh_vec_match = nothing
+        a.m_sh_vec_match = nothing
+        a.m_hh_vec_match = nothing
+        
+        a.muh_target = nothing
+        a.a = nothing
+        a.t0 = nothing
+        return a 
+    end
+end
+mutable struct MatInfo
+    ensinfo::EnsInfo
+    mu::Vector{Float64}
+    mat_list::Vector{Matrix{uwreal}}
+    y0::Int64
+    function MatInfo(_ensinfo::EnsInfo, _mat_list::Array{Array{T,2} where T,1}, _mu::Vector{Float64}, _y0::Int64)
+        a = new()
+        a.ensinfo = _ensinfo
+        a.mu = _mu
+        a.mat_list = _mat_list
+        a.y0 = _y0
+        return a 
+    end
+end
+