@doc raw"""
meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=false, 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, savepl::Bool=false,
filename::String="",plot_kw...)
meff(corr::Corr, plat::Vector{Int64}; kwargs..)
meff(corr::Vector{Corr},plat::Vector{Vector{Int64}}; kwargs...)
Computes the mass associated to given correlator `corr` at a given plateau `plat`.
Correlator can be passed as an `Corr` struct or `Vector{uwreal}`.
In case of Vector{Corr} it is better to not broadcast, especially in case of `data=true` or `pl=true`, since the the overloaded function for Vector{Corr} takes care of the extra outputs and organize them into a `Tuple` of Vector
# Arguments
- `pl::Bool=false` and `savepl::Bool=false`: show and save the plots. The plots are made through the functions `plot_data` and `plot_func`. It is possible to pass keyword arguments for these two functions.
When `savepl = true`, the plots are saved and closed. When `pl = true` the plot are generated, displayed and the last two returns are always `fig,ax` object to allow more costumability.
If both are `true`, `pl` is ignored.
See also: [`plot_data`](@ref), [`plot_func`](@ref).
- `data::Bool=false`: if set to `true`, the second return is a vector with the effective mass.
- `filename::String=""`: if given, overrides the automatically generated filename. `mkpath(dirname(filename))` is called before saving. If not given, a folder "meff/" is made and the correlator `μ` or `κ` are used to name the figure
# Examples
- with one correlator
```@example
data = read_mesons(path, "G5", "G5")[1]
corr_pp = corr_obs(data)
m = meff(corr_pp, [50, 60]) # return only meff
```
- with more correlator, making plots and returning data vector too
```@example
data = read_mesons(path, "G5", "G5")
corr_pp = corr_obs.(data)
m,fig,ax = meff(corr_pp, [[50, 60] for eachindex(corr_pp), pl=true)
# return meff, a vector of fig obj and a vector of ax obj
m,data = meff(corr_pp[1], [50, 60], pl=true,data=true, savepl=true)
# since pl and savepl are both true, pl is ignored, the plot is generated and saved
# the second arguments is given by the data vector
```
- with more correlator, making plots and returning effective PCAC masses
```@example
data_pp = read_mesons(path, "G5", "G5")
corr_pp = corr_obs.(data_pp)
m, data,fig,ax = meff(corr_a0p, corr_pp, [[50, 60] for eachindex(corr_pp), pl=true,data=true)
#the return is a Vector{uwreal} with the pseudoscalar masses, a Vector{Vector{uwreal}} with the effective ,
#a Vector{Figure} with the figure objects and a Vector{PyPlot.PyCall.PyObject} with the axes objects
```
"""
function meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=false, 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, savepl::Bool=false,
filename::String="",plot_kw...)
dim = length(corr)
aux = 0.5 .* log.((corr[2:dim-2] ./ corr[3:dim-1]).^2)
mass = plat_av(aux, plat, wpm)
if pl || savepl
abs_mass = abs2(mass);
isnothing(wpm) ? uwerr(abs_mass) : uwerr(abs_mass, wpm)
v,e = value(abs_mass), err(abs_mass);
x = [_ for _ in eachindex(aux)]
fig,ax = plot_data(abs2.(aux),x;plot_kw...)
fig,ax = plot_func(x->abs_mass,x[plat[1]:plat[2]],figs=(fig,ax) ;plot_kw...)
if !haskey(plot_kw,:title)
temp = isnothing(kappa) ? L"$\mu_1 = %$(mu[1])$ $\mu_2 = %$(mu[2])$" : L"$\kappa_1 = %$(kappa[1])$ $\kappa_2 = %$(kappa[2])$"
fig.suptitle(temp);
end
if !haskey(plot_kw,:ylabel)
ax.set_ylabel(L"$m_\mathrm{eff}$")
end
if !haskey(plot_kw,:xlabel)
ax.set_xlabel(L"$x_0$")
end
_, max_idx = findmax(value.(corr))
ax.set_xlim(left=max_idx)
ax.set_ylim(v-20*e,v+20*e)
if savepl
fig.tight_layout();
if filename ==""
filename = "meff/"
filename *= isnothing(kappa) ? "" : "kappa_1_$(kappa[1])_kappa_2_$(kappa[2]).png"
filename *= isnothing(mu) ? "" : "mu_1_$(mu[1])_mu_2_$(mu[2]).png"
end
mkpath(dirname(filename))
fig.savefig(filename)
close(fig)
elseif pl
display(fig);
end
end
if pl && !savepl
return data ? (mass,aux,fig,ax) : (mass,fig,ax)
else
return data ? (mass,aux) : mass
end
end
function meff(corr::Corr, plat::Vector{Int64}; kwargs...)
if corr.mu == [0.0, 0.0]
return meff(corr.obs, plat; kappa=corr.kappa, mu=nothing,kwargs...)
else
return meff(corr.obs, plat; mu=corr.mu, kappa=nothing, kwargs...)
end
end
function meff(corr::Vector{Corr},plat::Vector{Vector{Int64}}; kwargs...)
a = meff.(corr,plat; kwargs...)
if length(a[1])==1
return a
end
# Due to the broadcast, a is a Vector{Tuple{...}} where each element of the vector is the output of the function
# The following lines of code reshape a into a Matrix such that each column
# correspond to one Tuple and each row correspond to one component of the Tuple. The it will return each row
# separately.
#
# Example:
# V = [(a1,a2,a3),(b1,b2,b3)]
# V = collect.(V) #[[a1,a2,a3],[b1,b2,b3]]
# V stack(V) #[a1 b1; a2 b2; a3 b3]
# return Tuple([[r...] for e in eachrow(r)]) #([a1,b1],[a1,b2],[a3,b3])
a = stack(collect.(a));
return Tuple([r...] for r in eachrow(a))
end
function meff(corr::Vector{Vector{uwreal}},plat::Vector{Vector{Int64}}; kwargs...)
a = meff.(corr,plat; kwargs...)
if length(a[1])==1
return a
end
# Due to the broadcast, a is a Vector{Tuple{...}} where each element of the vector is the output of the function
# The following lines of code reshape a into a Matrix such that each column
# correspond to one Tuple and each row correspond to one component of the Tuple. The it will return each row
# separately.
#
# Example:
# V = [(a1,a2,a3),(b1,b2,b3)]
# V = collect.(V) #[[a1,a2,a3],[b1,b2,b3]]
# V stack(V) #[a1 b1; a2 b2; a3 b3]
# return Tuple([[r...] for e in eachrow(r)]) #([a1,b1],[a1,b2],[a3,b3])
a = stack(collect.(a));
return Tuple([r...] for r in eachrow(a))
end
@doc raw"""
mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64};
ca::Float64=0.0, pl::Bool=false, data::Bool=false, kappa::Union{Vector{Float64}, Nothing}=nothing,
mu::Union{Vector{Float64}, Nothing}=nothing, savepl::Bool = false,filename::String="",
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,
plot_kw...)
mpcac(a0p::Corr, pp::Corr, plat::Vector{Int64}; kwargs...)
mpcac(a0p::Vector{Corr},pp::Vector{Corr},plat::Vector{Vector{Int64}}; kwargs...)
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}`.
In case of Vector{Corr} it is better to not broadcast, especially in case of `data=true` or `pl=true`, since the the overloaded function for Vector{Corr} takes care of the extra outputs and organize them into a `Tuple` of Vector
The return structure is: mpcac[,data][,figs,axs]
# Arguments
- `pl::Bool=false` and `savepl::Bool=false`: show and save the plots. The plots are made through the functions `plot_data` and `plot_func`. It is possible to pass keyword arguments for these two functions.
When `savepl = true`, the plots are saved and closed. When `pl = true` the plot are generated, displayed and the last two returns are always `fig,ax` object to allow more costumability
If both are `true`, `pl` is ignored.
See also: [`plot_data`](@ref), [`plot_func`](@ref).
- `data::Bool=false`: if set to `true`, the second return is a vector with the effective PCAC mass.
- `ca::Float64=0.0`: if given, computes the PCAC mass with the improved axial current.
- `filename::String=""`: if given, overrides the automatically generated filename. `mkpath(dirname(filename))` is called before saving. If not given, a folder "mpcac/" is made and the correlator `μ` or `κ` are used to name the figure
# Examples
- with one correlator
```@example
data_pp = read_mesons(path, "G5", "G5")[1] #takes only the first component
data_a0p = read_mesons(path, "G5", "G0G5")[1] #takes only the first component
corr_pp = corr_obs(data_pp)
corr_a0p = corr_obs(data_a0p)
m12 = mpcac(corr_a0p, corr_pp, [50, 60]) # return is only mpcac
```
- with more correlator, `ca` improved, making plots and returning effective PCAC masses
```@example
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)
m12,fig,ax = mpcac(corr_a0p, corr_pp, [[50, 60] for eachindex(corr_pp), pl=true,)
# return mpcac, a vector of fig obj and a vector of ax obj
p0 = 9.2056
p1 = -13.9847
g2 = 1.73410
ca = -0.006033 * g2 *( 1 + exp(p0 + p1/g2))
m12,data = mpcac(corr_a0p[1], corr_pp[1], [50, 60], pl=true, ca=ca,data=true, savepl=true, filename = "mpcac/with_ca_improvement.png")
# since pl and savepl are both true, pl is ignored, the plot is generated and saved in mpcac/with_ca_improvement.png
# the second arguments is given by
```
- with more correlator, making plots and returning effective PCAC masses
```@example
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)
m12, data,fig,ax = mpcac(corr_a0p, corr_pp, [[50, 60] for eachindex(corr_pp), pl=true,data=true)
#the return is a Vector{uwreal} with the PCAC masses, a Vector{Vector{uwreal}} with the effective PCAC,
#a Vector{Figure} with the figure objects and a Vector{PyPlot.PyCall.PyObject} with the axes objects
```
"""
function mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=false, 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,filename::String="",plot_kw...)
corr_a0p = -a0p[1:end]
corr_pp = pp[1:end]
if ca != 0.0
der_a0p = (corr_a0p[3:end] .- corr_a0p[1:end-2]) / 2
der2_pp = (corr_pp[1:end-4] - 2*corr_pp[3:end-2] + corr_pp[5:end])/4
der_a0p = der_a0p[2:end-1] + ca * der2_pp
else
der_a0p = (corr_a0p[4:end-1] .- corr_a0p[2:end-3]) / 2
end
aux = der_a0p ./ ( 2. .* corr_pp[3:end-2])
mass = plat_av(aux, plat, wpm)
if pl|| savepl
x = [_ for _ in eachindex(aux)]
fig,ax = plot_data(aux,x,wpm=wpm,ylabel =L"$m_\mathrm{PCAC}$",xlabel=L"$x_0$";plot_kw...)
isnothing(wpm) ? uwerr(mass) : uwerr(mass,wpm)
v,e = value(mass), err(mass);
fig,ax = plot_func(x->mass,[plat[1]:plat[2]...],figs=(fig,ax);plot_kw...)
if !haskey(plot_kw,:title)
temp = isnothing(kappa) ? L"$\mu_1 = %$(mu[1])$ $\mu_2 = %$(mu[2])$" : L"$\kappa_1 = %$(kappa[1])$ $\kappa_2 = %$(kappa[2])$"
fig.suptitle(temp);
end
if !haskey(plot_kw,:ylabel)
ax.set_ylabel(L"$m_\mathrm{PCAC}$")
end
if !haskey(plot_kw,:xlabel)
ax.set_xlabel(L"$x_0$")
end
_, max_idx = findmax(value.(pp))
ax.set_xlim(left=max_idx)
ax.set_ylim(v-20*e,v+20*e)
if savepl
fig.tight_layout();
if filename ==""
filename = "mpcac/"
filename *= isnothing(kappa) ? "" : "kappa_1_$(kappa[1])_kappa_2_$(kappa[2]).png"
filename *= isnothing(mu) ? "" : "mu_1_$(mu[1])_mu_2_$(mu[2]).png"
end
mkpath(dirname(filename))
fig.savefig(filename)
close(fig)
elseif pl
display(fig);
end
end
if pl && !savepl
return data ? (mass, aux, fig, ax) : (mass, fig, ax)
else
return data ? (mass, aux) : mass
end
end
function mpcac(a0p::Corr, pp::Corr, plat::Vector{Int64}; kwargs... )
if a0p.kappa == pp.kappa || a0p.mu == pp.mu
if a0p.mu == [0.0, 0.0]
return mpcac(a0p.obs, pp.obs, plat; kappa=a0p.kappa, mu=nothing, kwargs...)
else
return mpcac(a0p.obs, pp.obs, plat; mu=a0p.mu, kappa=nothing, kwargs...)
end
else
error("mu or kappa values does not match")
end
end
function mpcac(a0p::Vector{Vector{uwreal}},pp::Vector{Vector{uwreal}},plat::Vector{Vector{Int64}}; kwargs...)
a = mpcac.(a0p,pp,plat; kwargs...)
if length(a[1])==1
return a
end
# Due to the broadcast, a is a Vector{Tuple{...}} where each element of the vector is the output of the function
# with one set of a0p, pp and plat. The following lines of code reshape a into a Matrix such that each column
# correspond to one Tuple and each row correspond to one component of the Tuple. The it will return each row
# separately.
#
# Example:
# V = [(a1,a2,a3),(b1,b2,b3)]
# V = collect.(V) #[[a1,a2,a3],[b1,b2,b3]]
# V stack(V) #[a1 b1; a2 b2; a3 b3]
# return Tuple([[r...] for e in eachrow(r)]) #([a1,b1],[a1,b2],[a3,b3])
a = stack(collect.(a));
return Tuple([r...] for r in eachrow(a))
end
function mpcac(a0p::Vector{Corr},pp::Vector{Corr},plat::Vector{Vector{Int64}}; kwargs...)
a = mpcac.(a0p,pp,plat; kwargs...)
if length(a[1])==1
return a
end
# Due to the broadcast, a is a Vector{Tuple{...}} where each element of the vector is the output of the function
# with one set of a0p, pp and plat. The following lines of code reshape a into a Matrix such that each column
# correspond to one Tuple and each row correspond to one component of the Tuple. The it will return each row
# separately.
#
# Example:
# V = [(a1,a2,a3),(b1,b2,b3)]
# V = collect.(V) #[[a1,a2,a3],[b1,b2,b3]]
# V stack(V) #[a1 b1; a2 b2; a3 b3]
# return Tuple([[r...] for e in eachrow(r)]) #([a1,b1],[a1,b2],[a3,b3])
a = stack(collect.(a));
return Tuple([r...] for r in eachrow(a))
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=false, data::Bool=false,
kappa::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,savepl::Bool=false,
filename::String="",plot_kw...)
dec_const(a0p::Corr, pp::Corr, plat::Vector{Int64}, m::uwreal; kwargs...)
dec_const(a0p::Vector{Corr},pp::Vector{Corr},plt::Vector{Vector{Int64}},m::Vector{uwreal}; kwargs...)
dec_const(a0pL::Vector{uwreal}, a0pR::Vector{uwreal}, ppL::Vector{uwreal}, ppR::Vector{uwreal}, plat::Vector{Int64},
m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=false, data::Bool=false, kappa::Union{Vector{Float64}, Nothing}=nothing,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,
savepl::Bool=false, filename::String="",plot_kw...)
dec_const(a0pL::Corr, a0pR::Corr, ppL::Corr, ppR::Corr, plat::Vector{Int64}, m::uwreal; kwargs...)
dec_const(a0pL::Vector{Corr},a0pR::Vector{Corr},ppL::Vector{Corr}, ppR::Vector{Corr}, plat::Vector{Vector{Int64}}, m::Vector{uwreal}; kwargs...)
dec_const(vv::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; pl::Bool=false, data::Bool=false,
kappa::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,
savepl::Bool=false,filename::String="",plot_kw...)
dec_const(vv::Corr, plat::Vector{Int64}, m::uwreal;kwargs...)
dec_const(vv::Vector{Corr}, plat::Vector{Vector{Int64}},m::Vector{uwreal}; kwargs...)
Computes the bare decay constant using `A_0P` and `PP` or the `VV` correlators and the associated masses `m`. 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, source position `y0` must be specified.
# Arguments
- `pl::Bool=false` and `savepl::Bool=false`: show and save the plots. The plots are made through the functions `plot_data` and `plot_func`. It is possible to pass keyword arguments for these two functions.
When `savepl = true`, the plots are saved and closed. When `pl = true` the plot are generated, displayed and the last two returns are always `fig,ax` objects to allow more costumability
If both are `true`, `pl` is ignored.
- `data::Bool=false`: if set to `true`, the second return is a vector with the ratio `R`(see below).
- `ca::Float64=0.0`: (ONLY IF `A_0P`) if given, computes the decay constant mass with the improved axial current.
- `filename::String=""`: if given, overrides the automatically generated filename. `mkpath(dirname(filename))` is called before saving. If not given, a folder "dec\_cont/" or "dec\_const\_vec" is made and the correlator's `κ` is used to name the figure
**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)}``
**If left and right correlators are included in the input. The result is computed with the following ratio**
`` R = \sqrt{f_A(x_0, y_0) * f_A(x_0, T - 1 - y_0) / f_P(T - 1 - y_0, y_0)}``
**If the vectorial correlator is used, ther results is computed with the following ratio**
`` R = \sqrt{\frac{2 f_V(x_0,y_0) * e^{m*(x0-y0)}{m}}}``
```@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)
corr_a0pL, corr_a0pR = [corr_a0p[1], corr_a0p[2]]
corr_ppL, corr_ppR = [corr_pp[1], corr_pp[2]]
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)
f_ratio = dec_const(corr_a0pL, corr_a0pR, corr_ppL, corr_ppR, [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=false, data::Bool=false,
kappa::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,savepl::Bool=false,
filename::String="",plot_kw...)
corr_a0p = -a0p[2:end-1]
corr_pp = pp[2:end-1]
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
T = length(corr_a0p)
aux = exp.((collect(1:T) .- (T)/2) .* [m for k = 1:T])
else
T = length(corr_a0p)
aux = exp.((collect(0:T-1) .- (T-1)/2) .* [m for k = 1:T])
end
R = corr_a0p .* aux ./ [((corr_pp[T-y0])^2)^(1/4) for k = 1:length(corr_a0p)]
R_av = plat_av(R, plat, wpm)
f = sqrt(2) * sqrt(R_av^2) / sqrt(m)
if pl || savepl
isnothing(wpm) ? uwerr(R_av) : uwerr(R_av, wpm)
isnothing(wpm) ? uwerr(f) : uwerr(f, wpm)
x = [1:length(R)...]
lbl = string(L"$af = $", sprint(show, f))
fig,ax = plot_data(R,x,label = lbl;plot_kw...)
fig,ax = plot_func(x->R_av,x[plat[1]:plat[2]],figs=(fig,ax),label = L"$R$";plot_kw...)
v,e = value(R_av),err(R_av)
ax.set_ylim(v-10*e, v+10*e)
if !haskey(plot_kw,:ylabel)
ax.set_ylabel(L"$R_\mathrm{av}$")
end
if !haskey(plot_kw,:xlabel)
ax.set_xlabel(L"$x_0$")
end
if !haskey(plot_kw,:title)
temp = isnothing(kappa) ? "" : L"$\kappa_1 = %$(kappa[1])$ $\kappa_2 = %$(kappa[2])$"
fig.suptitle(temp);
end
if savepl
if filename==""
filename = "dec_const/"
filename *= isnothing(kappa) ? "" : "kappa_1_$(kappa[1])_kappa_2_$(kappa[2]).png"
end
fig.tight_layout();
mkpath(dirname(filename));
fig.savefig(filename);
close(fig)
elseif pl
display(fig)
end
end
if pl &&!savepl
return data ? (f,R,fig,ax) : (f,fig,ax)
else
return data ? (f, R) : f
end
end
function dec_const(a0p::Corr, pp::Corr, plat::Vector{Int64}, m::uwreal; kwargs...)
if (a0p.y0 == pp.y0) && (a0p.kappa == pp.kappa)
return dec_const(a0p.obs, pp.obs, plat, m, a0p.y0,kappa=a0p.kappa; kwargs...)
else
error("y0 or kappa values does not match")
end
end
function dec_const(a0p::Vector{Corr},pp::Vector{Corr},plt::Vector{Vector{Int64}},m::Vector{uwreal}; kwargs...)
a = dec_const.(a0p,pp,plt,m;kwargs...)
if length(a[1])==1
return a;
end
# Due to the broadcast, a is a Vector{Tuple{...}} where each element of the vector is the output of the function
# . The following lines of code reshape a into a Matrix such that each column
# correspond to one Tuple and each row correspond to one component of the Tuple. The it will return each row
# separately.
#
# Example:
# V = [(a1,a2,a3),(b1,b2,b3)]
# V = collect.(V) #[[a1,a2,a3],[b1,b2,b3]]
# V stack(V) #[a1 b1; a2 b2; a3 b3]
# return Tuple([r...] for e in eachrow(r)) #([a1,b1],[a1,b2],[a3,b3])
a = stack(collect.(a));
return Tuple([r...] for r in eachrow(a))
end
function dec_const(a0pL::Vector{uwreal}, a0pR::Vector{uwreal}, ppL::Vector{uwreal}, ppR::Vector{uwreal}, plat::Vector{Int64},
m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=false, data::Bool=false, kappa::Union{Vector{Float64}, Nothing}=nothing,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,
savepl::Bool=false, filename::String="",plot_kw...)
corr_pp = (ppL + ppR[end:-1:1]) / 2
T = length(corr_pp)
if ca != 0.0
der_ppL = (ppL[3:end] - ppL[1:end-2]) / 2
der_ppR = (ppR[3:end] - ppR[1:end-2]) / 2
corr_a0pL = -a0pL[2:end-1] + ca * der_ppL
corr_a0pR = -a0pR[2:end-1] + ca * der_ppR
else
corr_a0pL = -a0pL[2:end-1]
corr_a0pR = -a0pR[2:end-1]
end
f1 = [corr_pp[T - y0] for k = 1:length(corr_a0pL)]
R = ((corr_a0pL .* corr_a0pR ./ f1).^2).^(1/4)
R_av = plat_av(R, plat, wpm)
f = sqrt(2) * sqrt(R_av^2) / sqrt(m)
if pl || savepl
isnothing(wpm) ? uwerr(R_av) : uwerr(R_av, wpm)
isnothing(wpm) ? uwerr(f) : uwerr(f, wpm)
x = [1:length(R)...]
lbl = string(L"$af = $", sprint(show, f))
fig,ax = plot_data(R,x,label = lbl;plot_kw...)
fig,ax = plot_func(x->R_av,x[plat[1],plat[2]],figs=(fig,ax),label = L"$R$";plot_kw...)
v,e = value(R_av),err(R_av)
ax.set_ylim(v-10*e, v+10*e)
if !haskey(plot_kw,:ylabel)
ax.set_ylabel(L"$R_\mathrm{av}$")
end
if !haskey(plot_kw,:xlabel)
ax.set_xlabel(L"$x_0$")
end
if !haskey(plot_kw,:title)
temp = isnothing(kappa) ? L"$\mu_1 = %$(mu[1])$ $\mu_2 = %$(mu[2])$" : L"$\kappa_1 = %$(kappa[1])$ $\kappa_2 = %$(kappa[2])$"
fig.suptitle(temp);
end
if savepl
if filename==""
filename = "dec_const/"
filename *= isnothing(kappa) ? "" : "kappa_1_$(kappa[1])_kappa_2_$(kappa[2]).png"
end
fig.tight_layout();
mkpath(dirname(filename));
fig.savefig(filename);
close(fig)
elseif pl
display(fig)
end
end
if pl&&!savepl
return data ? (f,R,fig,ax) : (f,fig,ax)
else
return data ? (f, R) : f
end
end
function dec_const(a0pL::Corr, a0pR::Corr, ppL::Corr, ppR::Corr, plat::Vector{Int64}, m::uwreal; kwargs...)
T = length(a0pL.obs)
if (a0pL.y0 == ppL.y0) && (a0pR.y0 == ppR.y0) && (a0pL.kappa == ppL.kappa) && (a0pR.kappa == ppR.kappa) && (a0pL.y0 == T - 1 - a0pR.y0)
return dec_const(a0pL.obs, a0pR.obs, ppL.obs, ppR.obs, plat, m, a0pL.y0; kwargs...)
else
error("y0 or kappa values does not match")
end
end
function dec_const(a0pL::Vector{Corr},a0pR::Vector{Corr},ppL::Vector{Corr}, ppR::Vector{Corr}, plat::Vector{Vector{Int64}}, m::Vector{uwreal}; kwargs...)
a = dec_const.(a0pL,a0pR,ppL,ppR,plat,m;kwargs...)
if length(a[1])==1
return a;
end
# Due to the broadcast, a is a Vector{Tuple{...}} where each element of the vector is the output of the function
# The following lines of code reshape a into a Matrix such that each column
# correspond to one Tuple and each row correspond to one component of the Tuple. The it will return each row
# separately.
#
# Example:
# V = [(a1,a2,a3),(b1,b2,b3)]
# V = collect.(V) #[[a1,a2,a3],[b1,b2,b3]]
# V = stack(V) #[a1 b1; a2 b2; a3 b3]
# return Tuple([r...] for e in eachrow(r)) #([a1,b1],[a1,b2],[a3,b3])
a = stack(collect.(a));
return Tuple([r...] for r in eachrow(a))
end
function dec_const(vv::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; pl::Bool=false, data::Bool=false,
kappa::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,
savepl::Bool=false,filename::String="",plot_kw...)
corr_vv = vv[2:end-1]
T = length(corr_vv)
aux = exp.((collect(1:T) .- y0 ) .* fill(m, T))
R = ((aux .* corr_vv).^2).^0.25
R_av = plat_av(R, plat, wpm)
f = sqrt(2 / m) * R_av
R .*= sqrt.(2 ./ [m for i in 1:length(R)])
if pl || savepl
uwerr.(R)
R_av *= sqrt(2/m)
isnothing(wpm) ? uwerr(R_av) : uwerr(R_av, wpm)
isnothing(wpm) ? uwerr(f) : uwerr(f, wpm)
x = [1:length(R)...]
lbl = string(L"$af = $", sprint(show, f))
fig,ax = plot_data(R,x,label=lbl;plot_kw...)
fig,ax = plot_func(x->R_av,x[plat[1]:plat[2]],label = L"$R$",figs=(fig,ax); plot_kw...)
v,e = value(R_av),err(R_av)
ax.set_xlim(left=y0)
ax.set_ylim(v-10*e, v+10*e)
if !haskey(plot_kw,:ylabel)
ax.set_ylabel(L"$R_\mathrm{av}$")
end
if !haskey(plot_kw,:xlabel)
ax.set_xlabel(L"$x_0$")
end
if !haskey(plot_kw,:title)
if !isnothing(kappa)
fig.suptitle(string( L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
end
ax.set_title(L"$f_{B^*}$")
end
if savepl
if filename==""
filename = "dec_const_vec/"
filename *= isnothing(kappa) ? "f_B_vec.png" : string( L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2],".png")
end
mkpath(filename)
fig.tight_layout()
fig.savepl(filename)
close(fig)
else
display(gcf())
end
end
if pl&&!savepl
return data ? (f,R,fig,ax) : (f,fig,ax)
else
return data ? (f,R) : f
end
end
function dec_const(vv::Corr, plat::Vector{Int64}, m::uwreal;kwargs...)
return dec_const(vv.obs, plat, m, vv.y0, kappa=vv.kappa; kwargs...)
end
function dec_const(vv::Vector{Corr}, plat::Vector{Vector{Int64}},m::Vector{uwreal}; kwargs...)
a= dec_const.(vv,plat,m;kwargs...)
if length(a[1])==1
return a;
end
# Due to the broadcast, a is a Vector{Tuple{...}} where each element of the vector is the output of the function
# The following lines of code reshape a into a Matrix such that each column
# correspond to one Tuple and each row correspond to one component of the Tuple. The it will return each row
# separately.
#
# Example:
# V = [(a1,a2,a3),(b1,b2,b3)]
# V = collect.(V) #[[a1,a2,a3],[b1,b2,b3]]
# V stack(V) #[a1 b1; a2 b2; a3 b3]
# return Tuple([r...] for e in eachrow(r)) #([a1,b1],[a1,b2],[a3,b3])
a = stack(collect.(a));
return Tuple([r...] for r in eachrow(a))
end
# this can be used also for pseudoscalar decay constants with wilson fermions
function dec_const_w(vv::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_vv = vv[2:end-1]
corr_pp = pp[2:end-1]
T = length(corr_vv)
if ca != 0.0
der_pp = (corr_pp[3:end] - corr_pp[1:end-2]) / 2
corr_vv = corr_vv[2:end-1] + ca* der_pp
aux = exp.((collect(1:T-2) .- y0 ) .* fill(m, T-2))
else
#corr_vv = corr_vv[2:end-1]
aux = exp.((collect(1:T) .- y0 ) .* fill(m/2, T))
end
R = (aux .* corr_vv ./ sqrt.(corr_pp))
R_av = plat_av(R, plat, wpm)
f = sqrt(2 / m) * R_av
R .*= sqrt.(2 ./ [m for i in 1:length(R)])
if pl || savepl
uwerr.(R)
R_av *= sqrt(2/m)
isnothing(wpm) ? uwerr(R_av) : uwerr(R_av, wpm)
isnothing(wpm) ? uwerr(f) : uwerr(f, wpm)
x = [1:length(R)...]
lbl = string(L"$af = $", sprint(show, f))
fig,ax = plot_data(R,x,label=lbl;plot_kw...)
fig,ax = plot_func(x->R_av,x[plat[1]:plat[2]],label = L"$R$",figs=(fig,ax); plot_kw...)
v,e = value(R_av),err(R_av)
ax.set_xlim(left=y0)
ax.set_ylim(v-10*e, v+10*e)
if !haskey(plot_kw,:ylabel)
ax.set_ylabel(L"$R_\mathrm{av}$")
end
if !haskey(plot_kw,:xlabel)
ax.set_xlabel(L"$x_0$")
end
if !haskey(plot_kw,:title)
if !isnothing(kappa)
fig.suptitle(string( L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
end
ax.set_title(L"$f_{B^*}$")
end
if savepl
if filename==""
filename = "dec_const_vec/"
filename *= isnothing(kappa) ? "f_B_vec.png" : string( L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2],".png")
end
mkpath(filename)
fig.tight_layout()
fig.savepl(filename)
close(fig)
else
display(gcf())
end
end
if pl&&!savepl
return data ? (f,R,fig,ax) : (f,fig,ax)
else
return data ? (f,R) : f
end
end
function dec_const_w(vv::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)
return dec_const(vv.obs, pp.obs, plat, m, vv.y0, kappa=vv.kappa, pl=pl, data=data, wpm=wpm)
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, 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, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
dec_const_pcvc(ppL::Vector{uwreal}, ppR::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(corrL::Corr, corrR::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`
must be specified.
The flags `pl` and `data` allow to show the plots and return data as an extra result.
**The method extract the matrix element assuming that the source is in the bulk. **
**If left and right correlators are included in the input. The result is computed with a ratio that cancels boundary effects:**
`` R = \sqrt{f_P(x_0, y_0) * f_P(x_0, T - 1 - y_0) / f_P(T - 1 - y_0, y_0)}``
```@example
data = read_mesons(path, "G5", "G5")
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)
#left and right correlators
f_ratio = dec_const_pcvc(ppL, ppR, [50, 60], m)
```
"""
function 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)
"""
compute the decay constant when the source is far from the boundaries
"""
corr_pp = corr[2:end-1]
dim = length(corr_pp)
aux = exp.((collect(1:dim) .- y0 ) .* [m for k in 1:dim])
R = ((aux .* corr_pp).^2).^0.25
R_av = plat_av(R, plat, wpm)
f = sqrt(2) * (mu[1] + mu[2]) *R_av / m^1.5
if pl
isnothing(wpm) ? uwerr(f) : uwerr(f, wpm)
isnothing(wpm) ? uwerr(R_av) : uwerr(R_av, wpm)
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(1:length(R), value.(R), err.(R), fmt="x", color="black", label=lbl)
ylabel(L"$R$")
xlabel(L"$x_0$")
xlim(y0,y0+plat[2])
ylim(v-10*e, v+10*e)
legend()
title(string(L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2]))
display(gcf())
end
if !data
return f
else
return (f, uwdot(sqrt(2) * (mu[1] + mu[2]) / m^1.5, R))
end
end
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
function dec_const_pcvc(ppL::Vector{uwreal}, ppR::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)
T = length(ppL)
f1 = [ppL[T - y0] for k = 1:T]
R = ((ppL .* ppR ./ f1).^2).^(1/4)
R = R[2:end-1]
R_av = plat_av(R, plat, wpm)
f = sqrt(2) * (mu[1] + mu[2]) * R_av / m^1.5
if pl
isnothing(wpm) ? uwerr(f) : uwerr(f, wpm)
isnothing(wpm) ? uwerr(R_av) : uwerr(R_av, wpm)
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(1:length(R), value.(R), err.(R), fmt="x", color="black", label=lbl)
ylabel(L"$R$")
xlabel(L"$x_0$")
legend()
title(string(L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2]))
display(gcf())
end
if !data
return f
else
return (f, uwdot(sqrt(2) * (mu[1] + mu[2]) / m^1.5, R))
end
end
function dec_const_pcvc(corrL::Corr, corrR::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
T = length(corrL.obs)
if (corrL.kappa == corrR.kappa) && (corrL.mu == corrR.mu) && (corrL.y0 == T - 1 - corrR.y0)
return dec_const_pcvc(corrL.obs, corrR.obs, plat, m, corrL.mu, corrL.y0, pl=pl, data=data, wpm=wpm)
else
error("y0, kappa or mu values does not match")
end
end
#t0
function get_model(x, p, n)
s = 0.0
for k = 1:n
s = s .+ p[k] .* x.^(k-1)
end
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, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, info::Bool=false)
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, info::Bool=false)
Computes `t0` using the energy density of the action `Ysl`(Yang-Mills action).
`t0` is computed in the plateau `plat`.
A polynomial interpolation in `t` is performed to find `t0`, where `npol` is the degree of the polynomial (linear fit by default)
The flag `pl` allows to show the plot.
The flag `info` provides extra output that contains information about the primary observables. The function returns the primary observables ``<WY>`` and ``<W>``
(it returns the observable <Y> if rw=nothing)
```@example
#Single replica
Y = read_ms(path)
rw = read_ms(path_rw)
t0, Yobs = comp_t0(Y, [38, 58], L=32, info=true)
t0_r, WYobs, Wobs = comp_t0(Y, [38, 58], L=32, rw=rw, info=true)
#Two replicas
Y1 = read_ms(path1)
Y2 = read_ms(path2)
rw1 = read_ms(path_rw1)
rw2 = read_ms(path_rw2)
t0 = comp_t0([Y1, Y2], [38, 58], L=32, pl=true)
t0_r = comp_t0(Y, [38, 58], L=32, rw=[rw1, rw2], pl=true)
```
"""
function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false,
rw::Union{Matrix{Float64}, Nothing}=nothing, npol::Int64=2, ws::ADerrors.wspace=ADerrors.wsg,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, info::Bool=false)
Ysl = Y.obs
t = Y.t
id = Y.id
replica = size.([Ysl], 1)
#Truncation
if id in keys(ADerrors.wsg.str2id)
n_ws = findfirst(x-> x == ws.str2id[id], ws.map_nob)
if !isnothing(n_ws)
ivrep_ws = ws.fluc[n_ws].ivrep
if length(ivrep_ws) != 1
error("Different number of replicas")
end
if replica[1] > ivrep_ws[1]
println("Automatic truncation in Ysl ", ivrep_ws[1], " / ", replica[1], ". R = 1")
Ysl = Ysl[1:ivrep_ws[1], :, :]
elseif replica[1] < ivrep_ws[1]
println(replica[1])
error("Automatic truncation failed. R = 1\nTry using truncate_data!")
end
end
end
xmax = size(Ysl, 2)
nt0 = t0_guess(t, Ysl, plat, L)
dt0 = iseven(npol) ? Int64(npol / 2) : Int64((npol+1)/ 2)
Y_aux = Matrix{uwreal}(undef, xmax, 2*dt0+1)
if !isnothing(rw)
Ysl_r, W = apply_rw(Ysl, rw)
W_obs = uwreal(W, id)
WY_aux = Matrix{uwreal}(undef, xmax, 2*dt0+1)
end
for i = 1:xmax
k = 1
for j = nt0-dt0:nt0+dt0
if isnothing(rw)
Y_aux[i, k] = uwreal(Ysl[:, i, j], id)
else
WY_aux[i, k] = uwreal(Ysl_r[:, i, j], id)
Y_aux[i, k] = WY_aux[i, k] / W_obs
end
k = k + 1
end
end
x = t[nt0-dt0:nt0+dt0]
t2E = [plat_av(Y_aux[:, j], plat, wpm) for j=1:2*dt0+1] .* x.^2 / L^3
model(x, p) = get_model(x, p, npol)
par, chi = fit_routine(model, x, t2E, npol)
fmin(x, p) = model(x, p) .- 0.3
t0 = root_error(fmin, t[nt0], par)
if pl
if isnothing(wpm)
uwerr(t0)
uwerr.(t2E)
else
uwerr(t0, wpm)
[uwerr(t2E_aux, wpm) for t2E_aux in t2E]
end
v = value.(t2E)
e = err.(t2E)
figure()
errorbar(x, v, e, fmt="d", capsize=2, mfc="none", color="black")
errorbar(value(t0), 0.3, xerr=err(t0), capsize=2, color="red", fmt="d")
xarr = Float64.(range(5.0, 5.2, length=100))
yarr = model(xarr, par); uwerr.(yarr)
fill_between(xarr, value.(yarr) .- err.(yarr), value.(yarr) .+ err.(yarr), color="royalblue", alpha=0.2 )
ylabel(L"$t^2\langle E(x_0, t)\rangle$")
xlabel(L"$t/a^2$")
tight_layout()
savefig("/Users/alessandroconigli/Desktop/t0_interp.pdf")
display(gcf())
figure()
t_pl = dt0 + 1
yyy = Y_aux[:, t_pl] * t[nt0]^2 / L^3 ; uwerr.(yyy)
# plot(collect(1:xmax), value.(Y_aux[:, t_pl]) * t[nt0]^2 / L^3, "x")
errorbar(collect(1:xmax), value.(yyy), err.(yyy), fmt="d", mfc="none", color="black", capsize=2)
fill_between(collect(plat[1]:plat[2]), v[t_pl]+e[t_pl], v[t_pl]-e[t_pl], alpha=0.75, color="royalblue")
ylabel(L"$t^2\langle E(x_0, t)\rangle$")
xlabel(L"$x_0/a$")
ylim(v[t_pl]-30*e[t_pl], v[t_pl]+50*e[t_pl])
title(string(L"$t/a^2 = $", t[nt0]))
tight_layout()
savefig("/Users/alessandroconigli/Desktop/t0_plateau.pdf")
display(gcf())
end
if info && !isnothing(rw)
return (t0, WY_aux, W_obs)
elseif info && isnothing(rw)
return (t0, Y_aux)
else
return t0
end
end
function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false,
rw::Union{Vector{Matrix{Float64}}, Nothing}=nothing, npol::Int64=2, ws::ADerrors.wspace=ADerrors.wsg,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, info::Bool=false)
nr = length(Y)
Ysl = getfield.(Y, :obs)
t = getfield.(Y, :t)
t = t[1]
id = getfield.(Y, :id)
replica = size.(Ysl, 1)
if !all(id .== id[1])
error("IDs are not equal")
end
id = id[1]
#Truncation
if id in keys(ADerrors.wsg.str2id)
n_ws = findfirst(x-> x == ws.str2id[id], ws.map_nob)
if !isnothing(n_ws)
ivrep_ws = ws.fluc[n_ws].ivrep
if length(replica) != length(ivrep_ws)
error("Different number of replicas")
end
for k = 1:length(replica)
if replica[k] > ivrep_ws[k]
println("Automatic truncation in Ysl ", ivrep_ws[k], " / ", replica[k], ". R = ", k)
Ysl[k] = Ysl[k][1:ivrep_ws[k], :, :]
elseif replica[k] < ivrep_ws[k]
error("Automatic truncation failed. R = ", replica[k], "\nTry using truncate_data!")
end
end
replica = size.(Ysl, 1)
end
end
tmp = Ysl[1]
[tmp = cat(tmp, Ysl[k], dims=1) for k=2:nr]
nt0 = t0_guess(t, tmp, plat, L)
xmax = size(tmp, 2)
dt0 = iseven(npol) ? Int64(npol / 2) : Int64((npol+1) / 2)
Y_aux = Matrix{uwreal}(undef, xmax, 2*dt0+1)
if !isnothing(rw)
Ysl_r, W = apply_rw(Ysl, rw)
tmp_r = Ysl_r[1]
tmp_W = W[1]
[tmp_r = cat(tmp_r, Ysl_r[k], dims=1) for k = 2:nr]
[tmp_W = cat(tmp_W, W[k], dims=1) for k = 2:nr]
W_obs = uwreal(tmp_W, id, replica)
WY_aux = Matrix{uwreal}(undef, xmax, 2*dt0+1)
end
for i = 1:xmax
k = 1
for j = nt0-dt0:nt0+dt0
if isnothing(rw)
Y_aux[i, k] = uwreal(tmp[:, i, j], id, replica)
else
WY_aux[i, k] = uwreal(tmp_r[:, i, j], id, replica)
Y_aux[i, k] = WY_aux[i, k] / W_obs
end
k = k + 1
end
end
x = t[nt0-dt0:nt0+dt0]
t2E = [plat_av(Y_aux[:, j], plat, wpm) for j=1:2*dt0+1] .* x.^2 / L^3
model(x, p) = get_model(x, p, npol)
par, chi = fit_routine(model, x, t2E, npol)
fmin(x, p) = model(x, p) .- 0.3
t0 = root_error(fmin, t[nt0], par)
if pl
if isnothing(wpm)
uwerr(t0)
uwerr.(t2E)
else
uwerr(t0, wpm)
[uwerr(t2E_aux, wpm) for t2E_aux in t2E]
end
v = value.(t2E)
e = err.(t2E)
figure()
errorbar(x, v, e, fmt="x")
errorbar(value(t0), 0.3, xerr=err(t0), fmt="x")
ylabel(L"$t^2E$")
xlabel(L"$t/a^2$")
display(gcf())
figure()
t_pl = dt0 + 1
plot(collect(1:xmax), value.(Y_aux[:, t_pl]) * t[nt0]^2 / L^3, "x")
fill_between(collect(plat[1]:plat[2]), v[t_pl]+e[t_pl], v[t_pl]-e[t_pl], alpha=0.75, color="green")
ylabel(L"$t^2E(x0, t)$")
xlabel(L"$x_0/a$")
title(string(L"$t/a^2 = $", t[nt0]))
display(gcf())
end
if info && !isnothing(rw)
return (t0, WY_aux, W_obs)
elseif info && isnothing(rw)
return (t0, Y_aux)
else
return t0
end
end
function t0_guess(t::Vector{Float64}, Ysl::Array{Float64, 3}, plat::Vector{Int64}, L::Int64)
t2E_ax = t.^2 .* mean(mean(Ysl[:, plat[1]:plat[2], :], dims=2), dims=1)[1, 1, :] / L^3
#look for values t2E: t2E > 0.3 and 0.0 < t2E_ax < 0.3
if !(any(t2E_ax .> 0.3) && any((t2E_ax .< 0.3) .& (t2E_ax .> 0.0)))
error("Error finding solutions. Check volume normalization!")
end
t0_aux = minimum(abs.(t2E_ax .- 0.3))
nt0 = findfirst(x-> abs(x - 0.3) == t0_aux, t2E_ax) #index closest value to t0
return nt0
end
function get_m(corr::juobs.Corr, ens::EnsInfo, PS::String;
tm::Union{Vector{Vector{Int64}}, Nothing}=nothing, tM::Union{Vector{Vector{Int64}}, Nothing}=nothing,
pl::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
corr_d = corr.obs
m_dat = 0.5 .* log.((corr_d[2:end-2] ./ corr_d[3:end-1]) .^ 2)
y0 = corr.y0
T = length(corr_d) - 1 - y0
isnothing(tm) ? tm = [[y0+10,y0+15,y0+20,y0+25,y0+30,y0+35,y0+40], [i for i in Int(round(T / 3)):Int(round(T / 3))+10]] : tm=tm
isnothing(tM) ? tM = [[T-10,T-15,T-20,T-25,T-30,T-35,T-40], [i for i in Int(round(2 * T / 3)):Int(round(2 * T / 3))+10]] : tM=tM
@.fit_exp(x,p) = p[1] + p[2] * exp(-p[3] * (x-y0)) + p[4] * exp(-p[5] * (T-x))
@.fit_const(x,p) = p[1] * x ^ 0
k1 = 5
k2 = 1
guess = value(m_dat[Int64(round(T / 2))])
m, syst, m_i, weight, pval = model_av([fit_exp, fit_const], m_dat, guess, tm=tm, tM=tM, k=[k1,k2], wpm=wpm)
if pl == true
isnothing(wpm) ? uwerr(m) : uwerr(m, wpm)
isnothing(wpm) ? uwerr.(m_dat) : [uwerr(m_dat[i], wpm) for i in 1:length(m_dat)]
isnothing(wpm) ? uwerr.(m_i) : [uwerr(m_i[i], wpm) for i in 1:length(m_i)]
v = value(m)
e = err(m)
fig = figure("eff")
errorbar(1:length(m_dat), value.(m_dat), err.(m_dat), fmt="x", color="black")
ylabel(L"$am_\mathrm{eff}$")
xlabel(L"$x_0/a$")
ylim(v-7*e, v+20*e)
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/m_",ens.id,"_",PS,"_plat.pdf"))
models = Array{String,1}()
for i in 1:length(tm[1])
for j in 1:length(tM[1])
push!(models, string("[",tm[1][i],",",tM[1][j],"]"))
end
end
for i in 1:length(tm[2])
for j in 1:length(tM[2])
push!(models, string("[",tm[2][i],",",tM[2][j],"]"))
end
end
fig = figure("model av")
subplot(411)
fill_between(1:length(m_dat), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(m_dat), value.(m_dat), err.(m_dat), fmt="x", color="black")
ylabel(L"$am_\mathrm{eff}$")
xlabel(L"$x_0/a$")
ylim(v-10*e, v+20*e)
ax = gca()
setp(ax.get_xticklabels(),visible=false)
subplot(412)
ylabel(L"$m_i$")
fill_between(1:length(m_i), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(m_i), value.(m_i), err.(m_i), fmt="x", color="black")
ax = gca()
ax[:set_xlim]([0, length(models)+1])
setp(ax.get_xticklabels(),visible=false)
subplot(413)
ylabel(L"$W_i$")
bar(1:length(m_i), weight, color="green")
ax = gca()
ax[:set_xlim]([0, length(models)+1])
setp(ax.get_xticklabels(),visible=false)
subplot(414)
ylabel("p-value")
bar(1:length(m_i), pval, color="green")
ax = gca()
ax[:set_xlim]([0, length(models)+1])
plt.xticks(collect(1:length(m_i)), models)
xticks(rotation=90)
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/m_",ens.id,"_",PS,".pdf"))
#close("all")
end
return m, syst, m_i, weight, pval
end
function get_m_pbc(corr::juobs.Corr, ens::EnsInfo, PS::String;
tm::Union{Vector{Int64}, Nothing}=nothing, tM::Union{Vector{Int64}, Nothing}=nothing,
pl::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,
method::String="cosh")
corr_d = corr.obs
m_dat = 0.5 .* log.((corr_d[2:end-2] ./ corr_d[3:end-1]) .^ 2)
if ens.id == "E250"
global T = 192
global T2 = 192
global Thalf = 96
elseif ens.id == "D450"
global T = 128
global T2 = 128
global Thalf = 64
end
guess = value(m_dat[Int64(round(T / 3))])
isnothing(tm) ? tm = [y0+10,y0+15,y0+20,y0+25,y0+30,y0+35,y0+40] : tm=tm
isnothing(tM) ? tM = [Thalf] : tM=tM
aux = [corr_d[i] / corr_d[i+1] for i in 2:length(corr_d)-1]
if ens.id == "E250"
global T = 96
elseif ens.id == "D450"
global T = 64
end
@.fit_fun(x,p) = cosh(p[1] * (x-T)) / cosh(p[1] * (x+1-T))
k1 = 1
aux2 = corr_d
@.fit_fun2(x,p) = p[2] * exp(-p[1] * (x-1)) + p[2] * exp(-p[1] * (T2+1-x))
k2 = 2
if method == "cosh"
m, syst, m_i, weight, pval = model_av(fit_fun, aux, guess, tm=tm, tM=tM, k=k1, wpm=wpm)
elseif method == "corr"
m, syst, m_i, weight, pval = model_av(fit_fun2, aux2, guess, tm=tm, tM=tM, k=k2, wpm=wpm)
end
if pl == true
##TODO
bla = 1
end
return m, syst, m_i, weight, pval
end
function get_mpcac(corr_pp::juobs.Corr, corr_ap::juobs.Corr, ens::EnsInfo, PS::String;
tm::Union{Vector{Vector{Int64}}, Nothing}=nothing, tM::Union{Vector{Vector{Int64}}, Nothing}=nothing,
impr::Bool=true, pl::Bool=false,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
ap_dat = -corr_ap.obs[2:end-1]
pp_dat = corr_pp.obs[2:end-1]
der_ap = (ap_dat[3:end] .- ap_dat[1:end-2]) / 2
if impr == true
ca = ens.ca
der2_pp = pp_dat[1:end-2] + pp_dat[3:end] - 2*pp_dat[2:end-1]
der_ap = der_ap + ca*der2_pp
end
mpcac_dat = der_ap ./ (2*pp_dat[2:end-1])
y0 = corr_pp.y0
T = length(pp_dat) - 1 - y0
isnothing(tm) ? tm = [[y0+10,y0+15,y0+20,y0+25,y0+30,y0+35,y0+40], [i for i in Int(round(T / 3)):Int(round(T / 3))+10]] : tm=tm
isnothing(tM) ? tM = [[T-10,T-15,T-20,T-25,T-30,T-35,T-40], [i for i in Int(round(2 * T / 3)):Int(round(2 * T / 3))+10]] : tM=tM
@.fit_exp(x,p) = p[1] + p[2] * exp(-p[3] * (x-y0)) + p[4] * exp(-p[5] * (T-x))
@.fit_const(x,p) = p[1] * x ^ 0
k1 = 5
k2 = 1
guess = value(mpcac_dat[Int64(round(T / 2))])
mpcac, syst, mpcac_i, weight, pval = model_av([fit_exp, fit_const], mpcac_dat, guess, tm=tm, tM=tM, k=[k1,k2], wpm=wpm)
if pl == true
isnothing(wpm) ? uwerr(mpcac) : uwerr(mpcac, wpm)
isnothing(wpm) ? uwerr.(mpcac_dat) : [uwerr(mpcac_dat[i], wpm) for i in 1:length(mpcac_dat)]
isnothing(wpm) ? uwerr.(mpcac_i) : [uwerr(mpcac_i[i], wpm) for i in 1:length(mpcac_i)]
v = value(mpcac)
e = err(mpcac)
fig = figure("eff")
errorbar(1:length(mpcac_dat), value.(mpcac_dat), err.(mpcac_dat), fmt="x", color="black")
ylabel(L"$am_\mathrm{pcac}$")
xlabel(L"$x_0/a$")
ylim(v-10*e, v+20*e)
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/mpcac_",ens.id,"_",PS,"_plat.pdf"))
fig = figure("model av")
subplot(411)
fill_between(1:length(mpcac_dat), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(mpcac_dat), value.(mpcac_dat), err.(mpcac_dat), fmt="x", color="black")
ylabel(L"$am_\mathrm{pcac}$")
xlabel(L"$x_0/a$")
ylim(v-10*e, v+10*e)
subplot(412)
ylabel(L"$mpcac_i$")
fill_between(1:length(mpcac_i), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(mpcac_i), value.(mpcac_i), err.(mpcac_i), fmt="x", color="black")
subplot(413)
ylabel(L"$W_i$")
bar(1:length(mpcac_i), weight, color="green")
subplot(414)
ylabel("p-value")
bar(1:length(mpcac_i), pval, color="green")
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/mpcac_",ens.id,"_",PS,".pdf"))
#close("all")
end
return mpcac, syst, mpcac_i, weight, pval
end
function get_f_wil(corr_pp::juobs.Corr, corr_ap::juobs.Corr, m::uwreal, ens::EnsInfo, PS::String;
tm::Union{Vector{Vector{Int64}}, Nothing}=nothing, tM::Union{Vector{Vector{Int64}}, Nothing}=nothing,
impr::Bool=true, pl::Bool=false,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
ap_dat = -corr_ap.obs
pp_dat = corr_pp.obs
T = length(pp_dat)
y0 = corr_pp.y0
if impr == true
ca = ens.ca
der_pp = (pp_dat[3:end] .- pp_dat[1:end-2]) / 2
ap_dat = ap_dat[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-1)/2) .* [m for k = 1:T])
end
R_dat = ap_dat .* aux ./ [((pp_dat[T-y0])^2)^(1/4) for k = 1:length(ap_dat)]
#f_dat = [sqrt(2) * sqrt(R_dat[i] ^ 2) / sqrt(m) for i in 1:length(R_dat)]
isnothing(tm) ? tm = [[y0+10,y0+15,y0+20,y0+25,y0+30,y0+35,y0+40], [i for i in Int(round(T / 3)):Int(round(T / 3))+10]] : tm=tm
isnothing(tM) ? tM = [[T-10,T-15,T-20,T-25,T-30,T-35,T-40], [i for i in Int(round(2 * T / 3)):Int(round(2 * T / 3))+10]] : tM=tM
@.fit_exp(x,p) = p[1] + p[2] * exp(-p[3] * (x-y0)) + p[4] * exp(-p[5] * (T-1-y0-x))
@.fit_const(x,p) = p[1] * x ^ 0
k1 = 5
k2 = 1
guess = value(R_dat[Int64(round(T / 2))])
R, syst, R_i, weight, pval = model_av([fit_exp, fit_const], R_dat, guess, tm=tm, tM=tM, k=[k1,k2], wpm=wpm)
f = sqrt(2) * sqrt(R^2) / sqrt(m)
if pl == true
isnothing(wpm) ? uwerr(R) : uwerr(R, wpm)
isnothing(wpm) ? uwerr.(R_dat) : [uwerr(R_dat[i], wpm) for i in 1:length(R_dat)]
isnothing(wpm) ? uwerr.(R_i) : [uwerr(R_i[i], wpm) for i in 1:length(R_i)]
v = value(R)
e = err(R)
fig = figure("eff")
errorbar(1:length(R_dat), value.(R_dat), err.(R_dat), fmt="x", color="black")
ylabel(L"$R_\mathrm{PS}$")
xlabel(L"$x_0/a$")
ylim(v-20*e, v+20*e)
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/R_",ens.id,"_",PS,"_plat.pdf"))
fig = figure("model av")
subplot(411)
fill_between(1:length(R_dat), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(R_dat), value.(R_dat), err.(R_dat), fmt="x", color="black")
ylabel(L"$R_\mathrm{PS}$")
xlabel(L"$x_0/a$")
ylim(v-10*e, v+10*e)
subplot(412)
ylabel(L"$R_i$")
fill_between(1:length(R_i), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(R_i), value.(R_i), err.(R_i), fmt="x", color="black")
subplot(413)
ylabel(L"$W_i$")
bar(1:length(R_i), weight, color="green")
subplot(414)
ylabel("p-value")
bar(1:length(R_i), pval, color="green")
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/R_",ens.id,"_",PS,".pdf"))
#close("all")
end
return f, syst, R_i, weight, pval
end
function get_f_wil(corr_ppL::juobs.Corr, corr_ppR::juobs.Corr, corr_apL::juobs.Corr, corr_apR::juobs.Corr,
m::uwreal, ens::EnsInfo, PS::String; tm::Union{Vector{Int64}, Nothing}=nothing,
tM::Union{Vector{Int64}, Nothing}=nothing, impr::Bool=true, pl::Bool=false,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
pp_dat = (corr_ppL.obs + corr_ppR.obs[end:-1:1]) / 2
T = length(corr_ppL.obs)
y0 = corr_ppL.y0
if impr == true
ca = ens.ca
der_ppL = (corr_ppL.obs[3:end] - corr_ppL.obs[1:end-2]) / 2
der_ppR = (corr_ppR.obs[3:end] - corr_ppR.obs[1:end-2]) / 2
apL_dat = -corr_apL.obs[2:end-1] + ca * der_ppL
apR_dat = -corr_apR.obs[2:end-1] + ca * der_ppR
else
apL_dat = -corr_apL.obs[2:end-1]
apR_dat = -corr_apR.obs[2:end-1]
end
f1 = [pp_dat[T - y0] for k = 1:length(apL_dat)]
R_dat = ((apL_dat .* apR_dat ./ f1).^2).^(1/4)
#f_dat = [sqrt(2) * sqrt(R_dat[i] ^ 2) / sqrt(m) for i in 1:length(R_dat)]
isnothing(tm) ? tm = [[y0+10,y0+15,y0+20,y0+25,y0+30,y0+35,y0+40], [i for i in Int(round(T / 3)):Int(round(T / 3))+10]] : tm=tm
isnothing(tM) ? tM = [[T-10,T-15,T-20,T-25,T-30,T-35,T-40], [i for i in Int(round(2 * T / 3)):Int(round(2 * T / 3))+10]] : tM=tM
@.fit_exp(x,p) = p[1] + p[2] * exp(-p[3] * (x-y0)) + p[4] * exp(-p[5] * (T-1-y0-x))
@.fit_const(x,p) = p[1] * x ^ 0
k1 = 5
k2 = 1
guess = value(R_dat[Int64(round(T / 2))])
R, syst, R_i, weight, pval = model_av([fit_exp, fit_const], R_dat, guess, tm=tm, tM=tM, k=[k1,k2], wpm=wpm)
f = sqrt(2) * sqrt(R^2) / sqrt(m)
if pl == true
isnothing(wpm) ? uwerr(R) : uwerr(R, wpm)
isnothing(wpm) ? uwerr.(R_dat) : [uwerr(R_dat[i], wpm) for i in 1:length(R_dat)]
isnothing(wpm) ? uwerr.(R_i) : [uwerr(R_i[i], wpm) for i in 1:length(R_i)]
v = value(R)
e = err(R)
fig = figure("eff")
errorbar(1:length(R_dat), value.(R_dat), err.(R_dat), fmt="x", color="black")
ylabel(L"$R_\mathrm{PS}$")
xlabel(L"$x_0/a$")
ylim(v-20*e, v+20*e)
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/R_",ens.id,"_",PS,"_plat.pdf"))
fig = figure("model av")
subplot(411)
fill_between(1:length(R_dat), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(R_dat), value.(R_dat), err.(R_dat), fmt="x", color="black")
ylabel(L"$R_\mathrm{PS}$")
xlabel(L"$x_0/a$")
ylim(v-10*e, v+10*e)
subplot(412)
ylabel(L"$R_i$")
fill_between(1:length(R_i), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(R_i), value.(R_i), err.(R_i), fmt="x", color="black")
subplot(413)
ylabel(L"$W_i$")
bar(1:length(R_i), weight, color="green")
subplot(414)
ylabel("p-value")
bar(1:length(R_i), pval, color="green")
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/R_",ens.id,"_",PS,".pdf"))
#close("all")
end
return f, syst, R_i, weight, pval
end
function get_f_wil_pbc(corr_pp::juobs.Corr, corr_ap::juobs.Corr, m::uwreal, ens::EnsInfo, PS::String;
tm::Union{Vector{Int64}, Nothing}=nothing, tM::Union{Vector{Int64}, Nothing}=nothing,
impr::Bool=true, pl::Bool=false, method::String="ratio",
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
ap_dat = -corr_ap.obs
pp_dat = corr_pp.obs
if ens.id == "E250"
global T = 192
global Thalf = 96
elseif ens.id == "D450"
global T = 128
global Thalf = 64
end
if impr == true
ca = ens.ca
der_pp = (pp_dat[3:end] .- pp_dat[1:end-2]) / 2
ap_dat = ap_dat[2:end-1] + ca * der_pp
pp_dat = pp_dat[2:end-1]
end
R_dat = ap_dat ./ sqrt.(pp_dat)
isnothing(tm) ? tm = [y0+10,y0+15,y0+20,y0+25,y0+30,y0+35,y0+40] : tm=tm
isnothing(tM) ? tM = [Thalf] : tM=tM
@.fit_fun(x,p) = p[1] * (-exp(-value(m) * (x-1)) + exp(-value(m) * (T+1-x))) / sqrt(exp(-value(m) * (x-1)) + exp(-value(m) * (T+1-x)))
@.fit_fun_ap(x,p) = p[1] * (-exp(-value(m) * (x-1)) + exp(-value(m) * (T+1-x)))
@.fit_fun_pp(x,p) = p[1] * (exp(-value(m) * (x-1)) + exp(-value(m) * (T+1-x)))
@.fit_const(x,p) = p[1] * x ^ 0
k = 1
function fit_corr_sim(x,p)
f = [p[1] * (-exp(-value(m) * x[i]) + exp(-value(m) * (T - x[i]))) for i in 1:div(length(x),2)]
g = [p[2] * (exp(-value(m) * x[i]) + exp(-value(m) * (T - x[i]))) for i in div(length(x),2)+1:length(x)]
return [f;g]
end
if method == "ratio"
R, syst, R_i, weight, pval = model_av(fit_fun, R_dat, .5, tm=tm, tM=tM, k=k, wpm=wpm)
f_i = sqrt.(2 ./ [m for i in 1:length(R_i)]) .* R_i
f = sqrt(2 / m) * R
elseif method == "corr"
cap, syst_ap, cap_i, weight_ap, pval_ap = model_av(fit_fun_ap, ap_dat, .5, tm=tm, tM=tM, k=k, wpm=wpm)
cpp, syst_pp, cpp_i, weight_pp, pval_pp = model_av(fit_fun_pp, pp_dat, .5, tm=tm, tM=tM, k=k, wpm=wpm)
f_i = sqrt.(2 ./ [m for i in 1:length(cap_i)]) .* cap_i ./ sqrt.(cpp_i)
f = sqrt(2 / m) * cap / sqrt(cpp)
syst = [syst_ap, syst_pp]
pval = [pval_ap, pval_pp]
weight = [weight_ap, weight_pp]
elseif method == "corr_sim"
pval = Array{Float64,1}()
TIC = Array{Float64,1}()
cap_i = Array{uwreal,1}()
cpp_i = Array{uwreal,1}()
for i in tm
y = [ap_dat[i:end]; pp_dat[i:end]]
x = collect(i:length(ap_dat))
x = [x; x]
uwerr.(y)
W = 1 ./ ADerrors.err.(y) .^ 2
uprm, chi2, chi_exp, pv = fit_alg(fit_corr_sim, x, y, 2, [.5, .5])
push!(cap_i, uprm[1])
push!(cpp_i, uprm[2])
push!(pval, pv)
push!(TIC, chi2 - 2 * chi_exp)
end
weight = exp.(-0.5 * TIC) ./ sum(exp.(-0.5 * TIC))
R_i = cap_i ./ sqrt.(cpp_i)
R_av = sum(R_i .* weight)
syst = sum(R_i .^ 2 .* weight) - R_av ^ 2
f_i = sqrt.(2 ./ [m for i in 1:length(cap_i)]) .* R_i
f = sqrt(2 / m) * R_av
elseif method == "pcac"
if impr == false
error("pcac method implemented only for improved axial current, set impr=true")
end
cpp, syst_pp, cap_i, weight, pval = model_av(fit_fun_pp, pp_dat, .5, tm=tm, tM=tM, k=k, wpm=wpm)
der_ap = (ap_dat[3:end] .- ap_dat[1:end-2]) / 2
der2_pp = pp_dat[1:end-2] + pp_dat[3:end] - 2*pp_dat[2:end-1]
der_ap = der_ap + ca*der2_pp
mpcac_dat = der_ap ./ (2*pp_dat[2:end-1])
f_dat = [sqrt(8 / m ^ 3 * cpp) for i in 1:length(mpcac_dat)] .* mpcac_dat
f, syst, f_i, weight, pval = model_av(fit_const, f_dat, .5, tm=tm, tM=tM, k=k, wpm=wpm)
end
if pl == true
##TODO
bla = 1
x = collect(1:length(R_dat))
R_dat = R_dat ./ [(-exp(-value(m) * (x[i]-1)) + exp(-value(m) * (T+1-x[i]))) / sqrt(exp(-value(m) * (x[i]-1)) + exp(-value(m) * (T+1-x[i]))) for i in 1:length(x)]
isnothing(wpm) ? uwerr.(R_dat) : [uwerr(R_dat[i], wpm) for i in 1:length(R_dat)]
v = value(R_dat[40])
e = err(R_dat[40])
figure()
errorbar(1:length(R_dat), value.(R_dat), err.(R_dat), fmt="x", color="black")
ylabel(L"$R_\mathrm{PS}$")
xlabel(L"$x_0/a$")
ylim(v-20*e, v+20*e)
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/R_",ens.id,"_",PS,"_plat.pdf"))
end
return f, syst, f_i, weight, pval
end
function get_f_tm(corr_pp::juobs.Corr, m::uwreal, ens::EnsInfo, PS::String;
tm::Union{Vector{Vector{Int64}}, Nothing}=nothing, tM::Union{Vector{Vector{Int64}}, Nothing}=nothing,
pl::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
T = length(corr_pp.obs)
y0 = corr_pp.y0
mu = corr_pp.mu
pp_dat = corr_pp.obs
aux = exp.((collect(0:T-1) .- (T-1)/2 ) .* [m for k in 1:T])
R_dat = pp_dat .* aux ./ [((pp_dat[T-y0])^2)^(1/4) for k = 1:T]
isnothing(tm) ? tm = [[y0+10,y0+15,y0+20,y0+25,y0+30,y0+35,y0+40], [i for i in Int(round(T / 3)):Int(round(T / 3))+10]] : tm=tm
isnothing(tM) ? tM = [[T-10,T-15,T-20,T-25,T-30,T-35,T-40], [i for i in Int(round(2 * T / 3)):Int(round(2 * T / 3))+10]] : tM=tM
@.fit_exp(x,p) = p[1] + p[2] * exp(-p[3] * (x-y0)) + p[4] * exp(-p[5] * (T-1-y0-x))
@.fit_const(x,p) = p[1] * x ^ 0
k1 = 5
k2 = 1
guess = value(R_dat[Int64(round(T / 2))])
R, syst, R_i, weight, pval = model_av([fit_exp, fit_const], R_dat, guess, tm=tm, tM=tM, k=[k1,k2], wpm=wpm)
f = sqrt(2) * (mu[1] + mu[2]) * R / m^1.5
if pl == true
isnothing(wpm) ? uwerr(R) : uwerr(R, wpm)
isnothing(wpm) ? uwerr.(R_dat) : [uwerr(R_dat[i], wpm) for i in 1:length(R_dat)]
isnothing(wpm) ? uwerr.(R_i) : [uwerr(R_i[i], wpm) for i in 1:length(R_i)]
v = value(R)
e = err(R)
fig = figure("eff")
errorbar(1:length(R_dat), value.(R_dat), err.(R_dat), fmt="x", color="black")
ylabel(L"$R_\mathrm{PS}$")
xlabel(L"$x_0/a$")
ylim(v-20*e, v+20*e)
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/R_",ens.id,"_",PS,"_plat.pdf"))
fig = figure("model av")
subplot(411)
fill_between(1:length(R_dat), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(R_dat), value.(R_dat), err.(R_dat), fmt="x", color="black")
ylabel(L"$R_\mathrm{PS}$")
xlabel(L"$x_0/a$")
ylim(v-10*e, v+10*e)
subplot(412)
ylabel(L"$R_i$")
fill_between(1:length(R_i), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(R_i), value.(R_i), err.(R_i), fmt="x", color="black")
subplot(413)
ylabel(L"$W_i$")
bar(1:length(R_i), weight, color="green")
subplot(414)
ylabel("p-value")
bar(1:length(R_i), pval, color="green")
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/R_",ens.id,"_",PS,".pdf"))
#close("all")
end
return f, syst, R_i, weight, pval
end
function get_f_tm(corr_ppL::juobs.Corr, corr_ppR::juobs.Corr, m::uwreal, ens::EnsInfo, PS::String;
tm::Union{Vector{Vector{Int64}}, Nothing}=nothing, tM::Union{Vector{Vector{Int64}}, Nothing}=nothing,
pl::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
T = length(corr_ppL.obs)
y0 = corr_ppL.y0
mu = corr_ppL.mu
ppL_dat = corr_ppL.obs
ppR_dat = corr_ppR.obs
f1 = [ppL_dat[T - y0] for k = 1:T]
R_dat = ((ppL_dat .* ppR_dat ./ f1).^2).^(1/4)
isnothing(tm) ? tm = [[y0+10,y0+15,y0+20,y0+25,y0+30,y0+35,y0+40], [i for i in Int(round(T / 3)):Int(round(T / 3))+10]] : tm=tm
isnothing(tM) ? tM = [[T-10,T-15,T-20,T-25,T-30,T-35,T-40], [i for i in Int(round(2 * T / 3)):Int(round(2 * T / 3))+10]] : tM=tM
@.fit_exp(x,p) = p[1] + p[2] * exp(-p[3] * (x-y0)) + p[4] * exp(-p[5] * (T-1-y0-x))
@.fit_const(x,p) = p[1] * x ^ 0
k1 = 5
k2 = 1
guess = value(R_dat[Int64(round(T / 2))])
R, syst, R_i, weight, pval = model_av([fit_exp, fit_const], R_dat, guess, tm=tm, tM=tM, k=[k1,k2], wpm=wpm)
f = sqrt(2) * (mu[1] + mu[2]) * R / m^1.5
if pl == true
isnothing(wpm) ? uwerr(R) : uwerr(R, wpm)
isnothing(wpm) ? uwerr.(R_dat) : [uwerr(R_dat[i], wpm) for i in 1:length(R_dat)]
isnothing(wpm) ? uwerr.(R_i) : [uwerr(R_i[i], wpm) for i in 1:length(R_i)]
v = value(R)
e = err(R)
fig = figure("eff")
errorbar(1:length(R_dat), value.(R_dat), err.(R_dat), fmt="x", color="black")
ylabel(L"$R_\mathrm{PS}$")
xlabel(L"$x_0/a$")
ylim(v-20*e, v+20*e)
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/R_",ens.id,"_",PS,"_plat.pdf"))
fig = figure("model av")
subplot(411)
fill_between(1:length(R_dat), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(R_dat), value.(R_dat), err.(R_dat), fmt="x", color="black")
ylabel(L"$R_\mathrm{PS}$")
xlabel(L"$x_0/a$")
ylim(v-10*e, v+10*e)
subplot(412)
ylabel(L"$R_i$")
fill_between(1:length(R_i), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(R_i), value.(R_i), err.(R_i), fmt="x", color="black")
subplot(413)
ylabel(L"$W_i$")
bar(1:length(R_i), weight, color="green")
subplot(414)
ylabel("p-value")
bar(1:length(R_i), pval, color="green")
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/R_",ens.id,"_",PS,".pdf"))
#close("all")
end
return f, syst, R_i, weight, pval
end
function get_f_tm_pbc(corr_pp::juobs.Corr, m::uwreal, ens::EnsInfo, PS::String;
tm::Union{Vector{Int64}, Nothing}=nothing, tM::Union{Vector{Int64}, Nothing}=nothing,
pl::Bool=false,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
pp_dat = corr_pp.obs[2:end]
mu = corr_pp.mu
if ens.id == "E250"
global T = 192
global Thalf = 96
elseif ens.id == "D450"
global T = 128
global Thalf = 64
end
isnothing(tm) ? tm = [y0+10,y0+15,y0+20,y0+25,y0+30,y0+35,y0+40] : tm=tm
isnothing(tM) ? tM = [Thalf] : tM=tM
@.fit_fun_pp(x,p) = p[1] * (exp(-value(m) * (x-1)) + exp(-value(m) * (T+1-x)))
k = 1
cpp, syst, cpp_i, weight, pval = model_av(fit_fun_pp, pp_dat, .5, tm=tm, tM=tM, k=k, wpm=wpm)
f_i = [sqrt(2) * (mu[1] + mu[2]) / m ^ 1.5 for i in 1:length(cpp_i)] .* sqrt.(cpp_i)
f = sqrt(2) * (mu[1] + mu[2]) / m ^ 1.5 * sqrt(cpp)
if pl == true
##TODO
bla = 1
x = collect(1:length(pp_dat))
R_dat = pp_dat #./ [exp(-m * (x[i]-1)) + exp(-m * (T+1-x[i])) for i in 1:length(x)]
vec = cpp .* [exp(-m * (x[i]-1)) + exp(-m * (T+1-x[i])) for i in 1:length(x)]
isnothing(wpm) ? uwerr(cpp) : uwerr(cpp, wpm)
isnothing(wpm) ? uwerr.(R_dat) : [uwerr(R_dat[i], wpm) for i in 1:length(R_dat)]
isnothing(wpm) ? uwerr.(vec) : [uwerr(vec[i], wpm) for i in 1:length(vec)]
v = value.(vec)
e = err.(vec)
#v = value(cpp)
#e = err(cpp)
figure()
errorbar(1:length(R_dat), value.(R_dat), err.(R_dat), fmt="x", color="black")
fill_between(x, v-e, v+e, color="gray", alpha=0.75)
ylabel(L"$R_\mathrm{PS}$")
xlabel(L"$x_0/a$")
#ylim(1e-8, 1e-7)
#xlim(0,length(R_dat))
yscale("log")
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/R_",ens.id,"_",PS,"_plat.pdf"))
#@gp x value.(R_dat) err.(R_dat) "w errorbars t ''"
#@gp:- x v-e v+e "w filledcurves t ''"
#@gp:- "set logscale y"
#save(string("/home/asaez/cls_ens/codes/lattA.jl/plots/R_",ens.id,"_",PS,"_plat.gp"))
end
return f, syst, f_i, weight, pval
end
function get_w0t0_plat(path::String, ens::EnsInfo, plat1::Vector{Int64}, plat2::Vector{Int64};
pl::Bool=false,
rw=false, npol::Int64=2, ws::ADerrors.wspace=ADerrors.wsg,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,
w0_guess::Union{Float64, Nothing}=nothing,
tau::Union{Float64, Nothing}=nothing)
y0 = 1 ## assumes this is the case, hardcoded, some ensembles will not fulfil !
println("WARNING!: make sure t_src is 1 in this ensemble")
t2YM, tdt2YM, W_obs, t, t_aux = get_YM_dYM(path, ens, rw=rw, w0=w0_guess, tau=tau)
dt0 = iseven(npol) ? Int64(npol / 2) : Int64((npol+1) / 2)
model(x, p) = get_model(x, p, npol)
fmin(x, p) = model(x, p) .- 0.3
function fit_defs(x,W) ## uncorrelated fit
chisq(p,d) = sum((d .- [p[1] for i in 1:length(x)]).^2 .* W)
return chisq
end
## w0
tdt2YM_guess = [plat_av(tdt2YM[:,i], plat1) for i in 1:length(tdt2YM[1,:])]
nt0 = findmin(abs.(value.(tdt2YM_guess[2:end-1]) .- 0.3))[2] + 1
x = t_aux[nt0-dt0:nt0+dt0]
xmax = size(tdt2YM, 1)
T = xmax - 1 - y0
tdt2E_i = Array{uwreal,1}()
syst_i = Array{uwreal,1}()
for j in 1:2*dt0+1
i = Int(round(dt0+0.5))
dat = tdt2YM[:,nt0-dt0-1+j]
if j == i
tdt2E_aux = plat_av(dat, plat1)
uwerr.(dat)
chisq = fit_defs(collect(plat1[1]:plat1[2]), 1 ./ err.(dat[plat1[1]:plat1[2]]) .^ 2)
chi2 = chisq(tdt2E_aux,dat[plat1[1]:plat1[2]])
isnothing(wpm) ? pval = pvalue(chisq,value(chi2),value.(tdt2E_aux),dat[plat1[1]:plat1[2]]) : pval = pvalue(chisq,value(chi2),value.(tdt2E_aux),dat[plat1[1]:plat1[2]],wpm=wpm)
if pl == true
isnothing(wpm) ? uwerr(tdt2E_aux) : uwerr(tdt2E_aux, wpm)
isnothing(wpm) ? uwerr.(dat) : [uwerr(dat[i], wpm) for i in 1:length(dat)]
v = value(tdt2E_aux)
e = err(tdt2E_aux)
fig = figure(string("pvalue=",pval))
errorbar(1:length(dat), value.(dat), err.(dat), fmt="x", color="black")
x_plot = collect(plat1[1]:plat1[2])
fill_between(x_plot, v-e, v+e, color="green", alpha=0.3)
ylabel(L"$t\frac{d}{dt}t^2E(x_0,t)$")
xlabel(L"$x_0$")
ylim(v-5*e, v+20*e)
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/tdt2E_",ens.id,"_plat.pdf"))
end
else
tdt2E_aux = plat_av(dat, plat1)
end
push!(tdt2E_i, tdt2E_aux)
end
if tdt2E_i[end] < 0.30 || tdt2E_i[1] > 0.30
println("WARNING!: extrapolating w0/a")
end
#uwerr.(tdt2E_i)
#uwerr.(syst_i)
#println("tdt2E_i = ", tdt2E_i)
#println("syst_i = ", syst_i)
par, aux = fit_alg(model, x[1:length(tdt2E_i)], tdt2E_i, npol, wpm=wpm)
w0 = root_error(fmin, t_aux[nt0], par)
if pl == true
v = value.(tdt2E_i)
e = err.(tdt2E_i)
uwerr(w0)
fig = figure("tdt2E_vs_t")
errorbar(x[1:length(tdt2E_i)], v, e, fmt="x")
errorbar(value(w0), 0.3, xerr=err(w0), fmt="x")
ylabel(L"$t\frac{d}{dt}t^2\left<E\right>$")
xlabel(L"$t/a^2$")
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/w0_",ens.id,".pdf"))
#close("all")
end
##
## t0
t2YM_guess = [plat_av(t2YM[:,i], plat2) for i in 1:length(t2YM[1,:])]
nt0 = findmin(abs.(value.(t2YM_guess[2:end-1]) .- 0.3))[2] + 1
x = t[nt0-dt0:nt0+dt0]
xmax = size(t2YM, 1)
T = xmax - 1 - y0
t2E_i = Array{uwreal,1}()
syst_i = Array{uwreal,1}()
for j in 1:2*dt0+1
i = Int(round(dt0+0.5))
dat = t2YM[:,nt0-dt0-1+j]
if j == i
t2E_aux = plat_av(dat, plat2)
uwerr.(dat)
chisq = fit_defs(collect(plat2[1]:plat2[2]), 1 ./ err.(dat[plat2[1]:plat2[2]]) .^ 2)
chi2 = chisq(t2E_aux,dat[plat2[1]:plat2[2]])
isnothing(wpm) ? pval = pvalue(chisq,value(chi2),value.(t2E_aux),dat[plat2[1]:plat2[2]]) : pval = pvalue(chisq,value(chi2),value.(t2E_aux),dat[plat2[1]:plat2[2]],wpm=wpm)
if pl == true
isnothing(wpm) ? uwerr(t2E_aux) : uwerr(t2E_aux, wpm)
isnothing(wpm) ? uwerr.(dat) : [uwerr(dat[i], wpm) for i in 1:length(dat)]
v = value(t2E_aux)
e = err(t2E_aux)
fig = figure(string("pvalue= ",pval))
errorbar(1:length(dat), value.(dat), err.(dat), fmt="x", color="black")
x_plot = collect(plat2[1]:plat2[2])
fill_between(x_plot, v-e, v+e, color="green", alpha=0.3)
ylabel(L"$t^2E(x_0,t)$")
xlabel(L"$x_0$")
ylim(v-5*e, v+20*e)
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/t2E_",ens.id,"_plat.pdf"))
end
else
t2E_aux = plat_av(dat, plat2)
end
push!(t2E_i, t2E_aux)
end
if t2E_i[end] < 0.30 || t2E_i[1] > 0.30
println("WARNING!: extrapolating t0/a")
end
#uwerr.(t2E_i)
#uwerr.(syst_i)
#println("t2E_i = ", t2E_i)
#println("syst_i = ", syst_i)
par, aux = fit_alg(model, x[1:length(t2E_i)], t2E_i, npol, wpm=wpm)
t0 = root_error(fmin, t[nt0], par)
if pl == true
v = value.(t2E_i)
e = err.(t2E_i)
uwerr(t0)
fig = figure("t2E_vs_t")
errorbar(x[1:length(t2E_i)], v, e, fmt="x")
errorbar(value(t0), 0.3, xerr=err(t0), fmt="x")
ylabel(L"$t^2\left<E\right>$")
xlabel(L"$t/a^2$")
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/t0_",ens.id,".pdf"))
#close("all")
end
##
return w0,t0
end