For 2025.5.22 release (#181)

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: AliceMazzoleni99 <alice.mazzoleni@psi.ch>
Co-authored-by: Mazzoleni Alice Francesca <mazzol_a@pc17378.psi.ch>
Co-authored-by: siebsi <sieb.patr@gmail.com>

.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:
@@ -24,13 +24,13 @@ jobs:
       - uses: actions/checkout@v4
 
       - name: Get conda
-        uses: conda-incubator/setup-miniconda@v3.0.4
+        uses: conda-incubator/setup-miniconda@v3
         with:
           python-version: ${{ matrix.python-version }}
+          environment-file: etc/dev-env.yml
+          miniforge-version: latest
           channels: conda-forge
-
-      - name: Prepare
-        run: conda install conda-build=24.9 conda-verify pytest anaconda-client
+          conda-remove-defaults: "true"
 
       - name: Enable upload
         run: conda config --set anaconda_upload yes

.github/workflows/build_conda.yml

@@ -24,14 +24,15 @@ jobs:
       - uses: actions/checkout@v4
 
       - name: Get conda
-        uses: conda-incubator/setup-miniconda@v3.0.4
+        uses: conda-incubator/setup-miniconda@v3
         with:
           python-version: ${{ matrix.python-version }}
+          environment-file: etc/dev-env.yml
+          miniforge-version: latest
           channels: conda-forge
+          conda-remove-defaults: "true"
 
-      - name: Prepare
-        run: conda install conda-build=24.9 conda-verify pytest anaconda-client
-
+  
       - name: Disable upload
         run: conda config --set anaconda_upload no
 

CMakeLists.txt

@@ -1,12 +1,17 @@
 cmake_minimum_required(VERSION 3.15)
 
 project(aare 
-    VERSION 1.0.0
     DESCRIPTION "Data processing library for PSI detectors"
     HOMEPAGE_URL "https://github.com/slsdetectorgroup/aare"
     LANGUAGES C CXX 
 )
 
+# Read VERSION file into project version
+set(VERSION_FILE "${CMAKE_CURRENT_SOURCE_DIR}/VERSION")
+file(READ "${VERSION_FILE}" VERSION_CONTENT)
+string(STRIP "${VERSION_CONTENT}" PROJECT_VERSION_STRING)
+set(PROJECT_VERSION ${PROJECT_VERSION_STRING})
+
 set(CMAKE_CXX_STANDARD 17) 
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_CXX_EXTENSIONS OFF)
@@ -74,6 +79,9 @@ 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)
@@ -85,6 +93,7 @@ if(AARE_BENCHMARKS)
 endif()
 
 
+
 set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
 
 if(AARE_FETCH_LMFIT)
@@ -389,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
 )
@@ -400,6 +408,9 @@ target_include_directories(aare_core PUBLIC
     "$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>" 
 )
 
