fixed conflicts

README.md

@@ -1,18 +1,28 @@
-# python_cpp_example
+# creader
 
-Minimal example building a C++ python extension. 
+Small python package to read cluster and raw files 
 
-Useful links:
+## Getting started
- * [Using NumPy C-API](https://numpy.org/doc/stable/user/c-info.html)
+
+Run
+```bash
+export PYTHONPATH=$PWD:$PYTHONPATH
+make 
+```
+And then have a look at the examples:
+
+- [Reading cluster files](examples/ClusterFile.ipynb)
+- [Reading Moench03/04 raw files](examples/RawFile.ipynb)
 
 
-### Build instructions
+## Build instructions
 
 **Simplified build using make** 
 ```bash
 $ make #build c extension inplace
 ```
 
+Check what is available
 ```
 $ make help
 clean                Remove the build folder and the shared library
@@ -54,3 +64,18 @@ typedef struct {
   int32_t data[9];
 } Cluster ;
 ```
+
+## Running tests
+
+```bash
+#Tell the program where the test data is located. 
+# Can change depending on how you mounted sls_det_storage
+$ export CREADER_TEST_DATA=/mnt/sls_det_storage/moench_data/cluster_reader_test/
+
+$ make test
+
+```
+
+## Useful links:
+ * [Using NumPy C-API](https://numpy.org/doc/stable/user/c-info.html)
+ * [Extending Python with C](https://docs.python.org/3/extending/extending.html)

creader/ClusterFile.py

@@ -2,5 +2,22 @@
 from . import ClusterFileReader
 
 class ClusterFile(ClusterFileReader):
-    def __init__(self, fname):
+    def __init__(self, *args, **kwargs):
-        super().__init__(fname)
+        super().__init__(*args, **kwargs)
+
+
+    def __iter__(self):
+        return self
+    
+    def __next__(self):
+        clusters = self.read()
+        if clusters.size == 0:
+            raise StopIteration
+        else:
+            return clusters
+        
+    def __enter__(self):
+        return self
+    
+    def __exit__(self, exc_type, exc_value, traceback):
+        pass

creader/RawFile.py

@@ -1,29 +1,45 @@
 
 import json
 from pathlib import Path
+import numpy as np
+import re
 from . import RawFileReader
-
+from .enums import DetectorType
 
 class RawFile:
     """
     Generic Raw File reader. Picks up settings from .json master file
     Currently supports: Moench03 =)
     """
-    def __init__(self, fname):
+    def __init__(self, fname, header = False):
         self.findex = 0
         self.fname = fname
+        self.json = False #Master file data comes from json
         fname = Path(fname)
-        if fname.suffix != '.json':
+        if fname.suffix == '.json':
-            raise ValueError("Need a master file in json format")
-        
             with open(fname) as f:
                 self.master = json.load(f)
+                self.json = True
+        elif fname.suffix == '.raw':
+            with open(fname) as f:
+                self._load_raw_master(f)
+        else:
+            raise ValueError(f'{fname.suffix} not supported')
 
         #Figure out which file to open
         if self.master['Detector Type'] == 'Moench' and self.master['Analog Samples'] == 5000:
             #TODO! pass settings to reader
             self._parse_fname()
-            self.reader = RawFileReader(self.data_fname(0,0))
+            self.reader = RawFileReader(self.data_fname(0,0), DetectorType.MOENCH_03, header = header)
+        elif self.master['Detector Type'] == 'ChipTestBoard' and self.master['Analog Samples'] == 5000:
+            self._parse_fname()
+
+            #Do we read analog or analog+digital
+            if self.master['Digital Flag'] == 1:
+                dt = DetectorType.MOENCH_04_AD
+            else: 
+                dt = DetectorType.MOENCH_04_AD
+            self.reader = RawFileReader(self.data_fname(0,0), dt, header = header)
         else:
             raise ValueError('unsupported file')
 
@@ -35,20 +51,128 @@ class RawFile:
         except:
             raise ValueError(f"Could not parse master file name: {self.fname}")
 
+    def _load_raw_master(self, f):
+        self.master = {}
+        lines = f.readlines()
+        it = iter(lines)
+
+        for line in it:
+            if line.startswith("#Frame"):
+                break
+            if line == "\n":
+                continue
+            if line.startswith("Scan Parameters"):
+                while not line.endswith("]\n"):
+                    line += next(it)
+
+            field, value = line.split(":", 1)
+            self.master[field.strip(" ")] = value.strip(" \n")
+
+        frame_header = {}
+        for line in it:
+            field, value = line.split(":", 1)
+            frame_header[field.strip()] = value.strip(" \n")
+
+        self.master["Frame Header"] = frame_header
+        self.master["Version"] = float(self.master["Version"])
+        self._parse_values()
+
+
+    def _parse_values(self):
+        int_fields = set(
+            (
+                "Analog Samples",
+                "Analog Flag",
+                "Digital Flag",
+                "Digital Samples",
+                "Max Frames Per File",
+                "Image Size",
+                "Frame Padding",
+                "Total Frames",
+                "Dynamic Range",
+                "Ten Giga",
+                "Quad",
+                "Number of Lines read out",
+                "Number of UDP Interfaces"
+            )
+        )
+        time_fields = set((
+            "Exptime", 
+            "Exptime1", 
+            "Exptime2",
+            "Exptime3",
+            "GateDelay1",
+            "GateDelay2",
+            "GateDelay3",
+            "SubExptime",#Eiger
+            "SubPeriod", #Eiger
+            "Period"
+
+        ))
+
+        #some fields might not exist for all detectors 
+        #hence using intersection
+        for field in time_fields.intersection(self.master.keys()):
+            self.master[field] = self.to_nanoseconds(self.master[field])
+
+        #Parse bothx .json and .raw master files
+        if self.json:
+            self.master['Image Size'] = self.master["Image Size in bytes"]
+            self.master['Pixels'] = (self.master['Pixels']['x'], self.master['Pixels']['y'])
+            self.master['nmod'] = int(self.master['Geometry']['x']*self.master['Geometry']['y'] )#ports not modules
+            if self.master['Detector Type'] == 'Eiger':
+                self.master['nmod'] = self.master['nmod'] // 2
+        else:
+            for field in int_fields.intersection(self.master.keys()):
+                self.master[field] = int(self.master[field].split()[0])
+            self.master["Pixels"] = tuple(
+                int(i) for i in self.master["Pixels"].strip("[]").split(",")
+            )
+
+        if "Rate Corrections" in self.master:
+            self.master["Rate Corrections"] = (
+                self.master["Rate Corrections"].strip("[]").split(",")
+            )
+            n = len(self.master["Rate Corrections"])
+            assert (
+                self.master["nmod"] == n
+            ), f'nmod from Rate Corrections {n} differs from nmod {self.master["nmod"]}'
+
+        #Parse threshold for Mythen3 (if needed)
+        if "Threshold Energies" in self.master.keys():
+            th = self.master["Threshold Energies"]
+            if isinstance(th, str):
+                th = [int(i) for i in th.strip('[]').split(',')]
+                self.master["Threshold Energies"] = th
+
+    @staticmethod
+    def to_nanoseconds(t):
+        nanoseconds = {"s": 1000 * 1000 * 1000, "ms": 1000 * 1000, "us": 1000, "ns": 1}
+        try:
+            value = re.match(r"(\d+(?:\.\d+)?)", t).group()
+            unit = t[len(value) :]
+            value = int(float(value) * nanoseconds[unit])
+            value = np.timedelta64(value, 'ns')
+        except:
+            raise ValueError(f"Could not convert: {t} to nanoseconds")
+        return value
+
     def data_fname(self, i, findex=0):
             return Path(f"{self.base}_d{i}_f{findex}_{self.run_id}.raw")
 
 
-    def read(self):
+    def read(self, *args):
-        return self.reader.read()
+        return self.reader.read(*args)
+    
+
     
     # Support iteration
     def __iter__(self):
         return self
     
     def __next__(self):
-        frame = self.reader.read()
+        res = self.reader.read()
-        if frame.shape[0] == 0:
+        if res.shape[0] == 0:
             raise StopIteration
     
     # Support with statement

creader/__init__.py

@@ -3,3 +3,5 @@ from _creader import *
 
 from .file_utils import open_file
 from .ClusterFile import ClusterFile
+from .enums import DetectorType
+from .RawFile import RawFile

creader/enums.py

@@ -0,0 +1,8 @@
+from enum import IntEnum
+
+class DetectorType(IntEnum):
+    #TODO! Move defenition to C to avoid mismatch
+    GENERIC = 0
+    MOENCH_03 = 3
+    MOENCH_04_A = 4
+    MOENCH_04_AD = 5

pyproject.toml

@@ -0,0 +1,3 @@
+[build-system]
+requires = ["setuptools", "oldest-supported-numpy/python"]
+build-backend = "setuptools.build_meta"

setup.py

@@ -21,7 +21,13 @@ c_ext = setuptools.Extension("_creader",
 c_ext.language = 'c'
 setuptools.setup(
     name= 'creader',
-    version = '2023.6.1',
+    version = '2023.10.17',
     description = 'Reading cluster files',
+    packages=setuptools.find_packages(exclude=[
+        'tests',
+    ]),
     ext_modules=[c_ext],
+     install_requires=[
+        'numpy',
+    ]
 )

src/ClusterReader.c

@@ -5,9 +5,9 @@
 
 //clang-format off
 typedef struct {
-    PyObject_HEAD 
+    PyObject_HEAD FILE *fp;
-    FILE *fp;
     int n_left;
+    Py_ssize_t chunk;
 } ClusterFileReader;
 //clang-format on
 
@@ -15,21 +15,40 @@ typedef struct {
 // raises python exception if something goes wrong
 // returned object should mean file is open and ready to read
 static int ClusterFileReader_init(ClusterFileReader *self, PyObject *args,
-                                  PyObject *Py_UNUSED(kwds)) {
+                                  PyObject *kwds) {
 
     // Parse file name, accepts string or pathlike objects
     const char *fname = NULL;
-    PyObject* buf;
+    self->n_left = 0;
+    self->chunk = 0;
+    PyObject *fname_obj = NULL;
+    PyObject *fname_bytes = NULL;
     Py_ssize_t len;
 
-    if (!PyArg_ParseTuple(args, "O&", PyUnicode_FSConverter, &buf))
+    static char *kwlist[] = {"fname", "chunk", NULL};
+
+    // if (!PyArg_ParseTuple(args, "O&", PyUnicode_FSConverter, &buf))
+    //     return -1;
+    // PyBytes_AsStringAndSize(buf, &fname, &len);
+
+    if (!PyArg_ParseTupleAndKeywords(args, kwds, "O|n", kwlist, &fname_obj,
+                                     &self->chunk)) {
         return -1;
-    PyBytes_AsStringAndSize(buf, &fname, &len);
+    }
+
+    if (fname_obj != Py_None)
+        if (!PyUnicode_FSConverter(fname_obj, &fname_bytes))
+            return -1;
+
+    PyBytes_AsStringAndSize(fname_bytes, &fname, &len);
+#ifdef CR_VERBOSE
+    printf("Opening: %s\n chunk: %lu\n", fname, self->chunk);
+#endif
 
     self->fp = fopen(fname, "rb");
-    self->n_left = 0;
     
-    //Keep the return code to not return before releasing buffer
+
+    // Keep the return code to not return before releasing buffer
     int rc = 0;
 
     // Raise python exception using information from errno
@@ -37,8 +56,8 @@ static int ClusterFileReader_init(ClusterFileReader *self, PyObject *args,
         PyErr_SetFromErrnoWithFilename(PyExc_OSError, fname);
         rc = -1;
     }
-    //Release buffer
+    // Release buffer
-    Py_DECREF(buf);
+    Py_DECREF(fname_bytes);
 
     // Success or fail
     return rc;
@@ -57,68 +76,76 @@ static void ClusterFileReader_dealloc(ClusterFileReader *self) {
 }
 
 // read method
-static PyObject *ClusterFileReader_read(ClusterFileReader *self, PyObject *args) {
+static PyObject *ClusterFileReader_read(ClusterFileReader *self,
+                                        PyObject *args) {
 
     const int ndim = 1;
     Py_ssize_t size = 0;
-   PyObject *noise_obj;
+    PyObject *noise_obj = NULL;
-   if (!PyArg_ParseTuple(args, "nO", &size,&noise_obj)) {
+    PyObject *noise_array = NULL;
-     PyErr_SetString(
+    if (!PyArg_ParseTuple(args, "|nO", &size, &noise_obj)) {
-		      PyExc_TypeError,
+        PyErr_SetString(PyExc_TypeError, "Could not parse args.");
-		      "Could not parse args.");
         return NULL;
     }
+
+    // Fall back on object default/config
+    if (size == 0)
+        size = self->chunk;
+
     npy_intp dims[] = {size};
 
-    // Create two numpy arrays from the passed objects, if possible numpy will
+    // If possible numpy will
+>>>>>>> e1f22fd014eb86d40e0ecb78bd65cb7a9bef71bc
     // use the underlying buffer, otherwise it will create a copy, for example
     // if data type is different or we pass in a list. The
     // NPY_ARRAY_C_CONTIGUOUS flag ensures that we have contiguous memory.
-    PyObject *noise_array = PyArray_FROM_OTF(noise_obj, NPY_DOUBLE, NPY_ARRAY_C_CONTIGUOUS);
-    int nx=0,ny=0;
-    double *noise_map=NULL;
 
+#ifdef CR_VERBOSE
+    printf("Getting ready to read: %lu clusters. Noise map: %p\n", size,
+           noise_obj);
+#endif
 
-    // If parsing of a or b fails we throw an exception in Python
+    // If the user passed a noise map we fetch a pointer to that array as well
-    if (noise_array ) {
+    int nx = 0, ny = 0;
+    double *noise_map = NULL;
+    if (noise_obj) {
+        noise_array =
+            PyArray_FROM_OTF(noise_obj, NPY_DOUBLE, NPY_ARRAY_C_CONTIGUOUS);
 
         int ndim_noise = PyArray_NDIM((PyArrayObject *)(noise_array));
         npy_intp *noise_shape = PyArray_SHAPE((PyArrayObject *)(noise_array));
 
-    
         // For the C++ function call we need pointers (or another C++ type/data
         // structure)
 
         noise_map = (double *)(PyArray_DATA((PyArrayObject *)(noise_array)));
 
 
-
         /* for (int i=0; i< ndim_noise; i++) { */
-    /*   printf("Dimension %d size %d pointer \n",i,noise_shape[i], noise_map); */
+        /*   printf("Dimension %d size %d pointer \n",i,noise_shape[i],
+         * noise_map); */
 
         /* } */
 
-      if (ndim_noise==2) {
 
-	nx=noise_shape[0];
+        if (ndim_noise == 2) {
-	ny=noise_shape[1];
 
-	//printf("Noise map found size %d %d %d\n",nx,ny,noise_map);
+            nx = noise_shape[0];
+            ny = noise_shape[1];
 
+            //	printf("Noise map found size %d %d %d\n",nx,ny,noise_map);
 
         } else {
-	nx=0;
+            nx = 0;
-	if (ndim_noise==1) 
+            if (ndim_noise == 1)
-	  nx=noise_shape[0];
+                nx = noise_shape[0];
-	ny=0;
+            ny = 0;
             noise_map = NULL;
-	//printf("NO Noise map found %d %d %d %d\n",ndim_noise,nx,ny,noise_map);
+            //	printf("NO Noise map found %d %d %d
+            //%d\n",ndim_noise,nx,ny,noise_map);
         }
-
     }
 
-
-
     // Create an uninitialized numpy array
     PyObject *clusters = PyArray_SimpleNewFromDescr(ndim,  dims, cluster_dt());
 
@@ -132,7 +159,8 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, PyObject *args)
     // Here goes the looping, removing frame numbers etc.
     int n_read = 0;
     if (noise_map)
-      read_clusters_with_cut(self->fp, size, buf, &self->n_left,noise_map, nx, ny);
+        n_read = read_clusters_with_cut(self->fp, size, buf, &self->n_left,
+                                        noise_map, nx, ny);
     else
        n_read = read_clusters(self->fp, size, buf, &self->n_left);
 

src/RawFileReader.c

@@ -1,7 +1,7 @@
 #include "RawFileReader.h"
+#include "arr_desc.h"
 #include "data_types.h"
 #include "raw_reader.h"
-#include "arr_desc.h"
 
 #include <stdbool.h>
 
@@ -11,6 +11,7 @@ typedef struct {
     // additional fields for size and decoder?
     int dtype;
     bool read_header;
+    int detector;
 } RawFileReader;
 //clang-format on
 
@@ -28,11 +29,12 @@ static int RawFileReader_init(RawFileReader *self, PyObject *args,
 
     // Should we read the header
     self->read_header = false;
+    self->detector = DT_GENERIC;
 
-    static char *kwlist[] = {"fname", "header", NULL};
+    static char *kwlist[] = {"fname", "detector_type", "header", NULL};
 
-    if (!PyArg_ParseTupleAndKeywords(args, kwds, "O|p", kwlist, &fname_obj,
+    if (!PyArg_ParseTupleAndKeywords(args, kwds, "Oi|p", kwlist, &fname_obj,
-                                     &self->read_header)) {
+                                     &self->detector, &self->read_header)) {
         return -1;
     }
 
@@ -41,11 +43,13 @@ static int RawFileReader_init(RawFileReader *self, PyObject *args,
             return -1;
 
     PyBytes_AsStringAndSize(fname_bytes, &fname, &len);
-    printf("%s\n read_header: %d\n", fname, self->read_header);
 
+#ifdef CR_VERBOSE
+    printf("fname: %s\n read_header: %d detector type: %d\n", fname,
+           self->read_header, self->detector);
+#endif
     self->fp = fopen(fname, "rb");
 
-
     // Keep the return code to not return before releasing buffer
     int rc = 0;
 
@@ -76,27 +80,46 @@ 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);
+    char *out_buf = PyArray_DATA((PyArrayObject *)frames);
+    PyObject *digital_frames = NULL;
+    char *digital_out = NULL;
+
     PyArray_FILLWBYTE((PyArrayObject *)frames, 0);
 
-    //Optional return the header
+    // Optional return the header
     PyObject *header = NULL;
-    char* header_out = NULL;
+    char *header_out = NULL;
-    if(self->read_header){
+    if (self->read_header) {
         header = PyArray_SimpleNewFromDescr(1, dims, frame_header_dt());
         header_out = PyArray_DATA((PyArrayObject *)header);
     }
 
-    const size_t frame_size = 400 * 400 * 2;
+    // Both analog and digital data
-    char *out_buf = PyArray_DATA((PyArrayObject *)frames);
+    if (self->detector == DT_MOENCH_04_AD) {
-    int64_t n_read =
+        digital_frames = PyArray_SimpleNew(ndim, dims, NPY_UINT8);
-        read_raw(self->fp, n_frames, frame_size, out_buf, header_out);
+        PyArray_FILLWBYTE((PyArrayObject *)digital_frames, 0);
+        digital_out = PyArray_DATA((PyArrayObject *)digital_frames);
+    }
+
+    int64_t n_read = 0;
+
+    switch (self->detector) {
+    case DT_MOENCH_03:
+        n_read = read_raw_m03(self->fp, n_frames, out_buf, header_out);
+        break;
+    case DT_MOENCH_04_A:
+    case DT_MOENCH_04_AD:
+        n_read =
+            read_raw_m04(self->fp, n_frames, out_buf, digital_out, header_out);
+        break;
+    default:
+        break;
+    }
 
     if (n_read != n_frames) {
         // resize the array to match the number of read photons
@@ -113,20 +136,41 @@ static PyObject *RawFileReader_read(RawFileReader *self, PyObject *args) {
         // resize the array to match the number of clusters read
         PyArray_Resize((PyArrayObject *)frames, &new_shape, 1, NPY_ANYORDER);
 
-        //if we also read header we need to reshape the header
+        if(digital_frames){
-        if(self->read_header){
+            PyArray_Resize((PyArrayObject *)digital_frames, &new_shape, 1, NPY_ANYORDER);
-            new_shape.len = 1;
-            PyArray_Resize((PyArrayObject *)header, &new_shape, 1, NPY_ANYORDER);
         }
 
+        // if we also read header we need to reshape the header
+        if (self->read_header) {
+            new_shape.len = 1;
+            PyArray_Resize((PyArrayObject *)header, &new_shape, 1,
+                           NPY_ANYORDER);
         }
-    PyObject *ret = frames;
+
-    if(self->read_header){
+        
-        ret = PyTuple_Pack(2, frames, header);
+    }
-        Py_DECREF(header);
+
+    // Build up a tuple with the return values
+    PyObject *ret = PyTuple_Pack(1, frames);
+    if (self->detector == DT_MOENCH_04_AD){
+        _PyTuple_Resize(&ret, 2);
+        PyTuple_SET_ITEM(ret, 1, digital_frames);
+    }
+
+    if (self->read_header){
+        Py_ssize_t old_size = PyTuple_GET_SIZE(ret);
+        _PyTuple_Resize(&ret, old_size+1);
+        PyTuple_SET_ITEM(ret, old_size, header);
+    }
+
+    // if we only have one item in the tuple lets return it instead of the tuple
+    if(PyTuple_GET_SIZE(ret) == 1){
+        Py_DECREF(ret);
+        return frames;
+    }else{
         Py_DECREF(frames);
-    }
         return ret;
+    }
     
 }
 

src/cluster_reader.c

@@ -26,7 +26,6 @@ int read_clusters(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_left) {
         // keep on reading frames and photons until reaching n_clusters
         while (fread(&iframe, sizeof(iframe), 1, fp)) {
 
-	  
             if (fread(&nph, sizeof(nph), 1, fp)) {
                 if (nph > n_clusters - nph_read)
                     nn = n_clusters - nph_read;
@@ -39,18 +38,14 @@ int read_clusters(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_left) {
             }
             if (nph_read >= n_clusters)
                 break;
-
-
-
         }
     }
     assert(nph_read <= n_clusters); // sanity check in debug mode
     return nph_read;
 }
 
-
+int read_clusters_with_cut(FILE *fp, int64_t n_clusters, Cluster *buf,
-
+                           int *n_left, double *noise_map, int nx, int ny) {
-int read_clusters_with_cut(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_left, double *noise_map, int nx, int ny) {
     int iframe = 0;
     int nph = *n_left;
 
@@ -58,11 +53,12 @@ int read_clusters_with_cut(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_le
 
     int32_t tot2[4], t2max, tot1;
     int32_t val, tot3;
-    Cluster *ptr=buf;
+    Cluster *ptr = buf;
-    int good=1;
+    int good = 1;
     double noise;
     // read photons left from previous frame
     if (nph) {
+<<<<<<< HEAD
 	for (int iph=0; iph<nph; iph++) {
 	  // read photons 1 by 1
 	  fread((void *)(ptr), sizeof(Cluster), 1, fp);
@@ -89,11 +85,41 @@ int read_clusters_with_cut(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_le
 	  if (nph_read >= n_clusters)
 	    break;
 	}
+=======
+        for (int iph = 0; iph < nph; iph++) {
+            // read photons 1 by 1
+            fread((void *)(ptr), sizeof(Cluster), 1, fp);
+            good = 1;
+            if (noise_map) {
+                if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 && ptr->y < ny) {
+                    tot1 = ptr->data[4];
+                    analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL, NULL, NULL,
+                                    NULL, NULL, NULL, NULL, NULL);
+                    noise = noise_map[ptr->y * nx + ptr->x];
+                    if (tot1 > noise && t2max > 2 * noise && tot3 > 3 * noise) {
+                        ;
+                    } else
+                        good = 0;
+                } else {
+                    printf("Bad pixel number %d %d\n", ptr->x, ptr->y);
+                    good = 0;
+                }
+            }
+            if (good) {
+                ptr++;
+                nph_read++;
+                (*n_left)--;
+            }
+            if (nph_read >= n_clusters)
+                break;
+        }
+>>>>>>> e1f22fd014eb86d40e0ecb78bd65cb7a9bef71bc
     }
 
     if (nph_read < n_clusters) {
         // keep on reading frames and photons until reaching n_clusters
         while (fread(&iframe, sizeof(iframe), 1, fp)) {
+<<<<<<< HEAD
 	  //printf("%d\n",nph_read);
 
 	  if (fread(&nph, sizeof(nph), 1, fp)) {
@@ -131,12 +157,55 @@ int read_clusters_with_cut(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_le
 	  if (nph_read >= n_clusters)
 	    break;
 	}
+=======
+            // printf("%d\n",nph_read);
+
+            if (fread(&nph, sizeof(nph), 1, fp)) {
+                // printf("** %d\n",nph);
+                *n_left = nph;
+                for (int iph = 0; iph < nph; iph++) {
+                    // read photons 1 by 1
+                    fread((void *)(ptr), sizeof(Cluster), 1, fp);
+                    good = 1;
+                    if (noise_map) {
+                        if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 &&
+                            ptr->y < ny) {
+                            tot1 = ptr->data[4];
+                            analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL,
+                                            NULL, NULL, NULL, NULL, NULL, NULL,
+                                            NULL);
+                            noise = noise_map[ptr->y * nx + ptr->x];
+                            if (tot1 > noise && t2max > 2 * noise &&
+                                tot3 > 3 * noise) {
+                                ;
+                            } else
+                                good = 0;
+                        } else {
+                            printf("Bad pixel number %d %d\n", ptr->x, ptr->y);
+                            good = 0;
                         }
-    //printf("%d\n",nph_read);
+                    }
+                    if (good) {
+                        ptr++;
+                        nph_read++;
+                        (*n_left)--;
+                    }
+                    if (nph_read >= n_clusters)
+                        break;
+                }
+            }
+
+            if (nph_read >= n_clusters)
+                break;
+        }
+>>>>>>> e1f22fd014eb86d40e0ecb78bd65cb7a9bef71bc
+    }
+    // printf("%d\n",nph_read);
     assert(nph_read <= n_clusters); // sanity check in debug mode
     return nph_read;
 }
 
+<<<<<<< HEAD
 
 
 
@@ -144,6 +213,9 @@ int read_clusters_with_cut(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_le
 
 
 int analyze_clusters(int64_t n_clusters, int32_t *cin, ClusterAnalysis *co, int csize) {
+=======
+int analyze_clusters(int64_t n_clusters, Cluster *cin, ClusterAnalysis *cout) {
+>>>>>>> e1f22fd014eb86d40e0ecb78bd65cb7a9bef71bc
 
     int32_t tot2[4], t2max;
     char quad;
@@ -153,6 +225,7 @@ int analyze_clusters(int64_t n_clusters, int32_t *cin, ClusterAnalysis *co, int
     //printf("csize is %d\n",csize);
     int ret;
     for (int ic = 0; ic < n_clusters; ic++) {
+<<<<<<< HEAD
         
       switch (csize) {
       case 2:
@@ -175,10 +248,22 @@ int analyze_clusters(int64_t n_clusters, int32_t *cin, ClusterAnalysis *co, int
       /* if (tot<=0) */
       /* printf("%d %d %d %d %d %d\n",ic,(cin+ic)->x, (cin+ic)->y, */
       /*   (cout+ic)->c, (cout+ic)->tot2, (cout+ic)->tot3); */
+=======
+
+        analyze_cluster(*(cin + ic), &t2max, &tot3, &quad, NULL, NULL, NULL,
+                        NULL, NULL, NULL, NULL, NULL);
+
+        (cout + ic)->c = quad;
+        (cout + ic)->tot2 = t2max;
+        (cout + ic)->tot3 = tot3;
+        // printf("%d %d %d %d %d %d\n",ic,(cin+ic)->x, (cin+ic)->y,
+        // (cout+ic)->c, (cout+ic)->tot2, (cout+ic)->tot3);
+>>>>>>> e1f22fd014eb86d40e0ecb78bd65cb7a9bef71bc
     }
     return nc;
 }
 
+<<<<<<< HEAD
 int analyze_cluster(Cluster cl, int32_t *t2, int32_t *t3, char *quad, double *eta2x, double *eta2y, double *eta3x, double *eta3y) {
 
   return analyze_data(cl.data, t2, t3, quad, eta2x, eta2y, eta3x, eta3y);
@@ -189,6 +274,14 @@ int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, double *et
 
 
   int ok=1;
+=======
+int analyze_cluster(Cluster cin, int32_t *t2, int32_t *t3, char *quad,
+                    double *eta2x, double *eta2y, double *eta3x, double *eta3y,
+                    double *eta2Lx, double *eta2Ly, double *eta3Xx,
+                    double *eta3Xy) {
+
+    int ok = 1;
+>>>>>>> e1f22fd014eb86d40e0ecb78bd65cb7a9bef71bc
 
     int32_t tot2[4], t2max;
     char c;
@@ -199,6 +292,7 @@ int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, double *et
         tot2[i] = 0;
     // t2max=0;
     for (int ix = 0; ix < 3; ix++) {
+<<<<<<< HEAD
       for (int iy = 0; iy < 3; iy++) {
 	val = data[iy * 3 + ix];
 	//	printf ("%d ",data[iy * 3 + ix]);
@@ -213,12 +307,25 @@ int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, double *et
 	  tot2[cTopRight] += val;
       }
       //	printf ("\n");
+=======
+        for (int iy = 0; iy < 3; iy++) {
+            val = cin.data[iy * 3 + ix];
+            tot3 += val;
+            if (ix <= 1 && iy <= 1)
+                tot2[0] += val;
+            if (ix >= 1 && iy <= 1)
+                tot2[1] += val;
+            if (ix <= 1 && iy >= 1)
+                tot2[2] += val;
+            if (ix >= 1 && iy >= 1)
+                tot2[3] += val;
+        }
+>>>>>>> e1f22fd014eb86d40e0ecb78bd65cb7a9bef71bc
     }
     //printf ("\n");
 
-
-
     if (t2 || quad) {
+<<<<<<< HEAD
       t2max = -1000;
       c = 0;
       for (int i = 0; i < 4; i++) {
@@ -227,12 +334,23 @@ int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, double *et
 	  c = i;
 	}
       }
+=======
+        t2max = tot2[0];
+        c = cBottomLeft;
+        for (int i = 1; i < 4; i++) {
+            if (tot2[i] > t2max) {
+                t2max = tot2[i];
+                c = i;
+            }
+        }
+>>>>>>> e1f22fd014eb86d40e0ecb78bd65cb7a9bef71bc
     }
     if (quad)
         *quad = c;
     if (t2)
         *t2 = t2max;
     if (t3)
+<<<<<<< HEAD
       *t3 = tot3;
     
     if (eta2x || eta2y) {
@@ -290,8 +408,9 @@ int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, double *et
     /* } */
       
 
+=======
+        *t3 = tot3;
+>>>>>>> e1f22fd014eb86d40e0ecb78bd65cb7a9bef71bc
 
     return ok;
-
-
 }

src/creader_module.c

@@ -7,8 +7,8 @@
 #include "RawFileReader.h"
 
 #include "arr_desc.h"
-#include "data_types.h"
 #include "cluster_reader.h"
+#include "data_types.h"
 
 /* static PyObject *clusterize(PyObject *Py_UNUSED(self), PyObject *args) {
 
@@ -172,13 +172,14 @@ static PyObject *clusterize(PyObject *Py_UNUSED(self), PyObject *args) {
 static PyObject *get_cluster_dt(PyObject *Py_UNUSED(self), PyObject *args) {
     if (!PyArg_ParseTuple(args, ""))
         return NULL;
-    return (PyObject*)cluster_dt();
+    return (PyObject *)cluster_dt();
 }
 
-static PyObject *get_frame_header_dt(PyObject *Py_UNUSED(self), PyObject *args) {
+static PyObject *get_frame_header_dt(PyObject *Py_UNUSED(self),
+                                     PyObject *args) {
     if (!PyArg_ParseTuple(args, ""))
         return NULL;
-    return (PyObject*)frame_header_dt();
+    return (PyObject *)frame_header_dt();
 }
 
 // Module docstring, shown as a part of help(creader)

src/data_types.h

@@ -35,7 +35,12 @@ typedef struct {
     double etay;
 } ClusterAnalysis;
 
-enum Decoder { MOENCH_03 = 3 };
+enum DetectorType {
+    DT_GENERIC = 0, 
+    DT_MOENCH_03 = 3,
+    DT_MOENCH_04_A = 4,
+    DT_MOENCH_04_AD = 5 
+};
 
 typedef struct
 {
@@ -56,4 +61,5 @@ typedef struct
 } Header;
 
 
+
 #endif

src/raw_reader.c

@@ -1,13 +1,14 @@
 #include "raw_reader.h"
 #include "data_types.h"
 
+#include <stdbool.h>
 #include <stdlib.h>
 #include <string.h>
 
-int64_t read_raw(FILE *fp, int64_t n_frames, size_t frame_size, char* out_buf, Header* header_out) {
+int64_t read_raw_m03(FILE *fp, int64_t n_frames, char *frame_out,
-
+                     Header *header_out) {
-    
 
+    const size_t frame_size = 400 * 400 * 2;
     char *tmp = malloc(frame_size);
     Header h;
 
@@ -17,7 +18,7 @@ int64_t read_raw(FILE *fp, int64_t n_frames, size_t frame_size, char* out_buf, H
         if (fread(&h, sizeof(Header), 1, fp) != 1) {
             break;
         }
-        if(header_out){
+        if (header_out) {
             memcpy(header_out, &h, sizeof(Header));
             header_out++;
         }
@@ -28,9 +29,9 @@ int64_t read_raw(FILE *fp, int64_t n_frames, size_t frame_size, char* out_buf, H
         }
 
         // Decode frame and write to numpy output array
-        decode_moench03((uint16_t *)tmp, (uint16_t *)out_buf);
+        decode_moench03((uint16_t *)tmp, (uint16_t *)frame_out);
 
-        out_buf += frame_size;
+        frame_out += frame_size;
         frames_read++;
     }
 
@@ -38,25 +39,70 @@ int64_t read_raw(FILE *fp, int64_t n_frames, size_t frame_size, char* out_buf, H
     return frames_read;
 }
 
-void decode_moench03(const uint16_t *buf, uint16_t *out_buf)
+int64_t read_raw_m04(FILE *fp, int64_t n_frames, char *frame_out,
-{
+                     char *digital_out, Header *header_out) {
-  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++)
+    const size_t frame_size = 400 * 400 * 2;
-    {
+    const size_t digital_frame_size = 5000 * 8;
+    char *tmp = malloc(frame_size);
+    char *digital_tmp = malloc(digital_frame_size);
+    Header h;
+
+    int64_t frames_read = 0;
+#ifdef CR_VERBOSE
+    printf("MOENCH04: digital_out: %p\n", digital_out);
+#endif
+    while (frames_read < n_frames) {
+        // Read header to temp buffer, return on fail
+        if (fread(&h, sizeof(Header), 1, fp) != 1) {
+            break;
+        }
+        if (header_out) {
+            memcpy(header_out, &h, sizeof(Header));
+            header_out++;
+        }
+
+        // Read frame to temporary buffer
+        if (fread(tmp, frame_size, 1, fp) != 1) {
+            break;
+        }
+
+        if (digital_out) {
+            if (fread(digital_tmp, digital_frame_size, 1, fp) != 1) {
+                break;
+            }
+        }
+
+        // Decode frame and write to numpy output array
+        // decode_moench03((uint16_t *)tmp, (uint16_t *)frame_out);
+        decode_moench04(tmp, digital_tmp, frame_out, digital_out);
+
+        frame_out += frame_size;
+        if (digital_out)
+            digital_out += 400*400; //one byte per pixel
+
+        frames_read++;
+    }
+
+    free(tmp);
+    free(digital_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)
+            if (i_sc / 4 % 2 == 0) {
-      {
                 row = 199 - (n_pixel / 25);
-      }
+            } else {
-      else
-      {
                 row = 200 + (n_pixel / 25);
             }
             int i_analog = n_pixel * 32 + i_sc;
@@ -64,3 +110,52 @@ void decode_moench03(const uint16_t *buf, uint16_t *out_buf)
         }
     }
 }
+
+void decode_moench04(const uint16_t *analog_data, const uint64_t *digital_data,
+                     uint16_t *analog_frame, uint8_t *digital_frame) {
+
+    int adc_numbers[32] = {9,  8,  11, 10, 13, 12, 15, 14, 1,  0,  3,
+                           2,  5,  4,  7,  6,  23, 22, 21, 20, 19, 18,
+                           17, 16, 31, 30, 29, 28, 27, 26, 25, 24};
+
+    int dbits_bit[32] = {-1, -1, -1, -1, -1, -1, 1,  3,  5,  7,  -1,
+                         -1, -1, -1, -1, -1, 62, 60, 58, 56, 54, 52,
+                         50, 48, 63, 61, 59, 57, 55, 53, 51, 49};
+
+    for (int n_pixel = 0; n_pixel < 5000; n_pixel++) {
+
+        uint64_t current_dbits = digital_data[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 % 16) * 25) + (n_pixel % 25);
+            int row = 0;
+            if (i_sc < 16) {
+                row = 199 - (n_pixel / 25);
+            } else {
+                row = 200 + (n_pixel / 25);
+            }
+            int i_analog = n_pixel * 32 + i_sc;
+            analog_frame[row * 400 + col] = analog_data[i_analog] & 0x3FFF;
+
+            // if we pass digital data to be decoded
+            if (digital_frame) {
+                int bit_loc = dbits_bit[adc_nr];
+
+                if (bit_loc < 0) {
+                    digital_frame[row * 400 + col] = 2;
+                } else {
+                    uint8_t dbit = (current_dbits >> bit_loc) & 1U;
+
+                    if ((6 <= adc_nr) & (adc_nr <= 9)) {
+                        digital_frame[row * 400 + col] = dbit;
+                    } else {
+                        digital_frame[row * 400 + col] = 1 - dbit;
+                    }
+                }
+            } else {
+                // digital_frame[row * 400 + col] = 2;
+            }
+        }
+    }
+}

src/raw_reader.h

@@ -1,8 +1,27 @@
 #pragma once
+#include <stdbool.h>
 #include <stdint.h>
 #include <stdio.h>
 #include "data_types.h"
 
-int64_t read_raw(FILE *fp, int64_t n_frames, size_t frame_size, char* out_buf, Header* header_out);
+int64_t read_raw_m03(
+    FILE *fp, 
+    int64_t n_frames, 
+    char* frame_out, 
+    Header* header_out
+);
+
+int64_t read_raw_m04(
+    FILE *fp, 
+    int64_t n_frames, 
+    char* frame_out, 
+    char* digital_out,
+    Header* header_out
+);
 
 void decode_moench03(const uint16_t *buf, uint16_t *out_buf);
+
+void decode_moench04(const uint16_t *analog_data,
+                     const uint64_t *digital_data,
+                     uint16_t *analog_frame,
+                     uint8_t *digital_frame);

tests/test_ClusterReader.py

@@ -85,10 +85,40 @@ def test_read_file_with_37_frames(data_path):
 def test_read_file_with_37_frames_in_chunks(data_path):
     #File shoud contain 37 frames with 5 clusters each
     #Full spec in utils/write_test_data.py
-    fname= (data_path/'37frames_with_5_clusters.clust').as_posix()
+    fname= data_path/'37frames_with_5_clusters.clust'
     r = ClusterFileReader(fname)
     total_clusters = 0
     while (clusters:=r.read(7)).size:
         total_clusters += clusters.size
     assert total_clusters == 185
     
+
+def test_read_file_with_noise_mask(data_path):
+    #No mask
+    fname= data_path/'noise_test.clust'
+    r = ClusterFileReader(fname)
+    cl = r.read(85) #file contains 70 clusters
+    assert cl.size == 70
+
+    #noise mask with zeros
+    noise_cut = np.zeros((400,400))
+    r = ClusterFileReader(fname)
+    cl = r.read(85, noise_cut)
+    assert cl.size == 70
+
+    #only pixel 80, 133 above noise
+    noise_cut[:] = 100
+
+    #TODO! Agree on orientation of noise mask!
+    # noise_cut[80,133] = 0
+    noise_cut[133,80] = 0
+    r = ClusterFileReader(fname)
+    cl = r.read(85, noise_cut)
+    assert cl.size == 10
+
+def test_chunk_config(data_path):
+    fname= data_path/'noise_test.clust'
+    #File contains total 70 clusters
+    r = ClusterFileReader(fname, chunk = 5)
+    assert r.read().size == 5
+    assert r.read(10).size == 10

tests/test_RawFileReader.py

@@ -1,31 +1,78 @@
 import pytest
 import os, sys
-from creader import RawFileReader
+from creader import RawFileReader, DetectorType
 from fixtures import data_path
 import numpy as np
 
+m03 = DetectorType.MOENCH_03
+m04a= DetectorType.MOENCH_04_A
+m04ad= DetectorType.MOENCH_04_AD
+
 def test_references_on_read(data_path):
-    fname= data_path/'test_d0_f0_0.raw'
+    fname= data_path/'Moench03_d0_f0_0.raw'
-    r = RawFileReader(fname)
+    r = RawFileReader(fname, m03)
     frames = r.read(10)
     assert sys.getrefcount(frames) == 2 #Over counts by one due to call by reference
 
 def test_references_on_read_with_header(data_path):
-    fname= data_path/'test_d0_f0_0.raw'
+    fname= data_path/'Moench03_d0_f0_0.raw'
-    r = RawFileReader(fname, header = True)
+    r = RawFileReader(fname, m03, header = True)
     frames, header = r.read(100)
     assert sys.getrefcount(frames) == 2 #Over counts by one due to call by reference
     assert sys.getrefcount(header) == 2 
 
 def test_reading_frame_numbers(data_path):
-    fname= data_path/'test_d0_f0_0.raw'
+    fname= data_path/'Moench03_d0_f0_0.raw'
-    r = RawFileReader(fname, header = True)
+    r = RawFileReader(fname, header = True, detector_type = m03)
     frames, header = r.read(1000)
     assert (header['Frame Number'] == np.arange(201,1201, dtype = np.uint64)).all()
 
 def test_reading_more_files_than_available(data_path):
-    fname= data_path/'test_d0_f0_0.raw'
+    fname= data_path/'Moench03_d0_f0_0.raw'
-    r = RawFileReader(fname, header = True)
+    r = RawFileReader(fname, m03, header = True)
     frames, header = r.read(1500)
     assert frames.shape == (1000,400, 400)
     assert header.size == 1000
+
+def test_references_on_m04_a(data_path):
+    fname= data_path/'Moench04_no_digital_d0_f0_0.raw'
+    r = RawFileReader(fname, m04a)
+    frames = r.read(3)
+    assert sys.getrefcount(frames) == 2 #Over counts by one due to call by reference
+    assert frames.shape == (3,400,400)
+
+def test_references_on_m04_a_with_heaer(data_path):
+    fname= data_path/'Moench04_no_digital_d0_f0_0.raw'
+    r = RawFileReader(fname, m04a, header = True)
+    frames, header = r.read(7)
+    assert sys.getrefcount(frames) == 2 #Over counts by one due to call by reference
+    assert frames.shape == (7,400,400)
+
+    assert sys.getrefcount(header) == 2 
+    assert header.size == 7
+
+def test_references_on_m04_ad(data_path):
+    fname= data_path/'Moench04_digital_d0_f0_0.raw'
+    r = RawFileReader(fname, m04ad)
+    analog, digital = r.read(3)
+    assert sys.getrefcount(analog) == 2 #Over counts by one due to call by reference
+    assert analog.shape == (3,400,400)
+    assert sys.getrefcount(digital) == 2
+
+def test_references_on_m04_ad_with_header(data_path):
+    fname= data_path/'Moench04_digital_d0_f0_0.raw'
+    r = RawFileReader(fname, m04ad, header = True)
+    analog, digital, header = r.read(3)
+    assert sys.getrefcount(analog) == 2 #Over counts by one due to call by reference
+    assert analog.shape == (3,400,400)
+    assert sys.getrefcount(digital) == 2
+    assert sys.getrefcount(header) == 2 
+
+def test_read_too_many_m04_frames(data_path):
+    fname= data_path/'Moench04_digital_d0_f0_0.raw'
+    r = RawFileReader(fname, m04ad, header = True)
+    analog, digital, header = r.read(150)
+    assert analog.shape == (100,400,400)
+    assert digital.shape == (100, 400, 400)
+    assert header.size == 100
+    assert sys.getrefcount(digital) == 2

utils/write_test_data.py

@@ -41,3 +41,20 @@ with open(path/'37frames_with_5_clusters.clust', 'wb') as f:
             print(data)
             data.tofile(f)
         
+#Writing out data to test noise cuts
+header[1] = 7
+
+with open(path/'noise_test.clust', 'wb') as f:
+    for i in range(10):
+        data['x'] = 50
+        data['y'] = 133
+        data['data'][:] = 50
+
+        header.tofile(f)
+        print(header)
+        header[0] += 1
+        
+        for j in range(7):
+            print(data)
+            data.tofile(f)
+            data['x'] += 10