@doc raw"""
    meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)     

    meff(corr::Corr, plat::Vector{Int64}; pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)

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, 
    mu::Union{Vector{Float64}, Nothing}=nothing, kappa::Union{Vector{Float64}, Nothing}=nothing,
    wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, savepl::Bool=false,
    label::Union{String, Nothing}=nothing)     
    dim = length(corr)                                                                                    
    aux = 0.5 .* log.((corr[2:dim-2] ./ corr[3:dim-1]).^2)     
    mass = plat_av(aux, plat, wpm)                                                                            
    if pl || savepl
        isnothing(wpm) ? uwerr(mass) : uwerr(mass, wpm)                       
        x = 1:length(aux)
        y = abs.(value.(aux))
        dy = err.(aux)
        v = abs(value(mass))
        e = err(mass)
        
        filename = "meff/";
        figure()
        fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75)
        l,=errorbar(x, y, dy, fmt="x", color="black")
        if !isnothing(label)
            l.set_label(label)
            legend()
            filename*=label*"_"
        end
        ylabel(L"$m_\mathrm{eff}$")
        xlabel(L"$x_0$")
        _, max_idx = findmax(value.(corr))
        ylim(v-80*e, v+80*e)
        xlim(left=max_idx)
        
    
        if !isnothing(kappa)
            title(string(L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
            filename*="kappa_1_$(kappa[1])_kappa_2_$(kappa[2]).png"
        end
        if !isnothing(mu)
            title(string(L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2]))
            filename*="mu_1_$(mu[1])_mu_2_$(mu[2]).png"
        end
        if savepl
            tight_layout()
            gcf().set_size_inches(10,5)
            savefig(filename)
            close(gcf())
        else
            display(gcf())
        end
        # close("all")
    end                                                  
    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, savepl::Bool=false,
    label::Union{String,Nothing}=nothing)
    
    if corr.mu == [0.0, 0.0]
        return meff(corr.obs, plat, pl=pl, data=data, kappa=corr.kappa, mu=nothing, wpm=wpm,savepl=savepl, label=label)
    else
        return meff(corr.obs, plat, pl=pl, data=data, mu=corr.mu, kappa=nothing, wpm=wpm,savepl=savepl, label=label)
    end
end
@doc raw"""
    mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, 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)

    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)

Computes the bare 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` variable 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, 
    kappa::Union{Vector{Float64}, Nothing}=nothing, mu::Union{Vector{Float64}, Nothing}=nothing, 
    wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,
    savepl::Bool = false, label::Union{String,Nothing}=nothing)     
    
    corr_a0p = -a0p[1:end] 
    corr_pp = pp[1:end]

    if ca != 0.0
        der_a0p = (corr_a0p[3:end] .- corr_a0p[1:end-2]) / 2 
        der2_pp = (corr_pp[1:end-4] - 2*corr_pp[3:end-2] + corr_pp[5:end])/4
        der_a0p = der_a0p[2:end-1] + ca * der2_pp
    else
        der_a0p = (corr_a0p[4:end-1] .- corr_a0p[2:end-3]) / 2 
    end

    #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[3:end-2])
    mass = plat_av(aux, plat, wpm)                                                                            
    if pl || savepl
        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)
        l, = errorbar(x, y, dy, fmt="x", color="black")
        filename = "mpcac/"
        if !isnothing(label)
            l.set_label(label)
            legend()
            filename *= label*"_"
        end
        ylabel(L"$m_\mathrm{PCAC}$")
        xlabel(L"$x_0$")
        _, max_idx = findmax(value.(pp))
        ylim(v-20*e, v+20*e)
        xlim(left=max_idx)
        
        if !isnothing(kappa)
            title(string(L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
            global filename *= "kappa_1_$(kappa[1])_kappa_2_$(kappa[2]).png"
        end
        if !isnothing(mu)
            title(string(L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2]))
            global filename *="mu_1_$(mu[1])_mu_2_$(mu[2]).png"
        end
        if savepl
            tight_layout()
            gcf().set_size_inches(10, 5)
            savefig(filename, dpi=100)
            savefig(filename)
            close(gcf())
        elseif pl
            display(gcf())
        end
    end                                                  
    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,
    savepl::Bool=false, label::Union{String,Nothing}=nothing)
    
    if a0p.kappa == pp.kappa || a0p.mu == pp.mu
        if a0p.mu == [0.0, 0.0]
            return mpcac(a0p.obs, pp.obs, plat, ca=ca, pl=pl, data=data, kappa=a0p.kappa, mu=nothing, wpm=wpm,savepl=davepl, label=label)
        else
            return mpcac(a0p.obs, pp.obs, plat, ca=ca, pl=pl, data=data, mu=a0p.mu, kappa=nothing, wpm=wpm,savepl=savepl,label=label)
        end
    else
        error("mu or kappa values does not match")
    end
end
#= 
@doc raw"""
    mpcac_ren(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}, za_zp::uwreal; ca::Float64=0.0,ba_bp_z::uwreal=0.0,pl::Bool=true, data::Bool=false,
            kappa::Union{Vector{Float64}, Nothing}=nothing, mu::Union{Vector{Float64}, Nothing}=nothing, 
            wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, savepl::Bool =false)
    
    mpcac_ren(a0p::Corr, pp::Corr, plat::Vector{Int64}, za_zp::uwreal; ca::Float64,ba_bp_z::uwreal, pl::Bool=true, data::Bool=false,
            wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,savepl::Bool =false)

Computes the renormalized 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`, `savepl` and `data` allow to show the plots, save it as a png file and return data as an extra result. The `ca` variable allows to compute `mpcac` using
the improved axial current. The `ba_bp_z` variable allows to compute `mpcac` using the renormalized improved axial and pseudoscalar currents. 
Without specifying `ba_bp_z` the results is equivalent to the output of `mpcac(a0p,pp,plat, ca=ca)` multiplied by `za_zp`

