Skip to content

J
juobs
Commit a4c6307a authored by Alessandro 's avatar Alessandro
Browse files
Options

gevp_tools folder removed

parent 3e20520f
Hide whitespace changes
Inline Side-by-side
Showing with 0 additions and 883 deletions
gevp tools/analysis.jl deleted 100644 → 0
View file @ 3e20520f
using juobs, ADerrors, DelimitedFiles, PyPlot, LaTeXStrings, Revise, LsqFit, DataFrames, CSV
include("/Users/ale/juobs/gevp tools/tools.jl")
include("/Users/ale/juobs/constants/juobs_const.jl")
##
const path_plot = "/Users/ale/Google Drive/phd/analysis/gevp analysis/plot"
const path_result = "/Users/ale/Google Drive/phd/analysis/gevp analysis"
const path_data = "/Users/ale/Desktop/data"
##
const ensembles = ["H400", "N200", "N203", "N300", "J303"]
const deg = [true, false, false, true, false]
const L = [32, 48, 48, 48, 64]
const beta = [3.46 , 3.55, 3.55, 3.70, 3.70]
const R = ["H400r001", ["N200r000", "N200r001"], ["N203r000", "N203r001"], "N300r002", "J303r003"]
mu_pp = Vector{Vector{Vector{Float64}}}(undef, length(ensembles))
mu_aa = Vector{Vector{Vector{Float64}}}(undef, length(ensembles))
m_lh = Vector{Vector{uwreal}}(undef, length(ensembles))
#### choose the ensembles and read data
iens =5
pp = read_dat(ensembles[iens], "G5","G5")
a1a1 = read_dat(ensembles[iens], "G1G5","G1G5")
a2a2 = read_dat(ensembles[iens], "G2G5","G2G5")
a3a3 = read_dat(ensembles[iens], "G3G5","G3G5")
pa1 = read_dat(ensembles[iens], "G5","G1G5")
pa2 = read_dat(ensembles[iens], "G5","G2G5")
pa3 = read_dat(ensembles[iens], "G5","G3G5")
pp_obs = corr_obs.(pp, L=L[iens])
a1a1_obs = corr_obs.(a1a1, L=L[iens])
a2a2_obs = corr_obs.(a2a2, L=L[iens])
a3a3_obs = corr_obs.(a3a3, L=L[iens])
pa1_obs = corr_obs.(pa1, L=L[iens])
pa2_obs = corr_obs.(pa2, L=L[iens])
pa3_obs = corr_obs.(pa3, L=L[iens])
mu_pp[iens] = getfield.(pp_obs, :mu)
##############################
## GEVP ANALYSIS #############
##############################
###################
## HH SECTOR
####################
#Initialize result container
##initialize partial results
mul, mus, muh = get_mu(mu_pp[iens], deg[iens])
m_lh = Vector{uwreal}(undef, length(muh))
m_sh = Vector{uwreal}(undef, length(muh))
m_lh_star = Vector{uwreal}(undef, length(muh))
m_sh_star = Vector{uwreal}(undef, length(muh))
m_hh = Vector{uwreal}(undef, length(muh))
m_hh_star = Vector{uwreal}(undef, length(muh))
f_lh = Vector{uwreal}(undef, length(muh))
f_sh = Vector{uwreal}(undef, length(muh))
f_lh_star = Vector{uwreal}(undef, length(muh))
f_sh_star = Vector{uwreal}(undef, length(muh))
f_hh = Vector{uwreal}(undef, length(muh))
f_hh_star = Vector{uwreal}(undef, length(muh))
##
hh11 = get_hh(mu_pp[iens], [pp_obs], deg[iens])
hh12 = get_hh(mu_pp[iens], [pa1_obs, pa2_obs, pa3_obs], deg[iens])
hh22 = get_hh(mu_pp[iens], [a1a1_obs, a2a2_obs, a3a3_obs], deg[iens])
gevp_hhmat = comp_mat(hh11, hh12, hh22)
@time hh_en, hh_evec = comp_energy(hh11,hh12,hh22, true)
## HH fitting procedure
mu_hh = 3
mu_hs = 3
mu_hl = 3
@. model(x,p)=p[1] + p[2] * exp(-(p[3]-p[1])*x)
ydata1 = hh_en[mu_hh].en_val[1][15:end-15]
ydata2 = hh_en[mu_hh].en_val[2][13:end-10]
xdata1 = collect(1:length(ydata1))
xdata2 = collect(1:length(ydata2))
wpm = Dict{Int64, Vector{Float64}}()
wpm[303] = [-1.0, 2.0,-1.0, -1.0]
##
hhpar1 = fit_routine(model, ydata1, wpm=wpm)
##
hhpar2 = fit_routine(model, ydata2, wpm=wpm)
##
fig = figure()
plt.errorbar(xdata1, value.(ydata1), yerr = err.(ydata1), ms=2, fmt="^", mfc="none", mec="cornflowerblue",ecolor="cornflowerblue", capsize=2, label=L"$am_{\eta_c}$") #am_{D_s}
plt.plot(xdata1, model(xdata1, value.(hhpar1)), lw=0.5, c="midnightblue" )
plt.errorbar(xdata2, value.(ydata2), yerr = err.(ydata2), ms=2, fmt="^", mfc="none", mec="tomato", ecolor="tomato" ,capsize=2, label=L"$am_{J/\psi}$") #am_{D_s^*}
plt.plot(xdata2, model(xdata2, value.(hhpar2)), lw=0.5, c="darkred" )
plt.xlim(1,50)
#plt.ylim(0.70,0.80)
plt.xlabel(L"$t$")
plt.ylabel(L"$E(t, t_0=5)$")
plt.title(string(ensembles[iens]))#," ", L"$\mu_c = $", fen[mu_sec].mu) )
plt.axvline(x=5, lw = 1, c ="green", ls="dashed", label=L"$t_0$")
plt.legend()
name = string(ensembles[iens], "_hh_gevp.pdf")
#plt.savefig(joinpath(path_plot, name), dpi=500)
display(fig)
## temporary container for hh masses
hhps_mass = hhpar1[1]
hhvec_mass = hhpar2[1]
m_hh[mu_hh] = hhpar1[1]
m_hh_star[mu_hh] = hhpar2[1]
## pseudo decay hh
ps_pneff = pseudo_mat_elem(hh_evec[mu_hh], hh11[1][mu_hh].obs, hhps_mass)
uwerr.(ps_pneff)
@. dec_mod(x,p) = p[1] + p[2] *x
ydat = ps_pneff[30:80]
ps_tofit = fit_routine(dec_mod, ydat, 2, wpm=wpm)
figu = figure()
plt.errorbar(collect(1:length(ydat)), value.(ydat), yerr = err.(ydat), fmt="d")
#plt.ylim(0.26,0.27)
plt.plot(collect(1:length(ydat)), (dec_mod(collect(1:length(ydat)),value.(ps_tofit) )))
display(figu)
hhps_dec = extract_dec_const(ps_tofit[1], hhps_mass, hh_en[mu_hh].mu)
uwerr(hhps_dec)
hhps_dec
f_hh[mu_hh] = hhps_dec
## vec decay hh
ak_corr = (hh22[1][mu_hh].obs .+ hh22[2][mu_hh].obs .+ hh22[3][mu_hh].obs)
aux_vec = vec_mat_elem(hh_evec[mu_hh], ak_corr, gevp_hhmat[mu_hh].mat, hhvec_mass)
uwerr.(aux_vec)
ydat = aux_vec[30:80]
vec_tofit = fit_routine(dec_mod, ydat, 2)
figu = figure()
plt.errorbar(collect(1:length(ydat)), value.(ydat), yerr = err.(ydat), fmt="d")
#plt.ylim(0.26,0.27)
plt.plot(collect(1:length(ydat)), (dec_mod(collect(1:length(ydat)),value.(vec_tofit) )))
display(figu)
hhvec_dec = extract_dec_const(vec_tofit[1], hhvec_mass, hh_en[mu_hh].mu)
uwerr(hhvec_dec)
hhvec_dec
f_hh_star[mu_hh] = hhvec_dec
#################
## HS SECTOR ####
#################
hs11 = get_sh(mu_pp[iens], [pp_obs], deg[iens])
hs12 = get_sh(mu_pp[iens], [pa1_obs, pa2_obs, pa3_obs], deg[iens])
hs22 = get_sh(mu_pp[iens], [a1a1_obs, a2a2_obs, a3a3_obs], deg[iens])
gevp_hsmat = comp_mat(hs11, hs12, hs22)
hs_en, hs_evec = comp_energy(hs11, hs12, hs22, true)
## HS fitting procedure
@. model(x,p)=p[1] + p[2] * exp(-(p[3]-p[1])*x)
ydata1 = hs_en[mu_hs].en_val[1][11:end-10]
ydata2 = hs_en[mu_hs].en_val[2][11:end-10]
xdata1 = collect(1:length(ydata1))
xdata2 = collect(1:length(ydata2))
wpm = Dict{Int64, Vector{Float64}}()
wpm[303] = [-1.0, 2.0,-1.0, -1.0]
println(hs_en[mu_hs].mu)
##
hspar1 = fit_routine(model, ydata1, wpm=wpm)
##
hspar2 = fit_routine(model, ydata2, wpm=wpm)
##
fig = figure()
plt.errorbar(xdata1, value.(ydata1), yerr = err.(ydata1), ms=2, fmt="^", mfc="none", mec="cornflowerblue",ecolor="cornflowerblue", capsize=2, label=L"$am_{\eta_c}$") #am_{D_s}
plt.plot(xdata1, model(xdata1, value.(hspar1)), lw=0.5, c="midnightblue" )
plt.errorbar(xdata2, value.(ydata2), yerr = err.(ydata2), ms=2, fmt="^", mfc="none", mec="tomato", ecolor="tomato" ,capsize=2, label=L"$am_{J/\psi}$") #am_{D_s^*}
plt.plot(xdata2, model(xdata2, value.(hspar2)), lw=0.5, c="darkred" )
plt.xlim(1,40)
#plt.ylim(0.40,0.7)
plt.xlabel(L"$t$")
plt.ylabel(L"$E(t, t_0=5)$")
plt.title(string(ensembles[iens]))#," ", L"$\mu_c = $", fen[mu_sec].mu) )
plt.axvline(x=5, lw = 1, c ="green", ls="dashed", label=L"$t_0$")
plt.legend()
name = string(ensembles[iens], "_hl_gevp.pdf")
#plt.savefig(joinpath(path_plot, name), dpi=500)
display(fig)
##temporary container for hs masses
hsps_mass = hspar1[1]
hsvec_mass = hspar2[1]
m_sh[mu_hs] = hspar1[1]
m_sh_star[mu_hs] = hspar2[1]
## pseudo decay hs
ps_pneff = pseudo_mat_elem(hs_evec[mu_hs], hs11[1][mu_hs].obs, hsps_mass)
uwerr.(ps_pneff)
@. dec_mod(x,p) = p[1] + p[2] *x
ydat = ps_pneff[30:60]
ps_tofit = fit_routine(dec_mod, ydat, 2)
figu = figure()
plt.errorbar(collect(1:length(ydat)), value.(ydat), yerr = err.(ydat), fmt="d")
#plt.ylim(0.26,0.27)
plt.plot(collect(1:length(ydat)), (dec_mod(collect(1:length(ydat)),value.(ps_tofit) )))
display(figu)
hsps_dec = extract_dec_const(ps_tofit[1], hsps_mass, hs_en[mu_hs].mu)
uwerr(hsps_dec)
hsps_dec
f_sh[mu_hs] = hsps_dec
## vec decay hs
ak_corr = (hs22[1][mu_hs].obs .+ hs22[2][mu_hs].obs .+ hs22[3][mu_hs].obs)
aux_vec = vec_mat_elem(hs_evec[mu_hs], ak_corr, gevp_hsmat[mu_hs].mat, hsvec_mass)
uwerr.(aux_vec)
ydat = aux_vec[25:50]
vec_tofit = fit_routine(dec_mod, ydat, 2)
figu = figure()
plt.errorbar(collect(1:length(ydat)), value.(ydat), yerr = err.(ydat), fmt="d")
#plt.ylim(0.26,0.27)
plt.plot(collect(1:length(ydat)), (dec_mod(collect(1:length(ydat)),value.(vec_tofit) )))
display(figu)
hsvec_dec = extract_dec_const(vec_tofit[1], hsvec_mass, hs_en[mu_hs].mu)
uwerr(hsvec_dec)
hsvec_dec
f_sh_star[mu_hs] = hsvec_dec
#################
## HL SECTOR ####
#################
hl11 = get_lh(mu_pp[iens], [pp_obs], deg[iens])
hl12 = get_lh(mu_pp[iens], [pa1_obs, pa2_obs, pa3_obs], deg[iens])
hl22 = get_lh(mu_pp[iens], [a1a1_obs, a2a2_obs, a3a3_obs], deg[iens])
gevp_hlmat = comp_mat(hl11, hl12, hl22)
hl_en, hl_evec = comp_energy(hl11, hl12, hl22, true)
## HL fitting procedure
@. model(x,p)=p[1] + p[2] * exp(-(p[3]-p[1])*x)
ydata1 = hl_en[mu_hl].en_val[1][10:end-10]
ydata2 = hl_en[mu_hl].en_val[2][10:end-10]
xdata1 = collect(1:length(ydata1))
xdata2 = collect(1:length(ydata2))
wpm = Dict{Int64, Vector{Float64}}()
wpm[303] = [-1.0, 2.0,-1.0, -1.0]
println(hl_en[mu_hl].mu)
##
hlpar1 = fit_routine(model, ydata1, wpm=wpm)
##
hlpar2 = fit_routine(model, ydata2, wpm=wpm)
##
fig = figure()
plt.errorbar(xdata1, value.(ydata1), yerr = err.(ydata1), ms=2, fmt="^", mfc="none", mec="cornflowerblue",ecolor="cornflowerblue", capsize=2, label=L"$am_{\eta_c}$") #am_{D_s}
plt.plot(xdata1, model(xdata1, value.(hlpar1)), lw=0.5, c="midnightblue" )
plt.errorbar(xdata2, value.(ydata2), yerr = err.(ydata2), ms=2, fmt="^", mfc="none", mec="tomato", ecolor="tomato" ,capsize=2, label=L"$am_{J/\psi}$") #am_{D_s^*}
plt.plot(xdata2, model(xdata2, value.(hlpar2)), lw=0.5, c="darkred" )
plt.xlim(1,40)
plt.ylim(0.40,0.7)
plt.xlabel(L"$t$")
plt.ylabel(L"$E(t, t_0=5)$")
plt.title(string(ensembles[iens]))#," ", L"$\mu_c = $", fen[mu_sec].mu) )
plt.axvline(x=5, lw = 1, c ="green", ls="dashed", label=L"$t_0$")
plt.legend()
name = string(ensembles[iens], "_hl_gevp.pdf")
#plt.savefig(joinpath(path_plot, name), dpi=500)
display(fig)
##temporary container for hs masses
hlps_mass = hlpar1[1]
hlvec_mass = hlpar2[1]
m_lh[mu_hl] = hlpar1[1]
m_lh_star[mu_hl] = hlpar2[1]
## pseudo decay hl
ps_pneff = pseudo_mat_elem(hl_evec[mu_hl], hl11[1][mu_hl].obs, hlps_mass)
uwerr.(ps_pneff)
@. dec_mod(x,p) = p[1] + p[2] *x
ydat = ps_pneff[20:50]
ps_tofit = fit_routine(dec_mod, ydat, 2)
figu = figure()
plt.errorbar(collect(1:length(ydat)), value.(ydat), yerr = err.(ydat), fmt="d")
#plt.ylim(0.26,0.27)
plt.plot(collect(1:length(ydat)), (dec_mod(collect(1:length(ydat)),value.(ps_tofit) )))
display(figu)
hlps_dec = extract_dec_const(ps_tofit[1], hlps_mass, hl_en[mu_hl].mu)
uwerr(hlps_dec)
hlps_dec
##store hl dec
f_lh[mu_hl] = hlps_dec
## vec decay hl
ak_corr = (hl22[1][mu_hl].obs .+ hl22[2][mu_hl].obs .+ hl22[3][mu_hl].obs)
aux_vec = vec_mat_elem(hl_evec[mu_hl], ak_corr, gevp_hlmat[mu_hl].mat, hlvec_mass)
uwerr.(aux_vec)
ydat = aux_vec[20:50]
vec_tofit = fit_routine(dec_mod, ydat, 2)
figu = figure()
plt.errorbar(collect(1:length(ydat)), value.(ydat), yerr = err.(ydat), fmt="d")
#plt.ylim(0.26,0.27)
plt.plot(collect(1:length(ydat)), (dec_mod(collect(1:length(ydat)),value.(vec_tofit) )))
display(figu)
hlvec_dec = extract_dec_const(vec_tofit[1], hlvec_mass, hl_en[mu_hl].mu)
uwerr(hlvec_dec)
hlvec_dec
##store hl_star dec
f_lh_star[mu_hl] = hlvec_dec
#############
## MATCHING
#############
mul, mus, muh = get_mu(mu_pp[iens], deg[iens])
##
target = a(beta[iens]) * (2*M[1] + 6*M[2] + M[3] + 3*M[4]) / (12*hc)
if !deg[iens]
muh_target = match_muc(muh, m_lh, m_sh, m_lh_star, m_sh_star, target)
else
muh_target = match_muc(muh, m_lh, m_lh_star, target)
end
uwerr(muh_target)
muh_target
##interpolate m_lh, m_lh_star, m_sh, m_sh_star, m_hh, m_hh_star
#m_lh
par, chi2exp = lin_fit(muh, m_lh)
m_lh_match = y_lin_fit(par, muh_target)
uwerr(m_lh_match)
#m_lh_star
par, chi2exp = lin_fit(muh, m_lh_star)
m_lh_star_match = y_lin_fit(par, muh_target)
uwerr(m_lh_star_match)
#m_sh
par, chi2exp = lin_fit(muh, m_sh)
m_sh_match = y_lin_fit(par, muh_target)
uwerr(m_sh_match)
#m_sh_star
par, chi2exp = lin_fit(muh, m_sh_star)
m_sh_star_match = y_lin_fit(par, muh_target)
uwerr(m_sh_star_match)
#m_hh
par, chi2exp = lin_fit(muh, m_hh)
m_hh_match = y_lin_fit(par, muh_target)
uwerr(m_hh_match)
#m_hh_star
par, chi2exp = lin_fit(muh, m_hh_star)
m_hh_star_match = y_lin_fit(par, muh_target)
uwerr(m_hh_star_match)
#f_lh
par, chi2exp = lin_fit(muh, f_lh)
f_lh_match = y_lin_fit(par, muh_target)
uwerr(f_lh_match)
#f_lh_star
par, chi2exp = lin_fit(muh, f_lh_star)
f_lh_star_match = y_lin_fit(par, muh_target)
uwerr(f_lh_star_match)
#f_sh
par, chi2exp = lin_fit(muh, f_sh)
f_sh_match = y_lin_fit(par, muh_target)
uwerr(f_sh_match)
#f_sh_star
par, chi2exp = lin_fit(muh, f_sh_star)
f_sh_star_match = y_lin_fit(par, muh_target)
uwerr(f_sh_star_match)
#f_hh
par, chi2exp = lin_fit(muh, f_hh)
f_hh_match = y_lin_fit(par, muh_target)
uwerr(f_hh_match)
#f_hh_star
par, chi2exp = lin_fit(muh, f_hh_star)
f_hh_star_match = y_lin_fit(par, muh_target)
uwerr(f_hh_star_match)
###################
## PLOTS #########
##################
## m_lh, m_lh_star
figure()
title(ensembles[iens])
xlabel(L"$a\mu_c$")
ylabel(L"$aM$")
errorbar(muh, value.(m_lh), err.(m_lh), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
errorbar(value(muh_target), value(m_lh_match), err(m_lh_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
errorbar(muh, value.(m_lh_star), err.(m_lh_star), fmt="*",ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
errorbar(value(muh_target), value(m_lh_star_match), err(m_lh_star_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="navy", lw=1)
legend([L"$m_D$", L"$\langle m_D \rangle$", L"$m_{D^*}$", L"$ \langle m_{D^*}\rangle $"])
display(gcf())
t = string(ensembles[iens], "_mD.pdf")
savefig(joinpath(path_plot,t))
close()
## m_sh, m_sh_star
figure()
title(ensembles[iens])
xlabel(L"$a\mu_c$")
ylabel(L"$aM$")
errorbar(muh, value.(m_sh), err.(m_sh), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
errorbar(value(muh_target), value(m_sh_match), err(m_sh_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
errorbar(muh, value.(m_sh_star), err.(m_sh_star), fmt="*",ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
errorbar(value(muh_target), value(m_sh_star_match), err(m_sh_star_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="navy", lw=1)
legend([L"$m_{D_s}$", L"$\langle m_{D_s} \rangle$", L"$m_{D^*_s}$", L"$ \langle m_{D^*_s}\rangle $"])
display(gcf())
t = string(ensembles[iens], "_mDs.pdf")
savefig(joinpath(path_plot,t))
close()
## m_hh, m_hh_star
figure()
title(ensembles[iens])
xlabel(L"$a\mu_c$")
ylabel(L"$aM$")
errorbar(muh, value.(m_hh), err.(m_hh), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
errorbar(value(muh_target), value(m_hh_match), err(m_hh_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
errorbar(muh, value.(m_hh_star), err.(m_hh_star), fmt="*",ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
errorbar(value(muh_target), value(m_hh_star_match), err(m_hh_star_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="navy", lw=1)
legend([L"$m_{\eta_c}$", L"$\langle m_{\eta_c} \rangle$", L"$m_{J/\psi}$", L"$ \langle m_{J/\psi}\rangle $"])
display(gcf())
t = string(ensembles[iens], "_m_eta.pdf")
savefig(joinpath(path_plot,t))
close()
## f_lh, f_sh
figure()
title(ensembles[iens])
xlabel(L"$a\mu_c$")
ylabel(L"$af$")
errorbar(muh, value.(f_lh), err.(f_lh), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
errorbar(value(muh_target), value(f_lh_match), err(f_lh_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
errorbar(muh, value.(f_sh), err.(f_sh), fmt="*",ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
errorbar(value(muh_target), value(f_sh_match), err(f_sh_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="navy", lw=1)
legend([L"$f_{D}$", L"$\langle f_{D} \rangle$", L"$f_{D_s}$", L"$ \langle f_{D_s}\rangle $"])
display(gcf())
t = string(ensembles[iens], "_fD.pdf")
savefig(joinpath(path_plot,t))
close()
## f_lh_star, f_sh_star
figure()
title(ensembles[iens])
xlabel(L"$a\mu_c$")
ylabel(L"$af$")
errorbar(muh, value.(f_lh_star), err.(f_lh_star), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
errorbar(value(muh_target), value(f_lh_star_match), err(f_lh_star_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
errorbar(muh, value.(f_sh_star), err.(f_sh_star), fmt="*",ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
errorbar(value(muh_target), value(f_sh_star_match), err(f_sh_star_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="navy", lw=1)
legend([L"$f_{D^*}$", L"$\langle f_{D^*} \rangle$", L"$f_{D_s^*}$", L"$ \langle f_{D_s^*}\rangle $"])
display(gcf())
t = string(ensembles[iens], "_fD*.pdf")
savefig(joinpath(path_plot,t))
close()
## f_hh, f_hh_star
figure()
title(ensembles[iens])
xlabel(L"$a\mu_c$")
ylabel(L"$af$")
errorbar(muh, value.(f_hh), err.(f_hh), fmt="*", ms=5, mfc="none", capsize=2, c="tomato", lw=1)
errorbar(value(muh_target), value(f_hh_match), err(f_hh_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="darkred", lw=1)
errorbar(muh, value.(f_hh_star), err.(f_hh_star), fmt="*",ms=5, mfc="none", capsize=2, c="cornflowerblue", lw=1)
errorbar(value(muh_target), value(f_hh_star_match), err(f_hh_star_match), err(muh_target), fmt="*", ms=5, mfc="none", capsize=2, c="navy", lw=1)
legend([L"$f_{\eta_c}$", L"$\langle f_{\eta_c} \rangle$", L"$f_{J/\psi}$", L"$ \langle f_{J/\psi}\rangle $"])
display(gcf())
t = string(ensembles[iens], "_f_eta.pdf")
savefig(joinpath(path_plot,t))
close()
###############
## RESULTS#####
###############
#Quark mass
muc = zm_tm(beta[iens]) * muh_target * sqrt(8 * t0(beta[iens]))
uwerr(muc)
#Meson masses
m_D = m_lh_match * sqrt(8 * t0(beta[iens]))
uwerr(m_D)
m_Ds = m_sh_match * sqrt(8 * t0(beta[iens]))
uwerr(m_Ds)
m_D_star = m_lh_star_match * sqrt(8 * t0(beta[iens]))
uwerr(m_D_star)
m_Ds_star = m_sh_star_match * sqrt(8 * t0(beta[iens]))
uwerr(m_Ds_star)
m_eta = m_hh_match * sqrt(8 * t0(beta[iens]))
uwerr(m_eta)
m_j = m_hh_star_match * sqrt(8 * t0(beta[iens]))
uwerr(m_j)
#Decay constant
f_D = f_lh_match * sqrt(8 * t0(beta[iens]))
uwerr(f_D)
f_Ds = f_sh_match * sqrt(8 * t0(beta[iens]))
uwerr(f_Ds)
f_D_star = f_lh_star_match * sqrt(8 * t0(beta[iens]))
uwerr(f_D_star)
f_Ds_star = f_sh_star_match * sqrt(8 * t0(beta[iens]))
uwerr(f_Ds_star)
f_eta = f_hh_match * sqrt(8 * t0(beta[iens]))
uwerr(f_eta)
f_j = f_hh_star_match * sqrt(8 * t0(beta[iens]))
uwerr(f_j)
println(" ", ensembles[iens])
println(L"Z^{tm}_M \mu_c \sqrt{8t_0} = ", muc)
println( L"M_D \sqrt{8t_0} = ", m_D)
println( L"M_Ds \sqrt{8t_0} = ", m_Ds)
println( L"M_D* \sqrt{8t_0} = ", m_D_star)
println( L"M_Ds* \sqrt{8t_0} = ", m_Ds_star)
println( L"M_eta* \sqrt{8t_0} = ", m_eta)
println( L"M_jpsi* \sqrt{8t_0} = ", m_j)
println( L"f_D \sqrt{8t_0} = ", f_D)
println( L"f_Ds \sqrt{8t_0} = ", f_Ds)
println( L"f_D* \sqrt{8t_0} = ", f_D_star)
println( L"f_Ds* \sqrt{8t_0} = ", f_Ds_star)
println( L"f_eta \sqrt{8t_0} = ", f_eta)
println( L"f_jpsi \sqrt{8t_0} = ", f_j)
## DataFrames with partial mass data
mass_df = DataFrame(ens=String[], muc=Float64[], mDs= Float64[], mDserr= Float64[], mDst= Float64[], mDsterr= Float64[], meta= Float64[], metaerr= Float64[], mJ= Float64[], mJerr=Float64[] )
##
push!(mass_df , ["N203", 0.173375, 0.61040, 35e-5, 0.6567, 10e-4, 0.93127, 14e-5, 0.96969, 25e-5])
push!(mass_df , ["N203", 0.18225, 0.62854, 34e-5, 0.6732, 10e-4, 0.96287, 14e-5, 1.00036, 24e-5])
push!(mass_df , ["N203", 0.191625, 0.64604, 41e-5, 0.6895, 10e-4, 0.99403, 13e-5, 1.03073, 22e-5])
push!(mass_df , ["N200", 0.17285, 0.61362, 39e-5, 0.66079, 65e-5,0.92917 , 14e-5 ,0.96765 ,31e-5])
push!(mass_df , ["N200", 0.1819, 0.63145, 44e-5, 0.67644, 66e-5, 0.96076, 13e-5 , 0.99831 , 30e-5])
push!(mass_df , ["N200", 0.190995, 0.64913, 41e-5, 0.69333, 66e-5, 0.99178, 13e-5 ,1.02860 ,29e-5 ])
##
path_result
##
CSV.write(joinpath(path_result,"uncomplete_result.csv"), mass_df)
gevp tools/results.rtf deleted 100644 → 0
View file @ 3e20520f
{\rtf1\ansi\ansicpg1252\cocoartf2509
\cocoatextscaling0\cocoaplatform0{\fonttbl\f0\fswiss\fcharset0 Helvetica;}
{\colortbl;\red255\green255\blue255;}
{\*\expandedcolortbl;;}
\paperw11900\paperh16840\margl1440\margr1440\vieww10800\viewh8400\viewkind0
\pard\tx566\tx1133\tx1700\tx2267\tx2834\tx3401\tx3968\tx4535\tx5102\tx5669\tx6236\tx6803\pardirnatural\partightenfactor0
\f0\fs24 \cf0 $H400\
$Z^\{tm\}_M \\mu_c \\sqrt\{8t_0\} = $3.030627780344985 +/- 0.07902594921440798\
$M_D \\sqrt\{8t_0\} = $3.981135432856687 +/- 0.05781090729513162\
$M_Ds \\sqrt\{8t_0\} = $3.981135432856687 +/- 0.05781090729513162\
$M_D* \\sqrt\{8t_0\} = $4.278807884262479 +/- 0.05002067830312485\
$M_Ds* \\sqrt\{8t_0\} = $4.278807884262479 +/- 0.05002067830312485\
$M_eta* \\sqrt\{8t_0\} = $6.096944664722402 +/- 0.09777085253831523\
$M_jpsi* \\sqrt\{8t_0\} = $6.3445407150455315 +/- 0.09568980579456061\
$f_D \\sqrt\{8t_0\} = $0.516849796654617 +/- 0.0030631272287505134\
$f_Ds \\sqrt\{8t_0\} = $0.516849796654617 +/- 0.0030631272287505134\
$f_D* \\sqrt\{8t_0\} = $0.21312484071325435 +/- 0.003968470514564074\
$f_Ds* \\sqrt\{8t_0\} = $0.21312484071325435 +/- 0.003968470514564074\
$f_eta \\sqrt\{8t_0\} = $0.9689427427736942 +/- 0.015398690942284236\
$f_jpsi \\sqrt\{8t_0\} = $0.4391703592723876 +/- 0.009322107275879348\
$\
\pard\tx566\tx1133\tx1700\tx2267\tx2834\tx3401\tx3968\tx4535\tx5102\tx5669\tx6236\tx6803\pardirnatural\partightenfactor0
\cf0 $N200\
$Z^\{tm\}_M \\mu_c \\sqrt\{8t_0\} = $3.095508084190539 +/- 0.07424875770953374\
$M_D \\sqrt\{8t_0\} = $3.958264727738654 +/- 0.05515406625525418\
$M_Ds \\sqrt\{8t_0\} = $4.072538504215807 +/- 0.05467443552294891\
$M_D* \\sqrt\{8t_0\} = $4.229272703354319 +/- 0.05006705711946871\
$M_Ds* \\sqrt\{8t_0\} = $4.361851498402342 +/- 0.04975162489041616\
$M_eta* \\sqrt\{8t_0\} = $6.198773176967105 +/- 0.0967669789483958\
$M_jpsi* \\sqrt\{8t_0\} = $6.439877922484902 +/- 0.0942314903150228\
$f_D \\sqrt\{8t_0\} = $0.47414683350372805 +/- 0.002821170657997782\
$f_Ds \\sqrt\{8t_0\} = $0.5204052951135102 +/- 0.002351052892691432\
$f_D* \\sqrt\{8t_0\} = $0.18581645952897405 +/- 0.0032030018851700124\
$f_Ds* \\sqrt\{8t_0\} = $0.20920280395715848 +/- 0.0024790842850004667\
$f_eta \\sqrt\{8t_0\} = $0.9266388913538548 +/- 0.011566795910907063\
$f_jpsi \\sqrt\{8t_0\} = $0.42904806660632133 +/- 0.007884932722702343\
$\
$N203\
$Z^\{tm\}_M \\mu_c \\sqrt\{8t_0\} = $3.07608073397371 +/- 0.07386718423714976\
$M_D \\sqrt\{8t_0\} = $3.966409130806326 +/- 0.054959061347465483\
$M_Ds \\sqrt\{8t_0\} = $4.030944118480624 +/- 0.05473249525087825\
$M_D* \\sqrt\{8t_0\} = $4.254701932623177 +/- 0.04997144163893854\
$M_Ds* \\sqrt\{8t_0\} = $4.320267491882463 +/- 0.049563464976048445\
$M_eta* \\sqrt\{8t_0\} = $6.175405773876267 +/- 0.09615293679832428\
$M_jpsi* \\sqrt\{8t_0\} = $6.418715171251667 +/- 0.09346001316883264\
$f_D \\sqrt\{8t_0\} = $0.48810439044765425 +/- 0.0026045102710287034\
$f_Ds \\sqrt\{8t_0\} = $0.513802783972585 +/- 0.002279763193762235\
$f_D* \\sqrt\{8t_0\} = $0.19370973720440873 +/- 0.0028186032755263607\
$f_Ds* \\sqrt\{8t_0\} = $0.20485862913653624 +/- 0.0024571277939203565\
$f_eta \\sqrt\{8t_0\} = $0.9262934445841321 +/- 0.011380030137209169\
$f_jpsi \\sqrt\{8t_0\} = $0.4287220496155749 +/- 0.007707612363034976\
$\
$N300\
$Z^\{tm\}_M \\mu_c \\sqrt\{8t_0\} = $2.996157019535518 +/- 0.07903034846063876\
$M_D \\sqrt\{8t_0\} = $3.9418982012432258 +/- 0.0596919625451062\
$M_Ds \\sqrt\{8t_0\} = $3.9418982012432258 +/- 0.0596919625451062\
$M_D* \\sqrt\{8t_0\} = $4.2918918432116255 +/- 0.04985913651457226\
$M_Ds* \\sqrt\{8t_0\} = $4.2918918432116255 +/- 0.04985913651457226\
$M_eta* \\sqrt\{8t_0\} = $6.077559366401819 +/- 0.10543357205926039\
$M_jpsi* \\sqrt\{8t_0\} = $6.322315861178219 +/- 0.1024859602694168\
$f_D \\sqrt\{8t_0\} = $0.502826744878361 +/- 0.0026201800440417583\
$f_Ds \\sqrt\{8t_0\} = $0.502826744878361 +/- 0.0026201800440417583\
$f_D* \\sqrt\{8t_0\} = $0.2068623601937379 +/- 0.003546907525867658\
$f_Ds* \\sqrt\{8t_0\} = $0.2068623601937379 +/- 0.003546907525867658\
$f_eta \\sqrt\{8t_0\} = $0.8669858354378843 +/- 0.010734940135830663\
$f_jpsi \\sqrt\{8t_0\} = $0.402925306634495 +/- 0.008184816863154212\
$\
$J303\
$Z^\{tm\}_M \\mu_c \\sqrt\{8t_0\} = $3.106805094851814 +/- 0.07428519131808012\
$M_D \\sqrt\{8t_0\} = $3.951817203146415 +/- 0.05634686437630071\
$M_Ds \\sqrt\{8t_0\} = $4.091122020600407 +/- 0.05584799044824401\
$M_D* \\sqrt\{8t_0\} = $4.220571557429657 +/- 0.050571226554394465\
$M_Ds* \\sqrt\{8t_0\} = $4.378430678774368 +/- 0.05075096037526913\
$M_eta* \\sqrt\{8t_0\} = $6.222653661097548 +/- 0.09980503048093914\
$M_jpsi* \\sqrt\{8t_0\} = $6.45871173311589 +/- 0.09688893871726818\
$f_D \\sqrt\{8t_0\} = $0.46790893343380957 +/- 0.0035987447382407898\
$f_Ds \\sqrt\{8t_0\} = $0.5204463638699105 +/- 0.0025343572218284698\
$f_D* \\sqrt\{8t_0\} = $0.17563003352045548 +/- 0.004603348336097794\
$f_Ds* \\sqrt\{8t_0\} = $0.20465609887290984 +/- 0.0030029574421112505\
$f_eta \\sqrt\{8t_0\} = $0.876239237941407 +/- 0.008803719053134774\
$f_jpsi \\sqrt\{8t_0\} = $0.4070411608133558 +/- 0.006650155013434349\
}
\ No newline at end of file
gevp tools/tools.jl deleted 100644 → 0
View file @ 3e20520f
function read_dat(ens::String, g1::String="G5", g2::String="G5")
path = joinpath(path_data, ens)
rep = readdir(path, join=true)
aux = read_mesons.(rep, g1, g2)
res = []
for j =1:length(aux[1])
aux2 = [aux[i][j] for i =1:length(aux)]
push!(res, aux2)
end
return res
end
function get_mu(mu_list::Vector{Vector{Float64}}, deg::Bool)
mu_sorted = unique(sort(minimum.(mu_list)))
mul = mu_sorted[1]
deg ? mus = 0.0 : mus = mu_sorted[2]
muh = unique(maximum.(mu_list))
muh = filter(x-> x > mul && x > mus, muh)
return mul, mus, muh
end
function get_lh(mu_list, obs::Array{juobs.Corr}, deg::Bool)
mul, mus, muh = get_mu(mu_list, deg)
obs_lh = Array{juobs.Corr}(undef, 0)
for k = 1:length(mu_list)
mu = mu_list[k]
if mul in mu && mu[1] != mu[2] && !(mus in mu) #l-h
push!(obs_lh, obs[k])
end
end
return obs_lh
end
function get_lh(mu_list, obs::Vector{Vector{juobs.Corr}}, deg::Bool)
mul, mus, muh = get_mu(mu_list, deg)
obs_lh = Array{Array{juobs.Corr}}(undef, length(obs), length(muh))
for n_obs = 1:length(obs)
obs_lh[n_obs] = get_lh(mu_list, obs[n_obs], deg)
end
return obs_lh
end
function get_sh(mu_list, obs::Array{juobs.Corr}, deg::Bool)
mul, mus, muh = get_mu(mu_list, deg)
obs_sh = Array{juobs.Corr}(undef, 0)
for k = 1:length(mu_list)
mu = mu_list[k]
if mus in mu && mu[1] != mu[2] && !(mul in mu)#s-h
push!(obs_sh, obs[k])
end
end
return obs_sh
end
function get_sh(mu_list, obs::Vector{Vector{juobs.Corr}}, deg::Bool)
mul, mus, muh = get_mu(mu_list, deg)
obs_sh = Array{Array{juobs.Corr}}(undef, length(obs), length(muh))
for n_obs = 1:length(obs)
obs_sh[n_obs] = get_sh(mu_list, obs[n_obs], deg)
end
return obs_sh
end
function get_hh(mu_list, obs::Array{juobs.Corr}, deg::Bool)
mul, mus, muh = get_mu(mu_list, deg)
obs_hh = Array{juobs.Corr}(undef,0)
for k = 1:length(mu_list)
mu = mu_list[k]
for i =1:length(muh)
if muh[i] in mu && mu[1] == mu[2]#h-h
push!(obs_hh, obs[k])
end
end
end
return obs_hh
end
function get_hh(mu_list, obs::Vector{Vector{juobs.Corr}}, deg::Bool)
mul, mus, muh = get_mu(mu_list, deg)
obs_hh = Array{Array{juobs.Corr}}(undef, length(obs), length(muh))
for n_obs = 1:length(obs)
obs_hh[n_obs] = get_hh(mu_list, obs[n_obs], deg)
end
return obs_hh
end
mutable struct infoMatt
mat::Vector{Matrix}
mu::Vector{Float64}
y0::Int64
function infoMatt(a11::Array{juobs.Corr}, a12::Array{juobs.Corr}, a22::Array{juobs.Corr})
mu = getfield(a11[1], :mu)
y0 = Int64(getfield(a11[1], :y0))
elem11 = (sum(a11[i].obs for i =1:length(a11)))
elem12 = (sum(a12[i].obs for i =1:length(a12)))
elem22 = (sum(a22[i].obs for i =1:length(a22)))
diag = [elem11, elem22]
subdiag =[ elem12 ]
#diag = [a11.obs, a22[1].obs.+a22[2].obs.+a22[3].obs]
#subdiag = [a12[1].obs.+a12[2].obs.+a12[3].obs]
mat = get_matrix(diag, subdiag)
return new(mat, mu, y0)
end
function infoMatt(a11::Array{juobs.Corr}, a12::Array{juobs.Corr}, a13::Array{juobs.Corr}, a22::Array{juobs.Corr}, a23::Array{juobs.Corr}, a33::Array{juobs.Corr})
mu = getfield(a11[1], :mu)
y0 = Int64(getfield(a11[1], :y0))
elem11 = sum(a11[i].obs for i =1:length(a11))
elem12 = sum(a12[i].obs for i =1:length(a12))
elem22 = sum(a22[i].obs for i =1:length(a22))
elem13 = sum(a13[i].obs for i =1:length(a13))
elem23 = sum(a23[i].obs for i =1:length(a23))
elem33 = sum(a33[i].obs for i =1:length(a33))
diag = [elem11, elem22, elem33]
subdiag =[ elem12, elem13, elem23 ]
#diag = [a11.obs, a22[1].obs.+a22[2].obs.+a22[3].obs]
#subdiag = [a12[1].obs.+a12[2].obs.+a12[3].obs]
mat = get_matrix(diag, subdiag)
return new(mat, mu, y0)
end
end
function comp_mat(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}})
mu = getfield.(tot11[1],:mu)
res = Array{infoMatt}(undef, length(mu))
for i = 1:length(mu)
a11 = [ tot11[j][i] for j = 1:size(tot11,1) ]
a12 = [ tot12[j][i] for j = 1:size(tot12,1) ]
a22 = [ tot22[j][i] for j = 1:size(tot22,1) ]
res[i] = infoMatt(a11, a12, a22)
#res[i]= infoMatt(tot11[i], [tot12[1][i],tot12[2][i],tot12[3][i]], [tot22[1][i],tot22[2][i],tot22[3][i]])
end
return res
end
function comp_mat(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr}}, tot13::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}}, tot33::Array{Array{juobs.Corr}}, tot23::Array{Array{juobs.Corr}})
mu = getfield.(tot11[1],:mu)
res = Array{infoMatt}(undef, length(mu))
for i = 1:length(mu)
a11 = [ tot11[j][i] for j = 1:size(tot11,1) ]
a12 = [ tot12[j][i] for j = 1:size(tot12,1) ]
a22 = [ tot22[j][i] for j = 1:size(tot22,1) ]
a33 = [ tot33[j][i] for j = 1:size(tot33,1) ]
a13 = [ tot13[j][i] for j = 1:size(tot13,1) ]
a23 = [ tot23[j][i] for j = 1:size(tot23,1) ]
res[i] = infoMatt(a11, a12, a13, a22, a23, a33)
#res[i]= infoMatt(tot11[i], [tot12[1][i],tot12[2][i],tot12[3][i]], [tot22[1][i],tot22[2][i],tot22[3][i]])
end
return res
end
function comp_energy(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}}, evec::Bool=false; tnot::Int64=3)
aux_mat = comp_mat(tot11, tot12, tot22)
y0 = getfield(aux_mat[1], :y0)
res_en = Array{infoEn}(undef, length(aux_mat))
if !evec
for i = 1:length(aux_mat)
mat_obs = getfield(aux_mat[i], :mat)[y0+1:end-1]
evalues = uwgevp_tot(mat_obs, tnot)
res_en[i] = infoEn(energies(evalues), aux_mat[i])
end
return res_en
else
evec = Array{Array{}}(undef, length(aux_mat))
for i = 1:length(aux_mat)
mat_obs = getfield(aux_mat[i], :mat)[y0+1:end-1]
evalues, eigvec = uwgevp_tot(mat_obs, tnot, evec=true, iter=50)
res_en[i] = infoEn(energies(evalues), aux_mat[i])
evec[i] = eigvec
end
return res_en, evec
end
end
function comp_energy(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr}}, tot13::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}}, tot33::Array{Array{juobs.Corr}}, tot23::Array{Array{juobs.Corr}}, evec::Bool=false; tnot::Int64=2)
aux_mat = comp_mat(tot11, tot12, tot13, tot22, tot23, tot33)
y0 = getfield(aux_mat[1], :y0)
res_en = Array{infoEn}(undef, length(aux_mat))
if !evec
for i = 1:length(aux_mat)
mat_obs = getfield(aux_mat[i], :mat)[y0+1:end-1]
evalues = uwgevp_tot(mat_obs, tnot, iter=50)
res_en[i] = infoEn(energies(evalues), aux_mat[i])
end
return res_en
else
evec = Array{Array{}}(undef, length(aux_mat))
for i = 1:length(aux_mat)
mat_obs = getfield(aux_mat[i], :mat)[y0+1:end-1]
evalues, eigvec = uwgevp_tot(mat_obs, tnot, evec=true)
println("in comp_energy, i is = ", i)
res_en[i] = infoEn(energies(evalues), aux_mat[i])
evec[i] = eigvec
end
return res_en, evec
end
end
mutable struct infoEn
en_val::Array{Array{uwreal}}
mu::Vector{Float64}
y0::Int64
function infoEn(energ::Array{Array{uwreal}},matr::infoMatt )
en_val = energ
mu = getfield(matr,:mu)
y0 = getfield(matr,:y0)
new(en_val, mu, y0)
end
end
function pseudo_mat_elem(evec::Array{Array{T,2} where T,1}, mass::uwreal, mu::Array{Float64})
Rn = [exp(mass * (t)/2) * evec[t][1] for t =2:length(evec)]
aux = sqrt(2)*sum(mu) / mass^1.5
return uwdot(aux,Rn)
end
function vec_mat_elem(evec::Array{Array{T,2} where T,1}, mass::uwreal, deg::Bool)
if !deg
Rn = [exp(mass * t/2) * evec[t][1] for t =1:length(evec)]
return uwdot(Rn , za(beta[iens]) * sqrt(2/ mass))
else
Rn = [exp(mass * t/2) * evec[t][4] for t =1:length(evec)]
println(za(beta[iens]))
println(mass)
return uwdot(Rn , za(beta[iens]) * sqrt(2/ mass))
#return uwdot(Rn , 1/mass)# * sqrt(2/ mass))
end
end
function extract_dec_const(mat_elem::uwreal, mass::uwreal, mu::Array{Float64})
return sqrt(2)*sum(mu)*mat_elem / mass^1.5
end
function match_muc(muh, m_lh, m_lh_star, target)
M = (m_lh .+ 3 .* m_lh_star) ./ 4
par, chi2exp = lin_fit(muh, M)
muh_target = x_lin_fit(par, target)
return muh_target
end
function match_muc(muh, m_lh, m_sh, m_lh_star, m_sh_star, target)
M = (2 .* m_lh .+ 6 .* m_lh_star .+ m_sh .+ 3 .* m_sh_star) ./ 12
par, chi2exp = lin_fit(muh, M)
muh_target = x_lin_fit(par, target)
return muh_target
end
function match_muc_flav(muh, m_lh, m_sh, target)
M = (2/3 .* m_lh .+ 1/3*m_sh )
#M = m_lh
par, chi2exp = lin_fit(muh, M)
muh_target = x_lin_fit(par, target)
return muh_target
end
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment