Skip to content

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

decay constants extraction added

parent 79a34cd1
Hide whitespace changes
Inline Side-by-side
Showing with 97 additions and 183 deletions
gevp.jl → src/gevp.jl
View file @ 43d0b8ef
using Revise, juobs, BDIO, ADerrors, DelimitedFiles, PyPlot, LaTeXStrings, LsqFit
using Revise, juobs, BDIO, DelimitedFiles, ADerrors, PyPlot, LaTeXStrings, LsqFit
#============= SET UP VARIABLES ===========#
const path_data = "/Users/ale/Desktop/data"
const path_plat = "/Users/ale/automation/plat.txt"
const path_results = "/Users/ale/Desktop/results_gevp"
const ensembles = [ "J303", "N300", "H400"]
const ensembles = [ "J303"]
const sector = Dict("ll"=>false, "ls"=>false, "lh"=>true, "ss"=>false, "sh"=>true, "hh"=>true )
const tau = 4
const _t0 = 2
......@@ -67,9 +67,11 @@ println("\n Computing observables \n")
println("\n Ensemble: ", ens.id)
en0_ps = gevp_mass(matinfo[i], t0=_t0, delta_t_in=range_t_fit, wpm=wpmm)
en0_vec = gevp_mass(matinfo_vec[i], t0=_t0, delta_t_in=range_t_fit, wpm=wpmm)
println( "\nground ps state: ", en0_ps)
println( "\nground vec state: ", en0_vec)
ensobs[i] = EnsObs(ens, gen_mulist[i], en0_ps, en0_vec)
println( "\nground ps mass: ", en0_ps)
println( "\nground vec mass: ", en0_vec)
dec0_ps = gevp_dec(matinfo[i], en0_ps, t0=_t0, delta_t_in=range_t_fit, wpm=wpmm)
dec0_vec = gevp_dec(matinfo_vec[i], en0_vec, t0=_t0, delta_t_in=range_t_fit, wpm=wpmm)
ensobs[i] = EnsObs(ens, gen_mulist[i], en0_ps, dec0_ps, en0_vec, dec0_vec)
if ens.deg
ensobs[i].m_ls = ensobs[i].m_ll
ensobs[i].m_ls_vec = ensobs[i].m_ll_vec
......@@ -77,6 +79,12 @@ println("\n Computing observables \n")
ensobs[i].m_ss_vec = ensobs[i].m_ll_vec
ensobs[i].m_sh = ensobs[i].m_lh
ensobs[i].m_sh_vec = ensobs[i].m_lh_vec
ensobs[i].f_ls = ensobs[i].f_ll
ensobs[i].f_ls_vec = ensobs[i].f_ll_vec
ensobs[i].f_ss = ensobs[i].f_ll
ensobs[i].f_ss_vec = ensobs[i].f_ll_vec
ensobs[i].f_sh = ensobs[i].f_lh
ensobs[i].f_sh_vec = ensobs[i].f_lh_vec
end
println("Time:")
end
......@@ -119,7 +127,7 @@ println("\n Matching Procedure \n")
par, chi2exp = lin_fit(muh, ens.m_lh_vec)
ens.m_lh_vec_match = y_lin_fit(par, value(ens.muh_target))
uwerr(ens.m_lh_vec_match)
#=
# f_lh
par, chi2exp = lin_fit(muh, ens.f_lh)
ens.f_lh_match = y_lin_fit(par, value(ens.muh_target))
......@@ -128,7 +136,7 @@ println("\n Matching Procedure \n")
par, chi2exp = lin_fit(muh, ens.f_lh_vec)
ens.f_lh_vec_match = y_lin_fit(par, value(ens.muh_target))
uwerr(ens.f_lh_vec_match)
=#
end
# interpolate sh
if sector["sh"]
......@@ -146,7 +154,7 @@ println("\n Matching Procedure \n")
par, chi2exp = lin_fit(muh, ens.m_sh_vec)
ens.m_sh_vec_match = y_lin_fit(par, value(ens.muh_target))
uwerr(ens.m_sh_vec_match)
#=
# f_sh
par, chi2exp = lin_fit(muh, ens.f_sh)
ens.f_sh_match = y_lin_fit(par, value(ens.muh_target))
......@@ -155,7 +163,7 @@ println("\n Matching Procedure \n")
par, chi2exp = lin_fit(muh, ens.f_sh_vec)
ens.f_sh_vec_match = y_lin_fit(par, value(ens.muh_target))
uwerr(ens.f_sh_vec_match)
=#
end
end
# interpolate hh
......@@ -168,7 +176,7 @@ println("\n Matching Procedure \n")
par, chi2exp = lin_fit(muh, ens.m_hh_vec)
ens.m_hh_vec_match = y_lin_fit(par, value(ens.muh_target))
uwerr(ens.m_hh_vec_match)
#=
# f_hh
par, chi2exp = lin_fit(muh, ens.f_hh)
ens.f_hh_match = y_lin_fit(par, value(ens.muh_target))
......@@ -177,7 +185,7 @@ println("\n Matching Procedure \n")
par, chi2exp = lin_fit(muh, ens.f_hh_vec)
ens.f_hh_vec_match = y_lin_fit(par, value(ens.muh_target))
uwerr(ens.f_hh_vec_match)
=#
end
end
println(" Total time:")
......@@ -212,7 +220,7 @@ for i in 1:length(ensobs)
name = string(ens.ensinfo.id,"_match_mD_mDvec.pdf")
savefig(joinpath(path_plot, name))
close()
#=
# decays
figure()
title(ens.ensinfo.id)
......@@ -229,7 +237,7 @@ for i in 1:length(ensobs)
name = string(ens.ensinfo.id,"_match_fD_fDvec.pdf")
savefig(joinpath(path_plot, name))
close()
=#
end
if sector["sh"] && !ens.ensinfo.deg
# masses
......@@ -248,7 +256,7 @@ for i in 1:length(ensobs)
name = string(ens.ensinfo.id,"_match_mDs_mDsvec.pdf")
savefig(joinpath(path_plot, name))
close()
#=
# decays
figure()
title(ens.ensinfo.id)
......@@ -265,7 +273,7 @@ for i in 1:length(ensobs)
name = string(ens.ensinfo.id,"_match_fDs_fDsvec.pdf")
savefig(joinpath(path_plot, name))
close()
=#
end
if sector["hh"]
# masses
......@@ -285,7 +293,7 @@ for i in 1:length(ensobs)
savefig(joinpath(path_plot, name))
close()
#=
# decays
figure()
title(ens.ensinfo.id)
......@@ -302,7 +310,7 @@ for i in 1:length(ensobs)
name = string(ens.ensinfo.id,"_match_fetac_fJpsi.pdf")
savefig(joinpath(path_plot, name))
close()
=#
end
end
......@@ -319,20 +327,20 @@ uwerr.(muc)
# ll
m_pi = Vector{uwreal}(undef, length(ensembles))
m_rho =Vector{uwreal}(undef, length(ensembles))
#=
f_pi = Vector{uwreal}(undef, length(ensembles))
f_rho =Vector{uwreal}(undef, length(ensembles))
=#
if sector["ll"]
try
m_pi .= getfield.(ensobs, :m_ll) .* sqrt.(8 * getfield.(ensobs,:t0))
m_rho .= getfield.(ensobs, :m_ll_vec) .* sqrt.(8 * getfield.(ensobs,:t0))
#f_pi .= getfield.(ensobs, :f_ll) .* sqrt.(8 * getfield.(ensobs,:t0))
#f_rho .= getfield.(ensobs, :f_ll_vec) .* sqrt.(8 * getfield.(ensobs,:t0))
f_pi .= getfield.(ensobs, :f_ll) .* sqrt.(8 * getfield.(ensobs,:t0))
f_rho .= getfield.(ensobs, :f_ll_vec) .* sqrt.(8 * getfield.(ensobs,:t0))
uwerr.(m_pi)
uwerr.(m_rho)
#uwerr.(f_pi)
#uwerr.(f_rho)
uwerr.(f_pi)
uwerr.(f_rho)
catch
error("Check ensobs[i].mu_list to see whether all the ensembles have a light-light sector.")
end
......@@ -340,18 +348,18 @@ end
# ls
m_k = Vector{uwreal}(undef, length(ensembles))
m_k_star =Vector{uwreal}(undef, length(ensembles))
#f_k= Vector{uwreal}(undef, length(ensembles))
#f_k_star =Vector{uwreal}(undef, length(ensembles))
f_k= Vector{uwreal}(undef, length(ensembles))
f_k_star =Vector{uwreal}(undef, length(ensembles))
if sector["ls"]
try
m_k .= getfield.(ensobs, :m_ls) .* sqrt.(8 * getfield.(ensobs,:t0))
m_k_star .= getfield.(ensobs, :m_ls_vec) .* sqrt.(8 * getfield.(ensobs,:t0))
#f_k .= getfield.(ensobs, :f_ls) .* sqrt.(8 * getfield.(ensobs,:t0))
#f_k_star .= getfield.(ensobs, :f_ls_vec) .* sqrt.(8 * getfield.(ensobs,:t0))
f_k .= getfield.(ensobs, :f_ls) .* sqrt.(8 * getfield.(ensobs,:t0))
f_k_star .= getfield.(ensobs, :f_ls_vec) .* sqrt.(8 * getfield.(ensobs,:t0))
uwerr.(m_k)
uwerr.(m_k_star)
#uwerr.(f_k)
#uwerr.(f_k_star)
uwerr.(f_k)
uwerr.(f_k_star)
catch
error("Check ensobs[i].mu_list to see whether all the ensembles have a light-stange sector.")
end
......@@ -359,18 +367,18 @@ end
# lh
m_D = Vector{uwreal}(undef, length(ensembles))
m_D_star =Vector{uwreal}(undef, length(ensembles))
#f_D= Vector{uwreal}(undef, length(ensembles))
#f_D_star =Vector{uwreal}(undef, length(ensembles))
f_D= Vector{uwreal}(undef, length(ensembles))
f_D_star =Vector{uwreal}(undef, length(ensembles))
if sector["lh"]
try
m_D .= getfield.(ensobs, :m_lh_match) .* sqrt.(8 * getfield.(ensobs,:t0))
m_D_star .= getfield.(ensobs, :m_lh_vec_match) .* sqrt.(8 * getfield.(ensobs,:t0))
#f_D .= getfield.(ensobs, :f_lh_match) .* sqrt.(8 * getfield.(ensobs,:t0))
#f_D_star .= getfield.(ensobs, :f_lh_vec_match) .* sqrt.(8 * getfield.(ensobs,:t0))
f_D .= getfield.(ensobs, :f_lh_match) .* sqrt.(8 * getfield.(ensobs,:t0))
f_D_star .= getfield.(ensobs, :f_lh_vec_match) .* sqrt.(8 * getfield.(ensobs,:t0))
uwerr.(m_D)
uwerr.(m_D_star)
#uwerr.(f_D)
#uwerr.(f_D_star)
uwerr.(f_D)
uwerr.(f_D_star)
catch
error("Check ensobs[i].mu_list to see whether all the ensembles have a light-heavy sector.")
end
......@@ -378,18 +386,18 @@ end
# ss
m_etaprime = Vector{uwreal}(undef, length(ensembles))
m_phi =Vector{uwreal}(undef, length(ensembles))
#f_etaprime= Vector{uwreal}(undef, length(ensembles))
#f_phi =Vector{uwreal}(undef, length(ensembles))
f_etaprime= Vector{uwreal}(undef, length(ensembles))
f_phi =Vector{uwreal}(undef, length(ensembles))
if sector["ss"]
try
m_etaprime .= getfield.(ensobs, :m_ss) .* sqrt.(8 * getfield.(ensobs,:t0))
m_phi .= getfield.(ensobs, :m_ss_vec) .* sqrt.(8 * getfield.(ensobs,:t0))
#f_etaprime .= getfield.(ensobs, :f_ss) .* sqrt.(8 * getfield.(ensobs,:t0))
#f_phi .= getfield.(ensobs, :f_ss_vec) .* sqrt.(8 * getfield.(ensobs,:t0))
f_etaprime .= getfield.(ensobs, :f_ss) .* sqrt.(8 * getfield.(ensobs,:t0))
f_phi .= getfield.(ensobs, :f_ss_vec) .* sqrt.(8 * getfield.(ensobs,:t0))
uwerr.(m_etaprime)
uwerr.(m_phi)
#uwerr.(f_etaprime)
#uwerr.(f_phi)
uwerr.(f_etaprime)
uwerr.(f_phi)
catch
error("Check ensobs[i].mu_list to see whether all the ensembles have a strange-strange sector.")
end
......@@ -397,18 +405,18 @@ end
# sh
m_Ds = Vector{uwreal}(undef, length(ensembles))
m_Ds_star =Vector{uwreal}(undef, length(ensembles))
#f_Ds= Vector{uwreal}(undef, length(ensembles))
#f_Ds_star =Vector{uwreal}(undef, length(ensembles))
f_Ds= Vector{uwreal}(undef, length(ensembles))
f_Ds_star =Vector{uwreal}(undef, length(ensembles))
if sector["sh"]
try
m_Ds .= getfield.(ensobs, :m_sh_match) .* sqrt.(8 * getfield.(ensobs,:t0))
m_Ds_star .= getfield.(ensobs, :m_sh_vec_match) .* sqrt.(8 * getfield.(ensobs,:t0))
#f_Ds .= getfield.(ensobs, :f_sh_match) .* sqrt.(8 * getfield.(ensobs,:t0))
#f_Ds_star .= getfield.(ensobs, :f_sh_vec_match) .* sqrt.(8 * getfield.(ensobs,:t0))
f_Ds .= getfield.(ensobs, :f_sh_match) .* sqrt.(8 * getfield.(ensobs,:t0))
f_Ds_star .= getfield.(ensobs, :f_sh_vec_match) .* sqrt.(8 * getfield.(ensobs,:t0))
uwerr.(m_Ds)
uwerr.(m_Ds_star)
#uwerr.(f_Ds)
#uwerr.(f_Ds_star)
uwerr.(f_Ds)
uwerr.(f_Ds_star)
catch
error("Check ensobs[i].mu_list to see whether all the ensembles have a strange-heavy sector.")
end
......@@ -416,18 +424,18 @@ end
# hh
m_etac = Vector{uwreal}(undef, length(ensembles))
m_jpsi =Vector{uwreal}(undef, length(ensembles))
#f_etac= Vector{uwreal}(undef, length(ensembles))
#f_jpsi =Vector{uwreal}(undef, length(ensembles))
f_etac= Vector{uwreal}(undef, length(ensembles))
f_jpsi =Vector{uwreal}(undef, length(ensembles))
if sector["hh"]
try
m_etac .= getfield.(ensobs, :m_hh_match) .* sqrt.(8 * getfield.(ensobs,:t0))
m_jpsi .= getfield.(ensobs, :m_hh_vec_match) .* sqrt.(8 * getfield.(ensobs,:t0))
#f_etac .= getfield.(ensobs, :f_hh_match) .* sqrt.(8 * getfield.(ensobs,:t0))
#f_jpsi .= getfield.(ensobs, :f_hh_vec_match) .* sqrt.(8 * getfield.(ensobs,:t0))
f_etac .= getfield.(ensobs, :f_hh_match) .* sqrt.(8 * getfield.(ensobs,:t0))
f_jpsi .= getfield.(ensobs, :f_hh_vec_match) .* sqrt.(8 * getfield.(ensobs,:t0))
uwerr.(m_etac)
uwerr.(m_jpsi)
#uwerr.(f_etac)
#uwerr.(f_jpsi)
uwerr.(f_etac)
uwerr.(f_jpsi)
catch
error("Check ensobs[i].mu_list to see whether all the ensembles have a heavy-heavy sector.")
end
......@@ -437,7 +445,7 @@ end
#===== CONTINUUM AND CHIRAL EXTRAPOLATION =======#
println("\n Continuum and Chiral Extrapolation \n")
#=
obs_db =Dict(
"muc"=> [muc],
"ll" => [m_pi, m_rho, f_pi, f_rho],
......@@ -447,7 +455,7 @@ obs_db =Dict(
"sh" => [m_Ds, m_Ds_star, f_Ds, f_Ds_star],
"hh" => [m_etac, m_jpsi, f_etac, f_jpsi]
)
=#
#=
obs_db =Dict(
"muc"=> [muc],
"ll" => [m_pi, m_rho],
......@@ -457,7 +465,7 @@ obs_db =Dict(
"sh" => [m_Ds, m_Ds_star],
"hh" => [m_etac, m_jpsi]
)
=#
#preparing observables of interest (data, name, ylabel)
obs = Vector{Vector{uwreal}}(undef,0)
......@@ -575,130 +583,3 @@ end
## ####################################
## NOT CONSIDER
data = read_mesons(path_dat, "G5", "G5")
corr_pp = corr_obs.(data, L=64)
##
hh_pp = corr_pp[6].obs
##
pp11 = hh_pp[97:end-7]
pp12 = hh_pp[98:end-6]
pp13 = hh_pp[99:end-5]
pp14 = hh_pp[100:end-4]
pp22 = hh_pp[99:end-5]
pp23 = hh_pp[100:end-4]
pp24 = hh_pp[101:end-3]
pp33 = hh_pp[101:end-3]
pp34 = hh_pp[102: end-2]
pp44 = hh_pp[103:end-1]
pp_sing = get_matrix([pp11, pp22, pp33, pp44], [pp12, pp13, pp14, pp23, pp24, pp34])
##
tau=4
pp11 = hh_pp[97:end-1-4*tau]
pp12 = hh_pp[97+tau:end-1-3*tau]
pp13 = hh_pp[97+2*tau:end-1-2*tau]
pp14 = hh_pp[100:end-4]
pp22 = hh_pp[97+2*tau:end-1-2*tau]
pp23 = hh_pp[97+3*tau:end-1-tau]
pp24 = hh_pp[101:end-1]
pp33 = hh_pp[97+4*tau:end-1]
pp34 = hh_pp[102: end-2]
pp44 = hh_pp[103:end-1]
pp_sing = get_matrix([pp11, pp22], [pp12])
##
##test lukas
evecs = uweigvecs(pp_sing[3], pp_sing[2])
proj_pp = Vector{uwreal}(undef, length(pp_sing))
for i = 1:length(pp_sing)
svec = evecs[:,1]
proj_pp[i] = uwdot(svec, uwdot(pp_sing[i], svec))[1]
uwerr(proj_pp[i])
end
uwerr.(pp11)
fig = figure()
errorbar(collect(1:length(proj_pp)), value.(proj_pp), err.(proj_pp), fmt="*",mec="blue", ecolor="blue")
errorbar(collect(1:length(pp11)), value.(pp11), err.(pp11), fmt="+",mec="red", ecolor="red")
yscale("log")
#ylim(-0.00000000001,0.0000000001)
display(fig )
##lukas eff mass
@.model(x,p) = p[1]*exp(-x*p[2]) + p[3]*exp(-x*p[4])
ydata = proj_pp[6:end-1]
uwerr.(ydata)
xdata = collect(1:length(ydata))
par = fit_routine(model, xdata, ydata,4)
##
mlmass = meff(proj_pp, [120-96, 155-96])
uwerr(lmass)
##
evals = getall_eigvals(pp_sing, 1)
evecs = getall_eigvecs(pp_sing, 1)
##
en = energies(evals)
fig = figure()
errorbar( collect(1:length(en[1])),value.(en[1]), yerr=err.(en[1]), fmt="+")
ylim(0.71,0.76)
display(fig)
## fitting the energy, very good
@.model(x,p) = p[1] + p[2]*exp(-x*(p[3]))
en1 = en[1][2:end-4]
xdata = collect(1:length(en1))
wpmm = Dict{Int64, Vector{Float64}}()
wpmm[303] = [-1.0, 1.5,-1.0, -1.0]
par = fit_routine(model, xdata, en1, 3, wpm=wpmm)
## try to extract matrix elements
function norm_evec(evec::Matrix{uwreal}, ctnot::Matrix{uwreal})
n=size(evec,1)
res_evec=Matrix{uwreal}(undef, n, n)
for i =1:n
aux_norm = (uwdot(evec[:,i], uwdot(ctnot, evec[:,i])).^2).^0.25
res_evec[i,:] = 1 ./aux_norm .* evec[:,i]
end
return res_evec
end
#1)normalise evecs and extract Ptilda for all times
norm = Array{Matrix{uwreal}}(undef, length(evecs))
for i=1:length(evecs)
norm[i] = norm_evec(evecs[i], pp_sing[4])
end
#2) construct matrix M(t,t_0)=c(t_0)Ptilda for all times
mtt0 = Array{Matrix{uwreal}}(undef, length(evecs))
for i=1:length(evecs)
mtt0[i] = uwdot(pp_sing[4], juobs.transpose(norm[i]))
end
## extract decays
function pseudo_mat_elem(evec::Array{Matrix{uwreal}}, mass::uwreal, mu::Array{Float64}, t0::Int64)
Rn = [exp(mass * (t0)/2) * evec[t][1] for t =1:length(evec)]
aux = sqrt(2)*sum(mu) / mass^1.5
return uwdot(aux,Rn)
end
test_dec = pseudo_mat_elem(mtt0, par[1], [0.14,0.14],4)
## Standard
mass = meff(hh_pp, [123, 155] )
f = dec_const_pcvc(hh_pp, [123,155], mass, [0.14,0.14],96)
## testing generalised
todiag = uwdot(juobs.invert(pp_sing[4]),pp_sing[10])
evals, evecs = uweigen(todiag)
##
#test = uwdot(juobs.invert(evecs), uwdot(uwdot(juobs.invert(pp_sing[4]),pp_sing[10]), evecs))
test = uwdot(evecs, uwdot(evals,juobs.invert(evecs)))
##
evals2, evecs2 = uweigen(pp_sing[10], pp_sing[4])
#test2 = uwdot(juobs.invert(evecs2), uwdot(uwdot(juobs.invert(pp_sing[4]),pp_sing[10]), evecs2))
test2 = uwdot(pp_sing[4],uwdot(evecs2, uwdot(evals2,juobs.invert(evecs2))))
\ No newline at end of file
tools.jl → src/tools.jl
View file @ 43d0b8ef
......@@ -164,6 +164,39 @@ function fit_mass(en_i::Array{uwreal}; t_in::Int64=3, wpm::Union{Nothing, Dict{I
println("\nPar[1] is ", par[1] )
return par[1]
end
function gevp_dec(mat_obs::Vector{MatInfo}, mass::Vector{uwreal}; 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_dec0 = Vector{uwreal}(undef, length(mat_obs))
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], 1) #ground and excited state energy
plat = select_plateau(mat_obs[i].ensinfo, mu_list)[end] .- y0s[i]
plat_dec0[i] = -dec(elem, plat, mass[i], mu_list[i])
end
return plat_dec0
end
function getall_eig(a::Vector{Matrix{uwreal}}, t0; iter=30)
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})
aux = plat_av(mat_elem, plat)
return sqrt(2/ (mass)^3) * sum(mu)* aux
end
function param_av(fit_masses::Vector{uwreal})
w = 1 ./ err.(fit_masses).^2
av = sum(w .* fit_masses) / sum(w)
......
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