Example of fitting process

This page will document the fitting process by taking the example of the high-statistics (time-smeared 0,16,32,48) run. The correlators are taken from /home/simone/qcd/onia/l2064f21b676m010m050-Vprob/Ds_k0.119_csw1.72/fits-avg-t0-t16-t32-t48/corrs. Interpreted, this is the 20^3x64 lattice with 2+1 flavors at beta=6.76 with m_u,d=.010 and m_s=.050. The results had the gamma-matrix bug at non-zero momentum so some quantities may not be trustworthy. The correlators were from the kappa=0.119 run with c_sw=1.72. The data is averaged from sources at t=0, 16, 32, and 48 for better statistics.

This document is organized as follows:

Notational notes
Initial selection of priors
Confirmation of priors
Other tests
The run
Summary of results

Notational notes

Jim generates data with filenames that look like, for example pi_d_2S_0.119_0.119_p000.bz2. Interpreted, this means this is a pi correlator, with local (delta) and 2S smearings. Both quarks are heavy quarks with kappa=0.119. It's the zero-momentum state (p000, and the data is compressed with bzip2.

My program generates a lot of plots, and therefore it's necessary for auto-generated names for those plots. I've come up with a naming scheme that encodes many of the parameters while keeping the name acceptably short. It typically follows the template: testy100_14p_p000_5_0,2_7_0,0.ro.8.boots. The name is organized as follows:

name (in this case, "test")
whether we're using Bayes statistics ("y" or "n")
number of bootstraps used (when bootstraps are done, in this case 100)
number of the dataset (14 refers to the run)
momentum (in this case 000)
model function (5 refers to a sum of coshes with coefficients constrained to be positive)
minimum correlator in fit (0 = local-local)
maximum correlator in fit (2 = 2S-2S)
SVD cut power for data (correlator) matrix (in this case, 1e-7)
tmin,tmax (0,0 means to use the default value as specified when selecting priors)
propagator (in this case, ro, but can also be pi, b1, a0, a1, hlpi, hlro)
number of energy levels in bootstrap result (in this case, 8)
suffix indicating it's a python bootstrap file (other suffixes include .ax and .jimboot)

Initial selection of priors

First we have to choose what priors (and which timeslices) to use. A helpful guide is from the effective mass plot. I check them using my eff_mass.py program. The resulting plots are accessible from the following links:

pi

ro

Summary

I choose as my fit ranges the values shown for each dataset. The energy priors are set as:
a0: E0=2.0+-.2, E1=2.5+-.2, Coefficients=exp{0+-2}
a1: E0=2.0+-.2, E1=2.5+-.2, Coefficients=exp{0+-2}
b1: E0=2.0+-.1, E1=2.5+-.2, Coefficients=exp{0+-2}
pi: E0=1.65+-.05, E1=2.15+-.2, Coefficients=exp{0+-2}, momentum-scale=.02
ro: E0=1.71+-.05, E1=2.21+-.2, Coefficients=exp{0+-2}, momentum-scale=.02

The ground state is meaningful as it is. The excited state is actually taken as a splitting, so in each case here it's internally taken to be 0.5 +- .2. Additionally, the model function actually wants the splitting as exp(E) so in reality it's not varying from .3 to .7, but rather from something like .34 to .74. Higher energy levels are calculated based on this splitting, with the width taken from the E1 width. Similarly, the coefficients have priors that vary on a logarithmic scale. What this means is that I specify 0+-2, and they actually have the central value exp(0)=1, and vary from exp(-2)=.14 to exp(2)=7.4. Finally, there is a momentum scale. This offsets the energy spectrum a little for each bit of momentum in the state. So the p000 data has no offset, p100 has an offset of momentum-scale, and p221 and p300 each have an offset of 9 times the momentum scale.

Confirmation of priors

Now we're ready for some fits. First, it's important to confirm that our priors are reasonable. It is sufficient to check our priors using only the local-local zero-momentum data. Again, it's time for some more plots. In the following, the number of states in the fit is on the x-axis. Additionally, the priors were varied (by factors of 2), and plotted in color. The red point is our default prior, as specified above. The black point to its left reduces the prior widths by half, while the green point to its right increases the prior widths by a factor of two. In all my plots, I vary the colors using the ordering red, green, blue, black, red... where the central value is always red.

a0

a0.E0-full.gif
a0.E0-last.gif
a0.E1-full.gif
a0.E1-last.gif

a1

a1.E0-full.gif
a1.E0-last.gif
a1.E1-full.gif
a1.E1-last.gif

b1

b1.E0-full.gif
b1.E0-last.gif
b1.E1-full.gif
b1.E1-last.gif

pi

pi.E0-full.gif
pi.E0-last.gif
pi.E1-full.gif
pi.E1-last.gif