+set(THREADS_PREFER_PTHREAD_FLAG ON)
+find_package(Threads REQUIRED)
+
 target_link_libraries(
     aare_core 
     PUBLIC 
@@ -408,6 +419,7 @@ target_link_libraries(
     ${STD_FS_LIB} # from helpers.cmake
     PRIVATE 
     aare_compiler_flags 
+    Threads::Threads
     $<BUILD_INTERFACE:lmfit>
 
 )
@@ -443,6 +455,7 @@ if(AARE_TESTS)
             ${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
 
         )

VERSION

@@ -0,0 +1 @@
+2025.5.22

conda-recipe/conda_build_config.yaml

@@ -1,28 +1,5 @@
 python:
     - 3.11
-    - 3.11
-    - 3.11
-    - 3.12
-    - 3.12
     - 3.12
     - 3.13
 
-
-    
-numpy: 
-    - 1.26
-    - 2.0
-    - 2.1
-    - 1.26
-    - 2.0
-    - 2.1
-    - 2.1
-
-
-zip_keys:
-    - python
-    - numpy
-
-pin_run_as_build:
-  numpy: x.x
-  python: x.x

conda-recipe/meta.yaml

@@ -1,11 +1,10 @@
+source:
+  path: ../
+
+{% set version = load_file_regex(load_file = 'VERSION', regex_pattern = '(\d+(?:\.\d+)*(?:[\+\w\.]+))').group(1) %}
 package:
   name: aare
-  version: 2025.4.22 #TODO! how to not duplicate this?
-
-
-
-
-
+  version: {{version}}
 
 source:
   path: ..
@@ -13,45 +12,39 @@ source:
 build:
   number: 0
   script:
-    - unset CMAKE_GENERATOR && {{ PYTHON }} -m pip install . -vv  # [not win]
-    - {{ PYTHON }} -m pip install . -vv  # [win]
+    - unset CMAKE_GENERATOR && {{ PYTHON }} -m pip install . -vv 
 
 requirements:
   build:
-    - python {{python}}
-    - numpy {{ numpy }}
     - {{ compiler('cxx') }}
-
-
-  host:
     - cmake
     - ninja
-    - python {{python}}
-    - numpy {{ numpy }}
+
+  host:
+    - python
     - pip
+    - numpy=2.1
     - scikit-build-core
     - pybind11 >=2.13.0
-    - fmt
-    - zeromq
-    - nlohmann_json
-    - catch2
+    - matplotlib # needed in host to solve the environment for run
 
   run:
-    - python {{python}}
-    - numpy {{ numpy }}
+    - python
+    - {{ pin_compatible('numpy') }}
     - matplotlib
+    
 
 
 test:
   imports:
     - aare
-  # requires:
-  #   - pytest
-  # source_files:
-  #   - tests
-  # commands:
-  #   - pytest tests
+  requires:
+    - pytest
+    - boost-histogram 
+  source_files:
+    - python/tests
+  commands:
+    - python -m pytest python/tests
 
 about:
-  summary: An example project built with pybind11 and scikit-build.
-  # license_file: LICENSE
+  summary: Data analysis library for hybrid pixel detectors from PSI 

etc/dev-env.yml

@@ -3,13 +3,11 @@ channels:
   - conda-forge
 dependencies:
   - anaconda-client
+  - conda-build
   - doxygen 
   - sphinx=7.1.2 
   - breathe 
-  - pybind11 
   - sphinx_rtd_theme 
   - furo 
-  - nlohmann_json 
   - zeromq 
-  - fmt 
-  - numpy
+

include/aare/ArrayExpr.hpp

@@ -1,22 +1,24 @@
 #pragma once
-#include <cstdint> //int64_t
-#include <cstddef> //size_t
+#include <cstdint> 
+#include <cstddef> 
 #include <array> 
-
 #include <cassert>
+#include "aare/defs.hpp"
+
+
 namespace aare {
 
-template <typename E, int64_t Ndim> class ArrayExpr {
+template <typename E, ssize_t Ndim> class ArrayExpr {
   public:
     static constexpr bool is_leaf = false;
 
     auto operator[](size_t i) const { return static_cast<E const &>(*this)[i]; }
     auto operator()(size_t i) const { return static_cast<E const &>(*this)[i]; }
     auto size() const { return static_cast<E const &>(*this).size(); }
-    std::array<int64_t, Ndim> shape() const { return static_cast<E const &>(*this).shape(); }
+    std::array<ssize_t, Ndim> shape() const { return static_cast<E const &>(*this).shape(); }
 };
 
-template <typename A, typename B, int64_t Ndim>
+template <typename A, typename B, ssize_t Ndim>
 class ArrayAdd : public ArrayExpr<ArrayAdd<A, B, Ndim>, Ndim> {
     const A &arr1_;
     const B &arr2_;
@@ -27,10 +29,10 @@ class ArrayAdd : public ArrayExpr<ArrayAdd<A, B, Ndim>, Ndim> {
     }
     auto operator[](int i) const { return arr1_[i] + arr2_[i]; }
     size_t size() const { return arr1_.size(); }
-    std::array<int64_t, Ndim> shape() const { return arr1_.shape(); }
+    std::array<ssize_t, Ndim> shape() const { return arr1_.shape(); }
 };
 
-template <typename A, typename B, int64_t Ndim>
+template <typename A, typename B, ssize_t Ndim>
 class ArraySub : public ArrayExpr<ArraySub<A, B, Ndim>, Ndim> {
     const A &arr1_;
     const B &arr2_;
@@ -41,10 +43,10 @@ class ArraySub : public ArrayExpr<ArraySub<A, B, Ndim>, Ndim> {
     }
     auto operator[](int i) const { return arr1_[i] - arr2_[i]; }
     size_t size() const { return arr1_.size(); }
-    std::array<int64_t, Ndim> shape() const { return arr1_.shape(); }
+    std::array<ssize_t, Ndim> shape() const { return arr1_.shape(); }
 };
 
-template <typename A, typename B, int64_t Ndim>
+template <typename A, typename B, ssize_t Ndim>
 class ArrayMul : public ArrayExpr<ArrayMul<A, B, Ndim>,Ndim> {
     const A &arr1_;
     const B &arr2_;
@@ -55,10 +57,10 @@ class ArrayMul : public ArrayExpr<ArrayMul<A, B, Ndim>,Ndim> {
     }
     auto operator[](int i) const { return arr1_[i] * arr2_[i]; }
     size_t size() const { return arr1_.size(); }
-    std::array<int64_t, Ndim> shape() const { return arr1_.shape(); }
+    std::array<ssize_t, Ndim> shape() const { return arr1_.shape(); }
 };
 
-template <typename A, typename B, int64_t Ndim>
+template <typename A, typename B, ssize_t Ndim>
 class ArrayDiv : public ArrayExpr<ArrayDiv<A, B, Ndim>, Ndim> {
     const A &arr1_;
     const B &arr2_;
@@ -69,27 +71,27 @@ class ArrayDiv : public ArrayExpr<ArrayDiv<A, B, Ndim>, Ndim> {
     }
     auto operator[](int i) const { return arr1_[i] / arr2_[i]; }
     size_t size() const { return arr1_.size(); }
-    std::array<int64_t, Ndim> shape() const { return arr1_.shape(); }
+    std::array<ssize_t, Ndim> shape() const { return arr1_.shape(); }
 };
 
 
 
-template <typename A, typename B, int64_t Ndim>
+template <typename A, typename B, ssize_t Ndim>
 auto operator+(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
     return ArrayAdd<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
 }
 
-template <typename A, typename B, int64_t Ndim>
+template <typename A, typename B, ssize_t Ndim>
 auto operator-(const ArrayExpr<A,Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
     return ArraySub<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
 }
 
-template <typename A, typename B, int64_t Ndim>
+template <typename A, typename B, ssize_t Ndim>
 auto operator*(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
     return ArrayMul<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
 }
 
-template <typename A, typename B, int64_t Ndim>
+template <typename A, typename B, ssize_t Ndim>
 auto operator/(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
     return ArrayDiv<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
 }

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

include/aare/FilePtr.hpp

@@ -18,8 +18,8 @@ class FilePtr {
     FilePtr(FilePtr &&other);
     FilePtr &operator=(FilePtr &&other);
     FILE *get();
-    int64_t tell();
-    void seek(int64_t offset, int whence = SEEK_SET) {
+    ssize_t tell();
+    void seek(ssize_t offset, int whence = SEEK_SET) {
         if (fseek(fp_, offset, whence) != 0)
             throw std::runtime_error("Error seeking in file");
     }

include/aare/Fit.hpp

@@ -15,6 +15,12 @@ NDArray<double, 1> gaus(NDView<double, 1> x, NDView<double, 1> par);
 double pol1(const double x, const double *par);
 NDArray<double, 1> pol1(NDView<double, 1> x, NDView<double, 1> par);
 
+double scurve(const double x, const double *par);
+NDArray<double, 1> scurve(NDView<double, 1> x, NDView<double, 1> par);
+
+double scurve2(const double x, const double *par);
+NDArray<double, 1> scurve2(NDView<double, 1> x, NDView<double, 1> par);
+
 } // namespace func
 
 
@@ -25,6 +31,9 @@ std::array<double, 3> gaus_init_par(const NDView<double, 1> x, const NDView<doub
 
 std::array<double, 2> pol1_init_par(const NDView<double, 1> x, const NDView<double, 1> y);
 
+std::array<double, 6> scurve_init_par(const NDView<double, 1> x, const NDView<double, 1> y); 
+std::array<double, 6> scurve2_init_par(const NDView<double, 1> x, const NDView<double, 1> y);
+
 static constexpr int DEFAULT_NUM_THREADS = 4;
 
 /**
@@ -38,7 +47,7 @@ NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y);
 /**
  * @brief Fit a 1D Gaussian to each pixel. Data layout [row, col, values]
  * @param x x values
- * @param y y vales, layout [row, col, values]
+ * @param y y values, layout [row, col, values]
  * @param n_threads number of threads to use
  */
 
@@ -51,7 +60,7 @@ NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y,
 /**
  * @brief Fit a 1D Gaussian with error estimates
  * @param x x values
- * @param y y vales, layout [row, col, values]
+ * @param y y values, layout [row, col, values]
  * @param y_err error in y, layout [row, col, values]
  * @param par_out output parameters
  * @param par_err_out output error parameters
@@ -64,7 +73,7 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
  * @brief Fit a 1D Gaussian to each pixel with error estimates. Data layout
  * [row, col, values]
  * @param x x values
- * @param y y vales, layout [row, col, values]
+ * @param y y values, layout [row, col, values]
  * @param y_err error in y, layout [row, col, values]
  * @param par_out output parameters, layout [row, col, values]
  * @param par_err_out output parameter errors, layout [row, col, values]
@@ -88,5 +97,19 @@ void fit_pol1(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
               NDView<double, 3> par_out, NDView<double, 3> par_err_out,NDView<double, 2> chi2_out,
               int n_threads = DEFAULT_NUM_THREADS);
 
+NDArray<double, 1> fit_scurve(NDView<double, 1> x, NDView<double, 1> y); 
+NDArray<double, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y, int n_threads); 
+void fit_scurve(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
+              NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2); 
+void fit_scurve(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
+              NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
+              int n_threads); 
 
+NDArray<double, 1> fit_scurve2(NDView<double, 1> x, NDView<double, 1> y); 
+NDArray<double, 3> fit_scurve2(NDView<double, 1> x, NDView<double, 3> y, int n_threads); 
+void fit_scurve2(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
+              NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2); 
+void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
+              NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
+              int n_threads); 
 } // namespace aare

include/aare/Frame.hpp

@@ -107,8 +107,8 @@ class Frame {
      * @return NDView<T, 2> 
      */
     template <typename T> NDView<T, 2> view() {
-        std::array<int64_t, 2> shape = {static_cast<int64_t>(m_rows),
-                                        static_cast<int64_t>(m_cols)};
+        std::array<ssize_t, 2> shape = {static_cast<ssize_t>(m_rows),
+                                        static_cast<ssize_t>(m_cols)};
         T *data = reinterpret_cast<T *>(m_data);
         return NDView<T, 2>(data, shape);
     }

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

include/aare/NDArray.hpp

@@ -22,10 +22,10 @@ TODO! Add expression templates for operators
 namespace aare {
 
 
-template <typename T, int64_t Ndim = 2>
+template <typename T, ssize_t Ndim = 2>
 class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
-    std::array<int64_t, Ndim> shape_;
-    std::array<int64_t, Ndim> strides_;
+    std::array<ssize_t, Ndim> shape_;
+    std::array<ssize_t, Ndim> strides_;
     size_t size_{};
     T *data_;
 
@@ -42,7 +42,7 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
      *
      * @param shape shape of the new NDArray
      */
-    explicit NDArray(std::array<int64_t, Ndim> shape)
+    explicit NDArray(std::array<ssize_t, Ndim> shape)
         : shape_(shape), strides_(c_strides<Ndim>(shape_)),
           size_(std::accumulate(shape_.begin(), shape_.end(), 1,
                                 std::multiplies<>())),
@@ -55,7 +55,7 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
      * @param shape shape of the new array
      * @param value value to initialize the array with
      */
-    NDArray(std::array<int64_t, Ndim> shape, T value) : NDArray(shape) {
+    NDArray(std::array<ssize_t, Ndim> shape, T value) : NDArray(shape) {
         this->operator=(value);
     }
 
@@ -186,22 +186,22 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
     }
 
     // TODO! is int the right type for index?
-    T &operator()(int64_t i) { return data_[i]; }
-    const T &operator()(int64_t i) const { return data_[i]; }
+    T &operator()(ssize_t i) { return data_[i]; }
+    const T &operator()(ssize_t i) const { return data_[i]; }
 
-    T &operator[](int64_t i) { return data_[i]; }
-    const T &operator[](int64_t i) const { return data_[i]; }
+    T &operator[](ssize_t i) { return data_[i]; }
+    const T &operator[](ssize_t i) const { return data_[i]; }
 
     T *data() { return data_; }
     std::byte *buffer() { return reinterpret_cast<std::byte *>(data_); }
     ssize_t size() const { return static_cast<ssize_t>(size_); }
     size_t total_bytes() const { return size_ * sizeof(T); }
-    std::array<int64_t, Ndim> shape() const noexcept { return shape_; }
-    int64_t shape(int64_t i) const noexcept { return shape_[i]; }
-    std::array<int64_t, Ndim> strides() const noexcept { return strides_; }
+    std::array<ssize_t, Ndim> shape() const noexcept { return shape_; }
+    ssize_t shape(ssize_t i) const noexcept { return shape_[i]; }
+    std::array<ssize_t, Ndim> strides() const noexcept { return strides_; }
     size_t bitdepth() const noexcept { return sizeof(T) * 8; }
 
-    std::array<int64_t, Ndim> byte_strides() const noexcept {
+    std::array<ssize_t, Ndim> byte_strides() const noexcept {
         auto byte_strides = strides_;
         for (auto &val : byte_strides)
             val *= sizeof(T);
@@ -228,7 +228,7 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
 };
 
 // Move assign
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 NDArray<T, Ndim> &
 NDArray<T, Ndim>::operator=(NDArray<T, Ndim> &&other) noexcept {
     if (this != &other) {
@@ -242,7 +242,7 @@ NDArray<T, Ndim>::operator=(NDArray<T, Ndim> &&other) noexcept {
     return *this;
 }
 
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const NDArray<T, Ndim> &other) {
     // check shape
     if (shape_ == other.shape_) {
@@ -254,7 +254,7 @@ NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const NDArray<T, Ndim> &other) {
     throw(std::runtime_error("Shape of ImageDatas must match"));
 }
 
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 NDArray<T, Ndim> &NDArray<T, Ndim>::operator-=(const NDArray<T, Ndim> &other) {
     // check shape
     if (shape_ == other.shape_) {
@@ -266,7 +266,7 @@ NDArray<T, Ndim> &NDArray<T, Ndim>::operator-=(const NDArray<T, Ndim> &other) {
     throw(std::runtime_error("Shape of ImageDatas must match"));
 }
 
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 NDArray<T, Ndim> &NDArray<T, Ndim>::operator*=(const NDArray<T, Ndim> &other) {
     // check shape
     if (shape_ == other.shape_) {
@@ -278,14 +278,14 @@ NDArray<T, Ndim> &NDArray<T, Ndim>::operator*=(const NDArray<T, Ndim> &other) {
     throw(std::runtime_error("Shape of ImageDatas must match"));
 }
 
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 NDArray<T, Ndim> &NDArray<T, Ndim>::operator&=(const T &mask) {
     for (auto it = begin(); it != end(); ++it)
         *it &= mask;
     return *this;
 }
 
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 NDArray<bool, Ndim> NDArray<T, Ndim>::operator>(const NDArray &other) {
     if (shape_ == other.shape_) {
         NDArray<bool, Ndim> result{shape_};
@@ -297,7 +297,7 @@ NDArray<bool, Ndim> NDArray<T, Ndim>::operator>(const NDArray &other) {
     throw(std::runtime_error("Shape of ImageDatas must match"));
 }
 
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 NDArray<T, Ndim> &NDArray<T, Ndim>::operator=(const NDArray<T, Ndim> &other) {
     if (this != &other) {
         delete[] data_;
@@ -310,7 +310,7 @@ NDArray<T, Ndim> &NDArray<T, Ndim>::operator=(const NDArray<T, Ndim> &other) {
     return *this;
 }
 
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 bool NDArray<T, Ndim>::operator==(const NDArray<T, Ndim> &other) const {
     if (shape_ != other.shape_)
         return false;
@@ -322,23 +322,23 @@ bool NDArray<T, Ndim>::operator==(const NDArray<T, Ndim> &other) const {
     return true;
 }
 
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 bool NDArray<T, Ndim>::operator!=(const NDArray<T, Ndim> &other) const {
     return !((*this) == other);
 }
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 NDArray<T, Ndim> &NDArray<T, Ndim>::operator++() {
     for (uint32_t i = 0; i < size_; ++i)
         data_[i] += 1;
     return *this;
 }
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 NDArray<T, Ndim> &NDArray<T, Ndim>::operator=(const T &value) {
     std::fill_n(data_, size_, value);
     return *this;
 }
 
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const T &value) {
     for (uint32_t i = 0; i < size_; ++i)
         data_[i] += value;
@@ -348,57 +348,57 @@ NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const T &value) {
 
 
 
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 NDArray<T, Ndim> NDArray<T, Ndim>::operator+(const T &value) {
     NDArray result = *this;
     result += value;
     return result;
 }
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 NDArray<T, Ndim> &NDArray<T, Ndim>::operator-=(const T &value) {
     for (uint32_t i = 0; i < size_; ++i)
         data_[i] -= value;
     return *this;
 }
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 NDArray<T, Ndim> NDArray<T, Ndim>::operator-(const T &value) {
     NDArray result = *this;
     result -= value;
     return result;
 }
 
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 NDArray<T, Ndim> &NDArray<T, Ndim>::operator/=(const T &value) {
     for (uint32_t i = 0; i < size_; ++i)
         data_[i] /= value;
     return *this;
 }
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 NDArray<T, Ndim> NDArray<T, Ndim>::operator/(const T &value) {
     NDArray result = *this;
     result /= value;
     return result;
 }
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 NDArray<T, Ndim> &NDArray<T, Ndim>::operator*=(const T &value) {
     for (uint32_t i = 0; i < size_; ++i)
         data_[i] *= value;
     return *this;
 }
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 NDArray<T, Ndim> NDArray<T, Ndim>::operator*(const T &value) {
     NDArray result = *this;
     result *= value;
     return result;
 }
-// template <typename T, int64_t Ndim> void NDArray<T, Ndim>::Print() {
+// template <typename T, ssize_t Ndim> void NDArray<T, Ndim>::Print() {
 //     if (shape_[0] < 20 && shape_[1] < 20)
 //         Print_all();
 //     else
 //         Print_some();
 // }
 
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 std::ostream &operator<<(std::ostream &os, const NDArray<T, Ndim> &arr) {
     for (auto row = 0; row < arr.shape(0); ++row) {
         for (auto col = 0; col < arr.shape(1); ++col) {
@@ -410,7 +410,7 @@ std::ostream &operator<<(std::ostream &os, const NDArray<T, Ndim> &arr) {
     return os;
 }
 
-template <typename T, int64_t Ndim> void NDArray<T, Ndim>::Print_all() {
+template <typename T, ssize_t Ndim> void NDArray<T, Ndim>::Print_all() {
     for (auto row = 0; row < shape_[0]; ++row) {
         for (auto col = 0; col < shape_[1]; ++col) {
             std::cout << std::setw(3);
@@ -419,7 +419,7 @@ template <typename T, int64_t Ndim> void NDArray<T, Ndim>::Print_all() {
         std::cout << "\n";
     }
 }
-template <typename T, int64_t Ndim> void NDArray<T, Ndim>::Print_some() {
+template <typename T, ssize_t Ndim> void NDArray<T, Ndim>::Print_some() {
     for (auto row = 0; row < 5; ++row) {
         for (auto col = 0; col < 5; ++col) {
             std::cout << std::setw(7);
@@ -429,7 +429,7 @@ template <typename T, int64_t Ndim> void NDArray<T, Ndim>::Print_some() {
     }
 }
 
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 void save(NDArray<T, Ndim> &img, std::string &pathname) {
     std::ofstream f;
     f.open(pathname, std::ios::binary);
@@ -437,9 +437,9 @@ void save(NDArray<T, Ndim> &img, std::string &pathname) {
     f.close();
 }
 
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 NDArray<T, Ndim> load(const std::string &pathname,
-                      std::array<int64_t, Ndim> shape) {
+                      std::array<ssize_t, Ndim> shape) {
     NDArray<T, Ndim> img{shape};
     std::ifstream f;
     f.open(pathname, std::ios::binary);

include/aare/NDView.hpp

@@ -14,10 +14,10 @@
 #include <vector>
 namespace aare {
 
-template <int64_t Ndim> using Shape = std::array<int64_t, Ndim>;
+template <ssize_t Ndim> using Shape = std::array<ssize_t, Ndim>;
 
 // TODO! fix mismatch between signed and unsigned
-template <int64_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;
@@ -25,41 +25,41 @@ template <int64_t Ndim> Shape<Ndim> make_shape(const std::vector<size_t> &shape)
     return arr;
 }
 
-template <int64_t Dim = 0, typename Strides> int64_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 <int64_t Dim = 0, typename Strides, typename... Ix>
-int64_t element_offset(const Strides &strides, int64_t i, Ix... index) {
+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 <int64_t Ndim> std::array<int64_t, Ndim> c_strides(const std::array<int64_t, Ndim> &shape) {
-    std::array<int64_t, Ndim> strides{};
+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 (int64_t i = Ndim - 1; i > 0; --i) {
+    for (ssize_t i = Ndim - 1; i > 0; --i) {
         strides[i - 1] = strides[i] * shape[i];
     }
     return strides;
 }
 
-template <int64_t Ndim> std::array<int64_t, Ndim> make_array(const std::vector<int64_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<int64_t, Ndim> arr{};
+    std::array<ssize_t, Ndim> arr{};
     std::copy_n(vec.begin(), Ndim, arr.begin());
     return arr;
 }
 
-template <typename T, int64_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;
     NDView(const NDView &) = default;
     NDView(NDView &&) = default;
 
-    NDView(T *buffer, std::array<int64_t, Ndim> shape)
+    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<>())) {}
 
-    // NDView(T *buffer, const std::vector<int64_t> &shape)
+    // 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<>())) {}
 
@@ -73,14 +73,14 @@ template <typename T, int64_t Ndim = 2> class NDView : public ArrayExpr<NDView<T
 
     ssize_t size() const { return static_cast<ssize_t>(size_); }
     size_t total_bytes() const { return size_ * sizeof(T); }
-    std::array<int64_t, Ndim> strides() const noexcept { return strides_; }
+    std::array<ssize_t, Ndim> strides() const noexcept { return strides_; }
 
     T *begin() { return buffer_; }
     T *end() { return buffer_ + size_; }
     T const *begin() const { return buffer_; }
     T const *end() const { return buffer_ + size_; }
-    T &operator()(int64_t i) const { return buffer_[i]; }
-    T &operator[](int64_t i) const { return buffer_[i]; }
+    T &operator()(ssize_t i) const { return buffer_[i]; }
+    T &operator[](ssize_t i) const { return buffer_[i]; }
 
     bool operator==(const NDView &other) const {
         if (size_ != other.size_)
@@ -136,15 +136,15 @@ template <typename T, int64_t Ndim = 2> class NDView : public ArrayExpr<NDView<T
     }
 
     auto &shape() const { return shape_; }
-    auto shape(int64_t i) const { return shape_[i]; }
+    auto shape(ssize_t i) const { return shape_[i]; }
 
     T *data() { return buffer_; }
     void print_all() const;
 
   private:
     T *buffer_{nullptr};
-    std::array<int64_t, Ndim> strides_{};
-    std::array<int64_t, Ndim> shape_{};
+    std::array<ssize_t, Ndim> strides_{};
+    std::array<ssize_t, Ndim> shape_{};
     uint64_t size_{};
 
     template <class BinaryOperation> NDView &elemenwise(T val, BinaryOperation op) {
@@ -160,7 +160,7 @@ template <typename T, int64_t Ndim = 2> class NDView : public ArrayExpr<NDView<T
         return *this;
     }
 };
-template <typename T, int64_t Ndim> void NDView<T, Ndim>::print_all() const {
+template <typename T, ssize_t Ndim> void NDView<T, Ndim>::print_all() const {
     for (auto row = 0; row < shape_[0]; ++row) {
         for (auto col = 0; col < shape_[1]; ++col) {
             std::cout << std::setw(3);
@@ -171,7 +171,7 @@ template <typename T, int64_t Ndim> void NDView<T, Ndim>::print_all() const {
 }
 
 
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 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) {
@@ -186,7 +186,7 @@ std::ostream& operator <<(std::ostream& os, const NDView<T, Ndim>& arr){
 
 template <typename T>
 NDView<T,1> make_view(std::vector<T>& vec){
-    return NDView<T,1>(vec.data(), {static_cast<int64_t>(vec.size())});
+    return NDView<T,1>(vec.data(), {static_cast<ssize_t>(vec.size())});
 }
 
 } // namespace aare

include/aare/NumpyFile.hpp

@@ -69,7 +69,7 @@ class NumpyFile : public FileInterface {
      */
     template <typename T, size_t NDim> NDArray<T, NDim> load() {
         NDArray<T, NDim> arr(make_shape<NDim>(m_header.shape));
-        if (fseek(fp, static_cast<int64_t>(header_size), SEEK_SET)) {
+        if (fseek(fp, static_cast<long>(header_size), SEEK_SET)) {
             throw std::runtime_error(LOCATION + "Error seeking to the start of the data");
         }
         size_t rc = fread(arr.data(), sizeof(T), arr.size(), fp);

include/aare/Pedestal.hpp

@@ -107,7 +107,7 @@ template <typename SUM_TYPE = double> class Pedestal {
         assert(frame.size() == m_rows * m_cols);
 
         // TODO! move away from m_rows, m_cols
-        if (frame.shape() != std::array<int64_t, 2>{m_rows, m_cols}) {
+        if (frame.shape() != std::array<ssize_t, 2>{m_rows, m_cols}) {
             throw std::runtime_error(
                 "Frame shape does not match pedestal shape");
         }
@@ -128,7 +128,7 @@ template <typename SUM_TYPE = double> class Pedestal {
         assert(frame.size() == m_rows * m_cols);
 
         // TODO! move away from m_rows, m_cols
-        if (frame.shape() != std::array<int64_t, 2>{m_rows, m_cols}) {
+        if (frame.shape() != std::array<ssize_t, 2>{m_rows, m_cols}) {
             throw std::runtime_error(
                 "Frame shape does not match pedestal shape");
         }

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();
 };

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;

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

include/aare/VarClusterFinder.hpp

@@ -28,7 +28,7 @@ template <typename T> class VarClusterFinder {
     };
 
   private:
-    const std::array<int64_t, 2> shape_;
+    const std::array<ssize_t, 2> shape_;
     NDView<T, 2> original_;
     NDArray<int, 2> labeled_;
     NDArray<int, 2> peripheral_labeled_;

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

include/aare/defs.hpp

@@ -204,23 +204,25 @@ 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{
-    int64_t xmin{};
-    int64_t xmax{};
-    int64_t ymin{};
-    int64_t ymax{};
+    ssize_t xmin{};
+    ssize_t xmax{};
+    ssize_t ymin{};
+    ssize_t ymax{};
   
-    int64_t height() const { return ymax - ymin; }
-    int64_t width() const { return xmax - xmin; }
-    bool contains(int64_t x, int64_t y) const {
+    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<int64_t>;
+using dynamic_shape = std::vector<ssize_t>;
 
 //TODO! Can we uniform enums between the libraries?
 

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,
+    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

pyproject.toml

@@ -1,10 +1,16 @@
+[tool.scikit-build.metadata.version]
+provider = "scikit_build_core.metadata.regex"
+input = "VERSION"
+regex = '^(?P<version>\d+(?:\.\d+)*(?:[\.\+\w]+)?)$'
+result = "{version}" 
+
 [build-system]
 requires = ["scikit-build-core>=0.10", "pybind11", "numpy"]
 build-backend = "scikit_build_core.build"
 
 [project]
 name = "aare"
-version = "2025.4.22"
+dynamic = ["version"]
 requires-python = ">=3.11"
 dependencies = [
     "numpy",

python/aare/__init__.py

@@ -15,7 +15,7 @@ from .ClusterFinder import ClusterFinder, ClusterCollector, ClusterFinderMT, Clu
 from .ClusterVector import ClusterVector
 
 
-from ._aare import fit_gaus, fit_pol1
+from ._aare import fit_gaus, fit_pol1, fit_scurve, fit_scurve2
 from ._aare import Interpolator
 from ._aare import calculate_eta2
 

python/aare/func.py

@@ -1 +1 @@
-from ._aare import gaus, pol1
+from ._aare import gaus, pol1, scurve, scurve2

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=}')

python/src/ctb_raw_file.hpp

@@ -34,7 +34,7 @@ m.def("adc_sar_05_decode64to16", [](py::array_t<uint8_t> input) {
     }
 
     //Create a 2D output array with the same shape as the input
-    std::vector<ssize_t> shape{input.shape(0), input.shape(1)/static_cast<int64_t>(bits_per_byte)};
+    std::vector<ssize_t> shape{input.shape(0), input.shape(1)/static_cast<ssize_t>(bits_per_byte)};
     py::array_t<uint16_t> output(shape);
 
     //Create a view of the input and output arrays
@@ -55,7 +55,7 @@ m.def("adc_sar_04_decode64to16", [](py::array_t<uint8_t> input) {
     }
 
     //Create a 2D output array with the same shape as the input
-    std::vector<ssize_t> shape{input.shape(0), input.shape(1)/static_cast<int64_t>(bits_per_byte)};
+    std::vector<ssize_t> shape{input.shape(0), input.shape(1)/static_cast<ssize_t>(bits_per_byte)};
     py::array_t<uint16_t> output(shape);
 
     //Create a view of the input and output arrays

python/src/file.hpp

@@ -198,7 +198,7 @@ void define_file_io_bindings(py::module &m) {
 
     py::class_<ROI>(m, "ROI")
         .def(py::init<>())
-        .def(py::init<int64_t, int64_t, int64_t, int64_t>(), py::arg("xmin"),
+        .def(py::init<ssize_t, ssize_t, ssize_t, ssize_t>(), py::arg("xmin"),
              py::arg("xmax"), py::arg("ymin"), py::arg("ymax"))
         .def_readwrite("xmin", &ROI::xmin)
         .def_readwrite("xmax", &ROI::xmax)

python/src/fit.hpp

@@ -55,6 +55,47 @@ void define_fit_bindings(py::module &m) {
         )",
         py::arg("x"), py::arg("par"));
 
+    m.def(
+        "scurve",
+        [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
+           py::array_t<double, py::array::c_style | py::array::forcecast> par) {
+            auto x_view = make_view_1d(x);
+            auto par_view = make_view_1d(par);
+            auto y = new NDArray<double, 1>{aare::func::scurve(x_view, par_view)};
+            return return_image_data(y);
+        },
+        R"(
+        Evaluate a 1D scurve function for all points in x using parameters par. 
+        
+        Parameters
+        ----------
+        x : array_like
+            The points at which to evaluate the scurve function.
+        par : array_like
+            The parameters of the scurve function. The first element is the background slope, the second element is the background intercept, the third element is the mean, the fourth element is the standard deviation, the fifth element is inflexion point count number, and the sixth element is C.     
+        )",
+        py::arg("x"), py::arg("par"));
+
+    m.def(
+        "scurve2",
+        [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
+           py::array_t<double, py::array::c_style | py::array::forcecast> par) {
+            auto x_view = make_view_1d(x);
+            auto par_view = make_view_1d(par);
+            auto y = new NDArray<double, 1>{aare::func::scurve2(x_view, par_view)};
+            return return_image_data(y);
+        },
+        R"(
+        Evaluate a 1D scurve2 function for all points in x using parameters par. 
+        
+        Parameters
+        ----------
+        x : array_like
+            The points at which to evaluate the scurve function.
+        par : array_like
+            The parameters of the scurve2 function. The first element is the background slope, the second element is the background intercept, the third element is the mean, the fourth element is the standard deviation, the fifth element is inflexion point count number, and the sixth element is C.     
+        )",
+        py::arg("x"), py::arg("par"));
 
     m.def(
         "fit_gaus",
@@ -235,6 +276,180 @@ n_threads : int, optional
         R"(
 Fit a 1D polynomial to data with error estimates.
 
+Parameters
+----------
+x : array_like
+    The x values.
+y : array_like
+    The y values.
+y_err : array_like
+    The error in the y values.
+n_threads : int, optional
+    The number of threads to use. Default is 4.
+)",
+        py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
+
+//=========
+    m.def(
+        "fit_scurve",
+        [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
+           py::array_t<double, py::array::c_style | py::array::forcecast> y,
+           int n_threads) {
+            if (y.ndim() == 3) {
+                auto par = new NDArray<double, 3>{};
+
+                auto x_view = make_view_1d(x);
+                auto y_view = make_view_3d(y);
+                *par = aare::fit_scurve(x_view, y_view, n_threads);
+                return return_image_data(par);
+            } else if (y.ndim() == 1) {
+                auto par = new NDArray<double, 1>{};
+                auto x_view = make_view_1d(x);
+                auto y_view = make_view_1d(y);
+                *par = aare::fit_scurve(x_view, y_view);
+                return return_image_data(par);
+            } else {
+                throw std::runtime_error("Data must be 1D or 3D");
+            }
+        },
+        py::arg("x"), py::arg("y"), py::arg("n_threads") = 4);
+
+    m.def(
+        "fit_scurve",
+        [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
+           py::array_t<double, py::array::c_style | py::array::forcecast> y,
+           py::array_t<double, py::array::c_style | py::array::forcecast> y_err,
+           int n_threads) {
+            if (y.ndim() == 3) {
+                auto par = new NDArray<double, 3>({y.shape(0), y.shape(1), 6});
+
+                auto par_err =
+                    new NDArray<double, 3>({y.shape(0), y.shape(1), 6});
+
+                auto y_view = make_view_3d(y);
+                auto y_view_err = make_view_3d(y_err);
+                auto x_view = make_view_1d(x);
+
+                auto chi2 = new NDArray<double, 2>({y.shape(0), y.shape(1)});
+
+                aare::fit_scurve(x_view, y_view, y_view_err, par->view(),
+                               par_err->view(), chi2->view(), n_threads);
+                return py::dict("par"_a = return_image_data(par),
+                                "par_err"_a = return_image_data(par_err), 
+                                "chi2"_a = return_image_data(chi2),
+                                "Ndf"_a = y.shape(2) - 2);
+
+
+            } else if (y.ndim() == 1) {
+                auto par = new NDArray<double, 1>({2});
+                auto par_err = new NDArray<double, 1>({2});
+
+                auto y_view = make_view_1d(y);
+                auto y_view_err = make_view_1d(y_err);
+                auto x_view = make_view_1d(x);
+
+                double chi2 = 0;
+
+                aare::fit_scurve(x_view, y_view, y_view_err, par->view(),
+                               par_err->view(), chi2);
+                return py::dict("par"_a = return_image_data(par),
+                                "par_err"_a = return_image_data(par_err),
+                                "chi2"_a = chi2, "Ndf"_a = y.size() - 2);
+
+            } else {
+                throw std::runtime_error("Data must be 1D or 3D");
+            }
+        },
+        R"(
+Fit a 1D polynomial to data with error estimates.
+
+Parameters
+----------
+x : array_like
+    The x values.
+y : array_like
+    The y values.
+y_err : array_like
+    The error in the y values.
+n_threads : int, optional
+    The number of threads to use. Default is 4.
+)",
+        py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
+    
+
+    m.def(
+        "fit_scurve2",
+        [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
+           py::array_t<double, py::array::c_style | py::array::forcecast> y,
+           int n_threads) {
+            if (y.ndim() == 3) {
+                auto par = new NDArray<double, 3>{};
+
+                auto x_view = make_view_1d(x);
+                auto y_view = make_view_3d(y);
+                *par = aare::fit_scurve2(x_view, y_view, n_threads);
+                return return_image_data(par);
+            } else if (y.ndim() == 1) {
+                auto par = new NDArray<double, 1>{};
+                auto x_view = make_view_1d(x);
+                auto y_view = make_view_1d(y);
+                *par = aare::fit_scurve2(x_view, y_view);
+                return return_image_data(par);
+            } else {
+                throw std::runtime_error("Data must be 1D or 3D");
+            }
+        },
+        py::arg("x"), py::arg("y"), py::arg("n_threads") = 4);
+
+    m.def(
+        "fit_scurve2",
+        [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
+           py::array_t<double, py::array::c_style | py::array::forcecast> y,
+           py::array_t<double, py::array::c_style | py::array::forcecast> y_err,
+           int n_threads) {
+            if (y.ndim() == 3) {
+                auto par = new NDArray<double, 3>({y.shape(0), y.shape(1), 6});
+
+                auto par_err =
+                    new NDArray<double, 3>({y.shape(0), y.shape(1), 6});
+
+                auto y_view = make_view_3d(y);
+                auto y_view_err = make_view_3d(y_err);
+                auto x_view = make_view_1d(x);
+
+                auto chi2 = new NDArray<double, 2>({y.shape(0), y.shape(1)});
+
+                aare::fit_scurve2(x_view, y_view, y_view_err, par->view(),
+                               par_err->view(), chi2->view(), n_threads);
+                return py::dict("par"_a = return_image_data(par),
+                                "par_err"_a = return_image_data(par_err), 
+                                "chi2"_a = return_image_data(chi2),
+                                "Ndf"_a = y.shape(2) - 2);
+
+
+            } else if (y.ndim() == 1) {
+                auto par = new NDArray<double, 1>({6});
+                auto par_err = new NDArray<double, 1>({6});
+
+                auto y_view = make_view_1d(y);
+                auto y_view_err = make_view_1d(y_err);
+                auto x_view = make_view_1d(x);
+
+                double chi2 = 0;
+
+                aare::fit_scurve2(x_view, y_view, y_view_err, par->view(),
+                               par_err->view(), chi2);
+                return py::dict("par"_a = return_image_data(par),
+                                "par_err"_a = return_image_data(par_err),
+                                "chi2"_a = chi2, "Ndf"_a = y.size() - 2);
+
+            } else {
+                throw std::runtime_error("Data must be 1D or 3D");
+            }
+        },
+        R"(
+Fit a 1D polynomial to data with error estimates.
+
 Parameters
 ----------
 x : array_like

python/src/np_helper.hpp

@@ -13,7 +13,7 @@ namespace py = pybind11;
 using namespace aare;
 
 // Pass image data back to python as a numpy array
-template <typename T, int64_t Ndim>
+template <typename T, ssize_t Ndim>
 py::array return_image_data(aare::NDArray<T, Ndim> *image) {
 
     py::capsule free_when_done(image, [](void *f) {

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);
 }

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

src/ClusterFile.cpp

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

src/FilePtr.cpp

@@ -21,7 +21,7 @@ FilePtr &FilePtr::operator=(FilePtr &&other) {
 
 FILE *FilePtr::get() { return fp_; }
 
-int64_t FilePtr::tell() {
+ssize_t FilePtr::tell() {
     auto pos = ftell(fp_);
     if (pos == -1)
         throw std::runtime_error(fmt::format("Error getting file position: {}", error_msg()));

src/Fit.cpp

@@ -34,6 +34,30 @@ NDArray<double, 1> pol1(NDView<double, 1> x, NDView<double, 1> par) {
     return y;
 }
 
+double scurve(const double x, const double * par) {
+    return (par[0] + par[1] * x) + 0.5 * (1 + erf((x - par[2]) / (sqrt(2) * par[3]))) * (par[4] + par[5] * (x - par[2]));
+}
+
+NDArray<double, 1> scurve(NDView<double, 1> x, NDView<double, 1> par) {
+    NDArray<double, 1> y({x.shape()}, 0);
+    for (ssize_t i = 0; i < x.size(); i++) {
+        y(i) = scurve(x(i), par.data());
+    }
+    return y;
+}
+
+double scurve2(const double x, const double * par) {
+    return (par[0] + par[1] * x) + 0.5 * (1 - erf((x - par[2]) / (sqrt(2) * par[3]))) * (par[4] + par[5] * (x - par[2]));
+}
+
+NDArray<double, 1> scurve2(NDView<double, 1> x, NDView<double, 1> par) {
+    NDArray<double, 1> y({x.shape()}, 0);
+    for (ssize_t i = 0; i < x.size(); i++) {
+        y(i) = scurve2(x(i), par.data());
+    }
+    return y;
+}
+
 } // namespace func
 
 NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y) {
@@ -81,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;
@@ -273,4 +297,229 @@ NDArray<double, 3> fit_pol1(NDView<double, 1> x, NDView<double, 3> y,
     return result;
 }
 
+// ~~ S-CURVES ~~
+
+// SCURVE --
+std::array<double, 6> scurve_init_par(const NDView<double, 1> x, const NDView<double, 1> y){
+        // Estimate the initial parameters for the fit
+        std::array<double, 6> start_par{0, 0, 0, 0, 0, 0};
+
+        auto ymax = std::max_element(y.begin(), y.end());
+        auto ymin = std::min_element(y.begin(), y.end());
+        start_par[4] = *ymin + (*ymax - *ymin) / 2;
+        
+        // Find the first x where the corresponding y value is above the threshold (start_par[4])
+        for (ssize_t i = 0; i < y.size(); ++i) {
+            if (y[i] >= start_par[4]) {
+                start_par[2] = x[i];
+                break; // Exit the loop after finding the first valid x
+            }
+        }
+
+        start_par[3] = 2 * sqrt(start_par[2]); 
+        start_par[0] = 100; 
+        start_par[1] = 0.25;
+        start_par[5] = 1;
+        return start_par;
+}
+
+// - No error
+NDArray<double, 1> fit_scurve(NDView<double, 1> x, NDView<double, 1> y) {
+    NDArray<double, 1> result = scurve_init_par(x, y);
+    lm_status_struct status;
+
+    lmcurve(result.size(), result.data(), x.size(), x.data(), y.data(),
+            aare::func::scurve, &lm_control_double, &status);
+
+    return result;
+}
+
+NDArray<double, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y, int n_threads) {
+    NDArray<double, 3> result({y.shape(0), y.shape(1), 6}, 0);
+
+    auto process = [&x, &y, &result](ssize_t first_row, ssize_t last_row) {
+        for (ssize_t row = first_row; row < last_row; row++) {
+            for (ssize_t col = 0; col < y.shape(1); col++) {
+                NDView<double, 1> values(&y(row, col, 0), {y.shape(2)});
+                auto res = fit_scurve(x, values);
+                result(row, col, 0) = res(0);
+                result(row, col, 1) = res(1);
+                result(row, col, 2) = res(2);
+                result(row, col, 3) = res(3);
+                result(row, col, 4) = res(4);
+                result(row, col, 5) = res(5);
+            }
+        }
+    };
+
+    auto tasks = split_task(0, y.shape(0), n_threads);
+    RunInParallel(process, tasks);
+    return result;
+}
+
+// - Error
+void fit_scurve(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
+              NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2) {
+
+    // Check that we have the correct sizes
+    if (y.size() != x.size() || y.size() != y_err.size() ||
+        par_out.size() != 6 || par_err_out.size() != 6) {
+        throw std::runtime_error("Data, x, data_err must have the same size "
+                                 "and par_out, par_err_out must have size 6");
+    }
+
+    lm_status_struct status;
+    par_out = scurve_init_par(x, y);
+    std::array<double, 36> cov = {0}; // size 6x6
+    // std::array<double, 4> cov{0, 0, 0, 0};
+
+    lmcurve2(par_out.size(), par_out.data(), par_err_out.data(), cov.data(),
+             x.size(), x.data(), y.data(), y_err.data(), aare::func::scurve,
+             &lm_control_double, &status);
+
+    // Calculate chi2
+    chi2 = 0;
+    for (ssize_t i = 0; i < y.size(); i++) {
+        chi2 += std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
+    }
+}
+
+void fit_scurve(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
+              NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
+              int n_threads) {
+
+    auto process = [&](ssize_t first_row, ssize_t last_row) {
+        for (ssize_t row = first_row; row < last_row; row++) {
+            for (ssize_t col = 0; col < y.shape(1); col++) {
+                NDView<double, 1> y_view(&y(row, col, 0), {y.shape(2)});
+                NDView<double, 1> y_err_view(&y_err(row, col, 0),
+                                             {y_err.shape(2)});
+                NDView<double, 1> par_out_view(&par_out(row, col, 0),
+                                               {par_out.shape(2)});
+                NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
+                                                   {par_err_out.shape(2)});
+
+                fit_scurve(x, y_view, y_err_view, par_out_view, par_err_out_view, chi2_out(row, col));
+
+            }
+        }
+    };
+
+    auto tasks = split_task(0, y.shape(0), n_threads);
+    RunInParallel(process, tasks);
+
+}
+
+// SCURVE2 ---
+
+std::array<double, 6> scurve2_init_par(const NDView<double, 1> x, const NDView<double, 1> y){
+        // Estimate the initial parameters for the fit
+        std::array<double, 6> start_par{0, 0, 0, 0, 0, 0};
+
+        auto ymax = std::max_element(y.begin(), y.end());
+        auto ymin = std::min_element(y.begin(), y.end());
+        start_par[4] = *ymin + (*ymax - *ymin) / 2;
+        
+        // Find the first x where the corresponding y value is above the threshold (start_par[4])
+        for (ssize_t i = 0; i < y.size(); ++i) {
+            if (y[i] <= start_par[4]) {
+                start_par[2] = x[i];
+                break; // Exit the loop after finding the first valid x
+            }
+        }
+
+        start_par[3] = 2 * sqrt(start_par[2]); 
+        start_par[0] = 100; 
+        start_par[1] = 0.25;
+        start_par[5] = -1;
+        return start_par;
+}
+
+// - No error
+NDArray<double, 1> fit_scurve2(NDView<double, 1> x, NDView<double, 1> y) {
+    NDArray<double, 1> result = scurve2_init_par(x, y);
+    lm_status_struct status;
+
+    lmcurve(result.size(), result.data(), x.size(), x.data(), y.data(),
+            aare::func::scurve2, &lm_control_double, &status);
+
+    return result;
+}
+
+NDArray<double, 3> fit_scurve2(NDView<double, 1> x, NDView<double, 3> y, int n_threads) {
+    NDArray<double, 3> result({y.shape(0), y.shape(1), 6}, 0);
+
+    auto process = [&x, &y, &result](ssize_t first_row, ssize_t last_row) {
+        for (ssize_t row = first_row; row < last_row; row++) {
+            for (ssize_t col = 0; col < y.shape(1); col++) {
+                NDView<double, 1> values(&y(row, col, 0), {y.shape(2)});
+                auto res = fit_scurve2(x, values);
+                result(row, col, 0) = res(0);
+                result(row, col, 1) = res(1);
+                result(row, col, 2) = res(2);
+                result(row, col, 3) = res(3);
+                result(row, col, 4) = res(4);
+                result(row, col, 5) = res(5);
+            }
+        }
+    };
+
+    auto tasks = split_task(0, y.shape(0), n_threads);
+    RunInParallel(process, tasks);
+    return result;
+}
+
+// - Error
+void fit_scurve2(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
+              NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2) {
+
+    // Check that we have the correct sizes
+    if (y.size() != x.size() || y.size() != y_err.size() ||
+        par_out.size() != 6 || par_err_out.size() != 6) {
+        throw std::runtime_error("Data, x, data_err must have the same size "
+                                 "and par_out, par_err_out must have size 6");
+    }
+
+    lm_status_struct status;
+    par_out = scurve2_init_par(x, y);
+    std::array<double, 36> cov = {0}; // size 6x6
+    // std::array<double, 4> cov{0, 0, 0, 0};
+
+    lmcurve2(par_out.size(), par_out.data(), par_err_out.data(), cov.data(),
+             x.size(), x.data(), y.data(), y_err.data(), aare::func::scurve2,
+             &lm_control_double, &status);
+
+    // Calculate chi2
+    chi2 = 0;
+    for (ssize_t i = 0; i < y.size(); i++) {
+        chi2 += std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
+    }
+}
+
+void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
+              NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
+              int n_threads) {
+
+    auto process = [&](ssize_t first_row, ssize_t last_row) {
+        for (ssize_t row = first_row; row < last_row; row++) {
+            for (ssize_t col = 0; col < y.shape(1); col++) {
+                NDView<double, 1> y_view(&y(row, col, 0), {y.shape(2)});
+                NDView<double, 1> y_err_view(&y_err(row, col, 0),
+                                             {y_err.shape(2)});
+                NDView<double, 1> par_out_view(&par_out(row, col, 0),
+                                               {par_out.shape(2)});
+                NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
+                                                   {par_err_out.shape(2)});
+
+                fit_scurve2(x, y_view, y_err_view, par_out_view, par_err_out_view, chi2_out(row, col));
+
+            }
+        }
+    };
+
+    auto tasks = split_task(0, y.shape(0), n_threads);
+    RunInParallel(process, tasks);
+
+}
+
 } // namespace aare

src/JungfrauDataFile.cpp

@@ -89,7 +89,7 @@ void JungfrauDataFile::seek(size_t frame_index) {
                             : frame_index;
     auto byte_offset = frame_offset * (m_bytes_per_frame + header_size);
     m_fp.seek(byte_offset);
-};
+}
 
 size_t JungfrauDataFile::tell() { return m_current_frame_index; }
 size_t JungfrauDataFile::total_frames() const { return m_total_frames; }
@@ -235,4 +235,4 @@ std::filesystem::path JungfrauDataFile::fpath(size_t file_index) const {
     return m_path / fname;
 }
 
-} // namespace aare
+} // namespace aare

src/NDArray.test.cpp

@@ -44,9 +44,9 @@ TEST_CASE("3D NDArray from NDView"){
     REQUIRE(image.size() == view.size());
     REQUIRE(image.data() != view.data());
 
-    for(int64_t i=0; i<image.shape(0); i++){
-        for(int64_t j=0; j<image.shape(1); j++){
-            for(int64_t k=0; k<image.shape(2); k++){
+    for(ssize_t i=0; i<image.shape(0); i++){
+        for(ssize_t j=0; j<image.shape(1); j++){
+            for(ssize_t k=0; k<image.shape(2); k++){
                 REQUIRE(image(i, j, k) == view(i, j, k));
             }
         }
@@ -54,7 +54,7 @@ TEST_CASE("3D NDArray from NDView"){
 }
 
 TEST_CASE("1D image") {
-    std::array<int64_t, 1> shape{{20}};
+    std::array<ssize_t, 1> shape{{20}};
     NDArray<short, 1> img(shape, 3);
     REQUIRE(img.size() == 20);
     REQUIRE(img(5) == 3);
@@ -71,7 +71,7 @@ TEST_CASE("Accessing a const object") {
 }
 
 TEST_CASE("Indexing of a 2D image") {
-    std::array<int64_t, 2> shape{{3, 7}};
+    std::array<ssize_t, 2> shape{{3, 7}};
     NDArray<long> img(shape, 5);
     for (uint32_t i = 0; i != img.size(); ++i) {
         REQUIRE(img(i) == 5);
@@ -114,7 +114,7 @@ TEST_CASE("Divide double by int") {
 }
 
 TEST_CASE("Elementwise multiplication of 3D image") {
-    std::array<int64_t, 3> shape{3, 4, 2};
+    std::array<ssize_t, 3> shape{3, 4, 2};
     NDArray<double, 3> a{shape};
     NDArray<double, 3> b{shape};
     for (uint32_t i = 0; i != a.size(); ++i) {
@@ -179,9 +179,9 @@ TEST_CASE("Compare two images") {
 }
 
 TEST_CASE("Size and shape matches") {
-    int64_t w = 15;
-    int64_t h = 75;
-    std::array<int64_t, 2> shape{w, h};
+    ssize_t w = 15;
+    ssize_t h = 75;
+    std::array<ssize_t, 2> shape{w, h};
     NDArray<double> a{shape};
     REQUIRE(a.size() == w * h);
     REQUIRE(a.shape() == shape);
@@ -224,7 +224,7 @@ TEST_CASE("Bitwise and on data") {
 
 
 TEST_CASE("Elementwise operations on images") {
-    std::array<int64_t, 2> shape{5, 5};
+    std::array<ssize_t, 2> shape{5, 5};
     double a_val = 3.0;
     double b_val = 8.0;
 

src/NDView.test.cpp

@@ -142,7 +142,7 @@ TEST_CASE("iterators") {
 //     for (int i = 0; i != 12; ++i) {
 //         vec.push_back(i);
 //     }
-//     std::vector<int64_t> shape{3, 4};
+//     std::vector<ssize_t> shape{3, 4};
 //     NDView<int, 2> data(vec.data(), shape);
 // }
 
@@ -151,8 +151,8 @@ TEST_CASE("divide with another span") {
     std::vector<int> vec1{3, 2, 1};
     std::vector<int> result{3, 6, 3};
 
-    NDView<int, 1> data0(vec0.data(), Shape<1>{static_cast<int64_t>(vec0.size())});
-    NDView<int, 1> data1(vec1.data(), Shape<1>{static_cast<int64_t>(vec1.size())});
+    NDView<int, 1> data0(vec0.data(), Shape<1>{static_cast<ssize_t>(vec0.size())});
+    NDView<int, 1> data1(vec1.data(), Shape<1>{static_cast<ssize_t>(vec1.size())});
 
     data0 /= data1;
 

src/NumpyFile.cpp

@@ -72,8 +72,8 @@ void NumpyFile::get_frame_into(size_t frame_number, std::byte *image_buf) {
     }
 }
 
-size_t NumpyFile::pixels_per_frame() { return m_pixels_per_frame; };
-size_t NumpyFile::bytes_per_frame() { return m_bytes_per_frame; };
+size_t NumpyFile::pixels_per_frame() { return m_pixels_per_frame; }
+size_t NumpyFile::bytes_per_frame() { return m_bytes_per_frame; }
 
 std::vector<Frame> NumpyFile::read_n(size_t n_frames) {
     // TODO: implement this in a more efficient way
@@ -197,4 +197,4 @@ void NumpyFile::load_metadata() {
     m_header = {dtype, fortran_order, shape};
 }
 
-} // namespace aare
+} // namespace aare

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,27 +16,18 @@ 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.");
     }
 }
 
-Frame RawFile::read_frame() { return get_frame(m_current_frame++); };
+Frame RawFile::read_frame() { return get_frame(m_current_frame++); }
 
 Frame RawFile::read_frame(size_t frame_number) {
     seek(frame_number);
@@ -52,13 +45,13 @@ void RawFile::read_into(std::byte *image_buf, size_t n_frames) {
 
 void RawFile::read_into(std::byte *image_buf) {
     return get_frame_into(m_current_frame++, image_buf);
-};
+}
 
 
 void RawFile::read_into(std::byte *image_buf, DetectorHeader *header) {
 
     return get_frame_into(m_current_frame++, image_buf, header);
-};
+}
 
 void RawFile::read_into(std::byte *image_buf, size_t n_frames, DetectorHeader *header) {
     // return get_frame_into(m_current_frame++, image_buf, header);
@@ -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() {
@@ -94,9 +87,9 @@ void RawFile::seek(size_t frame_index) {
                         frame_index, total_frames()));
     }
     m_current_frame = frame_index;
-};
+}
 
-size_t RawFile::tell() { return m_current_frame; };
+size_t RawFile::tell() { return m_current_frame; }
 
 size_t RawFile::total_frames() const { return m_master.frames_in_file(); }
 size_t RawFile::rows() const { return m_geometry.pixels_y; }
@@ -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());
-}
-
-RawFile::~RawFile() {
-
-    // TODO! Fix this, for file closing
-    for (auto &vec : subfiles) {
-        for (auto *subfile : vec) {
-            delete subfile;
-        }
-    }
+    return m_subfiles[0]->frame_number(frame_index);
 }
 
 
-
-
-
-} // namespace aare
+} // namespace aare

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>()));

src/RawMasterFile.cpp

@@ -87,7 +87,7 @@ int ScanParameters::start() const { return m_start; }
 int ScanParameters::stop() const { return m_stop; }
 void ScanParameters::increment_stop(){
     m_stop += 1;
-};
+}
 int ScanParameters::step() const { return m_step; }
 const std::string &ScanParameters::dac() const { return m_dac; }
 bool ScanParameters::enabled() const { return m_enabled; }
@@ -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
@@ -417,4 +421,4 @@ void RawMasterFile::parse_raw(const std::filesystem::path &fpath) {
     if(m_frames_in_file==0)
         m_frames_in_file = m_total_frames_expected;
 }
-} // namespace aare
+} // namespace aare

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

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>()));
+    }
+}

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);
+}

src/decode.cpp

@@ -26,8 +26,8 @@ void adc_sar_05_decode64to16(NDView<uint64_t, 2> input, NDView<uint16_t,2> outpu
         throw std::invalid_argument(LOCATION + " input and output shapes must match");
     }
 
-    for(int64_t i = 0; i < input.shape(0); i++){
-        for(int64_t j = 0; j < input.shape(1); j++){
+    for(ssize_t i = 0; i < input.shape(0); i++){
+        for(ssize_t j = 0; j < input.shape(1); j++){
             output(i,j) = adc_sar_05_decode64to16(input(i,j));
         }
     }
@@ -56,8 +56,8 @@ void adc_sar_04_decode64to16(NDView<uint64_t, 2> input, NDView<uint16_t,2> outpu
     if(input.shape() != output.shape()){
         throw std::invalid_argument(LOCATION + " input and output shapes must match");
     }
-    for(int64_t i = 0; i < input.shape(0); i++){
-        for(int64_t j = 0; j < input.shape(1); j++){
+    for(ssize_t i = 0; i < input.shape(0); i++){
+        for(ssize_t j = 0; j < input.shape(1); j++){
             output(i,j) = adc_sar_04_decode64to16(input(i,j));
         }
     }

update_version.py

@@ -0,0 +1,57 @@
+# SPDX-License-Identifier: LGPL-3.0-or-other
+# Copyright (C) 2021 Contributors to the Aare Package
+"""
+Script to update VERSION file with semantic versioning if provided as an argument, or with 0.0.0 if no argument is provided.
+"""
+
+import sys
+import os
+import re
+
+from packaging.version import Version, InvalidVersion
+
+
+SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
+
+def is_integer(value): 
+    try:
+        int(value)
+    except ValueError:
+        return False
+    else:
+        return True
+
+
+def get_version():
+
+    # Check at least one argument is passed
+    if len(sys.argv) < 2:
+        return "0.0.0"
+
+    version = sys.argv[1]
+
+    try:
+        v = Version(version)  # normalize check if version follows PEP 440 specification
+        
+        version_normalized = version.replace("-", ".")
+
+        version_normalized =  re.sub(r'0*(\d+)', lambda m : str(int(m.group(0))), version_normalized) #remove leading zeros 
+
+        return version_normalized 
+
+    except InvalidVersion as e:
+        print(f"Invalid version {version}. Version format must follow semantic versioning format of python PEP 440 version identification specification.")
+        sys.exit(1)
+    
+
+def write_version_to_file(version):
+    version_file_path = os.path.join(SCRIPT_DIR, "VERSION")
+    with open(version_file_path, "w") as version_file:
+        version_file.write(version)
+    print(f"Version {version} written to VERSION file.")
+
+# Main script
+if __name__ == "__main__":
+
+    version = get_version()
+    write_version_to_file(version)