Unanswered questions

• One way to prevent degenerate energy levels would be to introduce a mandatory minimum spacing. Currently we force an ordering by specifying delta_E = eps^2. We could introduce a minimum spacing as delta_E = minimum + eps^2. Any thoughts on whether this idea is good or bad?
• Is there a "right" way to decide which fit parameters to select? Obviously we want something that is stable under changes in priors, number of energy levels, etc. But that's hard to quantify.
• For a spectrum plot, we use the eta_c/Jpsi splitting to get a zero offset, and then h_c to set the scale. But these parameters have errors themselves. So the zero offset and scale are unknown to some degree. Do we need to reflect this fact in error bars elsewhere? How??
• I'm assuming it's reasonable to expect the 1S-1S correlator will be unable to see the 2S state. What does this imply about seeing the h_c? (The b1 correlator is the 1P state, right?)
• I need to start working on my poster for Lattice'02. Any suggestions on format (or a program to use to generate it) would be greatly appreciated. On a related note, is there any way to make axis generate thicker lines?
• Bugra and I would like to come up to Fermilab on either Monday or Tuesday. Which day works best for everyone?
• Would generating only the diagonal p000 states be faster? We don't have much time for debugging our correlator data, and the other states aren't needed for the fits anyway.

Log of charmonium spectrum fit results

In all of the following, colors are done so that:
   Red = 1 correlator only (delta-delta)
   Green = 2x2 correlation matrix (d and 1S)
   Blue = 3x3 matrix (d, 1S, 2S)
   blacK = 4x4 matrix (d, 1S, 2S, 3S)
   Cyan, where used, is to highlight problems
The x-axis shows the number of energy-levels I included in my fit.
In all plots, error bars are the result of 25 bootstraps.

I've looked at literally hundreds of plots, but the following show some of the most important qualitative results.

• pi_n_statecontamination.gif
  This is a non-bayesian fit to a pi propagator. It demonstrates the possibility of contamination from excited states in the sense of having E0 and E1 be degenerate, with E2 replacing E1 in the fit results. For example, when including 5 or 6 energy levels, we see a two-fold degeneracy for the ground state. Meanwhile, E3 has dropped to the expected value of E2.
  We can also see from this graph how quickly the ground state stabilizes depending on how many correlators are used in the fit. More correlators causes faster convergence to the final value, though it's probably not important since we'll use at least 3 energies in all our fits anyway.
  The excited state seems to be dropping as we add more energy levels, but it's all within errors.
• ro_n_E0.gif
  This is a zoom-in of the ground state for the ro propagator. It shows in detail the convergence to a final value, as well as showing an example of a degeneracy in the energy levels.
• a0_y_E0_drift.gif
  An anomaly was found when using 3 correlators (blue dataset) for the a0 propagator. After the fit had plateaued to a ground state, it suddenly drops further. My guess for an explanation is that with 7 energy levels, some of the higher energy levels were controlling the effects, resulting in the ground state being less important to the fit. This needs to be looked at further, but a possible solution would be to reduce the energy spacing in our priors. In the mean-time, we can limit ourselves to 5 energy levels at most.
• bestfit.gif
  The diamonds show my results for a bayesian fit to all four correlators, using 5 energy levels. The circles are (bayesian?) fit results from Jim (no error bars). The squares are from quenched data.

I'm working on getting more tools automated, such as adjusting the priors to see their effect on fit results. Unfortunately testing takes an excessive amount of computer time.

My program ran for 8 hours before crashing. Unfortunately that was only 1/10 of what I wanted it to do. I'll rerun it (on a smaller set of data) tonight.

The plotting program axis reliably dumps core with comment lines with excessive length. Grepping out comments before plotting is a quick workaround.

I figured out the cause of some of the drift I was seeing with large numbers of energy levels. bigtest_0y50_1.b1_E0.gif shows the drift for the ground state of the b1 state. Since we're using this state to set our scale, it is critical to get an accurate fit. But notice the ground state energy stabilizes, then drifts.
Looking at the bigger picture shows the cause. bigtest_0y50_1.b1_En.gif uses colors just to distinguish between energy levels: red is for E0, E2, ... while blue is for E1, E3, ... Notice that the energy maximum is 12. Apparently above that the computer can't handle the fit function. (Keep in mind it would have to compute exp(-E(t-Tmax)) = exp(-12(1-32)) = 3.6e161 and then balance that with a Z-factor of O(1e-80).) As we add more energies, they therefore get spaced closer together. This can be seen in the chi^2 as well, shown in bigtest_0y50_1.b1_chi2.gif. The chi2 starts out poor (225 for 30 dof (but the covariance matrix used ignored off-diagonal elements, so this chi2 is artificially small)) because only one energy level was used in the fit. Using additional energy levels brings the chi2 to a reasonable level (22/28 for n=2, for example). But when the fit breaks down due to machine limitations, the chi2 rapidly grows to reflect the bad fit.

