documentation added

(juobs_obs, juobs_reader and juobs_tools)

src/juobs.jl

 module juobs
-using ADerrors, PyPlot, Statistics, LaTeXStrings
+using ADerrors, PyPlot, Statistics, LaTeXStrings, LinearAlgebra
 
 include("juobs_types.jl")
 include("juobs_linalg.jl")
@@ -8,7 +8,7 @@ include("juobs_tools.jl")
 include("juobs_obs.jl")
 
 export read_mesons, read_ms1
-export apply_rw, corr_obs, plat_av, lin_fit, x_lin_fit, y_lin_fit
+export corr_obs, plat_av, lin_fit, x_lin_fit, y_lin_fit
 export meff, dec_const_pcvc 
 
 end # module

src/juobs_obs.jl

-function meff(obs::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bool=false )     
-    dim = length(obs)                                                                                    
-    aux = log.(obs[2:dim-2] ./ obs[3:dim-1])     
+@doc raw"""
+meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bool=false )     
+
+meff(corr::Corr, plat::Vector{Int64}; pl::Bool=true, data::Bool=false)
+
+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)     
+    dim = length(corr)                                                                                    
+    aux = log.(corr[2:dim-2] ./ corr[3:dim-1])     
     mass = plat_av(aux, plat)                                                                            
     uwerr(mass)                       
     if pl == true
@@ -27,12 +43,29 @@ meff(corr::Corr, plat::Vector{Int64}; pl::Bool=true, data::Bool=false) =
     meff(corr.obs, plat, pl=pl, data=data)
 
 ## Decay constants                                                                                       
-#TODO: test 
-function dec_const_pcvc(obs::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false)
+@doc raw"""
+dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false)meff(corr::Corr, plat::Vector{Int64}; pl::Bool=true, data::Bool=false)
+
+dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false)
+
+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.
+
+```@example
+data = read_mesons(path, "G5", "G5")
+corr_pp = corr_obs.(data)
+m = meff(corr_pp[1], [50, 60], pl=false)
+f = dec_const_pcvc(corr_pp[1], [50, 60], m, pl=false)
+```
+"""
+function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false)
     """
     compute the decay constant when the source is far from the boundaries
     """
-    corr_pp = obs[2:end-1]
+    corr_pp = corr[2:end-1]
     dim = length(corr_pp)
                                                                  
     aux = exp.((collect(1:dim) .- y0 ) .* [m for k in 1:dim])

src/juobs_reader.jl

@@ -69,8 +69,23 @@ function read_CHeader(path::String)
     return a
 end
 
-
-function read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int, Nothing}=nothing)
+@doc raw"""
+read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int64, Nothing}=nothing)
+
+This faction read a mesons dat file at a given path and returns a vector of CData structures for different masses and Dirac structures.
+Dirac structures g1 and/or g2 can be passed as string arguments in order to filter correaltors.
+ADerrors id can be specified as argument. If is not specified, the id is fixed according to the ensemble name (example: "H400"-> id = 400)
+
+Examples:
+```@example
+read_mesons(path)
+read_mesons(path, "G5")
+read_mesons(path, nothing, "G5")
+read_mesons(path, "G5", "G5")
+read_mesons(path, "G5", "G5", id=1)
+```
+"""
+function read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int64, Nothing}=nothing)
     
     isnothing(g1) ? t1=nothing : t1 = findfirst(x-> x==g1, gamma_name) - 1
     isnothing(g2) ? t2=nothing : t2 = findfirst(x-> x==g2, gamma_name) - 1
@@ -141,8 +156,21 @@ function read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union
     return res
 end
 
-
-function read_ms1(path::String, v::String="1.2")
+@doc raw"""
+read_ms1(path::String; v::String="1.2")
+
+Reads openQCD ms1 dat files at a given path. This method returns a matrix W[irw, icfg] that contains the reweighting factors, where
+irw is the rwf index and icfg the configuration number.
+The function is compatible with the output files of openQCD v=1.2, 1.4 and 1.6. Version can be specified as argument.
+
+Examples:
+```@example    
+read_ms1(path)
+read_ms1(path, v="1.4")
+read_ms1(path, v="1.6")
+```
+"""
+function read_ms1(path::String; v::String="1.2")
     data = open(path, "r")
     nrw = read(data, Int32)
 

src/juobs_tools.jl

@@ -18,7 +18,37 @@ function apply_rw(data::Array{Float64}, W::Array{Float64, 2})
     data_r = data .* W1 .* W2 / mean(W1 .* W2)
     return data_r
 end
+@doc raw"""
+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.
+
+@example```
+#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_r1 = read_mesons(path_r1, "G5", "G5")
+data_r2 = read_mesons(path_r2, "G5", "G5")
+rw1 = read_ms1(path_rw1)
+rw2 = read_ms1(path_rw2)
+
+cdata = [[data_r1[k], data_r2[k]] for k=1:length(data_r1)]
+
+corr_pp = corr_obs.(cdata)
+corr_pp_r = corr_obs.(cdata, rw=[rw1, rw2])
+```
+"""
 function corr_obs(cdata::CData; real::Bool=true, rw::Union{Array{Float64, 2}, Nothing}=nothing, L::Int64=1)
 
     real ? data = cdata.re_data ./ L^3 : data = cdata.im_data ./ L^3
@@ -77,7 +107,16 @@ function lin_fit(x::Vector{Float64}, v::Vector{Float64}, e::Vector{Float64})
     #C = [[Sxx/delta, -Sx/delta], [-Sx/delta,  S/delta]]
     return par
 end
+@doc raw"""
+lin_fit(x::Vector{Float64}, y::Vector{uwreal})
+
+Computes a linear fit of uwreal data points y. This method return uwreal fit parameters and chisqexpected.
 
+```@example
+fitp, csqexp = lin_fit(phi2, m2)
+m2_phys = fitp[1] + fitp[2] * phi2_phys
+```
+"""
 function lin_fit(x::Vector{Float64}, y::Vector{uwreal})
     uwerr.(y)
     par = lin_fit(x, value.(y), err.(y))
@@ -88,7 +127,17 @@ function lin_fit(x::Vector{Float64}, y::Vector{uwreal})
     return (fitp, csqexp)
 end
 
+@doc raw"""
+x_lin_fit(par::Vector{uwreal}, y::Union{uwreal, Float64})
+
+Computes the results of a linear interpolation/extrapolation in the x axis
+"""
 x_lin_fit(par::Vector{uwreal}, y::Union{uwreal, Float64}) = (y - par[1]) / par[2]
+@doc raw"""
+y_lin_fit(par::Vector{uwreal}, y::Union{uwreal, Float64})
+
+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