implementation of FlatField and MythenDetector class

include/aare/AngleCalibration.hpp

@@ -3,45 +3,232 @@
 #include <fstream>
 #include <iostream>
 #include <math.h>
+#include <memory>
 #include <sstream>
 #include <string>
 #include <vector>
 
-// function check connected that reads from a file if modules are connected or
+#include "NDArray.hpp"
-// not - which module though?
-
-// function read flatfield store in ff_corr probably correlation field
 
 namespace aare {
 
 using parameters =
     std::tuple<std::vector<double>, std::vector<double>, std::vector<double>>;
 
-struct MythenSpecifications {
+// TODO: some of these should be configurable - read them from a config file
+class MythenDetectorSpecifications {
 
-    static constexpr int32_t max_modules = 48;
+  public:
-    static constexpr int32_t strips_per_module = 1280;
+    // TODO: constructor that reads from a config file
-    static constexpr double pitch = 0.05; // strip width [mm] ?? TODO: not sure
+
-    static constexpr double ttmin = -180.0; // what is this?
+    MythenDetectorSpecifications() {
-    static constexpr float ttmax = 180.0;
+        num_strips_ = max_modules_ * strips_per_module_;
-    static constexpr float ttstep =
+
+        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>{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_; }
+
+    static constexpr size_t max_modules() { return max_modules_; }
+
+    static constexpr double exposure_time() { return exposure_time_; }
+
+    ssize_t num_strips() { return num_strips_; }
+
+  private:
+    static constexpr size_t max_modules_ = 48;
+    static constexpr size_t strips_per_module_ = 1280;
+    static constexpr double pitch_ = 0.05; // strip width [mm] ?? TODO: not sure
+    static constexpr double ttmin_ = -180.0; // what is this?
+    static constexpr float ttmax_ = 180.0;
+    static constexpr float ttstep_ =
         0.0036; // probably here to calculate bin size, what is this?
 
-    static constexpr double bloffset =
+    static constexpr double bloffset_ =
         1.532; // what is this? detector offset relative to what?
+
+    static constexpr double exposure_time_ =
+        5.0; // TODO: could read from acquired file but maybe should be
+             // configurable
+
+    static constexpr double dtt0_ = 0.0; // No idea what this is
+
+    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()});
+    }
+
+    void read_flatfield_from_file(const std::string &filename) {
+
+        std::string word;
+        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 >> word) {
+
+                module_number = std::stoi(word);
+
+                file_buffer >> word;
+                if (!mythen_detector->get_bad_channels()[module_number])
+                    flat_field[module_number] = std::stod(word);
+            }
+
+            file.close();
+        } catch (const std::exception &e) {
+            std::cerr << "Error: " << e.what() << std::endl;
+        }
+    }
+
+    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 < tolerance ? 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 mean / element; });
+
+        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; });
+
+        return normalized_flatfield;
+    }
+
+    // TODO: update is bad channels
+  private:
+    NDArray<uint32_t, 1> flat_field; // TODO: should be 2d
+    std::shared_ptr<MythenDetectorSpecifications> mythen_detector;
 };
 