For those of you who think my plots have too much information, be grateful I don't throw this at you!

I got three streams of my program to go on multiple correlators. All three streams crashed after about 7 hours of running. Since the core-file reported Signal 11, I suspect qcdhome has memory problems.

I realized a possible solution for the problem of N>5 energy levels hitting machine limits is to choose a smaller energy spacing in the priors. Unfortunately the solution to our problem of getting degenerate energy levels is to choose a larger energy spacing in the priors. I'm currently trying to dream up a way of having a small energy spacing, but with a limit to how small it can get. One option would be to do it as (delta E) = (min_spacing) + (eps^2). The minimum spacing just provides a zero-offset to our ordered energies. The only danger is if we encounter energy differences smaller than our prescribed minimum.

A summary of the analysis for different numbers of correlators in our matrix is bigtest_0y50_all_simple.gif. Here, RGBK colors correspond to a 1x1, 2x2, 3x3, or 4x4 correlation matrix in the fit. The horizontal axis comes from adding additional energy levels to the fits. Note that when including only two energy levels, the excited state is poorly determined, though using a larger correlation matrix helps. Including three energy levels gives a much better guess at the first excited state. My 4x4 fit had crashed at this point (hardware errors?), but I'll try running a 4-level fit with that data tonight.

I've been associating a0 with the chi_c0 state, and a1 with the chi_c1 state. Could someone confirm that this is correct?

bigtest_0y50_1_simple.a0.gif
bigtest_0y50_1_simple.a1.gif
bigtest_0y50_1_simple.b1.gif
bigtest_0y50_1_simple.pi.gif
bigtest_0y50_1_simple.ro.gif
bigtest_0y50_2_simple.a0.gif
bigtest_0y50_2_simple.a1.gif
bigtest_0y50_2_simple.b1.gif
bigtest_0y50_2_simple.pi.gif
bigtest_0y50_2_simple.ro.gif
bigtest_0y50_3_simple.a0.gif
bigtest_0y50_3_simple.a1.gif
bigtest_0y50_3_simple.b1.gif
bigtest_0y50_3_simple.pi.gif
bigtest_0y50_3_simple.ro.gif
bigtest_0y50_4_simple.ro.gif

cc_spectrum_orig.gif

Lepage states that chi^2 is of secondary importance in constrained curve fitting. I've seen cases where a really good chi^2 resulted in incorrect parameters, and therefore what we would call a "bad" fit. Goodness of fit comes from the parameters being stable under changes in priors, number of energy levels, etc. Unfortunately, this means we can't really automate any of this, since stability is difficult to quantify. Also, how are we to pick out the "best" fit, when we have several data-points that have plateaued. A possible solution is to fit to the plateau, using an algorithm that is robust to outlier points. A less rigorous method would be to do it by eye.

When producing a spectrum plot, I'm using the etc_c/Jpsi splitting to get a zero offset, and then h_c to set the scale. But these parameters have errors, themselves. So that means the zero offset and scale are unknown to some degree. Do I need to reflect this fact in error bars elsewhere? If so, how?

Taking the plateau energy fits for the 1x1 and 2x2 correlation matrices gives bestfit_q0y50.gif. The errors seem to be much larger for the 2x2 case. The 3x3 case looked to be even worse. This was using bayes statistics on the quenched configurations. I'll submit a run to investigate this without using bayes. It may also help (or not) to use the full correlator matrix (rho=1 instead of rho=0).

Programs crashed after running 5 hours. Again, the core reports signal 11.

I've now written a script to take all the fit result data and generate a list of plots on this webpage. Notation is as follows:
Color denotes prior width to make it easy to see if the priors are over-constraining the fit:

• 1*prior -- Red
• 2*prior -- Green
• 4*prior -- Blue
• 8*prior -- Black

