Commit 5bc786ac authored by Alberto Ramos's avatar Alberto Ramos

Documentation finished

parent f965e999
...@@ -259,8 +259,10 @@ function uwerror(a::uwreal, ws::wspace, wpm::Dict{Int64,Vector{Float64}}) ...@@ -259,8 +259,10 @@ function uwerror(a::uwreal, ws::wspace, wpm::Dict{Int64,Vector{Float64}})
gwin = view(a.cfd[j].gamm, 2:iw) gwin = view(a.cfd[j].gamm, 2:iw)
dbias = a.cfd[j].gamm[1] + 2.0*sum(gwin) dbias = a.cfd[j].gamm[1] + 2.0*sum(gwin)
a.cfd[j].gamm .= a.cfd[j].gamm .+ dbias/nd_eff if (dbias > 0.0)
a.cfd[j].gamm .= a.cfd[j].gamm .+ dbias/nd_eff
end
a.cfd[j].drho = zeros(nt) a.cfd[j].drho = zeros(nt)
if (a.cfd[j].gamm[1] != 0.0) if (a.cfd[j].gamm[1] != 0.0)
@inbounds for i in 1:nt @inbounds for i in 1:nt
...@@ -290,7 +292,7 @@ function uwerror(a::uwreal, ws::wspace, wpm::Dict{Int64,Vector{Float64}}) ...@@ -290,7 +292,7 @@ function uwerror(a::uwreal, ws::wspace, wpm::Dict{Int64,Vector{Float64}})
end end
end end
id_neg_taui = Vector{Int64}()
a.err = 0.0 a.err = 0.0
a.derr = 0.0 a.derr = 0.0
for j in 1:nid for j in 1:nid
...@@ -312,50 +314,66 @@ function uwerror(a::uwreal, ws::wspace, wpm::Dict{Int64,Vector{Float64}}) ...@@ -312,50 +314,66 @@ function uwerror(a::uwreal, ws::wspace, wpm::Dict{Int64,Vector{Float64}})
iw = 0 iw = 0
else else
wp = zeros(4) wp = zeros(4)
wp = get(wpm, a.ids[j], [-1.0,-1.0,-1.0,-1.0]) if haskey(wpm, a.ids[j])
if (wp[1] > 0.0) wp = get(wpm, a.ids[j], [-1.0,-1.0,-1.0,-1.0])
a.cfd[j].iw = round(wp[1]) if (wp[1] > 0.0)
elseif (wp[2] > 0.0) a.cfd[j].iw = round(wp[1])
a.cfd[j].iw = wopt_ulli(nd_eff, wp[2], a.cfd[j].gamm) elseif (wp[2] > 0.0)
end a.cfd[j].iw = wopt_ulli(nd_eff, wp[2], a.cfd[j].gamm)
end
if (wp[3] > 0.0)
iw = 1 if (wp[3] > 0.0)
for k in 2:nt iw = 1
if (a.cfd[j].drho[k]*wp[3] > a.cfd[j].gamm[k]/a.cfd[j].gamm[1]) for k in 2:nt
iw = k-1 if (a.cfd[j].drho[k]*wp[3] > a.cfd[j].gamm[k]/a.cfd[j].gamm[1])
break iw = k-1
break
end
end end
a.cfd[j].iw = iw
end
if (wp[4] > 0.0)
texp = wp[4]/ibn
else
texp = 0.0
end end
a.cfd[j].iw = iw
end
if (wp[4] > 0.0)
texp = wp[4]/ibn
else else
texp = 0.0 texp = 0.0
a.cfd[j].iw = wopt_ulli(nd_eff, DEFAULT_STAU, a.cfd[j].gamm)
end end
iw = a.cfd[j].iw iw = a.cfd[j].iw
gwin = view(a.cfd[j].gamm, 2:iw) gwin = view(a.cfd[j].gamm, 2:iw)
vti = 0.5 + sum(gwin)/a.cfd[j].gamm[1] vti = 0.5 + sum(gwin)/a.cfd[j].gamm[1]
a.cfd[j].dtaui = sqrt(vti^2 * (4.0*iw-2.0*vti+2.0)/nd_eff) a.cfd[j].dtaui = sqrt(vti^2 * (4.0*iw-2.0*vti+2.0)/nd_eff)
a.cfd[j].taui = vti + texp*a.cfd[j].gamm[iw+1]/a.cfd[j].gamm[1] a.cfd[j].taui = vti + texp*a.cfd[j].gamm[iw+1]/a.cfd[j].gamm[1]
end end
a.cfd[j].var = a.cfd[j].gamm[1] * 2.0*a.cfd[j].taui/nd_eff if (a.cfd[j].taui > 0.0)
a.cfd[j].var = a.cfd[j].gamm[1] * 2.0*a.cfd[j].taui/nd_eff
else
push!(id_neg_taui, a.ids[j])
end
end end
a.err = a.err + a.cfd[j].var a.err = a.err + a.cfd[j].var
if (iw > 1) if (iw > 1)
a.derr = a.derr + vti*(a.cfd[j].iw-0.5)/nd_eff if vti*(a.cfd[j].iw-0.5)/nd_eff > 0.0
a.derr = a.derr + vti*(a.cfd[j].iw-0.5)/nd_eff
end
end end
end end
a.err = sqrt(a.err)
a.derr = sqrt(a.derr)
return a.err if (length(id_neg_taui) == 0)
a.err = sqrt(a.err)
a.derr = sqrt(a.derr)
else
println(stderr, "ID's with negative tau_int: ", id_neg_taui)
error("Error analysis failed for some ID's. Choose your window more carefully")
end
return nothing
end end
function unique_ids_multi(a::Vector{uwreal}, ws::wspace) function unique_ids_multi(a::Vector{uwreal}, ws::wspace)
...@@ -682,12 +700,12 @@ uwreal(x::Float64) = ADerrors.uwreal(x, 0.0, 0.0, ...@@ -682,12 +700,12 @@ uwreal(x::Float64) = ADerrors.uwreal(x, 0.0, 0.0,
uwreal(data::Vector{Float64}, id::Int64) = ADerrors.uwcls(data::Vector{Float64}, id::Int64, wsg, [length(data)]) uwreal(data::Vector{Float64}, id::Int64) = ADerrors.uwcls(data::Vector{Float64}, id::Int64, wsg, [length(data)])
uwreal(data::Vector{Float64}, id::Int64, iv::Vector{Int64}) = ADerrors.uwcls(data::Vector{Float64}, id::Int64, wsg, iv) uwreal(data::Vector{Float64}, id::Int64, iv::Vector{Int64}) = ADerrors.uwcls(data::Vector{Float64}, id::Int64, wsg, iv)
uwreal(data::Vector{Float64}, id::Int64, idm::Vector{Int64}, nms::Int64) = ADerrors.uwcls_gaps(data::Vector{Float64}, uwreal(data::Vector{Float64}, id::Int64, idm::Vector{Int64}, nms::Int64) = ADerrors.uwcls_gaps(data::Vector{Float64},
id::Int64, ws::wspace, id::Int64, wsg,
[nms], [nms],
idm::Vector{Int64}, idm::Vector{Int64},
nms::Int64) nms::Int64)
uwreal(data::Vector{Float64}, id::Int64, iv::Vector{Int64}, idm::Vector{Int64}, nms::Int64) = ADerrors.uwcls_gaps(data::Vector{Float64}, uwreal(data::Vector{Float64}, id::Int64, iv::Vector{Int64}, idm::Vector{Int64}, nms::Int64) = ADerrors.uwcls_gaps(data::Vector{Float64},
id::Int64, ws::wspace, id::Int64, wsg,
iv::Vector{Int64}, iv::Vector{Int64},
idm::Vector{Int64}, idm::Vector{Int64},
nms::Int64) nms::Int64)
......
...@@ -85,7 +85,6 @@ uwerr(a) ...@@ -85,7 +85,6 @@ uwerr(a)
println("Error analysis result: ", a, " (tauint = ", taui(a, 666), ")") println("Error analysis result: ", a, " (tauint = ", taui(a, 666), ")")
``` ```
""" """
function taui(a::uwreal, mcid::Int64) function taui(a::uwreal, mcid::Int64)
idx = find_mcid(a, mcid) idx = find_mcid(a, mcid)
if (idx == nothing) if (idx == nothing)
...@@ -118,7 +117,6 @@ println("Error analysis result: ", a, ...@@ -118,7 +117,6 @@ println("Error analysis result: ", a,
" (tauint = ", taui(a, 666), " +/- ", dtaui(a, 666), ")") " (tauint = ", taui(a, 666), " +/- ", dtaui(a, 666), ")")
``` ```
""" """
function dtaui(a::uwreal, mcid::Int64) function dtaui(a::uwreal, mcid::Int64)
idx = find_mcid(a, mcid) idx = find_mcid(a, mcid)
if (idx == nothing) if (idx == nothing)
...@@ -221,6 +219,7 @@ end ...@@ -221,6 +219,7 @@ end
``` ```
""" """
function drho(a::uwreal, mcid::Int64) function drho(a::uwreal, mcid::Int64)
idx = find_mcid(a, mcid)
if (idx == nothing) if (idx == nothing)
error("No error available... maybe run uwerr") error("No error available... maybe run uwerr")
else else
......
...@@ -149,8 +149,6 @@ xp = [sum(value.(dt) ./ dx)/sum(1.0 ./ dx)] ...@@ -149,8 +149,6 @@ xp = [sum(value.(dt) ./ dx)/sum(1.0 ./ dx)]
# Compare chi^2 and expected chi^2 # Compare chi^2 and expected chi^2
println("chi^2 / chi_exp^2: ", chisq(xp, value.(dt)), " / ", chiexp(chisq, xp, dt)) println("chi^2 / chi_exp^2: ", chisq(xp, value.(dt)), " / ", chiexp(chisq, xp, dt))
```
""" """
function chiexp(chisq::Function, function chiexp(chisq::Function,
xp::Vector{Float64}, xp::Vector{Float64},
...@@ -253,6 +251,7 @@ xp = [sum(value.(dt) ./ dx)/sum(1.0 ./ dx)] ...@@ -253,6 +251,7 @@ xp = [sum(value.(dt) ./ dx)/sum(1.0 ./ dx)]
uwerr.(fitp) uwerr.(fitp)
println("Fit parameter: ", fitp[1]) println("Fit parameter: ", fitp[1])
println("chi^2 / chi_exp^2: ", chisq(xp, value.(dt)), " / ", csqexp) println("chi^2 / chi_exp^2: ", chisq(xp, value.(dt)), " / ", csqexp)
```
""" """
function fit_error(chisq::Function, function fit_error(chisq::Function,
xp::Vector{Float64}, xp::Vector{Float64},
......
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