mpcac added

src/juobs.jl

@@ -11,6 +11,6 @@ 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 meff, dec_const_pcvc, comp_t0
+export meff, mpcac, dec_const_pcvc, comp_t0
 
 end # module

src/juobs_obs.jl

@@ -47,6 +47,72 @@ end
 meff(corr::Corr, plat::Vector{Int64}; pl::Bool=true, data::Bool=false) = 
     meff(corr.obs, plat, pl=pl, data=data, mu=corr.mu)
 
+@doc raw"""
+    mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false, mu::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+
+    mpcac(a0p::Corr, pp::Corr, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false)
+
+Computes PCAC mass for a given correlator `a0p` and `pp` 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. The `ca` allows to compute `mpcac` using
+the improved axial current.
+
+```@example
+data_pp = read_mesons(path, "G5", "G5")
+data_a0p = read_mesons(path, "G5", "G0G5")
+corr_pp = corr_obs.(data_pp)
+corr_a0p = corr_obs.(data_a0p)
+m12 = mpcac(corr_a0p, corr_pp, [50, 60], pl=false)
+
+p0 = 9.2056
+p1 = -13.9847
+g2 = 1.73410
+ca = -0.006033 * g2 *( 1 + exp(p0 + p1/g2))
+
+m12 = mpcac(corr_a0p, corr_pp, [50, 60], pl=false, ca=ca)
+```
+"""
+function mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false, 
+    mu::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)     
+    corr_a0p = -a0p[2:end-1] 
+    corr_pp = pp[2:end-1]
+
+    der_a0p = (corr_a0p[3:end] .- corr_a0p[1:end-2]) / 2 
+    if ca != 0.0
+        der2_pp = corr_pp[1:end-2] + corr_pp[3:end] - 2 * corr_pp[2:end-1]
+        der_a0p = der_a0p + ca * der2_pp
+    end
+
+    aux = der_a0p ./ (2 .* corr_pp[2:end-1])
+    mass = plat_av(aux, plat, wpm)                                                                            
+    if pl == true
+        isnothing(wpm) ? uwerr(mass) : uwerr(mass, wpm)                       
+        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{PCAC}$")
+        xlabel(L"$x_0$")
+    
+        if !isnothing(mu)
+            title(string(L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2]))
+        end
+        display(gcf())
+    end                                                  
+    if data == false                   
+        return mass                                        
+    else               
+        return (mass, aux)     
+    end                    
+end   
+mpcac(a0p::Corr, pp::Corr, plat::Vector{Int64}; ca::Float64=0.0, pl::Bool=true, data::Bool=false) =
+    mpcac(a0p.obs, pp.obs, plat, ca=ca, pl=pl, data=data, mu=a0p.mu)
 ## Decay constants                                                                                       
 @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)