• S -- 1S-1S only
• 1 -- delta-delta only
• 2 -- d-d and 1S-1S
• 3 -- d-d, 1S-1S, 2S-2S

I name my plots in the following format:
runname_rqb#_pp_ee.gif
where

• runname = an identifying string
• r = the "rho" parameter
• q = "q"uenched or "u"nquenched
• b = number of bootstraps (0 = chi^2 errors)
• pp = propagator shown
• ee = energy level

In the following, I'm using only diagonal elements of the correlator matrix, so I constrain the Z-factors to be positive. On the first attempt, this introduced a bug.

The following plots are wrong! (AND their errors aren't plotted correctly)
study_0qy0_a0_E0.gif looks bad
study_0qy0_a0_E1.gif
study_0qy0_a1_E0.gif looks bad
study_0qy0_a1_E1.gif
study_0qy0_b1_E0.gif looks bad
study_0qy0_b1_E1.gif
study_0qy0_pi_E0.gif different states don't agree
study_0qy0_pi_E1.gif
study_0qy0_ro_E0.gif different states don't agree
study_0qy0_ro_E1.gif
study_0uy0_a0_E0.gif
study_0uy0_a0_E1.gif
study_0uy0_a1_E0.gif
study_0uy0_a1_E1.gif
study_0uy0_b1_E0.gif different states don't agree, S seems to be drifting
study_0uy0_b1_E1.gif
study_0uy0_pi_E0.gif
study_0uy0_pi_E1.gif
study_0uy0_ro_E0.gif
study_0uy0_ro_E1.gif

My attempt at a fix: (corrected error bars and z-factor bug)
alpha_0uy0_a0_E0.gif
alpha_0qy0_a0_E0.gif
alpha_0uy0_a0_E1.gif
alpha_0qy0_a0_E1.gif
alpha_0uy0_a1_E0.gif
alpha_0qy0_a1_E0.gif
alpha_0uy0_a1_E1.gif
alpha_0qy0_a1_E1.gif
alpha_0uy0_b1_E0.gif looks good, but should the errors be that large? Maybe so, the effective mass plot oscillatees wildly in that range.
alpha_0qy0_b1_E0.gif A few disturbing artifacts, but mostly within errors. Agrees with eff mass plot.
alpha_0uy0_b1_E1.gif
alpha_0qy0_b1_E1.gif
alpha_0uy0_ro_E0.gif
alpha_0qy0_ro_E0.gifEffective mass plot agrees with d-d, but other fits are way too high. This is very bad.
alpha_0uy0_ro_E1.gif
alpha_0qy0_ro_E1.gif
alpha_0uy0_pi_E0.gif S has drifted down, which is somewhat disturbing, even though it's within errors. The effective mass plot (for the S) agrees with the other data (not with the S). Plotting the S-state here with its excited states suggests a possible cause: apparently several of the excited states are degenerate. It seems they are forcing the ground state lower. alpha_0usy0_pi_E01.gif shows this.
alpha_0qy0_pi_E0.gif A very wide discrepancy between the states. This time, it's the d-d that is singled out. An effective mass plot shows it is the d-d that is correct (E0 should be around 1.42). This was a case where my E0 prior may have been a poor choice. I set E0=1.65+-0.1 because that was close for the b1 state. Apparently it's a poor choice for the pi state. I'll have to see if setting different priors for different states helps. Still, the fact that increasing the prior width didn't help is disturbing. And it selected a value above the prior anyway!
alpha_0uy0_pi_E1.gif
alpha_0qy0_pi_E1.gif

I think the next step is to select better priors for each of the ground states and repeat the tests. I'd like to try auto-selecting the prior by performing a fit to the effective mass plot, but that seems messy, so for now I guess I'll just do it the old-fashioned way and guesstimate it from the eff-mass plot.

The following is done with my guessed priors. I was unable to get a decent guess for a prior (from the eff-mass plot) for the new (unquenched) data with the a0, a1, and b1 correlators. Jim notes that the effective mass plots for the ro and pi are correlated. He claims this is due to a problem with gauge fixing (Coulomb fixing parameter isn't strong enough?), and is starting a new run.

beta_0uy0_a0_E0.gif
beta_0qy0_a0_E0.gif
beta_0uy0_a0_E1.gif
beta_0qy0_a0_E1.gif
beta_0uy0_a1_E0.gif
beta_0qy0_a1_E0.gif
beta_0uy0_a1_E1.gif
beta_0qy0_a1_E1.gif
beta_0uy0_b1_E0.gif Pretty much same as before
beta_0qy0_b1_E0.gif No improvement
beta_0uy0_b1_E1.gif
beta_0qy0_b1_E1.gif
beta_0uy0_ro_E0.gif smaller error bars than before
beta_0qy0_ro_E0.gif significantly better, but still not perfect
beta_0uy0_ro_E1.gif
beta_0qy0_ro_E1.gif
beta_0uy0_pi_E0.gif same as before
beta_0qy0_pi_E0.gif a lot better than before, but still not very good. Needs some close attention. -- It was getting a terrible chi^2, apparently because of a poor choice for the Z-factor. (I'd been setting Z-factors to 1+-2, but setting this to .5+-1 gives a great fit. I'm going to see about specifying my Z-factors separately for each data-set. They could be auto-guessed from the value of the correlator at t=0.
beta_0uy0_pi_E1.gif
beta_0qy0_pi_E1.gif

I ran with my Z-factors being auto-guessed, but saw some bad behavior in the 1S-1S pi ground state. Looking at the bigger picture shows the cause. This problem of running into machine limitations has gone on long enough.

I traced the limitation to be one of numeric underflow. I'm fitting to the function (Z^2)^2*(exp(-Et)+exp(-E(T-t))). But the unquenched data has T=64, and E reaches a maximum value of 11.827511. Note that at t=1, this means exp(-E(T-t))=2.47E-324. But the smallest (in magnitude) number that can be represented in double precision arithmetic is about 5E-324 (see this site for an explanation).

Fortunately, there's a workaround. We can rewrite our function as (Z^2)^2*exp(-ET/2)*2*cosh(E(T/2-t)). Now, the cosh protects us from the underflow. (For a while, anyway. We'll still underflow, but not until E=23.6. And that's enough room that we should be all right for a while.) I wrote and debugged these changes (only for the case of a diagonal correlator matrix, so far). (Actually, I'm not convinced it's debugged. It seems really slow, and if it's slow due to taking a bad path to convergence, then that might be a bug.) The results follow (errors are chi^2 errors):
delta_0uy0_a0_E0.gif Not very good
delta_0qy0_a0_E0.gif Smaller range than before, but states don't agree within errors!
delta_0uy0_a0_E1.gif
delta_0qy0_a0_E1.gif
delta_0uy0_a1_E0.gif States don't agree
delta_0qy0_a1_E0.gif States don't agree (though it's better than before, at least)
delta_0uy0_a1_E1.gif
delta_0qy0_a1_E1.gif
delta_0uy0_b1_E0.gif No real change
delta_0qy0_b1_E0.gif I think this needs to be looked at. States don't agree within errors!
delta_0uy0_b1_E1.gif
delta_0qy0_b1_E1.gif
delta_0uy0_ro_E0.gif Just another beautiful plot that shows off our 1S smearing wavefunction.
delta_0qy0_ro_E0.gif This should probably be looked at further.
delta_0uy0_ro_E1.gif S-state used to be centered here, now it's unable to locate the state. Strange....
delta_0qy0_ro_E1.gif Dramatic improvement, though still not great
delta_0uy0_pi_E0.gif Note that the 1S-1S state finds the ground level with only one energy. Must be a good smearing function. ;) The convergence of the rest of them is so impressive I had to make a zoomed-in plot!
delta_0qy0_pi_E0.gif No improvement, but the error bars increased to encompass the data.
delta_0uy0_pi_E1.gif It's interesting to note that here the 1S-1S can't find the first excited state, despite having found the ground state so perfectly. Again, I take credit for the wavefunction used in the smearing.
delta_0qy0_pi_E1.gif Dramatic improvement. Now we're within errors.

