Skip to content

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

major update to the twisted mass charm analysis code

parent 6c956856
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 2227 additions and 352 deletions
plat/plat_tm_dps.txt
View file @ 9fdef14c
......@@ -59,11 +59,11 @@ sh 83 110
hh 90 104
t0 20 110
#D200
ll 85 105
ll 80 87
ls 82 100
lh 80 88
lh 86 92
ss 85 100
sh 82 88
sh 86 92
hh 87 93
t0 20 110
#J303
......
plat/plat_tm_mps.txt
View file @ 9fdef14c
......@@ -61,11 +61,19 @@ t0 20 110
#D200
ll 79 86
ls 85 105
lh 84 87
lh 78 88
ss 85 105
sh 85 87
hh 90 93
sh 80 85
hh 80 90
t0 25 100
#N302
ll 87 103
ls 85 105
lh 88 100
ss 85 105
sh 88 110
hh 88 110
t0 20 110
#J303
ll 120 148
ls 123 140
......@@ -74,4 +82,28 @@ ss 123 140
sh 130 145
hh 133 175
t0 25 170
#E300
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
lh 125 145
ss 123 140
sh 130 145
hh 133 175
t0 25 170
#J501
ll 120 148
ls 123 140
lh 125 145
ss 123 140
sh 130 145
hh 133 175
t0 25 170
#end
src/aux obs computation/md.jl
View file @ 9fdef14c
......@@ -24,11 +24,11 @@ const path_plat_w = "/Users/alessandroconigli/Lattice/data/aux_obs_data/plat_wil
const path_rw = "/Users/alessandroconigli/Lattice/data/aux_obs_data/rwf"
const path_md = "/Users/alessandroconigli/Lattice/data/aux_obs_data/md"
const ens_list = ["D200"] # ["H101", "H102r002", "H400", "N202", "N200", "N203", "N300", "J303"]
const ens_list = ["H101", "H102r002", "H105", "H400", "N202", "N203", "N200", "D200", "N300", "J303"] # ["H101", "H102r002", "H400", "N202", "N200", "N203", "N300", "J303"]
const path_bdio_md = "/Users/alessandroconigli/Lattice/data/bdio_charm/md"
const path_bdio_w = "/Users/alessandroconigli/Lattice/data/bdio_charm/wilson"
const phi4_target = uwreal([1.117, 0.012], "phi4_target")
const phi4_target = uwreal([1.098, 0.010], "phi4_target")
const rwf = true
#========= INCLUDES ==========#
......@@ -49,15 +49,16 @@ for i in 1:length(ens_list)
end
#t0 [2.86, 3.659, 5.164, 8.595]
wpmm = Dict{String, Vector{Float64}}()
wpmm["H101"] = [-1.0, 4.0, -1.0, 14.0*2.86]
wpmm["H102r002"] = [-1.0, 4.0, -1.0, 14.0*2.86]
wpmm["H105"] = [-1.0, 4.0, -1.0, 14.0*2.86]
wpmm["H400"] = [-1.0, 4.0, -1.0, 14.0*3.659]
wpmm["N202"] = [-1.0, 4.0, -1.0, 14.0*5.164]
wpmm["N200"] = [-1.0, 4.0, -1.0, 14.0*5.164]
wpmm["N203"] = [-1.0, 4.0, -1.0, 14.0*5.164]
wpmm["N300"] = [-1.0, 4.0, -1.0, 14.0*8.595]
wpmm["J303"] = [-1.0, 4.0, -1.0, 14.0*8.595]
wpmm["H101"] = [-1.0, 2.0, -1.0, 14.0*2.86]
wpmm["H102r002"] = [-1.0, 2.0, -1.0, 14.0*2.86]
wpmm["H105"] = [-1.0, 2.0, -1.0, 14.0*2.86]
wpmm["H400"] = [-1.0, 2.0, -1.0, 14.0*3.659]
wpmm["N202"] = [-1.0, 2.0, -1.0, 14.0*5.164]
wpmm["N200"] = [-1.0, 2.0, -1.0, 14.0*5.164]
wpmm["N203"] = [-1.0, 2.0, -1.0, 14.0*5.164]
wpmm["D200"] = [-1.0, 2.0, -1.0, 14.0*5.164]
wpmm["N300"] = [-1.0, 2.0, -1.0, 14.0*8.595]
wpmm["J303"] = [-1.0, 2.0, -1.0, 14.0*8.595]
##
#======== ANALYSIS ==========#
gen_mulist = Vector{Vector{Vector{Float64}}}(undef, length(ens_list))
......@@ -68,6 +69,7 @@ obs_tm_shifted = Vector(undef, length(ensinfo))
#obs
obs_w = Vector(undef, length(ensinfo))
for (k, ens) in enumerate(ensinfo)
println(ens)
# t0
t0_ens, YW, W = get_t0(path_data_w, ens, rw=rwf, path_rw=path_rw, pl=true)
......
src/aux obs computation/md_old.jl 0 → 100644
View file @ 9fdef14c
#=
# The code produces BDIO files for md_info (see OrderedDict md_info) and Wilson observables (see OrderedDict obs_w)
# The BDIO files are saved in path_bdio_md and path_bdio_w respectivily. Each BDIO file contains information about a given ensemble.
# Plateaux for different ensembles are defined in txt files (path_plat_w)
# The dat files must be located in the following way
path_data_w/ensemble -> *mesons.dat (Wilson)
path_data_w/ensemble -> *ms.dat
path_rw/ensemble -> *.ms1.dat
path_md/ensemble -> *.pbp.dat
# Example: path_rw/H400 contains H400r001.ms1.dat and H400r002.ms1.dat
# The target value of phi4 can be modified
# Masses and decay constants are in terms of the reference scale t0
# See read_bdio.jl for more information
# WARNING
# here there is a possible bug: symmetric point ensembles are not shifted at a fix value of phi_2 as it should be
# please use md.jl
=#
using Base: String
using LaTeXStrings: length
using OrderedCollections
using juobs, BDIO, DelimitedFiles, ADerrors, LaTeXStrings
#================== PATHS + INFO ==================#
const path_data_w = "/Users/alessandroconigli/Lattice/data/aux_obs_data/wilson"
const path_plat_w = "/Users/alessandroconigli/Lattice/data/aux_obs_data/plat_wilson.txt"
const path_rw = "/Users/alessandroconigli/Lattice/data/aux_obs_data/rwf"
const path_md = "/Users/alessandroconigli/Lattice/data/aux_obs_data/md"
const ens_list = ["N300"] # ["H101", "H102r002", "H400", "N202", "N200", "N203", "N300", "J303"]
const path_bdio_md = "/Users/alessandroconigli/Lattice/data/bdio_charm/md"
const path_bdio_w = "/Users/alessandroconigli/Lattice/data/bdio_charm/wilson"
const phi4_target = uwreal([1.098, 0.010], "phi4_target")
const rwf = true
#========= INCLUDES ==========#
include("./types.jl")
include("./tools.jl")
include("./const.jl")
#======== GET ENSEMBLE INFORMATION 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
#t0 [2.86, 3.659, 5.164, 8.595]
wpmm = Dict{String, Vector{Float64}}()
wpmm["H101"] = [-1.0, 4.0, -1.0, 14.0*2.86]
wpmm["H102r002"] = [-1.0, 4.0, -1.0, 14.0*2.86]
wpmm["H105"] = [-1.0, 4.0, -1.0, 14.0*2.86]
wpmm["H400"] = [-1.0, 4.0, -1.0, 14.0*3.659]
wpmm["N202"] = [-1.0, 4.0, -1.0, 14.0*5.164]
wpmm["N200"] = [-1.0, 4.0, -1.0, 14.0*5.164]
wpmm["N203"] = [-1.0, 4.0, -1.0, 14.0*5.164]
wpmm["N300"] = [-1.0, 4.0, -1.0, 14.0*8.595]
wpmm["J303"] = [-1.0, 4.0, -1.0, 14.0*8.595]
##
#======== ANALYSIS ==========#
gen_mulist = Vector{Vector{Vector{Float64}}}(undef, length(ens_list))
md_info = Vector(undef, length(ensinfo))
obs_tm_shifted = Vector(undef, length(ensinfo))
#obs
obs_w = Vector(undef, length(ensinfo))
for (k, ens) in enumerate(ensinfo)
# t0
t0_ens, YW, W = get_t0(path_data_w, ens, rw=rwf, path_rw=path_rw, pl=true)
# WILSON correlators
pp_w, pp_wW, _ = get_corr(path_data_w, ens, "G5", "G5", rw=rwf, path_rw=path_rw, legacy=true)
pp_w_d1, _, _ = get_corr(path_data_w, ens, "G5_d1", "G5_d1", rw=rwf, path_rw=path_rw, legacy=true)
pp_w_d2, _, _ = get_corr(path_data_w, ens, "G5_d2", "G5_d2", rw=rwf, path_rw=path_rw, legacy=true)
# WILSON masses
m_w = get_meff_w(pp_w, ens, pl=true)
kappa_list = getfield.(pp_w, :kappa)
mll_w = get_ll_w(kappa_list, m_w, ens.deg)
mls_w = ens.deg ? mll_w : get_ls_w(kappa_list, m_w, ens.deg)
phi2 = 8 * t0_ens * mll_w^2
phi4 = ens.deg ? 8 * t0_ens * 1.5 * mll_w^2 : 8 * t0_ens * (0.5 * mll_w^2 + mls_w^2)
# MASS SHIFTS
dSdm = read_pbp(path_md, ens.id)
ppll = get_ll_w(kappa_list, pp_w, ens.deg)
ppllW = get_ll_w(kappa_list, pp_wW, ens.deg)
ppll_d1 = get_ll_w(kappa_list, pp_w_d1, ens.deg)
ppll_d2 = get_ll_w(kappa_list, pp_w_d2, ens.deg)
dphi4_s1, dphi4_s2 = md_sea(phi4, dSdm, YW, W) .+ md_sea(phi4, dSdm, ppllW, W)
dphi2_s1, dphi2_s2 = md_sea(phi2, dSdm, YW, W) .+ md_sea(phi2, dSdm, ppllW, W)
dphi4_v1, dphi4_v2 = md_val(phi4, ppll, [ppll_d1, ppll_d2])
dphi2_v1, dphi2_v2 = md_val(phi2, ppll, [ppll_d1, ppll_d2])
if !ens.deg
ppls = get_ls_w(kappa_list, pp_w, ens.deg)
pplsW = get_ls_w(kappa_list, pp_wW, ens.deg)
ppls_d1 = get_ls_w(kappa_list, pp_w_d1, ens.deg)
ppls_d2 = get_ls_w(kappa_list, pp_w_d2, ens.deg)
aux1, aux2 = md_val(phi4, ppls, [ppls_d1, ppls_d2])
dphi4_v1 = dphi4_v1 + aux1
dphi4_v2 = dphi4_v2 + aux2
dphi4_s1, dphi4_s2 = (dphi4_s1, dphi4_s2) .+ md_sea(phi4, dSdm, pplsW, W)
end
# Dict: Deltam, phi2, phi4, t0 (shifted)
deltam = (phi4_target - phi4) / (2*dphi4_s1 + dphi4_s2 + dphi4_v1 + dphi4_v2)
dt0_s1, dt0_s2 = md_sea(t0_ens, dSdm, YW, W)
t0_ens_shifted = t0_ens + deltam * (2 * dt0_s1 + dt0_s2)
phi2_shifted = phi2 + deltam * (2*dphi2_s1 + dphi2_s2 + dphi2_v1 + dphi2_v2)
phi4_shifted = phi4 + deltam * (2*dphi4_s1 + dphi4_s2 + dphi4_v1 + dphi4_v2)
md_info[k] = OrderedDict(
"deltam" => deltam,
"t0_shifted" => t0_ens_shifted,
"phi2_shifted" => phi2_shifted,
"phi4_shifted" => phi4_shifted
)
obs_w[k] = OrderedDict(
"mll" => sqrt(8 * t0_ens) * mll_w,
"mls" => sqrt(8 * t0_ens) * mls_w,
"t0" => t0_ens,
"phi2" => phi2,
"phi4" => phi4
)
end
## WRITE BDIO
# MD INFO
for (k, obs) in enumerate(md_info) #loop over ensembles
ens = ensinfo[k].id
p = joinpath(path_bdio_md, string(ens, ".bdio"))
fb = BDIO_open(p, "w", ens)
uinfo = 0
for (key, value) in obs #loop over observables
write_uwreal(value, fb, uinfo)
uinfo += 1
end
BDIO_close!(fb)
end
# WILSON
for (k, obs) in enumerate(obs_w) #loop over ensembles
ens = ensinfo[k].id
p = joinpath(path_bdio_w, string(ens, ".bdio"))
fb = BDIO_open(p, "w", ens)
uinfo = 0
for (key, value) in obs #loop over observables
write_uwreal(value, fb, uinfo)
uinfo += 1
end
BDIO_close!(fb)
end
\ No newline at end of file
src/aux obs computation/tools.jl
View file @ 9fdef14c
......@@ -23,7 +23,8 @@ function read_rw(path::String, ens::String)
try
length(rep) == 1 ? (return read_ms1(rep[1])) : (return read_ms1.(rep))
catch
length(rep) == 1 && length(rep)!=0 ? (return read_ms1(rep[1], v="1.4")) : (return read_ms1.(rep, v="1.4"))
ens == "H400" ? vv = "1.6" : vv="1.4"
length(rep) == 1 && length(rep)!=0 ? (return read_ms1(rep[1], v=vv)) : (return read_ms1.(rep, v=vv))
end
else
error("ms1.dat file not found for ensemble ", ens, " in path ", p)
......
src/gevp analysis/D200_gevp_test.jl 0 → 100644
View file @ 9fdef14c
This diff is collapsed.
src/gevp analysis/D200_test_noise_sources_ncnfg.jl 0 → 100644
View file @ 9fdef14c
using juobs, BDIO, DelimitedFiles, ADerrors, LaTeXStrings, PyPlot
using LaTeXStrings: length
using LsqFit
rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["text.usetex"] = false
rcParams["mathtext.fontset"] = "cm"
rcParams["font.size"] =10
rcParams["axes.labelsize"] =22
rcParams["axes.titlesize"] = 18
plt.rc("text", usetex=false) # set to true once you install latex
path_corr = "/Users/alessandroconigli/Lattice/data/charm_full_Dirac/N302/N302r001_correction_k0.137211_mup0.00292_mup0.00682_mup0.1320_mup0.1460.mesons.dat"
path_sl = "/Users/alessandroconigli/Lattice/data/charm_full_Dirac/N302/N302r001_sloppy_k0.137211_mup0.00292_mup0.00682_mup0.1320_mup0.1460.mesons.dat"
path_rw = "/Users/alessandroconigli/Lattice/data/aux_obs_data/rwf/N302/N302r001.ms1.dat"
ens = "N302"
##
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["H400"] = [-1.0, 1.5, -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["N300"] = [-1.0, 2.0, -1.0, -1.0]
wpmm["N302"] = [-1.0, 2.0, -1.0, -1.0]
wpmm["J303"] = [-1.0, 2.0, -1.0, -1.0]
wpmm["D200"] = [-1.0, 4.0, -1.0, -1.0]
##
# D200 param
# noise_trunc = [nothing, 1, 2, 4, 8]
# ncfg = [2000, 2000, 1000, 500, 250]
# idmsl = [collect(1:2000), collect(1:2000), collect(1:2:2000), collect(1:4:2000), collect(1:8:2000)]
noise_trunc = [16, 2, 4, 8]
ncfg = [60, 60, 30, 15]
idmsl = [collect(1:60), collect(1:60), collect(1:2:60), collect(1:4:60)]
## here i want to test, with full statistics, the TSM, sloppy, exact and exact-sloppy,corr observables
rwf = juobs.read_ms1(path_rw)
idm_corr= collect(1:20)
nms = 60
data_corr = juobs.read_mesons_correction(path_corr,"G5", "G5") # correction data combined
data_ex = juobs.read_mesons(path_corr, "G5", "G5" )[1:9] # sloppy data only
data_sl = juobs.read_mesons(path_sl, "G5", "G5") # exact data only
pp_corr = juobs.corr_obs.(data_corr, idm=idm_corr, nms=nms) # correction correlator
pp_ex = juobs.corr_obs.(data_ex, idm=idm_corr, nms=nms) # exact correlator
pp_sloppy = juobs.corr_obs.(data_sl, idm=idmsl[1], nms=nms) # sloppy correlator
pp_TSM = juobs.corr_obs_TSM.(data_sl, data_corr, idm_corr=idm_corr, idm_sl=idmsl[1], nms=nms) # exact correlator
for midx in [1,4,5,6]
TT = 75
tmp_corr = pp_corr[midx].obs[TT] ; uwerr(tmp_corr, wpmm)
tmp_ex = pp_ex[midx].obs[TT] ; uwerr(tmp_ex, wpmm)
tmp_sloppy = pp_sloppy[midx].obs[TT] ; uwerr(tmp_sloppy, wpmm)
tmp_TSM = pp_TSM[midx].obs[TT] ; uwerr(tmp_TSM, wpmm)
println("#############")
println("Midx = $(midx)")
println("TSM = $(tmp_TSM)")
println("Sloppy = $(tmp_sloppy)")
println("Exact = $(tmp_ex)")
println("Correc = $(tmp_corr)")
println("#############")
end
## here i am studing the tradeoff between noise sources and number of configurations.
ll = Vector()
pp = Vector{Vector{juobs.Corr}}(undef, 0)
for k in eachindex(noise_trunc)
println(k)
data_sl = juobs.read_mesons(path_sl, "G5", "G5", nnoise_trunc=noise_trunc[k])
ppaux = juobs.corr_obs_TSM.(data_sl, data_corr, idm_corr=idm_corr, idm_sl=idmsl[1], nms=nms )
push!(pp, ppaux )
llaux = string("cnfg:",ncfg[k], " nnoise:", noise_trunc[k])
push!(ll, llaux)
end
## pp corr vs cnfg and nnoise
PPfix = Vector{Vector{uwreal}}(undef, length(pp))
PPfixex = Vector{Vector{uwreal}}(undef, length(pp))
for k in eachindex(pp)
PPfix[k] = Vector{uwreal}(undef, 0)
PPfixex[k] = Vector{uwreal}(undef, 0)
for midx in [1,4,5,6]
yy = pp[k][midx].obs[75]
yyex = pp_ex[midx].obs[75]
push!(PPfix[k], yy)
push!(PPfixex[k], yyex)
end
end
tt = ["PPll", "PPlh", "PPsh", "PPhh"]
for k in eachindex(PPfix[1])
fig = figure(figsize=(10,7.5))
yy = getindex.(PPfix, k)
yyex = getindex.(PPfixex, k)
[uwerr(yy[ii], wpmm) for ii in eachindex(yy)]
[uwerr(yyex[ii], wpmm) for ii in eachindex(yyex)]
xx = collect(1:length(yy))
errorbar(xx, value.(yy), err.(yy), fmt="*", capsize=2, label="TSM" )
errorbar(xx.+0.1, value.(yyex), err.(yyex), fmt="*", capsize=2, label="Exact" )
xticks(xx, string.(ll), rotation=15)
title(tt[k])
legend()
display(fig)
# savefig(joinpath("/Users/alessandroconigli/Lattice",string(tt[k],"_corrppNoiseCnfg.pdf")))
close()
end
## print corr values at all nnoise and cnfg for all masses
for k in eachindex(PPfix[1])
yy = getindex.(PPfix, k)
[uwerr(yy[ii], wpmm) for ii in eachindex(yy)]
for elem in yy
println(elem)
end
println("\n")
end
## meff vs cnfg and nnoise
M = Vector{Vector{uwreal}}(undef, length(pp))
for k in eachindex(pp)
M[k] = Vector{uwreal}(undef, 0)
cc = 1
for midx in [1,4,5,6]
if cc == 1
tmin =collect(75:85)
tmax = collect(105:115)
tmax = tmax[1:2:end]
plat = [85,100]
else
tmin = cc==4 ? collect(86:88) : collect(78:84)
tmax = cc==4 ? collect(90:92) : collect(95:100)
plat = cc==4 ? [88,94] : [85,95]
end
@. fitt(x,p) = p[1]
maux = juobs.meff(pp[k][midx], plat, pl=true )
# maux, syst, all_res, FitDict = juobs.bayesian_av(fitt, data, tmin, tmax, 2, pl=true, data=true, plt_title="ciao", path_plt=nothing, label=L"$m$")
push!(M[k], maux)
cc+=1
end
end
##
# plot
tt = ["mll", "mlh", "msh", "mhh"]
for k in eachindex(M[1])
fig = figure(figsize=(10,7.5))
yy = getindex.(M, k)
xx = collect(1:length(yy))
errorbar(xx, value.(yy), err.(yy), fmt="*", capsize=2 )
xticks(xx, string.(ll), rotation=15)
title(tt[k])
display(fig)
savefig(joinpath("/Users/alessandroconigli/Lattice",string(tt[k],"_testNoiseCnfg.pdf")))
close()
end
##
idxm = 4
t0= meff(pp[1][idxm], [85,94])
t1 = meff(pp[2][idxm], [85,95])
t2 = meff(pp[3][idxm], [85,95])
t3 = meff(pp[4][idxm], [85,95])
t4 = meff(pp[5][idxm], [85,95])
\ No newline at end of file
src/gevp analysis/const.jl
View file @ 9fdef14c
#========= ENSEMBLE DATABASE ========#
ens_db = Dict(
#"ens_id"=>[L, beta, is_deg?, m_pi, dtr]
"H102r002" => [32, 3.4, false, 0.15306, 2],
#"ens_id"=>[L, beta, is_deg?, m_pi, dtr, idm_corr]
"H102r002" => [32, 3.4, false, 0.15306, 2, nothing],
#"H102r001" => [32, 3.4, false, 0.15306, 2],
"H105" => [32, 3.4, false, 0.12151, 2],
"H101" => [32, 3.4, true, 0.17979, 2],
"H400" => [32, 3.46, true, 0.16345, 1],
"N200" => [48, 3.55, false, 0.09222, 1],
"N202" => [48, 3.55, true, 0.13407, 2],
"N203" => [48, 3.55, false, 0.11224, 1],
"D200" => [64, 3.55, false, 0.06611, 2],
"N300" => [48, 3.70, true, 0.10630, 1],
"J303" => [64, 3.70, false, 0.06514, 2]
)
"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],
"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],
"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],
"J500" => [64, 3.85, true, 0.00000, 2, nothing],
"J501" => [64, 3.85, false, 0.00000, 2, nothing],
)
trunc_db = Dict(
"H102r002" => nothing,
#"H102r001" => nothing,
"H105" => nothing,
"H101" => nothing,
"H105r002" => nothing,
"H101" => [1001,1009],
"H400" => nothing,
"N200" => nothing,
"N202" => nothing,
"N203" => nothing,
"D200" => 1000,
"N300" => nothing,
"J303" => nothing
"D200" => 645, #1000
"N300" => 1279,
"N302" => nothing,
"J303" => 721,
"E300" => nothing,
"J500" => [745,649], #[751, 655],
"J501" => nothing
)
# total number of cnfg for each ens
ens_nms = Dict(
#"H102r001" => 1029,
"H102r002" => 1008,
"H105r002" => 1042,
"H101" => 2016,
"H400" => 1045,
"N200" => 1712,
"N202" => 899,
"N203" => 1543,
"D200" => 2001,
"N300" => 1540,
"N302" => 2201,
"J303" => 1073,
"E300" => 1139
)
#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)
......@@ -39,8 +65,8 @@ const ZM_error = [35, 27, 33, 38, 45] .* 1e-4
const ZM_tm_data = [2.6047, 2.6181, 2.6312, 2.6339, 2.6127]
const ZM_tm_error = [42, 33, 42, 48, 55] .* 1e-4
#1608.08900
const t0_data = [2.86, 3.659, 5.164, 8.595]
const t0_error = [11, 16, 18, 29] .* 1e-3
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
# 1808.09236
......@@ -49,7 +75,7 @@ const ZA_err = [72, 93, 43, 47, 47] .*1e-5
#Covariance matrices (Uncorrelated)
const C1 = zeros(5, 5)
const C2 = zeros(5, 5)
const C3 = zeros(4, 4)
const C3 = zeros(5, 5)
const C4 = zeros(6,6)
const C5 = zeros(5, 5)
for i = 1:6
......@@ -58,9 +84,7 @@ for i = 1:6
C1[i, i] = ZM_error[i] ^ 2
C2[i, i] = ZM_tm_error[i] ^ 2
C5[i, i] = ZA_err[i]^2
if i<=4
C3[i,i] = t0_error[i] ^ 2
end
C3[i,i] = t0_error[i] ^ 2
end
end
......@@ -84,6 +108,6 @@ Mrat = uwreal([0.9148, 0.0088], "M/mhad")
# Renormalization constants, t0 and a as functions (f(beta))
zm(beta::Float64) = ZM[b_values .== beta][1]
zm_tm(beta::Float64) = Mrat / ZP(beta)
t0(beta::Float64) = t0_[b_values2 .== beta][1]
a(beta::Float64) = a_[b_values2 .== beta][1]
t0(beta::Float64) = t0_[b_values .== beta][1]
a(beta::Float64) = a_[b_values .== beta][1]
za(beta::Float64) = Za[b_values .== beta][1]
\ No newline at end of file
src/gevp analysis/func_comb.jl
View file @ 9fdef14c
......@@ -183,4 +183,197 @@ for n = 2:f_npar+1
end
f_npar_comb = 2 * f_npar
f_n_param_comb = collect(5:2:f_npar_comb+5)
f_n_param_comb_clog = collect(6:2:f_npar_comb+6)
\ No newline at end of file
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})
# 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]) .*
((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
((4. .* p[1]) ./ x[:,5]) .* (3 .* phi2_sym .- 3 .* x[:,2])
# higher order terms HQET - careful, explicit p[1]=\phi_f^2 dependece is discouraged
r_chiphi2_hqet(x) = (4 ./ x[:,5].^2) .* (3 .* phi2_sym .- 3 .* x[:,2]) .* x[:,2]
r_chiphih_hqet(x) = (4 ./ x[:,5]) .* (3 .* phi2_sym .- 3 .* x[:,2]) ./ x[:,3]
# base model lattice spacing dependece
r_basemodel_continuum(x,p) = 1. .+ p[2] .* x[:,1] .* (3 .* phi2_sym .- 3 .* x[:,2])
# 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_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_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_labels = [["hqet"]]
for extrapar = 1:r_hqet_extra_par
idx_f_expar = filter(x->sum(x)==extrapar, r_hqet_func_map)
# build and store temporary hqet dependence
tmp_hqet_func = Vector{Function}()
tmp_hqet_param = []
tmp_hqet_label = []
for ii in idx_f_expar
aux_f = (x,p) -> r_basemodel_hqet(x,p) .+ sum( p[r_hqet_base_par+k] .* r_hqet_func_list[ii][k](x) for k in eachindex(r_hqet_func_list[ii]))
push!(tmp_hqet_param, r_hqet_base_par + extrapar)
push!(tmp_hqet_func, aux_f)
push!(tmp_hqet_label, r_hqet_func_label[ii])
end
# add continuum dependence - already code in view of adding more continuum higher order terms
for (jj, f_hqet) in enumerate(tmp_hqet_func)
aux1 = (x,p) -> f_hqet(x,p) .* r_basemodel_continuum(x,p)
push!(r_hqet_models, aux1)
push!(r_hqet_param, tmp_hqet_param[jj])
push!(r_hqet_labels, tmp_hqet_label[jj])
for contpar = 1:r_cont_extra_par
idx_cont_expar = filter(x->sum(x)==contpar, r_cont_func_map)
for ii in idx_cont_expar
aux_f_cont = (x,p) -> f_hqet(x,p) .* (r_basemodel_continuum(x,p) .+ sum(p[r_hqet_base_par+extrapar+k] .* r_cont_func_list[ii][k](x) for k in eachindex(r_cont_func_list[ii]) ) )
push!(r_hqet_models, aux_f_cont)
push!(r_hqet_param, tmp_hqet_param[jj]+contpar)
push!(r_hqet_labels, vcat(tmp_hqet_label[jj], r_cont_func_label[ii] ))
end
end
end
end
# at the end of the above loop all the HQET combinations are taken into acount and store in r_hqet_models, r_hqet_param, r_hqet_labels
# now ration functional forms for Taylor based expressions
# base model taylor
r_basemodel_taylor(x,p) = 1. .+ 3 .*(phi2_sym .- x[:,2]) .*
(p[1] .+ p[2] .* ( 1. ./ x[:,3] .- 1. ./ x[:,4]) )
#higher order models taylor
r_taylor_nlo_phi2(x) = 3 .* (phi2_sym .- x[:,2]) .* x[:,2]
# this redefinition is required since in hqet LO parameter in continuum dependece is p[2] while here is p[3]
r_basemodel_continuum_taylor(x,p) = 1. .+ p[3] .* x[:,1] .* (3 .* phi2_sym .- 3 .* x[:,2])
r_taylor_func_list = [r_taylor_nlo_phi2]
r_taylor_func_label = ["pm(2)"]
r_taylor_func_map = [Bool.([i]) for i=0:1]
r_taylor_base_par = 3 # 2 from taylor base model, 1 from LO lattice spacing dependence
r_taylor_extra_par = length(r_taylor_func_list)
r_taylor_models = Vector{Function}(undef, 1)# r_basemodel_taylor + extra param
r_taylor_models[1] = (x,p) -> r_basemodel_taylor(x,p) .* r_basemodel_continuum_taylor(x,p)
r_taylor_param = [3]
r_taylor_labels = [["taylor"]]
for extrapar in 1:r_taylor_extra_par
idx_f_expar = filter(x-> sum(x)==extrapar, r_taylor_func_map)
# build and store temporary taylor dependence
tmp_taylor_func = Vector{Function}()
tmp_taylor_param = []
tmp_taylor_label = []
for ii in idx_f_expar
aux_f = (x,p) -> r_basemodel_taylor(x,p) .+ sum(p[r_taylor_base_par+k] .* r_taylor_func_list[ii][k](x) for k in eachindex(r_taylor_func_list[ii]))
push!(tmp_taylor_param, r_taylor_base_par +extrapar)
push!(tmp_taylor_func, aux_f)
push!(tmp_taylor_label, r_taylor_func_label[ii])
end
# add continuum dependence - already code in view of adding more continuum higher order terms
for (jj, f_taylor) in enumerate(tmp_taylor_func)
aux1 = (x,p) -> f_taylor(x,p) .* r_basemodel_continuum_taylor(x,p)
push!(r_taylor_models, aux1)
push!(r_taylor_param, tmp_taylor_param[jj])
push!(r_taylor_labels, tmp_taylor_label[jj])
for contpar = 1:r_cont_extra_par
idx_cont_expar = filter(x->sum(x)==contpar, r_cont_func_map)
for ii in idx_cont_expar
aux_f_cont = (x,p) -> f_taylor(x,p) .* (r_basemodel_continuum_taylor(x,p) .+ sum(p[r_taylor_base_par+extrapar+k] .* r_cont_func_list[ii][k](x) for k in eachindex(r_cont_func_list[ii]) ))
push!(r_taylor_models, aux_f_cont)
push!(r_taylor_param, tmp_taylor_param[jj]+contpar)
push!(r_taylor_labels, vcat(tmp_taylor_label[jj], r_cont_func_label[ii]))
end
end
end
end
# at the end this second loop all the TAYLOR combinations are taken into acount and store in r_taylor_models, r_taylor_param, r_taylor_labels
# final ration models are
r_tot_models = vcat(r_hqet_models, r_taylor_models)
r_tot_param = vcat(r_hqet_param, r_taylor_param)
r_tot_labels = vcat(r_hqet_labels, r_taylor_labels)
##
# Operator overloading to call functions with uwreal parameters
function Base.length(uw::uwreal)
return length(value(uw))
end
# function Base.iterate(uw::uwreal)
# uwerr(uw)
# return iterate(value(uw), 1)
# out = iterate(value(uw), args...; kwargs...)
# return (uw, 1)
# isnothing(out2) ? (return nothing) : (return (uw, out2)) #(uw, out2) #(uw, getindex(out, 2))
# return iterate(uw, args...)
# return uwreal([iterate(value(uw), 1), err(uw)], " " )
# end
# function Base.iterate(uw::uwreal, i::Int64)
# uwerr(uw)
# return iterate(value(uw),i)
# end
function Base.iterate(uw::uwreal)
return iterate(uw, 1)
end
function Base.iterate(uw::uwreal, state::Int64)
if state > length(uw)
return nothing
else
return uw[state], state +1
end
end
function Base.iterate(uw::Vector{uwreal})
return iterate(uw, 1)
end
function Base.iterate(uw::Vector{uwreal}, state::Int64)
if state > length(uw)
return nothing
else
return uw[state], state +1
end
end
function Base.abs2(uw::uwreal)
return (uw^2)^0.5
end
function Base.getindex(uw::uwreal, ii::Int64)
idx = getindex(value(uw), ii)
return uw # uwreal([getindex(value(uw), kwargs...), err(uw)], " ")
end
# function Base.getindex(uw::Vector{uwreal}, ii::Int64)
# idx = getindex(value.(uw), ii)
# return uw[ii] # uwreal([getindex(value(uw), kwargs...), err(uw)], " ")
# end
aaa = iterate(uwreal([5., 0.05], "c"))
aaa = iterate([[uwreal([5., 0.05], "c")]])
\ No newline at end of file
src/gevp analysis/gevp_fit.jl
View file @ 9fdef14c
This diff is collapsed.
src/gevp analysis/gevp_plotter.jl
View file @ 9fdef14c
This diff is collapsed.
src/gevp analysis/gevp_tm.jl
View file @ 9fdef14c
......@@ -12,6 +12,7 @@
# The analysis is performed by solving a GEVP
##################################
using Base: String
using Base: @kwdef
using LaTeXStrings: length
using OrderedCollections
using juobs, BDIO, DelimitedFiles, ADerrors, LaTeXStrings, PyPlot
......@@ -46,14 +47,14 @@ const path_bdio_tm_shifted = "/Users/alessandroconigli/Lattice/data/bdio_charm/t
const path_plot = "/Users/alessandroconigli/MyDrive/phd/analysis/charm_gevp_2022/plots"
# ensembles to analyse
const ens_list = ["H102r002"]
const ens_list = ["D200"]
# reweighting and mass shift flags
const rwf = true
const mass_shift = true
const tau = 3 # gevp shift parameter
const tau = 3 # gevp shift parameter
const TSM = false # set whether TSM is used or not
@warning("\nTSM FLAG MODE: FALSE\n")
#@warning("\nTSM FLAG MODE: $(TSM)\n")
#=============== INCLUDES ===============#
......@@ -81,9 +82,12 @@ wpmm["H400"] = [-1.0, 4.0, -1.0, 14.0*3.659]
wpmm["N202"] = [-1.0, 4.0, -1.0, 14.0*5.164]
wpmm["N200"] = [-1.0, 4.0, -1.0, 14.0*5.164]
wpmm["N203"] = [-1.0, 4.0, -1.0, 14.0*5.164]
wpmm["D200"] = [-1.0, 4.0, -1.0, 14.0*5.164]
#wpmm["D200"] = [-1.0, 4.0, -1.0, 14.0*5.164]
wpmm["D200"] = [-1.0, 4.0, -1.0, -1.0]
wpmm["N300"] = [-1.0, 4.0, -1.0, 14.0*8.595]
wpmm["J303"] = [-1.0, 4.0, -1.0, 14.0*8.595]
wpmm["D200_corr"] = [1.0, -1.0, -1.0, -1.0]
##
#=============== ANALYSIS ===============#
......@@ -100,16 +104,18 @@ end
for (k,ens) in enumerate(ensinfo)
# t0
t0_ens, YW, W = get_t0(path_data_w, ens, rw=rwf, path_rw=path_rw)
#t0_ens = t0(ens.beta)
# tm correlators
pp_tm, pp_tmW, _ = get_corr(path_data_tm, ens,"G5", "G5", rw=rwf, path_rw=path_rw, tsm=TSM)
a1a1_tm, _, _ = get_corr(path_data_tm, ens, "G1G5", "G1G5", rw=rwf, path_rw=path_rw, tsm=TSM)
a2a2_tm, _, _ = get_corr(path_data_tm, ens, "G2G5", "G2G5", rw=rwf, path_rw=path_rw, tsm=TSM)
a3a3_tm, _, _ = get_corr(path_data_tm, ens, "G3G5", "G3G5", rw=rwf, path_rw=path_rw, tsm=TSM)
pp_tm, pp_tmW = get_corr(path_data_tm, ens,"G5", "G5", rw=rwf, path_rw=path_rw, tsm=TSM)
pa1_tm, _, _ = get_corr(path_data_tm, ens, "G5", "G1G5", rw=rwf, path_rw=path_rw, tsm=TSM)
pa2_tm, _, _ = get_corr(path_data_tm, ens, "G5", "G2G5", rw=rwf, path_rw=path_rw, tsm=TSM)
pa3_tm, _, _ = get_corr(path_data_tm, ens, "G5", "G3G5", rw=rwf, path_rw=path_rw, tsm=TSM)
a1a1_tm, _ = get_corr(path_data_tm, ens, "G1G5", "G1G5", rw=rwf, path_rw=path_rw, tsm=TSM)
a2a2_tm, _ = get_corr(path_data_tm, ens, "G2G5", "G2G5", rw=rwf, path_rw=path_rw, tsm=TSM)
a3a3_tm, _ = get_corr(path_data_tm, ens, "G3G5", "G3G5", rw=rwf, path_rw=path_rw, tsm=TSM)
pa1_tm, _ = get_corr(path_data_tm, ens, "G5", "G1G5", rw=rwf, path_rw=path_rw, tsm=TSM)
pa2_tm, _ = get_corr(path_data_tm, ens, "G5", "G2G5", rw=rwf, path_rw=path_rw, tsm=TSM)
pa3_tm, _ = get_corr(path_data_tm, ens, "G5", "G3G5", rw=rwf, path_rw=path_rw, tsm=TSM)
mu_list = gen_mulist[k] = getfield.(pp_tm, :mu)
......@@ -119,15 +125,15 @@ for (k,ens) in enumerate(ensinfo)
mat_tm2 = comp_mat_multigamma(ens, pp_tm, pa2_tm, a2a2_tm)
mat_tm3 = comp_mat_multigamma(ens, pp_tm, pa3_tm, a3a3_tm)
# gevp matrices for ps and vec decays
mat_dec_ps = comp_mat(ens, pp_tm, tau)
mat_dec_vec1 = comp_mat(ens, a1a1_tm, tau)
mat_dec_vec2 = comp_mat(ens, a2a2_tm, tau)
mat_dec_vec3 = comp_mat(ens, a3a3_tm, tau)
mat_dec_ps = comp_mat_tau_shift(ens, pp_tm, tau)
mat_dec_vec1 = comp_mat_tau_shift(ens, a1a1_tm, tau)
mat_dec_vec2 = comp_mat_tau_shift(ens, a2a2_tm, tau)
mat_dec_vec3 = comp_mat_tau_shift(ens, a3a3_tm, tau)
# gevp meson masses
m_tm_ps, m_tm_vec1 = gevp_mass(mat_tm1, path_plat_tm_mps, path_plat_tm_mvec, t0=3, pl=true)
_, m_tm_vec2 = gevp_mass(mat_tm2, path_plat_tm_mps, path_plat_tm_mvec, t0=3, pl=false)
_, m_tm_vec3 = gevp_mass(mat_tm3, path_plat_tm_mps, path_plat_tm_mvec, t0=3, pl=false)
m_tm_ps, m_tm_vec1 = gevp_mass(mat_tm1, path_plat_tm_mps, path_plat_tm_mvec, t0=3, wpm=wpmm, pl=true)
_, m_tm_vec2 = gevp_mass(mat_tm2, path_plat_tm_mps, path_plat_tm_mvec, t0=3, wpm=wpmm, pl=false)
_, m_tm_vec3 = gevp_mass(mat_tm3, path_plat_tm_mps, path_plat_tm_mvec, t0=3, wpm=wpmm, pl=false)
m_tm_vec = (m_tm_vec1 .+ m_tm_vec2 .+ m_tm_vec3 ) ./ 3.
# pion mass
pion_plat = select_plateau(ens, mu_list, path_plat_tm_mps)[1]
......@@ -143,10 +149,10 @@ for (k,ens) in enumerate(ensinfo)
# dec_vec = (dec_v1 .+ dec_v2 .+ dec_v3) ./ 3.
# decay constant from ratio
dec_ps = get_f_tm(pp_tm, m_tm_ps, ens, path_plat_tm_dps, pl=true )
dec_v1 = get_fstar_tm(a1a1_tm, m_tm_vec1, ens, path_plat_tm_dvec, pl=true)
dec_v2 = get_fstar_tm(a2a2_tm, m_tm_vec2, ens, path_plat_tm_dvec, pl=false)
dec_v3 = get_fstar_tm(a3a3_tm, m_tm_vec3, ens, path_plat_tm_dvec, pl=false)
dec_ps = get_f_tm(pp_tm, m_tm_ps, ens, path_plat_tm_dps, pl=true, wpm=wpmm )
dec_v1 = get_fstar_tm(a1a1_tm, m_tm_vec1, ens, path_plat_tm_dvec, pl=true, wpm=wpmm)
dec_v2 = get_fstar_tm(a2a2_tm, m_tm_vec2, ens, path_plat_tm_dvec, pl=false, wpm=wpmm)
dec_v3 = get_fstar_tm(a3a3_tm, m_tm_vec3, ens, path_plat_tm_dvec, pl=false, wpm=wpmm)
dec_vec = (dec_v1 .+ dec_v2 .+ dec_v3) ./ 3.
close("all")
......
src/gevp analysis/gevp_tm_BMA.jl 0 → 100644
View file @ 9fdef14c
This diff is collapsed.
src/gevp analysis/read_bdio.jl
View file @ 9fdef14c
......@@ -42,7 +42,8 @@ function read_BDIO(path::String, type::String, obs::String)
"flh_star" => 11,
"fsh_star" => 12,
"fhh_star" => 13,
"muh" => 14
"muh" => 14,
"fpik" => 15
)
dict_md = Dict(
"deltam" => 0,
......@@ -50,8 +51,15 @@ function read_BDIO(path::String, type::String, obs::String)
"phi2_shifted" => 2,
"phi4_shifted" => 3
)
dict2dict = Dict("w" => dict_w, "tm" => dict_tm, "md" => dict_md)
dict_ligth_obs = Dict(
"t0" => 1,
"mpi" => 2,
"mk" => 3,
"fpi" => 4,
"fk" => 5,
"dm" => 15
)
dict2dict = Dict("w" => dict_w, "tm" => dict_tm, "md" => dict_md, "ll_obs" => dict_ligth_obs)
if !(type in keys(dict2dict))
error("Incorrect type.\ntype = $(keys(dict2dict))")
......
src/gevp analysis/tools.jl
View file @ 9fdef14c
This diff is collapsed.
src/gevp analysis/types.jl
View file @ 9fdef14c
mutable struct EnsInfo
@kwdef mutable struct EnsInfo
id::String
L::Int64
beta::Float64
deg::Bool
mpi::Float64
dtr::Int64
idm_corr::Union{Int64, Nothing}
trunc::Union{Int64, Vector{Int64}, Nothing}
function EnsInfo(ens_id::String, info::Vector{Float64})
id = ens_id
L = info[1]
beta = info[2]
deg = info[3]
mpi = info[4]
dtr = info[5]
return new(id, L, beta, deg, mpi, dtr, nothing)
function EnsInfo(ens_id::String, info)
id = ens_id
L = info[1]
beta = info[2]
deg = info[3]
mpi = info[4]
dtr = info[5]
idm_corr = info[6]
return new(id, L, beta, deg, mpi, dtr, idm_corr, nothing)
end
function EnsInfo(ens_id::String, info::Vector{Float64}, trunc::Union{Nothing, Int64, Vector{Int64}})
id = ens_id
L = info[1]
beta = info[2]
deg = info[3]
mpi = info[4]
dtr = info[5]
return new(id, L, beta, deg, mpi, dtr, trunc)
function EnsInfo(ens_id::String, info, trunc::Union{Nothing, Int64, Vector{Int64}},)
id = ens_id
L = info[1]
beta = info[2]
deg = info[3]
mpi = info[4]
dtr = info[5]
idm_corr = info[6]
return new(id, L, beta, deg, mpi, dtr, idm_corr, trunc)
end
end
mutable struct EnsObs
......
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