ro

ro.E0-full.gif
ro.E0-last.gif
ro.E1-full.gif
ro.E1-last.gif

Summary

We see from the "full" plots that the results of interest are stable as we add energy levels. We see from the "last" plots that our choice of priors doesn't affect the results for the quantities of interest.

Other tests

We could investigate other things here, such as varying tmin/tmax (I selected values from the effective mass plots, as shown above), prior values (also selected from the effective mass plots), which datasets to include, and SVD cuts (1e-7 for both the data correlation matrix and for alpha-matrix inversion (in Marquardt algorithm)), but due to time constraints, I'm skipping this step. I'll come back to it later if I get a chance.

The complete run

Now we need to do a more complete run, including smeared correlators, non-zero momentum states, and bootstrap errorbars. Because of the wide variety of parameters, and due to the instability of qcdhome (hardware issues), it's a good idea to do this as several separate runs. The final results are given as plots and bootstrap files. Here are the plots:

a0

a0y100_14p_p000_5_0,0_7_0,0_0.a0.gif

a1

a1y100_14p_p000_5_0,0_7_0,0_0.a1.gif

b1

b1y100_14p_p000_5_0,0_7_0,0_0.b1.gif
b1y100_14p_p000_5_0,1_7_0,0_0.b1.gif
b1y100_14p_p000_5_0,2_7_0,0_0.b1.gif

pi

ro

Summary

The plots all indicate stability as we add additional energy levels to the fit. Also the bootstrap errorbars (the error on the n_states=8 point) are compatible with the chi^2 errors. The kinpi and kinro datasets can be used to find the M2 for the pi and ro states. Doing so gives (with confidence level chi^2/dof = 3.9/6 q=0.688):
a[ 0] = [2.8553 +0.0014 -0.0010]e+00 a[ 1] = [0.7493 +0.0073 -0.0036]e+00 a[ 2] = -[0.5446 +0.1192 -0.1090]e-01 a[ 3] = -[0.3197 +0.1048 -0.0727]e-01 m2 = [2.2550 +0.0104 -0.0214]e+00 Here are links to the dispersion relation bootstraps and a plot.

Summary of results

Finally we can see what the end results are. I feed the bootstraps into an analysis program, and find the following:

Files

a0: a0y100_14p_p000_5_0,0_7_0,0.a0.8.boots
a1: a1y100_14p_p000_5_0,0_7_0,0.a1.8.boots
b1: b1y100_14p_p000_5_0,2_7_0,0.b1.8.boots
pi: piy100_14p_p000_5_0,2_7_0,0.pi.8.boots
ro: roy100_14p_p000_5_0,2_7_0,0.ro.8.boots
kin1S: kin1Sy100_14p_5_0,0_7_0,0.1S.8.boots
Ds:
kinDs:

Quarkonia

M1
pi(1S): 1.642329 +0.000434 -0.000222
pi(2S): 2.043957 +0.026990 -0.019248
ro(1S): 1.705695 +0.000733 -0.000450
ro(2S): 2.102872 +0.017351 -0.019250
b1(1P): 1.978191 +0.005351 -0.007050
a0(1P): 1.942058 +0.008834 -0.012713
a1(1P): 1.999965 +0.010456 -0.016519
M2
kin1S: 2.254990 +0.010592 -0.021433
Splittings
(1P-1S) splitting: 0.28834 + 0.00494 - 0.00746 (a^(-1) = 1.59037 + 0.04227 - 0.02677 GeV)
(2S-1S) splitting: 0.39829 + 0.01982 - 0.01711 (a^(-1) = 1.49487 + 0.06709 - 0.07087 GeV)
hyperfine splitting: 0.10078 + 0.00308 - 0.00147 GeV
aM(hyperfine split): 0.06337 + 0.00049 - 0.00037
2S hyperfine splitting: 0.09370 + 0.02511 - 0.03451 GeV
aM(2S hyperfine split): 0.06337 + 0.00049 - 0.00037
R = ((2S-1S)/(1P-1S))_lat / ((2S-1S)/(1P-1S))_exp = a^{-1}_{1P-1S} / a^{-1}_{2S-1S}
R keeping correlations: 1.06388 + 0.06053 - 0.04428
R assuming no correlation: 1.06388 +- 0.05425
R_hyp = ((psi' - eta_c')/(psi-eta_c)_lat / ((psi'-eta_c')/(psi-eta_c))_exp R_hyp keeping correlations: 1.18464 + 0.28317 - 0.45451
R_hyp assuming no correlation: 1.18464 +- 0.37789

Heavy-light

M1
M2

Combined quantities

M2(onia1S)/M1(onia1P-1S) = 7.82066 + 0.22806 - 0.16502