additional methods for meff and dec_const_pcvc

dec_const_pcvc is completed, new methods, Corr->Obs, Obs constructors added,  L normalization

src/juobs_obs.jl

@@ -4,8 +4,7 @@ function plat_av(obs::Vector{uwreal}, plat::Vector{Int64})
     av = sum(w .* obs[plat[1]:plat[2]]) / sum(w)
     return av 
 end
-#TODO: mu1, y0, gamma -> return?
-function meff(obs::Vector{uwreal}, plat::Vector{Int64}; pl=true, data=false )     
+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])     
     mass = plat_av(aux, plat)                                                                            
@@ -30,30 +29,31 @@ function meff(obs::Vector{uwreal}, plat::Vector{Int64}; pl=true, data=false )
         return (mass, aux)     
     end                    
 end   
-meff(corr::Corr, plat::Vector{Int64}; pl=true, data=false) = meff(corr.obs, plat, pl=pl, data=data)  
+meff(corr::Obs, plat::Vector{Int64}; pl::Bool=true) = 
+    Obs(meff(corr.obs, plat, pl=pl, data=false), corr)  
 
-                                                                                                         
 ## Decay constants                                                                                       
-#TODO: mu1, y0, gamma -> return?
-#TODO: test
-function dec_const_pcvc(obs::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl=true, data=false)
+#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)
     """
     compute the decay constant when the source is far from the boundaries
     """
     corr_pp = obs[2:end-1]
     dim = length(corr_pp)
-
-    #m_array = Vector{uwreal}()
-    #[push!(m_array, m) for k in 1:dim]                                                                  
-    aux = exp.( (collect(1:dim) .- y0 ) .* [m for k in 1:dim] )
+                                                                 
+    aux = exp.((collect(1:dim) .- y0 ) .* [m for k in 1:dim])
     R = (aux .* corr_pp).^0.5
     R_av = plat_av(R, plat)
     f = sqrt(2) * (mu[1] + mu[2]) *R_av / m^1.5
     uwerr(f)
-    if pl == true                                
+    if pl == true
+        uwerr(R_av)
+        v = value(R_av)
+        e = err(R_av)                                
         uwerr.(R)
-        figure()                                            
-        errorbar(1:length(R), value.(R), err.(R), fmt="x")
+        figure()
+        fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75)                                            
+        errorbar(1:length(R), value.(R), err.(R), fmt="x", color="black")
         ylabel(L"$R$")
         xlabel(L"$x_0$")
         display(gcf())          
@@ -64,5 +64,9 @@ function dec_const_pcvc(obs::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu:
         return (f, R)                   
     end 
 end 
-dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl=true, data=false) =
-    dec_const_pcvc(corr.obs, plat, m, corr.mu, corr.y0, pl=pl, data=data)
\ No newline at end of file
+dec_const_pcvc(corr::Obs, plat::Vector{Int64}, m::uwreal; pl::Bool=true) =
+    Obs(dec_const_pcvc(corr.obs, plat, m, 
+    corr.mu, corr.y0, pl=pl, data=false), corr)
+dec_const_pcvc(corr::Obs, plat::Vector{Int64}, m::Obs; pl::Bool=true) =
+    Obs(dec_const_pcvc(corr.obs, plat, m.obs, 
+    corr.mu, corr.y0, pl=pl, data=false), corr)
\ No newline at end of file

src/juobs_reader.jl

@@ -70,7 +70,6 @@ function read_CHeader(path::String)
 end
 
 
-
 function read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int, Nothing}=nothing)
     
     isnothing(g1) ? t1=nothing : t1 = findfirst(x-> x==g1, gamma_name) - 1

src/juobs_tools.jl

@@ -19,20 +19,20 @@ function apply_rw(data::Array{Float64}, W::Array{Float64, 2})
     return data_r
 end
 
-function corr_obs(cdata::CData; real::Bool=true, rw::Union{Array{Float64, 2}, Nothing}=nothing)
+function corr_obs(cdata::CData; real::Bool=true, rw::Union{Array{Float64, 2}, Nothing}=nothing, L::Int64=1)
 
-    real ? data = cdata.re_data : data = cdata.im_data
+    real ? data = cdata.re_data ./ L^3 : data = cdata.im_data ./ L^3
     isnothing(rw) ? data_r = data : data_r = apply_rw(data, rw)
 
     nt = size(data)[2]
     obs = Vector{uwreal}(undef, nt)
     [obs[x0] = uwreal(data_r[:, x0], cdata.id) for x0 = 1:nt]
