Skip to content

J
juobs
Commit 20ec1516 authored by Alessandro 's avatar Alessandro
Browse files
Options

read_mesons and read_mesons_TSM now support automatic detection of missing measurements 

t0_phys changed from Bruno et al to the value from internal analysis
parent f3ec800f
Hide whitespace changes
Inline Side-by-side
Showing with 179 additions and 41 deletions
constants/juobs_const.jl
View file @ 20ec1516
......@@ -21,8 +21,8 @@ const ZS_over_ZP_err = [83, 64, 48, 23, 25 ] .* 1e-4
# const b_values2 = [3.40, 3.46, 3.55, 3.70]
const t0_data = [2.86, 3.659, 5.164, 8.595, 14.040]
const t0_error = [11, 16, 18, 29, 49] .* 1e-3
const t0_ph_value = [0.413]
const t0_ph_error = ones(1,1) .* 5e-3
const t0_ph_value = [0.4118]
const t0_ph_error = ones(1,1) .* 25e-4
#2101.10969
const RM_data = [2.335, 1.869, 1.523, 1.267, 1.149]
const RM_err = [31, 19, 14, 16, 18] .* 1e-3
......
src/juobs_linalg.jl
View file @ 20ec1516
......@@ -183,7 +183,7 @@ Julia>
function getall_eigvals(a::Vector{Matrix{uwreal}}, t0::Int64; iter=5 )
n = length(a)
res = Vector{Vector{uwreal}}(undef, n)
[res[i] = uweigvals(a[i], a[t0]) for i=1:n]
[res[i] = uweigvals(a[i], a[t0], iter=iter) for i=1:n]
return res
end
......@@ -782,6 +782,77 @@ function more_penrose_inv(A::Matrix)
end
########################
# OPERATOR OVERLOADING
########################
function Base.:*(x::uwreal, y::Array{Any})
N = size(y, 1)
return fill(x, N) .* y
end
Base.:*(x::Array{Any}, y::uwreal) = Base.:*(y,x)
function Base.:*(x::uwreal, y::Array{Float64})
N = size(y, 1)
return fill(x, N) .* y
end
Base.:*(x::Array{Float64}, y::uwreal) = Base.:*(y,x)
function Base.:*(x::uwreal, y::Array{uwreal})
N = size(y, 1)
return fill(x, N) .* y
end
Base.:*(x::Array{uwreal}, y::uwreal) = Base.:*(y,x)
function Base.:+(x::uwreal, y::Vector{uwreal})
N = size(y, 1)
return fill(x, N) .+ y
end
function Base.:+(x::Float64, y::Vector{Float64})
N = size(y, 1)
return fill(x, N) .+ y
end
function Base.:+(x::Float64, y::Vector{uwreal})
N = size(y, 1)
return fill(x, N) .+ y
end
function Base.:-(x::Float64, y::Vector{Any})
N = size(y, 1)
return fill(x, N) .- y
end
function Base.:-(x::Float64, y::Vector{Float64})
N = size(y, 1)
return fill(x, N) .- y
end
function Base.:-(x::Float64, y::Vector{uwreal})
N = size(y, 1)
return fill(x, N) .- y
end
function Base.length(uw::uwreal)
return length(value(uw))
end
function Base.iterate(uw::uwreal)
return iterate(uw, 1)
end
function Base.iterate(uw::uwreal, state::Int64)
if state > length(uw)
return nothing
else
return uw[state], state +1
end
end
function Base.iterate(uw::Vector{uwreal})
return iterate(uw, 1)
end
function Base.iterate(uw::Vector{uwreal}, state::Int64)
if state > length(uw)
return nothing
else
return uw[state], state +1
end
end
function Base.getindex(uw::uwreal, ii::Int64)
idx = getindex(value(uw), ii)
return uw # uwreal([getindex(value(uw), kwargs...), err(uw)], " ")
end
"""
Base.:*(x::uwreal, y::Matrix{uwreal})
......
src/juobs_obs.jl
View file @ 20ec1516
......@@ -33,9 +33,12 @@ function meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bo
errorbar(x, y, dy, fmt="x", color="black")
ylabel(L"$m_\mathrm{eff}$")
xlabel(L"$x_0$")
ylim(v-30*e, v+30*e)
# ylim(0.05, 0.075)
xlim(64,length(y))
_, max_idx = findmax(value.(corr))
ylim(v-20*e, v+20*e)
xlim(left=max_idx)
# ylim(v-30*e, v+30*e)
# ylim(0.15, 0.18)
#xlim(64,length(y))
if !isnothing(kappa)
title(string(L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
......@@ -51,7 +54,6 @@ function meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bo
return (mass, aux)
end
end
function meff(corr::Corr, plat::Vector{Int64}; pl::Bool=true, data::Bool=false,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
......@@ -123,7 +125,9 @@ function mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca:
errorbar(x, y, dy, fmt="x", color="black")
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]))
......@@ -352,7 +356,7 @@ function dec_const(vv::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64
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()
xlim(y0, y0 +plat[2] +15)
xlim(left=y0)
ylim(v-10*e, v+10*e)
ylabel(L"$af_{B}$")
#ylabel(L"$R_\mathrm{av}$")
......@@ -617,6 +621,7 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false,
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
......
src/juobs_reader.jl
View file @ 20ec1516
......@@ -412,6 +412,8 @@ function read_ms(path::String; id::Union{String, Nothing}=nothing, dtr::Int64=1
ntr = Int32((fsize - 3*4 - 8) / datsize)
if mod(ntr, dtr) != 0
println("ntr = ", ntr)
println("dtr = ", dtr)
error("ntr / dtr must be exact")
end
......
src/juobs_tools.jl
View file @ 20ec1516
......@@ -95,11 +95,16 @@ function corr_obs(cdata::CData; real::Bool=true, rw::Union{Array{Float64, 2}, No
real ? data = cdata.re_data ./ L^3 : data = cdata.im_data ./ L^3
nt = size(data)[2]
nms_vec = cdata.vcfg
delta_cnfg = nms_vec[2] - nms_vec[1]
idm = isnothing(idm) ? collect(1:delta_cnfg:nms) : idm
if isnothing(rw)
idm = isnothing(idm) ? collect(1:nms) : idm
# idm = isnothing(idm) ? collect(1:nms) : idm
obs = [uwreal(data[:, x0], cdata.id, idm, nms) for x0 = 1:nt]
else
idm = isnothing(idm) ? collect(1:nms) : idm
# idm = isnothing(idm) ? collect(1:nms) : idm
data_r, W = apply_rw(data, rw)
ow = [uwreal(data_r[:, x0], cdata.id, idm, nms) for x0 = 1:nt]
W_obs = uwreal(W, cdata.id, idm, nms)
......@@ -116,11 +121,22 @@ function corr_obs(cdata::CData; real::Bool=true, rw::Union{Array{Float64, 2}, No
end
#function corr_obs for R != 1
#TODO: vcfg with gaps
function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array{Float64, 2}, 1}, Nothing}=nothing, L::Int64=1, info::Bool=false)
function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array{Float64, 2}, 1}, Nothing}=nothing, L::Int64=1, info::Bool=false, idm::Union{Vector{Int64},Nothing}=nothing, nms::Union{Int64, Nothing}=nothing)
nr = length(cdata)
id = getfield.(cdata, :id)
vcfg = getfield.(cdata, :vcfg)
replica = Int64.(maximum.(vcfg))
nms = isnothing(nms) ? sum(replica) : nms
if isnothing(idm)
delta_vcfg = vcfg[1][2] - vcfg[1][1]
idm = collect(1:delta_vcfg:replica[1])
for ii in eachindex(cdata)[2:end]
delta_vcfg = vcfg[ii][2] - vcfg[ii][1]
aux = collect(1:delta_vcfg:replica[ii]) .+ replica[ii-1]
append!(idm, aux)
end
end
if !all(id .== id[1])
error("IDs are not equal")
......@@ -133,7 +149,7 @@ function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array
if isnothing(rw)
tmp = data[1]
[tmp = cat(tmp, data[k], dims=1) for k = 2:nr]
obs = [uwreal(tmp[:, x0], id[1], replica) for x0 = 1:nt]
obs = [uwreal(tmp[:, x0], id[1], replica, idm, nms) for x0 = 1:nt]
else
data_r, W = apply_rw(data, rw)
tmp = data_r[1]
......@@ -141,8 +157,8 @@ function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array
[tmp = cat(tmp, data_r[k], dims=1) for k = 2:nr]
[tmp_W = cat(tmp_W, W[k], dims=1) for k = 2:nr]
ow = [uwreal(tmp[:, x0], id[1], replica) for x0 = 1:nt]
W_obs = uwreal(tmp_W, id[1], replica)
ow = [uwreal(tmp[:, x0], id[1], replica, idm, nms) for x0 = 1:nt]
W_obs = uwreal(tmp_W, id[1], replica, idm, nms)
obs = [ow[x0] / W_obs for x0 = 1:nt]
end
if info && !isnothing(rw)
......@@ -204,7 +220,7 @@ corr_pp_r = corr_obs.(data, rw=[rw1, rw2])
"""
function corr_obs_TSM(cdata_sl::CData, cdata_corr::CData; real::Bool=true, rw::Union{Array{Float64, 2}, Nothing}=nothing, L::Int64=1, info::Bool=false, idm_sl::Union{Vector{Int64},Nothing}=nothing, idm_corr::Union{Vector{Int64},Nothing}=nothing, nms::Int64=Int64(maximum(cdata_sl.vcfg)))
# cdata1 is sloppy, cdata2 is correction
if cdata_sl.id != cdata_corr.id[1:4]
if cdata_sl.id != cdata_corr.id
error("Error: cdata_sl id != cdata_corr id")
end
if cdata_sl.header != cdata_corr.header # Base.:(==) and Base.:(!=) are redifined in juobs_types.jl
......@@ -218,7 +234,16 @@ function corr_obs_TSM(cdata_sl::CData, cdata_corr::CData; real::Bool=true, rw::U
nt2 = size(data2, 2)
id_sl = cdata_sl.id
id_corr = cdata_corr.id
# println(id_sl, " ", id_corr)
# added automatic idsm detection - to be tested
nms_vec_sl = cdata_sl.vcfg
nms_vec_corr = cdata_corr.vcfg
delta_cnfg_sl = nms_vec_sl[2] - nms_vec_sl[1]
delta_cnfg_corr = nms_vec_corr[2] - nms_vec_corr[1]
idm_sl = isnothing(idm_sl) ? collect(1:delta_cnfg_sl:nms_vec_sl[end]) : idm_sl
idm_corr = isnothing(idm_corr) ? collect(1:delta_cnfg_corr:nms_vec_corr[end]) : idm_corr
# end section to be tested
if isnothing(rw)
idm_sl = isnothing(idm_sl) ? collect(1:nms) : idm_sl
idm_corr = isnothing(idm_corr) ? collect(1:nms) : idm_corr
......@@ -250,7 +275,7 @@ function corr_obs_TSM(cdata_sl::CData, cdata_corr::CData; real::Bool=true, rw::U
end
end
function corr_obs_TSM(cdata1::Array{CData, 1}, cdata2::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array{Float64, 2}, 1}, Nothing}=nothing, L::Int64=1, info::Bool=false, idm::Union{Vector{Int64},Nothing}=nothing, nms::Union{Int64, Nothing}=nothing)
function corr_obs_TSM(cdata1::Array{CData, 1}, cdata2::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array{Float64, 2}, 1}, Nothing}=nothing, L::Int64=1, info::Bool=false, idm_sl::Union{Vector{Int64},Nothing}=nothing, idm_corr::Union{Vector{Int64},Nothing}=nothing, nms::Union{Int64, Nothing}=nothing)
if any(getfield.(cdata1, :id) .!= getfield.(cdata2, :id))
error("Error: cdata1 id != cdata2 id")
end
......@@ -259,13 +284,36 @@ function corr_obs_TSM(cdata1::Array{CData, 1}, cdata2::Array{CData, 1}; real::Bo
end
id = getfield.(cdata1, :id)
vcfg = getfield.(cdata1, :vcfg)
replica = Int64.(maximum.(vcfg))
if !all(id .== id[1])
error("IDs are not equal")
end
vcfg_sl = getfield.(cdata1, :vcfg)
replica_sl = Int64.(maximum.(vcfg_sl))
vcfg_corr = getfield.(cdata2, :vcfg)
replica_corr = Int64.(maximum.(vcfg_corr))
nms = isnothing(nms) ? sum(replica_sl) : nms # assuming vcfg_sl >= vcfg_corr
if isnothing(idm_sl)
delta_vcfg = vcfg_sl[1][2] - vcfg_sl[1][1]
idm_sl = collect(1:delta_vcfg:replica_sl[1])
for ii in eachindex(cdata1)[2:end]
delta_vcfg = vcfg_sl[ii][2] - vcfg_sl[ii][1]
aux = collect(1:delta_vcfg:replica_sl[ii]) .+ replica_sl[ii-1]
append!(idm_sl, aux)
end
end
if isnothing(idm_corr)
delta_vcfg = vcfg_corr[1][2] - vcfg_corr[1][1]
idm_corr = collect(1:delta_vcfg:replica_corr[1])
for ii in eachindex(cdata2)[2:end]
delta_vcfg = vcfg_corr[ii][2] - vcfg_corr[ii][1]
aux = collect(1:delta_vcfg:replica_corr[ii]) .+ replica_corr[ii-1]
append!(idm_corr, aux)
end
end
data1 = real ? getfield.(cdata1, :re_data) ./ L^3 : getfield.(cdata1, :im_data) ./ L^3
data2 = real ? getfield.(cdata2, :re_data) ./ L^3 : getfield.(cdata2, :im_data) ./ L^3
......@@ -275,8 +323,8 @@ function corr_obs_TSM(cdata1::Array{CData, 1}, cdata2::Array{CData, 1}; real::Bo
tmp1 = cat(data1..., dims=1)
tmp2 = cat(data2..., dims=1)
obs1 = [uwreal(tmp1[:, x0], id[1], replica) for x0 = 1:nt]
obs2 = [uwreal(tmp2[:, x0], id[1], replica) for x0 = 1:nt]
obs1 = [uwreal(tmp1[:, x0], id[1], replica_sl, idm_sl, nms) for x0 = 1:nt]
obs2 = [uwreal(tmp2[:, x0], id[1], replica_corr, idm_corr, nms) for x0 = 1:nt]
else
data1_r, W = apply_rw(data1, rw)
......@@ -286,10 +334,10 @@ function corr_obs_TSM(cdata1::Array{CData, 1}, cdata2::Array{CData, 1}; real::Bo
tmp2 = cat(data2_r..., dims=1)
tmp_W = cat(W..., dims=1)
ow1 = [uwreal(tmp1[:, x0], id[1], replica) for x0 = 1:nt]
ow2 = [uwreal(tmp2[:, x0], id[1], replica) for x0 = 1:nt]
ow1 = [uwreal(tmp1[:, x0], id[1], replica_sl, idm_sl, nms) for x0 = 1:nt]
ow2 = [uwreal(tmp2[:, x0], id[1], replica_corr, idm_corr, nms) for x0 = 1:nt]
W_obs = uwreal(tmp_W, id[1], replica)
W_obs = uwreal(tmp_W, id[1], replica_sl, idm_sl, nms )
obs1 = [ow1[x0] / W_obs for x0 = 1:nt]
obs2 = [ow2[x0] / W_obs for x0 = 1:nt]
......@@ -401,7 +449,7 @@ function md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ow::uwreal, w::Union{uwr
ivrep = ivrep[1]
if length(md) != length(ivrep)
warning("Nr obs != Nr md")
@info("Nr obs != Nr md")
end
#md_aux as a Matrix + Automatic truncation
......@@ -597,7 +645,7 @@ function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64},
for INDEX in tmin_array ## vary tmin
for j in tmax_array ## vary tmax
try
try
x = [i for i in INDEX+1:1:j]
yy = y[INDEX+1:1:j]
Ncut = total - length(x)
......@@ -614,9 +662,13 @@ function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64},
# if chi2 / dof(fit) > 5. || abs(value(up[1])) > 2. || err(up[1])/value(up[1])*100 >= 10.
# error()
# end
# println(up[1], " ", up[2], " ", up[3])
push!(Pval, Q)
push!(AIC, chi2 + 2*k + 2*Ncut)
push!(chi2chi2exp, chi2 / dof(fit))
# push!(AIC, chi2 + 2*k + 2*Ncut) # AIC criteria
push!(AIC, chi2 - 2*chi_exp) # TIC criteria
push!(chi2chi2exp, chi2 / dof(fit)) #substitute chi_exp with dof(fit) for chi^2/dof
push!(p1, up[1])
push!(mods,string("[", INDEX+1, ",", j, "]"))
# println("chi2_vs_exp " , chi2 / dof(fit), " chi2/dof ", sum(fit.resid.^2)/dof(fit))
......@@ -653,7 +705,9 @@ function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64},
p1_mean = sum(p1 .* weight_model)
isnothing(wpm) ? uwerr(p1_mean) : uwerr(p1_mean, wpm)
systematic_err = sqrt(sum(p1 .^ 2 .* weight_model) - (sum(p1 .* weight_model)) ^ 2)
# println(sum(p1 .^ 2 .* weight_model))
# println((sum(p1 .* weight_model)) ^ 2)
systematic_err = 0.0 #sqrt(sum(p1 .^ 2 .* weight_model) - (sum(p1 .* weight_model)) ^ 2)
if pl
fig = figure(figsize=(10,7.5))
......@@ -672,7 +726,7 @@ function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64},
fill_between(x, v-e, v+e, color="royalblue", alpha=0.2)
errorbar(x, y, dy, fmt="^", mfc="none", color="navy", capsize=2)
setp(ax1.get_xticklabels(),visible=false) # Disable x tick labels
isnothing(label) ? ylabel(L"$p_1$") : ylabel(label)
isnothing(label) ? ylabel(L"$p_1$") : ylabel(label)
subplot(412)
ax2=gca()
......@@ -692,7 +746,7 @@ function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64},
subplot(414)
bar(x, chi2chi2exp, alpha=0.4, color="forestgreen", edgecolor="darkgreen", linewidth=1.5 )
# ylabel(L"$\chi^2/\chi^2_{\mathrm{exp}}$")
ylabel(L"$\chi^2/dof$")
ylabel(L"$\chi^2/\mathrm{d.o.f.}$")
xticks(x, mods, rotation=45)
xlabel(L"$Models$")
......@@ -1054,11 +1108,12 @@ function fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}
(upar, chi_exp) = isnothing(wpm) ? fit_error(chisq, coef(fit), ydata, W=covar_inv) : fit_error(chisq, coef(fit), ydata, wpm, W=covar_inv)
end
chi2_fit_res = sum(fit.resid.^2 )
println("\n")
for i in 1:length(fit.resid)
println((fit.resid[i])^2)
end
println("\n")
# compute and print single point contribution to chi2
# println("\n")
# for i in 1:length(fit.resid)
# println((fit.resid[i])^2)
# end
# println("\n")
println("chi2 from fit residual = ", chi2_fit_res)
println("chi2 from chi2 function = ", chisq(coef(fit), ydata))
......@@ -1078,7 +1133,8 @@ function fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}
#println("Chisq / chiexp: ", chisq(coef(fit), ydata), " / ", chi_exp, " (dof: ", length(yval) - param,")")
println("Chisq / chiexp: ", chi2_fit_res, " / ", chi_exp, " (dof: ", length(yval) - param,")")
println("Chisq corrected: ", chis2_corrected)
return upar, chi2_fit_res / chi_exp
println("Return: params, chi2, chi2exp")
return upar, chi2_fit_res, chi_exp
end
function fit_routine(model::Vector{Function}, xdata::Vector{Array{Float64, N}} where N, ydata::Vector{Array{uwreal, N}} where N, param::Int64; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,
......@@ -1097,7 +1153,11 @@ function fit_routine(model::Vector{Function}, xdata::Vector{Array{Float64, N}} w
if isnothing(wpm)
uwerr.(ydata[i])
else
[uwerr(yaux, wpm) for yaux in ydata[i]]
# println(yaux for yaux in ydata[i])
for yaux in ydata[i]
[uwerr(yaux[k], wpm) for k in 1:length(yaux)]
end
# [[uwerr(yaux[k], wpm) for k in eachindex(yaux)] for yaux in ydata[i]]
end
e[i] = err.(ydata[i])
......@@ -1148,8 +1208,8 @@ function fit_routine(model::Vector{Function}, xdata::Vector{Array{Float64, N}} w
chis2_corrected = (length(dat) - param) * min_fun(sol.minimizer) / chi2_exp
println("Chisq corrected: ", chis2_corrected)
@info("fit_routine is returning par, chi2 / chi2_exp\n")
return upar, min_fun(sol.minimizer) / chi2_exp
@info("Returns: params, chi2, chi2_exp\n")
return upar, min_fun(sol.minimizer) , chi2_exp
end
function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}, param::Int64=3;
......
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