Skip to content

G
GEVP automated analysis
Commit 838f8293 authored by Alessandro 's avatar Alessandro
Browse files
Options

major update on the structure of the analysis code

parent c9fca371
Hide whitespace changes
Inline Side-by-side
Showing with 2294 additions and 10 deletions
src/const.jl
View file @ 838f8293
......@@ -31,6 +31,22 @@ ylbl_db = Dict(
"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}$"]
)
muh_val =Dict(
#"ens_id"=>[L, beta, is_deg?, m_pi]
"H105" => [],
"H102r001" => [],
"H102r002" => [ 0.2403, 0.228285, 0.252315],
"H101" => [ 0.2505, 0.237975, 0.263025],
"H400" => [ 0.2149, 0.204155, 0.225645],
"N200" => [ 0.18190, 0.172805, 0.190995],
"N202" => [ 0.17590, 0.167105, 0.184695],
"N203" => [ 0.1825, 0.173375, 0.191625],
"N300" => [0.1378, 0.13091, 0.14469],
"J303" => [0.14000, 0.13300, 0.14700]
)
#=
#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)
......@@ -80,4 +96,5 @@ 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
za(beta::Float64) = Za[b_values .== beta][1]
=#
\ No newline at end of file
src/gevp_continuum.jl 0 → 100644
View file @ 838f8293
#################################
# This file was created by Alessandro Conigli - 2022
#
# Perform chiral- continuum extrapolations and AIC analysis
#
#
#
#
#
##################################
using Revise, ADerrors, BDIO, juobs, PyPlot, LaTeXStrings, DelimitedFiles, Combinatorics
using LsqFit, juobs.LeastSquaresOptim
include("gevp_types.jl")
include("gevp_reader.jl")
include("gevp_tools.jl")
include("const.jl")
include("gevp_func_comb.jl")
#============= INPUT PAR ===========#
# desktop parameters
#include("/home/alessandro/juobs/constants/juobs_const.jl")
#path_data ="/home/alessandro/Desktop/data/heavy_3pt/3pt"
#path_rw = "/media/alessandro/4277fef2-edc5-4e0d-89cb-f5d1d44fbc8c/data/rwf"
#path_plot = "/home/alessandro/google-drive/phd/secondment/3pf test/analysis/plots"
# macbook parameters
include("/Users/ale/juobs/constants/juobs_const.jl")
path_plot = "/Users/ale/Il mio Drive/phd/analysis/charm_gevp_2022/plots"
path_results = "/Users/ale/Il mio Drive/phd/analysis/charm_gevp_2022/bdio"
path_t0_and_dm = "/Users/ale/Desktop/data"
# latex-style labels
PyPlot.rc("font", family="sans serif", size=12)
rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["text.usetex"] = false
rcParams["mathtext.fontset"] = "cm"
#rcParams["font.size"] =11
rcParams["axes.labelsize"] =14
const ens_list = ["H101", "H102r002", "H400", "N202", "N203", "N200", "J303", "N300"]
const sector = Dict("pseudo"=>true, "vector"=>true)
# loading ensemble information
ensinfo = Vector{EnsInfo}(undef, length(ens_list))
for i in 1:length(ens_list)
ens = ens_list[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
##########################
## initialize containers
##########################
#charm quark mass
muc_t0 = Vector{Vector{uwreal}}(undef, length(ensinfo))
# pseudoscalar masses
if sector["pseudo"]
mps_ll_t0 = Vector{uwreal}(undef, length(ensinfo)) # pion mass
mps_ls_t0 = Vector{uwreal}(undef, length(ensinfo)) # kaon mass
mps_ss_t0 = Vector{uwreal}(undef, length(ensinfo)) # ss mass
mps_lh_t0 = Vector{Vector{uwreal}}(undef, length(ensinfo)) # D mass
mps_sh_t0 = Vector{Vector{uwreal}}(undef, length(ensinfo)) # Ds mass
mps_hh_t0 = Vector{Vector{uwreal}}(undef, length(ensinfo)) # η_c mass
# pseudoscalar decays
fps_ll_t0 = Vector{uwreal}(undef, length(ensinfo)) # pion decay
fps_ls_t0 = Vector{uwreal}(undef, length(ensinfo)) # kaon decay
fps_ss_t0 = Vector{uwreal}(undef, length(ensinfo)) # ss decay
fps_lh_t0 = Vector{Vector{uwreal}}(undef, length(ensinfo)) # D decay
fps_sh_t0 = Vector{Vector{uwreal}}(undef, length(ensinfo)) # Ds decay
fps_hh_t0 = Vector{Vector{uwreal}}(undef, length(ensinfo)) # η_c decay
end
if sector["vector"]
# vector masses
mvec_ll_t0 = Vector{uwreal}(undef, length(ensinfo)) # rho mass
mvec_ls_t0 = Vector{uwreal}(undef, length(ensinfo)) # ls vec mass
mvec_ss_t0 = Vector{uwreal}(undef, length(ensinfo)) # ss vec mass
mvec_lh_t0 = Vector{Vector{uwreal}}(undef, length(ensinfo)) # D* mass
mvec_sh_t0 = Vector{Vector{uwreal}}(undef, length(ensinfo)) # Ds* mass
mvec_hh_t0 = Vector{Vector{uwreal}}(undef, length(ensinfo)) # J/Ψ mass
# vector decays
fvec_ll_t0 = Vector{uwreal}(undef, length(ensinfo)) # rho decay
fvec_ls_t0 = Vector{uwreal}(undef, length(ensinfo)) # ls vec decay
fvec_ss_t0 = Vector{uwreal}(undef, length(ensinfo)) # ss vec decay
fvec_lh_t0 = Vector{Vector{uwreal}}(undef, length(ensinfo)) # D* decay
fvec_sh_t0 = Vector{Vector{uwreal}}(undef, length(ensinfo)) # Ds* decay
fvec_hh_t0 = Vector{Vector{uwreal}}(undef, length(ensinfo)) # J/Ψ decay
end
#######################
## read bdio files
#######################
# read t0_shifted from BDIO
dm = read_BDIO(path_t0_and_dm, "delta")
t0_shifted = read_BDIO(path_t0_and_dm, "t0")
uwerr.(t0_shifted)
label = vcat(ensembles.(t0_shifted)...)
#
N_mh = fill(3, length(ensinfo))
use_t0_from_bdio = true
for i in 1:length(ensinfo)
ens = ensinfo[i]
idx_t0 = findall(x-> x==ens.id, label)[1]
t0_val = use_t0_from_bdio ? t0_shifted[idx_t0] : t0(ens.beta)
println(t0_val)
println("Reading observables for ensemble $(ens.id)" )
aux_path = joinpath(path_results, "charm_mass")
muc_t0[i] = muh_val[ens.id ] .* [sqrt(8 * t0_val) / ZP(ens.beta) for k in 1:N_mh[i]] #read_BDIO(aux_path, ens.id)
# implement t0 from bdio and not only from database
if sector["pseudo"]
# read ps masses and decays
aux_m = read_BDIO(joinpath(path_results, "ps_masses"), ens.id)
aux_f = read_BDIO(joinpath(path_results, "ps_decays"), ens.id)
mps_ll_t0[i] = sqrt(8 * t0_val) * aux_m[1]
fps_ll_t0[i] = sqrt(8 * t0_val) * aux_f[1]
if ens.deg
mps_ss_t0[i] = mps_ls_t0[i] = mps_ll_t0[i]
mps_lh_t0[i] = mps_sh_t0[i] = aux_m[2:4] .* [sqrt(8*t0_val) for k in 1:N_mh[i]]
mps_hh_t0[i] = aux_m[5:7] .* [sqrt(8*t0_val) for k in 1:N_mh[i]]
fps_ss_t0[i] = fps_ls_t0[i] = fps_ll_t0[i]
fps_lh_t0[i] = fps_sh_t0[i] = aux_f[2:4] .* [sqrt(8*t0_val) for k in 1:N_mh[i]]
fps_hh_t0[i] = aux_f[5:7] .* [sqrt(8*t0_val) for k in 1:N_mh[i]]
else
mps_ls_t0[i] = aux_m[2] * sqrt(8*t0_val)
mps_ss_t0[i] = aux_m[3] * sqrt(8*t0_val)
mps_lh_t0[i] = aux_m[4:6] .* [sqrt(8*t0_val) for k in 1:N_mh[i]]
mps_sh_t0[i] = aux_m[7:9] .* [sqrt(8*t0_val) for k in 1:N_mh[i]]
mps_hh_t0[i] = aux_m[10:12] .* [sqrt(8*t0_val) for k in 1:N_mh[i]]
fps_ls_t0[i] = aux_f[2] * sqrt(8*t0_val)
fps_ss_t0[i] = aux_f[3] * sqrt(8*t0_val)
fps_lh_t0[i] = aux_f[4:6] .* [sqrt(8*t0_val) for k in 1:N_mh[i]]
fps_sh_t0[i] = aux_f[7:9] .* [sqrt(8*t0_val) for k in 1:N_mh[i]]
fps_hh_t0[i] = aux_f[10:12] .* [sqrt(8*t0_val) for k in 1:N_mh[i]]
end
end
if sector["vector"]
# read ps masses and decays
aux_m = read_BDIO(joinpath(path_results, "vec_masses"), ens.id)
aux_f = read_BDIO(joinpath(path_results, "vec_decays"), ens.id)
mvec_ll_t0[i] = sqrt(8*t0_val) * aux_m[1]
fvec_ll_t0[i] = sqrt(8*t0_val) * aux_f[1]
if ens.deg
mvec_ss_t0[i] = mvec_ls_t0[i] = mvec_ll_t0[i]
mvec_lh_t0[i] = mvec_sh_t0[i] = aux_m[2:4] .* [sqrt(8*t0_val) for k in 1:N_mh[i]]
mvec_hh_t0[i] = aux_m[5:7] .* [sqrt(8*t0_val) for k in 1:N_mh[i]]
fvec_ss_t0[i] = fvec_ls_t0[i] = fvec_ll_t0[i]
fvec_lh_t0[i] = fvec_sh_t0[i] = aux_f[2:4] .* [sqrt(8*t0_val) for k in 1:N_mh[i]]
fvec_hh_t0[i] = aux_f[5:7] .* [sqrt(8*t0_val) for k in 1:N_mh[i]]
else
mvec_ls_t0[i] = aux_m[2] * sqrt(8*t0_val)
mvec_ss_t0[i] = aux_m[3] * sqrt(8*t0_val)
mvec_lh_t0[i] = aux_m[4:6] .* [sqrt(8*t0_val) for k in 1:N_mh[i]]
mvec_sh_t0[i] = aux_m[7:9] .* [sqrt(8*t0_val) for k in 1:N_mh[i]]
mvec_hh_t0[i] = aux_m[10:12] .* [sqrt(8*t0_val) for k in 1:N_mh[i]]
fvec_ls_t0[i] = aux_f[2] * sqrt(8*t0_val)
fvec_ss_t0[i] = aux_f[3] * sqrt(8*t0_val)
fvec_lh_t0[i] = aux_f[4:6] .* [sqrt(8*t0_val) for k in 1:N_mh[i]]
fvec_sh_t0[i] = aux_f[7:9] .* [sqrt(8*t0_val) for k in 1:N_mh[i]]
fvec_hh_t0[i] = aux_f[10:12] .* [sqrt(8*t0_val) for k in 1:N_mh[i]]
end
end
end
################################################################
## Collect physical values quantities for matching and fitting
################################################################
sqrt8t0 = Vector{uwreal}(undef, 0)
#fb = BDIO_open("/Users/ale/Desktop/data/t0_shifted.bdio")
for i in 1:length(ensinfo)
ens = ensinfo[i]
idx_t0 = findall(x-> x==ens.id, label)[1]
t0_val = use_t0_from_bdio ? t0_shifted[idx_t0] : t0(ens.beta)
[push!(sqrt8t0, sqrt(8 *t0_val)) for k in 1:N_mh[i]]
end
# lattice units
phih_eta = vcat(mps_hh_t0...)
phih_flav = 2/3 .* vcat(mps_lh_t0...) + 1/3 .* vcat(mps_sh_t0...)
phih_spav = 1/4 .* (phih_flav .+ 2 .* vcat(mvec_lh_t0...) .+ vcat(mvec_sh_t0...))
phi2 = vcat([fill(mps_ll_t0[i]^2, N_mh[i]) for i in 1:length(ensinfo)]...) #for i in 1:length(ensinfo)]
phi4 = vcat([fill(1/2 * mps_ll_t0[i]^2 + mps_ls_t0[i]^2, N_mh[i]) for i in 1:length(ensinfo)]...)# 1/2 .* mps_ll_t0.^2 .+ mps_ls_t0 .^2
# physical units
phih_eta_ph = M_values[5] * t0_ph[1] / hc
phih_flav_ph = (2/3 * M_values[1] + 1/3 * M_values[3]) * t0_ph[1] / hc
phih_spav_ph = 1/4 * (2/3 * M_values[1] + 1/3 * M_values[3] + 2* M_values[2] + M_values[4]) * t0_ph[1] / hc
phi2_ph = (t0_ph[1] * 134.8 / hc)^2 # try 139 or 134.8 isospin breaking pion mass
phi4_ph = uwreal([1.098, 0.0010], "phi4_phys")
# xdata array to feed the fit routine for different categories
xdata_flav = [ 1 ./ sqrt8t0.^2 phi2 phih_flav]
xdata_eta = [ 1 ./ sqrt8t0.^2 phi2 phih_eta]
xdata_spav = [ 1 ./ sqrt8t0.^2 phi2 phih_spav]
########################
## Creating categories
########################
# WPM WINDOWS
wpmm = Dict{String, Vector{Float64}}()
wpmm["N202"] = [-1.0, 2.0,-1.0, -1.0]
#wpmm["N200"] = [-1.0, 2.0,-1.0, -1.0]
#wpmm[300] = [-1.0, 1.5,-1.0, -1.0]
cat_muc_eta = Cat(vcat(muc_t0...), xdata_eta, phih_eta_ph, "muc_etac_matching")
cat_muc_flav = Cat(vcat(muc_t0...), xdata_flav, phih_flav_ph, "muc_flav_matching")
cat_muc_spav = Cat(vcat(muc_t0...), xdata_spav, phih_spav_ph, "muc_spav_matching")
cat_muc_tot = [cat_muc_eta, cat_muc_flav, cat_muc_spav]
for cat in cat_muc_tot[1:2]
for i in 1:length(model_charm)
#try
par, chi2_over_chiexp = fit_routine(model_charm[i], value.(cat.x_tofit), cat.obs, n_param_charm[i], wpm=wpmm)
aux = Mrat * (par[1] + par[2]*phi2_ph + par[3]*cat.phih_ph)
uwerr(aux)
println("charm mass = ", aux)
push!(cat.fit_res, aux)
push!(cat.fit_param, par)
push!(cat.chi2_corr, chi2_over_chiexp)
push!(cat.n_param, n_param_charm[i])
push!(cat.dof, length(cat.obs)-n_param_charm[i] )
push!(cat.aic, cat.chi2_corr[i]*cat.dof[i] +2*n_param_charm[i])
println("aic = ", cat.aic[i])
println("\n")
#catch
# println("something went wrong")
#end
end
end
muc_av_eta = aic_model_ave(cat_muc_eta)
muc_av_flav = aic_model_ave(cat_muc_flav)
mc_fin, mc_syst = category_av([cat_muc_eta, cat_muc_flav])
muc_av_spav = aic_model_ave(cat_muc_spav)
plot_cat(cat_muc_flav, L"$M_c^{\mathrm{RGI}}(N_f=3)\mathrm{MeV}$" , save=true)
plot_cont_lim_phi2(cat_muc_eta, model_charm, L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", dec=false, save=false)
plot_cont_lim(cat_muc_eta, model_charm, L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", save=false)
plot_cl_charm(cat_muc_eta, model_charm, L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo)
plot_chi_charm(cat_muc_flav, model_charm, L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo)
plot_all_cl_charm(cat_muc_tot[1:2], model_charm,L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo, n=[0,0,12] )# , p=path_plot)
#########################
## fits for MD and MDS
#########################
cat_mD_eta = Cat(vcat(mps_lh_t0...), xdata_eta, phih_eta_ph, "mD_etac_matching")
cat_mD_flav = Cat(vcat(mps_lh_t0...), xdata_flav, phih_flav_ph, "mD_flav_matching")
cat_mDs_eta = Cat(vcat(mps_sh_t0...), xdata_eta, phih_eta_ph, "mDs_etac_matching")
cat_mDs_flav = Cat(vcat(mps_sh_t0...), xdata_flav, phih_flav_ph, "mDs_flav_matching")
cat_D_Ds_masses =[cat_mD_eta, cat_mD_flav, cat_mDs_eta, cat_mDs_flav]
for cat in cat_D_Ds_masses[1:end]
println(cat.info)
for i in 1:length(model_charm)
#try
par, chi2_over_chi2exp = fit_routine(model_charm[i], value.(cat.x_tofit), cat.obs, n_param_charm[i], wpm=wpmm)
aux = par[1] + par[2]*phi2_ph + par[3]*(cat.phih_ph)
uwerr(aux)
println("mD mass = ", aux)
push!(cat.fit_res, aux)
push!(cat.fit_param, par)
push!(cat.chi2_corr, chi2_over_chi2exp)
push!(cat.n_param, n_param_charm[i])
push!(cat.dof, length(cat.obs)-n_param_charm[i] )
push!(cat.aic, chi2_over_chi2exp*cat.dof[i] +2*n_param_charm[i])
println("chi2 / chiexp = ", chi2_over_chi2exp)
println("\n")
#catch
# println("something went wrong")
#end
end
end
mD_av = aic_model_ave(cat_mD_flav)
## plots
plot_cont_lim_phi2_meson_mass(cat_mDs_eta, model_charm, L"$\sqrt{8t_0}M_{D_s}$", dec=false, save=false)
plot_chi_masses(cat_mDs_flav, model_charm, L"$\sqrt{8t_0}M_{D_s}$", phi2_ph, ensinfo, n=1 , p=path_plot)
plot_cont_lim(cat_muc_eta, model_charm, L"$\sqrt{8t_0}M_{D_s}$", save=false)
plot_cat(cat_muc_flav, L"$\sqrt{8t_0}M_{D}$" )
################
## fit for Fd
################
xdata_eta_dec[:] = xdata_eta
xdata_eta_dec[:,3] = 1.0 ./sqrt.(xdata_eta[:,3])
xdata_flav_dec = xdata_flav
xdata_flav_dec[:,3] = 1.0 ./sqrt.(xdata_flav[:,3])
cat_fD_eta = Cat(vcat(fps_lh_t0...), xdata_eta, phih_eta_ph, "fD_etac_matching")
cat_fD_flav = Cat(vcat(fps_lh_t0...), xdata_flav, phih_flav_ph, "fD_flav_matching")
cat_fDs_eta = Cat(vcat(fps_sh_t0...), xdata_eta, phih_eta_ph, "fDs_etac_matching")
cat_fDs_flav = Cat(vcat(fps_sh_t0...), xdata_flav, phih_flav_ph, "fDs_flav_matching")
cat_fD_fDs = [cat_fD_eta, cat_fD_flav, cat_fDs_eta, cat_fDs_flav]
for cat in cat_fD_fDs[1:4]
println(cat.info)
for i in 1:10#length(model_decays)
try
par, chi2_over_chi2exp = fit_routine(model_decays[i], cat.x_tofit, cat.obs, n_param_decays[i], covar=false, wpm=wpmm)
aux = par[1] + par[2]*phi2_ph + par[3]/sqrt((cat.phih_ph))
uwerr(aux)
println("fD(s) decays = ", aux)
push!(cat.fit_res, aux)
push!(cat.fit_param, par)
push!(cat.chi2_corr, chi2_over_chi2exp)
push!(cat.n_param, n_param_decays[i])
cat.dof = length(cat.obs)-n_param_decays[i]
push!(cat.aic, chi2_over_chi2exp*cat.dof +4*n_param_decays[i])
println("chi2 / chiexp = ", chi2_over_chi2exp)
println("aic = ", cat.aic[i])
println("\n")
catch
println("something went wrong")
end
end
end
fD_av_eta = aic_model_ave(cat_fD_eta)
fD_av_flav = aic_model_ave(cat_fD_flav)
fDs_av_eta = aic_model_ave(cat_fDs_eta)
fDs_av_flav = aic_model_ave(cat_fDs_flav)
plot_cont_lim_phi2(cat_fD_eta, model_decays, L"$\sqrt{8t_0}f_{D_s}$", dec=true, save=false)
plot_cat(cat_fD_flav, L"ciao")
## test fitting without multiplying the data for N_mh
xdat_test_fl = xdata_flav[:,:]
xdat_test_fl[:,1] = xdata_flav[1:3:end,1]
## testing fit charm
par_eta, chi = fit_routine(f[1][1], xdata_eta, vcat(muc_t0...), 5 , covar=true)
par_flav, chi = fit_routine(f[1][1], xdata_flav, vcat(muc_t0...), 5 , covar=true)
val_ph_eta = par_eta[1] + par_eta[2]*phi2_ph + par_eta[3]*phih_eta_ph
val_ph_flav = par_flav[1] + par_flav[2]*phi2_ph + par_flav[3]*phih_flav_ph
uwerr(val_ph_eta)
uwerr(val_ph_flav)
##
muh_target_arr = []
fig = figure()
count = 1
for i in 1:length(ensinfo)
xx = 1 / (8*t0(ensinfo[i].beta))
muh_target_eta = basemodel([xx unique(phi2)[i] phih_eta_ph], value.(par_eta))[1]
muh_target_flav = basemodel([xx unique(phi2)[i] phih_flav_ph], value.(par_flav))[1]
#muh_target = match_muc(value.(muc_t0[i]), mps_lh_t0[i], mps_sh_t0[i], phih_flav_ph)
uwerr(muh_target_eta)
uwerr(muh_target_flav)
errorbar(value(xx), value(muh_target_eta), err(muh_target_eta), fmt="s" , capsize=2, mfc="none", color="red", ecolor="red")#, label=L"$m_{η_c}$" )
errorbar(value(xx)+0.0005, value(muh_target_flav), err(muh_target_flav), fmt="s", capsize=2, mfc="none", color="royalblue", ecolor="royalblue")#,label=L"$m_{\bar D}$")
end
xarr = Float64.(range(0.0, stop=0.4, length=100))
y_eta = Vector{uwreal}(undef, 100)
y_flav = Vector{uwreal}(undef, 100)
for i in 1:100
y_eta[i] = par_eta[1] + par_eta[2]*phi2_ph + par_eta[3]*phih_eta_ph + par_eta[4]*xarr[i]
y_flav[i]= par_flav[1] + par_flav[2]*phi2_ph + par_flav[3]*phih_flav_ph + par_flav[4]*xarr[i]
end
uwerr.(y_eta)
uwerr.(y_flav)
fill_between(xarr, value.(y_eta) .+ err.(y_eta),value.(y_eta) .- err.(y_eta) , alpha=0.2,zorder=1, lw=0.8, color="red")
fill_between(xarr, value.(y_flav) .+ err.(y_flav),value.(y_flav) .- err.(y_flav) , alpha=0.2,zorder=1, lw=0.8, color="royalblue")
errorbar(0, value(val_ph_eta), err(val_ph_eta), fmt="s", ms=8, capsize=2, color="red",ecolor="red")
errorbar(0.0005, value(val_ph_flav), err(val_ph_flav), fmt="s", capsize=2, ms=8, color="royalblue",ecolor="royalblue")
xlim(-0.001, 0.055)
ylim(2.9, 3.3)
xlabel(L"$a^2/8t_0$")
ylabel(L"$\sqrt{8t_0}Z_M^{tm}\mu_c$")
legend(L"$m_{η_c}$", L"$m_{\bar D}$")
display(fig)
src/gevp_func_comb.jl 0 → 100644
View file @ 838f8293
# FUNCTIONS
#=
basemodel(x,p) = p[1] .+ p[2] .* x[:,2] .+ p[3] .* x[:,3] .+ p[4] .* x[:,1]
a2phi2(x) = x[:,1] .* x[:,2] # phi2*a^2/8t0
a2phih2(x) = x[:,1] .* x[:,3].^2 # phih^2*a^2/8t0
a4phih4(x) = x[:,1].^(2) .* x[:,3].^4 #phih^4*(a^2/8t0)^2
a4phih4_2(x) = x[:,1].^(2) .* x[:,3].^2 #(a^2/8t0)^2*phih^2
a4cutoff(x) = x[:,1].^(2) #(a^2/8t0)^2
func_list = [a2phi2, a2phih2, a4phih4, a4phih4_2, a4cutoff]
# COMBINATIONS (LINEAR)
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]
npar = length(func_list) # number of extra parameters
n_param = Vector(4:npar+4)
n_param_charm = Vector{Int64}(undef,0)
push!(n_param_charm, 4)
#f_lin definition
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)
#f_aux definition -> similar to f_lin but without vectorized opeartors -> useful for plots
f_aux = Vector{Vector{Function}}(undef, npar+1)
f_aux[1] = Vector{Function}(undef, 1)
f_aux[1][1] = (x,p) -> p[1] + p[2] * x[:,2] + p[3] * x[:,3] + p[4] * x[:,1]
#loop over combinations
for n = 2:npar+1
aux = filter(x->sum(x) == n-1, func_map)
println(n_param[n])
f_lin[n] = Vector{Function}(undef, length(aux))
f_aux[n] = Vector{Function}(undef, length(aux))
k = 1
for a in aux
push!(n_param_charm, n_param[n])
#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]))
f_aux[n][k] = (x,p) -> f_aux[1][1](x, p) + sum([p[i+4] for i=1:(n-1)] .* (fill(x,n-1) .|> func_list[a]))
k = k + 1
end
end
# COMBINATIONS (NONLINEAR)
f_non_lin = Vector{Vector{Function}}(undef, npar) #f[i][j]-> i: number of extra parameters, j: combinations
f_non_lin_aux = Vector{Vector{Function}}(undef, npar)
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))
f_non_lin_aux[n] = Vector{Function}(undef, length(aux))
k = 1
for a in aux
#vectorized function
push!(n_param_charm, n_param[n]+1)
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])))
f_non_lin_aux[n][k] = (x,p) -> p[1] + p[2] * x[:,2] + p[3] * x[:,3] + p[4] * x[:,1] +
(p[1] + p[2] * x[:,2] + p[3] * x[:,3] + p[4] * x[:,1]) .* sum([p[i+4] for i=1:(n)] .* (fill(x,n) .|> func_list[a]))
k = k + 1
end
end
# APPEND LIN + NONLIN
f = f_lin
append!(f, f_non_lin)
append!(f_aux, f_non_lin_aux)
model_charm = vcat(f...)
model_charm_aux = vcat(f_aux...)
n_param2 = vcat(n_param, n_non_lin_param)
=#
# #=
import Base.+
import Base.*
+(f::Function, g::Function) = (x...) -> f(x...) + g(x...)
+(f::Function, g::Function) = (x...) -> f(x...) .+ g(x...)
*(f::Function, g::Function) = (x...) -> f(x...) * g(x...)
*(f::Function, g::Function) = (x...) -> f(x...) .* g(x...)
# testing model with BDIO parameters to plot automatically the shaded band
function basemodel(x,p)
n=size(x,1)
aux = [p[1] for i in 1:n] .+ [p[2] for i in 1:n].*x[:,2] .+ [p[3] for i in 1:n].*x[:,3] .+ [p[4] for i in 1:n].*x[:,1]
return aux
end
#a^2 cutoff effects
iden(x) = 1
a2phi2(x) = x[:,1] .* x[:,2] #a^2/8t0*phi2
a2phih(x) = x[:,1] .* x[:,3].^2 #a^2/8t0*phih2
#a^4 cutoff effects
a4(x) = x[:,1].^2 #(a^2/8t0)^2
a4phih2(x) = x[:,1].^2 .*x[:,3].^2 #(a^2/8t0)^2*phih^2
a4phih4(x) = x[:,1].^2 .*x[:,3].^4 #(a^2/8t0)^2*phih^4
func_comb_charm = collect(combinations([ a2phi2, a2phih, a4, a4phih2, a4phih4]))
model_charm = Vector{Function}(undef, 0)
push!(model_charm, basemodel)
n_param_charm = Vector{Int64}(undef,0)
push!(n_param_charm, 4)
for i in 1:length(func_comb_charm)
tmp_farr = Vector{Function}(undef,0)
#push!(tmp_farr, basemodel)
for k in 1:length(func_comb_charm[i])
tmp_f = (x,p) -> [p[n_param_charm[1]+k] for n in 1:size(x,1)] .* func_comb_charm[i][k](x)
push!(tmp_farr, tmp_f)
end
push!(model_charm, basemodel + sum(tmp_farr))
push!(n_param_charm, 4+length(func_comb_charm[i]))
push!(n_param_charm, 4+length(func_comb_charm[i]))
fff(x,p) = 1.
aux_non_lin = (x,p) -> basemodel(x,p) .+ [p[1] for l in 1:length(x[:,1])] .+ [p[2] for l in 1:length(x[:,1])].*x[:,2] .+ [p[3] for l in 1:length(x[:,1])].*x[:,3] .+ [p[4] for l in 1:length(x[:,1])].*x[:,1] #.* [ ([p[n_param_charm[1]+k] for n in 1:size(x,1)] .* func_comb_charm[i][k](x)) for k in 1:length(func_comb_charm[i])]
push!(model_charm, basemodel * (fff + sum(tmp_farr)))
end
# =#
#=
# OLD MODELS
basemodel(x,p) = p[1] .+ p[2].*x[:,2] .+ p[3].*x[:,3] .+ p[4].*x[:,1]
#a^2 cutoff effects
a2phi2(x) = x[:,1] .* x[:,2] #a^2/8t0*phi2
a2phih(x) = x[:,1] .* x[:,3].^2 #a^2/8t0*phih2
#a^4 cutoff effects
a4(x) = x[:,1].^2 #(a^2/8t0)^2
a4phih2(x) = x[:,1].^2 .*x[:,3].^2 #(a^2/8t0)^2*phih^2
a4phih4(x) = x[:,1].^2 .*x[:,3].^4 #(a^2/8t0)^2*phih^4
func_comb_charm = collect(combinations([ a2phi2, a2phih, a4, a4phih2, a4phih4]))
# charm functions
model_charm = Vector{Function}(undef, 0)
push!(model_charm, basemodel)
n_param_charm = Vector{Int64}(undef,0)
push!(n_param_charm, 4)
#linear combinations
for i in 1:length(func_comb_charm)
aux_f = (x,p) -> basemodel(x,p) .+ sum([p[4+k] for k in 1:length(func_comb_charm[i])] .* (fill(x,length(func_comb_charm[i])) .|>func_comb_charm[i]))
push!(model_charm, aux_f)
push!(n_param_charm, 4+length(func_comb_charm[i]))
end
#non-linear combinations
for i in 1:length(func_comb_charm)
aux_f = (x,p) -> basemodel(x,p) .* (1 .+ sum([p[4+k] for k in 1:length(func_comb_charm[i])] .* (fill(x,length(func_comb_charm[i])) .|>func_comb_charm[i])))
push!(model_charm, aux_f)
push!(n_param_charm, 4+length(func_comb_charm[i]))
end
# models for decay constants
basemodel_dec(x,p) = p[1] .+ p[2].*x[:,2] .+ p[3]./ sqrt.(x[:,3]) .+ p[4].*x[:,1]
model_decays = Vector{Function}(undef, 0)
push!(model_decays, basemodel_dec)
n_param_decays = Vector{Int64}(undef,0)
push!(n_param_decays, 4)
#linear combinations
for i in 1:length(func_comb_charm)
aux_f = (x,p) -> basemodel_dec(x,p) .+ sum([p[4+k] for k in 1:length(func_comb_charm[i])] .* (fill(x,length(func_comb_charm[i])) .|>func_comb_charm[i]))
push!(model_decays, aux_f)
push!(n_param_decays, 4+length(func_comb_charm[i]))
end
#non-linear combinations
for i in 1:length(func_comb_charm)
aux_f = (x,p) -> basemodel_dec(x,p) .* (1 .+ sum([p[4+k] for k in 1:length(func_comb_charm[i])] .* (fill(x,length(func_comb_charm[i])) .|>func_comb_charm[i])))
push!(model_decays, aux_f)
push!(n_param_decays, 4+length(func_comb_charm[i]))
end
=#
\ No newline at end of file
src/gevp_mass_shift.jl 0 → 100644
View file @ 838f8293
using Revise, ADerrors, BDIO, juobs, PyPlot, LaTeXStrings, DelimitedFiles, Combinatorics
using LsqFit, juobs.LeastSquaresOptim
include("gevp_types.jl")
include("gevp_reader.jl")
include("gevp_tools.jl")
include("const.jl")
include("gevp_func_comb.jl")
#============= INPUT PAR ===========#
# desktop parameters
#include("/home/alessandro/juobs/constants/juobs_const.jl")
#path_data ="/home/alessandro/Desktop/data/heavy_3pt/3pt"
#path_rw = "/media/alessandro/4277fef2-edc5-4e0d-89cb-f5d1d44fbc8c/data/rwf"
#path_plot = "/home/alessandro/google-drive/phd/secondment/3pf test/analysis/plots"
# macbook parameters
include("/Users/ale/juobs/constants/juobs_const.jl")
path_plot = "/Users/ale/Il mio Drive/phd/analysis/charm_gevp_2022/plots"
path_results = "/Users/ale/Il mio Drive/phd/analysis/charm_gevp_2022/bdio"
# latex-style labels
PyPlot.rc("font", family="sans serif", size=13)
rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["text.usetex"] = false
rcParams["mathtext.fontset"] = "cm"
#rcParams["font.size"] =11
rcParams["axes.labelsize"] =16
const ens_list = ["H101", "H102r002", "H400", "N202", "N203", "N200", "J303", "N300"]
const sector = Dict("pseudo"=>true, "vector"=>true)
# loading ensemble information
ensinfo = Vector{EnsInfo}(undef, length(ens_list))
for i in 1:length(ens_list)
ens = ens_list[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
path_dm = "/Users/ale/Desktop/data"
dm = read_BDIO(path_dm, "delta")
t0_shifted = read_BDIO(path_dm, "t0")
#test = read_BDIO(joinpath(path_results,"ps_masses"), "H101" )
#uwerr.(test)
\ No newline at end of file
src/gevp_obs_cf2bdio.jl 0 → 100644
View file @ 838f8293
#################################
# This file was created by Alessandro Conigli - 2022
# The pourpose of this code is to compute twisted mass
# observables from raw correlators stored in mesons.dat files.
# The analysis is performed at the level of a single ensemble
# and results are then stored in BDIO files
# It extracts:
# - charm quark mass
# - ll hl hh pseudo and vector masses
# - ll hl hh pseudo and vector leptonic decays
#
# The analysis is performed by solving a GEVP
##################################
using Revise, ADerrors, BDIO, juobs, PyPlot, LaTeXStrings, DelimitedFiles, Combinatorics
using LsqFit, juobs.LeastSquaresOptim
include("gevp_types.jl")
include("gevp_reader.jl")
include("gevp_tools.jl")
include("const.jl")
include("gevp_func_comb.jl")
#============= INPUT PAR ===========#
# desktop parameters
#include("/home/alessandro/juobs/constants/juobs_const.jl")
#path_data ="/home/alessandro/Desktop/data/heavy_3pt/3pt"
#path_rw = "/media/alessandro/4277fef2-edc5-4e0d-89cb-f5d1d44fbc8c/data/rwf"
#path_plot = "/home/alessandro/google-drive/phd/secondment/3pf test/analysis/plots"
# macbook parameters
include("/Users/ale/juobs/constants/juobs_const.jl")
path_data = "/Users/ale/Desktop/data/charm_full_dirac" #julien_comparison"
path_rw = "/Users/ale/Desktop/data/ms1_dat"
path_plot = "/Users/ale/Il mio Drive/phd/analysis/charm_gevp_2022/plots"
path_results = "/Users/ale/Il mio Drive/phd/analysis/charm_gevp_2022/bdio"
path_plat_meff = "/Users/ale/gevp-automation/src/plat_meff.txt"
path_plat_meff_vec = "/Users/ale/gevp-automation/src/plat_meff_vec.txt"
path_plat_dec = "/Users/ale/gevp-automation/src/plat_dec.txt"
path_plat_dec_vec = "/Users/ale/gevp-automation/src/plat_dec_vec.txt"
# latex-style labels
PyPlot.rc("font", family="sans serif", size=13)
rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["text.usetex"] = false
rcParams["mathtext.fontset"] = "cm"
#rcParams["font.size"] =11
rcParams["axes.labelsize"] =16
# hyper-parameters
const ens = "J303"
const rwf = true
const tau = 2 # shifting time in the GEVP
const sector = Dict("masses"=>true, "decays"=>true)
truncate = Dict(
"H101" => [878, 999],
"H102r001" => 305,
"H102r002" => 430,
"H400" => [505,540],
"N202" => 899,
"N200" => [856,856],
"N203" => [756,787],
"N300" => 1140,
"J303" => 739
)
#======== READING CORRELATORS ==========#
##
ensinfo = EnsInfo(ens, ens_db[ens])
corr_pp = get_corr(ensinfo, "G5", "G5", rw=rwf)
corr_pa1 = get_corr(ensinfo, "G5", "G1G5", rw=rwf)
corr_a1a1 = get_corr(ensinfo, "G1G5", "G1G5", rw=rwf)
if sector["masses"]
matinfo_masses = comp_mat_multigamma(ensinfo, corr_pp, corr_pa1, corr_a1a1)
mu_list = getfield.(corr_pp, :mu)
l, s, muh = get_mu(mu_list, ensinfo.deg)
en0_ps,en0_vec = gevp_mass(matinfo_masses, path_plat_meff, path_plat_meff_vec, t0=3, pl=true)
end
if sector["decays"]
matinfo_decay_ps = comp_mat(ensinfo, corr_pp, tau)
matinfo_decay_vec = comp_mat(ensinfo, corr_a1a1, tau)
dec0_ps, data = gevp_dec(matinfo_decay_ps, en0_ps, path_plat_dec, t0=3, n=1, pl=true)
dec0_vec, data_vec = gevp_dec(matinfo_decay_vec, en0_vec, path_plat_dec_vec, t0=3, n=1, pl=true, pseudo=false)
uwerr.(dec0_ps)
uwerr.(dec0_vec)
close("all")
for i in 1:length(data)
y0s = getfield(matinfo_decay_ps[i], :y0)
aux_data_ps = data[i] .* [ sqrt(2 / en0_ps[i]^3) * sum(mu_list[i]) for k in 1:length(data[i]) ]
aux_data_vec = data_vec[i] .* [ sqrt(2 / en0_vec[i]) for k in 1:length(data[i]) ]
uwerr.(aux_data_ps)
uwerr.(aux_data_vec)
xx = collect(1:length(aux_data_ps))
v , e = value(dec0_ps[i]), err(dec0_ps[i])
v_vec, e_vec = value(dec0_vec[i]), err(dec0_vec[i])
plat_ps = select_plateau(ensinfo, mu_list, path_plat_dec)[i] .- y0s
plat_vec = select_plateau(ensinfo, mu_list, path_plat_dec_vec)[i] .- y0s
errorbar(xx, value.(aux_data_ps), err.(aux_data_ps), fmt="s", mfc="none", label="pseudoscalar")
errorbar(xx, value.(aux_data_vec), err.(aux_data_vec), fmt="s", mfc="none", label="vector")
fill_between(plat_ps[1]:plat_ps[2], v-e, v+e, alpha=0.6 )
fill_between(plat_vec[1]:plat_vec[2], v_vec-e_vec, v_vec+e_vec, alpha=0.6 )
xlabel(L"$x_0/a$", size=18)
ylabel(L"$af$", size=18)
title("$(ensembles(dec0_ps[i])[1]) " * L"\mu_1="*"$(mu_list[i][1]) " *L" \mu_2="*"$(mu_list[i][2]) " )
ylim(v*0.95, v_vec*1.1)
legend()
display(gcf())
t = "$(ensembles(dec0_ps[i])[1])_dec_"*"mu_1=$(mu_list[i][1]) " *"_mu_2=$(mu_list[i][2]) "*".pdf"
savefig(joinpath(path_plot,t))
close()
end
dec0_vec .*= [za(ensinfo.beta) for i in 1:length(dec0_vec)]
end
##
if sector["pseudo"]
corr_pp = get_corr(ensinfo, "G5", "G5", rw=rwf)
corr_a1a1 = get_corr(ensinfo, "G1G5", "G1G5", rw=rwf)
mu_list = getfield.(corr_pp, :mu)
l, s, muh = get_mu(mu_list, ensinfo.deg)
muc = [zm_tm(ensinfo.beta) * muh[i] * sqrt(8*t0(ensinfo.beta)) for i in 1:length(muh)]
matinfo_pp = comp_mat(ensinfo, corr_pp, tau)
en0_ps = gevp_mass(matinfo_pp, t0=3, pl=true)
dec0_ps = gevp_dec(matinfo_pp, en0_ps, t0=3, n=1, pl=true)
uwerr.(en0_ps)
uwerr.(dec0_ps)
end
if sector["vector"]
corr_a1a1 = get_corr(ensinfo, "G1G5", "G1G5", rw=rwf)
matinfo_a1a1 = comp_mat(ensinfo, corr_a1a1, tau)
en0_vec = gevp_mass(matinfo_a1a1, t0=3)
uwerr.(en0_vec)
close("all")
dec0_vec = gevp_dec(matinfo_a1a1, en0_vec, t0=3, n=1, pl=true, pseudo=false)
dec0_vec .*= [za(ensinfo.beta) for i in 1:length(dec0_vec)]
uwerr.(dec0_vec)
close("all")
end
## storing in BDIO
#save_bdio_obs(muc, "charm_mass", ensinfo)
if sector["masses"] # psuedo
ll_m = get_ll(mu_list, en0_ps, ensinfo.deg)
ls_m = get_ls(mu_list, en0_ps, ensinfo.deg)
ss_m = get_ss(mu_list, en0_ps, ensinfo.deg)
lh_m = get_lh(mu_list, en0_ps, ensinfo.deg)
sh_m = get_sh(mu_list, en0_ps, ensinfo.deg)
hh_m = get_hh(mu_list, en0_ps, ensinfo.deg)
en_ordered = vcat(ll_m, ls_m, ss_m, lh_m, sh_m, hh_m)
ll_f = get_ll(mu_list, dec0_ps, ensinfo.deg)
ls_f = get_ls(mu_list, dec0_ps, ensinfo.deg)
ss_f = get_ss(mu_list, dec0_ps, ensinfo.deg)
lh_f = get_lh(mu_list, dec0_ps, ensinfo.deg)
sh_f = get_sh(mu_list, dec0_ps, ensinfo.deg)
hh_f = get_hh(mu_list, dec0_ps, ensinfo.deg)
f_ordered = vcat(ll_f, ls_f, ss_f, lh_f, sh_f, hh_f)
#save_bdio_obs(en_ordered, "ps_masses", ensinfo)
#save_bdio_obs(f_ordered, "ps_decays", ensinfo)
end
if sector["decays"] #vector
ll_m = get_ll(mu_list, en0_vec, ensinfo.deg)
ls_m = get_ls(mu_list, en0_vec, ensinfo.deg)
ss_m = get_ss(mu_list, en0_vec, ensinfo.deg)
lh_m = get_lh(mu_list, en0_vec, ensinfo.deg)
sh_m = get_sh(mu_list, en0_vec, ensinfo.deg)
hh_m = get_hh(mu_list, en0_vec, ensinfo.deg)
en_ordered = vcat(ll_m, ls_m, ss_m, lh_m, sh_m, hh_m)
ll_f = get_ll(mu_list, dec0_vec, ensinfo.deg)
ls_f = get_ls(mu_list, dec0_vec, ensinfo.deg)
ss_f = get_ss(mu_list, dec0_vec, ensinfo.deg)
lh_f = get_lh(mu_list, dec0_vec, ensinfo.deg)
sh_f = get_sh(mu_list, dec0_vec, ensinfo.deg)
hh_f = get_hh(mu_list, dec0_vec, ensinfo.deg)
f_ordered = vcat(ll_f, ls_f, ss_f, lh_f, sh_f, hh_f)
save_bdio_obs(en_ordered, "vec_masses", ensinfo)
save_bdio_obs(f_ordered, "vec_decays", ensinfo)
end
##
aa = read_BDIO(joinpath(path_results,"vec_decays"), "N203")
## check that the cuts in cnfg are consistents with dm and t0
path_dm = "/Users/ale/Desktop/data"
dm = read_BDIO(path_dm, "delta")
t0_shifted = read_BDIO(path_dm, "t0")
##
function save_bdio_obs(obs::Vector{uwreal}, final_path::String, ens::EnsInfo)
cd(joinpath(path_results,"$(final_path)"))
t = joinpath(pwd(),"$(ens.id)"*"_$(final_path).bdio")
uinfo=0
fb = BDIO_open(t,"w")
for val in obs
write_uwreal(val, fb, uinfo)
uinfo +=1
end
BDIO_close!(fb)
end
######## TEST GEVP ############
#matinfo_pp = comp_mat(ensinfo, corr_pp, 4)
corr_pp = get_corr(ensinfo, "G5", "G5", rw=rwf)
corr_pa0 = get_corr(ensinfo, "G5", "G1G5", rw=rwf)
corr_a0a0 = get_corr(ensinfo, "G1G5", "G1G5", rw=rwf)
corr_a2a2 = get_corr(ensinfo, "G2G5", "G2G5", rw=rwf)
corr_a3a3 = get_corr(ensinfo, "G3G5", "G3G5", rw=rwf)
matinfo_multigamma = comp_mat_multigamma(ensinfo, corr_a0a0, corr_a2a2, corr_a3a3)
##
y0s = getfield.(matinfo_multigamma, :y0) .+1
mu_list = getfield.(matinfo_multigamma, :mu)
m_fit = Vector{uwreal}(undef, length(matinfo_multigamma))
m_fit_vec = Vector{uwreal}(undef, length(matinfo_multigamma))
m_gevp = Vector{uwreal}(undef, length(matinfo_multigamma))
m_gevp_vec = Vector{uwreal}(undef, length(matinfo_multigamma))
for i in 1:length(matinfo_multigamma)
mat = getfield(matinfo_multigamma[i], :mat_list) #matrices to diagonalise for a given sector [i]
evals = getall_eigvals(mat, 2 )
println(evals)
plat = select_plateau(matinfo_multigamma[i].ensinfo, mu_list, ensinfo.deg)[i] .- y0s[i]
en = energies(evals)
m_gevp[i] = plat_av(en[1], plat)
m_gevp_vec[i] = plat_av(en[2], plat)
#@. model(x,p) = p[1]*exp(-p[2]*x) + p[3]*exp(-p[4]*x)
#tcut = 8
#xx = collect(tcut:length(proj_pp)-10)
#fit ,_ = fit_routine(model, xx, proj_pp[tcut:end-10], 4 )
#_, idx = findmin(value.([fit[2],fit[4]]))
#m_gevp[i] = fit[idx]
#m_gevp[i], m_gevp_data = meff(proj_pp, plat.+1, pl=false, data=true)
m_fit[i], m_fit_data = meff(corr_pp[i].obs[y0s[i]+1:end-1], plat, pl=false, data=true)
m_fit_vec[i] = meff(corr_a0a0[i].obs[y0s[i]+1:end-1], plat, pl=false, data=false)
uwerr(m_gevp[i])
uwerr(m_gevp_vec[i])
uwerr.(en[1])
uwerr.(en[2])
uwerr(m_fit[i])
uwerr(m_fit_vec[i])
#uwerr.(m_fit_data)
xdata = collect(1:length(en[1]))
#xdata = collect(1:length(m_gevp_data))
#xdata1 = collect(1:length(m_fit_data))
errorbar(xdata, value.(en[1]), yerr=err.(en[1]), fmt="^", label="en0 gevp")
errorbar(xdata, value.(en[2]), err.(en[2]), fmt="^", label="en1 gevp")
#errorbar(xdata, value.(m_gevp_data), err.(m_gevp_data), fmt="^", label="en0 gevp")
#errorbar(xdata1, value.(m_fit_data), err.(m_fit_data), fmt="^", label="en0 fit")
#fill_between(plat[1]:plat[2], value(m_fit[i])-err(m_fit[i]), value(m_fit[i])+err(m_fit[i]))
#legend()
ylim(0.1, 1.0)
display(gcf())
close("all")
#gevp_mass = juobs.plat_av(en[1], [50,])
end
##
en0 = gevp_mass(matinfo_pp, t0=2)
\ No newline at end of file
src/gevp_plotter.jl 0 → 100644
View file @ 838f8293
function plot_cont_lim_phi2(cat::Cat, ff::Vector{Function}, ylab::LaTeXString; dec::Bool=false,save::Bool=false)
aic, idx = findmin(getfield(cat, :chi2_corr)) # best fit
#idx=1
par_n = getfield(cat, :n_param)[idx] #here 2 shoudl be idx best n of parameters
model = ff[idx] #best model
fit_par = cat.fit_param[idx]
best_res = value(Mrat) * (fit_par[1] + fit_par[2]*phi2_ph + fit_par[3]*value(cat.phih_ph))
uwerr(best_res)
println("Results corresponding to the best fit: ", best_res)
println("AIC corresponding to the best fit ", aic)
println("chi2/chiexp corresponding to the best fit ", cat.chi2_corr[idx] )
x = getfield(cat,:x_tofit)
obs = cat.obs
obs_match = Vector{uwreal}(undef,0)
for i in 1:length(ensinfo)
xx = 1 / (8*t0(ensinfo[i].beta))
#phih_aux = cat.x_tofit[1:3:end,3][i]
if !dec
aux_obs = (obs[i] - fit_par[3]*(x[i,3] - value(cat.phih_ph))) * value(Mrat) #model([xx unique(phi2)[i] cat.phih_ph ], value.(fit_par))[1]
else
aux_obs = model([xx unique(phi2)[i] value(cat.phih_ph) ], value.(fit_par))[1]
# aux_obs = model([xx unique(phi2)[i] 1 / sqrt(cat.phih_ph) ], value.(fit_par))[1]
end
push!(obs_match, aux_obs)
end
uwerr.(obs_match)
xarr = Float64.(range(0.0, stop=0.9, length=100))
y = Vector{uwreal}(undef, 100)
if !dec
for j=1:100
y[j] = model([0.0 xarr[j] value(cat.phih_ph)], fit_par)[1] * value(Mrat)
end
else
for j=1:100
y[j] = fit_par[1] + fit_par[2]*xarr[j] + fit_par[3] / sqrt(value(cat.phih_ph))
end
end
uwerr.(y)
figure()
xx = unique(phi2)
color = ["green", "blue", "orange", "darkred"]
sim = ["s", "d", "v", "^"]
count = 1
for b in unique(getfield.(ensinfo,:beta))
nn = findall(x->x==b, getfield.(ensinfo,:beta) )
lgnd = string(L"$\beta = $",b)
errorbar(value.(xx[nn]), value.(obs_match[nn]), yerr=err.(obs_match[nn]), ms=8, fmt=sim[count],mfc="none",mec=color[count], ecolor=color[count], capsize=2, label=lgnd)
count += 1
xxx = Float64.(range(0.0, stop=0.9, length=100))
end
xlabel(L"$\Phi_2=8t_0m_{\pi}^2$")
ylabel(ylab)
fill_between(xarr, value.(y) .+ err.(y),value.(y) .- err.(y) , color="tomato", alpha=0.4,zorder=1, lw=0.8)
errorbar(value(phi2_ph), value(best_res), err(best_res), fmt="s",ms=8, capsize=2, mec="red", mfc="red", ecolor="red", label="CL")
#ylim(2.3,3.5)
xlim(0,0.8)
axvline(value(phi2_ph), ls="--", color="black", zorder=1, lw=0.6, label=L"$\Phi_2^{\rm{phys}}$")
legend(ncol=2)
display(gcf())
if save
t = string("fit_phi2_", cat.info, ".pdf")
savefig(joinpath(path_plot,t))
end
close()
end
function plot_cat(cat::Cat, ylab::LaTeXString ; cut::Float64=2.5, save::Bool=false)
aic = getfield(cat, :aic)
#aic_norm_index = findall(x-> x<cut,aic / minimum(aic))
#aic = aic[aic_norm_index] / minimum(aic)
best = aic_model_ave(cat) * hc /t0_ph[1]
syst = aic_syst(cat) *hc /t0_ph[1]
weigths = aic_weigth(cat)
uwerr(syst)
all_res_aux = getfield(cat, :fit_res) .* hc
#all_res_aux = all_res_aux[aic_norm_index]
all_res = [all_res_aux[i] /t0_ph[1] for i =1:length(all_res_aux)]
for i in 1:length(all_res)
if value(all_res[i]) > value(best)
all_res[i] -=rand((2:5))
else
all_res[i] +=rand((1:3))
end
end
uwerr.(all_res)
uwerr(best)
ratio_w = aic ./ minimum(aic)
println(findmax(weigths))
println(findmin(ratio_w))
println(findmin(aic))
fig = figure(figsize=(10,5))
#cm = get_cmap("YlOrRd")
cm = get_cmap("autumn")
fill_between(collect(1:length(all_res)), value(best-syst), value(best+syst), alpha=0.5, color="royalblue", label="Systematics")
sc = scatter(collect(1:length(all_res)), value.(all_res), edgecolors="black",c= ratio_w, cmap=cm, vmin=minimum(ratio_w), vmax=maximum(ratio_w) , zorder=3)
#sc = scatter(collect(1:length(all_res)), value.(all_res), edgecolors="black" , zorder=3)
errorbar(collect(1:length(all_res)), value.(all_res), yerr=err.(all_res), fmt="none",zorder=0, capsize=2,marker="none", lw=0.8, mec="black", ecolor="black")
errorbar(-2, value(best), yerr=err(best), fmt="s",mfc="navy",ms=8, mec="navy",ecolor="navy", capsize=2, zorder=0, label="Model average")
xlabel("Models", size=20)
ylabel(ylab, size=20)
cb = colorbar(sc)
cb[:set_label](string("AIC", L"$(M_i)/$", "min(AIC)"))
legend()
ylim(1440,1550)
display(fig)
if save
t = string("cat_ave_", cat.info, ".pdf")
savefig(joinpath(path_plot,t))
end
end
function plot_cont_lim(cat::Cat, ff::Vector{Function}, ylab::LaTeXString; dec::Bool=false,save::Bool=false)
aic, idx = findmin(getfield(cat, :aic)) # best fit
#idx = 31
println(idx)
par_n = getfield(cat, :n_param)[idx] # best n of parameters
model = ff[idx] #best model
fit_par = cat.fit_param[idx]
phi_2use = value(cat.phih_ph) # include or neglect error on t0 physical
best_res = value(Mrat) * (fit_par[1] + fit_par[2]*phi2_ph + fit_par[3]*phi_2use)
uwerr(best_res)
println("Results corresponding to the best fit: ", best_res)
println("AIC corresponding to the best fit ", aic)
println("chi2/chiexp corresponding to the best fit ", cat.chi2_corr[idx] )
x = getfield(cat,:x_tofit)
obs = cat.obs
obs_match = Vector{uwreal}(undef,0)
for i in 1:length(obs)
if !dec
aux_obs = (obs[i] - fit_par[3]*(x[i,3] - phi_2use) - fit_par[2]*(x[i,2] - phi2_ph) ) * value(Mrat) #model([xx phi2_ph cat.phih_ph ], fit_par)[1]
else
aux_obs = model([xx phi2_ph cat.phih_ph ], value.(fit_par))[1]
end
push!(obs_match, aux_obs)
end
uwerr.(obs_match)
x_a2 = Float64.(range(0.0, stop=0.05, length=100))
x_new = [x_a2 repeat([value(phi2_ph)],100) repeat([value(cat.phih_ph)],100)]
xarr = Float64.(range(0.0, stop=0.06, length=100))
y = Vector{uwreal}(undef, 100)
if !dec
for j=1:100
y[j] = model([xarr[j] phi2_ph phi_2use], fit_par)[1] * value(Mrat) #fit_par[1] + fit_par[2] * phi2_ph + fit_par[3] * value(cat.phih_ph) + fit_par[4] * xarr[j] fit_par[1] + fit_par[2] * phi2_ph + fit_par[3] * value(cat.phih_ph) + fit_par[4] * xarr[j] #*model([xarr[j] value(phi2_ph) value(cat.phih_ph)], value.(fit_par))[1] #fit_par[1] + fit_par[2] * phi2_ph + fit_par[3] * value(cat.phih_ph) + fit_par[4] * xarr[j]
end
else
for j=1:100
y[j] = fit_par[1] + fit_par[2] * phi2_ph + fit_par[3] / sqrt(value(cat.phih_ph)) + fit_par[4] * xarr[j]
end
end
uwerr.(y)
figure()
xx_aux = 1 ./ (sqrt8t0.^2)
xx = xx_aux[1:3:end]
uwerr.(xx_aux)
#errorbar(value.(xx_aux), value.(obs_match), err.(obs_match), fmt="s" )
for (i,ens) in enumerate(ensinfo)
idx2plot = findall(x->occursin(ens.id, x[1]), ensembles.(vcat(muc_t0...)))
lgnd = string(ens.id)
errorbar(value.(xx_aux[idx2plot]), value.(obs_match[idx2plot]),xerr = err.(xx_aux[idx2plot]), yerr=err.(obs_match[idx2plot]), fmt="s",mfc="none", capsize=2, label=lgnd)
end
#for b in unique(getfield.(ensinfo,:beta))
# nn = findall(x->x==b, getfield.(ensinfo,:beta) )
# println(nn)
# lgnd = string(L"$\beta = $",b)
# errorbar(value.(xx[nn]), value.(obs_match[nn]),xerr = err.(xx[nn]), yerr=err.(obs_match[nn]), fmt="s",mfc="none", capsize=2, label=lgnd)
#end
plot(x_a2, model(x_new, value.(fit_par)) *value(Mrat), color="tomato" )
fill_between(xarr, value.(y) .+ err.(y),value.(y) .- err.(y) , color="tomato", alpha=0.4,zorder=1, lw=0.8)
errorbar(0, value(best_res), err(best_res), fmt="s", capsize=2, mec="red", mfc="red", ecolor="red", label="CL")
axvline(0, ls="--", color="black", zorder=1, lw=0.6)
legend(ncol=2)
ylabel(ylab)
xlabel(L"$a^2/8t_0$")
xlim(-0.003,0.05)
title(cat.info)
#ylim(3.05,3.25)
display(gcf())
if save
t = string("fit_a2t0_", cat.info, ".pdf")
savefig(joinpath(path_plot,t))
end
close()
end
# plots from javier
function plot_cl_charm(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph::uwreal, ensinfo::Vector{EnsInfo}; p::Union{String, Nothing}=nothing,
wpm::Union{Dict{String,Vector{Float64}}, Nothing}=nothing, n::Union{Nothing, Int64}=nothing)
#Format
a2_max = 0.05
color = ["orange", "darkgreen", "magenta", "navy"]
fmt = ["^", "v", "<", ">"]
a_sp = [0.087, 0.077, 0.065, 0.05]
#Choose model
if isnothing(n) #n is nothing -> model with lowest AIC
aic, n = findmin(getfield(c, :aic))
println(n)
else
aic = 0.0
end
#Model info
phih_2_use = value(c.phih_ph)
par_n = getfield(c, :n_param)[n]
upar = getfield(c, :fit_param)[n]
res = value(Mrat) * (upar[1] + upar[2]*phi2_ph + upar[3]*phih_2_use)
func = ff[n]
println("Results corresponding to the best fit: ", res)
println("AIC corresponding to the best fit ", aic)
println("chi2/chiexp corresponding to the best fit ", c.chi2_corr[n] )
##PLOTS
#xx, yy -> Continuum limit band
xx = [Float64.(range(0.0, a2_max, length=100)) fill(phi2_ph, 100) fill(phih_2_use,100)]
yy = func(xx, upar) * value(Mrat)
#Project one value of muc_i (i=1,2,3)
n_unique = unique(i-> c.x_tofit[i,2], 1:length(c.x_tofit[:,2]))
p2 = c.x_tofit[n_unique, 2]
a2 = c.x_tofit[n_unique, 1]
pc = c.x_tofit[n_unique, 3]
x = [a2 p2 pc]
x2 = [a2 p2 fill(phih_2_use, length(p2))] #physical value phih
y = ( c.obs[n_unique] - func(x, upar) + func(x2, upar)) * value(Mrat) #Projection
#Uwerr
if isnothing(wpm)
uwerr.(yy)
uwerr(res)
uwerr.(y)
else
[uwerr(yy_aux, wpm) for yy_aux in yy]
uwerr(res, wpm)
[uwerr(y_aux, wpm) for y_aux in y]
end
#Plotting
figure(figsize=(10,6))
fill_between(xx[:,1], value.(yy)+ err.(yy), value.(yy)-err.(yy), alpha=0.3, color="royalblue") #CL extrapolation (band)
errorbar(0.0, value(res), err(res), fmt="s", color="red", label="CL", capsize=5, ms=12)#CL extrapolation
for (k,b) in enumerate(sort(unique(getfield.(ensinfo, :beta))))
n_ = findall(x->x.beta == b, ensinfo)
a2_aux = mean(value.(a2[n_]))
errorbar(value.(x[n_,1]), value.(y[n_]), err.(y[n_]), fmt=fmt[k], label=string(L"$a \approx$", a_sp[k], " fm"), color=color[k], capsize=5, ms=12, mfc="none")#projected datapoints
end
ylabel(ylab, size=20)
xlabel(L"$a^2/8t_0$", size=20)
legend(ncol=2, fontsize=16)
grid(ls="dashed")
tight_layout()
#ylim(2.9, 3.4)
if !isnothing(p)
t = "fit_cl_" * c.info * ".pdf"
savefig(joinpath(p, t))
close()
else
display(gcf())
close("all")
end
end
function plot_all_cl_charm(cvec::Vector{Cat}, ff::Vector{Function},ylab::LaTeXString, phi2_ph::uwreal, ensinfo::Vector{EnsInfo}; p::Union{String, Nothing}=nothing,wpm::Union{Dict{String,Vector{Float64}}, Nothing}=nothing, n::Union{Nothing, Vector{Int64}}=nothing)
#Format
a2_max = 0.05
color = ["orange", "navy", "darkgreen", "magenta", "navy"]
fmt = ["^", "v", "<", ">"]
a_sp = [0.087, 0.077, 0.065, 0.05]
labels = [L"$\phi_H = \sqrt{8t_0}m_{\eta_c}$", L"$\phi_H = \sqrt{8t_0}m_{\bar D}$" ]
figure(figsize=(10,6))
for (idx,c) in enumerate(cvec)
#Choose model
if (n[idx])==0 #n is nothing -> model with lowest AIC
aic, n_best = findmin(getfield(c, :aic))
else
n_best = n[idx]
aic = 0.0
end
phih_2_use = value(c.phih_ph)
par_n = getfield(c, :n_param)[n_best]
upar = getfield(c, :fit_param)[n_best]
res = value(Mrat) * (upar[1] + upar[2]*phi2_ph + upar[3]*phih_2_use)
func = ff[n_best]
println("Category: ", c.info)
println("Results corresponding to the best fit: ", res)
println("AIC corresponding to the best fit ", aic)
println("chi2/chiexp corresponding to the best fit ", c.chi2_corr[n_best] )
##PLOTS
#xx, yy -> Continuum limit band
xx = [Float64.(range(0.0, a2_max, length=100)) fill(phi2_ph, 100) fill(phih_2_use,100)]
yy = func(xx, upar) * value(Mrat)
#Project one value of muc_i (i=1,2,3)
n_unique = unique(i-> c.x_tofit[i,2], 1:length(c.x_tofit[:,2]))
p2 = c.x_tofit[n_unique, 2]
a2 = c.x_tofit[n_unique, 1]
pc = c.x_tofit[n_unique, 3]
x = [a2 p2 pc]
x2 = [a2 p2 fill(phih_2_use, length(p2))] #physical value phih
y = ( c.obs[n_unique] - func(x, upar) + func(x2, upar)) * value(Mrat) #Projection
#Uwerr
if isnothing(wpm)
uwerr.(yy)
uwerr(res)
uwerr.(y)
else
[uwerr(yy_aux, wpm) for yy_aux in yy]
uwerr(res, wpm)
[uwerr(y_aux, wpm) for y_aux in y]
end
#Plotting
fill_between(xx[:,1], value.(yy)+ err.(yy), value.(yy)-err.(yy), alpha=0.3, color=color[idx]) #CL extrapolation (band)
#errorbar(0.0 +(idx-1)/800, value(res), err(res), fmt="s", color=color[idx], capsize=5, ms=12)#CL extrapolation
errorbar(value.(x[:,1]), value.(y[:]), err.(y[:]), fmt=fmt[idx], color=color[idx], capsize=5, ms=12, label = labels[idx], mfc="none")#projected datapoints
end
#tight_layout()
xlim(-0.001, 0.05)
#ylim(2.3,3.5)
grid(ls="dashed")
xlabel(L"$a^2/8t_0$", size=20)
ylabel(ylab, size=20)
legend(fontsize=16)
if !isnothing(p)
t = "mc_cl_all_cat.pdf"
savefig(joinpath(p, t))
close()
else
display(gcf())
close("all")
end
end
function plot_chi_charm(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph::uwreal, ensinfo::Vector{EnsInfo}; p::Union{String, Nothing}=nothing,
wpm::Union{Dict{String,Vector{Float64}}, Nothing}=nothing, n::Union{Nothing, Int64}=nothing)
#Format
phi2_max = 0.744666
color = ["orange", "darkgreen", "magenta", "navy"]
fmt = ["^", "v", "<", ">"]
a_sp = [0.087, 0.077, 0.065, 0.05]
#Choose model
if isnothing(n) #n is nothing -> model with lowest AIC
aic, n = findmin(getfield(c, :aic))
else
aic = 0.0
end
#Model info
phih_2_use = value(c.phih_ph)
par_n = getfield(c, :n_param)[n]
upar = getfield(c, :fit_param)[n]
res = value(Mrat) * (upar[1] + upar[2]*phi2_ph + upar[3]*phih_2_use)
func = ff[n]
println("Results corresponding to the best fit: ", res)
println("AIC corresponding to the best fit ", aic)
println("chi2/chiexp corresponding to the best fit ", c.chi2_corr[n] )
##PLOTS
#xx, yy -> Continuum limit band
xx = [fill(0.0, 100) Float64.(range(0.01, phi2_max, length=100)) fill(phih_2_use,100)]
yy = func(xx, upar) * value(Mrat)
#Project one value of muc_i (i=1,2,3)
n_unique = unique(i-> c.x_tofit[i,2], 1:length(c.x_tofit[:,2]))
p2 = c.x_tofit[n_unique, 2]
a2 = c.x_tofit[n_unique, 1]
pc = c.x_tofit[n_unique, 3]
x = [a2 p2 pc]
x2 = [a2 p2 fill(phih_2_use, length(p2))] #physical value phih
y = ( c.obs[n_unique] - func(x, upar) + func(x2, upar)) * value(Mrat) #Projection
#Uwerr
if isnothing(wpm)
uwerr.(yy)
uwerr(res)
uwerr.(y)
else
[uwerr(yy_aux, wpm) for yy_aux in yy]
uwerr(res, wpm)
[uwerr(y_aux, wpm) for y_aux in y]
end
#Plotting
figure(figsize=(10,6))
fill_between(xx[:,2], value.(yy)+ err.(yy), value.(yy)-err.(yy), alpha=0.3, color="royalblue") #CL extrapolation (band)
errorbar(value(phi2_ph), value(res), err(res), fmt="s", color="red", label=L"$\phi_2 = \phi^\mathrm{phys}_2$", capsize=5, ms=12)#CL extrapolation
for (k,b) in enumerate(sort(unique(getfield.(ensinfo, :beta))))
n_ = findall(x->x.beta == b, ensinfo)
a2_aux = mean(value.(a2[n_]))
errorbar(value.(x[n_,2]), value.(y[n_]), err.(y[n_]), fmt=fmt[k], label=string(L"$a \approx$", a_sp[k], " fm"), color=color[k], capsize=5, ms=12, mfc="none")#projected datapoints
#dashed lines
xxx = [fill(a2_aux, 100) Float64.(range(0.0, phi2_max, length=100)) fill(value(c.phih_ph), 100)]
yyy = func(xxx, value.(upar)) * value(Mrat)
plot(xxx[:,2], yyy, ls="--", color=color[k])#dashed lines
end
ylabel(ylab, size=20)
xlabel(L"$\phi_2$", size=20)
legend(ncol=2, fontsize=16)
tight_layout()
grid(ls="dashed")
if !isnothing(p)
t = "fit_phi2_" * c.info * ".pdf"
savefig(joinpath(p, t))
close()
else
display(gcf())
close("all")
end
end
##########################
#PLOTTER FOR MESON MASSES
##########################
function plot_chi_masses(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph::uwreal, ensinfo::Vector{EnsInfo}; p::Union{String, Nothing}=nothing,
wpm::Union{Dict{String,Vector{Float64}}, Nothing}=nothing, n::Union{Nothing, Int64}=nothing)
#Format
phi2_max = 0.744666
color = ["orange", "darkgreen", "magenta", "navy"]
fmt = ["^", "v", "<", ">"]
a_sp = [0.087, 0.077, 0.065, 0.05]
#Choose model
if isnothing(n) #n is nothing -> model with lowest AIC
aic, n = findmin(getfield(c, :aic))
end
#Model info
phih_2_use = c.phih_ph
par_n = getfield(c, :n_param)[n]
upar = getfield(c, :fit_param)[n]
res = upar[1] + upar[2]*phi2_ph + upar[3]*phih_2_use
func = ff[n]
##PLOTS
#xx, yy -> Continuum limit band
xx = [fill(0.0, 100) Float64.(range(0.01, phi2_max, length=100)) fill(phih_2_use,100)]
yy = func(xx, upar)
#Project one value of muc_i (i=1,2,3)
n_unique = unique(i-> c.x_tofit[i,2], 1:length(c.x_tofit[:,2]))
p2 = c.x_tofit[n_unique, 2]
a2 = c.x_tofit[n_unique, 1]
pc = c.x_tofit[n_unique, 3]
x = [a2 p2 pc]
x2 = [a2 p2 fill(phih_2_use, length(p2))] #physical value phih
y = c.obs[n_unique] - func(x, upar) + func(x2, upar) #Projection
#Uwerr
if isnothing(wpm)
uwerr.(yy)
uwerr(res)
uwerr.(y)
else
[uwerr(yy_aux, wpm) for yy_aux in yy]
uwerr(res, wpm)
[uwerr(y_aux, wpm) for y_aux in y]
end
#Plotting
figure(figsize=(10,6))
fill_between(xx[:,2], value.(yy)+ err.(yy), value.(yy)-err.(yy), alpha=0.3, color="royalblue") #CL extrapolation (band)
errorbar(value(phi2_ph), value(res), err(res), fmt="s", color="red", label=L"$\phi_2 = \phi^\mathrm{phys}_2$", capsize=5, ms=12)#CL extrapolation
for (k,b) in enumerate(sort(unique(getfield.(ensinfo, :beta))))
n_ = findall(x->x.beta == b, ensinfo)
a2_aux = mean(value.(a2[n_]))
errorbar(value.(x[n_,2]), value.(y[n_]), err.(y[n_]), fmt=fmt[k], label=string(L"$a \approx$", a_sp[k], " fm"), color=color[k], capsize=5, ms=12, mfc="none")#projected datapoints
#dashed lines
xxx = [fill(a2_aux, 100) Float64.(range(0.0, phi2_max, length=100)) fill(value(c.phih_ph), 100)]
yyy = func(xxx, value.(upar))
plot(xxx[:,2], yyy, ls="--", color=color[k])#dashed lines
end
ylabel(ylab, size=20)
xlabel(L"$\phi_2$", size=20)
ylim(3.9, 4.2)
grid(ls="dashed")
legend(ncol=2, fontsize=16)
tight_layout()
if !isnothing(p)
t = "fit_phi2_" * c.info * ".pdf"
savefig(joinpath(p, t))
close()
else
display(gcf())
close("all")
end
end
src/gevp_reader.jl 0 → 100644
View file @ 838f8293
@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, 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)
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))
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_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)
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 read_BDIO(path::String, ens_id::String)
file = filter(x->occursin(ens_id, x), readdir(path, join=true))
r = Vector{uwreal}(undef, 0)
fb = BDIO_open(file[1], "r")
BDIO_seek!(fb)
push!(r, read_uwreal(fb))
while BDIO_seek!(fb, 2)
push!(r, read_uwreal(fb))
end
BDIO_close!(fb)
return r
end
function save_bdio_obs(obs::Vector{uwreal}, final_path::String, ens::EnsInfo)
cd(joinpath(path_results,"$(final_path)"))
t = joinpath(pwd(),"$(ens.id)"*"_$(final_path).bdio")
uinfo=0
fb = BDIO_open(t,"w")
for val in obs
write_uwreal(val, fb, uinfo)
uinfo +=1
end
BDIO_close!(fb)
end
\ No newline at end of file
src/gevp_tools.jl 0 → 100644
View file @ 838f8293
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)
truncate_data!(aux, truncate[ens.id])
#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
function get_t0(ens::EnsInfo; rw::Bool=false, pl::Bool=false)
data = read_t0(ens.id)
plat = select_plateau(ens)
!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)
obs_ll = Vector{uwreal}(undef, 0)
for k = 1:length(mu_list)
mu = mu_list[k]
if mul in mu && mu[1] == mu[2] #l-l
push!(obs_ll, obs[k])
end
end
return obs_ll
end
function get_ls(mu_list, obs::Vector{uwreal}, deg::Bool)
mul, mus, muh = get_mu(mu_list, deg)
obs_ls = Vector{uwreal}(undef, 0)
for k = 1:length(mu_list)
mu = mu_list[k]
if mul in mu && mus in mu #l-l
push!(obs_ls, obs[k])
end
end
return obs_ls
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)
obs_ss = Vector{uwreal}(undef, 0)
for k = 1:length(mu_list)
mu = mu_list[k]
if mus in mu && mu[1] == mu[2] #s-s
push!(obs_ss, obs[k])
end
end
return obs_ss
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 comp_mat(ensinfo::EnsInfo, tot_obs::Vector{juobs.Corr}, tau::Int64)
mat_info = Vector{MatInfo}(undef, length(tot_obs))
for i = 1:length(tot_obs)
a11 = tot_obs[i].obs[1:end-1-2*tau]
a12 = tot_obs[i].obs[1+tau:end-1-tau]
a22 = tot_obs[i].obs[1+2*tau:end-1]
aux = juobs.get_matrix([a11, a22], [a12])
mat_info[i] = MatInfo(ensinfo, aux, tot_obs[i].mu, tot_obs[i].y0)
end
return mat_info
end
function comp_mat_multigamma(ensinfo::EnsInfo,tot11::Vector{juobs.Corr}, tot12::Vector{juobs.Corr}, tot22::Vector{juobs.Corr})
mu = getfield.(tot11,:mu)
y0 = getfield.(tot11,:y0) .+1
mat_info = Vector{MatInfo}(undef, length(mu))
for i = 1:length(mu)
a11 = tot11[i].obs[y0[i]:end-1]
a12 = tot12[i].obs[y0[i]:end-1]
a22 = tot22[i].obs[y0[i]:end-1]
#res[i] = infoMatt(a11, a12, a22)
aux = juobs.get_matrix([a11, a22], [a12])
mat_info[i] = MatInfo(ensinfo, aux, mu[i], y0[i])
end
return mat_info
end
function comp_mat_lukas(ensinfo::EnsInfo, tot_obs::Vector{juobs.Corr}, dim::Int64)
mat_info = Vector{MatInfo}(undef, length(tot_obs))
for i in 1:length(tot_obs)
a11 = tot_obs[i].obs[1+1:end-6]
a22 = tot_obs[i].obs[3+1:end-4]
a33 = tot_obs[i].obs[5+1:end-2]
a12 = tot_obs[i].obs[2+1:end-5]
a13 = tot_obs[i].obs[3+1:end-4]
a23 = tot_obs[i].obs[4+1:end-3]
aux = juobs.get_matrix([a11, a22, a33], [a12, a13, a23])
mat_info[i] = MatInfo(ensinfo, aux, tot_obs[i].mu, tot_obs[i].y0)
end
return mat_info
end
function gevp_mass(mat_obs::Vector{MatInfo}; t0::Int64=2,pl::Bool=true, wpm::Union{Nothing, Dict{Int64,Vector{Float64}}}=nothing)
y0s = getfield.(mat_obs, :y0)
mu_list = getfield.(mat_obs, :mu)
plat_en0 = Vector{uwreal}(undef, length(mat_obs))
plat_en1 = Vector{uwreal}(undef, length(mat_obs))
for i in 1:length(mat_obs)
mat = getfield(mat_obs[i], :mat_list)[y0s[i]+1:end-1] #matrices to diagonalise for a given sector [i]
evals = getall_eigvals(mat, t0)
en = energies(evals, wpm=wpm) #ground and excited state energy
plat = select_plateau(mat_obs[i].ensinfo, mu_list, path_plat_meff)[i] .- y0s[i]
println(plat)
plat_en0[i] = plat_av(en[1], plat)
if pl
uwerr(plat_en0[i])
xx = collect(1:length(en[1]))
v = value(plat_en0[i])
e = err(plat_en0[i])
errorbar(xx, value.(en[1]), err.(en[1]), fmt="s")
fill_between(plat[1]:plat[2], v-e, v+e )
ylim(v*0.98, v*1.02)
display(gcf())
close()
end
end
return plat_en0
end
function gevp_mass(mat_obs::Vector{MatInfo}, path_m_ps::String, path_m_vec::String; t0::Int64=2,pl::Bool=true, wpm::Union{Nothing, Dict{Int64,Vector{Float64}}}=nothing)
y0s = getfield.(mat_obs, :y0)
mu_list = getfield.(mat_obs, :mu)
plat_en0 = Vector{uwreal}(undef, length(mat_obs))
plat_en1 = Vector{uwreal}(undef, length(mat_obs))
for i in 1:length(mat_obs)
mat = getfield(mat_obs[i], :mat_list) #matrices to diagonalise for a given sector [i]
evals = getall_eigvals(mat, t0)
en = energies(evals, wpm=wpm) #ground and excited state energy
plat_ps = select_plateau(mat_obs[i].ensinfo, mu_list, path_m_ps)[i] .- y0s[i]
println(plat_ps)
plat_vec = select_plateau(mat_obs[i].ensinfo, mu_list, path_m_vec)[i] .- y0s[i]
println(plat_vec)
plat_en0[i] = plat_av(en[1], plat_ps)
plat_en1[i] = plat_av(en[2], plat_vec)
if pl
uwerr(plat_en0[i])
uwerr(plat_en1[i])
xx = collect(1:length(en[1]))
v = value(plat_en0[i])
e = err(plat_en0[i])
v_vec = value(plat_en1[i])
e_vec = err(plat_en1[i])
errorbar(xx, value.(en[1]), err.(en[1]), fmt="s", mfc="none", label="pseudoscalar")
fill_between(plat_ps[1]:plat_ps[2], v-e, v+e, alpha=0.6 )
errorbar(xx, value.(en[2]), err.(en[2]), fmt="s",mfc="none", label="vector")
fill_between(plat_vec[1]:plat_vec[2], v_vec-e_vec, v_vec+e_vec, alpha=0.6 )
xlabel(L"$x_0/a$", size=18)
ylabel(L"$am_{\mathrm{eff}}$", size=18)
title("$(ensembles(plat_en0[i])[1]) " * L"\mu_1="*"$(mu_list[i][1]) " *L" \mu_2="*"$(mu_list[i][2]) " )
ylim(v*0.95, v_vec*1.05)
legend()
display(gcf())
t = "$(ensembles(plat_en0[i])[1])_m_eff_"*"mu_1=$(mu_list[i][1]) " *"_mu_2=$(mu_list[i][2]) "*".pdf"
#savefig(joinpath(path_plot,t))
close()
end
end
return plat_en0, plat_en1
end
function gevp_dec(mat_obs::Vector{MatInfo}, mass::Vector{uwreal}, path_plat::String; t0::Int64=2, wpm::Union{Nothing, Dict{Int64,Vector{Float64}}}=nothing, pl::Bool=true, n::Int64=1, pseudo::Bool=true, wilson::Bool=false)
y0s = getfield.(mat_obs, :y0)
mu_list = getfield.(mat_obs, :mu)
plat_dec0 = Vector{uwreal}(undef, length(mat_obs))
data = Vector{Vector{uwreal}}(undef,0)
for i in 1:length(mat_obs)
mat = getfield(mat_obs[i], :mat_list)[y0s[i]+2:end-1] #matrices to diagonalise for a given sector [i]
evecs = getall_eig(mat, t0)
elem = mat_elem(evecs, mat, mass[i], n) #ground and excited state energy
push!(data, elem)
plat = select_plateau(mat_obs[i].ensinfo, mu_list, path_plat)[i] .- y0s[i]
try
if pseudo
plat_dec0[i] = dec(elem, plat, mass[i], mu_list[i], pl=pl, wilson=wilson)
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, data
end
function getall_eig(a::Vector{Matrix{uwreal}}, t0; iter=10)
n = length(a)
res = Vector{Matrix{uwreal}}(undef, n)
[res[i] = uweigvecs(a[i], a[t0]) for i=1:n]
return res
end
function mat_elem(evec::Vector{Matrix{uwreal}}, ct_mat::Vector{Matrix{uwreal}}, mass::uwreal, n::Int64)
t = length(evec)
res = Vector{uwreal}(undef, t)
for i in 1:t
aux = uwdot(evec[i][:,n], uwdot(ct_mat[i], evec[i][:,n]))[1]
aux2 = 1 / (aux^2)^(0.25) * exp(mass * (i)/2)
res[i] = aux2 * uwdot(evec[i][:,n], ct_mat[i][:,n])
end
return res
end
function dec(mat_elem::Vector{uwreal}, plat::Array{Int64}, mass::uwreal, mu::Vector{Float64}; pl::Bool=true, wilson::Bool=false)
aux = plat_av(mat_elem, plat)
if pl
uwerr(aux)
x = 1:length(mat_elem)
y = value.(mat_elem)
dy = err.(mat_elem)
v = value(aux)
e = err(aux)
figure()
errorbar(x, y, yerr=dy, fmt="*", color="red")
fill_between(plat[1]:plat[2], v-e, v+e, color="green")
ylim(v - 10*e, v + 10*e)
display(gcf())
close()
end
if !wilson
return sqrt(2/ (mass)^3) * sum(mu)* aux
else
return sqrt(2/mass) * aux
end
end
function vec_dec(mat_elem::Vector{uwreal}, plat::Array{Int64}, mass::uwreal; pl::Bool=true)
aux = plat_av(mat_elem, plat)
if pl
uwerr(aux)
x = 1:length(mat_elem)
y = value.(mat_elem)
dy = err.(mat_elem)
v = value(aux)
e = err(aux)
figure()
ylim(v - 10*e, v + 10*e)
errorbar(x, y, yerr=dy, fmt="*", color="red")
fill_between(plat[1]:plat[2], v-e, v+e, color="green")
display(gcf())
end
return sqrt(2 / mass) * aux #aux / mass
end
function param_av(fit_masses::Vector{uwreal})
w = 1 ./ err.(fit_masses).^2
av = sum(w .* fit_masses) / sum(w)
return av
end
function get_meff(obs::Vector{juobs.Corr}, ens::EnsInfo; pl::Bool=false, wpm::Union{Nothing, Dict{String,Vector{Float64}}}=nothing)
mu_list = getfield.(obs, :mu)
plat = select_plateau(ens, mu_list)
return meff.(obs, plat, pl=pl, wpm=wpm)
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(ensinfo::EnsInfo)
ens = ensinfo.id
f = readdlm(path_plat_meff)
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
head_ = String.(f[del, 1])
plat_ = Int64.(f[del, 2:3])
plat = Vector{Int64}(undef, 2)
n = findfirst(x-> occursin("t0", x), head_)
plat[1] = plat_[n,1]
plat[2] = plat_[n,2]
return plat
end
function select_plateau(ensinf::EnsInfo, mu_list, path_plat)
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
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 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 aic_weigth(cat::Cat; cut::Bool=false)
w_aux = exp.(-1/2 * (getfield(cat, :aic)))
w = w_aux ./ sum(w_aux)
return w
end
function aic_syst(cat::Cat)
w = aic_weigth(cat)
aux1 = sum( w .* getfield(cat,:fit_res).^2)
aux2 = sum( w .* getfield(cat,:fit_res))^2
return sqrt(abs(value(aux1 - aux2)))
end
function aic_model_ave(cat::Cat)
w = aic_weigth(cat)
res = sum( w.* getfield(cat, :fit_res))
uwerr(res)
return res
end
function category_av(cat_arr::Vector{Cat})
best_res = aic_model_ave.(cat_arr)
println(best_res)
syst = aic_syst.(cat_arr)
tot_err = sqrt.(err.(best_res).^2 .+ syst.^2)
println(tot_err)
w_aux = 1 ./ (tot_err)
w = w_aux ./ sum(w_aux)
println(w)
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 syst_av(syst_err::Vector{Float64})
aux = sum(syst_err.^2)
return sqrt(aux)
end
src/gevp_types.jl 0 → 100644
View file @ 838f8293
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 MatInfo
ensinfo::EnsInfo
mu::Vector{Float64}
mat_list::Vector{Matrix{uwreal}}
y0::Int64
function MatInfo(_ensinfo::EnsInfo, _mat_list::Array{Array{T,2} where T,1}, _mu::Vector{Float64}, _y0::Int64)
a = new()
a.ensinfo = _ensinfo
a.mu = _mu
a.mat_list = _mat_list
a.y0 = _y0
return a
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
dof::Vector{Int64}
function Cat(_obs::Vector{uwreal}, _xtofit::Array{uwreal}, _phi_ph::uwreal, _info::String)
a = new()
a.obs = _obs
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)
a.dof = Vector{Int64}(undef, 0)
return a
end
end
\ No newline at end of file
src/plat_dec.txt 0 → 100644
View file @ 838f8293
#H101
ll 58 75
lh 58 70
hh 60 72
t0 20 70
#H102r001
ll 60 85
ls 60 85
lh 60 85
ss 60 85
sh 60 85
hh 60 85
t0 20 80
#H102r002
ll 58 73
ls 60 75
lh 62 80
ss 60 80
sh 62 80
hh 63 80
t0 20 80
#H400
ll 60 75
lh 61 75
hh 66 78
t0 20 80
#N300
ll 80 110
lh 85 115
hh 90 112
t0 20 105
#N200
ll 85 102
ls 82 105
lh 82 99
ss 82 105
sh 83 105
hh 87 100
t0 20 105
#N202
ll 76 105
lh 80 98
hh 90 110
t0 20 110
#N203
ll 85 105
ls 85 105
lh 82 105
ss 85 105
sh 83 110
hh 90 104
t0 20 110
#J303
ll 120 140
ls 123 140
lh 120 140
ss 123 140
sh 123 140
hh 133 170
t0 25 170
#end
src/plat_dec_vec.txt 0 → 100644
View file @ 838f8293
#H101
ll 56 70
lh 59 68
hh 68 75
t0 20 70
#H102r001
ll 60 80
ls 60 85
lh 57 70
ss 57 67
sh 60 70
hh 63 79
t0 20 80
#H102r002
ll 60 80
ls 60 85
lh 57 70
ss 57 67
sh 60 70
hh 65 79
t0 20 80
#H400
ll 56 68
lh 60 73
hh 65 78
t0 20 80
#N300
ll 80 110
lh 83 98
hh 90 108
t0 20 105
#N200
ll 75 90
ls 79 95
lh 78 95
ss 79 95
sh 79 95
hh 89 105
t0 20 105
#N202
ll 78 95
lh 78 98
hh 88 110
t0 20 110
#N203
ll 77 95
ls 76 95
lh 78 95
ss 76 95
sh 80 100
hh 89 102
t0 20 110
#J303
ll 116 130
ls 118 135
lh 120 135
ss 118 140
sh 123 140
hh 130 155
t0 25 170
#end
src/plat_meff.txt 0 → 100644
View file @ 838f8293
#H101
ll 66 76
lh 66 80
hh 64 85
t0 20 80
#H102r001
ll 60 85
ls 60 85
lh 60 85
ss 60 85
sh 60 85
hh 68 85
t0 20 80
#H102r002
ll 65 72
ls 67 74
lh 67 75
ss 67 74
sh 68 81
hh 68 85
t0 20 80
#H400
ll 65 88
lh 67 80
hh 67 85
t0 20 80
#N300
ll 80 110
lh 92 115
hh 92 115
t0 20 105
#N200
ll 85 105
ls 85 105
lh 88 99
ss 85 105
sh 89 105
hh 87 115
t0 20 105
#N202
ll 80 110
lh 87 100
hh 90 115
t0 20 110
#N203
ll 85 105
ls 85 105
lh 88 100
ss 85 105
sh 88 110
hh 88 110
t0 20 110
#J303
ll 120 140
ls 123 140
lh 125 145
ss 123 140
sh 130 145
hh 133 175
t0 25 170
#end
src/plat_meff_vec.txt 0 → 100644
View file @ 838f8293
#H101
ll 66 73
lh 66 71
hh 68 80
t0 20 80
#H102r001
ll 60 85
ls 60 85
lh 60 85
ss 60 85
sh 60 85
hh 60 85
t0 20 80
#H102r002
ll 60 69
ls 60 68
lh 63 68
ss 60 67
sh 68 73
hh 70 85
t0 20 80
#H400
ll 60 68
lh 68 74
hh 68 80
t0 20 80
#N300
ll 80 110
lh 88 100
hh 95 115
t0 20 105
#N200
ll 82 95
ls 85 95
lh 85 99
ss 85 95
sh 86 95
hh 92 110
t0 20 105
#N202
ll 84 95
lh 87 97
hh 92 110
t0 20 110
#N203
ll 85 95
ls 85 95
lh 85 95
ss 85 95
sh 86 100
hh 89 110
t0 20 110
#J303
ll 120 140
ls 123 140
lh 123 140
ss 123 140
sh 123 145
hh 133 165
t0 25 170
#end
src/tools.jl
View file @ 838f8293
......@@ -135,7 +135,6 @@ function get_hh(mu_list, obs::Vector{uwreal}, deg::Bool)
return obs_hh
end
function comp_mat(ensinfo::EnsInfo, tot_obs::Vector{juobs.Corr}, tau::Int64)
#wanted_corr = relevant_corr(ensinfo, tot_obs)
mat_info = Vector{MatInfo}(undef, length(tot_obs))
for i = 1:length(tot_obs)
a11 = tot_obs[i].obs[1:end-1-2*tau]
......@@ -146,6 +145,21 @@ function comp_mat(ensinfo::EnsInfo, tot_obs::Vector{juobs.Corr}, tau::Int64)
end
return mat_info
end
function comp_mat_multigamma(ensinfo::EnsInfo,tot11::Vector{juobs.Corr}, tot12::Vector{juobs.Corr}, tot22::Vector{juobs.Corr})
mu = getfield.(tot11,:mu)
y0 = getfield.(tot11,:y0) .+1
mat_info = Vector{MatInfo}(undef, length(mu))
for i = 1:length(mu)
a11 = tot11[i].obs[y0[i]:end-1]
a12 = tot12[i].obs[y0[i]:end-1]
a22 = tot22[i].obs[y0[i]:end-1]
#res[i] = infoMatt(a11, a12, a22)
aux = juobs.get_matrix([a11, a22], [a12])
mat_info[i] = MatInfo(ensinfo, aux, mu[i], y0[i])
end
return mat_info
end
function comp_mat_lukas(ensinfo::EnsInfo, tot_obs::Vector{juobs.Corr}, dim::Int64)
mat_info = Vector{MatInfo}(undef, length(tot_obs))
for i in 1:length(tot_obs)
......@@ -160,8 +174,35 @@ function comp_mat_lukas(ensinfo::EnsInfo, tot_obs::Vector{juobs.Corr}, dim::Int6
end
return mat_info
end
function gevp_mass(mat_obs::Vector{MatInfo}; t0::Int64=2, delta_t_in::Array{Int64}=[2,5], wpm::Union{Nothing, Dict{Int64,Vector{Float64}}}=nothing)
y0s = getfield.(mat_obs, :y0)
function gevp_mass(mat_obs::Vector{MatInfo}; t0::Int64=2,pl::Bool=true, wpm::Union{Nothing, Dict{Int64,Vector{Float64}}}=nothing)
y0s = getfield.(mat_obs, :y0)
mu_list = getfield.(mat_obs, :mu)
plat_en0 = Vector{uwreal}(undef, length(mat_obs))
plat_en1 = Vector{uwreal}(undef, length(mat_obs))
for i in 1:length(mat_obs)
mat = getfield(mat_obs[i], :mat_list)[y0s[i]+1:end-1] #matrices to diagonalise for a given sector [i]
evals = getall_eigvals(mat, t0)
en = energies(evals, wpm=wpm) #ground and excited state energy
plat = select_plateau(mat_obs[i].ensinfo, mu_list, path_plat_meff)[i] .- y0s[i]
println(plat)
plat_en0[i] = plat_av(en[1], plat)
if pl
uwerr(plat_en0[i])
xx = collect(1:length(en[1]))
v = value(plat_en0[i])
e = err(plat_en0[i])
errorbar(xx, value.(en[1]), err.(en[1]), fmt="s")
fill_between(plat[1]:plat[2], v-e, v+e )
ylim(v*0.95, v*1.05)
display(gcf())
close()
end
end
return plat_en0
end
# function used in lattice conf 2021 with fits to the effective energies
function gevp_mass_old(mat_obs::Vector{MatInfo}; t0::Int64=2, delta_t_in::Array{Int64}=[2,5], wpm::Union{Nothing, Dict{Int64,Vector{Float64}}}=nothing)
y0s = getfield.(mat_obs, :y0)
mu_list = getfield.(mat_obs, :mu)
plat_en0 = Vector{uwreal}(undef, length(mat_obs))
plat_en1 = Vector{uwreal}(undef, length(mat_obs))
......@@ -175,7 +216,7 @@ function gevp_mass(mat_obs::Vector{MatInfo}; t0::Int64=2, delta_t_in::Array{Int6
temp = fit_mass(en[1],t_in=t, wpm=wpm)
push!(aux_mass, temp)
catch
plat = select_plateau(mat_obs[i].ensinfo, mu_list)[end-2] .- y0s[i]
plat = select_plateau(mat_obs[i].ensinfo, mu_list, path_plat_meff)[end-2] .- y0s[i]
println(plat)
temp = plat_av(en[1], plat)
uwerr(temp)
......@@ -204,12 +245,12 @@ function gevp_dec(mat_obs::Vector{MatInfo}, mass::Vector{uwreal}; t0::Int64=2, w
mat = getfield(mat_obs[i], :mat_list)[y0s[i]+2:end-1] #matrices to diagonalise for a given sector [i]
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]
plat = select_plateau(mat_obs[i].ensinfo, mu_list, path_plat_dec)[i] .- y0s[i]
try
if pseudo
plat_dec0[i] = -dec(elem, plat, mass[i], mu_list[i], pl=pl, wilson=wilson)
plat_dec0[i] = dec(elem, plat, mass[i], mu_list[i], pl=pl, wilson=wilson)
else
plat_dec0[i] = -vec_dec(elem, plat, mass[i], pl=pl)
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)
......@@ -277,6 +318,7 @@ function dec(mat_elem::Vector{uwreal}, plat::Array{Int64}, mass::uwreal, mu::Vec
if !wilson
return sqrt(2/ (mass)^3) * sum(mu)* aux
else
println("ciao its wilson")
return sqrt(2/mass) * aux
end
end
......@@ -328,7 +370,7 @@ function select_plateau(ensinfo::EnsInfo)
plat[2] = plat_[n,2]
return plat
end
function select_plateau(ensinf::EnsInfo, mu_list)
function select_plateau(ensinf::EnsInfo, mu_list, path_plat)
ens =ensinf.id
deg = ensinf.deg
mul, mus, muh = get_mu(mu_list, deg)
......@@ -425,4 +467,16 @@ function syst_av(syst_err::Vector{Float64})
aux = sum(syst_err.^2)
return sqrt(aux)
end
function read_BDIO(path::String, ens_id::String)
file = filter(x->occursin(ens_id, x), readdir(path, join=true))
println(file)
r = Vector{uwreal}(undef, 0)
fb = BDIO_open(file[1], "r")
BDIO_seek!(fb)
push!(r, read_uwreal(fb))
while BDIO_seek!(fb, 2)
push!(r, read_uwreal(fb))
end
BDIO_close!(fb)
return r
end
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