-    return Corr(obs, cdata)
+    return Obs(obs, cdata)
 end
 
 #function corr_obs for R != 1
 #TODO: vcfg with gaps
-function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array{Float64, 2}, 1}, Nothing}=nothing)
+function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array{Float64, 2}, 1}, Nothing}=nothing, L::Int64=1)
     nr = length(cdata)
     id = getfield.(cdata, :id)
     vcfg = getfield.(cdata, :vcfg)
@@ -43,7 +43,7 @@ function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array
         return nothing
     end
 
-    real ? data = getfield.(cdata, :re_data) : data = getfield.(cdata, :im_data)
+    real ? data = getfield.(cdata, :re_data) ./ L^3 : data = getfield.(cdata, :im_data) ./ L^3
     isnothing(rw) ? data_r = data : data_r = apply_rw.(data, rw)
     
 
@@ -55,5 +55,5 @@ function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array
 
     [obs[x0] = uwreal(tmp[:, x0], id[1], replica) for x0 = 1:nt]
 
-    return Corr(obs, cdata)
+    return Obs(obs, cdata)
 end    
\ No newline at end of file

src/juobs_types.jl

@@ -110,33 +110,36 @@ mutable struct CData
 end
 Base.copy(a::CData) = CData(a.header, a.vcfg, a.re_data, a.im_data, a.id)
 
-mutable struct Corr
-    obs::Vector{uwreal}
+mutable struct Obs
+    obs::Union{Vector{uwreal}, uwreal}
     mu::Vector{Float64}
     gamma::Vector{String}
     y0::Int64
-        function Corr(a::Vector{uwreal}, b::CData)
-            h = getfield(b, :header)
-            mu = [h.mu1, h.mu2]
-            gamma = [gamma_name[h.type1+1], gamma_name[h.type2+1]]
-            y0 = Int64(h.x0)
-            return new(a, mu, gamma, y0)
-        end 
-        function Corr(a::Vector{uwreal}, b::Vector{CData})
-            sym = [:mu1, :mu2, :type1, :type2, :x0]
-            h = getfield.(b, :header)
-            for s in sym
-                if !all(getfield.(h, s) .== getfield(h[1], s))
-                    println("Corr: Parameter mismatch")
-                    return nothing
-                end
+    function Obs(a::Vector{uwreal}, b::CData)
+        h = getfield(b, :header)
+        mu = [h.mu1, h.mu2]
+        gamma = [gamma_name[h.type1+1], gamma_name[h.type2+1]]
+        y0 = Int64(h.x0)
+        return new(a, mu, gamma, y0)
+    end 
+    function Obs(a::Vector{uwreal}, b::Vector{CData})
+        sym = [:mu1, :mu2, :type1, :type2, :x0]
+        h = getfield.(b, :header)
+        for s in sym
+            if !all(getfield.(h, s) .== getfield(h[1], s))
+                println("Obs: Parameter mismatch")
+                return nothing
             end
-            mu = [h[1].mu1, h[1].mu2]
-            gamma = [gamma_name[h[1].type1+1], gamma_name[h[1].type2+1]]
-            y0 = Int64(h[1].x0)
-            return new(a, mu, gamma, y0)
-        end 
+        end
+        mu = [h[1].mu1, h[1].mu2]
+        gamma = [gamma_name[h[1].type1+1], gamma_name[h[1].type2+1]]
+        y0 = Int64(h[1].x0)
+        return new(a, mu, gamma, y0)
+    end 
+    Obs(a::uwreal, b::Obs) = new(a, b.mu, b.gamma, b.y0)
+    Obs(a::Vector{uwreal}, b::Obs) = new(a, b.mu, b.gamma, b.y0)
 end
+
 function Base.show(io::IO, a::GHeader)
     print(io, "ncorr = ", a.ncorr, "\t")
     print(io, "nnoise = ", a.nnoise, "\t")
@@ -145,6 +148,7 @@ function Base.show(io::IO, a::GHeader)
     print(io, "\n")
 
 end
+
 function Base.show(io::IO, a::CHeader)
 
     fnames = fieldnames(CHeader)