Commit 17ce4bb1 authored by Javier's avatar Javier

uwLinAlg deleted and find_xmin commented

find_xmin to be tested (issues with double exponential fits)
parent 6a9891c7
......@@ -141,3 +141,52 @@ Computes the results of a linear interpolation/extrapolation in the y axis
y_lin_fit(par::Vector{uwreal}, x::Union{uwreal, Float64}) = par[1] + par[2] * x
#TODO: add combined fits
#=
using LsqFit
@doc raw"""
find_xmin(obs::Vector{uwreal}, y0::Int64; pl::Bool=false)
find_xmin(corr::Corr; pl::Bool=false)
Find the platau starting point for a given correlator.
The starting point xmin is determined
|C_2|^2 * exp(-M_2 * (xmin-y0)) / (2 * M_2) < dC(xmin, y0) / 4
where C_2 and M_2 are the matrix element and mass of the first excited state and dC is
the statistical error of the correlator
"""
function find_xmin(obs::Vector{uwreal}, y0::Int64; pl::Bool=false)
p0 = [0.5, 0.5, 1.0, 1.0]
@. model(t, p) = p[1] * exp(-p[2] * (t-y0)) + p[3] * exp(-p[4] * (t-y0))
x = Int64.(0:length(obs)-1)
uwerr.(obs)
y = value.(obs)
dy = err.(obs)
wt = 1 ./ dy.^2
fit = curve_fit(model, Float64.(x[y0+1:end]), y[y0+1:end], p0)
fit = curve_fit(model, Float64.(x[y0+1:end]), y[y0+1:end], wt[y0+1:end], fit.param)
par = fit.param
#println(par)
if par[2] > par[4]
p1 = par[1]
p2 = par[2]
else
p1 = par[3]
p2 = par[4]
end
@. f(t) = p1^2 * exp(-p2 * (t-y0)) / (2 * p2)
xmin = findfirst(t-> f(t) < 0.25*dy[t + 1], x)
if pl
errorbar(x, y, dy, fmt="x")
t = Int64.(y0:length(obs))
plot(t, model(t, par))
display(gcf())
end
return xmin
end
find_xmin(corr::Corr; pl::Bool=false) = find_xmin(corr.obs, corr.y0, pl=pl)
=#
\ No newline at end of file
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