Overall, it seems like this was a pretty good improvement. The pi and ro ground states are impressively well determined. Even their excited states are determined well, provided you give more information than the 1S-1S (which excludes the excited state). But there are still some problems with the b1, a0, and a1 states. Those problems might be related to the gauge fixing problem Jim mentioned. Given that I was unable to get reasonable effective-mass plots for those states, I'm guessing the fault doesn't lie in the fitter.

I'm not really sure what my next step is. Looks like I need data from the new run Jim is doing to see how that effects the b1, a0, and a1 states. In the mean time, I think my focus will be on speeding up my program. The fact that it is running so slowly (taking up to several minutes for a single fit) suggests there may be a bug in the choice it takes for varying parameters. Unfortunately, that's not something I can automate, so it will have to wait until I've gotten some sleep.

I realized my energies seem to be spreading apart further than I'd expect the priors to allow. I take this to mean that my priors aren't allowing Z-factors to get small enough, and aren't constraining the energy spacings enough. Regarding the energy spacings, I'd been reading in E+-s, and setting epsilon = sqrt(E), and sigma = sqrt(s). But in reality, the energy is (epsilon+-sigma)^2 = epsilon^2 +- 2*epsilon*sigma + sigma^2. So while epsilon = sqrt(E), I should set sigma=s/(2*epsilon) to approximate the expected prior. I could do the same for the Z-factors, but since they're fairly unknown in the first place, I'll just guess their value and then let the prior vary them by 100% of that.

