analysis + results updated

analysis/analysis.jl

-#TODO include rw, different plateaux depending on obs, print chi2, scaling+chiral plots
+#TODO include rw, different plateaux depending on obs, print chi2, compute t0, compute mpi
 using juobs, ADerrors, DelimitedFiles, PyPlot, LaTeXStrings
 const path = "/home/javier/Lattice/charm/production_2"
 const path_plat = "/home/javier/Lattice/juobs/analysis/plat.txt"
@@ -11,6 +11,7 @@ const beta = [3.46 , 3.55, 3.55, 3.70, 3.70]
 const R = ["H400r001", ["N200r000", "N200r001"], ["N203r000", "N203r001"], "N300r002", "J303r003"]
 include("/home/javier/Lattice/juobs/constants/juobs_const.jl")
 include("/home/javier/Lattice/juobs/analysis/functions.jl")
+const phi2 = 8 .* t0.(beta) .* [0.16345, 0.09222, 0.11224, 0.10630, 0.06514].^2 #8t0 m_pi^2
 
 
 m_lh = Vector{Vector{uwreal}}(undef, length(ensembles))
@@ -211,3 +212,52 @@ for iens=1:length(ensembles)
     println("(", ensembles[iens], ")", L"$f_D \sqrt{8t_0}= $", f_D[iens])
     println("(", ensembles[iens], ")", L"$f_Ds \sqrt{8t_0}= $", f_Ds[iens])
 end
+
+###########################
+###### FITS ###############
+###########################
+x = [1 ./(8 .* t0.(beta)) phi2] #x1 = a^2 / (8t0), x2 = 8t0 mpi^2
+uwerr.(x)
+phi2_ph = (t0_ph[1] * 139.57039 / hc)^2
+uwerr(phi2_ph)
+#f(a2/8t0,phi2) = p[1]+ p[2](a2/8t0) + (p[3] + p[4](a2/8t0)) * phi2 
+@. model(x, p) = (p[1] + p[2] * x[:, 1]) + (p[3] + p[4] * x[:, 1]) * x[:, 2] #linear fits
+obs = [muc, m_D, m_Ds, m_D_star, m_Ds_star, f_D, f_Ds] #sqrt(8t0) obs
+ttl_obs = ["muc", "m_D", "m_Ds", "m_D_star", "m_Ds_star", "f_D", "f_Ds"]
+ylbl = [L"$Z^{tm}_M \mu_c \sqrt{8t_0}$", L"$M_D \sqrt{8t_0}$", L"$M_{D_s} \sqrt{8t_0}$",
+L"$M_{D^*} \sqrt{8t_0}$", L"$M_{D^*_s} \sqrt{8t_0}$", L"$f_D \sqrt{8t_0}$", L"$f_{D_s} \sqrt{8t_0}$"]
+xlbl = L"$\phi_2$"
+obs_t0 = Vector{uwreal}(undef, length(obs)) #sqrt(8t0)obs @ CL & phi2_phys
+for k = 1:length(obs)
+    println("OBS ",ttl_obs[k])
+    par = fit_routine(model, value.(x), obs[k], 4)
+    obs_t0[k] = par[1] + par[3] * phi2_ph
+    uwerr(obs_t0[k])
+    figure()
+    for b in unique(beta) #select point with same beta
+        nn = findall(x-> x == b, beta)
+        lgnd = string(L"$\beta = $", b)
+        errorbar(value.(x[nn,2]), value.(obs[k][nn]), err.(obs[k][nn]), err.(x[nn,2]), fmt="x", label=lgnd)
+    end
+    lgnd=L"$\mathrm{CL}$"
+    errorbar(value(phi2_ph), value(obs_t0[k]), err(obs_t0[k]), err(phi2_ph), fmt="x", zorder=2, label=lgnd)
+    axvline(value(phi2_ph), ls="--", color="black", zorder=1, lw=0.6, label="")
+    xlabel(xlbl)
+    ylabel(ylbl[k])
+
+    legend()
+    display(gcf())
+    t = string("fit_", ttl_obs[k], ".pdf")
+    savefig(joinpath(path_plot, t))
+    close()
+
+end
+
+obs_ph = Vector{uwreal}(undef, length(obs))
+for k = 1:length(obs)
+    println(ylbl[k], " = ", obs_t0[k])
+    #phys
+    obs_ph[k] = obs_t0[k] * hc / t0_ph[1]
+    uwerr(obs_ph[k])
+    println(ttl_obs[k], "(MeV) = ", obs_ph[k])
+end

analysis/results.txt

@@ -34,3 +34,17 @@
 (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