some refactoring

- Fixing broken build of tests on RH8 by linking pthreads
- Removed extra ; causing warnings with -Wpedantic

.github/workflows/build_and_deploy_conda.yml

@@ -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)
@@ -57,6 +62,7 @@ option(AARE_FETCH_CATCH "Use FetchContent to download catch2" ON)
 option(AARE_FETCH_JSON "Use FetchContent to download nlohmann::json" ON)
 option(AARE_FETCH_ZMQ "Use FetchContent to download libzmq" ON)
 option(AARE_FETCH_LMFIT "Use FetchContent to download lmfit" ON)
+option(AARE_FETCH_HDF5 "Use FetchContent to download hdf5-devel" OFF)
 
 
 #Convenience option to use system libraries only (no FetchContent)
@@ -68,6 +74,7 @@ if(AARE_SYSTEM_LIBRARIES)
     set(AARE_FETCH_CATCH OFF CACHE BOOL "Disabled FetchContent for catch2" FORCE)
     set(AARE_FETCH_JSON OFF CACHE BOOL "Disabled FetchContent for nlohmann::json" FORCE)
     set(AARE_FETCH_ZMQ OFF CACHE BOOL "Disabled FetchContent for libzmq" FORCE)
+    set(AARE_FETCH_HDF5 OFF CACHE BOOL "Disabled FetchContent for hdf5" FORCE)
     # Still fetch lmfit when setting AARE_SYSTEM_LIBRARIES since this is not available
     # on conda-forge
 endif()
@@ -214,6 +221,23 @@ else()
     find_package(nlohmann_json 3.11.3 REQUIRED)
 endif()
 
