update fit to fit each event twice with different quad quantiles

This commit updates the fit program to fit each event and particle hypothesis
twice, once using the normal quad implementation and the other by cutting on
the 10% quantile of times. The first way is much much better when the event is
fully contained since quad will return a really good starting point, and the
second is much better for muons where we want to seed the fit near the entry
point of the muon.

Ideally we would only need a single call and I have an idea of how to update
QUAD to maybe return reasonable guesses in both cases. The idea is to take the
cloud of quad points and find the position and time that has the smallest time
such that it is only a certain Mahalabonis distance from the distribution. This
(I think) corresponds roughly to what I would do by eye where you look at the
distribution of quad points in the cloud and see that it forms a track, and
pick a point at the start of the track.

I started working on this second idea but haven't successfully tested it out
yet.

1 files changed, 30 insertions, 6 deletions

diff --git a/src/fit.c b/src/fit.c
index 6c6b66d..2fb11b9 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -5395,7 +5395,7 @@ double guess_energy(event *ev, double *pos, double *dir, int id)
     return T0;
 }
 
-int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double maxtime)
+int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double maxtime, double f)
 {
     /* Fit the event `ev` under the hypothesis that it's from `n` particles
      * with ids `id`. The best fit parameters are stored in `xopt` upon success
@@ -5446,7 +5446,7 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double
     nlopt_set_min_objective(opt,nlopt_nll2,&fpars);
 
     /* Guess the position and t0 of the event using the QUAD fitter. */
-    status = quad(ev,pos,&t0,10000,0.1);
+    status = quad(ev,pos,&t0,10000,f);
 
     if (status || t0 > GTVALID || NORM(pos) > PSUP_RADIUS) {
         /* If the QUAD fitter fails or returns something outside the PSUP or
@@ -6152,9 +6152,10 @@ size_t particle_combo_to_array(int particle_combo, int *id)
  *     maxtime    - maximum time in seconds to perform the minimization (after the
  *                  "quick" minimization phase)
  *     fmin_best  - the return value from nll_best()
+ *     f          - the quantile to pass to quad()
  *
  * Returns the return value from fit_event2(). */
-int do_fit(event *ev, int *id, size_t n, HDF5Fit *fit_result, double maxtime, double fmin_best)
+int do_fit(event *ev, int *id, size_t n, HDF5Fit *fit_result, double maxtime, double fmin_best, double f)
 {
     int rv;
     struct timeval tv_start, tv_stop;
@@ -6164,7 +6165,7 @@ int do_fit(event *ev, int *id, size_t n, HDF5Fit *fit_result, double maxtime, do
 
     /* Fit the event. */
     gettimeofday(&tv_start, NULL);
-    rv = fit_event2(ev,xopt,&fmin,id,n,maxtime);
+    rv = fit_event2(ev,xopt,&fmin,id,n,maxtime,f);
     gettimeofday(&tv_stop, NULL);
 
     if (rv == NLOPT_FORCED_STOP) return rv;
@@ -6604,7 +6605,27 @@ skip_mc:
             if (particle_combo) {
                 i = particle_combo_to_array(particle_combo,id);
 
-                if (do_fit(&ev,id,i,hdf5_fits+nfits,maxtime,fmin_best) != NLOPT_FORCED_STOP)
+                /* We call do_fit() twice here because we want to fit both when
+                 * quad is called normally, and when we cut on the first 10%
+                 * quantile of times. The first way is much much better when
+                 * the event is fully contained since quad will return a really
+                 * good starting point, and the second is much better for muons
+                 * where we want to seed the fit near the entry point of the
+                 * muon.
+                 *
+                 * FIXME: Ideally we would only need a single call and I have
+                 * an idea of how to update QUAD to maybe return reasonable
+                 * guesses in both cases. The idea is to take the cloud of quad
+                 * points and find the position and time that has the smallest
+                 * time such that it is only a certain Mahalabonis distance
+                 * from the distribution. This (I think) corresponds roughly to
+                 * what I would do by eye where you look at the distribution of
+                 * quad points in the cloud and see that it forms a track, and
+                 * pick a point at the start of the track. */
+                if (do_fit(&ev,id,i,hdf5_fits+nfits,maxtime,fmin_best,1.0) != NLOPT_FORCED_STOP)
+                    nfits++;
+
+                if (do_fit(&ev,id,i,hdf5_fits+nfits,maxtime,fmin_best,0.1) != NLOPT_FORCED_STOP)
                     nfits++;
 
                 if (quit) {
@@ -6625,7 +6646,10 @@ skip_mc:
                             id[k] = particles[combos[j*i+k]];
                         }
 
-                        if (do_fit(&ev,id,i,hdf5_fits+nfits,maxtime,fmin_best) != NLOPT_FORCED_STOP)
+                        if (do_fit(&ev,id,i,hdf5_fits+nfits,maxtime,fmin_best,1.0) != NLOPT_FORCED_STOP)
+                            nfits++;
+
+                        if (do_fit(&ev,id,i,hdf5_fits+nfits,maxtime,fmin_best,0.1) != NLOPT_FORCED_STOP)
                             nfits++;
 
                         if (quit) {