some refactoring

CMakeLists.txt

@@ -374,11 +374,13 @@ 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/AngleCalibration.hpp

@@ -11,6 +11,8 @@
 #include <string>
 #include <vector>
 
+#include "FlatField.hpp"
+#include "MythenDetectorSpecifications.hpp"
 #include "MythenFileReader.hpp"
 #include "NDArray.hpp"
 
@@ -19,304 +21,62 @@ namespace aare {
 using parameters =
     std::tuple<std::vector<double>, std::vector<double>, std::vector<double>>;
 
-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() { return bad_channels.view(); }
-
-    NDView<bool, 1> get_connected_modules() { 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() { return max_modules_; }
-
-    double exposure_time() { return exposure_time_; }
-
-    double bloffset() { return bloffset_; }
-
-    double dtt0() { return dtt0_; }
-
-    static constexpr double min_angle() { return min_angle_; }
-
-    static constexpr double max_angle() { return max_angle_; }
-
-    ssize_t num_strips() { 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
-};
-
-// 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() { return flat_field.view(); }
-
-    double calculate_mean(double tolerance = 0.001) {
-        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) {
-        double mean = calculate_mean(tolerance);
-
-        NDArray<double, 1> inverse_normalized_flatfield(flat_field.shape());
-
-        /*
-        std::transform(flat_field.begin(), flat_field.end(),
-                       inverse_normalized_flatfield.begin(),
-                       [&mean](const auto element) {
-                           return element == 0 ? 0.0 : mean / element;
-                       });
-        */
-
-        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) {
-        double mean = calculate_mean(tolerance);
-
-        NDArray<double, 1> normalized_flatfield(flat_field.shape());
-
-        /*
-        std::transform(flat_field.begin(), flat_field.end(),
-                       normalized_flatfield.begin(),
-                       [&mean](const auto element) { return element / mean; });
-        */
-
-        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;
-};
-
 class AngleCalibration {
 
   public:
     AngleCalibration(
         std::shared_ptr<MythenDetectorSpecifications> mythen_detector_,
         std::shared_ptr<FlatField> flat_field_,
-        std::shared_ptr<MythenFileReader> mythen_file_reader_)
+        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());
-
-        exposure_rate = 1. / mythen_detector->exposure_time();
-
-        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
-    }
 
     /** set the histogram bin width [degrees] */
-    void set_histogram_bin_width(double bin_width) {
+    void set_histogram_bin_width(double bin_width);
-        histogram_bin_width = bin_width;
 
-        num_bins =
+    double get_histogram_bin_width() const;
-            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 get_histogram_bin_width() { return histogram_bin_width; }
+    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() { return centers; }
+    std::vector<double> get_centers() const;
+    std::vector<double> get_conversions() const;
 
-    std::vector<double> get_conversions() { return conversions; }
+    std::vector<double> get_offsets() const;
 
-    std::vector<double> get_offsets() { return offsets; }
+    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();
+    parameters convert_to_EE_parameters() const;
 
     std::tuple<double, double, double>
-    convert_to_EE_parameters(const size_t module_index);
+    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 double offset) const;
 
     /** converts DG parameters to BC parameters e.g. best computing
      * parameters */
-    parameters convert_to_BC_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)
@@ -325,7 +85,7 @@ class AngleCalibration {
     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 double distance_to_strip = 0) const;
 
     /** calculates diffraction angle from EE module parameters (used in
      * Beer's Law)
@@ -338,37 +98,23 @@ class AngleCalibration {
      */
     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 size_t strip_index, const double distance_to_strip = 0) const;
 
     /** calculated the strip width expressed as angle [degrees]
-     * @param strip_index gloabl strip index of detector
+     * @param strip_index local strip index of module
      */
