added support for moench 03

creader/__init__.py

@@ -1 +1,2 @@
+#Make everything from the C extension available
 from _creader import * 

setup.py

@@ -8,7 +8,9 @@ c_ext = setuptools.Extension("_creader",
                         "src/creader_module.c", 
                         "src/cluster_reader.c", 
                         "src/ClusterReader.c",
-                        "src/arr_desc.c"
+                        "src/RawFileReader.c",
+                        "src/arr_desc.c",
+                        "src/raw_reader.c"
                         ],
                     include_dirs=[
                             np.get_include(),"src/"

src/RawFileReader.c

@@ -0,0 +1,125 @@
+#include "RawFileReader.h"
+#include "data_types.h"
+#include "raw_reader.h"
+
+//clang-format off
+typedef struct {
+    PyObject_HEAD FILE *fp;
+    // additional fields for size and decoder?
+    int dtype;
+} RawFileReader;
+//clang-format on
+
+// Constructor: sets the fp to NULL then tries to open the file
+// raises python exception if something goes wrong
+// returned object should mean file is open and ready to read
+static int RawFileReader_init(RawFileReader *self, PyObject *args,
+                              PyObject *Py_UNUSED(kwds)) {
+
+    // Parse file name, accepts string or pathlike objects
+    const char *fname = NULL;
+    PyObject *buf;
+    Py_ssize_t len;
+
+    if (!PyArg_ParseTuple(args, "O&", PyUnicode_FSConverter, &buf))
+        return -1;
+    PyBytes_AsStringAndSize(buf, &fname, &len);
+
+    self->fp = fopen(fname, "rb");
+    // self->n_left = 0;
+
+    // Keep the return code to not return before releasing buffer
+    int rc = 0;
+
+    // Raise python exception using information from errno
+    if (self->fp == NULL) {
+        PyErr_SetFromErrnoWithFilename(PyExc_OSError, fname);
+        rc = -1;
+    }
+    // Release buffer
+    Py_DECREF(buf);
+
+    // Success or fail
+    return rc;
+}
+
+// Custom destructor to make sure we close the file
+static void RawFileReader_dealloc(RawFileReader *self) {
+    if (self->fp) {
+        fclose(self->fp);
+        self->fp = NULL;
+    }
+    Py_TYPE(self)->tp_free((PyObject *)self);
+}
+
+// read method
+static PyObject *RawFileReader_read(RawFileReader *self, PyObject *args) {
+
+    self->dtype = NPY_UINT16;
+    const int ndim = 3;
+    Py_ssize_t n_frames = 1; //default number of frames to read
+    // Py_ssize_t moench_version = 3;
+    // Py_ssize_t analog_digital = 0;
+    if (!PyArg_ParseTuple(args, "|n", &n_frames))
+        return NULL;
+
+    npy_intp dims[3] = {n_frames, 400, 400};
+    PyObject *frames = PyArray_SimpleNew(ndim, dims, self->dtype);
+    PyArray_FILLWBYTE((PyArrayObject *)frames, 0);
+
+    const size_t frame_size = 400 * 400 * 2;
+    char *out_buf = PyArray_DATA((PyArrayObject *)frames);
+    int64_t n_read = read_raw(self->fp, n_frames, frame_size, out_buf);
+
+    if (n_read != n_frames) {
+        // resize the array to match the number of read photons
+        // this will reallocate memory
+
+        // create a new_shape struct on the stack
+        PyArray_Dims new_shape;
+
+        // reuse dims for the shape
+        dims[0] = n_read;
+        new_shape.ptr = dims;
+        new_shape.len = 3;
+
+        // resize the array to match the number of clusters read
+        PyArray_Resize((PyArrayObject *)frames, &new_shape, 3, NPY_ANYORDER);
+    }
+    return frames;
+}
+
+// List all methods in our ClusterFileReader class
+static PyMethodDef RawFileReader_methods[] = {
+    {"read", (PyCFunction)RawFileReader_read, METH_VARARGS,
+     "Read and decode frames"},
+    {NULL, NULL, 0, NULL} /* Sentinel */
+};
+
+// Class defenition
+static PyTypeObject RawFileReaderType = {
+    PyVarObject_HEAD_INIT(NULL, 0).tp_name = "_creader.RawFileReader",
+    .tp_doc = PyDoc_STR("RawFileReader implemented in C"),
+    .tp_basicsize = sizeof(RawFileReader),
+    .tp_itemsize = 0,
+    .tp_flags = Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE,
+    .tp_new = PyType_GenericNew,
+    .tp_dealloc = (destructor)RawFileReader_dealloc,
+    .tp_init = (initproc)RawFileReader_init,
+    .tp_methods = RawFileReader_methods,
+};
+
+PyObject *init_RawFileReader(PyObject *m) {
+    import_array();
+    if (PyType_Ready(&RawFileReaderType) < 0)
+        return NULL;
+
+    Py_INCREF(&RawFileReaderType);
+    if (PyModule_AddObject(m, "RawFileReader", (PyObject *)&RawFileReaderType) <
+        0) {
+        Py_DECREF(&RawFileReaderType);
+        Py_DECREF(m);
+        return NULL;
+    }
+    return m;
+}

src/RawFileReader.h

@@ -0,0 +1,7 @@
+#pragma once
+#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
+#define PY_SSIZE_T_CLEAN
+#include <Python.h>
+#include <numpy/arrayobject.h>
+
+PyObject* init_RawFileReader(PyObject *m);

src/creader_module.c

@@ -4,6 +4,8 @@
 #include <numpy/arrayobject.h>
 
 #include "ClusterReader.h"
+#include "RawFileReader.h"
+
 #include "arr_desc.h"
 #include "data_types.h"
 #include "cluster_reader.h"
@@ -89,5 +91,6 @@ PyMODINIT_FUNC PyInit__creader(void) {
     import_array();
 
     m = init_ClusterFileReader(m);
+    m = init_RawFileReader(m);
     return m;
 }

src/data_types.h

@@ -7,36 +7,52 @@ typedef struct {
     int16_t x;
     int16_t y;
     int32_t data[9];
-} Cluster ;
+} Cluster;
 
 typedef enum {
-  cBottomLeft=0,
+    cBottomLeft = 0,
-  cBottomRight=1,
+    cBottomRight = 1,
-  cTopLeft=2,
+    cTopLeft = 2,
-  cTopRight=3
+    cTopRight = 3
 } corner;
 
 typedef enum {
-  pBottomLeft=0,
+    pBottomLeft = 0,
-  pBottom=1,
+    pBottom = 1,
-  pBottomRight=2,
+    pBottomRight = 2,
-  pLeft=3,
+    pLeft = 3,
-  pCenter=4,
+    pCenter = 4,
-  pRight=5,
+    pRight = 5,
-  pTopLeft=6,
+    pTopLeft = 6,
-  pTop=7,
+    pTop = 7,
-  pTopRight=8
+    pTopRight = 8
 } pixel;
 
