Skip to content

J
juobs
Commit 9308ef32 authored by Javier's avatar Javier
Browse files
Options

dec_const added + juobs_obs wrappers updated

parent 64259011
Hide whitespace changes
Inline Side-by-side
Showing with 140 additions and 22 deletions
src/juobs.jl
View file @ 9308ef32
......@@ -11,6 +11,6 @@ include("juobs_obs.jl")
export read_mesons, read_ms1, read_ms, read_md, truncate_data!
export get_matrix, energies, uwdot, uweigvals, uweigvecs, uweigen, invert, getall_eigvals, getall_eigvecs
export corr_obs, md_sea, md_val, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine
export meff, mpcac, dec_const_pcvc, comp_t0
export meff, mpcac, dec_const, dec_const_pcvc, comp_t0
end # module
src/juobs_obs.jl
View file @ 9308ef32
@doc raw"""
meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bool=false )
meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
meff(corr::Corr, plat::Vector{Int64}; pl::Bool=true, data::Bool=false)
meff(corr::Corr, plat::Vector{Int64}; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
Computes effective mass for a given correlator corr at a given plateau `plat`.
Correlator can be passed as an `Corr` struct or `Vector{uwreal}`.
......@@ -14,7 +14,8 @@ corr_pp = corr_obs.(data)
m = meff(corr_pp[1], [50, 60], pl=false)
```
"""
function meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bool=false, mu::Union{Vector{Float64}, Nothing}=nothing,
function meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bool=false,
mu::Union{Vector{Float64}, Nothing}=nothing, kappa::Union{Vector{Float64}, Nothing}=nothing,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
dim = length(corr)
aux = 0.5 .* log.((corr[2:dim-2] ./ corr[3:dim-1]).^2)
......@@ -33,6 +34,9 @@ function meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bo
ylabel(L"$m_\mathrm{eff}$")
xlabel(L"$x_0$")
if !isnothing(kappa)
title(string(L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
end
if !isnothing(mu)
title(string(L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2]))
end
......@@ -44,18 +48,24 @@ function meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bo
return (mass, aux)
end
end
meff(corr::Corr, plat::Vector{Int64}; pl::Bool=true, data::Bool=false) =
meff(corr.obs, plat, pl=pl, data=data, mu=corr.mu)
function meff(corr::Corr, plat::Vector{Int64}; pl::Bool=true, data::Bool=false,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
if corr.mu == [0.0, 0.0]
return meff(corr.obs, plat, pl=pl, data=data, kappa=corr.kappa, mu=nothing, wpm=wpm)
else
return meff(corr.obs, plat, pl=pl, data=data, mu=corr.mu, kappa=nothing, wpm=wpm)
end
end
@doc raw"""
mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false, mu::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
mpcac(a0p::Corr, pp::Corr, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false)
mpcac(a0p::Corr, pp::Corr, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
Computes PCAC mass for a given correlator `a0p` and `pp` at a given plateau `plat`.
Computes the bare 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` and `data` allow to show the plots and return data as an extra result. The `ca` allows to compute `mpcac` using
The flags `pl` and `data` allow to show the plots and return data as an extra result. The `ca` variable allows to compute `mpcac` using
the improved axial current.
```@example
......@@ -74,7 +84,9 @@ 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,
mu::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
kappa::Union{Vector{Float64}, Nothing}=nothing, mu::Union{Vector{Float64}, Nothing}=nothing,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
corr_a0p = -a0p[2:end-1]
corr_pp = pp[2:end-1]
......@@ -100,6 +112,9 @@ function mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca:
ylabel(L"$m_\mathrm{PCAC}$")
xlabel(L"$x_0$")
if !isnothing(kappa)
title(string(L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
end
if !isnothing(mu)
title(string(L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2]))
end
......@@ -111,13 +126,112 @@ function mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca:
return (mass, aux)
end
end
mpcac(a0p::Corr, pp::Corr, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false) =
mpcac(a0p.obs, pp.obs, plat, ca=ca, pl=pl, data=data, mu=a0p.mu)
## Decay constants
function mpcac(a0p::Corr, pp::Corr, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
if a0p.kappa == pp.kappa || a0p.mu == pp.mu
if a0p.mu == [0.0, 0.0]
return mpcac(a0p.obs, pp.obs, plat, ca=ca, pl=pl, data=data, kappa=a0p.kappa, mu=nothing, wpm=wpm)
else
return mpcac(a0p.obs, pp.obs, plat, ca=ca, pl=pl, data=data, mu=a0p.mu, kappa=nothing, wpm=wpm)
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)
dec_const(a0p::Corr, pp::Corr, plat::Vector{Int64}, m::uwreal; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
Computes the bare decay constant using ``A_0P`` and ``PP`` correlators . The decay constant is computed in the plateau `plat`.
Correlator can be passed as an `Corr` struct or `Vector{uwreal}`. If it is passed as a uwreal vector, effective mass `m` and source position `y0`
must be specified.
The flags `pl` and `data` allow to show the plots and return data as an extra result. The `ca` variable allows to compute `dec_const` using
the improved axial current.
**The method assumes that the source is close to the boundary.** It takes the following ratio to cancel boundary effects.
``R = \frac{f_A(x_0, y_0)}{\sqrt{f_P(T-y_0, y_0)}} * e^{m (x_0 - T/2)}``
```@example
data_pp = read_mesons(path, "G5", "G5", legacy=true)
data_a0p = read_mesons(path, "G5", "G0G5", legacy=true)
corr_pp = corr_obs.(data_pp, L=32)
corr_a0p = corr_obs.(data_a0p, L=32)
m = meff(corr_pp[1], [50, 60], pl=false)
beta = 3.46
p0 = 9.2056
p1 = -13.9847
g2 = 6 / beta
ca = -0.006033 * g2 *( 1 + exp(p0 + p1/g2))
f = dec_const(corr_a0p[1], corr_pp[1], [50, 60], m, pl=true, ca=ca)
```
"""
function 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,
kappa::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
corr_a0p = -a0p
corr_pp = pp
T = length(corr_a0p)
if ca != 0.0
der_pp = (corr_pp[3:end] - corr_pp[1:end-2]) / 2
corr_a0p = corr_a0p[2:end-1] + ca * der_pp
aux = exp.((collect(1:T-2) .- (T-1)/2) .* [m for k = 1:T-2])
else
aux = exp.((collect(0:T-1) .- T/2) .* [m for k = 1:T])
end
R = corr_a0p .* aux ./ [sqrt(corr_pp[T-y0]) for k = 1:length(corr_a0p)]
R_av = plat_av(R, plat, wpm)
f = sqrt(2) * sqrt(R_av^2) / sqrt(m)
if pl == true
isnothing(wpm) ? uwerr(R_av) : uwerr(R_av, wpm)
isnothing(wpm) ? uwerr(f) : uwerr(f, wpm)
x = 1:length(R)
y = value.(R)
dy = err.(R)
v = value(R_av)
e = err(R_av)
figure()
lbl = string(L"$af = $", sprint(show, f))
fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75, label=L"$R$")
errorbar(x, y, dy, fmt="x", color="black", label=lbl)
legend()
ylabel(L"$R_\mathrm{av}$")
xlabel(L"$x_0$")
if !isnothing(kappa)
title(string(L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
end
display(gcf())
end
if !data
return f
else
return (f, R)
end
end
function dec_const(a0p::Corr, pp::Corr, plat::Vector{Int64}, m::uwreal; ca::Float64=0.0, pl::Bool=true, data::Bool=false,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
if (a0p.y0 == pp.y0) && (a0p.kappa == pp.kappa)
return dec_const(a0p.obs, pp.obs, plat, m, a0p.y0, ca=ca, kappa=a0p.kappa, pl=pl, data=data, wpm=wpm)
else
error("y0 or kappa values does not match")
end
end
@doc raw"""
dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false)meff(corr::Corr, plat::Vector{Int64}; pl::Bool=true, data::Bool=false)
dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false)
dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
Computes decay constant using the PCVC relation for twisted mass fermions. The decay constant is computed in the plateau `plat`.
Correlator can be passed as an `Corr` struct or `Vector{uwreal}`. If it is passed as a uwreal vector, vector of twisted masses `mu` and source position `y0`
......@@ -125,9 +239,10 @@ must be specified.
The flags `pl` and `data` allow to show the plots and return data as an extra result.
**The method assumes that the source is in the bulk.**
```@example
data = read_mesons(path, "G5", "G5")
corr_pp = corr_obs.(data)
corr_pp = corr_obs.(data, L=32)
m = meff(corr_pp[1], [50, 60], pl=false)
f = dec_const_pcvc(corr_pp[1], [50, 60], m, pl=false)
```
......@@ -171,8 +286,11 @@ function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu
return (f, uwdot(sqrt(2) * (mu[1] + mu[2]) / m^1.5, R))
end
end
dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false) =
dec_const_pcvc(corr.obs, plat, m, corr.mu, corr.y0, pl=pl, data=data)
function dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
return dec_const_pcvc(corr.obs, plat, m, corr.mu, corr.y0, pl=pl, data=data, wpm=wpm)
end
#t0
function get_model(x, p, n)
......@@ -183,9 +301,9 @@ function get_model(x, p, n)
return s
end
@doc raw"""
comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Matrix{Float64}, Nothing}=nothing, npol::Int64=2)
comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Matrix{Float64}, Nothing}=nothing, npol::Int64=2, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Vector{Matrix{Float64}}, Nothing}=nothing, npol::Int64=2)
comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Vector{Matrix{Float64}}, Nothing}=nothing, npol::Int64=2, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
Computes `t0` using the energy density of the action `Ysl`(Yang-Mills action).
`t0` is computed in the plateau `plat`.
......
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