Skip to content

J
juobs
Commit 934ec290 authored by Antonino D'Anna's avatar Antonino D'Anna
Browse files
Options

Reverted previouse changes, added read_mesons_chunks and mpcac_ren

parent fa605a1b
Hide whitespace changes
Inline Side-by-side
Showing with 163 additions and 27 deletions
src/juobs.jl
View file @ 934ec290
......@@ -9,9 +9,9 @@ include("juobs_reader.jl")
include("juobs_tools.jl")
include("juobs_obs.jl")
export read_mesons, read_mesons_correction, read_ms1, read_ms, read_md, truncate_data!, concat_data!
export read_mesons, read_mesons_correction, read_ms1, read_ms, read_md, truncate_data!
export get_matrix, energies, uwdot, uweigvals, uweigvecs, uweigen, invert, getall_eigvals, getall_eigvecs, hess_reduce, uwcholesky, transpose, tridiag_reduction, make_positive_def, invert_covar_matrix
export corr_obs, corr_obs_TSM, corr_sym, md_sea, md_val, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine, global_fit_routine, bayesian_av
export meff, mpcac, dec_const, dec_const_w, dec_const_pcvc, comp_t0
export meff, mpcac, mpcac_ren, dec_const, dec_const_w, dec_const_pcvc, comp_t0
end # module
src/juobs_obs.jl
View file @ 934ec290
......@@ -92,7 +92,7 @@ m12 = mpcac(corr_a0p, corr_pp, [50, 60], pl=false, ca=ca)
"""
function mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false,
kappa::Union{Vector{Float64}, Nothing}=nothing, mu::Union{Vector{Float64}, Nothing}=nothing,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, savepl::Bool = false)
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
corr_a0p = -a0p[1:end]
corr_pp = pp[1:end]
......@@ -129,21 +129,15 @@ function mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca:
_, max_idx = findmax(value.(pp))
ylim(v-20*e, v+20*e)
xlim(left=max_idx)
filename=""
if !isnothing(kappa)
title(string(L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
global filename = string("mpcac/kappa_1_$(kappa[1])_kappa_2_$(kappa[2]).png")
end
if !isnothing(mu)
title(string(L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2]))
global filename =string("mpcac/mu_1_$(mu[1])_mu_2_$(mu[2]).png")
end
if savepl
savefig(filename)
else
display(gcf())
end
end
display(gcf())
end
if !data
return mass
else
......@@ -164,6 +158,121 @@ function mpcac(a0p::Corr, pp::Corr, plat::Vector{Int64}; ca::Float64=0.0, pl::Bo
error("mu or kappa values does not match")
end
end
@doc raw"""
mpcac_ren(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}, za_zp::uwreal; ca::Float64=0.0,ba_bp_z::uwreal=0.0,pl::Bool=true, data::Bool=false,
kappa::Union{Vector{Float64}, Nothing}=nothing, mu::Union{Vector{Float64}, Nothing}=nothing,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, savepl::Bool =false)
mpcac_ren(a0p::Corr, pp::Corr, plat::Vector{Int64}, za_zp::uwreal; ca::Float64,ba_bp_z::uwreal, pl::Bool=true, data::Bool=false,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,savepl::Bool =false)
Computes the renormalized PCAC mass for a given correlator `a0p` and `pp` at a given plateau `plat`.
Correlator can be passed as an `Corr` struct or `Vector{uwreal}`.
The flags `pl`, `savepl` and `data` allow to show the plots, save it as a png file and return data as an extra result. The `ca` variable allows to compute `mpcac` using
the improved axial current. The `ba_bp_z` variable allows to compute `mpcac` using the renormalized improved axial and pseudoscalar currents.
Without specifying `ba_bp_z` the results is equivalent to the output of `mpcac(a0p,pp,plat, ca=ca)` multiplied by `za_zp`
```@example
include("juobs_const.jl")
data_pp = read_mesons(path, "G5", "G5")
data_a0p = read_mesons(path, "G5", "G0G5")
corr_pp = corr_obs.(data_pp)
corr_a0p = corr_obs.(data_a0p)
za_zp = za(beta)/ZP(beta)
m12 = mpcac_ren(corr_a0p, corr_pp, [50, 60], za_zp, pl=false)
#equivalent to m12 = za_zp.*mpcac(corr_a0p,corr_pp,[50,60],pl=false)
p0 = 9.2056
p1 = -13.9847
g2 = 1.73410
ca = -0.006033 * g2 *( 1 + exp(p0 + p1/g2))
m12 = mpcac_ren(corr_a0p, corr_pp, [50, 60],za_zp, pl=false, ca=ca)
#equivalent to m12 = za_zp.*mpcac(corr_a0p,corr_pp,[50,60],pl=false,ca=ca)
ba_bp_z = ba_bp(beta)*Z(beta)
m12 = mpcac_ren = (corr_a0p,corr_pp,[50,60],za_zp,pl=false, ca=ca, ba_bp_z=ba_bp_z)
```
"""
function mpcac_ren(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}, za_zp::uwreal; ca::Float64=0.0,ba_bp_z::uwreal=0.0,
pl::Bool=true, data::Bool=false,kappa::Union{Vector{Float64}, Nothing}=nothing, mu::Union{Vector{Float64}, Nothing}=nothing,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,savepl::Bool =false)
corr_a0p = -a0p[1:end]
corr_pp = pp[1:end]
der_a0p = (corr_a0p[4:end-1] .- corr_a0p[2:end-3]) / 2 # \de_0 a0p
m_pcac = der_a0p ./(2 * corr_pp[3:end-2]) # Unimproved, unrenormalized mpcac
if ba_bp_z !=0.0
der_a0p *= 1.0 .+ ba_bp_z*m_pcac #improvement with ba_bp
end
if ca !=0.0
der2_pp = (corr_pp[1:end-4] .- 2*corr_pp[3:end-2] .+ corr_pp[5:end])/4
der_a0p += ca*der2_pp
end
m_ren = za_zp.*(der_a0p./(2*corr_pp[3:end-2]))
mass = plat_av(m_ren,plat,wpm)
if pl || savepl
isnothing(wpm) ? uwerr(mass) : uwerr(mass, wpm)
x = 1:length(m_ren)
y = value.(m_ren)
dy = err.(m_ren)
v = value(mass)
e = err(mass)
plot= figure()
fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75)
errorbar(x, y, dy, fmt="x", color="black")
ylabel(L"$m^R_\mathrm{PCAC}$")
xlabel(L"$x_0$")
_, max_idx = findmax(value.(pp))
# ylim(v-20*e, v+20*e)
xlim(left=max_idx)
filename=""
if !isnothing(kappa)
title(string(L"$m_\mathrm{PCAC}^R$ $\mathcal{O(a)} improved","\n",L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
global filename = string("mpcac/kappa_1_$(kappa[1])_kappa_2_$(kappa[2])_renormalize_improved.png")
end
if !isnothing(mu)
title(string(L"$m_\mathrm{PCAC}^R$ $\mathcal{O(a)} improved","\n",L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2]))
global filename =string("mpcac/mu_1_$(mu[1])_mu_2_$(mu[2])_renomalized_improved.png")
end
if savepl
plot.tight_layout()
savefig(filename)
end
if pl
display(gcf())
else
close(plot)
end
end
if !data
return mass
end
return (mass, m_ren)
end
function mpcac_ren(a0p::Corr, pp::Corr, plat::Vector{Int64}, za_zp::uwreal; ca::Float64,ba_bp_z::uwreal, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,savepl::Bool =false)
if a0p.kappa == pp.kappa || a0p.mu == pp.mu
if a0p.mu == [0.0, 0.0]
return mpcac_ren(a0p.obs, pp.obs, plat, za_zp, ca=ca,ba_bp_z=ba_cp_z, pl=pl, data=data, kappa=a0p.kappa, mu=nothing, wpm=wpm, savepl=savepl)
else
return mpcac_ren(a0p.obs, pp.obs, plat, za_zp, ca=ca,ba_bp_z=ba_cp_z, pl=pl, data=data, mu=a0p.mu, kappa=nothing, wpm=wpm, savepl=savepl)
end
else
error("mu or kappa values does not match")
end
end
## Decay constants
@doc raw"""
dec_const(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
......
src/juobs_reader.jl
View file @ 934ec290
......@@ -284,6 +284,37 @@ function read_mesons_correction(path::Vector{String}, g1::Union{String, Nothing}
return cdata
end
@doc raw"""
read_mesons_chunks(paths::Vector{String}, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{String, Nothing}=nothing, legacy::Bool=false,
nnoise_trunc::Union{Int64, Nothing}=nothing)
This function reads mesons dat files stored in chunks and returns a single vector of `CData` structures for different masses and Dirac structures.
the paths to each chunks is given through the variable `paths::Vector{String}`
Dirac structures `g1` and/or `g2` can be passed as string arguments in order to filter correlators.
ADerrors id can be specified as argument. If is not specified, the `id` is fixed according to the ensemble name (example: "H400"-> id = "H400")
*For the old version (without smearing, distance preconditioning and theta) set legacy=true.
Examples:
```@example
read_mesons([path_chunk1,path_chunk2,path_chunk3])
read_mesons([path_chunk1,path_chunk2,path_chunk3], "G5")
read_mesons([path_chunk1,path_chunk2,path_chunk3], nothing, "G5")
read_mesons([path_chunk1,path_chunk2,path_chunk3], "G5", "G5")
read_mesons([path_chunk1,path_chunk2,path_chunk3], "G5", "G5", id="H100")
read_mesons([path_chunk1,path_chunk2,path_chunk3], "G5_d2", "G5_d2", legacy=true)
```
"""
function read_mesons_chunks(paths::Vector{String}, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{String, Nothing}=nothing, legacy::Bool=false,
nnoise_trunc::Union{Int64, Nothing}=nothing)
data = read_mesons(paths[1],g1,g2,id=id,legacy=legacy, nnoise_trunc=nnoise_trunc)
for i in 2:length(paths)
temp = read_mesons(paths[i],g1,g2 ,id=id,legacy=legacy, nnoise_trunc=nnoise_trunc)
concat_data!(data,temp)
end
return data
end
function read_rw(path::String; v::String="1.2")
data = open(path, "r")
nrw = read(data, Int32)
......
src/juobs_tools.jl
View file @ 934ec290
......@@ -1095,7 +1095,7 @@ fit_routine(model, xdata, ydata, param=3, covar=true)
```
"""
function fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}, param::Int64=3; info::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,
inv_cov::Union{Matrix{Float64}, Nothing}=nothing, verbose::Bool = true)
inv_cov::Union{Matrix{Float64}, Nothing}=nothing)
isnothing(wpm) ? uwerr.(ydata) : [uwerr(yaux, wpm) for yaux in ydata]
......@@ -1125,30 +1125,26 @@ function fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}
# println((fit.resid[i])^2)
# end
# println("\n")
if verbose
println("chi2 from fit residual = ", chi2_fit_res)
println("chi2 from chi2 function = ", chisq(coef(fit), ydata))
end
println("chi2 from fit residual = ", chi2_fit_res)
println("chi2 from chi2 function = ", chisq(coef(fit), ydata))
#Info
if info
for i = 1:length(upar)
isnothing(wpm) ? uwerr(upar[i]) : uwerr(upar[i], wpm)
if verbose
print("\n Fit parameter: ", i, ": ")
details(upar[i])
end
print("\n Fit parameter: ", i, ": ")
details(upar[i])
end
end
#chis2_corrected = (length(yval) - param) * chisq(coef(fit), ydata) / chi_exp
chis2_corrected = (length(yval) - param) * chi2_fit_res / chi_exp
if verbose
#println("Chisq / chiexp: ", chisq(coef(fit), ydata), " / ", chi_exp, " (dof: ", length(yval) - param,")")
println("Chisq / chiexp: ", chi2_fit_res, " / ", chi_exp, " (dof: ", length(yval) - param,")")
println("Chisq corrected: ", chis2_corrected)
println("Return: params, chi2, chi2exp")
end
#println("Chisq / chiexp: ", chisq(coef(fit), ydata), " / ", chi_exp, " (dof: ", length(yval) - param,")")
println("Chisq / chiexp: ", chi2_fit_res, " / ", chi_exp, " (dof: ", length(yval) - param,")")
println("Chisq corrected: ", chis2_corrected)
println("Return: params, chi2, chi2exp")
return upar, chi2_fit_res, chi_exp
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