-
 typedef struct {
     int32_t tot2;
     int32_t tot3;
     uint32_t c;
-} ClusterAnalysis ;
+} ClusterAnalysis;
-
 
+enum Decoder { MOENCH_03 = 3 };
 
+typedef struct
+{
+  int64_t frame_number;
+  int32_t subframe_number;
+  int32_t packet_number;
+  int64_t bunch_id;
+  int64_t timestamp;
+  int16_t module_id;
+  int16_t row;
+  int16_t column;
+  int16_t reserved;
+  int32_t debug;
+  int16_t round_robin_nb;
+  int8_t detector_type;
+  int8_t header_version;
+  int64_t packets_caught_mask[8];
+} Header;
 
 
 #endif

src/raw_reader.c

@@ -0,0 +1,61 @@
+#include "raw_reader.h"
+#include "data_types.h"
+
+#include <stdlib.h>
+
+int64_t read_raw(FILE *fp, int64_t n_frames, size_t frame_size, char* out_buf) {
+
+    
+
+    char *tmp = malloc(frame_size);
+    Header h;
+
+    int64_t frames_read = 0;
+    while (frames_read < n_frames) {
+        // Read header, return on fail
+        if (fread(&h, sizeof(Header), 1, fp) != 1) {
+            break;
+        }
+
+        // Read frame to temporary buffer
+        if (fread(tmp, frame_size, 1, fp) != 1) {
+            break;
+        }
+
+        // Decode frame and write to numpy output array
+        decode_moench03((uint16_t *)tmp, (uint16_t *)out_buf);
+
+        out_buf += frame_size;
+        frames_read++;
+    }
+
+    free(tmp);
+    return frames_read;
+}
+
+void decode_moench03(const uint16_t *buf, uint16_t *out_buf)
+{
+  int adc_numbers[32] = {12, 13, 14, 15, 12, 13, 14, 15, 8, 9, 10, 11, 8, 9, 10, 11,
+                         4, 5, 6, 7, 4, 5, 6, 7, 0, 1, 2, 3, 0, 1, 2, 3};
+  for (int n_pixel = 0; n_pixel < 5000; n_pixel++)
+  {
+
+    for (int i_sc = 0; i_sc < 32; i_sc++)
+    {
+      // analog part
+      int adc_nr = adc_numbers[i_sc];
+      int col = ((adc_nr * 25) + (n_pixel % 25));
+      int row = 0;
+      if (i_sc / 4 % 2 == 0)
+      {
+        row = 199 - (n_pixel / 25);
+      }
+      else
+      {
+        row = 200 + (n_pixel / 25);
+      }
+      int i_analog = n_pixel * 32 + i_sc;
+      out_buf[row * 400 + col] = buf[i_analog] & 0x3FFF;
+    }
+  }
+}

src/raw_reader.h

@@ -0,0 +1,6 @@
+#pragma once
+#include <stdint.h>
+#include <stdio.h>
+int64_t read_raw(FILE *fp, int64_t n_frames, size_t frame_size, char* out_buf);
+
+void decode_moench03(const uint16_t *buf, uint16_t *out_buf);

tests/test_RawFileReader.py

@@ -0,0 +1,11 @@
+import pytest
+import os, sys
+from creader import RawFileReader
+from fixtures import data_path
+import numpy as np
+
+def test_references_on_read(data_path):
+    fname= data_path/'test_d0_f0_0.raw'
+    r = RawFileReader(fname)
+    frames = r.read(10)
+    assert sys.getrefcount(frames) == 2 #Over counts by one due to call by reference