comp_t0 polynomial fit

comp_t0 is generalized for polynomials of degree npol

src/juobs_obs.jl

@@ -100,15 +100,22 @@ dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::
     dec_const_pcvc(corr.obs, plat, m, corr.mu, corr.y0, pl=pl, data=data)
 
 #t0
+function get_model(x, p, n)
+    sum = 0.0
+    for k = 1:n
+        sum = sum .+ p[k] .* x.^(k-1)
+    end
+    return sum
+end
 @doc raw"""
 
-comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Matrix{Float64}, Nothing}=nothing)
+comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Matrix{Float64}, Nothing}=nothing, npol::Int64=2)
 
-comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Vector{Matrix{Float64}}, Nothing}=nothing)
+comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Vector{Matrix{Float64}}, Nothing}=nothing, npol::Int64=2)
 
 Computes t0 using the energy density of the action Ysl(Yang-Mills action).
 t0 is computed in the plateau plat.
-A (linear) interpolation in t is performed to find t0.
+A polynomial interpolation in t is performed to find t0, where npol is the degree of the polynomial (linear fit by default)
 
 The flag pl allows to show the plot.
 
@@ -131,7 +138,7 @@ 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, rw::Union{Matrix{Float64}, Nothing}=nothing)
+function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Matrix{Float64}, Nothing}=nothing, npol::Int64=2)
     Ysl = Y.Ysl
     t = Y.t
     id = Y.id
@@ -140,19 +147,21 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Un
     xmax = size(Ysl, 2)
     nt0 = t0_guess(t, Ysl, plat, L)
     
-    Y_aux = Matrix{uwreal}(undef, xmax, 5)
+    dt0 = iseven(npol) ? Int64(npol / 2) : Int64((npol+1)/ 2)
+    Y_aux = Matrix{uwreal}(undef, xmax, 2*dt0+1)
     for i = 1:xmax
         k = 1
-        for j = nt0-2:nt0+2
+        for j = nt0-dt0:nt0+dt0
             Y_aux[i, k] = uwreal(Ysl[:, i, j], id)
             k = k + 1
         end
     end
-    x = t[nt0-2:nt0+2]
-    t2E = [plat_av(Y_aux[:, j], plat) for j=1:5] .* x.^2 / L^3
+    x = t[nt0-dt0:nt0+dt0]
+    t2E = [plat_av(Y_aux[:, j], plat) for j=1:2*dt0+1] .* x.^2 / L^3
     
-    @. model(x, p) = p[1]*x^3 + p[2]*x^2 + p[3]*x + p[4]
-    par = fit_routine(model, x, t2E, 4)
+    model(x, p) = get_model(x, p, npol)
+
+    par = fit_routine(model, x, t2E, npol)
     fmin(x, p) = model(x, p) .- 0.3
     t0 = root_error(fmin, t[nt0], par)
     if pl
@@ -170,8 +179,9 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Un
         display(gcf())
 
         figure()
-        plot(collect(1:xmax), value.(Y_aux[:, 3]) * t[nt0]^2 / L^3, "x") 
-        fill_between(collect(plat[1]:plat[2]), v[3]+e[3], v[3]-e[3], alpha=0.75, color="green")
+        t_pl = dt0 + 1
+        plot(collect(1:xmax), value.(Y_aux[:, t_pl]) * t[nt0]^2 / L^3, "x") 
+        fill_between(collect(plat[1]:plat[2]), v[t_pl]+e[t_pl], v[t_pl]-e[t_pl], alpha=0.75, color="green")
         ylabel(L"$t^2E(x0, t)$")
         xlabel(L"$x_0/a$")
         title(string(L"$t/a^2 = $", t[nt0]))
@@ -180,7 +190,7 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Un
     return t0
 end
 
-function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Vector{Matrix{Float64}}, Nothing}=nothing)
+function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Vector{Matrix{Float64}}, Nothing}=nothing, npol::Int64=2)
     nr = length(Y)
     Ysl = getfield.(Y, :Ysl)
     t = getfield.(Y, :t)
@@ -200,19 +210,21 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false
     nt0 = t0_guess(t, tmp, plat, L)
     xmax = size(tmp, 2)
 
-    Y_aux = Matrix{uwreal}(undef, xmax, 5)
+    dt0 = iseven(npol) ? Int64(npol / 2) : Int64((npol+1) / 2)
+    Y_aux = Matrix{uwreal}(undef, xmax, 2*dt0+1)
     for i = 1:xmax
         k = 1
-        for j = nt0-2:nt0+2
+        for j = nt0-dt0:nt0+dt0
             Y_aux[i, k] = uwreal(tmp[:, i, j], id[1], replica)
             k = k + 1
         end
     end
-    x = t[nt0-2:nt0+2]
-    t2E = [plat_av(Y_aux[:, j], plat) for j=1:5] .* x.^2 / L^3
+    x = t[nt0-dt0:nt0+dt0]
+    t2E = [plat_av(Y_aux[:, j], plat) for j=1:2*dt0+1] .* x.^2 / L^3
+    
+    model(x, p) = get_model(x, p, npol)
 
-    @. model(x, p) = p[1]*x^3 + p[2]*x^2 + p[3]*x + p[4]
-    par = fit_routine(model, x, t2E, 4)
+    par = fit_routine(model, x, t2E, npol)
     fmin(x, p) = model(x, p) .- 0.3
     t0 = root_error(fmin, t[nt0], par)
     if pl
@@ -230,8 +242,9 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false
         display(gcf())
 
         figure()
-        plot(collect(1:xmax), value.(Y_aux[:, 3]) * t[nt0]^2 / L^3, "x")
-        fill_between(collect(plat[1]:plat[2]), v[3]+e[3], v[3]-e[3], alpha=0.75, color="green")
+        t_pl = dt0 + 1
+        plot(collect(1:xmax), value.(Y_aux[:, t_pl]) * t[nt0]^2 / L^3, "x")
+        fill_between(collect(plat[1]:plat[2]), v[t_pl]+e[t_pl], v[t_pl]-e[t_pl], alpha=0.75, color="green")
         ylabel(L"$t^2E(x0, t)$")
         xlabel(L"$x_0/a$")
         title(string(L"$t/a^2 = $", t[nt0]))