now the matching of the charm quark mass is performed as a simultaneous fit...

now the matching of the charm quark mass is performed as a simultaneous fit and not ensemble by ensemble as it was previously

src/.vscode/settings.json

+{}
\ No newline at end of file

src/const.jl

@@ -44,8 +44,8 @@ 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
+const t0_ph_value = [0.415]
+const t0_ph_error = ones(1,1) .* 4e-3
 # 1808.09236
 const ZA_data = [0.75642, 0.76169, 0.76979, 0.78378, 0.79667]
 const ZA_err = [72, 93, 43, 47, 47] .*1e-5

src/func_comb.jl

+
+func_map = [Bool.([i,j,k,m,n]) for i=0:1 for j=0:1 for k=0:1 for m=0:1 for n=0:1]
+
+basemodel(x,p) = p[1] .+ p[2] .* x[:,2] .+ p[3] .* x[:,3]  .+ p[4] .* x[:,1]
+a2phi2_2(x) = x[:,1] .* x[:,2].^2  #a^2/8t0*phi2^2 
+a2phi2(x) = x[:,1] .* x[:,2]      # phi2*a^2/8t0
+a2phih2(x) = x[:,1] .* x[:,3].^2  # phih^2*a^2/8t0          
+a3phih3(x) = x[:,1].^(4) .* x[:,3].^4  #phih^3*(a^2/8t0)^3/2          
+a3cutoff(x) = x[:,1].^(4)   #(a^2/8t0)^4
+
+func_list = [a2phi2_2, a2phi2, a2phih2, a3phih3, a3cutoff]
+npar = 5 # number of extra parameters 
+n_lin_param = Vector(4:npar+4) #param for each functions
+f_lin = Vector{Vector{Function}}(undef, npar+1) #f[i][j]-> i: number of extra parameters, j: combinations
+f_lin[1] = Vector{Function}(undef, 1)
+f_lin[1][1] = (x,p) -> basemodel(x,p)
+for n = 2:npar+1
+    aux = filter(x->sum(x) == n-1, func_map)
+    f_lin[n] = Vector{Function}(undef, length(aux))
+    k=1
+    for a in aux 
+        #vectorized function 
+        f_lin[n][k] = (x,p) -> basemodel(x,p) .+ sum([p[i+4] for i=1:(n-1)] .* (fill(x,n-1) .|> func_list[a]))
+        k = k+1
+    end
+end
+
+f_non_lin = Vector{Vector{Function}}(undef, npar) #f[i][j]-> i: number of extra parameters, j: combinations
+n_non_lin_param = Vector(5:npar+4) #param for each non-linear functions
+for n = 1:npar
+    aux = filter(x->sum(x) == n, func_map)
+    f_non_lin[n] = Vector{Function}(undef, length(aux))
+    k=1
+    for a in aux 
+        #vectorized function 
+        f_non_lin[n][k] = (x,p) -> basemodel(x,p) .* (1 .+ sum([p[i+4] for i=1:(n)] .* (fill(x,n) .|> func_list[a])))
+        k = k+1
+    end
+end
+
+
+
+f_basemodel(x,p) = p[1] .+ p[2] .* x[:,2] .+ p[3] ./ sqrt.(x[:,3]) .+ p[4] .* x[:,1]
+f_a2phi2_2(x) = x[:,1] .* x[:,2].^2              #a^2/8t0 
+f_a2phi2(x) = x[:,1] .* x[:,2]      # phi2*a^2/8t0
+f_a2phih2(x) = x[:,1] ./ x[:,3]  # 1/phih*a^2/8t0          
+f_a3phih3(x) = x[:,1].^(4) ./ x[:,3].^2  #1/phih^3/2*(a^2/8t0)^3/2          
+f_a3cutoff(x) = x[:,1].^(3/2)   #(a^2/8t0)^3/2
+
+f_func_list = [f_a2phi2_2, f_a2phi2, f_a2phih2, f_a3phih3, f_a3cutoff]
+f_npar = 5 # number of extra parameters 
+f_f_lin = Vector{Vector{Function}}(undef, npar+1) #f[i][j]-> i: number of extra parameters, j: combinations
+f_f_lin[1] = Vector{Function}(undef, 1)
+f_f_lin[1][1] = (x,p) -> basemodel(x,p)
+for n = 2:npar+1
+    aux = filter(x->sum(x) == n-1, func_map)
+    f_f_lin[n] = Vector{Function}(undef, length(aux))
+    k=1
+    for a in aux 
+        #vectorized function 
+        f_f_lin[n][k] = (x,p) -> f_basemodel(x,p) .+ sum([p[i+4] for i=1:(n-1)] .* (fill(x,n-1) .|> f_func_list[a]))
+        k = k+1
+    end
+end
+
+f_f_non_lin = Vector{Vector{Function}}(undef, npar) #f[i][j]-> i: number of extra parameters, j: combinations
+n_non_lin_param = Vector(5:npar+4) #param for each non-linear functions
+for n = 1:npar
+    aux = filter(x->sum(x) == n, func_map)
+    f_f_non_lin[n] = Vector{Function}(undef, length(aux))
+    k=1
+    for a in aux 
+        #vectorized function 
+        f_f_non_lin[n][k] = (x,p) -> f_basemodel(x,p) .* (1 .+ sum([p[i+4] for i=1:(n)] .* (fill(x,n) .|> f_func_list[a])))
+        k = k+1
+    end
+end
+
+n_param = n_lin_param
+append!(n_param, n_non_lin_param)
+f = f_lin
+append!(f, f_non_lin)
+f_dec = f_f_lin 
+append!(f_dec, f_f_non_lin)
\ No newline at end of file

