analysis

.gitignore

@@ -2,3 +2,5 @@
 analysis/plat\.txt
 
 *.pdf
+
+analysis/*.txt

analysis/analysis_pcac.jl

+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
\ No newline at end of file

analysis/functions.jl

@@ -8,14 +8,19 @@ function read_dat(rep::Vector{String}, g1::String="G5", g2::String="G5")
     p = joinpath.(path, rep, "info")
     f = [filter(x-> contains(x, ".dat"), readdir(p[k]))[1] for k = 1:length(p)] 
     f = joinpath.(p, f)
-    aux = read_mesons.(f, g1, g2)
-    res = []
-    for j = 1:length(aux[1])
-        aux2 = [aux[i][j] for i = 1:length(aux)]
-        push!(res, aux2)
+    return read_mesons(f, g1, g2)
+end
+function read_rew(rep::String)
+    p = joinpath(path, rep, "rew")
+    f = filter(x-> contains(x, ".dat"), readdir(p))
+    f = joinpath(p, f[1])
+    try
+        return read_ms1(f)
+    catch
+        return read_ms1(f, v="1.4")
     end
-    return res
 end
+read_rew(rep::Vector{String}) = read_rew.(rep)
 
 function get_mu(mu_list::Vector{Vector{Float64}}, deg::Bool)
     mu_sorted = unique(sort(minimum.(mu_list)))
@@ -48,6 +53,16 @@ function get_ls(mu_list, obs::Vector{uwreal}, deg::Bool)
     end
 end
 
+function get_ss(mu_list, obs::Vector{uwreal}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mus in mu && mu[1] == mu[2] #s-s
+            return obs[k]
+        end
+    end
+end
+
 function get_lh(mu_list, obs::Vector{uwreal}, deg::Bool)
     mul, mus, muh = get_mu(mu_list, deg)
     obs_lh = Vector{uwreal}(undef, 0)
@@ -72,6 +87,17 @@ function get_sh(mu_list, obs::Vector{uwreal}, deg::Bool)
     return obs_sh
 end
 
+function get_hh(mu_list, obs::Vector{uwreal}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    obs_hh = Vector{uwreal}(undef, 0)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mu[1] == mu[2] && !(mul in mu) && !(mus in mu)
+            push!(obs_hh, obs[k])
+        end
+    end
+    return obs_hh
+end
 function select_plateau(ens::String, mu_list, deg::Bool)
     mul, mus, muh = get_mu(mu_list, deg)
     f = readdlm(path_plat)
@@ -116,6 +142,11 @@ function comp_meff(pp_obs::Vector{juobs.Corr}, deg::Bool, ens::String; pl::Bool=
     plat = select_plateau(ens, mu_list, deg)
     return meff.(pp_obs, plat, pl=pl)
 end
+function comp_pcac(a0p_obs::Vector{juobs.Corr}, pp_obs::Vector{juobs.Corr}, deg::Bool, ens::String; pl::Bool=false)
+    mu_list = getfield.(pp_obs, :mu)
+    plat = select_plateau(ens, mu_list, deg)
+    return mpcac.(a0p_obs, pp_obs, plat, pl=pl)
+end
 function comp_f(pp_obs::Vector{juobs.Corr}, m::Vector{uwreal}, deg::Bool, ens::String; pl::Bool=false)
     mu_list = getfield.(pp_obs, :mu)
     plat = select_plateau(ens, mu_list, deg)