1 files changed, 32 insertions, 4 deletions

diff --git a/src/fit.c b/src/fit.c
index e4fb97f..8a68105 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -43,6 +43,7 @@ typedef struct fitParams {
     int fast;
     int print;
     int id;
+    int charge_only;
 } fitParams;
 
 int particles[] = {
@@ -4975,7 +4976,7 @@ static struct startingParameters {
     { 800.0, 800.0, 800.0},
 };
 
-static double nopt_nll(unsigned int n, const double *x, double *grad, void *params)
+static double nlopt_nll(unsigned int n, const double *x, double *grad, void *params)
 {
     fitParams *fpars = (fitParams *) params;
     double theta, phi;
@@ -5007,7 +5008,7 @@ static double nopt_nll(unsigned int n, const double *x, double *grad, void *para
     v.n = 1;
 
     gettimeofday(&tv_start, NULL);
-    fval = nll(fpars->ev, &v, 1, fpars->dx, fpars->dx_shower, fpars->fast);
+    fval = nll(fpars->ev, &v, 1, fpars->dx, fpars->dx_shower, fpars->fast, fpars->charge_only);
     gettimeofday(&tv_stop, NULL);
 
     long long elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000;
@@ -5118,7 +5119,7 @@ int fit_event(event *ev, double *xopt, double *fmin, int id)
     struct timeval tv_start, tv_stop;
 
     opt = nlopt_create(NLOPT_LN_BOBYQA, 9);
-    nlopt_set_min_objective(opt,nopt_nll,&fpars);
+    nlopt_set_min_objective(opt,nlopt_nll,&fpars);
 
     /* Make a guess as to the energy. Right now we just use a simple
      * approximation that the muon produces approximately 6 photons/MeV.
@@ -5219,6 +5220,7 @@ int fit_event(event *ev, double *xopt, double *fmin, int id)
     fpars.dx_shower = 10.0;
     fpars.fast = 1;
     fpars.print = 0;
+    fpars.charge_only = 0;
     nlopt_set_ftol_abs(opt, 1.0);
     nlopt_set_maxeval(opt, 20);
 
@@ -5313,6 +5315,7 @@ int fit_event(event *ev, double *xopt, double *fmin, int id)
     fpars.dx_shower = 10.0;
     fpars.fast = 0;
     fpars.print = 1;
+    fpars.charge_only = 0;
     nlopt_set_ftol_abs(opt, 1e-5);
     nlopt_set_maxeval(opt, 1000);
 
@@ -5408,6 +5411,25 @@ end:
     return rv;
 }
 
+size_t get_nhit(event *ev)
+{
+    /* Returns the number of PMT hits in event `ev`.
+     *
+     * Note: Only hits on normal PMTs which aren't flagged are counted. */
+    size_t i, nhit;
+
+    nhit = 0;
+    for (i = 0; i < MAX_PMTS; i++) {
+        if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue;
+
+        if (!ev->pmt_hits[i].hit) continue;
+
+        nhit++;
+    }
+
+    return nhit;
+}
+
 int main(int argc, char **argv)
 {
     int i;
@@ -5425,6 +5447,8 @@ int main(int argc, char **argv)
     FILE *fout = NULL;
     int skip_second_event = 0;
     struct timeval tv_start, tv_stop;
+    long long elapsed;
+    double fmin_best;
 
     for (i = 1; i < argc; i++) {
         if (strlen(argv[i]) >= 2 && !strncmp(argv[i], "--", 2)) {
@@ -5581,6 +5605,7 @@ int main(int argc, char **argv)
             if (fout) {
                 if (first_ev) fprintf(fout, "    ev:\n");
                 fprintf(fout, "      - gtid: %i\n", ev.gtid);
+                fprintf(fout, "        nhit: %zu\n", get_nhit(&ev));
                 fprintf(fout, "        fit:\n");
             }
 
@@ -5592,7 +5617,9 @@ int main(int argc, char **argv)
                 }
                 gettimeofday(&tv_stop, NULL);
 
-                long long elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000;
+                elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000;
+
+                fmin_best = nll_best(&ev);
 
                 if (fout) {
                     fprintf(fout, "          %i:\n", particles[i]);
@@ -5607,6 +5634,7 @@ int main(int argc, char **argv)
                     fprintf(fout, "            z2: %.2f\n", xopt[8]);
                     fprintf(fout, "            fmin: %.2f\n", fmin);
                     fprintf(fout, "            time: %lld\n", elapsed);
+                    fprintf(fout, "            psi: %.2f\n", fmin-fmin_best);
                     fflush(fout);
                 }
             }