Merge branch 'javier'

src/juobs.jl

@@ -10,7 +10,7 @@ include("juobs_obs.jl")
 
 export read_mesons, read_ms1, read_ms, read_md, truncate_data!
 export get_matrix, energies, uwdot, uweigvals, uweigvecs, uweigen, invert, getall_eigvals, getall_eigvecs
-export corr_obs, md_sea, md_val, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine
+export corr_obs, corr_sym, md_sea, md_val, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine
 export meff, mpcac, dec_const, dec_const_pcvc, comp_t0
 
 end # module

src/juobs_obs.jl

@@ -20,7 +20,7 @@ function meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bo
     dim = length(corr)                                                                                    
     aux = 0.5 .* log.((corr[2:dim-2] ./ corr[3:dim-1]).^2)     
     mass = plat_av(aux, plat, wpm)                                                                            
-    if pl == true
+    if pl
         isnothing(wpm) ? uwerr(mass) : uwerr(mass, wpm)                       
         x = 1:length(aux)
         y = value.(aux)
@@ -42,12 +42,13 @@ function meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bo
         end
         display(gcf())
     end                                                  
-    if data == false                   
+    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)
     
@@ -98,7 +99,7 @@ function mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca:
 
     aux = der_a0p ./ (2 .* corr_pp[2:end-1])
     mass = plat_av(aux, plat, wpm)                                                                            
-    if pl == true
+    if pl
         isnothing(wpm) ? uwerr(mass) : uwerr(mass, wpm)                       
         x = 1:length(aux)
         y = value.(aux)
@@ -120,12 +121,13 @@ function mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca:
         end
         display(gcf())
     end                                                  
-    if data == false                   
+    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)
     
@@ -145,6 +147,11 @@ end
 
     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.
@@ -154,6 +161,10 @@ 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)
@@ -161,6 +172,9 @@ 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
@@ -170,6 +184,8 @@ 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,
@@ -191,7 +207,7 @@ function dec_const(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64},
     R_av = plat_av(R, plat, wpm)
 
     f = sqrt(2) * sqrt(R_av^2) / sqrt(m)
-    if pl == true
+    if pl
         isnothing(wpm) ? uwerr(R_av) : uwerr(R_av, wpm)
         isnothing(wpm) ? uwerr(f) : uwerr(f, wpm)
         x = 1:length(R)
@@ -219,6 +235,7 @@ function dec_const(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64},
         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) 
     
@@ -228,6 +245,69 @@ function dec_const(a0p::Corr, pp::Corr, plat::Vector{Int64}, m::uwreal; ca::Floa
         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
 @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)
 
@@ -259,7 +339,7 @@ function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu
     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 == true
+    if pl 
         if isnothing(wpm)
             uwerr(f)
             uwerr(R_av)
@@ -280,12 +360,13 @@ function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu
         title(string(L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2]))
         display(gcf())          
     end 
-    if data == false
+    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) 
     
@@ -497,6 +578,7 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false
     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
     t0_aux = minimum(abs.(t2E_ax .- 0.3))

src/juobs_tools.jl

@@ -125,7 +125,35 @@ function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array
 
     return Corr(obs, cdata)
 end    
+@doc raw"""
+    corr_sym(corrL::Corr, corrR::Corr, parity::Int64=1)
+
+Computes the symmetrized correlator using the left correlador `corrL` and the right correlator `corrR`. The source position
+of `corrR` must be `T - 1 - y0`, where `y0` is the source position of `corrL`. 
 
+```@example
+pp_sym = corr_sym(ppL, ppR, +1)
+a0p_sym = corr_sym(a0pL, a0pR, -1)
+```
+"""
+function corr_sym(corrL::Corr, corrR::Corr, parity::Int64=1)
+    T = length(corrL.obs)
+    sym = [:kappa, :mu, :gamma]
+    if corrL.y0 != T - 1 - corrR.y0
+        error("Corr: Parameter mismatch")
+    end
+    for s in sym
+        if getfield(corrL, s) != getfield(corrR, s)
+            error("Corr: Parameter mismatch")
+        end
+    end
+    if abs(parity) != 1
+        error("incorrect value of parity (+- 1)")
+    end
+
+    res = (corrL.obs[1:end] + parity * corrR.obs[end:-1:1]) / 2
+    return Corr(res, corrL.kappa, corrL.mu, corrL.gamma, corrL.y0)
+end
 #TODO: VECTORIZE, uwreal?
 @doc raw"""
     md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ws::ADerrors.wspace=ADerrors.wsg)

src/juobs_types.jl

@@ -167,6 +167,7 @@ mutable struct Corr
         y0 = Int64(h[1].x0)
         return new(a, kappa, mu, gamma, y0)
     end 
+    Corr(o::Vector{uwreal}, k::Vector{Float64}, m::Vector{Float64}, g::Vector{String}, y::Int64) = new(o, k, m, g, y)
 end
 
 mutable struct YData