diff options
Diffstat (limited to 'src/fit.c')
| -rw-r--r-- | src/fit.c | 36 |
1 files changed, 30 insertions, 6 deletions
@@ -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) { |