using juobs, ADerrors, DelimitedFiles, PyPlot, LaTeXStrings
const path = "/home/javier/Lattice/charm/production_2"
const path_plat = "/home/javier/Lattice/juobs/analysis/plat.txt"

const ensembles = ["H400", "N202", "N200", "N203", "N300", "J303"]
const deg = [true, true, false, false, true, false]
const L = [32, 48, 48, 48, 48, 64]
const beta = [3.46 , 3.55, 3.55, 3.55, 3.70, 3.70]
const R = [["H400r001", "H400r002"], "N202r001", ["N200r000", "N200r001"], ["N203r000", "N203r001"], "N300r002", "J303r003"]
include("/home/javier/Lattice/juobs/analysis/functions.jl")
include("/home/javier/Lattice/juobs/constants/juobs_const.jl")



m_ll = Vector{uwreal}(undef, length(ensembles))
m_ss = Vector{uwreal}(undef, length(ensembles))
m_hh = Vector{Vector{uwreal}}(undef, length(ensembles))
mu_pp = Vector{Vector{Vector{Float64}}}(undef, length(ensembles))

for iens = 1:length(ensembles)
    pp = read_dat(R[iens], "G5", "G5")
    a0p = read_dat(R[iens], "G5", "G0G5")
    rew = read_rew(R[iens])
    pp_obs = corr_obs.(pp, L=L[iens], rw=rew)
    a0p_obs = corr_obs.(a0p, L=L[iens], rw=rew)

    mu_pp[iens] = getfield.(pp_obs, :mu)

    m = comp_pcac(a0p_obs, pp_obs, deg[iens], ensembles[iens])

    m_ll[iens] = get_ll(mu_pp[iens], m, deg[iens])
    m_ss[iens] = deg[iens] ? m_ll[iens] : get_ss(mu_pp[iens], m, deg[iens])
    m_hh[iens] = get_hh(mu_pp[iens], m, deg[iens])
end

mm = get_mu.(mu_pp, deg)
mul_pp = getindex.(mm, 1)
mus_pp = getindex.(mm, 2)
muh_pp = getindex.(mm, 3)

cot_ll = Vector{uwreal}(undef, length(ensembles))
cot_ss = Vector{uwreal}(undef, length(ensembles))
cot_hh = Vector{Vector{uwreal}}(undef, length(ensembles))

for iens = 1:length(ensembles)
    cot_ll[iens] = za(beta[iens]) * m_ll[iens] / mul_pp[iens]
    cot_ss[iens] = za(beta[iens]) * m_ss[iens] / mus_pp[iens]
    cot_hh[iens] = fill(za(beta[iens]), 3) .* m_hh[iens] ./ muh_pp[iens]
end