Skip to content

J
juobs
Commit 0e6b14b7 authored by Antonino D'Anna's avatar Antonino D'Anna
Browse files
Options

Updated mpcac function using plot_data and plot_func, added overload for... 

Updated mpcac function using plot_data and plot_func, added overload for Vector{Corr}, added optional return to plot obj, updated documentation
parent 9dab5fc5
Hide whitespace changes
Inline Side-by-side
Showing with 222 additions and 270 deletions
src/juobs_obs.jl
View file @ 0e6b14b7
......@@ -17,7 +17,7 @@ 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, kappa::Union{Vector{Float64}, Nothing}=nothing,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, savepl::Bool=false,
label::Union{String, Nothing}=nothing)
label::Union{String, Nothing}=nothing,)
dim = length(corr)
aux = 0.5 .* log.((corr[2:dim-2] ./ corr[3:dim-1]).^2)
mass = plat_av(aux, plat, wpm)
......@@ -30,19 +30,19 @@ function meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bo
e = err(mass)
filename = "meff/";
figure()
fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75)
l,=errorbar(x, y, dy, fmt="x", color="black")
fig,ax = subplots(1)
ax.fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75)
l,=ax.errorbar(x, y, dy, fmt="x", color="black")
if !isnothing(label)
l.set_label(label)
legend()
ax.legend()
filename*=label*"_"
end
ylabel(L"$m_\mathrm{eff}$")
xlabel(L"$x_0$")
ax.set_ylabel(L"$m_\mathrm{eff}$")
ax.set_xlabel(L"$x_0$")
_, max_idx = findmax(value.(corr))
ylim(v-80*e, v+80*e)
xlim(left=max_idx)
ax.set_ylim(v-20*e, v+20*e)
ax.set_xlim(left=max_idx)
if !isnothing(kappa)
......@@ -54,19 +54,20 @@ function meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bo
filename*="mu_1_$(mu[1])_mu_2_$(mu[2]).png"
end
if savepl
tight_layout()
gcf().set_size_inches(10,5)
savefig(filename)
close(gcf())
fig.tight_layout()
fig.set_size_inches(10,5)
mkpath(dirname(filename));
fig.savefig(filename)
close(fig)
else
display(gcf())
end
# close("all")
end
if !data
return mass
else
return (mass, aux)
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...)
......@@ -78,113 +79,146 @@ function meff(corr::Corr, plat::Vector{Int64}; kwargs...)
end
end
function meff(corrs::Vector{Corr}, plats::Vector{Vector{Int64}}; data::Bool=false, kwargs...)
if !data
return meff.(corrs,plats;kwargs)
end
aux = meff.(corrs,plats; data=data, kwargs...)
return getfield.(aux,1), getfield.(aux,2);
function meff(pp::Vector{Corr},plat::Vector{Vector{Int64}}; kwargs...)
a = meff.(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
@doc raw"""
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,savepl::Bool=false, label::Union{String,Nothing}=nothing)
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::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,savepl::Bool=false, label::Union{String,Nothing}=nothing)
mpcac(a0p::Vector{Corr},pp::Vector{Corr},plat::Vector{Vector{Int64}}; kwargs...)
mpcac(a0p::Vector{Corr},pp::Vector{Corr},plat::Vector{Vector{int}}; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, savepl::Bool=false, label::Union{String,Nothing}=nothing)
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` variable allows to compute `mpcac` using
the improved axial current. The flag `savepl` allow to save automatically the plot, the variable `label` if given is used as lalel for the plot
and in the filename.
Correlator can be passed as an `Corr` struct or `Vector{uwreal}`.
In case of Vector{Corr} it is better to not broadcast, espectially 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
```@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 = mpcac(corr_a0p, corr_pp, [50, 60], pl=false)
The return structure is: mpcac[,data][,figs,axs]
p0 = 9.2056
p1 = -13.9847
g2 = 1.73410
ca = -0.006033 * g2 *( 1 + exp(p0 + p1/g2))
# 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_funct`](@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
m12 = mpcac(corr_a0p, corr_pp, [50, 60], pl=false, ca=ca, savepl=true, label="with_ca_improvement")
# 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=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, label::Union{String,Nothing}=nothing)
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
#der_a0p = (corr_a0p[3:end] .- corr_a0p[1:end-2]) / 2
#if ca != 0.0
# der2_pp = corr_pp[1:end-2] + corr_pp[3:end] - 2 * corr_pp[2:end-1]
# der_a0p = der_a0p + ca * der2_pp
#end
aux = der_a0p ./ ( 2. .* corr_pp[3:end-2])
mass = plat_av(aux, plat, wpm)
if pl || savepl
isnothing(wpm) ? uwerr(mass) : uwerr(mass, wpm)
x = 1:length(aux)
y = value.(aux)
dy = err.(aux)
v = value(mass)
e = err(mass)
figure()
fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75)
l, = errorbar(x, y, dy, fmt="x", color="black")
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...)
fig,ax = plot_func(x->mass,[plat[1]:plat[2]...],figs=(fig,ax);plot_kw...)
if !isnothing(kappa)
fig.suptitle(string(L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
end
if !isnothing(mu)
fig.suptitle(string(L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2]))
end
_, max_idx = findmax(value.(pp))
isnothing(wpm) ? uwerr(mass) : uwerr(mass,wpm)
v,e = value(mass), err(mass)
ax.set_ylim(v-20*e, v+20*e)
ax.set_xlim(left=max_idx)
if savepl
fig.tight_layout();
if filename ==""
filename = "mpcac/"
if !isnothing(label)
l.set_label(label)
legend()
filename *= label*"_"
end
ylabel(L"$m_\mathrm{PCAC}$")
xlabel(L"$x_0$")
_, max_idx = findmax(value.(pp))
ylim(v-20*e, v+20*e)
xlim(left=max_idx)
if !isnothing(kappa)
title(string(L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
global filename *= "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 *="mu_1_$(mu[1])_mu_2_$(mu[2]).png"
end
if savepl
tight_layout()
gcf().set_size_inches(10, 5)
savefig(filename, dpi=100)
savefig(filename)
close(gcf())
elseif pl
display(gcf())
end
end
if !data
return mass
else
return (mass, aux)
end
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... )
......@@ -200,137 +234,29 @@ function mpcac(a0p::Corr, pp::Corr, plat::Vector{Int64}; kwargs... )
end
end
function mpcac(a0p::Vector{Corr},pp::Vector{Corr},plat::Vector{Vector{Int64}}; data::Bool=false, kwargs...)
if !data
return mpcac.(a0p,pp,plat;kwargs...)
else
aux = mpcac.(a0p,pp,plat; data=true, kwargs...)
return (getfield.(aux,1),getfield.(aux,2))
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
#=
@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}; ca::Float64=0.0,ba_bp_z::uwreal=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, label::Union{String,Nothing}=nothing)
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 = (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()
filename="mpcac/"
fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75)
l,= errorbar(x, y, dy, fmt="x", color="black")
if !isnothing(label)
l.set_label(label)
legend()
filename*=label*"_"
end
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)
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]))
filename *="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]))
filename *="mu_1_$(mu[1])_mu_2_$(mu[2])_renomalized_improved.png"
end
if savepl
tight_layout()
gcf().set_size_inches(10,5)
savefig(filename)
end
if pl && !savepl
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}; ca::Float64=0.0,ba_bp_z::uwreal=uwreal(0.0), pl::Bool=true, data::Bool=false,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,savepl::Bool =false,label::Union{String,Nothing}=nothing)
if a0p.kappa == pp.kappa || a0p.mu == pp.mu
if a0p.mu == [0.0, 0.0]
return mpcac_ren(a0p.obs, pp.obs, plat, ca=ca,ba_bp_z=ba_bp_z, pl=pl, data=data, kappa=a0p.kappa, mu=nothing, wpm=wpm, savepl=savepl,label=label)
else
return mpcac_ren(a0p.obs, pp.obs, plat, ca=ca,ba_bp_z=ba_bp_z, pl=pl, data=data, mu=a0p.mu, kappa=nothing, wpm=wpm, savepl=savepl,label=label)
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)
......@@ -379,7 +305,7 @@ f_ratio = dec_const(corr_a0pL, corr_a0pR, corr_ppL, corr_ppR, [50, 60], m, pl=tr
```
"""
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)
kappa::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,savepl::Bool=false)
corr_a0p = -a0p[2:end-1]
corr_pp = pp[2:end-1]
......@@ -398,47 +324,73 @@ function dec_const(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64},
R_av = plat_av(R, plat, wpm)
f = sqrt(2) * sqrt(R_av^2) / sqrt(m)
if pl
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()
ylim(v-10*e, v+10*e)
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())
close()
end
if !data
return f
if pl || savepl
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)
fig,ax = subplots(1)
lbl = string(L"$af = $", sprint(show, f))
ax.fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75, label=L"$R$")
ax.errorbar(x, y, dy, fmt="x", color="black", label=lbl)
ax.legend()
ax.set_ylim(v-10*e, v+10*e)
ax.set_ylabel(L"$R_\mathrm{av}$")
ax.set_xlabel(L"$x_0$")
filename = "dec_const/"
if !isnothing(kappa)
fig.suptitle(string(L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
filename*=string(L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2],".png")
end
if savepl
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 (f, R)
return data ? (f, R) : f
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)
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, ca=ca, kappa=a0p.kappa, pl=pl, data=data, wpm=wpm)
return dec_const(a0p.obs, pp.obs, plat, m, a0p.y0; 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
# 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 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=true, data::Bool=false, kappa::Union{Vector{Float64}, Nothing}=nothing,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
......
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