-    double angular_strip_width_from_DG_parameters(const size_t strip_index);
+    double angular_strip_width_from_DG_parameters(
-
+        const double center, const double conversion, const double offset,
-    double angular_strip_width_from_DG_parameters(const double center,
+        const size_t local_strip_index) const;
-                                                  const double conversion,
-                                                  const double offset,
-                                                  const size_t strip_index);
-
-    double angular_strip_width_from_EE_parameters(const size_t strip_index);
 
     double angular_strip_width_from_EE_parameters(
         const double module_center_distance, const double normal_distance,
-        const double angle, const size_t strip_index);
+        const double angle, const size_t local_strip_index) const;
 
+  protected:
     /** converts global strip index to local strip index of that module */
-    size_t
+    size_t global_to_local_strip_index_conversion(
-    global_to_local_strip_index_conversion(const size_t global_strip_index);
+        const size_t global_strip_index) 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);
 
     /**
      * redistributes photon counts with of histogram using one bin per strip
@@ -381,11 +127,9 @@ class AngleCalibration {
     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);
+        NDView<double, 1> inverse_nromalized_flatfield) const;
 
-    void write_to_file(const std::string &filename);
+  private:
-
-  protected:
     // TODO: Design maybe have a struct with three vectors, store all three
     // sets of parameters as member variables
 
@@ -413,8 +157,6 @@ class AngleCalibration {
 
     ssize_t num_bins{};
 
-    double exposure_rate{};
-
     std::shared_ptr<MythenFileReader>
         mythen_file_reader; // TODO replace by FileInterface ptr
 };

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

@@ -75,21 +75,22 @@ class HDF5Dataset {
         dataspace.getSimpleExtentDims(shape.data(), nullptr);
     }
 
-    hsize_t get_shape(ssize_t index) { return shape[index]; }
+    hsize_t get_shape(ssize_t index) const { return shape[index]; }
 
-    std::vector<hsize_t> get_shape() { return shape; }
+    std::vector<hsize_t> get_shape() const { return shape; }
 
-    H5::DataType get_datatype() { return datatype; }
+    H5::DataType get_datatype() const { return datatype; }
 
-    const std::type_info *get_cpp_type() { return cpp_type; }
+    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,
+    void
-                          std::optional<const Subset> subset = std::nullopt) {
+    read_into_buffer(std::byte *buffer,
+                     std::optional<const Subset> subset = std::nullopt) const {
 
         if (subset) {
             // TODO treat scalar cases
@@ -119,7 +120,7 @@ class HDF5Dataset {
         }
     }
 
-    Frame store_as_frame() {
+    Frame store_as_frame() const {
         uint32_t rows{}, cols{};
         if (rank == 1) {
             rows = 1;
@@ -139,7 +140,8 @@ class HDF5Dataset {
         return frame;
     }
 
-    template <typename T, ssize_t NDim> NDArray<T, NDim> store_as_ndarray() {
+    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"
@@ -190,7 +192,7 @@ class HDF5FileReader {
 
     void close_file() { file.close(); }
 
-    HDF5Dataset get_dataset(const std::string &dataset_name) {
+    HDF5Dataset get_dataset(const std::string &dataset_name) const {
         H5::DataSet dataset;
         try {
             dataset = file.openDataSet(dataset_name);

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

@@ -43,7 +43,7 @@ class MythenFileReader : public HDF5FileReader {
         frame.photon_counts =
             dataset_photon_count.store_as_ndarray<uint32_t, 1>();
 
-        ++frame.photon_counts;
+        ++frame.photon_counts; // Why though?
 
         auto dataset_detector_angle =
             get_dataset("/entry/instrument/NDAttributes/DetectorAngle");

include/aare/NDArray.hpp

@@ -139,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;
 
@@ -438,4 +441,33 @@ NDArray<T, Ndim> load(const std::string &pathname,
     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>
@@ -17,7 +18,8 @@ namespace aare {
 template <ssize_t Ndim> using Shape = std::array<ssize_t, Ndim>;
 
 // TODO! fix mismatch between signed and unsigned
-template <ssize_t Ndim> Shape<Ndim> make_shape(const std::vector<size_t> &shape) {
+template <ssize_t Ndim>
+Shape<Ndim> make_shape(const std::vector<size_t> &shape) {
     if (shape.size() != Ndim)
         throw std::runtime_error("Shape size mismatch");
     Shape<Ndim> arr;
@@ -25,14 +27,18 @@ template <ssize_t Ndim> Shape<Ndim> make_shape(const std::vector<size_t> &shape)
     return arr;
 }
 
-template <ssize_t Dim = 0, typename Strides> ssize_t element_offset(const Strides & /*unused*/) { return 0; }
+template <ssize_t Dim = 0, typename Strides>
+ssize_t element_offset(const Strides & /*unused*/) {
+    return 0;
+}
 
 template <ssize_t Dim = 0, typename Strides, typename... Ix>
 ssize_t element_offset(const Strides &strides, ssize_t i, Ix... index) {
     return i * strides[Dim] + element_offset<Dim + 1>(strides, index...);
 }
 
