Initial commit

.vscode/settings.json

+{}
\ No newline at end of file

plat.txt

+#H400
+ll	40	83
+ls	40	83
+lh	40	83
+ss	40	83
+sh	40	83
+hh	40	83
+#N300
+ll	58	95
+lh	58	95
+hh	58	95
+#N200
+ll	85	105
+ls	85	105
+lh	85	105
+ss	85	105
+sh	83	110
+hh	83	110 
+#N203
+ll	85	105
+ls	85	105
+lh	85	105
+ss	85	105
+sh	83	110
+hh	81	95
+#J303
+ll	120	140
+ls	123	140
+lh	123	140
+ss	123	140
+sh	123	140
+hh	123	155
+#end

results_example/results.txt

+Data analysis results
+ 
+#================================# 
+ 
+Setup:
+Ensembles analysed: ["J303", "H400", "N300"]
+Data storage: /Users/ale/Desktop/data
+Plateaux file: /Users/ale/Google Drive/phd/analysis/automation/plat.txt
+Results storage: /Users/ale/Desktop/results
+Sector analysed: ["sh", "hh", "lh"]
+Reweighting factor: false
+Mass shift: false
+t0 values: extracted from const.jl 
+Phi2 = 8t0m_pi^2: m_pi extracted from ens_db in const.jl 
+ 
+#================================# 
+ 
+
+ Results at finite lattice spacings 
+ 
+Ensemble J303
+$Z^{tm}_M \mu_c \sqrt{8t_0}$ = 3.087441457058999 +/- 0.06439682032069077
+$m_{D_s} \sqrt{8t_0}$ = 4.078057960579361 +/- 0.007655108297108985
+$m_{D_s^*} \sqrt{8t_0}$ = 4.366977977493505 +/- 0.016814805398813727
+$f_{D_s} \sqrt{8t_0}$ = 0.5214236581218687 +/- 0.0030340134610223573
+$f_{D_s^*} \sqrt{8t_0}$ = 0.2032742010074851 +/- 0.003960549693487588
+$m_etac \sqrt{8t_0}$ = 6.199438233112716 +/- 0.010541567771375073
+$m_{J/ψ} \sqrt{8t_0}$ = 6.434900979492144 +/- 0.011365930108534554
+$f_{η_c} \sqrt{8t_0}$ = 0.8796445828751452 +/- 0.0021849815789922664
+$f_{J/ψ} \sqrt{8t_0}$ = 0.40718773063280017 +/- 0.0020634759520040164
+$m_D \sqrt{8t_0}$ = 3.933940114645234 +/- 0.010469403775963788
+$m_{D^*} \sqrt{8t_0}$ = 4.16768216402017 +/- 0.0515103684503361
+$f_D \sqrt{8t_0}$ = 0.4644351651742726 +/- 0.004867644723937044
+$f_{D^*} \sqrt{8t_0}$ = 0.1669527140863556 +/- 0.009789813214623628
+
+Ensemble H400
+$Z^{tm}_M \mu_c \sqrt{8t_0}$ = 3.032580394380745 +/- 0.06662773022505761
+$m_{D_s} \sqrt{8t_0}$ = 3.981978398192009 +/- 0.010337458568118262
+$m_{D_s^*} \sqrt{8t_0}$ = 4.280043260057022 +/- 0.024078918141912688
+$f_{D_s} \sqrt{8t_0}$ = 0.5171746574232967 +/- 0.002516809833076454
+$f_{D_s^*} \sqrt{8t_0}$ = 0.21320483338707588 +/- 0.0038707041022579716
+$m_etac \sqrt{8t_0}$ = 6.101838242920741 +/- 0.013458129500559598
+$m_{J/ψ} \sqrt{8t_0}$ = 6.3555552025914634 +/- 0.014190406675632187
+$f_{η_c} \sqrt{8t_0}$ = 0.9780867087401416 +/- 0.0051276448486433585
+$f_{J/ψ} \sqrt{8t_0}$ = 0.4511410830392943 +/- 0.001641511389275644
+$m_D \sqrt{8t_0}$ = 3.981978398192009 +/- 0.010337458568118262
+$m_{D^*} \sqrt{8t_0}$ = 4.280043260057022 +/- 0.024078918141912688
+$f_D \sqrt{8t_0}$ = 0.5171746574232967 +/- 0.002516809833076454
+$f_{D^*} \sqrt{8t_0}$ = 0.21320483338707588 +/- 0.0038707041022579716
+
+Ensemble N300
+$Z^{tm}_M \mu_c \sqrt{8t_0}$ = 3.0482237405271984 +/- 0.06522961439540911
+$m_{D_s} \sqrt{8t_0}$ = 3.981978398192391 +/- 0.011986238131587756
+$m_{D_s^*} \sqrt{8t_0}$ = 4.328632905699455 +/- 0.027254375933845283
+$f_{D_s} \sqrt{8t_0}$ = 0.5037122885948799 +/- 0.0037756921455479706
+$f_{D_s^*} \sqrt{8t_0}$ = 0.2085589430073183 +/- 0.0035401003493775417
+$m_etac \sqrt{8t_0}$ = 6.155963678990474 +/- 0.010568925225646663
+$m_{J/ψ} \sqrt{8t_0}$ = 6.401012264543191 +/- 0.011576539802341586
+$f_{η_c} \sqrt{8t_0}$ = 0.8900266271330511 +/- 0.003491761532756627
+$f_{J/ψ} \sqrt{8t_0}$ = 0.4165899334662931 +/- 0.0022199719238052123
+$m_D \sqrt{8t_0}$ = 3.981978398192391 +/- 0.011986238131587756
+$m_{D^*} \sqrt{8t_0}$ = 4.328632905699455 +/- 0.027254375933845283
+$f_D \sqrt{8t_0}$ = 0.5037122885948799 +/- 0.0037756921455479706
+$f_{D^*} \sqrt{8t_0}$ = 0.2085589430073183 +/- 0.0035401003493775417
+
+ 
+#================================# 
+ 
+
+ Results in the continuum limit 
+ 
+$Z^{tm}_M \mu_c \sqrt{8t_0}$ = 3.1154197662181335 +/- 0.06434746612878896
+mu_c(MeV) = 1488.514225535831 +/- 14.442336011452733
+
+$m_{D_s} \sqrt{8t_0}$ = 4.118192443077841 +/- 0.015300857438356928
+m_Ds(MeV) = 1967.6282796577461 +/- 25.10390083093398
+
+$m_{D_s^*} \sqrt{8t_0}$ = 4.419014477026606 +/- 0.03184398137612806
+m_Ds_star(MeV) = 2111.357828318512 +/- 29.81362983034594
+
+$f_{D_s} \sqrt{8t_0}$ = 0.5188425758435204 +/- 0.004915964933987481
+f_Ds(MeV) = 247.89743049433383 +/- 3.839298874911889
+
+$f_{D_s^*} \sqrt{8t_0}$ = 0.19762270647855212 +/- 0.006417471690963253
+f_Ds_star(MeV) = 94.42201435324021 +/- 3.268597172693712
+
+$m_etac \sqrt{8t_0}$ = 6.257721044563254 +/- 0.020838906788024733
+m_etac(MeV) = 2989.872149901212 +/- 37.62658949775549
+
+$m_{J/ψ} \sqrt{8t_0}$ = 6.482753815478666 +/- 0.02217903123519892
+m_jpsi(MeV) = 3097.3903997215107 +/- 39.033522687614926
+
+$f_{η_c} \sqrt{8t_0}$ = 0.810029856386968 +/- 0.005231302539793005
+f_etac(MeV) = 387.0235970815643 +/- 5.291825151882144
+
+$f_{J/ψ} \sqrt{8t_0}$ = 0.3776478603806127 +/- 0.003316526386863686
+f_jpsi(MeV) = 180.43610640054987 +/- 2.683138885393678
+
+$m_D \sqrt{8t_0}$ = 3.9138734989952284 +/- 0.018303398304057722
+m_D(MeV) = 1870.0068746351353 +/- 24.178312527760543
+
+$m_{D^*} \sqrt{8t_0}$ = 4.1364685514681 +/- 0.07600617878699434
+m_D_star(MeV) = 1976.3604086701266 +/- 43.30819593423517
+
+$f_D \sqrt{8t_0}$ = 0.4380487738792999 +/- 0.0072962610338359055
+f_D(MeV) = 209.29501650732362 +/- 4.262403623126786
+
+$f_{D^*} \sqrt{8t_0}$ = 0.14612895976638235 +/- 0.014223229623081625
+f_D_star(MeV) = 69.81885321959642 +/- 6.83761644174843
+