```@example
include("juobs_const.jl")
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)
za_zp = za(beta)/ZP(beta) 

m12 = mpcac_ren(corr_a0p, corr_pp, [50, 60], za_zp, pl=false)
#equivalent to m12 = za_zp.*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_ren(corr_a0p, corr_pp, [50, 60],za_zp, pl=false, ca=ca)
#equivalent to m12 = za_zp.*mpcac(corr_a0p,corr_pp,[50,60],pl=false,ca=ca)


ba_bp_z = ba_bp(beta)*Z(beta)
m12 = mpcac_ren = (corr_a0p,corr_pp,[50,60],za_zp,pl=false, ca=ca, ba_bp_z=ba_bp_z)

```
"""
function mpcac_ren(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca::Float64=0.0,ba_bp_z::uwreal=uwreal(0.0),
    pl::Bool=true, data::Bool=false,kappa::Union{Vector{Float64}, Nothing}=nothing, mu::Union{Vector{Float64}, Nothing}=nothing, 
    wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,savepl::Bool =false, label::Union{String,Nothing}=nothing)
    corr_a0p = -a0p[1:end] 
    corr_pp = pp[1:end]
    der_a0p = (corr_a0p[4:end-1] .- corr_a0p[2:end-3]) / 2  # \de_0 a0p
    m_pcac = der_a0p ./(2 * corr_pp[3:end-2]) # Unimproved, unrenormalized mpcac

    if ba_bp_z !=0.0 
        der_a0p .*= 1.0 .+ ba_bp_z*m_pcac  #improvement with ba_bp
    end

    if ca !=0.0
        der2_pp = (corr_pp[1:end-4] .- 2*corr_pp[3:end-2] .+ corr_pp[5:end])/4 
        der_a0p += ca*der2_pp  
    end
    
    m_ren = (der_a0p./(2*corr_pp[3:end-2]))
    mass = plat_av(m_ren,plat,wpm)

    if pl || savepl
        isnothing(wpm) ? uwerr(mass) : uwerr(mass, wpm)                       
        x = 1:length(m_ren)
        y = value.(m_ren)
        dy = err.(m_ren)
        v = value(mass)
        e = err(mass)
        
        plot= figure()
        filename="mpcac/"
        fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75)
        l,= errorbar(x, y, dy, fmt="x", color="black")
        if !isnothing(label)   
            l.set_label(label)
            legend()
            filename*=label*"_"
        end
        ylabel(L"$m^R_\mathrm{PCAC}$")
        xlabel(L"$x_0$")
        _, max_idx = findmax(value.(pp))
        ylim(v-20*e, v+20*e)
        xlim(left=max_idx)
        if !isnothing(kappa)
            title(string(L"$m_\mathrm{PCAC}^R$ $\mathcal{O(a)}$ improved","\n",L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
            filename *="kappa_1_$(kappa[1])_kappa_2_$(kappa[2])_renormalize_improved.png"
            
        end
        if !isnothing(mu)
            title(string(L"$m_\mathrm{PCAC}^R$ $\mathcal{O(a)}$ improved","\n",L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2]))
            filename *="mu_1_$(mu[1])_mu_2_$(mu[2])_renomalized_improved.png"
        end

        if savepl
            tight_layout()
            gcf().set_size_inches(10,5)
            savefig(filename)
        end
        if pl && !savepl
            display(gcf())
        else
            close(plot)
        end
    end  
    if !data                   
        return mass                                        
    end                
    return (mass, m_ren)  

end

function mpcac_ren(a0p::Corr, pp::Corr, plat::Vector{Int64}; ca::Float64=0.0,ba_bp_z::uwreal=uwreal(0.0), pl::Bool=true, data::Bool=false, 
    wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,savepl::Bool =false,label::Union{String,Nothing}=nothing)
    if a0p.kappa == pp.kappa || a0p.mu == pp.mu
        if a0p.mu == [0.0, 0.0]
            return mpcac_ren(a0p.obs, pp.obs, plat, ca=ca,ba_bp_z=ba_bp_z, pl=pl, data=data, kappa=a0p.kappa, mu=nothing, wpm=wpm, savepl=savepl,label=label)       
        else
            return mpcac_ren(a0p.obs, pp.obs, plat, ca=ca,ba_bp_z=ba_bp_z, pl=pl, data=data, mu=a0p.mu, kappa=nothing, wpm=wpm, savepl=savepl,label=label)
        end 
    else
        error("mu or kappa values does not match")
    end
end
 =#
## Decay constants
@doc raw"""
    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, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)

    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.