+if(AARE_FETCH_HDF5)
+    message(FATAL_ERROR "Fetching HDF5 via FetchContent is not supported here. Please install it via your system. 
+    For Ubuntu: sudo apt install libhdf5-dev
+    For Red Hat: sudo dnf install hdf5-devel
+    For MacOS: brew install hdf5")
+else()
+    find_package(HDF5 QUIET COMPONENTS CXX)
+    if (HDF5_FOUND)
+        message(STATUS "Found HDF5: ${HDF5_INCLUDE_DIRS}")
+    else()
+        message(FATAL_ERROR "HDF5 was NOT found! Please install it via your system  
+            For Ubuntu: sudo apt install libhdf5-dev
+            For Red Hat: sudo dnf install hdf5-devel
+            For MacOS: brew install hdf5")
+    endif()
+endif()
+
 include(GNUInstallDirs)
 
 # If conda build, always set lib dir to 'lib'
@@ -326,10 +350,6 @@ if(AARE_ASAN)
     )
 endif()
 
-
-
-
-
 if(AARE_TESTS)
     enable_testing()
     add_subdirectory(tests)
@@ -339,6 +359,7 @@ endif()
 ###------------------------------------------------------------------------------------------
 
 set(PUBLICHEADERS
+    include/aare/AngleCalibration.hpp
     include/aare/ArrayExpr.hpp
     include/aare/CalculateEta.hpp
     include/aare/Cluster.hpp
@@ -353,10 +374,14 @@ set(PUBLICHEADERS
     include/aare/Fit.hpp
     include/aare/FileInterface.hpp
     include/aare/FilePtr.hpp
+    include/aare/FlatField.hpp
     include/aare/Frame.hpp
     include/aare/GainMap.hpp
     include/aare/geo_helpers.hpp
+    include/aare/Hdf5FileReader.hpp
     include/aare/JungfrauDataFile.hpp
+    include/aare/MythenDetectorSpecifications.hpp
+    include/aare/MythenFileReader.hpp
     include/aare/NDArray.hpp
     include/aare/NDView.hpp
     include/aare/NumpyFile.hpp
@@ -372,6 +397,7 @@ set(PUBLICHEADERS
 
 
 set(SourceFiles
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/AngleCalibration.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/CtbRawFile.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/defs.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.cpp
@@ -400,14 +426,19 @@ 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 
     fmt::fmt 
     nlohmann_json::nlohmann_json 
+    HDF5::HDF5
     ${STD_FS_LIB} # from helpers.cmake
     PRIVATE 
     aare_compiler_flags 
+    Threads::Threads
     $<BUILD_INTERFACE:lmfit>
 
 )
@@ -424,12 +455,14 @@ endif()
 if(AARE_TESTS)
     set(TestSources
             ${CMAKE_CURRENT_SOURCE_DIR}/src/algorithm.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/AngleCalibration.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/defs.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/decode.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/geo_helpers.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/RawMasterFile.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/Hdf5FileReader.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/NDArray.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/NDView.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFinder.test.cpp
@@ -440,6 +473,7 @@ if(AARE_TESTS)
             ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFinderMT.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/Pedestal.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/JungfrauDataFile.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/MythenFileReader.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.test.cpp

VERSION

@@ -0,0 +1 @@
+0.0.0

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/AngleCalibration.hpp

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

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/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/FlatField.hpp

@@ -0,0 +1,112 @@
+
+/**
+ * stores flatfield for angle calibration
+ */
+
+#pragma once
+#include <cmath>
+#include <cstdint>
+
+#include <fstream>
+#include <iostream>
+#include <memory>
+#include <sstream>
+#include <string>
+
+#include "MythenDetectorSpecifications.hpp"
+#include "NDArray.hpp"
+
+namespace aare {
+// TODO maybe template now its uint32
+class FlatField {
+
+  public:
+    FlatField(std::shared_ptr<MythenDetectorSpecifications> mythen_detector_)
+        : mythen_detector(mythen_detector_) {
+
+        flat_field = NDArray<uint32_t, 1>(
+            std::array<ssize_t, 1>{mythen_detector->num_strips()}, 0);
+    }
+
+    void read_flatfield_from_file(const std::string &filename) {
+
+        std::string word;
+        uint32_t strip_number{};
+
+        try {
+            std::ifstream file(filename, std::ios_base::in);
+            if (!file.good()) {
+                throw std::logic_error("file does not exist");
+            }
+
+            std::stringstream file_buffer;
+            file_buffer << file.rdbuf();
+
+            while (file_buffer >> word) {
+
+                strip_number = std::stoi(word);
+
+                file_buffer >> word;
+                if (!mythen_detector->get_bad_channels()[strip_number])
+                    flat_field[strip_number] = std::stod(word);
+            }
+
+            file.close();
+        } catch (const std::exception &e) {
+            std::cerr << "Error: " << e.what() << std::endl;
+        }
+    }
+
+    NDView<uint32_t, 1> get_flatfield() const { return flat_field.view(); }
+
+    double calculate_mean(double tolerance = 0.001) const {
+        auto [sum, count] = std::accumulate(
+            flat_field.begin(), flat_field.end(),
+            std::make_pair<double, ssize_t>(0.0, 0),
+            [&tolerance](std::pair<double, ssize_t> acc, const auto &element) {
+                return element == 0 ? acc
+                                    : std::make_pair(acc.first + element,
+                                                     acc.second + 1);
+            });
+
+        return sum / count;
+    }
+
+    NDArray<double, 1>
+    inverse_normalized_flatfield(double tolerance = 0.001) const {
+        double mean = calculate_mean(tolerance);
+
+        NDArray<double, 1> inverse_normalized_flatfield(flat_field.shape());
+
+        for (ssize_t i = 0; i < flat_field.size(); ++i) {
+            inverse_normalized_flatfield[i] =
+                (flat_field[i] <= tolerance ? 0.0 : mean / flat_field[i]);
+            if (inverse_normalized_flatfield[i] < tolerance)
+                mythen_detector->get_bad_channels()[i] = true;
+        }
+
+        return inverse_normalized_flatfield; // TODO: better to have a copy in
+                                             // this context but unneccessary
+                                             // for angle calibration code
+        // maybe provide inplace and copy option
+        // maybe store as member variable access with view
+    }
+
+    NDArray<double, 1> normalized_flatfield(double tolerance = 0.001) const {
+        double mean = calculate_mean(tolerance);
+
+        NDArray<double, 1> normalized_flatfield(flat_field.shape());
+
+        for (ssize_t i = 0; i < flat_field.size(); ++i) {
+            normalized_flatfield[i] = (flat_field[i] == flat_field[i] / mean);
+            if (normalized_flatfield[i] < tolerance)
+                mythen_detector->get_bad_channels()[i] = true;
+        }
+        return normalized_flatfield;
+    }
+
+  private:
+    NDArray<uint32_t, 1> flat_field; // TODO: should be 2d
+    std::shared_ptr<MythenDetectorSpecifications> mythen_detector;
+};
+} // namespace aare

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/Hdf5FileReader.hpp

@@ -0,0 +1,212 @@
+/************************************************
+ * @file HD5FFileReader.hpp
+ * @short HDF5FileReader based on H5File object
+ ***********************************************/
+
+#pragma once
+
+#include "Frame.hpp"
+#include "NDArray.hpp"
+#include <H5Cpp.h>
+#include <array>
+#include <cxxabi.h>
+#include <optional>
+
+namespace aare {
+
+// return std::type_info
+inline const std::type_info &deduce_cpp_type(const H5::DataType datatype) {
+    if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_UINT8.getId())) {
+        return typeid(uint8_t);
+    } else if (H5Tequal(datatype.getId(),
+                        H5::PredType::NATIVE_UINT16.getId())) {
+        return typeid(uint16_t);
+    } else if (H5Tequal(datatype.getId(),
+                        H5::PredType::NATIVE_UINT32.getId())) {
+        return typeid(uint32_t);
+    } else if (H5Tequal(datatype.getId(),
+                        H5::PredType::NATIVE_UINT64.getId())) {
+        return typeid(uint64_t);
+    } else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT8.getId())) {
+        return typeid(int8_t);
+    } else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT16.getId())) {
+        return typeid(int16_t);
+    } else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT32.getId())) {
+        return typeid(int32_t);
+    } else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT64.getId())) {
+        return typeid(int64_t);
+    } else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT.getId())) {
+        return typeid(int);
+    } else if (H5Tequal(datatype.getId(), H5::PredType::IEEE_F64LE.getId())) {
+        return typeid(double);
+    } else if (H5Tequal(datatype.getId(), H5::PredType::IEEE_F32LE.getId())) {
+        return typeid(float);
+    } else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_FLOAT.getId())) {
+        return typeid(float);
+    } else if (H5Tequal(datatype.getId(),
+                        H5::PredType::NATIVE_DOUBLE.getId())) {
+        return typeid(float);
+    } else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_CHAR.getId()) &&
+               datatype.getId() == H5::PredType::NATIVE_CHAR.getId()) {
+        return typeid(char);
+    } else {
+        throw std::runtime_error("c++ type cannot be deduced");
+    }
+}
+
+struct Subset {
+    std::vector<hsize_t> shape;
+    std::vector<hsize_t> offset; // index where data subset should start
+};
+
+class HDF5Dataset {
+
+  public:
+    HDF5Dataset(const std::string &datasetname_, const H5::DataSet dataset_)
+        : datasetname(datasetname_), dataset(dataset_) {
+        datatype = dataset.getDataType();
+
+        cpp_type = &deduce_cpp_type(datatype);
+
+        dataspace = dataset.getSpace();
+        rank = dataspace.getSimpleExtentNdims(); // number of dimensions
+
+        shape.resize(rank);
+        dataspace.getSimpleExtentDims(shape.data(), nullptr);
+    }
+
+    hsize_t get_shape(ssize_t index) const { return shape[index]; }
+
+    std::vector<hsize_t> get_shape() const { return shape; }
+
+    H5::DataType get_datatype() const { return datatype; }
+
+    const std::type_info *get_cpp_type() const { return cpp_type; }
+
+    /**
+     * Reads subset of dataset into the buffer
+     * e.g. to read one 2d frame pass Subset({shape[1], shape[2]}, {frame_index,
+     * 0,0})
+     */
+    void
+    read_into_buffer(std::byte *buffer,
+                     std::optional<const Subset> subset = std::nullopt) const {
+
+        if (subset) {
+            // TODO treat scalar cases
+            if (static_cast<ssize_t>(subset->offset.size()) != rank) {
+                throw std::runtime_error("provide an offset for" +
+                                         std::to_string(rank) + "dimensions");
+            }
+            for (ssize_t i = 0; i < rank; ++i) {
+                hsize_t size =
+                    i < rank - static_cast<ssize_t>(subset->shape.size())
+                        ? 0
+                        : subset->shape[i - (rank - subset->shape.size())];
+                if ((size + subset->offset[i]) > shape[i]) {
+                    throw std::runtime_error(
+                        "subset is too large or offset is out of bounds");
+                }
+            }
+
+            H5::DataSpace memspace(static_cast<int>(subset->shape.size()),
+                                   subset->shape.data());
+
+            dataspace.selectHyperslab(H5S_SELECT_SET, subset->shape.data(),
+                                      subset->offset.data());
+            dataset.read(buffer, datatype, memspace, dataspace);
+        } else {
+            dataset.read(buffer, datatype);
+        }
+    }
+
+    Frame store_as_frame() const {
+        uint32_t rows{}, cols{};
+        if (rank == 1) {
+            rows = 1;
+            // TODO overflow
+            cols = static_cast<uint32_t>(shape[0]);
+        } else if (rank == 2) {
+            rows = static_cast<uint32_t>(shape[0]);
+            cols = static_cast<uint32_t>(shape[1]);
+        } else {
+            throw std::runtime_error("Frame only supports 2d images");
+        }
+
+        Frame frame(rows, cols, Dtype(*cpp_type));
+
+        read_into_buffer(frame.data());
+
+        return frame;
+    }
+
+    template <typename T, ssize_t NDim>
+    NDArray<T, NDim> store_as_ndarray() const {
+        if (NDim != rank) {
+            std::cout
+                << "Warning: dataset dimension and array dimension mismatch"
+                << std::endl; // TODO: replace with log - still valid if we
+                              // want subset
+        }
+        if (typeid(T) != *cpp_type) {
+            std::cout << "Warning: dataset and array type mismatch"
+                      << std::endl;
+        }
+        std::array<ssize_t, NDim> array_shape{};
+        std::transform(
+            shape.begin(), shape.end(), array_shape.begin(),
+            [](const auto dim) { return static_cast<ssize_t>(dim); });
+
+        aare::NDArray<T, NDim> dataset_array(array_shape);
+
+        read_into_buffer(reinterpret_cast<std::byte *>(dataset_array.data()));
+
+        return dataset_array;
+    }
+
+    // getMemDataSize()
+
+  private:
+    std::string datasetname{};
+    H5::DataSet dataset;
+    H5::DataSpace dataspace;
+    H5::DataType datatype;
+    const std::type_info *cpp_type;
+    ssize_t rank{};
+    std::vector<hsize_t> shape{};
+};
+
+class HDF5FileReader {
+
+  public:
+    HDF5FileReader() = default;
+
+    void open_file(const std::string &filename_) {
+        filename = filename_;
+        try {
+            file = H5::H5File(filename, H5F_ACC_RDONLY);
+        } catch (H5::Exception &e) {
+            std::cerr << "Error: " << e.getDetailMsg() << std::endl;
+        }
+    }
+
+    void close_file() { file.close(); }
+
+    HDF5Dataset get_dataset(const std::string &dataset_name) const {
+        H5::DataSet dataset;
+        try {
+            dataset = file.openDataSet(dataset_name);
+        } catch (H5::Exception &e) {
+            std::cerr << "Error: " << e.getDetailMsg() << std::endl;
+        }
+
+        // TODO use optional to handle error
+        return HDF5Dataset(dataset_name, dataset);
+    }
+
+  private:
+    std::string filename{};
+    H5::H5File file;
+};
+
+} // namespace aare