-template <ssize_t Ndim> std::array<ssize_t, Ndim> c_strides(const std::array<ssize_t, Ndim> &shape) {
+template <ssize_t Ndim>
+std::array<ssize_t, Ndim> c_strides(const std::array<ssize_t, Ndim> &shape) {
     std::array<ssize_t, Ndim> strides{};
     std::fill(strides.begin(), strides.end(), 1);
     for (ssize_t i = Ndim - 1; i > 0; --i) {
@@ -41,14 +47,16 @@ template <ssize_t Ndim> std::array<ssize_t, Ndim> c_strides(const std::array<ssi
     return strides;
 }
 
-template <ssize_t Ndim> std::array<ssize_t, Ndim> make_array(const std::vector<ssize_t> &vec) {
+template <ssize_t Ndim>
+std::array<ssize_t, Ndim> make_array(const std::vector<ssize_t> &vec) {
     assert(vec.size() == Ndim);
     std::array<ssize_t, Ndim> arr{};
     std::copy_n(vec.begin(), Ndim, arr.begin());
     return arr;
 }
 
-template <typename T, ssize_t Ndim = 2> class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
+template <typename T, ssize_t Ndim = 2>
+class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
   public:
     NDView() = default;
     ~NDView() = default;
@@ -57,17 +65,23 @@ template <typename T, ssize_t Ndim = 2> class NDView : public ArrayExpr<NDView<T
 
     NDView(T *buffer, std::array<ssize_t, Ndim> shape)
         : buffer_(buffer), strides_(c_strides<Ndim>(shape)), shape_(shape),
-          size_(std::accumulate(std::begin(shape), std::end(shape), 1, std::multiplies<>())) {}
+          size_(std::accumulate(std::begin(shape), std::end(shape), 1,
+                                std::multiplies<>())) {}
 
     // NDView(T *buffer, const std::vector<ssize_t> &shape)
-    //     : buffer_(buffer), strides_(c_strides<Ndim>(make_array<Ndim>(shape))), shape_(make_array<Ndim>(shape)),
+    //     : buffer_(buffer),
-    //       size_(std::accumulate(std::begin(shape), std::end(shape), 1, std::multiplies<>())) {}
+    //     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...)];
     }
 
@@ -92,18 +106,37 @@ template <typename T, ssize_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) {
-    NDView &operator/=(const T val) { return elemenwise(val, std::divides<T>()); }
+        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) {
-    template<size_t Size>
+        if (size() != static_cast<ssize_t>(arr.size()))
-    NDView& operator=(const std::array<T, Size> &arr) {
+            throw std::runtime_error(LOCATION +
-        if(size() != static_cast<ssize_t>(arr.size()))
+                                     "Array and NDView size mismatch");
-            throw std::runtime_error(LOCATION + "Array and NDView size mismatch");
         std::copy(arr.begin(), arr.end(), begin());
         return *this;
     }
@@ -147,13 +180,15 @@ template <typename T, ssize_t Ndim = 2> class NDView : public ArrayExpr<NDView<T
     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]);
         }
