@doc raw""" meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) meff(corr::Corr, plat::Vector{Int64}; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) Computes effective mass for a given correlator corr at a given plateau `plat`. Correlator can be passed as an `Corr` struct or `Vector{uwreal}`. The flags `pl` and `data` allow to show the plots and return data as an extra result. ```@example data = read_mesons(path, "G5", "G5") corr_pp = corr_obs.(data) 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, kappa::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) dim = length(corr) aux = 0.5 .* log.((corr[2:dim-2] ./ corr[3:dim-1]).^2) mass = plat_av(aux, plat, wpm) if pl isnothing(wpm) ? uwerr(mass) : uwerr(mass, wpm) x = 1:length(aux) y = value.(aux) dy = err.(aux) v = value(mass) e = err(mass) figure() fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75) errorbar(x, y, dy, fmt="x", color="black") ylabel(L"$m_\mathrm{eff}$") xlabel(L"$x_0$") ylim(v-100*e, v+100*e) if !isnothing(kappa) title(string(L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2])) end if !isnothing(mu) title(string(L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2])) end display(gcf()) end if !data return mass else 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) if corr.mu == [0.0, 0.0] return meff(corr.obs, plat, pl=pl, data=data, kappa=corr.kappa, mu=nothing, wpm=wpm) else return meff(corr.obs, plat, pl=pl, data=data, mu=corr.mu, kappa=nothing, wpm=wpm) end end @doc raw""" mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) mpcac(a0p::Corr, pp::Corr, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) Computes the bare PCAC mass for a given correlator `a0p` and `pp` at a given plateau `plat`. Correlator can be passed as an `Corr` struct or `Vector{uwreal}`. The flags `pl` and `data` allow to show the plots and return data as an extra result. The `ca` variable allows to compute `mpcac` using the improved axial current. ```@example data_pp = read_mesons(path, "G5", "G5") data_a0p = read_mesons(path, "G5", "G0G5") corr_pp = corr_obs.(data_pp) corr_a0p = corr_obs.(data_a0p) m12 = mpcac(corr_a0p, corr_pp, [50, 60], pl=false) p0 = 9.2056 p1 = -13.9847 g2 = 1.73410 ca = -0.006033 * g2 *( 1 + exp(p0 + p1/g2)) m12 = mpcac(corr_a0p, corr_pp, [50, 60], pl=false, ca=ca) ``` """ function mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false, kappa::Union{Vector{Float64}, Nothing}=nothing, mu::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) corr_a0p = -a0p[1:end] corr_pp = pp[1:end] der_a0p = (corr_a0p[3:end] .- corr_a0p[1:end-2]) / 2 if ca != 0.0 der2_pp = (corr_pp[1:end-4] - 2*corr_pp[3:end-2] + corr_pp[5:end])/4 der_a0p = der_a0p[2:end-1] + ca * der2_pp end #der_a0p = (corr_a0p[3:end] .- corr_a0p[1:end-2]) / 2 #if ca != 0.0 # der2_pp = corr_pp[1:end-2] + corr_pp[3:end] - 2 * corr_pp[2:end-1] # der_a0p = der_a0p + ca * der2_pp #end aux = der_a0p ./ ( corr_pp[3:end-2]) mass = plat_av(aux, plat, wpm) if pl isnothing(wpm) ? uwerr(mass) : uwerr(mass, wpm) x = 1:length(aux) y = value.(aux) dy = err.(aux) v = value(mass) e = err(mass) figure() fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75) errorbar(x, y, dy, fmt="x", color="black") ylabel(L"$m_\mathrm{PCAC}$") xlabel(L"$x_0$") ylim(v-100*e, v+100*e) if !isnothing(kappa) title(string(L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2])) end if !isnothing(mu) title(string(L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2])) end display(gcf()) end if !data return mass else return (mass, aux) end end function mpcac(a0p::Corr, pp::Corr, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) if a0p.kappa == pp.kappa || a0p.mu == pp.mu if a0p.mu == [0.0, 0.0] return mpcac(a0p.obs, pp.obs, plat, ca=ca, pl=pl, data=data, kappa=a0p.kappa, mu=nothing, wpm=wpm) else return mpcac(a0p.obs, pp.obs, plat, ca=ca, pl=pl, data=data, mu=a0p.mu, kappa=nothing, wpm=wpm) end else error("mu or kappa values does not match") end end ## Decay constants @doc raw""" dec_const(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) dec_const(a0p::Corr, pp::Corr, plat::Vector{Int64}, m::uwreal; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) dec_const(a0pL::Vector{uwreal}, a0pR::Vector{uwreal}, ppL::Vector{uwreal}, ppR::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) dec_const(a0pL::Corr, a0pR::Corr, ppL::Corr, ppR::Corr, plat::Vector{Int64}, m::uwreal; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) Computes the bare decay constant using ``A_0P`` and ``PP`` correlators . The decay constant is computed in the plateau `plat`. Correlator can be passed as an `Corr` struct or `Vector{uwreal}`. If it is passed as a uwreal vector, effective mass `m` and source position `y0` must be specified. The flags `pl` and `data` allow to show the plots and return data as an extra result. The `ca` variable allows to compute `dec_const` using the improved axial current. **The method assumes that the source is close to the boundary.** It takes the following ratio to cancel boundary effects. ``R = \frac{f_A(x_0, y_0)}{\sqrt{f_P(T-y_0, y_0)}} * e^{m (x_0 - T/2)}`` **If left and right correlators are included in the input. The result is computed with the following ratio** `` R = \sqrt{f_A(x_0, y_0) * f_A(x_0, T - 1 - y_0) / f_P(T - 1 - y_0, y_0)}`` ```@example data_pp = read_mesons(path, "G5", "G5", legacy=true) data_a0p = read_mesons(path, "G5", "G0G5", legacy=true) corr_pp = corr_obs.(data_pp, L=32) corr_a0p = corr_obs.(data_a0p, L=32) corr_a0pL, corr_a0pR = [corr_a0p[1], corr_a0p[2]] corr_ppL, corr_ppR = [corr_pp[1], corr_pp[2]] m = meff(corr_pp[1], [50, 60], pl=false) beta = 3.46 p0 = 9.2056 p1 = -13.9847 g2 = 6 / beta ca = -0.006033 * g2 *( 1 + exp(p0 + p1/g2)) f = dec_const(corr_a0p[1], corr_pp[1], [50, 60], m, pl=true, ca=ca) f_ratio = dec_const(corr_a0pL, corr_a0pR, corr_ppL, corr_ppR, [50, 60], m, pl=true, ca=ca) ``` """ function dec_const(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=true, data::Bool=false, kappa::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) corr_a0p = -a0p corr_pp = pp T = length(corr_a0p) if ca != 0.0 der_pp = (corr_pp[3:end] - corr_pp[1:end-2]) / 2 corr_a0p = corr_a0p[2:end-1] + ca * der_pp aux = exp.((collect(1:T-2) .- (T-1)/2) .* [m for k = 1:T-2]) else aux = exp.((collect(0:T-1) .- T/2) .* [m for k = 1:T]) end R = corr_a0p .* aux ./ [((corr_pp[T-y0])^2)^(1/4) for k = 1:length(corr_a0p)] R_av = plat_av(R, plat, wpm) f = sqrt(2) * sqrt(R_av^2) / sqrt(m) if pl isnothing(wpm) ? uwerr(R_av) : uwerr(R_av, wpm) isnothing(wpm) ? uwerr(f) : uwerr(f, wpm) x = 1:length(R) y = value.(R) dy = err.(R) v = value(R_av) e = err(R_av) figure() lbl = string(L"$af = $", sprint(show, f)) 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() ylim(v-10*e, v+10*e) ylabel(L"$R_\mathrm{av}$") xlabel(L"$x_0$") if !isnothing(kappa) title(string(L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2])) end display(gcf()) close() end if !data return f else return (f, R) end end function dec_const(a0p::Corr, pp::Corr, plat::Vector{Int64}, m::uwreal; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) if (a0p.y0 == pp.y0) && (a0p.kappa == pp.kappa) return dec_const(a0p.obs, pp.obs, plat, m, a0p.y0, ca=ca, kappa=a0p.kappa, pl=pl, data=data, wpm=wpm) else error("y0 or kappa values does not match") end end function dec_const(a0pL::Vector{uwreal}, a0pR::Vector{uwreal}, ppL::Vector{uwreal}, ppR::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=true, data::Bool=false, kappa::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) corr_pp = (ppL + ppR[end:-1:1]) / 2 T = length(corr_pp) if ca != 0.0 der_ppL = (ppL[3:end] - ppL[1:end-2]) / 2 der_ppR = (ppR[3:end] - ppR[1:end-2]) / 2 corr_a0pL = -a0pL[2:end-1] + ca * der_ppL corr_a0pR = -a0pR[2:end-1] + ca * der_ppR else corr_a0pL = -a0pL[2:end-1] corr_a0pR = -a0pR[2:end-1] end f1 = [corr_pp[T - y0] for k = 1:length(corr_a0pL)] R = ((corr_a0pL .* corr_a0pR ./ f1).^2).^(1/4) R_av = plat_av(R, plat, wpm) f = sqrt(2) * sqrt(R_av^2) / sqrt(m) if pl isnothing(wpm) ? uwerr(R_av) : uwerr(R_av, wpm) isnothing(wpm) ? uwerr(f) : uwerr(f, wpm) x = 1:length(R) y = value.(R) dy = err.(R) v = value(R_av) e = err(R_av) figure() lbl = string(L"$af = $", sprint(show, f)) 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() ylabel(L"$R_\mathrm{av}$") xlabel(L"$x_0$") if !isnothing(kappa) title(string(L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2])) end display(gcf()) end if !data return f else return (f, R) end end function dec_const(a0pL::Corr, a0pR::Corr, ppL::Corr, ppR::Corr, plat::Vector{Int64}, m::uwreal; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) T = length(a0pL.obs) if (a0pL.y0 == ppL.y0) && (a0pR.y0 == ppR.y0) && (a0pL.kappa == ppL.kappa) && (a0pR.kappa == ppR.kappa) && (a0pL.y0 == T - 1 - a0pR.y0) return dec_const(a0pL.obs, a0pR.obs, ppL.obs, ppR.obs, plat, m, a0pL.y0, ca=ca, kappa=a0pL.kappa, pl=pl, data=data, wpm=wpm) else error("y0 or kappa values does not match") end end ## ADD DOCUMENTATION FOR VECTOR DECAY function dec_const(vv::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; pl::Bool=true, data::Bool=false, kappa::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) corr_vv = vv[2:end-1] T = length(corr_vv) aux = exp.((collect(1:T) .- y0 ) .* fill(m, T)) R = ((aux .* corr_vv).^2).^0.25 R_av = plat_av(R, plat, wpm) f = sqrt(2 / m) * R_av if pl R .*= sqrt.(2 ./ [m for i in 1:length(R)]) uwerr.(R) R_av *= sqrt(2/m) isnothing(wpm) ? uwerr(R_av) : uwerr(R_av, wpm) isnothing(wpm) ? uwerr(f) : uwerr(f, wpm) x = 1:length(R) y = value.(R) dy = err.(R) v = value(R_av) e = err(R_av) figure() lbl = string(L"$af = $", sprint(show, f)) 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) ylim(v-10*e, v+10*e) ylabel(L"$af_{B}$") #ylabel(L"$R_\mathrm{av}$") xlabel(L"$x_0$") #title(L"$f_{B^*}$") if !isnothing(kappa) title(string( L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2])) end display(gcf()) #t = "ps_decay_$(kappa[1])_$(kappa[2]).pdf" #savefig(joinpath("/Users/ale/Il mio Drive/phd/secondment/3pf test/analysis/plots",t)) end if !data return f else return f, R end end function dec_const(vv::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) return dec_const(vv.obs, plat, m, vv.y0+1, kappa=vv.kappa, pl=pl, data=data, wpm=wpm) end @doc raw""" 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) dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) dec_const_pcvc(ppL::Vector{uwreal}, ppR::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) dec_const_pcvc(corrL::Corr, corrR::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) Computes decay constant using the PCVC relation for twisted mass fermions. The decay constant is computed in the plateau `plat`. Correlator can be passed as an `Corr` struct or `Vector{uwreal}`. If it is passed as a uwreal vector, vector of twisted masses `mu` and source position `y0` must be specified. The flags `pl` and `data` allow to show the plots and return data as an extra result. **The method extract the matrix element assuming that the source is in the bulk. ** **If left and right correlators are included in the input. The result is computed with a ratio that cancels boundary effects:** `` R = \sqrt{f_P(x_0, y_0) * f_P(x_0, T - 1 - y_0) / f_P(T - 1 - y_0, y_0)}`` ```@example data = read_mesons(path, "G5", "G5") corr_pp = corr_obs.(data, L=32) m = meff(corr_pp[1], [50, 60], pl=false) f = dec_const_pcvc(corr_pp[1], [50, 60], m, pl=false) #left and right correlators f_ratio = dec_const_pcvc(ppL, ppR, [50, 60], m) ``` """ 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 """ corr_pp = corr[2:end-1] dim = length(corr_pp) aux = exp.((collect(1:dim) .- y0 ) .* [m for k in 1:dim]) R = ((aux .* corr_pp).^2).^0.25 R_av = plat_av(R, plat, wpm) f = sqrt(2) * (mu[1] + mu[2]) *R_av / m^1.5 if pl isnothing(wpm) ? uwerr(f) : uwerr(f, wpm) isnothing(wpm) ? uwerr(R_av) : uwerr(R_av, wpm) v = value(R_av) e = err(R_av) figure() lbl = string(L"$af = $", sprint(show, f)) fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75, label=L"$R$") errorbar(1:length(R), value.(R), err.(R), fmt="x", color="black", label=lbl) ylabel(L"$R$") xlabel(L"$x_0$") legend() title(string(L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2])) display(gcf()) end if !data return f else return (f, uwdot(sqrt(2) * (mu[1] + mu[2]) / m^1.5, R)) end end function dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) return dec_const_pcvc(corr.obs, plat, m, corr.mu, corr.y0, pl=pl, data=data, wpm=wpm) end function dec_const_pcvc(ppL::Vector{uwreal}, ppR::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) T = length(ppL) f1 = [ppL[T - y0] for k = 1:T] R = ((ppL .* ppR ./ f1).^2).^(1/4) R = R[2:end-1] R_av = plat_av(R, plat, wpm) f = sqrt(2) * (mu[1] + mu[2]) * R_av / m^1.5 if pl isnothing(wpm) ? uwerr(f) : uwerr(f, wpm) isnothing(wpm) ? uwerr(R_av) : uwerr(R_av, wpm) v = value(R_av) e = err(R_av) figure() lbl = string(L"$af = $", sprint(show, f)) fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75, label=L"$R$") errorbar(1:length(R), value.(R), err.(R), fmt="x", color="black", label=lbl) ylabel(L"$R$") xlabel(L"$x_0$") legend() title(string(L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2])) display(gcf()) end if !data return f else return (f, uwdot(sqrt(2) * (mu[1] + mu[2]) / m^1.5, R)) end end function dec_const_pcvc(corrL::Corr, corrR::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) T = length(corrL.obs) if (corrL.kappa == corrR.kappa) && (corrL.mu == corrR.mu) && (corrL.y0 == T - 1 - corrR.y0) return dec_const_pcvc(corrL.obs, corrR.obs, plat, m, corrL.mu, corrL.y0, pl=pl, data=data, wpm=wpm) else error("y0, kappa or mu values does not match") end end #t0 function get_model(x, p, n) s = 0.0 for k = 1:n s = s .+ p[k] .* x.^(k-1) end return s end @doc raw""" comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Matrix{Float64}, Nothing}=nothing, npol::Int64=2, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{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, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) Computes `t0` using the energy density of the action `Ysl`(Yang-Mills action). `t0` is computed in the plateau `plat`. 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. ```@example #Single replica Y = read_ms(path) rw = read_ms(path_rw) t0 = comp_t0(Y, [38, 58], L=32) t0_r = comp_t0(Y, [38, 58], L=32, rw=rw) #Two replicas Y1 = read_ms(path1) Y2 = read_ms(path2) rw1 = read_ms(path_rw1) rw2 = read_ms(path_rw2) t0 = comp_t0([Y1, Y2], [38, 58], L=32, pl=true) 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, npol::Int64=2, ws::ADerrors.wspace=ADerrors.wsg, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) Ysl = Y.obs t = Y.t id = Y.id replica = size.([Ysl], 1) #Truncation if id in keys(ADerrors.wsg.str2id) n_ws = findfirst(x-> x == ws.str2id[id], ws.map_nob) if !isnothing(n_ws) ivrep_ws = ws.fluc[n_ws].ivrep if length(ivrep_ws) != 1 error("Different number of replicas") end if replica[1] > ivrep_ws[1] println("Automatic truncation in Ysl ", ivrep_ws[1], " / ", replica[1], ". R = 1") Ysl = Ysl[1:ivrep_ws[1], :, :] elseif replica[1] < ivrep_ws[1] error("Automatic truncation failed. R = 1\nTry using truncate_data!") end end end Ysl = isnothing(rw) ? Ysl : apply_rw(Ysl, rw) xmax = size(Ysl, 2) nt0 = t0_guess(t, Ysl, plat, L) 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-dt0:nt0+dt0 Y_aux[i, k] = uwreal(Ysl[:, i, j], id) k = k + 1 end end x = t[nt0-dt0:nt0+dt0] 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) par = fit_routine(model, x, t2E, npol) fmin(x, p) = model(x, p) .- 0.3 t0 = root_error(fmin, t[nt0], par) if pl if isnothing(wpm) uwerr(t0) uwerr.(t2E) else uwerr(t0, wpm) [uwerr(t2E_aux, wpm) for t2E_aux in t2E] end v = value.(t2E) e = err.(t2E) figure() errorbar(x, v, e, fmt="x") errorbar(value(t0), 0.3, xerr=err(t0), fmt="x") ylabel(L"$t^2E$") xlabel(L"$t/a^2$") display(gcf()) figure() 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])) display(gcf()) end return t0 end 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, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) nr = length(Y) Ysl = getfield.(Y, :obs) t = getfield.(Y, :t) t = t[1] id = getfield.(Y, :id) replica = size.(Ysl, 1) if !all(id .== id[1]) error("IDs are not equal") end id = id[1] #Truncation if id in keys(ADerrors.wsg.str2id) n_ws = findfirst(x-> x == ws.str2id[id], ws.map_nob) if !isnothing(n_ws) ivrep_ws = ws.fluc[n_ws].ivrep if length(replica) != length(ivrep_ws) error("Different number of replicas") end for k = 1:length(replica) if replica[k] > ivrep_ws[k] println("Automatic truncation in Ysl ", ivrep_ws[k], " / ", replica[k], ". R = ", k) Ysl[k] = Ysl[k][1:ivrep_ws[k], :, :] elseif replica[k] < ivrep_ws[k] error("Automatic truncation failed. R = ", replica[k], "\nTry using truncate_data!") end end replica = size.(Ysl, 1) end end Ysl = isnothing(rw) ? Ysl : apply_rw(Ysl, rw) tmp = Ysl[1] [tmp = cat(tmp, Ysl[k], dims=1) for k = 2:nr] nt0 = t0_guess(t, tmp, plat, L) xmax = size(tmp, 2) 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-dt0:nt0+dt0 Y_aux[i, k] = uwreal(tmp[:, i, j], id, replica) k = k + 1 end end x = t[nt0-dt0:nt0+dt0] 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) par = fit_routine(model, x, t2E, npol) fmin(x, p) = model(x, p) .- 0.3 t0 = root_error(fmin, t[nt0], par) if pl if isnothing(wpm) uwerr(t0) uwerr.(t2E) else uwerr(t0, wpm) [uwerr(t2E_aux, wpm) for t2E_aux in t2E] end v = value.(t2E) e = err.(t2E) figure() errorbar(x, v, e, fmt="x") errorbar(value(t0), 0.3, xerr=err(t0), fmt="x") ylabel(L"$t^2E$") xlabel(L"$t/a^2$") display(gcf()) figure() 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])) display(gcf()) end return t0 end function t0_guess(t::Vector{Float64}, Ysl::Array{Float64, 3}, plat::Vector{Int64}, L::Int64) t2E_ax = t.^2 .* mean(mean(Ysl[:, plat[1]:plat[2], :], dims=2), dims=1)[1, 1, :] / L^3 #look for values t2E: t2E > 0.3 and 0.0 < t2E_ax < 0.3 if !(any(t2E_ax .> 0.3) && any((t2E_ax .< 0.3) .& (t2E_ax .> 0.0))) error("Error finding solutions. Check volume normalization!") end t0_aux = minimum(abs.(t2E_ax .- 0.3)) nt0 = findfirst(x-> abs(x - 0.3) == t0_aux, t2E_ax) #index closest value to t0 return nt0 end