Developer (#164)

- State before merging the new cluster vector API

---------

Co-authored-by: Patrick <patrick.sieberer@psi.ch>
Co-authored-by: JulianHeymes <julian.heymes@psi.ch>
Co-authored-by: Dhanya Thattil <dhanya.thattil@psi.ch>
Co-authored-by: Xiangyu Xie <45243914+xiangyuxie@users.noreply.github.com>
Co-authored-by: xiangyu.xie <xiangyu.xie@psi.ch>
Co-authored-by: siebsi <sieb.patr@gmail.com>

.gitea/workflows/cmake_build.yml

@@ -2,9 +2,8 @@ name: Build the package using cmake then documentation
 
 on:
   workflow_dispatch:
-  push:
-    
 
+    
 
 permissions:
   contents: read
@@ -16,12 +15,12 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        platform: [ubuntu-latest, ] # macos-12, windows-2019]
-        python-version: ["3.12",]
+        platform: [ubuntu-latest, ] 
+        python-version: ["3.12", ] 
 
     runs-on: ${{ matrix.platform }}
 
-    # The setup-miniconda action needs this to activate miniconda
+
     defaults:
       run:
         shell: "bash -l {0}"
@@ -35,13 +34,13 @@ jobs:
             sudo apt-get -y install cmake gcc g++
 
       - name: Get conda
-        uses: conda-incubator/setup-miniconda@v3.0.4
+        uses: conda-incubator/setup-miniconda@v3
         with:
           python-version: ${{ matrix.python-version }}
+          environment-file: etc/dev-env.yml
+          miniforge-version: latest
           channels: conda-forge
-
-      - name: Prepare
-        run: conda install doxygen sphinx=7.1.2 breathe pybind11 sphinx_rtd_theme furo nlohmann_json zeromq fmt numpy
+          conda-remove-defaults: "true"
 
       - name: Build library 
         run: |
@@ -56,3 +55,4 @@ jobs:
 
 
 
+

.gitea/workflows/rh8-native.yml

@@ -0,0 +1,36 @@
+name: Build on RHEL8
+
+on:
+  push:
+  workflow_dispatch:
+
+permissions:
+  contents: read
+
+jobs:
+  build:
+    runs-on: "ubuntu-latest"
+    container:
+      image: gitea.psi.ch/images/rhel8-developer-gitea-actions
+    steps:
+      # workaround until actions/checkout@v4 is available for RH8
+      # - uses: actions/checkout@v4
+      - name: Clone repository
+        run: |
+          echo Cloning ${{ github.ref_name }}
+          git clone https://${{secrets.GITHUB_TOKEN}}@gitea.psi.ch/${{ github.repository }}.git --branch=${{ github.ref_name }} .
+
+
+      - name: Install dependencies
+        run: |
+          dnf install -y cmake python3.12 python3.12-devel python3.12-pip
+
+      - name: Build library 
+        run: |
+            mkdir build && cd build
+            cmake .. -DAARE_PYTHON_BINDINGS=ON -DAARE_TESTS=ON -DPython_FIND_VIRTUALENV=FIRST
+            make -j 2 
+
+      - name: C++ unit tests
+        working-directory: ${{gitea.workspace}}/build
+        run: ctest

.gitea/workflows/rh9-native.yml

@@ -0,0 +1,31 @@
+name: Build on RHEL9
+
+on:
+  push:
+  workflow_dispatch:
+
+permissions:
+  contents: read
+
+jobs:
+  build:
+    runs-on: "ubuntu-latest"
+    container:
+      image: gitea.psi.ch/images/rhel9-developer-gitea-actions
+    steps:
+      - uses: actions/checkout@v4
+
+
+      - name: Install dependencies
+        run: |
+          dnf install -y cmake python3.12 python3.12-devel python3.12-pip
+
+      - name: Build library 
+        run: |
+            mkdir build && cd build
+            cmake .. -DAARE_PYTHON_BINDINGS=ON -DAARE_TESTS=ON
+            make -j 2 
+
+      - name: C++ unit tests
+        working-directory: ${{gitea.workspace}}/build
+        run: ctest

.github/workflows/build_docs.yml

@@ -5,7 +5,6 @@ on:
   push:
     
 
-
 permissions:
   contents: read
   pages: write
@@ -16,12 +15,11 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        platform: [ubuntu-latest, ] # macos-12, windows-2019]
+        platform: [ubuntu-latest, ] 
         python-version: ["3.12",]
 
     runs-on: ${{ matrix.platform }}
 
-    # The setup-miniconda action needs this to activate miniconda
     defaults:
       run:
         shell: "bash -l {0}"
@@ -30,13 +28,13 @@ jobs:
       - uses: actions/checkout@v4
 
       - name: Get conda
-        uses: conda-incubator/setup-miniconda@v3.0.4
+        uses: conda-incubator/setup-miniconda@v3
         with:
           python-version: ${{ matrix.python-version }}
+          environment-file: etc/dev-env.yml
+          miniforge-version: latest
           channels: conda-forge
-
-      - name: Prepare
-        run: conda install doxygen sphinx=7.1.2 breathe pybind11 sphinx_rtd_theme furo nlohmann_json zeromq fmt numpy
+          conda-remove-defaults: "true"
 
       - name: Build library 
         run: |

.github/workflows/build_wheel.yml

@@ -0,0 +1,64 @@
+name: Build wheel
+
+on:
+    workflow_dispatch:
+    pull_request:
+    push:
+      branches:
+        - main
+    release:
+      types:
+        - published
+  
+
+jobs:
+  build_wheels:
+    name: Build wheels on ${{ matrix.os }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        os: [ubuntu-latest,]
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Build wheels
+        run: pipx run cibuildwheel==2.23.0
+
+      - uses: actions/upload-artifact@v4
+        with:
+          name: cibw-wheels-${{ matrix.os }}-${{ strategy.job-index }}
+          path: ./wheelhouse/*.whl
+    
+  build_sdist:
+    name: Build source distribution
+    runs-on: ubuntu-latest
+    steps:
+    - uses: actions/checkout@v4
+
+    - name: Build sdist
+      run: pipx run build --sdist
+
+    - uses: actions/upload-artifact@v4
+      with:
+        name: cibw-sdist
+        path: dist/*.tar.gz
+
+  upload_pypi:
+    needs: [build_wheels, build_sdist]
+    runs-on: ubuntu-latest
+    environment: pypi
+    permissions:
+        id-token: write
+    if: github.event_name == 'release' && github.event.action == 'published'
+    # or, alternatively, upload to PyPI on every tag starting with 'v' (remove on: release above to use this)
+    # if: github.event_name == 'push' && startsWith(github.ref, 'refs/tags/v')
+    steps:
+        - uses: actions/download-artifact@v4
+          with:
+            # unpacks all CIBW artifacts into dist/
+            pattern: cibw-*
+            path: dist
+            merge-multiple: true
+
+        - uses: pypa/gh-action-pypi-publish@release/v1

.gitignore

@@ -17,7 +17,8 @@ Testing/
 ctbDict.cpp
 ctbDict.h
 
-
+wheelhouse/
+dist/
 
 *.pyc
 */__pycache__/*

CMakeLists.txt

@@ -1,4 +1,4 @@
-cmake_minimum_required(VERSION 3.14)
+cmake_minimum_required(VERSION 3.15)
 
 project(aare 
     VERSION 1.0.0
@@ -11,6 +11,14 @@ set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_CXX_EXTENSIONS OFF)
 
+execute_process(
+        COMMAND git log -1 --format=%h
+        WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}
+        OUTPUT_VARIABLE GIT_HASH
+        OUTPUT_STRIP_TRAILING_WHITESPACE
+        )
+message(STATUS "Building from git hash: ${GIT_HASH}")
+
 if (${CMAKE_VERSION} VERSION_GREATER "3.24")
     cmake_policy(SET CMP0135 NEW) #Fetch content download timestamp
 endif()
@@ -31,7 +39,7 @@ set(CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH})
 
 
 # General options
-option(AARE_PYTHON_BINDINGS "Build python bindings" OFF)
+option(AARE_PYTHON_BINDINGS "Build python bindings" ON)
 option(AARE_TESTS "Build tests" OFF)
 option(AARE_BENCHMARKS "Build benchmarks" OFF)
 option(AARE_EXAMPLES "Build examples" OFF)
@@ -104,6 +112,7 @@ if(AARE_FETCH_LMFIT)
         )
     endif()
     
+
     #Disable what we don't need from lmfit
     set(BUILD_TESTING OFF CACHE BOOL "")
     set(LMFIT_CPPTEST OFF CACHE BOOL "")
@@ -331,8 +340,6 @@ endif()
 
 set(PUBLICHEADERS
     include/aare/ArrayExpr.hpp
-    include/aare/CalculateEta.hpp
-    include/aare/Cluster.hpp
     include/aare/ClusterFinder.hpp
     include/aare/ClusterFile.hpp
     include/aare/CtbRawFile.hpp
@@ -343,9 +350,10 @@ set(PUBLICHEADERS
     include/aare/File.hpp
     include/aare/Fit.hpp
     include/aare/FileInterface.hpp
+    include/aare/FilePtr.hpp
     include/aare/Frame.hpp
-    include/aare/GainMap.hpp
     include/aare/geo_helpers.hpp
+    include/aare/JungfrauDataFile.hpp
     include/aare/NDArray.hpp
     include/aare/NDView.hpp
     include/aare/NumpyFile.hpp
@@ -357,18 +365,22 @@ set(PUBLICHEADERS
     include/aare/RawSubFile.hpp
     include/aare/VarClusterFinder.hpp
     include/aare/utils/task.hpp
+
 )
 
 
 set(SourceFiles
     ${CMAKE_CURRENT_SOURCE_DIR}/src/CtbRawFile.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFile.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/defs.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/decode.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/File.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/FilePtr.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/Fit.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/geo_helpers.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/JungfrauDataFile.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/Interpolator.cpp
@@ -376,16 +388,20 @@ set(SourceFiles
     ${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/RawSubFile.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/RawMasterFile.cpp
+
     ${CMAKE_CURRENT_SOURCE_DIR}/src/utils/task.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/utils/ifstream_helpers.cpp
 )
 
 
 add_library(aare_core STATIC ${SourceFiles})
 target_include_directories(aare_core PUBLIC 
     "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>" 
-    "$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>" 
+    "$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>"
 )
 
+
+
 target_link_libraries(
     aare_core 
     PUBLIC 
@@ -394,7 +410,7 @@ target_link_libraries(
     ${STD_FS_LIB} # from helpers.cmake
     PRIVATE 
     aare_compiler_flags 
-    $<BUILD_INTERFACE:lmfit>
+    "$<BUILD_INTERFACE:lmfit>"
 
 )
 
@@ -411,6 +427,7 @@ if(AARE_TESTS)
     set(TestSources
             ${CMAKE_CURRENT_SOURCE_DIR}/src/algorithm.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/defs.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/decode.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/geo_helpers.test.cpp
@@ -419,10 +436,9 @@ if(AARE_TESTS)
             ${CMAKE_CURRENT_SOURCE_DIR}/src/NDView.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFinder.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterVector.test.cpp
-            ${CMAKE_CURRENT_SOURCE_DIR}/src/Cluster.test.cpp
-            ${CMAKE_CURRENT_SOURCE_DIR}/src/CalculateEta.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFile.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/Pedestal.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/JungfrauDataFile.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.test.cpp

benchmarks/CMakeLists.txt

@@ -1,27 +1,11 @@
+find_package(benchmark REQUIRED)
 
-include(FetchContent)
+add_executable(ndarray_benchmark ndarray_benchmark.cpp)
 
+target_link_libraries(ndarray_benchmark benchmark::benchmark aare_core aare_compiler_flags)
+# target_link_libraries(tests PRIVATE aare_core aare_compiler_flags)
 
-FetchContent_Declare(
-    benchmark
-    GIT_REPOSITORY https://github.com/google/benchmark.git
-    GIT_TAG        v1.8.3 # Change to the latest version if needed
-)
-
-# Ensure Google Benchmark is built correctly
-set(BENCHMARK_ENABLE_TESTING OFF CACHE BOOL "" FORCE)
-
-FetchContent_MakeAvailable(benchmark)
-
-add_executable(benchmarks)
-
-target_sources(benchmarks PRIVATE ndarray_benchmark.cpp calculateeta_benchmark.cpp)
-
-# Link Google Benchmark and other necessary libraries
-target_link_libraries(benchmarks PRIVATE benchmark::benchmark aare_core aare_compiler_flags)
-
-# Set output properties
-set_target_properties(benchmarks PROPERTIES
-    RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}
-    OUTPUT_NAME run_benchmarks
+set_target_properties(ndarray_benchmark PROPERTIES
+        RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}
+        # OUTPUT_NAME run_tests
 )

benchmarks/calculateeta_benchmark.cpp

@@ -1,70 +0,0 @@
-#include "aare/CalculateEta.hpp"
-#include "aare/ClusterFile.hpp"
-#include <benchmark/benchmark.h>
-
-using namespace aare;
-
-class ClusterFixture : public benchmark::Fixture {
-  public:
-    Cluster<int, 2, 2> cluster_2x2{};
-    Cluster<int, 3, 3> cluster_3x3{};
-
-  private:
-    using benchmark::Fixture::SetUp;
-
-    void SetUp([[maybe_unused]] const benchmark::State &state) override {
-        int temp_data[4] = {1, 2, 3, 1};
-        std::copy(std::begin(temp_data), std::end(temp_data),
-                  std::begin(cluster_2x2.data));
-
-        cluster_2x2.x = 0;
-        cluster_2x2.y = 0;
-
-        int temp_data2[9] = {1, 2, 3, 1, 3, 4, 5, 1, 20};
-        std::copy(std::begin(temp_data2), std::end(temp_data2),
-                  std::begin(cluster_3x3.data));
-
-        cluster_3x3.x = 0;
-        cluster_3x3.y = 0;
-    }
-
-    // void TearDown(::benchmark::State& state) {
-    // }
-};
-
-BENCHMARK_F(ClusterFixture, Calculate2x2Eta)(benchmark::State &st) {
-    for (auto _ : st) {
-        // This code gets timed
-        Eta2 eta = calculate_eta2(cluster_2x2);
-        benchmark::DoNotOptimize(eta);
-    }
-}
-
-// almost takes double the time
-BENCHMARK_F(ClusterFixture,
-            CalculateGeneralEtaFor2x2Cluster)(benchmark::State &st) {
-    for (auto _ : st) {
-        // This code gets timed
-        Eta2 eta = calculate_eta2<int, 2, 2>(cluster_2x2);
-        benchmark::DoNotOptimize(eta);
-    }
-}
-
-BENCHMARK_F(ClusterFixture, Calculate3x3Eta)(benchmark::State &st) {
-    for (auto _ : st) {
-        // This code gets timed
-        Eta2 eta = calculate_eta2(cluster_3x3);
-        benchmark::DoNotOptimize(eta);
-    }
-}
-
-// almost takes double the time
-BENCHMARK_F(ClusterFixture,
-            CalculateGeneralEtaFor3x3Cluster)(benchmark::State &st) {
-    for (auto _ : st) {
-        // This code gets timed
-        Eta2 eta = calculate_eta2<int, 3, 3>(cluster_3x3);
-        benchmark::DoNotOptimize(eta);
-    }
-}
-// BENCHMARK_MAIN();

conda-recipe/meta.yaml

@@ -1,6 +1,8 @@
 package:
   name: aare
-  version: 2025.4.1 #TODO! how to not duplicate this?
+  version: 2025.4.22 #TODO! how to not duplicate this?
+
+
 
 
 
@@ -37,6 +39,7 @@ requirements:
   run:
     - python {{python}}
     - numpy {{ numpy }}
+    - matplotlib
 
 
 test:

docs/src/JungfrauDataFile.rst

@@ -0,0 +1,25 @@
+JungfrauDataFile
+==================
+
+JungfrauDataFile is a class to read the .dat files that are produced by Aldo's receiver. 
+It is mostly used for calibration. 
+
+The structure of the file is:
+
+* JungfrauDataHeader
+* Binary data (256x256, 256x1024 or 512x1024)
+* JungfrauDataHeader
+* ...
+
+There is no metadata indicating number of frames or the size of the image, but this
+will be infered by this reader. 
+
+.. doxygenstruct:: aare::JungfrauDataHeader
+   :members:
+   :undoc-members:
+   :private-members:
+
+.. doxygenclass:: aare::JungfrauDataFile
+   :members:
+   :undoc-members:
+   :private-members:

docs/src/Tests.rst

@@ -0,0 +1,47 @@
+****************
+Tests
+****************
+
+We test the code both from the C++ and Python API. By default only tests that does not require image data is run. 
+
+C++
+~~~~~~~~~~~~~~~~~~
+
+.. code-block:: bash
+
+    mkdir build
+    cd build
+    cmake .. -DAARE_TESTS=ON
+    make -j 4
+
+    export AARE_TEST_DATA=/path/to/test/data
+    ./run_test [.files] #or using ctest, [.files] is the option to include tests needing data
+
+
+
+Python
+~~~~~~~~~~~~~~~~~~
+
+.. code-block:: bash
+
+    #From the root dir of the library
+    python -m pytest python/tests --files # passing --files will run the tests needing data
+
+
+
+Getting the test data
+~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. attention ::
+
+    The tests needing the test data are not run by default. To make the data available, you need to set the environment variable
+    AARE_TEST_DATA to the path of the test data directory. Then pass either [.files] for the C++ tests or --files for Python
+
+The image files needed for the test are large and are not included in the repository. They are stored
+using GIT LFS in a separate repository. To get the test data, you need to clone the repository.
+To do this, you need to have GIT LFS installed. You can find instructions on how to install it here: https://git-lfs.github.com/
+Once you have GIT LFS installed, you can clone the repository like any normal repo using:
+
+.. code-block:: bash
+    
+    git clone https://gitea.psi.ch/detectors/aare-test-data.git

docs/src/algorithm.rst

@@ -0,0 +1,5 @@
+algorithm
+=============
+
+.. doxygenfile:: algorithm.hpp
+

docs/src/index.rst

@@ -20,9 +20,6 @@ AARE
     Requirements
     Consume
 
-
-
-
 .. toctree::
     :caption: Python API
     :maxdepth: 1
@@ -31,6 +28,7 @@ AARE
     pyCtbRawFile
     pyClusterFile
     pyClusterVector
+    pyJungfrauDataFile
     pyRawFile
     pyRawMasterFile
     pyVarClusterFinder
@@ -42,6 +40,7 @@ AARE
     :caption: C++ API
     :maxdepth: 1
     
+    algorithm
     NDArray
     NDView
     Frame
@@ -51,6 +50,7 @@ AARE
     ClusterFinderMT
     ClusterFile
     ClusterVector
+    JungfrauDataFile
     Pedestal
     RawFile
     RawSubFile
@@ -59,4 +59,8 @@ AARE
 
 
 
-
+.. toctree::
+    :caption: Developer
+    :maxdepth: 3
+    
+    Tests

docs/src/pyJungfrauDataFile.rst

@@ -0,0 +1,10 @@
+JungfrauDataFile
+===================
+
+.. py:currentmodule:: aare
+
+.. autoclass:: JungfrauDataFile
+    :members:
+    :undoc-members:
+    :show-inheritance:
+    :inherited-members:

etc/dev-env.yml

@@ -0,0 +1,15 @@
+name: dev-environment
+channels:
+  - conda-forge
+dependencies:
+  - anaconda-client
+  - doxygen 
+  - sphinx=7.1.2 
+  - breathe 
+  - pybind11 
+  - sphinx_rtd_theme 
+  - furo 
+  - nlohmann_json 
+  - zeromq 
+  - fmt 
+  - numpy

include/aare/CalculateEta.hpp

@@ -1,122 +0,0 @@
-#pragma once
-
-#include "aare/Cluster.hpp"
-#include "aare/ClusterVector.hpp"
-#include "aare/NDArray.hpp"
-
-namespace aare {
-
-typedef enum {
-    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
-} pixel;
-
-template <typename T> struct Eta2 {
-    double x;
-    double y;
-    int c;
-    T sum;
-};
-
-/**
- * @brief Calculate the eta2 values for all clusters in a Clsutervector
- */
-template <typename ClusterType,
-          typename = std::enable_if_t<is_cluster_v<ClusterType>>>
-NDArray<double, 2> calculate_eta2(const ClusterVector<ClusterType> &clusters) {
-    NDArray<double, 2> eta2({static_cast<int64_t>(clusters.size()), 2});
-
-    for (size_t i = 0; i < clusters.size(); i++) {
-        auto e = calculate_eta2(clusters.at(i));
-        eta2(i, 0) = e.x;
-        eta2(i, 1) = e.y;
-    }
-
-    return eta2;
-}
-
-/**
- * @brief Calculate the eta2 values for a generic sized cluster and return them
- * in a Eta2 struct containing etay, etax and the index of the respective 2x2
- * subcluster.
- */
-template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType>
-Eta2<T>
-calculate_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
-    Eta2<T> eta{};
-
-    auto max_sum = cl.max_sum_2x2();
-    eta.sum = max_sum.first;
-    auto c = max_sum.second;
-
-    size_t index_bottom_left_max_2x2_subcluster =
-        (int(c / (ClusterSizeX - 1))) * ClusterSizeX + c % (ClusterSizeX - 1);
-
-    if ((cl.data[index_bottom_left_max_2x2_subcluster] +
-         cl.data[index_bottom_left_max_2x2_subcluster + 1]) != 0)
-        eta.x = static_cast<double>(
-                    cl.data[index_bottom_left_max_2x2_subcluster + 1]) /
-                static_cast<double>(
-                    (cl.data[index_bottom_left_max_2x2_subcluster] +
-                     cl.data[index_bottom_left_max_2x2_subcluster + 1]));
-
-    if ((cl.data[index_bottom_left_max_2x2_subcluster] +
-         cl.data[index_bottom_left_max_2x2_subcluster + ClusterSizeX]) != 0)
-        eta.y =
-            static_cast<double>(
-                cl.data[index_bottom_left_max_2x2_subcluster + ClusterSizeX]) /
-            static_cast<double>(
-                (cl.data[index_bottom_left_max_2x2_subcluster] +
-                 cl.data[index_bottom_left_max_2x2_subcluster + ClusterSizeX]));
-
-    eta.c = c; // TODO only supported for 2x2 and 3x3 clusters -> at least no
-               // underyling enum class
-    return eta;
-}
-
-// calculates Eta3 for 3x3 cluster based on code from analyze_cluster
-// TODO only supported for 3x3 Clusters
-template <typename T> Eta2<T> calculate_eta3(const Cluster<T, 3, 3> &cl) {
-
-    Eta2<T> eta{};
-
-    T sum = 0;
-
-    std::for_each(std::begin(cl.data), std::end(cl.data),
-                  [&sum](T x) { sum += x; });
-
-    eta.sum = sum;
-
-    eta.c = corner::cBottomLeft;
-
-    if ((cl.data[3] + cl.data[4] + cl.data[5]) != 0)
-
-        eta.x = static_cast<double>(-cl.data[3] + cl.data[3 + 2]) /
-
-                (cl.data[3] + cl.data[4] + cl.data[5]);
-
-    if ((cl.data[1] + cl.data[4] + cl.data[7]) != 0)
-
-        eta.y = static_cast<double>(-cl.data[1] + cl.data[2 * 3 + 1]) /
-
-                (cl.data[1] + cl.data[4] + cl.data[7]);
-
-    return eta;
-}
-
-} // namespace aare

include/aare/Cluster.hpp

@@ -1,121 +1,36 @@
-
-/************************************************
- * @file Cluster.hpp
- * @short definition of cluster, where CoordType (x,y) give
- * the cluster center coordinates and data the actual cluster data
- * cluster size is given as template parameters
- ***********************************************/
-
 #pragma once
 
 #include <algorithm>
 #include <array>
+#include <cstddef>
 #include <cstdint>
 #include <numeric>
