change output file format to YAML

1 files changed, 125 insertions, 47 deletions

diff --git a/src/fit.c b/src/fit.c
index e29e375..293c19c 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -24,7 +24,9 @@
 static int stop = 0;
 static nlopt_opt opt;
 
-#define EV_RECORD    0x45562020 // 'EV  '  (as written to ZDAB file)
+#define EV_RECORD      0x45562020
+#define MCTK_RECORD    0x4d43544b
+#define MCVX_RECORD    0x4d435658
 
 static size_t iter;
 
@@ -5316,17 +5318,64 @@ void sigint_handler(int dummy)
     stop = 1;
 }
 
+int get_event(zebraFile *f, event *ev, bank *b)
+{
+    /* Read all the PMT banks from the zebra file and update `ev`.
+     *
+     * Returns the last return value from get_bank(). */
+    int i, rv;
+    PMTBank bpmt;
+    int id, crate, card, channel;
+
+    for (i = 0; i < MAX_PMTS; i++) {
+        ev->pmt_hits[i].hit = 0;
+    }
+
+    while (1) {
+        rv = next_bank(f, b);
+
+        if (rv == -1) {
+            break;
+        } else if (rv == 1) {
+            /* EOF */
+            break;
+        }
+
+        switch (b->name) {
+        case PMT_RECORD:
+            unpack_pmt(b->data, &bpmt);
+            card = bpmt.pin/1024;
+            crate = (bpmt.pin % 1024)/32;
+            channel = bpmt.pin % 32;
+            id = crate*512 + card*32 + channel;
+            ev->pmt_hits[id].hit = 1;
+            ev->pmt_hits[id].t = bpmt.pt;
+            ev->pmt_hits[id].qhl = bpmt.phl;
+            ev->pmt_hits[id].qhs = bpmt.phs;
+            ev->pmt_hits[id].qlx = bpmt.plx;
+            //ev->pmt_hits[id].flags |= bpmt.pn & KPF_DIS;
+            break;
+        case EV_RECORD:
+        case MAST_RECORD:
+            goto end;
+        }
+    }
+
+end:
+
+    return rv;
+}
+
 int main(int argc, char **argv)
 {
     int i;
     zebraFile *f;
     bank b;
     int rv;
-    PMTBank bpmt;
     EVBank bev;
+    MCTKBank bmctk;
+    MCVXBank bmcvx;
     event ev = {0};
-    int crate, card, channel;
-    int id;
     double fmin;
     double xopt[9];
     char *filename = NULL;
@@ -5405,8 +5454,15 @@ int main(int argc, char **argv)
         exit(1);
     }
 
