old test for matrix elements

src/testing_mat_elem.jl

+using Revise, juobs, BDIO, ADerrors, DelimitedFiles, PyPlot, LaTeXStrings, LsqFit
+
+
+#============= SET UP VARIABLES ===========#
+const path_data = "/Users/ale/Desktop/data"
+const path_plat = "/Users/ale/automation/plat.txt"
+const path_results = "/Users/ale/Desktop/results_gevp"
+const ensembles = [ "J303", "N300"]
+const sector = Dict("ll"=>false, "ls"=>false, "lh"=>true, "ss"=>false, "sh"=>true, "hh"=>true )
+const tau = 4
+const _t0 = 2
+const range_t_fit = [2,5]
+const rwf = false
+const compute_t0 = false
+const mass_shift = false #not implemented yet
+#go to line 443 to change CL+chiral extrapolation fit model
+
+#========= DO NOT CHANGE AFTER THIS LINE ==========#
+include("types.jl")
+include("tools.jl")
+include("const.jl")
+
+#======== GET ENSEMBLE INFORMATION FROM DATABASE ==========#
+ensinfo = Vector{EnsInfo}(undef, length(ensembles))
+for i in 1:length(ensembles)
+    ens = ensembles[i]
+    try
+        ensinfo[i]= EnsInfo(ens, ens_db[ens])
+    catch
+        error("The ensemble id ", ens, " was not found in the const.jl ens_db database.
+        Please check the ensemble id or update the database")
+    end
+end
+gen_mulist = Vector{Vector{Vector{Float64}}}(undef, length(ensembles))
+println("\n Building matrices \n")
+@time begin
+    matinfo = Vector{Vector{MatInfo}}(undef,length(ensembles))
+    matinfo_vec = Vector{Vector{MatInfo}}(undef,length(ensembles))
+    for i in 1:length(ensinfo)
+        @time begin
+            ens = ensinfo[i]
+            println("\n Ensemble: ", ens.id)
+            pp = get_corr(ens,"G5", "G5", rw=rwf)
+            gen_mulist[i] = getfield.(pp, :mu)
+            a1a1 = get_corr(ens, "G1G5", "G1G5", rw=rwf)
+            matinfo[i] = comp_mat(ens, pp, tau)
+            matinfo_vec[i] = comp_mat(ens, a1a1, tau)
+            println("Time:")
+        end
+    end
+    println("Total time:")
+end
+
+## TESTING DECAYS
+
+## compute ps mass 
+mat = matinfo[1][6].mat_list[96:end-1]
+mu = matinfo[1][6].mu
+
+evals = getall_eigvals(mat, 2)
+en_eff = energies(evals)
+psmass = plat_av(en_eff[1], [25,70])
+
+## extract eigenvectors 
+function getall_eig(a::Vector{Matrix{uwreal}}, t0; iter=30)
+    n = length(a)
+    res = Vector{Matrix{uwreal}}(undef, n)
+    [res[i] = uweigvecs(a[i], a[t0]) for i=1:n]
+    return res 
+end
+##
+evec = getall_eig(mat, 2)
+
+##
+function rn(evec::Vector{Matrix{uwreal}}, ct_mat::Vector{Matrix{uwreal}}, mass::uwreal,  n::Int64)
+    t = length(evec)
+    res = Vector{uwreal}(undef, t)
+    for i in 1:t
+        aux = uwdot(evec[i][:,n], uwdot(ct_mat[i], evec[i][:,n]))[1]
+        aux2 = 1 / (aux)^(0.5)   * exp(mass * (i-2)/2)
+        res[i] = aux2 * uwdot(evec[i][:,n], ct_mat[i][:,n])
+    end
+    return res
+end
+
+r1 = rn(evec, mat, psmass, 1)
+uwerr.(r1)
+
+fig = figure()
+    x = collect(1:length(r1))
+    errorbar(x, value.(r1), yerr=err.(r1), fmt="*" )
+    ylim(-0.2,-0.15)
+display(fig)
+
+##
+function dec(mat_elem::Vector{uwreal}, plat::Vector{Int64}, mass::uwreal, mu::Vector{Float64})
+    return sqrt(2/ (mass)^3) * sum(mu)* aux 
+end
+
+fps = - dec(r1, [25,70], psmass, mu )
+uwerr(fps)
+##
+path = "/Users/ale/Desktop/data/N300/N300r002_k0.1372069_mu0.004119_mu0.13775_mu0.15225_mu0.145.mesons.dat"
+dat = read_mesons(path, "G5", "G5")
+
+pp_corr  = corr_obs.(dat, L=48)
+dec_const_pcvc(pp_corr[6], [60,110], psmass)
\ No newline at end of file