For a set of test parameters, my program takes about 100 seconds to perform 100 bootstraps. After several failed attempts to speed it up, I finally got some improvement by caching results of an energy-level computation. (I got the idea from earlier caching the results of a matrix inverse, which had sped up my program by a significant amount.) That got it to do the computation in about 90 seconds. Also caching the results of source_sink lookups dropped it to 85 seconds. Adding a variable called Tmid_t to replace many occurrances of T/2-t dropped it to 80 seconds. I noticed that one (fairly slow) routine, model.d_dp(), was being called 8500 times. I changed the calling routine to call it only 3/4 as often, and it still calls it 8500 times. Yes, I'm confused too. (Turns out that section of code was never called. Oops.) Moving a simple function (jackset) inline to the function that called it saved 30000 function calls and 1 second. I did several other minor optimizations, and eventually got it down to 75 seconds.

I've got a few lingering questions about the 1S-1S data. I assume we can expect it to do a poor job of determining the 2S state. How does this work for the b1 data (since that's P-wave)?

Tonight I'm running again, now that my priors are somewhat narrower. We'll see what happens.

epsilon_0uy0_a0_E0.gif
epsilon_0qy0_a0_E0.gif
epsilon_0uy0_a0_E1.gif
epsilon_0qy0_a0_E1.gif
epsilon_0uy0_a1_E0.gif
epsilon_0qy0_a1_E0.gif
epsilon_0uy0_a1_E1.gif
epsilon_0qy0_a1_E1.gif
epsilon_0uy0_b1_E0.gif
epsilon_0qy0_b1_E0.gif
epsilon_0uy0_b1_E1.gif
epsilon_0qy0_b1_E1.gif
epsilon_0uy0_ro_E0.gif
epsilon_0qy0_ro_E0.gif
epsilon_0uy0_ro_E1.gif
epsilon_0qy0_ro_E1.gif
epsilon_0uy0_pi_E0.gif
epsilon_0qy0_pi_E0.gif
epsilon_0uy0_pi_E1.gif
epsilon_0qy0_pi_E1.gif

The following is with the new run (improved gauge fixing):
zeta_0uy0_a0_E0.gif
zeta_0uy0_a0_E1.gif
zeta_0uy0_a1_E0.gif
zeta_0uy0_a1_E1.gif
zeta_0uy0_b1_E0.gif
zeta_0uy0_b1_E1.gif
zeta_0uy0_ro_E0.gif
zeta_0uy0_ro_E1.gif
zeta_0uy0_pi_E0.gif
zeta_0uy0_pi_E1.gif

Jim notes that our effective mass plots look good for d-d and 1S-d, but not for d-1S or 1S-1S. He recommends fitting to d-d and 1S-d only. It may also help to use a tmax. Even though the fits should ignore the data at large t if it has large errors, it will speed up the fitting program to have fewer points to include.

Since I'm apparently going to need to allow for Z to be positive or negative, I'm looking at improving the speed of that section of the code (similar to what I did above). Adding some of the more critical speed improvements resulted in a 25% speed increase.

I think introducing a minimum energy spacing is necessary if we are to automate our results at all. Degenerate energy levels are causing too much trouble.

Currently I only have expcosh implemented for positive definite Z-factors. I need to allow negative Z-factors there also.

Was at Fermilab most of the day. At night did some confirmation checks on my program.

See this comparison between my results (left) and Jim's results (right).