+    int first_ev = 1;
+    int first_mctk = 1;
+    int first_mcvx = 1;
+    int skip = 0;
     while (1) {
-        rv = next_bank(f, &b);
+        if (!skip)
+            rv = next_bank(f, &b);
+
+        skip = 0;
 
         if (rv == -1) {
             fprintf(stderr, "error getting bank: %s\n", zebra_err);
@@ -5416,56 +5472,78 @@ int main(int argc, char **argv)
             break;
         }
 
-        if (b.name == PMT_RECORD) {
-            unpack_pmt(b.data, &bpmt);
-            card = bpmt.pin/1024;
-            crate = (bpmt.pin % 1024)/32;
-            channel = bpmt.pin % 32;
-            id = crate*512 + card*32 + channel;
-            ev.pmt_hits[id].hit = 1;
-            ev.pmt_hits[id].t = bpmt.pt;
-            ev.pmt_hits[id].qhl = bpmt.phl;
-            ev.pmt_hits[id].qhs = bpmt.phs;
-            ev.pmt_hits[id].qlx = bpmt.plx;
-            //ev.pmt_hits[id].flags |= bpmt.pn & KPF_DIS;
-        } else if (b.name == EV_RECORD) {
-            unpack_ev(b.data, &bev);
-            if (ev.run != -1) {
-                if (ev.run != bev.run || ev.gtid != bev.gtr_id) {
-                    /* New event, so we need to fit the old event. */
-                    rv = fit_event(&ev,xopt,&fmin);
-
-                    if (rv == NLOPT_FORCED_STOP) {
-                        printf("ctrl-c caught. quitting...\n");
-                        break;
-                    }
-
-                    if (fout) {
-                        fprintf(fout, "%" PRIu32 " %10.2f %7.2f %7.2f %7.2f %5.2f %5.2f %6.2f %5.2f %5.2f f() = %7.3e\n",
-                                ev.gtid,
-                                xopt[0],
-                                xopt[1],
-                                xopt[2],
-                                xopt[3],
-                                xopt[4],
-                                xopt[5],
-                                xopt[6],
-                                xopt[7],
-                                xopt[8],
-                                fmin);
-                        fflush(fout);
-                    }
-                }
+        switch (b.name) {
+        case MAST_RECORD:
+            /* New event. */
+            if (fout) fprintf(fout, "-\n");
+            first_ev = 1;
+            first_mctk = 1;
+            first_mcvx = 1;
+            break;
+        case MCVX_RECORD:
+            /* New MC vertex. */
+            unpack_mcvx(b.data, &bmcvx);
+            if (fout) {
+                if (first_mcvx) fprintf(fout, "  mcvx:\n");
+                fprintf(fout, "    -\n");
+                fprintf(fout, "      posx: %.2f\n", bmcvx.x);
+                fprintf(fout, "      posy: %.2f\n", bmcvx.y);
+                fprintf(fout, "      posz: %.2f\n", bmcvx.z);
+            }
+            first_mcvx = 0;
+            break;
+        case MCTK_RECORD:
+            /* New MC track. */
+            unpack_mctk(b.data, &bmctk);
+            if (fout) {
+                if (first_mctk) fprintf(fout, "  mctk:\n");
+                fprintf(fout, "    -\n");
+                fprintf(fout, "      energy: %.2f\n", bmctk.ene);
+                fprintf(fout, "      dirx: %.4f\n", bmctk.drx);
+                fprintf(fout, "      diry: %.4f\n", bmctk.dry);
+                fprintf(fout, "      dirz: %.4f\n", bmctk.drz);
             }
+            first_mctk = 0;
+            break;
+        case EV_RECORD:
+            unpack_ev(b.data, &bev);
             ev.run = bev.run;
             ev.gtid = bev.gtr_id;
+            rv = get_event(f,&ev,&b);
 
-            for (i = 0; i < MAX_PMTS; i++) {
-                ev.pmt_hits[i].hit = 0;
+            if (fit_event(&ev,xopt,&fmin) == NLOPT_FORCED_STOP) {
+                printf("ctrl-c caught. quitting...\n");
+                goto end;
             }
+
+            if (fout) {
+                if (first_ev) fprintf(fout, "  ev:\n");
+                fprintf(fout, "    - gtid: %i\n", ev.gtid);
+                fprintf(fout, "      fit:\n");
+                fprintf(fout, "        -\n");
+                fprintf(fout, "          energy: %.2f\n", xopt[0]);
+                fprintf(fout, "          posx: %.2f\n", xopt[1]);
+                fprintf(fout, "          posy: %.2f\n", xopt[2]);
+                fprintf(fout, "          posz: %.2f\n", xopt[3]);
+                fprintf(fout, "          theta: %.4f\n", xopt[4]);
+                fprintf(fout, "          phi: %.4f\n", xopt[5]);
+                fprintf(fout, "          t0: %.2f\n", xopt[6]);
+                fprintf(fout, "          z1: %.2f\n", xopt[7]);
+                fprintf(fout, "          z2: %.2f\n", xopt[8]);
+                fprintf(fout, "          fmin: %.2f\n", fmin);
+                fflush(fout);
+            }
+
+            first_ev = 0;
+
+            /* Skip reading in the next bank since get_event() already read in
+             * the next bank. */
+            skip = 1;
         }
     }
 
+end:
+
     free_interpolation();
     free_charge();
     pmt_response_free();