include/aare/MythenDetectorSpecifications.hpp

@@ -0,0 +1,150 @@
+#pragma once
+#include <cmath>
+#include <cstdint>
+
+#include <fstream>
+#include <iostream>
+#include <sstream>
+#include <string>
+
+#include "NDArray.hpp"
+
+namespace aare {
+
+class MythenDetectorSpecifications {
+
+  public:
+    // TODO: constructor that reads from a config file
+
+    MythenDetectorSpecifications() {
+        num_strips_ = max_modules_ * strips_per_module_;
+
+        num_connected_modules_ = max_modules_;
+
+        bad_channels =
+            NDArray<bool, 1>(std::array<ssize_t, 1>{num_strips_}, false);
+
+        connected_modules = NDArray<bool, 1>(
+            std::array<ssize_t, 1>{static_cast<ssize_t>(max_modules_)}, true);
+    }
+
+    MythenDetectorSpecifications(const size_t max_modules,
+                                 const double exposure_time,
+                                 const double bloffset)
+        : max_modules_(max_modules), exposure_time_(exposure_time),
+          bloffset_(bloffset) {
+        num_strips_ = max_modules_ * strips_per_module_;
+
+        num_connected_modules_ = max_modules_;
+
+        bad_channels =
+            NDArray<bool, 1>(std::array<ssize_t, 1>{num_strips_}, false);
+
+        connected_modules = NDArray<bool, 1>(
+            std::array<ssize_t, 1>{static_cast<ssize_t>(max_modules_)}, true);
+    }
+
+    void read_bad_channels_from_file(const std::string &filename) {
+        std::string line;
+
+        try {
+            std::ifstream file(filename, std::ios_base::in);
+            if (!file.good()) {
+                throw std::logic_error("file does not exist");
+            }
+
+            while (std::getline(file, line)) {
+                std::size_t pos = line.find("-");
+
+                if (pos == std::string::npos) {
+                    bad_channels(std::stoi(line)) = true;
+                } else {
+                    size_t line_size = line.size();
+                    for (int i = std::stoi(line.substr(0, pos));
+                         i <= std::stoi(line.substr(pos + 1, line_size - pos));
+                         ++i)
+                        bad_channels(i) = true;
+                }
+            }
+
+            file.close();
+        } catch (const std::exception &e) {
+            std::cerr << "Error: " << e.what() << std::endl;
+        }
+    }
+
+    void read_unconnected_modules_from_file(const std::string &filename) {
+        std::string line;
+
+        try {
+            std::ifstream file(filename, std::ios_base::in);
+            if (!file.good()) {
+                throw std::logic_error("file does not exist");
+            }
+
+            std::stringstream file_buffer;
+            file_buffer << file.rdbuf();
+
+            file_buffer >> line;
+            num_connected_modules_ -= std::stoi(line);
+
+            while (file_buffer >> line) {
+                size_t module = std::stoi(line);
+                connected_modules[module] = false;
+                for (size_t i = module * strips_per_module_;
+                     i < (module + 1) * strips_per_module_; ++i)
+                    bad_channels[i] = true;
+            }
+        } catch (const std::exception &e) {
+            std::cerr << "Error: " << e.what() << std::endl;
+        }
+    }
+
+    NDView<bool, 1> get_bad_channels() const { return bad_channels.view(); }
+
+    NDView<bool, 1> get_connected_modules() const {
+        return connected_modules.view();
+    }
+
+    static constexpr double pitch() { return pitch_; }
+
+    static constexpr size_t strips_per_module() { return strips_per_module_; }
+
+    size_t max_modules() const { return max_modules_; }
+
+    double exposure_time() const { return exposure_time_; }
+
+    double bloffset() const { return bloffset_; }
+
+    double dtt0() const { return dtt0_; }
+
+    static constexpr double min_angle() { return min_angle_; }
+
+    static constexpr double max_angle() { return max_angle_; }
+
+    ssize_t num_strips() const { return num_strips_; }
+
+  private:
+    static constexpr size_t strips_per_module_ = 1280;
+    static constexpr double pitch_ = 0.05; // strip width [mm]
+    static constexpr double min_angle_ =
+        -180.0; // maybe shoudnt be static but configurable
+    static constexpr double max_angle_ = 180.0;
+    static constexpr double dtt0_ =
+        0.0; // No idea what this is - probably configurable
+
+    size_t max_modules_ = 48;
+
+    double exposure_time_ = 5.0; // TODO: could read from acquired file but
+                                 // maybe should be configurable
+    double bloffset_ = 1.532; // what is this? detector offset relative to what?
+
+    size_t num_connected_modules_{};
+
+    ssize_t num_strips_{};
+
+    NDArray<bool, 1> bad_channels;
+    NDArray<bool, 1> connected_modules; // connected modules
+};
+
+} // namespace aare

include/aare/MythenFileReader.hpp

@@ -0,0 +1,82 @@
+/************************************************
+ * @file MythenFileReader.hpp
+ * @short minimal file reader to read mythen files
+ ***********************************************/
+
+#include <bitset>
+#include <filesystem>
+#include <string>
+
+#include "Hdf5FileReader.hpp"
+#include "NDArray.hpp"
+
+namespace aare {
+
+struct MythenFrame {
+    NDArray<uint32_t, 1> photon_counts;
+    double detector_angle{};
+    // double reference_intensity{}; not needed
+    std::array<uint8_t, 3> channel_mask{};
+};
+
+/** minimal version for a mythen file reader */
+class MythenFileReader : public HDF5FileReader {
+
+  public:
+    MythenFileReader(const std::filesystem::path &file_path_,
+                     const std::string &file_prefix_)
+        : m_base_path(file_path_), file_prefix(file_prefix_) {};
+
+    MythenFrame read_frame(ssize_t frame_index) {
+        // TODO not a good design fixed number of digits in file name for frame
+        // number -> pad with zeros
+        //  not even sure if files have the same name
+        std::string current_file_name =
+            m_base_path / (file_prefix + std::to_string(frame_index) + ".h5");
+
+        MythenFrame frame;
+        open_file(current_file_name);
+
+        auto dataset_photon_count =
+            get_dataset("/entry/instrument/detector/data");
+
+        frame.photon_counts =
+            dataset_photon_count.store_as_ndarray<uint32_t, 1>();
+
+        ++frame.photon_counts; // Why though?
+
+        auto dataset_detector_angle =
+            get_dataset("/entry/instrument/NDAttributes/DetectorAngle");
+
+        dataset_detector_angle.read_into_buffer(
+            reinterpret_cast<std::byte *>(&frame.detector_angle));
+
+        auto dataset_channel_number =
+            get_dataset("/entry/instrument/NDAttributes/CounterMask");
+
+        uint8_t channel_number;
+
+        dataset_channel_number.read_into_buffer(
+            reinterpret_cast<std::byte *>(&channel_number));
+
+        std::bitset<3> binary_channel_numbers(channel_number); // 1 0 0
+
+        // binary_channel_numbers.flip(); // TODO not sure where most
+        // significant
+        //  bit is ask Anna again
+
+        frame.channel_mask = std::array<uint8_t, 3>{binary_channel_numbers[0],
+                                                    binary_channel_numbers[1],
+                                                    binary_channel_numbers[2]};
+
+        close_file();
+
+        return frame;
+    }
+
+  private:
+    std::filesystem::path m_base_path{};
+    std::string file_prefix{};
+};
+
+} // namespace aare

include/aare/NDArray.hpp