src/tools.jl

@@ -7,14 +7,13 @@ function read_data(ens::String, g1::String="G5", g2::String="G5"; legacy::Bool=f
     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))
+        length(rep)==1 ? (return read_mesons(rep[1], g1, g2, id=ens, legacy=legacy)) : (return read_mesons(rep, g1, g2, id=ens, 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))
+    rep = filter(x->occursin(ens, x), readdir(path_rw, join=true))
     if length(rep)!=0
         try
             length(rep) == 1 ? (return read_ms1(rep[1])) : (return read_ms1.(rep)) 
@@ -22,7 +21,19 @@ function read_rw(ens::String)
             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)
+        error("ms1.dat file not found for ensemble ", ens, " in path ", path_rw)
+    end
+end
+function read_mass_dev(ens::String)
+    rep = filter(x->occursin(ens,x), readdir(path_md, join=true))
+    try
+        if length(rep) == 1
+            return [read_md(rep[1])]
+        else
+            return [read_md(rep[1]), read_md(rep[2])]
+        end
+    catch
+        error("pbp.dat file not found for ensemble ", ens, " in path ", path_md)
     end
 end
 function read_t0(ens::String)
@@ -36,6 +47,12 @@ function read_t0(ens::String)
 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)
+    if ens.id =="J303"
+        truncate_data!(aux, 739)
+    end
+    if ens.id =="H101"
+        truncate_data!(aux, [878, 999])
+    end
     !rw ? obs = corr_obs.(aux, L=ens.L) : obs = corr_obs.(aux, L=ens.L, rw=read_rw(ens.id))
     return obs
 end
