Skip to content

J
juobs
Commit 1e233c9a authored by Javier's avatar Javier
Browse files
Options

wpm added 

functions that calls uwerr now includes wpm
parent 9d02c460
Hide whitespace changes
Inline Side-by-side
Showing with 55 additions and 26 deletions
src/juobs_obs.jl
View file @ 1e233c9a
...@@ -14,12 +14,13 @@ corr_pp = corr_obs.(data) ...@@ -14,12 +14,13 @@ corr_pp = corr_obs.(data)
m = meff(corr_pp[1], [50, 60], pl=false) 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,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
dim = length(corr) dim = length(corr)
aux = 0.5 .* log.((corr[2:dim-2] ./ corr[3:dim-1]).^2) aux = 0.5 .* log.((corr[2:dim-2] ./ corr[3:dim-1]).^2)
mass = plat_av(aux, plat) mass = plat_av(aux, plat, wpm)
uwerr(mass)
if pl == true if pl == true
isnothing(wpm) ? uwerr(mass) : uwerr(mass, wpm)
x = 1:length(aux) x = 1:length(aux)
y = value.(aux) y = value.(aux)
dy = err.(aux) dy = err.(aux)
...@@ -65,7 +66,8 @@ m = meff(corr_pp[1], [50, 60], pl=false) ...@@ -65,7 +66,8 @@ m = meff(corr_pp[1], [50, 60], pl=false)
f = dec_const_pcvc(corr_pp[1], [50, 60], m, pl=false) f = dec_const_pcvc(corr_pp[1], [50, 60], m, pl=false)
``` ```
""" """
function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false) 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 compute the decay constant when the source is far from the boundaries
""" """
...@@ -74,14 +76,21 @@ function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu ...@@ -74,14 +76,21 @@ function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu
aux = exp.((collect(1:dim) .- y0 ) .* [m for k in 1:dim]) aux = exp.((collect(1:dim) .- y0 ) .* [m for k in 1:dim])
R = ((aux .* corr_pp).^2).^0.25 R = ((aux .* corr_pp).^2).^0.25
R_av = plat_av(R, plat) R_av = plat_av(R, plat, wpm)
f = sqrt(2) * (mu[1] + mu[2]) *R_av / m^1.5 f = sqrt(2) * (mu[1] + mu[2]) *R_av / m^1.5
uwerr(f)
if pl == true if pl == true
uwerr(R_av) if isnothing(wpm)
uwerr(f)
uwerr(R_av)
uwerr.(R)
else
uwerr(f, wpm)
uwerr(R_av, wpm)
uwerr.(R, wpm)
end
v = value(R_av) v = value(R_av)
e = err(R_av) e = err(R_av)
uwerr.(R)
figure() figure()
fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75) fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75)
errorbar(1:length(R), value.(R), err.(R), fmt="x", color="black") errorbar(1:length(R), value.(R), err.(R), fmt="x", color="black")
...@@ -139,7 +148,8 @@ t0_r = comp_t0(Y, [38, 58], L=32, rw=[rw1, rw2], pl=true) ...@@ -139,7 +148,8 @@ 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, 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) 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)
Ysl = Y.Ysl Ysl = Y.Ysl
t = Y.t t = Y.t
...@@ -177,7 +187,7 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, ...@@ -177,7 +187,7 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false,
end end
end end
x = t[nt0-dt0:nt0+dt0] x = t[nt0-dt0:nt0+dt0]
t2E = [plat_av(Y_aux[:, j], plat) for j=1:2*dt0+1] .* x.^2 / L^3 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) model(x, p) = get_model(x, p, npol)
...@@ -185,8 +195,14 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, ...@@ -185,8 +195,14 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false,
fmin(x, p) = model(x, p) .- 0.3 fmin(x, p) = model(x, p) .- 0.3
t0 = root_error(fmin, t[nt0], par) t0 = root_error(fmin, t[nt0], par)
if pl if pl
uwerr(t0) if isnothing(wpm)
uwerr.(t2E) uwerr(t0)
uwerr.(t2E)
else
uwerr(t0, wpm)
uwerr.(t2E, wpm)
end
v = value.(t2E) v = value.(t2E)
e = err.(t2E) e = err.(t2E)
...@@ -211,7 +227,8 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, ...@@ -211,7 +227,8 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false,
end end
function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false, 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) 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)
nr = length(Y) nr = length(Y)
Ysl = getfield.(Y, :Ysl) Ysl = getfield.(Y, :Ysl)
...@@ -258,7 +275,7 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false ...@@ -258,7 +275,7 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false
end end
end end
x = t[nt0-dt0:nt0+dt0] x = t[nt0-dt0:nt0+dt0]
t2E = [plat_av(Y_aux[:, j], plat) for j=1:2*dt0+1] .* x.^2 / L^3 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) model(x, p) = get_model(x, p, npol)
...@@ -266,8 +283,14 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false ...@@ -266,8 +283,14 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false
fmin(x, p) = model(x, p) .- 0.3 fmin(x, p) = model(x, p) .- 0.3
t0 = root_error(fmin, t[nt0], par) t0 = root_error(fmin, t[nt0], par)
if pl if pl
uwerr(t0) if isnothing(wpm)
uwerr.(t2E) uwerr(t0)
uwerr.(t2E)
else
uwerr(t0, wpm)
uwerr.(t2E, wpm)
end
v = value.(t2E) v = value.(t2E)
e = err.(t2E) e = err.(t2E)
......
src/juobs_tools.jl
View file @ 1e233c9a
...@@ -183,8 +183,8 @@ function md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ws::ADerrors.wspace=ADer ...@@ -183,8 +183,8 @@ function md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ws::ADerrors.wspace=ADer
end end
end end
function plat_av(obs::Vector{uwreal}, plat::Vector{Int64}) function plat_av(obs::Vector{uwreal}, plat::Vector{Int64}, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing})
uwerr.(obs) isnothing(wpm) ? uwerr.(obs) : uwerr.(obs, wpm)
w = 1 ./ err.(obs)[plat[1]:plat[2]].^2 w = 1 ./ err.(obs)[plat[1]:plat[2]].^2
av = sum(w .* obs[plat[1]:plat[2]]) / sum(w) av = sum(w .* obs[plat[1]:plat[2]]) / sum(w)
return av return av
...@@ -212,13 +212,13 @@ fitp, csqexp = lin_fit(phi2, m2) ...@@ -212,13 +212,13 @@ fitp, csqexp = lin_fit(phi2, m2)
m2_phys = fitp[1] + fitp[2] * phi2_phys m2_phys = fitp[1] + fitp[2] * phi2_phys
``` ```
""" """
function lin_fit(x::Vector{<:Real}, y::Vector{uwreal}) function lin_fit(x::Vector{<:Real}, y::Vector{uwreal}; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
uwerr.(y) isnothing(wpm) ? uwerr.(y) : uwerr.(y, wpm)
par = lin_fit(x, value.(y), err.(y)) par = lin_fit(x, value.(y), err.(y))
chisq(p, d) = sum((d .- p[1] .- p[2].*x).^2 ./ err.(y) .^2) chisq(p, d) = sum((d .- p[1] .- p[2].*x).^2 ./ err.(y) .^2)
(fitp, csqexp) = fit_error(chisq, par, y) (fitp, csqexp) = fit_error(chisq, par, y)
for i = 1:length(fitp) for i = 1:length(fitp)
uwerr(fitp[i]) isnothing(wpm) ? uwerr(fitp[i]) : uwerr(fitp[i], wpm)
print("\n Fit parameter: ", i, ": ") print("\n Fit parameter: ", i, ": ")
details(fitp[i]) details(fitp[i])
end end
...@@ -252,14 +252,14 @@ The method return an array upar with the best fit parameters with their errors. ...@@ -252,14 +252,14 @@ The method return an array upar with the best fit parameters with their errors.
fit_routine(model, xdata, ydata, param=3) fit_routine(model, xdata, ydata, param=3)
""" """
function fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}, param::Int64=3; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) function fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}, param::Int64=3; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
uwerr.(ydata) isnothing(wpm) ? uwerr.(ydata) : uwerr.(ydata, wpm)
yval = value.(ydata) yval = value.(ydata)
yer = err.(ydata) yer = err.(ydata)
chisq = gen_chisq(model, xdata, yer) chisq = gen_chisq(model, xdata, yer)
fit = curve_fit(model, xdata, yval, 1.0 ./ yer.^2, fill(0.5, param)) fit = curve_fit(model, xdata, yval, 1.0 ./ yer.^2, fill(0.5, param))
(upar, chi_exp) = isnothing(wpm) ? fit_error(chisq, coef(fit), ydata) : fit_error(chisq, coef(fit), ydata, wpm) (upar, chi_exp) = isnothing(wpm) ? fit_error(chisq, coef(fit), ydata) : fit_error(chisq, coef(fit), ydata, wpm)
for i = 1:length(upar) for i = 1:length(upar)
uwerr(upar[i]) isnothing(wpm) ? uwerr(upar[i]) : uwerr(upar[i], wpm)
print("\n Fit parameter: ", i, ": ") print("\n Fit parameter: ", i, ": ")
details(upar[i]) details(upar[i])
end end
...@@ -275,10 +275,16 @@ end ...@@ -275,10 +275,16 @@ end
function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}, param::Int64=3; function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}, param::Int64=3;
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, covar::Bool=false) wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, covar::Bool=false)
uwerr.(ydata) if isnothing(wpm)
uwerr.(ydata)
uwerr.(xdata)
else
uwerr.(ydata, wpm)
uwerr.(xdata, wpm)
end
yval = value.(ydata) yval = value.(ydata)
yer = err.(ydata) yer = err.(ydata)
uwerr.(xdata)
xval = value.(xdata) xval = value.(xdata)
xer = err.(xdata) xer = err.(xdata)
...@@ -307,7 +313,7 @@ function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal} ...@@ -307,7 +313,7 @@ function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}
#### chisq_full, min_fun out of conditional -> #### chisq_full, min_fun out of conditional ->
#### COMPILER WARNING ** incremental compilation may be fatally broken for this module ** #### COMPILER WARNING ** incremental compilation may be fatally broken for this module **
for i = 1:length(upar) for i = 1:length(upar)
uwerr(upar[i]) isnothing(wpm) ? uwerr(upar[i]) : uwerr(upar[i], wpm)
print("\n Fit parameter: ", i, ": ") print("\n Fit parameter: ", i, ": ")
details(upar[i]) details(upar[i])
end end
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment