Binding 5x5, 7x7 and 9x9 clusters in python (#188 )

- New binding code with macros to bind all cluster templates - Simplified factory function on the python side - 5x5, 7x7 and 9x9 bindings in python
Refactor cluster bindings (#185 )
2025-06-14 08:17:13 +02:00 · 2025-06-05 08:57:59 +02:00 · 2025-06-03 08:43:40 +02:00 · 2025-05-22 11:23:57 +02:00 · 2025-05-22 11:10:23 +02:00 · 2025-05-22 11:00:03 +02:00
48 changed files with 1441 additions and 2173 deletions
--- a/.github/workflows/build_and_deploy_conda.yml
+++ b/.github/workflows/build_and_deploy_conda.yml
@ -1,9 +1,9 @@
 name: Build pkgs and deploy if on main

 on:
-  push:
-    branches:
-     - main
+  release:
+      types:
+        - published

 jobs:
  build:
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@ -62,7 +62,6 @@ option(AARE_FETCH_CATCH "Use FetchContent to download catch2" ON)
 option(AARE_FETCH_JSON "Use FetchContent to download nlohmann::json" ON)
 option(AARE_FETCH_ZMQ "Use FetchContent to download libzmq" ON)
 option(AARE_FETCH_LMFIT "Use FetchContent to download lmfit" ON)
-option(AARE_FETCH_HDF5 "Use FetchContent to download hdf5-devel" OFF)


 #Convenience option to use system libraries only (no FetchContent)
@ -74,13 +73,15 @@ if(AARE_SYSTEM_LIBRARIES)
    set(AARE_FETCH_CATCH OFF CACHE BOOL "Disabled FetchContent for catch2" FORCE)
    set(AARE_FETCH_JSON OFF CACHE BOOL "Disabled FetchContent for nlohmann::json" FORCE)
    set(AARE_FETCH_ZMQ OFF CACHE BOOL "Disabled FetchContent for libzmq" FORCE)
-    set(AARE_FETCH_HDF5 OFF CACHE BOOL "Disabled FetchContent for hdf5" FORCE)
    # Still fetch lmfit when setting AARE_SYSTEM_LIBRARIES since this is not available
    # on conda-forge
 endif()

 if(AARE_VERBOSE)
    add_compile_definitions(AARE_VERBOSE)
+    add_compile_definitions(AARE_LOG_LEVEL=aare::logDEBUG5)
+else()
+    add_compile_definitions(AARE_LOG_LEVEL=aare::logERROR)
 endif()

 if(AARE_CUSTOM_ASSERT)
@ -92,6 +93,7 @@ if(AARE_BENCHMARKS)
 endif()


+
 set(CMAKE_EXPORT_COMPILE_COMMANDS ON)

 if(AARE_FETCH_LMFIT)
@ -221,23 +223,6 @@ else()
    find_package(nlohmann_json 3.11.3 REQUIRED)
 endif()

-if(AARE_FETCH_HDF5)
-    message(FATAL_ERROR "Fetching HDF5 via FetchContent is not supported here. Please install it via your system. 
-    For Ubuntu: sudo apt install libhdf5-dev
-    For Red Hat: sudo dnf install hdf5-devel
-    For MacOS: brew install hdf5")
-else()
-    find_package(HDF5 QUIET COMPONENTS CXX)
-    if (HDF5_FOUND)
-        message(STATUS "Found HDF5: ${HDF5_INCLUDE_DIRS}")
-    else()
-        message(FATAL_ERROR "HDF5 was NOT found! Please install it via your system  
-            For Ubuntu: sudo apt install libhdf5-dev
-            For Red Hat: sudo dnf install hdf5-devel
-            For MacOS: brew install hdf5")
-    endif()
-endif()
-
 include(GNUInstallDirs)

 # If conda build, always set lib dir to 'lib'
@ -350,6 +335,10 @@ if(AARE_ASAN)
    )
 endif()

+
+
+
+
 if(AARE_TESTS)
    enable_testing()
    add_subdirectory(tests)
@ -359,7 +348,6 @@ endif()
 ###------------------------------------------------------------------------------------------

 set(PUBLICHEADERS
-    include/aare/AngleCalibration.hpp
    include/aare/ArrayExpr.hpp
    include/aare/CalculateEta.hpp
    include/aare/Cluster.hpp
@ -374,14 +362,10 @@ set(PUBLICHEADERS
    include/aare/Fit.hpp
    include/aare/FileInterface.hpp
    include/aare/FilePtr.hpp
-    include/aare/FlatField.hpp
    include/aare/Frame.hpp
    include/aare/GainMap.hpp
    include/aare/geo_helpers.hpp
-    include/aare/Hdf5FileReader.hpp
    include/aare/JungfrauDataFile.hpp
-    include/aare/MythenDetectorSpecifications.hpp
-    include/aare/MythenFileReader.hpp
    include/aare/NDArray.hpp
    include/aare/NDView.hpp
    include/aare/NumpyFile.hpp
@ -397,7 +381,6 @@ set(PUBLICHEADERS


 set(SourceFiles
-    ${CMAKE_CURRENT_SOURCE_DIR}/src/AngleCalibration.cpp
    ${CMAKE_CURRENT_SOURCE_DIR}/src/CtbRawFile.cpp
    ${CMAKE_CURRENT_SOURCE_DIR}/src/defs.cpp
    ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.cpp
@ -415,7 +398,6 @@ 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
 )
@ -434,7 +416,6 @@ target_link_libraries(
    PUBLIC 
    fmt::fmt 
    nlohmann_json::nlohmann_json 
-    HDF5::HDF5
    ${STD_FS_LIB} # from helpers.cmake
    PRIVATE 
    aare_compiler_flags 
@ -455,14 +436,12 @@ endif()
 if(AARE_TESTS)
    set(TestSources
            ${CMAKE_CURRENT_SOURCE_DIR}/src/algorithm.test.cpp
-            ${CMAKE_CURRENT_SOURCE_DIR}/src/AngleCalibration.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
            ${CMAKE_CURRENT_SOURCE_DIR}/src/RawMasterFile.test.cpp
-            ${CMAKE_CURRENT_SOURCE_DIR}/src/Hdf5FileReader.test.cpp
            ${CMAKE_CURRENT_SOURCE_DIR}/src/NDArray.test.cpp
            ${CMAKE_CURRENT_SOURCE_DIR}/src/NDView.test.cpp
            ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFinder.test.cpp
@ -473,10 +452,10 @@ if(AARE_TESTS)
            ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFinderMT.test.cpp
            ${CMAKE_CURRENT_SOURCE_DIR}/src/Pedestal.test.cpp
            ${CMAKE_CURRENT_SOURCE_DIR}/src/JungfrauDataFile.test.cpp
-            ${CMAKE_CURRENT_SOURCE_DIR}/src/MythenFileReader.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
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/RawSubFile.test.cpp
            ${CMAKE_CURRENT_SOURCE_DIR}/src/utils/task.test.cpp

        )
--- a/2
+++ b/2
@ -1 +1 @@
-0.0.0
+2025.5.22
--- a/include/aare/AngleCalibration.hpp
+++ b/include/aare/AngleCalibration.hpp
@ -1,164 +0,0 @@
-#pragma once
-#include <algorithm>
-#include <cmath>
-#include <cstdint>
-
-#include <fstream>
-#include <iomanip>
-#include <iostream>
-#include <memory>
-#include <sstream>
-#include <string>
-#include <vector>
-
-#include "FlatField.hpp"
-#include "MythenDetectorSpecifications.hpp"
-#include "MythenFileReader.hpp"
-#include "NDArray.hpp"
-
-namespace aare {
-
-using parameters =
-    std::tuple<std::vector<double>, std::vector<double>, std::vector<double>>;
-
-class AngleCalibration {
-
-  public:
-    AngleCalibration(
-        std::shared_ptr<MythenDetectorSpecifications> mythen_detector_,
-        std::shared_ptr<FlatField> flat_field_,
-        std::shared_ptr<MythenFileReader> mythen_file_reader_);
-
-    /** set the histogram bin width [degrees] */
-    void set_histogram_bin_width(double bin_width);
-
-    double get_histogram_bin_width() const;
-
-    ssize_t get_new_num_bins() const;
-
-    /** reads the historical Detector Group (DG) parameters from file **/
-    void read_initial_calibration_from_file(const std::string &filename);
-
-    std::vector<double> get_centers() const;
-    std::vector<double> get_conversions() const;
-
-    std::vector<double> get_offsets() const;
-
-    NDView<double, 1> get_new_photon_counts() const;
-
-    NDView<double, 1> get_new_statistical_errors() const;
-
-    /** converts DG parameters to easy EE parameters e.g.geometric
-     * parameters */
-    parameters convert_to_EE_parameters() const;
-
-    std::tuple<double, double, double>
-    convert_to_EE_parameters(const size_t module_index) const;
-
-    std::tuple<double, double, double>
-    convert_to_EE_parameters(const double center, const double conversion,
-                             const double offset) const;
-
-    /** converts DG parameters to BC parameters e.g. best computing
-     * parameters */
-    parameters convert_to_BC_parameters() const;
-
-    /**
-     * calculates new histogram with fixed sized angle bins
-     * for several acquisitions at different detector angles for given frame
-     * indices
-     * @param start_frame_index, end_frame_index gives range of frames
-     */
-    void
-    calculate_fixed_bin_angle_width_histogram(const size_t start_frame_index,
-                                              const size_t end_frame_index);
-
-    void write_to_file(const std::string &filename,
-                       const bool store_nonzero_bins = false,
-                       const std::filesystem::path &filepath =
-                           std::filesystem::current_path()) const;
-
-    /** calculates diffraction angle from EE module parameters (used in
-     * Beer's Law)
-     * @param strip_index local strip index of module
-     */
-    double diffraction_angle_from_EE_parameters(
-        const double module_center_distance, const double normal_distance,
-        const double angle, const size_t strip_index,
-        const double distance_to_strip = 0) const;
-
-    /** calculates diffraction angle from EE module parameters (used in
-     * Beer's Law)
-     * @param center module center
-     * @param conversion module conversion
-     * @param offset module offset
-     * @param strip_index local strip index of module
-     * @param distance_to_strip distance to strip given by strip_index and
-     * module -> note needs to be small enough to be in the respective module
-     */
-    double diffraction_angle_from_DG_parameters(
-        const double center, const double conversion, const double offset,
-        const size_t strip_index, const double distance_to_strip = 0) const;
-
-    /** calculated the strip width expressed as angle [degrees]
-     * @param strip_index local strip index of module
-     */
-    double angular_strip_width_from_DG_parameters(
-        const double center, const double conversion, const double offset,
-        const size_t local_strip_index) const;
-
-    double angular_strip_width_from_EE_parameters(
-        const double module_center_distance, const double normal_distance,
-        const double angle, const size_t local_strip_index) const;
-
-  protected:
-    /** converts global strip index to local strip index of that module */
-    size_t global_to_local_strip_index_conversion(
-        const size_t global_strip_index) const;
-
-    /**
-     * redistributes photon counts with of histogram using one bin per strip
-     * to histogram with fixed size angle bins
-     * @param frame MythenFrame storing data from image
-     * @param bin_counts accumulate new photon counts
-     * @param new_statistical_weights accumulate new statistical weights
-     * @param new_errors accumulate new_errors
-     */
-    void redistribute_photon_counts_to_fixed_angle_bins(
-        const MythenFrame &frame, NDView<double, 1> bin_counts,
-        NDView<double, 1> new_statistical_weights, NDView<double, 1> new_errors,
-        NDView<double, 1> inverse_nromalized_flatfield) const;
-
-  private:
-    // TODO: Design maybe have a struct with three vectors, store all three
-    // sets of parameters as member variables
-
-    // TODO: check if interpretation and units are correct
-    // historical DG parameters
-    // TODO change to NDArray
-    std::vector<double> centers;     // orthogonal projection of sample onto
-                                     // detector (given in strip number) [mm]
-                                     // D/pitch
-    std::vector<double> conversions; // pitch/(normal distance from sample
-                                     // to detector (R)) [mm]
-                                     // //used for easy conversion
-    std::vector<double>
-        offsets; // position of strip zero relative to sample [degrees] phi
-                 // - 180/pi*D/R TODO: expected an arcsin(D/R)?
-
-    std::shared_ptr<MythenDetectorSpecifications> mythen_detector;
-
-    std::shared_ptr<FlatField> flat_field;
-
-    NDArray<double, 1> new_photon_counts;
-    NDArray<double, 1> new_photon_count_errors;
-
-    double histogram_bin_width = 0.0036; // [degrees]
-
-    ssize_t num_bins{};
-
-    std::shared_ptr<MythenFileReader>
-        mythen_file_reader; // TODO replace by FileInterface ptr
-};
-
-} // namespace aare
--- a/include/aare/ClusterFile.hpp
+++ b/include/aare/ClusterFile.hpp
@ -5,6 +5,8 @@
 #include "aare/GainMap.hpp"
 #include "aare/NDArray.hpp"
 #include "aare/defs.hpp"
+#include "aare/logger.hpp"
+
 #include <filesystem>
 #include <fstream>
 #include <optional>
@ -369,11 +371,15 @@ ClusterFile<ClusterType, Enable>::read_frame_without_cut() {
                                 "Could not read number of clusters");
    }

+    LOG(logDEBUG1) << "Reading " << n_clusters
+                  << " clusters from frame " << frame_number;
+
    ClusterVector<ClusterType> clusters(n_clusters);
    clusters.set_frame_number(frame_number);
-
    clusters.resize(n_clusters);

+    LOG(logDEBUG1) << "clusters.item_size(): " << clusters.item_size();
+
    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");
--- a/include/aare/ClusterFileSink.hpp
+++ b/include/aare/ClusterFileSink.hpp
@ -21,7 +21,7 @@ class ClusterFileSink {

    void process() {
        m_stopped = false;
-        fmt::print("ClusterFileSink started\n");
+        LOG(logDEBUG) << "ClusterFileSink started";
        while (!m_stop_requested || !m_source->isEmpty()) {
            if (ClusterVector<ClusterType> *clusters = m_source->frontPtr();
                clusters != nullptr) {
@ -41,13 +41,16 @@ class ClusterFileSink {
                std::this_thread::sleep_for(m_default_wait);
            }
        }
-        fmt::print("ClusterFileSink stopped\n");
+        LOG(logDEBUG) << "ClusterFileSink stopped";
        m_stopped = true;
    }

  public:
    ClusterFileSink(ClusterFinderMT<ClusterType, uint16_t, double> *source,
                    const std::filesystem::path &fname) {
+        LOG(logDEBUG) << "ClusterFileSink: "
+                       << "source: " << source->sink()
+                       << ", file: " << fname.string();
        m_source = source->sink();
        m_thread = std::thread(&ClusterFileSink::process, this);
        m_file.open(fname, std::ios::binary);
--- a/include/aare/ClusterFinder.hpp
+++ b/include/aare/ClusterFinder.hpp
@ -38,7 +38,12 @@ class ClusterFinder {
        : 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) {};
+          m_pedestal(image_size[0], image_size[1]), m_clusters(capacity) {
+        LOG(logDEBUG ) << "ClusterFinder: "
+                       << "image_size: " << image_size[0] << "x" << image_size[1]
+                       << ", nSigma: " << nSigma
+                       << ", capacity: " << capacity;
+    }

    void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) {
        m_pedestal.push(frame);
--- a/include/aare/ClusterFinderMT.hpp
+++ b/include/aare/ClusterFinderMT.hpp
@ -7,6 +7,7 @@

 #include "aare/ClusterFinder.hpp"
 #include "aare/NDArray.hpp"
+#include "aare/logger.hpp"
 #include "aare/ProducerConsumerQueue.hpp"

 namespace aare {
@ -123,6 +124,12 @@ class ClusterFinderMT {
                    size_t capacity = 2000, size_t n_threads = 3)
        : m_n_threads(n_threads) {

+        LOG(logDEBUG1) << "ClusterFinderMT: "
+            << "image_size: " << image_size[0] << "x" << image_size[1]
+            << ", nSigma: " << nSigma
+            << ", capacity: " << capacity
+            << ", n_threads: " << n_threads;
+
        for (size_t i = 0; i < n_threads; i++) {
            m_cluster_finders.push_back(
                std::make_unique<
--- a/include/aare/ClusterVector.hpp
+++ b/include/aare/ClusterVector.hpp
@ -133,9 +133,9 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
     */
    size_t capacity() const { return m_data.capacity(); }

-    const auto begin() const { return m_data.begin(); }
+    auto begin() const { return m_data.begin(); }

-    const auto end() const { return m_data.end(); }
+    auto end() const { return m_data.end(); }

    /**
     * @brief Return the size in bytes of a single cluster
--- a/include/aare/FlatField.hpp
+++ b/include/aare/FlatField.hpp
@ -1,112 +0,0 @@
-
-/**
- * stores flatfield for angle calibration
- */
-
-#pragma once
-#include <cmath>
-#include <cstdint>
-
-#include <fstream>
-#include <iostream>
-#include <memory>
-#include <sstream>
-#include <string>
-
-#include "MythenDetectorSpecifications.hpp"
-#include "NDArray.hpp"
-
-namespace aare {
-// TODO maybe template now its uint32
-class FlatField {
-
-  public:
-    FlatField(std::shared_ptr<MythenDetectorSpecifications> mythen_detector_)
-        : mythen_detector(mythen_detector_) {
-
-        flat_field = NDArray<uint32_t, 1>(
-            std::array<ssize_t, 1>{mythen_detector->num_strips()}, 0);
-    }
-
-    void read_flatfield_from_file(const std::string &filename) {
-
-        std::string word;
-        uint32_t strip_number{};
-
-        try {
-            std::ifstream file(filename, std::ios_base::in);
-            if (!file.good()) {
-                throw std::logic_error("file does not exist");
-            }
-
-            std::stringstream file_buffer;
-            file_buffer << file.rdbuf();
-
-            while (file_buffer >> word) {
-
-                strip_number = std::stoi(word);
-
-                file_buffer >> word;
-                if (!mythen_detector->get_bad_channels()[strip_number])
-                    flat_field[strip_number] = std::stod(word);
-            }
-
-            file.close();
-        } catch (const std::exception &e) {
-            std::cerr << "Error: " << e.what() << std::endl;
-        }
-    }
-
-    NDView<uint32_t, 1> get_flatfield() const { return flat_field.view(); }
-
-    double calculate_mean(double tolerance = 0.001) const {
-        auto [sum, count] = std::accumulate(
-            flat_field.begin(), flat_field.end(),
-            std::make_pair<double, ssize_t>(0.0, 0),
-            [&tolerance](std::pair<double, ssize_t> acc, const auto &element) {
-                return element == 0 ? acc
-                                    : std::make_pair(acc.first + element,
-                                                     acc.second + 1);
-            });
-
-        return sum / count;
-    }
-
-    NDArray<double, 1>
-    inverse_normalized_flatfield(double tolerance = 0.001) const {
-        double mean = calculate_mean(tolerance);
-
-        NDArray<double, 1> inverse_normalized_flatfield(flat_field.shape());
-
-        for (ssize_t i = 0; i < flat_field.size(); ++i) {
-            inverse_normalized_flatfield[i] =
-                (flat_field[i] <= tolerance ? 0.0 : mean / flat_field[i]);
-            if (inverse_normalized_flatfield[i] < tolerance)
-                mythen_detector->get_bad_channels()[i] = true;
-        }
-
-        return inverse_normalized_flatfield; // TODO: better to have a copy in
-                                             // this context but unneccessary
-                                             // for angle calibration code
-        // maybe provide inplace and copy option
-        // maybe store as member variable access with view
-    }
-
-    NDArray<double, 1> normalized_flatfield(double tolerance = 0.001) const {
-        double mean = calculate_mean(tolerance);
-
-        NDArray<double, 1> normalized_flatfield(flat_field.shape());
-
-        for (ssize_t i = 0; i < flat_field.size(); ++i) {
-            normalized_flatfield[i] = (flat_field[i] == flat_field[i] / mean);
-            if (normalized_flatfield[i] < tolerance)
-                mythen_detector->get_bad_channels()[i] = true;
-        }
-        return normalized_flatfield;
-    }
-
-  private:
-    NDArray<uint32_t, 1> flat_field; // TODO: should be 2d
-    std::shared_ptr<MythenDetectorSpecifications> mythen_detector;
-};
-} // namespace aare
--- a/include/aare/Hdf5FileReader.hpp
+++ b/include/aare/Hdf5FileReader.hpp
@ -1,212 +0,0 @@
-/************************************************
- * @file HD5FFileReader.hpp
- * @short HDF5FileReader based on H5File object
- ***********************************************/
-
-#pragma once
-
-#include "Frame.hpp"
-#include "NDArray.hpp"
-#include <H5Cpp.h>
-#include <array>
-#include <cxxabi.h>
-#include <optional>
-
-namespace aare {
-
-// return std::type_info
-inline const std::type_info &deduce_cpp_type(const H5::DataType datatype) {
-    if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_UINT8.getId())) {
-        return typeid(uint8_t);
-    } else if (H5Tequal(datatype.getId(),
-                        H5::PredType::NATIVE_UINT16.getId())) {
-        return typeid(uint16_t);
-    } else if (H5Tequal(datatype.getId(),
-                        H5::PredType::NATIVE_UINT32.getId())) {
-        return typeid(uint32_t);
-    } else if (H5Tequal(datatype.getId(),
-                        H5::PredType::NATIVE_UINT64.getId())) {
-        return typeid(uint64_t);
-    } else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT8.getId())) {
-        return typeid(int8_t);
-    } else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT16.getId())) {
-        return typeid(int16_t);
-    } else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT32.getId())) {
-        return typeid(int32_t);
-    } else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT64.getId())) {
-        return typeid(int64_t);
-    } else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT.getId())) {
-        return typeid(int);
-    } else if (H5Tequal(datatype.getId(), H5::PredType::IEEE_F64LE.getId())) {
-        return typeid(double);
-    } else if (H5Tequal(datatype.getId(), H5::PredType::IEEE_F32LE.getId())) {
-        return typeid(float);
-    } else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_FLOAT.getId())) {
-        return typeid(float);
-    } else if (H5Tequal(datatype.getId(),
-                        H5::PredType::NATIVE_DOUBLE.getId())) {
-        return typeid(float);
-    } else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_CHAR.getId()) &&
-               datatype.getId() == H5::PredType::NATIVE_CHAR.getId()) {
-        return typeid(char);
-    } else {
-        throw std::runtime_error("c++ type cannot be deduced");
-    }
-}
-
-struct Subset {
-    std::vector<hsize_t> shape;
-    std::vector<hsize_t> offset; // index where data subset should start
-};
-
-class HDF5Dataset {
-
-  public:
-    HDF5Dataset(const std::string &datasetname_, const H5::DataSet dataset_)
-        : datasetname(datasetname_), dataset(dataset_) {
-        datatype = dataset.getDataType();
-
-        cpp_type = &deduce_cpp_type(datatype);
-
-        dataspace = dataset.getSpace();
-        rank = dataspace.getSimpleExtentNdims(); // number of dimensions
-
-        shape.resize(rank);
-        dataspace.getSimpleExtentDims(shape.data(), nullptr);
-    }
-
-    hsize_t get_shape(ssize_t index) const { return shape[index]; }
-
-    std::vector<hsize_t> get_shape() const { return shape; }
-
-    H5::DataType get_datatype() const { return datatype; }
-
-    const std::type_info *get_cpp_type() const { return cpp_type; }
-
-    /**
-     * Reads subset of dataset into the buffer
-     * e.g. to read one 2d frame pass Subset({shape[1], shape[2]}, {frame_index,
-     * 0,0})
-     */
-    void
-    read_into_buffer(std::byte *buffer,
-                     std::optional<const Subset> subset = std::nullopt) const {
-
-        if (subset) {
-            // TODO treat scalar cases
-            if (static_cast<ssize_t>(subset->offset.size()) != rank) {
-                throw std::runtime_error("provide an offset for" +
-                                         std::to_string(rank) + "dimensions");
-            }
-            for (ssize_t i = 0; i < rank; ++i) {
-                hsize_t size =
-                    i < rank - static_cast<ssize_t>(subset->shape.size())
-                        ? 0
-                        : subset->shape[i - (rank - subset->shape.size())];
-                if ((size + subset->offset[i]) > shape[i]) {
-                    throw std::runtime_error(
-                        "subset is too large or offset is out of bounds");
-                }
-            }
-
-            H5::DataSpace memspace(static_cast<int>(subset->shape.size()),
-                                   subset->shape.data());
-
-            dataspace.selectHyperslab(H5S_SELECT_SET, subset->shape.data(),
-                                      subset->offset.data());
-            dataset.read(buffer, datatype, memspace, dataspace);
-        } else {
-            dataset.read(buffer, datatype);
-        }
-    }
-
-    Frame store_as_frame() const {
-        uint32_t rows{}, cols{};
-        if (rank == 1) {
-            rows = 1;
-            // TODO overflow
-            cols = static_cast<uint32_t>(shape[0]);
-        } else if (rank == 2) {
-            rows = static_cast<uint32_t>(shape[0]);
-            cols = static_cast<uint32_t>(shape[1]);
-        } else {
-            throw std::runtime_error("Frame only supports 2d images");
-        }
-
-        Frame frame(rows, cols, Dtype(*cpp_type));
-
-        read_into_buffer(frame.data());
-
-        return frame;
-    }
-
-    template <typename T, ssize_t NDim>
-    NDArray<T, NDim> store_as_ndarray() const {
-        if (NDim != rank) {
-            std::cout
-                << "Warning: dataset dimension and array dimension mismatch"
-                << std::endl; // TODO: replace with log - still valid if we
-                              // want subset
-        }
-        if (typeid(T) != *cpp_type) {
-            std::cout << "Warning: dataset and array type mismatch"
-                      << std::endl;
-        }
-        std::array<ssize_t, NDim> array_shape{};
-        std::transform(
-            shape.begin(), shape.end(), array_shape.begin(),
-            [](const auto dim) { return static_cast<ssize_t>(dim); });
-
-        aare::NDArray<T, NDim> dataset_array(array_shape);
-
-        read_into_buffer(reinterpret_cast<std::byte *>(dataset_array.data()));
-
-        return dataset_array;
-    }
-
-    // getMemDataSize()
-
-  private:
-    std::string datasetname{};
-    H5::DataSet dataset;
-    H5::DataSpace dataspace;
-    H5::DataType datatype;
-    const std::type_info *cpp_type;
-    ssize_t rank{};
-    std::vector<hsize_t> shape{};
-};
-
-class HDF5FileReader {
-
-  public:
-    HDF5FileReader() = default;
-
-    void open_file(const std::string &filename_) {
-        filename = filename_;
-        try {
-            file = H5::H5File(filename, H5F_ACC_RDONLY);
-        } catch (H5::Exception &e) {
-            std::cerr << "Error: " << e.getDetailMsg() << std::endl;
-        }
-    }
-
-    void close_file() { file.close(); }
-
-    HDF5Dataset get_dataset(const std::string &dataset_name) const {
-        H5::DataSet dataset;
-        try {
-            dataset = file.openDataSet(dataset_name);
-        } catch (H5::Exception &e) {
-            std::cerr << "Error: " << e.getDetailMsg() << std::endl;
-        }
-
-        // TODO use optional to handle error
-        return HDF5Dataset(dataset_name, dataset);
-    }
-
-  private:
-    std::string filename{};
-    H5::H5File file;
-};
-
-} // namespace aare
--- a/include/aare/Interpolator.hpp
+++ b/include/aare/Interpolator.hpp
@ -51,7 +51,7 @@ Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) {
            Photon photon;
            photon.x = cluster.x;
            photon.y = cluster.y;
-            photon.energy = eta.sum;
+            photon.energy = static_cast<decltype(photon.energy)>(eta.sum);

            // auto ie = nearest_index(m_energy_bins, photon.energy)-1;
            // auto ix = nearest_index(m_etabinsx, eta.x)-1;
@ -99,7 +99,7 @@ Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) {
            Photon photon;
            photon.x = cluster.x;
            photon.y = cluster.y;
-            photon.energy = eta.sum;
+            photon.energy = static_cast<decltype(photon.energy)>(eta.sum);

            // Now do some actual interpolation.
            // Find which energy bin the cluster is in
--- a/include/aare/MythenDetectorSpecifications.hpp
+++ b/include/aare/MythenDetectorSpecifications.hpp
@ -1,150 +0,0 @@
-#pragma once
-#include <cmath>
-#include <cstdint>
-
-#include <fstream>
-#include <iostream>
-#include <sstream>
-#include <string>
-
-#include "NDArray.hpp"
-
-namespace aare {
-
-class MythenDetectorSpecifications {
-
-  public:
-    // TODO: constructor that reads from a config file
-
-    MythenDetectorSpecifications() {
-        num_strips_ = max_modules_ * strips_per_module_;
-
-        num_connected_modules_ = max_modules_;
-
-        bad_channels =
-            NDArray<bool, 1>(std::array<ssize_t, 1>{num_strips_}, false);
-
-        connected_modules = NDArray<bool, 1>(
-            std::array<ssize_t, 1>{static_cast<ssize_t>(max_modules_)}, true);
-    }
-
-    MythenDetectorSpecifications(const size_t max_modules,
-                                 const double exposure_time,
-                                 const double bloffset)
-        : max_modules_(max_modules), exposure_time_(exposure_time),
-          bloffset_(bloffset) {
-        num_strips_ = max_modules_ * strips_per_module_;
-
-        num_connected_modules_ = max_modules_;
-
-        bad_channels =
-            NDArray<bool, 1>(std::array<ssize_t, 1>{num_strips_}, false);
-
-        connected_modules = NDArray<bool, 1>(
-            std::array<ssize_t, 1>{static_cast<ssize_t>(max_modules_)}, true);
-    }
-
-    void read_bad_channels_from_file(const std::string &filename) {
-        std::string line;
-
-        try {
-            std::ifstream file(filename, std::ios_base::in);
-            if (!file.good()) {
-                throw std::logic_error("file does not exist");
-            }
-
-            while (std::getline(file, line)) {
-                std::size_t pos = line.find("-");
-
-                if (pos == std::string::npos) {
-                    bad_channels(std::stoi(line)) = true;
-                } else {
-                    size_t line_size = line.size();
-                    for (int i = std::stoi(line.substr(0, pos));
-                         i <= std::stoi(line.substr(pos + 1, line_size - pos));
-                         ++i)
-                        bad_channels(i) = true;
-                }
-            }
-
-            file.close();
-        } catch (const std::exception &e) {
-            std::cerr << "Error: " << e.what() << std::endl;
-        }
-    }
-
-    void read_unconnected_modules_from_file(const std::string &filename) {
-        std::string line;
-
-        try {
-            std::ifstream file(filename, std::ios_base::in);
-            if (!file.good()) {
-                throw std::logic_error("file does not exist");
-            }
-
-            std::stringstream file_buffer;
-            file_buffer << file.rdbuf();
-
-            file_buffer >> line;
-            num_connected_modules_ -= std::stoi(line);
-
-            while (file_buffer >> line) {
-                size_t module = std::stoi(line);
-                connected_modules[module] = false;
-                for (size_t i = module * strips_per_module_;
-                     i < (module + 1) * strips_per_module_; ++i)
-                    bad_channels[i] = true;
-            }
-        } catch (const std::exception &e) {
-            std::cerr << "Error: " << e.what() << std::endl;
-        }
-    }
-
-    NDView<bool, 1> get_bad_channels() const { return bad_channels.view(); }
-
-    NDView<bool, 1> get_connected_modules() const {
-        return connected_modules.view();
-    }
-
-    static constexpr double pitch() { return pitch_; }
-
-    static constexpr size_t strips_per_module() { return strips_per_module_; }
-
-    size_t max_modules() const { return max_modules_; }
-
-    double exposure_time() const { return exposure_time_; }
-
-    double bloffset() const { return bloffset_; }
-
-    double dtt0() const { return dtt0_; }
-
-    static constexpr double min_angle() { return min_angle_; }
-
-    static constexpr double max_angle() { return max_angle_; }
-
-    ssize_t num_strips() const { return num_strips_; }
-
-  private:
-    static constexpr size_t strips_per_module_ = 1280;
-    static constexpr double pitch_ = 0.05; // strip width [mm]
-    static constexpr double min_angle_ =
-        -180.0; // maybe shoudnt be static but configurable
-    static constexpr double max_angle_ = 180.0;
-    static constexpr double dtt0_ =
-        0.0; // No idea what this is - probably configurable
-
-    size_t max_modules_ = 48;
-
-    double exposure_time_ = 5.0; // TODO: could read from acquired file but
-                                 // maybe should be configurable
-    double bloffset_ = 1.532; // what is this? detector offset relative to what?
-
-    size_t num_connected_modules_{};
-
-    ssize_t num_strips_{};
-
-    NDArray<bool, 1> bad_channels;
-    NDArray<bool, 1> connected_modules; // connected modules
-};
-
-} // namespace aare
--- a/include/aare/MythenFileReader.hpp
+++ b/include/aare/MythenFileReader.hpp
@ -1,82 +0,0 @@
-/************************************************
- * @file MythenFileReader.hpp
- * @short minimal file reader to read mythen files
- ***********************************************/
-
-#include <bitset>
-#include <filesystem>
-#include <string>
-
-#include "Hdf5FileReader.hpp"
-#include "NDArray.hpp"
-
-namespace aare {
-
-struct MythenFrame {
-    NDArray<uint32_t, 1> photon_counts;
-    double detector_angle{};
-    // double reference_intensity{}; not needed
-    std::array<uint8_t, 3> channel_mask{};
-};
-
-/** minimal version for a mythen file reader */
-class MythenFileReader : public HDF5FileReader {
-
-  public:
-    MythenFileReader(const std::filesystem::path &file_path_,
-                     const std::string &file_prefix_)
-        : m_base_path(file_path_), file_prefix(file_prefix_) {};
-
-    MythenFrame read_frame(ssize_t frame_index) {
-        // TODO not a good design fixed number of digits in file name for frame
-        // number -> pad with zeros
-        //  not even sure if files have the same name
-        std::string current_file_name =
-            m_base_path / (file_prefix + std::to_string(frame_index) + ".h5");
-
-        MythenFrame frame;
-        open_file(current_file_name);
-
-        auto dataset_photon_count =
-            get_dataset("/entry/instrument/detector/data");
-
-        frame.photon_counts =
-            dataset_photon_count.store_as_ndarray<uint32_t, 1>();
-
-        ++frame.photon_counts; // Why though?
-
-        auto dataset_detector_angle =
-            get_dataset("/entry/instrument/NDAttributes/DetectorAngle");
-
-        dataset_detector_angle.read_into_buffer(
-            reinterpret_cast<std::byte *>(&frame.detector_angle));
-
-        auto dataset_channel_number =
-            get_dataset("/entry/instrument/NDAttributes/CounterMask");
-
-        uint8_t channel_number;
-
-        dataset_channel_number.read_into_buffer(
-            reinterpret_cast<std::byte *>(&channel_number));
-
-        std::bitset<3> binary_channel_numbers(channel_number); // 1 0 0
-
-        // binary_channel_numbers.flip(); // TODO not sure where most
-        // significant
-        //  bit is ask Anna again
-
-        frame.channel_mask = std::array<uint8_t, 3>{binary_channel_numbers[0],
-                                                    binary_channel_numbers[1],
-                                                    binary_channel_numbers[2]};
-
-        close_file();
-
-        return frame;
-    }
-
-  private:
-    std::filesystem::path m_base_path{};
-    std::string file_prefix{};
-};
-
-} // namespace aare
--- a/include/aare/NDArray.hpp
+++ b/include/aare/NDArray.hpp
@ -21,6 +21,7 @@ TODO! Add expression templates for operators

 namespace aare {

+
 template <typename T, ssize_t Ndim = 2>
 class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
    std::array<ssize_t, Ndim> shape_;
@ -47,6 +48,7 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
                                std::multiplies<>())),
          data_(new T[size_]) {}

+
    /**
     * @brief Construct a new NDArray object with a shape and value.
     *
@ -67,8 +69,8 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
        std::copy(v.begin(), v.end(), begin());
    }

-    template <size_t Size>
-    NDArray(const std::array<T, Size> &arr) : NDArray<T, 1>({Size}) {
+    template<size_t Size>
+    NDArray(const std::array<T, Size>& arr) : NDArray<T,1>({Size}) {
        std::copy(arr.begin(), arr.end(), begin());
    }

@ -77,6 +79,7 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
        : shape_(other.shape_), strides_(c_strides<Ndim>(shape_)),
          size_(other.size_), data_(other.data_) {
        other.reset(); // TODO! is this necessary?
+
    }

    // Copy constructor
@ -110,10 +113,10 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
    NDArray &operator-=(const NDArray &other);
    NDArray &operator*=(const NDArray &other);

-    // Write directly to the data array, or create a new one
-    template <size_t Size>
-    NDArray<T, 1> &operator=(const std::array<T, Size> &other) {
-        if (Size != size_) {
+    //Write directly to the data array, or create a new one
+    template<size_t Size>
+    NDArray<T,1>& operator=(const std::array<T,Size> &other){
+        if(Size != size_){
            delete[] data_;
            size_ = Size;
            data_ = new T[size_];
@ -139,9 +142,6 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {

    NDArray<bool, Ndim> operator>(const NDArray &other);

-    bool equals(const NDArray<T, Ndim> &other,
-                const T tolerance = std::numeric_limits<T>::epsilon()) const;
-
    bool operator==(const NDArray &other) const;
    bool operator!=(const NDArray &other) const;

@ -157,6 +157,11 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {

    NDArray &operator&=(const T & /*mask*/);

+  
+
+
+    
+
    void sqrt() {
        for (int i = 0; i < size_; ++i) {
            data_[i] = std::sqrt(data_[i]);
@ -340,6 +345,9 @@ NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const T &value) {
    return *this;
 }

+
+
+
 template <typename T, ssize_t Ndim>
 NDArray<T, Ndim> NDArray<T, Ndim>::operator+(const T &value) {
    NDArray result = *this;
@ -383,7 +391,6 @@ NDArray<T, Ndim> NDArray<T, Ndim>::operator*(const T &value) {
    result *= value;
    return result;
 }
-
 // template <typename T, ssize_t Ndim> void NDArray<T, Ndim>::Print() {
 //     if (shape_[0] < 20 && shape_[1] < 20)
 //         Print_all();
@ -426,7 +433,7 @@ template <typename T, ssize_t Ndim>
 void save(NDArray<T, Ndim> &img, std::string &pathname) {
    std::ofstream f;
    f.open(pathname, std::ios::binary);
-    f.write(reinterpret_cast<char *>(img.buffer()), img.size() * sizeof(T));
+    f.write(img.buffer(), img.size() * sizeof(T));
    f.close();
 }

@ -436,38 +443,11 @@ NDArray<T, Ndim> load(const std::string &pathname,
    NDArray<T, Ndim> img{shape};
    std::ifstream f;
    f.open(pathname, std::ios::binary);
-    f.read(reinterpret_cast<char *>(img.buffer()), img.size() * sizeof(T));
+    f.read(img.buffer(), img.size() * sizeof(T));
    f.close();
    return img;
 }

-template <typename T, ssize_t Ndim = 1>
-NDArray<T, Ndim> load_non_binary_file(const std::string &filename,
-                                      const std::array<ssize_t, Ndim> shape) {
-    std::string word;
-    NDArray<T, Ndim> array(shape);
-    try {
-        std::ifstream file(filename, std::ios_base::in);
-        if (!file.good()) {
-            throw std::logic_error("file does not exist");
-        }

-        std::stringstream file_buffer;
-        file_buffer << file.rdbuf();
-
-        ssize_t counter = 0;
-        while (file_buffer >> word && counter < size) {
-            array[counter] = static_cast<T>(
-                std::stod(word)); // TODO change for different Types
-            ++counter;
-        }
-
-        file.close();
-    } catch (const std::exception &e) {
-        std::cerr << "Error: " << e.what() << std::endl;
-    }
-
-    return array;
-}

 } // namespace aare
--- a/include/aare/NDView.hpp
+++ b/include/aare/NDView.hpp
@ -1,12 +1,11 @@
 #pragma once
-#include "aare/ArrayExpr.hpp"
 #include "aare/defs.hpp"
+#include "aare/ArrayExpr.hpp"

 #include <algorithm>
 #include <array>
 #include <cassert>
 #include <cstdint>
-#include <fstream>
 #include <functional>
 #include <iomanip>
 #include <iostream>
@ -18,8 +17,7 @@ namespace aare {
 template <ssize_t Ndim> using Shape = std::array<ssize_t, Ndim>;

 // TODO! fix mismatch between signed and unsigned
-template <ssize_t Ndim>
-Shape<Ndim> make_shape(const std::vector<size_t> &shape) {
+template <ssize_t Ndim> Shape<Ndim> make_shape(const std::vector<size_t> &shape) {
    if (shape.size() != Ndim)
        throw std::runtime_error("Shape size mismatch");
    Shape<Ndim> arr;
@ -27,18 +25,14 @@ Shape<Ndim> make_shape(const std::vector<size_t> &shape) {
    return arr;
 }

-template <ssize_t Dim = 0, typename Strides>
-ssize_t element_offset(const Strides & /*unused*/) {
-    return 0;
-}
+template <ssize_t Dim = 0, typename Strides> ssize_t element_offset(const Strides & /*unused*/) { return 0; }

 template <ssize_t Dim = 0, typename Strides, typename... Ix>
 ssize_t element_offset(const Strides &strides, ssize_t i, Ix... index) {
    return i * strides[Dim] + element_offset<Dim + 1>(strides, index...);
 }

-template <ssize_t Ndim>
-std::array<ssize_t, Ndim> c_strides(const std::array<ssize_t, Ndim> &shape) {
+template <ssize_t Ndim> std::array<ssize_t, Ndim> c_strides(const std::array<ssize_t, Ndim> &shape) {
    std::array<ssize_t, Ndim> strides{};
    std::fill(strides.begin(), strides.end(), 1);
    for (ssize_t i = Ndim - 1; i > 0; --i) {
@ -47,16 +41,14 @@ std::array<ssize_t, Ndim> c_strides(const std::array<ssize_t, Ndim> &shape) {
    return strides;
 }

-template <ssize_t Ndim>
-std::array<ssize_t, Ndim> make_array(const std::vector<ssize_t> &vec) {
+template <ssize_t Ndim> std::array<ssize_t, Ndim> make_array(const std::vector<ssize_t> &vec) {
    assert(vec.size() == Ndim);
    std::array<ssize_t, Ndim> arr{};
    std::copy_n(vec.begin(), Ndim, arr.begin());
    return arr;
 }

-template <typename T, ssize_t Ndim = 2>
-class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
+template <typename T, ssize_t Ndim = 2> class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
  public:
    NDView() = default;
    ~NDView() = default;
@ -65,23 +57,17 @@ class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {

    NDView(T *buffer, std::array<ssize_t, Ndim> shape)
        : buffer_(buffer), strides_(c_strides<Ndim>(shape)), shape_(shape),
-          size_(std::accumulate(std::begin(shape), std::end(shape), 1,
-                                std::multiplies<>())) {}
+          size_(std::accumulate(std::begin(shape), std::end(shape), 1, std::multiplies<>())) {}

    // NDView(T *buffer, const std::vector<ssize_t> &shape)
-    //     : buffer_(buffer),
-    //     strides_(c_strides<Ndim>(make_array<Ndim>(shape))),
-    //     shape_(make_array<Ndim>(shape)),
-    //       size_(std::accumulate(std::begin(shape), std::end(shape), 1,
-    //       std::multiplies<>())) {}
+    //     : buffer_(buffer), strides_(c_strides<Ndim>(make_array<Ndim>(shape))), shape_(make_array<Ndim>(shape)),
+    //       size_(std::accumulate(std::begin(shape), std::end(shape), 1, std::multiplies<>())) {}

-    template <typename... Ix>
-    std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) {
+    template <typename... Ix> std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) {
        return buffer_[element_offset(strides_, index...)];
    }

-    template <typename... Ix>
-    std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) const {
+    template <typename... Ix> std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) const {
        return buffer_[element_offset(strides_, index...)];
    }

@ -106,37 +92,18 @@ class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
        return true;
    }

-    bool equals(const NDView<T, Ndim> &other, const T tolerance) const {
-        if (shape_ != other.shape_)
-            return false;
-
-        using SignedT = typename make_signed<T>::type;
-
-        for (uint32_t i = 0; i != size_; ++i)
-            if (std::abs(static_cast<SignedT>(buffer_[i]) -
-                         static_cast<SignedT>(other.buffer_[i])) > tolerance)
-                return false;
-
-        return true;
-    }
-
    NDView &operator+=(const T val) { return elemenwise(val, std::plus<T>()); }
    NDView &operator-=(const T val) { return elemenwise(val, std::minus<T>()); }
-    NDView &operator*=(const T val) {
-        return elemenwise(val, std::multiplies<T>());
-    }
-    NDView &operator/=(const T val) {
-        return elemenwise(val, std::divides<T>());
-    }
+    NDView &operator*=(const T val) { return elemenwise(val, std::multiplies<T>()); }
+    NDView &operator/=(const T val) { return elemenwise(val, std::divides<T>()); }

-    NDView &operator/=(const NDView &other) {
-        return elemenwise(other, std::divides<T>());
-    }
+    NDView &operator/=(const NDView &other) { return elemenwise(other, std::divides<T>()); }

-    template <size_t Size> NDView &operator=(const std::array<T, Size> &arr) {
-        if (size() != static_cast<ssize_t>(arr.size()))
-            throw std::runtime_error(LOCATION +
-                                     "Array and NDView size mismatch");
+
+    template<size_t Size>
+    NDView& operator=(const std::array<T, Size> &arr) {
+        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;
    }
@ -180,15 +147,13 @@ class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
    std::array<ssize_t, Ndim> shape_{};
    uint64_t size_{};

-    template <class BinaryOperation>
-    NDView &elemenwise(T val, BinaryOperation op) {
+    template <class BinaryOperation> NDView &elemenwise(T val, BinaryOperation op) {
        for (uint64_t i = 0; i != size_; ++i) {
            buffer_[i] = op(buffer_[i], val);
        }
        return *this;
    }
-    template <class BinaryOperation>
-    NDView &elemenwise(const NDView &other, BinaryOperation op) {
+    template <class BinaryOperation> NDView &elemenwise(const NDView &other, BinaryOperation op) {
        for (uint64_t i = 0; i != size_; ++i) {
            buffer_[i] = op(buffer_[i], other.buffer_[i]);
        }
@ -205,8 +170,9 @@ template <typename T, ssize_t Ndim> void NDView<T, Ndim>::print_all() const {
    }
 }

+
 template <typename T, ssize_t Ndim>
-std::ostream &operator<<(std::ostream &os, const NDView<T, Ndim> &arr) {
+std::ostream& operator <<(std::ostream& os, const NDView<T, Ndim>& arr){
    for (auto row = 0; row < arr.shape(0); ++row) {
        for (auto col = 0; col < arr.shape(1); ++col) {
            os << std::setw(3);
@ -217,16 +183,10 @@ std::ostream &operator<<(std::ostream &os, const NDView<T, Ndim> &arr) {
    return os;
 }

-template <typename T> NDView<T, 1> make_view(std::vector<T> &vec) {
-    return NDView<T, 1>(vec.data(), {static_cast<ssize_t>(vec.size())});
-}

-template <typename T, ssize_t Ndim>
-void save(NDView<T, Ndim> img, const std::string &pathname) {
-    std::ofstream f;
-    f.open(pathname, std::ios::binary);
-    f.write(reinterpret_cast<char *>(img.data()), img.size() * sizeof(T));
-    f.close();
+template <typename T>
+NDView<T,1> make_view(std::vector<T>& vec){
+    return NDView<T,1>(vec.data(), {static_cast<ssize_t>(vec.size())});
 }

 } // namespace aare
--- a/include/aare/RawFile.hpp
+++ b/include/aare/RawFile.hpp
@ -30,22 +30,11 @@ struct ModuleConfig {
 * Consider using that unless you need raw file specific functionality.
 */
 class RawFile : public FileInterface {
-    size_t n_subfiles{}; //f0,f1...fn
-    size_t n_subfile_parts{}; // d0,d1...dn
-    //TODO! move to vector of SubFile instead of pointers
-    std::vector<std::vector<RawSubFile *>> subfiles; //subfiles[f0,f1...fn][d0,d1...dn]
-    // std::vector<xy> positions;
-    
+    std::vector<std::unique_ptr<RawSubFile>> m_subfiles;
    ModuleConfig cfg{0, 0};
-
    RawMasterFile m_master;
-
    size_t m_current_frame{};
-
-    // std::vector<ModuleGeometry> m_module_pixel_0;
-    // size_t m_rows{};
-    // size_t m_cols{};
-
+    size_t m_current_subfile{};
    DetectorGeometry m_geometry;

  public:
@ -56,7 +45,7 @@ class RawFile : public FileInterface {

     */
    RawFile(const std::filesystem::path &fname, const std::string &mode = "r");
-    virtual ~RawFile() override;
+    virtual ~RawFile() override = default;

    Frame read_frame() override;
    Frame read_frame(size_t frame_number) override;
@ -80,7 +69,7 @@ class RawFile : public FileInterface {
    size_t cols() const override;
    size_t bitdepth() const override;
    xy geometry();
-    size_t n_mod() const;
+    size_t n_modules() const;
    
    RawMasterFile master() const;

@ -115,9 +104,6 @@ class RawFile : public FileInterface {
     */
    static DetectorHeader read_header(const std::filesystem::path &fname);

-    // void update_geometry_with_roi();
-    int find_number_of_subfiles();
-
    void open_subfiles();
    void find_geometry();
 };
--- a/include/aare/RawMasterFile.hpp
+++ b/include/aare/RawMasterFile.hpp
@ -121,6 +121,7 @@ class RawMasterFile {

    size_t total_frames_expected() const;
    xy geometry() const;
+    size_t n_modules() const;

    std::optional<size_t> analog_samples() const;
    std::optional<size_t> digital_samples() const;
--- a/include/aare/RawSubFile.hpp
+++ b/include/aare/RawSubFile.hpp
@ -18,11 +18,20 @@ class RawSubFile {
    std::ifstream m_file;
    DetectorType m_detector_type;
    size_t m_bitdepth;
-    std::filesystem::path m_fname;
+    std::filesystem::path m_path;     //!< path to the subfile
+    std::string m_base_name;          //!< base name used for formatting file names
+    size_t m_offset{}; //!< file index of the first file, allow starting at non zero file
+    size_t m_total_frames{};    //!< total number of frames in the series of files
    size_t m_rows{};
    size_t m_cols{};
    size_t m_bytes_per_frame{};
-    size_t m_num_frames{};
+
+
+    int m_module_index{};
+    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
+
    uint32_t m_pos_row{};
    uint32_t m_pos_col{};
 
@ -67,12 +76,17 @@ 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; }
+    size_t frames_in_file() const { return m_total_frames; }

 private:
  template <typename T>
  void read_with_map(std::byte *image_buf);

+  void parse_fname(const std::filesystem::path &fname);
+  void scan_files();
+  void open_file(size_t file_index);
+  std::filesystem::path fpath(size_t file_index) const;
+
 };

 } // namespace aare
--- a/include/aare/algorithm.hpp
+++ b/include/aare/algorithm.hpp
@ -107,5 +107,16 @@ std::vector<T> cumsum(const std::vector<T>& vec) {
 }


+template <typename Container> bool all_equal(const Container &c) {
+    if (!c.empty() &&
+        std::all_of(begin(c), end(c),
+                    [c](const typename Container::value_type &element) {
+                        return element == c.front();
+                    }))
+        return true;
+    return false;
+}
+
+

 } // namespace aare
--- a/include/aare/defs.hpp
+++ b/include/aare/defs.hpp
@ -3,16 +3,16 @@
 #include "aare/Dtype.hpp"

 #include <array>
+#include <stdexcept>
 #include <cassert>
 #include <cstdint>
 #include <cstring>
-#include <stdexcept>
 #include <string>
 #include <string_view>
-#include <type_traits>
 #include <variant>
 #include <vector>

+
 /**
 * @brief LOCATION macro to get the current location in the code
 */
@ -20,24 +20,28 @@
    std::string(__FILE__) + std::string(":") + std::to_string(__LINE__) +      \
        ":" + std::string(__func__) + ":"

+
+
 #ifdef AARE_CUSTOM_ASSERT
-#define AARE_ASSERT(expr)                                                      \
-    if (expr) {                                                                \
-    } else                                                                     \
+#define AARE_ASSERT(expr)\
+    if (expr)\
+        {}\
+    else\
        aare::assert_failed(LOCATION + " Assertion failed: " + #expr + "\n");
 #else
-#define AARE_ASSERT(cond)                                                      \
-    do {                                                                       \
-        (void)sizeof(cond);                                                    \
-    } while (0)
+#define AARE_ASSERT(cond)\
+    do { (void)sizeof(cond); } while(0)
 #endif

+
 namespace aare {

 inline constexpr size_t bits_per_byte = 8;

 void assert_failed(const std::string &msg);

+
+
 class DynamicCluster {
  public:
    int cluster_sizeX;
@ -51,7 +55,7 @@ class DynamicCluster {

  public:
    DynamicCluster(int cluster_sizeX_, int cluster_sizeY_,
-                   Dtype dt_ = Dtype(typeid(int32_t)))
+            Dtype dt_ = Dtype(typeid(int32_t)))
        : cluster_sizeX(cluster_sizeX_), cluster_sizeY(cluster_sizeY_),
          dt(dt_) {
        m_data = new std::byte[cluster_sizeX * cluster_sizeY * dt.bytes()]{};
@ -175,24 +179,24 @@ template <typename T> struct t_xy {
 };
 using xy = t_xy<uint32_t>;

+
 /**
- * @brief Class to hold the geometry of a module. Where pixel 0 is located and
- * the size of the module
+ * @brief Class to hold the geometry of a module. Where pixel 0 is located and the size of the module
 */
-struct ModuleGeometry {
+struct ModuleGeometry{
    int origin_x{};
    int origin_y{};
    int height{};
    int width{};
    int row_index{};
-    int col_index{};
+    int col_index{}; 
 };

 /**
- * @brief Class to hold the geometry of a detector. Number of modules, their
- * size and where pixel 0 for each module is located
+ * @brief Class to hold the geometry of a detector. Number of modules, their size and where pixel 0 
+ * for each module is located
 */
-struct DetectorGeometry {
+struct DetectorGeometry{
    int modules_x{};
    int modules_y{};
    int pixels_x{};
@ -200,32 +204,35 @@ struct DetectorGeometry {
    int module_gap_row{};
    int module_gap_col{};
    std::vector<ModuleGeometry> module_pixel_0;
+    
+    auto size() const { return module_pixel_0.size(); }
 };

-struct ROI {
+struct ROI{
    ssize_t xmin{};
    ssize_t xmax{};
    ssize_t ymin{};
    ssize_t ymax{};
-
+  
    ssize_t height() const { return ymax - ymin; }
    ssize_t width() const { return xmax - xmin; }
    bool contains(ssize_t x, ssize_t y) const {
        return x >= xmin && x < xmax && y >= ymin && y < ymax;
    }
-};
+  };
+

 using dynamic_shape = std::vector<ssize_t>;

-// TODO! Can we uniform enums between the libraries?
+//TODO! Can we uniform enums between the libraries?

 /**
- * @brief Enum class to identify different detectors.
+ * @brief Enum class to identify different detectors. 
 * The values are the same as in slsDetectorPackage
 * Different spelling to avoid confusion with the slsDetectorPackage
 */
 enum class DetectorType {
-    // Standard detectors match the enum values from slsDetectorPackage
+    //Standard detectors match the enum values from slsDetectorPackage
    Generic,
    Eiger,
    Gotthard,
@ -236,9 +243,8 @@ enum class DetectorType {
    Gotthard2,
    Xilinx_ChipTestBoard,

-    // Additional detectors used for defining processing. Variants of the
-    // standard ones.
-    Moench03 = 100,
+    //Additional detectors used for defining processing. Variants of the standard ones.
+    Moench03=100,
    Moench03_old,
    Unknown
 };
@ -259,12 +265,4 @@ template <> FrameDiscardPolicy StringTo(const std::string & /*mode*/);

 using DataTypeVariants = std::variant<uint16_t, uint32_t>;

-template <typename T, bool = std::is_integral_v<T>> struct make_signed {
-    using type = T;
-};
-
-template <typename T> struct make_signed<T, true> {
-    using type = std::make_signed_t<T>;
-};
-
 } // namespace aare
--- a/include/aare/logger.hpp
+++ b/include/aare/logger.hpp
@ -0,0 +1,139 @@
+#pragma once
+/*Utility to log to console*/
+
+
+#include <iostream>
+#include <sstream>
+#include <sys/time.h>
+
+namespace aare {
+
+#define RED      "\x1b[31m"
+#define GREEN    "\x1b[32m"
+#define YELLOW   "\x1b[33m"
+#define BLUE     "\x1b[34m"
+#define MAGENTA  "\x1b[35m"
+#define CYAN     "\x1b[36m"
+#define GRAY     "\x1b[37m"
+#define DARKGRAY "\x1b[30m"
+
+#define BG_BLACK   "\x1b[48;5;232m"
+#define BG_RED     "\x1b[41m"
+#define BG_GREEN   "\x1b[42m"
+#define BG_YELLOW  "\x1b[43m"
+#define BG_BLUE    "\x1b[44m"
+#define BG_MAGENTA "\x1b[45m"
+#define BG_CYAN    "\x1b[46m"
+#define RESET      "\x1b[0m"
+#define BOLD       "\x1b[1m"
+
+
+enum TLogLevel {
+    logERROR,
+    logWARNING,
+    logINFOBLUE,
+    logINFOGREEN,
+    logINFORED,
+    logINFOCYAN,
+    logINFOMAGENTA,
+    logINFO,
+    logDEBUG,  // constructors, destructors etc. should still give too much output
+    logDEBUG1,
+    logDEBUG2,
+    logDEBUG3,
+    logDEBUG4,
+    logDEBUG5
+};
+
+// Compiler should optimize away anything below this value
+#ifndef AARE_LOG_LEVEL
+#define AARE_LOG_LEVEL "LOG LEVEL NOT SET IN CMAKE" //This is configured in the main CMakeLists.txt
+#endif
+
+#define __AT__                                                                 \
+    std::string(__FILE__) + std::string("::") + std::string(__func__) +        \
+        std::string("(): ")
+#define __SHORT_FORM_OF_FILE__                                                 \
+    (strrchr(__FILE__, '/') ? strrchr(__FILE__, '/') + 1 : __FILE__)
+#define __SHORT_AT__                                                           \
+    std::string(__SHORT_FORM_OF_FILE__) + std::string("::") +                  \
+        std::string(__func__) + std::string("(): ")
+
+class Logger {
+    std::ostringstream os;
+    TLogLevel m_level = AARE_LOG_LEVEL;
+
+  public:
+    Logger() = default;
+    explicit Logger(TLogLevel level) : m_level(level){};
+    ~Logger() {
+        // output in the destructor to allow for << syntax
+        os << RESET << '\n';
+        std::clog << os.str() << std::flush; // Single write
+    }
+
+    static TLogLevel &ReportingLevel() { // singelton eeh TODO! Do we need a runtime option?
+        static TLogLevel reportingLevel = logDEBUG5;
+        return reportingLevel;
+    }
+
+    // Danger this buffer need as many elements as TLogLevel
+    static const char *Color(TLogLevel level) noexcept {
+        static const char *const colors[] = {
+            RED BOLD, YELLOW BOLD, BLUE,  GREEN, RED,   CYAN,  MAGENTA,
+            RESET,    RESET,       RESET, RESET, RESET, RESET, RESET};
+        // out of bounds
+        if (level < 0 || level >= sizeof(colors) / sizeof(colors[0])) {
+            return RESET;
+        }
+        return colors[level];
+    }
+
+    // Danger this buffer need as many elements as TLogLevel
+    static std::string ToString(TLogLevel level) {
+        static const char *const buffer[] = {
+            "ERROR",  "WARNING", "INFO",   "INFO",  "INFO",
+            "INFO",   "INFO",    "INFO",   "DEBUG", "DEBUG1",
+            "DEBUG2", "DEBUG3",  "DEBUG4", "DEBUG5"};
+        // out of bounds
+        if (level < 0 || level >= sizeof(buffer) / sizeof(buffer[0])) {
+            return "UNKNOWN";
+        }
+        return buffer[level];
+    }
+
+    std::ostringstream &Get() {
+        os << Color(m_level) << "- " << Timestamp() << " " << ToString(m_level)
+           << ": ";
+        return os;
+    }
+
+    static std::string Timestamp() {
+        constexpr size_t buffer_len = 12;
+        char buffer[buffer_len];
+        time_t t;
+        ::time(&t);
+        tm r;
+        strftime(buffer, buffer_len, "%X", localtime_r(&t, &r));
+        buffer[buffer_len - 1] = '\0';
+        struct timeval tv;
+        gettimeofday(&tv, nullptr);
+        constexpr size_t result_len = 100;
+        char result[result_len];
+        snprintf(result, result_len, "%s.%03ld", buffer,
+                 static_cast<long>(tv.tv_usec) / 1000);
+        result[result_len - 1] = '\0';
+        return result;
+    }
+};
+
+// TODO! Do we need to keep the runtime option?
+#define LOG(level)                                                             \
+    if (level > AARE_LOG_LEVEL)                                                \
+        ;                                                                      \
+    else if (level > aare::Logger::ReportingLevel())                           \
+        ;                                                                      \
+    else                                                                       \
+        aare::Logger(level).Get()
+
+} // namespace aare
--- a/python/aare/ClusterFinder.py
+++ b/python/aare/ClusterFinder.py
@ -1,22 +1,47 @@

-from ._aare import ClusterFinder_Cluster3x3i, ClusterFinder_Cluster2x2i, ClusterFinderMT_Cluster3x3i, ClusterFinderMT_Cluster2x2i, ClusterCollector_Cluster3x3i, ClusterCollector_Cluster2x2i
+# from ._aare import ClusterFinder_Cluster3x3i, ClusterFinder_Cluster2x2i, ClusterFinderMT_Cluster3x3i, ClusterFinderMT_Cluster2x2i, ClusterCollector_Cluster3x3i, ClusterCollector_Cluster2x2i


-from ._aare import ClusterFileSink_Cluster3x3i, ClusterFileSink_Cluster2x2i
+# from ._aare import ClusterFileSink_Cluster3x3i, ClusterFileSink_Cluster2x2i
+
+from . import _aare
 import numpy as np

+_supported_cluster_sizes = [(2,2), (3,3), (5,5), (7,7), (9,9),]
+
+# def _get_class()
+
+def _type_to_char(dtype):
+    if dtype == np.int32:
+        return 'i'
+    elif dtype == np.float32:
+        return 'f'
+    elif dtype == np.float64:
+        return 'd'
+    else:
+        raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32, np.float32, and np.float64 are supported.")
+
+def _get_class(name, cluster_size, dtype):
+    """
+    Helper function to get the class based on the name, cluster size, and dtype.
+    """
+    try:
+        class_name = f"{name}_Cluster{cluster_size[0]}x{cluster_size[1]}{_type_to_char(dtype)}"
+        cls = getattr(_aare, class_name)
+    except AttributeError:
+        raise ValueError(f"Unsupported combination of type and cluster size: {dtype}/{cluster_size} when requesting {class_name}")
+    return cls
+
+
+
 def ClusterFinder(image_size, cluster_size, n_sigma=5, dtype = np.int32, capacity = 1024):
    """
    Factory function to create a ClusterFinder object. Provides a cleaner syntax for 
    the templated ClusterFinder in C++.
    """
-    if dtype == np.int32 and cluster_size == (3,3):
-        return ClusterFinder_Cluster3x3i(image_size, n_sigma = n_sigma, capacity=capacity)
-    elif dtype == np.int32 and cluster_size == (2,2): 
-        return ClusterFinder_Cluster2x2i(image_size, n_sigma = n_sigma, capacity=capacity)
-    else:
-        #TODO! add the other formats
-        raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32 is supported.")
+    cls = _get_class("ClusterFinder", cluster_size, dtype)
+    return cls(image_size, n_sigma=n_sigma, capacity=capacity)
+


 def ClusterFinderMT(image_size, cluster_size = (3,3), dtype=np.int32, n_sigma=5, capacity = 1024, n_threads = 3): 
@ -25,15 +50,9 @@ def ClusterFinderMT(image_size, cluster_size = (3,3), dtype=np.int32, n_sigma=5,
    the templated ClusterFinderMT in C++.
    """

-    if dtype == np.int32 and cluster_size == (3,3):
-        return ClusterFinderMT_Cluster3x3i(image_size, n_sigma = n_sigma,
-                    capacity = capacity, n_threads = n_threads)
-    elif dtype == np.int32 and cluster_size == (2,2):
-        return ClusterFinderMT_Cluster2x2i(image_size, n_sigma = n_sigma,
-                    capacity = capacity, n_threads = n_threads)
-    else:
-        #TODO! add the other formats
-        raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32 is supported.")
+    cls = _get_class("ClusterFinderMT", cluster_size, dtype)
+    return cls(image_size, n_sigma=n_sigma, capacity=capacity, n_threads=n_threads)
+


 def ClusterCollector(clusterfindermt, cluster_size = (3,3), dtype=np.int32): 
@ -42,14 +61,8 @@ def ClusterCollector(clusterfindermt, cluster_size = (3,3), dtype=np.int32):
    the templated ClusterCollector in C++.
    """

-    if dtype == np.int32 and cluster_size == (3,3): 
-        return ClusterCollector_Cluster3x3i(clusterfindermt)
-    elif dtype == np.int32 and cluster_size == (2,2):
-        return ClusterCollector_Cluster2x2i(clusterfindermt)
-     
-    else: 
-        #TODO! add the other formats
-        raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32 is supported.")
+    cls = _get_class("ClusterCollector", cluster_size, dtype)
+    return cls(clusterfindermt)

 def ClusterFileSink(clusterfindermt, cluster_file, dtype=np.int32): 
    """ 
@ -57,11 +70,15 @@ def ClusterFileSink(clusterfindermt, cluster_file, dtype=np.int32):
    the templated ClusterCollector in C++.
    """

-    if dtype == np.int32 and clusterfindermt.cluster_size == (3,3): 
-        return ClusterFileSink_Cluster3x3i(clusterfindermt, cluster_file)
-    elif dtype == np.int32 and clusterfindermt.cluster_size == (2,2):
-        return ClusterFileSink_Cluster2x2i(clusterfindermt, cluster_file)
-     
-    else: 
-        #TODO! add the other formats
-        raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32 is supported.")
+    cls = _get_class("ClusterFileSink", clusterfindermt.cluster_size, dtype)
+    return cls(clusterfindermt, cluster_file)
+
+
+def ClusterFile(fname, cluster_size=(3,3), dtype=np.int32):
+    """
+    Factory function to create a ClusterFile object. Provides a cleaner syntax for
+    the templated ClusterFile in C++.
+    """
+
+    cls = _get_class("ClusterFile", cluster_size, dtype)
+    return cls(fname)
--- a/python/aare/init.py
+++ b/python/aare/init.py
@ -5,13 +5,12 @@ from . import _aare
 from ._aare import File, RawMasterFile, RawSubFile, JungfrauDataFile
 from ._aare import Pedestal_d, Pedestal_f, ClusterFinder_Cluster3x3i, VarClusterFinder
 from ._aare import DetectorType
-from ._aare import ClusterFile_Cluster3x3i as ClusterFile
 from ._aare import hitmap
 from ._aare import ROI

 # from ._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink, ClusterVector_i

-from .ClusterFinder import ClusterFinder, ClusterCollector, ClusterFinderMT, ClusterFileSink
+from .ClusterFinder import ClusterFinder, ClusterCollector, ClusterFinderMT, ClusterFileSink, ClusterFile
 from .ClusterVector import ClusterVector


--- a/python/examples/play.py
+++ b/python/examples/play.py
@ -1,79 +1,89 @@
 import sys
 sys.path.append('/home/l_msdetect/erik/aare/build')

-from aare._aare import ClusterVector_i, Interpolator

-import pickle 
-import numpy as np
-import matplotlib.pyplot as plt
-import boost_histogram as bh
-import torch
-import math
-import time
+from aare import RawSubFile, DetectorType, RawFile
+
+from pathlib import Path
+path = Path("/home/l_msdetect/erik/data/aare-test-data/raw/jungfrau/")
+f = RawSubFile(path/"jungfrau_single_d0_f0_0.raw", DetectorType.Jungfrau, 512, 1024, 16)
+
+# f = RawFile(path/"jungfrau_single_master_0.json")
+
+
+# from aare._aare import ClusterVector_i, Interpolator
+
+# import pickle 
+# import numpy as np
+# import matplotlib.pyplot as plt
+# import boost_histogram as bh
+# import torch
+# import math
+# import time



-def gaussian_2d(mx, my, sigma = 1, res=100, grid_size = 2):
-    """
-    Generate a 2D gaussian as position mx, my, with sigma=sigma. 
-    The gaussian is placed on a 2x2 pixel matrix with resolution 
-    res in one dimesion.
-    """
-    x = torch.linspace(0, pixel_size*grid_size, res)
-    x,y = torch.meshgrid(x,x, indexing="ij")
-    return 1 / (2*math.pi*sigma**2) * \
-      torch.exp(-((x - my)**2 / (2*sigma**2) + (y - mx)**2 / (2*sigma**2)))
+# def gaussian_2d(mx, my, sigma = 1, res=100, grid_size = 2):
+#     """
+#     Generate a 2D gaussian as position mx, my, with sigma=sigma. 
+#     The gaussian is placed on a 2x2 pixel matrix with resolution 
+#     res in one dimesion.
+#     """
+#     x = torch.linspace(0, pixel_size*grid_size, res)
+#     x,y = torch.meshgrid(x,x, indexing="ij")
+#     return 1 / (2*math.pi*sigma**2) * \
+#       torch.exp(-((x - my)**2 / (2*sigma**2) + (y - mx)**2 / (2*sigma**2)))

-scale = 1000 #Scale factor when converting to integer
-pixel_size = 25 #um
-grid = 2
-resolution = 100
-sigma_um = 10
-xa = np.linspace(0,grid*pixel_size,resolution)
-ticks = [0, 25, 50]
+# scale = 1000 #Scale factor when converting to integer
+# pixel_size = 25 #um
+# grid = 2
+# resolution = 100
+# sigma_um = 10
+# xa = np.linspace(0,grid*pixel_size,resolution)
+# ticks = [0, 25, 50]

-hit = np.array((20,20))
-etahist_fname = "/home/l_msdetect/erik/tmp/test_hist.pkl"
+# hit = np.array((20,20))
+# etahist_fname = "/home/l_msdetect/erik/tmp/test_hist.pkl"

-local_resolution = 99
-grid_size = 3
-xaxis = np.linspace(0,grid_size*pixel_size, local_resolution)
-t = gaussian_2d(hit[0],hit[1], grid_size = grid_size, sigma = 10, res = local_resolution)
-pixels = t.reshape(grid_size, t.shape[0] // grid_size, grid_size, t.shape[1] // grid_size).sum(axis = 3).sum(axis = 1)
-pixels = pixels.numpy()
-pixels = (pixels*scale).astype(np.int32)
-v = ClusterVector_i(3,3)
-v.push_back(1,1, pixels)
+# local_resolution = 99
+# grid_size = 3
+# xaxis = np.linspace(0,grid_size*pixel_size, local_resolution)
+# t = gaussian_2d(hit[0],hit[1], grid_size = grid_size, sigma = 10, res = local_resolution)
+# pixels = t.reshape(grid_size, t.shape[0] // grid_size, grid_size, t.shape[1] // grid_size).sum(axis = 3).sum(axis = 1)
+# pixels = pixels.numpy()
+# pixels = (pixels*scale).astype(np.int32)
+# v = ClusterVector_i(3,3)
+# v.push_back(1,1, pixels)

-with open(etahist_fname, "rb") as f:
-        hist = pickle.load(f)
-eta = hist.view().copy()
-etabinsx = np.array(hist.axes.edges.T[0].flat)
-etabinsy = np.array(hist.axes.edges.T[1].flat)
-ebins = np.array(hist.axes.edges.T[2].flat)
-p = Interpolator(eta, etabinsx[0:-1], etabinsy[0:-1], ebins[0:-1])
+# with open(etahist_fname, "rb") as f:
+#         hist = pickle.load(f)
+# eta = hist.view().copy()
+# etabinsx = np.array(hist.axes.edges.T[0].flat)
+# etabinsy = np.array(hist.axes.edges.T[1].flat)
+# ebins = np.array(hist.axes.edges.T[2].flat)
+# p = Interpolator(eta, etabinsx[0:-1], etabinsy[0:-1], ebins[0:-1])




-#Generate the hit
+# #Generate the hit




-tmp = p.interpolate(v)
-print(f'tmp:{tmp}')
-pos = np.array((tmp['x'], tmp['y']))*25
+# tmp = p.interpolate(v)
+# print(f'tmp:{tmp}')
+# pos = np.array((tmp['x'], tmp['y']))*25


-print(pixels)
-fig, ax = plt.subplots(figsize = (7,7))
-ax.pcolormesh(xaxis, xaxis, t)
-ax.plot(*pos, 'o')
-ax.set_xticks([0,25,50,75])
-ax.set_yticks([0,25,50,75])
-ax.set_xlim(0,75)
-ax.set_ylim(0,75)
-ax.grid()
-print(f'{hit=}')
-print(f'{pos=}')
+# print(pixels)
+# fig, ax = plt.subplots(figsize = (7,7))
+# ax.pcolormesh(xaxis, xaxis, t)
+# ax.plot(*pos, 'o')
+# ax.set_xticks([0,25,50,75])
+# ax.set_yticks([0,25,50,75])
+# ax.set_xlim(0,75)
+# ax.set_ylim(0,75)
+# ax.grid()
+# print(f'{hit=}')
+# print(f'{pos=}')
--- a/python/src/bind_Cluster.hpp
+++ b/python/src/bind_Cluster.hpp
@ -0,0 +1,64 @@
+#include "aare/Cluster.hpp"
+
+#include <cstdint>
+#include <fmt/format.h>
+#include <filesystem>
+#include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
+#include <pybind11/stl.h>
+#include <pybind11/stl_bind.h>
+
+namespace py = pybind11;
+using pd_type = double;
+
+using namespace aare;
+
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wunused-parameter"
+
+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);
+
+    py::class_<Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>>(
+        m, class_name.c_str(), py::buffer_protocol())
+
+        .def(py::init([](uint8_t x, uint8_t y, py::array_t<Type> data) {
+            py::buffer_info buf_info = data.request();
+            Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType> cluster;
+            cluster.x = x;
+            cluster.y = y;
+            auto r = data.template unchecked<1>(); // no bounds checks
+            for (py::ssize_t i = 0; i < data.size(); ++i) {
+                cluster.data[i] = r(i);
+            }
+            return cluster;
+        }));
+
+    /*
+    //TODO! Review if to keep or not
+    .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!!!
+
+        });
+    */
+}
+
+#pragma GCC diagnostic pop
--- a/python/src/bind_ClusterCollector.hpp
+++ b/python/src/bind_ClusterCollector.hpp
@ -0,0 +1,46 @@
+#include "aare/ClusterCollector.hpp"
+#include "aare/ClusterFileSink.hpp"
+#include "aare/ClusterFinder.hpp"
+#include "aare/ClusterFinderMT.hpp"
+#include "aare/ClusterVector.hpp"
+#include "aare/NDView.hpp"
+#include "aare/Pedestal.hpp"
+#include "np_helper.hpp"
+
+#include <cstdint>
+#include <filesystem>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/stl_bind.h>
+
+namespace py = pybind11;
+using pd_type = double;
+
+using namespace aare;
+
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wunused-parameter"
+
+
+template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = uint16_t>
+void define_ClusterCollector(py::module &m,
+                                       const std::string &typestr) {
+    auto class_name = fmt::format("ClusterCollector_{}", typestr);
+
+    using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
+
+    py::class_<ClusterCollector<ClusterType>>(m, class_name.c_str())
+        .def(py::init<ClusterFinderMT<ClusterType, uint16_t, double> *>())
+        .def("stop", &ClusterCollector<ClusterType>::stop)
+        .def(
+            "steal_clusters",
+            [](ClusterCollector<ClusterType> &self) {
+                auto v = new std::vector<ClusterVector<ClusterType>>(
+                    self.steal_clusters());
+                return v; // TODO change!!!
+            },
+            py::return_value_policy::take_ownership);
+}
+
+#pragma GCC diagnostic pop
--- a/python/src/bind_ClusterFile.hpp
+++ b/python/src/bind_ClusterFile.hpp
@ -21,7 +21,7 @@ using namespace ::aare;

 template <typename Type, uint8_t CoordSizeX, uint8_t CoordSizeY,
          typename CoordType = uint16_t>
-void define_cluster_file_io_bindings(py::module &m,
+void define_ClusterFile(py::module &m,
                                     const std::string &typestr) {

    using ClusterType = Cluster<Type, CoordSizeX, CoordSizeY, CoordType>;
--- a/python/src/bind_ClusterFileSink.hpp
+++ b/python/src/bind_ClusterFileSink.hpp
@ -0,0 +1,44 @@
+#include "aare/ClusterCollector.hpp"
+#include "aare/ClusterFileSink.hpp"
+#include "aare/ClusterFinder.hpp"
+#include "aare/ClusterFinderMT.hpp"
+#include "aare/ClusterVector.hpp"
+#include "aare/NDView.hpp"
+#include "aare/Pedestal.hpp"
+#include "np_helper.hpp"
+
+#include <cstdint>
+#include <filesystem>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/stl_bind.h>
+
+namespace py = pybind11;
+using pd_type = double;
+
+using namespace aare;
+
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wunused-parameter"
+
+
+
+
+
+
+template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = uint16_t>
+void define_ClusterFileSink(py::module &m,
+                                       const std::string &typestr) {
+    auto class_name = fmt::format("ClusterFileSink_{}", typestr);
+
+    using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
+
+    py::class_<ClusterFileSink<ClusterType>>(m, class_name.c_str())
+        .def(py::init<ClusterFinderMT<ClusterType, uint16_t, double> *,
+                      const std::filesystem::path &>())
+        .def("stop", &ClusterFileSink<ClusterType>::stop);
+}
+
+
+#pragma GCC diagnostic pop
--- a/python/src/bind_ClusterFinder.hpp
+++ b/python/src/bind_ClusterFinder.hpp
@ -0,0 +1,77 @@
+#include "aare/ClusterCollector.hpp"
+#include "aare/ClusterFileSink.hpp"
+#include "aare/ClusterFinder.hpp"
+#include "aare/ClusterFinderMT.hpp"
+#include "aare/ClusterVector.hpp"
+#include "aare/NDView.hpp"
+#include "aare/Pedestal.hpp"
+#include "np_helper.hpp"
+
+#include <cstdint>
+#include <filesystem>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/stl_bind.h>
+
+namespace py = pybind11;
+using pd_type = double;
+
+using namespace aare;
+
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wunused-parameter"
+
+template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = uint16_t>
+void define_ClusterFinder(py::module &m, const std::string &typestr) {
+    auto class_name = fmt::format("ClusterFinder_{}", typestr);
+
+    using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
+
+    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"),
+             py::arg("n_sigma") = 5.0, py::arg("capacity") = 1'000'000)
+        .def("push_pedestal_frame",
+             [](ClusterFinder<ClusterType, 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(
+            "steal_clusters",
+            [](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
+               bool realloc_same_capacity) {
+                ClusterVector<ClusterType> clusters =
+                    self.steal_clusters(realloc_same_capacity);
+                return clusters;
+            },
+            py::arg("realloc_same_capacity") = false)
+        .def(
+            "find_clusters",
+            [](ClusterFinder<ClusterType, 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);
+}
+
+#pragma GCC diagnostic pop
--- a/python/src/bind_ClusterFinderMT.hpp
+++ b/python/src/bind_ClusterFinderMT.hpp
@ -0,0 +1,81 @@
+#include "aare/ClusterCollector.hpp"
+#include "aare/ClusterFileSink.hpp"
+#include "aare/ClusterFinder.hpp"
+#include "aare/ClusterFinderMT.hpp"
+#include "aare/ClusterVector.hpp"
+#include "aare/NDView.hpp"
+#include "aare/Pedestal.hpp"
+#include "np_helper.hpp"
+
+#include <cstdint>
+#include <filesystem>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/stl_bind.h>
+
+namespace py = pybind11;
+using pd_type = double;
+
+using namespace aare;
+
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wunused-parameter"
+
+template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = uint16_t>
+void define_ClusterFinderMT(py::module &m,
+                                       const std::string &typestr) {
+    auto class_name = fmt::format("ClusterFinderMT_{}", typestr);
+
+    using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
+
+    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)
+        .def("push_pedestal_frame",
+             [](ClusterFinderMT<ClusterType, 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,
+               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_property_readonly("cluster_size", [](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self){
+            return py::make_tuple(ClusterSizeX, ClusterSizeY);
+        })
+        .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);
+}
+
+
+#pragma GCC diagnostic pop
--- a/python/src/bind_ClusterVector.hpp
+++ b/python/src/bind_ClusterVector.hpp
@ -101,4 +101,6 @@ void define_ClusterVector(py::module &m, const std::string &typestr) {

              return hitmap;
          });
-}
+}
+
+#pragma GCC diagnostic pop
--- a/python/src/cluster.hpp
+++ b/python/src/cluster.hpp
@ -1,211 +0,0 @@
-#include "aare/ClusterCollector.hpp"
-#include "aare/ClusterFileSink.hpp"
-#include "aare/ClusterFinder.hpp"
-#include "aare/ClusterFinderMT.hpp"
-#include "aare/ClusterVector.hpp"
-#include "aare/NDView.hpp"
-#include "aare/Pedestal.hpp"
-#include "np_helper.hpp"
-
-#include <cstdint>
-#include <filesystem>
-#include <pybind11/pybind11.h>
-#include <pybind11/stl.h>
-#include <pybind11/stl_bind.h>
-
-namespace py = pybind11;
-using pd_type = double;
-
-using namespace aare;
-
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wunused-parameter"
-
-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);
-
-    py::class_<Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>>(
-        m, class_name.c_str(), py::buffer_protocol())
-
-        .def(py::init([](uint8_t x, uint8_t y, py::array_t<Type> data) {
-            py::buffer_info buf_info = data.request();
-            Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType> cluster;
-            cluster.x = x;
-            cluster.y = y;
-            auto r = data.template unchecked<1>(); // no bounds checks
-            for (py::ssize_t i = 0; i < data.size(); ++i) {
-                cluster.data[i] = r(i);
-            }
-            return cluster;
-        }));
-
-    /*
-    .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 T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType = uint16_t>
-void define_cluster_finder_mt_bindings(py::module &m,
-                                       const std::string &typestr) {
-    auto class_name = fmt::format("ClusterFinderMT_{}", typestr);
-
-    using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
-
-    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)
-        .def("push_pedestal_frame",
-             [](ClusterFinderMT<ClusterType, 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,
-               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_property_readonly("cluster_size", [](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self){
-            return py::make_tuple(ClusterSizeX, ClusterSizeY);
-        })
-        .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);
-}
-
-template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType = uint16_t>
-void define_cluster_collector_bindings(py::module &m,
-                                       const std::string &typestr) {
-    auto class_name = fmt::format("ClusterCollector_{}", typestr);
-
-    using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
-
-    py::class_<ClusterCollector<ClusterType>>(m, class_name.c_str())
-        .def(py::init<ClusterFinderMT<ClusterType, uint16_t, double> *>())
-        .def("stop", &ClusterCollector<ClusterType>::stop)
-        .def(
-            "steal_clusters",
-            [](ClusterCollector<ClusterType> &self) {
-                auto v = new std::vector<ClusterVector<ClusterType>>(
-                    self.steal_clusters());
-                return v; // TODO change!!!
-            },
-            py::return_value_policy::take_ownership);
-}
-
-template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType = uint16_t>
-void define_cluster_file_sink_bindings(py::module &m,
-                                       const std::string &typestr) {
-    auto class_name = fmt::format("ClusterFileSink_{}", typestr);
-
-    using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
-
-    py::class_<ClusterFileSink<ClusterType>>(m, class_name.c_str())
-        .def(py::init<ClusterFinderMT<ClusterType, uint16_t, double> *,
-                      const std::filesystem::path &>())
-        .def("stop", &ClusterFileSink<ClusterType>::stop);
-}
-
-template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType = uint16_t>
-void define_cluster_finder_bindings(py::module &m, const std::string &typestr) {
-    auto class_name = fmt::format("ClusterFinder_{}", typestr);
-
-    using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
-
-    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"),
-             py::arg("n_sigma") = 5.0, py::arg("capacity") = 1'000'000)
-        .def("push_pedestal_frame",
-             [](ClusterFinder<ClusterType, 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(
-            "steal_clusters",
-            [](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
-               bool realloc_same_capacity) {
-                ClusterVector<ClusterType> clusters =
-                    self.steal_clusters(realloc_same_capacity);
-                return clusters;
-            },
-            py::arg("realloc_same_capacity") = false)
-        .def(
-            "find_clusters",
-            [](ClusterFinder<ClusterType, 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);
-}
-#pragma GCC diagnostic pop
--- a/python/src/module.cpp
+++ b/python/src/module.cpp
@ -1,11 +1,15 @@
 // Files with bindings to the different classes

 //New style file naming
+#include "bind_Cluster.hpp"
+#include "bind_ClusterCollector.hpp"
+#include "bind_ClusterFinder.hpp"
+#include "bind_ClusterFinderMT.hpp"
+#include "bind_ClusterFile.hpp"
+#include "bind_ClusterFileSink.hpp"
 #include "bind_ClusterVector.hpp"

 //TODO! migrate the other names
-#include "cluster.hpp"
-#include "cluster_file.hpp"
 #include "ctb_raw_file.hpp"
 #include "file.hpp"
 #include "fit.hpp"
@ -24,6 +28,25 @@

 namespace py = pybind11;

+/* MACRO that defines Cluster bindings for a specific size and type
+
+T - Storage type of the cluster data (int, float, double)
+N - Number of rows in the cluster
+M - Number of columns in the cluster
+U - Type of the pixel data (e.g., uint16_t)
+TYPE_CODE - A character representing the type code (e.g., 'i' for int, 'd' for double, 'f' for float)
+
+*/
+#define DEFINE_CLUSTER_BINDINGS(T, N, M, U, TYPE_CODE) \
+    define_ClusterFile<T, N, M, U>(m, "Cluster" #N "x" #M #TYPE_CODE); \
+    define_ClusterVector<T, N, M, U>(m, "Cluster" #N "x" #M #TYPE_CODE); \
+    define_ClusterFinder<T, N, M, U>(m, "Cluster" #N "x" #M #TYPE_CODE); \
+    define_ClusterFinderMT<T, N, M, U>(m, "Cluster" #N "x" #M #TYPE_CODE); \
+    define_ClusterFileSink<T, N, M, U>(m, "Cluster" #N "x" #M #TYPE_CODE); \
+    define_ClusterCollector<T, N, M, U>(m, "Cluster" #N "x" #M #TYPE_CODE); \
+    define_Cluster<T, N, M, U>(m, #N "x" #M #TYPE_CODE); \
+    register_calculate_eta<T, N, M, U>(m);
+
 PYBIND11_MODULE(_aare, m) {
    define_file_io_bindings(m);
    define_raw_file_io_bindings(m);
@ -38,59 +61,23 @@ PYBIND11_MODULE(_aare, m) {
    define_interpolation_bindings(m);
    define_jungfrau_data_file_io_bindings(m);

-    define_cluster_file_io_bindings<int, 3, 3, uint16_t>(m, "Cluster3x3i");
-    define_cluster_file_io_bindings<double, 3, 3, uint16_t>(m, "Cluster3x3d");
-    define_cluster_file_io_bindings<float, 3, 3, uint16_t>(m, "Cluster3x3f");
-    define_cluster_file_io_bindings<int, 2, 2, uint16_t>(m, "Cluster2x2i");
-    define_cluster_file_io_bindings<float, 2, 2, uint16_t>(m, "Cluster2x2f");
-    define_cluster_file_io_bindings<double, 2, 2, uint16_t>(m, "Cluster2x2d");
+    DEFINE_CLUSTER_BINDINGS(int, 3, 3, uint16_t, i);
+    DEFINE_CLUSTER_BINDINGS(double, 3, 3, uint16_t, d);
+    DEFINE_CLUSTER_BINDINGS(float, 3, 3, uint16_t, f);

-    define_ClusterVector<int, 3, 3, uint16_t>(m, "Cluster3x3i");
-    define_ClusterVector<double, 3, 3, uint16_t>(m, "Cluster3x3d");
-    define_ClusterVector<float, 3, 3, uint16_t>(m, "Cluster3x3f");
-    define_ClusterVector<int, 2, 2, uint16_t>(m, "Cluster2x2i");
-    define_ClusterVector<double, 2, 2, uint16_t>(m, "Cluster2x2d");
-    define_ClusterVector<float, 2, 2, uint16_t>(m, "Cluster2x2f");
+    DEFINE_CLUSTER_BINDINGS(int, 2, 2, uint16_t, i);
+    DEFINE_CLUSTER_BINDINGS(double, 2, 2, uint16_t, d);
+    DEFINE_CLUSTER_BINDINGS(float, 2, 2, uint16_t, f);

-    define_cluster_finder_bindings<int, 3, 3, uint16_t>(m, "Cluster3x3i");
-    define_cluster_finder_bindings<double, 3, 3, uint16_t>(m, "Cluster3x3d");
-    define_cluster_finder_bindings<float, 3, 3, uint16_t>(m, "Cluster3x3f");
-    define_cluster_finder_bindings<int, 2, 2, uint16_t>(m, "Cluster2x2i");
-    define_cluster_finder_bindings<double, 2, 2, uint16_t>(m, "Cluster2x2d");
-    define_cluster_finder_bindings<float, 2, 2, uint16_t>(m, "Cluster2x2f");
+    DEFINE_CLUSTER_BINDINGS(int, 5, 5, uint16_t, i);
+    DEFINE_CLUSTER_BINDINGS(double, 5, 5, uint16_t, d);
+    DEFINE_CLUSTER_BINDINGS(float, 5, 5, uint16_t, f);

-    define_cluster_finder_mt_bindings<int, 3, 3, uint16_t>(m, "Cluster3x3i");
-    define_cluster_finder_mt_bindings<double, 3, 3, uint16_t>(m, "Cluster3x3d");
-    define_cluster_finder_mt_bindings<float, 3, 3, uint16_t>(m, "Cluster3x3f");
-    define_cluster_finder_mt_bindings<int, 2, 2, uint16_t>(m, "Cluster2x2i");
-    define_cluster_finder_mt_bindings<double, 2, 2, uint16_t>(m, "Cluster2x2d");
-    define_cluster_finder_mt_bindings<float, 2, 2, uint16_t>(m, "Cluster2x2f");
+    DEFINE_CLUSTER_BINDINGS(int, 7, 7, uint16_t, i);
+    DEFINE_CLUSTER_BINDINGS(double, 7, 7, uint16_t, d);
+    DEFINE_CLUSTER_BINDINGS(float, 7, 7, uint16_t, f);

-    define_cluster_file_sink_bindings<int, 3, 3, uint16_t>(m, "Cluster3x3i");
-    define_cluster_file_sink_bindings<double, 3, 3, uint16_t>(m, "Cluster3x3d");
-    define_cluster_file_sink_bindings<float, 3, 3, uint16_t>(m, "Cluster3x3f");
-    define_cluster_file_sink_bindings<int, 2, 2, uint16_t>(m, "Cluster2x2i");
-    define_cluster_file_sink_bindings<double, 2, 2, uint16_t>(m, "Cluster2x2d");
-    define_cluster_file_sink_bindings<float, 2, 2, uint16_t>(m, "Cluster2x2f");
-
-    define_cluster_collector_bindings<int, 3, 3, uint16_t>(m, "Cluster3x3i");
-    define_cluster_collector_bindings<double, 3, 3, uint16_t>(m, "Cluster3x3f");
-    define_cluster_collector_bindings<float, 3, 3, uint16_t>(m, "Cluster3x3d");
-    define_cluster_collector_bindings<int, 2, 2, uint16_t>(m, "Cluster2x2i");
-    define_cluster_collector_bindings<double, 2, 2, uint16_t>(m, "Cluster2x2f");
-    define_cluster_collector_bindings<float, 2, 2, uint16_t>(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");
-
-    register_calculate_eta<int, 3, 3, uint16_t>(m);
-    register_calculate_eta<float, 3, 3, uint16_t>(m);
-    register_calculate_eta<double, 3, 3, uint16_t>(m);
-    register_calculate_eta<int, 2, 2, uint16_t>(m);
-    register_calculate_eta<float, 2, 2, uint16_t>(m);
-    register_calculate_eta<double, 2, 2, uint16_t>(m);
+    DEFINE_CLUSTER_BINDINGS(int, 9, 9, uint16_t, i);
+    DEFINE_CLUSTER_BINDINGS(double, 9, 9, uint16_t, d);
+    DEFINE_CLUSTER_BINDINGS(float, 9, 9, uint16_t, f);
 }
--- a/python/src/raw_file.hpp
+++ b/python/src/raw_file.hpp
@ -32,7 +32,7 @@ void define_raw_file_io_bindings(py::module &m) {
                 shape.push_back(self.cols());

                 // return headers from all subfiles
-                 py::array_t<DetectorHeader> header(self.n_mod());
+                 py::array_t<DetectorHeader> header(self.n_modules());

                 const uint8_t item_size = self.bytes_per_pixel();
                 if (item_size == 1) {
@ -61,10 +61,10 @@ void define_raw_file_io_bindings(py::module &m) {
            
                // return headers from all subfiles
                py::array_t<DetectorHeader> header;
-                if (self.n_mod() == 1) {
+                if (self.n_modules() == 1) {
                    header = py::array_t<DetectorHeader>(n_frames);
                } else {
-                    header = py::array_t<DetectorHeader>({self.n_mod(), n_frames});
+                    header = py::array_t<DetectorHeader>({self.n_modules(), n_frames});
                }
                // py::array_t<DetectorHeader> header({self.n_mod(), n_frames});

@ -100,7 +100,7 @@ void define_raw_file_io_bindings(py::module &m) {
        .def_property_readonly("cols", &RawFile::cols)
        .def_property_readonly("bitdepth", &RawFile::bitdepth)
        .def_property_readonly("geometry", &RawFile::geometry)
-        .def_property_readonly("n_mod", &RawFile::n_mod)
+        .def_property_readonly("n_modules", &RawFile::n_modules)
        .def_property_readonly("detector_type", &RawFile::detector_type)
        .def_property_readonly("master", &RawFile::master);
 }
--- a/python/tests/test_RawSubFile.py
+++ b/python/tests/test_RawSubFile.py
@ -5,32 +5,35 @@ from aare import RawSubFile, DetectorType

@pytest.mark.files
 def test_read_a_jungfrau_RawSubFile(test_data_path):
+
+    # Starting with f1 there is now 7 frames left in the series of files
    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
+        assert f.frames_in_file == 7

        headers, frames = f.read()

-    assert headers.size == 3
-    assert frames.shape == (3, 512, 1024)
+    assert headers.size == 7
+    assert frames.shape == (7, 512, 1024)
    
-    # Frame numbers in this file should be 4, 5, 6
-    for i,h in zip(range(4,7,1), headers):
+
+    for i,h in zip(range(4,11,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)
+    assert np.all(data[3:] == 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")

+    # Given the first subfile in a series we can read all frames from f0, f1, f2...fN
    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
+        assert i == 10
+        assert header["frameNumber"] == 10
--- a/src/AngleCalibration.cpp
+++ b/src/AngleCalibration.cpp
@ -1,377 +0,0 @@
-#include "aare/AngleCalibration.hpp"
-
-namespace aare {
-
-AngleCalibration::AngleCalibration(
-    std::shared_ptr<MythenDetectorSpecifications> mythen_detector_,
-    std::shared_ptr<FlatField> flat_field_,
-    std::shared_ptr<MythenFileReader> mythen_file_reader_)
-    : mythen_detector(mythen_detector_), flat_field(flat_field_),
-      mythen_file_reader(mythen_file_reader_) {
-    centers.reserve(mythen_detector->max_modules());
-    conversions.reserve(mythen_detector->max_modules());
-    offsets.reserve(mythen_detector->max_modules());
-
-    num_bins = std::floor(mythen_detector->max_angle() / histogram_bin_width) -
-               std::floor(mythen_detector->min_angle() / histogram_bin_width) +
-               1; // TODO only works if negative
-                  // and positive angle
-}
-
-void AngleCalibration::set_histogram_bin_width(double bin_width) {
-    histogram_bin_width = bin_width;
-
-    num_bins = std::floor(mythen_detector->max_angle() / histogram_bin_width) -
-               std::floor(mythen_detector->min_angle() / histogram_bin_width) +
-               1; // TODO only works if negative
-                  // and positive angle
-}
-
-double AngleCalibration::get_histogram_bin_width() const {
-    return histogram_bin_width;
-}
-
-ssize_t AngleCalibration::get_new_num_bins() const { return num_bins; }
-
-std::vector<double> AngleCalibration::get_centers() const { return centers; }
-
-std::vector<double> AngleCalibration::get_conversions() const {
-    return conversions;
-}
-
-std::vector<double> AngleCalibration::get_offsets() const { return offsets; }
-
-NDView<double, 1> AngleCalibration::get_new_photon_counts() const {
-    return new_photon_counts.view();
-}
-
-NDView<double, 1> AngleCalibration::get_new_statistical_errors() const {
-    return new_photon_count_errors.view();
-}
-
-void AngleCalibration::read_initial_calibration_from_file(
-    const std::string &filename) {
-
-    std::string line;
-    uint32_t module_number{};
-
-    try {
-        std::ifstream file(filename, std::ios_base::in);
-        if (!file.good()) {
-            throw std::logic_error("file does not exist");
-        }
-
-        std::stringstream file_buffer;
-        file_buffer << file.rdbuf();
-
-        while (file_buffer >> line) {
-            if (line == "module") {
-                file_buffer >> line;
-                module_number = std::stoi(line);
-            }
-            if (line == "center") {
-                file_buffer >> line;
-                centers.insert(centers.begin() + module_number,
-                               std::stod(line));
-            }
-            if (line == "conversion") {
-                file_buffer >> line;
-                conversions.insert(conversions.begin() + module_number,
-                                   std::stod(line));
-            }
-            if (line == "offset") {
-                file_buffer >> line;
-                offsets.insert(offsets.begin() + module_number,
-                               std::stod(line));
-            }
-        }
-
-        file.close();
-    } catch (const std::exception &e) {
-        std::cerr << "Error: " << e.what() << std::endl;
-    }
-}
-
-parameters AngleCalibration::convert_to_EE_parameters() const {
-
-    // normal distance between sample and detector (R)
-    std::vector<double> normal_distances(centers.size());
-    // distances between intersection point of sample normal and module origin
-    // (D)
-    std::vector<double> module_center_distances(centers.size());
-    // angles between undiffracted beam and orthogonal sample projection on
-    // detector (phi)
-    std::vector<double> angles(centers.size());
-
-    for (size_t i = 0; i < centers.size(); ++i) {
-        auto [module_center_distance, normal_distance, angle] =
-            convert_to_EE_parameters(i);
-        normal_distances[i] = normal_distance;
-        module_center_distances[i] = module_center_distance;
-        angles[i] = angle;
-    }
-
-    return std::make_tuple(module_center_distances, normal_distances, angles);
-}
-
-std::tuple<double, double, double>
-AngleCalibration::convert_to_EE_parameters(const size_t module_index) const {
-    return convert_to_EE_parameters(centers[module_index],
-                                    conversions[module_index],
-                                    offsets[module_index]);
-}
-
-std::tuple<double, double, double> AngleCalibration::convert_to_EE_parameters(
-    const double center, const double conversion, const double offset) const {
-    const double module_center_distance =
-        center * MythenDetectorSpecifications::pitch();
-    const double normal_distance =
-        MythenDetectorSpecifications::pitch() / std::abs(conversion);
-    const double angle = offset + 180.0 / M_PI * center * std::abs(conversion);
-
-    return std::make_tuple(module_center_distance, normal_distance, angle);
-}
-
-size_t AngleCalibration::global_to_local_strip_index_conversion(
-    const size_t global_strip_index) const {
-    const size_t module_index =
-        global_strip_index / MythenDetectorSpecifications::strips_per_module();
-    // local strip index in module
-    size_t local_strip_index =
-        global_strip_index -
-        module_index * MythenDetectorSpecifications::strips_per_module();
-    // switch if indexing is in clock-wise direction
-    local_strip_index =
-        std::signbit(conversions[module_index])
-            ? MythenDetectorSpecifications::strips_per_module() - 1 -
-                  local_strip_index
-            : local_strip_index;
-
-    return local_strip_index;
-}
-
-/*
-parameters
-AngleCalibration::convert_to_BC_parameters() {}
-*/
-
-double AngleCalibration::diffraction_angle_from_DG_parameters(
-    const double center, const double conversion, const double offset,
-    const size_t strip_index, const double distance_to_strip) const {
-
-    return offset + 180.0 / M_PI *
-                        (center * std::abs(conversion) -
-                         atan((center - (strip_index + distance_to_strip)) *
-                              std::abs(conversion)));
-}
-
-double AngleCalibration::diffraction_angle_from_EE_parameters(
-    const double module_center_distance, const double normal_distance,
-    const double angle, const size_t strip_index,
-    const double distance_to_strip) const {
-
-    return angle - 180.0 / M_PI *
-                       atan((module_center_distance -
-                             MythenDetectorSpecifications::pitch() *
-                                 (strip_index + distance_to_strip)) /
-                            normal_distance); // TODO: why is it minus
-                                              // is it defined counter
-                                              // clockwise? thought
-                                              // should have a flipped
-                                              // sign
-}
-
-double AngleCalibration::angular_strip_width_from_DG_parameters(
-    const double center, const double conversion, const double offset,
-    const size_t local_strip_index) const {
-
-    return std::abs(diffraction_angle_from_DG_parameters(
-                        center, conversion, offset, local_strip_index, -0.5) -
-                    diffraction_angle_from_DG_parameters(
-                        center, conversion, offset, local_strip_index, 0.5));
-}
-
-double AngleCalibration::angular_strip_width_from_EE_parameters(
-    const double module_center_distance, const double normal_distance,
-    const double angle, const size_t local_strip_index) const {
-
-    return std::abs(diffraction_angle_from_EE_parameters(
-                        module_center_distance, normal_distance, angle,
-                        local_strip_index, -0.5) -
-                    diffraction_angle_from_EE_parameters(
-                        module_center_distance, normal_distance, angle,
-                        local_strip_index, 0.5));
-
-    // TODO: again not sure about division order - taking abs anyway
-}
-
-void AngleCalibration::calculate_fixed_bin_angle_width_histogram(
-    const size_t start_frame_index, const size_t end_frame_index) {
-
-    new_photon_counts = NDArray<double, 1>(std::array<ssize_t, 1>{num_bins});
-
-    new_photon_count_errors =
-        NDArray<double, 1>(std::array<ssize_t, 1>{num_bins});
-
-    // TODO: maybe group these into a 2d array - better cache usage
-    NDArray<double, 1> bin_counts(std::array<ssize_t, 1>{num_bins}, 0.0);
-    NDArray<double, 1> new_statistical_weights(std::array<ssize_t, 1>{num_bins},
-                                               0.0);
-    NDArray<double, 1> new_errors(std::array<ssize_t, 1>{num_bins}, 0.0);
-
-    NDArray<double, 1> inverse_normalized_flatfield =
-        flat_field->inverse_normalized_flatfield();
-
-    for (size_t frame_index = start_frame_index; frame_index < end_frame_index;
-         ++frame_index) {
-        MythenFrame frame = mythen_file_reader->read_frame(frame_index);
-        redistribute_photon_counts_to_fixed_angle_bins(
-            frame, bin_counts.view(), new_statistical_weights.view(),
-            new_errors.view(), inverse_normalized_flatfield.view());
-    }
-
-    for (ssize_t i = 0; i < new_photon_counts.size(); ++i) {
-        new_photon_counts[i] = (new_statistical_weights[i] <=
-                                std::numeric_limits<double>::epsilon())
-                                   ? 0.
-                                   : bin_counts[i] / new_statistical_weights[i];
-        new_photon_count_errors[i] =
-            (bin_counts[i] <= std::numeric_limits<double>::epsilon())
-                ? 0.
-                : 1.0 / std::sqrt(bin_counts[i]);
-    }
-}
-
-void AngleCalibration::redistribute_photon_counts_to_fixed_angle_bins(
-    const MythenFrame &frame, NDView<double, 1> bin_counts,
-    NDView<double, 1> new_statistical_weights, NDView<double, 1> new_errors,
-    NDView<double, 1> inverse_normalized_flatfield) const {
-
-    ssize_t channel = 0; // TODO handle mask - FlatField still 1d
-
-    if (frame.photon_counts.shape()[0] != mythen_detector->num_strips()) {
-        throw std::runtime_error("wrong number of strips read");
-    }
-
-    ssize_t num_bins1 = mythen_detector->min_angle() / histogram_bin_width;
-    ssize_t num_bins2 = mythen_detector->max_angle() / histogram_bin_width;
-
-    std::cout << "position: " << frame.detector_angle
-              << std::endl; // replace with log
-
-    double exposure_rate = 1. / mythen_detector->exposure_time();
-
-    for (ssize_t strip_index = 0; strip_index < mythen_detector->num_strips();
-         ++strip_index) {
-
-        size_t module_index =
-            strip_index / MythenDetectorSpecifications::strips_per_module();
-
-        if (mythen_detector->get_bad_channels()[strip_index] ||
-            !mythen_detector->get_connected_modules()[module_index])
-            continue;
-
-        double poisson_error = std::sqrt(frame.photon_counts(strip_index)) *
-                               inverse_normalized_flatfield(strip_index) *
-                               exposure_rate;
-
-        double corrected_photon_count =
-            frame.photon_counts(strip_index) *
-            inverse_normalized_flatfield(strip_index) * exposure_rate;
-
-        size_t local_strip_index =
-            global_to_local_strip_index_conversion(strip_index);
-
-        double diffraction_angle = diffraction_angle_from_DG_parameters(
-            centers[module_index], conversions[module_index],
-            offsets[module_index], local_strip_index);
-
-        diffraction_angle += (frame.detector_angle + mythen_detector->dtt0() +
-                              mythen_detector->bloffset());
-
-        if (diffraction_angle < mythen_detector->min_angle() ||
-            diffraction_angle > mythen_detector->max_angle())
-            continue;
-
-        double angle_covered_by_strip = angular_strip_width_from_DG_parameters(
-            centers[module_index], conversions[module_index],
-            offsets[module_index], local_strip_index);
-
-        double photon_count_per_bin = histogram_bin_width *
-                                      corrected_photon_count /
-                                      angle_covered_by_strip;
-        double error_photon_count_per_bin =
-            histogram_bin_width * poisson_error / angle_covered_by_strip;
-
-        double statistical_weights =
-            1.0 / std::pow(error_photon_count_per_bin, 2); // 1./sigma²
-
-        double strip_boundary_left =
-            diffraction_angle - 0.5 * angle_covered_by_strip;
-        double strip_boundary_right =
-            diffraction_angle + 0.5 * angle_covered_by_strip;
-
-        ssize_t left_bin_index = std::max(
-            num_bins1,
-            static_cast<ssize_t>(
-                std::floor(strip_boundary_left / histogram_bin_width) - 1));
-        ssize_t right_bin_index = std::min(
-            num_bins2,
-            static_cast<ssize_t>(
-                std::ceil(strip_boundary_right / histogram_bin_width) + 1));
-
-        // TODO should it be < or <=
-        for (ssize_t bin = left_bin_index; bin <= right_bin_index; ++bin) {
-            double bin_coverage = std::min(strip_boundary_right,
-                                           (bin + 0.5) * histogram_bin_width) -
-                                  std::max(strip_boundary_left,
-                                           (bin - 0.5) * histogram_bin_width);
-
-            double bin_coverage_factor = bin_coverage / histogram_bin_width;
-
-            ssize_t bin_index = bin - num_bins1;
-            // TODO: maybe have this threshold configurable
-            if (bin_coverage >= 0.0001) {
-                new_statistical_weights(bin_index) +=
-                    statistical_weights * bin_coverage_factor;
-                bin_counts(bin_index) += statistical_weights *
-                                         bin_coverage_factor *
-                                         photon_count_per_bin;
-                new_errors(bin_index) += statistical_weights *
-                                         bin_coverage_factor *
-                                         std::pow(photon_count_per_bin, 2);
-            }
-        }
-    }
-}
-
-void AngleCalibration::write_to_file(
-    const std::string &filename, const bool store_nonzero_bins,
-    const std::filesystem::path &filepath) const {
-
-    std::ofstream output_file(filepath / filename);
-
-    if (!output_file) {
-        std::cerr << "Error opening file!"
-                  << std::endl; // TODO: replace with log
-    }
-
-    output_file << std::fixed << std::setprecision(15);
-
-    for (ssize_t i = 0; i < num_bins; ++i) {
-        if (new_photon_counts[i] <= std::numeric_limits<double>::epsilon() &&
-            store_nonzero_bins) {
-            continue;
-        }
-
-        output_file << std::floor(mythen_detector->min_angle() /
-                                  histogram_bin_width) *
-                               histogram_bin_width +
-                           i * histogram_bin_width
-                    << " " << new_photon_counts[i] << " "
-                    << new_photon_count_errors[i] << std::endl;
-    }
-    output_file.close();
-}
-
-} // namespace aare
--- a/src/AngleCalibration.test.cpp
+++ b/src/AngleCalibration.test.cpp
@ -1,234 +0,0 @@
-/************************************************
- * @file AngleCalibration.test.cpp
- * @short test case for angle calibration class
- ***********************************************/
-
-#include "aare/AngleCalibration.hpp"
-
-#include <filesystem>
-
-#include "test_config.hpp"
-
-#include <iomanip>
-#include <type_traits>
-
-#include <catch2/catch_all.hpp>
-#include <catch2/catch_test_macros.hpp>
-#include <catch2/matchers/catch_matchers_floating_point.hpp>
-
-using namespace aare;
-
-TEST_CASE("read initial angle calibration file",
-          "[.anglecalibration] [.files]") {
-
-    std::shared_ptr<MythenDetectorSpecifications> mythen_detector_ptr =
-        std::make_shared<MythenDetectorSpecifications>();
-
-    AngleCalibration anglecalibration(mythen_detector_ptr,
-                                      std::shared_ptr<FlatField>{},
-                                      std::shared_ptr<MythenFileReader>{});
-
-    std::string filename = test_data_path() / "AngleCalibration_Test_Data" /
-                           "Angcal_2E_Feb2023_P29.off";
-
-    REQUIRE(std::filesystem::exists(filename));
-
-    anglecalibration.read_initial_calibration_from_file(filename);
-
-    auto centers = anglecalibration.get_centers();
-    auto conversions = anglecalibration.get_conversions();
-    auto offsets = anglecalibration.get_offsets();
-
-    std::cout.precision(17);
-
-    CHECK(centers.size() == 48);
-    CHECK(conversions.size() == 48);
-    CHECK(offsets.size() == 48);
-
-    CHECK(centers[9] == Catch::Approx(660.342326));
-    CHECK(offsets[47] == Catch::Approx(5.8053312));
-    CHECK(conversions[27] == Catch::Approx(-0.6581179125e-4));
-}
-
-TEST_CASE("read bad channels",
-          "[.anglecalibration][.mythenspecifications][.files]") {
-
-    MythenDetectorSpecifications mythen_detector;
-
-    std::string bad_channels_filename = test_data_path() /
-                                        "AngleCalibration_Test_Data" /
-                                        "bc2023_003_RING.chans";
-
-    REQUIRE(std::filesystem::exists(bad_channels_filename));
-
-    mythen_detector.read_bad_channels_from_file(bad_channels_filename);
-
-    CHECK(mythen_detector.get_bad_channels().size() == 61440);
-
-    CHECK(mythen_detector.get_bad_channels()[61437] == true);
-    CHECK(std::all_of(mythen_detector.get_bad_channels().begin() + 30720,
-                      mythen_detector.get_bad_channels().begin() + 61439,
-                      [](const bool element) { return element; }));
-}
-
-TEST_CASE("read unconnected modules",
-          "[.anglecalibration][.mythenspecifications][.files]") {
-
-    MythenDetectorSpecifications mythen_detector;
-
-    std::string unconnected_modules_filename =
-        test_data_path() / "AngleCalibration_Test_Data" / "ModOut.txt";
-
-    REQUIRE(std::filesystem::exists(unconnected_modules_filename));
-
-    mythen_detector.read_unconnected_modules_from_file(
-        unconnected_modules_filename);
-
-    CHECK(mythen_detector.get_connected_modules().size() == 48);
-
-    CHECK(std::all_of(mythen_detector.get_connected_modules().begin(),
-                      mythen_detector.get_connected_modules().end(),
-                      [](const bool element) { return element; }));
-}
-
-TEST_CASE("read flatfield", "[.anglecalibration][.flatfield][.files]") {
-
-    std::shared_ptr<MythenDetectorSpecifications> mythen_detector_ptr =
-        std::make_shared<MythenDetectorSpecifications>();
-
-    FlatField flatfield(mythen_detector_ptr);
-
-    std::string flatfield_filename =
-        test_data_path() / "AngleCalibration_Test_Data" /
-        "Flatfield_E22p0keV_T11000eV_up_48M_a_LONG_Feb2023_open_WS_SUMC.raw";
-
-    REQUIRE(std::filesystem::exists(flatfield_filename));
-
-    flatfield.read_flatfield_from_file(flatfield_filename);
-
-    auto flatfield_data = flatfield.get_flatfield();
-
-    CHECK(flatfield_data.size() == 61440);
-
-    CHECK(flatfield_data[0] == 0);
-    CHECK(flatfield_data[21] == 4234186);
-}
-
-TEST_CASE("compare result with python code", "[.anglecalibration] [.files]") {
-
-    auto fpath = test_data_path() / "AngleCalibration_Test_Data";
-
-    REQUIRE(std::filesystem::exists(fpath));
-
-    std::shared_ptr<MythenDetectorSpecifications> mythen_detector_ptr =
-        std::make_shared<MythenDetectorSpecifications>();
-
-    std::string bad_channels_filename = fpath / "bc2023_003_RING.chans";
-
-    REQUIRE(std::filesystem::exists(bad_channels_filename));
-
-    mythen_detector_ptr->read_bad_channels_from_file(bad_channels_filename);
-
-    std::string unconnected_modules_filename = fpath / "ModOut.txt";
-
-    REQUIRE(std::filesystem::exists(unconnected_modules_filename));
-
-    mythen_detector_ptr->read_unconnected_modules_from_file(
-        unconnected_modules_filename);
-
-    std::shared_ptr<FlatField> flat_field_ptr =
-        std::make_shared<FlatField>(mythen_detector_ptr);
-
-    std::string flatfield_filename =
-        fpath /
-        "Flatfield_E22p0keV_T11000eV_up_48M_a_LONG_Feb2023_open_WS_SUMC.raw";
-
-    REQUIRE(std::filesystem::exists(flatfield_filename));
-
-    flat_field_ptr->read_flatfield_from_file(flatfield_filename);
-
-    std::shared_ptr<MythenFileReader> mythen_file_reader_ptr =
-        std::make_shared<MythenFileReader>(fpath,
-                                           "ang1up_22keV_LaB60p3mm_48M_a_0");
-
-    AngleCalibration anglecalibration(mythen_detector_ptr, flat_field_ptr,
-                                      mythen_file_reader_ptr);
-
-    std::string initial_angles_filename = fpath / "Angcal_2E_Feb2023_P29.off";
-
-    REQUIRE(std::filesystem::exists(initial_angles_filename));
-
-    anglecalibration.read_initial_calibration_from_file(
-        initial_angles_filename);
-
-    anglecalibration.calculate_fixed_bin_angle_width_histogram(320, 340);
-
-    // anglecalibration.write_to_file("cpp_new_photon_counts.xye");
-
-    auto expected_filename_photons =
-        test_data_path() / "AngleCalibration_Test_Data" / "new_photons.bin";
-
-    REQUIRE(std::filesystem::exists(expected_filename_photons));
-
-    auto expected_filename_errors =
-        test_data_path() / "AngleCalibration_Test_Data" / "new_errors.bin";
-
-    REQUIRE(std::filesystem::exists(expected_filename_errors));
-
-    ssize_t new_num_bins = anglecalibration.get_new_num_bins();
-
-    auto python_output_errors = load<double, 1>(
-        expected_filename_errors, std::array<ssize_t, 1>{new_num_bins});
-
-    auto python_output_photons = load<double, 1>(
-        expected_filename_photons, std::array<ssize_t, 1>{new_num_bins});
-
-    CHECK(anglecalibration.get_new_photon_counts().equals(
-        python_output_photons.view(),
-        1e-8)); // not sure about precision does not exactly match to all
-                // decimal digits
-
-    CHECK(anglecalibration.get_new_statistical_errors().equals(
-        python_output_errors.view(),
-        1e-8)); //
-}
-
-TEST_CASE("check conversion from DG to EE parameters", "[.anglecalibration]") {
-
-    std::shared_ptr<MythenDetectorSpecifications> mythen_detector_ptr =
-        std::make_shared<MythenDetectorSpecifications>();
-
-    AngleCalibration anglecalibration(mythen_detector_ptr,
-                                      std::shared_ptr<FlatField>{},
-                                      std::shared_ptr<MythenFileReader>{});
-
-    // DG test parameters
-    const double center = 642.197591224993;
-    const double conversion = 0.657694036246975e-4;
-    const double offset = 5.004892881251670;
-    const ssize_t local_strip_index = 1;
-
-    double diffraction_angle_DG_param =
-        anglecalibration.diffraction_angle_from_DG_parameters(
-            center, conversion, offset, local_strip_index);
-
-    auto [distance_center, normal_distance, angle] =
-        anglecalibration.convert_to_EE_parameters(center, conversion, offset);
-
-    double diffraction_angle_EE_param =
-        anglecalibration.diffraction_angle_from_EE_parameters(
-            distance_center, normal_distance, angle, local_strip_index);
-
-    CHECK(diffraction_angle_EE_param ==
-          Catch::Approx(diffraction_angle_DG_param));
-
-    double strip_width_DG_param =
-        anglecalibration.angular_strip_width_from_DG_parameters(
-            center, conversion, offset, local_strip_index);
-
-    double strip_width_EE_param =
-        anglecalibration.angular_strip_width_from_EE_parameters(
-            distance_center, normal_distance, angle, local_strip_index);
-
-    CHECK(strip_width_DG_param == Catch::Approx(strip_width_EE_param));
-}
--- a/src/ClusterFile.cpp
+++ b/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
--- a/src/Fit.cpp
+++ b/src/Fit.cpp
@ -105,7 +105,7 @@ std::array<double, 3> gaus_init_par(const NDView<double, 1> x, const NDView<doub
    auto delta = x[1] - x[0];
    start_par[2] =
        std::count_if(y.begin(), y.end(),
-                      [e, delta](double val) { return val > *e / 2; }) *
+                      [e](double val) { return val > *e / 2; }) *
        delta / 2.35;

    return start_par;
--- a/src/Hdf5FileReader.test.cpp
+++ b/src/Hdf5FileReader.test.cpp
@ -1,110 +0,0 @@
-/************************************************
- * @file Hdf5FileReader.test.cpp
- * @short test case for reading hdf5 files
- ***********************************************/
-
-#include <filesystem>
-
-#include "test_config.hpp"
-
-#include <H5Cpp.h>
-
-#include "aare/Hdf5FileReader.hpp"
-#include "aare/NDArray.hpp"
-
-#include <catch2/catch_all.hpp>
-#include <catch2/catch_test_macros.hpp>
-#include <catch2/matchers/catch_matchers_floating_point.hpp>
-
-using namespace aare;
-
-TEST_CASE("read hdf5 file", "[.hdf5file][.files]") {
-
-    // TODO generalize datasetpath
-    std::string filename = test_data_path() / "AngleCalibration_Test_Data" /
-                           "ang1up_22keV_LaB60p3mm_48M_a_0320.h5";
-
-    REQUIRE(std::filesystem::exists(filename));
-
-    HDF5FileReader file_reader;
-
-    file_reader.open_file(filename);
-
-    auto dataset = file_reader.get_dataset("/entry/data/data");
-
-    auto shape = dataset.get_shape();
-
-    CHECK(shape[0] == 61440);
-
-    auto type = dataset.get_datatype();
-
-    const std::type_info *type_info = dataset.get_cpp_type();
-
-    CHECK(*type_info == typeid(uint32_t));
-
-    SECTION("read dataset into NDArray") {
-
-        NDArray<uint32_t, 1> dataset_array =
-            dataset.store_as_ndarray<uint32_t, 1>();
-
-        CHECK(dataset_array(0) == 866);
-        CHECK(dataset_array(61439) == 1436);
-    }
-
-    SECTION("read dataset into Frame") {
-        Frame frame = dataset.store_as_frame();
-        CHECK(*(reinterpret_cast<uint32_t *>(frame.pixel_ptr(0, 0))) == 866);
-        CHECK(*(reinterpret_cast<uint32_t *>(frame.pixel_ptr(0, 61439))) ==
-              1436);
-    }
-    SECTION("read subset of dataset") {
-        Frame frame(1, 10, Dtype(typeid(uint32_t)));
-
-        Subset subset{std::vector<hsize_t>{10}, std::vector<hsize_t>{10}};
-
-        dataset.read_into_buffer(frame.data(), subset);
-
-        CHECK(*(reinterpret_cast<uint32_t *>(frame.pixel_ptr(0, 0))) == 664);
-        CHECK(*(reinterpret_cast<uint32_t *>(frame.pixel_ptr(0, 9))) == 654);
-    }
-    /*
-    SECTION("read scalar") {
-    }
-    */
-}
-
-TEST_CASE("test datatypes", "[.hdf5file]") {
-
-    auto [dtype, expected_type_info] = GENERATE(
-        std::make_tuple(H5::DataType(H5::PredType::NATIVE_INT), &typeid(int)),
-        std::make_tuple(H5::DataType(H5::PredType::NATIVE_INT8),
-                        &typeid(int8_t)),
-        std::make_tuple(H5::DataType(H5::PredType::NATIVE_UINT16),
-                        &typeid(uint16_t)),
-        std::make_tuple(H5::DataType(H5::PredType::NATIVE_INT16),
-                        &typeid(int16_t)),
-        std::make_tuple(H5::DataType(H5::PredType::STD_U32LE),
-                        &typeid(uint32_t)),
-        std::make_tuple(H5::DataType(H5::PredType::STD_I32LE),
-                        &typeid(int32_t)),
-        std::make_tuple(H5::DataType(H5::PredType::NATIVE_INT32),
-                        &typeid(int32_t)),
-        std::make_tuple(H5::DataType(H5::PredType::IEEE_F64LE),
-                        &typeid(double)),
-        std::make_tuple(H5::DataType(H5::PredType::IEEE_F32LE), &typeid(float)),
-        std::make_tuple(H5::DataType(H5::PredType::NATIVE_FLOAT),
-                        &typeid(float)),
-        std::make_tuple(H5::DataType(H5::PredType::NATIVE_DOUBLE),
-                        &typeid(double)),
-        std::make_tuple(H5::DataType(H5::PredType::NATIVE_CHAR),
-                        &typeid(int8_t)));
-
-    const std::type_info &type_info = deduce_cpp_type(dtype);
-
-    CHECK(type_info == *expected_type_info);
-
-    // TODO: handle bit swapping
-    REQUIRE_THROWS(deduce_cpp_type(
-        H5::DataType(H5::PredType::IEEE_F32BE))); // does not convert from big
-                                                  // to little endian
-}
--- a/src/MythenFileReader.test.cpp
+++ b/src/MythenFileReader.test.cpp
@ -1,33 +0,0 @@
-/************************************************
- * @file MythenFileReader.test.cpp
- * @short test case for angle calibration class
- ***********************************************/
-
-#include "aare/MythenFileReader.hpp"
-
-#include <filesystem>
-
-#include "test_config.hpp"
-
-#include <catch2/catch_all.hpp>
-#include <catch2/catch_test_macros.hpp>
-#include <catch2/matchers/catch_matchers_floating_point.hpp>
-
-using namespace aare;
-
-TEST_CASE("test mythenfile_reader", "[.mythenfilereader][.files]") {
-
-    auto fpath = test_data_path() / "AngleCalibration_Test_Data";
-
-    REQUIRE(std::filesystem::exists(fpath));
-
-        MythenFileReader file_reader(fpath, "ang1up_22keV_LaB60p3mm_48M_a_0");
-
-    auto frame = file_reader.read_frame(320);
-
-    CHECK(frame.detector_angle == 0.99955);
-
-    CHECK(frame.channel_mask == std::array<uint8_t, 3>{0, 0, 1});
-
-    CHECK(frame.photon_counts.size() == 61440);
-}
--- a/src/RawFile.cpp
+++ b/src/RawFile.cpp
@ -1,6 +1,8 @@
 #include "aare/RawFile.hpp"
+#include "aare/algorithm.hpp"
 #include "aare/PixelMap.hpp"
 #include "aare/defs.hpp"
+#include "aare/logger.hpp"
 #include "aare/geo_helpers.hpp"

 #include <fmt/format.h>
@ -14,23 +16,14 @@ RawFile::RawFile(const std::filesystem::path &fname, const std::string &mode)
    : m_master(fname) {
    m_mode = mode;
    if (mode == "r") {
-
-        n_subfiles = find_number_of_subfiles(); // f0,f1...fn
-        n_subfile_parts =
-            m_master.geometry().col * m_master.geometry().row; // d0,d1...dn
-
-        
-
        find_geometry();
-
        if (m_master.roi()){
            m_geometry = update_geometry_with_roi(m_geometry, m_master.roi().value());
        }
-        
        open_subfiles();
    } else {
        throw std::runtime_error(LOCATION +
-                                 "Unsupported mode. Can only read RawFiles.");
+                                 " Unsupported mode. Can only read RawFiles.");
    }
 }

@ -67,12 +60,12 @@ void RawFile::read_into(std::byte *image_buf, size_t n_frames, DetectorHeader *h
        this->get_frame_into(m_current_frame++, image_buf, header);
        image_buf += bytes_per_frame();
        if(header) 
-            header+=n_mod();
+            header+=n_modules();
    }

 }

-size_t RawFile::n_mod() const { return n_subfile_parts; }
+size_t RawFile::n_modules() const { return m_master.n_modules(); }


 size_t RawFile::bytes_per_frame() {
@ -106,17 +99,11 @@ xy RawFile::geometry() { return m_master.geometry(); }

 void RawFile::open_subfiles() {
    if (m_mode == "r")
-        for (size_t i = 0; i != n_subfiles; ++i) {
-            auto v = std::vector<RawSubFile *>(n_subfile_parts);
-            for (size_t j = 0; j != n_subfile_parts; ++j) {
-                auto pos = m_geometry.module_pixel_0[j];
-                v[j] = new RawSubFile(m_master.data_fname(j, i),
-                                      m_master.detector_type(), pos.height,
-                                      pos.width, m_master.bitdepth(),
-                                      pos.row_index, pos.col_index);
-
-            }
-            subfiles.push_back(v);
+        for (size_t i = 0; i != n_modules(); ++i) {
+            auto pos = m_geometry.module_pixel_0[i];
+            m_subfiles.emplace_back(std::make_unique<RawSubFile>(
+                m_master.data_fname(i, 0), m_master.detector_type(), pos.height,
+                pos.width, m_master.bitdepth(), pos.row_index, pos.col_index));
        }
    else {
        throw std::runtime_error(LOCATION +
@ -141,18 +128,6 @@ DetectorHeader RawFile::read_header(const std::filesystem::path &fname) {
    return h;
 }

-int RawFile::find_number_of_subfiles() {
-    int n_files = 0;
-    // f0,f1...fn How many files is the data split into?
-    while (std::filesystem::exists(m_master.data_fname(0, n_files)))
-        n_files++; // increment after test
-
-#ifdef AARE_VERBOSE
-    fmt::print("Found: {} subfiles\n", n_files);
-#endif
-    return n_files;
-
-}

 RawMasterFile RawFile::master() const { return m_master; }

@ -168,7 +143,7 @@ void RawFile::find_geometry() {
    uint16_t c{};


-    for (size_t i = 0; i < n_subfile_parts; i++) {
+    for (size_t i = 0; i < n_modules(); i++) {
        auto h = read_header(m_master.data_fname(i, 0));
        r = std::max(r, h.row);
        c = std::max(c, h.column);
@ -210,70 +185,58 @@ size_t RawFile::bytes_per_pixel() const {
 }

 void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, DetectorHeader *header) {
+    LOG(logDEBUG) << "RawFile::get_frame_into(" << frame_index << ")";
    if (frame_index >= total_frames()) {
        throw std::runtime_error(LOCATION + "Frame number out of range");
    }
-    std::vector<size_t> frame_numbers(n_subfile_parts);
-    std::vector<size_t> frame_indices(n_subfile_parts, frame_index);
+    std::vector<size_t> frame_numbers(n_modules());
+    std::vector<size_t> frame_indices(n_modules(), frame_index);


    // sync the frame numbers

-    if (n_subfile_parts != 1) {
-        for (size_t part_idx = 0; part_idx != n_subfile_parts; ++part_idx) {
-            auto subfile_id = frame_index / m_master.max_frames_per_file();
-            if (subfile_id >= subfiles.size()) {
-                throw std::runtime_error(LOCATION +
-                                         " Subfile out of range. Possible missing data.");
-            }
-            frame_numbers[part_idx] =
-                subfiles[subfile_id][part_idx]->frame_number(
-                    frame_index % m_master.max_frames_per_file());
+    if (n_modules() != 1) { //if we have more than one module
+        for (size_t part_idx = 0; part_idx != n_modules(); ++part_idx) {
+            frame_numbers[part_idx] = m_subfiles[part_idx]->frame_number(frame_index); 
        }
+
        // 1. if frame number vector is the same break
-        while (std::adjacent_find(frame_numbers.begin(), frame_numbers.end(),
-                                  std::not_equal_to<>()) !=
-               frame_numbers.end()) {
+        while (!all_equal(frame_numbers)) {
+
            // 2. find the index of the minimum frame number,
            auto min_frame_idx = std::distance(
                frame_numbers.begin(),
                std::min_element(frame_numbers.begin(), frame_numbers.end()));
+
            // 3. increase its index and update its respective frame number
            frame_indices[min_frame_idx]++;
+
            // 4. if we can't increase its index => throw error
            if (frame_indices[min_frame_idx] >= total_frames()) {
                throw std::runtime_error(LOCATION +
                                         "Frame number out of range");
            }
-            auto subfile_id =
-                frame_indices[min_frame_idx] / m_master.max_frames_per_file();
+
            frame_numbers[min_frame_idx] =
-                subfiles[subfile_id][min_frame_idx]->frame_number(
-                    frame_indices[min_frame_idx] %
-                    m_master.max_frames_per_file());
+                m_subfiles[min_frame_idx]->frame_number(frame_indices[min_frame_idx]);
        }
    }

    if (m_master.geometry().col == 1) {
        // get the part from each subfile and copy it to the frame
-        for (size_t part_idx = 0; part_idx != n_subfile_parts; ++part_idx) {
+        for (size_t part_idx = 0; part_idx != n_modules(); ++part_idx) {
            auto corrected_idx = frame_indices[part_idx];
-            auto subfile_id = corrected_idx / m_master.max_frames_per_file();
-            if (subfile_id >= subfiles.size()) {
-                throw std::runtime_error(LOCATION +
-                                         " Subfile out of range. Possible missing data.");
-            }
-
+      
            // This is where we start writing
            auto offset = (m_geometry.module_pixel_0[part_idx].origin_y * m_geometry.pixels_x +
                m_geometry.module_pixel_0[part_idx].origin_x)*m_master.bitdepth()/8;

            if (m_geometry.module_pixel_0[part_idx].origin_x!=0)
-                throw std::runtime_error(LOCATION + "Implementation error. x pos not 0.");
+                throw std::runtime_error(LOCATION + " Implementation error. x pos not 0.");
            
-            //TODO! Risk for out of range access
-            subfiles[subfile_id][part_idx]->seek(corrected_idx % m_master.max_frames_per_file());
-            subfiles[subfile_id][part_idx]->read_into(frame_buffer + offset, header);
+            //TODO! What if the files don't match?
+            m_subfiles[part_idx]->seek(corrected_idx);
+            m_subfiles[part_idx]->read_into(frame_buffer + offset, header);
            if (header)
                ++header;
        }
@ -282,26 +245,21 @@ void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, Detect
        //TODO! should we read row by row?

        // create a buffer large enough to hold a full module
-
        auto bytes_per_part = m_master.pixels_y() * m_master.pixels_x() *
                              m_master.bitdepth() /
                              8; // TODO! replace with image_size_in_bytes
+
        auto *part_buffer = new std::byte[bytes_per_part];

        // TODO! if we have many submodules we should reorder them on the module
        // level

-        for (size_t part_idx = 0; part_idx != n_subfile_parts; ++part_idx) {
+        for (size_t part_idx = 0; part_idx != n_modules(); ++part_idx) {
            auto pos = m_geometry.module_pixel_0[part_idx];
            auto corrected_idx = frame_indices[part_idx];
-            auto subfile_id = corrected_idx / m_master.max_frames_per_file();
-            if (subfile_id >= subfiles.size()) {
-                throw std::runtime_error(LOCATION +
-                                         " Subfile out of range. Possible missing data.");
-            }

-            subfiles[subfile_id][part_idx]->seek(corrected_idx % m_master.max_frames_per_file());
-            subfiles[subfile_id][part_idx]->read_into(part_buffer, header);
+            m_subfiles[part_idx]->seek(corrected_idx);
+            m_subfiles[part_idx]->read_into(part_buffer, header);
            if(header)
                ++header;

@ -321,6 +279,7 @@ void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, Detect
        }
        delete[] part_buffer;
    }
+        
 }

 std::vector<Frame> RawFile::read_n(size_t n_frames) {
@ -337,27 +296,8 @@ size_t RawFile::frame_number(size_t frame_index) {
    if (frame_index >= m_master.frames_in_file()) {
        throw std::runtime_error(LOCATION + " Frame number out of range");
    }
-    size_t subfile_id = frame_index / m_master.max_frames_per_file();
-    if (subfile_id >= subfiles.size()) {
-        throw std::runtime_error(
-            LOCATION + " Subfile out of range. Possible missing data.");
-    }
-    return subfiles[subfile_id][0]->frame_number(
-        frame_index % m_master.max_frames_per_file());
+    return m_subfiles[0]->frame_number(frame_index);
 }

-RawFile::~RawFile() {
-
-    // TODO! Fix this, for file closing
-    for (auto &vec : subfiles) {
-        for (auto *subfile : vec) {
-            delete subfile;
-        }
-    }
-}
-
-
-
-

 } // namespace aare
--- a/src/RawFile.test.cpp
+++ b/src/RawFile.test.cpp
@ -99,11 +99,11 @@ TEST_CASE("Read frame numbers from a raw file", "[.integration]") {
    }
 }

-TEST_CASE("Compare reading from a numpy file with a raw file", "[.integration]") {
-    auto fpath_raw = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json";
+TEST_CASE("Compare reading from a numpy file with a raw file", "[.files]") {
+    auto fpath_raw = test_data_path() / "raw/jungfrau" / "jungfrau_single_master_0.json";
    REQUIRE(std::filesystem::exists(fpath_raw));

-    auto fpath_npy = test_data_path() / "jungfrau" / "jungfrau_single_0.npy";
+    auto fpath_npy = test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy";
    REQUIRE(std::filesystem::exists(fpath_npy));

    File raw(fpath_raw, "r");
@ -113,6 +113,7 @@ TEST_CASE("Compare reading from a numpy file with a raw file", "[.integration]")
    CHECK(npy.total_frames() == 10);

    for (size_t i = 0; i < 10; ++i) {
+        CHECK(raw.tell() == i);
        auto raw_frame = raw.read_frame();
        auto npy_frame = npy.read_frame();
        CHECK((raw_frame.view<uint16_t>() == npy_frame.view<uint16_t>()));
--- a/src/RawMasterFile.cpp
+++ b/src/RawMasterFile.cpp
@ -140,6 +140,10 @@ std::optional<size_t> RawMasterFile::number_of_rows() const {

 xy RawMasterFile::geometry() const { return m_geometry; }

+size_t RawMasterFile::n_modules() const {
+    return m_geometry.row * m_geometry.col;
+}
+
 std::optional<uint8_t> RawMasterFile::quad() const { return m_quad; }

 // optional values, these may or may not be present in the master file
--- a/src/RawSubFile.cpp
+++ b/src/RawSubFile.cpp
@ -1,9 +1,15 @@
 #include "aare/RawSubFile.hpp"
 #include "aare/PixelMap.hpp"
+#include "aare/algorithm.hpp"
 #include "aare/utils/ifstream_helpers.hpp"
+#include "aare/logger.hpp"
+
+
 #include <cstring> // memcpy
 #include <fmt/core.h>
 #include <iostream>
+#include <regex>
+



@ -12,51 +18,51 @@ namespace aare {
 RawSubFile::RawSubFile(const std::filesystem::path &fname,
                       DetectorType detector, size_t rows, size_t cols,
                       size_t bitdepth, uint32_t pos_row, uint32_t pos_col)
-    : m_detector_type(detector), m_bitdepth(bitdepth), m_fname(fname),
+    : m_detector_type(detector), m_bitdepth(bitdepth),
      m_rows(rows), m_cols(cols),
      m_bytes_per_frame((m_bitdepth / 8) * m_rows * m_cols), m_pos_row(pos_row),
      m_pos_col(pos_col) {
+
+    LOG(logDEBUG) << "RawSubFile::RawSubFile()";  
    if (m_detector_type == DetectorType::Moench03_old) {
        m_pixel_map = GenerateMoench03PixelMap();
    } else if (m_detector_type == DetectorType::Eiger && m_pos_row % 2 == 0) {
        m_pixel_map = GenerateEigerFlipRowsPixelMap();
    }

-    if (std::filesystem::exists(fname)) {
-        m_num_frames = std::filesystem::file_size(fname) /
-                   (sizeof(DetectorHeader) + rows * cols * bitdepth / 8);
-    } else {
-        throw std::runtime_error(
-            LOCATION + fmt::format("File {} does not exist", m_fname.string()));
-    }

-    // fp = fopen(m_fname.string().c_str(), "rb");
-    m_file.open(m_fname, std::ios::binary);
-    if (!m_file.is_open()) {
-        throw std::runtime_error(
-            LOCATION + fmt::format("Could not open file {}", m_fname.string()));
-    }
-
-#ifdef AARE_VERBOSE
-    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));
-#endif
+    parse_fname(fname);
+    scan_files();
+    open_file(m_current_file_index); // open the first file
 }

 void RawSubFile::seek(size_t frame_index) {
-    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));
+    LOG(logDEBUG) << "RawSubFile::seek(" << frame_index << ")";
+    if (frame_index >= m_total_frames) {
+        throw std::runtime_error(LOCATION + " Frame index out of range: " +
+                                 std::to_string(frame_index));
    }
-    m_file.seekg((sizeof(DetectorHeader) + bytes_per_frame()) * 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 + sizeof(DetectorHeader));
+    m_file.seekg(byte_offset);
 }

 size_t RawSubFile::tell() {
-    return m_file.tellg() / (sizeof(DetectorHeader) + bytes_per_frame());
+    LOG(logDEBUG) << "RawSubFile::tell():" << m_current_frame_index;
+    return m_current_frame_index;
 }

 void RawSubFile::read_into(std::byte *image_buf, DetectorHeader *header) {
+    LOG(logDEBUG) << "RawSubFile::read_into()";
+
    if (header) {
        m_file.read(reinterpret_cast<char *>(header), sizeof(DetectorHeader));
    } else {
@ -90,6 +96,13 @@ void RawSubFile::read_into(std::byte *image_buf, DetectorHeader *header) {
    if (m_file.fail()){
        throw std::runtime_error(LOCATION + ifstream_error_msg(m_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 RawSubFile::read_into(std::byte *image_buf, size_t n_frames, DetectorHeader *header) {
@ -130,4 +143,69 @@ size_t RawSubFile::frame_number(size_t frame_index) {
    return h.frameNumber;
 }

+void RawSubFile::parse_fname(const std::filesystem::path &fname) {
+    LOG(logDEBUG) << "RawSubFile::parse_fname()";
+    // data has the format: /path/too/data/jungfrau_single_d0_f1_0.raw
+    // d0 is the module index, will not change for this file
+    // f1 is the file index - thi is the one we need
+    // 0 is the measurement index, will not change
+    m_path = fname.parent_path();
+    m_base_name = fname.filename();
+
+    // Regex to extract numbers after 'd' and 'f'
+    std::regex pattern(R"(^(.*_d)(\d+)(_f)(\d+)(_\d+\.raw)$)");
+    std::smatch match;
+
+    if (std::regex_match(m_base_name, match, pattern)) {
+        m_offset = std::stoi(match[4].str()); // find the first file index in case of a truncated series
+        m_base_name = match[1].str() + match[2].str() + match[3].str() + "{}" + match[5].str();
+        LOG(logDEBUG) << "Base name: " << m_base_name;
+        LOG(logDEBUG) << "Offset: " << m_offset;
+        LOG(logDEBUG) << "Path: " << m_path.string();
+    } else {
+        throw std::runtime_error(
+            LOCATION + fmt::format("Could not parse file name {}", fname.string()));
+    }
+}
+
+std::filesystem::path RawSubFile::fpath(size_t file_index) const {
+    auto fname = fmt::format(m_base_name, file_index);
+    return m_path / fname;
+}
+
+void RawSubFile::open_file(size_t file_index) {
+    m_file.close();
+    auto fname = fpath(file_index+m_offset);
+    LOG(logDEBUG) << "RawSubFile::open_file():  " << fname.string();
+    m_file.open(fname, std::ios::binary);
+    if (!m_file.is_open()) {
+        throw std::runtime_error(
+            LOCATION + fmt::format("Could not open file {}", fpath(file_index).string()));
+    }
+    m_current_file_index = file_index;
+}
+
+void RawSubFile::scan_files() {
+    LOG(logDEBUG) << "RawSubFile::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 + sizeof(DetectorHeader));
+        m_last_frame_in_file.push_back(n_frames);
+        LOG(logDEBUG) << "Found: " << n_frames << " frames in file: " << fpath(file_index).string();
+        ++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);
+    if(m_last_frame_in_file.empty()){
+        m_total_frames = 0;
+    }else{
+        m_total_frames = m_last_frame_in_file.back();
+    }
+}
+
 } // namespace aare
--- a/src/RawSubFile.test.cpp
+++ b/src/RawSubFile.test.cpp
@ -0,0 +1,76 @@
+#include "aare/RawSubFile.hpp"
+#include "aare/File.hpp"
+#include "aare/NDArray.hpp"
+#include <catch2/catch_test_macros.hpp>
+#include "test_config.hpp"
+
+using namespace aare;
+
+TEST_CASE("Read frames directly from a RawSubFile", "[.files]"){
+    auto fpath_raw = test_data_path() / "raw/jungfrau" / "jungfrau_single_d0_f0_0.raw";
+    REQUIRE(std::filesystem::exists(fpath_raw));
+
+    RawSubFile f(fpath_raw, DetectorType::Jungfrau, 512, 1024, 16);
+    REQUIRE(f.rows() == 512);
+    REQUIRE(f.cols() == 1024);
+    REQUIRE(f.pixels_per_frame() == 512 * 1024);
+    REQUIRE(f.bytes_per_frame() == 512 * 1024 * 2);
+    REQUIRE(f.bytes_per_pixel() == 2);
+
+
+    auto fpath_npy = test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy";
+    REQUIRE(std::filesystem::exists(fpath_npy));
+
+    //Numpy file with the same data to use as reference
+    File npy(fpath_npy, "r");
+
+    CHECK(f.frames_in_file() == 10);
+    CHECK(npy.total_frames() == 10);
+
+
+    DetectorHeader header{};
+    NDArray<uint16_t, 2> image({static_cast<ssize_t>(f.rows()), static_cast<ssize_t>(f.cols())});
+    for (size_t i = 0; i < 10; ++i) {
+        CHECK(f.tell() == i);
+        f.read_into(image.buffer(), &header);
+        auto npy_frame = npy.read_frame();
+        CHECK((image.view() == npy_frame.view<uint16_t>()));
+    }
+}
+
+TEST_CASE("Read frames directly from a RawSubFile starting at the second file", "[.files]"){
+    // we know this file has 10 frames with frame numbers 1 to 10
+    // f0 1,2,3
+    // f1 4,5,6 <-- starting here
+    // f2 7,8,9
+    // f3 10
+    
+    auto fpath_raw = test_data_path() / "raw/jungfrau" / "jungfrau_single_d0_f1_0.raw";
+    REQUIRE(std::filesystem::exists(fpath_raw));
+
+    RawSubFile f(fpath_raw, DetectorType::Jungfrau, 512, 1024, 16);
+
+
+    auto fpath_npy = test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy";
+    REQUIRE(std::filesystem::exists(fpath_npy));
+
+    //Numpy file with the same data to use as reference
+    File npy(fpath_npy, "r");
+    npy.seek(3);
+
+    CHECK(f.frames_in_file() == 7);
+    CHECK(npy.total_frames() == 10);
+
+
+    DetectorHeader header{};
+    NDArray<uint16_t, 2> image({static_cast<ssize_t>(f.rows()), static_cast<ssize_t>(f.cols())});
+    for (size_t i = 0; i < 7; ++i) {
+        CHECK(f.tell() == i);
+        f.read_into(image.buffer(), &header);
+        // frame numbers start at 1 frame index at 0 
+        // adding 3 + 1 to verify the frame number
+        CHECK(header.frameNumber == i + 4); 
+        auto npy_frame = npy.read_frame();
+        CHECK((image.view() == npy_frame.view<uint16_t>()));
+    }
+}
--- a/src/algorithm.test.cpp
+++ b/src/algorithm.test.cpp
@ -160,3 +160,36 @@ TEST_CASE("cumsum works with negative numbers", "[algorithm]") {
    REQUIRE(result[3] == -6);
    REQUIRE(result[4] == -10);
 }
+
+
+TEST_CASE("cumsum on an empty vector", "[algorithm]") {
+    std::vector<double> vec = {};
+    auto result = aare::cumsum(vec);
+    REQUIRE(result.size() == 0);
+
+}
+
+TEST_CASE("All equal on an empty vector is false", "[algorithm]") {
+    std::vector<int> vec = {};
+    REQUIRE(aare::all_equal(vec) == false);
+}
+
+TEST_CASE("All equal on a vector with 1 element is true", "[algorithm]") {
+    std::vector<int> vec = {1};
+    REQUIRE(aare::all_equal(vec) == true);
+}
+
+TEST_CASE("All equal on a vector with 2 elements is true", "[algorithm]") {
+    std::vector<int> vec = {1, 1};
+    REQUIRE(aare::all_equal(vec) == true);
+}
+
+TEST_CASE("All equal on a vector with two different  elements is false", "[algorithm]") {
+    std::vector<int> vec = {1, 2};
+    REQUIRE(aare::all_equal(vec) == false);
+}
+
+TEST_CASE("Last element is different", "[algorithm]") {
+    std::vector<int> vec = {1, 1, 1, 1, 2};
+    REQUIRE(aare::all_equal(vec) == false);
+}
Author	SHA1	Message	Date
Erik Fröjdh	1bc2fd770a	Binding 5x5, 7x7 and 9x9 clusters in python (#188 ) All checks were successful Build on RHEL8 / build (push) Successful in 2m55s Details Build on RHEL9 / build (push) Successful in 2m58s Details - New binding code with macros to bind all cluster templates - Simplified factory function on the python side - 5x5, 7x7 and 9x9 bindings in python	2025-06-05 08:57:59 +02:00
Erik Fröjdh	69964e08d5	Refactor cluster bindings (#185 ) All checks were successful Build on RHEL9 / build (push) Successful in 2m19s Details Build on RHEL8 / build (push) Successful in 2m34s Details - Split up the file for cluster bindings - new file names according to bind_ClassName.hpp	2025-06-03 08:43:40 +02:00
froejdh_e	9ecf4f4b44	merge All checks were successful Build on RHEL9 / build (push) Successful in 2m22s Details Build on RHEL8 / build (push) Successful in 2m30s Details	2025-05-22 11:23:57 +02:00
froejdh_e	f2a024644b	bumped version upload on release	2025-05-22 11:10:23 +02:00
Erik Fröjdh	9e1b8731b0	RawSubFile support multi file access (#173 ) This PR is a fix/improvement to a problem that Jonathan had. (#156) The original implementation opened all subfiles at once witch works for normal sized datasets but fails at a certain point (thousands of files). - This solution uses RawSubFile to manage the different file indicies and only opens the file we need - Added logger.h from slsDetectorPackage for debug printing (in production no messages should be visible)	2025-05-22 11:00:03 +02:00
Erik Fröjdh	a6eebbe9bd	removed extra const on return type, added cast (#177 ) All checks were successful Build on RHEL9 / build (push) Successful in 2m31s Details Build on RHEL8 / build (push) Successful in 2m34s Details Fixed warnings on apple clang: - removed extra const on return type - added cast to suppress a float to double conversion warning	2025-05-20 15:27:38 +02:00
Erik Fröjdh	fd0196f2fd	Developer (#164 ) All checks were successful Build on RHEL9 / build (push) Successful in 1m58s Details Build on RHEL8 / build (push) Successful in 2m22s Details - 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>	2025-04-22 16:41:48 +02:00
 @ -1 +1 @@
 .0.0
 .5.22