Merge branch 'master' of gitlab.ift.uam-csic.es:Conigli/gevp-automated-analysis

src/const.jl

 #========= ENSEMBLE DATABASE ========#
 ens_db = Dict(
     #"ens_id"=>[L, beta, is_deg?, m_pi]
+    "H102r002" => [32, 3.4, false, 0.17979],
+    "H102r001" => [32, 3.4, false, 0.17979],
+    "H101"     => [32, 3.4, true, 0.18302],
     "H400"     => [32, 3.46, true, 0.16345],
     "N200"     => [48, 3.55, false, 0.09222],
+    "N202"     => [48, 3.55, true, 0.13407],
     "N203"     => [48, 3.55, false, 0.11224],
     "N300"     => [48, 3.70, true, 0.10630],
     "J303"     => [64, 3.70, false, 0.06514]
@@ -42,16 +46,21 @@ 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
+# 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
 #Covariance matrices
 const C1 = zeros(5, 5)
 const C2 = zeros(5, 5)
 const C3 = zeros(4, 4)
 const C4 = zeros(6,6)
+const C5 = zeros(5, 5)
 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
+        C5[i, i] = ZA_err[i]^2
         if i<=4
             C3[i,i] = t0_error[i] ^ 2
         end
@@ -64,8 +73,10 @@ 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_)
+const Za = cobs(ZA_data, C5, "ZA")
 
 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]
+za(beta::Float64) = Za[b_values .== beta][1]
\ No newline at end of file

src/gevp.jl

@@ -2,10 +2,11 @@ using Revise, juobs, BDIO, DelimitedFiles, ADerrors, PyPlot, LaTeXStrings, LsqFi
 
 
 #============= SET UP VARIABLES ===========#
-const path_data = "/Users/ale/Desktop/data"
-const path_plat = "/Users/ale/automation/plat.txt"
-const path_results = "/Users/ale/Desktop/results_gevp"
-const ensembles = [ "J303","N300", "H400"]
+const path_data = "/media/alessandro/4277fef2-edc5-4e0d-89cb-f5d1d44fbc8c/data"
+const path_plat = "/home/alessandro/automated-analysis/plat.txt"
+const path_results = "/home/alessandro/Desktop/results_gevp"
+#some problem with the pion mass in N203
+const ensembles = ["J303", "H400", "N200" , "N202", "N300", "H101", "N203"]
 const sector = Dict("ll"=>false, "ls"=>false, "lh"=>true, "ss"=>false, "sh"=>true, "hh"=>true )
 const tau = 3
 const _t0 = 2
@@ -112,11 +113,14 @@ println("\n Matching Procedure \n")
             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
+        #target = ens.a * (2/3*M[1] + 1/3* M[3])/ hc
+        target = ens.a * M[5] / hc
         if !ens.ensinfo.deg
-            ens.muh_target = match_muc(muh, ens.m_lh, ens.m_sh, target)
+            #ens.muh_target = match_muc(muh, ens.m_lh, ens.m_sh, target)
+            ens.muh_target = match_muc_heavymass(muh, ens.m_hh, target)
         else
-            ens.muh_target = match_muc(muh, ens.m_lh, ens.m_sh, target)
+            #ens.muh_target = match_muc(muh, ens.m_lh, ens.m_sh, target)
+            ens.muh_target = match_muc_heavymass(muh, ens.m_hh, target)
         end
         uwerr(ens.muh_target)
         println(t0, "\n", a, "\n", ens.muh_target)
@@ -435,7 +439,7 @@ if sector["hh"]
         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))        
+        f_jpsi .= za.(getfield.(ensinfo,:beta)) .* getfield.(ensobs, :f_hh_vec_match) .* sqrt.(8 * getfield.(ensobs,:t0))        
         uwerr.(m_etac)
         uwerr.(m_jpsi)
         uwerr.(f_etac)
@@ -444,7 +448,7 @@ catch
         error("Check ensobs[i].mu_list to see whether all the ensembles have a heavy-heavy sector.")
     end
 end
-
+##
 
 #===== CONTINUUM AND CHIRAL EXTRAPOLATION =======#
 
@@ -487,6 +491,7 @@ for lab in label
         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
@@ -496,8 +501,9 @@ uwerr(phi2_ph)
 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
+@. cl_chir_model(x, p) = (p[1] + p[2] * x[:, 1]) + (p[3] + p[4] * x[:,1]) * x[:, 2] #linear fits
 
+## FIT PROJECTING TO PHI2=0
 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)
@@ -521,11 +527,43 @@ for k = 1:length(obs)
     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)
