Updated juobs_uwreal_const to use Ref{uwreal}

d3e648ea · Antonino D'Anna · 8fe484e9 · d3e648ea · d3e648ea · d3e648ea
Commit d3e648ea authored Feb 03, 2025 by Antonino D'Anna
39 changed files
--- a/.gitignore
+++ b/.gitignore
 analysis/plat\.txt
 *.pdf
 analysis/*.txt
--- a/Manifest.toml
+++ b/Manifest.toml
--- a/README.md
+++ b/README.md
 # juobs
 Analysis code (Julia)
\ No newline at end of file
--- a/analysis/analysis.jl
+++ b/analysis/analysis.jl
 #TODO include rw, different plateaux depending on obs, print chi2, compute t0, compute mpi
 using juobs, ADerrors, DelimitedFiles, PyPlot, LaTeXStrings
 const path = "/home/javier/Lattice/charm/production_2"
 const path_plat = "/home/javier/Lattice/juobs/analysis/plat.txt"
 const path_plot = "/home/javier/Lattice/juobs/analysis/plots"
 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"]
 include("/home/javier/Lattice/juobs/constants/juobs_const.jl")
 include("/home/javier/Lattice/juobs/analysis/functions.jl")
 const phi2 = 8 .* t0.(beta) .* [0.16345, 0.09222, 0.11224, 0.10630, 0.06514].^2 #8t0 m_pi^2
 m_lh = Vector{Vector{uwreal}}(undef, length(ensembles))
 m_sh = Vector{Vector{uwreal}}(undef, length(ensembles))
 m_lh_star = Vector{Vector{uwreal}}(undef, length(ensembles))
 m_sh_star = Vector{Vector{uwreal}}(undef, length(ensembles))
 f_lh = Vector{Vector{uwreal}}(undef, length(ensembles))
 f_sh = Vector{Vector{uwreal}}(undef, length(ensembles))
 mu_pp = Vector{Vector{Vector{Float64}}}(undef, length(ensembles))
 mu_aa = Vector{Vector{Vector{Float64}}}(undef, length(ensembles))
 ###########################
 ###### COMPUTATION ########
 ###########################
 for iens = 1:length(ensembles)
    pp = read_dat(R[iens], "G5", "G5")
    aa1 = read_dat(R[iens], "G1G5", "G1G5")
    pp_obs = corr_obs.(pp, L=L[iens])
    aa1_obs = corr_obs.(aa1, L=L[iens])
    mu_pp[iens] = getfield.(pp_obs, :mu)
    mu_aa[iens] = getfield.(aa1_obs, :mu)
    m = comp_meff(pp_obs, deg[iens], ensembles[iens])
    m_star = comp_meff(aa1_obs, deg[iens], ensembles[iens])
    f = comp_f(pp_obs, m, deg[iens], ensembles[iens])
    m_lh[iens] = get_lh(mu_pp[iens], m, deg[iens])
    deg[iens] ? m_sh[iens] = m_lh[iens] : m_sh[iens] = get_sh(mu_pp[iens], m, deg[iens])
    m_lh_star[iens] = get_lh(mu_aa[iens], m_star, deg[iens])
    deg[iens] ? m_sh_star[iens] = m_lh_star[iens] : m_sh_star[iens] = get_sh(mu_aa[iens], m_star, deg[iens])
    f_lh[iens] = get_lh(mu_pp[iens], f, deg[iens])
    deg[iens] ? f_sh[iens] = f_lh[iens] : f_sh[iens] = get_sh(mu_pp[iens], f, deg[iens])
 end
 mm = get_mu.(mu_pp, deg)
 mul_pp = getindex.(mm, 1)
 mus_pp = getindex.(mm, 2)
 muh_pp = getindex.(mm, 3)
 mm = get_mu.(mu_aa, deg)
 mul_aa = getindex.(mm, 1)
 mus_aa = getindex.(mm, 2)
 muh_aa = getindex.(mm, 3)
 m_lh_match = Vector{uwreal}(undef, length(ensembles))
 m_sh_match = Vector{uwreal}(undef, length(ensembles))
 m_lh_v_match = Vector{uwreal}(undef, length(ensembles))
 m_sh_v_match = Vector{uwreal}(undef, length(ensembles))
 f_lh_match = Vector{uwreal}(undef, length(ensembles))
 f_sh_match = Vector{uwreal}(undef, length(ensembles))
 muh_target = Vector{uwreal}(undef, length(ensembles))
 ###########################
 ###### MATCHING ###########
 ###########################
 for iens = 1:length(ensembles)
    target = a(beta[iens]) * (2*M[1] + 6*M[2] + M[3] + 3*M[4]) / (12*hc)
    if !deg[iens] 
        muh_target[iens] = match_muc(muh_pp[iens], m_lh[iens], m_sh[iens], m_lh_star[iens], m_sh_star[iens], target)
    else
        muh_target[iens] = match_muc(muh_pp[iens], m_lh[iens], m_lh_star[iens], target)
    end
    uwerr(muh_target[iens])
    #Interpolate m_lh m_lh_star, m_sh, m_sh_tar
    par, chi2exp = lin_fit(muh_pp[iens], m_lh[iens])
    m_lh_match[iens] = y_lin_fit(par, muh_target[iens])
    par, chi2exp = lin_fit(muh_aa[iens], m_lh_star[iens])
    m_lh_v_match[iens] = y_lin_fit(par, muh_target[iens])
    uwerr.(m_lh[iens])
    uwerr.(m_lh_star[iens])
    uwerr(m_lh_match[iens])
    uwerr(m_lh_v_match[iens])
    if !deg[iens]
        par, chi2exp = lin_fit(muh_pp[iens], m_sh[iens])
        m_sh_match[iens] = y_lin_fit(par, muh_target[iens])
        par, chi2exp = lin_fit(muh_aa[iens], m_sh_star[iens])
        m_sh_v_match[iens] = y_lin_fit(par, muh_target[iens])
        uwerr.(m_sh[iens])
        uwerr.(m_sh_star[iens])
        uwerr(m_sh_match[iens])
        uwerr(m_sh_v_match[iens])
    else
        m_sh_match[iens] = m_lh_match[iens]
        m_sh_v_match[iens] = m_lh_v_match[iens]
    end
    #Interpolate f_lh, f_sh
    par, chi2exp = lin_fit(muh_pp[iens], f_lh[iens])
    f_lh_match[iens] = y_lin_fit(par, muh_target[iens])
    uwerr.(f_lh[iens])
    uwerr(f_lh_match[iens])
    if !deg[iens]
        par, chi2exp = lin_fit(muh_pp[iens], f_sh[iens])
        f_sh_match[iens] = y_lin_fit(par, muh_target[iens])
        uwerr.(f_sh[iens])
        uwerr(f_sh_match[iens])
    else
        f_sh_match[iens] = f_lh_match[iens]
    end
 end
 ###########################
 ###### PLOTS ##############
 ###########################
 for iens = 1:length(ensembles)
    #m_lh m_lh_star
    figure()
    title(ensembles[iens])
    xlabel(L"$a\mu$")
    ylabel(L"$aM$")
    errorbar(muh_pp[iens], value.(m_lh[iens]), err.(m_lh[iens]), fmt="x")
    errorbar(value(muh_target[iens]), value(m_lh_match[iens]), err(m_lh_match[iens]), err(muh_target[iens]), fmt="x")
    errorbar(muh_aa[iens], value.(m_lh_star[iens]), err.(m_lh_star[iens]), fmt="x")
    errorbar(value(muh_target[iens]), value(m_lh_v_match[iens]), err(m_lh_v_match[iens]), err(muh_target[iens]), fmt="x")  
    legend([L"$m_D$", L"$m_D (\mathrm{int})$", L"$m_{D^*}$", L"$m_{D^*} (\mathrm{int})$"])  
    display(gcf())
    t = string(ensembles[iens], "_mD.pdf")
    savefig(joinpath(path_plot, t))
    close()
    #m_sh m_sh_star
    if !deg[iens]
        figure()
        title(ensembles[iens])
        xlabel(L"$a\mu$")
        ylabel(L"$aM$")
        errorbar(muh_pp[iens], value.(m_sh[iens]), err.(m_sh[iens]), fmt="x")
        errorbar(value(muh_target[iens]), value(m_sh_match[iens]), err(m_sh_match[iens]), err(muh_target[iens]), fmt="x")
        errorbar(muh_aa[iens], value.(m_sh_star[iens]), err.(m_sh_star[iens]), fmt="x")
        errorbar(value(muh_target[iens]), value(m_sh_v_match[iens]), err(m_sh_v_match[iens]), err(muh_target[iens]), fmt="x")   
        legend([L"$m_{D_s}$", L"$m_{D_s} (\mathrm{int})$", L"$m_{D^*_s}$", L"$m_{D^*_s} (\mathrm{int})$"])   
        display(gcf())
        t = string(ensembles[iens], "_mDs.pdf")
        savefig(joinpath(path_plot, t))
        close()
    end
    #f_lh f_sh
    figure()
    title(ensembles[iens])
    xlabel(L"$a\mu$")
    ylabel(L"$af$")
    errorbar(muh_pp[iens], value.(f_lh[iens]), err.(f_lh[iens]), fmt="x")
    errorbar(value(muh_target[iens]), value(f_lh_match[iens]), err(f_lh_match[iens]), err(muh_target[iens]), fmt="x")
    l = [L"$f_D$", L"$f_D (\mathrm{int})$"]
    if !deg[iens]
        errorbar(muh_pp[iens], value.(f_sh[iens]), err.(f_sh[iens]), fmt="x")
        errorbar(value(muh_target[iens]), value(f_sh_match[iens]), err(f_sh_match[iens]), err(muh_target[iens]), fmt="x")
        push!(l, L"$f_{D_s}$")
        push!(l, L"$f_{D_s} (\mathrm{int})$")
    end
    legend(l)
    display(gcf())
    t = string(ensembles[iens], "_fD.pdf")
    savefig(joinpath(path_plot, t))
    close()
 end
 ###########################
 ###### RESULTS ###########
 ###########################
 #Quark mass
 muc = zm_tm.(beta) .* muh_target .* sqrt.(8 * t0.(beta))
 uwerr.(muc)
 #Meson masses
 m_D = m_lh_match .* sqrt.(8 * t0.(beta))
 m_Ds = m_sh_match .* sqrt.(8 * t0.(beta))
 m_D_star = m_lh_v_match .* sqrt.(8 * t0.(beta))
 m_Ds_star = m_sh_v_match .* sqrt.(8 * t0.(beta))
 uwerr.(m_D)
 uwerr.(m_Ds)
 uwerr.(m_D_star)
 uwerr.(m_Ds_star)
 #Decay const
 f_D = f_lh_match .* sqrt.(8 * t0.(beta))
 f_Ds = f_sh_match .* sqrt.(8 * t0.(beta))
 uwerr.(f_D)
 uwerr.(f_Ds)
 for iens=1:length(ensembles)
    println("(", ensembles[iens], ")", L"$Z^{tm}_M \mu_c \sqrt{8t_0} = $", muc[iens])
    println("(", ensembles[iens], ")", L"$M_D \sqrt{8t_0}= $", m_D[iens])
    println("(", ensembles[iens], ")", L"$M_Ds \sqrt{8t_0}= $", m_Ds[iens])
    println("(", ensembles[iens], ")", L"$M_D* \sqrt{8t_0}= $", m_D_star[iens])
    println("(", ensembles[iens], ")", L"$M_Ds* \sqrt{8t_0}= $", m_Ds_star[iens])
    println("(", ensembles[iens], ")", L"$f_D \sqrt{8t_0}= $", f_D[iens])
    println("(", ensembles[iens], ")", L"$f_Ds \sqrt{8t_0}= $", f_Ds[iens])
 end
 ###########################
 ###### FITS ###############
 ###########################
 x = [1 ./(8 .* t0.(beta)) phi2] #x1 = a^2 / (8t0), x2 = 8t0 mpi^2
 uwerr.(x)
 phi2_ph = (t0_ph[1] * 139.57039 / hc)^2
 uwerr(phi2_ph)
 #f(a2/8t0,phi2) = p[1]+ p[2](a2/8t0) + (p[3] + p[4](a2/8t0)) * phi2 
 @. model(x, p) = (p[1] + p[2] * x[:, 1]) + (p[3] + p[4] * x[:, 1]) * x[:, 2] #linear fits
 obs = [muc, m_D, m_Ds, m_D_star, m_Ds_star, f_D, f_Ds] #sqrt(8t0) obs
 ttl_obs = ["muc", "m_D", "m_Ds", "m_D_star", "m_Ds_star", "f_D", "f_Ds"]
 ylbl = [L"$Z^{tm}_M \mu_c \sqrt{8t_0}$", L"$M_D \sqrt{8t_0}$", L"$M_{D_s} \sqrt{8t_0}$",
 L"$M_{D^*} \sqrt{8t_0}$", L"$M_{D^*_s} \sqrt{8t_0}$", L"$f_D \sqrt{8t_0}$", L"$f_{D_s} \sqrt{8t_0}$"]
 xlbl = L"$\phi_2$"
 obs_t0 = Vector{uwreal}(undef, length(obs)) #sqrt(8t0)obs @ CL & phi2_phys
 for k = 1:length(obs)
    println("OBS ",ttl_obs[k])
    par = fit_routine(model, value.(x), obs[k], 4)
    obs_t0[k] = par[1] + par[3] * phi2_ph
    uwerr(obs_t0[k])
    figure()
    for b in unique(beta) #select point with same beta
        nn = findall(x-> x == b, beta)
        lgnd = string(L"$\beta = $", b)
        errorbar(value.(x[nn,2]), value.(obs[k][nn]), err.(obs[k][nn]), err.(x[nn,2]), fmt="x", label=lgnd)
    end
    lgnd=L"$\mathrm{CL}$"
    errorbar(value(phi2_ph), value(obs_t0[k]), err(obs_t0[k]), err(phi2_ph), fmt="x", zorder=2, label=lgnd)
    axvline(value(phi2_ph), ls="--", color="black", zorder=1, lw=0.6, label="")
    xlabel(xlbl)
    ylabel(ylbl[k])
    legend()
    display(gcf())
    t = string("fit_", ttl_obs[k], ".pdf")
    savefig(joinpath(path_plot, 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
--- a/analysis/analysis_pcac.jl
+++ b/analysis/analysis_pcac.jl
 using juobs, ADerrors, DelimitedFiles, PyPlot, LaTeXStrings
 const path = "/home/javier/Lattice/charm/production_2"
 const path_plat = "/home/javier/Lattice/juobs/analysis/plat.txt"
 const ensembles = ["H400", "N202", "N200", "N203", "N300", "J303"]
 const deg = [true, true, false, false, true, false]
 const L = [32, 48, 48, 48, 48, 64]
 const beta = [3.46 , 3.55, 3.55, 3.55, 3.70, 3.70]
 const R = [["H400r001", "H400r002"], "N202r001", ["N200r000", "N200r001"], ["N203r000", "N203r001"], "N300r002", "J303r003"]
 include("/home/javier/Lattice/juobs/analysis/functions.jl")
 include("/home/javier/Lattice/juobs/constants/juobs_const.jl")
 m_ll = Vector{uwreal}(undef, length(ensembles))
 m_ss = Vector{uwreal}(undef, length(ensembles))
 m_hh = Vector{Vector{uwreal}}(undef, length(ensembles))
 mu_pp = Vector{Vector{Vector{Float64}}}(undef, length(ensembles))
 for iens = 1:length(ensembles)
    pp = read_dat(R[iens], "G5", "G5")
    a0p = read_dat(R[iens], "G5", "G0G5")
    rew = read_rew(R[iens])
    pp_obs = corr_obs.(pp, L=L[iens], rw=rew)
    a0p_obs = corr_obs.(a0p, L=L[iens], rw=rew)
    mu_pp[iens] = getfield.(pp_obs, :mu)
    m = comp_pcac(a0p_obs, pp_obs, deg[iens], ensembles[iens])
    m_ll[iens] = get_ll(mu_pp[iens], m, deg[iens])
    m_ss[iens] = deg[iens] ? m_ll[iens] : get_ss(mu_pp[iens], m, deg[iens])
    m_hh[iens] = get_hh(mu_pp[iens], m, deg[iens])
 end
 mm = get_mu.(mu_pp, deg)
 mul_pp = getindex.(mm, 1)
 mus_pp = getindex.(mm, 2)
 muh_pp = getindex.(mm, 3)
 cot_ll = Vector{uwreal}(undef, length(ensembles))
 cot_ss = Vector{uwreal}(undef, length(ensembles))
 cot_hh = Vector{Vector{uwreal}}(undef, length(ensembles))
 for iens = 1:length(ensembles)
    cot_ll[iens] = za(beta[iens]) * m_ll[iens] / mul_pp[iens]
    cot_ss[iens] = za(beta[iens]) * m_ss[iens] / mus_pp[iens]
    cot_hh[iens] = fill(za(beta[iens]), 3) .* m_hh[iens] ./ muh_pp[iens]
 end
\ No newline at end of file
--- a/analysis/fits.log
+++ b/analysis/fits.log
 OBS muc
 Fit parameter: 1: 3.3279451286618453 +/- 0.1872060642719553
 ## Number of error sources: 21
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  74.20   216
  #                       sqrt(8 t0) (fm) 00000001  15.37            -
  #                                            200   3.99   350,856
  #                                            400   1.19   505
  #                          ZM_tm values 00000004   1.07            -
  #                                            300   0.94   245
  #                                            203   0.91   756,622
  #                                    t0 00000003   0.76            -
  #                                    t0 00000004   0.73            -
  #                          ZM_tm values 00000003   0.71            -
  #                                    t0 00000002   0.09            -
  #                          ZM_tm values 00000002   0.03            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                          ZM_tm values 00000001   0.00            -
  #                          ZM_tm values 00000005   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Fit parameter: 2: -7.701331946021094 +/- 7.646190205835605
 ## Number of error sources: 21
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  75.70   216
  #                                            200  11.97   350,856
  #                                            400   3.22   505
  #                                            203   2.18   756,622
  #                                    t0 00000003   2.12            -
  #                          ZM_tm values 00000003   1.96            -
  #                          ZM_tm values 00000004   1.46            -
  #                                    t0 00000004   1.03            -
  #                                    t0 00000002   0.26            -
  #                          ZM_tm values 00000002   0.08            -
  #                                            300   0.02   245
  #                       sqrt(8 t0) (fm) 00000001   0.01            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                          ZM_tm values 00000001   0.00            -
  #                          ZM_tm values 00000005   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Fit parameter: 3: -0.4540918027823863 +/- 0.2388566093430024
 ## Number of error sources: 21
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  75.19   216
  #                                            300  13.16   245
  #                                            400   4.41   505
  #                                            200   3.62   350,856
  #                                            203   1.17   756,622
  #                                    t0 00000003   0.80            -
  #                          ZM_tm values 00000003   0.74            -
  #                                    t0 00000002   0.35            -
  #                          ZM_tm values 00000004   0.29            -
  #                                    t0 00000004   0.15            -
  #                          ZM_tm values 00000002   0.11            -
  #                       sqrt(8 t0) (fm) 00000001   0.01            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                          ZM_tm values 00000001   0.00            -
  #                          ZM_tm values 00000005   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Fit parameter: 4: 12.11810047771023 +/- 10.384037332178496
 ## Number of error sources: 21
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  67.54   216
  #                                            200  10.58   350,856
  #                                            400  10.38   505
  #                                            300   3.78   245
  #                                            203   1.95   756,622
  #                                    t0 00000003   1.88            -
  #                          ZM_tm values 00000003   1.74            -
  #                                    t0 00000002   0.83            -
  #                          ZM_tm values 00000004   0.63            -
  #                                    t0 00000004   0.40            -
  #                          ZM_tm values 00000002   0.26            -
  #                       sqrt(8 t0) (fm) 00000001   0.04            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                          ZM_tm values 00000001   0.00            -
  #                          ZM_tm values 00000005   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Chisq / chiexp: 0.037007412036443685 (Error not available... maybe run uwerr) / 0.050909328579011835 (dof: 1)
 OBS m_D
 Fit parameter: 1: 4.048595157522108 +/- 0.13295625736285707
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  74.71   216
  #                       sqrt(8 t0) (fm) 00000001  17.76            -
  #                                            200   4.01   350,856
  #                                            400   1.38   505
  #                                            300   1.05   245
  #                                            203   0.96   756,622
  #                                    t0 00000004   0.06            -
  #                                    t0 00000003   0.06            -
  #                                    t0 00000002   0.01            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Fit parameter: 2: -4.367427212140693 +/- 5.281775769423441
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  80.57   216
  #                                            200  12.73   350,856
  #                                            400   3.95   505
  #                                            203   2.43   756,622
  #                                    t0 00000003   0.17            -
  #                                    t0 00000004   0.09            -
  #                                            300   0.02   245
  #                       sqrt(8 t0) (fm) 00000001   0.02            -
  #                                    t0 00000002   0.02            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Fit parameter: 3: -0.16797828889715144 +/- 0.1700443403112658
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  75.35   216
  #                                            300  14.61   245
  #                                            400   5.10   505
  #                                            200   3.60   350,856
  #                                            203   1.25   756,622
  #                                    t0 00000003   0.06            -
  #                                    t0 00000002   0.02            -
  #                                    t0 00000004   0.01            -
  #                       sqrt(8 t0) (fm) 00000001   0.00            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Fit parameter: 4: 8.006344213712698 +/- 7.269508682875713
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  69.99   216
  #                                            400  12.40   505
  #                                            200  10.91   350,856
  #                                            300   4.32   245
  #                                            203   2.15   756,622
  #                                    t0 00000003   0.14            -
  #                                    t0 00000002   0.05            -
  #                                    t0 00000004   0.03            -
  #                       sqrt(8 t0) (fm) 00000001   0.00            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Chisq / chiexp: 0.01494132594664716 (Error not available... maybe run uwerr) / 0.046712860766424225 (dof: 1)
 OBS m_Ds
 Fit parameter: 1: 4.31362022727613 +/- 0.13685549095602995
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  76.09   216
  #                       sqrt(8 t0) (fm) 00000001  16.72            -
  #                                            200   3.95   350,856
  #                                            400   1.30   505
  #                                            300   1.00   245
  #                                            203   0.89   756,622
  #                                    t0 00000003   0.03            -
  #                                    t0 00000004   0.02            -
  #                                    t0 00000002   0.01            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Fit parameter: 2: -6.726552717896041 +/- 5.457720362216015
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  81.42   216
  #                                            200  12.45   350,856
  #                                            400   3.69   505
  #                                            203   2.23   756,622
  #                                    t0 00000003   0.09            -
  #                       sqrt(8 t0) (fm) 00000001   0.05            -
  #                                    t0 00000004   0.04            -
  #                                            300   0.02   245
  #                                    t0 00000002   0.01            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Fit parameter: 3: -0.5111388042444888 +/- 0.17512221322135785
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  76.66   216
  #                                            300  13.78   245
  #                                            400   4.80   505
  #                                            200   3.54   350,856
  #                                            203   1.16   756,622
  #                                    t0 00000003   0.03            -
  #                                    t0 00000002   0.02            -
  #                                    t0 00000004   0.00            -
  #                       sqrt(8 t0) (fm) 00000001   0.00            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Fit parameter: 4: 11.125384531033184 +/- 7.4804808444962925
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  71.33   216
  #                                            400  11.68   505
  #                                            200  10.75   350,856
  #                                            300   4.10   245
  #                                            203   1.99   756,622
  #                                    t0 00000003   0.08            -
  #                                    t0 00000002   0.05            -
  #                       sqrt(8 t0) (fm) 00000001   0.02            -
  #                                    t0 00000004   0.01            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Chisq / chiexp: 0.0034879355163730487 (Error not available... maybe run uwerr) / 0.04783649734148854 (dof: 1)
 OBS m_D_star
 Fit parameter: 1: 4.078728153874206 +/- 0.11043414144470609
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  76.50   216
  #                       sqrt(8 t0) (fm) 00000001  18.62            -
  #                                            200   3.87   350,856
  #                                            203   0.62   756,622
  #                                            400   0.22   505
  #                                            300   0.17   245
  #                                    t0 00000003   0.00            -
  #                                    t0 00000004   0.00            -
  #                                    t0 00000002   0.00            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Fit parameter: 2: 4.546752906721969 +/- 4.320814627740389
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  85.04   216
  #                                            200  12.57   350,856
  #                                            203   1.66   756,622
  #                                            400   0.66   505
  #                       sqrt(8 t0) (fm) 00000001   0.05            -
  #                                    t0 00000003   0.01            -
  #                                    t0 00000004   0.01            -
  #                                            300   0.00   245
  #                                    t0 00000002   0.00            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Fit parameter: 3: 0.28567507471029296 +/- 0.1299619002885983
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  91.13   216
  #                                            200   4.18   350,856
  #                                            300   2.79   245
  #                                            400   0.96   505
  #                                            203   0.91   756,622
  #                       sqrt(8 t0) (fm) 00000001   0.02            -
  #                                    t0 00000003   0.01            -
  #                                    t0 00000002   0.00            -
  #                                    t0 00000004   0.00            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Fit parameter: 4: -6.680949348261337 +/- 5.609252379522268
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  83.04   216
  #                                            200  12.16   350,856
  #                                            400   2.32   505
  #                                            203   1.62   756,622
  #                                            300   0.80   245
  #                       sqrt(8 t0) (fm) 00000001   0.03            -
  #                                    t0 00000003   0.01            -
  #                                    t0 00000002   0.01            -
  #                                    t0 00000004   0.00            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Chisq / chiexp: 2.4765381445760555e-5 (Error not available... maybe run uwerr) / 0.026775165015589318 (dof: 1)
 OBS m_Ds_star
 Fit parameter: 1: 4.523261012567903 +/- 0.10344559001693072
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  70.42   216
  #                       sqrt(8 t0) (fm) 00000001  25.59            -
  #                                            200   3.18   350,856
  #                                            203   0.32   756,622
  #                                            400   0.25   505
  #                                            300   0.19   245
  #                                    t0 00000004   0.02            -
  #                                    t0 00000003   0.02            -
  #                                    t0 00000002   0.00            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Fit parameter: 2: -3.9473716737053306 +/- 3.848265028792789
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  86.60   216
  #                                            200  11.45   350,856
  #                                            203   0.95   756,622
  #                                            400   0.83   505
  #                       sqrt(8 t0) (fm) 00000001   0.07            -
  #                                    t0 00000003   0.06            -
  #                                    t0 00000004   0.03            -
  #                                            300   0.00   245
  #                                    t0 00000002   0.00            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Fit parameter: 3: -0.2888313108867589 +/- 0.11690482054266071
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  90.96   216
  #                                            200   3.72   350,856
  #                                            300   3.45   245
  #                                            400   1.19   505
  #                                            203   0.52   756,622
  #                       sqrt(8 t0) (fm) 00000001   0.12            -
  #                                    t0 00000003   0.02            -
  #                                    t0 00000004   0.01            -
  #                                    t0 00000002   0.00            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Fit parameter: 4: 4.3198464617638015 +/- 5.012696506358427
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  83.99   216
  #                                            200  11.01   350,856
  #                                            400   2.90   505
  #                                            300   1.00   245
  #                                            203   0.92   756,622
  #                       sqrt(8 t0) (fm) 00000001   0.09            -
  #                                    t0 00000003   0.06            -
  #                                    t0 00000004   0.02            -
  #                                    t0 00000002   0.01            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Chisq / chiexp: 0.015177074710628055 (Error not available... maybe run uwerr) / 0.016118934044216784 (dof: 1)
 OBS f_D
 Fit parameter: 1: 0.429406512463374 +/- 0.02014695034107166
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  89.06   216
  #                                            200   4.46   350,856
  #                                            203   1.80   756,622
  #                                    t0 00000004   1.28            -
  #                                            300   1.26   245
  #                                    t0 00000003   1.12            -
  #                                            400   0.92   505
  #                                    t0 00000002   0.11            -
  #                       sqrt(8 t0) (fm) 00000001   0.00            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Fit parameter: 2: 0.8840051508330702 +/- 0.8877169695716438
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  78.08   216
  #                                            200  12.01   350,856
  #                                            203   3.20   756,622
  #                                    t0 00000003   2.59            -
  #                                            400   2.07   505
  #                                    t0 00000004   1.51            -
  #                       sqrt(8 t0) (fm) 00000001   0.26            -
  #                                    t0 00000002   0.25            -
  #                                            300   0.04   245
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Fit parameter: 3: 0.08110781249666073 +/- 0.028073841349955025
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  75.65   216
  #                                            300  14.09   245
  #                                            400   3.02   505
  #                                            200   2.99   350,856
  #                                            203   2.56   756,622
  #                                    t0 00000003   1.03            -
  #                                    t0 00000002   0.36            -
  #                                    t0 00000004   0.25            -
  #                       sqrt(8 t0) (fm) 00000001   0.05            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Fit parameter: 4: -0.23531500117442072 +/- 1.194580435504181
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  70.95   216
  #                                            200  10.38   350,856
  #                                            400   7.06   505
  #                                            300   4.41   245
  #                                            203   3.34   756,622
  #                                    t0 00000003   2.39            -
  #                                    t0 00000002   0.84            -
  #                                    t0 00000004   0.63            -
  #                       sqrt(8 t0) (fm) 00000001   0.00            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Chisq / chiexp: 0.47283542288170355 (Error not available... maybe run uwerr) / 0.7452145868659791 (dof: 1)
 OBS f_Ds
 Fit parameter: 1: 0.5306749798655852 +/- 0.012769894088102694
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  81.45   216
  #                                    t0 00000004   3.84            -
  #                                            200   3.75   350,856
  #                                            300   3.28   245
  #                                    t0 00000003   3.19            -
  #                                            400   2.22   505
  #                                            203   1.73   756,622
  #                       sqrt(8 t0) (fm) 00000001   0.27            -
  #                                    t0 00000002   0.27            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Fit parameter: 2: 0.25158377584140645 +/- 0.5678900623551665
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  70.12   216
  #                                            200  10.37   350,856
  #                                    t0 00000003   7.18            -
  #                                            400   4.68   505
  #                                    t0 00000004   4.27            -
  #                                            203   2.47   756,622
  #                                    t0 00000002   0.56            -
  #                                            300   0.19   245
  #                       sqrt(8 t0) (fm) 00000001   0.17            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Fit parameter: 3: -0.05116934552402741 +/- 0.01962706259221214
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  56.86   216
  #                                            300  28.30   245
  #                                            400   6.33   505
  #                                            203   2.63   756,622
  #                                    t0 00000003   2.49            -
  #                                            200   1.84   350,856
  #                                    t0 00000004   0.78            -
  #                                    t0 00000002   0.76            -
  #                       sqrt(8 t0) (fm) 00000001   0.02            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Fit parameter: 4: 0.6338434751229951 +/- 0.8165786251257829
 ## Number of error sources: 16
 ## Number of MC ids       : 5
 ## Contribution to error  :               Ensemble  [%]     [MC length]
  #                                            303  55.80   216
  #                                            400  14.68   505
  #                                            300   9.99   245
  #                                            200   7.69   350,856
  #                                    t0 00000003   5.83            -
  #                                            203   2.51   756,622
  #                                    t0 00000002   1.76            -
  #                                    t0 00000004   1.67            -
  #                       sqrt(8 t0) (fm) 00000001   0.07            -
  #                  charmed meson masses 00000004   0.00            -
  #                  charmed meson masses 00000002   0.00            -
  #                  charmed meson masses 00000001   0.00            -
  #                  charmed meson masses 00000003   0.00            -
  #                                    t0 00000001   0.00            -
  #                  charmed meson masses 00000005   0.00            -
  #                  charmed meson masses 00000006   0.00            -
 Chisq / chiexp: 1.8860007864688775 (Error not available... maybe run uwerr) / 0.5314364650138338 (dof: 1)
--- a/analysis/functions.jl
+++ b/analysis/functions.jl
 function read_dat(rep::String, g1::String="G5", g2::String="G5")
    p = joinpath(path, rep, "info")
    f = filter(x-> contains(x, ".dat"), readdir(p))
    f = joinpath(p, f[1])
    return read_mesons(f, g1, g2)
 end
 function read_dat(rep::Vector{String}, g1::String="G5", g2::String="G5")
    p = joinpath.(path, rep, "info")
    f = [filter(x-> contains(x, ".dat"), readdir(p[k]))[1] for k = 1:length(p)] 
    f = joinpath.(p, f)
    return read_mesons(f, g1, g2)
 end
 function read_rew(rep::String)
    p = joinpath(path, rep, "rew")
    f = filter(x-> contains(x, ".dat"), readdir(p))
    f = joinpath(p, f[1])
    try
        return read_ms1(f)
    catch
        return read_ms1(f, v="1.4")
    end
 end
 read_rew(rep::Vector{String}) = read_rew.(rep)
 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 && mu[1] != mu[2] #l-s
            return obs[k]
        end
    end
 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_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_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]
        if mu[1] == mu[2] && !(mul in mu) && !(mus in mu)
            push!(obs_hh, obs[k])
        end
    end
    return obs_hh
 end
 function select_plateau(ens::String, mu_list, deg::Bool)
    mul, mus, muh = get_mu(mu_list, deg)
    f = readdlm(path_plat)
    head = String.(f[:, 1])
    delim = findall(x-> contains(x, "#"), head)
    edelim = findfirst(x-> contains(x, ens), head)
    cdelim = findfirst(x-> x.== edelim, delim)
    del = delim[cdelim]+1 : delim[cdelim+1]-1
    plat = Vector{Vector{Int64}}(undef, 0)
    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-> contains(x, "lh"), head_)
            push!(plat, [plat_[n, 1], plat_[n, 2]])
        elseif mul in mu && mu[1] == mu[2] #light-light
            n = findfirst(x-> contains(x, "ll"), head_)
            push!(plat, [plat_[n, 1], plat_[n, 2]])
        elseif mul in mu && mus in mu#strange-light
            n = findfirst(x-> contains(x, "ls"), 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-> contains(x, "sh"), 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-> contains(x, "ss"), head_)
            push!(plat, [plat_[n, 1], plat_[n, 2]])
        elseif !(mul in mu) && !(mus in mu) #heavy-heavy
            n = findfirst(x-> contains(x, "hh"), head_)
            push!(plat, [plat_[n, 1], plat_[n, 2]])
        end
    end
    return plat
 end
 function comp_meff(pp_obs::Vector{juobs.Corr}, deg::Bool, ens::String; pl::Bool=false)
    mu_list = getfield.(pp_obs, :mu)
    plat = select_plateau(ens, mu_list, deg)
    return meff.(pp_obs, plat, pl=pl)
 end
 function comp_pcac(a0p_obs::Vector{juobs.Corr}, pp_obs::Vector{juobs.Corr}, deg::Bool, ens::String; pl::Bool=false)
    mu_list = getfield.(pp_obs, :mu)
    plat = select_plateau(ens, mu_list, deg)
    return mpcac.(a0p_obs, pp_obs, plat, pl=pl)
 end
 function comp_f(pp_obs::Vector{juobs.Corr}, m::Vector{uwreal}, deg::Bool, ens::String; pl::Bool=false)
    mu_list = getfield.(pp_obs, :mu)
    plat = select_plateau(ens, mu_list, deg)
    return dec_const_pcvc.(pp_obs, plat, m, pl=pl)
 end
 function match_muc(muh, m_lh, m_lh_star, target)
    M = (m_lh .+ 3 .* m_lh_star) ./ 4
    par, chi2exp = lin_fit(muh, M)
    muh_target = x_lin_fit(par, target)
    return muh_target
 end
 function match_muc(muh, m_lh, m_sh, m_lh_star, m_sh_star, target)
    M = (2 .* m_lh .+ 6 .* m_lh_star .+ m_sh .+ 3 .* m_sh_star) ./ 12
    par, chi2exp = lin_fit(muh, M)
    muh_target = x_lin_fit(par, target)
    return muh_target
 end
\ No newline at end of file
--- a/analysis/results.txt
+++ b/analysis/results.txt
 (H400)$Z^{tm}_M \mu_c \sqrt{8t_0} = $3.030909671357194 +/- 0.07903701818131038
 (H400)$M_D \sqrt{8t_0}= $3.9807556606105736 +/- 0.05781703939545818
 (H400)$M_Ds \sqrt{8t_0}= $3.9807556606105736 +/- 0.05781703939545818
 (H400)$M_D* \sqrt{8t_0}= $4.278934408455874 +/- 0.05001764960908615
 (H400)$M_Ds* \sqrt{8t_0}= $4.278934408455874 +/- 0.05001764960908615
 (H400)$f_D \sqrt{8t_0}= $0.5171340154200049 +/- 0.0030551745985989815
 (H400)$f_Ds \sqrt{8t_0}= $0.5171340154200049 +/- 0.0030551745985989815
 (N200)$Z^{tm}_M \mu_c \sqrt{8t_0} = $3.0777083812680863 +/- 0.07529964071401364
 (N200)$M_D \sqrt{8t_0}= $3.9484624535288577 +/- 0.056300551744400275
 (N200)$M_Ds \sqrt{8t_0}= $4.0641686659916445 +/- 0.05571406258479428
 (N200)$M_D* \sqrt{8t_0}= $4.2322082817565505 +/- 0.05089413966187937
 (N200)$M_Ds* \sqrt{8t_0}= $4.366150683510803 +/- 0.050567123757276136
 (N200)$f_D \sqrt{8t_0}= $0.4781332278412574 +/- 0.002533317431472261
 (N200)$f_Ds \sqrt{8t_0}= $0.5253557430045498 +/- 0.0019233408904497612
 (N203)$Z^{tm}_M \mu_c \sqrt{8t_0} = $3.069739432147581 +/- 0.07467450492462731
 (N203)$M_D \sqrt{8t_0}= $3.9619597653172356 +/- 0.05581411314938938
 (N203)$M_Ds \sqrt{8t_0}= $4.027646770204338 +/- 0.055472228755110334
 (N203)$M_D* \sqrt{8t_0}= $4.253505640330456 +/- 0.050187214812839276
 (N203)$M_Ds* \sqrt{8t_0}= $4.326727158575913 +/- 0.049836834652424916
 (N203)$f_D \sqrt{8t_0}= $0.488543655769974 +/- 0.0026562881901046583
 (N203)$f_Ds \sqrt{8t_0}= $0.5157212901463102 +/- 0.0021384519222027373
 (N300)$Z^{tm}_M \mu_c \sqrt{8t_0} = $2.997384760490272 +/- 0.07895636913539895
 (N300)$M_D \sqrt{8t_0}= $3.94373962398667 +/- 0.059571174171691506
 (N300)$M_Ds \sqrt{8t_0}= $3.94373962398667 +/- 0.059571174171691506
 (N300)$M_D* \sqrt{8t_0}= $4.291279319094004 +/- 0.04988459701774591
 (N300)$M_Ds* \sqrt{8t_0}= $4.291279319094004 +/- 0.04988459701774591
 (N300)$f_D \sqrt{8t_0}= $0.5032296543520727 +/- 0.003757345679706973
 (N300)$f_Ds \sqrt{8t_0}= $0.5032296543520727 +/- 0.003757345679706973
 (J303)$Z^{tm}_M \mu_c \sqrt{8t_0} = $3.13490498080823 +/- 0.08328602171541571
 (J303)$M_D \sqrt{8t_0}= $3.970055446897706 +/- 0.062204612871266905
 (J303)$M_Ds \sqrt{8t_0}= $4.113876464562304 +/- 0.06232345158893755
 (J303)$M_D* \sqrt{8t_0}= $4.199854539735192 +/- 0.05475423487017079
 (J303)$M_Ds* \sqrt{8t_0}= $4.399899824892723 +/- 0.05477252163421222
 (J303)$f_D \sqrt{8t_0}= $0.4649182529771295 +/- 0.004855221613601069
 (J303)$f_Ds \sqrt{8t_0}= $0.5220832840476 +/- 0.0031391370517194982
 $Z^{tm}_M \mu_c \sqrt{8t_0}$ = 3.289196404345117 +/- 0.16908049991604063
 muc(MeV) = 1571.5428436124926 +/- 74.65754166879076
 $M_D \sqrt{8t_0}$ = 4.034261174024082 +/- 0.12045227374006602
 m_D(MeV) = 1927.5268174700511 +/- 51.13372706582599
 $M_Ds \sqrt{8t_0}$ = 4.270003547803183 +/- 0.12354622671856343
 m_Ds(MeV) = 2040.1620009329038 +/- 52.898332267818304
 $M_D* \sqrt{8t_0}$ = 4.103105481644236 +/- 0.10100061920706542
 m_D_star(MeV) = 1960.4198909335223 +/- 42.35821470980299
 $M_Ds* \sqrt{8t_0}$ = 4.498614355722275 +/- 0.09456340196942646
 m_Ds_star(MeV) = 2149.3897985442322 +/- 37.956898903502925
 $f_D \sqrt{8t_0}$ = 0.4363276333625155 +/- 0.017871060655147757
 f_D(MeV) = 208.47267398669123 +/- 8.883408580389085
 $f_Ds \sqrt{8t_0}$ = 0.5263085789505523 +/- 0.011256710293275131
 f_Ds(MeV) = 251.46460688360165 +/- 6.054227498288841
--- a/constants/juobs_const.jl
+++ b/constants/juobs_const.jl
 const hc = 197.3269804 #MeV fm
 const ens_db = Dict(
    #"ens_id"=>[L, T, beta, dtr, vrw, trunc_cnfg]
    "H101"     => [32, 96, 3.4, 2, "1.2", [1001,1009]],
    "H102r001" => [32, 96, 3.4, 2, "2.0", [997]],
    "H102r002" => [32, 96, 3.4, 2, "2.0", [1008]],
    "H105"     => [32, 96, 3.4, 2, "2.0", [947,1042]],
    "H105r005" => [32, 96, 3.4, 1, "1.2", [837]],
    "H400"     => [32, 96, 3.46, 1, "1.2", [505,540]],
    "D450"     => [64, 128, 3.46, 1, ["2.0"], [1000]],
    "N202"     => [48, 128, 3.55, 2, "1.2", [899]],
    "N203"     => [48, 128, 3.55, 1, "2.0", [756,787]],
    "N200"     => [48, 128, 3.55, 1, "2.0", [856,856]],
    "D200"     => [64, 128, 3.55, 2, "2.0", [2001]],
    "E250"     => [96, 192, 3.55, 1, ["2.0"], [1009]],
    "N300"     => [48, 128, 3.70, 1, "1.2", [1521]],
    "N302"     => [48, 128, 3.70, 1, "1.2", [2201]],
    "J303"     => [64, 192, 3.70, 2, "2.0", [1073]],
    "E300"     => [96, 192, 3.70, 1, ["1.4"], [1139]],
    "J500"     => [64, 192, 3.85, 2, ["1.2", "1.4", "2.0"], [789,655,431]],
    "J501"     => [64, 192, 3.85, 1, ["1.2", "1.4", "2.0"], [1635,1142,1150]]
 )
 const db = Dict(
    "J501" => ["ts001_chunk1", "ts001_chunk3", "ts190_chunk1", "ts190_chunk2"],
    "J500" => ["ts001_chunk1", "ts001_chunk2", "ts190_chunk1", "ts190_chunk2"],
    "E250" => ["ts001", "ts097"],
    "E300" => ["ts001_chunk2", "ts001_chunk3", "ts190_chunk1"],#, "ts001_chunk1"],
    "D450" => ["ts001", "ts065_1", "ts065_2"]
 )
 const db_c = Dict(
    "J501" => ["ts001", "ts190"],
    "J500" => ["ts001_chunk1", "ts190_chunk2"],
    "E250" => ["ts001"],
    "E300" => ["ts001_chunk2"],#, "ts001_chunk1"]
    "D450" => ["ts001_1", "ts001_2"]
 )
 const sym_bool = Dict(
    "J501" => true,
    "J500" => true,
    "E250" => false,
    "E300" => false,
    "D450" => false
 )
 const flag_s = Dict(
    "J303r003" => [324,325,326],
    "H105r001" => [100,105,106],
    "H105r002" => [1],
    "H105r005" => [254, 255, 256, 257, 259, 260, 261, 264, 265, 266, 269, 280, 282, 283, 284, 285, 286, 287, 288, 289, 291, 299, 301, 313, 314, 315, 316, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341, 342]
 )
 const ensemble_inv = Dict(
    "H101"     => 1,
    "H102r001" => 2,
    "H102r002" => 3,
    "H105"     => 4,
    "H105r005" => 5,
    "H400"     => 6,
    "D450"     => 7,
    "N202"     => 8,
    "N203"     => 9,
    "N200"     => 10,
    "D200"     => 11,
    "E250"     => 12,
    "N300"     => 13,
    "N302"     => 14,
    "J303"     => 15,
    "E300"     => 16,
    "J500"     => 17,
    "J501"     => 18
 )
\ No newline at end of file
--- a/constants/juobs_uwreal_const.jl
+++ b/constants/juobs_uwreal_const.jl
 import ADerrors
-const b_values = [3.40, 3.46, 3.55, 3.70, 3.85]
+const beta_to_idx = Dict{Float64,Int64}(3.40 =>1, 3.46=>2, 3.55=>3, 3.70=>4, 3.85=>5)
-const M = let C = zeros(8, 8)
+let 
-#PDG
+  C = zeros(8, 8)
  M_values = [1869.65, 2010.26, 1968.34, 2112.2, 2980.3, 3096.916, 5279.65, 5366.3] #MD, MD*, MDs, MDs*, \eta_c, J/\psi , MB, MBs(MeV)
  M_error = [0.05, 0.05, 0.07, 0.4, 1.2, 0.011, 0.12, 0.6]
-  for i = 1:8
+  M = Ref{Vector{uwreal}}()
-    C[i, i] = M_error[i] ^ 2
+  for i = 1:8
-  end
+    C[i, i] = M_error[i] ^ 2
-  ADerrors.cobs(M_values, C, "charmed meson masses")
+  end
-end
+  global function M()
+    if !isassigned(M) 
-#1802.05243 taking into account correlations in zm_tm
+      M.x = ADerrors.cobs(M_values, C, "charmed meson masses")
-let C = [[0.375369e-6, 0.429197e-6, -0.186896e-5] [0.429197e-6, 0.268393e-4, 0.686776e-4] [-0.186896e-5, 0.686776e-4, 0.212386e-3]]
+    end
+    return M.x;
-  z =  ADerrors.cobs([0.348629, 0.020921, 0.070613], C, "z")
+  end
+end
-  global ZP(b) = z[1] + z[2] * (b - 3.79) + z[3] * (b - 3.79)^2
-  CC = [[0.164635e-4, 0.215658e-4, -0.754203e-4] [0.215658e-4, 0.121072e-2, 0.308890e-2] [-0.754203e-4, 0.308890e-2, 0.953843e-2]]
+let 
-  zz = ADerrors.cobs([2.270073, 0.121644, -0.464575], CC, "zm")
+  #1802.05243 taking into account correlations in zm_tm
+  C = [0.375369e-6, 0.429197e-6, -0.186896e-5; 0.429197e-6, 0.268393e-4, 0.686776e-4; -0.186896e-5, 0.686776e-4, 0.212386e-3]
-  # same for zm wilson
+  z = Ref{Vector{uwreal}}();
-  Mrat = ADerrors.uwreal([.9148, 0.0088], "M/mhad")
+  global function ZP(beta::Float64)
+    if !isassigned(z)
-  global z_covar(b) = Mrat.mean*(zz[1] + zz[2]*(b-3.79) + zz[3]*(b-3.79)^2)
+      z.x =  ADerrors.cobs([0.348629, 0.020921, 0.070613], C, "z")
-  global zm_tm_covar(b) = Mrat / ZP(b)
+    end
-end 
+    return z.x[1] + z.x[2] *(beta-3.79) + z.x[3]*(beta-3.79)^2
+  end
-## taking into account correlations in r_m
+  CC = [0.164635e-4, 0.215658e-4, -0.754203e-4; 0.215658e-4, 0.121072e-2, 0.308890e-2; -0.754203e-4, 0.308890e-2, 0.953843e-2]
-#Crm = [[3.699760, -2.198586, -1.476913e-3] [-2.198586, 1.306512, 8.776569e-4] [-1.476913e-3, 8.776569e-4, 5.895710e-7] ]
+  zm = Ref{Vector{uwreal}}()
-#c_i = ADerrors.cobs([-7.86078, 5.49175, -0.54078], Crm, "c_i in rm")
+  Mrat_ = Ref{uwreal}()
-#rm(b::Float64) = 1.0 +  0.004630*(6/b)^4 * ((1. +c_i[1]*(6/b)^2 + c_i[2]*(6/b)^4) / (1. + c_i[3]*(6/b)^2))
+  global function Mrat_()
-let C = zeros(5, 5)
+    if !isassigned(Mrat_)
-  #1802.05243 
+      Mrat_.x =ADerrors.uwreal([.9148, 0.0088], "M/mhad")
-  ZM_error = [35, 27, 33, 38, 45] .* 1e-4
+    end
-  ZM_data = [1.9684, 1.9935, 2.0253, 2.0630, 2.0814]
+    return Mrat_.x
+  end
-  for i in 1:5
-    C[i, i] = ZM_error[i] ^ 2
+  global function z_covar(beta::Float64)
-  end
+    if !isassigned(zz)
-  ZM = ADerrors.cobs(ZM_data, C, "ZM values")
+      zm.x = ADerrors.cobs([2.270073, 0.121644, -0.464575], CC, "zm")
-  global zm(beta::Float64)         = ZM[b_values .== beta][1]
+    end
-end
+    return Mrat().mean*(zz.x[1] + zz.x[2]*(beta-3.79)+zz.x[3]*(b-3.79)^2)
+  end
-let C = zeros(5, 5)
-  #1802.05243 
+  global  zm_tm_covar(beta::Float64) = Mrat() / ZP(beta)
-  ZM_tm_data = [2.6047, 2.6181, 2.6312, 2.6339, 2.6127]
+end
-  ZM_tm_error = [42, 33, 42, 48, 55] .* 1e-4
-  for i in 1:5
+## taking into account correlations in r_m
-    C[i, i] = ZM_tm_error[i] ^ 2
+#Crm = [[3.699760, -2.198586, -1.476913e-3] [-2.198586, 1.306512, 8.776569e-4] [-1.476913e-3, 8.776569e-4, 5.895710e-7] ]
-  end
+#c_i = ADerrors.cobs([-7.86078, 5.49175, -0.54078], Crm, "c_i in rm")
-  ZM_tm = ADerrors.cobs(ZM_tm_data, C, "ZM_tm values")
+#rm(b::Float64) = 1.0 +  0.004630*(6/b)^4 * ((1. +c_i[1]*(6/b)^2 + c_i[2]*(6/b)^4) / (1. + c_i[3]*(6/b)^2))
-  global zm_tm(beta::Float64)      = ZM_tm[b_values .== beta][1]
-end
+let C = zeros(5, 5)
+  #1802.05243 
-let C = zeros(5, 5)
+  ZM_error = [35, 27, 33, 38, 45] .* 1e-4
-  #2101.02694
+  ZM_data = [1.9684, 1.9935, 2.0253, 2.0630, 2.0814]
-  t0_data = [2.846, 3.634, 5.140, 8.564, 14.064]
-  t0_error = [8, 13, 21, 24, 63] .* 1e-3
+  for i in 1:5
+    C[i, i] = ZM_error[i] ^ 2
+  end
-  for i in 1:5
+  ZM = Ref{Vector{uwreal}}()
-    C[i, i] = t0_error[i] ^ 2
+  global function zm(beta::Float64)
-  end
+    if !isassigned(ZM)
+      ZM.x =  ADerrors.cobs(ZM_data, C, "ZM values")
-  t0_ = ADerrors.cobs(t0_data, C, "t0")
+    end
-  global t0(beta::Float64)         = t0_[b_values .== beta][1]
+    return ZM[beta_to_idx[beta]]
+  end
-  # from AS analysis
+end
-  t0_ph =ADerrors.uwreal([0.4118,0.25e-4],"sqrt(8 t0) (fm)") 
-  a_ = [t0_ph] ./ sqrt.(8 .* t0_)
+let C = zeros(5, 5)
-  global a(beta::Float64)          = a_[b_values .== beta][1]
+  #1802.05243 
-end
+  ZM_tm_data = [2.6047, 2.6181, 2.6312, 2.6339, 2.6127]
+  ZM_tm_error = [42, 33, 42, 48, 55] .* 1e-4
-let C = zeros(5, 5)
+  for i in 1:5
-  #1808.09236
+    C[i, i] = ZM_tm_error[i] ^ 2
-  ZA_data = [0.75642, 0.76169, 0.76979, 0.78378, 0.79667]
+  end
-  ZA_err = [72, 93, 43, 47, 47] .*1e-5
+  ZM_tm = Ref{Vector{uwreal}}()
-  for i in 1:5
+  global function zm_tm(beta::Float64)
-    C[i, i] = ZA_err[i] ^ 2
+    if !isassigned(ZM_tm)
-  end
+      ZM_tm.x = ADerrors.cobs(ZM_tm_data, C, "ZM_tm values")
+    end
-  Za = ADerrors.cobs(ZA_data, C, "ZA")
+    return ZM_tm.x[beta_to_idx[beta]]
-  global za(beta::Float64)         = Za[b_values .== beta][1]
+  end
 end
 let C = zeros(5, 5)
-  #1808.09236
+  #2101.02694
-  ZV_data = [0.72259, 0.72845, 0.73886, 0.75574, 0.77074]
+  t0_data = [2.846, 3.634, 5.140, 8.564, 14.064]
-  ZV_err = [74, 89, 46, 48, 49] .* 1e-5
+  t0_error = [8, 13, 21, 24, 63] .* 1e-3
  for i in 1:5
-    C[i, i] = ZV_err[i] ^ 2 
+    C[i, i] = t0_error[i] ^ 2
  end
-  Zv = ADerrors.cobs(ZV_data, C, "ZV")
-  global zv(beta::Float64) = Zv[b_values .== beta][1]
+  t0_ = Ref{Vector{uwreal}}()
-end
+  global function t0(beta::Float64)
+    if !isassigned(beta)
-let C = zeros(5, 5)
+      t0_.x = ADerrors.cobs(t0_data, C, "t0")
-  #2005.01352
+    end
-  ZS_over_ZP_data = [1.3497, 1.2914, 1.2317, 1.1709, 1.1343]
+    return t0_.x[beta_to_idx[beta]];
-  ZS_over_ZP_err = [83, 64, 48, 23, 25 ] .* 1e-4
+  end
-  C = zeros(5,5)
+  t0_ph_ = Ref{uwreal}()
-  for i in 1:5
+  global function t0_ph()
-    C[i, i] = ZS_over_ZP_err[i] ^ 2
+    if !isassigned(t0_ph)
-  end
+      t0_ph_.x = ADerrors.uwreal([0.4118,0.25e-4],"sqrt(8 t0) (fm)") 
-  ZS_over_ZP = ADerrors.cobs(ZS_over_ZP_data, C, "ZS/ZP")
+    end
-  global zs_over_zp(beta::Float64) = ZS_over_ZP[b_values .== beta][1]
+    return t0_ph_.x
-end
+  end
-let C = zeros(5, 5)
+  global function a(beta::Float64)
-  #2101.10969
+    return t0_ph() / sqrt(8*t0_.x[beta_to_idx[beta]]);
-  RM_err  = [31, 19, 14, 16, 18] .* 1e-3
+  end
-  RM_data = [2.335, 1.869, 1.523, 1.267, 1.149]
+end
-  C=zeros(5,5)
-  for i in 1:5
-    C[i, i] = RM_err[i] ^2 
+let C = zeros(5, 5)
-  end
+  #1808.09236
-  RM = ADerrors.cobs(RM_data, C, "rm")
+  ZA_data = [0.75642, 0.76169, 0.76979, 0.78378, 0.79667]
-  global rm(beta::Float64)         = RM[b_values .== beta][1]
+  ZA_err = [72, 93, 43, 47, 47] .*1e-5
-end
+  for i in 1:5
-let C = zeros(5, 5)
+    C[i, i] = ZA_err[i] ^ 2
-  #1906.03445 LCP1 for heavy quarks
+  end
-  BA_MINUS_BP_data = [-0.324, -0.265, -0.196, -0.119, -0.073]
+  Za = Ref{Vector{uwreal}}()
-  BA_MINUS_BP_err  = [17, 14, 14, 14, 12] .*1e-3 
-  C=zeros(5,5)
+  global function za(beta::Float64)
-  for i in 1:5
+    if !isassigned(Za)
-    C[i,i]  = BA_MINUS_BP_err[i] ^2
+      Za.x = ADerrors.cobs(ZA_data, C, "ZA")
-  end
+    end
+    return za.x[beta_to_idx[beta]];
-  BA_BP = ADerrors.cobs(BA_MINUS_BP_data, C, "ba-bp")
+  end
-  global ba_bp(beta::Float64)      = BA_BP[b_values .== beta][1]
+end
-end
+let C = zeros(5, 5)
-let C = zeros(5, 5)
+  #1808.09236
-  #1906.03445 LCP1 for heavy quarks 
+  ZV_data = [0.72259, 0.72845, 0.73886, 0.75574, 0.77074]
-  ZDATA = [0.8981, 0.9468, 1.0015, 1.0612, 1.0971]
+  ZV_err = [74, 89, 46, 48, 49] .* 1e-5
-  ZERR  = [35, 35, 30, 17, 18 ] .* 1e-4
  for i in 1:5
-    C[i,i]  = ZERR[i] ^2
+    C[i, i] = ZV_err[i] ^ 2 
  end
-  ZZ = ADerrors.cobs(ZDATA, C, "Z")
+  Zv = Ref{Vector{uwreal}}()
-  global Z(beta::Float64)          = ZZ[b_values .== beta][1]
-end
+  global function zv(beta::Float64)
+    if !isassigned(Zv)
+      Zv.x = ADerrors.cobs(ZV_data, C, "ZV")
+    end
+    return zv.x[beta_to_idx[beta]];
+  end
+end
+let C = zeros(5, 5)
+  #2005.01352
+  ZS_over_ZP_data = [1.3497, 1.2914, 1.2317, 1.1709, 1.1343]
+  ZS_over_ZP_err = [83, 64, 48, 23, 25 ] .* 1e-4
+  C = zeros(5,5)
+  for i in 1:5
+    C[i, i] = ZS_over_ZP_err[i] ^ 2
+  end
+  ZS_over_ZP = Ref{Vector{uwreal}}()
+  global function zs_over_zp(beta::Float64)
+    if !isassigned(ZS_over_ZP_data)
+      ZS_over_ZP.x = ADerrors.cobs(ZS_over_ZP_data, C, "ZS/ZP")
+    end
+    return ZS_over_ZP.x[beta_to_idx[beta]];
+  end 
+end
+let C = zeros(5, 5)
+  #2101.10969
+  RM_err  = [31, 19, 14, 16, 18] .* 1e-3
+  RM_data = [2.335, 1.869, 1.523, 1.267, 1.149]
+  C=zeros(5,5)
+  for i in 1:5
+    C[i, i] = RM_err[i] ^2 
+  end
+  RM  = Ref{Vector{uwreal}}()
+  global function rm(beta::Float64)
+    if !isassigned(rm)
+      RM.x = ADerrors.cobs(RM_data, C, "rm")
+    end
+    return rm.x[beta_to_idx[beta]];
+  end  
+end
+let C = zeros(5, 5)
+  #1906.03445 LCP1 for heavy quarks
+  BA_MINUS_BP_data = [-0.324, -0.265, -0.196, -0.119, -0.073]
+  BA_MINUS_BP_err  = [17, 14, 14, 14, 12] .*1e-3 
+  C=zeros(5,5)
+  for i in 1:5
+    C[i,i]  = BA_MINUS_BP_err[i] ^2
+  end
+  BA_BP = Ref{Vector{uwreal}}()
+  global function ba_bp(beta::Float64)
+    if !isassigned(BA_BP)
+      BA_BP.x =  ADerrors.cobs(BA_MINUS_BP_data, C, "ba-bp")
+    end
+    return  BA_BP.x[beta_to_idx[beta]];
+  end
+end
+let C = zeros(5, 5)
+  #1906.03445 LCP1 for heavy quarks 
+  ZDATA = [0.8981, 0.9468, 1.0015, 1.0612, 1.0971]
+  ZERR  = [35, 35, 30, 17, 18 ] .* 1e-4
+  for i in 1:5
+    C[i,i]  = ZERR[i] ^2
+  end
+  ZZ = Ref{Vector{uwreal}}();
+  global function Z(beta::Float64)
+    if !isassinged(ZZ)
+      ZZ.x =ADerrors.cobs(ZDATA, C, "Z")
+    end
+    return ZZ[beta_to_idx[beta]];
+  end
+end
--- a/constants/path_csv.jl
+++ b/constants/path_csv.jl
--- a/docs/build/assets/documenter.js
+++ b/docs/build/assets/documenter.js
 // Generated by Documenter.jl
 requirejs.config({
  paths: {
    'highlight-julia': 'https://cdnjs.cloudflare.com/ajax/libs/highlight.js/11.0.1/languages/julia.min',
    'headroom': 'https://cdnjs.cloudflare.com/ajax/libs/headroom/0.12.0/headroom.min',
    'jqueryui': 'https://cdnjs.cloudflare.com/ajax/libs/jqueryui/1.12.1/jquery-ui.min',
    'katex-auto-render': 'https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.13.11/contrib/auto-render.min',
    'jquery': 'https://cdnjs.cloudflare.com/ajax/libs/jquery/3.6.0/jquery.min',
    'headroom-jquery': 'https://cdnjs.cloudflare.com/ajax/libs/headroom/0.12.0/jQuery.headroom.min',
    'katex': 'https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.13.11/katex.min',
    'highlight': 'https://cdnjs.cloudflare.com/ajax/libs/highlight.js/11.0.1/highlight.min',
    'highlight-julia-repl': 'https://cdnjs.cloudflare.com/ajax/libs/highlight.js/11.0.1/languages/julia-repl.min',
  },
  shim: {
  "highlight-julia": {
    "deps": [
      "highlight"
    ]
  },
  "katex-auto-render": {
    "deps": [
      "katex"
    ]
  },
  "headroom-jquery": {
    "deps": [
      "jquery",
      "headroom"
    ]
  },
  "highlight-julia-repl": {
    "deps": [
      "highlight"
    ]
  }
 }
 });
 ////////////////////////////////////////////////////////////////////////////////
 require(['jquery', 'katex', 'katex-auto-render'], function($, katex, renderMathInElement) {
 $(document).ready(function() {
  renderMathInElement(
    document.body,
    {
  "delimiters": [
    {
      "left": "$",
      "right": "$",
      "display": false
    },
    {
      "left": "$$",
      "right": "$$",
      "display": true
    },
    {
      "left": "\\[",
      "right": "\\]",
      "display": true
    }
  ]
 }
  );
 })
 })
 ////////////////////////////////////////////////////////////////////////////////
 require(['jquery', 'highlight', 'highlight-julia', 'highlight-julia-repl'], function($) {
 $(document).ready(function() {
    hljs.highlightAll();
 })
 })
 ////////////////////////////////////////////////////////////////////////////////
 require([], function() {
 function addCopyButtonCallbacks() {
  for (const el of document.getElementsByTagName("pre")) {
    const button = document.createElement("button");
    button.classList.add("copy-button", "fas", "fa-copy");
    el.appendChild(button);
    const success = function () {
      button.classList.add("success", "fa-check");
      button.classList.remove("fa-copy");
    };
    const failure = function () {
      button.classList.add("error", "fa-times");
      button.classList.remove("fa-copy");
    };
    button.addEventListener("click", function () {
      copyToClipboard(el.innerText).then(success, failure);
      setTimeout(function () {
        button.classList.add("fa-copy");
        button.classList.remove("success", "fa-check", "fa-times");
      }, 5000);
    });
  }
 }
 function copyToClipboard(text) {
  // clipboard API is only available in secure contexts
  if (window.navigator && window.navigator.clipboard) {
    return window.navigator.clipboard.writeText(text);
  } else {
    return new Promise(function (resolve, reject) {
      try {
        const el = document.createElement("textarea");
        el.textContent = text;
        el.style.position = "fixed";
        el.style.opacity = 0;
        document.body.appendChild(el);
        el.select();
        document.execCommand("copy");
        resolve();
      } catch (err) {
        reject(err);
      } finally {
        document.body.removeChild(el);
      }
    });
  }
 }
 if (document.readyState === "loading") {
  document.addEventListener("DOMContentLoaded", addCopyButtonCallbacks);
 } else {
  addCopyButtonCallbacks();
 }
 })
 ////////////////////////////////////////////////////////////////////////////////
 require(['jquery', 'headroom', 'headroom-jquery'], function($, Headroom) {
 // Manages the top navigation bar (hides it when the user starts scrolling down on the
 // mobile).
 window.Headroom = Headroom; // work around buggy module loading?
 $(document).ready(function() {
  $('#documenter .docs-navbar').headroom({
    "tolerance": {"up": 10, "down": 10},
  });
 })
 })
 ////////////////////////////////////////////////////////////////////////////////
 require(['jquery'], function($) {
 // Modal settings dialog
 $(document).ready(function() {
  var settings = $('#documenter-settings');
  $('#documenter-settings-button').click(function(){
    settings.toggleClass('is-active');
  });
  // Close the dialog if X is clicked
  $('#documenter-settings button.delete').click(function(){
    settings.removeClass('is-active');
  });
  // Close dialog if ESC is pressed
  $(document).keyup(function(e) {
    if (e.keyCode == 27) settings.removeClass('is-active');
  });
 });
 })
 ////////////////////////////////////////////////////////////////////////////////
 require(['jquery'], function($) {
 // Manages the showing and hiding of the sidebar.
 $(document).ready(function() {
  var sidebar = $("#documenter > .docs-sidebar");
  var sidebar_button = $("#documenter-sidebar-button")
  sidebar_button.click(function(ev) {
    ev.preventDefault();
    sidebar.toggleClass('visible');
    if (sidebar.hasClass('visible')) {
      // Makes sure that the current menu item is visible in the sidebar.
      $("#documenter .docs-menu a.is-active").focus();
    }
  });
  $("#documenter > .docs-main").bind('click', function(ev) {
    if ($(ev.target).is(sidebar_button)) {
      return;
    }
    if (sidebar.hasClass('visible')) {
      sidebar.removeClass('visible');
    }
  });
 })
 // Resizes the package name / sitename in the sidebar if it is too wide.
 // Inspired by: https://github.com/davatron5000/FitText.js
 $(document).ready(function() {
  e = $("#documenter .docs-autofit");
  function resize() {
    var L = parseInt(e.css('max-width'), 10);
    var L0 = e.width();
    if(L0 > L) {
      var h0 = parseInt(e.css('font-size'), 10);
      e.css('font-size', L * h0 / L0);
      // TODO: make sure it survives resizes?
    }
  }
  // call once and then register events
  resize();
  $(window).resize(resize);
  $(window).on('orientationchange', resize);
 });
 // Scroll the navigation bar to the currently selected menu item
 $(document).ready(function() {
  var sidebar = $("#documenter .docs-menu").get(0);
  var active = $("#documenter .docs-menu .is-active").get(0);
  if(typeof active !== 'undefined') {
    sidebar.scrollTop = active.offsetTop - sidebar.offsetTop - 15;
  }
 })
 })
 ////////////////////////////////////////////////////////////////////////////////
 require(['jquery'], function($) {
 function set_theme(theme) {
  var active = null;
  var disabled = [];
  for (var i = 0; i < document.styleSheets.length; i++) {
    var ss = document.styleSheets[i];
    var themename = ss.ownerNode.getAttribute("data-theme-name");
    if(themename === null) continue; // ignore non-theme stylesheets
    // Find the active theme
    if(themename === theme) active = ss;
    else disabled.push(ss);
  }
  if(active !== null) {
    active.disabled = false;
    if(active.ownerNode.getAttribute("data-theme-primary") === null) {
      document.getElementsByTagName('html')[0].className = "theme--" + theme;
    } else {
      document.getElementsByTagName('html')[0].className = "";
    }
    disabled.forEach(function(ss){
      ss.disabled = true;
    });
  }
  // Store the theme in localStorage
  if(typeof(window.localStorage) !== "undefined") {
    window.localStorage.setItem("documenter-theme", theme);
  } else {
    console.error("Browser does not support window.localStorage");
  }
 }
 // Theme picker setup
 $(document).ready(function() {
  // onchange callback
  $('#documenter-themepicker').change(function themepick_callback(ev){
    var themename = $('#documenter-themepicker option:selected').attr('value');
    set_theme(themename);
  });
  // Make sure that the themepicker displays the correct theme when the theme is retrieved
  // from localStorage
  if(typeof(window.localStorage) !== "undefined") {
    var theme =  window.localStorage.getItem("documenter-theme");
    if(theme !== null) {
      $('#documenter-themepicker option').each(function(i,e) {
        e.selected = (e.value === theme);
      })
    } else {
      $('#documenter-themepicker option').each(function(i,e) {
        e.selected = $("html").hasClass(`theme--${e.value}`);
      })
    }
  }
 })
 })
 ////////////////////////////////////////////////////////////////////////////////
 require(['jquery'], function($) {
 // update the version selector with info from the siteinfo.js and ../versions.js files
 $(document).ready(function() {
  var version_selector = $("#documenter .docs-version-selector");
  var version_selector_select = $("#documenter .docs-version-selector select");
  version_selector_select.change(function(x) {
    target_href = version_selector_select.children("option:selected").get(0).value;
    window.location.href = target_href;
  });
  // add the current version to the selector based on siteinfo.js, but only if the selector is empty
  if (typeof DOCUMENTER_CURRENT_VERSION !== 'undefined' && $('#version-selector > option').length == 0) {
    var option = $("<option value='#' selected='selected'>" + DOCUMENTER_CURRENT_VERSION + "</option>");
    version_selector_select.append(option);
  }
  if (typeof DOC_VERSIONS !== 'undefined') {
    var existing_versions = version_selector_select.children("option");
    var existing_versions_texts = existing_versions.map(function(i,x){return x.text});
    DOC_VERSIONS.forEach(function(each) {
      var version_url = documenterBaseURL + "/../" + each;
      var existing_id = $.inArray(each, existing_versions_texts);
      // if not already in the version selector, add it as a new option,
      // otherwise update the old option with the URL and enable it
      if (existing_id == -1) {
        var option = $("<option value='" + version_url + "'>" + each + "</option>");
        version_selector_select.append(option);
      } else {
        var option = existing_versions[existing_id];
        option.value = version_url;
        option.disabled = false;
      }
    });
  }
  // only show the version selector if the selector has been populated
  if (version_selector_select.children("option").length > 0) {
    version_selector.toggleClass("visible");
  }
 })
 })
--- a/docs/build/assets/search.js
+++ b/docs/build/assets/search.js
 // Generated by Documenter.jl
 requirejs.config({
  paths: {
    'lunr': 'https://cdnjs.cloudflare.com/ajax/libs/lunr.js/2.3.9/lunr.min',
    'lodash': 'https://cdnjs.cloudflare.com/ajax/libs/lodash.js/4.17.21/lodash.min',
    'jquery': 'https://cdnjs.cloudflare.com/ajax/libs/jquery/3.6.0/jquery.min',
  }
 });
 ////////////////////////////////////////////////////////////////////////////////
 require(['jquery', 'lunr', 'lodash'], function($, lunr, _) {
 $(document).ready(function() {
  // parseUri 1.2.2
  // (c) Steven Levithan <stevenlevithan.com>
  // MIT License
  function parseUri (str) {
    var  o   = parseUri.options,
      m   = o.parser[o.strictMode ? "strict" : "loose"].exec(str),
      uri = {},
      i   = 14;
    while (i--) uri[o.key[i]] = m[i] || "";
    uri[o.q.name] = {};
    uri[o.key[12]].replace(o.q.parser, function ($0, $1, $2) {
      if ($1) uri[o.q.name][$1] = $2;
    });
    return uri;
  };
  parseUri.options = {
    strictMode: false,
    key: ["source","protocol","authority","userInfo","user","password","host","port","relative","path","directory","file","query","anchor"],
    q:   {
      name:   "queryKey",
      parser: /(?:^|&)([^&=]*)=?([^&]*)/g
    },
    parser: {
      strict: /^(?:([^:\/?#]+):)?(?:\/\/((?:(([^:@]*)(?::([^:@]*))?)?@)?([^:\/?#]*)(?::(\d*))?))?((((?:[^?#\/]*\/)*)([^?#]*))(?:\?([^#]*))?(?:#(.*))?)/,
      loose:  /^(?:(?![^:@]+:[^:@\/]*@)([^:\/?#.]+):)?(?:\/\/)?((?:(([^:@]*)(?::([^:@]*))?)?@)?([^:\/?#]*)(?::(\d*))?)(((\/(?:[^?#](?![^?#\/]*\.[^?#\/.]+(?:[?#]|$)))*\/?)?([^?#\/]*))(?:\?([^#]*))?(?:#(.*))?)/
    }
  };
  $("#search-form").submit(function(e) {
    e.preventDefault()
  })
  // list below is the lunr 2.1.3 list minus the intersect with names(Base)
  // (all, any, get, in, is, only, which) and (do, else, for, let, where, while, with)
  // ideally we'd just filter the original list but it's not available as a variable
  lunr.stopWordFilter = lunr.generateStopWordFilter([
    'a',
    'able',
    'about',
    'across',
    'after',
    'almost',
    'also',
    'am',
    'among',
    'an',
    'and',
    'are',
    'as',
    'at',
    'be',
    'because',
    'been',
    'but',
    'by',
    'can',
    'cannot',
    'could',
    'dear',
    'did',
    'does',
    'either',
    'ever',
    'every',
    'from',
    'got',
    'had',
    'has',
    'have',
    'he',
    'her',
    'hers',
    'him',
    'his',
    'how',
    'however',
    'i',
    'if',
    'into',
    'it',
    'its',
    'just',
    'least',
    'like',
    'likely',
    'may',
    'me',
    'might',
    'most',
    'must',
    'my',
    'neither',
    'no',
    'nor',
    'not',
    'of',
    'off',
    'often',
    'on',
    'or',
    'other',
    'our',
    'own',
    'rather',
    'said',
    'say',
    'says',
    'she',
    'should',
    'since',
    'so',
    'some',
    'than',
    'that',
    'the',
    'their',
    'them',
    'then',
    'there',
    'these',
    'they',
    'this',
    'tis',
    'to',
    'too',
    'twas',
    'us',
    'wants',
    'was',
    'we',
    'were',
    'what',
    'when',
    'who',
    'whom',
    'why',
    'will',
    'would',
    'yet',
    'you',
    'your'
  ])
  // add . as a separator, because otherwise "title": "Documenter.Anchors.add!"
  // would not find anything if searching for "add!", only for the entire qualification
  lunr.tokenizer.separator = /[\s\-\.]+/
  // custom trimmer that doesn't strip @ and !, which are used in julia macro and function names
  lunr.trimmer = function (token) {
    return token.update(function (s) {
      return s.replace(/^[^a-zA-Z0-9@!]+/, '').replace(/[^a-zA-Z0-9@!]+$/, '')
    })
  }
  lunr.Pipeline.registerFunction(lunr.stopWordFilter, 'juliaStopWordFilter')
  lunr.Pipeline.registerFunction(lunr.trimmer, 'juliaTrimmer')
  var index = lunr(function () {
    this.ref('location')
    this.field('title',{boost: 100})
    this.field('text')
    documenterSearchIndex['docs'].forEach(function(e) {
        this.add(e)
    }, this)
  })
  var store = {}
  documenterSearchIndex['docs'].forEach(function(e) {
      store[e.location] = {title: e.title, category: e.category, page: e.page}
  })
  $(function(){
    searchresults = $('#documenter-search-results');
    searchinfo = $('#documenter-search-info');
    searchbox = $('#documenter-search-query');
    function update_search(querystring) {
      tokens = lunr.tokenizer(querystring)
      results = index.query(function (q) {
        tokens.forEach(function (t) {
          q.term(t.toString(), {
            fields: ["title"],
            boost: 100,
            usePipeline: true,
            editDistance: 0,
            wildcard: lunr.Query.wildcard.NONE
          })
          q.term(t.toString(), {
            fields: ["title"],
            boost: 10,
            usePipeline: true,
            editDistance: 2,
            wildcard: lunr.Query.wildcard.NONE
          })
          q.term(t.toString(), {
            fields: ["text"],
            boost: 1,
            usePipeline: true,
            editDistance: 0,
            wildcard: lunr.Query.wildcard.NONE
          })
        })
      })
      searchinfo.text("Number of results: " + results.length)
      searchresults.empty()
      results.forEach(function(result) {
        data = store[result.ref]
        link = $('<a class="docs-label">'+data.title+'</a>')
        link.attr('href', documenterBaseURL+'/'+result.ref)
        if (data.category != "page"){
          cat = $('<span class="docs-category">('+data.category+', '+data.page+')</span>')
        } else {
          cat = $('<span class="docs-category">('+data.category+')</span>')
        }
        li = $('<li>').append(link).append(" ").append(cat)
        searchresults.append(li)
      })
    }
    function update_search_box() {
      querystring = searchbox.val()
      update_search(querystring)
    }
    searchbox.keyup(_.debounce(update_search_box, 250))
    searchbox.change(update_search_box)
    search_query_uri = parseUri(window.location).queryKey["q"]
    if(search_query_uri !== undefined) {
      search_query = decodeURIComponent(search_query_uri.replace(/\+/g, '%20'))
      searchbox.val(search_query)
    }
    update_search_box();
  })
 })
 })
--- a/docs/build/assets/themes/documenter-dark.css
+++ b/docs/build/assets/themes/documenter-dark.css
--- a/docs/build/assets/themes/documenter-light.css
+++ b/docs/build/assets/themes/documenter-light.css
--- a/docs/build/assets/themeswap.js
+++ b/docs/build/assets/themeswap.js
 // Small function to quickly swap out themes. Gets put into the <head> tag..
 function set_theme_from_local_storage() {
  // Intialize the theme to null, which means default
  var theme = null;
  // If the browser supports the localstorage and is not disabled then try to get the
  // documenter theme
  if(window.localStorage != null) {
    // Get the user-picked theme from localStorage. May be `null`, which means the default
    // theme.
    theme =  window.localStorage.getItem("documenter-theme");
  }
  // Check if the browser supports user color preference
  var darkPreference = false;
  // Check if the users preference is for dark color scheme
  if(window.matchMedia('(prefers-color-scheme: dark)').matches === true) {
    darkPreference = true;
  }
  // Initialize a few variables for the loop:
  //
  //  - active: will contain the index of the theme that should be active. Note that there
  //    is no guarantee that localStorage contains sane values. If `active` stays `null`
  //    we either could not find the theme or it is the default (primary) theme anyway.
  //    Either way, we then need to stick to the primary theme.
  //
  //  - disabled: style sheets that should be disabled (i.e. all the theme style sheets
  //    that are not the currently active theme)
  var active = null; var disabled = []; var darkTheme = null;
  for (var i = 0; i < document.styleSheets.length; i++) {
    var ss = document.styleSheets[i];
    // The <link> tag of each style sheet is expected to have a data-theme-name attribute
    // which must contain the name of the theme. The names in localStorage much match this.
    var themename = ss.ownerNode.getAttribute("data-theme-name");
    // attribute not set => non-theme stylesheet => ignore
    if(themename === null) continue;
    // To distinguish the default (primary) theme, it needs to have the data-theme-primary
    // attribute set.
    var isprimary = (ss.ownerNode.getAttribute("data-theme-primary") !== null);
    // Check if the theme is primary dark theme
    var isDarkTheme = (ss.ownerNode.getAttribute("data-theme-primary-dark") !== null);
    // If ss is for dark theme then set the value of darkTheme to the name of the theme
    if(isDarkTheme) darkTheme = themename;
    // If we find a matching theme (and it's not the default), we'll set active to non-null
    if(themename === theme) active = i;
    // Store the style sheets of inactive themes so that we could disable them
    if(themename !== theme) disabled.push(ss);
  }
  if(active !== null) {
    // If we did find an active theme, we'll (1) add the theme--$(theme) class to <html>
    document.getElementsByTagName('html')[0].className = "theme--" + theme;
    // and (2) disable all the other theme stylesheets
    disabled.forEach(function(ss){
      ss.disabled = true;
    });
  }
  else if(darkTheme !== null && darkPreference === true) {
    // If we did find an active theme, we'll (1) add the theme--$(theme) class to <html>
    document.getElementsByTagName('html')[0].className = "theme--" + darkTheme;
    // and (2) disable all the other theme stylesheets
    disabled.forEach(function(ss){
      if (ss.ownerNode.getAttribute("data-theme-name") !== darkTheme) {
        ss.disabled = true;
      }
    });
  }
 }
 set_theme_from_local_storage();
--- a/docs/build/assets/warner.js
+++ b/docs/build/assets/warner.js
 function maybeAddWarning () {
    // DOCUMENTER_NEWEST is defined in versions.js, DOCUMENTER_CURRENT_VERSION and DOCUMENTER_STABLE
    // in siteinfo.js.
    // If either of these are undefined something went horribly wrong, so we abort.
    if (
            window.DOCUMENTER_NEWEST === undefined ||
            window.DOCUMENTER_CURRENT_VERSION === undefined ||
            window.DOCUMENTER_STABLE === undefined
       ) {
        return
    };
    // Current version is not a version number, so we can't tell if it's the newest version. Abort.
    if (!/v(\d+\.)*\d+/.test(window.DOCUMENTER_CURRENT_VERSION)) {
        return
    };
    // Current version is newest version, so no need to add a warning.
    if (window.DOCUMENTER_NEWEST === window.DOCUMENTER_CURRENT_VERSION) {
        return
    };
    // Add a noindex meta tag (unless one exists) so that search engines don't index this version of the docs.
    if (document.body.querySelector('meta[name="robots"]') === null) {
        const meta = document.createElement('meta');
        meta.name = 'robots';
        meta.content = 'noindex';
        document.getElementsByTagName('head')[0].appendChild(meta);
    };
    const div = document.createElement('div');
    div.classList.add('outdated-warning-overlay');
    const closer = document.createElement('button');
    closer.classList.add('outdated-warning-closer', 'delete');
    closer.addEventListener('click', function () {
        document.body.removeChild(div);
    });
    const href = window.documenterBaseURL + '/../' + window.DOCUMENTER_STABLE;
    div.innerHTML = 'This documentation is not for the latest version. <br> <a href="' + href + '">Go to the latest documentation</a>.';
    div.appendChild(closer);
    document.body.appendChild(div);
 };
 if (document.readyState === 'loading') {
    document.addEventListener('DOMContentLoaded', maybeAddWarning);
 } else {
    maybeAddWarning();
 };
--- a/docs/build/index.html
+++ b/docs/build/index.html
 <!DOCTYPE html>
 <html lang="en"><head><meta charset="UTF-8"/><meta name="viewport" content="width=device-width, initial-scale=1.0"/><title>Home · juobs Documentation</title><script data-outdated-warner src="assets/warner.js"></script><link href="https://cdnjs.cloudflare.com/ajax/libs/lato-font/3.0.0/css/lato-font.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/juliamono/0.039/juliamono-regular.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.3/css/fontawesome.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.3/css/solid.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.3/css/brands.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.13.11/katex.min.css" rel="stylesheet" type="text/css"/><script>documenterBaseURL="."</script><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.6/require.min.js" data-main="assets/documenter.js"></script><script src="siteinfo.js"></script><script src="../versions.js"></script><link class="docs-theme-link" rel="stylesheet" type="text/css" href="assets/themes/documenter-dark.css" data-theme-name="documenter-dark" data-theme-primary-dark/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="assets/themes/documenter-light.css" data-theme-name="documenter-light" data-theme-primary/><script src="assets/themeswap.js"></script></head><body><div id="documenter"><nav class="docs-sidebar"><div class="docs-package-name"><span class="docs-autofit"><a href="index.html">juobs Documentation</a></span></div><form class="docs-search" action="search.html"><input class="docs-search-query" id="documenter-search-query" name="q" type="text" placeholder="Search docs"/></form><ul class="docs-menu"><li class="is-active"><a class="tocitem" href="index.html">Home</a><ul class="internal"><li><a class="tocitem" href="#Contents"><span>Contents</span></a></li></ul></li><li><a class="tocitem" href="reader.html">Reader</a></li><li><a class="tocitem" href="tools.html">Tools</a></li><li><a class="tocitem" href="obs.html">Observables</a></li><li><a class="tocitem" href="linalg.html">Linear Algebra</a></li></ul><div class="docs-version-selector field has-addons"><div class="control"><span class="docs-label button is-static is-size-7">Version</span></div><div class="docs-selector control is-expanded"><div class="select is-fullwidth is-size-7"><select id="documenter-version-selector"></select></div></div></div></nav><div class="docs-main"><header class="docs-navbar"><nav class="breadcrumb"><ul class="is-hidden-mobile"><li class="is-active"><a href="index.html">Home</a></li></ul><ul class="is-hidden-tablet"><li class="is-active"><a href="index.html">Home</a></li></ul></nav><div class="docs-right"><a class="docs-edit-link" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs" title="Edit on GitLab"><span class="docs-icon fab"></span><span class="docs-label is-hidden-touch">Edit on GitLab</span></a><a class="docs-settings-button fas fa-cog" id="documenter-settings-button" href="#" title="Settings"></a><a class="docs-sidebar-button fa fa-bars is-hidden-desktop" id="documenter-sidebar-button" href="#"></a></div></header><article class="content" id="documenter-page"><h1 id="DOCUMENTATION"><a class="docs-heading-anchor" href="#DOCUMENTATION">DOCUMENTATION</a><a id="DOCUMENTATION-1"></a><a class="docs-heading-anchor-permalink" href="#DOCUMENTATION" title="Permalink"></a></h1><h2 id="Contents"><a class="docs-heading-anchor" href="#Contents">Contents</a><a id="Contents-1"></a><a class="docs-heading-anchor-permalink" href="#Contents" title="Permalink"></a></h2><ul><li><a href="reader.html#Reader">Reader</a></li><li><a href="tools.html#Tools">Tools</a></li><li><a href="obs.html#Observables">Observables</a></li><li><a href="linalg.html#Linear-Algebra">Linear Algebra</a></li></ul></article><nav class="docs-footer"><a class="docs-footer-nextpage" href="reader.html">Reader »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 0.27.15 on <span class="colophon-date" title="Friday 1 July 2022 12:33">Friday 1 July 2022</span>. Using Julia version 1.6.5.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>
--- a/docs/build/linalg.html
+++ b/docs/build/linalg.html
 <!DOCTYPE html>
 <html lang="en"><head><meta charset="UTF-8"/><meta name="viewport" content="width=device-width, initial-scale=1.0"/><title>Linear Algebra · juobs Documentation</title><script data-outdated-warner src="assets/warner.js"></script><link href="https://cdnjs.cloudflare.com/ajax/libs/lato-font/3.0.0/css/lato-font.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/juliamono/0.039/juliamono-regular.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.3/css/fontawesome.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.3/css/solid.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.3/css/brands.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.13.11/katex.min.css" rel="stylesheet" type="text/css"/><script>documenterBaseURL="."</script><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.6/require.min.js" data-main="assets/documenter.js"></script><script src="siteinfo.js"></script><script src="../versions.js"></script><link class="docs-theme-link" rel="stylesheet" type="text/css" href="assets/themes/documenter-dark.css" data-theme-name="documenter-dark" data-theme-primary-dark/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="assets/themes/documenter-light.css" data-theme-name="documenter-light" data-theme-primary/><script src="assets/themeswap.js"></script></head><body><div id="documenter"><nav class="docs-sidebar"><div class="docs-package-name"><span class="docs-autofit"><a href="index.html">juobs Documentation</a></span></div><form class="docs-search" action="search.html"><input class="docs-search-query" id="documenter-search-query" name="q" type="text" placeholder="Search docs"/></form><ul class="docs-menu"><li><a class="tocitem" href="index.html">Home</a></li><li><a class="tocitem" href="reader.html">Reader</a></li><li><a class="tocitem" href="tools.html">Tools</a></li><li><a class="tocitem" href="obs.html">Observables</a></li><li class="is-active"><a class="tocitem" href="linalg.html">Linear Algebra</a></li></ul><div class="docs-version-selector field has-addons"><div class="control"><span class="docs-label button is-static is-size-7">Version</span></div><div class="docs-selector control is-expanded"><div class="select is-fullwidth is-size-7"><select id="documenter-version-selector"></select></div></div></div></nav><div class="docs-main"><header class="docs-navbar"><nav class="breadcrumb"><ul class="is-hidden-mobile"><li class="is-active"><a href="linalg.html">Linear Algebra</a></li></ul><ul class="is-hidden-tablet"><li class="is-active"><a href="linalg.html">Linear Algebra</a></li></ul></nav><div class="docs-right"><a class="docs-edit-link" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs" title="Edit on GitLab"><span class="docs-icon fab"></span><span class="docs-label is-hidden-touch">Edit on GitLab</span></a><a class="docs-settings-button fas fa-cog" id="documenter-settings-button" href="#" title="Settings"></a><a class="docs-sidebar-button fa fa-bars is-hidden-desktop" id="documenter-sidebar-button" href="#"></a></div></header><article class="content" id="documenter-page"><h1 id="Linear-Algebra"><a class="docs-heading-anchor" href="#Linear-Algebra">Linear Algebra</a><a id="Linear-Algebra-1"></a><a class="docs-heading-anchor-permalink" href="#Linear-Algebra" title="Permalink"></a></h1><article class="docstring"><header><a class="docstring-binding" id="juobs.uweigvals" href="#juobs.uweigvals"><code>juobs.uweigvals</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">uweigvals(a::Matrix{uwreal}; iter = 30)
 uweigvals(a::Matrix{uwreal}, b::Matrix{uwreal};  iter = 30)</code></pre><p>This function computes the eigenvalues of a matrix of uwreal objects. If a second matrix b is given as input, it returns the generalised eigenvalues instead. It takes as input:</p><ul><li><p><code>a::Matrix{uwreal}</code> : a matrix of uwreal</p></li><li><p><code>b::Matrix{uwreal}</code> : a matrix of uwreal, optional</p></li><li><p><code>iter=30</code>: optional flag to set the iterations of the qr algorithm used to solve the eigenvalue problem</p></li></ul><p>It returns:</p><ul><li><code>res = Vector{uwreal}</code>: a vector where each elements is an eigenvalue </li></ul><pre><code class="language- hljs">a = Matrix{uwreal}(nothing, n,n) ## n*n matrix of uwreal with nothing entries
 b = Matrix{uwreal}(nothing, n,n) ## n*n matrix of uwreal with nothing entries
 res = uweigvals(a)  ##eigenvalues
 res1 = uweigvals(a,b) ## generalised eigenvalues</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article><article class="docstring"><header><a class="docstring-binding" id="juobs.uweigvecs" href="#juobs.uweigvecs"><code>juobs.uweigvecs</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">uweigvecs(a::Matrix{uwreal}; iter = 30)
 uweigvecs(a::Matrix{uwreal}, b::Matrix{uwreal};  iter = 30)</code></pre><p>This function computes the eigenvectors of a matrix of uwreal objects. If a second matrix b is given as input, it returns the generalised eigenvectors instead. It takes as input:</p><ul><li><p><code>a::Matrix{uwreal}</code> : a matrix of uwreal</p></li><li><p><code>b::Matrix{uwreal}</code> : a matrix of uwreal, optional</p></li><li><p><code>iter=30</code> : the number of iterations of the qr algorithm used to extract the eigenvalues </p></li></ul><p>It returns:</p><ul><li><code>res = Matrix{uwreal}</code>: a matrix where each column is an eigenvector </li></ul><p>Examples:</p><pre><code class="language- hljs">a = Matrix{uwreal}(nothing, n,n) ## n*n matrix of uwreal with nothing entries
 b = Matrix{uwreal}(nothing, n,n) ## n*n matrix of uwreal with nothing entries
 res = uweigvecs(a)  ##eigenvectors in column 
 res1 = uweigvecs(a,b) ## generalised eigenvectors in column </code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article><article class="docstring"><header><a class="docstring-binding" id="juobs.uweigen" href="#juobs.uweigen"><code>juobs.uweigen</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">uweigen(a::Matrix{uwreal}; iter = 30)
 uweigen(a::Matrix{uwreal}, b::Matrix{uwreal};  iter = 30)</code></pre><p>This function computes the eigenvalues and the eigenvectors  of a matrix of uwreal objects. If a second matrix b is given as input, it returns the generalised eigenvalues and eigenvectors instead. It takes as input:</p><ul><li><p><code>a::Matrix{uwreal}</code> : a matrix of uwreal</p></li><li><p><code>b::Matrix{uwreal}</code> : a matrix of uwreal, optional</p></li><li><p><code>iter=30</code> : the number of iterations of the qr algorithm used to extract the eigenvalues </p></li></ul><p>It returns:</p><ul><li><p><code>evals = Vector{uwreal}</code>: a vector where each elements is an eigenvalue </p></li><li><p><code>evecs  = Matrix{uwreal}</code>: a matrix where the i-th column is the eigenvector of the i-th eigenvalue</p></li></ul><p>Examples:</p><pre><code class="language- hljs">a = Matrix{uwreal}(nothing, n,n) ## n*n matrix of uwreal with nothing entries
 b = Matrix{uwreal}(nothing, n,n) ## n*n matrix of uwreal with nothing entries
 eval, evec = uweigen(a)   
 eval1, evec1 = uweigvecs(a,b) </code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article><article class="docstring"><header><a class="docstring-binding" id="juobs.get_matrix" href="#juobs.get_matrix"><code>juobs.get_matrix</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">get_matrix(diag::Vector{Array}, upper::Vector{Array} )</code></pre><p>This function allows the user to build an array of matrices, where each matrix is a  symmetric matrix of correlators at a given timeslice.   </p><p>It takes as input:</p><ul><li><p><code>diag::Vector{Array}</code>:  vector of correlators. Each correlator will constitute a diagonal element of the matrices A[i] where i runs over the timeslices, i.e.                <code>A[i][1,1] = diag[1], .... A[i][n,n] = diag[n]</code>                <code>Given n=length(diag)</code>, the matrices will have dimension n*n </p></li><li><p><code>upper::Vector{Array}</code>:  vector of correlators liying on the upper diagonal position.               <code>A[i][1,2] = upper[1], .. , A[i][1,n] = upper[n-1],</code>               <code>A[i][2,3] = upper[n], .. , A[i][n-1,n] = upper[n(n-1)/2]</code>                Given n, <code>length(upper)=n(n-1)/2</code></p></li></ul><p>The method returns an array of symmetric matrices of dimension n for each timeslice  </p><p>Examples:</p><pre><code class="language- hljs">## load data 
 pp_data = read_mesons(path, &quot;G5&quot;, &quot;G5&quot;)
 pa_data = read_mesons(path, &quot;G5&quot;, &quot;G0G5&quot;)
 aa_data = read_mesons(path, &quot;G0G5&quot;, &quot;G0G5&quot;)
 ## create Corr struct
 corr_pp = corr_obs.(pp_data)
 corr_pa = corr_obs.(pa_data)
 corr_aa = corr_obs.(aa_data) # array of correlators for different \mu_q combinations
 ## set up matrices
 corr_diag = [corr_pp[1], corr_aa[1]] 
 corr_upper = [corr_pa[1]]
 matrices = [corr_diag, corr_upper]
 Julia&gt; matrices[i]
       pp[i]  ap[i]
       pa[i]  aa[i]
 ## where i runs over the timeslices</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article><article class="docstring"><header><a class="docstring-binding" id="juobs.energies" href="#juobs.energies"><code>juobs.energies</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">energies(evals::Vector{Array}; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)</code></pre><p>This method computes the energy level from the eigenvalues according to:</p><p><span>$E_i(t) = \log(λ(t) / λ(t+1))$</span></p><p>where <code>i=1,..,n</code> with <code>n=length(evals[1])</code> and <code>t=1,..,T</code>  total time slices. It returns a vector array en where each entry en[i][t] contains the i-th states energy at time t </p><p>Examples:</p><pre><code class="language- hljs">## load data
 pp_data = read_mesons(path, &quot;G5&quot;, &quot;G5&quot;)
 pa_data = read_mesons(path, &quot;G5&quot;, &quot;G0G5&quot;)
 aa_data = read_mesons(path, &quot;G0G5&quot;, &quot;G0G5&quot;)
 ## create Corr struct
 corr_pp = corr_obs.(pp_data)
 corr_pa = corr_obs.(pa_data)
 corr_aa = corr_obs.(aa_data) # array of correlators for different \mu_q combinations
 ## set up matrices 
 corr_diag = [corr_pp[1], corr_aa[1]] 
 corr_upper = [corr_pa[1]]
 matrices = [corr_diag, corr_upper]
 ## solve the GEVP
 evals = getall_eigvals(matrices, 5) #where t_0=5
 en = energies(evals)
 Julia&gt; en[i] # i-th state energy at each timeslice</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article><article class="docstring"><header><a class="docstring-binding" id="juobs.getall_eigvals" href="#juobs.getall_eigvals"><code>juobs.getall_eigvals</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">getall_eigvals(a::Vector{Matrix}, t0; iter=30 )</code></pre><p>This function solves a GEVP problem, returning the eigenvalues, for a list of matrices, taking as generalised matrix the one  at index t0, i.e:</p><p><span>$C(t_i)v_i = λ_i C(t_0) v_i$</span>,   with   i=1:lenght(a)</p><p>It takes as input:</p><ul><li><p><code>a::Vector{Matrix}</code> : a vector of matrices</p></li><li><p><code>t0::Int64</code> : idex value at which the fixed matrix is taken</p></li><li><p><code>iter=30</code> : the number of iterations of the qr algorithm used to extract the eigenvalues </p></li></ul><p>It returns:</p><ul><li><code>res</code> = Vector{Vector{uwreal}}</li></ul><p>where <code>res[i]</code> are the generalised eigenvalues of the i-th matrix of the input array. </p><p>Examples:</p><pre><code class="language- hljs">## load data
 pp_data = read_mesons(path, &quot;G5&quot;, &quot;G5&quot;)
 pa_data = read_mesons(path, &quot;G5&quot;, &quot;G0G5&quot;)
 aa_data = read_mesons(path, &quot;G0G5&quot;, &quot;G0G5&quot;)
 ## create Corr struct
 corr_pp = corr_obs.(pp_data)
 corr_pa = corr_obs.(pa_data)
 corr_aa = corr_obs.(aa_data) # array of correlators for different \mu_q combinations
 ## set up matrices 
 corr_diag = [corr_pp[1], corr_aa[1]] 
 corr_upper = [corr_pa[1]]
 matrices = [corr_diag, corr_upper]
 ## solve the GEVP
 #evals = getall_eigvals(matrices, 5) #where t_0=5
 Julia&gt;</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article><article class="docstring"><header><a class="docstring-binding" id="juobs.getall_eigvecs" href="#juobs.getall_eigvecs"><code>juobs.getall_eigvecs</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">getall_eigvecs(a::Vector{Matrix}, delta_t; iter=30 )</code></pre><p>This function solves a GEVP problem, returning the eigenvectors, for a list of matrices.</p><p><span>$C(t_i)v_i = λ_i C(t_i-\delta_t) v_i$</span>, with i=1:lenght(a)</p><p>Here <code>delta_t</code> is the time shift within the two matrices of the problem, and is kept fixed. It takes as input:</p><ul><li><p><code>a::Vector{Matrix}</code> : a vector of matrices</p></li><li><p><code>delta_t::Int64</code> : the fixed time shift t-t_0</p></li><li><p><code>iter=30</code> : the number of iterations of the qr algorithm used to extract the eigenvalues </p></li></ul><p>It returns:</p><ul><li><code>res = Vector{Matrix{uwreal}}</code></li></ul><p>where each <code>res[i]</code> is a matrix  with the eigenvectors as  columns Examples:</p><pre><code class="language- hljs">mat_array = get_matrix(diag, upper_diag)
 evecs = getall_eigvecs(mat_array, 5)</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="obs.html">« Observables</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 0.27.15 on <span class="colophon-date" title="Friday 1 July 2022 12:33">Friday 1 July 2022</span>. Using Julia version 1.6.5.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>
--- a/docs/build/obs.html
+++ b/docs/build/obs.html
 <!DOCTYPE html>
 <html lang="en"><head><meta charset="UTF-8"/><meta name="viewport" content="width=device-width, initial-scale=1.0"/><title>Observables · juobs Documentation</title><script data-outdated-warner src="assets/warner.js"></script><link href="https://cdnjs.cloudflare.com/ajax/libs/lato-font/3.0.0/css/lato-font.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/juliamono/0.039/juliamono-regular.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.3/css/fontawesome.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.3/css/solid.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.3/css/brands.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.13.11/katex.min.css" rel="stylesheet" type="text/css"/><script>documenterBaseURL="."</script><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.6/require.min.js" data-main="assets/documenter.js"></script><script src="siteinfo.js"></script><script src="../versions.js"></script><link class="docs-theme-link" rel="stylesheet" type="text/css" href="assets/themes/documenter-dark.css" data-theme-name="documenter-dark" data-theme-primary-dark/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="assets/themes/documenter-light.css" data-theme-name="documenter-light" data-theme-primary/><script src="assets/themeswap.js"></script></head><body><div id="documenter"><nav class="docs-sidebar"><div class="docs-package-name"><span class="docs-autofit"><a href="index.html">juobs Documentation</a></span></div><form class="docs-search" action="search.html"><input class="docs-search-query" id="documenter-search-query" name="q" type="text" placeholder="Search docs"/></form><ul class="docs-menu"><li><a class="tocitem" href="index.html">Home</a></li><li><a class="tocitem" href="reader.html">Reader</a></li><li><a class="tocitem" href="tools.html">Tools</a></li><li class="is-active"><a class="tocitem" href="obs.html">Observables</a></li><li><a class="tocitem" href="linalg.html">Linear Algebra</a></li></ul><div class="docs-version-selector field has-addons"><div class="control"><span class="docs-label button is-static is-size-7">Version</span></div><div class="docs-selector control is-expanded"><div class="select is-fullwidth is-size-7"><select id="documenter-version-selector"></select></div></div></div></nav><div class="docs-main"><header class="docs-navbar"><nav class="breadcrumb"><ul class="is-hidden-mobile"><li class="is-active"><a href="obs.html">Observables</a></li></ul><ul class="is-hidden-tablet"><li class="is-active"><a href="obs.html">Observables</a></li></ul></nav><div class="docs-right"><a class="docs-edit-link" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs" title="Edit on GitLab"><span class="docs-icon fab"></span><span class="docs-label is-hidden-touch">Edit on GitLab</span></a><a class="docs-settings-button fas fa-cog" id="documenter-settings-button" href="#" title="Settings"></a><a class="docs-sidebar-button fa fa-bars is-hidden-desktop" id="documenter-sidebar-button" href="#"></a></div></header><article class="content" id="documenter-page"><h1 id="Observables"><a class="docs-heading-anchor" href="#Observables">Observables</a><a id="Observables-1"></a><a class="docs-heading-anchor-permalink" href="#Observables" title="Permalink"></a></h1><article class="docstring"><header><a class="docstring-binding" id="juobs.meff" href="#juobs.meff"><code>juobs.meff</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)     
 meff(corr::Corr, plat::Vector{Int64}; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)</code></pre><p>Computes effective mass for a given correlator corr at a given plateau <code>plat</code>. Correlator can be passed as an <code>Corr</code> struct or <code>Vector{uwreal}</code>.</p><p>The flags <code>pl</code> and <code>data</code> allow to show the plots and return data as an extra result.</p><pre><code class="language- hljs">data = read_mesons(path, &quot;G5&quot;, &quot;G5&quot;)
 corr_pp = corr_obs.(data)
 m = meff(corr_pp[1], [50, 60], pl=false)</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article><article class="docstring"><header><a class="docstring-binding" id="juobs.mpcac" href="#juobs.mpcac"><code>juobs.mpcac</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
 mpcac(a0p::Corr, pp::Corr, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)</code></pre><p>Computes the bare PCAC mass for a given correlator <code>a0p</code> and <code>pp</code> at a given plateau <code>plat</code>. Correlator can be passed as an <code>Corr</code> struct or <code>Vector{uwreal}</code>.</p><p>The flags <code>pl</code> and <code>data</code> allow to show the plots and return data as an extra result. The <code>ca</code> variable allows to compute <code>mpcac</code> using the improved axial current.</p><pre><code class="language- hljs">data_pp = read_mesons(path, &quot;G5&quot;, &quot;G5&quot;)
 data_a0p = read_mesons(path, &quot;G5&quot;, &quot;G0G5&quot;)
 corr_pp = corr_obs.(data_pp)
 corr_a0p = corr_obs.(data_a0p)
 m12 = mpcac(corr_a0p, corr_pp, [50, 60], pl=false)
 p0 = 9.2056
 p1 = -13.9847
 g2 = 1.73410
 ca = -0.006033 * g2 *( 1 + exp(p0 + p1/g2))
 m12 = mpcac(corr_a0p, corr_pp, [50, 60], pl=false, ca=ca)</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article><article class="docstring"><header><a class="docstring-binding" id="juobs.dec_const" href="#juobs.dec_const"><code>juobs.dec_const</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">dec_const(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
 dec_const(a0p::Corr, pp::Corr, plat::Vector{Int64}, m::uwreal; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
 dec_const(a0pL::Vector{uwreal}, a0pR::Vector{uwreal}, ppL::Vector{uwreal}, ppR::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
 dec_const(a0pL::Corr, a0pR::Corr, ppL::Corr, ppR::Corr, plat::Vector{Int64}, m::uwreal; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)</code></pre><p>Computes the bare decay constant using <span>$A_0P$</span> and <span>$PP$</span> correlators . The decay constant is computed in the plateau <code>plat</code>. Correlator can be passed as an <code>Corr</code> struct or <code>Vector{uwreal}</code>. If it is passed as a uwreal vector, effective mass <code>m</code> and source position <code>y0</code> must be specified.</p><p>The flags <code>pl</code> and <code>data</code> allow to show the plots and return data as an extra result. The <code>ca</code> variable allows to compute <code>dec_const</code> using the improved axial current.</p><p><strong>The method assumes that the source is close to the boundary.</strong> It takes the following ratio to cancel boundary effects. <span>$R = \frac{f_A(x_0, y_0)}{\sqrt{f_P(T-y_0, y_0)}} * e^{m (x_0 - T/2)}$</span></p><p><strong>If left and right correlators are included in the input. The result is computed with the following ratio</strong> <span>$R = \sqrt{f_A(x_0, y_0) * f_A(x_0, T - 1 - y_0) / f_P(T - 1 - y_0, y_0)}$</span></p><pre><code class="language- hljs">data_pp = read_mesons(path, &quot;G5&quot;, &quot;G5&quot;, legacy=true)
 data_a0p = read_mesons(path, &quot;G5&quot;, &quot;G0G5&quot;, legacy=true)
 corr_pp = corr_obs.(data_pp, L=32)
 corr_a0p = corr_obs.(data_a0p, L=32)
 corr_a0pL, corr_a0pR = [corr_a0p[1], corr_a0p[2]]
 corr_ppL, corr_ppR = [corr_pp[1], corr_pp[2]]
 m = meff(corr_pp[1], [50, 60], pl=false)
 beta = 3.46
 p0 = 9.2056
 p1 = -13.9847
 g2 = 6 / beta
 ca = -0.006033 * g2 *( 1 + exp(p0 + p1/g2))
 f = dec_const(corr_a0p[1], corr_pp[1], [50, 60], m, pl=true, ca=ca)
 f_ratio = dec_const(corr_a0pL, corr_a0pR, corr_ppL, corr_ppR, [50, 60], m, pl=true, ca=ca)</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article><article class="docstring"><header><a class="docstring-binding" id="juobs.dec_const_pcvc" href="#juobs.dec_const_pcvc"><code>juobs.dec_const_pcvc</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
 dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
 dec_const_pcvc(ppL::Vector{uwreal}, ppR::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
 dec_const_pcvc(corrL::Corr, corrR::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)</code></pre><p>Computes decay constant using the PCVC relation for twisted mass fermions. The decay constant is computed in the plateau <code>plat</code>. Correlator can be passed as an <code>Corr</code> struct or <code>Vector{uwreal}</code>. If it is passed as a uwreal vector, vector of twisted masses <code>mu</code> and source position <code>y0</code> must be specified.</p><p>The flags <code>pl</code> and <code>data</code> allow to show the plots and return data as an extra result.</p><p><strong>The method extract the matrix element assuming that the source is in the bulk. ** **If left and right correlators are included in the input. The result is computed with a ratio that cancels boundary effects:</strong> <span>$R = \sqrt{f_P(x_0, y_0) * f_P(x_0, T - 1 - y_0) / f_P(T - 1 - y_0, y_0)}$</span></p><pre><code class="language- hljs">data = read_mesons(path, &quot;G5&quot;, &quot;G5&quot;)
 corr_pp = corr_obs.(data, L=32)
 m = meff(corr_pp[1], [50, 60], pl=false)
 f = dec_const_pcvc(corr_pp[1], [50, 60], m, pl=false)
 #left and right correlators
 f_ratio = dec_const_pcvc(ppL, ppR, [50, 60], m)</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article><article class="docstring"><header><a class="docstring-binding" id="juobs.comp_t0" href="#juobs.comp_t0"><code>juobs.comp_t0</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Matrix{Float64}, Nothing}=nothing, npol::Int64=2, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, info::Bool=false)
 comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Vector{Matrix{Float64}}, Nothing}=nothing, npol::Int64=2, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, info::Bool=false)</code></pre><p>Computes <code>t0</code> using the energy density of the action <code>Ysl</code>(Yang-Mills action). <code>t0</code> is computed in the plateau <code>plat</code>. A polynomial interpolation in <code>t</code> is performed to find <code>t0</code>, where <code>npol</code> is the degree of the polynomial (linear fit by default)</p><p>The flag <code>pl</code> allows to show the plot. </p><p>The flag <code>info</code> provides extra output that contains information about the primary observables. The function returns the primary observables <span>$&lt;WY&gt;$</span> and <span>$&lt;W&gt;$</span> (it returns the observable &lt;Y&gt; if rw=nothing)</p><pre><code class="language- hljs">#Single replica
 Y = read_ms(path)
 rw = read_ms(path_rw)
 t0, Yobs = comp_t0(Y, [38, 58], L=32, info=true)
 t0_r, WYobs, Wobs = comp_t0(Y, [38, 58], L=32, rw=rw, info=true)
 #Two replicas
 Y1 = read_ms(path1)
 Y2 = read_ms(path2)
 rw1 = read_ms(path_rw1)
 rw2 = read_ms(path_rw2)
 t0 = comp_t0([Y1, Y2], [38, 58], L=32, pl=true)
 t0_r = comp_t0(Y, [38, 58], L=32, rw=[rw1, rw2], pl=true)
 </code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="tools.html">« Tools</a><a class="docs-footer-nextpage" href="linalg.html">Linear Algebra »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 0.27.15 on <span class="colophon-date" title="Friday 1 July 2022 12:33">Friday 1 July 2022</span>. Using Julia version 1.6.5.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>
--- a/docs/build/reader.html
+++ b/docs/build/reader.html
 <!DOCTYPE html>
 <html lang="en"><head><meta charset="UTF-8"/><meta name="viewport" content="width=device-width, initial-scale=1.0"/><title>Reader · juobs Documentation</title><script data-outdated-warner src="assets/warner.js"></script><link href="https://cdnjs.cloudflare.com/ajax/libs/lato-font/3.0.0/css/lato-font.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/juliamono/0.039/juliamono-regular.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.3/css/fontawesome.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.3/css/solid.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.3/css/brands.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.13.11/katex.min.css" rel="stylesheet" type="text/css"/><script>documenterBaseURL="."</script><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.6/require.min.js" data-main="assets/documenter.js"></script><script src="siteinfo.js"></script><script src="../versions.js"></script><link class="docs-theme-link" rel="stylesheet" type="text/css" href="assets/themes/documenter-dark.css" data-theme-name="documenter-dark" data-theme-primary-dark/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="assets/themes/documenter-light.css" data-theme-name="documenter-light" data-theme-primary/><script src="assets/themeswap.js"></script></head><body><div id="documenter"><nav class="docs-sidebar"><div class="docs-package-name"><span class="docs-autofit"><a href="index.html">juobs Documentation</a></span></div><form class="docs-search" action="search.html"><input class="docs-search-query" id="documenter-search-query" name="q" type="text" placeholder="Search docs"/></form><ul class="docs-menu"><li><a class="tocitem" href="index.html">Home</a></li><li class="is-active"><a class="tocitem" href="reader.html">Reader</a></li><li><a class="tocitem" href="tools.html">Tools</a></li><li><a class="tocitem" href="obs.html">Observables</a></li><li><a class="tocitem" href="linalg.html">Linear Algebra</a></li></ul><div class="docs-version-selector field has-addons"><div class="control"><span class="docs-label button is-static is-size-7">Version</span></div><div class="docs-selector control is-expanded"><div class="select is-fullwidth is-size-7"><select id="documenter-version-selector"></select></div></div></div></nav><div class="docs-main"><header class="docs-navbar"><nav class="breadcrumb"><ul class="is-hidden-mobile"><li class="is-active"><a href="reader.html">Reader</a></li></ul><ul class="is-hidden-tablet"><li class="is-active"><a href="reader.html">Reader</a></li></ul></nav><div class="docs-right"><a class="docs-edit-link" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs" title="Edit on GitLab"><span class="docs-icon fab"></span><span class="docs-label is-hidden-touch">Edit on GitLab</span></a><a class="docs-settings-button fas fa-cog" id="documenter-settings-button" href="#" title="Settings"></a><a class="docs-sidebar-button fa fa-bars is-hidden-desktop" id="documenter-sidebar-button" href="#"></a></div></header><article class="content" id="documenter-page"><h1 id="Reader"><a class="docs-heading-anchor" href="#Reader">Reader</a><a id="Reader-1"></a><a class="docs-heading-anchor-permalink" href="#Reader" title="Permalink"></a></h1><article class="docstring"><header><a class="docstring-binding" id="juobs.read_mesons" href="#juobs.read_mesons"><code>juobs.read_mesons</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{String, Nothing}=nothing, legacy::Bool=false)
 read_mesons(path::Vector{String}, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{String, Nothing}=nothing, legacy::Bool=false)</code></pre><p>This function read a mesons dat file at a given path and returns a vector of <code>CData</code> structures for different masses and Dirac structures. Dirac structures <code>g1</code> and/or <code>g2</code> can be passed as string arguments in order to filter correaltors. ADerrors id can be specified as argument. If is not specified, the <code>id</code> is fixed according to the ensemble name (example: &quot;H400&quot;-&gt; id = &quot;H400&quot;)</p><p>*For the old version (without smearing, distance preconditioning and theta) set legacy=true.</p><p>Examples:</p><pre><code class="language- hljs">read_mesons(path)
 read_mesons(path, &quot;G5&quot;)
 read_mesons(path, nothing, &quot;G5&quot;)
 read_mesons(path, &quot;G5&quot;, &quot;G5&quot;)
 read_mesons(path, &quot;G5&quot;, &quot;G5&quot;, id=&quot;H100&quot;)
 read_mesons(path, &quot;G5_d2&quot;, &quot;G5_d2&quot;, legacy=true)</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article><article class="docstring"><header><a class="docstring-binding" id="juobs.read_ms" href="#juobs.read_ms"><code>juobs.read_ms</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">read_ms(path::String; id::Union{String, Nothing}=nothing, dtr::Int64=1, obs::String=&quot;Y&quot;)</code></pre><p>Reads openQCD ms dat files at a given path. This method return YData: </p><ul><li><p><code>t(t)</code>: flow time values</p></li><li><p><code>obs(icfg, x0, t)</code>: the time-slice sums of the densities of the observable (Wsl, Ysl or Qsl)</p></li><li><p><code>vtr</code>: vector that contains trajectory number</p></li><li><p><code>id</code>: ensmble id</p></li></ul><p><code>dtr</code> = <code>dtr_cnfg</code> / <code>dtr_ms</code>, where <code>dtr_cnfg</code> is the number of trajectories computed before saving the configuration. <code>dtr_ms</code> is the same but applied to the ms.dat file.</p><p>Examples:</p><pre><code class="language- hljs">Y = read_ms(path)</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article><article class="docstring"><header><a class="docstring-binding" id="juobs.read_ms1" href="#juobs.read_ms1"><code>juobs.read_ms1</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">read_ms1(path::String; v::String=&quot;1.2&quot;)</code></pre><p>Reads openQCD ms1 dat files at a given path. This method returns a matrix <code>W[irw, icfg]</code> that contains the reweighting factors, where <code>irw</code> is the <code>rwf</code> index and icfg the configuration number. The function is compatible with the output files of openQCD v=1.2, 1.4 and 1.6. Version can be specified as argument.</p><p>Examples:</p><pre><code class="language- hljs">read_ms1(path)
 read_ms1(path, v=&quot;1.4&quot;)
 read_ms1(path, v=&quot;1.6&quot;)</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article><article class="docstring"><header><a class="docstring-binding" id="juobs.read_md" href="#juobs.read_md"><code>juobs.read_md</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">read_md(path::String)</code></pre><p>Reads openQCD  pbp.dat files at a given path. This method returns a matrix <code>md[irw, icfg]</code> that contains the derivatives <span>$dS/dm$</span>, where <span>$md[irw=1] = dS/dm_l$</span> and <span>$md[irw=2] = dS/dm_s$</span></p><p><span>$Seff = -tr(log(D+m))$</span></p><p><span>$dSeff/ dm = -tr((D+m)^-1)$</span></p><p>Examples:</p><pre><code class="language- hljs">md = read_md(path)</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article><article class="docstring"><header><a class="docstring-binding" id="juobs.truncate_data!" href="#juobs.truncate_data!"><code>juobs.truncate_data!</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">truncate_data!(data::YData, nc::Int64)
 truncate_data!(data::Vector{YData}, nc::Vector{Int64})
 truncate_data!(data::Vector{CData}, nc::Int64)
 truncate_data!(data::Vector{Vector{CData}}, nc::Vector{Int64})</code></pre><p>Truncates the output of <code>read_mesons</code> and <code>read_ms</code> taking the first <code>nc</code> configurations.</p><p>Examples:</p><pre><code class="language- hljs">#Single replica
 dat = read_mesons(path, &quot;G5&quot;, &quot;G5&quot;)
 Y = read_ms(path)
 truncate_data!(dat, nc)
 truncate_data!(Y, nc)
 #Two replicas
 dat = read_mesons([path1, path2], &quot;G5&quot;, &quot;G5&quot;)
 Y = read_ms.([path1_ms, path2_ms])
 truncate_data!(dat, [nc1, nc2])
 truncate_data!(Y, [nc1, nc2])</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="index.html">« Home</a><a class="docs-footer-nextpage" href="tools.html">Tools »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 0.27.15 on <span class="colophon-date" title="Friday 1 July 2022 12:33">Friday 1 July 2022</span>. Using Julia version 1.6.5.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>
--- a/docs/build/search.html
+++ b/docs/build/search.html
 <!DOCTYPE html>
 <html lang="en"><head><meta charset="UTF-8"/><meta name="viewport" content="width=device-width, initial-scale=1.0"/><title>Search · juobs Documentation</title><script data-outdated-warner src="assets/warner.js"></script><link href="https://cdnjs.cloudflare.com/ajax/libs/lato-font/3.0.0/css/lato-font.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/juliamono/0.039/juliamono-regular.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.3/css/fontawesome.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.3/css/solid.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.3/css/brands.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.13.11/katex.min.css" rel="stylesheet" type="text/css"/><script>documenterBaseURL="."</script><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.6/require.min.js" data-main="assets/documenter.js"></script><script src="siteinfo.js"></script><script src="../versions.js"></script><link class="docs-theme-link" rel="stylesheet" type="text/css" href="assets/themes/documenter-dark.css" data-theme-name="documenter-dark" data-theme-primary-dark/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="assets/themes/documenter-light.css" data-theme-name="documenter-light" data-theme-primary/><script src="assets/themeswap.js"></script></head><body><div id="documenter"><nav class="docs-sidebar"><div class="docs-package-name"><span class="docs-autofit"><a href="index.html">juobs Documentation</a></span></div><form class="docs-search" action="search.html"><input class="docs-search-query" id="documenter-search-query" name="q" type="text" placeholder="Search docs"/></form><ul class="docs-menu"><li><a class="tocitem" href="index.html">Home</a></li><li><a class="tocitem" href="reader.html">Reader</a></li><li><a class="tocitem" href="tools.html">Tools</a></li><li><a class="tocitem" href="obs.html">Observables</a></li><li><a class="tocitem" href="linalg.html">Linear Algebra</a></li></ul><div class="docs-version-selector field has-addons"><div class="control"><span class="docs-label button is-static is-size-7">Version</span></div><div class="docs-selector control is-expanded"><div class="select is-fullwidth is-size-7"><select id="documenter-version-selector"></select></div></div></div></nav><div class="docs-main"><header class="docs-navbar"><nav class="breadcrumb"><ul class="is-hidden-mobile"><li class="is-active"><a href="search.html">Search</a></li></ul><ul class="is-hidden-tablet"><li class="is-active"><a href="search.html">Search</a></li></ul></nav><div class="docs-right"><a class="docs-settings-button fas fa-cog" id="documenter-settings-button" href="#" title="Settings"></a><a class="docs-sidebar-button fa fa-bars is-hidden-desktop" id="documenter-sidebar-button" href="#"></a></div></header><article><p id="documenter-search-info">Loading search...</p><ul id="documenter-search-results"></ul></article><nav class="docs-footer"><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 0.27.15 on <span class="colophon-date" title="Friday 1 July 2022 12:33">Friday 1 July 2022</span>. Using Julia version 1.6.5.</p></section><footer class="modal-card-foot"></footer></div></div></div></body><script src="search_index.js"></script><script src="assets/search.js"></script></html>
--- a/docs/build/search_index.js
+++ b/docs/build/search_index.js
 var documenterSearchIndex = {"docs":
 [{"location":"reader.html#Reader","page":"Reader","title":"Reader","text":"","category":"section"},{"location":"reader.html","page":"Reader","title":"Reader","text":"read_mesons\nread_ms\nread_ms1\nread_md\ntruncate_data!","category":"page"},{"location":"reader.html#juobs.read_mesons","page":"Reader","title":"juobs.read_mesons","text":"read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{String, Nothing}=nothing, legacy::Bool=false)\n\nread_mesons(path::Vector{String}, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{String, Nothing}=nothing, legacy::Bool=false)\n\nThis function read a mesons dat file at a given path and returns a vector of CData structures for different masses and Dirac structures. Dirac structures g1 and/or g2 can be passed as string arguments in order to filter correaltors. ADerrors id can be specified as argument. If is not specified, the id is fixed according to the ensemble name (example: \"H400\"-> id = \"H400\")\n\n*For the old version (without smearing, distance preconditioning and theta) set legacy=true.\n\nExamples:\n\nread_mesons(path)\nread_mesons(path, \"G5\")\nread_mesons(path, nothing, \"G5\")\nread_mesons(path, \"G5\", \"G5\")\nread_mesons(path, \"G5\", \"G5\", id=\"H100\")\nread_mesons(path, \"G5_d2\", \"G5_d2\", legacy=true)\n\n\n\n\n\n","category":"function"},{"location":"reader.html#juobs.read_ms","page":"Reader","title":"juobs.read_ms","text":"read_ms(path::String; id::Union{String, Nothing}=nothing, dtr::Int64=1, obs::String=\"Y\")\n\nReads openQCD ms dat files at a given path. This method return YData: \n\nt(t): flow time values\nobs(icfg, x0, t): the time-slice sums of the densities of the observable (Wsl, Ysl or Qsl)\nvtr: vector that contains trajectory number\nid: ensmble id\n\ndtr = dtr_cnfg / dtr_ms, where dtr_cnfg is the number of trajectories computed before saving the configuration. dtr_ms is the same but applied to the ms.dat file.\n\nExamples:\n\nY = read_ms(path)\n\n\n\n\n\n","category":"function"},{"location":"reader.html#juobs.read_ms1","page":"Reader","title":"juobs.read_ms1","text":"read_ms1(path::String; v::String=\"1.2\")\n\nReads openQCD ms1 dat files at a given path. This method returns a matrix W[irw, icfg] that contains the reweighting factors, where irw is the rwf index and icfg the configuration number. The function is compatible with the output files of openQCD v=1.2, 1.4 and 1.6. Version can be specified as argument.\n\nExamples:\n\nread_ms1(path)\nread_ms1(path, v=\"1.4\")\nread_ms1(path, v=\"1.6\")\n\n\n\n\n\n","category":"function"},{"location":"reader.html#juobs.read_md","page":"Reader","title":"juobs.read_md","text":"read_md(path::String)\n\nReads openQCD  pbp.dat files at a given path. This method returns a matrix md[irw, icfg] that contains the derivatives dSdm, where mdirw=1 = dSdm_l and mdirw=2 = dSdm_s\n\nSeff = -tr(log(D+m))\n\ndSeff dm = -tr((D+m)^-1)\n\nExamples:\n\nmd = read_md(path)\n\n\n\n\n\n","category":"function"},{"location":"reader.html#juobs.truncate_data!","page":"Reader","title":"juobs.truncate_data!","text":"truncate_data!(data::YData, nc::Int64)\n\ntruncate_data!(data::Vector{YData}, nc::Vector{Int64})\n\ntruncate_data!(data::Vector{CData}, nc::Int64)\n\ntruncate_data!(data::Vector{Vector{CData}}, nc::Vector{Int64})\n\nTruncates the output of read_mesons and read_ms taking the first nc configurations.\n\nExamples:\n\n#Single replica\ndat = read_mesons(path, \"G5\", \"G5\")\nY = read_ms(path)\ntruncate_data!(dat, nc)\ntruncate_data!(Y, nc)\n\n#Two replicas\ndat = read_mesons([path1, path2], \"G5\", \"G5\")\nY = read_ms.([path1_ms, path2_ms])\ntruncate_data!(dat, [nc1, nc2])\ntruncate_data!(Y, [nc1, nc2])\n\n\n\n\n\n","category":"function"},{"location":"linalg.html#Linear-Algebra","page":"Linear Algebra","title":"Linear Algebra","text":"","category":"section"},{"location":"linalg.html","page":"Linear Algebra","title":"Linear Algebra","text":"uweigvals\nuweigvecs\nuweigen\nget_matrix\nenergies\ngetall_eigvals\ngetall_eigvecs","category":"page"},{"location":"linalg.html#juobs.uweigvals","page":"Linear Algebra","title":"juobs.uweigvals","text":"uweigvals(a::Matrix{uwreal}; iter = 30)\n\nuweigvals(a::Matrix{uwreal}, b::Matrix{uwreal};  iter = 30)\n\nThis function computes the eigenvalues of a matrix of uwreal objects. If a second matrix b is given as input, it returns the generalised eigenvalues instead. It takes as input:\n\na::Matrix{uwreal} : a matrix of uwreal\nb::Matrix{uwreal} : a matrix of uwreal, optional\niter=30: optional flag to set the iterations of the qr algorithm used to solve the eigenvalue problem\n\nIt returns:\n\nres = Vector{uwreal}: a vector where each elements is an eigenvalue \n\na = Matrix{uwreal}(nothing, n,n) ## n*n matrix of uwreal with nothing entries\nb = Matrix{uwreal}(nothing, n,n) ## n*n matrix of uwreal with nothing entries\n\nres = uweigvals(a)  ##eigenvalues\nres1 = uweigvals(a,b) ## generalised eigenvalues\n\n\n\n\n\n","category":"function"},{"location":"linalg.html#juobs.uweigvecs","page":"Linear Algebra","title":"juobs.uweigvecs","text":"uweigvecs(a::Matrix{uwreal}; iter = 30)\n\nuweigvecs(a::Matrix{uwreal}, b::Matrix{uwreal};  iter = 30)\n\nThis function computes the eigenvectors of a matrix of uwreal objects. If a second matrix b is given as input, it returns the generalised eigenvectors instead. It takes as input:\n\na::Matrix{uwreal} : a matrix of uwreal\nb::Matrix{uwreal} : a matrix of uwreal, optional\niter=30 : the number of iterations of the qr algorithm used to extract the eigenvalues \n\nIt returns:\n\nres = Matrix{uwreal}: a matrix where each column is an eigenvector \n\nExamples:\n\na = Matrix{uwreal}(nothing, n,n) ## n*n matrix of uwreal with nothing entries\nb = Matrix{uwreal}(nothing, n,n) ## n*n matrix of uwreal with nothing entries\n\nres = uweigvecs(a)  ##eigenvectors in column \nres1 = uweigvecs(a,b) ## generalised eigenvectors in column \n\n\n\n\n\n","category":"function"},{"location":"linalg.html#juobs.uweigen","page":"Linear Algebra","title":"juobs.uweigen","text":"uweigen(a::Matrix{uwreal}; iter = 30)\n\nuweigen(a::Matrix{uwreal}, b::Matrix{uwreal};  iter = 30)\n\nThis function computes the eigenvalues and the eigenvectors  of a matrix of uwreal objects. If a second matrix b is given as input, it returns the generalised eigenvalues and eigenvectors instead. It takes as input:\n\na::Matrix{uwreal} : a matrix of uwreal\nb::Matrix{uwreal} : a matrix of uwreal, optional\niter=30 : the number of iterations of the qr algorithm used to extract the eigenvalues \n\nIt returns:\n\nevals = Vector{uwreal}: a vector where each elements is an eigenvalue \nevecs  = Matrix{uwreal}: a matrix where the i-th column is the eigenvector of the i-th eigenvalue\n\nExamples:\n\na = Matrix{uwreal}(nothing, n,n) ## n*n matrix of uwreal with nothing entries\nb = Matrix{uwreal}(nothing, n,n) ## n*n matrix of uwreal with nothing entries\n\neval, evec = uweigen(a)   \neval1, evec1 = uweigvecs(a,b) \n\n\n\n\n\n","category":"function"},{"location":"linalg.html#juobs.get_matrix","page":"Linear Algebra","title":"juobs.get_matrix","text":"get_matrix(diag::Vector{Array}, upper::Vector{Array} )\n\nThis function allows the user to build an array of matrices, where each matrix is a  symmetric matrix of correlators at a given timeslice.   \n\nIt takes as input:\n\ndiag::Vector{Array}:  vector of correlators. Each correlator will constitute a diagonal element of the matrices A[i] where i runs over the timeslices, i.e.                A[i][1,1] = diag[1], .... A[i][n,n] = diag[n]                Given n=length(diag), the matrices will have dimension n*n \nupper::Vector{Array}:  vector of correlators liying on the upper diagonal position.               A[i][1,2] = upper[1], .. , A[i][1,n] = upper[n-1],               A[i][2,3] = upper[n], .. , A[i][n-1,n] = upper[n(n-1)/2]                Given n, length(upper)=n(n-1)/2\n\nThe method returns an array of symmetric matrices of dimension n for each timeslice  \n\nExamples:\n\n## load data \npp_data = read_mesons(path, \"G5\", \"G5\")\npa_data = read_mesons(path, \"G5\", \"G0G5\")\naa_data = read_mesons(path, \"G0G5\", \"G0G5\")\n\n## create Corr struct\ncorr_pp = corr_obs.(pp_data)\ncorr_pa = corr_obs.(pa_data)\ncorr_aa = corr_obs.(aa_data) # array of correlators for different \\mu_q combinations\n\n## set up matrices\ncorr_diag = [corr_pp[1], corr_aa[1]] \ncorr_upper = [corr_pa[1]]\n\nmatrices = [corr_diag, corr_upper]\n\nJulia> matrices[i]\n       pp[i]  ap[i]\n       pa[i]  aa[i]\n\n## where i runs over the timeslices\n\n\n\n\n\n","category":"function"},{"location":"linalg.html#juobs.energies","page":"Linear Algebra","title":"juobs.energies","text":"energies(evals::Vector{Array}; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)\n\nThis method computes the energy level from the eigenvalues according to:\n\nE_i(t) = log(λ(t)  λ(t+1))\n\nwhere i=1,..,n with n=length(evals[1]) and t=1,..,T  total time slices. It returns a vector array en where each entry en[i][t] contains the i-th states energy at time t \n\nExamples:\n\n## load data\npp_data = read_mesons(path, \"G5\", \"G5\")\npa_data = read_mesons(path, \"G5\", \"G0G5\")\naa_data = read_mesons(path, \"G0G5\", \"G0G5\")\n\n## create Corr struct\ncorr_pp = corr_obs.(pp_data)\ncorr_pa = corr_obs.(pa_data)\ncorr_aa = corr_obs.(aa_data) # array of correlators for different \\mu_q combinations\n\n## set up matrices \ncorr_diag = [corr_pp[1], corr_aa[1]] \ncorr_upper = [corr_pa[1]]\n\nmatrices = [corr_diag, corr_upper]\n\n## solve the GEVP\nevals = getall_eigvals(matrices, 5) #where t_0=5\nen = energies(evals)\n\nJulia> en[i] # i-th state energy at each timeslice\n\n\n\n\n\n","category":"function"},{"location":"linalg.html#juobs.getall_eigvals","page":"Linear Algebra","title":"juobs.getall_eigvals","text":"getall_eigvals(a::Vector{Matrix}, t0; iter=30 )\n\nThis function solves a GEVP problem, returning the eigenvalues, for a list of matrices, taking as generalised matrix the one  at index t0, i.e:\n\nC(t_i)v_i = λ_i C(t_0) v_i,   with   i=1:lenght(a)\n\nIt takes as input:\n\na::Vector{Matrix} : a vector of matrices\nt0::Int64 : idex value at which the fixed matrix is taken\niter=30 : the number of iterations of the qr algorithm used to extract the eigenvalues \n\nIt returns:\n\nres = Vector{Vector{uwreal}}\n\nwhere res[i] are the generalised eigenvalues of the i-th matrix of the input array. \n\nExamples:\n\n## load data\npp_data = read_mesons(path, \"G5\", \"G5\")\npa_data = read_mesons(path, \"G5\", \"G0G5\")\naa_data = read_mesons(path, \"G0G5\", \"G0G5\")\n\n## create Corr struct\ncorr_pp = corr_obs.(pp_data)\ncorr_pa = corr_obs.(pa_data)\ncorr_aa = corr_obs.(aa_data) # array of correlators for different \\mu_q combinations\n\n## set up matrices \ncorr_diag = [corr_pp[1], corr_aa[1]] \ncorr_upper = [corr_pa[1]]\n\nmatrices = [corr_diag, corr_upper]\n\n## solve the GEVP\n#evals = getall_eigvals(matrices, 5) #where t_0=5\n\n\nJulia>\n\n\n\n\n\n","category":"function"},{"location":"linalg.html#juobs.getall_eigvecs","page":"Linear Algebra","title":"juobs.getall_eigvecs","text":"getall_eigvecs(a::Vector{Matrix}, delta_t; iter=30 )\n\nThis function solves a GEVP problem, returning the eigenvectors, for a list of matrices.\n\nC(t_i)v_i = λ_i C(t_i-delta_t) v_i, with i=1:lenght(a)\n\nHere delta_t is the time shift within the two matrices of the problem, and is kept fixed. It takes as input:\n\na::Vector{Matrix} : a vector of matrices\ndelta_t::Int64 : the fixed time shift t-t_0\niter=30 : the number of iterations of the qr algorithm used to extract the eigenvalues \n\nIt returns:\n\nres = Vector{Matrix{uwreal}}\n\nwhere each res[i] is a matrix  with the eigenvectors as  columns Examples:\n\nmat_array = get_matrix(diag, upper_diag)\nevecs = getall_eigvecs(mat_array, 5)\n\n\n\n\n\n","category":"function"},{"location":"obs.html#Observables","page":"Observables","title":"Observables","text":"","category":"section"},{"location":"obs.html","page":"Observables","title":"Observables","text":"meff\nmpcac\ndec_const\ndec_const_pcvc\ncomp_t0","category":"page"},{"location":"obs.html#juobs.meff","page":"Observables","title":"juobs.meff","text":"meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)     \n\nmeff(corr::Corr, plat::Vector{Int64}; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)\n\nComputes effective mass for a given correlator corr at a given plateau plat. Correlator can be passed as an Corr struct or Vector{uwreal}.\n\nThe flags pl and data allow to show the plots and return data as an extra result.\n\ndata = read_mesons(path, \"G5\", \"G5\")\ncorr_pp = corr_obs.(data)\nm = meff(corr_pp[1], [50, 60], pl=false)\n\n\n\n\n\n","category":"function"},{"location":"obs.html#juobs.mpcac","page":"Observables","title":"juobs.mpcac","text":"mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)\n\nmpcac(a0p::Corr, pp::Corr, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)\n\nComputes the bare PCAC mass for a given correlator a0p and pp at a given plateau plat. Correlator can be passed as an Corr struct or Vector{uwreal}.\n\nThe flags pl and data allow to show the plots and return data as an extra result. The ca variable allows to compute mpcac using the improved axial current.\n\ndata_pp = read_mesons(path, \"G5\", \"G5\")\ndata_a0p = read_mesons(path, \"G5\", \"G0G5\")\ncorr_pp = corr_obs.(data_pp)\ncorr_a0p = corr_obs.(data_a0p)\nm12 = mpcac(corr_a0p, corr_pp, [50, 60], pl=false)\n\np0 = 9.2056\np1 = -13.9847\ng2 = 1.73410\nca = -0.006033 * g2 *( 1 + exp(p0 + p1/g2))\n\nm12 = mpcac(corr_a0p, corr_pp, [50, 60], pl=false, ca=ca)\n\n\n\n\n\n","category":"function"},{"location":"obs.html#juobs.dec_const","page":"Observables","title":"juobs.dec_const","text":"dec_const(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)\n\ndec_const(a0p::Corr, pp::Corr, plat::Vector{Int64}, m::uwreal; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)\n\ndec_const(a0pL::Vector{uwreal}, a0pR::Vector{uwreal}, ppL::Vector{uwreal}, ppR::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)\n\ndec_const(a0pL::Corr, a0pR::Corr, ppL::Corr, ppR::Corr, plat::Vector{Int64}, m::uwreal; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)\n\nComputes the bare decay constant using A_0P and PP correlators . The decay constant is computed in the plateau plat. Correlator can be passed as an Corr struct or Vector{uwreal}. If it is passed as a uwreal vector, effective mass m and source position y0 must be specified.\n\nThe flags pl and data allow to show the plots and return data as an extra result. The ca variable allows to compute dec_const using the improved axial current.\n\nThe method assumes that the source is close to the boundary. It takes the following ratio to cancel boundary effects. R = fracf_A(x_0 y_0)sqrtf_P(T-y_0 y_0) * e^m (x_0 - T2)\n\nIf left and right correlators are included in the input. The result is computed with the following ratio R = sqrtf_A(x_0 y_0) * f_A(x_0 T - 1 - y_0)  f_P(T - 1 - y_0 y_0)\n\ndata_pp = read_mesons(path, \"G5\", \"G5\", legacy=true)\ndata_a0p = read_mesons(path, \"G5\", \"G0G5\", legacy=true)\n\ncorr_pp = corr_obs.(data_pp, L=32)\ncorr_a0p = corr_obs.(data_a0p, L=32)\n\ncorr_a0pL, corr_a0pR = [corr_a0p[1], corr_a0p[2]]\ncorr_ppL, corr_ppR = [corr_pp[1], corr_pp[2]]\n\nm = meff(corr_pp[1], [50, 60], pl=false)\n\nbeta = 3.46\np0 = 9.2056\np1 = -13.9847\ng2 = 6 / beta\nca = -0.006033 * g2 *( 1 + exp(p0 + p1/g2))\n\nf = dec_const(corr_a0p[1], corr_pp[1], [50, 60], m, pl=true, ca=ca)\n\nf_ratio = dec_const(corr_a0pL, corr_a0pR, corr_ppL, corr_ppR, [50, 60], m, pl=true, ca=ca)\n\n\n\n\n\n","category":"function"},{"location":"obs.html#juobs.dec_const_pcvc","page":"Observables","title":"juobs.dec_const_pcvc","text":"dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)\n\ndec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)\n\ndec_const_pcvc(ppL::Vector{uwreal}, ppR::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)\n\ndec_const_pcvc(corrL::Corr, corrR::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)\n\nComputes decay constant using the PCVC relation for twisted mass fermions. The decay constant is computed in the plateau plat. Correlator can be passed as an Corr struct or Vector{uwreal}. If it is passed as a uwreal vector, vector of twisted masses mu and source position y0 must be specified.\n\nThe flags pl and data allow to show the plots and return data as an extra result.\n\nThe method extract the matrix element assuming that the source is in the bulk. ** **If left and right correlators are included in the input. The result is computed with a ratio that cancels boundary effects: R = sqrtf_P(x_0 y_0) * f_P(x_0 T - 1 - y_0)  f_P(T - 1 - y_0 y_0)\n\ndata = read_mesons(path, \"G5\", \"G5\")\ncorr_pp = corr_obs.(data, L=32)\nm = meff(corr_pp[1], [50, 60], pl=false)\nf = dec_const_pcvc(corr_pp[1], [50, 60], m, pl=false)\n\n#left and right correlators\nf_ratio = dec_const_pcvc(ppL, ppR, [50, 60], m)\n\n\n\n\n\n","category":"function"},{"location":"obs.html#juobs.comp_t0","page":"Observables","title":"juobs.comp_t0","text":"comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Matrix{Float64}, Nothing}=nothing, npol::Int64=2, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, info::Bool=false)\n\ncomp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Vector{Matrix{Float64}}, Nothing}=nothing, npol::Int64=2, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, info::Bool=false)\n\nComputes t0 using the energy density of the action Ysl(Yang-Mills action). t0 is computed in the plateau plat. A polynomial interpolation in t is performed to find t0, where npol is the degree of the polynomial (linear fit by default)\n\nThe flag pl allows to show the plot. \n\nThe flag info provides extra output that contains information about the primary observables. The function returns the primary observables WY and W (it returns the observable <Y> if rw=nothing)\n\n#Single replica\nY = read_ms(path)\nrw = read_ms(path_rw)\n\nt0, Yobs = comp_t0(Y, [38, 58], L=32, info=true)\nt0_r, WYobs, Wobs = comp_t0(Y, [38, 58], L=32, rw=rw, info=true)\n\n#Two replicas\nY1 = read_ms(path1)\nY2 = read_ms(path2)\nrw1 = read_ms(path_rw1)\nrw2 = read_ms(path_rw2)\n\nt0 = comp_t0([Y1, Y2], [38, 58], L=32, pl=true)\nt0_r = comp_t0(Y, [38, 58], L=32, rw=[rw1, rw2], pl=true)\n\n\n\n\n\n\n","category":"function"},{"location":"index.html#DOCUMENTATION","page":"Home","title":"DOCUMENTATION","text":"","category":"section"},{"location":"index.html#Contents","page":"Home","title":"Contents","text":"","category":"section"},{"location":"index.html","page":"Home","title":"Home","text":"Pages = [\"reader.md\", \"tools.md\", \"obs.md\", \"linalg.md\"]\nDepth = 3","category":"page"},{"location":"tools.html#Tools","page":"Tools","title":"Tools","text":"","category":"section"},{"location":"tools.html","page":"Tools","title":"Tools","text":"corr_obs\ncorr_sym\nmd_sea\nmd_val\nlin_fit\nfit_routine\nbayesian_av","category":"page"},{"location":"tools.html#juobs.corr_obs","page":"Tools","title":"juobs.corr_obs","text":"corr_obs(cdata::CData; real::Bool=true, rw::Union{Array{Float64, 2}, Nothing}=nothing, L::Int64=1, info::Bool=false)\n\ncorr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array{Float64, 2}, 1}, Nothing}=nothing, L::Int64=1, info::Bool=false)\n\nCreates a Corr struct with the given CData struct cdata (read_mesons) for a single replica. An array of CData can be passed as argument for multiple replicas.\n\nThe flag real select the real or imaginary part of the correlator. If rw is specified, the method applies reweighting. rw is passed as a matrix of Float64 (read_ms1) The correlator can be normalized with the volume factor if L is fixed.\n\nThe flag info provides extra output that contains information about the primary observables. The function returns the primary observables WO and W (it returns the observable <O> if rw=nothing)\n\n#Single replica\ndata = read_mesons(path, \"G5\", \"G5\")\nrw = read_ms1(path_rw)\ncorr_pp = corr_obs.(data)\ncorr_pp_r = corr_obs.(data, rw=rw)\n\n#Single replica + Info\ndata = read_mesons(path, \"G5\", \"G5\")\nrw = read_ms1(path_rw)\ncorr_pp, O = corr_obs(data[1], info=true)\ncorr_pp_r, WO, W = corr_obs(data[1], rw=rw, info=true)\n\n#Two replicas\ndata = read_mesons([path_r1, path_r2], \"G5\", \"G5\")\nrw1 = read_ms1(path_rw1)\nrw2 = read_ms1(path_rw2)\n\ncorr_pp = corr_obs.(data)\ncorr_pp_r = corr_obs.(data, rw=[rw1, rw2])\n\n\n\n\n\n","category":"function"},{"location":"tools.html#juobs.corr_sym","page":"Tools","title":"juobs.corr_sym","text":"corr_sym(corrL::Corr, corrR::Corr, parity::Int64=1)\n\nComputes the symmetrized correlator using the left correlador corrL and the right correlator corrR. The source position of corrR must be T - 1 - y0, where y0 is the source position of corrL. \n\npp_sym = corr_sym(ppL, ppR, +1)\na0p_sym = corr_sym(a0pL, a0pR, -1)\n\n\n\n\n\n","category":"function"},{"location":"tools.html#juobs.md_sea","page":"Tools","title":"juobs.md_sea","text":"md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ow::uwreal, w::Union{uwreal, Nothing}=nothing, ws::ADerrors.wspace=ADerrors.wsg)\n\nmd_sea(a::uwreal, md::Vector{Matrix{Float64}}, ow::Array{uwreal}, w::Union{uwreal, Nothing}=nothing, ws::ADerrors.wspace=ADerrors.wsg)\n\nComputes the derivative of an observable A with respect to the sea quark masses.\n\nfracd Adm(sea) = sum_i fracpartial Apartial O_i  fracd O_id m(sea)\n\nfracd O_idm(sea) = O_i fracpartial Spartial m - O_i fracpartial Spartial m\n\nwhere O_i are primary observables. The function returns  fracpartial Apartial O_i  fracd O_id m(sea), where the primary observable O_i is specified as an input argument.\n\now is a primary observable or a vector of primary observables that contains OW (O if reweighting is not applied) and  w is a primary obserable that contains the reweighting factors W (w=nothing if reweighting is not applied).\n\nmd is a vector that contains the derivative of the action S with respect to the sea quark masses for each replica. md[irep][irw, icfg]\n\nmd_sea returns a tuple of uwreal observables (dAdm_l dAdm_s)_sea,  where m_l and m_s are the light and strange quark masses.\n\n#Single replica\ndata = read_mesons(path, \"G5\", \"G5\")\nmd = read_md(path_md)\nrw = read_ms1(path_rw)\n\ncorr_pp, wpp, w = corr_obs.(data[1], rw=rw, info=true)\nm = meff(corr_pp[1], plat)\nm_mdl, m_mds = md_sea(m, [md], wpp, w, ADerrors.wsg)\nm_shifted = m + 2 * dml * m_mdl + dms * m_mds\n\n#Two replicas\ndata = read_mesons([path_r1, path_r2], \"G5\", \"G5\")\nmd1 = read_md(path_md1)\nmd2 = read_md(path_md2)\n\ncorr_pp, pp = corr_obs(data[1], info=true)\nm = meff(corr_pp[1], plat)\nm_mdl, m_mds = md_sea(m, [md1, md2], pp, ADerrors.wsg)\nm_shifted = m + 2 * dml * m_mdl + dms * m_mds\n\n\n\n\n\n","category":"function"},{"location":"tools.html#juobs.md_val","page":"Tools","title":"juobs.md_val","text":"md_val(a::uwreal, obs::Corr, derm::Vector{Corr})\n\nComputes the derivative of an observable A with respect to the valence quark masses.\n\nfracd Adm(val) = sum_i fracpartial Apartial O_i  fracd O_id m(val)\n\nfracd O_idm(val) = fracpartial O_ipartial m(val)\n\nwhere O_i are primary observables \n\nmd is a vector that contains the derivative of the action S with respect to the sea quark masses for each replica. md[irep][irw, icfg]\n\nmd_val returns a tuple of uwreal observables (dAdm_1 dAdm_2)_val,  where m_1 and m_2 are the correlator masses.\n\ndata = read_mesons(path, \"G5\", \"G5\", legacy=true)\ndata_d1 = read_mesons(path, \"G5_d1\", \"G5_d1\", legacy=true)\ndata_d2 = read_mesons(path, \"G5_d2\", \"G5_d2\", legacy=true)\n\nrw = read_ms1(path_rw)\n\ncorr_pp = corr_obs.(data, rw=rw)\ncorr_pp_d1 = corr_obs.(data_d1, rw=rw)\ncorr_pp_d2 = corr_obs.(data_d2, rw=rw)\nderm = [[corr_pp_d1[k], corr_pp_d2[k]] for k = 1:length(pp_d1)]\n\nm = meff(corr_pp[1], plat)\nm_md1, m_md2 = md_val(m, corr_pp[1], derm[1])\nm_shifted = m + dm1 * m_md1 + dm2 * m_md2\n\n\n\n\n\n","category":"function"},{"location":"tools.html#juobs.lin_fit","page":"Tools","title":"juobs.lin_fit","text":"lin_fit(x::Vector{<:Real}, y::Vector{uwreal})\n\nComputes a linear fit of uwreal data points y. This method return uwreal fit parameters and chisqexpected.\n\nfitp, csqexp = lin_fit(phi2, m2)\nm2_phys = fitp[1] + fitp[2] * phi2_phys\n\n\n\n\n\n","category":"function"},{"location":"tools.html#juobs.fit_routine","page":"Tools","title":"juobs.fit_routine","text":"fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}, param::Int64=3; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)\n\nfit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}, param::Int64=3; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, covar::Bool=false)\n\nGiven a model function with a number param of parameters and an array of uwreal, this function fit ydata with the given model and print fit information The method return an array upar with the best fit parameters with their errors. The flag wpm is an optional array of Float64 of lenght 4. The first three paramenters specify the criteria to determine the summation windows:\n\nvp[1]: The autocorrelation function is summed up to t = round(vp1).\nvp[2]: The sumation window is determined using U. Wolff poposal with S_tau = wpm2\nvp[3]: The autocorrelation function Gamma(t) is summed up a point where its error deltaGamma(t) is a factor vp[3] times larger than the signal.\n\nAn additional fourth parameter vp[4], tells ADerrors to add a tail to the error with tau_exp = wpm4. Negative values of wpm[1:4] are ignored and only one component of wpm[1:3] needs to be positive. If the flag covaris set to true, fit_routine takes into account covariances between x and y for each data point.\n\n@. model(x,p) = p[1] + p[2] * exp(-(p[3]-p[1])*x)\n@. model2(x,p) = p[1] + p[2] * x[:, 1] + (p[3] + p[4] * x[:, 1]) * x[:, 2]\nfit_routine(model, xdata, ydata, param=3)\nfit_routine(model, xdata, ydata, param=3, covar=true)\n\n\n\n\n\n","category":"function"},{"location":"tools.html#juobs.bayesian_av","page":"Tools","title":"juobs.bayesian_av","text":"bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64}, tmax_array::Array{Int64}, k::Int64, pl::Bool, data::Bool; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)\n\nbayesian_av(fun1::Function, fun2::Function, y::Array{uwreal}, tmin_array::Array{Int64}, tmax_array::Array{Int64}, k1::Int64, k2::Int64, pl::Bool, data::Bool; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)\n\nbayesian_av(fun::Array{Function}, y::Array{uwreal}, tmin_array::Array{Int64}, tmax_array::Array{Int64}, k::Array{Int64}, pl::Bool, data::Bool; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)\n\nComputes bayesian average of data. For a given fit function, it explores choices of fit intervals, assigning each of them a weight. The function saves the first fit parameter of your function, and then it does the weighted average of it and assigns a systematic. See https://arxiv.org/abs/2008.01069 \n\nThe function takes as input the fit intervals to explore. \n\ntmin_array is an array of integers with the lower bounds on the fit intervals to explore, ***ordered from lower to higher***.\n\ntmax_array is an array of integers with the upper bounds on the fit intervals to explore, ***ordered from lower to higher***.\n\nk is the number of parameters of the fit function to use.\n\nYou can also use as input two fit functions, and two values of k, one for each function. Then, for each fit interval choice, the function explores the two fit functions. This means that for each fit interval choice you get two results: one for the first fit funcction, and another for the second. You can also use a vector of functions and a vector of k (numer of parameters of each funtion) to apply the bayesian averaging method to multiple functions.\n\nThe method returns two objects: first, the weighted average as an uwreal object, with mean value and statistichal error. The second object returned is the systematic error coming from the fit interval variation. If data is true, then returns 4 objects: weighted average, systematic error, a vector with the results of the fit for each fit interval choice, and a vector with the weights associated to each fit.\n\n@.fun(x,p) = p[1] * x ^0\nk = 1\ntmin_array = [10,11,12,13,14,15]\ntmax_array = [80,81,82,83,84,85]\n(average, systematics, data, weights) = bayesian_av(fun,x,tmin_array,tmax_array,k,pl=true,data=true)\n\n@.fun1(x,p) = p[1] * x ^0\n@.fun2(x,p) = p[1] + p[2] * exp( - p[3] * (x))\nk1 = 1\nk2 = 3\ntmin_array = [10,11,12,13,14,15]\ntmax_array = [80,81,82,83,84,85]\n(average, systematics) = bayesian_av(fun1,fun2,x,tmin_array,tmax_array,k1,k2)\n\n@.fun1(x,p) = p[1] * x ^0\n@.fun2(x,p) = p[1] + p[2] * exp( - p[3] * (x))\nk1 = 1\nk2 = 3\ntmin_array = [10,11,12,13,14,15]\ntmax_array = [80,81,82,83,84,85]\n(average, systematics) = bayesian_av([fun1,fun2],x,tmin_array,tmax_array,[k1,k2])\n\n\n\n\n\n","category":"function"}]
 }
--- a/docs/build/tools.html
+++ b/docs/build/tools.html
 <!DOCTYPE html>
 <html lang="en"><head><meta charset="UTF-8"/><meta name="viewport" content="width=device-width, initial-scale=1.0"/><title>Tools · juobs Documentation</title><script data-outdated-warner src="assets/warner.js"></script><link href="https://cdnjs.cloudflare.com/ajax/libs/lato-font/3.0.0/css/lato-font.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/juliamono/0.039/juliamono-regular.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.3/css/fontawesome.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.3/css/solid.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.3/css/brands.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.13.11/katex.min.css" rel="stylesheet" type="text/css"/><script>documenterBaseURL="."</script><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.6/require.min.js" data-main="assets/documenter.js"></script><script src="siteinfo.js"></script><script src="../versions.js"></script><link class="docs-theme-link" rel="stylesheet" type="text/css" href="assets/themes/documenter-dark.css" data-theme-name="documenter-dark" data-theme-primary-dark/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="assets/themes/documenter-light.css" data-theme-name="documenter-light" data-theme-primary/><script src="assets/themeswap.js"></script></head><body><div id="documenter"><nav class="docs-sidebar"><div class="docs-package-name"><span class="docs-autofit"><a href="index.html">juobs Documentation</a></span></div><form class="docs-search" action="search.html"><input class="docs-search-query" id="documenter-search-query" name="q" type="text" placeholder="Search docs"/></form><ul class="docs-menu"><li><a class="tocitem" href="index.html">Home</a></li><li><a class="tocitem" href="reader.html">Reader</a></li><li class="is-active"><a class="tocitem" href="tools.html">Tools</a></li><li><a class="tocitem" href="obs.html">Observables</a></li><li><a class="tocitem" href="linalg.html">Linear Algebra</a></li></ul><div class="docs-version-selector field has-addons"><div class="control"><span class="docs-label button is-static is-size-7">Version</span></div><div class="docs-selector control is-expanded"><div class="select is-fullwidth is-size-7"><select id="documenter-version-selector"></select></div></div></div></nav><div class="docs-main"><header class="docs-navbar"><nav class="breadcrumb"><ul class="is-hidden-mobile"><li class="is-active"><a href="tools.html">Tools</a></li></ul><ul class="is-hidden-tablet"><li class="is-active"><a href="tools.html">Tools</a></li></ul></nav><div class="docs-right"><a class="docs-edit-link" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs" title="Edit on GitLab"><span class="docs-icon fab"></span><span class="docs-label is-hidden-touch">Edit on GitLab</span></a><a class="docs-settings-button fas fa-cog" id="documenter-settings-button" href="#" title="Settings"></a><a class="docs-sidebar-button fa fa-bars is-hidden-desktop" id="documenter-sidebar-button" href="#"></a></div></header><article class="content" id="documenter-page"><h1 id="Tools"><a class="docs-heading-anchor" href="#Tools">Tools</a><a id="Tools-1"></a><a class="docs-heading-anchor-permalink" href="#Tools" title="Permalink"></a></h1><article class="docstring"><header><a class="docstring-binding" id="juobs.corr_obs" href="#juobs.corr_obs"><code>juobs.corr_obs</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">corr_obs(cdata::CData; real::Bool=true, rw::Union{Array{Float64, 2}, Nothing}=nothing, L::Int64=1, info::Bool=false)
 corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array{Float64, 2}, 1}, Nothing}=nothing, L::Int64=1, info::Bool=false)</code></pre><p>Creates a <code>Corr</code> struct with the given <code>CData</code> struct <code>cdata</code> (<code>read_mesons</code>) for a single replica. An array of <code>CData</code> can be passed as argument for multiple replicas.</p><p>The flag <code>real</code> select the real or imaginary part of the correlator. If <code>rw</code> is specified, the method applies reweighting. <code>rw</code> is passed as a matrix of Float64 (<code>read_ms1</code>) The correlator can be normalized with the volume factor if <code>L</code> is fixed.</p><p>The flag <code>info</code> provides extra output that contains information about the primary observables. The function returns the primary observables <span>$&lt;WO&gt;$</span> and <span>$&lt;W&gt;$</span> (it returns the observable &lt;O&gt; if rw=nothing)</p><pre><code class="language- hljs">#Single replica
 data = read_mesons(path, &quot;G5&quot;, &quot;G5&quot;)
 rw = read_ms1(path_rw)
 corr_pp = corr_obs.(data)
 corr_pp_r = corr_obs.(data, rw=rw)
 #Single replica + Info
 data = read_mesons(path, &quot;G5&quot;, &quot;G5&quot;)
 rw = read_ms1(path_rw)
 corr_pp, O = corr_obs(data[1], info=true)
 corr_pp_r, WO, W = corr_obs(data[1], rw=rw, info=true)
 #Two replicas
 data = read_mesons([path_r1, path_r2], &quot;G5&quot;, &quot;G5&quot;)
 rw1 = read_ms1(path_rw1)
 rw2 = read_ms1(path_rw2)
 corr_pp = corr_obs.(data)
 corr_pp_r = corr_obs.(data, rw=[rw1, rw2])</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article><article class="docstring"><header><a class="docstring-binding" id="juobs.corr_sym" href="#juobs.corr_sym"><code>juobs.corr_sym</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">corr_sym(corrL::Corr, corrR::Corr, parity::Int64=1)</code></pre><p>Computes the symmetrized correlator using the left correlador <code>corrL</code> and the right correlator <code>corrR</code>. The source position of <code>corrR</code> must be <code>T - 1 - y0</code>, where <code>y0</code> is the source position of <code>corrL</code>. </p><pre><code class="language- hljs">pp_sym = corr_sym(ppL, ppR, +1)
 a0p_sym = corr_sym(a0pL, a0pR, -1)</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article><article class="docstring"><header><a class="docstring-binding" id="juobs.md_sea" href="#juobs.md_sea"><code>juobs.md_sea</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ow::uwreal, w::Union{uwreal, Nothing}=nothing, ws::ADerrors.wspace=ADerrors.wsg)
 md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ow::Array{uwreal}, w::Union{uwreal, Nothing}=nothing, ws::ADerrors.wspace=ADerrors.wsg)</code></pre><p>Computes the derivative of an observable A with respect to the sea quark masses.</p><p><span>$\frac{d &lt;A&gt;}{dm(sea)} = \sum_i \frac{\partial &lt;A&gt;}{\partial &lt;O_i&gt;}  \frac{d &lt;O_i&gt;}{d m(sea)}$</span></p><p><span>$\frac{d &lt;O_i&gt;}{dm(sea)} = &lt;O_i&gt; &lt;\frac{\partial S}{\partial m}&gt; - &lt;O_i \frac{\partial S}{\partial m}&gt;$</span></p><p>where <span>$O_i$</span> are primary observables. The function returns  <span>$\frac{\partial &lt;A&gt;}{\partial &lt;O_i&gt;}  \frac{d &lt;O_i&gt;}{d m(sea)}$</span>, where the primary observable <span>$O_i$</span> is specified as an input argument.</p><p><code>ow</code> is a primary observable or a vector of primary observables that contains <span>$&lt;OW&gt;$</span> (<span>$&lt;O&gt;$</span> if reweighting is not applied) and  <code>w</code> is a primary obserable that contains the reweighting factors <span>$&lt;W&gt;$</span> (<code>w</code>=nothing if reweighting is not applied).</p><p><code>md</code> is a vector that contains the derivative of the action <span>$S$</span> with respect to the sea quark masses for each replica. <code>md[irep][irw, icfg]</code></p><p><code>md_sea</code> returns a tuple of uwreal observables <span>$(dA/dm_l, dA/dm_s)|_{sea}$</span>,  where <span>$m_l$</span> and <span>$m_s$</span> are the light and strange quark masses.</p><pre><code class="language- hljs">#Single replica
 data = read_mesons(path, &quot;G5&quot;, &quot;G5&quot;)
 md = read_md(path_md)
 rw = read_ms1(path_rw)
 corr_pp, wpp, w = corr_obs.(data[1], rw=rw, info=true)
 m = meff(corr_pp[1], plat)
 m_mdl, m_mds = md_sea(m, [md], wpp, w, ADerrors.wsg)
 m_shifted = m + 2 * dml * m_mdl + dms * m_mds
 #Two replicas
 data = read_mesons([path_r1, path_r2], &quot;G5&quot;, &quot;G5&quot;)
 md1 = read_md(path_md1)
 md2 = read_md(path_md2)
 corr_pp, pp = corr_obs(data[1], info=true)
 m = meff(corr_pp[1], plat)
 m_mdl, m_mds = md_sea(m, [md1, md2], pp, ADerrors.wsg)
 m_shifted = m + 2 * dml * m_mdl + dms * m_mds</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article><article class="docstring"><header><a class="docstring-binding" id="juobs.md_val" href="#juobs.md_val"><code>juobs.md_val</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">md_val(a::uwreal, obs::Corr, derm::Vector{Corr})</code></pre><p>Computes the derivative of an observable A with respect to the valence quark masses.</p><p><span>$\frac{d &lt;A&gt;}{dm(val)} = \sum_i \frac{\partial &lt;A&gt;}{\partial &lt;O_i&gt;}  \frac{d &lt;O_i&gt;}{d m(val)}$</span></p><p><span>$\frac{d &lt;O_i&gt;}{dm(val)} = &lt;\frac{\partial O_i}{\partial m(val)}&gt;$</span></p><p>where <span>$O_i$</span> are primary observables </p><p><code>md</code> is a vector that contains the derivative of the action <span>$S$</span> with respect to the sea quark masses for each replica. <code>md[irep][irw, icfg]</code></p><p><code>md_val</code> returns a tuple of <code>uwreal</code> observables <span>$(dA/dm_1, dA/dm_2)|_{val}$</span>,  where <span>$m_1$</span> and <span>$m_2$</span> are the correlator masses.</p><pre><code class="language- hljs">data = read_mesons(path, &quot;G5&quot;, &quot;G5&quot;, legacy=true)
 data_d1 = read_mesons(path, &quot;G5_d1&quot;, &quot;G5_d1&quot;, legacy=true)
 data_d2 = read_mesons(path, &quot;G5_d2&quot;, &quot;G5_d2&quot;, legacy=true)
 rw = read_ms1(path_rw)
 corr_pp = corr_obs.(data, rw=rw)
 corr_pp_d1 = corr_obs.(data_d1, rw=rw)
 corr_pp_d2 = corr_obs.(data_d2, rw=rw)
 derm = [[corr_pp_d1[k], corr_pp_d2[k]] for k = 1:length(pp_d1)]
 m = meff(corr_pp[1], plat)
 m_md1, m_md2 = md_val(m, corr_pp[1], derm[1])
 m_shifted = m + dm1 * m_md1 + dm2 * m_md2</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article><article class="docstring"><header><a class="docstring-binding" id="juobs.lin_fit" href="#juobs.lin_fit"><code>juobs.lin_fit</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">lin_fit(x::Vector{&lt;:Real}, y::Vector{uwreal})</code></pre><p>Computes a linear fit of uwreal data points y. This method return uwreal fit parameters and chisqexpected.</p><pre><code class="language- hljs">fitp, csqexp = lin_fit(phi2, m2)
 m2_phys = fitp[1] + fitp[2] * phi2_phys</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article><article class="docstring"><header><a class="docstring-binding" id="juobs.fit_routine" href="#juobs.fit_routine"><code>juobs.fit_routine</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">fit_routine(model::Function, xdata::Array{&lt;:Real}, ydata::Array{uwreal}, param::Int64=3; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
 fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}, param::Int64=3; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, covar::Bool=false)</code></pre><p>Given a model function with a number param of parameters and an array of <code>uwreal</code>, this function fit ydata with the given <code>model</code> and print fit information The method return an array <code>upar</code> with the best fit parameters with their errors. The flag <code>wpm</code> is an optional array of Float64 of lenght 4. The first three paramenters specify the criteria to determine the summation windows:</p><ul><li><p><code>vp[1]</code>: The autocorrelation function is summed up to <span>$t = round(vp[1])$</span>.</p></li><li><p><code>vp[2]</code>: The sumation window is determined using U. Wolff poposal with <span>$S_\tau = wpm[2]$</span></p></li><li><p><code>vp[3]</code>: The autocorrelation function <span>$\Gamma(t)$</span> is summed up a point where its error <span>$\delta\Gamma(t)$</span> is a factor <code>vp[3]</code> times larger than the signal.</p></li></ul><p>An additional fourth parameter <code>vp[4]</code>, tells ADerrors to add a tail to the error with <span>$\tau_{exp} = wpm[4]$</span>. Negative values of <code>wpm[1:4]</code> are ignored and only one component of <code>wpm[1:3]</code> needs to be positive. If the flag <code>covar</code>is set to true, <code>fit_routine</code> takes into account covariances between x and y for each data point.</p><pre><code class="language- hljs">@. model(x,p) = p[1] + p[2] * exp(-(p[3]-p[1])*x)
 @. model2(x,p) = p[1] + p[2] * x[:, 1] + (p[3] + p[4] * x[:, 1]) * x[:, 2]
 fit_routine(model, xdata, ydata, param=3)
 fit_routine(model, xdata, ydata, param=3, covar=true)</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article><article class="docstring"><header><a class="docstring-binding" id="juobs.bayesian_av" href="#juobs.bayesian_av"><code>juobs.bayesian_av</code></a> — <span class="docstring-category">Function</span></header><section><div><pre><code class="language-julia hljs">bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64}, tmax_array::Array{Int64}, k::Int64, pl::Bool, data::Bool; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
 bayesian_av(fun1::Function, fun2::Function, y::Array{uwreal}, tmin_array::Array{Int64}, tmax_array::Array{Int64}, k1::Int64, k2::Int64, pl::Bool, data::Bool; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
 bayesian_av(fun::Array{Function}, y::Array{uwreal}, tmin_array::Array{Int64}, tmax_array::Array{Int64}, k::Array{Int64}, pl::Bool, data::Bool; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)</code></pre><p>Computes bayesian average of data. For a given fit function, it explores choices of fit intervals, assigning each of them a weight. The function saves the <code>first</code> fit parameter of your function, and then it does the weighted average of it and assigns a systematic. See https://arxiv.org/abs/2008.01069 </p><p>The function takes as input the fit intervals to explore. </p><p><code>tmin_array</code> is an array of integers with the lower bounds on the fit intervals to explore, ***ordered from lower to higher***.</p><p><code>tmax_array</code> is an array of integers with the upper bounds on the fit intervals to explore, ***ordered from lower to higher***.</p><p><code>k</code> is the number of parameters of the fit function to use.</p><p>You can also use as input two fit functions, and two values of <code>k</code>, one for each function. Then, for each fit interval choice, the function explores the two fit functions. This means that for each fit interval choice you get two results: one for the first fit funcction, and another for the second. You can also use a vector of functions and a vector of k (numer of parameters of each funtion) to apply the bayesian averaging method to multiple functions.</p><p>The method returns two objects: first, the weighted average as an uwreal object, with mean value and statistichal error. The second object returned is the systematic error coming from the fit interval variation. If <code>data</code> is <code>true</code>, then returns 4 objects: weighted average, systematic error, a vector with the results of the fit for each fit interval choice, and a vector with the weights associated to each fit.</p><pre><code class="language- hljs">@.fun(x,p) = p[1] * x ^0
 k = 1
 tmin_array = [10,11,12,13,14,15]
 tmax_array = [80,81,82,83,84,85]
 (average, systematics, data, weights) = bayesian_av(fun,x,tmin_array,tmax_array,k,pl=true,data=true)
 @.fun1(x,p) = p[1] * x ^0
 @.fun2(x,p) = p[1] + p[2] * exp( - p[3] * (x))
 k1 = 1
 k2 = 3
 tmin_array = [10,11,12,13,14,15]
 tmax_array = [80,81,82,83,84,85]
 (average, systematics) = bayesian_av(fun1,fun2,x,tmin_array,tmax_array,k1,k2)
 @.fun1(x,p) = p[1] * x ^0
 @.fun2(x,p) = p[1] + p[2] * exp( - p[3] * (x))
 k1 = 1
 k2 = 3
 tmin_array = [10,11,12,13,14,15]
 tmax_array = [80,81,82,83,84,85]
 (average, systematics) = bayesian_av([fun1,fun2],x,tmin_array,tmax_array,[k1,k2])</code></pre></div><a class="docs-sourcelink" target="_blank" href="https://gitlab.ift.uam-csic.es/jugarrio/juobs">source</a></section></article></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="reader.html">« Reader</a><a class="docs-footer-nextpage" href="obs.html">Observables »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 0.27.15 on <span class="colophon-date" title="Friday 1 July 2022 12:33">Friday 1 July 2022</span>. Using Julia version 1.6.5.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>
--- a/docs/make.jl
+++ b/docs/make.jl
 using Documenter, juobs
 push!(LOAD_PATH, "../src")
 makedocs(sitename = "juobs Documentation", doctest=true,
 		  repo = "https://gitlab.ift.uam-csic.es/jugarrio/juobs",
 		  pages = [
 			   "Home" => "index.md",
 			   "Reader" => "reader.md",
 			   "Tools" => "tools.md",
 			   "Observables" =>  "obs.md",
 			   "Linear Algebra" => "linalg.md"],
 		  format = Documenter.HTML(prettyurls=false))
--- a/docs/src/index.md
+++ b/docs/src/index.md
 # DOCUMENTATION
 ## Contents
 ```@contents
 Pages = ["reader.md", "tools.md", "obs.md", "linalg.md"]
 Depth = 3
 ```
--- a/docs/src/linalg.md
+++ b/docs/src/linalg.md
 # Linear Algebra
 ```@docs
 uweigvals
 uweigvecs
 uweigen
 get_matrix
 energies
 getall_eigvals
 getall_eigvecs
 ```
--- a/docs/src/obs.md
+++ b/docs/src/obs.md
 # Observables
 ```@docs
 meff
 mpcac
 dec_const
 dec_const_pcvc
 comp_t0
 ```
--- a/docs/src/reader.md
+++ b/docs/src/reader.md
 # Reader
 ```@docs
 read_mesons
 read_ms
 read_ms1
 read_md
 truncate_data!
 ```
--- a/docs/src/tools.md
+++ b/docs/src/tools.md
 # Tools
 ```@docs
 corr_obs
 corr_sym
 md_sea
 md_val
 lin_fit
 fit_routine
 bayesian_av
 ```
--- a/gevp tools/tools.jl
+++ b/gevp tools/tools.jl
 function read_dat(ens::String, g1::String="G5", g2::String="G5")
    path = joinpath(path_data, ens)
    rep = readdir(path, join=true)
    aux = read_mesons.(rep, g1, g2)
    res = []
    for j =1:length(aux[1])
        aux2 = [aux[i][j] for i =1:length(aux)]
        push!(res, aux2)
    end
    return res
 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_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_lh(mu_list, obs::Vector{Vector{juobs.Corr}}, deg::Bool)
    mul, mus, muh = get_mu(mu_list, deg)
    obs_lh = Array{Array{juobs.Corr}}(undef, length(obs), length(muh))
    for n_obs = 1:length(obs)
        obs_lh[n_obs] = get_lh(mu_list, obs[n_obs], deg)
    end
    return obs_lh
 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_sh(mu_list, obs::Vector{Vector{juobs.Corr}}, deg::Bool)
    mul, mus, muh = get_mu(mu_list, deg)
    obs_sh = Array{Array{juobs.Corr}}(undef, length(obs), length(muh))
    for n_obs = 1:length(obs)
        obs_sh[n_obs] = get_sh(mu_list, obs[n_obs], deg)
    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
 function get_hh(mu_list, obs::Vector{Vector{juobs.Corr}}, deg::Bool)
    mul, mus, muh = get_mu(mu_list, deg)
    obs_hh = Array{Array{juobs.Corr}}(undef, length(obs), length(muh))
    for n_obs = 1:length(obs)
        obs_hh[n_obs] = get_hh(mu_list, obs[n_obs], deg)
    end
    return obs_hh
 end
 mutable struct  infoMatt
    mat::Vector{Matrix}
    mu::Vector{Float64}
    y0::Int64
    function infoMatt(a11::Array{juobs.Corr}, a12::Array{juobs.Corr}, a22::Array{juobs.Corr})
        mu = getfield(a11[1], :mu)
        y0 = Int64(getfield(a11[1], :y0))
        elem11 = (sum(a11[i].obs for i =1:length(a11)))
        elem12 = (sum(a12[i].obs for i =1:length(a12)))
        elem22 = (sum(a22[i].obs for i =1:length(a22)))
        diag = [elem11, elem22]
        subdiag =[ elem12 ]
        #diag = [a11.obs, a22[1].obs.+a22[2].obs.+a22[3].obs]
        #subdiag = [a12[1].obs.+a12[2].obs.+a12[3].obs]
        mat = get_matrix(diag, subdiag)
        return new(mat, mu, y0)
    end
 end
 function comp_mat(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}})
    mu = getfield.(tot11[1],:mu)
    res = Array{infoMatt}(undef, length(mu))
    for i = 1:length(mu)
        a11 = [ tot11[j][i] for j = 1:size(tot11,1) ]
        a12 = [ tot12[j][i] for j = 1:size(tot12,1) ]
        a22 = [ tot22[j][i] for j = 1:size(tot22,1) ]
        res[i] = infoMatt(a11, a12, a22)
        #res[i]=  infoMatt(tot11[i], [tot12[1][i],tot12[2][i],tot12[3][i]], [tot22[1][i],tot22[2][i],tot22[3][i]])
    end
    return res
 end
 function comp_mat(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr}}, tot13::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}}, tot33::Array{Array{juobs.Corr}}, tot23::Array{Array{juobs.Corr}})
    mu = getfield.(tot11[1],:mu)
    res = Array{infoMatt}(undef, length(mu))
    for i = 1:length(mu)
        a11 = [ tot11[j][i] for j = 1:size(tot11,1) ]
        a12 = [ tot12[j][i] for j = 1:size(tot12,1) ]
        a22 = [ tot22[j][i] for j = 1:size(tot22,1) ]
        a33 = [ tot33[j][i] for j = 1:size(tot33,1) ]
        a13 = [ tot13[j][i] for j = 1:size(tot13,1) ]
        a23 = [ tot23[j][i] for j = 1:size(tot23,1) ]
        res[i] = infoMatt(a11, a12, a13, a22, a23, a33)
        #res[i]=  infoMatt(tot11[i], [tot12[1][i],tot12[2][i],tot12[3][i]], [tot22[1][i],tot22[2][i],tot22[3][i]])
    end
    return res
 end
 function comp_energy(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}}, evec::Bool=false; tnot::Int64=3)
    aux_mat = comp_mat(tot11, tot12, tot22)
    y0 = getfield(aux_mat[1], :y0)
    res_en = Array{infoEn}(undef, length(aux_mat))
    if !evec
        for i = 1:length(aux_mat)
            mat_obs = getfield(aux_mat[i], :mat)[y0+1:end-1]
            evalues = uwgevp_tot(mat_obs, tnot)
            res_en[i] = infoEn(energies(evalues), aux_mat[i])
        end
        return res_en
    else
        evec = Array{Array{}}(undef, length(aux_mat))
        for i = 1:length(aux_mat)
            mat_obs = getfield(aux_mat[i], :mat)[y0+1:end-1]
            eigvec = uwgevp_tot(mat_obs, delta_t=3)
            #res_en[i] = infoEn(energies(evalues), aux_mat[i])
            evec[i] = eigvec
        end
        return  evec
    end
 end
 function comp_energy(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr}}, tot13::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}}, tot33::Array{Array{juobs.Corr}}, tot23::Array{Array{juobs.Corr}}, evec::Bool=false; tnot::Int64=2)
    aux_mat = comp_mat(tot11, tot12, tot13, tot22, tot23, tot33)
    y0 = getfield(aux_mat[1], :y0)
    res_en = Array{infoEn}(undef, length(aux_mat))
    if !evec
        for i = 1:length(aux_mat)
            mat_obs = getfield(aux_mat[i], :mat)[y0+1:end-1]
            evalues = uwgevp_tot(mat_obs, tnot, iter=50)
            res_en[i] = infoEn(energies(evalues), aux_mat[i])
        end
        return res_en
    else
        evec = Array{Array{}}(undef, length(aux_mat))
        for i = 1:length(aux_mat)
            mat_obs = getfield(aux_mat[i], :mat)[y0+1:end-1]
            eigvec = uwgevp_tot(mat_obs, delta_t=3)
            println("in comp_energy, i is = ", i)
            #res_en[i] = infoEn(energies(evalues), aux_mat[i])
            evec[i] = eigvec
        end
        return  evec
    end
 end
 mutable struct infoEn
    en_val::Array{Array{uwreal}}
    mu::Vector{Float64}
    y0::Int64
    function infoEn(energ::Array{Array{uwreal}},matr::infoMatt )
        en_val = energ
        mu = getfield(matr,:mu)
        y0 = getfield(matr,:y0)
        new(en_val, mu, y0)
    end
 end
 function pseudo_mat_elem(evec::Array{Array{T,2} where T,1}, mass::uwreal, mu::Array{Float64})
    Rn = [exp(mass * (t)/2) * evec[t][1] for t =1:length(evec)]
    aux = sqrt(2)*sum(mu)  / mass^1.5 
    return uwdot(aux,Rn)
 end
 function vec_mat_elem(evec::Array{Array{T,2} where T,1}, mass::uwreal, deg::Bool)
    if !deg
        Rn = [exp(mass * t/2) * evec[t][1] for t =1:length(evec)]
        return uwdot(Rn , za(beta[iens]) * sqrt(2/ mass))
    else
        Rn = [exp(mass * t/2) * evec[t][4] for t =1:length(evec)]
        println(za(beta[iens]))
        println(mass)
        return uwdot(Rn , za(beta[iens]) * sqrt(2/ mass))
        #return uwdot(Rn , 1/mass)# * sqrt(2/ mass))
 end
 end
 function extract_dec_const(mat_elem::uwreal, mass::uwreal, mu::Array{Float64})
    return sqrt(2)*sum(mu)*mat_elem / mass^1.5
 end
 function match_muc(muh, m_lh, m_lh_star, target)
    M = (m_lh .+ 3 .* m_lh_star) ./ 4
    par, chi2exp = lin_fit(muh, M)
    muh_target = x_lin_fit(par, target)
    return muh_target
 end
 function match_muc(muh, m_lh, m_sh, m_lh_star, m_sh_star, target)
    M = (2 .* m_lh .+ 6 .* m_lh_star .+ m_sh .+ 3 .* m_sh_star) ./ 12
    par, chi2exp = lin_fit(muh, M)
    muh_target = x_lin_fit(par, target)
    return muh_target
 end
\ No newline at end of file
--- a/src/Changelog.md
+++ b/src/Changelog.md
 This file should keep track of all the changes made on a brach since the last merging with master.
 On the master brach, it should keep track of all the major, potentially code-codebreaking, changes perform.
 This file should not be merged into master, but only kept as a reference within its own branch
 ### Antonino Branch
  - Changes before Dic 8 2024:
    - Added juobs_plot.jl with the plot routine: plot_data to plot uwreal and plot_func to plot a function.
    - Updated some documentation in juobs_reader.jl
    - Updated documentation in juobs_tools.jl
    - Updated documentation in juobs_type.jl 
  - Dic 8 2024: added functions plat_av to the existing one. This change should not be code-breaking, since it only add different version of the exising one. Now it is possible to pass a Vector or a Matrix W containing the weight to use in the plateau average. If W is a Matrix, a correlated fit is perform. If W is not given, then the old function is called. 
  - Dic 10 2024: mpcac function now follow the convention more clearly. documentation updated juobs_tools. 
                 pvalue function for uncorrelated data now accept a W vector. default value is Vector{Float64}(), and if W is not provided, it is computed as the inverse of the error squared of the data. 
\ No newline at end of file
--- a/src/juobs.jl
+++ b/src/juobs.jl
 module juobs
 using ADerrors, PyPlot, LaTeXStrings, LinearAlgebra, LsqFit, Optim, ForwardDiff
 using Base: Float64
 import Statistics: mean
 include("juobs_types.jl")
 include("juobs_linalg.jl")
 include("juobs_reader.jl")
 include("juobs_tools.jl")
 include("juobs_plots.jl")
 include("juobs_obs.jl")
 include("../constants/juobs_const.jl")
 include("../constants/path_csv.jl")
 export EnsInfo
 export ens_db, db, db_c, flag_s, sym_bool, ensemble_inv
 export read_mesons, read_mesons_correction, read_mesons_multichunks, read_mesons_correction_multichunks, read_ens_tm, read_ens_tm_sym, read_ens_TSM, read_ens_wil, read_ens_csv, get_corr_tm, get_corr_TSM, get_corr_TSM_multichunks, get_corr_wil, read_ms1, read_ms, read_md, get_YM, get_YM_dYM, truncate_data!, concat_data!
 export get_matrix, energies, uwdot, uweigvals, uweigvecs, uweigen, invert, getall_eigvals, getall_eigvecs, hess_reduce, uwcholesky, transpose, tridiag_reduction, make_positive_def, invert_covar_matrix
 export corr_obs, corr_obs_TSM, corr_sym, corr_sym_E250, corr_sym_D450, md_sea, md_val, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine, global_fit_routine, bayesian_av, pvalue
 export plot_data,plot_func, pvalue, fve, model_av, fit_alg
 export meff, mpcac, dec_const, dec_const_w, dec_const_w, dec_const_pcvc, comp_t0, get_m, get_m_pbc, get_mpcac, get_f_wil, get_f_wil_pbc, get_f_tm, get_f_tm_pbc, get_w0t0_plat
 end # module
--- a/src/juobs_linalg.jl
+++ b/src/juobs_linalg.jl
 #=
 using ForwardDiff, LinearAlgebra
 for op in (:eigvals, :eigvecs)
    @eval function LinearAlgebra.$op(a::Matrix{uwreal})
        function fvec(x::Matrix)
            return LinearAlgebra.$op(x)
        end
        return fvec(value.(a))
        #return uwreal(LinearAlgebra.$op(a.mean), a.prop, ForwardDiff.derivative($op, a.mean)*a.der)
    end
 end
 =#
 @doc raw"""
    get_matrix(diag::Vector{Array}, upper::Vector{Array} )
 This function allows the user to build an array of matrices, where each matrix is a  symmetric matrix of correlators at a given timeslice.   
 It takes as input:
 - `diag::Vector{Array}`:  vector of correlators. Each correlator will constitute a diagonal element of the matrices A[i] where i runs over the timeslices, i.e.
                 `A[i][1,1] = diag[1], .... A[i][n,n] = diag[n]`
                 `Given n=length(diag)`, the matrices will have dimension n*n 
 - `upper::Vector{Array}`:  vector of correlators liying on the upper diagonal position.
                `A[i][1,2] = upper[1], .. , A[i][1,n] = upper[n-1], `
                `A[i][2,3] = upper[n], .. , A[i][n-1,n] = upper[n(n-1)/2]`
                 Given n, `length(upper)=n(n-1)/2`
 The method returns an array of symmetric matrices of dimension n for each timeslice  
 Examples:
 ```@example 
 ## load data 
 pp_data = read_mesons(path, "G5", "G5")
 pa_data = read_mesons(path, "G5", "G0G5")
 aa_data = read_mesons(path, "G0G5", "G0G5")
 ## create Corr struct
 corr_pp = corr_obs.(pp_data)
 corr_pa = corr_obs.(pa_data)
 corr_aa = corr_obs.(aa_data) # array of correlators for different \mu_q combinations
 ## set up matrices
 corr_diag = [corr_pp[1], corr_aa[1]] 
 corr_upper = [corr_pa[1]]
 matrices = [corr_diag, corr_upper]
 Julia> matrices[i]
       pp[i]  ap[i]
       pa[i]  aa[i]
 ## where i runs over the timeslices
 ```
 """
 function get_matrix(corr_diag::Vector{Vector{uwreal}}, corr_upper::Vector{Vector{uwreal}})
    time = length(corr_diag[1]) # total time slices 
    n = length(corr_diag) # matrix dimension
    d = length(corr_upper) # upper elements
    if d != (n*(n-1) / 2)
        error("Error get_matrix")
    end
    res = Vector{Matrix{uwreal}}(undef, time) # array with all the matrices at each time step.
    aux = Matrix{uwreal}(undef, n, n)
    for t in 1:time  
        count=0
        for i in range(n-1,1, step=-1)
            for j in range(n, i+1, step=-1)
                aux[i,j] = corr_upper[d-count][t]
                aux[j,i] = corr_upper[d-count][t]
                count +=1
            end
            aux[i,i] = corr_diag[i][t]
        end
        aux[n,n] = corr_diag[n][t]
        res[t] = copy(aux)
    end
    return res
 end
 get_matrix(corr_diag::Vector{Corr}, corr_upper::Vector{Corr}) = get_matrix(getfield.(corr_diag, :obs), getfield.(corr_upper, :obs))
 @doc raw"""
    energies(evals::Vector{Array}; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
 This method computes the energy level from the eigenvalues according to:
 ``E_i(t) = \log(λ(t) / λ(t+1))``
 where `i=1,..,n` with `n=length(evals[1])` and `t=1,..,T`  total time slices.
 It returns a vector array en where each entry en[i][t] contains the i-th states energy at time t 
 Examples:
 ```@example 
 ## load data
 pp_data = read_mesons(path, "G5", "G5")
 pa_data = read_mesons(path, "G5", "G0G5")
 aa_data = read_mesons(path, "G0G5", "G0G5")
 ## create Corr struct
 corr_pp = corr_obs.(pp_data)
 corr_pa = corr_obs.(pa_data)
 corr_aa = corr_obs.(aa_data) # array of correlators for different \mu_q combinations
 ## set up matrices 
 corr_diag = [corr_pp[1], corr_aa[1]] 
 corr_upper = [corr_pa[1]]
 matrices = [corr_diag, corr_upper]
 ## solve the GEVP
 evals = getall_eigvals(matrices, 5) #where t_0=5
 en = energies(evals)
 Julia> en[i] # i-th state energy at each timeslice
 ```
 """
 function energies(evals::Union{Vector{Vector{uwreal}},Array{Array{uwreal}} }; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)     
    time = length(evals) 
    n = length(evals[1])
    eff_en = Array{Array{uwreal}}(undef, n)
    aux_en = Vector{uwreal}(undef, time-1)
    for i in 1:n
        for t in 1:time-1
            ratio = evals[t][i] / evals[t+1][i]
            aux_en[t] = 0.5*log(ratio * ratio)
        end
        #uwerr.(aux_en)
        #isnothing(wpm) ? uwerr.(aux_en) : [uwerr(a, wpm) for a in aux_en]
        eff_en[i] = copy(aux_en)
    end
    return eff_en
 end
 @doc raw"""
    getall_eigvals(a::Vector{Matrix}, t0; iter=30 )
 This function solves a GEVP problem, returning the eigenvalues, for a list of matrices, taking as generalised matrix the one 
 at index t0, i.e:
 ``C(t_i)v_i = λ_i C(t_0) v_i``,   with   i=1:lenght(a)
 It takes as input:
 - `a::Vector{Matrix}` : a vector of matrices
 - `t0::Int64` : idex value at which the fixed matrix is taken
 - `iter=30` : the number of iterations of the qr algorithm used to extract the eigenvalues 
 It returns:
 - `res` = Vector{Vector{uwreal}}
 where `res[i]` are the generalised eigenvalues of the i-th matrix of the input array. 
 Examples:
 ```@example 
 ## load data
 pp_data = read_mesons(path, "G5", "G5")
 pa_data = read_mesons(path, "G5", "G0G5")
 aa_data = read_mesons(path, "G0G5", "G0G5")
 ## create Corr struct
 corr_pp = corr_obs.(pp_data)
 corr_pa = corr_obs.(pa_data)
 corr_aa = corr_obs.(aa_data) # array of correlators for different \mu_q combinations
 ## set up matrices 
 corr_diag = [corr_pp[1], corr_aa[1]] 
 corr_upper = [corr_pa[1]]
 matrices = [corr_diag, corr_upper]
 ## solve the GEVP
 #evals = getall_eigvals(matrices, 5) #where t_0=5
 Julia>
 ```
 """
 function getall_eigvals(a::Vector{Matrix{uwreal}}, t0::Int64; iter=5 )
    n = length(a)
    res = Vector{Vector{uwreal}}(undef, n)
    [res[i] = uweigvals(a[i], a[t0], iter=iter) for i=1:n]
    return res
 end
 @doc raw"""
    uweigvals(a::Matrix{uwreal}; iter = 30)
    uweigvals(a::Matrix{uwreal}, b::Matrix{uwreal};  iter = 30)
 This function computes the eigenvalues of a matrix of uwreal objects.
 If a second matrix b is given as input, it returns the generalised eigenvalues instead.
 It takes as input:
 - `a::Matrix{uwreal}` : a matrix of uwreal
 - `b::Matrix{uwreal}` : a matrix of uwreal, optional
 - `iter=30`: optional flag to set the iterations of the qr algorithm used to solve the eigenvalue problem
 It returns:
 - `res = Vector{uwreal}`: a vector where each elements is an eigenvalue 
 ```@example
 a = Matrix{uwreal}(nothing, n,n) ## n*n matrix of uwreal with nothing entries
 b = Matrix{uwreal}(nothing, n,n) ## n*n matrix of uwreal with nothing entries
 res = uweigvals(a)  ##eigenvalues
 res1 = uweigvals(a,b) ## generalised eigenvalues
 ```
 """
 function uweigvals(a::Matrix{uwreal}; iter::Int64=5, diag::Bool=true)
    n = size(a,1)
    if n == 2 
        return uweigvals_dim_geq3(a, iter=iter, diag=diag) #uweigvals_dim2(a) 
    else
        return uweigvals_dim_geq3(a, iter=iter, diag=diag)
    end
 end
 function uweigvals(a::Matrix{uwreal}, b::Matrix{uwreal};  iter = 5, diag::Bool=true)
    c = uwdot(invert(b),a)
    n = size(c,1)
    if n == 2
        return uweigvals_dim_geq3(c, iter=iter, diag=diag) #uweigvals_dim2(c) 
    else
        return uweigvals_dim_geq3(c, iter=iter, diag=diag)
    end
    #for k in 1:iter
    #    q_k, r_k = qr(c)
    #    c = uwdot(r_k, q_k)
    #end
 end
 function uweigvals_dim_geq3(a::Matrix{uwreal}; iter = 5, diag::Bool=true)
    n = size(a,1)
    #a = tridiag_reduction(a)
    for k in 1:iter
        q_k, r_k = qr(a)
        a = uwdot(r_k, q_k)
    end
    diag == true ? (return get_diag(a)) : (return a)
 end
 function uweigvals_dim2(a::Matrix{uwreal})
    res = Vector{uwreal}(undef,2)
    aux_sum = a[1,1] + a[2,2]
    delta = ADerrors.sqrt((aux_sum)^2 - 4*(a[1,1]*a[2,2]- a[1,2]*a[2,1]))
    res[1] = (aux_sum + delta) / 2
    res[2] = (aux_sum - delta) / 2
    return res
 end
 @doc raw"""
    getall_eigvecs(a::Vector{Matrix}, delta_t; iter=30 )
 This function solves a GEVP problem, returning the eigenvectors, for a list of matrices.
 ``C(t_i)v_i = λ_i C(t_i-\delta_t) v_i  ``, with i=1:lenght(a)
 Here `delta_t` is the time shift within the two matrices of the problem, and is kept fixed.
 It takes as input:
 - `a::Vector{Matrix}` : a vector of matrices
 - `delta_t::Int64` : the fixed time shift t-t_0
 - `iter=30` : the number of iterations of the qr algorithm used to extract the eigenvalues 
 It returns:
 - `res = Vector{Matrix{uwreal}}`
 where each `res[i]` is a matrix  with the eigenvectors as  columns
 Examples:
 ```@example
 mat_array = get_matrix(diag, upper_diag)
 evecs = getall_eigvecs(mat_array, 5)
 ```
 """
 function getall_eigvecs(a::Vector{Matrix{uwreal}}, delta_t; iter=5)
    n = length(a)-delta_t
    res = Vector{Matrix{uwreal}}(undef, n)
    [res[i] = uweigvecs(a[i+delta_t], a[i]) for i=1:n]
    return res 
 end
 @doc raw"""
    uweigvecs(a::Matrix{uwreal}; iter = 30)
    uweigvecs(a::Matrix{uwreal}, b::Matrix{uwreal};  iter = 30)
 This function computes the eigenvectors of a matrix of uwreal objects.
 If a second matrix b is given as input, it returns the generalised eigenvectors instead.
 It takes as input:
 - `a::Matrix{uwreal}` : a matrix of uwreal
 - `b::Matrix{uwreal}` : a matrix of uwreal, optional
 - `iter=30` : the number of iterations of the qr algorithm used to extract the eigenvalues 
 It returns:
 - `res = Matrix{uwreal}`: a matrix where each column is an eigenvector 
 Examples:
 ```@example
 a = Matrix{uwreal}(nothing, n,n) ## n*n matrix of uwreal with nothing entries
 b = Matrix{uwreal}(nothing, n,n) ## n*n matrix of uwreal with nothing entries
 res = uweigvecs(a)  ##eigenvectors in column 
 res1 = uweigvecs(a,b) ## generalised eigenvectors in column 
 ```
 """
 function uweigvecs(a::Matrix{uwreal}; iter = 10)
    n = size(a,1)
    # if n > 2
        # a, q = tridiag_reduction(a, q_mat=true)
    # end
    u = idty(n)
    for k in 1:iter
        q_k, r_k = qr(a)
        a = uwdot(r_k, q_k)
        u = uwdot(u, q_k)
    end
    #n > 2 ? (return  uwdot(q,u)) : (return u)
    return u
 end
 function uweigvecs(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 10)
    c  = uwdot(invert(b), a)
    #return uweigvecs(c, iter=iter)
    n = size(c,1)
    u = idty(n)
    for k in 1:iter
        q_k, r_k = qr(c)
        c = uwdot(r_k, q_k)
        u = uwdot(u, q_k)
    end
    return  u
 end
 @doc raw"""
    uweigen(a::Matrix{uwreal}; iter = 30)
    uweigen(a::Matrix{uwreal}, b::Matrix{uwreal};  iter = 30)
 This function computes the eigenvalues and the eigenvectors  of a matrix of uwreal objects.
 If a second matrix b is given as input, it returns the generalised eigenvalues and eigenvectors instead.
 It takes as input:
 - `a::Matrix{uwreal}` : a matrix of uwreal
 - `b::Matrix{uwreal}` : a matrix of uwreal, optional
 - `iter=30` : the number of iterations of the qr algorithm used to extract the eigenvalues 
 It returns:
 - `evals = Vector{uwreal}`: a vector where each elements is an eigenvalue 
 - `evecs  = Matrix{uwreal}`: a matrix where the i-th column is the eigenvector of the i-th eigenvalue
 Examples:
 ```@example
 a = Matrix{uwreal}(nothing, n,n) ## n*n matrix of uwreal with nothing entries
 b = Matrix{uwreal}(nothing, n,n) ## n*n matrix of uwreal with nothing entries
 eval, evec = uweigen(a)   
 eval1, evec1 = uweigvecs(a,b) 
 ```
 """ 
 function uweigen(a::Matrix{uwreal}; iter = 5, diag::Bool=true)
    return uweigvals(a, iter=iter,diag=diag), uweigvecs(a, iter=iter)
 end
 function uweigen(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 5, diag::Bool=true)
    return uweigvals(a, b, iter=iter, diag=diag), uweigvecs(a, b, iter=iter)
 end
 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
 function get_diag(a::Matrix{uwreal})
    n = size(a,1)
    res = [a[k, k] for k = 1:n]
    return res
 end
 function get_diag(v::Vector{uwreal})
    n = length(v)
    res = idty(n)
    for i =1:n
        res[i,i] = v[i]
    end
    return res
 end
 """
 Given an input matrix of uwreals, this function returns the inverse of that matrix using QR decomposition and
 backward substitution.
 """
 function invert(a::Matrix{uwreal})
    q,r = qr(a)
    r_inv = upper_inversion(r)
    return uwdot(r_inv, transpose(q))
 end
 """
 Performs a Cholesky decomposition of A, which must 
    be a symmetric and positive definite matrix. The function
    returns the lower variant triangular matrix, L.
 """
 function uwcholesky(A::Matrix{uwreal})
    n = size(A,1)
    L = fill(uwreal(0), n, n)
    for i in 1:n 
        for k in 1:i
            tmp_sum = sum(L[i,j] * L[k,j] for j in 1:k)
            if i==k 
                L[i,k] = sqrt(((A[i,i] - tmp_sum)))
            else
                L[i,k] = (1.0 / L[k,k] * (A[i,k]-tmp_sum))
            end
        end
    end
    return L 
 end
 """
 qr(A::Matrix{uwreal})
 qr(A::Matrix{Float64})
 This methods returns the QR decomposition of a matrix A of either uwreal or Float64 variables
 """
 function qr(A::Matrix{uwreal})
    m,n = size(A)
    Q = idty(m)
    for i in 1:(n-(m==n))
        H = idty(m)
        H[i:m,i:m] = make_householder(A[i:m,i])
        Q = uwdot(Q,H)
        A = uwdot(H,A)
    end
    for i in 2:n
        for j in 1:i-1
            A[i,j] = 0.0
        end
    end
    return Q,A
 end
 function qr(A::Matrix{Float64})
    m,n = size(A)
    Q = Matrix{Float64}(I,m,m)
    for i in 1:(n-(m==n))
        H =  Matrix{Float64}(I,m,m)
        H[i:m,i:m] = make_householder(A[i:m,i])
        X = similar(Q)
        Y = similar(A)
        Q = mul!(X,Q,H)
        A = mul!(Y,H,A)
    end
    return Q,A
 end
 """
 make_householder(a::Vector{uwreal})
 make_householder(a::Vector{Float64})
 This method returns the householder matrix used in qr decomposition to get an upper triangular matrix
 It takes as input either a vector of uwreal or a vector of Float64
 """
 function make_householder(a::Vector{uwreal})
    n = length(a)
    v = Vector{uwreal}(undef, n)
    for i in 1:n
        v[i] = a[i]  / (a[1] + uwcopysign(uwnorm(a), a[1]))
    end
    v[1] = uwreal(1.0) 
    H = idty(n)
    H = H -  (2.0 / uwdot(v,v)) * uwdot(reshape(v, :,1), transpose(v)) #this last term is a vector product  
    return H
 end
 function make_householder(a::Vector{Float64})
    n = length(a)
    v = a / (a[1] + copysign(norm(a), a[1])) 
    v[1] = 1.0
    H = Matrix{Float64}(I, n, n)
    println(typeof(v))
    y = similar(H)
    println(typeof(y))
    H = H .-  (2.0 / dot(v,v)) * mul!(y, v, transpose(v))
    return H
 end
 """
 This function computes a vector u that generates a Householder reflection H = I - uu* satisfying Ha=nu*e1.
 Accumulating this transformations any matrix can
 be reduced to upper Hessenberg form.
 """
 function house_gen(a::Vector{uwreal})
    u = deepcopy(a)
    nu = uwnorm(a)
    if nu == 0 
        u[1] = sqrt(2)
    end
    if u[1] != 0 
        rho = u[1] / sqrt(u[1]*u[1])
    else
        rho = 1
    end
    u = uwdot(rho/nu, u) 
    u[1] += 1
    u =  uwdot(u , 1/sqrt(u[1]))
    nu = -(rho) * nu
    return u, nu 
 end
 """
 hess_red(a::Matrix{uwreal})
 Given a matrix of uwreal variables, this function returns the Hessenberg form of the matrix by applying householder transformations.
 """
 function hess_reduce(a::Matrix{uwreal}; q_mat::Bool=false)
    n = size(a,1)
    H = deepcopy(a)
    Q = idty(n)
    for k in 1:n-2
        u, H[k+1,k] = house_gen(H[k+1:n,k])
        Q[k+1:n,k] = u
        v = uwdot(u,H[k+1:n, k+1:n])
        H[k+1:n, k+1:n] .= H[k+1:n, k+1:n] .- uwdot(reshape(u,length(u),1),v)
        H[k+2:n,k] .= 0.0 
        v = uwdot(H[1:n,k+1:n], u)
        H[1:n,k+1:n] .= H[1:n,k+1:n] .- uwdot(v,reshape(u,1,length(u)))
    end
    if q_mat 
        Id = idty(n)
        for k  in  reverse(1:n-2)
            u = Q[k+1:n, k]
            v = uwdot(u, Q[k+1:n,k+1:n])
            Q[k+1:n,k+1:n] .= Q[k+1:n,k+1:n] .- uwdot(reshape(u,length(u),1),v)
            Q[:,k] = Id[:,k]
        end
        return H,Q
    else
        return H
    end
 end
 """
 """
 function tridiag_reduction(a::Matrix{uwreal}; q_mat::Bool=false)
    n = size(a,1)
    H = deepcopy(a)
    Q = idty(n)
    for k in 1:n-2
        u, H[k+1,k] = house_gen(H[k+1:n,k])
        Q[k+1:n,k] = u
        v = uwdot(u,H[k+1:n, k+1:n])
        H[k+1:n, k+1:n] .= H[k+1:n, k+1:n] .- uwdot(reshape(u,length(u),1),v)
        H[k+2:n,k] .= 0.0 
        v = uwdot(H[1:n,k+1:n], u)
        H[1:n,k+1:n] .= H[1:n,k+1:n] .- uwdot(v,reshape(u,1,length(u)))
        H[k,k+2:n] .= 0.0 
    end
    if q_mat 
        Id = idty(n)
        for k  in  reverse(1:n-2)
            u = Q[k+1:n, k]
            v = uwdot(u, Q[k+1:n,k+1:n])
            Q[k+1:n,k+1:n] .= Q[k+1:n,k+1:n] .- uwdot(reshape(u,length(u),1),v)
            Q[:,k] = Id[:,k]
        end
        return H,Q
    else
        return H
    end
 end
 """
 upper_inversion(u::Matrix{uwreal})
 This method uses backward propagation to compute the inverse of an upper triangular matrix u. Note that the inverse of u must be upper triangular as well.
 Tests executed so far always confirmed this property.
 """
 function upper_inversion(u::Matrix{uwreal})
    n = size(u,1)
    u_inv = Matrix{uwreal}(undef, n, n)
    for i in 1:n
        e_i = canonical_basis(n, i)
        x = backward_sub(u, e_i)
        u_inv[:,i] = x
    end
    return u_inv
 end
 """
 backward_sub(u::Matrix{uwreal}, y::Vector{uwreal})
 This method operates backward substitution to solve a liner system Ux = y where U is an upper triangular uwreal matrix 
 and y is a uwreal vector. It returns a vector of uwreal x that might be used to inverte the upper triangular matrix U.
 """
 function backward_sub(u::Matrix{uwreal}, y::Vector{uwreal})
    n = length(y)
    x =  Vector{uwreal}(undef, n)
    x[n] = y[n] / u[n,n]
    for i in range(n-1,1, step=-1)
        temp = 0.0
        for j in i+1:n
            temp += u[i,j] * x[j]
        end
        x[i] = (y[i] - temp)  /u[i,i]
        if x[i] !=0
        end
    end
    return x
 end
 """
 canonical_basis(n::Int64, k::Int64)
 Given the dimension n and the position k, canonical_basis returns a vector of uwreal with all entries set to zero
 except for the k-th entry which is set to 1.0.
 This method is used as input vector in backward_sub to compute the inverse of an upper triangular matrix.
 """
 function canonical_basis(n::Int64, k::Int64)
    e = zeros(n)
    e[k] = 1.0
    return [uwreal(e[i]) for i=1:n]
 end
 """
 uwdot(a::Vector{uwreal}, b::Vector{uwreal})
 uwdot(a::Matrix{uwreal}, b::Matrix{uwreal})
 Depending on the arguments, uwdot returns the vector vector or the matrix matrix product 
 of uwreal data type.
 """
 uwdot(a::Vector{uwreal}, b::Vector{uwreal}) = sum(a .* b)
 uwdot(a::Vector{uwreal}, b::uwreal) = [a[i] * b for i=1:length(a)]
 uwdot(a::uwreal, b::Vector{uwreal}) = uwdot(b,a) 
 uwdot(a::Matrix{uwreal}, b::Vector{uwreal}) = uwdot(a, reshape(b,(length(b),1)))
 uwdot(a::Vector{uwreal}, b::Matrix{uwreal}) = uwdot(reshape(a,(1,length(a))), b)
 function uwdot(a::Matrix{uwreal}, b::Matrix{uwreal})
    n,m = size(a)
    k,p = size(b)
    c = Matrix{uwreal}(undef, n, p)
    if m != k
        error("Error uwdot: in a*b the size of a is ", size(a), " while the size of b is ", size(b) )
    end
    for i = 1:n
        for j = 1:p
            c[i, j] = sum(a[i, :] .* b[:, j])
        end
    end
    return c
 end
 """
 uwnorm(a::Vector{uwreal})
 It returns the norm of a vector of uwreal
 """
 function uwnorm(a::Vector{uwreal}; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
    norm = sqrt(uwdot(a,a))
    #isnothing(wpm) ? uwerr(norm) : uwerr(norm, wpm)
    return norm 
 end
 """
 Return the transpose of a matrix or vector of uwreals
 """
 function Base.transpose(a::Matrix{uwreal})
    res = similar(a)
    n = size(a, 1)
    if size(a,1) == 1
        return reshape(a, :,1)
    elseif size(a,2) == 1
        return reshape(a, 1,:)    
    end
    for i in 1:n-1
        for j in i+1:n
            temp = a[i,j]
            res[i,j] = a[j,i]
            res[j,i] = temp
        end
        res[i,i] = a[i,i]
    end
    res[n,n] = a[n,n]
    return res
 end
 function Base.transpose(vec::Vector{uwreal}) 
    res = deepcopy(vec)
    return reshape(res,1,:)
 end
 """
 uwcopysign(a::uwreal, b::uwreal)
 It implements  the copysign function for uwreal variables
 """
 function uwcopysign(a::uwreal, b::uwreal)
    c = deepcopy(a)
    aux =  copysign(value(a), value(b))
    c.mean = aux
    return c
 end
 """
 idty(n::Int64)
 It returns an indetity matrix of uwreal variables given the size n 
 """
 function idty(n::Int64)
    res = Matrix{uwreal}(undef,n,n)
    for i in 1:n
        for j in 1:n
            if i==j
                res[i,j]=1.0
            else
                res[i,j]=0.0
            end
        end
    end
    return res
 end
 function make_positive_def(A::Matrix; eps::Float64=10^(-14))
    vals, vecs = eigen(Symmetric(A))
    idx = findall(x-> x<0.0, vals)
    vals[idx] .= eps
    return  Symmetric(vecs * Diagonal(vals) * vecs')
    #B = 0.5 * (A + A')
    #U, S, V = svd(B)
    #H = V * Diagonal(S) * V'
    #A_hat = 0.5 * (B + H)
    #A_hat = 0.5 * ( A_hat + A_hat')
    #k = 0
    #while !isposdef(A_hat)
    #    mineig = minimum(eigvals(A_hat))
    #    eps =  2.220446049250313e-16
    #    A_hat = A_hat + (-mineig*k^2 .+ eps*Matrix(I, size(A)))
    #    k+=1
    #end
    #return A_hat
 end
 function invert_covar_matrix(A::Matrix; eps::Float64=10^(-9))
    F = svd(A)
    s_diag = F.S
    for i in 1:length(s_diag)
        if s_diag[i] < eps
            s_diag[i] = 0.0
        else
            s_diag[i] = 1 / s_diag[i]
        end
    end
    cov_inv = F.V * Diagonal(s_diag) * F.U'
    return cov_inv  
 end
 function more_penrose_inv(A::Matrix)
    F = svd(A)
 end
 ########################
 # OPERATOR OVERLOADING
 ########################
 function Base.:*(x::uwreal, y::Array{Any})
    N = size(y, 1)
    return fill(x, N) .* y
 end  
 Base.:*(x::Array{Any}, y::uwreal) = Base.:*(y,x)
 function Base.:*(x::uwreal, y::Array{Float64})
    N = size(y, 1)
    return fill(x, N) .* y
 end
 Base.:*(x::Array{Float64}, y::uwreal) = Base.:*(y,x)
 function Base.:*(x::uwreal, y::Array{uwreal})
    N = size(y, 1)
    return fill(x, N) .* y
 end
 Base.:*(x::Array{uwreal}, y::uwreal) = Base.:*(y,x)
 function Base.:+(x::uwreal, y::Vector{uwreal})
    N = size(y, 1)
    return fill(x, N) .+ y
 end
 function Base.:+(x::Float64, y::Vector{Float64})
    N = size(y, 1)
    return fill(x, N) .+ y
 end
 function Base.:+(x::Float64, y::Vector{uwreal})
    N = size(y, 1)
    return fill(x, N) .+ y
 end
 function Base.:-(x::Float64, y::Vector{Any})
    N = size(y, 1)
    return fill(x, N) .- y
 end
 function Base.:-(x::Float64, y::Vector{Float64})
    N = size(y, 1)
    return fill(x, N) .- y
 end
 function Base.:-(x::Float64, y::Vector{uwreal})
    N = size(y, 1)
    return fill(x, N) .- y
 end
 function Base.length(uw::uwreal)
    return length(value(uw))
 end
 function Base.iterate(uw::uwreal)
    return iterate(uw, 1)
 end
 function Base.iterate(uw::uwreal, state::Int64)
    if state > length(uw)
        return nothing
    else
        return uw[state], state +1
    end
 end
 function Base.iterate(uw::Vector{uwreal})
    return iterate(uw, 1)
 end
 function Base.iterate(uw::Vector{uwreal}, state::Int64)
    if state > length(uw)
        return nothing
    else
        return uw[state], state +1
    end
 end
 function Base.getindex(uw::uwreal, ii::Int64)
    idx = getindex(value(uw), ii)
    return  uw # uwreal([getindex(value(uw), kwargs...), err(uw)], " ")
 end
 function Base.abs2(uw::uwreal)
    return (uw^2)^0.5
 end
 """
 Base.:*(x::uwreal, y::Matrix{uwreal})
 Multiplication operator overloading between uwreal variable and matrix of uwreal 
 """
 function Base.:*(x::uwreal, y::Matrix{uwreal})
    n,m = size(y)
    res = Matrix{uwreal}(undef, n,m)
    for i in 1:n
        for j in 1:m
            res[i,j] = x * y[i,j]
            if res[i,j] !=0 && res[i,j] !=1
            end
        end        
    end 
    return res   
 end
--- a/src/juobs_obs.jl
+++ b/src/juobs_obs.jl
--- a/src/juobs_plots.jl
+++ b/src/juobs_plots.jl
 @doc raw"""
    function plot_data(data::Vector{uwreal},xdata::Union{Vector{uwreal},Vector{T} where T<:real};
      wpm::Union{Dict{String, Vector{Float64}},Nothing} = nothing,
      c::String = "k",
      fmt::String = "x",
      label::AbstractString = "",
      title::AbstractString = "",
      xlabel::AbstractString = "x",
      ylabel::AbstractString = "y",
      figs::Union{Tuple{Figure, PyPlot.PyCall.PyObject},Nothing}=nothing)
 Generate a plot with `data` as function of `xdata`. If `xdata` is a `Vector{uwreal}`
 it plots the error bar also in the x-axis. 
 It return a tuple `fig,ax` to easily allow any modification to the plot. 
 # Arguments  
 - `wpm::Union{Dict{String,Vector{Float64}},Dict{Int64,Vector{Float64}},Nothing}=nothing`
  windows parameter used by `ADerrors` to estimate the integration window
 - `label::AbstractString = ""`: set the label of `data`. When set, it prints the legend
 - `figs::Union{Tuple{Figure, PyPlot.PyCall.PyObject},Nothing}=nothing`: 
  A tuple with an exiting `fig` and `ax`. If passed, the fuction prints the data on
  the existing figure
 # Examples
 ```@example
  using ADerror, juobs
  xdata = collect(1:100)
  ydata = [uwreal([sin(x*0.1*π),0.01],"sin") for x in xdata]
  fig,ax = data_plot(ydata,xdata,label=L"$\sin(\frac{x*\pi}{10}$")
  ydata2 = [uwreal([cos(x*0.1*π),0.01],"cos") for x in xdata]
  xdata2 = [uwreal([x,0.1], "xdata") for x in xdata];
  fig,ax = data_plot(ydata2,xdata2,label=L"$\cos(\frac{x*\pi}{10}$",figs = (fig,ax))
 ```
 See also: [`plot_func`](@ref)
 """
 function plot_data(data::Vector{uwreal},xdata::Vector{T} where T <: Real; 
  wpm::Union{Dict{String, Vector{Float64}},Nothing} = nothing,
  c = "k",
  fmt = "x",
  label::AbstractString = "",
  title::AbstractString = "",
  xlabel::AbstractString = "x",
  ylabel::AbstractString = "y",
  figs::Union{Tuple{Figure, PyPlot.PyCall.PyObject},Nothing}=nothing)
  isnothing(wpm) ? uwerr.(data) : [uwerr(d,wpm) for d in data];
  fig, ax = isnothing(figs) ?  subplots(1) : figs;
  line, = ax.errorbar(xdata,value.(data), err.(data), c=c,fmt=fmt)
  if length(label)!=0
    line.set_label(label)
    ax.legend();
  end
  fig.suptitle(title);
  ax.set_ylabel(ylabel);
  ax.set_xlabel(xlabel);
  return fig,ax;
 end
 function plot_data(data::Vector{uwreal},xdata::Vector{uwreal}; 
  wpm::Union{Dict{String, Vector{Float64}},Nothing} = nothing,
  c = "k",
  fmt = "x",
  label::AbstractString = "",
  title::AbstractString = "",
  xlabel::AbstractString = "x",
  ylabel::AbstractString = "y",
  figs::Union{Tuple{Figure, PyPlot.PyCall.PyObject},Nothing}=nothing)
  isnothing(wpm) ? uwerr.(data)  : [uwerr(d,wpm) for d in data];
  isnothing(wpm) ? uwerr.(xdata) : [uwerr(x,wpm) for x in xdata];
  fig, ax = isnothing(figs) ?  subplots(1) : figs;
  line, = ax.errorbar(value.(xdata),value.(data), ADerrors.err.(data), ADerrors.err.(xdata),c=c,fmt=fmt)
  if length(label)!=0
    line.set_label(label)
    ax.legend();
  end
  fig.suptitle(title);
  ax.set_ylabel(ylabel);
  ax.set_xlabel(xlabel);
  return fig,ax;
 end
 @doc raw"""
    function plot_func(f::Function,xdata::Vector{T} where T<:Real;
      wpm::Union{Dict{String,Vector{Float64}},Nothing} = nothing,
      c::String ="green",
      alpha::Float64=0.25,
      label::AbstractString = "",
      figs::Union{Tuple{Figure, PyPlot.PyCall.PyObject},Nothing}=nothing,
 It plots the function `f` evaluated at `xdata`through broadcasting. 
 The function assumes that `f` return type is `uwreal` and is called as `f(x)`.
 It returns a tuple with `fig,ax` to easiy edit the plot
 # Arguments  
 - `wpm::Union{Dict{String,Vector{Float64}},Dict{Int64,Vector{Float64}},Nothing}=nothing`
  windows parameter used by `ADerrors` to estimate the integration window
 - `label::AbstractString = ""`: set the label of `data`. When set, it prints the legend
 - `figs::Union{Tuple{Figure, PyPlot.PyCall.PyObject},Nothing}=nothing`: 
  A tuple with an exiting `fig` and `ax`. If passed, the fuction prints the data on
  the existing figure
 # Example
 ```@example
  using ADerror, juobs
  xdata = collect(1:100)
  f(x) = uwreal([sin(x*0.1*π),0.1],"sin")
  fig,ax = plot_func(f,xdata,label=L"\sin(\frac{x*\pi}{10})")
 ```
 See also: [`plot_data`](@ref)
 """
 function plot_func(f::Function, xdata::Vector{T} where T<:Real; 
  wpm::Union{Dict{String,Vector{Float64}},Nothing} = nothing,
  c::String ="green",
  alpha::Float64=0.25,
  figs::Union{Tuple{Figure, PyPlot.PyCall.PyObject},Nothing}=nothing,
  label::AbstractString = "",
  kwargs...) 
  fig,ax = isnothing(figs) ? subplots(1) : figs
  y = f.(xdata);
  isnothing(wpm) ? uwerr.(y) : [uwerr(a,wpm) for a in y];
  v,e = value.(y), err.(y);
  l = ax.fill_between(xdata,v.-e,v.+e, color=c,alpha=alpha);
  if label!=""
    l.set_label(label);
    ax.legend();
  end
  v,i = findmax(x->value(x), y);
  e = err(y[i])
  ax.set_ylim(top=v+20*e)
  return fig,ax;
 end
\ No newline at end of file
--- a/src/juobs_reader.jl
+++ b/src/juobs_reader.jl
--- a/src/juobs_tools.jl
+++ b/src/juobs_tools.jl
--- a/src/juobs_types.jl
+++ b/src/juobs_types.jl
 ##############UTILITY NAMES##############
  #=  STARTING AT 0
    noise_name=['Z2', 'GAUSS', 'U1', 'Z4']
    gamma_name = ['G0', 'G1', 'G2', 'G3',
 'invalid', 'G5', '1', 'G0G1', 'G0G2', 'G0G3',
 'G0G5', 'G1G2', 'G1G3', 'G1G5', 'G2G3', 'G2G5', 'G3G5']
    # quark and gauge Smearing types
    qs_name=['Local', 'Wuppertal', '3D Gradient Flow', 'Gradient Flow']
    gs_name=['Local', 'APE', '3D Wilson Flow', 'Quark 3D Gradient Flow', 'Quark Gradient Flow']
 =#
 const noise_name=["Z2", "GAUSS", "U1", "Z4"]
 const gamma_name = ["G0", "G1", "G2", "G3",
 "invalid", "G5", "1", "G0G1", "G0G2", "G0G3",
 "G0G5", "G1G2", "G1G3", "G1G5", "G2G3", "G2G5", "G3G5", 
 "G0_d1", "G1_d1", "G2_d1", "G3_d1",
 "invalid", "G5_d1", "1_d1", "G0G1_d1", "G0G2_d1", "G0G3_d1",
 "G0G5_d1", "G1G2_d1", "G1G3_d1", "G1G5_d1", "G2G3_d1", "G2G5_d1", "G3G5_d1",
 "G0_d2", "G1_d2", "G2_d2", "G3_d2",
 "invalid", "G5_d2", "1_d2", "G0G1_d2", "G0G2_d2", "G0G3_d2",
 "G0G5_d2", "G1G2_d2", "G1G3_d2", "G1G5_d2", "G2G3_d2", "G2G5_d2", "G3G5_d2"]
 const qs_name=["Local", "Wuppertal", "3D Gradient Flow", "Gradient Flow"]
 const gs_name=["Local", "APE", "3D Wilson Flow", "Quark 3D Gradient Flow", "Quark Gradient Flow"]
 mutable struct GHeader
    ncorr::Int32
    nnoise::Int32
    tvals::Int32
    noisetype::Int32
    hsize::Int32 #headersize
    GHeader(a, b, c, d) = new(a, b, c, d, 4*4)
    function GHeader(x::Vector{Int32})
        a = new(x[1], x[2], x[3], x[4], 4*4)
        return a
    end
 end
 mutable struct Sm
    type::Int32
    niter::Int32
    eps::Float64
    qg::Int32 #1->fermionic 2->gluonic
    Sm(t,q) = new(t, 0, 0.0, q)
    Sm(t, n, e, q) = new(t, n, e, q)
 end
 mutable struct CHeader
    k1::Float64
    k2::Float64
    mu1::Float64
    mu2::Float64
    dp1::Float64
    dp2::Float64
    theta1::Vector{Float64}
    theta2::Vector{Float64}
    q1::Sm
    q2::Sm
    g1::Sm
    g2::Sm
    type1::Int32
    type2::Int32
    x0::Int32
    is_real::Int32
    hsize::Int32 #headersize
    dsize::Int32 #datasize / (nnoise * T * ncfg)
    function CHeader(aux_f::Vector{Float64}, aux_i::Vector{Int32}, theta::Vector{Float64}, sm_par::Vector{Sm})
        a = new()
        a.k1 = aux_f[1]
        a.k2 = aux_f[2]
        a.mu1 = aux_f[3]
        a.mu2 = aux_f[4]
        a.dp1 = aux_f[5]
        a.dp2 = aux_f[6]
        a.type1 = aux_i[1]
        a.type2 = aux_i[2]
        a.x0 = aux_i[3]
        a.is_real = aux_i[4]
        a.theta1 = theta[1:3]
        a.theta2 = theta[4:6]
        a.q1 = sm_par[1]
        a.q2 = sm_par[2]
        a.g1 = sm_par[3]
        a.g2 = sm_par[4]
        a.hsize = 8*12 + 4*8 #without smearing parameters niter, neps
        if a.q1.type != 0
            a.hsize += 8+4
        end
        if a.q2.type != 0
            a.hsize += 8+4
        end
        if a.g1.type != 0 && a.g1.type != 3 && a.g1.type != 4
            a.hsize += 8+4
        end
        if a.g2.type != 0 && a.g2.type != 3 && a.g2.type != 4
            a.hsize += 8+4
        end
        a.dsize = 16 - 8* a.is_real 
        return a
    end
    function CHeader(aux_f::Vector{Float64}, aux_i::Vector{Int32})
        a = new()
        a.k1 = aux_f[1]
        a.k2 = aux_f[2]
        a.mu1 = aux_f[3]
        a.mu2 = aux_f[4]
        a.dp1 = 0.0
        a.dp2 = 0.0
        a.type1 = aux_i[1]
        a.type2 = aux_i[2]
        a.x0 = aux_i[3]
        a.is_real = aux_i[4]
        a.theta1 = zeros(3)
        a.theta2 = zeros(3)
        a.q1 = Sm(0, 1)
        a.q2 = Sm(0, 1)
        a.g1 = Sm(0, 2)
        a.g2 = Sm(0, 2)
        a.hsize = 8*4 + 4*4
        a.dsize = 16 - 8* a.is_real 
        return a
    end
 end
 function Base.:(==)(a::CHeader, b::CHeader)
    for s in [:k1, :k2, :mu1, :mu2, :dp1, :dp2, :type1, :type2, :x0, :is_real, :theta1, :theta2]
        if getfield(a, s) != getfield(b, s)
            return false
        end
    end
    return true
 end
 #Base.:(!=)(a::CHeader, b::CHeader) = !(a == b)
 mutable struct CData
    header::CHeader
    vcfg::Array{Int32}
    re_data::Array{Float64}
    im_data::Array{Float64}
    id::String
    CData(a, b, c, d, e) = new(a, b, c, d, e)
 end
 Base.copy(a::CData) = CData(a.header, a.vcfg, a.re_data, a.im_data, a.id)
 mutable struct Corr
    obs::Vector{uwreal}
    kappa::Vector{Float64}
    mu::Vector{Float64}
    gamma::Vector{String}
    y0::Int64
    theta1::Vector{Float64}
    theta2::Vector{Float64}
    function Corr(a::Vector{uwreal}, b::CData)
        h = getfield(b, :header)
        kappa = [h.k1, h.k2]
        mu = [h.mu1, h.mu2]
        gamma = [gamma_name[h.type1+1], gamma_name[h.type2+1]]
        y0 = Int64(h.x0)
        theta1 = h.theta1
        theta2 = h.theta2
        return new(a, kappa, mu, gamma, y0, theta1, theta2)
    end 
    function Corr(a::Vector{uwreal}, b::Vector{CData})
        sym = [:k1, :k2, :mu1, :mu2, :type1, :type2, :x0]
        h = getfield.(b, :header)
        for s in sym
            if !all(getfield.(h, s) .== getfield(h[1], s))
                error("Corr: Parameter mismatch")
            end
        end
        kappa = [h[1].k1, h[1].k2]
        mu = [h[1].mu1, h[1].mu2]
        gamma = [gamma_name[h[1].type1+1], gamma_name[h[1].type2+1]]
        y0 = Int64(h[1].x0)
        theta1 = h[1].theta1
        theta2 = h[1].theta2
        return new(a, kappa, mu, gamma, y0, theta1, theta2)
    end 
    Corr(o::Vector{uwreal}, k::Vector{Float64}, m::Vector{Float64}, g::Vector{String}, y::Int64, th1::Vector{Float64}, th2::Vector{Float64} ) = new(o, k, m, g, y, th1, th2)
 end
 mutable struct YData
    vtr::Vector{Int32}
    t::Vector{Float64}
    obs::Array{Float64, 3}
    id::String
    YData(a, b, c, d) = new(a, b, c, d)
 end
 @doc raw"""
    EnsInfo(ens_id::String, info::Vector{Any})
 Stores information about the ensamble `ens_id`
 `info::Vector{Any}` is a vector with 6 components that are:
  -`L::Int64` = number of site in the spatial direction
  -`T::Int64` = number of site in the temporal direction
  -`beta::Float64`
  -`ca::Float64` = c_A improvement parameters of the axial current
  -`dtr::Int64` = ??
  -`vrw::Union{String,Vector{string}}` = version of the reweighting factor. 
  -`cnfg::Vector{Int64}`= number of configuration per replica
 The fields in `EnsInfo` are the same as before, plus the `id::String` of the ensemble
 Examples:
 ```@example    
 id = "H101"
 ens = EnsInfo(id, ens_db[id])
 ```
 """
 mutable struct EnsInfo
    id::String
    L::Int64
    T::Int64
    beta::Float64
    ca::Float64
    dtr::Int64
    vrw::Union{String, Vector{String}}
    cnfg::Vector{Int64}
    function EnsInfo(ens_id::String, info::Vector{Any})
        id = ens_id
        L = info[1]
        T = info[2]
        beta = info[3]
        dtr = info[4]
        vrw = info[5]
 	    cnfg = info[6]
        p0 = 9.2056
        p1 = -13.9847
        g2 = 6 ./ beta 
        ca = - 0.006033 .* g2 .*( 1 .+exp.(p0 .+ p1./g2)) 
        return new(id, L, T, beta, ca, dtr, vrw, cnfg)
    end
    function EnsInfo(ens_id::String)
      if haskey(ens_db,ens_id)
        return EnsInfo(ens_id,ens_db[ens_id])
      else
        error("$ens_id is not in database")
      end
    end
 end
 function Base.show(io::IO, a::GHeader)
    print(io, "ncorr = ", a.ncorr, "\t")
    print(io, "nnoise = ", a.nnoise, "\t")
    print(io, "tvals = ", a.tvals, "\t")
    print(io, "noisetype = ", noise_name[a.noisetype + 1], "\t")
    print(io, "\n")
 end
 function Base.show(io::IO, a::CHeader)
    fnames = fieldnames(CHeader)
    for k in fnames
        f = getfield(a, k)
        if k != :type1 && k!= :type2
            print(io, "$k = $f\t")
        else
            print(io, "$k = ", gamma_name[f + 1], "\t")   
        end
    end
    print(io, "\n")
 end
 function Base.show(io::IO, a::Sm)
    print(io, "\n")
    if a.qg == 1
        print(io, "type = ", qs_name[a.type+1], "\t")
    elseif a.qg == 2
        print(io, "type = ", gs_name[a.type+1], "\t")
    end
    print(io, "niter = ", a.niter, "\t")
    print(io, "eps = ", a.eps, "\t")
    print(io, "\n")
 end
 function Base.show(io::IO, a::CData)
    print(io, a.header)
 end
 function Base.show(io::IO, a::Corr)
    fnames = fieldnames(Corr)
    for k in fnames
        f = getfield(a, k)
        if k != :obs
            print(io, "$k = $f\t")
        end
    end
 end
\ No newline at end of file