The flags `pl` and `data` allow to show the plots and return data as an extra result. The `ca` variable allows to compute `dec_const` using
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)

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
p0 = 9.2056
p1 = -13.9847
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,
    kappa::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]
    
    if ca != 0.0
        der_pp = (corr_pp[3:end] - corr_pp[1:end-2]) / 2
        corr_a0p = corr_a0p[2:end-1] + ca * der_pp
        T = length(corr_a0p)
        aux = exp.((collect(1:T) .- (T)/2) .* [m for k = 1:T])
    else
        T = length(corr_a0p)
        aux = exp.((collect(1:T) .- T/2) .* [m for k = 1:T])
    end

    R = corr_a0p .* aux ./ [((corr_pp[T-y0])^2)^(1/4) for k = 1:length(corr_a0p)]
    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()
        ylim(v-10*e, v+10*e)
        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())
        close()
    end
    if !data
        return f
    else
        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) 
    
    if (a0p.y0 == pp.y0) && (a0p.kappa == pp.kappa)
        return dec_const(a0p.obs, pp.obs, plat, m, a0p.y0, ca=ca, kappa=a0p.kappa, pl=pl, data=data, wpm=wpm)
    else 
        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
## ADD DOCUMENTATION FOR VECTOR DECAY 
function dec_const(vv::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; 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_vv = vv[2:end-1]
    T = length(corr_vv)

    aux = exp.((collect(1:T) .- y0  ) .* fill(m, T))

    R = ((aux .* corr_vv).^2).^0.25
    R_av = plat_av(R, plat, wpm)
    f = sqrt(2 / m)  * R_av
    R .*= sqrt.(2 ./ [m for i in 1:length(R)])
    if pl
        uwerr.(R)
        R_av *= sqrt(2/m)
        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()
        xlim(left=y0)
        ylim(v-10*e, v+10*e)
        ylabel(L"$af_{B}$")
        #ylabel(L"$R_\mathrm{av}$")
        xlabel(L"$x_0$")
        #title(L"$f_{B^*}$")
        if !isnothing(kappa)
            title(string( L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
        end
        display(gcf())
        #t = "ps_decay_$(kappa[1])_$(kappa[2]).pdf"
        #savefig(joinpath("/Users/ale/Il mio Drive/phd/secondment/3pf test/analysis/plots",t))
    end
    if !data
        return f
    else
        return f, R
    end
end
function dec_const(vv::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false,
    wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) 
    
    return dec_const(vv.obs, plat, m, vv.y0, kappa=vv.kappa, pl=pl, data=data, wpm=wpm)
end
# this can be used also for pseudoscalar decay constants with wilson fermions 
function dec_const_w(vv::Vector{uwreal}, pp::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_vv = vv[2:end-1]
    corr_pp = pp[2:end-1]
    T = length(corr_vv)
    if ca != 0.0
        der_pp = (corr_pp[3:end] - corr_pp[1:end-2]) / 2
        corr_vv = corr_vv[2:end-1] + ca* der_pp
        aux = exp.((collect(1:T-2) .- y0   ) .* fill(m, T-2))
    else
        #corr_vv = corr_vv[2:end-1]
        aux = exp.((collect(1:T) .- y0  ) .* fill(m/2, T))
    end
  

    R = (aux .* corr_vv ./ sqrt.(corr_pp))
    R_av = plat_av(R, plat, wpm)
    f = sqrt(2 / m)  * R_av
    if pl
        R .*= sqrt.(2 ./ [m for i in 1:length(R)])
        uwerr.(R)
        R_av *= sqrt(2/m)
        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()
        xlim(y0, y0 +plat[2] )
        ylim(v-10*e, v+10*e)
        ylabel(L"$af_{vec}$")
        #ylabel(L"$R_\mathrm{av}$")
        xlabel(L"$x_0$")
        #title(L"$f_{B^*}$")
        if !isnothing(kappa)
            title(string( L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
        end
        display(gcf())
        #t = "ps_decay_$(kappa[1])_$(kappa[2]).pdf"
        #savefig(joinpath("/Users/ale/Il mio Drive/phd/secondment/3pf test/analysis/plots",t))
    end
    if !data
        return f
    else
        return f, R
    end
end
function dec_const_w(vv::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) 
    
    return dec_const_w(vv.obs, pp.obs, plat, m, vv.y0, ca=ca, kappa=vv.mu, pl=pl, data=data, wpm=wpm)
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)

    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 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,
    wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
    """
    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).^2).^0.25
    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$")
        xlim(y0,y0+plat[2])
        ylim(v-10*e, v+10*e)
        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(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) 
    
    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
    for k = 1:n
        s = s .+ p[k] .* x.^(k-1)
    end
    return s
end
@doc raw"""
    comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Matrix{Float64}, Nothing}=nothing, npol::Int64=2, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, info::Bool=false)

    comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Vector{Matrix{Float64}}, Nothing}=nothing, npol::Int64=2, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, info::Bool=false)

Computes `t0` using the energy density of the action `Ysl`(Yang-Mills action).
`t0` is computed in the plateau `plat`.
A polynomial interpolation in `t` is performed to find `t0`, where `npol` is the degree of the polynomial (linear fit by default)

The flag `pl` allows to show the plot. 

The flag `info` provides extra output that contains information about the primary observables. The function returns the primary observables ``<WY>`` and ``<W>``
(it returns the observable <Y> if rw=nothing)

```@example
#Single replica
Y = read_ms(path)
rw = read_ms(path_rw)

t0, Yobs = comp_t0(Y, [38, 58], L=32, info=true)
t0_r, WYobs, Wobs = comp_t0(Y, [38, 58], L=32, rw=rw, info=true)

#Two replicas
Y1 = read_ms(path1)
Y2 = read_ms(path2)
rw1 = read_ms(path_rw1)
rw2 = read_ms(path_rw2)

t0 = comp_t0([Y1, Y2], [38, 58], L=32, pl=true)
t0_r = comp_t0(Y, [38, 58], L=32, rw=[rw1, rw2], pl=true)

```
"""
function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, 
    rw::Union{Matrix{Float64}, Nothing}=nothing, npol::Int64=2, ws::ADerrors.wspace=ADerrors.wsg,
    wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, info::Bool=false)

    Ysl = Y.obs
    t = Y.t
    id = Y.id
    replica = size.([Ysl], 1)

    #Truncation
    if id in keys(ADerrors.wsg.str2id)
        n_ws = findfirst(x-> x == ws.str2id[id], ws.map_nob)
        if !isnothing(n_ws)
            ivrep_ws = ws.fluc[n_ws].ivrep

            if length(ivrep_ws) != 1
                error("Different number of replicas")
            end

            if replica[1] > ivrep_ws[1]
                println("Automatic truncation in Ysl ", ivrep_ws[1], " / ", replica[1], ". R = 1")
                Ysl = Ysl[1:ivrep_ws[1], :, :]
            elseif replica[1] < ivrep_ws[1]
                println(replica[1])
                error("Automatic truncation failed. R = 1\nTry using truncate_data!")
            end
        end
    end

    xmax = size(Ysl, 2)
    nt0 = t0_guess(t, Ysl, plat, L)
    
    dt0 = iseven(npol) ? Int64(npol / 2) : Int64((npol+1)/ 2)
    Y_aux = Matrix{uwreal}(undef, xmax, 2*dt0+1)
    
    if !isnothing(rw)
        Ysl_r, W = apply_rw(Ysl, rw)
        W_obs = uwreal(W, id)
        WY_aux = Matrix{uwreal}(undef, xmax, 2*dt0+1)
    end

    for i = 1:xmax
        k = 1
        for j = nt0-dt0:nt0+dt0
            if isnothing(rw)
                Y_aux[i, k] = uwreal(Ysl[:, i, j], id)
            else
                WY_aux[i, k] = uwreal(Ysl_r[:, i, j], id)
                Y_aux[i, k] = WY_aux[i, k] / W_obs
            end
            k = k + 1
        end
    end
    x = t[nt0-dt0:nt0+dt0]
    t2E = [plat_av(Y_aux[:, j], plat, wpm) for j=1:2*dt0+1] .* x.^2 / L^3
    
    model(x, p) = get_model(x, p, npol)

    par, chi  = fit_routine(model, x, t2E, npol)
    fmin(x, p) = model(x, p) .- 0.3
    t0 = root_error(fmin, t[nt0], par)
    if pl
        if isnothing(wpm)
            uwerr(t0)
            uwerr.(t2E)
        else
            uwerr(t0, wpm)
            [uwerr(t2E_aux, wpm) for t2E_aux in t2E]
        end

        v = value.(t2E)
        e = err.(t2E)

        figure()
        errorbar(x, v, e, fmt="d", capsize=2, mfc="none", color="black")
        errorbar(value(t0), 0.3, xerr=err(t0), capsize=2, color="red", fmt="d")
        xarr = Float64.(range(5.0, 5.2, length=100))
        yarr = model(xarr, par); uwerr.(yarr)
        fill_between(xarr, value.(yarr) .- err.(yarr), value.(yarr) .+ err.(yarr), color="royalblue", alpha=0.2 )
        ylabel(L"$t^2\langle E(x_0, t)\rangle$")
        xlabel(L"$t/a^2$")
        tight_layout()
        savefig("/Users/alessandroconigli/Desktop/t0_interp.pdf")
        display(gcf())

        figure()
        t_pl = dt0 + 1
        yyy = Y_aux[:, t_pl] * t[nt0]^2 / L^3 ; uwerr.(yyy)
        # plot(collect(1:xmax), value.(Y_aux[:, t_pl]) * t[nt0]^2 / L^3, "x") 
        errorbar(collect(1:xmax), value.(yyy), err.(yyy), fmt="d", mfc="none", color="black", capsize=2) 
        fill_between(collect(plat[1]:plat[2]), v[t_pl]+e[t_pl], v[t_pl]-e[t_pl], alpha=0.75, color="royalblue")
        ylabel(L"$t^2\langle E(x_0, t)\rangle$")
        xlabel(L"$x_0/a$")
        ylim(v[t_pl]-30*e[t_pl], v[t_pl]+50*e[t_pl])
        title(string(L"$t/a^2 = $", t[nt0]))
        tight_layout()
        savefig("/Users/alessandroconigli/Desktop/t0_plateau.pdf")
        display(gcf())
    end
    if info && !isnothing(rw)
        return (t0, WY_aux, W_obs)
    elseif info && isnothing(rw)
        return (t0, Y_aux)
    else
        return t0
    end
end

function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false, 
    rw::Union{Vector{Matrix{Float64}}, Nothing}=nothing, npol::Int64=2, ws::ADerrors.wspace=ADerrors.wsg,
    wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, info::Bool=false)

    nr = length(Y)
    Ysl = getfield.(Y, :obs)
    t = getfield.(Y, :t)
    t = t[1]
    id = getfield.(Y, :id)
    replica = size.(Ysl, 1)

    if !all(id .== id[1])
        error("IDs are not equal")
    end
    id = id[1]
    #Truncation
    if id in keys(ADerrors.wsg.str2id)
        n_ws = findfirst(x-> x == ws.str2id[id], ws.map_nob)
        if !isnothing(n_ws)
            ivrep_ws = ws.fluc[n_ws].ivrep

            if length(replica) != length(ivrep_ws)
                error("Different number of replicas")
            end

            for k = 1:length(replica)
                if replica[k] > ivrep_ws[k]
                    println("Automatic truncation in Ysl ", ivrep_ws[k], " / ", replica[k], ". R = ", k)
                    Ysl[k] = Ysl[k][1:ivrep_ws[k], :, :]
                elseif replica[k] < ivrep_ws[k]
                    error("Automatic truncation failed. R = ", replica[k], "\nTry using truncate_data!")
                end
            end
            replica = size.(Ysl, 1)
        end
    end

    tmp = Ysl[1]
    [tmp = cat(tmp, Ysl[k], dims=1) for k=2:nr]
    nt0 = t0_guess(t, tmp, plat, L)
    xmax = size(tmp, 2)

    dt0 = iseven(npol) ? Int64(npol / 2) : Int64((npol+1) / 2)
    Y_aux = Matrix{uwreal}(undef, xmax, 2*dt0+1)

    if !isnothing(rw)
        Ysl_r, W = apply_rw(Ysl, rw)
        tmp_r = Ysl_r[1]
        tmp_W = W[1]
        [tmp_r = cat(tmp_r, Ysl_r[k], dims=1) for k = 2:nr]
        [tmp_W = cat(tmp_W, W[k], dims=1) for k = 2:nr]
        W_obs = uwreal(tmp_W, id, replica)
        WY_aux = Matrix{uwreal}(undef, xmax, 2*dt0+1)
    end
    for i = 1:xmax
        k = 1
        for j = nt0-dt0:nt0+dt0
            if isnothing(rw)
                Y_aux[i, k] = uwreal(tmp[:, i, j], id, replica)
            else
                WY_aux[i, k] = uwreal(tmp_r[:, i, j], id, replica)
                Y_aux[i, k] = WY_aux[i, k] / W_obs
            end
            k = k + 1
        end
    end
    x = t[nt0-dt0:nt0+dt0]
    t2E = [plat_av(Y_aux[:, j], plat, wpm) for j=1:2*dt0+1] .* x.^2 / L^3
    
    model(x, p) = get_model(x, p, npol)

    par, chi = fit_routine(model, x, t2E, npol)
    fmin(x, p) = model(x, p) .- 0.3
    t0 = root_error(fmin, t[nt0], par)
    if pl
        if isnothing(wpm)
            uwerr(t0)
            uwerr.(t2E)
        else
            uwerr(t0, wpm)
            [uwerr(t2E_aux, wpm) for t2E_aux in t2E]
        end

        v = value.(t2E)
        e = err.(t2E)

        figure()
        errorbar(x, v, e, fmt="x")
        
        errorbar(value(t0), 0.3, xerr=err(t0), fmt="x")
        ylabel(L"$t^2E$")
        xlabel(L"$t/a^2$")
        display(gcf())

        figure()
        t_pl = dt0 + 1
        plot(collect(1:xmax), value.(Y_aux[:, t_pl]) * t[nt0]^2 / L^3, "x")
        fill_between(collect(plat[1]:plat[2]), v[t_pl]+e[t_pl], v[t_pl]-e[t_pl], alpha=0.75, color="green")
        ylabel(L"$t^2E(x0, t)$")
        xlabel(L"$x_0/a$")
        title(string(L"$t/a^2 = $", t[nt0]))
        display(gcf())
    end
    if info && !isnothing(rw)
        return (t0, WY_aux, W_obs)
    elseif info && isnothing(rw)
        return (t0, Y_aux)
    else
        return t0
    end
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
    #look for values t2E: t2E > 0.3 and 0.0 < t2E_ax < 0.3
    if !(any(t2E_ax .> 0.3) && any((t2E_ax .< 0.3) .& (t2E_ax .> 0.0))) 
        error("Error finding solutions. Check volume normalization!")
    end
    t0_aux = minimum(abs.(t2E_ax .- 0.3))
    nt0 = findfirst(x-> abs(x - 0.3) == t0_aux, t2E_ax) #index closest value to t0
    return nt0
end