Skip to content

J
juobs
Commit 7f62176c authored by Javier's avatar Javier
Browse files
Options

first analysis added 

Computes matching in the charm sector 
Scaling + chiral extrapolations are missing
parent 5d427866
Hide whitespace changes
Inline Side-by-side
Showing with 397 additions and 0 deletions
analysis/analysis.jl 0 → 100644
View file @ 7f62176c
#TODO include rw, different plateaux depending on obs, print chi2, scaling+chiral plots
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")
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
analysis/functions.jl 0 → 100644
View file @ 7f62176c
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)
aux = read_mesons.(f, 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::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 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_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
analysis/plat.txt 0 → 100644
View file @ 7f62176c
#H400
ll 40 83
lh 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 12 140
ls 123 140
lh 123 140
ss 123 140
sh 123 140
hh 123 155
#end
analysis/results.txt 0 → 100644
View file @ 7f62176c
(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
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