Merge branch 'javier'

.gitignore

+analysis/plots/*
+
+analysis/plat\.txt

Project.toml

@@ -8,5 +8,6 @@ ADerrors = "5e92007d-7bf1-471c-8ceb-4591b8b567a9"
 BDIO = "375f315e-f2c4-11e9-2ef9-134f02f79e27"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

analysis/functions.jl

+function read_dat(rep::String, g1::String="G5", g2::String="G5")
+    p = joinpath(path, rep, "info")
+    f = filter(x-> contains(x, ".dat"), readdir(p))
+    f = joinpath(p, f[1])
+    return read_mesons(f, g1, g2)
+end
+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)
+    end
+    return res
+end
+
+function get_mu(mu_list::Vector{Vector{Float64}}, deg::Bool)
+    mu_sorted = unique(sort(minimum.(mu_list)))
+
+    mul = mu_sorted[1]
+    deg ? mus = 0.0 : mus = mu_sorted[2]
+    muh = unique(maximum.(mu_list))
+    muh = filter(x-> x > mul && x > mus, muh)
+
+    return mul, mus, muh
+end
+
+function get_ll(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 mul in mu && mu[1] == mu[2] #l-l
+            return obs[k]
+        end
+    end
+end
+
+function get_ls(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 mul in mu && mus in mu && mu[1] != mu[2] #l-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)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mul in mu && mu[1] != mu[2] && !(mus in mu) #l-h
+            push!(obs_lh, obs[k])
+        end
+    end
+    return obs_lh
+end
+
+function get_sh(mu_list, obs::Vector{uwreal}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    obs_sh = Vector{uwreal}(undef, 0)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mus in mu && mu[1] != mu[2] && !(mul in mu)#s-h
+            push!(obs_sh, obs[k])
+        end
+    end
+    return obs_sh
+end
+
+function select_plateau(ens::String, mu_list, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    f = readdlm(path_plat)
+    head = String.(f[:, 1])
+
+    delim = findall(x-> contains(x, "#"), head)
+    edelim = findfirst(x-> contains(x, ens), head)
+    cdelim = findfirst(x-> x.== edelim, delim)
+
+    del = delim[cdelim]+1 : delim[cdelim+1]-1
+
+    plat = Vector{Vector{Int64}}(undef, 0)
+    plat_ = Int64.(f[del, 2:3])
+    head_ = String.(f[del, 1])
+    
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mul in mu && mu[1] != mu[2] && !(mus in mu)  #heavy-light
+            n = findfirst(x-> contains(x, "lh"), head_)
+            push!(plat, [plat_[n, 1], plat_[n, 2]])
+        elseif mul in mu && mu[1] == mu[2] #light-light
+            n = findfirst(x-> contains(x, "ll"), head_)
+            push!(plat, [plat_[n, 1], plat_[n, 2]])
+        elseif mul in mu && mus in mu#strange-light
+            n = findfirst(x-> contains(x, "ls"), head_)
+            push!(plat, [plat_[n, 1], plat_[n, 2]])
+        elseif mus in mu && mu[1] != mu[2] && !(mul in mu)#heavy-strange
+            n = findfirst(x-> contains(x, "sh"), head_)
+            push!(plat, [plat_[n, 1], plat_[n, 2]])
+        elseif mus in mu && mu[1] == mu[2] && !(mul in mu)#strange-strange
+            n = findfirst(x-> contains(x, "ss"), head_)
+            push!(plat, [plat_[n, 1], plat_[n, 2]])
+        elseif !(mul in mu) && !(mus in mu) #heavy-heavy
+            n = findfirst(x-> contains(x, "hh"), head_)
+            push!(plat, [plat_[n, 1], plat_[n, 2]])
+        end
+    end
+    return plat
+end
+function comp_meff(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 meff.(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)
+    return dec_const_pcvc.(pp_obs, plat, m, pl=pl)
+end
+
+function match_muc(muh, m_lh, m_lh_star, target)
+    M = (m_lh .+ 3 .* m_lh_star) ./ 4
+
+    par, chi2exp = lin_fit(muh, M)
+    muh_target = x_lin_fit(par, target)
+    return muh_target
+end
+function match_muc(muh, m_lh, m_sh, m_lh_star, m_sh_star, target)
+    M = (2 .* m_lh .+ 6 .* m_lh_star .+ m_sh .+ 3 .* m_sh_star) ./ 12
+
+    par, chi2exp = lin_fit(muh, M)
+    muh_target = x_lin_fit(par, target)
+    return muh_target
+end
\ No newline at end of file

analysis/results.txt

+(H400)$Z^{tm}_M \mu_c \sqrt{8t_0} = $3.030909671357194 +/- 0.07903701818131038
+(H400)$M_D \sqrt{8t_0}= $3.9807556606105736 +/- 0.05781703939545818
+(H400)$M_Ds \sqrt{8t_0}= $3.9807556606105736 +/- 0.05781703939545818
+(H400)$M_D* \sqrt{8t_0}= $4.278934408455874 +/- 0.05001764960908615
+(H400)$M_Ds* \sqrt{8t_0}= $4.278934408455874 +/- 0.05001764960908615
+(H400)$f_D \sqrt{8t_0}= $0.5171340154200049 +/- 0.0030551745985989815
+(H400)$f_Ds \sqrt{8t_0}= $0.5171340154200049 +/- 0.0030551745985989815
+(N200)$Z^{tm}_M \mu_c \sqrt{8t_0} = $3.0777083812680863 +/- 0.07529964071401364
+(N200)$M_D \sqrt{8t_0}= $3.9484624535288577 +/- 0.056300551744400275
+(N200)$M_Ds \sqrt{8t_0}= $4.0641686659916445 +/- 0.05571406258479428
+(N200)$M_D* \sqrt{8t_0}= $4.2322082817565505 +/- 0.05089413966187937
+(N200)$M_Ds* \sqrt{8t_0}= $4.366150683510803 +/- 0.050567123757276136
+(N200)$f_D \sqrt{8t_0}= $0.4781332278412574 +/- 0.002533317431472261
+(N200)$f_Ds \sqrt{8t_0}= $0.5253557430045498 +/- 0.0019233408904497612
+(N203)$Z^{tm}_M \mu_c \sqrt{8t_0} = $3.069739432147581 +/- 0.07467450492462731
+(N203)$M_D \sqrt{8t_0}= $3.9619597653172356 +/- 0.05581411314938938
+(N203)$M_Ds \sqrt{8t_0}= $4.027646770204338 +/- 0.055472228755110334
+(N203)$M_D* \sqrt{8t_0}= $4.253505640330456 +/- 0.050187214812839276
+(N203)$M_Ds* \sqrt{8t_0}= $4.326727158575913 +/- 0.049836834652424916
+(N203)$f_D \sqrt{8t_0}= $0.488543655769974 +/- 0.0026562881901046583
+(N203)$f_Ds \sqrt{8t_0}= $0.5157212901463102 +/- 0.0021384519222027373
+(N300)$Z^{tm}_M \mu_c \sqrt{8t_0} = $2.997384760490272 +/- 0.07895636913539895
+(N300)$M_D \sqrt{8t_0}= $3.94373962398667 +/- 0.059571174171691506
+(N300)$M_Ds \sqrt{8t_0}= $3.94373962398667 +/- 0.059571174171691506
+(N300)$M_D* \sqrt{8t_0}= $4.291279319094004 +/- 0.04988459701774591
+(N300)$M_Ds* \sqrt{8t_0}= $4.291279319094004 +/- 0.04988459701774591
+(N300)$f_D \sqrt{8t_0}= $0.5032296543520727 +/- 0.003757345679706973
+(N300)$f_Ds \sqrt{8t_0}= $0.5032296543520727 +/- 0.003757345679706973
+(J303)$Z^{tm}_M \mu_c \sqrt{8t_0} = $3.13490498080823 +/- 0.08328602171541571
+(J303)$M_D \sqrt{8t_0}= $3.970055446897706 +/- 0.062204612871266905
+(J303)$M_Ds \sqrt{8t_0}= $4.113876464562304 +/- 0.06232345158893755
+(J303)$M_D* \sqrt{8t_0}= $4.199854539735192 +/- 0.05475423487017079
+(J303)$M_Ds* \sqrt{8t_0}= $4.399899824892723 +/- 0.05477252163421222
+(J303)$f_D \sqrt{8t_0}= $0.4649182529771295 +/- 0.004855221613601069
+(J303)$f_Ds \sqrt{8t_0}= $0.5220832840476 +/- 0.0031391370517194982
+
+$Z^{tm}_M \mu_c \sqrt{8t_0}$ = 3.289196404345117 +/- 0.16908049991604063
+muc(MeV) = 1571.5428436124926 +/- 74.65754166879076
+$M_D \sqrt{8t_0}$ = 4.034261174024082 +/- 0.12045227374006602
+m_D(MeV) = 1927.5268174700511 +/- 51.13372706582599
+$M_Ds \sqrt{8t_0}$ = 4.270003547803183 +/- 0.12354622671856343
+m_Ds(MeV) = 2040.1620009329038 +/- 52.898332267818304
+$M_D* \sqrt{8t_0}$ = 4.103105481644236 +/- 0.10100061920706542
+m_D_star(MeV) = 1960.4198909335223 +/- 42.35821470980299
+$M_Ds* \sqrt{8t_0}$ = 4.498614355722275 +/- 0.09456340196942646
+m_Ds_star(MeV) = 2149.3897985442322 +/- 37.956898903502925
+$f_D \sqrt{8t_0}$ = 0.4363276333625155 +/- 0.017871060655147757
+f_D(MeV) = 208.47267398669123 +/- 8.883408580389085
+$f_Ds \sqrt{8t_0}$ = 0.5263085789505523 +/- 0.011256710293275131
+f_Ds(MeV) = 251.46460688360165 +/- 6.054227498288841

constants/juobs_const.jl

 using ADerrors
 #PDG
 const hc = 197.3269804 #MeV fm
-const M_values = [1869.65, 2010.26, 1968.34, 2112.2] #MD, MD*, MDs, MDs* (MeV)
-const M_error = [0.05, 0.05, 0.07, 0.4]
+const M_values = [1869.65, 2010.26, 1968.34, 2112.2, 2980.3, 3096.916] #MD, MD*, MDs, MDs*, \eta_c, J/\psi (MeV)
+const M_error = [0.05, 0.05, 0.07, 0.4, 1.2, 0.011]
 #1802.05243 
 const b_values = [3.40, 3.46, 3.55, 3.70, 3.85]
 const ZM_data = [1.9684, 1.9935, 2.0253, 2.0630, 2.0814]
@@ -20,13 +20,15 @@ const t0_ph_error = ones(1,1) .* 5e-3
 const C1 = zeros(5, 5)
 const C2 = zeros(5, 5)
 const C3 = zeros(4, 4)
-const C4 = zeros(4, 4)
-for i = 1:5
-    C1[i, i] = ZM_error[i] ^ 2
-    C2[i, i] = ZM_tm_error[i] ^ 2
-    if i <= 4
-        C3[i, i] = t0_error[i] ^ 2
-        C4[i, i] = M_error[i] ^ 2
+const C4 = zeros(6, 6)
+for i = 1:6
+    C4[i, i] = M_error[i] ^ 2
+    if i<=5
+        C1[i, i] = ZM_error[i] ^ 2
+        C2[i, i] = ZM_tm_error[i] ^ 2
+        if i <= 4
+            C3[i, i] = t0_error[i] ^ 2
+        end
     end
 
 end

gevp tools/tools.jl

+function read_dat(ens::String, g1::String="G5", g2::String="G5")
+    path = joinpath(path_data, ens)
+    rep = readdir(path, join=true)
+    aux = read_mesons.(rep, g1, g2)
+    res = []
+    for j =1:length(aux[1])
+        aux2 = [aux[i][j] for i =1:length(aux)]
+        push!(res, aux2)
+    end
+    return res
+end
+function get_mu(mu_list::Vector{Vector{Float64}}, deg::Bool)
+    mu_sorted = unique(sort(minimum.(mu_list)))
+
+    mul = mu_sorted[1]
+    deg ? mus = 0.0 : mus = mu_sorted[2]
+    muh = unique(maximum.(mu_list))
+    muh = filter(x-> x > mul && x > mus, muh)
+
+    return mul, mus, muh
+end
+function get_lh(mu_list, obs::Array{juobs.Corr}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    obs_lh = Array{juobs.Corr}(undef, 0)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mul in mu && mu[1] != mu[2] && !(mus in mu) #l-h
+            push!(obs_lh, obs[k])
+        end
+    end
+    return obs_lh
+end
+function get_lh(mu_list, obs::Vector{Vector{juobs.Corr}}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    obs_lh = Array{Array{juobs.Corr}}(undef, length(obs), length(muh))
+    for n_obs = 1:length(obs)
+        obs_lh[n_obs] = get_lh(mu_list, obs[n_obs], deg)
+    end
+    return obs_lh
+end
+function get_sh(mu_list, obs::Array{juobs.Corr}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    obs_sh = Array{juobs.Corr}(undef, 0)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        if mus in mu && mu[1] != mu[2] && !(mul in mu)#s-h
+            push!(obs_sh, obs[k])
+        end
+    end
+    return obs_sh
+end
+function get_sh(mu_list, obs::Vector{Vector{juobs.Corr}}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    obs_sh = Array{Array{juobs.Corr}}(undef, length(obs), length(muh))
+    for n_obs = 1:length(obs)
+        obs_sh[n_obs] = get_sh(mu_list, obs[n_obs], deg)
+    end
+    return obs_sh
+end
+function get_hh(mu_list, obs::Array{juobs.Corr}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    obs_hh = Array{juobs.Corr}(undef,0)
+    for k = 1:length(mu_list)
+        mu = mu_list[k]
+        for i =1:length(muh)
+            if muh[i] in mu && mu[1] == mu[2]#h-h
+                push!(obs_hh, obs[k])
+            end
+        end
+    end
+    return obs_hh
+end
+function get_hh(mu_list, obs::Vector{Vector{juobs.Corr}}, deg::Bool)
+    mul, mus, muh = get_mu(mu_list, deg)
+    obs_hh = Array{Array{juobs.Corr}}(undef, length(obs), length(muh))
+    for n_obs = 1:length(obs)
+        obs_hh[n_obs] = get_hh(mu_list, obs[n_obs], deg)
+    end
+    return obs_hh
+end
+mutable struct  infoMatt
+    mat::Vector{Matrix}
+    mu::Vector{Float64}
+    y0::Int64
+    function infoMatt(a11::Array{juobs.Corr}, a12::Array{juobs.Corr}, a22::Array{juobs.Corr})
+        mu = getfield(a11[1], :mu)
+        y0 = Int64(getfield(a11[1], :y0))
+        elem11 = (sum(a11[i].obs for i =1:length(a11)))
+        elem12 = (sum(a12[i].obs for i =1:length(a12)))
+        elem22 = (sum(a22[i].obs for i =1:length(a22)))
+        diag = [elem11, elem22]
+        subdiag =[ elem12 ]
+        #diag = [a11.obs, a22[1].obs.+a22[2].obs.+a22[3].obs]
+        #subdiag = [a12[1].obs.+a12[2].obs.+a12[3].obs]
+        mat = get_matrix(diag, subdiag)
+        return new(mat, mu, y0)
+    end
+end
+function comp_mat(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}})
+    mu = getfield.(tot11[1],:mu)
+    res = Array{infoMatt}(undef, length(mu))
+    for i = 1:length(mu)
+        a11 = [ tot11[j][i] for j = 1:size(tot11,1) ]
+        a12 = [ tot12[j][i] for j = 1:size(tot12,1) ]
+        a22 = [ tot22[j][i] for j = 1:size(tot22,1) ]
+        res[i] = infoMatt(a11, a12, a22)
+        #res[i]=  infoMatt(tot11[i], [tot12[1][i],tot12[2][i],tot12[3][i]], [tot22[1][i],tot22[2][i],tot22[3][i]])
+    end
+    return res
+end
+function comp_mat(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr}}, tot13::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}}, tot33::Array{Array{juobs.Corr}}, tot23::Array{Array{juobs.Corr}})
+    mu = getfield.(tot11[1],:mu)
+    res = Array{infoMatt}(undef, length(mu))
+    for i = 1:length(mu)
+        a11 = [ tot11[j][i] for j = 1:size(tot11,1) ]
+        a12 = [ tot12[j][i] for j = 1:size(tot12,1) ]
+        a22 = [ tot22[j][i] for j = 1:size(tot22,1) ]
+        a33 = [ tot33[j][i] for j = 1:size(tot33,1) ]
+        a13 = [ tot13[j][i] for j = 1:size(tot13,1) ]
+        a23 = [ tot23[j][i] for j = 1:size(tot23,1) ]
+        res[i] = infoMatt(a11, a12, a13, a22, a23, a33)
+        #res[i]=  infoMatt(tot11[i], [tot12[1][i],tot12[2][i],tot12[3][i]], [tot22[1][i],tot22[2][i],tot22[3][i]])
+    end
+    return res
+end
+function comp_energy(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}}, evec::Bool=false; tnot::Int64=3)
+    aux_mat = comp_mat(tot11, tot12, tot22)
+    y0 = getfield(aux_mat[1], :y0)
+    res_en = Array{infoEn}(undef, length(aux_mat))
+    if !evec
+        for i = 1:length(aux_mat)
+            mat_obs = getfield(aux_mat[i], :mat)[y0+1:end-1]
+            evalues = uwgevp_tot(mat_obs, tnot)
+            res_en[i] = infoEn(energies(evalues), aux_mat[i])
+        end
+        return res_en
+    else
+        evec = Array{Array{}}(undef, length(aux_mat))
+        for i = 1:length(aux_mat)
+            mat_obs = getfield(aux_mat[i], :mat)[y0+1:end-1]
+            evalues, eigvec = uwgevp_tot(mat_obs, tnot, evec=true, iter=50)
+            res_en[i] = infoEn(energies(evalues), aux_mat[i])
+            evec[i] = eigvec
+        end
+        return res_en, evec
+    end
+end
+function comp_energy(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr}}, tot13::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}}, tot33::Array{Array{juobs.Corr}}, tot23::Array{Array{juobs.Corr}}, evec::Bool=false; tnot::Int64=2)
+    aux_mat = comp_mat(tot11, tot12, tot13, tot22, tot23, tot33)
+    y0 = getfield(aux_mat[1], :y0)
+    res_en = Array{infoEn}(undef, length(aux_mat))
+    if !evec
+        for i = 1:length(aux_mat)
+            mat_obs = getfield(aux_mat[i], :mat)[y0+1:end-1]
+            evalues = uwgevp_tot(mat_obs, tnot, iter=50)
+            res_en[i] = infoEn(energies(evalues), aux_mat[i])
+        end
+        return res_en
+    else
+        evec = Array{Array{}}(undef, length(aux_mat))
+        for i = 1:length(aux_mat)
+            mat_obs = getfield(aux_mat[i], :mat)[y0+1:end-1]
+            evalues, eigvec = uwgevp_tot(mat_obs, tnot, evec=true)
+            println("in comp_energy, i is = ", i)
+            res_en[i] = infoEn(energies(evalues), aux_mat[i])
+            evec[i] = eigvec
+        end
+        return res_en, evec
+    end
+end
+mutable struct infoEn
+    en_val::Array{Array{uwreal}}
+    mu::Vector{Float64}
+    y0::Int64
+    function infoEn(energ::Array{Array{uwreal}},matr::infoMatt )
+        en_val = energ
+        mu = getfield(matr,:mu)
+        y0 = getfield(matr,:y0)
+        new(en_val, mu, y0)
+    end
+end
+function pseudo_mat_elem(evec::Array{Array{T,2} where T,1}, mass::uwreal, mu::Array{Float64})
+    Rn = [exp(mass * (t)/2) * evec[t][1] for t =2:length(evec)]
+    aux = sqrt(2)*sum(mu)  / mass^1.5 
+    return uwdot(aux,Rn)
+end
+function vec_mat_elem(evec::Array{Array{T,2} where T,1}, mass::uwreal)
+    Rn = [exp(mass * t/2) * evec[t][1] for t =1:length(evec)]
+    return uwdot(Rn , 1/ mass)
+end
+function extract_dec_const(mat_elem::uwreal, mass::uwreal, mu::Array{Float64})
+    return sqrt(2)*sum(mu)*mat_elem / mass^1.5
+end
+function match_muc(muh, m_lh, m_lh_star, target)
+    M = (m_lh .+ 3 .* m_lh_star) ./ 4
+
+    par, chi2exp = lin_fit(muh, M)
+    muh_target = x_lin_fit(par, target)
+    return muh_target
+end
+function match_muc(muh, m_lh, m_sh, m_lh_star, m_sh_star, target)
+    M = (2 .* m_lh .+ 6 .* m_lh_star .+ m_sh .+ 3 .* m_sh_star) ./ 12
+
+    par, chi2exp = lin_fit(muh, M)
+    muh_target = x_lin_fit(par, target)
+    return muh_target
+end
\ No newline at end of file

src/juobs.jl

 module juobs
-using ADerrors, PyPlot, Statistics, LaTeXStrings, LinearAlgebra
+using ADerrors, PyPlot, Statistics, LaTeXStrings, LinearAlgebra, LsqFit
 
 include("juobs_types.jl")
 include("juobs_linalg.jl")
@@ -7,9 +7,9 @@ include("juobs_reader.jl")
 include("juobs_tools.jl")
 include("juobs_obs.jl")
 
-export read_mesons, read_ms1
-export get_matrix, uwgevp_tot, energies
-export corr_obs, plat_av, lin_fit, x_lin_fit, y_lin_fit
-export meff, dec_const_pcvc 
+export read_mesons, read_ms1, read_ms, read_md
+export get_matrix, uwgevp_tot, energies, uwdot
+export corr_obs, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine
+export meff, dec_const_pcvc, comp_t0
 
 end # module

src/juobs_linalg.jl

@@ -85,16 +85,16 @@ Given a vector where each entry evals[t] is a uwreal array of eigenvalues, this
 The index t here runs from 1:T=lenght(evals), while the index i stands for the number of energy levels computed: i = length(evals[t])
 It returns a vector array eff_en where each entry  eff_en[t] contains the first N states energies as uwreal objects 
 """
-function energies(evals::Vector{Vector{uwreal}})
+function energies(evals::Union{Vector{Vector{uwreal}},Array{Array{uwreal}} })
     time = length(evals) 
     n = length(evals[1])
-    eff_en = Vector{Vector}(undef, n)
+    eff_en = Array{Array{uwreal}}(undef, n)
     aux_en = Vector{uwreal}(undef, time-1)
     for i in 1:n
         #aux_en = Vector{uwreal}(undef, time-1)
         for t in 1:time-1
             ratio = evals[t][i] / evals[t+1][i]
-            aux_en[t] = log(sqrt(ratio * ratio))
+            aux_en[t] = 0.5*log(ratio * ratio)
         end
         uwerr.(aux_en)
         eff_en[i] = copy(aux_en)
@@ -117,21 +117,35 @@ The columns of such matrix are the eigenvectors associated with eigenvalues eval
 The keyword iter, set by default to 30, selects the number of iterations of the qr algorithm before stopping.
 """
 function uwgevp_tot(mat_list::Vector{Matrix}, tnot::Int64; iter::Int64 = 30, evec::Bool = false)
+    n = length(mat_list)
+    evals = Array{Array{uwreal}}(undef, 0 )
     if !evec
-        aux_evals = uwgevp.(mat_list, c_tnot=mat_list[tnot], iter = iter)
-        evals = get_diag.(aux_evals)
+        println(tnot)
+        c_inv = invert(mat_list[tnot])
+        for i =1:n
+            aux_evals = uwgevp(mat_list[i], c_inv, iter = iter)
+            push!(evals , get_diag(aux_evals))
+        end
         return evals
     else
-        aux_evals, aux_evec =  uwgevp.(mat_list, c_tnot=mat_list[tnot], iter = iter, evec= true)
-        evals = get_diag.(aux_evals)
-        evecs = aux_evec
+        evecs = Array{Matrix}(undef, 0)
+        for i = 2:n-4
+            c_inv = invert(mat_list[i])
+            aux_evals, aux_evec =  uwgevp(mat_list[i+3], c_inv,  evec=evec, iter=iter)
+            push!(evals,  get_diag(aux_evals))
+            ptilda = norm_evec(aux_evec, mat_list[i])
+            m = uwdot(mat_list[i], ptilda)
+            push!(evecs, m)
+            println( "\n t is ", i+3, " tnot is ", i)            
+            #if  2*tnot < i && i<n-5
+            #    tnot+=5
+            #end
+        end
         return evals, evecs
     end
 end
-
-
 """
-uwgevp(c_t:: Matrix{uwreal}, c_tnot::Matrix{uwreal}; iter = 30, evec = false)
+uwgevp(c_t:: Matrix{uwreal}, c_tnot_inv::Matrix{uwreal}; iter = 30, evec = false)
 
 Given two matrices of uwreal C(t) and C(t0), this method solves the generalized eigenvalue problem
 C(t)V = C(t_0)D V, where V is the eigenvector matrix and D is the diagonal matrix of eigenvalues.
@@ -139,8 +153,7 @@ C(t)V = C(t_0)D V, where V is the eigenvector matrix and D is the diagonal matri
 If evec = true also the eigenvectors are returned as columns of the orthogonal matrix u.
 The keyword iter, set by default to 30, selects the number of iterations of the qr algorithm before stopping.
 """
-function uwgevp(c_t:: Matrix{uwreal}; c_tnot::Matrix{uwreal}, iter::Int64 = 30, evec::Bool = false)
-    c_tnot_inv = invert(c_tnot)
+function uwgevp(c_t:: Matrix{uwreal}, c_tnot_inv::Matrix{uwreal}; iter::Int64 = 30, evec::Bool = false)
     c = uwdot(c_tnot_inv, c_t) #matrix to diagonalize
     return uwevp(c, iter = iter, evec = evec)
 end
@@ -171,7 +184,18 @@ function uwevp(a::Matrix{uwreal}; iter = 30, evec = false)
     end
 end
 
-
+function norm_evec(evec::Matrix{uwreal}, ctnot::Matrix{uwreal})
+    n=size(evec,1)
+    res_evec=Matrix{uwreal}(undef, n, n)
+    for i =1:n
+        aux_norm = (uwdot(evec[:,i], uwdot(ctnot, evec[:,i])).^2).^0.25
+        res_evec[:,i] = 1 ./aux_norm .* evec[:,i]
+    end
+    #println("res_evec is ", res_evec)
+    #println("evec is ", evec)
+    #println(res_evec)
+    return res_evec
+end
 function get_diag(a::Matrix{uwreal})
     n = size(a,1)
     res = [a[k, k] for k = 1:n]
@@ -354,6 +378,10 @@ of uwreal data type.
 
 """
 uwdot(a::Vector{uwreal}, b::Vector{uwreal}) = sum(a .* b)
+uwdot(a::Vector{uwreal}, b::uwreal) = [a[i] * b for i=1:length(a)]
+uwdot(a::uwreal, b::Vector{uwreal}) = uwdot(b,a) 
+uwdot(a::Matrix{uwreal}, b::Vector{uwreal}) = uwdot(a, reshape(b,(length(b),1)))
+uwdot(a::Vector{uwreal}, b::Matrix{uwreal}) = uwdot(reshape(a,(1,length(a))), b)
 
 function uwdot(a::Matrix{uwreal}, b::Matrix{uwreal})
     n,m = size(a)

src/juobs_obs.jl

@@ -98,3 +98,135 @@ function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu
 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)
+
+#t0
+@doc raw"""
+
+comp_t0(y::YData, plat::Vector{Int64}; pl::Bool=false, L::Int64=1)
+
+comp_t0(Y::Vector{YData}, plat::Vector{Int64}; pl::Bool=false, L::Int64=1)
+
+Computes t0 using the energy density of the action Ysl(Yang-Mills action).
+t0 is computed in the plateau plat.
+A (linear) interpolation in t is performed to find t0.
+
+The flag pl allows to show the plot.
+
+```@example
+#Single replica
+Y = read_ms(path)
+rw = read_ms(path_rw)
+
+t0 = comp_t0(Y, [38, 58], L=32)
+t0_r = comp_t0(Y, [38, 58], L=32, rw=rw)
+
+#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}; pl::Bool=false, L::Int64=1, rw::Union{Matrix{Float64}, Nothing}=nothing)
+    Ysl = Y.Ysl
+    t = Y.t
+    id = Y.id
+
+    Ysl = isnothing(rw) ? Ysl : apply_rw(Ysl, rw)
+    
+    xmax = size(mean(Ysl, dims=3))[1]
+    nt0 = t0_guess(t, Ysl, plat, L)
+    
+    Y = Matrix{uwreal}(undef, xmax, 3)
+    for i = 1:xmax
+        k = 1
+        for j = nt0-1:nt0+1
+            Y[i, k] = uwreal(Ysl[i, j, :], id)
+            k = k + 1
+        end
+    end
+    x = t[nt0-1:nt0+1]
+    t2E = [plat_av(Y[:, j], plat) for j=1:3] .* x.^2 / L^3
+    
+    par, chi2exp = lin_fit(x, t2E)
+
+    t0 = x_lin_fit(par, 0.3)
+    if pl
+        uwerr(t0)
+        uwerr.(t2E)
+        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())        
+        end
+    return t0
+end
+
+function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; pl::Bool=false, L::Int64=1, rw::Union{Vector{Matrix{Float64}}, Nothing}=nothing)
+    nr = length(Y)
+    Ysl = getfield.(Y, :Ysl)
+    t = getfield.(Y, :t)
+    id = getfield.(Y, :id)
+    vtr = getfield.(Y, :vtr)
+    replica = length.(vtr)
+    
+
+    if !all(id .== id[1])
+        println("IDs are not equal")
+        return nothing
+    end
+
+    Ysl = isnothing(rw) ? Ysl : apply_rw.(Ysl, rw)
+
+    xmax = size(mean(Ysl[1], dims=3))[1]
+
+    tmp = Ysl[1]
+    [tmp = cat(tmp, Ysl[k], dims=3) for k = 2:nr]
+    nt0 = t0_guess(t[1], tmp, plat, L)
+
+    Y = Matrix{uwreal}(undef, xmax, 3)
+    for i = 1:xmax
+        k = 1
+        for j = nt0-1:nt0+1
+            Y[i, k] = uwreal(tmp[i, j, :], id[1], replica)
+            k = k + 1
+        end
+    end
+    x = t[1][nt0-1:nt0+1]
+    t2E = [plat_av(Y[:, j], plat) for j=1:3] .* x.^2 / L^3
+    println(t2E)
+    par, chi2exp = lin_fit(x, t2E)
+
+    t0 = x_lin_fit(par, 0.3)
+    if pl
+        uwerr(t0)
+        uwerr.(t2E)
+        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())        
+    end
+    return t0
+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=1), dims=3)[1, :, 1] / L^3
+    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
\ No newline at end of file

src/juobs_reader.jl

@@ -156,21 +156,7 @@ function read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union
     return res
 end
 
-@doc raw"""
-read_ms1(path::String; v::String="1.2")
-
-Reads openQCD ms1 dat files at a given path. This method returns a matrix W[irw, icfg] that contains the reweighting factors, where
-irw is the rwf index and icfg the configuration number.
-The function is compatible with the output files of openQCD v=1.2, 1.4 and 1.6. Version can be specified as argument.
-
-Examples:
-```@example    
-read_ms1(path)
-read_ms1(path, v="1.4")
-read_ms1(path, v="1.6")
-```
-"""
-function read_ms1(path::String; v::String="1.2")
+function read_rw(path::String; v::String="1.2")
     data = open(path, "r")
     nrw = read(data, Int32)
 
@@ -208,9 +194,109 @@ function read_ms1(path::String; v::String="1.2")
             end
         end
     end
-    
+    close(data)
+    return r_data
+end
+
+@doc raw"""
+read_ms1(path::String; v::String="1.2")
+
+Reads openQCD ms1 dat files at a given path. This method returns a matrix W[irw, icfg] that contains the reweighting factors, where
+irw is the rwf index and icfg the configuration number.
+The function is compatible with the output files of openQCD v=1.2, 1.4 and 1.6. Version can be specified as argument.
+
+Examples:
+```@example    
+read_ms1(path)
+read_ms1(path, v="1.4")
+read_ms1(path, v="1.6")
+```
+"""
+function read_ms1(path::String; v::String="1.2")
+    r_data = read_rw(path, v=v)
+    nrw = length(r_data)
+    ncnfg = size(r_data[1])[3]
     W = zeros(Float64, nrw, ncnfg)
     [W[k, :] = prod(mean(exp.(.-r_data[k]), dims=2), dims=1) for k = 1:nrw]
-    close(data)
     return W
-end
\ No newline at end of file
+end
+@doc raw"""
+read_md(path::String)
+
+Reads openQCD  pbp.dat files at a given path. This method returns a matrix md[irw, icfg] that contains the derivatives dS/dm, where
+md[irw=1] = dS/dm_l and md[irw=2] = dS/dm_s
+
+Examples:
+```@example    
+md = read_md(path)
+```
+"""
+function read_md(path::String)
+    r_data = read_rw(path, v="1.4")
+    nrw = length(r_data)
+    ncnfg = size(r_data[1])[3]
+    md = zeros(Float64, nrw, ncnfg)
+    [md[k, :] = prod(mean(r_data[k], dims=2), dims=1) for k = 1:nrw]
+    return md
+end
+
+@doc raw"""
+read_ms(path::String)
+
+Reads openQCD ms dat files at a given path. This method return: 
+    t(t): flow time values
+    Wsl(x0, t, icfg): the time-slice sums of the densities of the Wilson plaquette action
+    Ysl(x0, t, icfg): the time-slice sums of the densities of the Yang-Mills action 
+    Qsl(x0, t, icfg): the time-slice sums of the densities of the topological charge
+
+Examples:
+```@example    
+t, W, Y, Q = read_ms(path)
+```
+"""
+function read_ms(path::String; id::Union{Int64, Nothing}=nothing)  
+    if isnothing(id)
+        bname = basename(path)
+        m = findfirst(r"[A-Z][0-9]{3}r[0-9]{3}", bname)
+        id = parse(Int64, bname[m[2:4]])
+    end  
+    data = open(path, "r")
+
+    dn = read(data, Int32)
+    nn = read(data, Int32)
+    tvals = read(data, Int32)
+    eps = read(data, Float64)
+
+    fsize = filesize(path)
+    datsize=4 + 3*8*(nn + 1) * tvals # measurement size of each cnfg
+
+    ntr = Int32((fsize - 3*4 - 8) / datsize)
+
+    vntr = Vector{Int32}(undef, ntr)
+    
+    # x0, t, cfg 
+    Wsl = Array{Float64}(undef, tvals, nn + 1, ntr) 
+    Ysl = Array{Float64}(undef, tvals, nn + 1, ntr)
+    Qsl = Array{Float64}(undef, tvals, nn + 1, ntr)
+
+    for itr = 1:ntr
+        vntr[itr] = read(data, Int32)
+        for iobs = 1:3
+            for inn = 0:nn
+                tmp = Vector{Float64}(undef, tvals)
+                read!(data, tmp)
+                if iobs == 1
+                    Wsl[:, inn + 1, itr] = tmp
+                elseif iobs == 2
+                    Ysl[:, inn + 1, itr] = tmp
+                elseif iobs == 3
+                    Qsl[:, inn + 1, itr] = tmp
+                end
+
+            end
+        end
+    end
+    close(data)
+    t = Float64.(0:nn) .* dn .* eps
+    return YData(vntr, t, Ysl, id)
+end

src/juobs_tools.jl

@@ -95,7 +95,7 @@ function plat_av(obs::Vector{uwreal}, plat::Vector{Int64})
     return av 
 end
 
-function lin_fit(x::Vector{Float64}, v::Vector{Float64}, e::Vector{Float64})
+function lin_fit(x::Vector{<:Real}, v::Vector{Float64}, e::Vector{Float64})
     sig2 = e .* e
     S = sum(1 ./ sig2)
     Sx = sum(x ./ sig2)
@@ -108,7 +108,7 @@ function lin_fit(x::Vector{Float64}, v::Vector{Float64}, e::Vector{Float64})
     return par
 end
 @doc raw"""
-lin_fit(x::Vector{Float64}, y::Vector{uwreal})
+lin_fit(x::Vector{<:Real}, y::Vector{uwreal})
 
 Computes a linear fit of uwreal data points y. This method return uwreal fit parameters and chisqexpected.
 
@@ -117,7 +117,7 @@ fitp, csqexp = lin_fit(phi2, m2)
 m2_phys = fitp[1] + fitp[2] * phi2_phys
 ```
 """
-function lin_fit(x::Vector{Float64}, y::Vector{uwreal})
+function lin_fit(x::Vector{<:Real}, y::Vector{uwreal})
     uwerr.(y)
     par = lin_fit(x, value.(y), err.(y))
 
@@ -140,6 +140,34 @@ Computes the results of a linear interpolation/extrapolation in the y axis
 """
 y_lin_fit(par::Vector{uwreal}, x::Union{uwreal, Float64}) = par[1] + par[2] * x
 
+@doc raw"""
+fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}, param::Int64=3)
+
+Given a model function with a number param of parameters and an array of uwreal,
+this function fit ydata with the given model and print  fit information
+The method return an array upar with the best fit parameters with their errors.
+'''@example
+@. model(x,p) = p[1] + p[2] * exp(-(p[3]-p[1])*x)
+fit_routine(model, ydata, param=3)
+"""
+function fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}, param::Int64=3; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing )
+    yval = value.(ydata)
+    yer = err.(ydata)
+    chisq = gen_chisq(model, xdata, yer)
+    fit = curve_fit(model, xdata, yval, 1.0 ./ yer.^2, fill(0.5, param))
+    (upar, chi_exp) = isnothing(wpm) ? fit_error(chisq, coef(fit), ydata) : fit_error(chisq, coef(fit), ydata, wpm)
+    for i = 1:length(upar)
+        uwerr(upar[i])
+        print("\n Fit parameter: ", i, ": ")
+        details(upar[i])
+    end
+    println("Chisq / chiexp: ", chisq(coef(fit), ydata), " / ", chi_exp, " (dof: ", dof(fit),")")
+    return upar
+end
+function gen_chisq(f::Function, x::Array{<:Real}, err::Vector{Float64})
+    chisq(par, dat) = sum((dat .- f(x,par)).^2 ./err.^2)
+    return chisq
+end
 #TODO: add combined fits
 #=
 using LsqFit

src/juobs_types.jl

@@ -138,6 +138,14 @@ mutable struct Corr
     end 
 end
 
+mutable struct YData
+    vtr::Vector{Int32}
+    t::Vector{Float64}
+    Ysl::Array{Float64, 3}
+    id::Int64
+    YData(a, b, c, d) = new(a, b, c, d)
+end
+
 function Base.show(io::IO, a::GHeader)
     print(io, "ncorr = ", a.ncorr, "\t")
     print(io, "nnoise = ", a.nnoise, "\t")