-#include <type_traits>
 
 namespace aare {
 
-template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType = int16_t>
-constexpr bool is_valid_cluster =
-    std::is_arithmetic_v<T> && std::is_integral_v<CoordType> &&
-    (ClusterSizeX > 0) && (ClusterSizeY > 0);
-
-// requires clause c++20 maybe update
-template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType = int16_t,
-          typename Enable = std::enable_if_t<
-              is_valid_cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>>
-struct Cluster {
-    CoordType x;
-    CoordType y;
-    T data[ClusterSizeX * ClusterSizeY];
-
-    T sum() const {
-        return std::accumulate(data, data + ClusterSizeX * ClusterSizeY, 0);
-    }
-
-    std::pair<T, int> max_sum_2x2() const {
-
-        constexpr size_t num_2x2_subclusters =
-            (ClusterSizeX - 1) * (ClusterSizeY - 1);
-
-        std::array<T, num_2x2_subclusters> sum_2x2_subcluster;
-        for (size_t i = 0; i < ClusterSizeY - 1; ++i) {
-            for (size_t j = 0; j < ClusterSizeX - 1; ++j)
-                sum_2x2_subcluster[i * (ClusterSizeX - 1) + j] =
-                    data[i * ClusterSizeX + j] +
-                    data[i * ClusterSizeX + j + 1] +
-                    data[(i + 1) * ClusterSizeX + j] +
-                    data[(i + 1) * ClusterSizeX + j + 1];
-        }
-
-        int index = std::max_element(sum_2x2_subcluster.begin(),
-                                     sum_2x2_subcluster.end()) -
-                    sum_2x2_subcluster.begin();
-        return std::make_pair(sum_2x2_subcluster[index], index);
-    }
-};
-
-// Specialization for 2x2 clusters (only one sum exists)
-template <typename T> struct Cluster<T, 2, 2, int16_t> {
+//TODO! Template this? 
+struct Cluster3x3 {
     int16_t x;
     int16_t y;
-    T data[4];
+    int32_t data[9];
 
-    T sum() const { return std::accumulate(data, data + 4, 0); }
+    int32_t sum_2x2() const{
+        std::array<int32_t, 4> total;
+        total[0] = data[0] + data[1] + data[3] + data[4];
+        total[1] = data[1] + data[2] + data[4] + data[5];
+        total[2] = data[3] + data[4] + data[6] + data[7];
+        total[3] = data[4] + data[5] + data[7] + data[8];
+        return *std::max_element(total.begin(), total.end());
+    }
 
-    std::pair<T, int> max_sum_2x2() const {
-        return std::make_pair(data[0] + data[1] + data[2] + data[3],
-                              0); // Only one possible 2x2 sum
+    int32_t sum() const{
+        return std::accumulate(data, data + 9, 0);
     }
 };
-
-// Specialization for 3x3 clusters
-template <typename T> struct Cluster<T, 3, 3, int16_t> {
+struct Cluster2x2 {
     int16_t x;
     int16_t y;
-    T data[9];
-
-    T sum() const { return std::accumulate(data, data + 9, 0); }
-
-    std::pair<T, int> max_sum_2x2() const {
-        std::array<T, 4> sum_2x2_subclusters;
-        sum_2x2_subclusters[0] = data[0] + data[1] + data[3] + data[4];
-        sum_2x2_subclusters[1] = data[1] + data[2] + data[4] + data[5];
-        sum_2x2_subclusters[2] = data[3] + data[4] + data[6] + data[7];
-        sum_2x2_subclusters[3] = data[4] + data[5] + data[7] + data[8];
-        int index = std::max_element(sum_2x2_subclusters.begin(),
-                                     sum_2x2_subclusters.end()) -
-                    sum_2x2_subclusters.begin();
-        return std::make_pair(sum_2x2_subclusters[index], index);
-    }
+    int32_t data[4];
 };
 
-// Type Traits for is_cluster_type
-template <typename T>
-struct is_cluster : std::false_type {}; // Default case: Not a Cluster
-
-template <typename T, uint8_t X, uint8_t Y, typename CoordType>
-struct is_cluster<Cluster<T, X, Y, CoordType>> : std::true_type {}; // Cluster
-
-template <typename T> constexpr bool is_cluster_v = is_cluster<T>::value;
-
-template <typename ClusterType,
-          typename = std::enable_if_t<is_cluster_v<ClusterType>>>
-struct extract_template_arguments; // Forward declaration
-
-// helper struct to extract template argument
-template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType>
-struct extract_template_arguments<
-    Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
-
-    using value_type = T;
-    static constexpr int cluster_size_x = ClusterSizeX;
-    static constexpr int cluster_size_y = ClusterSizeY;
-    using coordtype = CoordType;
-};
-
-} // namespace aare
+} // namespace aare

include/aare/ClusterCollector.hpp

@@ -2,31 +2,29 @@
 #include <atomic>
 #include <thread>
 
-#include "aare/ClusterFinderMT.hpp"
-#include "aare/ClusterVector.hpp"
 #include "aare/ProducerConsumerQueue.hpp"
+#include "aare/ClusterVector.hpp"
+#include "aare/ClusterFinderMT.hpp"
 
 namespace aare {
 
-template <typename ClusterType,
-          typename = std::enable_if_t<is_cluster_v<ClusterType>>>
-class ClusterCollector {
-    ProducerConsumerQueue<ClusterVector<ClusterType>> *m_source;
-    std::atomic<bool> m_stop_requested{false};
-    std::atomic<bool> m_stopped{true};
-    std::chrono::milliseconds m_default_wait{1};
-    std::thread m_thread;
-    std::vector<ClusterVector<ClusterType>> m_clusters;
+class ClusterCollector{
+        ProducerConsumerQueue<ClusterVector<int>>* m_source;
+        std::atomic<bool> m_stop_requested{false};
+        std::atomic<bool> m_stopped{true};
+        std::chrono::milliseconds m_default_wait{1};
+        std::thread m_thread;
+        std::vector<ClusterVector<int>> m_clusters;
 
-    void process() {
+    void process(){
         m_stopped = false;
         fmt::print("ClusterCollector started\n");
-        while (!m_stop_requested || !m_source->isEmpty()) {
-            if (ClusterVector<ClusterType> *clusters = m_source->frontPtr();
+        while (!m_stop_requested || !m_source->isEmpty()) { 
+            if (ClusterVector<int> *clusters = m_source->frontPtr();
                 clusters != nullptr) {
                 m_clusters.push_back(std::move(*clusters));
                 m_source->popFront();
-            } else {
+            }else{
                 std::this_thread::sleep_for(m_default_wait);
             }
         }
@@ -34,21 +32,21 @@ class ClusterCollector {
         m_stopped = true;
     }
 
-  public:
-    ClusterCollector(ClusterFinderMT<ClusterType, uint16_t, double> *source) {
-        m_source = source->sink();
-        m_thread = std::thread(&ClusterCollector::process, this);
-    }
-    void stop() {
-        m_stop_requested = true;
-        m_thread.join();
-    }
-    std::vector<ClusterVector<ClusterType>> steal_clusters() {
-        if (!m_stopped) {
-            throw std::runtime_error("ClusterCollector is still running");
+    public:
+        ClusterCollector(ClusterFinderMT<uint16_t, double, int32_t>* source){
+            m_source = source->sink();
+            m_thread = std::thread(&ClusterCollector::process, this);
+        }
+        void stop(){
+            m_stop_requested = true;
+            m_thread.join();
+        }
+        std::vector<ClusterVector<int>> steal_clusters(){
+            if(!m_stopped){
+                throw std::runtime_error("ClusterCollector is still running");
+            }
+            return std::move(m_clusters);
         }
-        return std::move(m_clusters);
-    }
 };
 
 } // namespace aare

include/aare/ClusterFile.hpp

@@ -2,7 +2,6 @@
 
 #include "aare/Cluster.hpp"
 #include "aare/ClusterVector.hpp"
-#include "aare/GainMap.hpp"
 #include "aare/NDArray.hpp"
 #include "aare/defs.hpp"
 #include <filesystem>
@@ -11,18 +10,43 @@
 
 namespace aare {
 
-/*
-Binary cluster file. Expects data to be layed out as:
-int32_t frame_number
-uint32_t number_of_clusters
-int16_t x, int16_t y, int32_t data[9] x number_of_clusters
-int32_t frame_number
-uint32_t number_of_clusters
-....
-*/
 
-// TODO: change to support any type of clusters, e.g. header line with
-// clsuter_size_x, cluster_size_y,
+//TODO! Legacy enums, migrate to enum class
+typedef enum {
+    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
+} pixel;
+
+struct Eta2 {
+    double x;
+    double y;
+    corner c;
+    int32_t sum;
+};
+
+struct ClusterAnalysis {
+    uint32_t c;
+    int32_t tot;
+    double etax;
+    double etay;
+};
+
+
+
 /**
  * @brief Class to read and write cluster files
  * Expects data to be laid out as:
@@ -35,18 +59,14 @@ uint32_t number_of_clusters
  *       uint32_t number_of_clusters
  *       etc.
  */
-template <typename ClusterType,
-          typename Enable = std::enable_if_t<is_cluster_v<ClusterType>>>
 class ClusterFile {
     FILE *fp{};
-    uint32_t m_num_left{};    /*Number of photons left in frame*/
-    size_t m_chunk_size{};    /*Number of clusters to read at a time*/
-    const std::string m_mode; /*Mode to open the file in*/
-    std::optional<ROI> m_roi; /*Region of interest, will be applied if set*/
-    std::optional<NDArray<int32_t, 2>>
-        m_noise_map; /*Noise map to cut photons, will be applied if set*/
-    std::optional<GainMap> m_gain_map; /*Gain map to apply to the clusters, will
-                                          be applied if set*/
+    uint32_t m_num_left{};            /*Number of photons left in frame*/
+    size_t m_chunk_size{};            /*Number of clusters to read at a time*/
+    const std::string m_mode;         /*Mode to open the file in*/
+    std::optional<ROI> m_roi;         /*Region of interest, will be applied if set*/
+    std::optional<NDArray<int32_t, 2>> m_noise_map; /*Noise map to cut photons, will be applied if set*/
+    std::optional<NDArray<double, 2>> m_gain_map; /*Gain map to apply to the clusters, will be applied if set*/
 
   public:
     /**
@@ -60,7 +80,8 @@ class ClusterFile {
      */
     ClusterFile(const std::filesystem::path &fname, size_t chunk_size = 1000,
                 const std::string &mode = "r");
-
+    
+    
     ~ClusterFile();
 
     /**
@@ -68,388 +89,64 @@ class ClusterFile {
      * If EOF is reached the returned vector will have less than n_clusters
      * clusters
      */
-    ClusterVector<ClusterType> read_clusters(size_t n_clusters);
+    ClusterVector<int32_t> read_clusters(size_t n_clusters);
+
+    ClusterVector<int32_t> read_clusters(size_t n_clusters, ROI roi);
 
     /**
      * @brief Read a single frame from the file and return the clusters. The
      * cluster vector will have the frame number set.
-     * @throws std::runtime_error if the file is not opened for reading or the
-     * file pointer not at the beginning of a frame
+     * @throws std::runtime_error if the file is not opened for reading or the file pointer not
+     * at the beginning of a frame
      */
-    ClusterVector<ClusterType> read_frame();
+    ClusterVector<int32_t> read_frame();
 
-    void write_frame(const ClusterVector<ClusterType> &clusters);
 
+    void write_frame(const ClusterVector<int32_t> &clusters);
+    
     /**
      * @brief Return the chunk size
      */
     size_t chunk_size() const { return m_chunk_size; }
 
     /**
-     * @brief Set the region of interest to use when reading clusters. If set
-     * only clusters within the ROI will be read.
+     * @brief Set the region of interest to use when reading clusters. If set only clusters within
+     * the ROI will be read.
      */
     void set_roi(ROI roi);
 
     /**
-     * @brief Set the noise map to use when reading clusters. If set clusters
-     * below the noise level will be discarded. Selection criteria one of:
-     * Central pixel above noise, highest 2x2 sum above 2 * noise, total sum
-     * above 3 * noise.
+     * @brief Set the noise map to use when reading clusters. If set clusters below the noise
+     * level will be discarded. Selection criteria one of: Central pixel above noise, highest
+     * 2x2 sum above 2 * noise, total sum above 3 * noise.
      */
     void set_noise_map(const NDView<int32_t, 2> noise_map);
 
     /**
-     * @brief Set the gain map to use when reading clusters. If set the gain map
-     * will be applied to the clusters that pass ROI and noise_map selection.
+     * @brief Set the gain map to use when reading clusters. If set the gain map will be applied
+     * to the clusters that pass ROI and noise_map selection. The gain map is expected to be in ADU/energy.
      */
     void set_gain_map(const NDView<double, 2> gain_map);
-
-    void set_gain_map(const GainMap &gain_map);
-
-    void set_gain_map(const GainMap &&gain_map);
-
+    
+    
     /**
-     * @brief Close the file. If not closed the file will be closed in the
-     * destructor
+     * @brief Close the file. If not closed the file will be closed in the destructor
      */
     void close();
 
-  private:
-    ClusterVector<ClusterType> read_clusters_with_cut(size_t n_clusters);
-    ClusterVector<ClusterType> read_clusters_without_cut(size_t n_clusters);
-    ClusterVector<ClusterType> read_frame_with_cut();
-    ClusterVector<ClusterType> read_frame_without_cut();
-    bool is_selected(ClusterType &cl);
-    ClusterType read_one_cluster();
+    private:
+        ClusterVector<int32_t> read_clusters_with_cut(size_t n_clusters);
+        ClusterVector<int32_t> read_clusters_without_cut(size_t n_clusters);
+        ClusterVector<int32_t> read_frame_with_cut();
+        ClusterVector<int32_t> read_frame_without_cut();
+        bool is_selected(Cluster3x3 &cl);
+        Cluster3x3 read_one_cluster();
 };
 
-template <typename ClusterType, typename Enable>
-ClusterFile<ClusterType, Enable>::ClusterFile(
-    const std::filesystem::path &fname, size_t chunk_size,
-    const std::string &mode)
-    : m_chunk_size(chunk_size), m_mode(mode) {
+//TODO! helper functions that doesn't really belong here
+NDArray<double, 2> calculate_eta2(ClusterVector<int> &clusters);
+Eta2 calculate_eta2(Cluster3x3 &cl);
+Eta2 calculate_eta2(Cluster2x2 &cl);
 
-    if (mode == "r") {
-        fp = fopen(fname.c_str(), "rb");
-        if (!fp) {
-            throw std::runtime_error("Could not open file for reading: " +
-                                     fname.string());
-        }
-    } else if (mode == "w") {
-        fp = fopen(fname.c_str(), "wb");
-        if (!fp) {
-            throw std::runtime_error("Could not open file for writing: " +
-                                     fname.string());
-        }
-    } else if (mode == "a") {
-        fp = fopen(fname.c_str(), "ab");
-        if (!fp) {
-            throw std::runtime_error("Could not open file for appending: " +
-                                     fname.string());
-        }
-    } else {
-        throw std::runtime_error("Unsupported mode: " + mode);
-    }
-}
-
-template <typename ClusterType, typename Enable>
-ClusterFile<ClusterType, Enable>::~ClusterFile() {
-    close();
-}
-
-template <typename ClusterType, typename Enable>
-void ClusterFile<ClusterType, Enable>::close() {
-    if (fp) {
-        fclose(fp);
-        fp = nullptr;
-    }
-}
-template <typename ClusterType, typename Enable>
-void ClusterFile<ClusterType, Enable>::set_roi(ROI roi) {
-    m_roi = roi;
-}
-template <typename ClusterType, typename Enable>
-void ClusterFile<ClusterType, Enable>::set_noise_map(
-    const NDView<int32_t, 2> noise_map) {
-    m_noise_map = NDArray<int32_t, 2>(noise_map);
-}
-template <typename ClusterType, typename Enable>
-void ClusterFile<ClusterType, Enable>::set_gain_map(
-    const NDView<double, 2> gain_map) {
-    m_gain_map = GainMap(gain_map);
-}
-
-template <typename ClusterType, typename Enable>
-void ClusterFile<ClusterType, Enable>::set_gain_map(const GainMap &gain_map) {
-    m_gain_map = gain_map;
-}
-
-template <typename ClusterType, typename Enable>
-void ClusterFile<ClusterType, Enable>::set_gain_map(const GainMap &&gain_map) {
-    m_gain_map = gain_map;
-}
-
-// TODO generally supported for all clsuter types
-template <typename ClusterType, typename Enable>
-void ClusterFile<ClusterType, Enable>::write_frame(
-    const ClusterVector<ClusterType> &clusters) {
-    if (m_mode != "w" && m_mode != "a") {
-        throw std::runtime_error("File not opened for writing");
-    }
-    if (!(clusters.cluster_size_x() == 3) &&
-        !(clusters.cluster_size_y() == 3)) {
-        throw std::runtime_error("Only 3x3 clusters are supported");
-    }
-    int32_t frame_number = clusters.frame_number();
-    fwrite(&frame_number, sizeof(frame_number), 1, fp);
-    uint32_t n_clusters = clusters.size();
-    fwrite(&n_clusters, sizeof(n_clusters), 1, fp);
-    fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp);
-}
-
-template <typename ClusterType, typename Enable>
-ClusterVector<ClusterType>
-ClusterFile<ClusterType, Enable>::read_clusters(size_t n_clusters) {
-    if (m_mode != "r") {
-        throw std::runtime_error("File not opened for reading");
-    }
-    if (m_noise_map || m_roi) {
-        return read_clusters_with_cut(n_clusters);
-    } else {
-        return read_clusters_without_cut(n_clusters);
-    }
-}
-
-template <typename ClusterType, typename Enable>
-ClusterVector<ClusterType>
-ClusterFile<ClusterType, Enable>::read_clusters_without_cut(size_t n_clusters) {
-    if (m_mode != "r") {
-        throw std::runtime_error("File not opened for reading");
-    }
-
-    ClusterVector<ClusterType> clusters(n_clusters);
-
-    int32_t iframe = 0; // frame number needs to be 4 bytes!
-    size_t nph_read = 0;
-    uint32_t nn = m_num_left;
-    uint32_t nph = m_num_left; // number of clusters in frame needs to be 4
-
-    // auto buf = reinterpret_cast<Cluster3x3 *>(clusters.data());
-    auto buf = clusters.data();
-    // if there are photons left from previous frame read them first
-    if (nph) {
-        if (nph > n_clusters) {
-            // if we have more photons left in the frame then photons to read we
-            // read directly the requested number
-            nn = n_clusters;
-        } else {
-            nn = nph;
-        }
-        nph_read += fread((buf + nph_read * clusters.item_size()),
-                          clusters.item_size(), nn, fp);
-        m_num_left = nph - nn; // write back the number of photons left
-    }
-
-    if (nph_read < n_clusters) {
-        // keep on reading frames and photons until reaching n_clusters
-        while (fread(&iframe, sizeof(iframe), 1, fp)) {
-            clusters.set_frame_number(iframe);
-            // read number of clusters in frame
-            if (fread(&nph, sizeof(nph), 1, fp)) {
-                if (nph > (n_clusters - nph_read))
-                    nn = n_clusters - nph_read;
-                else
-                    nn = nph;
-
-                nph_read += fread((buf + nph_read * clusters.item_size()),
-                                  clusters.item_size(), nn, fp);
-                m_num_left = nph - nn;
-            }
-            if (nph_read >= n_clusters)
-                break;
-        }
-    }
-
-    // Resize the vector to the number of clusters.
-    // No new allocation, only change bounds.
-    clusters.resize(nph_read);
-    if (m_gain_map)
-        m_gain_map->apply_gain_map(clusters);
-    return clusters;
-}
-
-template <typename ClusterType, typename Enable>
-ClusterVector<ClusterType>
-ClusterFile<ClusterType, Enable>::read_clusters_with_cut(size_t n_clusters) {
-    ClusterVector<ClusterType> clusters;
-    clusters.reserve(n_clusters);
-
-    // if there are photons left from previous frame read them first
-    if (m_num_left) {
-        while (m_num_left && clusters.size() < n_clusters) {
-            ClusterType c = read_one_cluster();
-            if (is_selected(c)) {
-                clusters.push_back(c);
-            }
-        }
-    }
-
-    // we did not have enough clusters left in the previous frame
-    // keep on reading frames until reaching n_clusters
-    if (clusters.size() < n_clusters) {
-        // sanity check
-        if (m_num_left) {
-            throw std::runtime_error(
-                LOCATION + "Entered second loop with clusters left\n");
-        }
-
-        int32_t frame_number = 0; // frame number needs to be 4 bytes!
-        while (fread(&frame_number, sizeof(frame_number), 1, fp)) {
-            if (fread(&m_num_left, sizeof(m_num_left), 1, fp)) {
-                clusters.set_frame_number(
-                    frame_number); // cluster vector will hold the last frame
-                                   // number
-                while (m_num_left && clusters.size() < n_clusters) {
-                    ClusterType c = read_one_cluster();
-                    if (is_selected(c)) {
-                        clusters.push_back(c);
-                    }
-                }
-            }
-
-            // we have enough clusters, break out of the outer while loop
-            if (clusters.size() >= n_clusters)
-                break;
-        }
-    }
-    if (m_gain_map)
-        m_gain_map->apply_gain_map(clusters);
-
-    return clusters;
-}
-
-template <typename ClusterType, typename Enable>
-ClusterType ClusterFile<ClusterType, Enable>::read_one_cluster() {
-    ClusterType c;
-    auto rc = fread(&c, sizeof(c), 1, fp);
-    if (rc != 1) {
-        throw std::runtime_error(LOCATION + "Could not read cluster");
-    }
-    --m_num_left;
-    return c;
-}
-
-template <typename ClusterType, typename Enable>
-ClusterVector<ClusterType> ClusterFile<ClusterType, Enable>::read_frame() {
-    if (m_mode != "r") {
-        throw std::runtime_error(LOCATION + "File not opened for reading");
-    }
-    if (m_noise_map || m_roi) {
-        return read_frame_with_cut();
-    } else {
-        return read_frame_without_cut();
-    }
-}
-
-template <typename ClusterType, typename Enable>
-ClusterVector<ClusterType>
-ClusterFile<ClusterType, Enable>::read_frame_without_cut() {
-    if (m_mode != "r") {
-        throw std::runtime_error("File not opened for reading");
-    }
-    if (m_num_left) {
-        throw std::runtime_error(
-            "There are still photons left in the last frame");
-    }
-    int32_t frame_number;
-    if (fread(&frame_number, sizeof(frame_number), 1, fp) != 1) {
-        throw std::runtime_error(LOCATION + "Could not read frame number");
-    }
-
-    int32_t n_clusters; // Saved as 32bit integer in the cluster file
-    if (fread(&n_clusters, sizeof(n_clusters), 1, fp) != 1) {
-        throw std::runtime_error(LOCATION +
-                                 "Could not read number of clusters");
-    }
-
-    ClusterVector<ClusterType> clusters(n_clusters);
-    clusters.set_frame_number(frame_number);
-
-    if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) !=
-        static_cast<size_t>(n_clusters)) {
-        throw std::runtime_error(LOCATION + "Could not read clusters");
-    }
-    clusters.resize(n_clusters);
-    if (m_gain_map)
-        m_gain_map->apply_gain_map(clusters);
-    return clusters;
-}
-
-template <typename ClusterType, typename Enable>
-ClusterVector<ClusterType>
-ClusterFile<ClusterType, Enable>::read_frame_with_cut() {
-    if (m_mode != "r") {
-        throw std::runtime_error("File not opened for reading");
-    }
-    if (m_num_left) {
-        throw std::runtime_error(
-            "There are still photons left in the last frame");
-    }
-    int32_t frame_number;
-    if (fread(&frame_number, sizeof(frame_number), 1, fp) != 1) {
-        throw std::runtime_error("Could not read frame number");
-    }
-
-    if (fread(&m_num_left, sizeof(m_num_left), 1, fp) != 1) {
-        throw std::runtime_error("Could not read number of clusters");
-    }
-
-    ClusterVector<ClusterType> clusters;
-    clusters.reserve(m_num_left);
-    clusters.set_frame_number(frame_number);
-    while (m_num_left) {
-        ClusterType c = read_one_cluster();
-        if (is_selected(c)) {
-            clusters.push_back(c);
-        }
-    }
-    if (m_gain_map)
-        m_gain_map->apply_gain_map(clusters);
-    return clusters;
-}
-
-template <typename ClusterType, typename Enable>
-bool ClusterFile<ClusterType, Enable>::is_selected(ClusterType &cl) {
-    // Should fail fast
-    if (m_roi) {
-        if (!(m_roi->contains(cl.x, cl.y))) {
-            return false;
-        }
-    }
-
-    auto cluster_size_x = extract_template_arguments<
-        std::remove_reference_t<decltype(cl)>>::cluster_size_x;
-    auto cluster_size_y = extract_template_arguments<
-        std::remove_reference_t<decltype(cl)>>::cluster_size_y;
-
-    size_t cluster_center_index =
-        (cluster_size_x / 2) + (cluster_size_y / 2) * cluster_size_x;
-
-    if (m_noise_map) {
-        auto sum_1x1 = cl.data[cluster_center_index]; // central pixel
-        auto sum_2x2 = cl.max_sum_2x2().first; // highest sum of 2x2 subclusters
-        auto total_sum = cl.sum();             // sum of all pixels
-
-        auto noise =
-            (*m_noise_map)(cl.y, cl.x); // TODO! check if this is correct
-        if (sum_1x1 <= noise || sum_2x2 <= 2 * noise ||
-            total_sum <= 3 * noise) {
-            return false;
-        }
-    }
-    // we passed all checks
-    return true;
-}
 
 } // namespace aare

include/aare/ClusterFileSink.hpp

@@ -3,41 +3,35 @@
 #include <filesystem>
 #include <thread>
 
-#include "aare/ClusterFinderMT.hpp"
-#include "aare/ClusterVector.hpp"
 #include "aare/ProducerConsumerQueue.hpp"
+#include "aare/ClusterVector.hpp"
+#include "aare/ClusterFinderMT.hpp"
 
-namespace aare {
+namespace aare{
 
-template <typename ClusterType,
-          typename = std::enable_if_t<is_cluster_v<ClusterType>>>
-class ClusterFileSink {
-    ProducerConsumerQueue<ClusterVector<ClusterType>> *m_source;
+class ClusterFileSink{
+    ProducerConsumerQueue<ClusterVector<int>>* m_source;
     std::atomic<bool> m_stop_requested{false};
     std::atomic<bool> m_stopped{true};
     std::chrono::milliseconds m_default_wait{1};
     std::thread m_thread;
     std::ofstream m_file;
 
-    void process() {
+
+    void process(){
         m_stopped = false;
         fmt::print("ClusterFileSink started\n");
-        while (!m_stop_requested || !m_source->isEmpty()) {
-            if (ClusterVector<ClusterType> *clusters = m_source->frontPtr();
+        while (!m_stop_requested || !m_source->isEmpty()) { 
+            if (ClusterVector<int> *clusters = m_source->frontPtr();
                 clusters != nullptr) {
                 // Write clusters to file
-                int32_t frame_number =
-                    clusters->frame_number(); // TODO! Should we store frame
-                                              // number already as int?
+                int32_t frame_number = clusters->frame_number(); //TODO! Should we store frame number already as int?
                 uint32_t num_clusters = clusters->size();
-                m_file.write(reinterpret_cast<const char *>(&frame_number),
-                             sizeof(frame_number));
-                m_file.write(reinterpret_cast<const char *>(&num_clusters),
-                             sizeof(num_clusters));
-                m_file.write(reinterpret_cast<const char *>(clusters->data()),
-                             clusters->size() * clusters->item_size());
+                m_file.write(reinterpret_cast<const char*>(&frame_number), sizeof(frame_number));
+                m_file.write(reinterpret_cast<const char*>(&num_clusters), sizeof(num_clusters));
+                m_file.write(reinterpret_cast<const char*>(clusters->data()), clusters->size() * clusters->item_size());
                 m_source->popFront();
-            } else {
+            }else{
                 std::this_thread::sleep_for(m_default_wait);
             }
         }
@@ -45,18 +39,18 @@ class ClusterFileSink {
         m_stopped = true;
     }
 
-  public:
-    ClusterFileSink(ClusterFinderMT<ClusterType, uint16_t, double> *source,
-                    const std::filesystem::path &fname) {
-        m_source = source->sink();
-        m_thread = std::thread(&ClusterFileSink::process, this);
-        m_file.open(fname, std::ios::binary);
-    }
-    void stop() {
-        m_stop_requested = true;
-        m_thread.join();
-        m_file.close();
-    }
+    public:
+        ClusterFileSink(ClusterFinderMT<uint16_t, double, int32_t>* source, const std::filesystem::path& fname){
+            m_source = source->sink();
+            m_thread = std::thread(&ClusterFileSink::process, this);
+            m_file.open(fname, std::ios::binary);
+        }
+        void stop(){
+            m_stop_requested = true;
+            m_thread.join();
+            m_file.close();
+        }
 };
 
+
 } // namespace aare

include/aare/ClusterFileV2.hpp

@@ -1,16 +1,15 @@
 #pragma once
 #include "aare/core/defs.hpp"
 #include <filesystem>
-#include <fmt/format.h>
 #include <string>
+#include <fmt/format.h>
 
 namespace aare {
 struct ClusterHeader {
     int32_t frame_number;
     int32_t n_clusters;
     std::string to_string() const {
-        return "frame_number: " + std::to_string(frame_number) +
-               ", n_clusters: " + std::to_string(n_clusters);
+        return "frame_number: " + std::to_string(frame_number) + ", n_clusters: " + std::to_string(n_clusters);
     }
 };
 
@@ -25,8 +24,7 @@ struct ClusterV2_ {
                 data_str += std::to_string(d) + ", ";
             }
             data_str += "]";
-            return "x: " + std::to_string(x) + ", y: " + std::to_string(y) +
-                   ", data: " + data_str;
+            return "x: " + std::to_string(x) + ", y: " + std::to_string(y) + ", data: " + data_str;
         }
         return "x: " + std::to_string(x) + ", y: " + std::to_string(y);
     }
@@ -36,31 +34,27 @@ struct ClusterV2 {
     ClusterV2_ cluster;
     int32_t frame_number;
     std::string to_string() const {
-        return "frame_number: " + std::to_string(frame_number) + ", " +
-               cluster.to_string();
+        return "frame_number: " + std::to_string(frame_number) + ", " + cluster.to_string();
     }
 };
 
 /**
  * @brief
- * important not: fp always points to the clusters header and does not point to
- * individual clusters
+ * important not: fp always points to the clusters header and does not point to individual clusters
  *
  */
 class ClusterFileV2 {
-    std::filesystem::path m_fpath;
+    std::filesystem::path m_fpath; 
     std::string m_mode;
     FILE *fp{nullptr};
 
-    void check_open() {
+    void check_open(){
         if (!fp)
-            throw std::runtime_error(
-                fmt::format("File: {} not open", m_fpath.string()));
+            throw std::runtime_error(fmt::format("File: {} not open", m_fpath.string()));
     }
 
   public:
-    ClusterFileV2(std::filesystem::path const &fpath, std::string const &mode)
-        : m_fpath(fpath), m_mode(mode) {
+    ClusterFileV2(std::filesystem::path const &fpath, std::string const &mode): m_fpath(fpath), m_mode(mode) {
         if (m_mode != "r" && m_mode != "w")
             throw std::invalid_argument("mode must be 'r' or 'w'");
         if (m_mode == "r" && !std::filesystem::exists(m_fpath))
@@ -83,7 +77,7 @@ class ClusterFileV2 {
         check_open();
 
         ClusterHeader header;
-        fread(&header, sizeof(ClusterHeader), 1, fp);
+        fread(&header, sizeof(ClusterHeader), 1, fp); 
         std::vector<ClusterV2_> clusters_(header.n_clusters);
         fread(clusters_.data(), sizeof(ClusterV2_), header.n_clusters, fp);
         std::vector<ClusterV2> clusters;
@@ -123,7 +117,7 @@ class ClusterFileV2 {
 
     size_t write(std::vector<std::vector<ClusterV2>> const &clusters) {
         check_open();
-        if (m_mode != "w")
+        if (m_mode != "w") 
             throw std::runtime_error("File not opened in write mode");
 
         size_t n_clusters = 0;

include/aare/ClusterFinder.hpp

@@ -10,21 +10,17 @@
 
 namespace aare {
 
-template <typename ClusterType = Cluster<int32_t, 3, 3>,
-          typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double>
+template <typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double,
+          typename CT = int32_t>
 class ClusterFinder {
     Shape<2> m_image_size;
+    const int m_cluster_sizeX;
+    const int m_cluster_sizeY;
     const PEDESTAL_TYPE m_nSigma;
     const PEDESTAL_TYPE c2;
     const PEDESTAL_TYPE c3;
     Pedestal<PEDESTAL_TYPE> m_pedestal;
-    ClusterVector<ClusterType> m_clusters;
-
-    static const uint8_t ClusterSizeX =
-        extract_template_arguments<ClusterType>::cluster_size_x;
-    static const uint8_t ClusterSizeY =
-        extract_template_arguments<ClusterType>::cluster_size_x;
-    using CT = typename extract_template_arguments<ClusterType>::value_type;
+    ClusterVector<CT> m_clusters;
 
   public:
     /**
@@ -35,12 +31,15 @@ class ClusterFinder {
      * @param capacity initial capacity of the cluster vector
      *
      */
-    ClusterFinder(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0,
-                  size_t capacity = 1000000)
-        : m_image_size(image_size), m_nSigma(nSigma),
-          c2(sqrt((ClusterSizeY + 1) / 2 * (ClusterSizeX + 1) / 2)),
-          c3(sqrt(ClusterSizeX * ClusterSizeY)),
-          m_pedestal(image_size[0], image_size[1]), m_clusters(capacity) {};
+    ClusterFinder(Shape<2> image_size, Shape<2> cluster_size,
+                  PEDESTAL_TYPE nSigma = 5.0, size_t capacity = 1000000)
+        : m_image_size(image_size), m_cluster_sizeX(cluster_size[0]),
+          m_cluster_sizeY(cluster_size[1]),
+          m_nSigma(nSigma),
+          c2(sqrt((m_cluster_sizeY + 1) / 2 * (m_cluster_sizeX + 1) / 2)),
+          c3(sqrt(m_cluster_sizeX * m_cluster_sizeY)),
+          m_pedestal(image_size[0], image_size[1]),
+          m_clusters(m_cluster_sizeX, m_cluster_sizeY, capacity) {};
 
     void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) {
         m_pedestal.push(frame);
@@ -57,29 +56,23 @@ class ClusterFinder {
      * same capacity as the old one
      *
      */
-    ClusterVector<ClusterType>
-    steal_clusters(bool realloc_same_capacity = false) {
-        ClusterVector<ClusterType> tmp = std::move(m_clusters);
+    ClusterVector<CT> steal_clusters(bool realloc_same_capacity = false) {
+        ClusterVector<CT> tmp = std::move(m_clusters);
         if (realloc_same_capacity)
-            m_clusters = ClusterVector<ClusterType>(tmp.capacity());
+            m_clusters = ClusterVector<CT>(m_cluster_sizeX, m_cluster_sizeY,
+                                           tmp.capacity());
         else
-            m_clusters = ClusterVector<ClusterType>{};
+            m_clusters = ClusterVector<CT>(m_cluster_sizeX, m_cluster_sizeY);
         return tmp;
     }
     void find_clusters(NDView<FRAME_TYPE, 2> frame, uint64_t frame_number = 0) {
         // // TODO! deal with even size clusters
         // // currently 3,3 -> +/- 1
         // //  4,4 -> +/- 2
-        int dy = ClusterSizeY / 2;
-        int dx = ClusterSizeX / 2;
-        int has_center_pixel_x =
-            ClusterSizeX %
-            2; // for even sized clusters there is no proper cluster center and
-               // even amount of pixels around the center
-        int has_center_pixel_y = ClusterSizeY % 2;
-
+        int dy = m_cluster_sizeY / 2;
+        int dx = m_cluster_sizeX / 2;
         m_clusters.set_frame_number(frame_number);
-        std::vector<CT> cluster_data(ClusterSizeX * ClusterSizeY);
+        std::vector<CT> cluster_data(m_cluster_sizeX * m_cluster_sizeY);
         for (int iy = 0; iy < frame.shape(0); iy++) {
             for (int ix = 0; ix < frame.shape(1); ix++) {
 
@@ -94,8 +87,8 @@ class ClusterFinder {
                     continue; // NEGATIVE_PEDESTAL go to next pixel
                               // TODO! No pedestal update???
 
-                for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) {
-                    for (int ic = -dx; ic < dx + has_center_pixel_x; ic++) {
+                for (int ir = -dy; ir < dy + 1; ir++) {
+                    for (int ic = -dx; ic < dx + 1; ic++) {
                         if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
                             iy + ir >= 0 && iy + ir < frame.shape(0)) {
                             PEDESTAL_TYPE val =
@@ -116,12 +109,8 @@ class ClusterFinder {
                     // pass
                 } else {
                     // m_pedestal.push(iy, ix, frame(iy, ix));   // Safe option
-                    m_pedestal.push_fast(
-                        iy, ix,
-                        frame(iy,
-                              ix)); // Assume we have reached n_samples in the
-                                    // pedestal, slight performance improvement
-                    continue;       // It was a pedestal value nothing to store
+                    m_pedestal.push_fast(iy, ix, frame(iy, ix)); // Assume we have reached n_samples in the pedestal, slight performance improvement
+                    continue; // It was a pedestal value nothing to store
                 }
 
                 // Store cluster
@@ -133,14 +122,13 @@ class ClusterFinder {
                     // It's worth redoing the look since most of the time we
                     // don't have a photon
                     int i = 0;
-                    for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) {
-                        for (int ic = -dx; ic < dx + has_center_pixel_y; ic++) {
+                    for (int ir = -dy; ir < dy + 1; ir++) {
+                        for (int ic = -dx; ic < dx + 1; ic++) {
                             if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
                                 iy + ir >= 0 && iy + ir < frame.shape(0)) {
                                 CT tmp =
                                     static_cast<CT>(frame(iy + ir, ix + ic)) -
-                                    static_cast<CT>(
-                                        m_pedestal.mean(iy + ir, ix + ic));
+                                    m_pedestal.mean(iy + ir, ix + ic);
                                 cluster_data[i] =
                                     tmp; // Watch for out of bounds access
                                 i++;
@@ -148,13 +136,10 @@ class ClusterFinder {
                         }
                     }
 
-                    ClusterType new_cluster{};
-                    new_cluster.x = ix;
-                    new_cluster.y = iy;
-                    std::copy(cluster_data.begin(), cluster_data.end(),
-                              new_cluster.data);
                     // Add the cluster to the output ClusterVector
-                    m_clusters.push_back(new_cluster);
+                    m_clusters.push_back(
+                        ix, iy,
+                        reinterpret_cast<std::byte *>(cluster_data.data()));
                 }
             }
         }

include/aare/ClusterFinderMT.hpp

@@ -30,16 +30,14 @@ struct FrameWrapper {
  * @tparam PEDESTAL_TYPE type of the pedestal data
  * @tparam CT type of the cluster data
  */
-template <typename ClusterType = Cluster<int32_t, 3, 3>,
-          typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double>
+template <typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double,
+          typename CT = int32_t>
 class ClusterFinderMT {
-
-    using CT = typename extract_template_arguments<ClusterType>::value_type;
     size_t m_current_thread{0};
     size_t m_n_threads{0};
-    using Finder = ClusterFinder<ClusterType, FRAME_TYPE, PEDESTAL_TYPE>;
+    using Finder = ClusterFinder<FRAME_TYPE, PEDESTAL_TYPE, CT>;
     using InputQueue = ProducerConsumerQueue<FrameWrapper>;
-    using OutputQueue = ProducerConsumerQueue<ClusterVector<ClusterType>>;
+    using OutputQueue = ProducerConsumerQueue<ClusterVector<int>>;
     std::vector<std::unique_ptr<InputQueue>> m_input_queues;
     std::vector<std::unique_ptr<OutputQueue>> m_output_queues;
 
@@ -68,8 +66,7 @@ class ClusterFinderMT {
                 switch (frame->type) {
                 case FrameType::DATA:
                     cf->find_clusters(frame->data.view(), frame->frame_number);
-                    m_output_queues[thread_id]->write(
-                        cf->steal_clusters(realloc_same_capacity));
+                    m_output_queues[thread_id]->write(cf->steal_clusters(realloc_same_capacity));
                     break;
 
                 case FrameType::PEDESTAL:
@@ -117,31 +114,28 @@ class ClusterFinderMT {
      * expected number of clusters in a frame per frame.
      * @param n_threads number of threads to use
      */
-    ClusterFinderMT(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0,
-                    size_t capacity = 2000, size_t n_threads = 3)
+    ClusterFinderMT(Shape<2> image_size, Shape<2> cluster_size,
+                    PEDESTAL_TYPE nSigma = 5.0, size_t capacity = 2000,
+                    size_t n_threads = 3)
         : m_n_threads(n_threads) {
         for (size_t i = 0; i < n_threads; i++) {
             m_cluster_finders.push_back(
-                std::make_unique<
-                    ClusterFinder<ClusterType, FRAME_TYPE, PEDESTAL_TYPE>>(
-                    image_size, nSigma, capacity));
+                std::make_unique<ClusterFinder<FRAME_TYPE, PEDESTAL_TYPE, CT>>(
+                    image_size, cluster_size, nSigma, capacity));
         }
         for (size_t i = 0; i < n_threads; i++) {
             m_input_queues.emplace_back(std::make_unique<InputQueue>(200));
             m_output_queues.emplace_back(std::make_unique<OutputQueue>(200));
         }
-        // TODO! Should we start automatically?
+        //TODO! Should we start automatically?
         start();
     }
 
     /**
      * @brief Return the sink queue where all the clusters are collected
-     * @warning You need to empty this queue otherwise the cluster finder will
-     * wait forever
+     * @warning You need to empty this queue otherwise the cluster finder will wait forever
      */
-    ProducerConsumerQueue<ClusterVector<ClusterType>> *sink() {
-        return &m_sink;
-    }
+    ProducerConsumerQueue<ClusterVector<int>> *sink() { return &m_sink; }
 
     /**
      * @brief Start all processing threads

include/aare/ClusterVector.hpp

@@ -1,5 +1,4 @@
 #pragma once
-#include "aare/Cluster.hpp" //TODO maybe store in seperate file !!!
 #include <algorithm>
 #include <array>
 #include <cstddef>
@@ -14,10 +13,6 @@
 
 namespace aare {
 
-template <typename ClusterType,
-          typename = std::enable_if_t<is_cluster_v<ClusterType>>>
-class ClusterVector; // Forward declaration
-
 /**
  * @brief ClusterVector is a container for clusters of various sizes. It uses a
  * contiguous memory buffer to store the clusters. It is templated on the data
@@ -29,16 +24,15 @@ class ClusterVector; // Forward declaration
  * @tparam CoordType data type of the x and y coordinates of the cluster
  * (normally int16_t)
  */
-#if 0
-template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType>
-class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
-
+template <typename T, typename CoordType = int16_t> class ClusterVector {
+    using value_type = T;
+    size_t m_cluster_size_x;
+    size_t m_cluster_size_y;
     std::byte *m_data{};
     size_t m_size{0};
     size_t m_capacity;
     uint64_t m_frame_number{0}; // TODO! Check frame number size and type
-    /** 
+    /*
     Format string used in the python bindings to create a numpy
     array from the buffer
     = - native byte order
@@ -49,26 +43,29 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
     constexpr static char m_fmt_base[] = "=h:x:\nh:y:\n({},{}){}:data:";
 
   public:
-    using value_type = T;
-    using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
-
     /**
      * @brief Construct a new ClusterVector object
+     * @param cluster_size_x size of the cluster in x direction
+     * @param cluster_size_y size of the cluster in y direction
      * @param capacity initial capacity of the buffer in number of clusters
      * @param frame_number frame number of the clusters. Default is 0, which is
      * also used to indicate that the clusters come from many frames
      */
-    ClusterVector(size_t capacity = 1024, uint64_t frame_number = 0)
-        : m_capacity(capacity), m_frame_number(frame_number) {
-        allocate_buffer(m_capacity);
+    ClusterVector(size_t cluster_size_x = 3, size_t cluster_size_y = 3,
+                  size_t capacity = 1024, uint64_t frame_number = 0)
+        : m_cluster_size_x(cluster_size_x), m_cluster_size_y(cluster_size_y),
+          m_capacity(capacity), m_frame_number(frame_number) {
+        allocate_buffer(capacity);
     }
 
     ~ClusterVector() { delete[] m_data; }
 
     // Move constructor
     ClusterVector(ClusterVector &&other) noexcept
-        : m_data(other.m_data), m_size(other.m_size),
-          m_capacity(other.m_capacity), m_frame_number(other.m_frame_number) {
+        : m_cluster_size_x(other.m_cluster_size_x),
+          m_cluster_size_y(other.m_cluster_size_y), m_data(other.m_data),
+          m_size(other.m_size), m_capacity(other.m_capacity),
+          m_frame_number(other.m_frame_number) {
         other.m_data = nullptr;
         other.m_size = 0;
         other.m_capacity = 0;
@@ -78,6 +75,8 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
     ClusterVector &operator=(ClusterVector &&other) noexcept {
         if (this != &other) {
             delete[] m_data;
+            m_cluster_size_x = other.m_cluster_size_x;
+            m_cluster_size_y = other.m_cluster_size_y;
             m_data = other.m_data;
             m_size = other.m_size;
             m_capacity = other.m_capacity;
@@ -104,22 +103,26 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
 
     /**
      * @brief Add a cluster to the vector
+     * @param x x-coordinate of the cluster
+     * @param y y-coordinate of the cluster
+     * @param data pointer to the data of the cluster
+     * @warning The data pointer must point to a buffer of size cluster_size_x *
+     * cluster_size_y * sizeof(T)
      */
-    void push_back(const ClusterType &cluster) {
+    void push_back(CoordType x, CoordType y, const std::byte *data) {
         if (m_size == m_capacity) {
             allocate_buffer(m_capacity * 2);
         }
         std::byte *ptr = element_ptr(m_size);
-        *reinterpret_cast<CoordType *>(ptr) = cluster.x;
+        *reinterpret_cast<CoordType *>(ptr) = x;
         ptr += sizeof(CoordType);
-        *reinterpret_cast<CoordType *>(ptr) = cluster.y;
+        *reinterpret_cast<CoordType *>(ptr) = y;
         ptr += sizeof(CoordType);
 
-        std::memcpy(ptr, cluster.data, ClusterSizeX * ClusterSizeY * sizeof(T));
-
+        std::copy(data, data + m_cluster_size_x * m_cluster_size_y * sizeof(T),
+                  ptr);
         m_size++;
     }
-
     ClusterVector &operator+=(const ClusterVector &other) {
         if (m_size + other.m_size > m_capacity) {
             allocate_buffer(m_capacity + other.m_size);
@@ -134,11 +137,10 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
      * @brief Sum the pixels in each cluster
      * @return std::vector<T> vector of sums for each cluster
      */
-    /*
     std::vector<T> sum() {
         std::vector<T> sums(m_size);
         const size_t stride = item_size();
-        const size_t n_pixels = ClusterSizeX * ClusterSizeY;
+        const size_t n_pixels = m_cluster_size_x * m_cluster_size_y;
         std::byte *ptr = m_data + 2 * sizeof(CoordType); // skip x and y
 
         for (size_t i = 0; i < m_size; i++) {
@@ -149,33 +151,43 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
         }
         return sums;
     }
-    */
 
     /**
-     * @brief Sum the pixels in the 2x2 subcluster with the biggest pixel sum in
-     * each cluster
+     * @brief Return the maximum sum of the 2x2 subclusters in each cluster
      * @return std::vector<T> vector of sums for each cluster
-     */ //TODO if underlying container is a vector use std::for_each
-    /*
+     * @throws std::runtime_error if the cluster size is not 3x3
+     * @warning Only 3x3 clusters are supported for the 2x2 sum.
+     */
     std::vector<T> sum_2x2() {
-        std::vector<T> sums_2x2(m_size);
+        std::vector<T> sums(m_size);
+        const size_t stride = item_size();
+
+        if (m_cluster_size_x != 3 || m_cluster_size_y != 3) {
+            throw std::runtime_error(
+                "Only 3x3 clusters are supported for the 2x2 sum.");
+        }
+        std::byte *ptr = m_data + 2 * sizeof(CoordType); // skip x and y
 
         for (size_t i = 0; i < m_size; i++) {
-            sums_2x2[i] = at(i).max_sum_2x2;
+            std::array<T, 4> total;
+            auto T_ptr = reinterpret_cast<T *>(ptr);
+            total[0] = T_ptr[0] + T_ptr[1] + T_ptr[3] + T_ptr[4];
+            total[1] = T_ptr[1] + T_ptr[2] + T_ptr[4] + T_ptr[5];
+            total[2] = T_ptr[3] + T_ptr[4] + T_ptr[6] + T_ptr[7];
+            total[3] = T_ptr[4] + T_ptr[5] + T_ptr[7] + T_ptr[8];
+
+            sums[i] = *std::max_element(total.begin(), total.end());
+            ptr += stride;
         }
-        return sums_2x2;
+
+        return sums;
     }
-    */
 
     /**
      * @brief Return the number of clusters in the vector
      */
     size_t size() const { return m_size; }
 
-    uint8_t cluster_size_x() const { return ClusterSizeX; }
-
-    uint8_t cluster_size_y() const { return ClusterSizeY; }
-
     /**
      * @brief Return the capacity of the buffer in number of clusters. This is
      * the number of clusters that can be stored in the current buffer without
@@ -187,7 +199,8 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
      * @brief Return the size in bytes of a single cluster
      */
     size_t item_size() const {
-        return 2 * sizeof(CoordType) + ClusterSizeX * ClusterSizeY * sizeof(T);
+        return 2 * sizeof(CoordType) +
+               m_cluster_size_x * m_cluster_size_y * sizeof(T);
     }
 
     /**
@@ -207,6 +220,9 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
         return m_data + element_offset(i);
     }
 
+    size_t cluster_size_x() const { return m_cluster_size_x; }
+    size_t cluster_size_y() const { return m_cluster_size_y; }
+
     std::byte *data() { return m_data; }
     std::byte const *data() const { return m_data; }
 
@@ -214,12 +230,8 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
      * @brief Return a reference to the i-th cluster casted to type V
      * @tparam V type of the cluster
      */
-    ClusterType &at(size_t i) {
-        return *reinterpret_cast<ClusterType *>(element_ptr(i));
-    }
-
-    const ClusterType &at(size_t i) const {
-        return *reinterpret_cast<const ClusterType *>(element_ptr(i));
+    template <typename V> V &at(size_t i) {
+        return *reinterpret_cast<V *>(element_ptr(i));
     }
 
     template <typename V> const V &at(size_t i) const {
@@ -231,7 +243,7 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
         return m_fmt_base;
     }
 
-    /** 
+    /**
      * @brief Return the frame number of the clusters. 0 is used to indicate
      * that the clusters come from many frames
      */
@@ -241,7 +253,7 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
         m_frame_number = frame_number;
     }
 
-    /** 
+    /**
      * @brief Resize the vector to contain new_size clusters. If new_size is
      * greater than the current capacity, a new buffer is allocated. If the size
      * is smaller no memory is freed, size is just updated.
@@ -256,6 +268,28 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
         m_size = new_size;
     }
 
+    void apply_gain_map(const NDView<double> gain_map){
+        //in principle we need to know the size of the image for this lookup
+        //TODO! check orientations
+        std::array<int64_t, 9> xcorr = {-1, 0, 1, -1, 0, 1, -1, 0, 1};
+        std::array<int64_t, 9> ycorr = {-1, -1, -1, 0, 0, 0, 1, 1, 1};
+        for (size_t i=0; i<m_size; i++){
+            auto& cl = at<Cluster3x3>(i);
+
+            if (cl.x > 0 && cl.y > 0 && cl.x < gain_map.shape(1)-1 && cl.y < gain_map.shape(0)-1){
+                for (size_t j=0; j<9; j++){
+                    size_t x = cl.x + xcorr[j];
+                    size_t y = cl.y + ycorr[j];
+                    cl.data[j] = static_cast<T>(cl.data[j] * gain_map(y, x));
+                }
+            }else{
+                memset(cl.data, 0, 9*sizeof(T)); //clear edge clusters
+            }
+
+            
+        }
+    }
+
   private:
     void allocate_buffer(size_t new_capacity) {
         size_t num_bytes = item_size() * new_capacity;
@@ -266,124 +300,5 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
         m_capacity = new_capacity;
     }
 };
-#endif
-
-/**
- * @brief ClusterVector is a container for clusters of various sizes. It
- * uses a contiguous memory buffer to store the clusters. It is templated on
- * the data type and the coordinate type of the clusters.
- * @note push_back can invalidate pointers to elements in the container
- * @warning ClusterVector is currently move only to catch unintended copies,
- * but this might change since there are probably use cases where copying is
- * needed.
- * @tparam T data type of the pixels in the cluster
- * @tparam CoordType data type of the x and y coordinates of the cluster
- * (normally int16_t)
- */
-template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType>
-class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
-
-    std::vector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> m_data{};
-    uint64_t m_frame_number{0}; // TODO! Check frame number size and type
-
-  public:
-    using value_type = T;
-    using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
-
-    /**
-     * @brief Construct a new ClusterVector object
-     * @param capacity initial capacity of the buffer in number of clusters
-     * @param frame_number frame number of the clusters. Default is 0, which is
-     * also used to indicate that the clusters come from many frames
-     */
-    ClusterVector(size_t capacity = 1024, uint64_t frame_number = 0)
-        : m_frame_number(frame_number) {
-        m_data.reserve(capacity);
-    }
-
-    // Move constructor
-    ClusterVector(ClusterVector &&other) noexcept
-        : m_data(other.m_data), m_frame_number(other.m_frame_number) {
-        other.m_data.clear();
-    }
-
-    // Move assignment operator
-    ClusterVector &operator=(ClusterVector &&other) noexcept {
-        if (this != &other) {
-            m_data = other.m_data;
-            m_frame_number = other.m_frame_number;
-            other.m_data.clear();
-            other.m_frame_number = 0;
-        }
-        return *this;
-    }
-
-    /**
-     * @brief Reserve space for at least capacity clusters
-     * @param capacity number of clusters to reserve space for
-     * @note If capacity is less than the current capacity, the function does
-     * nothing.
-     */
-    void reserve(size_t capacity) { m_data.reserve(capacity); }
-
-    void resize(size_t size) { m_data.resize(size); }
-
-    void push_back(const ClusterType &cluster) { m_data.push_back(cluster); }
-
-    ClusterVector &operator+=(const ClusterVector &other) {
-        m_data.insert(m_data.end(), other.begin(), other.end());
-
-        return *this;
-    }
-
-    /**
-     * @brief Return the number of clusters in the vector
-     */
-    size_t size() const { return m_data.size(); }
-
-    uint8_t cluster_size_x() const { return ClusterSizeX; }
-
-    uint8_t cluster_size_y() const { return ClusterSizeY; }
-
-    /**
-     * @brief Return the capacity of the buffer in number of clusters. This is
-     * the number of clusters that can be stored in the current buffer without
-     * reallocation.
-     */
-    size_t capacity() const { return m_data.capacity(); }
-
-    const auto begin() const { return m_data.begin(); }
-
-    const auto end() const { return m_data.end(); }
-
-    /**
-     * @brief Return the size in bytes of a single cluster
-     */
-    size_t item_size() const {
-        return 2 * sizeof(CoordType) + ClusterSizeX * ClusterSizeY * sizeof(T);
-    }
-
-    ClusterType *data() { return m_data.data(); }
-    ClusterType const *data() const { return m_data.data(); }
-
-    /**
-     * @brief Return a reference to the i-th cluster casted to type V
-     * @tparam V type of the cluster
-     */
-    ClusterType &at(size_t i) { return m_data[i]; }
-
-    const ClusterType &at(size_t i) const { return m_data[i]; }
-
-    /**
-     * @brief Return the frame number of the clusters. 0 is used to indicate
-     * that the clusters come from many frames
-     */
-    uint64_t frame_number() const { return m_frame_number; }
-
-    void set_frame_number(uint64_t frame_number) {
-        m_frame_number = frame_number;
-    }
-};
 
 } // namespace aare

include/aare/FilePtr.hpp

@@ -0,0 +1,30 @@
+#pragma once
+#include <cstdio>
+#include <filesystem>
+
+namespace aare {
+
+/**
+ * \brief RAII wrapper for FILE pointer
+ */
+class FilePtr {
+    FILE *fp_{nullptr};
+
+  public:
+    FilePtr() = default;
+    FilePtr(const std::filesystem::path& fname, const std::string& mode);
+    FilePtr(const FilePtr &) = delete;            // we don't want a copy
+    FilePtr &operator=(const FilePtr &) = delete; // since we handle a resource
+    FilePtr(FilePtr &&other);
+    FilePtr &operator=(FilePtr &&other);
+    FILE *get();
+    int64_t tell();
+    void seek(int64_t offset, int whence = SEEK_SET) {
+        if (fseek(fp_, offset, whence) != 0)
+            throw std::runtime_error("Error seeking in file");
+    }
+    std::string error_msg();
+    ~FilePtr();
+};
+
+} // namespace aare

include/aare/GainMap.hpp

@@ -1,58 +0,0 @@
-/************************************************
- * @file ApplyGainMap.hpp
- * @short function to apply gain map of image size to a vector of clusters
- ***********************************************/
-
-#pragma once
-#include "aare/Cluster.hpp"
-#include "aare/ClusterVector.hpp"
-#include "aare/NDArray.hpp"
-#include "aare/NDView.hpp"
-#include <memory>
-
-namespace aare {
-
-class GainMap {
-
-  public:
-    explicit GainMap(const NDArray<double, 2> &gain_map)
-        : m_gain_map(gain_map) {};
-
-    explicit GainMap(const NDView<double, 2> gain_map) {
-        m_gain_map = NDArray<double, 2>(gain_map);
-    }
-
-    template <typename ClusterType,
-              typename = std::enable_if_t<is_cluster_v<ClusterType>>>
-    void apply_gain_map(ClusterVector<ClusterType> &clustervec) {
-        // in principle we need to know the size of the image for this lookup
-        size_t ClusterSizeX = clustervec.cluster_size_x();
-        size_t ClusterSizeY = clustervec.cluster_size_y();
-
-        using T = typename ClusterVector<ClusterType>::value_type;
-
-        int64_t index_cluster_center_x = ClusterSizeX / 2;
-        int64_t index_cluster_center_y = ClusterSizeY / 2;
-        for (size_t i = 0; i < clustervec.size(); i++) {
-            auto &cl = clustervec.at(i);
-
-            if (cl.x > 0 && cl.y > 0 && cl.x < m_gain_map.shape(1) - 1 &&
-                cl.y < m_gain_map.shape(0) - 1) {
-                for (size_t j = 0; j < ClusterSizeX * ClusterSizeY; j++) {
-                    size_t x = cl.x + j % ClusterSizeX - index_cluster_center_x;
-                    size_t y = cl.y + j / ClusterSizeX - index_cluster_center_y;
-                    cl.data[j] = cl.data[j] * static_cast<T>(m_gain_map(y, x));
-                }
-            } else {
-                memset(cl.data, 0,
-                       ClusterSizeX * ClusterSizeY *
-                           sizeof(T)); // clear edge clusters
-            }
-        }
-    }
-
-  private:
-    NDArray<double, 2> m_gain_map{};
-};
-
-} // end of namespace aare

include/aare/Interpolator.hpp

@@ -1,132 +1,29 @@
 #pragma once
-
-#include "aare/CalculateEta.hpp"
-#include "aare/Cluster.hpp"
-#include "aare/ClusterFile.hpp" //Cluster_3x3
-#include "aare/ClusterVector.hpp"
 #include "aare/NDArray.hpp"
 #include "aare/NDView.hpp"
-#include "aare/algorithm.hpp"
+#include "aare/ClusterVector.hpp"
+#include "aare/ClusterFile.hpp" //Cluster_3x3
+namespace aare{
 
-namespace aare {
-
-struct Photon {
+struct Photon{
     double x;
     double y;
     double energy;
 };
 
-class Interpolator {
+class Interpolator{
     NDArray<double, 3> m_ietax;
     NDArray<double, 3> m_ietay;
 
     NDArray<double, 1> m_etabinsx;
     NDArray<double, 1> m_etabinsy;
     NDArray<double, 1> m_energy_bins;
+    public:
+        Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins, NDView<double, 1> ybins, NDView<double, 1> ebins);
+        NDArray<double, 3> get_ietax(){return m_ietax;}
+        NDArray<double, 3> get_ietay(){return m_ietay;}
 
-  public:
-    Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
-                 NDView<double, 1> ybins, NDView<double, 1> ebins);
-    NDArray<double, 3> get_ietax() { return m_ietax; }
-    NDArray<double, 3> get_ietay() { return m_ietay; }
-
-    template <typename ClusterType,
-              typename Eanble = std::enable_if_t<is_cluster_v<ClusterType>>>
-    std::vector<Photon> interpolate(const ClusterVector<ClusterType> &clusters);
+        std::vector<Photon> interpolate(const ClusterVector<int32_t>& clusters);
 };
 
-// TODO: generalize to support any clustertype!!! otherwise add std::enable_if_t
-// to only take Cluster2x2 and Cluster3x3
-template <typename ClusterType, typename Enable>
-std::vector<Photon>
-Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) {
-    std::vector<Photon> photons;
-    photons.reserve(clusters.size());
-
-    if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) {
-        for (size_t i = 0; i < clusters.size(); i++) {
-
-            auto cluster = clusters.at(i);
-            auto eta = calculate_eta2(cluster);
-
-            Photon photon;
-            photon.x = cluster.x;
-            photon.y = cluster.y;
-            photon.energy = eta.sum;
-
-            // auto ie = nearest_index(m_energy_bins, photon.energy)-1;
-            // auto ix = nearest_index(m_etabinsx, eta.x)-1;
-            // auto iy = nearest_index(m_etabinsy, eta.y)-1;
-            // Finding the index of the last element that is smaller
-            // should work fine as long as we have many bins
-            auto ie = last_smaller(m_energy_bins, photon.energy);
-            auto ix = last_smaller(m_etabinsx, eta.x);
-            auto iy = last_smaller(m_etabinsy, eta.y);
-
-            // fmt::print("ex: {}, ix: {}, iy: {}\n", ie, ix, iy);
-
-            double dX, dY;
-            // cBottomLeft = 0,
-            // cBottomRight = 1,
-            // cTopLeft = 2,
-            // cTopRight = 3
-            switch (eta.c) {
-            case cTopLeft:
-                dX = -1.;
-                dY = 0;
-                break;
-            case cTopRight:;
-                dX = 0;
-                dY = 0;
-                break;
-            case cBottomLeft:
-                dX = -1.;
-                dY = -1.;
-                break;
-            case cBottomRight:
-                dX = 0.;
-                dY = -1.;
-                break;
-            }
-            photon.x += m_ietax(ix, iy, ie) * 2 + dX;
-            photon.y += m_ietay(ix, iy, ie) * 2 + dY;
-            photons.push_back(photon);
-        }
-    } else if (clusters.cluster_size_x() == 2 ||
-               clusters.cluster_size_y() == 2) {
-        for (size_t i = 0; i < clusters.size(); i++) {
-            auto cluster = clusters.at(i);
-            auto eta = calculate_eta2(cluster);
-
-            Photon photon;
-            photon.x = cluster.x;
-            photon.y = cluster.y;
-            photon.energy = eta.sum;
-
-            // Now do some actual interpolation.
-            // Find which energy bin the cluster is in
-            //  auto ie = nearest_index(m_energy_bins, photon.energy)-1;
-            //  auto ix = nearest_index(m_etabinsx, eta.x)-1;
-            //  auto iy = nearest_index(m_etabinsy, eta.y)-1;
-            // Finding the index of the last element that is smaller
-            // should work fine as long as we have many bins
-            auto ie = last_smaller(m_energy_bins, photon.energy);
-            auto ix = last_smaller(m_etabinsx, eta.x);
-            auto iy = last_smaller(m_etabinsy, eta.y);
-
-            photon.x += m_ietax(ix, iy, ie) *
-                        2; // eta goes between 0 and 1 but we could move the hit
-                           // anywhere in the 2x2
-            photon.y += m_ietay(ix, iy, ie) * 2;
-            photons.push_back(photon);
-        }
-
-    } else {
-        throw std::runtime_error(
-            "Only 3x3 and 2x2 clusters are supported for interpolation");
-    }
-
-    return photons;
-}
-
 } // namespace aare

include/aare/JungfrauDataFile.hpp

@@ -0,0 +1,106 @@
+#pragma once
+#include <cstdint>
+#include <filesystem>
+#include <vector>
+
+#include "aare/FilePtr.hpp"
+#include "aare/defs.hpp"
+#include "aare/NDArray.hpp"
+#include "aare/FileInterface.hpp"
+namespace aare {
+
+
+struct JungfrauDataHeader{
+    uint64_t framenum;
+    uint64_t bunchid;
+};
+
+class JungfrauDataFile : public FileInterface {
+
+    size_t m_rows{};     //!< number of rows in the image, from find_frame_size();
+    size_t m_cols{};     //!< number of columns in the image, from find_frame_size();
+    size_t m_bytes_per_frame{}; //!< number of bytes per frame excluding header
+    size_t m_total_frames{};    //!< total number of frames in the series of files
+    size_t m_offset{}; //!< file index of the first file, allow starting at non zero file
+    size_t m_current_file_index{};   //!< The index of the open file
+    size_t m_current_frame_index{};  //!< The index of the current frame (with reference to all files)
+
+    std::vector<size_t> m_last_frame_in_file{}; //!< Used for seeking to the correct file
+    std::filesystem::path m_path; //!< path to the files
+    std::string m_base_name;      //!< base name used for formatting file names
+
+    FilePtr m_fp; //!< RAII wrapper for a FILE* 
+
+    
+    using pixel_type = uint16_t;
+    static constexpr size_t header_size = sizeof(JungfrauDataHeader); 
+    static constexpr size_t n_digits_in_file_index = 6; //!< to format file names
+
+  public:
+    JungfrauDataFile(const std::filesystem::path &fname);
+
+    std::string base_name() const;      //!< get the base name of the file (without path and extension)
+    size_t bytes_per_frame() override;           
+    size_t pixels_per_frame() override;    
+    size_t bytes_per_pixel() const;  
+    size_t bitdepth() const override;
+    void seek(size_t frame_index) override;  //!< seek to the given frame index (note not byte offset)
+    size_t tell() override;            //!< get the frame index of the file pointer
+    size_t total_frames() const override;
+    size_t rows() const override;
+    size_t cols() const override;
+    std::array<ssize_t,2> shape() const;
+    size_t n_files() const; //!< get the number of files in the series. 
+
+    // Extra functions needed for FileInterface
+    Frame read_frame() override;
+    Frame read_frame(size_t frame_number) override;
+    std::vector<Frame> read_n(size_t n_frames=0) override;
+    void read_into(std::byte *image_buf) override;
+    void read_into(std::byte *image_buf, size_t n_frames) override;
+    size_t frame_number(size_t frame_index) override;
+    DetectorType detector_type() const override;
+
+    /**
+     * @brief Read a single frame from the file into the given buffer.  
+     * @param image_buf buffer to read the frame into. (Note the caller is responsible for allocating the buffer)
+     * @param header pointer to a JungfrauDataHeader or nullptr to skip header)
+     */
+    void read_into(std::byte *image_buf, JungfrauDataHeader *header = nullptr);
+
+    /**
+     * @brief Read a multiple frames from the file into the given buffer.  
+     * @param image_buf buffer to read the frame into. (Note the caller is responsible for allocating the buffer)
+     * @param n_frames number of frames to read
+     * @param header pointer to a JungfrauDataHeader or nullptr to skip header)
+     */
+    void read_into(std::byte *image_buf, size_t n_frames, JungfrauDataHeader *header = nullptr);
+    
+    /** 
+     * @brief Read a single frame from the file into the given NDArray
+     * @param image NDArray to read the frame into.
+     */
+    void read_into(NDArray<uint16_t>* image, JungfrauDataHeader* header = nullptr);
+
+    JungfrauDataHeader read_header();
+    std::filesystem::path current_file() const { return fpath(m_current_file_index+m_offset); }
+
+    
+    private:
+    /**
+     * @brief Find the size of the frame in the file. (256x256, 256x1024, 512x1024)
+     * @param fname path to the file
+     * @throws std::runtime_error if the file is empty or the size cannot be determined
+     */
+     void find_frame_size(const std::filesystem::path &fname);
+
+
+     void parse_fname(const std::filesystem::path &fname);
+     void scan_files();
+     void open_file(size_t file_index);
+     std::filesystem::path fpath(size_t frame_index) const;
+
+
+  };
+
+} // namespace aare

include/aare/NDArray.hpp

@@ -194,7 +194,7 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
 
     T *data() { return data_; }
     std::byte *buffer() { return reinterpret_cast<std::byte *>(data_); }
-    size_t size() const { return size_; }
+    ssize_t size() const { return static_cast<ssize_t>(size_); }
     size_t total_bytes() const { return size_ * sizeof(T); }
     std::array<int64_t, Ndim> shape() const noexcept { return shape_; }
     int64_t shape(int64_t i) const noexcept { return shape_[i]; }

include/aare/NDView.hpp

@@ -71,7 +71,7 @@ template <typename T, int64_t Ndim = 2> class NDView : public ArrayExpr<NDView<T
         return buffer_[element_offset(strides_, index...)];
     }
 
-    size_t size() const { return size_; }
+    ssize_t size() const { return static_cast<ssize_t>(size_); }
     size_t total_bytes() const { return size_ * sizeof(T); }
     std::array<int64_t, Ndim> strides() const noexcept { return strides_; }
 
@@ -102,7 +102,7 @@ template <typename T, int64_t Ndim = 2> class NDView : public ArrayExpr<NDView<T
 
     template<size_t Size>
     NDView& operator=(const std::array<T, Size> &arr) {
-        if(size() != arr.size())
+        if(size() != static_cast<ssize_t>(arr.size()))
             throw std::runtime_error(LOCATION + "Array and NDView size mismatch");
         std::copy(arr.begin(), arr.end(), begin());
         return *this;
@@ -184,4 +184,9 @@ std::ostream& operator <<(std::ostream& os, const NDView<T, Ndim>& arr){
 }
 
 
+template <typename T>
+NDView<T,1> make_view(std::vector<T>& vec){
+    return NDView<T,1>(vec.data(), {static_cast<int64_t>(vec.size())});
+}
+
 } // namespace aare

include/aare/RawSubFile.hpp

@@ -22,7 +22,7 @@ class RawSubFile {
     size_t m_rows{};
     size_t m_cols{};
     size_t m_bytes_per_frame{};
-    size_t n_frames{};
+    size_t m_num_frames{};
     uint32_t m_pos_row{};
     uint32_t m_pos_col{};
  
@@ -53,6 +53,7 @@ class RawSubFile {
     size_t tell();
 
     void read_into(std::byte *image_buf, DetectorHeader *header = nullptr);
+    void read_into(std::byte *image_buf, size_t n_frames, DetectorHeader *header= nullptr);
     void get_part(std::byte *buffer, size_t frame_index);
     
     void read_header(DetectorHeader *header);
@@ -66,6 +67,8 @@ class RawSubFile {
     size_t pixels_per_frame() const { return m_rows * m_cols; }
     size_t bytes_per_pixel() const { return m_bitdepth / bits_per_byte; }
 
+    size_t frames_in_file() const { return m_num_frames; }
+
 private:
   template <typename T>
   void read_with_map(std::byte *image_buf);

include/aare/VarClusterFinder.hpp

@@ -226,7 +226,7 @@ template <typename T> void VarClusterFinder<T>::single_pass(NDView<T, 2> img) {
 
 template <typename T> void VarClusterFinder<T>::first_pass() {
 
-    for (size_t i = 0; i < original_.size(); ++i) {
+    for (ssize_t i = 0; i < original_.size(); ++i) {
         if (use_noise_map)
             threshold_ = 5 * noiseMap(i);
         binary_(i) = (original_(i) > threshold_);
@@ -250,7 +250,7 @@ template <typename T> void VarClusterFinder<T>::first_pass() {
 
 template <typename T> void VarClusterFinder<T>::second_pass() {
 
-    for (size_t i = 0; i != labeled_.size(); ++i) {
+    for (ssize_t i = 0; i != labeled_.size(); ++i) {
         auto cl = labeled_(i);
         if (cl != 0) {
             auto it = child.find(cl);

include/aare/algorithm.hpp

@@ -7,13 +7,20 @@
 
 namespace aare {
 /**
- * @brief Find the index of the last element smaller than val
- * assume a sorted array
+ * @brief Index of the last element that is smaller than val.
+ * Requires a sorted array. Uses >= for ordering. If all elements
+ * are smaller it returns the last element and if all elements are
+ * larger it returns the first element.
+ * @param first iterator to the first element
+ * @param last iterator to the last element
+ * @param val value to compare
+ * @return index of the last element that is smaller than val
+ * 
  */
 template <typename T>
 size_t last_smaller(const T* first, const T* last, T val) {
     for (auto iter = first+1; iter != last; ++iter) {
-        if (*iter > val) {
+        if (*iter >= val) {
             return std::distance(first, iter-1);
         }
     }
@@ -25,7 +32,49 @@ size_t last_smaller(const NDArray<T, 1>& arr, T val) {
     return last_smaller(arr.begin(), arr.end(), val);
 }
 
+template <typename T>
+size_t last_smaller(const std::vector<T>& vec, T val) {
+    return last_smaller(vec.data(), vec.data()+vec.size(), val);
+}
 
+/**
+ * @brief Index of the first element that is larger than val.
+ * Requires a sorted array. Uses > for ordering. If all elements
+ * are larger it returns the first element and if all elements are
+ * smaller it returns the last element.
+ * @param first iterator to the first element
+ * @param last iterator to the last element
+ * @param val value to compare
+ * @return index of the first element that is larger than val
+ */
+template <typename T>
+size_t first_larger(const T* first, const T* last, T val) {
+    for (auto iter = first; iter != last; ++iter) {
+        if (*iter > val) {
+            return std::distance(first, iter);
+        }
+    }
+    return std::distance(first, last-1);
+}
+
+template <typename T>
+size_t first_larger(const NDArray<T, 1>& arr, T val) {
+    return first_larger(arr.begin(), arr.end(), val);
+}
+
+template <typename T>
+size_t first_larger(const std::vector<T>& vec, T val) {
+    return first_larger(vec.data(), vec.data()+vec.size(), val);
+}
+
+/**
+ * @brief Index of the nearest element to val.
+ * Requires a sorted array. If there is no difference it takes the first element.
+ * @param first iterator to the first element
+ * @param last iterator to the last element
+ * @param val value to compare
+ * @return index of the nearest element
+ */
 template <typename T>
 size_t nearest_index(const T* first, const T* last, T val) {
     auto iter = std::min_element(first, last,
@@ -50,6 +99,13 @@ size_t nearest_index(const std::array<T,N>& arr, T val) {
     return nearest_index(arr.data(), arr.data()+arr.size(), val);
 }
 
+template <typename T>
+std::vector<T> cumsum(const std::vector<T>& vec) {
+    std::vector<T> result(vec.size());
+    std::partial_sum(vec.begin(), vec.end(), result.begin());
+    return result;
+}
+
 
 
 } // namespace aare

include/aare/decode.hpp

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <cstdint>
+#include <vector>
 #include <aare/NDView.hpp>
 namespace aare {
 
@@ -10,4 +11,16 @@ uint16_t adc_sar_04_decode64to16(uint64_t input);
 void adc_sar_05_decode64to16(NDView<uint64_t, 2> input, NDView<uint16_t,2> output);
 void adc_sar_04_decode64to16(NDView<uint64_t, 2> input, NDView<uint16_t,2> output);
 
-} // namespace aare
+
+/**
+ * @brief Apply custom weights to a 16-bit input value. Will sum up weights[i]**i
+ * for each bit i that is set in the input value.
+ * @throws std::out_of_range if weights.size() < 16
+ * @param input 16-bit input value
+ * @param weights vector of weights, size must be less than or equal to 16
+ */
+double apply_custom_weights(uint16_t input, const NDView<double, 1> weights);
+
+void apply_custom_weights(NDView<uint16_t, 1> input, NDView<double, 1> output, const NDView<double, 1> weights);
+
+} // namespace aare

include/aare/utils/ifstream_helpers.hpp

@@ -0,0 +1,12 @@
+#pragma once
+
+#include <fstream>
+#include <string>
+namespace aare {
+
+/** 
+ * @brief Get the error message from an ifstream object
+*/
+std::string ifstream_error_msg(std::ifstream &ifs);
+
+} // namespace aare

pyproject.toml

@@ -4,14 +4,32 @@ build-backend = "scikit_build_core.build"
 
 [project]
 name = "aare"
-version = "2025.4.1"
+version = "2025.4.22"
+requires-python = ">=3.11"
+dependencies = [
+    "numpy",
+    "matplotlib",
+]
+
+
+[tool.cibuildwheel]
+
+build = "cp{311,312,313}-manylinux_x86_64"
+
 
 
 
 [tool.scikit-build]
-cmake.verbose = true
+build.verbose = true
+cmake.build-type = "Release"
+install.components = ["python"]
 
 [tool.scikit-build.cmake.define]
 AARE_PYTHON_BINDINGS = "ON"
-AARE_SYSTEM_LIBRARIES = "ON"
-AARE_INSTALL_PYTHONEXT = "ON"
+AARE_INSTALL_PYTHONEXT = "ON"
+
+
+[tool.pytest.ini_options]
+markers = [
+    "files: marks tests that need additional data (deselect with '-m \"not files\"')",
+]

python/CMakeLists.txt

@@ -1,12 +1,13 @@
 
-find_package (Python 3.10 COMPONENTS Interpreter Development REQUIRED)
+find_package (Python 3.10 COMPONENTS Interpreter Development.Module REQUIRED)
+set(PYBIND11_FINDPYTHON ON) # Needed for RH8
 
 # Download or find pybind11 depending on configuration
 if(AARE_FETCH_PYBIND11)
     FetchContent_Declare(
         pybind11
         GIT_REPOSITORY https://github.com/pybind/pybind11
-        GIT_TAG        v2.13.0
+        GIT_TAG        v2.13.6
     )
     FetchContent_MakeAvailable(pybind11)
 else()
@@ -58,10 +59,16 @@ endforeach(FILE ${PYTHON_EXAMPLES})
 
 
 if(AARE_INSTALL_PYTHONEXT)
-    install(TARGETS _aare
+    install(
+        TARGETS _aare
         EXPORT "${TARGETS_EXPORT_NAME}"
         LIBRARY DESTINATION aare
+        COMPONENT python
     )
 
-    install(FILES ${PYTHON_FILES} DESTINATION aare)
+    install(
+        FILES ${PYTHON_FILES} 
+        DESTINATION aare
+        COMPONENT python
+        )
 endif()

python/aare/__init__.py

@@ -2,22 +2,26 @@
 from . import _aare
 
 
-# from ._aare import File, RawMasterFile, RawSubFile
-# from ._aare import Pedestal_d, Pedestal_f, ClusterFinder, VarClusterFinder
+from ._aare import File, RawMasterFile, RawSubFile, JungfrauDataFile
+from ._aare import Pedestal_d, Pedestal_f, ClusterFinder, VarClusterFinder
 from ._aare import DetectorType
-from ._aare import ClusterFile_Cluster3x3i as ClusterFile
+from ._aare import ClusterFile
 from ._aare import hitmap
 from ._aare import ROI
 
-# from ._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink, ClusterVector_i
+from ._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink, ClusterVector_i
 
 from ._aare import fit_gaus, fit_pol1
 from ._aare import Interpolator
+
+
+from ._aare import apply_custom_weights
+
 from .CtbRawFile import CtbRawFile
 from .RawFile import RawFile
 from .ScanParameters import ScanParameters
 
-from .utils import random_pixels, random_pixel, flat_list
+from .utils import random_pixels, random_pixel, flat_list, add_colorbar
 
 
 #make functions available in the top level API

python/aare/utils.py

@@ -1,4 +1,6 @@
 import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import make_axes_locatable
 
 def random_pixels(n_pixels, xmin=0, xmax=512, ymin=0, ymax=1024):
     """Return a list of random pixels.
@@ -24,4 +26,11 @@ def random_pixel(xmin=0, xmax=512, ymin=0, ymax=1024):
 
 def flat_list(xss):
     """Flatten a list of lists."""
-    return [x for xs in xss for x in xs]
+    return [x for xs in xss for x in xs]
+
+def add_colorbar(ax, im, size="5%", pad=0.05):
+    """Add a colorbar with the same height as the image."""
+    divider = make_axes_locatable(ax)
+    cax = divider.append_axes("right", size=size, pad=pad)
+    plt.colorbar(im, cax=cax)
+    return ax, im, cax

python/src/cluster.hpp

@@ -16,246 +16,149 @@
 namespace py = pybind11;
 using pd_type = double;
 
-using namespace aare;
-
-template <typename Type, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType>
-void define_cluster(py::module &m, const std::string &typestr) {
-    auto class_name = fmt::format("Cluster{}", typestr);
-
-    using ClusterType =
-        Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType, void>;
-    py::class_<Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType, void>>(
-        m, class_name.c_str())
-
-        .def(py::init([](uint8_t x, uint8_t y, py::array_t<Type> data) {
-            py::buffer_info buf_info = data.request();
-            Type *ptr = static_cast<Type *>(buf_info.ptr);
-            Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType, void> cluster;
-            cluster.x = x;
-            cluster.y = y;
-            std::copy(ptr, ptr + ClusterSizeX * ClusterSizeY,
-                      cluster.data); // Copy array contents
-            return cluster;
-        }))
-
-        //.def(py::init<>())
-        .def_readwrite("x", &ClusterType::x)
-        .def_readwrite("y", &ClusterType::y)
-        .def_property(
-            "data",
-            [](ClusterType &c) -> py::array {
-                return py::array(py::buffer_info(
-                    c.data, sizeof(Type),
-                    py::format_descriptor<Type>::format(), // Type
-                                                           // format
-                    1, // Number of dimensions
-                    {static_cast<ssize_t>(ClusterSizeX *
-                                          ClusterSizeY)}, // Shape (flattened)
-                    {sizeof(Type)} // Stride (step size between elements)
-                    ));
-            },
-            [](ClusterType &c, py::array_t<Type> arr) {
-                py::buffer_info buf_info = arr.request();
-                Type *ptr = static_cast<Type *>(buf_info.ptr);
-                std::copy(ptr, ptr + ClusterSizeX * ClusterSizeY,
-                          c.data); // TODO dont iterate over centers!!!
-            });
-}
-
-template <typename Type, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType = uint16_t>
+template <typename T>
 void define_cluster_vector(py::module &m, const std::string &typestr) {
-    using ClusterType =
-        Cluster<Type, ClusterSizeX, ClusterSizeY, uint16_t, void>;
     auto class_name = fmt::format("ClusterVector_{}", typestr);
-
-    py::class_<ClusterVector<ClusterType>>(m, class_name.c_str(),
-                                           py::buffer_protocol())
-
-        .def(py::init()) // TODO change!!!
-        /*
+    py::class_<ClusterVector<T>>(m, class_name.c_str(), py::buffer_protocol())
+        .def(py::init<int, int>(),
+    py::arg("cluster_size_x") = 3, py::arg("cluster_size_y") = 3)
         .def("push_back",
-             [](ClusterVector<ClusterType> &self, ClusterType &cl) {
-                 //  auto view = make_view_2d(data);
-                 self.push_back(cl);
+             [](ClusterVector<T> &self, int x, int y, py::array_t<T> data) {
+                //  auto view = make_view_2d(data);
+                 self.push_back(x, y, reinterpret_cast<const std::byte*>(data.data()));
              })
-        */
-        /*
-        .def(
-            "push_back",
-            [](ClusterVector<ClusterType> &self, py::object obj) {
-                ClusterType &cl = py::cast<ClusterType &>(obj);
-                self.push_back(cl);
-            },
-            py::arg("cluster"))
-        */
-
-        .def("push_back",
-             [](ClusterVector<ClusterType> &self, const ClusterType &cluster) {
-                 self.push_back(cluster);
-             })
-
-        //.def("push_back", &ClusterVector<ClusterType>::push_back) //TODO
-        // implement push_back
-        .def_property_readonly("size", &ClusterVector<ClusterType>::size)
-        .def("item_size", &ClusterVector<ClusterType>::item_size)
+        .def_property_readonly("size", &ClusterVector<T>::size)
+        .def("item_size", &ClusterVector<T>::item_size)
         .def_property_readonly("fmt",
-                               [typestr]() { return fmt_format<ClusterType>; })
-        /*
+                               [typestr](ClusterVector<T> &self) {
+                                   return fmt::format(
+                                       self.fmt_base(), self.cluster_size_x(),
+                                       self.cluster_size_y(), typestr);
+                               })
         .def("sum",
-             [](ClusterVector<ClusterType> &self) {
+             [](ClusterVector<T> &self) {
                  auto *vec = new std::vector<T>(self.sum());
                  return return_vector(vec);
              })
-        .def("sum_2x2",
-             [](ClusterVector<ClusterType> &self) {
-                 auto *vec = new std::vector<T>(self.sum_2x2());
-                 return return_vector(vec);
-             })
-        */
-        .def_property_readonly("cluster_size_x",
-                               &ClusterVector<ClusterType>::cluster_size_x)
-        .def_property_readonly("cluster_size_y",
-                               &ClusterVector<ClusterType>::cluster_size_y)
-        .def_property_readonly("capacity",
-                               &ClusterVector<ClusterType>::capacity)
-        .def_property("frame_number", &ClusterVector<ClusterType>::frame_number,
-                      &ClusterVector<ClusterType>::set_frame_number)
-        .def_buffer(
-            [typestr](ClusterVector<ClusterType> &self) -> py::buffer_info {
-                return py::buffer_info(
-                    self.data(),             /* Pointer to buffer */
-                    self.item_size(),        /* Size of one scalar */
-                    fmt_format<ClusterType>, /* Format descriptor */
-                    1,                       /* Number of dimensions */
-                    {self.size()},           /* Buffer dimensions */
-                    {self.item_size()} /* Strides (in bytes) for each index */
-                );
-            });
+        .def("sum_2x2", [](ClusterVector<T> &self) {
+            auto *vec = new std::vector<T>(self.sum_2x2());
+            return return_vector(vec);
+        })
+        .def_property_readonly("cluster_size_x", &ClusterVector<T>::cluster_size_x)
+        .def_property_readonly("cluster_size_y", &ClusterVector<T>::cluster_size_y)
+        .def_property_readonly("capacity", &ClusterVector<T>::capacity)
+        .def_property("frame_number", &ClusterVector<T>::frame_number,
+                      &ClusterVector<T>::set_frame_number)
+        .def_buffer([typestr](ClusterVector<T> &self) -> py::buffer_info {
+            return py::buffer_info(
+                self.data(),      /* Pointer to buffer */
+                self.item_size(), /* Size of one scalar */
+                fmt::format(self.fmt_base(), self.cluster_size_x(),
+                            self.cluster_size_y(),
+                            typestr), /* Format descriptor */
+                1,                    /* Number of dimensions */
+                {self.size()},        /* Buffer dimensions */
+                {self.item_size()}    /* Strides (in bytes) for each index */
+            );
+        });
 }
 
-template <typename ClusterType>
-void define_cluster_finder_mt_bindings(py::module &m,
-                                       const std::string &typestr) {
-    auto class_name = fmt::format("ClusterFinderMT_{}", typestr);
-
-    py::class_<ClusterFinderMT<ClusterType, uint16_t, pd_type>>(
-        m, class_name.c_str())
-        .def(py::init<Shape<2>, pd_type, size_t, size_t>(),
-             py::arg("image_size"), py::arg("n_sigma") = 5.0,
-             py::arg("capacity") = 2048, py::arg("n_threads") = 3)
+void define_cluster_finder_mt_bindings(py::module &m) {
+    py::class_<ClusterFinderMT<uint16_t, pd_type>>(m, "ClusterFinderMT")
+        .def(py::init<Shape<2>, Shape<2>, pd_type, size_t, size_t>(),
+             py::arg("image_size"), py::arg("cluster_size"),
+             py::arg("n_sigma") = 5.0, py::arg("capacity") = 2048,
+             py::arg("n_threads") = 3)
         .def("push_pedestal_frame",
-             [](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
+             [](ClusterFinderMT<uint16_t, pd_type> &self,
                 py::array_t<uint16_t> frame) {
                  auto view = make_view_2d(frame);
                  self.push_pedestal_frame(view);
              })
         .def(
             "find_clusters",
-            [](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
+            [](ClusterFinderMT<uint16_t, pd_type> &self,
                py::array_t<uint16_t> frame, uint64_t frame_number) {
                 auto view = make_view_2d(frame);
                 self.find_clusters(view, frame_number);
                 return;
             },
             py::arg(), py::arg("frame_number") = 0)
-        .def("clear_pedestal",
-             &ClusterFinderMT<ClusterType, uint16_t, pd_type>::clear_pedestal)
-        .def("sync", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::sync)
-        .def("stop", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::stop)
-        .def("start", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::start)
-        .def(
-            "pedestal",
-            [](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
-               size_t thread_index) {
-                auto pd = new NDArray<pd_type, 2>{};
-                *pd = self.pedestal(thread_index);
-                return return_image_data(pd);
-            },
-            py::arg("thread_index") = 0)
-        .def(
-            "noise",
-            [](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
-               size_t thread_index) {
-                auto arr = new NDArray<pd_type, 2>{};
-                *arr = self.noise(thread_index);
-                return return_image_data(arr);
-            },
-            py::arg("thread_index") = 0);
+        .def("clear_pedestal", &ClusterFinderMT<uint16_t, pd_type>::clear_pedestal)
+        .def("sync", &ClusterFinderMT<uint16_t, pd_type>::sync)
+        .def("stop", &ClusterFinderMT<uint16_t, pd_type>::stop)
+        .def("start", &ClusterFinderMT<uint16_t, pd_type>::start)
+        .def("pedestal",
+             [](ClusterFinderMT<uint16_t, pd_type> &self, size_t thread_index) {
+                 auto pd = new NDArray<pd_type, 2>{};
+                 *pd = self.pedestal(thread_index);
+                 return return_image_data(pd);
+             },py::arg("thread_index") = 0)
+        .def("noise",
+             [](ClusterFinderMT<uint16_t, pd_type> &self, size_t thread_index) {
+                 auto arr = new NDArray<pd_type, 2>{};
+                 *arr = self.noise(thread_index);
+                 return return_image_data(arr);
+             },py::arg("thread_index") = 0);
 }
 
-template <typename ClusterType>
-void define_cluster_collector_bindings(py::module &m,
-                                       const std::string &typestr) {
-    auto class_name = fmt::format("ClusterCollector_{}", typestr);
-
-    py::class_<ClusterCollector<ClusterType>>(m, class_name.c_str())
-        .def(py::init<ClusterFinderMT<ClusterType, uint16_t, double> *>())
-        .def("stop", &ClusterCollector<ClusterType>::stop)
+void define_cluster_collector_bindings(py::module &m) {
+    py::class_<ClusterCollector>(m, "ClusterCollector")
+        .def(py::init<ClusterFinderMT<uint16_t, double, int32_t> *>())
+        .def("stop", &ClusterCollector::stop)
         .def(
             "steal_clusters",
-            [](ClusterCollector<ClusterType> &self) {
-                auto v = new std::vector<ClusterVector<ClusterType>>(
-                    self.steal_clusters());
+            [](ClusterCollector &self) {
+                auto v =
+                    new std::vector<ClusterVector<int>>(self.steal_clusters());
                 return v;
             },
             py::return_value_policy::take_ownership);
 }
 
-template <typename ClusterType>
-void define_cluster_file_sink_bindings(py::module &m,
-                                       const std::string &typestr) {
-    auto class_name = fmt::format("ClusterFileSink_{}", typestr);
-
-    py::class_<ClusterFileSink<ClusterType>>(m, class_name.c_str())
-        .def(py::init<ClusterFinderMT<ClusterType, uint16_t, double> *,
+void define_cluster_file_sink_bindings(py::module &m) {
+    py::class_<ClusterFileSink>(m, "ClusterFileSink")
+        .def(py::init<ClusterFinderMT<uint16_t, double, int32_t> *,
                       const std::filesystem::path &>())
-        .def("stop", &ClusterFileSink<ClusterType>::stop);
+        .def("stop", &ClusterFileSink::stop);
 }
 
-template <typename ClusterType>
-void define_cluster_finder_bindings(py::module &m, const std::string &typestr) {
-    auto class_name = fmt::format("ClusterFinder_{}", typestr);
-
-    py::class_<ClusterFinder<ClusterType, uint16_t, pd_type>>(
-        m, class_name.c_str())
-        .def(py::init<Shape<2>, pd_type, size_t>(), py::arg("image_size"),
+void define_cluster_finder_bindings(py::module &m) {
+    py::class_<ClusterFinder<uint16_t, pd_type>>(m, "ClusterFinder")
+        .def(py::init<Shape<2>, Shape<2>, pd_type, size_t>(),
+             py::arg("image_size"), py::arg("cluster_size"),
              py::arg("n_sigma") = 5.0, py::arg("capacity") = 1'000'000)
         .def("push_pedestal_frame",
-             [](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
+             [](ClusterFinder<uint16_t, pd_type> &self,
                 py::array_t<uint16_t> frame) {
                  auto view = make_view_2d(frame);
                  self.push_pedestal_frame(view);
              })
-        .def("clear_pedestal",
-             &ClusterFinder<ClusterType, uint16_t, pd_type>::clear_pedestal)
-        .def_property_readonly(
-            "pedestal",
-            [](ClusterFinder<ClusterType, uint16_t, pd_type> &self) {
-                auto pd = new NDArray<pd_type, 2>{};
-                *pd = self.pedestal();
-                return return_image_data(pd);
-            })
-        .def_property_readonly(
-            "noise",
-            [](ClusterFinder<ClusterType, uint16_t, pd_type> &self) {
-                auto arr = new NDArray<pd_type, 2>{};
-                *arr = self.noise();
-                return return_image_data(arr);
-            })
+        .def("clear_pedestal", &ClusterFinder<uint16_t, pd_type>::clear_pedestal)
+        .def_property_readonly("pedestal",
+                               [](ClusterFinder<uint16_t, pd_type> &self) {
+                                   auto pd = new NDArray<pd_type, 2>{};
+                                   *pd = self.pedestal();
+                                   return return_image_data(pd);
+                               })
+        .def_property_readonly("noise",
+                               [](ClusterFinder<uint16_t, pd_type> &self) {
+                                   auto arr = new NDArray<pd_type, 2>{};
+                                   *arr = self.noise();
+                                   return return_image_data(arr);
+                               })
         .def(
             "steal_clusters",
-            [](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
+            [](ClusterFinder<uint16_t, pd_type> &self,
                bool realloc_same_capacity) {
-                auto v = new ClusterVector<ClusterType>(
+                auto v = new ClusterVector<int>(
                     self.steal_clusters(realloc_same_capacity));
                 return v;
             },
             py::arg("realloc_same_capacity") = false)
         .def(
             "find_clusters",
-            [](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
+            [](ClusterFinder<uint16_t, pd_type> &self,
                py::array_t<uint16_t> frame, uint64_t frame_number) {
                 auto view = make_view_2d(frame);
                 self.find_clusters(view, frame_number);
@@ -264,7 +167,7 @@ void define_cluster_finder_bindings(py::module &m, const std::string &typestr) {
             py::arg(), py::arg("frame_number") = 0);
 
     m.def("hitmap",
-          [](std::array<size_t, 2> image_size, ClusterVector<ClusterType> &cv) {
+          [](std::array<size_t, 2> image_size, ClusterVector<int32_t> &cv) {
               py::array_t<int32_t> hitmap(image_size);
               auto r = hitmap.mutable_unchecked<2>();
 
@@ -283,4 +186,24 @@ void define_cluster_finder_bindings(py::module &m, const std::string &typestr) {
               }
               return hitmap;
           });
-}
+    define_cluster_vector<int>(m, "i");
+    define_cluster_vector<double>(m, "d");
+    define_cluster_vector<float>(m, "f");
+
+    py::class_<DynamicCluster>(m, "DynamicCluster", py::buffer_protocol())
+        .def(py::init<int, int, Dtype>())
+        .def("size", &DynamicCluster::size)
+        .def("begin", &DynamicCluster::begin)
+        .def("end", &DynamicCluster::end)
+        .def_readwrite("x", &DynamicCluster::x)
+        .def_readwrite("y", &DynamicCluster::y)
+        .def_buffer([](DynamicCluster &c) -> py::buffer_info {
+            return py::buffer_info(c.data(), c.dt.bytes(), c.dt.format_descr(),
+                                   1, {c.size()}, {c.dt.bytes()});
+        })
+
+        .def("__repr__", [](const DynamicCluster &a) {
+            return "<DynamicCluster: x: " + std::to_string(a.x) +
+                   ", y: " + std::to_string(a.y) + ">";
+        });
+}

python/src/cluster_file.hpp

@@ -1,4 +1,3 @@
-#include "aare/CalculateEta.hpp"
 #include "aare/ClusterFile.hpp"
 #include "aare/defs.hpp"
 
@@ -11,81 +10,64 @@
 #include <pybind11/stl/filesystem.h>
 #include <string>
 
-// Disable warnings for unused parameters, as we ignore some
-// in the __exit__ method
+//Disable warnings for unused parameters, as we ignore some
+//in the __exit__ method
 #pragma GCC diagnostic push
 #pragma GCC diagnostic ignored "-Wunused-parameter"
 
+
 namespace py = pybind11;
 using namespace ::aare;
 
-template <typename ClusterType>
-void define_cluster_file_io_bindings(py::module &m,
-                                     const std::string &typestr) {
-    // PYBIND11_NUMPY_DTYPE(Cluster<int, 3, 3>, x, y,
-    // data); // is this used - maybe use as cluster type
+void define_cluster_file_io_bindings(py::module &m) {
+    PYBIND11_NUMPY_DTYPE(Cluster3x3, x, y, data);
 
-    auto class_name = fmt::format("ClusterFile_{}", typestr);
-
-    py::class_<ClusterFile<ClusterType>>(m, class_name.c_str())
+    py::class_<ClusterFile>(m, "ClusterFile")
         .def(py::init<const std::filesystem::path &, size_t,
                       const std::string &>(),
              py::arg(), py::arg("chunk_size") = 1000, py::arg("mode") = "r")
-        .def(
-            "read_clusters",
-            [](ClusterFile<ClusterType> &self, size_t n_clusters) {
-                auto v = new ClusterVector<ClusterType>(
-                    self.read_clusters(n_clusters));
+        .def("read_clusters",
+             [](ClusterFile &self, size_t n_clusters) {
+                auto v = new ClusterVector<int32_t>(self.read_clusters(n_clusters));
                 return v;
-            },
-            py::return_value_policy::take_ownership)
+             },py::return_value_policy::take_ownership)
         .def("read_frame",
-             [](ClusterFile<ClusterType> &self) {
-                 auto v = new ClusterVector<ClusterType>(self.read_frame());
-                 return v;
+             [](ClusterFile &self) {
+                auto v = new ClusterVector<int32_t>(self.read_frame());
+                return v;
              })
-        .def("set_roi", &ClusterFile<ClusterType>::set_roi)
-        .def(
-            "set_noise_map",
-            [](ClusterFile<ClusterType> &self, py::array_t<int32_t> noise_map) {
-                auto view = make_view_2d(noise_map);
-                self.set_noise_map(view);
-            })
-
-        .def("set_gain_map",
-             [](ClusterFile<ClusterType> &self, py::array_t<double> gain_map) {
-                 auto view = make_view_2d(gain_map);
-                 self.set_gain_map(view);
-             })
-
-        // void set_gain_map(const GainMap &gain_map); //TODO do i need a
-        // gainmap constructor?
-
-        .def("close", &ClusterFile<ClusterType>::close)
-        .def("write_frame", &ClusterFile<ClusterType>::write_frame)
-        .def("__enter__", [](ClusterFile<ClusterType> &self) { return &self; })
+        .def("set_roi", &ClusterFile::set_roi)
+        .def("set_noise_map", [](ClusterFile &self, py::array_t<int32_t> noise_map) {
+            auto view = make_view_2d(noise_map);
+            self.set_noise_map(view);
+        })
+        .def("set_gain_map", [](ClusterFile &self, py::array_t<double> gain_map) {
+            auto view = make_view_2d(gain_map);
+            self.set_gain_map(view);
+        })
+        .def("close", &ClusterFile::close)
+        .def("write_frame", &ClusterFile::write_frame)
+        .def("__enter__", [](ClusterFile &self) { return &self; })
         .def("__exit__",
-             [](ClusterFile<ClusterType> &self,
+             [](ClusterFile &self,
                 const std::optional<pybind11::type> &exc_type,
                 const std::optional<pybind11::object> &exc_value,
                 const std::optional<pybind11::object> &traceback) {
                  self.close();
              })
-        .def("__iter__", [](ClusterFile<ClusterType> &self) { return &self; })
-        .def("__next__", [](ClusterFile<ClusterType> &self) {
-            auto v = new ClusterVector<ClusterType>(
-                self.read_clusters(self.chunk_size()));
+        .def("__iter__", [](ClusterFile &self) { return &self; })
+        .def("__next__", [](ClusterFile &self) {
+            auto v = new ClusterVector<int32_t>(self.read_clusters(self.chunk_size()));
             if (v->size() == 0) {
                 throw py::stop_iteration();
             }
             return v;
         });
 
-    m.def("calculate_eta2",
-          [](const aare::ClusterVector<ClusterType> &clusters) {
-              auto eta2 = new NDArray<double, 2>(calculate_eta2(clusters));
-              return return_image_data(eta2);
-          });
+        m.def("calculate_eta2", []( aare::ClusterVector<int32_t> &clusters) {
+            auto eta2 = new NDArray<double, 2>(calculate_eta2(clusters));
+            return return_image_data(eta2);
+        });
 }
 
 #pragma GCC diagnostic pop

python/src/ctb_raw_file.hpp

@@ -10,6 +10,8 @@
 #include "aare/decode.hpp"
 // #include "aare/fClusterFileV2.hpp"
 
+#include "np_helper.hpp"
+
 #include <cstdint>
 #include <filesystem>
 #include <pybind11/iostream.h>
@@ -65,35 +67,54 @@ m.def("adc_sar_04_decode64to16", [](py::array_t<uint8_t> input) {
     return output;
 });
 
- py::class_<CtbRawFile>(m, "CtbRawFile")
-        .def(py::init<const std::filesystem::path &>())
-        .def("read_frame",
-             [](CtbRawFile &self) {
-                 size_t image_size = self.image_size_in_bytes();
-                 py::array image;
-                 std::vector<ssize_t> shape;
-                 shape.reserve(2);
-                 shape.push_back(1);
-                 shape.push_back(image_size);
+m.def(
+    "apply_custom_weights",
+    [](py::array_t<uint16_t, py::array::c_style | py::array::forcecast> &input,
+       py::array_t<double, py::array::c_style | py::array::forcecast>
+           &weights) {
+        
 
-                 py::array_t<DetectorHeader> header(1);
+        // Create new array with same shape as the input array (uninitialized values)
+        py::buffer_info buf = input.request();
+        py::array_t<double> output(buf.shape);
 
-                 // always read bytes
-                 image = py::array_t<uint8_t>(shape);
+        // Use NDViews to call into the C++ library
+        auto weights_view = make_view_1d(weights);
+        NDView<uint16_t, 1> input_view(input.mutable_data(), {input.size()});
+        NDView<double, 1> output_view(output.mutable_data(), {output.size()});
 
-                 self.read_into(
-                     reinterpret_cast<std::byte *>(image.mutable_data()),
-                     header.mutable_data());
+        apply_custom_weights(input_view, output_view, weights_view);
+        return output;
+    });
 
-                 return py::make_tuple(header, image);
-             })
-        .def("seek", &CtbRawFile::seek)
-        .def("tell", &CtbRawFile::tell)
-        .def("master", &CtbRawFile::master)
+py::class_<CtbRawFile>(m, "CtbRawFile")
+    .def(py::init<const std::filesystem::path &>())
+    .def("read_frame",
+         [](CtbRawFile &self) {
+             size_t image_size = self.image_size_in_bytes();
+             py::array image;
+             std::vector<ssize_t> shape;
+             shape.reserve(2);
+             shape.push_back(1);
+             shape.push_back(image_size);
 
-        .def_property_readonly("image_size_in_bytes",
-                               &CtbRawFile::image_size_in_bytes)
+             py::array_t<DetectorHeader> header(1);
 
-        .def_property_readonly("frames_in_file", &CtbRawFile::frames_in_file);
+             // always read bytes
+             image = py::array_t<uint8_t>(shape);
 
-}
+             self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()),
+                            header.mutable_data());
+
+             return py::make_tuple(header, image);
+         })
+    .def("seek", &CtbRawFile::seek)
+    .def("tell", &CtbRawFile::tell)
+    .def("master", &CtbRawFile::master)
+
+    .def_property_readonly("image_size_in_bytes",
+                           &CtbRawFile::image_size_in_bytes)
+
+    .def_property_readonly("frames_in_file", &CtbRawFile::frames_in_file);
+
+}

python/src/file.hpp

@@ -20,6 +20,9 @@
 namespace py = pybind11;
 using namespace ::aare;
 
+
+
+
 //Disable warnings for unused parameters, as we ignore some
 //in the __exit__ method
 #pragma GCC diagnostic push
@@ -214,36 +217,9 @@ void define_file_io_bindings(py::module &m) {
     
 
 
-    py::class_<RawSubFile>(m, "RawSubFile")
-        .def(py::init<const std::filesystem::path &, DetectorType, size_t,
-                      size_t, size_t>())
-        .def_property_readonly("bytes_per_frame", &RawSubFile::bytes_per_frame)
-        .def_property_readonly("pixels_per_frame",
-                               &RawSubFile::pixels_per_frame)
-        .def("seek", &RawSubFile::seek)
-        .def("tell", &RawSubFile::tell)
-        .def_property_readonly("rows", &RawSubFile::rows)
-        .def_property_readonly("cols", &RawSubFile::cols)
-        .def("read_frame",
-             [](RawSubFile &self) {
-                 const uint8_t item_size = self.bytes_per_pixel();
-                 py::array image;
-                 std::vector<ssize_t> shape;
-                 shape.reserve(2);
-                 shape.push_back(self.rows());
-                 shape.push_back(self.cols());
-                 if (item_size == 1) {
-                     image = py::array_t<uint8_t>(shape);
-                 } else if (item_size == 2) {
-                     image = py::array_t<uint16_t>(shape);
-                 } else if (item_size == 4) {
-                     image = py::array_t<uint32_t>(shape);
-                 }
-                 fmt::print("item_size: {} rows: {} cols: {}\n", item_size, self.rows(), self.cols());
-                 self.read_into(
-                     reinterpret_cast<std::byte *>(image.mutable_data()));
-                 return image;
-             });
+
+
+
 
 #pragma GCC diagnostic pop
     // py::class_<ClusterHeader>(m, "ClusterHeader")

python/src/interpolation.hpp

@@ -8,54 +8,31 @@
 #include <pybind11/stl.h>
 
 namespace py = pybind11;
-
-template <typename ClusterType>
-void register_interpolate(py::class_<aare::Interpolator> &interpolator,
-                          const std::string &typestr) {
-    auto name = fmt::format("interpolate_{}", typestr);
-
-    interpolator.def(name.c_str(),
-                     [](aare::Interpolator &self,
-                        const ClusterVector<ClusterType> &clusters) {
-                         auto photons = self.interpolate<ClusterType>(clusters);
-                         auto *ptr = new std::vector<Photon>{photons};
-                         return return_vector(ptr);
-                     });
-}
-
 void define_interpolation_bindings(py::module &m) {
 
-    PYBIND11_NUMPY_DTYPE(aare::Photon, x, y, energy);
+    PYBIND11_NUMPY_DTYPE(aare::Photon, x,y,energy);
 
-    auto interpolator =
-        py::class_<aare::Interpolator>(m, "Interpolator")
-            .def(py::init([](py::array_t<double, py::array::c_style |
-                                                     py::array::forcecast>
-                                 etacube,
-                             py::array_t<double> xbins,
-                             py::array_t<double> ybins,
-                             py::array_t<double> ebins) {
-                return Interpolator(make_view_3d(etacube), make_view_1d(xbins),
-                                    make_view_1d(ybins), make_view_1d(ebins));
-            }))
-            .def("get_ietax",
-                 [](Interpolator &self) {
-                     auto *ptr = new NDArray<double, 3>{};
-                     *ptr = self.get_ietax();
-                     return return_image_data(ptr);
-                 })
-            .def("get_ietay", [](Interpolator &self) {
-                auto *ptr = new NDArray<double, 3>{};
-                *ptr = self.get_ietay();
-                return return_image_data(ptr);
-            });
-
-    register_interpolate<Cluster<int, 3, 3>>(interpolator, "Cluster3x3i");
-    register_interpolate<Cluster<float, 3, 3>>(interpolator, "Cluster3x3f");
-    register_interpolate<Cluster<double, 3, 3>>(interpolator, "Cluster3x3d");
-    register_interpolate<Cluster<int, 2, 2>>(interpolator, "Cluster2x2i");
-    register_interpolate<Cluster<float, 2, 2>>(interpolator, "Cluster2x2f");
-    register_interpolate<Cluster<double, 2, 2>>(interpolator, "Cluster2x2d");
+    py::class_<aare::Interpolator>(m, "Interpolator")
+        .def(py::init([](py::array_t<double, py::array::c_style | py::array::forcecast> etacube, py::array_t<double> xbins,
+                         py::array_t<double> ybins, py::array_t<double> ebins) {
+            return Interpolator(make_view_3d(etacube), make_view_1d(xbins),
+                                make_view_1d(ybins), make_view_1d(ebins));
+        }))
+        .def("get_ietax", [](Interpolator& self){
+            auto*ptr = new NDArray<double,3>{};
+            *ptr = self.get_ietax();
+            return return_image_data(ptr);
+        })
+        .def("get_ietay", [](Interpolator& self){
+            auto*ptr = new NDArray<double,3>{};
+            *ptr = self.get_ietay();
+            return return_image_data(ptr);
+        })
+        .def("interpolate", [](Interpolator& self, const ClusterVector<int32_t>& clusters){
+            auto photons = self.interpolate(clusters);
+            auto* ptr = new std::vector<Photon>{photons};
+            return return_vector(ptr);
+        });
 
     // TODO! Evaluate without converting to double
     m.def(

python/src/jungfrau_data_file.hpp

@@ -0,0 +1,116 @@
+
+#include "aare/JungfrauDataFile.hpp"
+#include "aare/defs.hpp"
+
+#include <cstdint>
+#include <filesystem>
+#include <pybind11/iostream.h>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/stl/filesystem.h>
+#include <string>
+
+namespace py = pybind11;
+using namespace ::aare;
+
+// Disable warnings for unused parameters, as we ignore some
+// in the __exit__ method
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wunused-parameter"
+
+auto read_dat_frame(JungfrauDataFile &self) {
+    py::array_t<JungfrauDataHeader> header(1);
+    py::array_t<uint16_t> image({
+        self.rows(), 
+        self.cols() 
+    });
+
+    self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()),
+                   header.mutable_data());
+
+    return py::make_tuple(header, image);
+}
+
+auto read_n_dat_frames(JungfrauDataFile &self, size_t n_frames) {
+    // adjust for actual frames left in the file
+    n_frames = std::min(n_frames, self.total_frames() - self.tell());
+    if (n_frames == 0) {
+        throw std::runtime_error("No frames left in file");
+    }
+
+    py::array_t<JungfrauDataHeader> header(n_frames);
+    py::array_t<uint16_t> image({
+        n_frames, self.rows(), 
+        self.cols()});
+
+    self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()),
+                   n_frames, header.mutable_data());
+
+    return py::make_tuple(header, image);
+}
+
+void define_jungfrau_data_file_io_bindings(py::module &m) {
+    // Make the JungfrauDataHeader usable from numpy
+    PYBIND11_NUMPY_DTYPE(JungfrauDataHeader, framenum, bunchid);
+
+    py::class_<JungfrauDataFile>(m, "JungfrauDataFile")
+        .def(py::init<const std::filesystem::path &>())
+        .def("seek", &JungfrauDataFile::seek,
+             R"(
+               Seek to the given frame index.
+               )")
+        .def("tell", &JungfrauDataFile::tell,
+             R"(
+               Get the current frame index.
+               )")
+        .def_property_readonly("rows", &JungfrauDataFile::rows)
+        .def_property_readonly("cols", &JungfrauDataFile::cols)
+        .def_property_readonly("base_name", &JungfrauDataFile::base_name)
+        .def_property_readonly("bytes_per_frame",
+                               &JungfrauDataFile::bytes_per_frame)
+        .def_property_readonly("pixels_per_frame",
+                               &JungfrauDataFile::pixels_per_frame)
+        .def_property_readonly("bytes_per_pixel",
+                               &JungfrauDataFile::bytes_per_pixel)
+        .def_property_readonly("bitdepth", &JungfrauDataFile::bitdepth)
+        .def_property_readonly("current_file", &JungfrauDataFile::current_file)
+        .def_property_readonly("total_frames", &JungfrauDataFile::total_frames)
+        .def_property_readonly("n_files", &JungfrauDataFile::n_files)
+        .def("read_frame", &read_dat_frame,
+             R"(
+               Read a single frame from the file.
+               )")
+        .def("read_n", &read_n_dat_frames,
+             R"(
+               Read maximum n_frames frames from the file.
+               )")
+        .def(
+            "read",
+            [](JungfrauDataFile &self) {
+                self.seek(0);
+                auto n_frames = self.total_frames();
+                return read_n_dat_frames(self, n_frames);
+            },
+            R"(
+              Read all frames from the file. Seeks to the beginning before reading.
+              )")
+        .def("__enter__", [](JungfrauDataFile &self) { return &self; })
+        .def("__exit__",
+             [](JungfrauDataFile &self,
+                const std::optional<pybind11::type> &exc_type,
+                const std::optional<pybind11::object> &exc_value,
+                const std::optional<pybind11::object> &traceback) {
+                 //  self.close();
+             })
+        .def("__iter__", [](JungfrauDataFile &self) { return &self; })
+        .def("__next__", [](JungfrauDataFile &self) {
+            try {
+                return read_dat_frame(self);
+            } catch (std::runtime_error &e) {
+                throw py::stop_iteration();
+            }
+        });
+}
+
+#pragma GCC diagnostic pop

python/src/module.cpp

@@ -1,17 +1,20 @@
-// Files with bindings to the different classes
-#include "cluster.hpp"
-#include "cluster_file.hpp"
-#include "ctb_raw_file.hpp"
+//Files with bindings to the different classes
 #include "file.hpp"
-#include "fit.hpp"
-#include "interpolation.hpp"
-#include "pedestal.hpp"
-#include "pixel_map.hpp"
 #include "raw_file.hpp"
+#include "ctb_raw_file.hpp"
 #include "raw_master_file.hpp"
 #include "var_cluster.hpp"
+#include "pixel_map.hpp"
+#include "pedestal.hpp"
+#include "cluster.hpp"
+#include "cluster_file.hpp"
+#include "fit.hpp"
+#include "interpolation.hpp"
+#include "raw_sub_file.hpp"
 
-// Pybind stuff
+#include "jungfrau_data_file.hpp"
+
+//Pybind stuff
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 
@@ -20,61 +23,20 @@ namespace py = pybind11;
 PYBIND11_MODULE(_aare, m) {
     define_file_io_bindings(m);
     define_raw_file_io_bindings(m);
+    define_raw_sub_file_io_bindings(m);
     define_ctb_raw_file_io_bindings(m);
     define_raw_master_file_bindings(m);
     define_var_cluster_finder_bindings(m);
     define_pixel_map_bindings(m);
     define_pedestal_bindings<double>(m, "Pedestal_d");
     define_pedestal_bindings<float>(m, "Pedestal_f");
+    define_cluster_finder_bindings(m);
+    define_cluster_finder_mt_bindings(m);
+    define_cluster_file_io_bindings(m);
+    define_cluster_collector_bindings(m);
+    define_cluster_file_sink_bindings(m);
     define_fit_bindings(m);
     define_interpolation_bindings(m);
+    define_jungfrau_data_file_io_bindings(m);
 
-    define_cluster_file_io_bindings<Cluster<int, 3, 3>>(m, "Cluster3x3i");
-    define_cluster_file_io_bindings<Cluster<double, 3, 3>>(m, "Cluster3x3d");
-    define_cluster_file_io_bindings<Cluster<float, 3, 3>>(m, "Cluster3x3f");
-    define_cluster_file_io_bindings<Cluster<int, 2, 2>>(m, "Cluster2x2i");
-    define_cluster_file_io_bindings<Cluster<float, 2, 2>>(m, "Cluster2x2f");
-    define_cluster_file_io_bindings<Cluster<double, 2, 2>>(m, "Cluster2x2d");
-
-    define_cluster_vector<int, 3, 3, uint16_t>(m, "Cluster3x3i");
-    define_cluster_vector<double, 3, 3, uint16_t>(m, "Cluster3x3d");
-    define_cluster_vector<float, 3, 3, uint16_t>(m, "Cluster3x3f");
-    define_cluster_vector<int, 2, 2, uint16_t>(m, "Cluster2x2i");
-    define_cluster_vector<double, 2, 2, uint16_t>(m, "Cluster2x2d");
-    define_cluster_vector<float, 2, 2, uint16_t>(m, "Cluster2x2f");
-
-    define_cluster_finder_bindings<Cluster<int, 3, 3>>(m, "Cluster3x3i");
-    define_cluster_finder_bindings<Cluster<double, 3, 3>>(m, "Cluster3x3d");
-    define_cluster_finder_bindings<Cluster<float, 3, 3>>(m, "Cluster3x3f");
-    define_cluster_finder_bindings<Cluster<int, 2, 2>>(m, "Cluster2x2i");
-    define_cluster_finder_bindings<Cluster<double, 2, 2>>(m, "Cluster2x2d");
-    define_cluster_finder_bindings<Cluster<float, 2, 2>>(m, "Cluster2x2f");
-
-    define_cluster_finder_mt_bindings<Cluster<int, 3, 3>>(m, "Cluster3x3i");
-    define_cluster_finder_mt_bindings<Cluster<double, 3, 3>>(m, "Cluster3x3d");
-    define_cluster_finder_mt_bindings<Cluster<float, 3, 3>>(m, "Cluster3x3f");
-    define_cluster_finder_mt_bindings<Cluster<int, 2, 2>>(m, "Cluster2x2i");
-    define_cluster_finder_mt_bindings<Cluster<double, 2, 2>>(m, "Cluster2x2d");
-    define_cluster_finder_mt_bindings<Cluster<float, 2, 2>>(m, "Cluster2x2f");
-
-    define_cluster_file_sink_bindings<Cluster<int, 3, 3>>(m, "Cluster3x3i");
-    define_cluster_file_sink_bindings<Cluster<double, 3, 3>>(m, "Cluster3x3d");
-    define_cluster_file_sink_bindings<Cluster<float, 3, 3>>(m, "Cluster3x3f");
-    define_cluster_file_sink_bindings<Cluster<int, 2, 2>>(m, "Cluster2x2i");
-    define_cluster_file_sink_bindings<Cluster<double, 2, 2>>(m, "Cluster2x2d");
-    define_cluster_file_sink_bindings<Cluster<float, 2, 2>>(m, "Cluster2x2f");
-
-    define_cluster_collector_bindings<Cluster<int, 3, 3>>(m, "Cluster3x3i");
-    define_cluster_collector_bindings<Cluster<double, 3, 3>>(m, "Cluster3x3f");
-    define_cluster_collector_bindings<Cluster<float, 3, 3>>(m, "Cluster3x3d");
-    define_cluster_collector_bindings<Cluster<int, 2, 2>>(m, "Cluster2x2i");
-    define_cluster_collector_bindings<Cluster<double, 2, 2>>(m, "Cluster2x2f");
-    define_cluster_collector_bindings<Cluster<float, 2, 2>>(m, "Cluster2x2d");
-
-    define_cluster<int, 3, 3, uint16_t>(m, "3x3i");
-    define_cluster<float, 3, 3, uint16_t>(m, "3x3f");
-    define_cluster<double, 3, 3, uint16_t>(m, "3x3d");
-    define_cluster<int, 2, 2, uint16_t>(m, "2x2i");
-    define_cluster<float, 2, 2, uint16_t>(m, "2x2f");
-    define_cluster<double, 2, 2, uint16_t>(m, "2x2d");
-}
+}

python/src/np_helper.hpp

@@ -10,7 +10,6 @@
 #include "aare/NDView.hpp"
 
 namespace py = pybind11;
-using namespace aare;
 
 // Pass image data back to python as a numpy array
 template <typename T, int64_t Ndim>
@@ -41,46 +40,25 @@ template <typename T> py::array return_vector(std::vector<T> *vec) {
 }
 
 // todo rewrite generic
-template <class T, int Flags>
-auto get_shape_3d(const py::array_t<T, Flags> &arr) {
+template <class T, int Flags> auto get_shape_3d(const py::array_t<T, Flags>& arr) {
     return aare::Shape<3>{arr.shape(0), arr.shape(1), arr.shape(2)};
 }
 
-template <class T, int Flags> auto make_view_3d(py::array_t<T, Flags> &arr) {
+template <class T, int Flags> auto make_view_3d(py::array_t<T, Flags>& arr) {
     return aare::NDView<T, 3>(arr.mutable_data(), get_shape_3d<T, Flags>(arr));
 }
 
-template <class T, int Flags>
-auto get_shape_2d(const py::array_t<T, Flags> &arr) {
+template <class T, int Flags> auto get_shape_2d(const py::array_t<T, Flags>& arr) {
     return aare::Shape<2>{arr.shape(0), arr.shape(1)};
 }
 
-template <class T, int Flags>
-auto get_shape_1d(const py::array_t<T, Flags> &arr) {
+template <class T, int Flags> auto get_shape_1d(const py::array_t<T, Flags>& arr) {
     return aare::Shape<1>{arr.shape(0)};
 }
 
-template <class T, int Flags> auto make_view_2d(py::array_t<T, Flags> &arr) {
+template <class T, int Flags> auto make_view_2d(py::array_t<T, Flags>& arr) {
     return aare::NDView<T, 2>(arr.mutable_data(), get_shape_2d<T, Flags>(arr));
 }
-template <class T, int Flags> auto make_view_1d(py::array_t<T, Flags> &arr) {
+template <class T, int Flags> auto make_view_1d(py::array_t<T, Flags>& arr) {
     return aare::NDView<T, 1>(arr.mutable_data(), get_shape_1d<T, Flags>(arr));
-}
-
-template <typename ClusterType> struct fmt_format_trait; // forward declaration
-
-template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType>
-struct fmt_format_trait<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
-
-    static std::string value() {
-        return fmt::format("T{{{}:x;{}:y;{}:data;}}",
-                           py::format_descriptor<CoordType>::format(),
-                           py::format_descriptor<CoordType>::format(),
-                           fmt::format("{}{}", ClusterSizeX * ClusterSizeY,
-                                       py::format_descriptor<T>::format()));
-    }
-};
-
-template <typename ClusterType>
-auto fmt_format = fmt_format_trait<ClusterType>::value();
+}

python/src/raw_sub_file.hpp

@@ -0,0 +1,110 @@
+#include "aare/CtbRawFile.hpp"
+#include "aare/File.hpp"
+#include "aare/Frame.hpp"
+#include "aare/RawFile.hpp"
+#include "aare/RawMasterFile.hpp"
+#include "aare/RawSubFile.hpp"
+
+#include "aare/defs.hpp"
+// #include "aare/fClusterFileV2.hpp"
+
+#include <cstdint>
+#include <filesystem>
+#include <pybind11/iostream.h>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/stl/filesystem.h>
+#include <string>
+
+namespace py = pybind11;
+using namespace ::aare;
+
+auto read_frame_from_RawSubFile(RawSubFile &self) {
+    py::array_t<DetectorHeader> header(1);
+    const uint8_t item_size = self.bytes_per_pixel();
+    std::vector<ssize_t> shape{static_cast<ssize_t>(self.rows()),
+        static_cast<ssize_t>(self.cols())};
+    
+    py::array image;
+    if (item_size == 1) {
+        image = py::array_t<uint8_t>(shape);
+    } else if (item_size == 2) {
+        image = py::array_t<uint16_t>(shape);
+    } else if (item_size == 4) {
+        image = py::array_t<uint32_t>(shape);
+    }
+    self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()),
+                   header.mutable_data());
+
+    return py::make_tuple(header, image);
+}
+
+auto read_n_frames_from_RawSubFile(RawSubFile &self, size_t n_frames) {
+    py::array_t<DetectorHeader> header(n_frames);
+    const uint8_t item_size = self.bytes_per_pixel();
+    std::vector<ssize_t> shape{
+        static_cast<ssize_t>(n_frames), 
+        static_cast<ssize_t>(self.rows()),
+        static_cast<ssize_t>(self.cols())
+    };
+    
+    py::array image;
+    if (item_size == 1) {
+        image = py::array_t<uint8_t>(shape);
+    } else if (item_size == 2) {
+        image = py::array_t<uint16_t>(shape);
+    } else if (item_size == 4) {
+        image = py::array_t<uint32_t>(shape);
+    }
+    self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()), n_frames,
+                   header.mutable_data());
+
+    return py::make_tuple(header, image);
+}
+
+
+//Disable warnings for unused parameters, as we ignore some
+//in the __exit__ method
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wunused-parameter"
+
+void define_raw_sub_file_io_bindings(py::module &m) {
+    py::class_<RawSubFile>(m, "RawSubFile")
+        .def(py::init<const std::filesystem::path &, DetectorType, size_t,
+                      size_t, size_t>())
+        .def_property_readonly("bytes_per_frame", &RawSubFile::bytes_per_frame)
+        .def_property_readonly("pixels_per_frame",
+                               &RawSubFile::pixels_per_frame)
+        .def_property_readonly("bytes_per_pixel", &RawSubFile::bytes_per_pixel)                    
+        .def("seek", &RawSubFile::seek)
+        .def("tell", &RawSubFile::tell)
+        .def_property_readonly("rows", &RawSubFile::rows)
+        .def_property_readonly("cols", &RawSubFile::cols)
+        .def_property_readonly("frames_in_file", &RawSubFile::frames_in_file)
+        .def("read_frame", &read_frame_from_RawSubFile)
+        .def("read_n", &read_n_frames_from_RawSubFile)
+        .def("read", [](RawSubFile &self){
+            self.seek(0);
+            auto n_frames = self.frames_in_file();
+            return read_n_frames_from_RawSubFile(self, n_frames);
+        })
+        .def("__enter__", [](RawSubFile &self) { return &self; })
+        .def("__exit__",
+             [](RawSubFile &self,
+                const std::optional<pybind11::type> &exc_type,
+                const std::optional<pybind11::object> &exc_value,
+                const std::optional<pybind11::object> &traceback) {
+             })
+        .def("__iter__", [](RawSubFile &self) { return &self; })
+        .def("__next__", [](RawSubFile &self) {
+            try {
+                return read_frame_from_RawSubFile(self);
+            } catch (std::runtime_error &e) {
+                throw py::stop_iteration();
+            }
+        });
+
+}
+
+#pragma GCC diagnostic pop

python/tests/conftest.py

@@ -0,0 +1,29 @@
+import os
+from pathlib import Path
+import pytest
+
+
+
+def pytest_addoption(parser):
+    parser.addoption(
+        "--files", action="store_true", default=False, help="run slow tests"
+    )
+
+
+def pytest_configure(config):
+    config.addinivalue_line("markers", "files: mark test as needing image files to run")
+
+
+def pytest_collection_modifyitems(config, items):
+    if config.getoption("--files"):
+        return
+    skip = pytest.mark.skip(reason="need --files option to run")
+    for item in items:
+        if "files" in item.keywords:
+            item.add_marker(skip)
+
+
+@pytest.fixture
+def test_data_path():
+    return Path(os.environ["AARE_TEST_DATA"])
+

python/tests/test_Cluster.py

@@ -1,64 +0,0 @@
-import pytest 
-import numpy as np
-
-from _aare import ClusterVector_Cluster3x3i, Interpolator, Cluster3x3i, ClusterFinder_Cluster3x3i
-
-def test_ClusterVector(): 
-    """Test ClusterVector""" 
-
-    clustervector = ClusterVector_Cluster3x3i()
-    assert clustervector.cluster_size_x == 3
-    assert clustervector.cluster_size_y == 3
-    assert clustervector.item_size() == 4+9*4
-    assert clustervector.frame_number == 0
-    assert clustervector.capacity == 1024
-    assert clustervector.size == 0
-
-    cluster = Cluster3x3i(0,0,np.ones(9, dtype=np.int32))
-    
-    clustervector.push_back(cluster) 
-    assert clustervector.size == 1
-
-    #push_back - check size
-
-
-
-def test_Interpolator(): 
-    """Test Interpolator"""
-
-    ebins = np.linspace(0,10, 20, dtype=np.float64)
-    xbins = np.linspace(0, 5, 30, dtype=np.float64)
-    ybins = np.linspace(0, 5, 30, dtype=np.float64)
-
-    etacube = np.zeros(shape=[30, 30, 20], dtype=np.float64)
-    interpolator = Interpolator(etacube, xbins, ybins, ebins)
-
-    assert interpolator.get_ietax().shape == (30,30,20)
-    assert interpolator.get_ietay().shape == (30,30,20)
-    clustervector = ClusterVector_Cluster3x3i()
-    
-    cluster = Cluster3x3i(0,0, np.ones(9, dtype=np.int32))
-    #clustervector.push_back(cluster) 
-    #num_clusters = 1; 
-
-    #assert interpolator.interpolate_Cluster3x3i(clustervector).shape == (num_clusters, 3)
-
-
-#def test_cluster_file(): 
-
-#def test_cluster_finder(): 
-    #"""Test ClusterFinder""" 
-
-    #clusterfinder = ClusterFinder_Cluster3x3i([100,100])
-
-    #clusterfinder.find_clusters()
-
-    #clusters = clusterfinder.steal_clusters()
-
-    #print("cluster size: ", clusters.size()) 
-
-
-
-
-
-

python/tests/test_RawSubFile.py

@@ -0,0 +1,36 @@
+import pytest
+import numpy as np
+from aare import RawSubFile, DetectorType
+
+
+@pytest.mark.files
+def test_read_a_jungfrau_RawSubFile(test_data_path):
+    with RawSubFile(test_data_path / "raw/jungfrau/jungfrau_single_d0_f1_0.raw", DetectorType.Jungfrau, 512, 1024, 16) as f:
+        assert f.frames_in_file == 3
+
+        headers, frames = f.read()
+
+    assert headers.size == 3
+    assert frames.shape == (3, 512, 1024)
+    
+    # Frame numbers in this file should be 4, 5, 6
+    for i,h in zip(range(4,7,1), headers):
+        assert h["frameNumber"] == i
+
+    # Compare to canned data using numpy
+    data = np.load(test_data_path / "raw/jungfrau/jungfrau_single_0.npy")
+    assert np.all(data[3:6] == frames)
+
+@pytest.mark.files
+def test_iterate_over_a_jungfrau_RawSubFile(test_data_path):
+
+    data = np.load(test_data_path / "raw/jungfrau/jungfrau_single_0.npy")
+
+    with RawSubFile(test_data_path / "raw/jungfrau/jungfrau_single_d0_f0_0.raw", DetectorType.Jungfrau, 512, 1024, 16) as f:
+        i = 0
+        for header, frame in f:
+            assert header["frameNumber"] == i+1
+            assert np.all(frame == data[i])
+            i += 1
+        assert i == 3
+        assert header["frameNumber"] == 3

python/tests/test_jungfrau_dat_files.py

@@ -0,0 +1,92 @@
+import pytest
+import numpy as np
+from aare import JungfrauDataFile
+
+@pytest.mark.files
+def test_jfungfrau_dat_read_number_of_frames(test_data_path):
+    with JungfrauDataFile(test_data_path / "dat/AldoJF500k_000000.dat") as dat_file:
+        assert dat_file.total_frames == 24
+
+    with JungfrauDataFile(test_data_path / "dat/AldoJF250k_000000.dat") as dat_file:
+        assert dat_file.total_frames == 53
+    
+    with JungfrauDataFile(test_data_path / "dat/AldoJF65k_000000.dat") as dat_file:
+        assert dat_file.total_frames == 113
+
+
+@pytest.mark.files
+def test_jfungfrau_dat_read_number_of_file(test_data_path):
+    with JungfrauDataFile(test_data_path / "dat/AldoJF500k_000000.dat") as dat_file:
+        assert dat_file.n_files == 4
+
+    with JungfrauDataFile(test_data_path / "dat/AldoJF250k_000000.dat") as dat_file:
+        assert dat_file.n_files == 7
+    
+    with JungfrauDataFile(test_data_path / "dat/AldoJF65k_000000.dat") as dat_file:
+        assert dat_file.n_files == 7
+
+
+@pytest.mark.files
+def test_read_module(test_data_path):
+    """
+    Read all frames from the series of .dat files. Compare to canned data in npz format. 
+    """
+
+    # Read all frames from the .dat file
+    with JungfrauDataFile(test_data_path / "dat/AldoJF500k_000000.dat") as f:
+        header, data = f.read()
+
+    #Sanity check
+    n_frames = 24
+    assert header.size == n_frames
+    assert data.shape == (n_frames, 512, 1024)
+
+    # Read reference data using numpy
+    with np.load(test_data_path / "dat/AldoJF500k.npz") as f:
+        ref_header = f["headers"]
+        ref_data = f["frames"]
+
+    # Check that the data is the same
+    assert np.all(ref_header == header)
+    assert np.all(ref_data == data)
+
+@pytest.mark.files
+def test_read_half_module(test_data_path):
+
+    # Read all frames from the .dat file
+    with JungfrauDataFile(test_data_path / "dat/AldoJF250k_000000.dat") as f:
+        header, data = f.read()
+
+    n_frames = 53
+    assert header.size == n_frames
+    assert data.shape == (n_frames, 256, 1024)
+
+    # Read reference data using numpy
+    with np.load(test_data_path / "dat/AldoJF250k.npz") as f:
+        ref_header = f["headers"]
+        ref_data = f["frames"]
+
+    # Check that the data is the same
+    assert np.all(ref_header == header)
+    assert np.all(ref_data == data)
+
+
+@pytest.mark.files
+def test_read_single_chip(test_data_path):
+
+    # Read all frames from the .dat file
+    with JungfrauDataFile(test_data_path / "dat/AldoJF65k_000000.dat") as f:
+        header, data = f.read()
+
+    n_frames = 113
+    assert header.size == n_frames
+    assert data.shape == (n_frames, 256, 256)
+
+    # Read reference data using numpy
+    with np.load(test_data_path / "dat/AldoJF65k.npz") as f:
+        ref_header = f["headers"]
+        ref_data = f["frames"]
+
+    # Check that the data is the same
+    assert np.all(ref_header == header)
+    assert np.all(ref_data == data)

src/CalculateEta.test.cpp

@@ -1,62 +0,0 @@
-/************************************************
- * @file CalculateEta.test.cpp
- * @short test case to calculate_eta2
- ***********************************************/
-
-#include "aare/CalculateEta.hpp"
-#include "aare/Cluster.hpp"
-#include "aare/ClusterFile.hpp"
-
-// #include "catch.hpp"
-#include <array>
-#include <catch2/catch_all.hpp>
-#include <catch2/catch_test_macros.hpp>
-
-using namespace aare;
-
-using ClusterTypes =
-    std::variant<Cluster<int, 2, 2>, Cluster<int, 3, 3>, Cluster<int, 5, 5>,
-                 Cluster<int, 4, 2>, Cluster<int, 2, 3>>;
-
-auto get_test_parameters() {
-    return GENERATE(
-        std::make_tuple(ClusterTypes{Cluster<int, 2, 2>{0, 0, {1, 2, 3, 1}}},
-                        Eta2<int>{2. / 3, 3. / 4, corner::cBottomLeft, 7}),
-        std::make_tuple(
-            ClusterTypes{Cluster<int, 3, 3>{0, 0, {1, 2, 3, 4, 5, 6, 1, 2, 7}}},
-            Eta2<int>{6. / 11, 2. / 7, corner::cTopRight, 20}),
-        std::make_tuple(ClusterTypes{Cluster<int, 5, 5>{
-                            0, 0, {1, 6, 7, 6, 5, 4, 3, 2, 1, 8, 8, 9, 2,
-                                   1, 4, 5, 6, 7, 8, 4, 1, 1, 1, 1, 1}}},
-                        Eta2<int>{9. / 17, 5. / 13, 8, 28}),
-        std::make_tuple(
-            ClusterTypes{Cluster<int, 4, 2>{0, 0, {1, 4, 7, 2, 5, 6, 4, 3}}},
-            Eta2<int>{7. / 11, 6. / 10, 1, 21}),
-        std::make_tuple(
-            ClusterTypes{Cluster<int, 2, 3>{0, 0, {1, 3, 2, 3, 4, 2}}},
-            Eta2<int>{3. / 5, 4. / 6, 1, 11}));
-}
-
-TEST_CASE("compute_largest_2x2_subcluster", "[.eta_calculation]") {
-    auto [cluster, expected_eta] = get_test_parameters();
-
-    auto [sum, index] = std::visit(
-        [](const auto &clustertype) { return clustertype.max_sum_2x2(); },
-        cluster);
-    CHECK(expected_eta.c == index);
-    CHECK(expected_eta.sum == sum);
-}
-
-TEST_CASE("calculate_eta2", "[.eta_calculation]") {
-
-    auto [cluster, expected_eta] = get_test_parameters();
-
-    auto eta = std::visit(
-        [](const auto &clustertype) { return calculate_eta2(clustertype); },
-        cluster);
-
-    CHECK(eta.x == expected_eta.x);
-    CHECK(eta.y == expected_eta.y);
-    CHECK(eta.c == expected_eta.c);
-    CHECK(eta.sum == expected_eta.sum);
-}

src/Cluster.test.cpp

@@ -1,28 +0,0 @@
-/************************************************
- * @file test-Cluster.cpp
- * @short test case for generic Cluster, ClusterVector, and calculate_eta2
- ***********************************************/
-
-#include "aare/Cluster.hpp"
-#include "aare/CalculateEta.hpp"
-#include "aare/ClusterFile.hpp"
-
-// #include "catch.hpp"
-#include <array>
-#include <catch2/catch_all.hpp>
-#include <catch2/catch_test_macros.hpp>
-
-using namespace aare;
-
-TEST_CASE("Correct Instantiation of Cluster and ClusterVector",
-          "[.cluster][.instantiation]") {
-
-    CHECK(is_valid_cluster<double, 3, 3>);
-    CHECK(is_valid_cluster<double, 3, 2>);
-    CHECK(not is_valid_cluster<int, 0, 0>);
-    CHECK(not is_valid_cluster<std::string, 2, 2>);
-    CHECK(not is_valid_cluster<int, 2, 2, double>);
-
-    CHECK(not is_cluster_v<int>);
-    CHECK(is_cluster_v<Cluster<int, 3, 3>>);
-}

src/ClusterFile.cpp

@@ -0,0 +1,402 @@
+#include "aare/ClusterFile.hpp"
+
+#include <algorithm>
+
+namespace aare {
+
+ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size,
+                         const std::string &mode)
+    : m_chunk_size(chunk_size), m_mode(mode) {
+
+    if (mode == "r") {
+        fp = fopen(fname.c_str(), "rb");
+        if (!fp) {
+            throw std::runtime_error("Could not open file for reading: " +
+                                     fname.string());
+        }
+    } else if (mode == "w") {
+        fp = fopen(fname.c_str(), "wb");
+        if (!fp) {
+            throw std::runtime_error("Could not open file for writing: " +
+                                     fname.string());
+        }
+    } else if (mode == "a") {
+        fp = fopen(fname.c_str(), "ab");
+        if (!fp) {
+            throw std::runtime_error("Could not open file for appending: " +
+                                     fname.string());
+        }
+    } else {
+        throw std::runtime_error("Unsupported mode: " + mode);
+    }
+}
+
+void ClusterFile::set_roi(ROI roi){
+    m_roi = roi;
+}
+
+void ClusterFile::set_noise_map(const NDView<int32_t, 2> noise_map){
+    m_noise_map = NDArray<int32_t, 2>(noise_map);
+}
+
+void ClusterFile::set_gain_map(const NDView<double, 2> gain_map){
+    m_gain_map = NDArray<double, 2>(gain_map);
+    
+    // Gain map is passed as ADU/keV to avoid dividing in when applying the gain
+    // map we invert it here
+    for (auto &item : m_gain_map->view()) {
+        item = 1.0 / item;
+    }
+}
+
+ClusterFile::~ClusterFile() { close(); }
+
+void ClusterFile::close() {
+    if (fp) {
+        fclose(fp);
+        fp = nullptr;
+    }
+}
+
+void ClusterFile::write_frame(const ClusterVector<int32_t> &clusters) {
+    if (m_mode != "w" && m_mode != "a") {
+        throw std::runtime_error("File not opened for writing");
+    }
+    if (!(clusters.cluster_size_x() == 3) &&
+        !(clusters.cluster_size_y() == 3)) {
+        throw std::runtime_error("Only 3x3 clusters are supported");
+    }
+    //First write the frame number - 4 bytes
+    int32_t frame_number = clusters.frame_number();
+    if(fwrite(&frame_number, sizeof(frame_number), 1, fp)!=1){
+        throw std::runtime_error(LOCATION + "Could not write frame number");
+    }
+
+    //Then write the number of clusters - 4 bytes
+    uint32_t n_clusters = clusters.size();
+    if(fwrite(&n_clusters, sizeof(n_clusters), 1, fp)!=1){
+        throw std::runtime_error(LOCATION + "Could not write number of clusters");
+    }
+
+    //Now write the clusters in the frame
+    if(fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp)!=clusters.size()){
+        throw std::runtime_error(LOCATION + "Could not write clusters");
+    }
+}
+
+
+ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters){
+    if (m_mode != "r") {
+        throw std::runtime_error("File not opened for reading");
+    }
+    if (m_noise_map || m_roi){
+        return read_clusters_with_cut(n_clusters);
+    }else{
+        return read_clusters_without_cut(n_clusters);
+    }
+}
+
+ClusterVector<int32_t> ClusterFile::read_clusters_without_cut(size_t n_clusters) {
+    if (m_mode != "r") {
+        throw std::runtime_error("File not opened for reading");
+    }
+    
+    ClusterVector<int32_t> clusters(3,3, n_clusters);
+
+    int32_t iframe = 0; // frame number needs to be 4 bytes!
+    size_t nph_read = 0;
+    uint32_t nn = m_num_left;
+    uint32_t nph = m_num_left; // number of clusters in frame needs to be 4
+
+    // auto buf = reinterpret_cast<Cluster3x3 *>(clusters.data());
+    auto buf = clusters.data();
+    // if there are photons left from previous frame read them first
+    if (nph) {
+        if (nph > n_clusters) {
+            // if we have more photons left in the frame then photons to read we
+            // read directly the requested number
+            nn = n_clusters;
+        } else {
+            nn = nph;
+        }
+        nph_read += fread((buf + nph_read*clusters.item_size()),
+                          clusters.item_size(), nn, fp);
+        m_num_left = nph - nn; // write back the number of photons left
+    }
+
+    if (nph_read < n_clusters) {
+        // keep on reading frames and photons until reaching n_clusters
+        while (fread(&iframe, sizeof(iframe), 1, fp)) {
+            clusters.set_frame_number(iframe);
+            // read number of clusters in frame
+            if (fread(&nph, sizeof(nph), 1, fp)) {
+                if (nph > (n_clusters - nph_read))
+                    nn = n_clusters - nph_read;
+                else
+                    nn = nph;
+
+                nph_read += fread((buf + nph_read*clusters.item_size()),
+                                  clusters.item_size(), nn, fp);
+                m_num_left = nph - nn;
+            }
+            if (nph_read >= n_clusters)
+                break;
+        }
+    }
+
+    // Resize the vector to the number of clusters.
+    // No new allocation, only change bounds.
+    clusters.resize(nph_read);
+    if(m_gain_map)
+        clusters.apply_gain_map(m_gain_map->view());
+    return clusters;
+}
+
+
+ClusterVector<int32_t> ClusterFile::read_clusters_with_cut(size_t n_clusters) {
+    ClusterVector<int32_t> clusters(3,3);
+    clusters.reserve(n_clusters);
+
+    // if there are photons left from previous frame read them first
+    if (m_num_left) {
+        while(m_num_left && clusters.size() < n_clusters){
+            Cluster3x3 c = read_one_cluster();
+            if(is_selected(c)){
+                clusters.push_back(c.x, c.y, reinterpret_cast<std::byte*>(c.data));
+            }
+        }
+    }
+
+    // we did not have enough clusters left in the previous frame
+    // keep on reading frames until reaching n_clusters
+    if (clusters.size() < n_clusters) {
+        // sanity check
+        if (m_num_left) {
+            throw std::runtime_error(LOCATION + "Entered second loop with clusters left\n");
+        }
+        
+        int32_t frame_number = 0; // frame number needs to be 4 bytes!
+        while (fread(&frame_number, sizeof(frame_number), 1, fp)) {
+            if (fread(&m_num_left, sizeof(m_num_left), 1, fp)) {
+                clusters.set_frame_number(frame_number); //cluster vector will hold the last frame number
+                while(m_num_left && clusters.size() < n_clusters){
+                    Cluster3x3 c = read_one_cluster();
+                    if(is_selected(c)){
+                        clusters.push_back(c.x, c.y, reinterpret_cast<std::byte*>(c.data));
+                    }
+                }
+            }
+
+            // we have enough clusters, break out of the outer while loop
+            if (clusters.size() >= n_clusters)
+                break;
+        }
+
+    }
+    if(m_gain_map)
+        clusters.apply_gain_map(m_gain_map->view());
+
+    return clusters;
+}
+
+Cluster3x3 ClusterFile::read_one_cluster(){
+    Cluster3x3 c;
+    auto rc = fread(&c, sizeof(c), 1, fp);
+    if (rc != 1) {
+        throw std::runtime_error(LOCATION + "Could not read cluster");
+    }
+    --m_num_left;
+    return c;
+}
+
+ClusterVector<int32_t> ClusterFile::read_frame(){
+    if (m_mode != "r") {
+        throw std::runtime_error(LOCATION + "File not opened for reading");
+    }
+    if (m_noise_map || m_roi){
+        return read_frame_with_cut();
+    }else{
+        return read_frame_without_cut();
+    }
+}
+
+ClusterVector<int32_t> ClusterFile::read_frame_without_cut() {
+    if (m_mode != "r") {
+        throw std::runtime_error("File not opened for reading");
+    }
+    if (m_num_left) {
+        throw std::runtime_error(
+            "There are still photons left in the last frame");
+    }
+    int32_t frame_number;
+    if (fread(&frame_number, sizeof(frame_number), 1, fp) != 1) {
+        throw std::runtime_error(LOCATION + "Could not read frame number");
+    }
+
+    int32_t n_clusters; // Saved as 32bit integer in the cluster file
+    if (fread(&n_clusters, sizeof(n_clusters), 1, fp) != 1) {
+        throw std::runtime_error(LOCATION + "Could not read number of clusters");
+    }
+
+    ClusterVector<int32_t> clusters(3, 3, n_clusters);
+    clusters.set_frame_number(frame_number);
+
+    if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) !=
+        static_cast<size_t>(n_clusters)) {
+        throw std::runtime_error(LOCATION + "Could not read clusters");
+    }
+    clusters.resize(n_clusters);
+    if (m_gain_map)
+        clusters.apply_gain_map(m_gain_map->view());
+    return clusters;
+}
+
+ClusterVector<int32_t> ClusterFile::read_frame_with_cut() {
+    if (m_mode != "r") {
+        throw std::runtime_error("File not opened for reading");
+    }
+    if (m_num_left) {
+        throw std::runtime_error(
+            "There are still photons left in the last frame");
+    }
+    int32_t frame_number;
+    if (fread(&frame_number, sizeof(frame_number), 1, fp) != 1) {
+        throw std::runtime_error("Could not read frame number");
+    }
+
+
+    if (fread(&m_num_left, sizeof(m_num_left), 1, fp) != 1) {
+        throw std::runtime_error("Could not read number of clusters");
+    }
+    
+    ClusterVector<int32_t> clusters(3, 3);
+    clusters.reserve(m_num_left);
+    clusters.set_frame_number(frame_number);
+    while(m_num_left){
+        Cluster3x3 c = read_one_cluster();
+        if(is_selected(c)){
+            clusters.push_back(c.x, c.y, reinterpret_cast<std::byte*>(c.data));
+        }
+    }
+    if (m_gain_map)
+        clusters.apply_gain_map(m_gain_map->view());
+    return clusters;
+}
+
+
+
+bool ClusterFile::is_selected(Cluster3x3 &cl) {
+    //Should fail fast
+    if (m_roi) {
+        if (!(m_roi->contains(cl.x, cl.y))) {
+            return false;
+        }
+    }
+    if (m_noise_map){
+        int32_t sum_1x1 = cl.data[4]; // central pixel
+        int32_t sum_2x2 = cl.sum_2x2(); // highest sum of 2x2 subclusters
+        int32_t sum_3x3 = cl.sum(); // sum of all pixels
+
+        auto noise = (*m_noise_map)(cl.y, cl.x); //TODO! check if this is correct
+        if (sum_1x1 <= noise || sum_2x2 <= 2 * noise || sum_3x3 <= 3 * noise) {
+            return false;
+        }
+    }
+    //we passed all checks
+    return true;
+}
+
+NDArray<double, 2> calculate_eta2(ClusterVector<int> &clusters) {
+    //TOTO! make work with 2x2 clusters
+    NDArray<double, 2> eta2({static_cast<int64_t>(clusters.size()), 2});
+    
+    if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) {
+        for (size_t i = 0; i < clusters.size(); i++) {
+            auto e = calculate_eta2(clusters.at<Cluster3x3>(i));
+            eta2(i, 0) = e.x;
+            eta2(i, 1) = e.y;
+        }
+    }else if(clusters.cluster_size_x() == 2 || clusters.cluster_size_y() == 2){
+        for (size_t i = 0; i < clusters.size(); i++) {
+            auto e = calculate_eta2(clusters.at<Cluster2x2>(i));
+            eta2(i, 0) = e.x;
+            eta2(i, 1) = e.y;
+        }
+    }else{
+        throw std::runtime_error("Only 3x3 and 2x2 clusters are supported");
+    }
+    
+    return eta2;
+}
+
+/** 
+ * @brief Calculate the eta2 values for a 3x3 cluster and return them in a Eta2 struct
+ * containing etay, etax and the corner of the cluster. 
+*/
+Eta2 calculate_eta2(Cluster3x3 &cl) {
+    Eta2 eta{};
+
+    std::array<int32_t, 4> tot2;
+    tot2[0] = cl.data[0] + cl.data[1] + cl.data[3] + cl.data[4];
+    tot2[1] = cl.data[1] + cl.data[2] + cl.data[4] + cl.data[5];
+    tot2[2] = cl.data[3] + cl.data[4] + cl.data[6] + cl.data[7];
+    tot2[3] = cl.data[4] + cl.data[5] + cl.data[7] + cl.data[8];
+
+    auto c = std::max_element(tot2.begin(), tot2.end()) - tot2.begin();
+    eta.sum = tot2[c];
+    switch (c) {
+    case cBottomLeft:
+        if ((cl.data[3] + cl.data[4]) != 0)
+            eta.x =
+                static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
+        if ((cl.data[1] + cl.data[4]) != 0)
+            eta.y =
+                static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
+        eta.c = cBottomLeft;
+        break;
+    case cBottomRight:
+        if ((cl.data[2] + cl.data[5]) != 0)
+            eta.x =
+                static_cast<double>(cl.data[5]) / (cl.data[4] + cl.data[5]);
+        if ((cl.data[1] + cl.data[4]) != 0)
+            eta.y =
+                static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
+        eta.c = cBottomRight;
+        break;
+    case cTopLeft:
+        if ((cl.data[7] + cl.data[4]) != 0)
+            eta.x =
+                static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
+        if ((cl.data[7] + cl.data[4]) != 0)
+            eta.y =
+                static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
+        eta.c = cTopLeft;
+        break;
+    case cTopRight:
+        if ((cl.data[5] + cl.data[4]) != 0)
+            eta.x =
+                static_cast<double>(cl.data[5]) / (cl.data[5] + cl.data[4]);
+        if ((cl.data[7] + cl.data[4]) != 0)
+            eta.y =
+                static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
+        eta.c = cTopRight;
+        break;
+    // no default to allow compiler to warn about missing cases
+    }
+    return eta;
+}
+
+
+Eta2 calculate_eta2(Cluster2x2 &cl) {
+    Eta2 eta{};
+    if ((cl.data[0] + cl.data[1]) != 0)
+        eta.x = static_cast<double>(cl.data[1]) / (cl.data[0] + cl.data[1]);
+    if ((cl.data[0] + cl.data[2]) != 0)
+        eta.y = static_cast<double>(cl.data[2]) / (cl.data[0] + cl.data[2]);
+    eta.sum = cl.data[0] + cl.data[1] + cl.data[2]+ cl.data[3];
+    eta.c = cBottomLeft; //TODO! This is not correct, but need to put something
+    return eta;
+}
+
+
+} // namespace aare

src/ClusterFile.test.cpp

@@ -1,32 +1,35 @@
 #include "aare/ClusterFile.hpp"
 #include "test_config.hpp"
 
+
 #include "aare/defs.hpp"
 #include <catch2/catch_test_macros.hpp>
 #include <filesystem>
 
-using aare::Cluster;
+
+
+
 using aare::ClusterFile;
 
-TEST_CASE("Read one frame from a a cluster file", "[.integration]") {
-    // We know that the frame has 97 clusters
-    auto fpath =
-        test_data_path() / "clusters" / "single_frame_97_clustrers.clust";
+
+TEST_CASE("Read one frame from a a cluster file", "[.files]") {
+    //We know that the frame has 97 clusters
+    auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust";
     REQUIRE(std::filesystem::exists(fpath));
 
-    ClusterFile<Cluster<int32_t, 3, 3>> f(fpath);
+    ClusterFile f(fpath);
     auto clusters = f.read_frame();
     REQUIRE(clusters.size() == 97);
     REQUIRE(clusters.frame_number() == 135);
 }
 
-TEST_CASE("Read one frame using ROI", "[.integration]") {
-    // We know that the frame has 97 clusters
-    auto fpath =
-        test_data_path() / "clusters" / "single_frame_97_clustrers.clust";
+
+TEST_CASE("Read one frame using ROI", "[.files]") {
+    //We know that the frame has 97 clusters
+    auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust";
     REQUIRE(std::filesystem::exists(fpath));
 
-    ClusterFile<Cluster<int32_t, 3, 3>> f(fpath);
+    ClusterFile f(fpath);
     aare::ROI roi;
     roi.xmin = 0;
     roi.xmax = 50;
@@ -37,39 +40,45 @@ TEST_CASE("Read one frame using ROI", "[.integration]") {
     REQUIRE(clusters.size() == 49);
     REQUIRE(clusters.frame_number() == 135);
 
-    // Check that all clusters are within the ROI
+    //Check that all clusters are within the ROI
     for (size_t i = 0; i < clusters.size(); i++) {
-        auto c = clusters.at(i);
+        auto c = clusters.at<aare::Cluster3x3>(i);
         REQUIRE(c.x >= roi.xmin);
         REQUIRE(c.x <= roi.xmax);
         REQUIRE(c.y >= roi.ymin);
         REQUIRE(c.y <= roi.ymax);
     }
+
 }
 
-TEST_CASE("Read clusters from single frame file", "[.integration]") {
 
-    auto fpath =
-        test_data_path() / "clusters" / "single_frame_97_clustrers.clust";
+
+TEST_CASE("Read clusters from single frame file", "[.files]") {
+
+    auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust";
+
     REQUIRE(std::filesystem::exists(fpath));
 
     SECTION("Read fewer clusters than available") {
-        ClusterFile<Cluster<int32_t, 3, 3>> f(fpath);
+        ClusterFile f(fpath);
         auto clusters = f.read_clusters(50);
         REQUIRE(clusters.size() == 50);
-        REQUIRE(clusters.frame_number() == 135);
+         REQUIRE(clusters.frame_number() == 135);
     }
     SECTION("Read more clusters than available") {
-        ClusterFile<Cluster<int32_t, 3, 3>> f(fpath);
+        ClusterFile f(fpath);
         // 100 is the maximum number of clusters read
         auto clusters = f.read_clusters(100);
         REQUIRE(clusters.size() == 97);
         REQUIRE(clusters.frame_number() == 135);
     }
     SECTION("Read all clusters") {
-        ClusterFile<Cluster<int32_t, 3, 3>> f(fpath);
+        ClusterFile f(fpath);
         auto clusters = f.read_clusters(97);
         REQUIRE(clusters.size() == 97);
         REQUIRE(clusters.frame_number() == 135);
     }
+
+
+    
 }

src/ClusterFinder.test.cpp

@@ -1,18 +1,19 @@
 #include "aare/ClusterFinder.hpp"
 #include "aare/Pedestal.hpp"
-#include <catch2/catch_test_macros.hpp>
 #include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <catch2/catch_test_macros.hpp>
 #include <chrono>
 #include <random>
 
 using namespace aare;
 
-// TODO! Find a way to test the cluster finder
+//TODO! Find a way to test the cluster finder
+
+
 
 // class ClusterFinderUnitTest : public ClusterFinder {
 //   public:
-//     ClusterFinderUnitTest(int cluster_sizeX, int cluster_sizeY, double nSigma
-//     = 5.0, double threshold = 0.0)
+//     ClusterFinderUnitTest(int cluster_sizeX, int cluster_sizeY, double nSigma = 5.0, double threshold = 0.0)
 //         : ClusterFinder(cluster_sizeX, cluster_sizeY, nSigma, threshold) {}
 //     double get_c2() { return c2; }
 //     double get_c3() { return c3; }
@@ -36,8 +37,8 @@ using namespace aare;
 //     REQUIRE_THAT(cf.get_c3(), Catch::Matchers::WithinRel(c3, 1e-9));
 // }
 
-TEST_CASE("Construct a cluster finder") {
-    ClusterFinder clusterFinder({400, 400});
+TEST_CASE("Construct a cluster finder"){
+    ClusterFinder clusterFinder({400,400}, {3,3});
     // REQUIRE(clusterFinder.get_cluster_sizeX() == 3);
     // REQUIRE(clusterFinder.get_cluster_sizeY() == 3);
     // REQUIRE(clusterFinder.get_threshold() == 1);
@@ -48,17 +49,16 @@ TEST_CASE("Construct a cluster finder") {
 //     aare::Pedestal pedestal(10, 10, 5);
 //     NDArray<double, 2> frame({10, 10});
 //     frame = 0;
-//     ClusterFinder clusterFinder(3, 3, 1, 1); // 3x3 cluster, 1 nSigma, 1
-//     threshold
+//     ClusterFinder clusterFinder(3, 3, 1, 1); // 3x3 cluster, 1 nSigma, 1 threshold
 
-//     auto clusters =
-//     clusterFinder.find_clusters_without_threshold(frame.span(), pedestal);
+//     auto clusters = clusterFinder.find_clusters_without_threshold(frame.span(), pedestal);
 
 //     REQUIRE(clusters.size() == 0);
 
 //     frame(5, 5) = 10;
-//     clusters = clusterFinder.find_clusters_without_threshold(frame.span(),
-//     pedestal); REQUIRE(clusters.size() == 1); REQUIRE(clusters[0].x == 5);
+//     clusters = clusterFinder.find_clusters_without_threshold(frame.span(), pedestal);
+//     REQUIRE(clusters.size() == 1);
+//     REQUIRE(clusters[0].x == 5);
 //     REQUIRE(clusters[0].y == 5);
 //     for (int i = 0; i < 3; i++) {
 //         for (int j = 0; j < 3; j++) {

src/ClusterVector.test.cpp

@@ -1,17 +1,21 @@
-#include "aare/ClusterVector.hpp"
 #include <cstdint>
+#include "aare/ClusterVector.hpp"
 
-#include <catch2/catch_all.hpp>
-#include <catch2/catch_test_macros.hpp>
 #include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <catch2/catch_test_macros.hpp>
 
-using aare::Cluster;
 using aare::ClusterVector;
 
-TEST_CASE("ClusterVector 2x2 int32_t capacity 4, push back then read",
-          "[.ClusterVector]") {
+struct Cluster_i2x2 {
+    int16_t x;
+    int16_t y;
+    int32_t data[4];
+};
 
-    ClusterVector<Cluster<int32_t, 2, 2>> cv(4);
+TEST_CASE("ClusterVector 2x2 int32_t capacity 4, push back then read") {
+    
+
+    ClusterVector<int32_t> cv(2, 2, 4);
     REQUIRE(cv.capacity() == 4);
     REQUIRE(cv.size() == 0);
     REQUIRE(cv.cluster_size_x() == 2);
@@ -19,102 +23,112 @@ TEST_CASE("ClusterVector 2x2 int32_t capacity 4, push back then read",
     // int16_t, int16_t, 2x2 int32_t = 20 bytes
     REQUIRE(cv.item_size() == 20);
 
-    // Create a cluster and push back into the vector
-    Cluster<int32_t, 2, 2> c1 = {1, 2, {3, 4, 5, 6}};
-    cv.push_back(c1);
+    //Create a cluster and push back into the vector
+    Cluster_i2x2 c1 = {1, 2, {3, 4, 5, 6}};
+    cv.push_back(c1.x, c1.y, reinterpret_cast<std::byte*>(&c1.data[0]));
     REQUIRE(cv.size() == 1);
     REQUIRE(cv.capacity() == 4);
 
-    auto c2 = cv.at(0);
+    //Read the cluster back out using copy. TODO! Can we improve the API?
+    Cluster_i2x2 c2;
+    std::byte *ptr = cv.element_ptr(0);
+    std::copy(ptr, ptr + cv.item_size(), reinterpret_cast<std::byte*>(&c2));
 
-    // Check that the data is the same
+    //Check that the data is the same
     REQUIRE(c1.x == c2.x);
     REQUIRE(c1.y == c2.y);
-    for (size_t i = 0; i < 4; i++) {
+    for(size_t i = 0; i < 4; i++) {
         REQUIRE(c1.data[i] == c2.data[i]);
     }
 }
 
-TEST_CASE("Summing 3x1 clusters of int64", "[.ClusterVector]") {
-    ClusterVector<Cluster<int32_t, 3, 1>> cv(2);
+TEST_CASE("Summing 3x1 clusters of int64"){
+    struct Cluster_l3x1{
+        int16_t x;
+        int16_t y;
+        int32_t data[3];
+    };
+
+    ClusterVector<int32_t> cv(3, 1, 2);
     REQUIRE(cv.capacity() == 2);
     REQUIRE(cv.size() == 0);
     REQUIRE(cv.cluster_size_x() == 3);
     REQUIRE(cv.cluster_size_y() == 1);
 
-    // Create a cluster and push back into the vector
-    Cluster<int32_t, 3, 1> c1 = {1, 2, {3, 4, 5}};
-    cv.push_back(c1);
+    //Create a cluster and push back into the vector
+    Cluster_l3x1 c1 = {1, 2, {3, 4, 5}};
+    cv.push_back(c1.x, c1.y, reinterpret_cast<std::byte*>(&c1.data[0]));
     REQUIRE(cv.capacity() == 2);
     REQUIRE(cv.size() == 1);
 
-    Cluster<int32_t, 3, 1> c2 = {6, 7, {8, 9, 10}};
-    cv.push_back(c2);
+    Cluster_l3x1 c2 = {6, 7, {8, 9, 10}};
+    cv.push_back(c2.x, c2.y, reinterpret_cast<std::byte*>(&c2.data[0]));
     REQUIRE(cv.capacity() == 2);
     REQUIRE(cv.size() == 2);
 
-    Cluster<int32_t, 3, 1> c3 = {11, 12, {13, 14, 15}};
-    cv.push_back(c3);
+    Cluster_l3x1 c3 = {11, 12, {13, 14, 15}};
+    cv.push_back(c3.x, c3.y, reinterpret_cast<std::byte*>(&c3.data[0]));
     REQUIRE(cv.capacity() == 4);
     REQUIRE(cv.size() == 3);
 
-    /*
     auto sums = cv.sum();
     REQUIRE(sums.size() == 3);
     REQUIRE(sums[0] == 12);
     REQUIRE(sums[1] == 27);
     REQUIRE(sums[2] == 42);
-    */
 }
 
-TEST_CASE("Storing floats", "[.ClusterVector]") {
-    ClusterVector<Cluster<float, 2, 4>> cv(10);
+TEST_CASE("Storing floats"){
+    struct Cluster_f4x2{
+        int16_t x;
+        int16_t y;
+        float data[8];
+    };
+
+    ClusterVector<float> cv(2, 4, 10);
     REQUIRE(cv.capacity() == 10);
     REQUIRE(cv.size() == 0);
     REQUIRE(cv.cluster_size_x() == 2);
     REQUIRE(cv.cluster_size_y() == 4);
 
-    // Create a cluster and push back into the vector
-    Cluster<float, 2, 4> c1 = {1, 2, {3.0, 4.0, 5.0, 6.0, 3.0, 4.0, 5.0, 6.0}};
-    cv.push_back(c1);
+    //Create a cluster and push back into the vector
+    Cluster_f4x2 c1 = {1, 2, {3.0, 4.0, 5.0, 6.0,3.0, 4.0, 5.0, 6.0}};
+    cv.push_back(c1.x, c1.y, reinterpret_cast<std::byte*>(&c1.data[0]));
     REQUIRE(cv.capacity() == 10);
     REQUIRE(cv.size() == 1);
+    
 
-    Cluster<float, 2, 4> c2 = {
-        6, 7, {8.0, 9.0, 10.0, 11.0, 8.0, 9.0, 10.0, 11.0}};
-    cv.push_back(c2);
+    Cluster_f4x2 c2 = {6, 7, {8.0, 9.0, 10.0, 11.0,8.0, 9.0, 10.0, 11.0}};
+    cv.push_back(c2.x, c2.y, reinterpret_cast<std::byte*>(&c2.data[0]));
     REQUIRE(cv.capacity() == 10);
     REQUIRE(cv.size() == 2);
 
-    /*
     auto sums = cv.sum();
     REQUIRE(sums.size() == 2);
     REQUIRE_THAT(sums[0], Catch::Matchers::WithinAbs(36.0, 1e-6));
     REQUIRE_THAT(sums[1], Catch::Matchers::WithinAbs(76.0, 1e-6));
-    */
 }
 
-TEST_CASE("Push back more than initial capacity", "[.ClusterVector]") {
-
-    ClusterVector<Cluster<int32_t, 2, 2>> cv(2);
+TEST_CASE("Push back more than initial capacity"){
+    
+    ClusterVector<int32_t> cv(2, 2, 2);
     auto initial_data = cv.data();
-    Cluster<int32_t, 2, 2> c1 = {1, 2, {3, 4, 5, 6}};
-    cv.push_back(c1);
+    Cluster_i2x2 c1 = {1, 2, {3, 4, 5, 6}};
+    cv.push_back(c1.x, c1.y, reinterpret_cast<std::byte*>(&c1.data[0]));
     REQUIRE(cv.size() == 1);
     REQUIRE(cv.capacity() == 2);
 
-    Cluster<int32_t, 2, 2> c2 = {6, 7, {8, 9, 10, 11}};
-    cv.push_back(c2);
+    Cluster_i2x2 c2 = {6, 7, {8, 9, 10, 11}};
+    cv.push_back(c2.x, c2.y, reinterpret_cast<std::byte*>(&c2.data[0]));
     REQUIRE(cv.size() == 2);
     REQUIRE(cv.capacity() == 2);
 
-    Cluster<int32_t, 2, 2> c3 = {11, 12, {13, 14, 15, 16}};
-    cv.push_back(c3);
-    REQUIRE(cv.size() == 3);
+    Cluster_i2x2 c3 = {11, 12, {13, 14, 15, 16}};
+    cv.push_back(c3.x, c3.y, reinterpret_cast<std::byte*>(&c3.data[0]));
+    REQUIRE(cv.size() == 3); 
     REQUIRE(cv.capacity() == 4);
 
-    Cluster<int32_t, 2, 2> *ptr =
-        reinterpret_cast<Cluster<int32_t, 2, 2> *>(cv.data());
+    Cluster_i2x2* ptr = reinterpret_cast<Cluster_i2x2*>(cv.data());
     REQUIRE(ptr[0].x == 1);
     REQUIRE(ptr[0].y == 2);
     REQUIRE(ptr[1].x == 6);
@@ -122,31 +136,29 @@ TEST_CASE("Push back more than initial capacity", "[.ClusterVector]") {
     REQUIRE(ptr[2].x == 11);
     REQUIRE(ptr[2].y == 12);
 
-    // We should have allocated a new buffer, since we outgrew the initial
-    // capacity
+    //We should have allocated a new buffer, since we outgrew the initial capacity
     REQUIRE(initial_data != cv.data());
+
 }
 
-TEST_CASE("Concatenate two cluster vectors where the first has enough capacity",
-          "[.ClusterVector]") {
-    ClusterVector<Cluster<int32_t, 2, 2>> cv1(12);
-    Cluster<int32_t, 2, 2> c1 = {1, 2, {3, 4, 5, 6}};
-    cv1.push_back(c1);
-    Cluster<int32_t, 2, 2> c2 = {6, 7, {8, 9, 10, 11}};
-    cv1.push_back(c2);
+TEST_CASE("Concatenate two cluster vectors where the first has enough capacity"){
+    ClusterVector<int32_t> cv1(2, 2, 12);
+    Cluster_i2x2 c1 = {1, 2, {3, 4, 5, 6}};
+    cv1.push_back(c1.x, c1.y, reinterpret_cast<std::byte*>(&c1.data[0]));
+    Cluster_i2x2 c2 = {6, 7, {8, 9, 10, 11}};
+    cv1.push_back(c2.x, c2.y, reinterpret_cast<std::byte*>(&c2.data[0]));
 
-    ClusterVector<Cluster<int32_t, 2, 2>> cv2(2);
-    Cluster<int32_t, 2, 2> c3 = {11, 12, {13, 14, 15, 16}};
-    cv2.push_back(c3);
-    Cluster<int32_t, 2, 2> c4 = {16, 17, {18, 19, 20, 21}};
-    cv2.push_back(c4);
+    ClusterVector<int32_t> cv2(2, 2, 2);
+    Cluster_i2x2 c3 = {11, 12, {13, 14, 15, 16}};
+    cv2.push_back(c3.x, c3.y, reinterpret_cast<std::byte*>(&c3.data[0]));
+    Cluster_i2x2 c4 = {16, 17, {18, 19, 20, 21}};
+    cv2.push_back(c4.x, c4.y, reinterpret_cast<std::byte*>(&c4.data[0]));
 
     cv1 += cv2;
     REQUIRE(cv1.size() == 4);
     REQUIRE(cv1.capacity() == 12);
 
-    Cluster<int32_t, 2, 2> *ptr =
-        reinterpret_cast<Cluster<int32_t, 2, 2> *>(cv1.data());
+    Cluster_i2x2* ptr = reinterpret_cast<Cluster_i2x2*>(cv1.data());
     REQUIRE(ptr[0].x == 1);
     REQUIRE(ptr[0].y == 2);
     REQUIRE(ptr[1].x == 6);
@@ -157,26 +169,24 @@ TEST_CASE("Concatenate two cluster vectors where the first has enough capacity",
     REQUIRE(ptr[3].y == 17);
 }
 
-TEST_CASE("Concatenate two cluster vectors where we need to allocate",
-          "[.ClusterVector]") {
-    ClusterVector<Cluster<int32_t, 2, 2>> cv1(2);
-    Cluster<int32_t, 2, 2> c1 = {1, 2, {3, 4, 5, 6}};
-    cv1.push_back(c1);
-    Cluster<int32_t, 2, 2> c2 = {6, 7, {8, 9, 10, 11}};
-    cv1.push_back(c2);
+TEST_CASE("Concatenate two cluster vectors where we need to allocate"){
+    ClusterVector<int32_t> cv1(2, 2, 2);
+    Cluster_i2x2 c1 = {1, 2, {3, 4, 5, 6}};
+    cv1.push_back(c1.x, c1.y, reinterpret_cast<std::byte*>(&c1.data[0]));
+    Cluster_i2x2 c2 = {6, 7, {8, 9, 10, 11}};
+    cv1.push_back(c2.x, c2.y, reinterpret_cast<std::byte*>(&c2.data[0]));
 
-    ClusterVector<Cluster<int32_t, 2, 2>> cv2(2);
-    Cluster<int32_t, 2, 2> c3 = {11, 12, {13, 14, 15, 16}};
-    cv2.push_back(c3);
-    Cluster<int32_t, 2, 2> c4 = {16, 17, {18, 19, 20, 21}};
-    cv2.push_back(c4);
+    ClusterVector<int32_t> cv2(2, 2, 2);
+    Cluster_i2x2 c3 = {11, 12, {13, 14, 15, 16}};
+    cv2.push_back(c3.x, c3.y, reinterpret_cast<std::byte*>(&c3.data[0]));
+    Cluster_i2x2 c4 = {16, 17, {18, 19, 20, 21}};
+    cv2.push_back(c4.x, c4.y, reinterpret_cast<std::byte*>(&c4.data[0]));
 
     cv1 += cv2;
     REQUIRE(cv1.size() == 4);
     REQUIRE(cv1.capacity() == 4);
 
-    Cluster<int32_t, 2, 2> *ptr =
-        reinterpret_cast<Cluster<int32_t, 2, 2> *>(cv1.data());
+    Cluster_i2x2* ptr = reinterpret_cast<Cluster_i2x2*>(cv1.data());
     REQUIRE(ptr[0].x == 1);
     REQUIRE(ptr[0].y == 2);
     REQUIRE(ptr[1].x == 6);
@@ -185,49 +195,4 @@ TEST_CASE("Concatenate two cluster vectors where we need to allocate",
     REQUIRE(ptr[2].y == 12);
     REQUIRE(ptr[3].x == 16);
     REQUIRE(ptr[3].y == 17);
-}
-
-struct ClusterTestData {
-    uint8_t ClusterSizeX;
-    uint8_t ClusterSizeY;
-    std::vector<int64_t> index_map_x;
-    std::vector<int64_t> index_map_y;
-};
-
-TEST_CASE("Gain Map Calculation Index Map", "[.ClusterVector][.gain_map]") {
-
-    auto clustertestdata = GENERATE(
-        ClusterTestData{3,
-                        3,
-                        {-1, 0, 1, -1, 0, 1, -1, 0, 1},
-                        {-1, -1, -1, 0, 0, 0, 1, 1, 1}},
-        ClusterTestData{
-            4,
-            4,
-            {-2, -1, 0, 1, -2, -1, 0, 1, -2, -1, 0, 1, -2, -1, 0, 1},
-            {-2, -2, -2, -2, -1, -1, -1, -1, 0, 0, 0, 0, 1, 1, 1, 1}},
-        ClusterTestData{2, 2, {-1, 0, -1, 0}, {-1, -1, 0, 0}},
-        ClusterTestData{5,
-                        5,
-                        {-2, -1, 0,  1,  2, -2, -1, 0,  1,  2, -2, -1, 0,
-                         1,  2,  -2, -1, 0, 1,  2,  -2, -1, 0, 1,  2},
-                        {-2, -2, -2, -2, -2, -1, -1, -1, -1, -1, 0, 0, 0,
-                         0,  0,  1,  1,  1,  1,  1,  2,  2,  2,  2, 2}});
-
-    uint8_t ClusterSizeX = clustertestdata.ClusterSizeX;
-    uint8_t ClusterSizeY = clustertestdata.ClusterSizeY;
-
-    std::vector<int64_t> index_map_x(ClusterSizeX * ClusterSizeY);
-    std::vector<int64_t> index_map_y(ClusterSizeX * ClusterSizeY);
-
-    int64_t index_cluster_center_x = ClusterSizeX / 2;
-    int64_t index_cluster_center_y = ClusterSizeY / 2;
-
-    for (size_t j = 0; j < ClusterSizeX * ClusterSizeY; j++) {
-        index_map_x[j] = j % ClusterSizeX - index_cluster_center_x;
-        index_map_y[j] = j / ClusterSizeX - index_cluster_center_y;
-    }
-
-    CHECK(index_map_x == clustertestdata.index_map_x);
-    CHECK(index_map_y == clustertestdata.index_map_y);
 }

src/File.cpp

@@ -1,4 +1,5 @@
 #include "aare/File.hpp"
+#include "aare/JungfrauDataFile.hpp"
 #include "aare/NumpyFile.hpp"
 #include "aare/RawFile.hpp"
 
@@ -27,6 +28,8 @@ File::File(const std::filesystem::path &fname, const std::string &mode,
     else if (fname.extension() == ".npy") {
         // file_impl = new NumpyFile(fname, mode, cfg);
         file_impl = std::make_unique<NumpyFile>(fname, mode, cfg);
+    }else if(fname.extension() == ".dat"){
+        file_impl = std::make_unique<JungfrauDataFile>(fname);
     } else {
         throw std::runtime_error("Unsupported file type");
     }

src/FilePtr.cpp

@@ -0,0 +1,44 @@
+
+#include "aare/FilePtr.hpp"
+#include <fmt/format.h>
+#include <stdexcept>
+#include <utility>
+
+namespace aare {
+
+FilePtr::FilePtr(const std::filesystem::path& fname, const std::string& mode = "rb") {
+    fp_ = fopen(fname.c_str(), mode.c_str());
+    if (!fp_)
+        throw std::runtime_error(fmt::format("Could not open: {}", fname.c_str()));
+}
+
+FilePtr::FilePtr(FilePtr &&other) { std::swap(fp_, other.fp_); }
+
+FilePtr &FilePtr::operator=(FilePtr &&other) {
+    std::swap(fp_, other.fp_);
+    return *this;
+}
+
+FILE *FilePtr::get() { return fp_; }
+
+int64_t FilePtr::tell() {
+    auto pos = ftell(fp_);
+    if (pos == -1)
+        throw std::runtime_error(fmt::format("Error getting file position: {}", error_msg()));
+    return pos;
+} 
+FilePtr::~FilePtr() {
+    if (fp_)
+        fclose(fp_); // check?
+}
+
+std::string FilePtr::error_msg(){
+    if (feof(fp_)) {
+        return "End of file reached";
+    }
+    if (ferror(fp_)) {
+        return fmt::format("Error reading file: {}", std::strerror(errno));
+    }
+    return "";
+}
+} // namespace aare

src/Fit.cpp

@@ -18,7 +18,7 @@ double gaus(const double x, const double *par) {
 
 NDArray<double, 1> gaus(NDView<double, 1> x, NDView<double, 1> par) {
     NDArray<double, 1> y({x.shape(0)}, 0);
-    for (size_t i = 0; i < x.size(); i++) {
+    for (ssize_t i = 0; i < x.size(); i++) {
         y(i) = gaus(x(i), par.data());
     }
     return y;
@@ -28,7 +28,7 @@ double pol1(const double x, const double *par) { return par[0] * x + par[1]; }
 
 NDArray<double, 1> pol1(NDView<double, 1> x, NDView<double, 1> par) {
     NDArray<double, 1> y({x.shape()}, 0);
-    for (size_t i = 0; i < x.size(); i++) {
+    for (ssize_t i = 0; i < x.size(); i++) {
         y(i) = pol1(x(i), par.data());
     }
     return y;
@@ -153,7 +153,7 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
 
     // Calculate chi2
     chi2 = 0;
-    for (size_t i = 0; i < y.size(); i++) {
+    for (ssize_t i = 0; i < y.size(); i++) {
         chi2 += std::pow((y(i) - func::gaus(x(i), par_out.data())) / y_err(i), 2);
     }
 }
@@ -205,7 +205,7 @@ void fit_pol1(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
 
     // Calculate chi2
     chi2 = 0;
-    for (size_t i = 0; i < y.size(); i++) {
+    for (ssize_t i = 0; i < y.size(); i++) {
         chi2 += std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
     }
 }

src/Interpolator.cpp

@@ -1,11 +1,11 @@
 #include "aare/Interpolator.hpp"
+#include "aare/algorithm.hpp"
 
 namespace aare {
 
 Interpolator::Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
                            NDView<double, 1> ybins, NDView<double, 1> ebins)
-    : m_ietax(etacube), m_ietay(etacube), m_etabinsx(xbins), m_etabinsy(ybins),
-      m_energy_bins(ebins) {
+    : m_ietax(etacube), m_ietay(etacube), m_etabinsx(xbins), m_etabinsy(ybins), m_energy_bins(ebins) {
     if (etacube.shape(0) != xbins.size() || etacube.shape(1) != ybins.size() ||
         etacube.shape(2) != ebins.size()) {
         throw std::invalid_argument(
@@ -53,4 +53,87 @@ Interpolator::Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
     }
 }
 
+std::vector<Photon> Interpolator::interpolate(const ClusterVector<int32_t>& clusters) {
+    std::vector<Photon> photons;
+    photons.reserve(clusters.size());
+
+    if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) {
+        for (size_t i = 0; i<clusters.size(); i++){
+
+            auto cluster = clusters.at<Cluster3x3>(i);
+            Eta2 eta= calculate_eta2(cluster);
+    
+            Photon photon;
+            photon.x = cluster.x;
+            photon.y = cluster.y;
+            photon.energy = eta.sum;
+    
+
+            //Finding the index of the last element that is smaller
+            //should work fine as long as we have many bins
+            auto ie = last_smaller(m_energy_bins, photon.energy);
+            auto ix = last_smaller(m_etabinsx, eta.x);
+            auto iy = last_smaller(m_etabinsy, eta.y); 
+            
+            double dX{}, dY{};
+            // cBottomLeft = 0,
+            // cBottomRight = 1,
+            // cTopLeft = 2,
+            // cTopRight = 3
+            switch (eta.c) {
+            case cTopLeft:
+                dX = -1.;
+                dY = 0.;
+                break;
+            case cTopRight:;
+                dX = 0.;
+                dY = 0.;
+                break;
+            case cBottomLeft:
+                dX = -1.;
+                dY = -1.;
+                break;
+            case cBottomRight:
+                dX = 0.; 
+                dY = -1.;
+                break;
+            }
+            photon.x += m_ietax(ix, iy, ie)*2 + dX;
+            photon.y += m_ietay(ix, iy, ie)*2 + dY;
+            photons.push_back(photon);
+        }
+    }else if(clusters.cluster_size_x() == 2 || clusters.cluster_size_y() == 2){
+        for (size_t i = 0; i<clusters.size(); i++){
+            auto cluster = clusters.at<Cluster2x2>(i);
+            Eta2 eta= calculate_eta2(cluster);
+    
+            Photon photon;
+            photon.x = cluster.x;
+            photon.y = cluster.y;
+            photon.energy = eta.sum;
+    
+            //Now do some actual interpolation. 
+            //Find which energy bin the cluster is in
+            // auto ie = nearest_index(m_energy_bins, photon.energy)-1;
+            // auto ix = nearest_index(m_etabinsx, eta.x)-1;
+            // auto iy = nearest_index(m_etabinsy, eta.y)-1;  
+            //Finding the index of the last element that is smaller
+            //should work fine as long as we have many bins
+            auto ie = last_smaller(m_energy_bins, photon.energy);
+            auto ix = last_smaller(m_etabinsx, eta.x);
+            auto iy = last_smaller(m_etabinsy, eta.y); 
+
+            photon.x += m_ietax(ix, iy, ie)*2; //eta goes between 0 and 1 but we could move the hit anywhere in the 2x2
+            photon.y += m_ietay(ix, iy, ie)*2;
+            photons.push_back(photon);
+        }
+        
+    }else{
+        throw std::runtime_error("Only 3x3 and 2x2 clusters are supported for interpolation");
+    }
+  
+
+    return photons;
+}
+
 } // namespace aare

src/JungfrauDataFile.cpp

@@ -0,0 +1,238 @@
+#include "aare/JungfrauDataFile.hpp"
+#include "aare/algorithm.hpp"
+#include "aare/defs.hpp"
+
+#include <cerrno>
+#include <fmt/format.h>
+
+namespace aare {
+
+JungfrauDataFile::JungfrauDataFile(const std::filesystem::path &fname) {
+
+    if (!std::filesystem::exists(fname)) {
+        throw std::runtime_error(LOCATION +
+                                 "File does not exist: " + fname.string());
+    }
+    find_frame_size(fname);
+    parse_fname(fname);
+    scan_files();
+    open_file(m_current_file_index);
+}
+
+
+// FileInterface
+
+Frame JungfrauDataFile::read_frame(){
+    Frame f(rows(), cols(), Dtype::UINT16);
+    read_into(reinterpret_cast<std::byte *>(f.data()), nullptr);
+    return f;
+}
+
+Frame JungfrauDataFile::read_frame(size_t frame_number){
+    seek(frame_number);
+    Frame f(rows(), cols(), Dtype::UINT16);
+    read_into(reinterpret_cast<std::byte *>(f.data()), nullptr);
+    return f;
+}
+
+std::vector<Frame> JungfrauDataFile::read_n(size_t n_frames) {
+    std::vector<Frame> frames;
+    for(size_t i = 0; i < n_frames; ++i){
+        frames.push_back(read_frame());
+    }
+    return frames;
+}
+
+void JungfrauDataFile::read_into(std::byte *image_buf) {
+    read_into(image_buf, nullptr);
+}
+void JungfrauDataFile::read_into(std::byte *image_buf, size_t n_frames) {
+    read_into(image_buf, n_frames, nullptr);
+}   
+
+size_t JungfrauDataFile::frame_number(size_t frame_index) {
+    seek(frame_index);
+    return read_header().framenum;
+}
+
+std::array<ssize_t, 2> JungfrauDataFile::shape() const {
+    return {static_cast<ssize_t>(rows()), static_cast<ssize_t>(cols())};
+}
+
+DetectorType JungfrauDataFile::detector_type() const { return DetectorType::Jungfrau; }
+
+std::string JungfrauDataFile::base_name() const { return m_base_name; }
+
+size_t JungfrauDataFile::bytes_per_frame() { return m_bytes_per_frame; }
+
+size_t JungfrauDataFile::pixels_per_frame() { return m_rows * m_cols; }
+
+size_t JungfrauDataFile::bytes_per_pixel() const { return sizeof(pixel_type); }
+
+size_t JungfrauDataFile::bitdepth() const {
+    return bytes_per_pixel() * bits_per_byte;
+}
+
+void JungfrauDataFile::seek(size_t frame_index) {
+    if (frame_index >= m_total_frames) {
+        throw std::runtime_error(LOCATION + "Frame index out of range: " +
+                                 std::to_string(frame_index));
+    }
+    m_current_frame_index = frame_index;
+    auto file_index = first_larger(m_last_frame_in_file, frame_index);
+
+    if (file_index != m_current_file_index)
+        open_file(file_index);
+
+    auto frame_offset = (file_index)
+                            ? frame_index - m_last_frame_in_file[file_index - 1]
+                            : frame_index;
+    auto byte_offset = frame_offset * (m_bytes_per_frame + header_size);
+    m_fp.seek(byte_offset);
+};
+
+size_t JungfrauDataFile::tell() { return m_current_frame_index; }
+size_t JungfrauDataFile::total_frames() const { return m_total_frames; }
+size_t JungfrauDataFile::rows() const { return m_rows; }
+size_t JungfrauDataFile::cols() const { return m_cols; }
+
+size_t JungfrauDataFile::n_files() const { return m_last_frame_in_file.size(); }
+
+void JungfrauDataFile::find_frame_size(const std::filesystem::path &fname) {
+
+    static constexpr size_t module_data_size =
+        header_size + sizeof(pixel_type) * 512 * 1024;
+    static constexpr size_t half_data_size =
+        header_size + sizeof(pixel_type) * 256 * 1024;
+    static constexpr size_t chip_data_size =
+        header_size + sizeof(pixel_type) * 256 * 256;
+
+    auto file_size = std::filesystem::file_size(fname);
+    if (file_size == 0) {
+        throw std::runtime_error(LOCATION +
+                                 "Cannot guess frame size: file is empty");
+    }
+
+    if (file_size % module_data_size == 0) {
+        m_rows = 512;
+        m_cols = 1024;
+        m_bytes_per_frame = module_data_size - header_size;
+    } else if (file_size % half_data_size == 0) {
+        m_rows = 256;
+        m_cols = 1024;
+        m_bytes_per_frame = half_data_size - header_size;
+    } else if (file_size % chip_data_size == 0) {
+        m_rows = 256;
+        m_cols = 256;
+        m_bytes_per_frame = chip_data_size - header_size;
+    } else {
+        throw std::runtime_error(LOCATION +
+                                 "Cannot find frame size: file size is not a "
+                                 "multiple of any known frame size");
+    }
+}
+
+void JungfrauDataFile::parse_fname(const std::filesystem::path &fname) {
+    m_path = fname.parent_path();
+    m_base_name = fname.stem();
+
+    // find file index, then remove if from the base name
+    if (auto pos = m_base_name.find_last_of('_'); pos != std::string::npos) {
+        m_offset = std::stoul(m_base_name.substr(pos + 1));
+        m_base_name.erase(pos);
+    }
+}
+
+void JungfrauDataFile::scan_files() {
+    // find how many files we have and the number of frames in each file
+    m_last_frame_in_file.clear();
+    size_t file_index = m_offset;
+    while (std::filesystem::exists(fpath(file_index))) {
+        auto n_frames = std::filesystem::file_size(fpath(file_index)) /
+                        (m_bytes_per_frame + header_size);
+        m_last_frame_in_file.push_back(n_frames);
+        ++file_index;
+    }
+
+    // find where we need to open the next file and total number of frames
+    m_last_frame_in_file = cumsum(m_last_frame_in_file);
+    m_total_frames = m_last_frame_in_file.back();
+}
+
+void JungfrauDataFile::read_into(std::byte *image_buf,
+                                 JungfrauDataHeader *header) {
+
+    // read header if not passed nullptr
+    if (header) {
+        if (auto rc = fread(header, sizeof(JungfrauDataHeader), 1, m_fp.get());
+            rc != 1) {
+            throw std::runtime_error(
+                LOCATION +
+                "Could not read header from file:" + m_fp.error_msg());
+        }
+    } else {
+        m_fp.seek(header_size, SEEK_CUR);
+    }
+
+    // read data
+    if (auto rc = fread(image_buf, 1, m_bytes_per_frame, m_fp.get());
+        rc != m_bytes_per_frame) {
+        throw std::runtime_error(LOCATION + "Could not read image from file" +
+                                 m_fp.error_msg());
+    }
+
+    // prepare for next read
+    // if we are at the end of the file, open the next file
+    ++m_current_frame_index;
+    if (m_current_frame_index >= m_last_frame_in_file[m_current_file_index] &&
+        (m_current_frame_index < m_total_frames)) {
+        ++m_current_file_index;
+        open_file(m_current_file_index);
+    }
+}
+
+void JungfrauDataFile::read_into(std::byte *image_buf, size_t n_frames,
+                                 JungfrauDataHeader *header) {
+    if (header) {
+        for (size_t i = 0; i < n_frames; ++i)
+            read_into(image_buf + i * m_bytes_per_frame, header + i); 
+    }else{
+        for (size_t i = 0; i < n_frames; ++i)
+            read_into(image_buf + i * m_bytes_per_frame, nullptr);
+    }
+}
+
+void JungfrauDataFile::read_into(NDArray<uint16_t>* image, JungfrauDataHeader* header) {
+    if(image->shape()!=shape()){
+        throw std::runtime_error(LOCATION +
+            "Image shape does not match file size: " + std::to_string(rows()) + "x" + std::to_string(cols()));
+    }
+    read_into(reinterpret_cast<std::byte *>(image->data()), header);
+}
+
+
+JungfrauDataHeader JungfrauDataFile::read_header() {
+    JungfrauDataHeader header;
+    if (auto rc = fread(&header, 1, sizeof(header), m_fp.get());
+        rc != sizeof(header)) {
+        throw std::runtime_error(LOCATION + "Could not read header from file" +
+                                 m_fp.error_msg());
+    }
+    m_fp.seek(-header_size, SEEK_CUR);
+    return header;
+}
+
+void JungfrauDataFile::open_file(size_t file_index) {
+    // fmt::print(stderr, "Opening file: {}\n",
+    // fpath(file_index+m_offset).string());
+    m_fp = FilePtr(fpath(file_index + m_offset), "rb");
+    m_current_file_index = file_index;
+}
+
+std::filesystem::path JungfrauDataFile::fpath(size_t file_index) const {
+    auto fname = fmt::format("{}_{:0{}}.dat", m_base_name, file_index,
+                             n_digits_in_file_index);
+    return m_path / fname;
+}
+
+} // namespace aare

src/JungfrauDataFile.test.cpp

@@ -0,0 +1,114 @@
+#include "aare/JungfrauDataFile.hpp"
+
+#include <catch2/catch_test_macros.hpp>
+#include "test_config.hpp"
+
+using aare::JungfrauDataFile;
+using aare::JungfrauDataHeader;
+TEST_CASE("Open a Jungfrau data file", "[.files]") {
+    //we know we have 4 files with 7, 7, 7, and 3 frames
+    //firs frame number if 1 and the bunch id is frame_number**2
+    //so we can check the header
+    auto fpath = test_data_path() / "dat" / "AldoJF500k_000000.dat";
+    REQUIRE(std::filesystem::exists(fpath));
+
+    JungfrauDataFile f(fpath);
+    REQUIRE(f.rows() == 512);
+    REQUIRE(f.cols() == 1024);
+    REQUIRE(f.bytes_per_frame() == 1048576);    
+    REQUIRE(f.pixels_per_frame() == 524288);
+    REQUIRE(f.bytes_per_pixel() == 2);
+    REQUIRE(f.bitdepth() == 16);
+    REQUIRE(f.base_name() == "AldoJF500k");
+    REQUIRE(f.n_files() == 4);
+    REQUIRE(f.tell() == 0);
+    REQUIRE(f.total_frames() == 24);
+    REQUIRE(f.current_file() == fpath);
+
+    //Check that the frame number and buch id is read correctly
+    for (size_t i = 0; i < 24; ++i) {
+        JungfrauDataHeader header;
+        aare::NDArray<uint16_t> image(f.shape());
+        f.read_into(&image, &header);
+        REQUIRE(header.framenum == i + 1);
+        REQUIRE(header.bunchid == (i + 1) * (i + 1));
+        REQUIRE(image.shape(0) == 512);
+        REQUIRE(image.shape(1) == 1024);
+    }
+}
+
+TEST_CASE("Seek in a JungfrauDataFile", "[.files]"){
+    auto fpath = test_data_path() / "dat" / "AldoJF65k_000000.dat";
+    REQUIRE(std::filesystem::exists(fpath));
+
+    JungfrauDataFile f(fpath);
+
+    //The file should have 113 frames
+    f.seek(19);
+    REQUIRE(f.tell() == 19);
+    auto h = f.read_header();
+    REQUIRE(h.framenum == 19+1);
+
+    //Reading again does not change the file pointer
+    auto h2 = f.read_header();
+    REQUIRE(h2.framenum == 19+1);
+
+    f.seek(59);
+    REQUIRE(f.tell() == 59);
+    auto h3 = f.read_header();
+    REQUIRE(h3.framenum == 59+1);
+
+    JungfrauDataHeader h4;
+    aare::NDArray<uint16_t> image(f.shape());
+    f.read_into(&image, &h4);
+    REQUIRE(h4.framenum == 59+1);
+
+    //now we should be on the next frame
+    REQUIRE(f.tell() == 60);
+    REQUIRE(f.read_header().framenum == 60+1);
+
+    REQUIRE_THROWS(f.seek(86356)); //out of range
+}
+
+TEST_CASE("Open a Jungfrau data file with non zero file index", "[.files]"){
+
+    auto fpath = test_data_path() / "dat" / "AldoJF65k_000003.dat";
+    REQUIRE(std::filesystem::exists(fpath));
+
+    JungfrauDataFile f(fpath);
+
+    //18 files per data file, opening the 3rd file we ignore the first 3
+    REQUIRE(f.total_frames() == 113-18*3);
+    REQUIRE(f.tell() == 0);
+
+    //Frame numbers start at 1 in the first file
+    REQUIRE(f.read_header().framenum == 18*3+1);
+
+    // moving relative to the third file
+    f.seek(5);
+    REQUIRE(f.read_header().framenum == 18*3+1+5);
+
+    // ignoring the first 3 files
+    REQUIRE(f.n_files() == 4);
+
+    REQUIRE(f.current_file().stem() == "AldoJF65k_000003");
+
+}
+
+TEST_CASE("Read into throws if size doesn't match", "[.files]"){
+    auto fpath = test_data_path() / "dat" / "AldoJF65k_000000.dat";
+    REQUIRE(std::filesystem::exists(fpath));
+
+    JungfrauDataFile f(fpath);
+
+    aare::NDArray<uint16_t> image({39, 85});
+    JungfrauDataHeader header;
+
+    REQUIRE_THROWS(f.read_into(&image, &header));
+    REQUIRE_THROWS(f.read_into(&image, nullptr));
+    REQUIRE_THROWS(f.read_into(&image));
+
+    REQUIRE(f.tell() == 0);
+
+
+}

src/NDArray.test.cpp

@@ -183,14 +183,14 @@ TEST_CASE("Size and shape matches") {
     int64_t h = 75;
     std::array<int64_t, 2> shape{w, h};
     NDArray<double> a{shape};
-    REQUIRE(a.size() == static_cast<uint64_t>(w * h));
+    REQUIRE(a.size() == w * h);
     REQUIRE(a.shape() == shape);
 }
 
 TEST_CASE("Initial value matches for all elements") {
     double v = 4.35;
     NDArray<double> a{{5, 5}, v};
-    for (uint32_t i = 0; i < a.size(); ++i) {
+    for (int i = 0; i < a.size(); ++i) {
         REQUIRE(a(i) == v);
     }
 }

src/NDView.test.cpp

@@ -3,6 +3,7 @@
 
 #include <iostream>
 #include <vector>
+#include <numeric>
 
 using aare::NDView;
 using aare::Shape;
@@ -21,10 +22,8 @@ TEST_CASE("Element reference 1D") {
 }
 
 TEST_CASE("Element reference 2D") {
-    std::vector<int> vec;
-    for (int i = 0; i != 12; ++i) {
-        vec.push_back(i);
-    }
+    std::vector<int> vec(12);
+    std::iota(vec.begin(), vec.end(), 0);
 
     NDView<int, 2> data(vec.data(), Shape<2>{3, 4});
     REQUIRE(vec.size() == static_cast<size_t>(data.size()));
@@ -58,10 +57,8 @@ TEST_CASE("Element reference 3D") {
 }
 
 TEST_CASE("Plus and miuns with single value") {
-    std::vector<int> vec;
-    for (int i = 0; i != 12; ++i) {
-        vec.push_back(i);
-    }
+    std::vector<int> vec(12);
+    std::iota(vec.begin(), vec.end(), 0);
     NDView<int, 2> data(vec.data(), Shape<2>{3, 4});
     data += 5;
     int i = 0;
@@ -116,10 +113,8 @@ TEST_CASE("elementwise assign") {
 }
 
 TEST_CASE("iterators") {
-    std::vector<int> vec;
-    for (int i = 0; i != 12; ++i) {
-        vec.push_back(i);
-    }
+    std::vector<int> vec(12);
+    std::iota(vec.begin(), vec.end(), 0);
     NDView<int, 1> data(vec.data(), Shape<1>{12});
     int i = 0;
     for (const auto item : data) {
@@ -167,27 +162,31 @@ TEST_CASE("divide with another span") {
 }
 
 TEST_CASE("Retrieve shape") {
-    std::vector<int> vec;
-    for (int i = 0; i != 12; ++i) {
-        vec.push_back(i);
-    }
+    std::vector<int> vec(12);
+    std::iota(vec.begin(), vec.end(), 0);
     NDView<int, 2> data(vec.data(), Shape<2>{3, 4});
     REQUIRE(data.shape()[0] == 3);
     REQUIRE(data.shape()[1] == 4);
 }
 
 TEST_CASE("compare two views") {
-    std::vector<int> vec1;
-    for (int i = 0; i != 12; ++i) {
-        vec1.push_back(i);
-    }
+    std::vector<int> vec1(12);
+    std::iota(vec1.begin(), vec1.end(), 0);
     NDView<int, 2> view1(vec1.data(), Shape<2>{3, 4});
 
-    std::vector<int> vec2;
-    for (int i = 0; i != 12; ++i) {
-        vec2.push_back(i);
-    }
+    std::vector<int> vec2(12);
+    std::iota(vec2.begin(), vec2.end(), 0);
     NDView<int, 2> view2(vec2.data(), Shape<2>{3, 4});
 
     REQUIRE((view1 == view2));
+}
+
+
+TEST_CASE("Create a view over a vector"){
+    std::vector<int> vec(12);
+    std::iota(vec.begin(), vec.end(), 0);
+    auto v = aare::make_view(vec);
+    REQUIRE(v.shape()[0] == 12);
+    REQUIRE(v[0] == 0);
+    REQUIRE(v[11] == 11);
 }

src/RawSubFile.cpp

@@ -1,9 +1,12 @@
 #include "aare/RawSubFile.hpp"
 #include "aare/PixelMap.hpp"
+#include "aare/utils/ifstream_helpers.hpp"
 #include <cstring> // memcpy
 #include <fmt/core.h>
 #include <iostream>
 
+
+
 namespace aare {
 
 RawSubFile::RawSubFile(const std::filesystem::path &fname,
@@ -20,7 +23,7 @@ RawSubFile::RawSubFile(const std::filesystem::path &fname,
     }
 
     if (std::filesystem::exists(fname)) {
-        n_frames = std::filesystem::file_size(fname) /
+        m_num_frames = std::filesystem::file_size(fname) /
                    (sizeof(DetectorHeader) + rows * cols * bitdepth / 8);
     } else {
         throw std::runtime_error(
@@ -35,7 +38,7 @@ RawSubFile::RawSubFile(const std::filesystem::path &fname,
     }
 
 #ifdef AARE_VERBOSE
-    fmt::print("Opened file: {} with {} frames\n", m_fname.string(), n_frames);
+    fmt::print("Opened file: {} with {} frames\n", m_fname.string(), m_num_frames);
     fmt::print("m_rows: {}, m_cols: {}, m_bitdepth: {}\n", m_rows, m_cols,
                m_bitdepth);
     fmt::print("file size: {}\n", std::filesystem::file_size(fname));
@@ -43,8 +46,8 @@ RawSubFile::RawSubFile(const std::filesystem::path &fname,
 }
 
 void RawSubFile::seek(size_t frame_index) {
-    if (frame_index >= n_frames) {
-        throw std::runtime_error(LOCATION + fmt::format("Frame index {} out of range in a file with {} frames", frame_index, n_frames));
+    if (frame_index >= m_num_frames) {
+        throw std::runtime_error(LOCATION + fmt::format("Frame index {} out of range in a file with {} frames", frame_index, m_num_frames));
     }
     m_file.seekg((sizeof(DetectorHeader) + bytes_per_frame()) * frame_index);
 }
@@ -60,6 +63,10 @@ void RawSubFile::read_into(std::byte *image_buf, DetectorHeader *header) {
         m_file.seekg(sizeof(DetectorHeader), std::ios::cur);
     }
 
+    if (m_file.fail()){
+        throw std::runtime_error(LOCATION + ifstream_error_msg(m_file));
+    }
+
     // TODO! expand support for different bitdepths
     if (m_pixel_map) {
         // read into a temporary buffer and then copy the data to the buffer
@@ -79,8 +86,24 @@ void RawSubFile::read_into(std::byte *image_buf, DetectorHeader *header) {
         // read directly into the buffer
         m_file.read(reinterpret_cast<char *>(image_buf), bytes_per_frame());
     }
+
+    if (m_file.fail()){
+        throw std::runtime_error(LOCATION + ifstream_error_msg(m_file));
+    }
 }
 
+void RawSubFile::read_into(std::byte *image_buf, size_t n_frames, DetectorHeader *header) {
+    for (size_t i = 0; i < n_frames; i++) {
+        read_into(image_buf, header);
+        image_buf += bytes_per_frame();
+        if (header) {
+            ++header;
+        }
+    }
+}
+
+
+
 template <typename T>
 void RawSubFile::read_with_map(std::byte *image_buf) {
     auto part_buffer = new std::byte[bytes_per_frame()];

src/algorithm.test.cpp

@@ -1,11 +1,12 @@
 
 
-#include <aare/algorithm.hpp>
 #include <catch2/catch_test_macros.hpp>
+#include <aare/algorithm.hpp>
+
 
 TEST_CASE("Find the closed index in a 1D array", "[algorithm]") {
     aare::NDArray<double, 1> arr({5});
-    for (size_t i = 0; i < arr.size(); i++) {
+    for (ssize_t i = 0; i < arr.size(); i++) {
         arr[i] = i;
     }
     // arr 0, 1, 2, 3, 4
@@ -16,9 +17,9 @@ TEST_CASE("Find the closed index in a 1D array", "[algorithm]") {
     REQUIRE(aare::nearest_index(arr, -1.0) == 0);
 }
 
-TEST_CASE("Passing integers to nearest_index works", "[algorithm]") {
+TEST_CASE("Passing integers to nearest_index works", "[algorithm]"){
     aare::NDArray<int, 1> arr({5});
-    for (size_t i = 0; i < arr.size(); i++) {
+    for (ssize_t i = 0; i < arr.size(); i++) {
         arr[i] = i;
     }
     // arr 0, 1, 2, 3, 4
@@ -29,7 +30,8 @@ TEST_CASE("Passing integers to nearest_index works", "[algorithm]") {
     REQUIRE(aare::nearest_index(arr, -1) == 0);
 }
 
-TEST_CASE("nearest_index works with std::vector", "[algorithm]") {
+
+TEST_CASE("nearest_index works with std::vector", "[algorithm]"){
     std::vector<double> vec = {0, 1, 2, 3, 4};
     REQUIRE(aare::nearest_index(vec, 2.123) == 2);
     REQUIRE(aare::nearest_index(vec, 2.66) == 3);
@@ -38,7 +40,7 @@ TEST_CASE("nearest_index works with std::vector", "[algorithm]") {
     REQUIRE(aare::nearest_index(vec, -10.0) == 0);
 }
 
-TEST_CASE("nearest index works with std::array", "[algorithm]") {
+TEST_CASE("nearest index works with std::array", "[algorithm]"){
     std::array<double, 5> arr = {0, 1, 2, 3, 4};
     REQUIRE(aare::nearest_index(arr, 2.123) == 2);
     REQUIRE(aare::nearest_index(arr, 2.501) == 3);
@@ -47,9 +49,20 @@ TEST_CASE("nearest index works with std::array", "[algorithm]") {
     REQUIRE(aare::nearest_index(arr, -10.0) == 0);
 }
 
-TEST_CASE("last smaller", "[algorithm]") {
+TEST_CASE("nearest index when there is no different uses the first element", "[algorithm]"){
+    std::vector<int> vec = {5, 5, 5, 5, 5};
+    REQUIRE(aare::nearest_index(vec, 5) == 0);
+}
+
+TEST_CASE("nearest index when there is no different uses the first element also when all smaller", "[algorithm]"){
+    std::vector<int> vec = {5, 5, 5, 5, 5};
+    REQUIRE(aare::nearest_index(vec, 10) == 0);
+}
+
+
+TEST_CASE("last smaller", "[algorithm]"){
     aare::NDArray<double, 1> arr({5});
-    for (size_t i = 0; i < arr.size(); i++) {
+    for (ssize_t i = 0; i < arr.size(); i++) {
         arr[i] = i;
     }
     // arr 0, 1, 2, 3, 4
@@ -59,11 +72,88 @@ TEST_CASE("last smaller", "[algorithm]") {
     REQUIRE(aare::last_smaller(arr, 253.) == 4);
 }
 
-TEST_CASE("returns last bin strictly smaller", "[algorithm]") {
+TEST_CASE("returns last bin strictly smaller", "[algorithm]"){
     aare::NDArray<double, 1> arr({5});
-    for (size_t i = 0; i < arr.size(); i++) {
+    for (ssize_t i = 0; i < arr.size(); i++) {
         arr[i] = i;
     }
     // arr 0, 1, 2, 3, 4
-    REQUIRE(aare::last_smaller(arr, 2.0) == 2);
-}
+    REQUIRE(aare::last_smaller(arr, 2.0) == 1);
+
+}
+
+TEST_CASE("last_smaller with all elements smaller returns last element", "[algorithm]"){
+    aare::NDArray<double, 1> arr({5});
+    for (ssize_t i = 0; i < arr.size(); i++) {
+        arr[i] = i;
+    }
+    // arr 0, 1, 2, 3, 4
+    REQUIRE(aare::last_smaller(arr, 50.) == 4);
+}
+
+TEST_CASE("last_smaller with all elements bigger returns first element", "[algorithm]"){
+    aare::NDArray<double, 1> arr({5});
+    for (ssize_t i = 0; i < arr.size(); i++) {
+        arr[i] = i;
+    }
+    // arr 0, 1, 2, 3, 4
+    REQUIRE(aare::last_smaller(arr, -50.) == 0);
+}
+
+TEST_CASE("last smaller with all elements equal returns the first element", "[algorithm]"){
+    std::vector<int> vec = {5,5,5,5,5,5,5};
+    REQUIRE(aare::last_smaller(vec, 5) == 0);
+}
+
+
+TEST_CASE("first_lager with vector", "[algorithm]"){
+    std::vector<double> vec = {0, 1, 2, 3, 4};
+    REQUIRE(aare::first_larger(vec, 2.5) == 3);
+}
+
+TEST_CASE("first_lager with all elements smaller returns last element", "[algorithm]"){
+    std::vector<double> vec = {0, 1, 2, 3, 4};
+    REQUIRE(aare::first_larger(vec, 50.) == 4);
+}
+
+TEST_CASE("first_lager with all elements bigger returns first element", "[algorithm]"){
+    std::vector<double> vec = {0, 1, 2, 3, 4};
+    REQUIRE(aare::first_larger(vec, -50.) == 0);
+}
+
+TEST_CASE("first_lager with all elements the same as the check returns last", "[algorithm]"){
+    std::vector<int> vec = {14, 14, 14, 14, 14};
+    REQUIRE(aare::first_larger(vec, 14) == 4);
+}
+
+TEST_CASE("first larger with the same element", "[algorithm]"){
+    std::vector<int> vec = {7,8,9,10,11};
+    REQUIRE(aare::first_larger(vec, 9) == 3);
+}
+
+TEST_CASE("cumsum works", "[algorithm]"){
+    std::vector<double> vec = {0, 1, 2, 3, 4};
+    auto result = aare::cumsum(vec);
+    REQUIRE(result.size() == vec.size());
+    REQUIRE(result[0] == 0);
+    REQUIRE(result[1] == 1);
+    REQUIRE(result[2] == 3);
+    REQUIRE(result[3] == 6);
+    REQUIRE(result[4] == 10);
+}
+TEST_CASE("cumsum works with empty vector", "[algorithm]"){
+    std::vector<double> vec = {};
+    auto result = aare::cumsum(vec);
+    REQUIRE(result.size() == 0);
+}
+TEST_CASE("cumsum works with negative numbers", "[algorithm]"){
+    std::vector<double> vec = {0, -1, -2, -3, -4};
+    auto result = aare::cumsum(vec);
+    REQUIRE(result.size() == vec.size());
+    REQUIRE(result[0] == 0);
+    REQUIRE(result[1] == -1);
+    REQUIRE(result[2] == -3);
+    REQUIRE(result[3] == -6);
+    REQUIRE(result[4] == -10);
+}
+

src/decode.cpp

@@ -1,5 +1,5 @@
 #include "aare/decode.hpp"
-
+#include <cmath>
 namespace aare {
 
 uint16_t adc_sar_05_decode64to16(uint64_t input){
@@ -22,6 +22,10 @@ uint16_t adc_sar_05_decode64to16(uint64_t input){
 }
 
 void adc_sar_05_decode64to16(NDView<uint64_t, 2> input, NDView<uint16_t,2> output){
+    if(input.shape() != output.shape()){
+        throw std::invalid_argument(LOCATION + " input and output shapes must match");
+    }
+
     for(int64_t i = 0; i < input.shape(0); i++){
         for(int64_t j = 0; j < input.shape(1); j++){
             output(i,j) = adc_sar_05_decode64to16(input(i,j));
@@ -49,6 +53,9 @@ uint16_t adc_sar_04_decode64to16(uint64_t input){
 }
 
 void adc_sar_04_decode64to16(NDView<uint64_t, 2> input, NDView<uint16_t,2> output){
+    if(input.shape() != output.shape()){
+        throw std::invalid_argument(LOCATION + " input and output shapes must match");
+    }
     for(int64_t i = 0; i < input.shape(0); i++){
         for(int64_t j = 0; j < input.shape(1); j++){
             output(i,j) = adc_sar_04_decode64to16(input(i,j));
@@ -56,6 +63,40 @@ void adc_sar_04_decode64to16(NDView<uint64_t, 2> input, NDView<uint16_t,2> outpu
     }
 }
 
+double apply_custom_weights(uint16_t input, const NDView<double, 1> weights) {
+    if(weights.size() > 16){
+        throw std::invalid_argument("weights size must be less than or equal to 16");
+    }
+
+    double result = 0.0;
+    for (ssize_t i = 0; i < weights.size(); ++i) {
+        result += ((input >> i) & 1) * std::pow(weights[i], i);
+    }
+    return result;
+
+}
+
+void apply_custom_weights(NDView<uint16_t, 1> input, NDView<double, 1> output, const NDView<double,1> weights) {
+    if(input.shape() != output.shape()){
+        throw std::invalid_argument(LOCATION + " input and output shapes must match");
+    }
+
+    //Calculate weights to avoid repeatedly calling std::pow
+    std::vector<double> weights_powers(weights.size());
+    for (ssize_t i = 0; i < weights.size(); ++i) {
+        weights_powers[i] = std::pow(weights[i], i);
+    }
+
+    // Apply custom weights to each element in the input array
+    for (ssize_t i = 0; i < input.shape(0); i++) {            
+        double result = 0.0;
+        for (size_t bit_index = 0; bit_index < weights_powers.size(); ++bit_index) {
+            result += ((input(i) >> bit_index) & 1) * weights_powers[bit_index];
+        }
+        output(i) = result;
+    }
+}
+
 
 
 } // namespace aare

src/decode.test.cpp

@@ -0,0 +1,80 @@
+#include "aare/decode.hpp"
+
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include "aare/NDArray.hpp"
+using Catch::Matchers::WithinAbs;
+#include <vector>
+
+TEST_CASE("test_adc_sar_05_decode64to16"){
+    uint64_t input = 0;
+    uint16_t output = aare::adc_sar_05_decode64to16(input);
+    CHECK(output == 0);
+
+
+    // bit 29 on th input is bit 0 on the output
+    input = 1UL << 29;
+    output = aare::adc_sar_05_decode64to16(input);
+    CHECK(output == 1);
+
+    // test all bits by iteratting through the bitlist
+    std::vector<int> bitlist = {29, 19, 28, 18, 31, 21, 27, 20, 24, 23, 25, 22};
+    for (size_t i = 0; i < bitlist.size(); i++) {
+        input = 1UL << bitlist[i];
+        output = aare::adc_sar_05_decode64to16(input);
+        CHECK(output == (1 << i));
+    }
+
+    
+    // test a few "random" values
+    input = 0;
+    input |= (1UL << 29);
+    input |= (1UL << 19);
+    input |= (1UL << 28);
+    output = aare::adc_sar_05_decode64to16(input);
+    CHECK(output == 7UL);
+
+
+    input = 0;
+    input |= (1UL << 18);
+    input |= (1UL << 27);
+    input |= (1UL << 25);
+    output = aare::adc_sar_05_decode64to16(input);
+    CHECK(output == 1096UL);
+
+    input = 0;
+    input |= (1UL << 25);
+    input |= (1UL << 22);
+    output = aare::adc_sar_05_decode64to16(input);
+    CHECK(output == 3072UL);
+ }
+
+
+ TEST_CASE("test_apply_custom_weights") {
+
+    uint16_t input = 1;
+    aare::NDArray<double, 1> weights_data({3}, 0.0);
+    weights_data(0) = 1.7;
+    weights_data(1) = 2.1;
+    weights_data(2) = 1.8;
+
+    auto weights = weights_data.view();
+
+
+    double output = aare::apply_custom_weights(input, weights);
+    CHECK_THAT(output, WithinAbs(1.0, 0.001));
+
+    input = 1 << 1;
+    output = aare::apply_custom_weights(input, weights);
+    CHECK_THAT(output, WithinAbs(2.1, 0.001));
+
+
+    input = 1 << 2;
+    output = aare::apply_custom_weights(input, weights);
+    CHECK_THAT(output, WithinAbs(3.24, 0.001));
+
+    input = 0b111;
+    output = aare::apply_custom_weights(input, weights);
+    CHECK_THAT(output, WithinAbs(6.34, 0.001));
+
+ }

src/utils/ifstream_helpers.cpp

@@ -0,0 +1,18 @@
+#include "aare/utils/ifstream_helpers.hpp"
+
+namespace aare {
+
+std::string ifstream_error_msg(std::ifstream &ifs) {
+    std::ios_base::iostate state = ifs.rdstate();
+    if (state & std::ios_base::eofbit) {
+        return " End of file reached";
+    } else if (state & std::ios_base::badbit) {
+        return " Bad file stream";
+    } else if (state & std::ios_base::failbit) {
+        return " File read failed";
+    }else{
+        return " Unknown/no error";
+    }
+}
+
+} // namespace aare