@@ -160,7 +177,7 @@ function fit_mass(en_i::Array{uwreal}; t_in::Int64=3, wpm::Union{Nothing, Dict{I
     @.model(x,p) = p[1] + p[2]*exp(-x*p[3])
     ydata = en_i[t_in:end-5]
     xdata = collect(1:length(ydata))
-    par = fit_routine(model, xdata, ydata, 3, wpm=wpm)
+    par, _ = fit_routine(model, xdata, ydata, 3, wpm=wpm)
     println("\nPar[1] is ", par[1] )
     return par[1]
 end
@@ -173,9 +190,6 @@ function gevp_dec(mat_obs::Vector{MatInfo}, mass::Vector{uwreal}; t0::Int64=2, w
         evecs = getall_eig(mat, t0)
         elem = mat_elem(evecs, mat, mass[i], n) #ground and excited state energy
         plat = select_plateau(mat_obs[i].ensinfo, mu_list)[end] .- y0s[i]
-        println("plateau = ", plat)
-        println("\n elem = ", elem)
-        println("\n mu_list = ", mu_list)
         try
             if pseudo
                 plat_dec0[i] = -dec(elem, plat, mass[i], mu_list[i], pl=pl)
@@ -323,66 +337,47 @@ function match_muc_heavymass(muh, m_hh, target)
     muh_target = x_lin_fit(par, target)
     return muh_target
 end
-
-function get_ll(mu_list, obs::Array{juobs.Corr}, deg::Bool)
-    mul, mus, muh = get_mu(mu_list, deg)
-    for k = 1:length(mu_list)
-        mu = mu_list[k]
-        if mul in mu && mu[1] == mu[2] #l-l
-            return obs[k]
-        end
-    end
-end
-function get_ls(mu_list, obs::Array{juobs.Corr}, deg::Bool)
-    mul, mus, muh = get_mu(mu_list, deg)
-    for k = 1:length(mu_list)
-        mu = mu_list[k]
-        if mul in mu && mus in mu  #l-l
-            return obs[k]
-        end
-    end
+function match_phih(muh,phih, phih_ph)
+    par, chi2exp = lin_fit(muh, phih)
+    muh_target = x_lin_fit(par, phih_ph)
+    
+    return muh_target
 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
+
+function aic_weigth(cat::Cat)
+    w_aux = exp.(-1/2 * getfield(cat, :aic))
+    w_aux1=filter(x->x<10000, w_aux)
+    index = findall(x->x<10000, w_aux)
+    norm = sum(w_aux1)
+    w = w_aux1 ./ norm
+    return w, index
+end
+function aic_syst(cat::Cat)
+    w,index = aic_weigth(cat)
+    aux1 = sum( w .* getfield(cat,:fit_res)[index].^2)
+    aux2 = sum( w .* getfield(cat,:fit_res)[index])^2
+    return sqrt(abs(value(aux1 - aux2)))
+end
+function aic_model_ave(cat::Cat)
+    w, index = aic_weigth(cat)
+    res = sum( w.* getfield(cat, :fit_res)[index])
+    uwerr(res)
+    return res
 end
-function get_ss(mu_list, obs::Array{juobs.Corr}, deg::Bool)
-    mul, mus, muh = get_mu(mu_list, deg)
-    for k = 1:length(mu_list)
-        mu = mu_list[k]
-        if mus in mu && mu[1] == mu[2] #s-s
-            return obs[k]
-        end
-    end
+function category_av(cat_arr::Vector{Cat})
+    best_res = aic_model_ave.(cat_arr)
+    syst = aic_syst.(cat_arr)
+    tot_err = err.(best_res).^2 .+ syst.^2
+    w_aux = 1 ./ (tot_err)
+    w = w_aux ./ sum(w_aux)
+    aux1 = sum( w .* best_res.^2)
+    aux2 = sum( w .* best_res)^2
+    av = sum( w .* best_res) 
+    return av, sqrt(abs(value(aux1 - aux2)))
 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
+
+function syst_av(syst_err::Vector{Float64})
+    aux = sum(syst_err.^2)
+    return sqrt(aux)
 end
-function get_hh(mu_list, obs::Array{juobs.Corr}, deg::Bool)
-    mul, mus, muh = get_mu(mu_list, deg)
-    obs_hh = Array{juobs.Corr}(undef,0)
-    for k = 1:length(mu_list)
-        mu = mu_list[k]
-        for i =1:length(muh)
-            if muh[i] in mu && mu[1] == mu[2]#h-h
-                push!(obs_hh, obs[k])
-            end
-        end
-    end
-    return obs_hh
-end
\ No newline at end of file
+

src/types.jl

@@ -61,6 +61,7 @@ mutable struct EnsObs
     muh_target::Union{Nothing,uwreal}
     a::Union{Nothing,uwreal}
     t0::Union{Nothing,uwreal}
+    deltam::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})
@@ -112,6 +113,7 @@ mutable struct EnsObs
         a.muh_target = nothing
         a.a = nothing
         a.t0 = nothing
+        a.deltam = nothing
         return a 
     end
     function EnsObs(ens::EnsInfo, mu::Vector{Vector{Float64}}, m_ps::Vector{uwreal}, m_vec::Vector{uwreal})
@@ -163,3 +165,30 @@ mutable struct MatInfo
     end
 end
 
+mutable struct Cat
+    obs::Vector{uwreal}
+    phih::Vector{uwreal}
+    phih_ph::uwreal
+    #models::Vector{Vector{Function}}
+    x_tofit::Array{uwreal}
+    info::String
+    fit_param::Vector{Vector{uwreal}} #store all fit parameters
+    fit_res::Vector{uwreal} # store all fit results
+    chi2_corr::Vector{Float64} # store all chi2 corrected
+    n_param::Vector{Int64} # store the number of parameters. Maybe this is redundant
+    aic::Vector{Float64} # store the AIC value 
+    function Cat(_obs::Vector{uwreal}, _xtofit::Array{uwreal}, _phi_ph::uwreal, _info::String)
+        a = new()
+        a.obs = _obs
+        #a.models = _models
+        a.info = _info
+        a.x_tofit = _xtofit
+        a.phih_ph = _phi_ph
+        a.fit_param = Vector{Vector{uwreal}}(undef, 0)
+        a.fit_res = Vector{uwreal}(undef, 0)
+        a.chi2_corr = Vector{Float64}(undef, 0)
+        a.n_param = Vector{Int64}(undef, 0)
+        a.aic = Vector{Float64}(undef, 0)
+        return a 
+    end
+end
\ No newline at end of file