change output file format to YAML

3 files changed, 187 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();
diff --git a/src/zdab_utils.c b/src/zdab_utils.c
index a689f92..1cde62a 100644
--- a/src/zdab_utils.c
+++ b/src/zdab_utils.c
@@ -49,6 +49,24 @@ void swap_int16(int16_t *val_pt, int count)
     return;
 }
 
+void unpack_mcvx(uint32_t *data, MCVXBank *b)
+{
+    unpack((uint8_t *) data,      "l",&b->cls);
+    unpack((uint8_t *) (data+1),  "l",&b->inc);
+    unpack((uint8_t *) (data+2),  "f",&b->x);
+    unpack((uint8_t *) (data+3),  "f",&b->y);
+    unpack((uint8_t *) (data+4),  "f",&b->z);
+    unpack((uint8_t *) (data+5),  "F",&b->tim);
+    unpack((uint8_t *) (data+7),  "l",&b->rgn);
+    unpack((uint8_t *) (data+8),  "l",&b->idm);
+    unpack((uint8_t *) (data+9),  "l",&b->rg2);
+    unpack((uint8_t *) (data+10), "l",&b->im2);
+    unpack((uint8_t *) (data+12), "f",&b->bnx);
+    unpack((uint8_t *) (data+13), "f",&b->bny);
+    unpack((uint8_t *) (data+14), "f",&b->bnz);
+    unpack((uint8_t *) (data+15), "l",&b->cer);
+}
+
 void unpack_mctk(uint32_t *data, MCTKBank *b)
 {
     unpack((uint8_t *) data,      "l",&b->idp);
diff --git a/src/zdab_utils.h b/src/zdab_utils.h
index 666d2d9..f70482a 100644
--- a/src/zdab_utils.h
+++ b/src/zdab_utils.h
@@ -18,6 +18,17 @@
 #define SWAP_INT16(a,b)
 #endif
 
+/* Status bits for the MC Vertex bank.
+ *
+ * See http://www.lip.pt/~barros/SNO/doc/snoman/companion_frames.html. */
+#define KMCVX_SRC 0x00000001  /* Source */
+#define KMCVX_BOU 0x00000002  /* Boundary */
+#define KMCVX_INT 0x00000004  /* Interaction */
+#define KMCVX_SNK 0x00000008  /* Sink */
+#define KMCVX_PSC 0x00000010  /* Pre-source */
+#define KMCVX_CRE 0x00000100  /* Creation i.e. vertex creates new particle (ignoring Cerenkov photon). */
+#define KMCVX_IPM 0x00000200  /* Indirect MCPM. The MCPM hit comes indirectly from the vertex. */
+
 /* Status bits for the MC Track bank.
  *
  * Bits 4-17 are used to store the "history" of Cerenkov photons. The bit
@@ -52,6 +63,38 @@
 #define KPF_NO_CAL   0x08000000
 #define KPF_BAD_CAL  0x10000000
 
+typedef struct MCVXBank {
+    /* Class:  = 1 Source, = 2 Boundary, = 3 Interaction, = 4 Sink
+	       = 5 Pre-source. */
+    uint32_t cls;
+    /* Interaction code. */
+    uint32_t inc;
+    /* Position X. */
+    float x;
+    /* Position Y. */
+    float y;
+    /* Position Z. */
+    float z;
+    /* Time in ns since MC Generation time (in MC bank). */
+    double tim;
+    /* First region code. */
+    uint32_t rgn;
+    /* First physical media code. */
+    uint32_t idm;
+    /* Second region code. */
+    uint32_t rg2;
+    /* Second physical media code. */
+    uint32_t im2;
+    /* Boundary normal X (only defined if boundary status bit set). */
+    float bnx;
+    /* Boundary normal Y (only defined if boundary status bit set). */
+    float bny;
+    /* Boundary normal Z (only defined if boundary status bit set). */
+    float bnz;
+    /* Number of Cerenkov photons. Only non-zero for CBV vertex. */
+    uint32_t cer;
+} MCVXBank;
+
 typedef struct MCTKBank {
     /* Particle id code. */
     uint32_t idp;
@@ -239,6 +282,7 @@ typedef struct PMTBank {
     uint32_t qrc;
 } PMTBank;
 
+void unpack_mcvx(uint32_t *data, MCVXBank *b);
 void unpack_mctk(uint32_t *data, MCTKBank *b);
 void unpack_ev(uint32_t *data, EVBank *b);
 void unpack_pmt(uint32_t *data, PMTBank *b);