@@ -21,11 +21,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,20 +41,19 @@ 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<>())),
           data_(new T[size_]) {}
 
-
     /**
      * @brief Construct a new NDArray object with a shape and value.
      *
      * @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);
     }
 
@@ -69,8 +67,8 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
         std::copy(v.begin(), v.end(), begin());
     }
 
-    template<size_t Size>
-    NDArray(const std::array<T, Size>& arr) : NDArray<T,1>({Size}) {
+    template <size_t Size>
+    NDArray(const std::array<T, Size> &arr) : NDArray<T, 1>({Size}) {
         std::copy(arr.begin(), arr.end(), begin());
     }
 
@@ -79,7 +77,6 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
         : shape_(other.shape_), strides_(c_strides<Ndim>(shape_)),
           size_(other.size_), data_(other.data_) {
         other.reset(); // TODO! is this necessary?
-
     }
 
     // Copy constructor
@@ -113,10 +110,10 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
     NDArray &operator-=(const NDArray &other);
     NDArray &operator*=(const NDArray &other);
 
-    //Write directly to the data array, or create a new one
-    template<size_t Size>
-    NDArray<T,1>& operator=(const std::array<T,Size> &other){
-        if(Size != size_){
+    // Write directly to the data array, or create a new one
+    template <size_t Size>
+    NDArray<T, 1> &operator=(const std::array<T, Size> &other) {
+        if (Size != size_) {
             delete[] data_;
             size_ = Size;
             data_ = new T[size_];
@@ -142,6 +139,9 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
 
     NDArray<bool, Ndim> operator>(const NDArray &other);
 
+    bool equals(const NDArray<T, Ndim> &other,
+                const T tolerance = std::numeric_limits<T>::epsilon()) const;
+
     bool operator==(const NDArray &other) const;
     bool operator!=(const NDArray &other) const;
 
@@ -157,11 +157,6 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
 
     NDArray &operator&=(const T & /*mask*/);
 
-  
-
-
-    
-
     void sqrt() {
         for (int i = 0; i < size_; ++i) {
             data_[i] = std::sqrt(data_[i]);
@@ -186,22 +181,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 +223,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 +237,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 +249,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 +261,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 +273,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 +292,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 +305,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,83 +317,81 @@ 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;
     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>
+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 +403,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 +412,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,25 +422,52 @@ 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);
-    f.write(img.buffer(), img.size() * sizeof(T));
+    f.write(reinterpret_cast<char *>(img.buffer()), img.size() * sizeof(T));
     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);
-    f.read(img.buffer(), img.size() * sizeof(T));
+    f.read(reinterpret_cast<char *>(img.buffer()), img.size() * sizeof(T));
     f.close();
     return img;
 }
 
+template <typename T, ssize_t Ndim = 1>
+NDArray<T, Ndim> load_non_binary_file(const std::string &filename,
+                                      const std::array<ssize_t, Ndim> shape) {
+    std::string word;
+    NDArray<T, Ndim> array(shape);
+    try {
+        std::ifstream file(filename, std::ios_base::in);
+        if (!file.good()) {
+            throw std::logic_error("file does not exist");
+        }
 
+        std::stringstream file_buffer;
+        file_buffer << file.rdbuf();
+
+        ssize_t counter = 0;
+        while (file_buffer >> word && counter < size) {
+            array[counter] = static_cast<T>(
+                std::stod(word)); // TODO change for different Types
+            ++counter;
+        }
+
+        file.close();
+    } catch (const std::exception &e) {
+        std::cerr << "Error: " << e.what() << std::endl;
+    }
+
+    return array;
+}
 
 } // namespace aare

include/aare/NDView.hpp

@@ -1,11 +1,12 @@
 #pragma once
-#include "aare/defs.hpp"
 #include "aare/ArrayExpr.hpp"
+#include "aare/defs.hpp"
 
 #include <algorithm>
 #include <array>
 #include <cassert>
 #include <cstdint>
+#include <fstream>
 #include <functional>
 #include <iomanip>
 #include <iostream>
@@ -14,10 +15,11 @@
 #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,62 +27,74 @@ 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<>())) {}
+          size_(std::accumulate(std::begin(shape), std::end(shape), 1,
+                                std::multiplies<>())) {}
 
-    // NDView(T *buffer, const std::vector<int64_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<>())) {}
+    // 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<>())) {}
 
-    template <typename... Ix> std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) {
+    template <typename... Ix>
+    std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) {
         return buffer_[element_offset(strides_, index...)];
     }
 
-    template <typename... Ix> std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) const {
+    template <typename... Ix>
+    std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) const {
         return buffer_[element_offset(strides_, index...)];
     }
 
     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_)
@@ -92,18 +106,37 @@ template <typename T, int64_t Ndim = 2> class NDView : public ArrayExpr<NDView<T
         return true;
     }
 
+    bool equals(const NDView<T, Ndim> &other, const T tolerance) const {
+        if (shape_ != other.shape_)
+            return false;
+
+        using SignedT = typename make_signed<T>::type;
+
+        for (uint32_t i = 0; i != size_; ++i)
+            if (std::abs(static_cast<SignedT>(buffer_[i]) -
+                         static_cast<SignedT>(other.buffer_[i])) > tolerance)
+                return false;
+
+        return true;
+    }
+
     NDView &operator+=(const T val) { return elemenwise(val, std::plus<T>()); }
     NDView &operator-=(const T val) { return elemenwise(val, std::minus<T>()); }
-    NDView &operator*=(const T val) { return elemenwise(val, std::multiplies<T>()); }
-    NDView &operator/=(const T val) { return elemenwise(val, std::divides<T>()); }
+    NDView &operator*=(const T val) {
+        return elemenwise(val, std::multiplies<T>());
+    }
+    NDView &operator/=(const T val) {
+        return elemenwise(val, std::divides<T>());
+    }
 
-    NDView &operator/=(const NDView &other) { return elemenwise(other, std::divides<T>()); }
+    NDView &operator/=(const NDView &other) {
+        return elemenwise(other, std::divides<T>());
+    }
 
