dec_const_pcvc new method

dec_const uses time reflected correlators

src/juobs_obs.jl

@@ -313,18 +313,27 @@ end
 
     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)
 
+    dec_const_pcvc(ppL::Vector{uwreal}, ppR::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)
+
+    dec_const_pcvc(corrL::Corr, corrR::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+
 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.
 
-**The method assumes that the source is in the bulk.**
+**The method extract the matrix element assuming that the source is in the bulk. **
+**If left and right correlators are included in the input. The result is computed with a ratio that cancels boundary effects:**
+`` R = \sqrt{f_P(x_0, y_0) * f_P(x_0, T - 1 - y_0) / f_P(T - 1 - y_0, y_0)}``
 ```@example
 data = read_mesons(path, "G5", "G5")
 corr_pp = corr_obs.(data, L=32)
 m = meff(corr_pp[1], [50, 60], pl=false)
 f = dec_const_pcvc(corr_pp[1], [50, 60], m, pl=false)
+
+#left and right correlators
+f_ratio = dec_const_pcvc(ppL, ppR, [50, 60], m)
 ```
 """
 function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false,
@@ -340,23 +349,18 @@ function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu
     R_av = plat_av(R, plat, wpm)
     f = sqrt(2) * (mu[1] + mu[2]) *R_av / m^1.5
     if pl 
-        if isnothing(wpm)
-            uwerr(f)
-            uwerr(R_av)
-            uwerr.(R)
-        else
-            uwerr(f, wpm)
-            uwerr(R_av, wpm)
-            [uwerr(Raux, wpm) for Raux in R]
-        end
-    
+        isnothing(wpm) ? uwerr(f) : uwerr(f, wpm)
+        isnothing(wpm) ? uwerr(R_av) : uwerr(R_av, wpm)
         v = value(R_av)
-        e = err(R_av)                                
+        e = err(R_av)  
+
         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")
+        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(1:length(R), value.(R), err.(R), fmt="x", color="black", label=lbl)
         ylabel(L"$R$")
         xlabel(L"$x_0$")
+        legend()
         title(string(L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2]))
         display(gcf())          
     end 
@@ -373,6 +377,50 @@ function dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=tru
     return dec_const_pcvc(corr.obs, plat, m, corr.mu, corr.y0, pl=pl, data=data, wpm=wpm)
 end
 
+function dec_const_pcvc(ppL::Vector{uwreal}, ppR::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)
+
+    T = length(ppL)
+    f1 = [ppL[T - y0] for k = 1:T]
+
+    R = ((ppL .* ppR ./ f1).^2).^(1/4)
+    R = R[2:end-1]
+    R_av = plat_av(R, plat, wpm)
+    f = sqrt(2) * (mu[1] + mu[2]) * R_av / m^1.5
+    if pl
+        isnothing(wpm) ? uwerr(f) : uwerr(f, wpm)
+        isnothing(wpm) ? uwerr(R_av) : uwerr(R_av, wpm)
+        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(1:length(R), value.(R), err.(R), fmt="x", color="black", label=lbl)
+        ylabel(L"$R$")
+        xlabel(L"$x_0$")
+        legend()
+        title(string(L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2]))
+        display(gcf())   
+    end
+
+    if !data
+        return f
+    else
+        return (f, uwdot(sqrt(2) * (mu[1] + mu[2]) / m^1.5, R))
+    end
+end
+
+function dec_const_pcvc(corrL::Corr, corrR::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false,
+    wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+
+    T = length(corrL.obs)
+    if (corrL.kappa == corrR.kappa) && (corrL.mu == corrR.mu) && (corrL.y0 == T - 1 - corrR.y0)
+        return dec_const_pcvc(corrL.obs, corrR.obs, plat, m, corrL.mu, corrL.y0, pl=pl, data=data, wpm=wpm)
+    else
+        error("y0, kappa or mu values does not match")
+    end
+end
 #t0
 function get_model(x, p, n)
     s = 0.0