From e43899cca88746d594985980b49a6da9af24f3e1 Mon Sep 17 00:00:00 2001
From: Erik Frojdh <erik.frojdh@gmail.com>
Date: Fri, 2 Jun 2023 11:29:34 +0200
Subject: [PATCH] added support for moench 03
---
creader/__init__.py | 3 +-
setup.py | 4 +-
src/RawFileReader.c | 125 ++++++++++++++++++++++++++++++++++++
src/RawFileReader.h | 7 ++
src/creader_module.c | 3 +
src/data_types.h | 62 +++++++++++-------
src/raw_reader.c | 61 ++++++++++++++++++
src/raw_reader.h | 6 ++
tests/test_RawFileReader.py | 11 ++++
9 files changed, 257 insertions(+), 25 deletions(-)
create mode 100644 src/RawFileReader.c
create mode 100644 src/RawFileReader.h
create mode 100644 src/raw_reader.c
create mode 100644 src/raw_reader.h
create mode 100644 tests/test_RawFileReader.py
diff --git a/creader/__init__.py b/creader/__init__.py
index 88a46e7..3efaa4f 100644
--- a/creader/__init__.py
+++ b/creader/__init__.py
@@ -1 +1,2 @@
-from _creader import *
\ No newline at end of file
+#Make everything from the C extension available
+from _creader import *
\ No newline at end of file
diff --git a/setup.py b/setup.py
index 8d19abc..bccd187 100644
--- a/setup.py
+++ b/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/"
diff --git a/src/RawFileReader.c b/src/RawFileReader.c
new file mode 100644
index 0000000..51fdc9f
--- /dev/null
+++ b/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;
+}
\ No newline at end of file
diff --git a/src/RawFileReader.h b/src/RawFileReader.h
new file mode 100644
index 0000000..21b2fbf
--- /dev/null
+++ b/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);
\ No newline at end of file
diff --git a/src/creader_module.c b/src/creader_module.c
index 7db53ce..3c082b9 100644
--- a/src/creader_module.c
+++ b/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;
}
diff --git a/src/data_types.h b/src/data_types.h
index 13a071b..71bc019 100644
--- a/src/data_types.h
+++ b/src/data_types.h
@@ -4,39 +4,55 @@
// Data type for the cluster
// has to be extended
typedef struct {
- int16_t x;
- int16_t y;
- int32_t data[9];
-} Cluster ;
+ int16_t x;
+ int16_t y;
+ int32_t data[9];
+} Cluster;
typedef enum {
- cBottomLeft=0,
- cBottomRight=1,
- cTopLeft=2,
- cTopRight=3
+ cBottomLeft = 0,
+ cBottomRight = 1,
+ cTopLeft = 2,
+ cTopRight = 3
} corner;
typedef enum {
- pBottomLeft=0,
- pBottom=1,
- pBottomRight=2,
- pLeft=3,
- pCenter=4,
- pRight=5,
- pTopLeft=6,
- pTop=7,
- pTopRight=8
+ pBottomLeft = 0,
+ pBottom = 1,
+ pBottomRight = 2,
+ pLeft = 3,
+ pCenter = 4,
+ pRight = 5,
+ pTopLeft = 6,
+ pTop = 7,
+ pTopRight = 8
} pixel;
-
typedef struct {
- int32_t tot2;
- int32_t tot3;
- uint32_t c;
-} ClusterAnalysis ;
-
+ int32_t tot2;
+ int32_t tot3;
+ uint32_t c;
+} 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
diff --git a/src/raw_reader.c b/src/raw_reader.c
new file mode 100644
index 0000000..b38ec5f
--- /dev/null
+++ b/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;
+ }
+ }
+}
diff --git a/src/raw_reader.h b/src/raw_reader.h
new file mode 100644
index 0000000..f79fb49
--- /dev/null
+++ b/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);
\ No newline at end of file
diff --git a/tests/test_RawFileReader.py b/tests/test_RawFileReader.py
new file mode 100644
index 0000000..fb6fe49
--- /dev/null
+++ b/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
\ No newline at end of file