src/analysis.jl

+#=
+Software for automated data analysis of Wilson_tm fermions.
+Go to SET UP VARIABLES  section to properly configure before running.
+You have to up set the following variables:
+- path_data where data are stored
+- path_plat where a plat.txt file is stored with plateaus for all ensembles
+- path_result where you want to save results.txt and plots 
+- ensembles to analyse 
+- sector you are interested in: light-light, light-heavy, heavy-heavy...
+- rwf set to true if you want to include reweighting factors
+- compute_t0 set to true if you want to extract t0 from ms.dat file. Otherwise, t0 is taken from 1608.08900 in const.jl
+- mass_shift set to true if you want to include the mass shift
+- go to line 477 to modify continuum limit + chiral extrapolatio fit model
+
+Data in path_data need to be stored as follow:
+data->
+    -ens_id_1->
+             -mesons.dat
+             -ms.dat
+             -ms1.dat
+             -pbp.dat
+    -ens_id_2->
+             -r0_mesons.dat
+             -r1_mesons.dat
+             -r0_ms.dat
+             -r1_ms.dat
+             -r0_ms1.dat
+             -r1_ms1.dat
+             -r0_pbp.dat
+             -r1_pbp.dat
+    -and so on
+=#
+
+using Revise, juobs, PyPlot, LaTeXStrings, ADerrors, DelimitedFiles
+
+#============= SET UP VARIABLES ===========#
+const path_data = "/Users/ale/Desktop/data"
+const path_plat = "/Users/ale/Google Drive/phd/analysis/automation/plat.txt"
+const path_results = "/Users/ale/Desktop/results"
+const ensembles = ["J303", "H400", "N300"]
+const sector = Dict("ll"=>false, "ls"=>false, "lh"=>true, "ss"=>false, "sh"=>true, "hh"=>true )
+const rwf = false
+const compute_t0 = false
+const mass_shift = false #not implemented yet
+#go to line 443 to change CL+chiral extrapolation fit model
+
+#========= DO NOT CHANGE AFTER THIS LINE ==========#
+include("types.jl")
+include("tools.jl")
+include("const.jl")
+
+#======== GET ENSEMBLE INFORMATION FROM DATABASE ==========#
+ensinfo = Vector{EnsInfo}(undef, length(ensembles))
+for i in 1:length(ensembles)
+    ens = ensembles[i]
+    try
+        ensinfo[i]= EnsInfo(ens, ens_db[ens])
+    catch
+        error("The ensemble id ", ens, " was not found in the const.jl ens_db database.
+        Please check the ensemble id or update the database")
+    end
+end
+
+#======== ANALYSIS ==========#
+
+println("\n Computing Observables \n")
+
+    
+@time begin
+    
+    ensobs = Vector{EnsObs}(undef,length(ensembles))
+    for i in 1:length(ensinfo)
+        @time begin
+            ens = ensinfo[i]
+            println("\n Ensemble: ", ens.id)
+            pp = get_corr(ens,"G5", "G5", rw=rwf)
+            a1a1 = get_corr(ens, "G1G5", "G1G5", rw=rwf)
+            mu_list = getfield.(pp, :mu)
+            m_ps = get_meff(pp, ens)
+            m_vec = get_meff(a1a1, ens)
+            f_ps = get_f(pp, m_ps, ens)
+            f_vec = get_f(a1a1, m_vec, ens)
+            ensobs[i] =  EnsObs(ens, mu_list, m_ps, f_ps, m_vec, f_vec)
+            if ens.deg
+                ensobs[i].m_ls = ensobs[i].m_ll
+                ensobs[i].m_ls_vec = ensobs[i].m_ll_vec
+                ensobs[i].m_ss = ensobs[i].m_ll
+                ensobs[i].m_ss_vec = ensobs[i].m_ll_vec
+                ensobs[i].m_sh = ensobs[i].m_lh
+                ensobs[i].m_sh_vec = ensobs[i].m_lh_vec
+                ensobs[i].f_ls = ensobs[i].f_ll
+                ensobs[i].f_ls_vec = ensobs[i].f_ll_vec
+                ensobs[i].f_ss = ensobs[i].f_ll
+                ensobs[i].f_ss_vec = ensobs[i].f_ll_vec
+                ensobs[i].f_sh = ensobs[i].f_lh
+                ensobs[i].f_sh_vec = ensobs[i].f_lh_vec
+            end
+            println(" Time:")
+        end
+    end
+    println("Total time:")
+end
+    
+
+#=========== MATCHING =============#
+
+println("\n Matching Procedure \n")
+
+@time begin
+    for i in 1:length(ensobs)
+        ens = ensobs[i]
+        mul, mus, muh = get_mu(ens.mu_list, ens.ensinfo.deg)
+        if compute_t0
+            est_t0 = get_t0(ens.ensinfo, rw=rwf)
+            ens.t0 = est_t0
+            est_a = t0_ph[1] / sqrt(8 * ens.t0)
+            ens.a = est_a
+        else
+            ens.t0 = t0(ens.ensinfo.beta)
+            ens.a = a(ens.ensinfo.beta)
+        end
+        target = ens.a * (2/3*M[1] + 1/3* M[3])/ hc
+        if !ens.ensinfo.deg
+            ens.muh_target = match_muc(muh, ens.m_lh, ens.m_sh, target)
+        else
+            ens.muh_target = match_muc(muh, ens.m_lh, ens.m_sh, target)
+        end
+        uwerr(ens.muh_target)
+        println(t0, "\n", a, "\n", ens.muh_target)
+        
+        # interpolate lh
+        if sector["lh"]
+            # m_lh
+            par, chi2exp = lin_fit(muh, ens.m_lh)
+            ens.m_lh_match = y_lin_fit(par, value(ens.muh_target))
+            uwerr(ens.m_lh_match)
+            #m_lh_star
+            par, chi2exp = lin_fit(muh, ens.m_lh_vec)
+            ens.m_lh_vec_match = y_lin_fit(par, value(ens.muh_target))
+            uwerr(ens.m_lh_vec_match)
+            # f_lh
+            par, chi2exp = lin_fit(muh, ens.f_lh)
+            ens.f_lh_match = y_lin_fit(par, value(ens.muh_target))
+            uwerr(ens.f_lh_match)
+            #f_lh_star
+            par, chi2exp = lin_fit(muh, ens.f_lh_vec)
+            ens.f_lh_vec_match = y_lin_fit(par, value(ens.muh_target))
+            uwerr(ens.f_lh_vec_match)
+        end
+        # interpolate sh
+        if sector["sh"]
+            if ens.ensinfo.deg
+                ens.m_sh_match = ens.m_lh_match
+                ens.m_sh_vec_match = ens.m_lh_vec_match
+                ens.f_sh_match = ens.f_lh_match
+                ens.f_sh_vec_match = ens.f_lh_vec_match
+            else
+                # m_sh
+                par, chi2exp = lin_fit(muh, ens.m_sh)     
+                ens.m_sh_match = y_lin_fit(par, value(ens.muh_target))
+                uwerr(ens.m_sh_match)
+                # m_sh_star
+                par, chi2exp = lin_fit(muh, ens.m_sh_vec)     
+                ens.m_sh_vec_match = y_lin_fit(par, value(ens.muh_target))
+                uwerr(ens.m_sh_vec_match)
+                # f_sh
+                par, chi2exp = lin_fit(muh, ens.f_sh)     
+                ens.f_sh_match = y_lin_fit(par, value(ens.muh_target))
+                uwerr(ens.f_sh_match)
+                # f_sh_star
+                par, chi2exp = lin_fit(muh, ens.f_sh_vec)     
+                ens.f_sh_vec_match = y_lin_fit(par, value(ens.muh_target))
+                uwerr(ens.f_sh_vec_match)
+            end
+        end
+        # interpolate hh
+        if sector["hh"]
+            # m_hh
+            par, chi2exp = lin_fit(muh, ens.m_hh)
+            ens.m_hh_match = y_lin_fit(par, value(ens.muh_target))
+            uwerr(ens.m_hh_match)
+            # m_hh_vec
+            par, chi2exp = lin_fit(muh, ens.m_hh_vec)
+            ens.m_hh_vec_match = y_lin_fit(par, value(ens.muh_target))
+            uwerr(ens.m_hh_vec_match)
+            # f_hh
+            par, chi2exp = lin_fit(muh, ens.f_hh)
+            ens.f_hh_match = y_lin_fit(par, value(ens.muh_target))
+            uwerr(ens.f_hh_match)
+            # f_hh_vec
+            par, chi2exp = lin_fit(muh, ens.f_hh_vec)
+            ens.f_hh_vec_match = y_lin_fit(par, value(ens.muh_target))
+            uwerr(ens.f_hh_vec_match)
+        end
+    end
+    println(" Total time:")
+end
+    
+    
+#========== MATCHING PLOTS ==========#
+
+println("\n Matchin Plots \n ")
+check_folder = filter(x-> occursin("plots", x), readdir(path_results))
+if !("plots" in check_folder)
+    mkdir(joinpath(path_results,"plots"))
+end
+path_plot = joinpath(path_results,"plots")
+for i in 1:length(ensobs)
+    ens = ensobs[i]
+    mul, mus, muh = get_mu(ens.mu_list, ens.ensinfo.deg)
+    if sector["lh"]
+        # masses
+        figure()
+        title(ens.ensinfo.id)
+        xlabel(L"$a\mu$")
+        ylabel(L"$aM$")
+        #m_lh
+        errorbar(muh, value.(ens.m_lh), yerr=err.(ens.m_lh), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.m_lh_match), yerr=err(ens.m_lh_match), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
+        #m_lh_star
+        errorbar(muh, value.(ens.m_lh_vec), yerr=err.(ens.m_lh_vec), fmt="*", ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.m_lh_vec_match), yerr=err(ens.m_lh_vec_match), fmt="*", ms=5, mfc="none", capsize=2, c="navy", lw=1)
+        legend([L"$m_D$", L"$\langle m_D \rangle$", L"$m_{D^*}$", L"$ \langle m_{D^*}\rangle $"])  
+        display(gcf())
+        name = string(ens.ensinfo.id,"_match_mD_mDvec.pdf")
+        savefig(joinpath(path_plot, name))
+        close()
+        
+        # decays
+        figure()
+        title(ens.ensinfo.id)
+        xlabel(L"$a\mu$")
+        ylabel(L"$af$")
+        #m_lh
+        errorbar(muh, value.(ens.f_lh), yerr=err.(ens.f_lh), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.f_lh_match), yerr=err(ens.f_lh_match), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
+        #m_lh_star
+        errorbar(muh, value.(ens.f_lh_vec), yerr=err.(ens.f_lh_vec), fmt="*", ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.f_lh_vec_match), yerr=err(ens.f_lh_vec_match), fmt="*", ms=5, mfc="none", capsize=2, c="navy", lw=1)
+        legend([L"$f_D$", L"$\langle f_D \rangle$", L"$f_{D^*}$", L"$ \langle f_{D^*}\rangle $"])  
+        display(gcf())
+        name = string(ens.ensinfo.id,"_match_fD_fDvec.pdf")
+        savefig(joinpath(path_plot, name))
+        close()
+    end
+    if sector["sh"] && !ens.ensinfo.deg
+        # masses
+        figure()
+        title(ens.ensinfo.id)
+        xlabel(L"$a\mu$")
+        ylabel(L"$aM$")
+        #m_sh
+        errorbar(muh, value.(ens.m_sh), yerr=err.(ens.m_sh), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.m_sh_match), yerr=err(ens.m_sh_match), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
+        #m_sh_star
+        errorbar(muh, value.(ens.m_sh_vec), yerr=err.(ens.m_sh_vec), fmt="*", ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.m_sh_vec_match), yerr=err(ens.m_sh_vec_match), fmt="*", ms=5, mfc="none", capsize=2, c="navy", lw=1)
+        legend([L"$m_{D_s}$", L"$\langle m_{D_s} \rangle$", L"$m_{D^*_s}$", L"$ \langle m_{D^*_s}\rangle $"])  
+        display(gcf())
+        name = string(ens.ensinfo.id,"_match_mDs_mDsvec.pdf")
+        savefig(joinpath(path_plot, name))
+        close()
+        
+        # decays
+        figure()
+        title(ens.ensinfo.id)
+        xlabel(L"$a\mu$")
+        ylabel(L"$af$")
+        #f_sh
+        errorbar(muh, value.(ens.f_sh), yerr=err.(ens.f_sh), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.f_sh_match), yerr=err(ens.f_sh_match), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
+        #fsh_star
+        errorbar(muh, value.(ens.f_sh_vec), yerr=err.(ens.f_sh_vec), fmt="*", ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.f_sh_vec_match), yerr=err(ens.f_sh_vec_match), fmt="*", ms=5, mfc="none", capsize=2, c="navy", lw=1)
+        legend([L"$f_{D_s}$", L"$\langle f_{D_s} \rangle$", L"$f_{D^*_s}$", L"$ \langle f_{D^*_s}\rangle $"])  
+        display(gcf())
+        name = string(ens.ensinfo.id,"_match_fDs_fDsvec.pdf")
+        savefig(joinpath(path_plot, name))
+        close()
+    end
+    if sector["hh"] 
+        # masses
+        figure()
+        title(ens.ensinfo.id)
+        xlabel(L"$a\mu$")
+        ylabel(L"$aM$")
+        #m_hh
+        errorbar(muh, value.(ens.m_hh), yerr=err.(ens.m_hh), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.m_hh_match), yerr=err(ens.m_hh_match), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
+        #m_hh_star
+        errorbar(muh, value.(ens.m_hh_vec), yerr=err.(ens.m_hh_vec), fmt="*", ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.m_hh_vec_match), yerr=err(ens.m_hh_vec_match), fmt="*", ms=5, mfc="none", capsize=2, c="navy", lw=1)
+        legend([L"$m_{\eta_c}$", L"$\langle m_{\eta_c} \rangle$", L"$m_{J/\psi}$", L"$ \langle m_{J/\psi}\rangle $"])  
+        display(gcf())
+        name = string(ens.ensinfo.id,"_match_metac_mJpsi.pdf")
+        savefig(joinpath(path_plot, name))
+        close()
+        
+        # decays
+        figure()
+        title(ens.ensinfo.id)
+        xlabel(L"$a\mu$")
+        ylabel(L"$af$")
+        #f_hh
+        errorbar(muh, value.(ens.f_hh), yerr=err.(ens.f_hh), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.f_hh_match), yerr=err(ens.f_hh_match), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
+        #f_hh_star
+        errorbar(muh, value.(ens.f_hh_vec), yerr=err.(ens.f_hh_vec), fmt="*", ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
+        errorbar(value(ens.muh_target), xerr=err(ens.muh_target), value(ens.f_hh_vec_match), yerr=err(ens.f_hh_vec_match), fmt="*", ms=5, mfc="none", capsize=2, c="navy", lw=1)
+        legend([L"$f_{\eta_c}$", L"$\langle f_{\eta_c} \rangle$", L"$f_{J/\psi}$", L"$ \langle f_{J/\psi}\rangle $"])  
+        display(gcf())
+        name = string(ens.ensinfo.id,"_match_fetac_fJpsi.pdf")
+        savefig(joinpath(path_plot, name))
+        close()
+    end
+end
+
+#======= RESULTS ==========#
+# charm quark mass 
+muc = zm_tm.(getfield.(ensinfo,:beta)) .* getfield.(ensobs,:muh_target) .* sqrt.(8 * getfield.(ensobs,:t0))
+uwerr.(muc)
+
+# ll
+m_pi = Vector{uwreal}(undef, length(ensembles))
+m_rho =Vector{uwreal}(undef, length(ensembles))
+f_pi = Vector{uwreal}(undef, length(ensembles))
+f_rho =Vector{uwreal}(undef, length(ensembles))
+if sector["ll"]
+    try 
+        m_pi .= getfield.(ensobs, :m_ll) .* sqrt.(8 * getfield.(ensobs,:t0))
+        m_rho .= getfield.(ensobs, :m_ll_vec) .* sqrt.(8 * getfield.(ensobs,:t0))
+        f_pi .= getfield.(ensobs, :f_ll) .* sqrt.(8 * getfield.(ensobs,:t0))
+        f_rho .= getfield.(ensobs, :f_ll_vec) .* sqrt.(8 * getfield.(ensobs,:t0))        
+        uwerr.(m_pi)
+        uwerr.(m_rho)
+        uwerr.(f_pi)
+        uwerr.(f_rho)
+    catch
+        error("Check ensobs[i].mu_list to see whether all the ensembles have a light-light sector.")
+    end
+end
+# ls
+m_k = Vector{uwreal}(undef, length(ensembles))
+m_k_star =Vector{uwreal}(undef, length(ensembles))
+f_k= Vector{uwreal}(undef, length(ensembles))
+f_k_star =Vector{uwreal}(undef, length(ensembles))
+if sector["ls"] 
+    try 
+        m_k .= getfield.(ensobs, :m_ls) .* sqrt.(8 * getfield.(ensobs,:t0))
+        m_k_star .= getfield.(ensobs, :m_ls_vec) .* sqrt.(8 * getfield.(ensobs,:t0))
+        f_k .= getfield.(ensobs, :f_ls) .* sqrt.(8 * getfield.(ensobs,:t0))
+        f_k_star .= getfield.(ensobs, :f_ls_vec) .* sqrt.(8 * getfield.(ensobs,:t0))        
+        uwerr.(m_k)
+        uwerr.(m_k_star)
+        uwerr.(f_k)
+        uwerr.(f_k_star)
+    catch
+        error("Check ensobs[i].mu_list to see whether all the ensembles have a light-stange sector.")
+    end
+end
+# lh
+m_D = Vector{uwreal}(undef, length(ensembles))
+m_D_star =Vector{uwreal}(undef, length(ensembles))
+f_D= Vector{uwreal}(undef, length(ensembles))
+f_D_star =Vector{uwreal}(undef, length(ensembles))
+if sector["lh"] 
+    try 
+        m_D .= getfield.(ensobs, :m_lh_match) .* sqrt.(8 * getfield.(ensobs,:t0))
+        m_D_star .= getfield.(ensobs, :m_lh_vec_match) .* sqrt.(8 * getfield.(ensobs,:t0))
+        f_D .= getfield.(ensobs, :f_lh_match) .* sqrt.(8 * getfield.(ensobs,:t0))
+        f_D_star .= getfield.(ensobs, :f_lh_vec_match) .* sqrt.(8 * getfield.(ensobs,:t0))        
+        uwerr.(m_D)
+        uwerr.(m_D_star)
+        uwerr.(f_D)
+        uwerr.(f_D_star)
+catch
+        error("Check ensobs[i].mu_list to see whether all the ensembles have a light-heavy sector.")
+    end
+end
+# ss
+m_etaprime = Vector{uwreal}(undef, length(ensembles))
+m_phi =Vector{uwreal}(undef, length(ensembles))
+f_etaprime= Vector{uwreal}(undef, length(ensembles))
+f_phi =Vector{uwreal}(undef, length(ensembles))
+if sector["ss"] 
+    try 
+        m_etaprime .= getfield.(ensobs, :m_ss) .* sqrt.(8 * getfield.(ensobs,:t0))
+        m_phi .= getfield.(ensobs, :m_ss_vec) .* sqrt.(8 * getfield.(ensobs,:t0))
+        f_etaprime .= getfield.(ensobs, :f_ss) .* sqrt.(8 * getfield.(ensobs,:t0))
+        f_phi .= getfield.(ensobs, :f_ss_vec) .* sqrt.(8 * getfield.(ensobs,:t0))        
+        uwerr.(m_etaprime)
+        uwerr.(m_phi)
+        uwerr.(f_etaprime)
+        uwerr.(f_phi)
+catch
+        error("Check ensobs[i].mu_list to see whether all the ensembles have a strange-strange sector.")
+    end
+end
+# sh 
+m_Ds = Vector{uwreal}(undef, length(ensembles))
+m_Ds_star =Vector{uwreal}(undef, length(ensembles))
+f_Ds= Vector{uwreal}(undef, length(ensembles))
+f_Ds_star =Vector{uwreal}(undef, length(ensembles))
+if sector["sh"] 
+    try 
+        m_Ds .= getfield.(ensobs, :m_sh_match) .* sqrt.(8 * getfield.(ensobs,:t0))
+        m_Ds_star .= getfield.(ensobs, :m_sh_vec_match) .* sqrt.(8 * getfield.(ensobs,:t0))
+        f_Ds .= getfield.(ensobs, :f_sh_match) .* sqrt.(8 * getfield.(ensobs,:t0))
+        f_Ds_star .= getfield.(ensobs, :f_sh_vec_match) .* sqrt.(8 * getfield.(ensobs,:t0))        
+        uwerr.(m_Ds)
+        uwerr.(m_Ds_star)
+        uwerr.(f_Ds)
+        uwerr.(f_Ds_star)
+catch
+        error("Check ensobs[i].mu_list to see whether all the ensembles have a strange-heavy sector.")
+    end
+end
+# hh 
+m_etac = Vector{uwreal}(undef, length(ensembles))
+m_jpsi =Vector{uwreal}(undef, length(ensembles))
+f_etac= Vector{uwreal}(undef, length(ensembles))
+f_jpsi =Vector{uwreal}(undef, length(ensembles))
+if sector["hh"] 
+    try 
+        m_etac .= getfield.(ensobs, :m_hh_match) .* sqrt.(8 * getfield.(ensobs,:t0))
+        m_jpsi .= getfield.(ensobs, :m_hh_vec_match) .* sqrt.(8 * getfield.(ensobs,:t0))
+        f_etac .= getfield.(ensobs, :f_hh_match) .* sqrt.(8 * getfield.(ensobs,:t0))
+        f_jpsi .= getfield.(ensobs, :f_hh_vec_match) .* sqrt.(8 * getfield.(ensobs,:t0))        
+        uwerr.(m_etac)
+        uwerr.(m_jpsi)
+        uwerr.(f_etac)
+        uwerr.(f_jpsi)
+catch
+        error("Check ensobs[i].mu_list to see whether all the ensembles have a heavy-heavy sector.")
+    end
+end
+
+#===== CONTINUUM AND CHIRAL EXTRAPOLATION =======#
+
+println("\n Continuum and Chiral Extrapolation \n")
+obs_db =Dict(
+    "muc"=> [muc],
+    "ll" => [m_pi, m_rho, f_pi, f_rho],
+    "ls" => [m_k, m_k_star, f_k, f_k_star],
+    "lh" => [m_D, m_D_star, f_D, f_D_star],
+    "ss" => [m_etaprime, m_phi, f_etaprime, f_phi],
+    "sh" => [m_Ds, m_Ds_star, f_Ds, f_Ds_star],
+    "hh" => [m_etac, m_jpsi, f_etac, f_jpsi]
+)
+
+
+#preparing observables of interest (data, name, ylabel)
+obs = Vector{Vector{uwreal}}(undef,0)
+ttl_obs = Vector{String}(undef,0)
+ylbl = Vector{String}(undef,0)
+push!(obs, obs_db["muc"][1])
+push!(ttl_obs, ttl_obs_db["muc"][1])
+push!(ylbl, ylbl_db["muc"][1])
+
+label = findall(sector)
+for lab in label
+    for i in 1:length(obs_db[lab])
+        push!(obs, obs_db[lab][i])
+        push!(ttl_obs, ttl_obs_db[lab][i])
+        push!(ylbl, ylbl_db[lab][i])
+    end
+end
+# estimate phi2 from database or form analysis if ll sector is computed 
+phi2 = Vector{uwreal}(undef,length(ensembles))
+sector["ll"] ? phi2 = 8 .* getfield.(ensobs,:t0) .* getfield.(ensobs,:m_ll).^2 : phi2 = 8 .* getfield.(ensobs,:t0) .* getfield.(getfield.(ensobs,:ensinfo),:mpi).^2
+phi2_ph = (t0_ph[1] * 139.57039 / hc)^2
+uwerr(phi2_ph)
+#fit routine 
+x = [1 ./(8 .* getfield.(ensobs,:t0)) phi2] #x[1] = a^2 / (8t0), x[2] = 8t0 mpi^2
+uwerr.(x)
+#modify CL+chiral fit  model
+@. cl_chir_model(x, p) = (p[1] + p[2] * x[:, 1]) + (p[3]) * x[:, 2] #linear fits
+
+obs_t0 = Vector{uwreal}(undef, length(obs)) #sqrt(8t0)obs @ CL & phi2_phys
+xarr = Float64.(range(0.0, stop=0.05, length=100))
+for k = 1:length(obs)
+    par = fit_routine(cl_chir_model, value.(x), obs[k], 4)
+    obs_t0[k] = par[1] + par[3] * phi2_ph
+    uwerr(obs_t0[k])    
+    y = Vector{uwreal}(undef, 100)
+    for j=1:100
+        y[j] = par[1] + par[2] * xarr[j] + par[3] * phi2_ph
+    end
+    uwerr.(y)
+    figure()
+    for b in unique(phi2) #select point with same beta
+        nn = findall(x-> x == b, phi2)
+        lgnd = string(L"$m_\pi = $", Int32(round(sqrt(value(b))*hc/t0_ph_value[1], sigdigits=3)), " MeV")
+        errorbar(value.(x[nn,1]), value.(obs[k][nn]), err.(obs[k][nn]), ms=5, fmt="s", mfc="none", mew=0.8, elinewidth=0.8, capsize=2, label=lgnd)
+    end
+    fill_between(xarr, value.(y) .+ err.(y),value.(y) .- err.(y) , color="tomato", alpha=0.4)
+    lgnd=L"$\mathrm{CL}$"
+    errorbar(0, value(obs_t0[k]), err(obs_t0[k]), fmt="s", zorder=2, ms=4,  mfc="red", c = "red", mew=0.8, elinewidth=1, capsize=2, label=lgnd)
+    axvline(0, ls="--", color="black", zorder=1, lw=0.6, label="")
+    xlabel(L"$a^2/8t_0$")
+    ylabel(ylbl[k])
+    xlim(-0.002,0.04)
+    #ylim(2.8,3.4)
+    legend(ncol=2)
+    display(gcf())
+    t = string("fit_", ttl_obs[k], ".pdf")
+    savefig(joinpath(path_plot, t))
+    #savefig(joinpath("/Users/ale/Desktop/plottest", t))
+    close()
+
+end
+
+obs_ph = Vector{uwreal}(undef, length(obs))
+for k = 1:length(obs)
+    println(ylbl[k], " = ", obs_t0[k])
+    #phys
+    obs_ph[k] = obs_t0[k] * hc / t0_ph[1]
+    uwerr(obs_ph[k])
+    println(ttl_obs[k], "(MeV) = ", obs_ph[k])
+end
+
+#======== STORE RESULTS ========#
+
+open(joinpath(path_results, "results.txt"), "w") do f 
+    print(f, "Data analysis results", "\n \n")
+    print(f, "#================================# \n \n")
+    print(f, "Setup:\n")
+    print(f, "Ensembles analysed: ", ensembles, "\n")
+    print(f, "Data storage: ", path_data, "\n")
+    print(f, "Plateaux file: ", path_plat, "\n")
+    print(f, "Results storage: ", path_results, "\n")
+    print(f, "Sector analysed: ", findall(sector), "\n")
+    print(f, "Reweighting factor: ", rwf, "\n")
+    print(f, "Mass shift: ", mass_shift, "\n")
+    if compute_t0
+        print(f, "t0 values: computed from the ms.dat files \n")
+    else
+        print(f, "t0 values: extracted from const.jl \n")
+    end
+    if sector["ll"]
+        print(f, "Phi2 = 8t0m_pi^2: m_pi computed from mesons.dat files \n")
+    else
+        print(f, "Phi2 = 8t0m_pi^2: m_pi extracted from ens_db in const.jl \n")
+    end
+    print(f, " \n#================================# \n \n")
+    print(f, "\n Results at finite lattice spacings \n \n")
+
+    lab = findall(sector)
+    for i in 1:length(ensembles)
+        print(f, "Ensemble ", ensobs[i].ensinfo.id, "\n")
+        for j in 1:length(obs)
+            print(f, ylbl[j], " = ", obs[j][i], "\n")
+        end
+        print(f, "\n")
+    end
+
+    print(f, " \n#================================# \n \n")
+    print(f, "\n Results in the continuum limit \n \n")
+    for k = 1:length(obs)
+        print(f, ylbl[k], " = ", obs_t0[k], "\n")
+        #phys
+        print(f, ttl_obs[k], "(MeV) = ", obs_ph[k], "\n\n")
+    end
+end
\ No newline at end of file

src/const.jl

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

src/tools.jl

+@doc raw"""
+    read_data(ens::String, g1::String="G5", g2::String="G5")
+Given the ensemble id and the Dirac structure, this method reads the data in path_data and
+returns a CData object. Replicas are automatically detected and loaded if present.
+"""
+function read_data(ens::String, g1::String="G5", g2::String="G5"; legacy::Bool=false)
+    path = joinpath(path_data, ens)
+    rep = filter(x->occursin("mesons.dat", x), readdir(path, join=true))
+    if length(rep)!=0    
+        length(rep)==1 ? (return read_mesons(rep[1], g1, g2, legacy=legacy)) : (return read_mesons(rep, g1, g2, legacy=legacy))
+    else
+        error("mesons.dat file not found for ensemble ", ens, " in path ", path)
+    end
+end
+function read_rw(ens::String)
+    path = joinpath(path_data, ens)
+    rep = filter(x->occursin("ms1.dat", x), readdir(path, join=true))
+    if length(rep)!=0
+        try
+            length(rep) == 1 ? (return read_ms1(rep[1])) : (return read_ms1.(rep)) 
+        catch
+            length(rep) == 1 && length(rep)!=0 ? (return read_ms1(rep[1], v="1.4")) : (return read_ms1.(rep, v="1.4")) 
+        end
+    else
+        error("ms1.dat file not found for ensemble ", ens, " in path ", path)
+    end
+end
+function read_t0(ens::String)
+    path = joinpath(path_data, ens)
+    rep = filter(x->occursin("ms.dat", x), readdir(path, join=true))
+    if length(rep)!=0
+        length(rep) == 1 ? (return read_ms(rep[1])) : (return read_ms.(rep)) 
+    else
+        error("ms.dat file not found for ensemble ", ens, " in path ", path)
+    end
+end
+function get_corr(ens::EnsInfo, g1::String="G5", g2::String="G5"; rw::Bool=false, legacy::Bool=false)
+    aux = read_data(ens.id, g1, g2, legacy=legacy)
+    !rw ? obs = corr_obs.(aux, L=ens.L) : obs = corr_obs.(aux, L=ens.L, rw=read_rw(ens.id))
+    return obs
+end
+function get_t0(ens::EnsInfo; rw::Bool=false, pl::Bool=false)
+    data = read_t0(ens.id)
+    plat = [40,160] 
+    !rw ? t0 = comp_t0(data, plat, L=ens.L, pl=pl) : t0 = comp_t0(data, plat, L=ens.L, pl=pl, rw=read_rw(ens.id))
+    return t0
+end
+function get_mu(mu_list::Vector{Vector{Float64}}, deg::Bool)
+    mu_sorted = unique(sort(minimum.(mu_list)))
+
+    mul = mu_sorted[1]
+    deg ? mus = 0.0 : mus = mu_sorted[2]
+    muh = unique(maximum.(mu_list))
+    muh = filter(x-> x > mul && x > mus, muh)
+
+    return mul, mus, muh
+end
+function get_ll(mu_list, obs::Vector{uwreal}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mul in mu && mu[1] == mu[2] #l-l
+            return obs[k]
+        end
+    end
+end
+function get_ls(mu_list, obs::Vector{uwreal}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mul in mu && mus in mu  #l-l
+            return obs[k]
+        end
+    end
+end
+function get_lh(mu_list, obs::Vector{uwreal}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    obs_lh = Vector{uwreal}(undef, 0)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mul in mu && mu[1] != mu[2] && !(mus in mu) #l-h
+            push!(obs_lh, obs[k])
+        end
+    end
+    return obs_lh
+end
+function get_ss(mu_list, obs::Vector{uwreal}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mus in mu && mu[1] == mu[2] #s-s
+            return obs[k]
+        end
+    end
+end
+function get_sh(mu_list, obs::Vector{uwreal}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    obs_sh = Vector{uwreal}(undef, 0)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mus in mu && mu[1] != mu[2] && !(mul in mu)#s-h
+            push!(obs_sh, obs[k])
+        end
+    end
+    return obs_sh
+end
+function get_hh(mu_list, obs::Vector{uwreal}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    obs_hh = Vector{uwreal}(undef,0)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        for i =1:length(muh)
+            if muh[i] in mu && mu[1] == mu[2]#h-h
+                push!(obs_hh, obs[k])
+            end
+        end
+    end
+    return obs_hh
+end
+function get_meff(obs::Vector{juobs.Corr}, ens::EnsInfo; pl::Bool=false)
+    mu_list  = getfield.(obs, :mu)
+    plat = select_plateau(ens, mu_list)
+    return meff.(obs, plat, pl=pl)
+end
+function get_f(obs::Vector{juobs.Corr}, m::Vector{uwreal}, ens::EnsInfo; pl::Bool=false)
+    mu_list = getfield.(obs, :mu)
+    plat = select_plateau(ens, mu_list)
+    return dec_const_pcvc.(obs, plat, m, pl=pl)
+end
+function select_plateau(ensinf::EnsInfo, mu_list)
+    ens =ensinf.id
+    deg = ensinf.deg
+    mul, mus, muh = get_mu(mu_list, deg)
+    f = readdlm(path_plat)
+    head = String.(f[:, 1])
+
+    delim = findall(x-> occursin("#", x ), head)
+    edelim = findfirst(x-> occursin(ens, x), head)
+    cdelim = findfirst(x-> x.== edelim, delim)
+
+    del = delim[cdelim]+1 : delim[cdelim+1]-1
+    #del = delim
+    plat = Vector{Vector{Int64}}(undef, 0)
+    #plat = Dict()
+    plat_ = Int64.(f[del, 2:3])
+    head_ = String.(f[del, 1])
+    
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mul in mu && mu[1] != mu[2] && !(mus in mu)  #heavy-light
+            n = findfirst(x-> occursin("lh", x), head_)
+            push!(plat,  [plat_[n, 1], plat_[n, 2]]  ) 
+        elseif mul in mu && mu[1] == mu[2] #light-light
+            n = findfirst(x-> occursin("ll", x), head_)
+            push!(plat,  [plat_[n, 1], plat_[n, 2]] )               
+        elseif mul in mu && mus in mu#strange-light
+            n = findfirst(x-> occursin("ls", x), head_)
+            push!(plat,  [plat_[n, 1], plat_[n, 2]] )
+        elseif mus in mu && mu[1] != mu[2] && !(mul in mu)#heavy-strange
+            n = findfirst(x-> occursin("sh", x ), head_)
+            push!(plat, [plat_[n, 1], plat_[n, 2]] )
+        elseif mus in mu && mu[1] == mu[2] && !(mul in mu)#strange-strange
+            n = findfirst(x-> occursin("ss", x), head_)
+            push!(plat, [plat_[n, 1], plat_[n, 2]] )
+        elseif !(mul in mu) && !(mus in mu) #heavy-heavy
+            n = findfirst(x-> occursin("hh", x), head_)
+            push!(plat, [plat_[n, 1], plat_[n, 2]] )
+        end
+    end
+    return plat
+end
+
+function match_muc(muh, m_lh, m_sh, target)
+    M = (2/3 .* m_lh .+ 1/3 .* m_sh) 
+    par, chi2exp = lin_fit(muh, M)
+    muh_target = x_lin_fit(par, target)
+    return muh_target
+end
\ No newline at end of file

src/types.jl

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