Commit d1cf0981 authored by Alessandro 's avatar Alessandro

changes for ift laptop, decay constant from gevp, renormalization constant za...

changes for ift laptop, decay constant from gevp, renormalization constant za added, exploring multiple continuum+chiral extrapolations
parent 515a5d52
#========= 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]
"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]
)
ttl_obs_db = Dict(
"muc"=> ["mu_c"],
......@@ -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
......@@ -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
......@@ -435,7 +436,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 +445,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 +488,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 +498,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 +524,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_", ttl_obs[k], ".pdf")
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_a=0", ttl_obs[k], ".pdf")
savefig(joinpath(path_plot, t))
#savefig(joinpath("/Users/ale/Desktop/plottest", t))
close()
......@@ -541,6 +576,67 @@ 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)) .* (2/3 .* getfield.(ensobs, :m_lh_match) .+ 2/3 .* getfield.(ensobs, :m_sh_match))
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[5] * x[:, 1].^2 ) + (p[3]) * x[:, 2] + p[4] * x[:,2].^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, 5)
test_charm = par[1] + par[3] * phi2_ph + par[4] * phi2_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 + par[4] * phi2_ph^2
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.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
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
......
......@@ -173,10 +173,19 @@ 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]
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)
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
......@@ -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
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment