juobs_obs added

juobs_obs added,   meff and plat_av moved to juobs_obs, Corr now includes y0,  syntax simplified(Symbol), Plots->PyPlot, LatexString support

Manifest.toml

@@ -72,6 +72,12 @@ git-tree-sha1 = "7c4f882c41faa72118841185afc58a2eb00ef612"
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
 version = "0.3.3+0"
 
+[[Conda]]
+deps = ["JSON", "VersionParsing"]
+git-tree-sha1 = "7a58bb32ce5d85f8bf7559aa7c2842f9aecf52fc"
+uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d"
+version = "1.4.1"
+
 [[Contour]]
 deps = ["StaticArrays"]
 git-tree-sha1 = "81685fee51fc5168898e3cbd8b0f01506cd9148e"
@@ -420,6 +426,18 @@ version = "1.6.3"
 deps = ["Unicode"]
 uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
+[[PyCall]]
+deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Serialization", "VersionParsing"]
+git-tree-sha1 = "3a3fdb9000d35958c9ba2323ca7c4958901f115d"
+uuid = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
+version = "1.91.4"
+
+[[PyPlot]]
+deps = ["Colors", "LaTeXStrings", "PyCall", "Sockets", "Test", "VersionParsing"]
+git-tree-sha1 = "67dde2482fe1a72ef62ed93f8c239f947638e5a2"
+uuid = "d330b81b-6aea-500a-939a-2ce795aea3ee"
+version = "2.9.0"
+
 [[QuadGK]]
 deps = ["DataStructures", "LinearAlgebra"]
 git-tree-sha1 = "12fbe86da16df6679be7521dfb39fbc861e1dc7b"
@@ -550,6 +568,11 @@ git-tree-sha1 = "af0c29913f108f649999e74098814c7ef0f644de"
 uuid = "b8865327-cd53-5732-bb35-84acbb429228"
 version = "1.2.0"
 
+[[VersionParsing]]
+git-tree-sha1 = "80229be1f670524750d905f8fc8148e5a8c4537f"
+uuid = "81def892-9a0e-5fdd-b105-ffc91e053289"
+version = "1.2.0"
+
 [[Zlib_jll]]
 deps = ["Libdl", "Pkg"]
 git-tree-sha1 = "fdd89e5ab270ea0f2a0174bd9093e557d06d4bfa"

Project.toml

@@ -6,5 +6,6 @@ version = "0.1.0"
 [deps]
 ADerrors = "5e92007d-7bf1-471c-8ceb-4591b8b567a9"
 BDIO = "375f315e-f2c4-11e9-2ef9-134f02f79e27"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
+PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

src/juobs.jl

 module juobs
-using ADerrors, Plots, Statistics
+using ADerrors, PyPlot, Statistics, LaTeXStrings
 
 include("juobs_types.jl")
 include("juobs_reader.jl")
 include("juobs_tools.jl")
+include("juobs_obs.jl")
 
 export read_mesons, read_ms1
-export meff, plat_av, apply_rw, corr_obs
+export apply_rw, corr_obs
+export plat_av, meff, dec_const_pcvc 
 
 end # module

src/juobs_obs.jl

