Tools

corr_obs(cdata::CData; real::Bool=true, rw::Union{Array{Float64, 2}, Nothing}=nothing, L::Int64=1)

corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array{Float64, 2}, 1}, Nothing}=nothing, L::Int64=1)

Creates a Corr struct with the given CData struct cdata (read_mesons) for a single replica. An array of CData can be passed as argument for multiple replicas.

The flag real select the real or imaginary part of the correlator. If rw is specified, the method applies reweighting. rw is passed as a matrix of Float64 (read_ms1) The correlator can be normalized with the volume factor if L is fixed.

#Single replica
data = read_mesons(path, "G5", "G5")
rw = read_ms1(path_rw)
corr_pp = corr_obs.(data)
corr_pp_r = corr_obs.(data, rw=rw)

#Two replicas
data = read_mesons([path_r1, path_r2], "G5", "G5")
rw1 = read_ms1(path_rw1)
rw2 = read_ms1(path_rw2)

corr_pp = corr_obs.(data)
corr_pp_r = corr_obs.(data, rw=[rw1, rw2])

md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ws::ADerrors.wspace=ADerrors.wsg)

Computes the derivative of an observable A with respect to the sea quark masses.

$\frac{d <A>}{dm(sea)} = \sum_i \frac{\partial <A>}{\partial <O_i>} \frac{d <O_i>}{d m(sea)}$

$\frac{d <O_i>}{dm(sea)} = <O_i> <\frac{\partial S}{\partial m}> - <O_i \frac{\partial S}{\partial m}> = - <(O_i - <O_i>) (\frac{\partial S}{\partial m} - <\frac{\partial S}{\partial m}>)>$

where $O_i$ are primary observables

md is a vector that contains the derivative of the action $S$ with respect to the sea quark masses for each replica. md[irep][irw, icfg]

md_sea returns a tuple of uwreal observables $(dA/dm_l, dA/dm_s)|_{sea}$, where $m_l$ and $m_s$ are the light and strange quark masses.

#Single replica
data = read_mesons(path, "G5", "G5")
md = read_md(path_md)
rw = read_ms1(path_rw)

corr_pp = corr_obs.(data, rw=rw)
m = meff(corr_pp[1], plat)
m_mdl, m_mds = md_sea(m, [md], ADerrors.wsg)
m_shifted = m + 2 * dml * m_mdl + dms * m_mds

#Two replicas
data = read_mesons([path_r1, path_r2], "G5", "G5")
md1 = read_md(path_md1)
md2 = read_md(path_md2)

corr_pp = corr_obs.(data)
m = meff(corr_pp[1], plat)
m_mdl, m_mds = md_sea(m, [md1, md2], ADerrors.wsg)
m_shifted = m + 2 * dml * m_mdl + dms * m_mds

md_val(a::uwreal, obs::Corr, derm::Vector{Corr})

Computes the derivative of an observable A with respect to the valence quark masses.

$\frac{d <A>}{dm(val)} = \sum_i \frac{\partial <A>}{\partial <O_i>} \frac{d <O_i>}{d m(val)}$

$\frac{d <O_i>}{dm(val)} = <\frac{\partial O_i}{\partial m(val)}>$

where $O_i$ are primary observables

md is a vector that contains the derivative of the action $S$ with respect to the sea quark masses for each replica. md[irep][irw, icfg]

md_val returns a tuple of uwreal observables $(dA/dm_1, dA/dm_2)|_{val}$, where $m_1$ and $m_2$ are the correlator masses.

data = read_mesons(path, "G5", "G5", legacy=true)
data_d1 = read_mesons(path, "G5_d1", "G5_d1", legacy=true)
data_d2 = read_mesons(path, "G5_d2", "G5_d2", legacy=true)

rw = read_ms1(path_rw)

corr_pp = corr_obs.(data, rw=rw)
corr_pp_d1 = corr_obs.(data_d1, rw=rw)
corr_pp_d2 = corr_obs.(data_d2, rw=rw)
derm = [[corr_pp_d1[k], corr_pp_d2[k]] for k = 1:length(pp_d1)]

m = meff(corr_pp[1], plat)
m_md1, m_md2 = md_val(m, corr_pp[1], derm[1])
m_shifted = m + 2 * dm1 * m_md1 + dm2 * m_md2

lin_fit(x::Vector{<:Real}, y::Vector{uwreal})

Computes a linear fit of uwreal data points y. This method return uwreal fit parameters and chisqexpected.

fitp, csqexp = lin_fit(phi2, m2)
m2_phys = fitp[1] + fitp[2] * phi2_phys

fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}, param::Int64=3; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)

fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}, param::Int64=3; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, covar::Bool=false)

Given a model function with a number param of parameters and an array of uwreal, this function fit ydata with the given model and print fit information The method return an array upar with the best fit parameters with their errors. The flag wpm is an optional array of Float64 of lenght 4. The first three paramenters specify the criteria to determine the summation windows:

• vp[1]: The autocorrelation function is summed up to $t = round(vp[1])$.

• vp[2]: The sumation window is determined using U. Wolff poposal with $S_\tau = wpm[2]$

• vp[3]: The autocorrelation function $\Gamma(t)$ is summed up a point where its error $\delta\Gamma(t)$ is a factor vp[3] times larger than the signal.

An additional fourth parameter vp[4], tells ADerrors to add a tail to the error with $\tau_{exp} = wpm[4]$. Negative values of wpm[1:4] are ignored and only one component of wpm[1:3] needs to be positive. If the flag covaris set to true, fit_routine takes into account covariances between x and y for each data point.

@. model(x,p) = p[1] + p[2] * exp(-(p[3]-p[1])*x)
@. model2(x,p) = p[1] + p[2] * x[:, 1] + (p[3] + p[4] * x[:, 1]) * x[:, 2]
fit_routine(model, xdata, ydata, param=3)
fit_routine(model, xdata, ydata, param=3, covar=true)