added moench04 support

README.md

@@ -12,6 +12,7 @@ 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

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)
-        with open(fname) as f:
+                self.json = True
-            self.master = json.load(f)
+        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

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,12 +43,13 @@ static int RawFileReader_init(RawFileReader *self, PyObject *args,
             return -1;
 
     PyBytes_AsStringAndSize(fname_bytes, &fname, &len);
+
 #ifdef CR_VERBOSE
-    printf("fname: %s\n read_header: %d\n", fname, self->read_header);
+    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;
 
@@ -77,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
@@ -114,20 +136,35 @@ 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 we also read header we need to reshape the header
-        if(self->read_header){
+        if (self->read_header) {
             new_shape.len = 1;
-            PyArray_Resize((PyArrayObject *)header, &new_shape, 1, NPY_ANYORDER);
+            PyArray_Resize((PyArrayObject *)header, &new_shape, 1,
+                           NPY_ANYORDER);
         }
+    }
 
+    // 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);
     }
-    PyObject *ret = frames;
+
-    if(self->read_header){
+    if (self->read_header){
-        ret = PyTuple_Pack(2, frames, header);
+        Py_ssize_t old_size = PyTuple_GET_SIZE(ret);
-        Py_DECREF(header);
+        _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;
     }
-   return ret;
     
 }
 

src/data_types.h

@@ -34,7 +34,12 @@ typedef struct {
     uint32_t c;
 } 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
 {
@@ -55,4 +60,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,9 +18,9 @@ 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));
+            memcpy(header_out, &h, sizeof(Header));
-          header_out++;
+            header_out++;
         }
 
         // Read frame to temporary buffer
@@ -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,29 +39,123 @@ 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;
-      // analog part
+    char *tmp = malloc(frame_size);
-      int adc_nr = adc_numbers[i_sc];
+    char *digital_tmp = malloc(digital_frame_size);
-      int col = ((adc_nr * 25) + (n_pixel % 25));
+    Header h;
-      int row = 0;
+
-      if (i_sc / 4 % 2 == 0)
+    int64_t frames_read = 0;
-      {
+#ifdef CR_VERBOSE
-        row = 199 - (n_pixel / 25);
+    printf("MOENCH04: digital_out: %p\n", digital_out);
-      }
+#endif
-      else
+    while (frames_read < n_frames) {
-      {
+        // Read header to temp buffer, return on fail
-        row = 200 + (n_pixel / 25);
+        if (fread(&h, sizeof(Header), 1, fp) != 1) {
-      }
+            break;
-      int i_analog = n_pixel * 32 + i_sc;
+        }
-      out_buf[row * 400 + col] = buf[i_analog] & 0x3FFF;
+        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) {
+                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;
+        }
+    }
+}
+
+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_RawFileReader.py

@@ -1,31 +1,69 @@
 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