+    xlim(-0.002,0.05)
+    #ylim(2.8,3.4)
+    legend(ncol=2)
+    display(gcf())
+    t = string("fit_phi2=0", ttl_obs[k], ".pdf")
+    savefig(joinpath(path_plot, t))
+    #savefig(joinpath("/Users/ale/Desktop/plottest", t))
+    close()
+
+end
+
+## FIT PROJECTING TO a=0
+xarr = Float64.(range(0.0, stop=0.9, length=100))
+for k = 1:length(obs)
+    par = fit_routine(cl_chir_model, value.(x), obs[k], 4)   
+    y = Vector{uwreal}(undef, 100)
+    for j=1:100
+        y[j] = par[1]  + par[3] * xarr[j]# phi2_ph
+    end
+    uwerr.(y)
+    figure()
+    for b in unique(getfield.(ensinfo,:beta)) #select point with same beta
+        nn = findall(x-> x == b, getfield.(ensinfo,: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
+    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"$\Phi_2$")
+    ylabel(ylbl[k])
+    xlim(-0.002,0.9)
     #ylim(2.8,3.4)
     legend(ncol=2)
     display(gcf())
-    t = string("fit_", ttl_obs[k], ".pdf")
+    t = string("fit_a=0", ttl_obs[k], ".pdf")
     savefig(joinpath(path_plot, t))
     #savefig(joinpath("/Users/ale/Desktop/plottest", t))
     close()
@@ -541,6 +579,71 @@ for k = 1:length(obs)
     println(ttl_obs[k], "(MeV) = ", obs_ph[k])
 end
 
+
+## TRYING TO ESTIMATE CHARM MASS WITH DIFFERENT MODELS TO ADDRESS SYSTEMATICS
+
+# 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
+phih = Vector{uwreal}(undef, length(ensembles))
+phih = sqrt.( 8 .* getfield.(ensobs, :t0)) .*getfield.(ensobs, :m_hh_match)  
+phih_ph =  M_values[5] * t0_ph[1] / hc
+#phih = sqrt.( 8 .* getfield.(ensobs, :t0)) .*(2/3 .* getfield.(ensobs, :m_lh_match) .+ 1/3 .* getfield.(ensobs, :m_sh_match))
+#phih_ph =  (2/3 * M_values[1] + 1/3 * M_values[3]) * t0_ph[1] / hc
+
+phi2_ph = (t0_ph[1] * 139.57039 / hc)^2
+uwerr(phi2_ph)
+#fit routine 
+x = [1 ./(8 .* getfield.(ensobs,:t0)) phi2 phih] #x[1] = a^2 / (8t0), x[2] = 8t0 mpi^2, x[3] = sqrt(8t0) m_flav_av
+uwerr.(x)
+#modify CL+chiral fit  model
+@. cl_chir_model(x, p) = (p[1] + p[2] * x[:, 1] ) + (p[3]) * x[:, 2] + p[4] * x[:,1] * x[:,3].^2 #linear fits
+
+test_charm = uwreal(0.0)
+xarr = Float64.(range(0.0, stop=0.05, length=100))
+for k = [1:1]
+    par = fit_routine(cl_chir_model, value.(x), muc, 4)
+    test_charm = par[1] + par[3] * phi2_ph # + par[4] * phih_ph^2
+    uwerr(test_charm)    
+    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.(muc[nn]), err.(muc[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(test_charm), err(test_charm), 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.05)
+    ylim(2.4,3.4)
+    legend(ncol=2)
+    display(gcf())
+    t = string("fit_phi2=0", ttl_obs[k], ".pdf")
+    #savefig(joinpath(path_plot, t))
+    #savefig(joinpath("/Users/ale/Desktop/plottest", t))
+    close()
+
+end
+
+test_charm_phys = test_charm * hc / t0_ph[1]
+uwerr(test_charm_phys)
+println("\n Results in the continuum limit \n \n")
+println( ylbl[1], " = ", test_charm, "\n")
+#phys
+println(ttl_obs[1], "(MeV) = ", test_charm_phys, "\n\n")
+    
+
+
+
+##
 #======== STORE RESULTS ========#
 
 open(joinpath(path_results, "results.txt"), "w") do f 

src/tools.jl

@@ -173,11 +173,20 @@ 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)
             else
                 plat_dec0[i] = -vec_dec(elem, plat, mass[i], pl=pl)
             end
+        catch
+            println("Decay constant for the ensemble ", mat_obs[i].ensinfo.id, " failed in the sector ", mu_list[i], " pseudo sector?", pseudo)
+            println("The associated effective mass is ", mass[i])
+            plat_dec0[i] = uwreal(0)
+        end
     end
     return plat_dec0
 end
@@ -227,7 +236,7 @@ function vec_dec(mat_elem::Vector{uwreal}, plat::Array{Int64}, mass::uwreal; pl:
         fill_between(plat[1]:plat[2], v-e, v+e, color="green")    
         display(gcf())
     end
-    return aux / mass
+    return sqrt(2 / mass) * aux #aux / mass
 end
 function param_av(fit_masses::Vector{uwreal})
     w = 1 ./ err.(fit_masses).^2
@@ -309,6 +318,11 @@ function match_muc(muh, m_lh, m_sh, target)
     muh_target = x_lin_fit(par, target)
     return muh_target
 end
+function match_muc_heavymass(muh, m_hh, target)
+    par, chi2exp = lin_fit(muh, m_hh)
+    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)