-
-    template<size_t Size>
-    NDView& operator=(const std::array<T, Size> &arr) {
-        if(size() != static_cast<ssize_t>(arr.size()))
-            throw std::runtime_error(LOCATION + "Array and NDView size mismatch");
+    template <size_t Size> NDView &operator=(const std::array<T, Size> &arr) {
+        if (size() != static_cast<ssize_t>(arr.size()))
+            throw std::runtime_error(LOCATION +
+                                     "Array and NDView size mismatch");
         std::copy(arr.begin(), arr.end(), begin());
         return *this;
     }
@@ -136,31 +169,33 @@ 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) {
+    template <class BinaryOperation>
+    NDView &elemenwise(T val, BinaryOperation op) {
         for (uint64_t i = 0; i != size_; ++i) {
             buffer_[i] = op(buffer_[i], val);
         }
         return *this;
     }
-    template <class BinaryOperation> NDView &elemenwise(const NDView &other, BinaryOperation op) {
+    template <class BinaryOperation>
+    NDView &elemenwise(const NDView &other, BinaryOperation op) {
         for (uint64_t i = 0; i != size_; ++i) {
             buffer_[i] = op(buffer_[i], other.buffer_[i]);
         }
         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);
@@ -170,9 +205,8 @@ template <typename T, int64_t Ndim> void NDView<T, Ndim>::print_all() const {
     }
 }
 
-
-template <typename T, int64_t Ndim>
-std::ostream& operator <<(std::ostream& os, const NDView<T, Ndim>& arr){
+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) {
             os << std::setw(3);
@@ -183,10 +217,16 @@ std::ostream& operator <<(std::ostream& os, const NDView<T, Ndim>& arr){
     return os;
 }
 
+template <typename T> NDView<T, 1> make_view(std::vector<T> &vec) {
+    return NDView<T, 1>(vec.data(), {static_cast<ssize_t>(vec.size())});
+}
 
-template <typename T>
-NDView<T,1> make_view(std::vector<T>& vec){
-    return NDView<T,1>(vec.data(), {static_cast<int64_t>(vec.size())});
+template <typename T, ssize_t Ndim>
+void save(NDView<T, Ndim> img, const std::string &pathname) {
+    std::ofstream f;
+    f.open(pathname, std::ios::binary);
+    f.write(reinterpret_cast<char *>(img.data()), img.size() * sizeof(T));
+    f.close();
 }
 
 } // 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/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/defs.hpp

@@ -3,16 +3,16 @@
 #include "aare/Dtype.hpp"
 
 #include <array>
-#include <stdexcept>
 #include <cassert>
 #include <cstdint>
 #include <cstring>
+#include <stdexcept>
 #include <string>
 #include <string_view>
+#include <type_traits>
 #include <variant>
 #include <vector>
 
-
 /**
  * @brief LOCATION macro to get the current location in the code
  */
@@ -20,28 +20,24 @@
     std::string(__FILE__) + std::string(":") + std::to_string(__LINE__) +      \
         ":" + std::string(__func__) + ":"
 
-
-
 #ifdef AARE_CUSTOM_ASSERT
-#define AARE_ASSERT(expr)\
-    if (expr)\
-        {}\
-    else\
+#define AARE_ASSERT(expr)                                                      \
+    if (expr) {                                                                \
+    } else                                                                     \
         aare::assert_failed(LOCATION + " Assertion failed: " + #expr + "\n");
 #else
-#define AARE_ASSERT(cond)\
-    do { (void)sizeof(cond); } while(0)
+#define AARE_ASSERT(cond)                                                      \
+    do {                                                                       \
+        (void)sizeof(cond);                                                    \
+    } while (0)
 #endif
 
-
 namespace aare {
 
 inline constexpr size_t bits_per_byte = 8;
 
 void assert_failed(const std::string &msg);
 
-
-
 class DynamicCluster {
   public:
     int cluster_sizeX;
@@ -55,7 +51,7 @@ class DynamicCluster {
 
   public:
     DynamicCluster(int cluster_sizeX_, int cluster_sizeY_,
-            Dtype dt_ = Dtype(typeid(int32_t)))
+                   Dtype dt_ = Dtype(typeid(int32_t)))
         : cluster_sizeX(cluster_sizeX_), cluster_sizeY(cluster_sizeY_),
           dt(dt_) {
         m_data = new std::byte[cluster_sizeX * cluster_sizeY * dt.bytes()]{};
@@ -179,24 +175,24 @@ template <typename T> struct t_xy {
 };
 using xy = t_xy<uint32_t>;
 
-
 /**
- * @brief Class to hold the geometry of a module. Where pixel 0 is located and the size of the module
+ * @brief Class to hold the geometry of a module. Where pixel 0 is located and
+ * the size of the module
  */
-struct ModuleGeometry{
+struct ModuleGeometry {
     int origin_x{};
     int origin_y{};
     int height{};
     int width{};
     int row_index{};
-    int col_index{}; 
+    int col_index{};
 };
 
 /**
- * @brief Class to hold the geometry of a detector. Number of modules, their size and where pixel 0 
- * for each module is located
+ * @brief Class to hold the geometry of a detector. Number of modules, their
+ * size and where pixel 0 for each module is located
  */
-struct DetectorGeometry{
+struct DetectorGeometry {
     int modules_x{};
     int modules_y{};
     int pixels_x{};
@@ -206,31 +202,30 @@ struct DetectorGeometry{
     std::vector<ModuleGeometry> module_pixel_0;
 };
 
-struct ROI{
-    int64_t xmin{};
-    int64_t xmax{};
-    int64_t ymin{};
-    int64_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 {
+struct ROI {
+    ssize_t xmin{};
+    ssize_t xmax{};
+    ssize_t ymin{};
+    ssize_t ymax{};
+
+    ssize_t height() const { return ymax - ymin; }
+    ssize_t width() const { return xmax - xmin; }
+    bool contains(ssize_t x, ssize_t y) const {
         return x >= xmin && x < xmax && y >= ymin && y < ymax;
     }
-  };
+};
 
+using dynamic_shape = std::vector<ssize_t>;
 
-using dynamic_shape = std::vector<int64_t>;
-
-//TODO! Can we uniform enums between the libraries?
+// TODO! Can we uniform enums between the libraries?
 
 /**
- * @brief Enum class to identify different detectors. 
+ * @brief Enum class to identify different detectors.
  * The values are the same as in slsDetectorPackage
  * Different spelling to avoid confusion with the slsDetectorPackage
  */
 enum class DetectorType {
-    //Standard detectors match the enum values from slsDetectorPackage
+    // Standard detectors match the enum values from slsDetectorPackage
     Generic,
     Eiger,
     Gotthard,
@@ -241,8 +236,9 @@ enum class DetectorType {
     Gotthard2,
     Xilinx_ChipTestBoard,
 
-    //Additional detectors used for defining processing. Variants of the standard ones.
-    Moench03=100,
+    // Additional detectors used for defining processing. Variants of the
+    // standard ones.
+    Moench03 = 100,
     Moench03_old,
     Unknown
 };
@@ -263,4 +259,12 @@ template <> FrameDiscardPolicy StringTo(const std::string & /*mode*/);
 
 using DataTypeVariants = std::variant<uint16_t, uint32_t>;
 
+template <typename T, bool = std::is_integral_v<T>> struct make_signed {
+    using type = T;
+};
+
+template <typename T> struct make_signed<T, true> {
+    using type = std::make_signed_t<T>;
+};
+
 } // namespace aare

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/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) {

src/AngleCalibration.cpp

@@ -0,0 +1,377 @@
+#include "aare/AngleCalibration.hpp"
+
+namespace aare {
+
+AngleCalibration::AngleCalibration(
+    std::shared_ptr<MythenDetectorSpecifications> mythen_detector_,
+    std::shared_ptr<FlatField> flat_field_,
+    std::shared_ptr<MythenFileReader> mythen_file_reader_)
+    : mythen_detector(mythen_detector_), flat_field(flat_field_),
+      mythen_file_reader(mythen_file_reader_) {
+    centers.reserve(mythen_detector->max_modules());
+    conversions.reserve(mythen_detector->max_modules());
+    offsets.reserve(mythen_detector->max_modules());
+
+    num_bins = std::floor(mythen_detector->max_angle() / histogram_bin_width) -
+               std::floor(mythen_detector->min_angle() / histogram_bin_width) +
+               1; // TODO only works if negative
+                  // and positive angle
+}
+
+void AngleCalibration::set_histogram_bin_width(double bin_width) {
+    histogram_bin_width = bin_width;
+
+    num_bins = std::floor(mythen_detector->max_angle() / histogram_bin_width) -
+               std::floor(mythen_detector->min_angle() / histogram_bin_width) +
+               1; // TODO only works if negative
+                  // and positive angle
+}
+
+double AngleCalibration::get_histogram_bin_width() const {
+    return histogram_bin_width;
+}
+
+ssize_t AngleCalibration::get_new_num_bins() const { return num_bins; }
+
+std::vector<double> AngleCalibration::get_centers() const { return centers; }
+
+std::vector<double> AngleCalibration::get_conversions() const {
+    return conversions;
+}
+
+std::vector<double> AngleCalibration::get_offsets() const { return offsets; }
+
+NDView<double, 1> AngleCalibration::get_new_photon_counts() const {
+    return new_photon_counts.view();
+}
+
+NDView<double, 1> AngleCalibration::get_new_statistical_errors() const {
+    return new_photon_count_errors.view();
+}
+
+void AngleCalibration::read_initial_calibration_from_file(
+    const std::string &filename) {
+
+    std::string line;
+    uint32_t module_number{};
+
+    try {
+        std::ifstream file(filename, std::ios_base::in);
+        if (!file.good()) {
+            throw std::logic_error("file does not exist");
+        }
+
+        std::stringstream file_buffer;
+        file_buffer << file.rdbuf();
+
+        while (file_buffer >> line) {
+            if (line == "module") {
+                file_buffer >> line;
+                module_number = std::stoi(line);
+            }
+            if (line == "center") {
+                file_buffer >> line;
+                centers.insert(centers.begin() + module_number,
+                               std::stod(line));
+            }
+            if (line == "conversion") {
+                file_buffer >> line;
+                conversions.insert(conversions.begin() + module_number,
+                                   std::stod(line));
+            }
+            if (line == "offset") {
+                file_buffer >> line;
+                offsets.insert(offsets.begin() + module_number,
+                               std::stod(line));
+            }
+        }
+
+        file.close();
+    } catch (const std::exception &e) {
+        std::cerr << "Error: " << e.what() << std::endl;
+    }
+}
+
+parameters AngleCalibration::convert_to_EE_parameters() const {
+
+    // normal distance between sample and detector (R)
+    std::vector<double> normal_distances(centers.size());
+    // distances between intersection point of sample normal and module origin
+    // (D)
+    std::vector<double> module_center_distances(centers.size());
+    // angles between undiffracted beam and orthogonal sample projection on
+    // detector (phi)
+    std::vector<double> angles(centers.size());
+
+    for (size_t i = 0; i < centers.size(); ++i) {
+        auto [module_center_distance, normal_distance, angle] =
+            convert_to_EE_parameters(i);
+        normal_distances[i] = normal_distance;
+        module_center_distances[i] = module_center_distance;
+        angles[i] = angle;
+    }
+
+    return std::make_tuple(module_center_distances, normal_distances, angles);
+}
+
+std::tuple<double, double, double>
+AngleCalibration::convert_to_EE_parameters(const size_t module_index) const {
+    return convert_to_EE_parameters(centers[module_index],
+                                    conversions[module_index],
+                                    offsets[module_index]);
+}
+
+std::tuple<double, double, double> AngleCalibration::convert_to_EE_parameters(
+    const double center, const double conversion, const double offset) const {
+    const double module_center_distance =
+        center * MythenDetectorSpecifications::pitch();
+    const double normal_distance =
+        MythenDetectorSpecifications::pitch() / std::abs(conversion);
+    const double angle = offset + 180.0 / M_PI * center * std::abs(conversion);
+
+    return std::make_tuple(module_center_distance, normal_distance, angle);
+}
+
+size_t AngleCalibration::global_to_local_strip_index_conversion(
+    const size_t global_strip_index) const {
+    const size_t module_index =
+        global_strip_index / MythenDetectorSpecifications::strips_per_module();
+    // local strip index in module
+    size_t local_strip_index =
+        global_strip_index -
+        module_index * MythenDetectorSpecifications::strips_per_module();
+    // switch if indexing is in clock-wise direction
+    local_strip_index =
+        std::signbit(conversions[module_index])
+            ? MythenDetectorSpecifications::strips_per_module() - 1 -
+                  local_strip_index
+            : local_strip_index;
+
+    return local_strip_index;
+}
+
+/*
+parameters
+AngleCalibration::convert_to_BC_parameters() {}
+*/
+
+double AngleCalibration::diffraction_angle_from_DG_parameters(
+    const double center, const double conversion, const double offset,
+    const size_t strip_index, const double distance_to_strip) const {
+
+    return offset + 180.0 / M_PI *
+                        (center * std::abs(conversion) -
+                         atan((center - (strip_index + distance_to_strip)) *
+                              std::abs(conversion)));
+}
+
+double AngleCalibration::diffraction_angle_from_EE_parameters(
+    const double module_center_distance, const double normal_distance,
+    const double angle, const size_t strip_index,
+    const double distance_to_strip) const {
+
+    return angle - 180.0 / M_PI *
+                       atan((module_center_distance -
+                             MythenDetectorSpecifications::pitch() *
+                                 (strip_index + distance_to_strip)) /
+                            normal_distance); // TODO: why is it minus
+                                              // is it defined counter
+                                              // clockwise? thought
+                                              // should have a flipped
+                                              // sign
+}
+
+double AngleCalibration::angular_strip_width_from_DG_parameters(
+    const double center, const double conversion, const double offset,
+    const size_t local_strip_index) const {
+
+    return std::abs(diffraction_angle_from_DG_parameters(
+                        center, conversion, offset, local_strip_index, -0.5) -
+                    diffraction_angle_from_DG_parameters(
+                        center, conversion, offset, local_strip_index, 0.5));
+}
+
+double AngleCalibration::angular_strip_width_from_EE_parameters(
+    const double module_center_distance, const double normal_distance,
+    const double angle, const size_t local_strip_index) const {
+
+    return std::abs(diffraction_angle_from_EE_parameters(
+                        module_center_distance, normal_distance, angle,
+                        local_strip_index, -0.5) -
+                    diffraction_angle_from_EE_parameters(
+                        module_center_distance, normal_distance, angle,
+                        local_strip_index, 0.5));
+
+    // TODO: again not sure about division order - taking abs anyway
+}
+
+void AngleCalibration::calculate_fixed_bin_angle_width_histogram(
+    const size_t start_frame_index, const size_t end_frame_index) {
+
+    new_photon_counts = NDArray<double, 1>(std::array<ssize_t, 1>{num_bins});
+
+    new_photon_count_errors =
+        NDArray<double, 1>(std::array<ssize_t, 1>{num_bins});
+
+    // TODO: maybe group these into a 2d array - better cache usage
+    NDArray<double, 1> bin_counts(std::array<ssize_t, 1>{num_bins}, 0.0);
+    NDArray<double, 1> new_statistical_weights(std::array<ssize_t, 1>{num_bins},
+                                               0.0);
+    NDArray<double, 1> new_errors(std::array<ssize_t, 1>{num_bins}, 0.0);
+
+    NDArray<double, 1> inverse_normalized_flatfield =
+        flat_field->inverse_normalized_flatfield();
+
+    for (size_t frame_index = start_frame_index; frame_index < end_frame_index;
+         ++frame_index) {
+        MythenFrame frame = mythen_file_reader->read_frame(frame_index);
+        redistribute_photon_counts_to_fixed_angle_bins(
+            frame, bin_counts.view(), new_statistical_weights.view(),
+            new_errors.view(), inverse_normalized_flatfield.view());
+    }
+
+    for (ssize_t i = 0; i < new_photon_counts.size(); ++i) {
+        new_photon_counts[i] = (new_statistical_weights[i] <=
+                                std::numeric_limits<double>::epsilon())
+                                   ? 0.
+                                   : bin_counts[i] / new_statistical_weights[i];
+        new_photon_count_errors[i] =
+            (bin_counts[i] <= std::numeric_limits<double>::epsilon())
+                ? 0.
+                : 1.0 / std::sqrt(bin_counts[i]);
+    }
+}
+
+void AngleCalibration::redistribute_photon_counts_to_fixed_angle_bins(
+    const MythenFrame &frame, NDView<double, 1> bin_counts,
+    NDView<double, 1> new_statistical_weights, NDView<double, 1> new_errors,
+    NDView<double, 1> inverse_normalized_flatfield) const {
+
+    ssize_t channel = 0; // TODO handle mask - FlatField still 1d
+
+    if (frame.photon_counts.shape()[0] != mythen_detector->num_strips()) {
+        throw std::runtime_error("wrong number of strips read");
+    }
+
+    ssize_t num_bins1 = mythen_detector->min_angle() / histogram_bin_width;
+    ssize_t num_bins2 = mythen_detector->max_angle() / histogram_bin_width;
+
+    std::cout << "position: " << frame.detector_angle
+              << std::endl; // replace with log
+
+    double exposure_rate = 1. / mythen_detector->exposure_time();
+
+    for (ssize_t strip_index = 0; strip_index < mythen_detector->num_strips();
+         ++strip_index) {
+
+        size_t module_index =
+            strip_index / MythenDetectorSpecifications::strips_per_module();
+
+        if (mythen_detector->get_bad_channels()[strip_index] ||
+            !mythen_detector->get_connected_modules()[module_index])
+            continue;
+
+        double poisson_error = std::sqrt(frame.photon_counts(strip_index)) *
+                               inverse_normalized_flatfield(strip_index) *
+                               exposure_rate;
+
+        double corrected_photon_count =
+            frame.photon_counts(strip_index) *
+            inverse_normalized_flatfield(strip_index) * exposure_rate;
+
+        size_t local_strip_index =
+            global_to_local_strip_index_conversion(strip_index);
+
+        double diffraction_angle = diffraction_angle_from_DG_parameters(
+            centers[module_index], conversions[module_index],
+            offsets[module_index], local_strip_index);
+
+        diffraction_angle += (frame.detector_angle + mythen_detector->dtt0() +
+                              mythen_detector->bloffset());
+
+        if (diffraction_angle < mythen_detector->min_angle() ||
+            diffraction_angle > mythen_detector->max_angle())
+            continue;
+
+        double angle_covered_by_strip = angular_strip_width_from_DG_parameters(
+            centers[module_index], conversions[module_index],
+            offsets[module_index], local_strip_index);
+
+        double photon_count_per_bin = histogram_bin_width *
+                                      corrected_photon_count /
+                                      angle_covered_by_strip;
+        double error_photon_count_per_bin =
+            histogram_bin_width * poisson_error / angle_covered_by_strip;
+
+        double statistical_weights =
+            1.0 / std::pow(error_photon_count_per_bin, 2); // 1./sigma²
+
+        double strip_boundary_left =
+            diffraction_angle - 0.5 * angle_covered_by_strip;
+        double strip_boundary_right =
+            diffraction_angle + 0.5 * angle_covered_by_strip;
+
+        ssize_t left_bin_index = std::max(
+            num_bins1,
+            static_cast<ssize_t>(
+                std::floor(strip_boundary_left / histogram_bin_width) - 1));
+        ssize_t right_bin_index = std::min(
+            num_bins2,
+            static_cast<ssize_t>(
+                std::ceil(strip_boundary_right / histogram_bin_width) + 1));
+
+        // TODO should it be < or <=
+        for (ssize_t bin = left_bin_index; bin <= right_bin_index; ++bin) {
+            double bin_coverage = std::min(strip_boundary_right,
+                                           (bin + 0.5) * histogram_bin_width) -
+                                  std::max(strip_boundary_left,
+                                           (bin - 0.5) * histogram_bin_width);
+
+            double bin_coverage_factor = bin_coverage / histogram_bin_width;
+
+            ssize_t bin_index = bin - num_bins1;
+            // TODO: maybe have this threshold configurable
+            if (bin_coverage >= 0.0001) {
+                new_statistical_weights(bin_index) +=
+                    statistical_weights * bin_coverage_factor;
+                bin_counts(bin_index) += statistical_weights *
+                                         bin_coverage_factor *
+                                         photon_count_per_bin;
+                new_errors(bin_index) += statistical_weights *
+                                         bin_coverage_factor *
+                                         std::pow(photon_count_per_bin, 2);
+            }
+        }
+    }
+}
+
+void AngleCalibration::write_to_file(
+    const std::string &filename, const bool store_nonzero_bins,
+    const std::filesystem::path &filepath) const {
+
+    std::ofstream output_file(filepath / filename);
+
+    if (!output_file) {
+        std::cerr << "Error opening file!"
+                  << std::endl; // TODO: replace with log
+    }
+
+    output_file << std::fixed << std::setprecision(15);
+
+    for (ssize_t i = 0; i < num_bins; ++i) {
+        if (new_photon_counts[i] <= std::numeric_limits<double>::epsilon() &&
+            store_nonzero_bins) {
+            continue;
+        }
+
+        output_file << std::floor(mythen_detector->min_angle() /
+                                  histogram_bin_width) *
+                               histogram_bin_width +
+                           i * histogram_bin_width
+                    << " " << new_photon_counts[i] << " "
+                    << new_photon_count_errors[i] << std::endl;
+    }
+    output_file.close();
+}
+
+} // namespace aare

src/AngleCalibration.test.cpp

@@ -0,0 +1,234 @@
+/************************************************
+ * @file AngleCalibration.test.cpp
+ * @short test case for angle calibration class
+ ***********************************************/
+
+#include "aare/AngleCalibration.hpp"
+
+#include <filesystem>
+
+#include "test_config.hpp"
+
+#include <iomanip>
+#include <type_traits>
+
+#include <catch2/catch_all.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+
+using namespace aare;
+
+TEST_CASE("read initial angle calibration file",
+          "[.anglecalibration] [.files]") {
+
+    std::shared_ptr<MythenDetectorSpecifications> mythen_detector_ptr =
+        std::make_shared<MythenDetectorSpecifications>();
+
+    AngleCalibration anglecalibration(mythen_detector_ptr,
+                                      std::shared_ptr<FlatField>{},
+                                      std::shared_ptr<MythenFileReader>{});
+
+    std::string filename = test_data_path() / "AngleCalibration_Test_Data" /
+                           "Angcal_2E_Feb2023_P29.off";
+
+    REQUIRE(std::filesystem::exists(filename));
+
+    anglecalibration.read_initial_calibration_from_file(filename);
+
+    auto centers = anglecalibration.get_centers();
+    auto conversions = anglecalibration.get_conversions();
+    auto offsets = anglecalibration.get_offsets();
+
+    std::cout.precision(17);
+
+    CHECK(centers.size() == 48);
+    CHECK(conversions.size() == 48);
+    CHECK(offsets.size() == 48);
+
+    CHECK(centers[9] == Catch::Approx(660.342326));
+    CHECK(offsets[47] == Catch::Approx(5.8053312));
+    CHECK(conversions[27] == Catch::Approx(-0.6581179125e-4));
+}
+
+TEST_CASE("read bad channels",
+          "[.anglecalibration][.mythenspecifications][.files]") {
+
+    MythenDetectorSpecifications mythen_detector;
+
+    std::string bad_channels_filename = test_data_path() /
+                                        "AngleCalibration_Test_Data" /
+                                        "bc2023_003_RING.chans";
+
+    REQUIRE(std::filesystem::exists(bad_channels_filename));
+
+    mythen_detector.read_bad_channels_from_file(bad_channels_filename);
+
+    CHECK(mythen_detector.get_bad_channels().size() == 61440);
+
+    CHECK(mythen_detector.get_bad_channels()[61437] == true);
+    CHECK(std::all_of(mythen_detector.get_bad_channels().begin() + 30720,
+                      mythen_detector.get_bad_channels().begin() + 61439,
+                      [](const bool element) { return element; }));
+}
+
+TEST_CASE("read unconnected modules",
+          "[.anglecalibration][.mythenspecifications][.files]") {
+
+    MythenDetectorSpecifications mythen_detector;
+
+    std::string unconnected_modules_filename =
+        test_data_path() / "AngleCalibration_Test_Data" / "ModOut.txt";
+
+    REQUIRE(std::filesystem::exists(unconnected_modules_filename));
+
+    mythen_detector.read_unconnected_modules_from_file(
+        unconnected_modules_filename);
+
+    CHECK(mythen_detector.get_connected_modules().size() == 48);
+
+    CHECK(std::all_of(mythen_detector.get_connected_modules().begin(),
+                      mythen_detector.get_connected_modules().end(),
+                      [](const bool element) { return element; }));
+}
+
+TEST_CASE("read flatfield", "[.anglecalibration][.flatfield][.files]") {
+
+    std::shared_ptr<MythenDetectorSpecifications> mythen_detector_ptr =
+        std::make_shared<MythenDetectorSpecifications>();
+
+    FlatField flatfield(mythen_detector_ptr);
+
+    std::string flatfield_filename =
+        test_data_path() / "AngleCalibration_Test_Data" /
+        "Flatfield_E22p0keV_T11000eV_up_48M_a_LONG_Feb2023_open_WS_SUMC.raw";
+
+    REQUIRE(std::filesystem::exists(flatfield_filename));
+
+    flatfield.read_flatfield_from_file(flatfield_filename);
+
+    auto flatfield_data = flatfield.get_flatfield();
+
+    CHECK(flatfield_data.size() == 61440);
+
+    CHECK(flatfield_data[0] == 0);
+    CHECK(flatfield_data[21] == 4234186);
+}
+
+TEST_CASE("compare result with python code", "[.anglecalibration] [.files]") {
+
+    auto fpath = test_data_path() / "AngleCalibration_Test_Data";
+
+    REQUIRE(std::filesystem::exists(fpath));
+
+    std::shared_ptr<MythenDetectorSpecifications> mythen_detector_ptr =
+        std::make_shared<MythenDetectorSpecifications>();
+
+    std::string bad_channels_filename = fpath / "bc2023_003_RING.chans";
+
+    REQUIRE(std::filesystem::exists(bad_channels_filename));
+
+    mythen_detector_ptr->read_bad_channels_from_file(bad_channels_filename);
+
+    std::string unconnected_modules_filename = fpath / "ModOut.txt";
+
+    REQUIRE(std::filesystem::exists(unconnected_modules_filename));
+
+    mythen_detector_ptr->read_unconnected_modules_from_file(
+        unconnected_modules_filename);
+
+    std::shared_ptr<FlatField> flat_field_ptr =
+        std::make_shared<FlatField>(mythen_detector_ptr);
+
+    std::string flatfield_filename =
+        fpath /
+        "Flatfield_E22p0keV_T11000eV_up_48M_a_LONG_Feb2023_open_WS_SUMC.raw";
+
+    REQUIRE(std::filesystem::exists(flatfield_filename));
+
+    flat_field_ptr->read_flatfield_from_file(flatfield_filename);
+
+    std::shared_ptr<MythenFileReader> mythen_file_reader_ptr =
+        std::make_shared<MythenFileReader>(fpath,
+                                           "ang1up_22keV_LaB60p3mm_48M_a_0");
+
+    AngleCalibration anglecalibration(mythen_detector_ptr, flat_field_ptr,
+                                      mythen_file_reader_ptr);
+
+    std::string initial_angles_filename = fpath / "Angcal_2E_Feb2023_P29.off";
+
+    REQUIRE(std::filesystem::exists(initial_angles_filename));
+
+    anglecalibration.read_initial_calibration_from_file(
+        initial_angles_filename);
+
+    anglecalibration.calculate_fixed_bin_angle_width_histogram(320, 340);
+
+    // anglecalibration.write_to_file("cpp_new_photon_counts.xye");
+
+    auto expected_filename_photons =
+        test_data_path() / "AngleCalibration_Test_Data" / "new_photons.bin";
+
+    REQUIRE(std::filesystem::exists(expected_filename_photons));
+
+    auto expected_filename_errors =
+        test_data_path() / "AngleCalibration_Test_Data" / "new_errors.bin";
+
+    REQUIRE(std::filesystem::exists(expected_filename_errors));
+
+    ssize_t new_num_bins = anglecalibration.get_new_num_bins();
+
+    auto python_output_errors = load<double, 1>(
+        expected_filename_errors, std::array<ssize_t, 1>{new_num_bins});
+
+    auto python_output_photons = load<double, 1>(
+        expected_filename_photons, std::array<ssize_t, 1>{new_num_bins});
+
+    CHECK(anglecalibration.get_new_photon_counts().equals(
+        python_output_photons.view(),
+        1e-8)); // not sure about precision does not exactly match to all
+                // decimal digits
+
+    CHECK(anglecalibration.get_new_statistical_errors().equals(
+        python_output_errors.view(),
+        1e-8)); //
+}
+
+TEST_CASE("check conversion from DG to EE parameters", "[.anglecalibration]") {
+
+    std::shared_ptr<MythenDetectorSpecifications> mythen_detector_ptr =
+        std::make_shared<MythenDetectorSpecifications>();
+
+    AngleCalibration anglecalibration(mythen_detector_ptr,
+                                      std::shared_ptr<FlatField>{},
+                                      std::shared_ptr<MythenFileReader>{});
+
+    // DG test parameters
+    const double center = 642.197591224993;
+    const double conversion = 0.657694036246975e-4;
+    const double offset = 5.004892881251670;
+    const ssize_t local_strip_index = 1;
+
+    double diffraction_angle_DG_param =
+        anglecalibration.diffraction_angle_from_DG_parameters(
+            center, conversion, offset, local_strip_index);
+
+    auto [distance_center, normal_distance, angle] =
+        anglecalibration.convert_to_EE_parameters(center, conversion, offset);
+
+    double diffraction_angle_EE_param =
+        anglecalibration.diffraction_angle_from_EE_parameters(
+            distance_center, normal_distance, angle, local_strip_index);
+
+    CHECK(diffraction_angle_EE_param ==
+          Catch::Approx(diffraction_angle_DG_param));
+
+    double strip_width_DG_param =
+        anglecalibration.angular_strip_width_from_DG_parameters(
+            center, conversion, offset, local_strip_index);
+
+    double strip_width_EE_param =
+        anglecalibration.angular_strip_width_from_EE_parameters(
+            distance_center, normal_distance, angle, local_strip_index);
+
+    CHECK(strip_width_DG_param == Catch::Approx(strip_width_EE_param));
+}

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) {
@@ -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/Hdf5FileReader.test.cpp

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

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/MythenFileReader.test.cpp

@@ -0,0 +1,33 @@
+/************************************************
+ * @file MythenFileReader.test.cpp
+ * @short test case for angle calibration class
+ ***********************************************/
+
+#include "aare/MythenFileReader.hpp"
+
+#include <filesystem>
+
+#include "test_config.hpp"
+
+#include <catch2/catch_all.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+
+using namespace aare;
+
+TEST_CASE("test mythenfile_reader", "[.mythenfilereader][.files]") {
+
+    auto fpath = test_data_path() / "AngleCalibration_Test_Data";
+
+    REQUIRE(std::filesystem::exists(fpath));
+
+        MythenFileReader file_reader(fpath, "ang1up_22keV_LaB60p3mm_48M_a_0");
+
+    auto frame = file_reader.read_frame(320);
+
+    CHECK(frame.detector_angle == 0.99955);
+
+    CHECK(frame.channel_mask == std::array<uint8_t, 3>{0, 0, 1});
+
+    CHECK(frame.photon_counts.size() == 61440);
+}

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

@@ -34,7 +34,7 @@ RawFile::RawFile(const std::filesystem::path &fname, const std::string &mode)
     }
 }
 
-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 +52,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);
@@ -70,7 +70,7 @@ void RawFile::read_into(std::byte *image_buf, size_t n_frames, DetectorHeader *h
             header+=n_mod();
     }
 
-};
+}
 
 size_t RawFile::n_mod() const { return n_subfile_parts; }
 
@@ -94,9 +94,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; }
@@ -360,4 +360,4 @@ RawFile::~RawFile() {
 
 
 
-} // namespace aare
+} // namespace aare

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; }
@@ -417,4 +417,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/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)