juobs_obs.jl

@doc raw"""
meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bool=false )     

meff(corr::Corr, plat::Vector{Int64}; pl::Bool=true, data::Bool=false)

Computes effective mass for a given correlator corr 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.

```@example
data = read_mesons(path, "G5", "G5")
corr_pp = corr_obs.(data)
m = meff(corr_pp[1], [50, 60], pl=false)
```
"""
function meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bool=false)     
    dim = length(corr)                                                                                    
    aux = log.(corr[2:dim-2] ./ corr[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::Bool=true, data::Bool=false) = 
    meff(corr.obs, plat, pl=pl, data=data)

## 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)

dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false)

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.

```@example
data = read_mesons(path, "G5", "G5")
corr_pp = corr_obs.(data)
m = meff(corr_pp[1], [50, 60], pl=false)
f = dec_const_pcvc(corr_pp[1], [50, 60], m, pl=false)
```
"""
function dec_const_pcvc(corr::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 = corr[2:end-1]
    dim = length(corr_pp)
                                                                 
    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_av)
        v = value(R_av)
        e = err(R_av)                                
        uwerr.(R)
        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())          
    end 
    if data == false
        return f
    else
        return (f, R)                   
    end 
end 
dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false) =
    dec_const_pcvc(corr.obs, plat, m, corr.mu, corr.y0, pl=pl, data=data)