@@ -170,9 +205,8 @@ template <typename T, ssize_t Ndim> void NDView<T, Ndim>::print_all() const {
     }
 }
 
-
 template <typename T, ssize_t Ndim>
-std::ostream& operator <<(std::ostream& os, const NDView<T, Ndim>& arr){
+std::ostream &operator<<(std::ostream &os, const NDView<T, Ndim> &arr) {
     for (auto row = 0; row < arr.shape(0); ++row) {
         for (auto col = 0; col < arr.shape(1); ++col) {
             os << std::setw(3);
@@ -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>
+template <typename T, ssize_t Ndim>
-NDView<T,1> make_view(std::vector<T>& vec){
+void save(NDView<T, Ndim> img, const std::string &pathname) {
-    return NDView<T,1>(vec.data(), {static_cast<ssize_t>(vec.size())});
+    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/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)\
+#define AARE_ASSERT(expr)                                                      \
-    if (expr)\
+    if (expr) {                                                                \
-        {}\
+    } else                                                                     \
-    else\
         aare::assert_failed(LOCATION + " Assertion failed: " + #expr + "\n");
 #else
-#define AARE_ASSERT(cond)\
+#define AARE_ASSERT(cond)                                                      \
-    do { (void)sizeof(cond); } while(0)
+    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,11 +175,11 @@ 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{};
@@ -193,10 +189,10 @@ struct ModuleGeometry{
 };
 
 /**
- * @brief Class to hold the geometry of a detector. Number of modules, their size and where pixel 0 
+ * @brief Class to hold the geometry of a detector. Number of modules, their
- * for each module is located
+ * size and where pixel 0 for each module is located
  */
-struct DetectorGeometry{
+struct DetectorGeometry {
     int modules_x{};
     int modules_y{};
     int pixels_x{};
@@ -206,7 +202,7 @@ struct DetectorGeometry{
     std::vector<ModuleGeometry> module_pixel_0;
 };
 
-struct ROI{
+struct ROI {
     ssize_t xmin{};
     ssize_t xmax{};
     ssize_t ymin{};
@@ -217,12 +213,11 @@ struct ROI{
     bool contains(ssize_t x, ssize_t y) const {
         return x >= xmin && x < xmax && y >= ymin && y < ymax;
     }
-  };
+};
-
 
 using dynamic_shape = std::vector<ssize_t>;
 
-//TODO! Can we uniform enums between the libraries?
+// TODO! Can we uniform enums between the libraries?
 
 /**
  * @brief Enum class to identify different detectors.
@@ -230,7 +225,7 @@ using dynamic_shape = std::vector<ssize_t>;
  * 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.
+    // Additional detectors used for defining processing. Variants of the
-    Moench03=100,
+    // 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

src/AngleCalibration.cpp

@@ -2,6 +2,53 @@
 
 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) {
 
@@ -45,7 +92,7 @@ void AngleCalibration::read_initial_calibration_from_file(
     }
 }
 
-parameters AngleCalibration::convert_to_EE_parameters() {
+parameters AngleCalibration::convert_to_EE_parameters() const {
 
     // normal distance between sample and detector (R)
     std::vector<double> normal_distances(centers.size());
@@ -68,14 +115,14 @@ parameters AngleCalibration::convert_to_EE_parameters() {
 }
 
 std::tuple<double, double, double>
-AngleCalibration::convert_to_EE_parameters(const size_t module_index) {
+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 double center, const double conversion, const double offset) const {
     const double module_center_distance =
         center * MythenDetectorSpecifications::pitch();
     const double normal_distance =
@@ -86,7 +133,7 @@ std::tuple<double, double, double> AngleCalibration::convert_to_EE_parameters(
 }
 
 size_t AngleCalibration::global_to_local_strip_index_conversion(
-    const size_t global_strip_index) {
+    const size_t global_strip_index) const {
     const size_t module_index =
         global_strip_index / MythenDetectorSpecifications::strips_per_module();
     // local strip index in module
@@ -110,7 +157,7 @@ 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 size_t strip_index, const double distance_to_strip) const {
 
     return offset + 180.0 / M_PI *
                         (center * std::abs(conversion) -
@@ -121,7 +168,7 @@ double AngleCalibration::diffraction_angle_from_DG_parameters(
 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 double distance_to_strip) const {
 
     return angle - 180.0 / M_PI *
                        atan((module_center_distance -
@@ -134,23 +181,9 @@ double AngleCalibration::diffraction_angle_from_EE_parameters(
                                               // sign
 }
 
-double AngleCalibration::angular_strip_width_from_DG_parameters(
-    const size_t strip_index) {
-
-    const size_t module_index =
-        strip_index / MythenDetectorSpecifications::strips_per_module();
-
-    return angular_strip_width_from_DG_parameters(
-        centers[module_index], conversions[module_index], offsets[module_index],
-        strip_index);
-}
-
 double AngleCalibration::angular_strip_width_from_DG_parameters(
     const double center, const double conversion, const double offset,
-    const size_t strip_index) {
+    const size_t local_strip_index) const {
-
-    const size_t local_strip_index =
-        global_to_local_strip_index_conversion(strip_index);
 
     return std::abs(diffraction_angle_from_DG_parameters(
                         center, conversion, offset, local_strip_index, -0.5) -
@@ -158,25 +191,9 @@ double AngleCalibration::angular_strip_width_from_DG_parameters(
                         center, conversion, offset, local_strip_index, 0.5));
 }
 
-double AngleCalibration::angular_strip_width_from_EE_parameters(
-    const size_t strip_index) {
-
-    const size_t module_index =
-        strip_index / MythenDetectorSpecifications::strips_per_module();
-
-    const auto [module_center_distance, normal_distance, angle] =
-        convert_to_EE_parameters(module_index);
-
-    return angular_strip_width_from_EE_parameters(
-        module_center_distance, normal_distance, angle, strip_index);
-}
-
 double AngleCalibration::angular_strip_width_from_EE_parameters(
     const double module_center_distance, const double normal_distance,
-    const double angle, const size_t strip_index) {
+    const double angle, const size_t local_strip_index) const {
-
-    const size_t local_strip_index =
-        global_to_local_strip_index_conversion(strip_index);
 
     return std::abs(diffraction_angle_from_EE_parameters(
                         module_center_distance, normal_distance, angle,
@@ -196,10 +213,10 @@ void AngleCalibration::calculate_fixed_bin_angle_width_histogram(
     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 =
@@ -228,7 +245,7 @@ void AngleCalibration::calculate_fixed_bin_angle_width_histogram(
 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) {
+    NDView<double, 1> inverse_normalized_flatfield) const {
 
     ssize_t channel = 0; // TODO handle mask - FlatField still 1d
 
@@ -242,6 +259,8 @@ void AngleCalibration::redistribute_photon_counts_to_fixed_angle_bins(
     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) {
 
@@ -254,7 +273,8 @@ void AngleCalibration::redistribute_photon_counts_to_fixed_angle_bins(
 
         double poisson_error = std::sqrt(frame.photon_counts(strip_index)) *
                                inverse_normalized_flatfield(strip_index) *
-                               exposure_rate; // not sure what this is
+                               exposure_rate;
+
         double corrected_photon_count =
             frame.photon_counts(strip_index) *
             inverse_normalized_flatfield(strip_index) * exposure_rate;
@@ -273,8 +293,9 @@ void AngleCalibration::redistribute_photon_counts_to_fixed_angle_bins(
             diffraction_angle > mythen_detector->max_angle())
             continue;
 
-        double angle_covered_by_strip =
+        double angle_covered_by_strip = angular_strip_width_from_DG_parameters(
-            angular_strip_width_from_DG_parameters(strip_index);
+            centers[module_index], conversions[module_index],
+            offsets[module_index], local_strip_index);
 
         double photon_count_per_bin = histogram_bin_width *
                                       corrected_photon_count /
@@ -324,8 +345,11 @@ void AngleCalibration::redistribute_photon_counts_to_fixed_angle_bins(
     }
 }
 
-void AngleCalibration::write_to_file(const std::string &filename) {
+void AngleCalibration::write_to_file(
-    std::ofstream output_file(filename);
+    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!"
@@ -335,7 +359,8 @@ void AngleCalibration::write_to_file(const std::string &filename) {
     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()) {
+        if (new_photon_counts[i] <= std::numeric_limits<double>::epsilon() &&
+            store_nonzero_bins) {
             continue;
         }
 

src/AngleCalibration.test.cpp

@@ -18,101 +18,18 @@
 
 using namespace aare;
 
-template <typename T, ssize_t Ndim = 1>
-NDArray<T, Ndim> read_into_array(const std::string &filename,
-                                 const std::array<ssize_t, Ndim> size) {
-    std::string word;
-    NDArray<T, Ndim> array(size);
-    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) {
-            array[counter] = 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;
-}
-
-template <typename T, bool = std::is_integral_v<T>> struct safe_make_signed {
-    using type = T;
-};
-
-template <typename T> struct safe_make_signed<T, true> {
-    using type = std::make_signed_t<T>;
-};
-
-template <typename T, ssize_t Ndim>
-bool check_equality_of_arrays(NDView<T, Ndim> array1, NDView<T, Ndim> array2) {
-    bool equal = true;
-
-    using SignedT = typename safe_make_signed<T>::type;
-
-    for (ssize_t i = 0; i < array1.size(); ++i) {
-        if (std::abs(static_cast<SignedT>(array1[i]) -
-                     static_cast<SignedT>(array2[i])) > 1e-6) {
-            std::cout << "index: " << i << std::endl;
-            std::cout << std::setprecision(15) << array1[i] << std::endl;
-            std::cout << std::setprecision(15) << array2[i] << std::endl;
-            equal = false;
-            break;
-        }
-    }
-
-    return equal;
-}
-
-class AngleCalibrationTestClass : public AngleCalibration {
-
-  public:
-    AngleCalibrationTestClass(
-        std::shared_ptr<MythenDetectorSpecifications> mythen_detector_,
-        std::shared_ptr<FlatField> flat_field_,
-        std::shared_ptr<MythenFileReader> mythen_file_reader_)
-        : AngleCalibration(mythen_detector_, flat_field_, mythen_file_reader_) {
-    }
-    ~AngleCalibrationTestClass() = default;
-
-    std::vector<double> get_centers() { return centers; }
-
-    std::vector<double> get_conversions() { return conversions; }
-
-    std::vector<double> get_offsets() { return offsets; }
-};
-
 TEST_CASE("read initial angle calibration file",
           "[.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::shared_ptr<FlatField> flat_field_ptr =
+    AngleCalibration anglecalibration(mythen_detector_ptr,
-        std::make_shared<FlatField>(mythen_detector_ptr);
+                                      std::shared_ptr<FlatField>{},
+                                      std::shared_ptr<MythenFileReader>{});
 
-    std::shared_ptr<MythenFileReader> mythen_file_reader_ptr =
+    std::string filename = test_data_path() / "AngleCalibration_Test_Data" /
-        std::make_shared<MythenFileReader>(fpath,
+                           "Angcal_2E_Feb2023_P29.off";
-                                           "ang1up_22keV_LaB60p3mm_48M_a_0");
-
-    AngleCalibrationTestClass anglecalibration(
-        mythen_detector_ptr, flat_field_ptr, mythen_file_reader_ptr);
-
-    std::string filename = fpath / "Angcal_2E_Feb2023_P29.off";
 
     REQUIRE(std::filesystem::exists(filename));
 
@@ -135,13 +52,12 @@ TEST_CASE("read initial angle calibration file",
 
 TEST_CASE("read bad channels",
           "[.anglecalibration][.mythenspecifications][.files]") {
-    auto fpath = test_data_path() / "AngleCalibration_Test_Data";
-
-    REQUIRE(std::filesystem::exists(fpath));
 
     MythenDetectorSpecifications mythen_detector;
 
-    std::string bad_channels_filename = fpath / "bc2023_003_RING.chans";
+    std::string bad_channels_filename = test_data_path() /
+                                        "AngleCalibration_Test_Data" /
+                                        "bc2023_003_RING.chans";
 
     REQUIRE(std::filesystem::exists(bad_channels_filename));
 
@@ -157,13 +73,11 @@ TEST_CASE("read bad channels",
 
 TEST_CASE("read unconnected modules",
           "[.anglecalibration][.mythenspecifications][.files]") {
-    auto fpath = test_data_path() / "AngleCalibration_Test_Data";
-
-    REQUIRE(std::filesystem::exists(fpath));
 
     MythenDetectorSpecifications mythen_detector;
 
-    std::string unconnected_modules_filename = fpath / "ModOut.txt";
+    std::string unconnected_modules_filename =
+        test_data_path() / "AngleCalibration_Test_Data" / "ModOut.txt";
 
     REQUIRE(std::filesystem::exists(unconnected_modules_filename));
 
@@ -178,9 +92,6 @@ TEST_CASE("read unconnected modules",
 }
 
 TEST_CASE("read flatfield", "[.anglecalibration][.flatfield][.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>();
@@ -188,7 +99,7 @@ TEST_CASE("read flatfield", "[.anglecalibration][.flatfield][.files]") {
     FlatField flatfield(mythen_detector_ptr);
 
     std::string flatfield_filename =
-        fpath /
+        test_data_path() / "AngleCalibration_Test_Data" /
         "Flatfield_E22p0keV_T11000eV_up_48M_a_LONG_Feb2023_open_WS_SUMC.raw";
 
     REQUIRE(std::filesystem::exists(flatfield_filename));
@@ -203,112 +114,7 @@ TEST_CASE("read flatfield", "[.anglecalibration][.flatfield][.files]") {
     CHECK(flatfield_data[21] == 4234186);
 }
 
-TEST_CASE("check flatfield values", "[.anglecalibration][.flatfield][.files]") {
+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);
-
-    FlatField 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));
-
-    flatfield.read_flatfield_from_file(flatfield_filename);
-
-    auto flatfield_data = flatfield.get_flatfield();
-
-    std::string expected_flatfield_filename = fpath / "flatfield.txt";
-
-    REQUIRE(std::filesystem::exists(expected_flatfield_filename));
-
-    NDArray<double, 1> expected_flatfield = read_into_array<double, 1>(
-        expected_flatfield_filename, flatfield_data.shape());
-
-    auto bad_channels = mythen_detector_ptr->get_bad_channels();
-
-    bool equal_flatfield = true;
-    for (ssize_t i = 0; i < flatfield_data.size(); ++i) {
-        if (!bad_channels[i] &&
-            std::abs(flatfield_data[i] - expected_flatfield[i]) >
-                std::numeric_limits<double>::epsilon()) {
-            std::cout << "index: " << i << std::endl;
-            std::cout << flatfield_data[i] << std::endl;
-            std::cout << expected_flatfield[i] << std::endl;
-            equal_flatfield = false;
-            break;
-        }
-    }
-    CHECK(equal_flatfield);
-}
-
-TEST_CASE("check inverse flatfield values",
-          "[.anglecalibration][.flatfield][.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);
-
-    FlatField 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));
-
-    flatfield.read_flatfield_from_file(flatfield_filename);
-
-    auto inverse_flatfield = flatfield.inverse_normalized_flatfield();
-
-    std::string expected_inverseflatfield_filename =
-        fpath / "inverseflatfield.txt";
-
-    REQUIRE(std::filesystem::exists(expected_inverseflatfield_filename));
-
-    NDArray<double, 1> expected_inverseflatfield = read_into_array<double, 1>(
-        expected_inverseflatfield_filename, inverse_flatfield.shape());
-
-    auto bad_channels = mythen_detector_ptr->get_bad_channels();
-
-    bool equal = true;
-    for (ssize_t i = 0; i < inverse_flatfield.size(); ++i) {
-        if (!bad_channels[i] &&
-            std::abs(inverse_flatfield[i] - expected_inverseflatfield[i]) >
-                1e-10) {
-            std::cout << "index: " << i << std::endl;
-            std::cout << std::setprecision(15) << inverse_flatfield[i]
-                      << std::endl;
-            std::cout << std::setprecision(15) << expected_inverseflatfield[i]
-                      << std::endl;
-            equal = false;
-            break;
-        }
-    }
-    CHECK(equal);
-}
-
-TEST_CASE("calculate new fixed angle width bins histogram",
-          "[.anglecalibration] [.files]") {
 
     auto fpath = test_data_path() / "AngleCalibration_Test_Data";
 
@@ -357,78 +163,72 @@ TEST_CASE("calculate new fixed angle width bins histogram",
 
     anglecalibration.calculate_fixed_bin_angle_width_histogram(320, 340);
 
-    anglecalibration.write_to_file(
+    // anglecalibration.write_to_file("cpp_new_photon_counts.xye");
-        "cpp_new_photon_counts.xye"); // TODO adjust output path
+
+    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("compare result with python code", "[.anglecalibration][.files]") {
+TEST_CASE("check conversion from DG to EE parameters", "[.anglecalibration]") {
-    auto expected_filename =
-        test_data_path() / "AngleCalibration_Test_Data" / "out2.xye";
-
-    REQUIRE(std::filesystem::exists(expected_filename));
-
-    auto expected_array = read_into_array<double, 2>(
-        expected_filename.string(), std::array<ssize_t, 2>{61440, 3});
-
-    auto filename =
-        std::filesystem::current_path() / "../build/cpp_new_photon_counts.xye";
-
-    // auto array = load<double, 2>(filename, std::array<ssize_t, 2>{61440, 3});
-    auto array = read_into_array<double, 2>(filename.string(),
-                                            std::array<ssize_t, 2>{61440, 3});
-
-    CHECK(check_equality_of_arrays(array.view(), expected_array.view()));
-}
-
-TEST_CASE("check conversion from DG to EE parameters",
-          "[.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::shared_ptr<FlatField> flat_field_ptr =
+    AngleCalibration anglecalibration(mythen_detector_ptr,
-        std::make_shared<FlatField>(mythen_detector_ptr);
+                                      std::shared_ptr<FlatField>{},
+                                      std::shared_ptr<MythenFileReader>{});
 
-    std::shared_ptr<MythenFileReader> mythen_file_reader_ptr =
+    // DG test parameters
-        std::make_shared<MythenFileReader>(fpath,
+    const double center = 642.197591224993;
-                                           "ang1up_22keV_LaB60p3mm_48M_a_0");
+    const double conversion = 0.657694036246975e-4;
-
+    const double offset = 5.004892881251670;
-    AngleCalibration anglecalibration(mythen_detector_ptr, flat_field_ptr,
+    const ssize_t local_strip_index = 1;
-                                      mythen_file_reader_ptr);
-
-    // DG parameters
-    double center = 642.197591224993;
-    double conversion = 0.657694036246975e-4;
-    double offset = 5.004892881251670;
-
-    double global_strip_index =
-        MythenDetectorSpecifications::strips_per_module() + 1;
 
     double diffraction_angle_DG_param =
         anglecalibration.diffraction_angle_from_DG_parameters(
-            center, conversion, offset, 1);
+            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, 1);
+            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, global_strip_index);
+            center, conversion, offset, local_strip_index);
 
     double strip_width_EE_param =
         anglecalibration.angular_strip_width_from_EE_parameters(
-            distance_center, normal_distance, angle, global_strip_index);
+            distance_center, normal_distance, angle, local_strip_index);
 
     CHECK(strip_width_DG_param == Catch::Approx(strip_width_EE_param));
 }