-// number_of_activated_modules
-// number_of_channles
-// number_of_dimension
-// is bad array keeping track of bad channels!!
 class AngleCalibration {
 
   public:
-    AngleCalibration() {
+    AngleCalibration(
-        centers.reserve(MythenSpecifications::max_modules);
+        std::shared_ptr<MythenDetectorSpecifications> mythen_detector_,
-        conversions.reserve(MythenSpecifications::max_modules);
+        std::shared_ptr<FlatField> flat_field_)
-        offsets.reserve(MythenSpecifications::max_modules);
+        : mythen_detector(mythen_detector_), flat_field(flat_field_) {
+        centers.reserve(MythenDetectorSpecifications::max_modules());
+        conversions.reserve(MythenDetectorSpecifications::max_modules());
+        offsets.reserve(MythenDetectorSpecifications::max_modules());
+
+        num_strips = MythenDetectorSpecifications::max_modules() *
+                     MythenDetectorSpecifications::strips_per_module();
     }
 
     /** reads the historical Detector Group (DG) parameters from file **/
@@ -65,6 +252,15 @@ class AngleCalibration {
         const double normal_distance, const double module_center_distance,
         const double angle, const size_t strip_index);
 
+    /** calculates diffraction angle from EE module parameters (used in Beer's
+     * Law)
+     * @param strip_index local strip index of module
+     */
+    double diffraction_angle_from_DG_parameters(const double center,
+                                                const double conversion,
+                                                const double offset,
+                                                const size_t strip_index);
+
     /** calculated the strip width expressed as angle [degrees]
      * @param strip_index gloabl strip index of detector
      */
@@ -74,26 +270,38 @@ class AngleCalibration {
     size_t
     global_to_local_strip_index_conversion(const size_t global_strip_index) {
         const size_t module_index =
-            global_strip_index / MythenSpecifications::strips_per_module;
+            global_strip_index /
+            MythenDetectorSpecifications::strips_per_module();
         // local strip index in module
         size_t local_strip_index =
             global_strip_index -
-            module_index * MythenSpecifications::strips_per_module;
+            module_index * MythenDetectorSpecifications::strips_per_module();
         // switch if indexing is in clock-wise direction
         local_strip_index =
             signbit(conversions[module_index])
-                ? MythenSpecifications::strips_per_module - local_strip_index
+                ? MythenDetectorSpecifications::strips_per_module() -
+                      local_strip_index
                 : local_strip_index;
 
         return local_strip_index;
     }
 
+    /**
+     * redistributes photon counts with of histogram using one bin per strip to
+     * histogram with fixed size angle bins
+     * @param filename where histogram is stored
+     */
+    // TODO: pass frame or filename?
+    void
+    redistribute_photon_counts_to_fixed_angle_bins(const std::string &filename);
+
   protected:
     // 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
@@ -103,11 +311,16 @@ class AngleCalibration {
     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;
+
+    ssize_t num_strips{};
+
+    double exposure_rate = 1. / MythenDetectorSpecifications::exposure_time();
 };
 
-// read hdf5 files - > do they store the histogram? what angles do they store?
-
-// TODO what kind of file does it need to support? - probably code a csv parser
 void AngleCalibration::read_initial_calibration_from_file(
     const std::string &filename) {
 
@@ -176,9 +389,10 @@ parameters AngleCalibration::convert_to_EE_parameters() {
 std::tuple<double, double, double>
 AngleCalibration::convert_to_EE_parameters(const size_t module_index) {
     const double normal_distance =
-        centers[module_index] * MythenSpecifications::pitch;
+        centers[module_index] * MythenDetectorSpecifications::pitch();
     const double module_center_distance =
-        MythenSpecifications::pitch / std::abs(conversions[module_index]);
+        MythenDetectorSpecifications::pitch() /
+        std::abs(conversions[module_index]);
     const double angle =
         offsets[module_index] +
         180.0 / M_PI * centers[module_index] *
@@ -193,24 +407,33 @@ 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) {
+    return offset + 180.0 / M_PI *
+                        (center * conversion -
+                         atan((center - strip_index) * conversion));
+}
+
 double AngleCalibration::diffraction_angle_from_EE_parameters(
     const double normal_distance, const double module_center_distance,
     const double angle, const size_t strip_index) {
 
-    return angle - 180.0 / M_PI *
+    return angle -
-                       atan((module_center_distance -
+           180.0 / M_PI *
-                             MythenSpecifications::pitch * strip_index) /
+               atan((module_center_distance -
-                            normal_distance); // TODO: why is it minus
+                     MythenDetectorSpecifications::pitch() * strip_index) /
-                                              // is it defined counter
+                    normal_distance); // TODO: why is it minus
-                                              // clockwise? thought
+                                      // is it defined counter
-                                              // should have a flipped
+                                      // clockwise? thought
-                                              // sign
+                                      // should have a flipped
+                                      // sign
 }
 
 double AngleCalibration::angular_strip_width(const size_t strip_index) {
 
     const size_t module_index =
-        strip_index / MythenSpecifications::strips_per_module;
+        strip_index / MythenDetectorSpecifications::strips_per_module();
 
     const auto [normal_distance, module_center_distance, angle] =
         convert_to_EE_parameters(module_index);
@@ -229,20 +452,49 @@ double AngleCalibration::angular_strip_width(const size_t strip_index) {
 }
 
 /*
-void AngleCalibration::Calibrate() {
+void AngleCalibration::redistribute_photon_counts_to_fixed_angle_bins(
+    const std::string &filename) {
 
-    // iterates over each strip
+    // TODO: do i want to have a MythenReader that reads everything e.g. angle
-    // skips the one which are not enabled or have a bad channel
+    // in read_frame()
-    // get module index
+    double detector_angle;      // read from file
-    // get index of strip in module - make sure to use function correctly based
+    double reference_intensity; // read from file - Whats the difference between
-    // on sign
+                                // this and flatfield
 
-    // then I read some images - probably the histogram - h5 file - need a
+    NDArray<uint32_t, 2> photon_counts; // read from file e.g read frame - maybe
-    // reader need to read the filesystem!! -
+                                        // keep it as a frame and not an NDArray
 
-    // we modify it with ff_corr computed by flatfield why?
+    ssize_t strategy = 0; // photon_counts is 2d do we have to loop over all
+                          // strategies as well - probably
 
+    if (photon_counts.shape()[0] != num_strips) {
+        throw std::runtime_error("wrong number of strips read");
+    }
 
+    NDArray<uint32_t, 2> inverse_normalized_flatfield =
+        flat_field->inverse_normalized_flatfield();
+
+    for (size_t strip_index = 0; strip_index < num_strips; ++strip_index) {
+
+        size_t module_index =
+            strip_index / MythenDetectorSpecifications::strips_per_module;
+
+        if (mythen_detector->bad_channels[strip_index] ||
+            !mythen_detector->enabled_modules[module_index])
+            continue;
+
+        double poisson_error =
+            sqrt(photon_counts[strip_index, strategy];) *
+            inverse_normalized_flatfield[strip_index, strategy] *
+            exposure_rate; // not sure what this is
+        double corrected_photon_count =
+            photon_counts[strip_index, strategy] *
+            inverse_normalized_flatfield[strip_index, strategy] *
+            exposure_rate; // not sure what this is
+
+        size_t local_strip_index =
+            global_to_local_strip_index_conversion(strip_index);
+    }
 }
 */
 

src/AngleCalibration.test.cpp

@@ -16,7 +16,10 @@ using namespace aare;
 class AngleCalibrationTestClass : public aare::AngleCalibration {
 
   public:
-    AngleCalibrationTestClass() = default;
+    AngleCalibrationTestClass(
+        std::shared_ptr<MythenDetectorSpecifications> mythen_detector_,
+        std::shared_ptr<FlatField> flat_field_)
+        : aare::AngleCalibration(mythen_detector_, flat_field_) {}
     ~AngleCalibrationTestClass() = default;
 
     std::vector<double> get_centers() { return centers; }
@@ -29,7 +32,14 @@ class AngleCalibrationTestClass : public aare::AngleCalibration {
 TEST_CASE("read initial angle calibration file",
           "[.anglecalibration][.fileread]") {
 
-    AngleCalibrationTestClass anglecalibration;
+    std::shared_ptr<MythenDetectorSpecifications> mythen_detector_ptr =
+        std::make_shared<MythenDetectorSpecifications>();
+
+    std::shared_ptr<FlatField> flat_field_ptr =
+        std::make_shared<FlatField>(mythen_detector_ptr);
+
+    AngleCalibrationTestClass anglecalibration(mythen_detector_ptr,
+                                               flat_field_ptr);
 
     std::string filename =
         "/home/mazzol_a/Documents/mythen3tools/beamline/"