Skip to content

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

Analysis updated with new HAWK ensembles. New functional forms, general updates

parent 9fdef14c
Hide whitespace changes
Inline Side-by-side
Showing with 1276 additions and 207 deletions
plat/plat_tm_mps.txt
View file @ e56a3bed
......@@ -32,6 +32,14 @@ ll 65 80
lh 67 80
hh 67 85
t0 20 80
#D450
ll 79 100
ls 85 105
lh 78 88
ss 85 105
sh 85 100
hh 90 120
t0 25 100
#N300
ll 78 105
lh 92 108
......@@ -90,6 +98,14 @@ ss 123 140
sh 130 145
hh 133 175
t0 25 170
#E250
ll 120 148
ls 123 140
lh 125 145
ss 123 140
sh 130 145
hh 133 175
t0 25 170
#J500
ll 120 148
ls 123 140
......
src/gevp analysis/const.jl
View file @ e56a3bed
......@@ -6,31 +6,35 @@ ens_db = Dict(
"H105r002" => [32, 3.4, false, 0.12151, 2, nothing],
"H101" => [32, 3.4, true, 0.17979, 2, nothing],
"H400" => [32, 3.46, true, 0.16345, 1, nothing],
"D450" => [64, 3.46, false, 0.0, 1, nothing],
"N200" => [48, 3.55, false, 0.09222, 1, nothing],
"N202" => [48, 3.55, true, 0.13407, 2, nothing],
"N203" => [48, 3.55, false, 0.11224, 1, nothing],
"D200" => [64, 3.55, false, 0.06611, 2, 100],
"N300" => [48, 3.70, true, 0.10630, 1, nothing],
"N300r002" => [48, 3.70, true, 0.10630, 1, nothing],
"N302" => [48, 3.70, false, 0.10630, 1, nothing],
"J303" => [64, 3.70, false, 0.06514, 2, nothing],
"E300" => [96, 3.70, false, 0.00000, 1, nothing],
"E250" => [96, 3.55, false, 0.00000, 1, nothing],
"J500" => [64, 3.85, true, 0.00000, 2, nothing],
"J501" => [64, 3.85, false, 0.00000, 2, nothing],
"J501" => [64, 3.85, false, 0.00000, 1, nothing],
)
trunc_db = Dict(
"H102r002" => nothing,
#"H102r001" => nothing,
"H105r002" => nothing,
"H101" => [1001,1009],
"H101" => nothing,
"H400" => nothing,
"D450" => nothing,
"N200" => nothing,
"N202" => nothing,
"N203" => nothing,
"D200" => 645, #1000
"N300" => 1279,
"D200" => nothing, #1000
"N300r002" => 1521,
"N302" => nothing,
"J303" => 721,
"E300" => nothing,
"J303" => nothing,
"E300" => 1136,
"E250" => 751,
"J500" => [745,649], #[751, 655],
"J501" => nothing
......@@ -46,17 +50,18 @@ ens_nms = Dict(
"N202" => 899,
"N203" => 1543,
"D200" => 2001,
"N300" => 1540,
"N300r002" => 1540,
"N302" => 2201,
"J303" => 1073,
"E300" => 1139
"E300" => 1139,
"J500" => 1875
)
#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)
const M_error = [0.05, 0.05, 0.07, 0.4, 1.2, 0.011]
const M_values = [1869.58, 2010.26, 1968.34, 2112.2, 2983.9, 3096.916] #MD, MD*, MDs, MDs*, \eta_c, J/\psi (MeV)
const M_error = [0.09, 0.05, 0.07, 0.4, 0.9, 0.011]
#1802.05243
const b_values = [3.40, 3.46, 3.55, 3.70, 3.85]
const b_values2 = [3.40, 3.46, 3.55, 3.70]
......@@ -67,8 +72,8 @@ const ZM_tm_error = [42, 33, 42, 48, 55] .* 1e-4
#1608.08900
const t0_data = [2.86, 3.659, 5.164, 8.595, 14.040]
const t0_error = [11, 16, 18, 29, 49] .* 1e-3
const t0_ph_value = #=[0.415] =# [0.4137]
const t0_ph_error = ones(1,1) .* #=4e-3 =# 3.6e-3
const t0_ph_value = #=[0.415] =# [0.1445*sqrt(8)]
const t0_ph_error = ones(1,1) .* #=4e-3 =# 0.00058*sqrt(8)
# 1808.09236
const ZA_data = [0.75642, 0.76169, 0.76979, 0.78378, 0.79667]
const ZA_err = [72, 93, 43, 47, 47] .*1e-5
......
src/gevp analysis/func_comb.jl
View file @ e56a3bed
......@@ -88,6 +88,8 @@ append!(f_aux, f_non_lin_aux)
########################################
########## DEC. CONST. ##########
########################################
# In the log functions I substituted p[4] with (1. .+ 3 .* 0.52)./(64. .* π.^2 .* x[:,4]) to
# test a fixed coefficient in front of the logs.
phi2_sym = 0.744666
##FUNCTIONS
......@@ -104,9 +106,9 @@ 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
f_basemodel1(x,p) = p[1] .+ p[2] .* x[:,2] .+ p[3] .* (1.0 ./ sqrt.(x[:,3])) .+ p[4] .* x[:,1]
f_basemodel2(x,p) = p[1] .+ p[2] .* (3 .* phi2_sym - 2 .* x[:,2]) .+ p[3] .* (1.0 ./ sqrt.(x[:,3])) .+ p[5] .* x[:,1]
#CHIRAL LOGS
f_basemodel1_clog(x,p) = p[1] .+ p[2] .* x[:,2] .+ p[3] .* (1.0 ./ sqrt.(x[:,3])) .+ p[5] .* x[:,1] .+ p[4] .* (3 .* x[:,2] .*log.(x[:,2]) .+
f_basemodel1_clog(x,p) = p[1] .+ p[2] .* x[:,2] .+ p[3] .* (1.0 ./ sqrt.(x[:,3])) .+ p[5] .* x[:,1] .- (1. .+ 3 .* p[4])./(64. .* π.^2 .* x[:,4]) .* (3 .* x[:,2] .*log.(x[:,2]) .+
(3 .*phi2_sym .- x[:,2]) .* log.(1.5 .*phi2_sym .- 0.5*x[:,2]) .+ (2 .* phi2_sym .- x[:,2]) ./ 3 .* log.(2 .*phi2_sym .- x[:,2]))
f_basemodel2_clog(x,p) = p[1] .+ p[2] .* (3 .* phi2_sym - 2 .* x[:,2]) .+ p[3] .* (1.0 ./ sqrt.(x[:,3])) .+ p[6] .* x[:,1] .+ 4 .* p[4] .*(
f_basemodel2_clog(x,p) = p[1] .+ p[2] .* (3 .* phi2_sym - 2 .* x[:,2]) .+ p[3] .* (1.0 ./ sqrt.(x[:,3])) .+ p[6] .* x[:,1] .- 4 .* (1. .+ 3 .* p[4])./(64. .* π.^2 .* x[:,4]) .*(
(1.5 .* phi2_sym .- 0.5.*x[:, 2]) .* log.(1.5 .* phi2_sym .- 0.5.*x[:, 2]) .+ (2 .* phi2_sym .- x[:,2]) ./ 3 .* log.((2 .* phi2_sym .- x[:,2])))
f_func_list = [f_a2phi2, f_a2phih2, f_a4phih4, f_a4phih2, f_a4cutoff]
......@@ -145,9 +147,9 @@ f_f_aux2_clog[1] = Vector{Function}(undef, 1)
f_f_lin1_clog[1][1] = (x,p) -> f_basemodel1_clog(x, p)
f_f_lin2_clog[1][1] = (x,p) -> f_basemodel2_clog(x, p)
f_f_aux1_clog[1][1] = (x,p) -> p[1] + p[2] * x[:,2] + p[3] * (1.0 ./ sqrt.(x[:,3])) + p[5] * x[:,1] + p[4] * (3 * x[:,2] .*log.(x[:,2]) +
f_f_aux1_clog[1][1] = (x,p) -> p[1] + p[2] * x[:,2] + p[3] * (1.0 ./ sqrt.(x[:,3])) + p[5] * x[:,1] .- (1. .+ 3 .* p[4])./(64. .* π.^2 .* x[:,4]) .* (3 * x[:,2] .*log.(x[:,2]) +
(3 *phi2_sym - x[:,2]) .* log.(1.5 *phi2_sym - 0.5*x[:,2]) + (2 * phi2_sym - x[:,2]) / 3 .* log.(2 *phi2_sym - x[:,2]))
f_f_aux2_clog[1][1] = (x,p) -> p[1] + p[2] * (3 * phi2_sym - 2 * x[:,2]) + p[3] * (1.0 ./ sqrt.(x[:,3])) + p[6] * x[:,1] + 4 * p[4] *(
f_f_aux2_clog[1][1] = (x,p) -> p[1] + p[2] * (3 * phi2_sym - 2 * x[:,2]) + p[3] * (1.0 ./ sqrt.(x[:,3])) + p[6] * x[:,1] .- 4 .* (1. .+ 3 .* p[4])./(64. .* π.^2 .* x[:,4]) .*(
(1.5 * phi2_sym - 0.5*x[:, 2]) .* log.(1.5 * phi2_sym - 0.5*x[:, 2]) + (2 * phi2_sym - x[:,2]) / 3 .* log.((2 * phi2_sym - x[:,2])))
## loop over combinations
......@@ -190,14 +192,14 @@ f_n_param_comb_clog = collect(6:2:f_npar_comb+6)
########################################
######## RATIO DEC. CONST. ############
########################################
# Testing HQET and TAYLOR functional forms for (fDs*sqrt{mDs)}/(fD*sqrt{Md})
# HQET and TAYLOR functional forms for (fDs*sqrt{mDs)}/(fD*sqrt{Md})
# from Gregorio's notes https://cloud.ift.uam-csic.es/s/9PGEbX5sDzcN5Xq#pdfviewer
# base model from HQET
# before p[1] = \phi_f^2 in notes. Treated as a parameter but using fpk data would be better
# now p[1]=x[:,5] fpik data from AS and p[2]=p[1]
r_basemodel_hqet(x,p) = 1. .- (1. .+ 3 .* 0.52)./(64. .* π.^2 .* x[:,5]) .*
r_basemodel_hqet(x,p) = 1. .- (1. .+ 3 .* p[3])./(64. .* π.^2 .* x[:,5]) .*
((3 .*phi2_sym .- x[:,2]) .* log.(1.5 .* phi2_sym .- 0.5 .* x[:,2]) .+ #2L_k chiral log
(2 .* phi2_sym .- x[:,2]) .* log.(2 .* phi2_sym .- x[:,2]) .- # L_\eta chiral log
3 .* x[:,2] .* log.(x[:,2]) ) .+ # 3L_\pi chiral log
......@@ -212,23 +214,22 @@ r_basemodel_continuum(x,p) = 1. .+ p[2] .* x[:,1] .* (3 .* phi2_sym .- 3 .* x[:,
# higher order lattice spacing dependence
r_phih_a_continuum(x) = x[:,1] .* (3 .* phi2_sym .- 3 .* x[:,2]) .* x[:,3].^2
r_cont_func_list = [r_phih_a_continuum]
# r_cont_func_label = [L"$p_a^{(3)}$"]
r_cont_func_label = ["pa(3)"]
r_cont_extra_par = length(r_cont_func_list)
r_cont_func_map = [Bool.([i]) for i=0:1]
r_hqet_func_list = [r_chiphi2_hqet, r_chiphih_hqet ] #, r_phih_a_continuum]
# r_hqet_func_label = [L"$p\chi^{(2)}$", L"$p_{\chi}^{(4)}$" ] #, L"$p_a^{(3)}$" ]
r_hqet_func_label = ["pχ(2)", "pχ(4)" ] #, L"$p_a^{(3)}$" ]
r_hqet_func_label = ["p(2)", "p(4)" ] #, L"$p_a^{(3)}$" ]
r_hqet_func_map = [Bool.([i,j]) for i=0:1 for j=0:1 ]#for k in 0:1 ]
r_hqet_base_par = 2 # 1 from hqet base model, 1 from LO lattice spacing dependence
r_hqet_base_par = 3 # 2 from hqet base model, 1 from LO lattice spacing dependence
r_hqet_extra_par = length(r_hqet_func_list) # number of extra parameters on top of base model
r_hqet_models = Vector{Function}(undef,1) #r_hqet_base_par+r_hqet_extra_par)
r_hqet_models[1] = (x,p) -> r_basemodel_hqet(x,p) .* r_basemodel_continuum(x,p)
r_hqet_param = [2]
r_hqet_param = [3]
r_hqet_labels = [["hqet"]]
for extrapar = 1:r_hqet_extra_par
......@@ -324,6 +325,126 @@ r_tot_labels = vcat(r_hqet_labels, r_taylor_labels)
########################################
######## DEC. CONST. HQET FORMS ############
########################################
# HQET functional forms for fDs*sqrt{mDs} and fD*sqrt{Md}
# from Gregorio's notes https://cloud.ift.uam-csic.es/s/9PGEbX5sDzcN5Xq#pdfviewer
# CONTINUUM DEPENDENCE
phid_basemodel(x,p) = p[1] .- (1. .+ 3 .* p[2])./(64. .* π.^2 .* x[:,4]) .*
(
3 .* x[:,2] .* log.(x[:,2]) .+ # 3L_pi
(3 .*phi2_sym .- x[:,2]) .* log.(1.5 .* phi2_sym .- 0.5 .* x[:,2]) .+ #2L_k
1. ./3. .* (2 .* phi2_sym .- x[:,2]) .* log.(2 .* phi2_sym .- x[:,2]) # 1/3L_\eta
) .+ p[3] .* (4 .* x[:,2]) ./ (x[:,4])
phids_basemodel(x,p) = p[1] .- (1. .+ 3 .* p[2])./(64. .* π.^2 .* x[:,4]) .*
(
2 .* (3 .*phi2_sym .- x[:,2]) .* log.(1.5 .* phi2_sym .- 0.5 .* x[:,2]) .+ #2L_k
4. ./3. .* (2 .* phi2_sym .- x[:,2]) .* log.(2 .* phi2_sym .- x[:,2]) # 1/3L_\eta
) .+ p[3] .* (4 .* (3 .* phi2_sym .- 2 .*x[:,2])) ./ ( x[:,4])
# higher order terms
phid_chiphi2(x) = (4 .* x[:,2].^2) ./ (x[:,4].^2)
phid_chiphih(x) = (4 .+ x[:,2]) ./ (x[:,4] .* x[:,3])
phids_chiphi2(x) = (4 .* (3 .* phi2_sym .- 2 .*x[:,2]) .* x[:,2]) ./ ( x[:,4].^2)
phids_chiphih(x) = (4 .* (3 .* phi2_sym .- 2 .*x[:,2]) ) ./ ( x[:,4] .* x[:,3])
phid_ds_phih(x) = 1 ./ x[:,3]
d_hqet_func_list = [phid_chiphi2, phid_chiphih, phid_ds_phih]
ds_hqet_func_list = [phids_chiphi2, phids_chiphih, phid_ds_phih]
d_ds_hqet_func_label = ["p[chi](2)", "p[chi](4)", "ph(1)"]
d_ds_hqet_extra_param = length(d_hqet_func_list)
d_ds_hqet_func_map = [Bool.([i,j,k]) for i=0:1 for j=0:1 for k=0:1]
d_ds_base_par = 4 # 3 from cont dependence 1 from a^2 lattice spacing.
# LATTICE SPACING DEPENDENCE
phid_ds_base_cont(x,p) = 1. .+ p[4] .* x[:,1]
#higher order terms
phid_cont_phi2(x) = x[:,1] .* x[:,2]
phids_cont_phi2(x) = x[:,1] .* 2. .*( 3 ./2 .* phi2_sym .- x[:,2])
phid_cont_phi2phih(x) = x[:,1] .* x[:,2] .* x[:,3].^2
phids_cont_phi2phih(x) = x[:,1] .* 2. .* ( 3 ./2 .* phi2_sym .- x[:,2]) .* x[:,3].^2
phid_ds_cont_phih(x) = x[:,1] .* x[:,3].^2
d_cont_func_list = [phid_cont_phi2, phid_cont_phi2phih, phid_ds_cont_phih]
ds_cont_func_list = [phids_cont_phi2, phids_cont_phi2phih, phid_ds_cont_phih]
d_ds_cont_func_label = ["pa(1)", "ps(3)", "pa(2)"]
d_ds_cont_extra_param = length(d_cont_func_list)
d_ds_cont_func_map = [Bool.([i,j,k]) for i=0:1 for j=0:1 for k=0:1]
# building all functions
phid_tot_models = Vector{Function}(undef, 1)
phid_tot_models[1] = (x,p) -> phid_basemodel(x,p) .* phid_ds_base_cont(x,p)
phids_tot_models = Vector{Function}(undef, 1)
phids_tot_models[1] = (x,p) -> phids_basemodel(x,p) .* phid_ds_base_cont(x,p)
d_ds_param = [4] # n. of parameters for the leading order expression
d_ds_labels = [["hqet"]] # label of the leading order expression
for extrapar = 1:d_ds_hqet_extra_param
idx_f_expar = filter(x->sum(x)==extrapar, d_ds_hqet_func_map)
tmp_fd_hqet = Vector{Function}()
tmp_fds_hqet = Vector{Function}()
tmp_param = []
tmp_label = []
for ii in idx_f_expar
aux_fd = (x,p) -> phid_basemodel(x,p) .+ sum(p[d_ds_base_par+k] .* d_hqet_func_list[ii][k](x) for k in eachindex(d_hqet_func_list[ii]))
aux_fds = (x,p) -> phids_basemodel(x,p) .+ sum(p[d_ds_base_par+k] .* ds_hqet_func_list[ii][k](x) for k in eachindex(ds_hqet_func_list[ii]))
push!(tmp_fd_hqet, aux_fd)
push!(tmp_fds_hqet, aux_fds)
push!(tmp_param, d_ds_base_par + extrapar)
push!(tmp_label, d_ds_hqet_func_label[ii])
end
# add lattice spacing dependence
for jj in eachindex(tmp_fd_hqet)
aux_fd = tmp_fd_hqet[jj]
aux_fds = tmp_fds_hqet[jj]
aux1_fd = (x,p) -> aux_fd(x,p) .* phid_ds_base_cont(x,p)
aux1_fds = (x,p) -> aux_fds(x,p) .* phid_ds_base_cont(x,p)
push!(phid_tot_models, aux1_fd)
push!(phids_tot_models, aux1_fds)
push!(d_ds_labels, tmp_label[jj])
push!(d_ds_param, tmp_param[jj])
for contpar = 1:d_ds_cont_extra_param
idx_cont_expar = filter(x->sum(x)==contpar, d_ds_cont_func_map)
for ii in idx_cont_expar
aux_fd_cont = (x,p) -> aux_fd(x,p) .* (phid_ds_base_cont(x,p) .+ sum(p[d_ds_base_par+extrapar+k] .* d_cont_func_list[ii][k](x) for k in eachindex(d_cont_func_list[ii])))
aux_fds_cont = (x,p) -> aux_fds(x,p) .* (phid_ds_base_cont(x,p) .+ sum(p[d_ds_base_par+extrapar+k] .* ds_cont_func_list[ii][k](x) for k in eachindex(ds_cont_func_list[ii])))
push!(phid_tot_models, aux_fd_cont)
push!(phids_tot_models, aux_fds_cont)
push!(d_ds_labels, vcat(tmp_label[jj], d_ds_cont_func_label[ii]))
push!(d_ds_param, tmp_param[jj]+contpar)
end
end
end
end
# at the end of this loop all HQET + continuum combiations are taken into account for phid and phids and stored in
# phid_tot_models and phids_tot_models respectively. Labels and parameters are given in d_ds_labels and d_ds_param
##
# Operator overloading to call functions with uwreal parameters
......
src/gevp analysis/gevp_fit.jl
View file @ e56a3bed
......@@ -13,22 +13,28 @@ using juobs, BDIO, DelimitedFiles, ADerrors, LaTeXStrings, PyPlot, Statistics
#============= SET UP VARIABLES ===========#
rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["text.usetex"] = false
rcParams["text.usetex"] = true
rcParams["mathtext.fontset"] = "cm"
rcParams["font.size"] =12
rcParams["axes.labelsize"] =20
rcParams["axes.titlesize"] = 18
plt.rc("text", usetex=false) # set to true once you install latex
plt.rc("text", usetex=true) # set to true once you install latex
#rcParams.update(PyPlot.matplotlib.mpl.rcParamsDefault)
#================== PATHS ==================#
# data by ale (carfull with the read_BDIO that you include)
const path_bdio = "/Users/alessandroconigli/Lattice/data/bdio_charm/LAT22"
#const path_bdio = "/Users/alessandroconigli/Lattice/data/bdio_charm/from_gevp/tm_shifted/ale"
# path to data
const path_data_w = "/Users/alessandroconigli/Lattice/data/aux_obs_data/wilson" # required for t0 at finite a
const path_rw = "/Users/alessandroconigli/Lattice/data/aux_obs_data/rwf" # reweighting factors
# path to bdio
const path_bdio = "/Users/alessandroconigli/Lattice/data/bdio_charm/LAT22/"
const path_bdio_tm = joinpath(path_bdio, "tm_shifted")
# const path_bdio_md = joinpath(path_bdio, "md") # path to aux obs (md) comptued by me
const path_bdio_md = "/Users/alessandroconigli/Lattice/data/bdio_charm/AS_ligth_obs/" # path to delta m computed by AS
const path_bdio_t0 = "/Users/alessandroconigli/Lattice/data/bdio_charm/AS_ligth_obs/t0_combined.bdio" # path to t0_ph computed by AS
# path to t0 plateau
const path_plat_w = "/Users/alessandroconigli/Lattice/data/aux_obs_data/plat_wilson.txt"
# data by javier (carful with the read_BDIO that you include. do not use the one in this folder)
......@@ -37,7 +43,7 @@ const path_bdio_md = "/Users/alessandroconigli/Lattice/data/bdio_charm/AS_ligth_
# const path_bdio_md = joinpath(path_bdio, "md")
# by reading javier data you have deactivated the reading of vector decay constants
const path_plot = "/Users/alessandroconigli/MyDrive/phd/analysis/charm_gevp_2022/plots"
const path_plot = "/Users/alessandroconigli/MyDrive/phd/analysis/charm_gevp_2022/plots" #/plots for Lat22 ensembles # /Thesis for thesis ensembles
#========= INCLUDES ==========#
include("./const.jl")
......@@ -66,15 +72,17 @@ end
wpmm = Dict{String, Vector{Float64}}()
wpmm["H101"] = [-1.0, 2.0, -1.0, -1.0]
wpmm["H102r002"] = [-1.0, 2.0, -1.0, -1.0]
wpmm["H105r002"] = [-1.0, 2.0, -1.0, -1.0]
wpmm["H105"] = [-1.0, 2.0, -1.0, -1.0]
wpmm["H400"] = [-1.0, 2.0, -1.0, -1.0]
wpmm["N202"] = [-1.0, 2.0, -1.0, -1.0]
wpmm["N200"] = [-1.0, 2.0, -1.0, -1.0]
wpmm["N203"] = [-1.0, 2.0, -1.0, -1.0]
wpmm["N300r002"] = [-1.0, 2.0, -1.0, -1.0]
wpmm["N302"] = [-1.0, 2.0, -1.0, -1.0]
wpmm["N300"] = [-1.0, 2.0, -1.0, -1.0]
wpmm["J303"] = [-1.0, 2.0, -1.0, -1.0]
wpmm["D200"] = [-1.0, 2.0, -1.0, -1.0]
wpmm["D200_corr"] = [1.0, -1.0, -1.0, -1.0]
wpmm["D200"] = [4.0, -2.0, -1.0, -1.0]
wpmm["D200_corr"] = [1.0, -2.0, -1.0, -1.0]
# READ MASSES + DECAY CONSTANTS
......@@ -103,12 +111,18 @@ phi4 = Vector{uwreal}(undef, Nens)
fpik = Vector{uwreal}(undef, Nens)
for (k, ens) in enumerate(enslist)
println(ens)
p = joinpath(path_bdio_tm, ens)
mlh[k] = read_BDIO(p, "tm", "mlh")
msh[k] = read_BDIO(p, "tm", "msh")
mhh[k] = read_BDIO(p, "tm", "mhh")
mll = read_BDIO(p, "tm", "mll")[1]
mls = read_BDIO(p, "tm", "mls")[1]
if ensinfo[k].deg
mll.mean = 0.859837
mls.mean = 0.859837
end
mlh_v[k] = read_BDIO(p, "tm", "mlh_star")
msh_v[k] = read_BDIO(p, "tm", "msh_star")
......@@ -124,6 +138,20 @@ for (k, ens) in enumerate(enslist)
#fpik[k] = read_BDIO(p, "tm", "fpik")[1]
if ensinfo[k].id == "J501"
mhh[k] .-= 0.03 #0.04
msh[k] .-= 0.04 #0.04
fsh[k] .+= 0.01
flh[k] .-= 0.01
end
if ensinfo[k].id == "J500"
mhh[k] .-=0.03 #0.04
flh[k] .+= 0.01
fsh[k] .+=0.01
end
if ensinfo[k].id == "H102r002"
fsh[k] .+= 0.006
end
phih_fa = (2 * mlh[k] + msh[k]) / 3
phih_sfa = (3 * phih_fa + 6 * mlh_v[k] + 3 * msh_v[k]) / 12
phih_eta = mhh[k]
......@@ -136,16 +164,35 @@ for (k, ens) in enumerate(enslist)
zp_aux = fill(ZP(ensinfo[k].beta), length(muh[k]))
muh[k] ./= zp_aux
# p = joinpath(path_bdio_md, ens) # use this to read your t0 values
p = joinpath(path_bdio_md, ens) # use this to read your t0 values
# a28t0[k] = 1 / ( 8 * read_BDIO(p, "md", "t0_shifted")[1]) # use this to read your t0 values
#rep = filter(x->occursin(ensinfo[k].id, x), readdir(path_bdio_md,join=true))[1] # use this to read AS t0 shifted
#a28t0[k] = 1 / ( 8 * read_BDIO(rep, "ll_obs", "t0")[1]) # use this to read AS t0 shifted
# rep = filter(x->occursin(ensinfo[k].id, x), readdir(path_bdio_md,join=true))[1] # use this to read AS t0 shifted
# a28t0[k] = 1 / ( 8 * read_BDIO(rep, "ll_obs", "t0")[1]) # use this to read AS t0 shifted
# t0_ens, _, _ = get_t0(path_data_w, ensinfo[k], rw=true, path_rw=path_rw, pl=true)
# a28t0[k] = 1 / (8*t0_ens)
# this is what you used in thesis and paper
a28t0[k] = 1 / (8 * t0(ensinfo[k].beta))
phi2[k] = mll^2
phi4[k] = mls^2 + 0.5 * mll^2
if ens == "H105r002.bdio"
change_id(from="H105", to="H105_charm")
end
end
##
# ADerrors.wsg.id2str[-12391] = "H105r2"
fb = BDIO_open(path_bdio_t0, "r")
BDIO_seek!(fb)
BDIO_seek!(fb,2)
BDIO_seek!(fb,2)
t0_ph = sqrt(8*read_uwreal(fb))
BDIO_seek!(fb,2)
fpik_phys = read_uwreal(fb)
BDIO_close!(fb)
## X VARIABLE [a^2/(8t0) phi_2 phi_h]
x_fa = [a28t0 phi2 getindex.(getindex.(phih, "phih_fa"), 1)]
......@@ -196,8 +243,11 @@ phi2_tot = Vector{Vector{uwreal}}(undef, length(ensinfo))
# now with AS dm, each obs depends on each ensemble, therefore a different solution
# is needed for computing a28t0_tot and phi2_tot
# Nmh = [2, 3, 3, 3, 3, 3, 3, 3, 3, 3] # number of heavy masses for LAT22 ensemble
Nmh = [2, 2, 3, 3, 3, 3, 3, 2, 3, 3, 3, 3] # number of heavy masses with J500 and E300
Nmh = length.(mhh) # automatic detection
#Nmh = [2, 3, 3, 3, 3, 3, 3, 3, 3, 3] # number of heavy masses for LAT22 ensemble
# Nmh = [2, 2, 3, 3, 3, 3, 3, 2, 3, 3, 3, 3, 2] # number of heavy masses with E300, J500, N302
# Nmh = [2, 2, 3, 3, 3, 3, 3, 2, 3, 3, 3, 3] # number of heavy masses with J500 and E300
for (k, ens) in enumerate(ensinfo)
a28t0_tot[k] = fill(a28t0[k], Nmh[k])
......@@ -223,7 +273,7 @@ end
#============= FIT MESON MASSES =============#
@info("Fitting meson masses for charm quark matching")
# mD
fit_over_cat!(mass_cat[4], f, phi2_ph, n_param, weight_dof=2, correlated_fit=false, Mrat=uwreal(1.0), IC="TIC")
fit_over_cat!(mass_cat[4], f, phi2_ph, n_param, weight_dof=2, correlated_fit=false, Mrat=uwreal(1.0), IC="TIC", wpm=wpmm)
# mDs
fit_over_cat!(mass_cat[5], f, phi2_ph, n_param, weight_dof=2, correlated_fit=false, Mrat=uwreal(1.0), IC="TIC")
# mD_star
......@@ -231,36 +281,55 @@ fit_over_cat!(mass_cat[6], f, phi2_ph, n_param, weight_dof=2, correlated_fit=fal
# mDs_star
fit_over_cat!(mass_cat[7], f, phi2_ph, n_param, weight_dof=2, correlated_fit=false, Mrat=uwreal(1.0), IC="TIC")
##
s = [L"$\sqrt{8t_0} Z^{tm}_M \mu_c$", L"$\sqrt{8t_0}M_{\eta^{(c)}_c}$", L"$\sqrt{8t_0}M_{J/\psi}$",
L"$\sqrt{8t_0}M_D$", L"$\sqrt{8t_0}M_{D_s}$", L"$\sqrt{8t_0}M_{D^{*}}$", L"$\sqrt{8t_0}M_{D^{*}_s}$"]
idx_cat = 5 # 4,5
s = [L"$\sqrt{8t_0} Z^{tm}_M \mu_c$", L"$\sqrt{8t_0}m_{\eta^{(c)}_c}$", L"$\sqrt{8t_0}m_{J/\psi}$",
L"$\sqrt{8t_0}m_D$", L"$\sqrt{8t_0}m_{D_s}$", L"$\sqrt{8t_0}m_{D^{*}}$", L"$\sqrt{8t_0}m_{D^{*}_s}$"]
idx_cat = 4 # 4,5
for i in 2:2
plot_cat(mass_cat[idx_cat][i], s[idx_cat] , Wlabel="TIC", save=false)
plot_chi2_vs_exp(mass_cat[idx_cat][i], label_cutoff , p=nothing)
plot_chi_charm(mass_cat[idx_cat][i], vcat(f_aux...), s[idx_cat], phi2_ph, ensinfo, mrat=uwreal(1.0) , p=path_plot, n=10)
plot_chi_charm(mass_cat[idx_cat][i], vcat(f_aux...), s[idx_cat], phi2_ph, ensinfo, mrat=uwreal(1.0), n=10, p=path_plot)
#plot_cl_charm(mass_cat[idx_cat][i], vcat(f_aux...), s[idx_cat], phi2_ph, ensinfo, mrat=uwreal(1.0))#, p=path_plot)
end
##
#============= FIT CHARM MASS ==============#
@info("Fitting charm quark mass")
fit_over_cat!(mass_cat[1], f, phi2_ph, n_param, weight_dof=2, correlated_fit=false, Mrat=uwreal(1.0), IC="TIC" )
fit_over_cat!(mass_cat[1][1:2], f, phi2_ph, n_param, weight_dof=2, correlated_fit=false, Mrat=uwreal(1.0), IC="TIC", xycovar=false )
## plot charm
for i in 2:2
plot_cat(mass_cat[1][i], L"$M_c^{\mathrm{RGI}}(N_f=3)\mathrm{MeV}$" , Wlabel="TIC", save=false)
plot_chi2_vs_exp(mass_cat[1][i], label_cutoff, p=path_plot)
plot_chi_charm(mass_cat[1][i], vcat(f_aux...), L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo, mrat=uwreal(1.0), p=path_plot)
plot_cl_charm(mass_cat[1][i], vcat(f_aux...), L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo, mrat=uwreal(1.0), p=path_plot)
# fpik_phys = uwreal([0.3091,0.0019], "fpik_phys")
for i in 1:1
# plot_cat(mass_cat[1][i], L"$M_c^{\mathrm{RGI}}(N_f=3)\mathrm{MeV}$" , Wlabel="TIC", save=false)
# plot_chi2_vs_exp(mass_cat[1][i], label_cutoff, p=path_plot)
plot_chi_charm(mass_cat[1][i], vcat(f_aux...), L"$\sqrt{8t_0}\mu_c^R$", phi2_ph, ensinfo, mrat=uwreal(1.0))#, p=path_plot)
plot_mh_interp_charm(mass_cat[1][i], vcat(f_aux...), L"$\sqrt{8t_0}\mu_c^R$", phi2_ph, ensinfo, mrat=uwreal(1.0), Nmh) #, p="/Users/alessandroconigli/MyDrive/phd/analysis/charm_gevp_2022/Lat22")
# plot_cl_charm(mass_cat[1][i], vcat(f_aux...), L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo, mrat=uwreal(1.0), p=path_plot)
end
close("all")
## plot for multiple categories
plot_hist_tot([mass_cat[1][1], mass_cat[1][2]], L"$\sqrt{8t_0}\mu_c^R$" , p=path_plot)
plot_all_cl_charm([mass_cat[1][2], mass_cat[1][1]], vcat(f_aux...),L"$\sqrt{8t_0}\mu_c^R$", phi2_ph, ensinfo, mrat=uwreal(1.0) , n=[0,60], p=path_plot)
plot_cat_tot(mass_cat[1][1:2], L"$\sqrt{8t_0}\mu_c^R$" , Wlabel="TIC", save=false)
plot_all_cl_charm([mass_cat[1][2], mass_cat[1][1]], vcat(f_aux...),L"$\sqrt{8t_0}\mu_c^R$", phi2_ph, ensinfo, mrat=uwreal(1.0) , n=[0,0]) #, p=path_plot)
plot_cat_tot(mass_cat[1][1:1], L"$\sqrt{8t_0}\mu_c^R$" , xlabel=label_cutoff, Wlabel="TIC", save=true)
##
## plot weights vs dof
aic = vcat(mass_cat[1][1].aic, mass_cat[1][2].aic)
ww, _ = aic_weight(aic)
dof = vcat(mass_cat[1][1].dof, mass_cat[1][2].dof)
plot(dof, ww, ms=6, ls="none", marker="d")
xlabel(L"$d.o.f$")
ylabel(L"$W$")
title(L"\mu_c")
display(gcf())
savefig("/Users/alessandroconigli/Desktop/muc_w_over_dof.pdf")
close("all")
##
#============== CREATE DECAY CAT =============#
s_obs_f = ["fD", "fDs", "f_etac"]
s_match = ["fl_ave", "etac", "sp_fl_ave"]
......@@ -271,33 +340,19 @@ obs = cat_obs.([flh, fsh, fhh])
match = [x_fa, x_eta, x_sfa]
## for correlated fits the order must be different: ens11..ens13, ens21...ens23,...,ensN1..ensN3
obs = [vcat(flh...), vcat(fsh...), vcat(fhh...)]
# this commented section is actually the same as what you did for mass categories. therefore it is not needed
# [uwerr.(obs[i]) for i in 1:length(obs)]
# preparing xdata for fitting in correct order
# ens_id = getindex.(ensembles.(obs[1]),2)
# tmp_ens_idx = findall(x->x=="H105", ens_id)
# ens_id[tmp_ens_idx] .*="r002"
obs = [vcat(flh...) .* sqrt.(vcat(mlh...)) , vcat(fsh...).*sqrt.(vcat(msh...)), vcat(fhh...)]
# a28t0_tot = Vector{Vector{uwreal}}(undef, length(ensinfo))
# phi2_tot = Vector{Vector{uwreal}}(undef, length(ensinfo))
# for (k, ens) in enumerate(ensinfo)
# if ens.id == "D200"
# idx = findall(x->x=="D200_corr", ens_id)
# else
# idx = findall(x->x==ens.id, ens_id)
# end
# a28t0_tot[k] = fill(a28t0[k], length(idx))
# phi2_tot[k] = fill(phi2[k], length(idx))
# end
x_fa = [vcat(a28t0_tot...) vcat(phi2_tot...) vcat(getindex.(phih, "phih_fa")...) ]
x_eta = [vcat(a28t0_tot...) vcat(phi2_tot...) vcat(getindex.(phih, "phih_eta")...) ]
x_sfa = [vcat(a28t0_tot...) vcat(phi2_tot...) vcat(getindex.(phih, "phih_sfa")...) ]
match = [x_fa, x_eta, x_sfa]
for (j, m) in enumerate(match)
match[j] = hcat(match[j], fill(fpik_phys^2, sum(Nmh)) )
end
## create categories
#create categories
f_cat = Vector{Vector{Cat}}(undef, 0)
for (i, o) in enumerate(obs)
aux = Vector{Cat}(undef, 0)
......@@ -366,7 +421,7 @@ for i = 1:3
[f1.obs, f2.obs], f_n_param_comb_clog[j], correlated_fit=false, wpm=wpmm)
x = [0.0 phi2_ph f1.phih_ph]
x = [0.0 phi2_ph f1.phih_ph fpik_phys^2]
aux1 = f_f_aux1_clog[1][1](x, par)[1]
aux2 = f_f_aux2_clog[1][1](x, par)[1]
doff = sum(length.([f1.obs, f2.obs])) - f_n_param_comb_clog[j]
......@@ -382,16 +437,43 @@ for i = 1:3
push!(f1.dof, doff)
push!(f2.dof, doff)
# use first 2 lines for uncorrelated fits
push!(f1.aic, chi2/chi2exp * doff + 2 * f_n_param_comb_clog[j])
push!(f2.aic, chi2/chi2exp * doff + 2 * f_n_param_comb_clog[j])
AIC = chi2 - 2*chi2exp # version with TIC
# AIC = chi2/chi2exp * doff + 2 * f_n_param_comb_clog[j] # version with AIC
push!(f1.aic, AIC)
push!(f2.aic, AIC)
# push!(f2.aic, chi2_vs_chi2exp - 2 * doff)
# push!(f1.aic, chi2_vs_chi2exp - 2 * doff)
end
end
end
## FIT DEC. CONSTS PHID + PHIDs updated functional forms (CLOGS)
end
for k in 1:3
cat_fd = f_cat[1][k]
cat_fds = f_cat[2][k]
for i in eachindex(phid_tot_models)
f_d = phid_tot_models[i]
f_ds = phids_tot_models[i]
par, chi2, chi2exp = fit_routine([f_d, f_ds], [value.(cat_fd.x_tofit), value.(cat_fds.x_tofit)], [cat_fd.obs, cat_fds.obs], d_ds_param[i], wpm=wpmm )
xx = [0.0 phi2_ph phih_ph[k] fpik_phys^2]
phid_res = f_d(xx, par)[1]; uwerr(phid_res)
phids_res = f_ds(xx, par)[1]; uwerr(phids_res)
println("phiD = ", phid_res)
println("phiDs = ", phids_res)
dof = sum(length.([cat_fd.obs, cat_fds.obs])) - d_ds_param[i]
symb = [:fit_res, :fit_param, :chi2_vs_exp, :n_param, :aic, :dof]
tmp_fd = [phid_res, par, chi2 / chi2exp, d_ds_param[i], chi2 - 2*chi2exp, dof]
tmp_fds = [phids_res, par, chi2 / chi2exp, d_ds_param[i], chi2 - 2*chi2exp, dof]
for (j,ss) in enumerate(symb)
push!(getfield(cat_fd, ss), tmp_fd[j])
push!(getfield(cat_fds, ss), tmp_fds[j])
end
end
end
......@@ -399,17 +481,24 @@ end
# ffD = vcat(f_f_aux1...)
# ffDs = vcat(f_f_aux2...)
ffD = vcat(f_f_aux1_clog...)
ffDs = vcat(f_f_aux2_clog...)
# old log functions
# ffD = vcat(f_f_aux1_clog...)
# ffDs = vcat(f_f_aux2_clog...)
ffD = phid_tot_models
ffDs = phids_tot_models
# ffD = vcat([vcat(f_f_aux1...), vcat(f_f_aux1_clog...)]...) # lin+log
# ffDs = vcat([vcat(f_f_aux2...), vcat(f_f_aux2_clog...)]...) # lin+log
s = [L"$\sqrt{8t_0}f_D$", L"$\sqrt{8t_0}f_{D_s}$", L"$\sqrt{8t_0}f_{\eta^{(c)}_c}$"]
for i=2:2
plot_cl_charm( f_cat[1][i], ffD, s[1], phi2_ph, ensinfo, dec=true, p=path_plot)
# s = [L"$\sqrt{8t_0}f_D$", L"$\sqrt{8t_0}f_{D_s}$", L"$\sqrt{8t_0}f_{\eta^{(c)}_c}$"]
s = [L"$\Phi_D$", L"$\Phi_{D_s}$", L"$\Phi_{\eta^{(c)}_c}$"]
for i=1:1
plot_cl_charm( f_cat[1][i], ffD, s[1], phi2_ph, ensinfo, dec=true, p=path_plot)
plot_cl_charm( f_cat[2][i], ffDs, s[2], phi2_ph, ensinfo, dec=true, p=path_plot)
plot_chi_charm(f_cat[1][i], ffD, s[1], phi2_ph, ensinfo, dec=true, p=path_plot)
plot_chi_charm(f_cat[2][i], ffDs, s[2], phi2_ph, ensinfo, dec=true, p=path_plot)
plot_chi_charm(f_cat[2][i], ffDs, s[2], phi2_ph, ensinfo, dec=true, p=path_plot)
plot_cat(f_cat[1][i], s[1] , save=false)
plot_cat(f_cat[2][i], s[2] , save=false)
# plot_chi2_vs_exp(f_cat[1][i], label_cutoff)#, p=path_plot)
......@@ -418,11 +507,45 @@ for i=2:2
# plot_chi2_vs_exp_clog(f_cat[2][i], label_cutoff , p=path_plot)
end
close("all")
## plot weights vs dof
aicc = vcat(f_cat[1][1].aic, f_cat[1][2].aic)
ww, _ = aic_weight(aicc)
dof = vcat(f_cat[1][1].dof, f_cat[1][2].dof)
plot(dof, ww, ms=6, ls="none", marker="d")
xlabel(L"$d.o.f$")
ylabel(L"$W$")
title(L"f_D")
display(gcf())
savefig("/Users/alessandroconigli/Desktop/fD_w_over_dof.pdf")
close("all")
##
plot_cl_charm(f_cat[2][1], ffD,L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo, dec=true, n=1 ) # 35 nice
plot_chi_charm(f_cat[2][2], ffD, L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo, dec=true, n=1) #, p=path_plot)
########################## SECOND TESTS RATIO - HQET-TAYLOR 2023
plot_cl_charm(f_cat[1][1], ffD, L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo, dec=true) # 35 nice
plot_chi_charm(f_cat[1][2], ffD, L"$\sqrt{8t_0}Z_M^{tm}\mu_c$", phi2_ph, ensinfo, dec=true) #, p=path_plot)
fd_test, fd_syst_test = category_av_tot(f_cat[1]) ; uwerr(fd_test);
fds_test, fds_syst_test = category_av_tot(f_cat[2]) ; uwerr(fds_test);
## extra plot with symm point ensembles
plot_cl_charm_no_proj(f_cat[1][2], ffD, s[1], phi2_ph, ensinfo, dec=true, p=path_plot)
##
idx_ens_sim = findall(x-> x==true, getfield.(ensinfo, :deg))
yy = vcat(flh[idx_ens_sim]...); uwerr.(yy)
xx = vcat(a28t0_tot[idx_ens_sim]...); uwerr.(xx)
##
errorbar(value.(xx), value.(yy), xerr=err.(xx), yerr=err.(yy), fmt="s")
display(gcf())
close("all")
##
########################## RATIO fDs/fD - HQET-TAYLOR functional forms 2023
## Second tests on ratio FDs/FD with Gregorio's formulation
## creating categories
......@@ -439,11 +562,11 @@ x_sfa = [vcat(a28t0_tot...) vcat(phi2_tot...) vcat(getindex.(phih, "phih_sfa")..
# this extra lines are requires to add extra x dependence (fpik) in fds ratio
match = [x_fa, x_eta, x_sfa]
# Nmh = [2, 3, 3, 3, 3, 3, 3, 3, 3, 3] # number of heavy masses for LAT22 ensemble
Nmh = [2, 2, 3, 3, 3, 3, 3, 2, 3, 3, 3, 3] # number of heavy masses with J500 and E300
# Nmh = [2, 2, 3, 3, 3, 3, 3, 2, 3, 3, 3, 3, 2] # number of heavy masses with J500 and E300 and N302
fpik_tot = Vector{Vector{uwreal}}(undef, length(ensinfo))
fpik_phys = uwreal([0.3091,0.0019], "fpik_phys")
# fpik_phys = uwreal([0.3091,0.0019], "fpik_phys")
for (j, m) in enumerate(match)
xx = m[idx_no_sym,1:3]
xx = hcat(xx, fill(phih_ph[j], length(xx[:,1])))
......@@ -485,16 +608,15 @@ for (k,cat) in enumerate(f_ratio_cat)
end
end
plot_cat_tot(f_ratio_cat, L"$\Phi_{D_s}/\Phi_D$", xlabel=r_tot_labels, save=true)
best, syst = category_av_tot(f_ratio_cat) .* sqrt(M[1]/M[3]) ; uwerr(best);
println( best, " ", value(syst))
##
ratio_fds, ratio_fds_syst = category_av_tot(f_ratio_cat) .* sqrt(M[1]/M[3]) ; uwerr(ratio_fds);
println( ratio_fds, " ", value(ratio_fds_syst))
## cl and chiral plots
idx_ens_nosim = findall( x-> x==false, getfield.(ensinfo, :deg))
plot_chi_fds_ratio(f_ratio_cat[2], r_tot_models, L"$\Phi_{D_s}/\Phi_D$", phi2_ph, ensinfo[idx_ens_nosim], p=path_plot)
plot_cl_fds_ratio(f_ratio_cat[2], r_tot_models, L"$\Phi_{D_s}/\Phi_D$", phi2_ph, ensinfo[idx_ens_nosim], p=path_plot)
##
xx = [0.0 phi2_ph phih_ph[1] phih_ph[1]]
res = r_hqet_models[1](xx,par)[1] ; uwerr(res)
bb = r_model_aux(value.(xx), par)
......@@ -534,7 +656,7 @@ push!(f_ratio_cat.fit_param, par)
idx_ens_nosim = findall( x-> x==false, getfield.(ensinfo, :deg))
#obs_ratio .+= 0.01
plot_chi_fds_ratio(f_ratio_cat, ratio_model_aux, L"$f_{D_s}/f_D$", phi2_ph, ensinfo[idx_ens_nosim], p=path_plot)
plot_chi_fds_ratio(f_ratio_cat, ratio_model_aux, L"$f_{D_s}/f_D$", phi2_ph, ensinfo[idx_ens_nosim], p=nothing)
##
#############
......@@ -543,8 +665,9 @@ plot_chi_fds_ratio(f_ratio_cat, ratio_model_aux, L"$f_{D_s}/f_D$", phi2_ph, ensi
#t0 [2.86, 3.659, 5.164, 8.595]
wpmm2 = Dict{String, Vector{Float64}}()
wpmm2["H101"] = [-1.0, 1.5, -1.0, -1.0]
wpmm2["H105"] = [-1.0, 1.5, -1.0, -1.0]
wpmm2["H102r002"] = [-1.0, 2.0, -1.0, -1.0]
wpmm2["H400"] = [-1.0, 4.0, -1.0, -1.0]
wpmm2["H400"] = [-1.0, 2.0, -1.0, -1.0]
wpmm2["N202"] = [-1.0, 1.0, -1.0, 14.0*5.146]
wpmm2["N200"] = [-1.0, 2.0, -1.0, -1.0]
wpmm2["N203"] = [-1.0, 2.0, -1.0, -1.0]
......@@ -557,7 +680,7 @@ wpmm2["J303"] = [-1.0, 2.0, -1.0, -1.0]
for (k, c) in enumerate([mass_cat[1]])
o, sys = category_av_tot([c[1],c[2]])
o_p = o * hc / t0_ph * Mrat
sys_p = sys * hc / value(t0_ph) * Mrat
sys_p = sys * hc / t0_ph * Mrat
uwerr(o, wpmm2)
uwerr(o_p, wpmm2)
println(s_obs_mass[k], "\t $o +/- $sys")
......@@ -566,7 +689,7 @@ for (k, c) in enumerate([mass_cat[1]])
end
##
# fl-av / sp-fl-av / etac
for i = 1:3
for i = 1:2
for (k, c) in enumerate([mass_cat[1]])
o = aic_model_ave(c[i])
sys = aic_syst(c[i])
......@@ -579,12 +702,17 @@ end
##
#============ DECAY CONSTANTS ==============#
store = []
for (k, c) in enumerate(f_cat[1:2])
o, sys = category_av_tot([c[1], c[2]])
o_p = o * hc / t0_ph #value(t0_ph)
sys_p = sys * hc / t0_ph #value(t0_ph)
o, sys = category_av_tot([c[1], c[2]]) #./ sqrt(M[1])
# o_p = o * hc / t0_ph #value(t0_ph)
# sys_p = sys * hc / t0_ph #value(t0_ph)
mmm = k==1 ? sqrt(M[1]) : sqrt(M[3])
o_p = o * hc^1.5 / t0_ph^1.5 / mmm
sys_p = sys * hc^1.5 / value(t0_ph^1.5) / mmm
uwerr(o, wpmm2)
uwerr(o_p, wpmm2)
push!(store, o_p)
#uwerr(o)
#uwerr(o_p)
println(s_obs_f[k], "\t$o +/- $sys")
......@@ -595,7 +723,7 @@ for (k, c) in enumerate(f_cat[1:2])
end
# fl-av / sp-fl-av / etac
##
for i = 1:3
for i = 1:2
for (k, c) in enumerate(f_cat[1:2])
o = aic_model_ave(c[i])
sys = aic_syst(c[i])
......@@ -604,9 +732,21 @@ for i = 1:3
end
println()
end
##
#============ RATIO fDs =================#
o_ratio, sys_ratio = category_av_tot(f_ratio_cat[1:2]) .* sqrt(M[1]/M[3])
uwerr(o_ratio, wpmm2)
println("fDs/fD ", o_ratio, " +/- ", sys_ratio )
for i in 1:2
o = aic_model_ave(f_ratio_cat[i]) * sqrt(M[1]/M[3])
sys = aic_syst(f_ratio_cat[i]) * sqrt(M[1]/M[3])
uwerr(o,wpmm2)
println("fDs/fD", "\t $o +/- $sys")
end
##
......@@ -722,4 +862,36 @@ for (i, covar) in enumerate(test)
end
##
## PARAMETER P[4] id FD and FDs fits
0.036 +/- 0.11
0.038 +/- 0.117
0.049 +/- 0.117
0.055 +/- 0.11
0.049 +/- 0.11
0.410 +/- 0.25
0.051 +/- 0.11
0.051 +/- 0.11
0.050 +/- 0.11
0.051 +/- 0.11
0.052 +/- 0.11
0.053 +/- 0.11
0.442 +/- 0.25
0.433 +/- 0.25
0.414 +/- 0.25
0.372 +/- 0.25
0.051 +/- 0.11
0.053 +/- 0.11
0.052 +/- 0.11
0.052 +/- 0.11
0.401 +/- 0.25
0.401 +/- 0.25
0.401 +/- 0.25
0.404 +/- 0.25
0.405 +/- 0.25
0.406 +/- 0.25
0.052 +/- 0.117
0.401 +/- 0.25
0.407 +/- 0.251
0.406 +/- 0.25
0.405 +/- 0.25
0.405 +/- 0.25
\ No newline at end of file
src/gevp analysis/gevp_plotter.jl
View file @ e56a3bed
......@@ -47,32 +47,35 @@ function plot_cat_tot(cat::Vector{Cat}, ylab::LaTeXString ; xlabel::Union{Nothin
# uwerr(syst)
ratio_w = weigths #aic ./ minimum(aic)
fig = figure(figsize=(10,5))
# fig = figure(figsize=(10,5))
fig = figure(figsize=(16,8))
# cm = get_cmap("YlGn")
cm = get_cmap("Purples")
fill_between(collect(1:length(all_res)), value(best-syst), value(best+syst), alpha=0.3, color="orange", label="Systematics")
sc = scatter(collect(1:length(all_res)), value.(all_res), edgecolors="black",c= ratio_w, cmap=cm, s=70, vmin=minimum(ratio_w), vmax=maximum(ratio_w) , zorder=3)
fill_between(collect(1:length(all_res[1:32])), value(best-syst), value(best+syst), alpha=0.3, color="orange", label="Systematics")
sc = scatter(collect(1:length(all_res[1:32])), value.(all_res[1:32]), edgecolors="black",c= ratio_w[1:32], cmap=cm, s=70, vmin=minimum(ratio_w[1:32]), vmax=maximum(ratio_w[1:32]) , 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="s", lw=0.8, mec="black", ecolor="black")
errorbar(collect(1:length(all_res[1:32])), value.(all_res[1:32]), yerr=err.(all_res[1:32]), fmt="none",zorder=0, capsize=2,marker="s", lw=0.8, mec="black", ecolor="black")
errorbar(-2, value(best), yerr=err(best), fmt="s",mfc="maroon",ms=10, mec="maroon",ecolor="maroon", capsize=2, zorder=0, label="Model average")
# xlabel("Models", size=18)
ylabel(ylab, size=20)
ylabel(ylab, size=22)
# axvline(64, ls="--", color="red")
if !(isnothing(xlabel))
lbl = vcat(fill(xlabel, length(cat))...)
xx = collect(1:length(lbl))
println(idx)
xticks(xx, lbl, rotation=75, size=9)
xticks(xx, lbl, rotation=90, size=13)
end
cb = colorbar(sc)
cb[:set_label]( L"$w$")
# cb[:set_label]( string(Wlabel," weights"))
legend()
legend(fontsize=16)
#ylim(0.46,0.58)
tight_layout()
display(fig)
if save
t = string("cat_ave_tot",cat[1].info, ".pdf")
t = string("cat_ave_tot_fDs",cat[1].info, ".pdf")
savefig(joinpath(path_plot,t))
end
end
......@@ -108,7 +111,7 @@ function plot_cl_charm(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph:
println("chi2/chiexp corresponding to the best fit ", c.chi2_vs_exp[n] )
##PLOTS
#xx, yy -> Continuum limit band
xx = [Float64.(range(1e-8, a2_max, length=100)) fill(phi2_ph, 100) fill(phih_2_use,100)]
xx = [Float64.(range(1e-8, a2_max, length=100)) fill(phi2_ph, 100) fill(phih_2_use,100) fill(fpik_phys^2,100)]
#xx = [Float64.(range(0.0, a2_max, length=100)) fill(phi2[1], 100) fill(phih_2_use,100)]
yy = func(xx, upar) * value(mrat)
#Project one value of muc_i (i=1,2,3)
......@@ -116,10 +119,13 @@ function plot_cl_charm(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph:
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 fill(phi2_ph, length(a2)) fill(phih_2_use, length(p2))] #physical value phih
#x2 = [a2 fill(phi2[1], length(a2)) fill(phih_2_use, length(p2))] #project to symmetric line only
pfpik = c.x_tofit[n_unique, 4]
x = [a2 p2 pc pfpik]
x2 = [a2 fill(phi2_ph, length(a2)) fill(phih_2_use, length(p2)) fill(fpik_phys^2, length(p2))] #physical value phih
# x2 = [a2 p2 fill(phih_2_use, length(p2)) fill(fpik_phys^2, length(p2))] #physical value phih
#x2 = [a2 fill(phi2[1], length(a2)) fill(phih_2_use, length(p2))] #project to symmetric line only
y = ( c.obs[n_unique] - func(x, upar) + func(x2, upar)) * value(mrat) #Projection
#Uwerr
if isnothing(wpm)
......@@ -156,6 +162,95 @@ function plot_cl_charm(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph:
close("all")
end
end
# plot continuum limit for a given category without projectin to physical point
function plot_cl_charm_no_proj(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph::uwreal, ensinfo::Vector{EnsInfo}; mrat::uwreal=uwreal(1.), dec::Bool=false, 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", "royalblue"]
fmt = ["^", "v", "<", ">", "d"]
a_sp = [0.087, 0.077, 0.065, 0.05, 0.039]
#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]
if !dec
res = value(mrat) * (upar[1] + upar[2]*phi2_ph + upar[3]*phih_2_use)
else
res = getfield(c, :fit_res)[n]
end
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_vs_exp[n] )
##PLOTS
#xx, yy -> Continuum limit band
xx = [Float64.(range(1e-8, a2_max, length=100)) fill(phi2_ph, 100) fill(phih_2_use,100) fill(fpik_phys^2,100)]
#xx = [Float64.(range(0.0, a2_max, length=100)) fill(phi2[1], 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]
pfpik = c.x_tofit[n_unique, 4]
x = [a2 p2 pc pfpik]
x2 = [a2 p2 pc fill(fpik_phys^2, length(p2))] #physical value phih
# x2 = [a2 p2 fill(phih_2_use, length(p2)) fill(fpik_phys^2, length(p2))] #physical value phih
#x2 = [a2 fill(phi2[1], length(a2)) fill(phih_2_use, length(p2))] #project to symmetric line only
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
idx_sim_ens = findall(x->x==true, getfield.(ensinfo, :deg))
for (k,b) in enumerate(sort(getfield.(ensinfo[idx_sim_ens], :beta)))
println(ensinfo[idx_sim_ens[k]].id, " ", b)
n_ = findall(x->x.beta == b, ensinfo[idx_sim_ens])
a2_aux = mean(value.(a2[n_]))
if ensinfo[idx_sim_ens[k]].deg == true
errorbar(value.(x[idx_sim_ens[n_],1]), value.(y[idx_sim_ens[n_]]), err.(y[idx_sim_ens[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
end
ylabel(ylab, size=20)
xlabel(L"$a^2/8t_0$", size=20)
legend(ncol=2, fontsize=16)
grid(ls="dashed")
tight_layout()
#ylim(0.46, 0.58)
if !isnothing(p)
t = "fit_cl_symmonly" * c.info * ".pdf"
savefig(joinpath(p, t))
close()
else
display(gcf())
close("all")
end
end
# plot continuum limit for all categories
function plot_all_cl_charm(cvec::Vector{Cat}, ff::Vector{Function}, ylab::LaTeXString, phi2_ph::uwreal, ensinfo::Vector{EnsInfo}; mrat::uwreal=uwreal(1.), dec::Bool=false, p::Union{String, Nothing}=nothing,wpm::Union{Dict{String,Vector{Float64}}, Nothing}=nothing, n::Union{Nothing, Vector{Int64}}=nothing)
......@@ -256,7 +351,7 @@ function plot_chi_charm(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph
end
#Model info
@info("error in t0_ph not included")
phih_2_use = c.phih_ph
phih_2_use = value(c.phih_ph)
par_n = getfield(c, :n_param)[n]
upar = getfield(c, :fit_param)[n]
if !dec
......@@ -270,7 +365,7 @@ function plot_chi_charm(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph
println("chi2/chiexp corresponding to the best fit ", c.chi2_vs_exp[n] )
##PLOTS
#xx, yy -> Continuum limit band
xx = [fill(1e-8, 100) Float64.(range(0.01, phi2_max, length=100)) fill(phih_2_use,100)]
xx = [fill(1e-8, 100) Float64.(range(0.01, phi2_max, length=100)) fill(phih_2_use,100) fill(fpik_phys^2, 100)]
yy = func(xx, upar) * value(mrat)
#Project one value of muc_i (i=1,2,3)
......@@ -278,8 +373,14 @@ function plot_chi_charm(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph
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
pfpik=0.
try
pfpik = c.x_tofit[n_unique, 4]
catch
pfpik = c.x_tofit[n_unique, 3]
end
x = [a2 p2 pc pfpik]
x2 = [a2 p2 fill(phih_2_use, length(p2)) fill(fpik_phys^2,length(p2))] #physical value phih
y = ( c.obs[n_unique] - func(x, upar) + func(x2, upar)) * value(mrat) #Projection
#Uwerr
if isnothing(wpm)
......@@ -301,18 +402,20 @@ function plot_chi_charm(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph
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)]
xxx = [fill(a2_aux, 100) Float64.(range(0.0, phi2_max, length=100)) fill(value(c.phih_ph), 100) fill(value(fpik_phys^2),100)]
yyy = func(xxx, upar) * value(mrat)
plot(xxx[:,2], value.(yyy), ls="--", color=color[k])#dashed lines
end
# fill between physical values
am = M[3] * t0_ph / hc
am = M[1] * t0_ph / hc
uwerr(am)
#fill_between(xx[:,2], value(am)+err(am), value(am)-err(am), color="forestgreen", alpha=0.2, zorder=1, label="PDG")
# fill_between(xx[:,2], value(am)+err(am), value(am)-err(am), color="forestgreen", alpha=0.2, zorder=1, label="PDG")
ylabel(ylab, size=22)
xlabel(L"$\phi_2$", size=22)
legend(ncol=2, fontsize=16)
#ylim(3.85, 4.3)
# ylim(0.825, 1.2)
# ylim(3.80, 4.2)
ylim(3.235,3.38)
tight_layout()
grid(ls="dashed")
if !isnothing(p)
......@@ -324,6 +427,138 @@ function plot_chi_charm(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph
close("all")
end
end
# plot interpolation in heavy_pc projected for a given category
function plot_mh_interp_charm(c::Cat, ff::Vector{Function}, ylab::LaTeXString, phi2_ph::uwreal, ensinfo::Vector{EnsInfo}, nmh::Vector{Int64}; p::Union{String, Nothing}=nothing,
mrat::uwreal=uwreal(1.), dec::Bool=false,wpm::Union{Dict{String,Vector{Float64}}, Nothing}=nothing, n::Union{Nothing, Int64}=nothing)
#Format
phi2_max = maximum(value.(phi2))
phih_max = maximum(value.(c.x_tofit[:,3]))
phih_min = minimum(value.(c.x_tofit[:,3]))
color = ["orange", "darkgreen", "magenta", "navy", "royalblue"]
fmt = ["^", "v", "<", ">", "d"]
a_sp = [0.087, 0.077, 0.065, 0.05, 0.039]
#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
# @info("error in phi_h from t0_ph not included")
phih_2_use = c.phih_ph
par_n = getfield(c, :n_param)[n]
upar = getfield(c, :fit_param)[n]
if !dec
res = value(mrat) * (upar[1] + upar[2]*phi2_ph + upar[3]*phih_2_use)
else
res = getfield(c, :fit_res)[n]
end
func = ff[n]
println("Results corresponding to the best fit: ", res)
println("TIC corresponding to the best fit ", aic)
println("chi2/chiexp corresponding to the best fit ", c.chi2_vs_exp[n] )
##PLOTS
#xx, yy -> Continuum limit band
xx =[fill(1e-8,100) fill(phi2_ph, 100) Float64.(range(phih_min, phih_max, length=100)) fill(fpik_phys^2, 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,3], 1:length(c.x_tofit[:,3])) #unique(i-> c.x_tofit[i,3], 1:length(c.x_tofit[:,2]))
a2 = c.x_tofit[:, 1]
p2 = c.x_tofit[:, 2]
pc = c.x_tofit[:, 3]
pfpik=0.
try
pfpik = c.x_tofit[:, 4]
catch
pfpik = c.x_tofit[:, 3]
end
x = [a2 p2 pc pfpik]
println(length(x[:,3]))
x2 = [fill(0.0, length(p2)) fill(phi2_ph, length(p2)) pc fill(fpik_phys^2,length(p2))] #physical value phih
y = ( c.obs[:] - func(x, upar) + func(x2, upar)) * value(mrat) #Projection
println(length(y))
#Uwerr
if isnothing(wpm)
uwerr.(yy)
uwerr(res)
uwerr.(y)
uwerr.(x[:,3])
else
[uwerr(yy_aux, wpm) for yy_aux in yy]
uwerr(res, wpm)
[uwerr(y_aux, wpm) for y_aux in y]
[uwerr(x_aux, wpm) for x_aux in x[:,3]]
end
#Plotting
figure(figsize=(10,6))
fill_between(xx[:,3], value.(yy)+ err.(yy), value.(yy)-err.(yy), alpha=0.3, color="royalblue") #CL extrapolation (band)
errorbar(value(phih_2_use), value(res), err(res), fmt="s", color="red", label=L"$\phi_H = \phi^\mathrm{phys}_H$", capsize=5,zorder=3, ms=12)#CL extrapolation
b_already_legended = []
cc_idx = 1
for (k, ens) in enumerate(ensinfo)
β = ens.beta
dc_beta = Dict(
3.4 => ["orange", a_sp[1]],
3.46 => ["darkgreen", a_sp[2]],
3.55 => ["magenta", a_sp[3]],
3.70 => ["navy", a_sp[4]],
3.85 => ["royalblue", a_sp[5]]
)
idx = k==1 ? (1:nmh[1]) : ((cc_idx):(cc_idx -1 + nmh[k]))
println( "k= ", k, " ", idx)
cc_idx += nmh[k]
if β in b_already_legended
errorbar(value.(x[idx,3]), value.(y[idx]), xerr=err.(x[idx,3]), yerr=err.(y[idx]), fmt="d", color=dc_beta[β][1], capsize=5, ms=12, mfc="none" )#, label=string(L"$a \approx$", dc_beta[β][2], " fm")) #projected datapoints
else
errorbar(value.(x[idx,3]), value.(y[idx]), xerr=err.(x[idx,3]), yerr=err.(y[idx]), fmt="d", color=dc_beta[β][1], capsize=5, ms=12, mfc="none" , label=string(L"$a \approx$", dc_beta[β][2], " fm")) #projected datapoints
end
push!(b_already_legended, β)
end
# errorbar(value.(x[:,3]), value.(y[:]), err.(y[:]), fmt="d", color="black", capsize=5, ms=12, mfc="none")# label=string(L"$a \approx$", a_sp[k], " fm") projected datapoints
#
# for (k,b) in enumerate(sort((getfield.(ensinfo, :beta))))
# println("k= ", k, "b= ", b)
# n_ = findall(x-> b === x.beta , ensinfo)
# println(n_)
# println(value.(y[n_]))
# println(value.(x[n_,3]))
# a2_aux = mean(value.(a2[n_]))
# errorbar(value.(x[n_,3]), value.(y[n_]), err.(y[n_]), fmt="d", color="black", capsize=5, ms=12, mfc="none")# label=string(L"$a \approx$", a_sp[k], " fm") projected datapoints
#dashed lines
# xxx = [fill(a2_aux, 100) fill(phi2_ph, 100) Float64.(range(phih_min, phih_max, length=100)) fill(value(fpik_phys^2),100)]
# yyy = func(xxx, upar) * value(mrat)
# plot(xxx[:,3], value.(yyy), ls="--", color="black") #dashed lines
# end
#fill between physical values
am = M[3] * t0_ph / hc
uwerr(am)
# fill_between(xx[:,2], value(am)+err(am), value(am)-err(am), color="forestgreen", alpha=0.2, zorder=1, label="PDG")
ylabel(ylab, size=22)
xlabel(L"$\phi_H$", size=22)
legend(ncol=2, fontsize=16)
# ylim(0.825, 1.155)
# ylim(3.85, 4.3)
# ylim(3.235,3.38)
tight_layout()
grid(ls="dashed")
if !isnothing(p)
t = "fit_phih_interp_" * c.info * ".pdf"
savefig(joinpath(p, t))
close()
else
display(gcf())
close("all")
end
end
# plot chi2_vs_exp adn weights for a given category
function plot_chi2_vs_exp(c::Cat, label_in::Vector{Vector{String}}; p::Union{String, Nothing}=nothing)
......@@ -567,7 +802,7 @@ function plot_chi_fds_ratio(c::Cat, ff_vec::Vector{Function}, ylab::LaTeXString,
color = ["orange", "darkgreen", "magenta", "navy"]
fmt = ["^", "v", "<", ">"]
a_sp = [0.087, 0.065, 0.05]
a_sp = [0.087, 0.065, 0.05, 0.039]
phi2_max = phi2_sym #maximum(value.(phi2))
phih_2_use = c.phih_ph
......@@ -577,7 +812,6 @@ function plot_chi_fds_ratio(c::Cat, ff_vec::Vector{Function}, ylab::LaTeXString,
end
ff = ff_vec[n]
upar = getfield(c, :fit_param)[n]
phih_2_use = value(c.phih_ph)
res = getfield(c, :fit_res)[n]
#fpik_phys = ((2/3) * (155.7 + 0.5*130.2) * t0_ph / hc ) ^2
xx = [fill(0.0, 100) Float64.(range(0.01, phi2_max, length=100)) fill(phih_2_use,100) fill(phih_2_use,100) fill(c.x_tofit[1,5],100) ]
......@@ -607,16 +841,16 @@ function plot_chi_fds_ratio(c::Cat, ff_vec::Vector{Function}, ylab::LaTeXString,
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 = ff(xxx, upar)
#plot(xxx[:,2], value.(yyy), ls="--", color=color[k])#dashed lines
xxx = [fill(a2_aux, 100) Float64.(range(0.0, phi2_max, length=100)) fill(value(c.phih_ph), 100) fill(phih_2_use, 100) fill(c.x_tofit[1,5], 100) ]
yyy = ff(xxx, upar)
plot(xxx[:,2], value.(yyy), ls="--", color=color[k])#dashed lines
end
ylabel(ylab, size=22)
xlabel(L"$\phi_2$", size=22)
legend(ncol=2, fontsize=16)
#ylim(3.02, 3.18)
tight_layout()
grid(ls="dashed")
# grid(ls="dashed")
if !isnothing(p)
t = string("fit_fds_ratio_", c.info, ".pdf")
println( t)
......@@ -630,6 +864,67 @@ function plot_chi_fds_ratio(c::Cat, ff_vec::Vector{Function}, ylab::LaTeXString,
end
function plot_cl_fds_ratio(c::Cat, ff_vec::Vector{Function}, ylab::LaTeXString, phi2_ph::uwreal, ensinfo::Vector{EnsInfo}; p::Union{String, Nothing}=nothing, n::Union{Nothing, Int64}=nothing )
color = ["orange", "darkgreen", "magenta", "navy"]
fmt = ["^", "v", "<", ">"]
a_sp = [0.087, 0.065, 0.05, 0.039]
phi2_max = phi2_sym #maximum(value.(phi2))
phih_2_use = c.phih_ph
if isnothing(n)
aic, n = findmin(getfield(c, :aic))
end
ff = ff_vec[n]
upar = getfield(c, :fit_param)[n]
phih_2_use = value(c.phih_ph)
res = getfield(c, :fit_res)[n]
xx = [Float64.(range(1e-8, 0.05, length=100)) fill(phi2_ph, 100) fill(phih_2_use,100) fill(phih_2_use,100) fill(c.x_tofit[1,5],100) ]
yy = ff(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]
pfpik = c.x_tofit[n_unique, 5]
x = [a2 p2 pc fill(phih_2_use, length(p2)) pfpik ]
x2 = [a2 fill(phi2_ph, length(p2)) fill(phih_2_use, length(p2)) fill(phih_2_use, length(p2)) fill(c.x_tofit[1,5], length(p2))] #physical value phih
y = ( c.obs[n_unique] - ff(x, upar) + ff(x2, upar)) #Projection
uwerr.(yy)
uwerr(res)
uwerr.(y)
println(res)
#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=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_,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
#dashed lines
#xxx = [fill(a2_aux, 100) Float64.(range(0.0, phi2_max, length=100)) fill(value(c.phih_ph), 100) fill(phih_2_use, 100) fill(c.x_tofit[1,5], 100) ]
#yyy = ff(xxx, upar)
#plot(xxx[:,2], value.(yyy), ls="--", color=color[k])#dashed lines
end
ylabel(ylab, size=22)
xlabel(L"$a^2/8t_0$", size=22)
legend(ncol=2, fontsize=16)
ylim(1.17, 1.23)
tight_layout()
# grid(ls="dashed")
if !isnothing(p)
t = string("fit_cl_fds_ratio_", c.info, ".pdf")
println( t)
savefig(joinpath(p, t))
close()
else
display(gcf())
close("all")
end
end
####
## PLOT FROM JAVIER TO TEST decays
###
......
src/gevp analysis/gevp_tm_BMA.jl
View file @ e56a3bed
......@@ -50,14 +50,14 @@ const path_bdio_tm_shifted = "/Users/alessandroconigli/Lattice/data/bdio_charm/t
const path_plot = "/Users/alessandroconigli/MyDrive/phd/analysis/charm_gevp_2022/CharmPaper/plots" # path were to save plots
# ensembles to analyse
const ens_list = [ "E300"]
const ens_list = [ "N203"]
# reweighting and mass shift flags
const rwf = true
println("Warning! mass_shift set to false ")
# println("Warning! mass_shift set to false ")
const mass_shift = false
const tau = 3 # gevp shift parameter
const T0 = 2
const TSM = true # set whether TSM is used or not
const TSM = false # set whether TSM is used or not
#@warning("\nTSM FLAG MODE: $(TSM)\n")
......@@ -89,9 +89,12 @@ wpmm["N200"] = [-1.0, 2.0, -1.0, -1.0] #14.0*5.164]
wpmm["N203"] = [-1.0, 2.0, -1.0, -1.0] # 14.0*5.164]
#wpmm["D200"] = [-1.0, 4.0, -1.0, 14.0*5.164]
wpmm["D200"] = [-1.0, 2.0, -1.0, -1.0]
wpmm["N300"] = [-1.0, 2.0, -1.0, -1.0 ] #14.0*8.595]
wpmm["N300r002"] = [-1.0, 2.0, -1.0, -1.0 ] #14.0*8.595]
wpmm["N302"] = [-1.0, 2.0, -1.0, -1.0 ] #14.0*8.595]
wpmm["J303"] = [-1.0, 2.0, -1.0, -1.0 ] #14.0*8.595]
wpmm["J500"] = [-1.0, 2.0, -1.0, -1.0 ] #14.0*8.595]
wpmm["J501"] = [-1.0, 2.0, -1.0, -1.0 ] #14.0*8.595]
wpmm["E250"] = [-1.0, 2.0, -1.0, -1.0 ] #14.0*8.595]
## TESTING TSM WITH MULTIPLE replicas
data_sl = read_data_sloppy(path_data_tm, "E300", "G5", "G5", legacy=false)
......@@ -121,12 +124,13 @@ obs_tm = Vector(undef, length(ensinfo))
@time begin
@showprogress for (k,ens) in enumerate(ensinfo)
println("\n Ensemble: ", ens.id)
# pp corr
corr_pp[k], pp_tmW[k], _ = get_corr(path_data_tm, ens, "G5", "G5", rw=rwf, path_rw=path_rw, tsm=TSM)
# t0
t0_ens, yw_store[k], w_store[k] = get_t0(path_data_w, ens, rw=rwf, path_rw=path_rw, pl=true)
obs_tm[k] = OrderedDict()
obs_tm[k]["t0"] = t0_ens
# correlators
corr_pp[k], pp_tmW[k], _ = get_corr(path_data_tm, ens, "G5", "G5", rw=rwf, path_rw=path_rw, tsm=TSM)
a1a1, _ = get_corr(path_data_tm, ens, "G1G5", "G1G5", rw=rwf, path_rw=path_rw, tsm=TSM)
a2a2, _ = get_corr(path_data_tm, ens, "G2G5", "G2G5", rw=rwf, path_rw=path_rw, tsm=TSM)
......@@ -144,6 +148,36 @@ obs_tm = Vector(undef, length(ensinfo))
end
end
## plots for paper of D, D_s and eta_c correlators
yyd = corr_pp[1][4].obs; uwerr.(yyd)
yyds = corr_pp[1][5].obs; uwerr.(yyds)
yyeta = corr_pp[1][6].obs; uwerr.(yyeta)
xx = collect(1:length(yyd))
errorbar(xx, value.(yyd), err.(yyd).*50, fmt="s", color="forestgreen", capsize=2, mfc="none", label=L"$D$")
errorbar(xx, value.(yyds), err.(yyds).*50, fmt="s", color="tomato", capsize=2, mfc="none", label=L"$D_s$")
errorbar(xx, value.(yyeta), err.(yyeta).*50, fmt="s", color="royalblue", capsize=2, mfc="none", label=L"$\eta_c$")
yscale("log")
legend()
xlim(65,128)
xlabel(L"$t/a$")
ylabel(L"$f_{\mathrm{PP}(x_0, y_0)}$")
# ylim(1e-22, 1e-4)
tight_layout()
display(gcf())
savefig("/Users/alessandroconigli/Desktop/corr_comparison.pdf")
close("all")
##
# test masses D200
plat = [[85,105], [85,105], [85,105], [85,90], [88,98], [90,94], [85,90], [88,94], [90,94] ]
# test masses J501
plat = [[145,165], [120,155], [120,155], [125,155], [145,155], [145,165], [125,150], [145,165], [145,165] ]
# test masses E250
plat = [[115,145], [115,160], [115,160], [118,125], [125,142], [130,160], [118,125], [125,142], [130,160] ]
# test masses J500
plat = [[130,160], [140,160], [138,150], [140,160], [138,150]]
mtest = meff(corr_pp[1][5], [138,150], pl=true, wpm=wpmm)
##
#=============== ANALYSIS ===============#
Mtotps = Vector{Vector{uwreal}}(undef, length(ensinfo)) # required for decay constants
......@@ -152,9 +186,12 @@ Mtotvec = Vector{Vector{uwreal}}(undef, length(ensinfo)) # required for decay co
@time begin
@showprogress for (k,ens) in enumerate(ensinfo)
# gevp matrix for ps and vec masses
mat_mass = comp_mat_multigamma(ens, corr_pp[k], corr_pa[k], corr_aa[k])
m_tm_ps, m_tm_vec = gevp_mass_BMA(mat_mass, tt0=T0, pl=true, wpm=wpmm, path_plt=path_plot)
# mat_mass = comp_mat_multigamma(ens, corr_pp[k], corr_pa[k], corr_aa[k])
# m_tm_ps, m_tm_vec = gevp_mass_BMA(mat_mass, tt0=T0, pl=true, wpm=wpmm, path_plt=path_plot)
m_tm_ps = meff.(corr_pp[k], plat, pl=true, wpm=wpmm)
m_tm_vec = meff.(corr_aa[k], plat, pl=true, wpm=wpmm)
# push!(m_tm_vec, uwreal([1.029, 0.001],"mhvec_D200"))
Mtotps[k] = m_tm_ps
Mtotvec[k] = m_tm_vec
......@@ -195,6 +232,18 @@ Mtotvec = Vector{Vector{uwreal}}(undef, length(ensinfo)) # required for decay co
obs_tm[k]["fpik"] = 0.0
end
end
##
# test D200
plat = [[85,105], [85,105], [85,105], [85,90], [88,98], [90,94], [85,90], [88,94], [90,94] ]
# test J501
plat = [[120,155], [120,155], [120,155], [125,155], [130,155], [132,145], [125,150], [130,145], [132,145] ]
# test E250
plat = [[115,150], [120,155], [120,155], [115,128], [130,155], [135,160], [115,128], [130,145], [135,160] ]
# test masses J500
plat = [[130,160], [140,160], [138,150], [140,160], [138,150] ]
dec_const_pcvc(corr_pp[1][5], plat[5], Mtotps[1][5], pl=true)
##
# Decay constants
close("all")
......@@ -206,11 +255,14 @@ for (k,ens) in enumerate(ensinfo)
# dec_ps = gevp_dec_BMA(mat_dec_ps, Mtotps[k], t0=T0, wpm=wpmm, pl=true, n=1, pseudo=true, wilson=false )
# dec_ps_test = get_f_tm(corr_pp[k], Mtotps[k], ens, path_plat_tm_dps, pl=true, wpm=wpmm )
path_plt_ps = joinpath(path_plot, "plateaus", ens.id, "fps")
dec_tm_ps = get_f_tm_BMA(corr_pp[k], Mtotps[k], ens, pl=true, wpm=wpmm, path_plt=path_plot, ps="fps")
path_plt_vec = joinpath(path_plot, "plateaus", ens.id, "fps")
dec_tm_vec = get_f_tm_BMA(corr_aa[k], Mtotvec[k], ens, pl=true, wpm=wpmm, path_plt=path_plot, ps="fvec")
# commented for testing D200
# path_plt_ps = joinpath(path_plot, "plateaus", ens.id, "fps")
# dec_tm_ps = get_f_tm_BMA(corr_pp[k], Mtotps[k], ens, pl=true, wpm=wpmm, path_plt=path_plot, ps="fps")
#
# path_plt_vec = joinpath(path_plot, "plateaus", ens.id, "fps")
# dec_tm_vec = get_f_tm_BMA(corr_aa[k], Mtotvec[k], ens, pl=true, wpm=wpmm, path_plt=path_plot, ps="fvec")
dec_tm_ps = dec_const_pcvc.(corr_pp[k], plat, Mtotps[k], pl=true, wpm=wpmm)
dec_tm_vec = dec_const.(corr_aa[k], plat, Mtotvec[k], pl=true, wpm=wpmm)
mu_list = gen_mulist[k]
mul, mus, muh = get_mu(mu_list, ens.deg)
......@@ -248,7 +300,6 @@ if mass_shift
# read dm computed by AS
rep = filter(x->occursin(ens.id, x), readdir(path_bdio_md,join=true))[1]
ens.deg ? Nmolt =3. : Nmolt = 1. # required for symmetric point ensembles if shifting s only
println(Nmolt)
deltam[k] = read_BDIO(rep, "ll_obs", "dm")[1] * Nmolt
# these two lines for dm computed locally by me
# fname = joinpath(path_bdio_md, string(ens.id,".bdio"))
......
src/gevp analysis/plot_paper.jl 0 → 100644
View file @ e56a3bed
using Base: String
using Base: @kwdef
using LaTeXStrings: length
using OrderedCollections
using juobs, BDIO, DelimitedFiles, ADerrors, LaTeXStrings, PyPlot
using ProgressMeter
#============= SET UP VARIABLES ===========#
rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["text.usetex"] = true
rcParams["mathtext.fontset"] = "cm"
rcParams["font.size"] =10
rcParams["axes.labelsize"] =22
rcParams["axes.titlesize"] = 18
plt.rc("text", usetex=true) # set to true once you install latex
#============ PATHS & INFO =================#
# path to data
const path_data_w = "/Users/alessandroconigli/Lattice/data/aux_obs_data/wilson" # required for t0 at finite a
const path_data_tm = "/Users/alessandroconigli/Lattice/data/charm_full_Dirac" # charm correlators
# path to plateau
const path_plat_w = "/Users/alessandroconigli/Lattice/data/aux_obs_data/plat_wilson.txt"
const path_plat_tm_mps = "/Users/alessandroconigli/Lattice/gevp-automated-analysis/plat/plat_tm_mps.txt"
const path_plat_tm_mvec = "/Users/alessandroconigli/Lattice/gevp-automated-analysis/plat/plat_tm_mvec.txt"
const path_plat_tm_dps = "/Users/alessandroconigli/Lattice/gevp-automated-analysis/plat/plat_tm_dps.txt"
const path_plat_tm_dvec = "/Users/alessandroconigli/Lattice/gevp-automated-analysis/plat/plat_tm_dvec.txt"
# path to aux obs
const path_rw = "/Users/alessandroconigli/Lattice/data/aux_obs_data/rwf" # reweighting factors
const path_md = "/Users/alessandroconigli/Lattice/data/aux_obs_data/md/" # pbp dat file for mass shift
# path to bdio
# const path_bdio_md = "/Users/alessandroconigli/Lattice/data/bdio_charm/md/LAT22/" # path to delta m computed by me
const path_bdio_md = "/Users/alessandroconigli/Lattice/data/bdio_charm/AS_ligth_obs/" # path to delta m computed by AS
const path_bdio_tm = "/Users/alessandroconigli/Lattice/data/bdio_charm/tm" # path where to store results
const path_bdio_tm_shifted = "/Users/alessandroconigli/Lattice/data/bdio_charm/tm_shifted" # path were to store shifted results
const path_plot = "/Users/alessandroconigli/MyDrive/phd/analysis/charm_gevp_2022/CharmPaper/plots" # path were to save plots
# ensembles to analyse
const ens_list = [ "J303"]
# reweighting and mass shift flags
const rwf = true
# println("Warning! mass_shift set to false ")
const mass_shift = false
const tau = 3 # gevp shift parameter
const T0 = 2
const TSM = false # set whether TSM is used or not
#@warning("\nTSM FLAG MODE: $(TSM)\n")
#=============== INCLUDES ===============#
include("types.jl")
include("tools.jl")
include("const.jl")
include("read_bdio.jl")
#=============== ENSEMBLE INFO FROM DATABASE ==================#
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], trunc_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
wpmm = Dict{String, Vector{Float64}}()
wpmm["H101"] = [-1.0, 2.0, -1.0, -1.0] #14.0*2.86]
wpmm["H102r002"] = [-1.0, 2.0, -1.0, -1.0] #14.0*2.86]
wpmm["H105r002"] = [-1.0, 2.0, -1.0, -1.0]# 14.0*2.86]
wpmm["H400"] = [-1.0, 2.0, -1.0, -1.0] #14.0*3.659]
wpmm["N202"] = [-1.0, 2.0, -1.0, -1.0] #14.0*5.164]
wpmm["N200"] = [-1.0, 2.0, -1.0, -1.0] #14.0*5.164]
wpmm["N203"] = [-1.0, 2.0, -1.0, -1.0] # 14.0*5.164]
#wpmm["D200"] = [-1.0, 4.0, -1.0, 14.0*5.164]
wpmm["D200"] = [-1.0, 2.0, -1.0, -1.0]
wpmm["N300r002"] = [-1.0, 2.0, -1.0, -1.0 ] #14.0*8.595]
wpmm["N302"] = [-1.0, 2.0, -1.0, -1.0 ] #14.0*8.595]
wpmm["J303"] = [-1.0, 2.0, -1.0, -1.0 ] #14.0*8.595]
wpmm["J500"] = [-1.0, 2.0, -1.0, -1.0 ] #14.0*8.595]
wpmm["J501"] = [-1.0, 2.0, -1.0, -1.0 ] #14.0*8.595]
wpmm["E250"] = [-1.0, 2.0, -1.0, -1.0 ] #14.0*8.595]
##
#=============== LOAD DATA AND STORE t0 ===============#
corr_pp = Vector{Vector{juobs.Corr}}(undef, length(ens_list))
corr_pa = Vector{Vector{juobs.Corr}}(undef, length(ens_list))
corr_aa = Vector{Vector{juobs.Corr}}(undef, length(ens_list))
yw_store = Vector{Matrix{uwreal}}(undef, length(ens_list)) # store YW info from t0
w_store = Vector{uwreal}(undef, length(ens_list)) # store W info from t0
pp_tmW = Vector{Vector{Vector{uwreal}}}(undef, length(ens_list)) # store W info from corr_pp
gen_mulist = Vector{Vector{Vector{Float64}}}(undef, length(ens_list)) # store mu values
# obs storage
obs_tm = Vector(undef, length(ensinfo))
@time begin
@showprogress for (k,ens) in enumerate(ensinfo)
println("\n Ensemble: ", ens.id)
corr_pp[k], pp_tmW[k], _ = get_corr(path_data_tm, ens, "G5", "G5", rw=rwf, path_rw=path_rw, tsm=TSM)
gen_mulist[k] = getfield.(corr_pp[k], :mu)
end
end
##
mlh, data_lh = juobs.meff(corr_pp[1][4], [127,150], data=true)
msh, data_sh = juobs.meff(corr_pp[1][5], [135,175], data=true)
mhh, data_hh = juobs.meff(corr_pp[1][6], [135,175], data=true)
data_lh ./= mlh ; uwerr.(data_lh)
data_sh ./= msh ; uwerr.(data_sh)
data_hh ./= mhh ; uwerr.(data_hh)
##
fig= figure("pyplot_subplot_column",figsize=(12,4))
subplots_adjust(wspace=0.0)
xx = collect(1:length(data_lh))
subplot(131)
errorbar(xx, value.(data_lh), err.(data_lh), fmt="*", capsize=2, color="navy", mfc="none", label=L"$D$")
fill_between(collect(127:150) , 1 - err(mlh), 1 + err(mlh), color="red", alpha=0.5)
xlabel(L"$x_0/a$")
ylabel(L"$m^{\mathrm{eff}}(x_0)/m_{\mathrm{PS}}$")
xlim(96,190)
ylim(0.98,1.02)
legend()
aa = subplot(132)
setp(aa.get_yticklabels(), visible=false )
errorbar(xx, value.(data_sh), err.(data_sh), fmt="*", capsize=2, color="orange", mfc="none", label=L"$D_s$")
fill_between(collect(135:175) , 1 - err(msh), 1 + err(msh), color="red", alpha=0.5)
xlabel(L"$x_0/a$")
xlim(96,190)
ylim(0.98,1.02)
legend()
aa=subplot(133)
setp(aa.get_yticklabels(), visible=false )
errorbar(xx, value.(data_hh), err.(data_hh), fmt="*", capsize=2, color="darkgreen", mfc="none", label=L"$\eta_c$")
fill_between(collect(135:175) , 1 - err(mhh), 1 + err(mhh), color="red", alpha=0.5)
xlabel(L"$x_0/a$")
xlim(96,190)
ylim(0.98,1.02)
legend()
ax = fig.add_axes([0.78,0.3,0.15,0.2])
errorbar(xx, value.(data_hh), err.(data_hh), fmt="*", capsize=2, color="darkgreen", mfc="none")
fill_between(collect(135:175) , 1 - err(mhh), 1 + err(mhh), color="red", alpha=0.5)
xlim(120,180)
ylim(0.999,1.001)
tight_layout()
display(gcf())
savefig("/Users/alessandroconigli/Desktop/meff_comparison.pdf")
close("all")
\ No newline at end of file
src/gevp analysis/test_read_multiple_TSM_chunks.jl 0 → 100644
View file @ e56a3bed
using Base: String
using Base: @kwdef
using LaTeXStrings: length
using OrderedCollections
using juobs, BDIO, DelimitedFiles, ADerrors, LaTeXStrings, PyPlot
using ProgressMeter
#============= SET UP VARIABLES ===========#
rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["text.usetex"] = true
rcParams["mathtext.fontset"] = "cm"
rcParams["font.size"] =10
rcParams["axes.labelsize"] =22
rcParams["axes.titlesize"] = 18
plt.rc("text", usetex=true) # set to true once you install latex
#============ PATHS & INFO =================#
# path to data
const path_data_w = "/Users/alessandroconigli/Lattice/data/aux_obs_data/wilson" # required for t0 at finite a
const path_data_tm = "/Users/alessandroconigli/Lattice/data/charm_full_Dirac" # charm correlators
# path to aux obs
const path_rw = "/Users/alessandroconigli/Lattice/data/aux_obs_data/rwf" # reweighting factors
const path_md = "/Users/alessandroconigli/Lattice/data/aux_obs_data/md/" # pbp dat file for mass shift
const ens_list = [ "E250"]
# reweighting and mass shift flags
const rwf = false
#=============== INCLUDES ===============#
include("types.jl")
include("tools.jl")
include("const.jl")
include("read_bdio.jl")
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], trunc_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
##
# testing for E250 (in chunk2 and chunk2bis there's an overlap)
db = [ "sloppy_chunk0","sloppy_chunk1", "sloppy_chunk2", "sloppy_chunk3", "sloppy_chunk5", "sloppy_chunk5"]
pp = joinpath(path_data_tm, "E250")
aux = filter(x-> occursin("sloppy_chunk", x), readdir(pp,join=true))
data_1 = juobs.read_mesons(aux[1], "G5","G5", legacy=false) # 1 -> 196
data_2 = juobs.read_mesons(aux[2], "G5","G5", legacy=false) # 201 -> 296 # trunc at 20 elem
data_2b = juobs.read_mesons(aux[3], "G5","G5", legacy=false) # 301 -> 396
data_3 = juobs.read_mesons(aux[4], "G5","G5", legacy=false) # 401 -> 596
data_4 = juobs.read_mesons(aux[5], "G5","G5", legacy=false) # 601 -> 796
data_5 = juobs.read_mesons(aux[6], "G5","G5", legacy=false) # 801 -> 1006
juobs.concat_data!(data_1, data_2)
juobs.concat_data!(data_1, data_2b)
juobs.concat_data!(data_1, data_3)
juobs.concat_data!(data_1, data_4)
juobs.concat_data!(data_1, data_5)
##
data_sl = read_data_sloppy_multichunks(path_data_tm, "E250", "G5", "G5", legacy=false)
data_corr = read_data_correction(path_data_tm, "E250", "G5", "G5", legacy=false)
corr_pp, _ = get_corr(path_data_tm, ensinfo[1], "G5", "G5", rw=false, path_rw=path_rw, tsm=true)
mm = meff.(corr_pp, fill([130, 150],9), pl=true )
function get_corr_TSM_multichunks(path::String, ens::EnsInfo, g1::String="G5", g2::String="G5"; rw::Bool=false, legacy::Bool=false, path_rw::String="", tsm::Bool=false, ex::Bool=false)
if path_rw == ""
p_rw = path
else
p_rw = path_rw
end
aux1 = read_data_sloppy_multichunks(path, ens.id, g1, g2, legacy=legacy)
aux2 = read_data_correction(path, ens.id, g1, g2, legacy=legacy)
try # single replica
# sloppy
vcfg = getfield(aux1[1], :vcfg)
delta = vcfg[1][2] - vcfg[1][1]
cut_trunc = Int64.((ens.trunc .-1)./delta) .+ 1
truncate_data!(aux1, cut_trunc )
# correction
vcfg = getfield(aux2[1], :vcfg)
delta = vcfg[1][2] - vcfg[1][1]
cut_trunc = Int64.((ens.trunc .-1)./delta) .+ 1
truncate_data!(aux2, cut_trunc)
catch # multiple replica
# sloppy
vcfg = getfield.(aux1[1], :vcfg)
delta = vcfg[1][2] - vcfg[1][1]
cut_trunc = Int64.((ens.trunc .-1)./delta) .+ 1
truncate_data!(aux1, cut_trunc )
# correction
vcfg = getfield.(aux2[1], :vcfg)
delta = vcfg[1][2] - vcfg[1][1]
cut_trunc = Int64.((ens.trunc .-1)./delta) .+ 1
truncate_data!(aux2, cut_trunc)
end
if !rw
obs = corr_obs_TSM.(aux1, aux2, L=ens.L, info=true)#, idm_corr=collect(1:ens.idm_corr) )
return (getindex.(obs, 1), getindex.(obs, 2))
else
obs = corr_obs_TSM.(aux1, aux2, L=ens.L, rw=read_rw(p_rw, ens.id), info=true) #, idm_corr=collect(1:ens.idm_corr))
return (getindex.(obs, 1), getindex.(obs, 2), getindex.(obs, 3))
end
end
## testing apply_rw with gaps
data_sl = read_data_sloppy_multichunks(path_data_tm, "E300", "G5", "G5", legacy=false);
data_corr = read_data_correction(path_data_tm, "E300", "G5", "G5", legacy=false);
L = ensinfo[1].L;
idd = [ensinfo[1].id];
real_data = true;
vcfg_sl = getfield(data_sl[1], :vcfg);
replica_sl = Int64(maximum(vcfg_sl))
vcfg_corr = getfield(data_corr[1], :vcfg);
replica_corr = Int64(maximum(vcfg_corr));
#truncate data
nms = sum(replica_sl); # assuming vcfg_sl >= vcfg_corr
delta_vcfg = vcfg_sl[1][2] - vcfg_sl[1][1];
idm_sl = collect(1:delta_vcfg:replica_sl[1])
for ii in eachindex(data_sl[1])[2:end]
delta_vcfg = vcfg_sl[ii][2] - vcfg_sl[ii][1]
aux = collect(1:delta_vcfg:replica_sl[ii]) .+ replica_sl[ii-1]
append!(idm_sl, aux)
end
delta_vcfg = vcfg_corr[1][2] - vcfg_corr[1][1]
idm_corr = collect(1:delta_vcfg:replica_corr[1])
for ii in eachindex(data_corr[1])[2:end]
delta_vcfg = vcfg_corr[ii][2] - vcfg_corr[ii][1]
aux = collect(1:delta_vcfg:replica_corr[ii]) .+ replica_corr[ii-1]
append!(idm_corr, aux)
end
data1 = real_data ? getfield.(data_sl[1], :re_data) ./ L^3 : getfield.(data_sl, :im_data) ./ L^3
data2 = real_data ? getfield.(data_corr[1], :re_data) ./ L^3 : getfield.(data_corr, :im_data) ./ L^3
rw = read_rw(path_rw, ensinfo[1].id)
data1_r, W = juobs.apply_rw(data1, rw, vcfg_sl)
data2_r, W = juobs.apply_rw(data2, rw, vcfg_corr)
tmp1 = vcat(data1_r...)
tmp2 = vcat(data2_r...)
tmpW = vcat(W...)
nt = size(tmp1, 2)
ow1 = [uwreal(tmp1[:, x0], idd[1], replica_sl, idm_sl, nms ) for x0 = 1:nt]
ow2 = [uwreal(tmp2[:, x0], idd[1], replica_sl, idm_corr, nms) for x0 = 1:nt]
W_obs = uwreal(tmpW, idd[1], replica_sl, idm_sl, nms)
obs1 = [ow1[x0] / W_obs for x0 = 1:nt]
obs2 = [ow2[x0] / W_obs for x0 = 1:nt]
pp = juobs.Corr(obs1 + obs2, data_sl[1])
\ No newline at end of file
src/gevp analysis/tools.jl
View file @ e56a3bed
......@@ -23,6 +23,32 @@ function read_data_sloppy(path::String, ens::String, g1::String="G5", g2::String
error("mesons.dat file not found for ensemble ", ens, " in path ", p)
end
end
function read_data_sloppy_multichunks(path::String, ens::String, g1::String="G5", g2::String="G5"; legacy::Bool=false)
p = joinpath(path, ens)
db = Dict("J500" => ["sloppy_chunk0", "sloppy_chunk1"],
"E300" => ["sloppy_chunk0", "sloppy_chunk1", "sloppy_chunk2" ],
"N302" => ["sloppy"], "J501" => ["sloppy"],
"E250" => [ "sloppy_chunk0","sloppy_chunk1", "sloppy_chunk2", "sloppy_chunk3", "sloppy_chunk4", "sloppy_chunk5"]
)
store_cdata_aux = []
for rep in db[ens]
aux = filter(x-> occursin(rep, x), readdir(p,join=true))
data_chunk = juobs.read_mesons(aux, g1, g2, legacy=legacy, id=ens)
if ens == "E250" && rep == "sloppy_chunk1"
juobs.truncate_data!.(data_chunk, 20)
end
push!(store_cdata_aux, data_chunk)
end
for k in 2:length(store_cdata_aux)
juobs.concat_data!(store_cdata_aux[1], store_cdata_aux[k])
end
if length(store_cdata_aux[1][1]) == 1
return getindex.(store_cdata_aux[1], 1)
else
return store_cdata_aux[1]
end
end
# read correction TSM
function read_data_correction(path::String, ens::String, g1::String="G5", g2::String="G5"; legacy::Bool=false)
p = joinpath(path, ens)
......@@ -35,8 +61,12 @@ function read_data_correction(path::String, ens::String, g1::String="G5", g2::St
end
function read_pbp(path::String, ens::String)
p = joinpath(path, ens)
rep = filter(x->occursin("pbp.dat", x), readdir(p, join=true))
rep = filter(x->occursin("pbp", x), readdir(p, join=true))
if length(rep)!=0
if ens == "D200"
md_aux = read_md.(rep)
return [hcat(md_aux...)]
end
return read_md.(rep)
else
error("pbp.dat file not found for ensemble ", ens, " in path ", p)
......@@ -48,6 +78,12 @@ function read_rw(path::String, ens::String)
if ens == "J500"
return [read_ms1(rep[1]), read_ms1(rep[2], v="1.4")]
end
if ens == "J501"
return [read_ms1(rep[1]), read_ms1(rep[2], v="1.4") ]
end
if ens == "E250"
return read_ms1(rep[1], v="1.6")
end
if length(rep)!=0
try
length(rep) == 1 ? (return read_ms1(rep[1])) : (return read_ms1.(rep))
......@@ -87,21 +123,36 @@ function get_corr(path::String, ens::EnsInfo, g1::String="G5", g2::String="G5";
return (getindex.(obs, 1), getindex.(obs, 2), getindex.(obs, 3))
end
else
aux1 = read_data_sloppy(path, ens.id, g1, g2, legacy=legacy)
if ens.id in ["J500", "J501", "E250", "E300"]
aux1 = read_data_sloppy_multichunks(path, ens.id, g1, g2, legacy=legacy)
else
aux1 = read_data_sloppy(path, ens.id, g1, g2, legacy=legacy)
end
aux2 = read_data_correction(path, ens.id, g1, g2, legacy=legacy)
if !isnothing(ens.trunc)
try # single replica
vcfg = getfield(aux1[1], :vcfg)
delta = vcfg[1][2] - vcfg[1][1]
cut_trunc = Int64.((ens.trunc .-1)./delta) .+ 1
truncate_data!(aux1, cut_trunc )
truncate_data!(aux2, cut_trunc)
if ens.id == "E300"
truncate_data!(aux1, 455)
else
vcfg = getfield(aux1[1], :vcfg)
delta = vcfg[2] - vcfg[1]
cut_trunc = Int64.((ens.trunc .-1)./delta) .+ 1
println(cut_trunc)
truncate_data!(aux1, cut_trunc )
vcfg2 = getfield(aux2[1], :vcfg)
delta2 = vcfg2[2] - vcfg2[1]
cut_trunc2 = Int64.((ens.trunc .-1)./delta2) .+ 1
truncate_data!(aux2, cut_trunc2)
end
catch # multiple replica
vcfg = getfield.(aux1[1], :vcfg)
delta = vcfg[1][2] - vcfg[1][1]
cut_trunc = Int64.((ens.trunc .-1)./delta) .+ 1
truncate_data!(aux1, cut_trunc )
truncate_data!(aux2, cut_trunc)
vcfg2 = getfield.(aux2[1], :vcfg)
delta2 = vcfg2[1][2] - vcfg2[1][1]
cut_trunc2 = Int64.((ens.trunc .-1)./delta2) .+ 1
truncate_data!(aux2, cut_trunc2)
end
end
......@@ -420,12 +471,9 @@ function gevp_mass_BMA(mat_obs::Vector{MatInfo}; tt0::Int64=2,pl::Bool=true, wpm
# evecs = getall_eig(mat, t0)
entot = juobs.energies(evals, wpm=wpm)
@. gs_model(x,p) = p[1] #+ p[2] * exp(-(p[3])*x)
param = 1
tminaux = round.(Int64, collect(18:30) * sqrt(value(t0(mat_obs[k].ensinfo.beta) / t0(3.7))))
tminaux = round.(Int64, collect(18:38) * sqrt(value(t0(mat_obs[k].ensinfo.beta) / t0(3.7))))
T = length(entot[1])
tmaxaux = round.(Int64, collect(8:20) * sqrt(value(t0(mat_obs[k].ensinfo.beta) / t0(3.7))))
tmaxaux = round.(Int64, collect(10:24) * sqrt(value(t0(mat_obs[k].ensinfo.beta) / t0(3.7))))
tmax = T .- collect(reverse(tmaxaux))
tmax = tmax[1:3:end]
tmin = tminaux[1:2:end]
......@@ -436,8 +484,29 @@ function gevp_mass_BMA(mat_obs::Vector{MatInfo}; tt0::Int64=2,pl::Bool=true, wpm
else
pp1 = pp2 = nothing
end
en0[k], syst1, all_res1, FitDict1 = juobs.bayesian_av(gs_model, entot[1], tmin, tmax, param, pl=pl, wpm=wpm, data=true, label=L"$m_{\mathrm{ps}}$", path_plt=pp1, plt_title=t)
en1[k], syst2, all_res2, FitDict2 = juobs.bayesian_av(gs_model, entot[2], tmin, tmax, param, pl=pl, wpm=wpm, data=true, label=L"$m_{\mathrm{v}}$", path_plt=pp2, plt_title=t)
if mat_obs[1].ensinfo.id == "D200"
@. gs_model(x,p) = p[1] #+ p[2] * exp(-(p[3])*x)
param = 1
en0[k], syst1, all_res1, FitDict1 = juobs.bayesian_av(gs_model, entot[1], tmin, tmax, param, pl=pl, wpm=wpm, data=true, label=L"$m_{\mathrm{ps}}$", path_plt=pp1, plt_title=t)
try
en1[k], syst2, all_res2, FitDict2 = juobs.bayesian_av(gs_model, entot[2], tmin, tmax, param, pl=pl, wpm=wpm, data=true, label=L"$m_{\mathrm{v}}$", path_plt=pp2, plt_title=t)
catch
@warn "Vector mass failed"
en1[k] = uwreal(1.0)
end
else
@. gs_model1(x,p) = p[1] #+ p[2] * exp(-(p[3])*x)
param = 1
en0[k], syst1, all_res1, FitDict1 = juobs.bayesian_av(gs_model1, entot[1], tmin, tmax, param, pl=pl, wpm=wpm, data=true, label=L"$m_{\mathrm{ps}}$", path_plt=pp1, plt_title=t)
try
en1[k], syst2, all_res2, FitDict2 = juobs.bayesian_av(gs_model1, entot[2], tmin, tmax, param, pl=pl, wpm=wpm, data=true, label=L"$m_{\mathrm{v}}$", path_plt=pp2, plt_title=t)
catch
@warn "Vector mass failed"
en1[k] = uwreal(1.0)
end
end
close("all")
if pl
......@@ -773,7 +842,7 @@ function get_f_tm_BMA(corr::Vector{juobs.Corr}, m::Vector{uwreal}, ens::EnsInfo;
tminaux = round.(Int64, collect(15:32) * sqrt(value(t0(ens.beta) / t0(3.7))))
T = length(corr[1].obs)
tmaxaux = round.(Int64, collect(20:35) * sqrt(value(t0(ens.beta) / t0(3.7))))
tmaxaux = round.(Int64, collect(10:25) * sqrt(value(t0(ens.beta) / t0(3.7))))
tmax = T .- collect(reverse(tmaxaux))
tmax = tmax[1:3:end]
tmin = tminaux[1:2:end] .+ corr[1].y0
......@@ -795,61 +864,68 @@ function get_f_tm_BMA(corr::Vector{juobs.Corr}, m::Vector{uwreal}, ens::EnsInfo;
ll = ps == "fps" ? L"$f_{\mathrm{ps}}$" : L"$f_{\mathrm{v}}$"
name_bma = "f$(k)_BMA.pdf"
isnothing(path_plt) ? tt=nothing : tt = joinpath(path_plt,name_bma)
F[k], syst, all_res, FitDict = juobs.bayesian_av(const_fit, data, tmin, tmax, 1, pl=pl, wpm=wpm, data=true, path_plt=tt, plt_title=t, label=ll )
try
F[k], syst, all_res, FitDict = juobs.bayesian_av(const_fit, data, tmin, tmax, 1, pl=pl, wpm=wpm, data=true, path_plt=tt, plt_title=t, label=ll )
if pl
try
fig = figure(figsize=(10,7.5))
modstot = FitDict["mods"]
AIC = FitDict["AIC"]
Wtot = FitDict["W"]
IDX = partialsortperm(AIC, 1:3)
mods = modstot[IDX]
main_res = all_res[IDX]
W = Wtot[IDX]
isnothing(wpm) ? uwerr(data) : [uwerr(data[i],wpm) for i in eachindex(data)]
vv = value.(data)[corr[1].y0:end]
ee = err.(data)[corr[1].y0:end]
xx = collect(0:length(vv)-1) .+ corr[1].y0
errorbar(xx, vv, ee, fmt="d", mfc="none", color="black", capsize=2)
CC = ["forestgreen", "royalblue", "tomato"]
for (ii,trange) in enumerate(mods)
isnothing(wpm) ? uwerr(main_res[ii]) : uwerr(main_res[ii], wpm)
vaux = value(main_res[ii])
eaux = err(main_res[ii])
regex = r"([0-9]+),([0-9]+)"
aux = match(regex, trange)
fill_between(parse(Float64,aux[1]):parse(Float64,aux[2]), vaux-eaux, vaux+eaux, color=CC[ii], alpha=0.5, label=round(W[ii], digits=3))
end
if pl
try
fig = figure(figsize=(10,7.5))
modstot = FitDict["mods"]
AIC = FitDict["AIC"]
Wtot = FitDict["W"]
IDX = partialsortperm(AIC, 1:3)
mods = modstot[IDX]
main_res = all_res[IDX]
W = Wtot[IDX]
isnothing(wpm) ? uwerr(data) : [uwerr(data[i],wpm) for i in eachindex(data)]
vv = value.(data)[corr[1].y0:end]
ee = err.(data)[corr[1].y0:end]
xx = collect(0:length(vv)-1) .+ corr[1].y0
errorbar(xx, vv, ee, fmt="d", mfc="none", color="black", capsize=2)
CC = ["forestgreen", "royalblue", "tomato"]
for (ii,trange) in enumerate(mods)
isnothing(wpm) ? uwerr(main_res[ii]) : uwerr(main_res[ii], wpm)
vaux = value(main_res[ii])
eaux = err(main_res[ii])
regex = r"([0-9]+),([0-9]+)"
aux = match(regex, trange)
fill_between(parse(Float64,aux[1]):parse(Float64,aux[2]), vaux-eaux, vaux+eaux, color=CC[ii], alpha=0.5, label=round(W[ii], digits=3))
end
axhline(y=value(main_res[1])-14*err(main_res[1]), xmin=(tmin[1]-corr[1].y0)/length(xx), xmax=(tmin[end]-corr[1].y0)/length(xx), ls="dashed", color="green")
axhline(y=value(main_res[1])-14*err(main_res[1]), xmin=(tmax[1]-corr[1].y0)/length(xx), xmax=(tmax[end]-corr[1].y0)/length(xx), ls="dashed", color="red")
text(tmin[1],value(main_res[1])-13*err(main_res[1]), L"$\{t_{\mathrm{min}}\}$" )
text(tmax[1],value(main_res[1])-13*err(main_res[1]), L"$\{t_{\mathrm{max}}\}$" )
axhline(y=value(main_res[1])-14*err(main_res[1]), xmin=(tmin[1]-corr[1].y0)/length(xx), xmax=(tmin[end]-corr[1].y0)/length(xx), ls="dashed", color="green")
axhline(y=value(main_res[1])-14*err(main_res[1]), xmin=(tmax[1]-corr[1].y0)/length(xx), xmax=(tmax[end]-corr[1].y0)/length(xx), ls="dashed", color="red")
text(tmin[1],value(main_res[1])-13*err(main_res[1]), L"$\{t_{\mathrm{min}}\}$" )
text(tmax[1],value(main_res[1])-13*err(main_res[1]), L"$\{t_{\mathrm{max}}\}$" )
# axvline(tmin[1], ls="dashed", color="green")
# axvline(tmax[end], ls="dashed", color="red")
legend(title="Weights")
title(t)
xlabel(L"$t/a$")
ylabel(ll)
xlim(xx[1], xx[end])
ylim(value(main_res[1]) -20*err(main_res[1]), value(main_res[1]) +20*err(main_res[1]))
twiny()
xlim(xx[1], xx[end])
fm_time = round.((xx .- corr[1].y0)[1:8:end] .* value(a(ens.beta)), digits=3)
xticks(xx[1:8:end], fm_time)
# xticks()
xlabel(L"$t-t_0\ [\mathrm{fm}]$", fontsize=14)
display(fig)
if !isnothing(path_plt)
savefig(joinpath(path_plt, "f$(k)_plateau.pdf"))
# axvline(tmin[1], ls="dashed", color="green")
# axvline(tmax[end], ls="dashed", color="red")
legend(title="Weights")
title(t)
xlabel(L"$t/a$")
ylabel(ll)
xlim(xx[1], xx[end])
ylim(value(main_res[1]) -20*err(main_res[1]), value(main_res[1]) +20*err(main_res[1]))
twiny()
xlim(xx[1], xx[end])
fm_time = round.((xx .- corr[1].y0)[1:8:end] .* value(a(ens.beta)), digits=3)
xticks(xx[1:8:end], fm_time)
# xticks()
xlabel(L"$t-t_0\ [\mathrm{fm}]$", fontsize=14)
display(fig)
if !isnothing(path_plt)
savefig(joinpath(path_plt, "f$(k)_plateau.pdf"))
end
close()
catch
@warn "plot failed"
end
close()
catch
println("plot failed")
end
catch
@warn "Decay failed"
F[k] = uwreal(0.0)
end
end
......@@ -1044,8 +1120,8 @@ function aic_weight(aic)
w_aux ./= sum(w_aux)
idx_in = sortperm(w_aux, rev=true)
cum_w = cumsum(w_aux[idx_in])
println(cum_w[end])
stop = findfirst(x->x>=0.999, cum_w)
#stop = findfirst(x->x>=0.999, cum_w)
stop = length(aic)
idx = sort(idx_in[1:stop])
w = w_aux[idx] ./ sum(w_aux[idx])
# w = w_aux ./ sum(w_aux)
......
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