+function plat_av(obs::Vector{uwreal}, plat::Vector{Int64})
+    uwerr.(obs)
+    w = 1 ./ err.(obs)[plat[1]:plat[2]].^2
+    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 )     
+    dim = length(obs)                                                                                    
+    aux = log.(obs[2:dim-2] ./ obs[3:dim-1])     
+    mass = plat_av(aux, plat)                                                                            
+    uwerr(mass)                       
+    if pl == true
+        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$")
+        display(gcf())
+    end                                                  
+    if data == false                   
+        return mass                                        
+    else               
+        return (mass, aux)     
+    end                    
+end   
+meff(corr::Corr, plat::Vector{Int64}; pl=true, data=false) = meff(corr.obs, plat, pl=pl, data=data)  
+
+                                                                                                         
+## 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)
+    """
+    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] )
+    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                                
+        uwerr.(R)
+        figure()                                            
+        errorbar(1:length(R), value.(R), err.(R), fmt="x")
+        ylabel(L"$R$")
+        xlabel(L"$x_0$")
+        display(gcf())          
+    end 
+    if data == false
+        return f
+    else
+        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

src/juobs_tools.jl

-function plat_av(obs::Vector{uwreal}, plat::Vector{Int64})
-    uwerr.(obs)
-    w = 1 ./ err.(obs)[plat[1]:plat[2]].^2
-    av = sum(w .* obs[plat[1]:plat[2]]) / sum(w)
-    return av 
-end
-
-function meff(obs::Vector{uwreal}, plat::Vector{Int64})
-    dim = length(obs)
-    aux = log.(obs[1:dim-1] ./ obs[2:dim])
-    mass = plat_av(aux, plat)
-    uwerr(mass)
-    return mass
-end
-
-function apply_rw(cdata::CData, W::Array{Float64})
-    W1 = W[1, :]
-    W2 = W[2, :]
+#TODO: apply_rw with gaps
+function apply_rw(cdata::CData, W::Array{Float64, 2})
+    nc = maximum(cdata.vcfg)
+    W1 = W[1, 1:nc]
+    W2 = W[2, 1:nc]
     cdata_r = copy(cdata)
 
     cdata_r.re_data = cdata.re_data .* W1 .* W2 / mean(W1 .* W2)
@@ -23,15 +10,16 @@ function apply_rw(cdata::CData, W::Array{Float64})
     return cdata_r
 end
 
-function apply_rw(data::Array{Float64}, W::Array{Float64})
-    W1 = W[1, :]
-    W2 = W[2, :]
+function apply_rw(data::Array{Float64}, W::Array{Float64, 2})
+    nc = size(data)[1]
+    W1 = W[1, 1:nc]
+    W2 = W[2, 1:nc]
 
     data_r = data .* W1 .* W2 / mean(W1 .* W2)
     return data_r
 end
 
-function corr_obs(cdata::CData; real::Bool=true, rw::Union{Array{Float64}, Nothing}=nothing)
+function corr_obs(cdata::CData; real::Bool=true, rw::Union{Array{Float64, 2}, Nothing}=nothing)
 
     real ? data = cdata.re_data : data = cdata.im_data
     isnothing(rw) ? data_r = data : data_r = apply_rw(data, rw)
@@ -46,8 +34,8 @@ end
 #TODO: vcfg with gaps
 function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array{Float64, 2}, 1}, Nothing}=nothing)
     nr = length(cdata)
-    id = getfield.(cdata, Symbol("id"))
-    vcfg = getfield.(cdata, Symbol("vcfg"))
+    id = getfield.(cdata, :id)
+    vcfg = getfield.(cdata, :vcfg)
     replica = Int64.(maximum.(vcfg))
     
     if !all(id .== id[1])
@@ -55,7 +43,7 @@ function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array
         return nothing
     end
 
-    real ? data = getfield.(cdata, Symbol("re_data")) : data = getfield.(cdata, Symbol("im_data"))
+    real ? data = getfield.(cdata, :re_data) : data = getfield.(cdata, :im_data)
     isnothing(rw) ? data_r = data : data_r = apply_rw.(data, rw)
     
 

src/juobs_types.jl

@@ -114,24 +114,27 @@ mutable struct Corr
     obs::Vector{uwreal}
     mu::Vector{Float64}
     gamma::Vector{String}
+    y0::Int64
         function Corr(a::Vector{uwreal}, b::CData)
-            h = getfield(b, Symbol("header"))
+            h = getfield(b, :header)
             mu = [h.mu1, h.mu2]
             gamma = [gamma_name[h.type1+1], gamma_name[h.type2+1]]
-            return new(a, mu, gamma)
+            y0 = Int64(h.x0)
+            return new(a, mu, gamma, y0)
         end 
         function Corr(a::Vector{uwreal}, b::Vector{CData})
-            sym = Symbol.(["mu1", "mu2", "type1", "type2"])
-            h = getfield.(b, Symbol("header"))
+            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 missmatch")
+                    println("Corr: Parameter mismatch")
                     return nothing
                 end
             end
             mu = [h[1].mu1, h[1].mu2]
             gamma = [gamma_name[h[1].type1+1], gamma_name[h[1].type2+1]]
-            return new(a, mu, gamma)
+            y0 = Int64(h[1].x0)
+            return new(a, mu, gamma, y0)
         end 
 end
 function Base.show(io::IO, a::GHeader)
@@ -147,7 +150,7 @@ function Base.show(io::IO, a::CHeader)
     fnames = fieldnames(CHeader)
     for k in fnames
         f = getfield(a, k)
-        if k!=Symbol("type1") && k!=Symbol("type2")
+        if k != :type1 && k!= :type2
             print(io, "$k = $f\t")
         else
             print(io, "$k = ", gamma_name[f + 1], "\t")