diff --git a/CMakeLists.txt b/CMakeLists.txt index 6a726fa..2180f2b 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/include/aare/AngleCalibration.hpp b/include/aare/AngleCalibration.hpp index 86a987d..a96b4b9 100644 --- a/include/aare/AngleCalibration.hpp +++ b/include/aare/AngleCalibration.hpp @@ -11,6 +11,8 @@ #include #include +#include "FlatField.hpp" +#include "MythenDetectorSpecifications.hpp" #include "MythenFileReader.hpp" #include "NDArray.hpp" @@ -19,304 +21,62 @@ namespace aare { using parameters = std::tuple, std::vector, std::vector>; -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(std::array{num_strips_}, false); - - connected_modules = NDArray( - std::array{static_cast(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(std::array{num_strips_}, false); - - connected_modules = NDArray( - std::array{static_cast(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 get_bad_channels() { return bad_channels.view(); } - - NDView 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 bad_channels; - NDArray connected_modules; // connected modules -}; - -// TODO maybe template now its uint32 -class FlatField { - - public: - FlatField(std::shared_ptr mythen_detector_) - : mythen_detector(mythen_detector_) { - - flat_field = NDArray( - std::array{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 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(0.0, 0), - [&tolerance](std::pair acc, const auto &element) { - return element == 0 ? acc - : std::make_pair(acc.first + element, - acc.second + 1); - }); - - return sum / count; - } - - NDArray inverse_normalized_flatfield(double tolerance = 0.001) { - double mean = calculate_mean(tolerance); - - NDArray 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 normalized_flatfield(double tolerance = 0.001) { - double mean = calculate_mean(tolerance); - - NDArray 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 flat_field; // TODO: should be 2d - std::shared_ptr mythen_detector; -}; - class AngleCalibration { public: AngleCalibration( std::shared_ptr mythen_detector_, std::shared_ptr flat_field_, - std::shared_ptr 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 - } + std::shared_ptr mythen_file_reader_); /** set the histogram bin width [degrees] */ - void set_histogram_bin_width(double bin_width) { - histogram_bin_width = bin_width; + void set_histogram_bin_width(double 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 get_histogram_bin_width() const; - 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 get_centers() { return centers; } + std::vector get_centers() const; + std::vector get_conversions() const; - std::vector get_conversions() { return conversions; } + std::vector get_offsets() const; - std::vector get_offsets() { return offsets; } + NDView get_new_photon_counts() const; + + NDView 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 - convert_to_EE_parameters(const size_t module_index); + convert_to_EE_parameters(const size_t module_index) const; std::tuple 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, - const size_t strip_index); - - double angular_strip_width_from_EE_parameters(const size_t strip_index); + 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 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 - global_to_local_strip_index_conversion(const size_t global_strip_index); - - /** - * 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); + 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 @@ -381,11 +127,9 @@ class AngleCalibration { void redistribute_photon_counts_to_fixed_angle_bins( const MythenFrame &frame, NDView bin_counts, NDView new_statistical_weights, NDView new_errors, - NDView inverse_nromalized_flatfield); + NDView inverse_nromalized_flatfield) const; - void write_to_file(const std::string &filename); - - protected: + private: // 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 mythen_file_reader; // TODO replace by FileInterface ptr }; diff --git a/include/aare/FlatField.hpp b/include/aare/FlatField.hpp new file mode 100644 index 0000000..cb85934 --- /dev/null +++ b/include/aare/FlatField.hpp @@ -0,0 +1,112 @@ + +/** + * stores flatfield for angle calibration + */ + +#pragma once +#include +#include + +#include +#include +#include +#include +#include + +#include "MythenDetectorSpecifications.hpp" +#include "NDArray.hpp" + +namespace aare { +// TODO maybe template now its uint32 +class FlatField { + + public: + FlatField(std::shared_ptr mythen_detector_) + : mythen_detector(mythen_detector_) { + + flat_field = NDArray( + std::array{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 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(0.0, 0), + [&tolerance](std::pair acc, const auto &element) { + return element == 0 ? acc + : std::make_pair(acc.first + element, + acc.second + 1); + }); + + return sum / count; + } + + NDArray + inverse_normalized_flatfield(double tolerance = 0.001) const { + double mean = calculate_mean(tolerance); + + NDArray 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 normalized_flatfield(double tolerance = 0.001) const { + double mean = calculate_mean(tolerance); + + NDArray 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 flat_field; // TODO: should be 2d + std::shared_ptr mythen_detector; +}; +} // namespace aare \ No newline at end of file diff --git a/include/aare/Hdf5FileReader.hpp b/include/aare/Hdf5FileReader.hpp index ad86b84..3b2fe5f 100644 --- a/include/aare/Hdf5FileReader.hpp +++ b/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 get_shape() { return shape; } + std::vector 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, - std::optional subset = std::nullopt) { + void + read_into_buffer(std::byte *buffer, + std::optional 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 NDArray store_as_ndarray() { + template + NDArray 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); diff --git a/include/aare/MythenDetectorSpecifications.hpp b/include/aare/MythenDetectorSpecifications.hpp new file mode 100644 index 0000000..f93be95 --- /dev/null +++ b/include/aare/MythenDetectorSpecifications.hpp @@ -0,0 +1,150 @@ +#pragma once +#include +#include + +#include +#include +#include +#include + +#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(std::array{num_strips_}, false); + + connected_modules = NDArray( + std::array{static_cast(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(std::array{num_strips_}, false); + + connected_modules = NDArray( + std::array{static_cast(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 get_bad_channels() const { return bad_channels.view(); } + + NDView 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 bad_channels; + NDArray connected_modules; // connected modules +}; + +} // namespace aare \ No newline at end of file diff --git a/include/aare/MythenFileReader.hpp b/include/aare/MythenFileReader.hpp index 0602520..6d78da3 100644 --- a/include/aare/MythenFileReader.hpp +++ b/include/aare/MythenFileReader.hpp @@ -43,7 +43,7 @@ class MythenFileReader : public HDF5FileReader { frame.photon_counts = dataset_photon_count.store_as_ndarray(); - ++frame.photon_counts; + ++frame.photon_counts; // Why though? auto dataset_detector_angle = get_dataset("/entry/instrument/NDAttributes/DetectorAngle"); diff --git a/include/aare/NDArray.hpp b/include/aare/NDArray.hpp index 31e15d8..bd73f94 100644 --- a/include/aare/NDArray.hpp +++ b/include/aare/NDArray.hpp @@ -139,6 +139,9 @@ class NDArray : public ArrayExpr, Ndim> { NDArray operator>(const NDArray &other); + bool equals(const NDArray &other, + const T tolerance = std::numeric_limits::epsilon()) const; + bool operator==(const NDArray &other) const; bool operator!=(const NDArray &other) const; @@ -438,4 +441,33 @@ NDArray load(const std::string &pathname, return img; } +template +NDArray load_non_binary_file(const std::string &filename, + const std::array shape) { + std::string word; + NDArray 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( + 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 \ No newline at end of file diff --git a/include/aare/NDView.hpp b/include/aare/NDView.hpp index 56054e2..629eee4 100644 --- a/include/aare/NDView.hpp +++ b/include/aare/NDView.hpp @@ -1,11 +1,12 @@ #pragma once -#include "aare/defs.hpp" #include "aare/ArrayExpr.hpp" +#include "aare/defs.hpp" #include #include #include #include +#include #include #include #include @@ -17,7 +18,8 @@ namespace aare { template using Shape = std::array; // TODO! fix mismatch between signed and unsigned -template Shape make_shape(const std::vector &shape) { +template +Shape make_shape(const std::vector &shape) { if (shape.size() != Ndim) throw std::runtime_error("Shape size mismatch"); Shape arr; @@ -25,14 +27,18 @@ template Shape make_shape(const std::vector &shape) return arr; } -template ssize_t element_offset(const Strides & /*unused*/) { return 0; } +template +ssize_t element_offset(const Strides & /*unused*/) { + return 0; +} template ssize_t element_offset(const Strides &strides, ssize_t i, Ix... index) { return i * strides[Dim] + element_offset(strides, index...); } -template std::array c_strides(const std::array &shape) { +template +std::array c_strides(const std::array &shape) { std::array strides{}; std::fill(strides.begin(), strides.end(), 1); for (ssize_t i = Ndim - 1; i > 0; --i) { @@ -41,14 +47,16 @@ template std::array c_strides(const std::array std::array make_array(const std::vector &vec) { +template +std::array make_array(const std::vector &vec) { assert(vec.size() == Ndim); std::array arr{}; std::copy_n(vec.begin(), Ndim, arr.begin()); return arr; } -template class NDView : public ArrayExpr, Ndim> { +template +class NDView : public ArrayExpr, Ndim> { public: NDView() = default; ~NDView() = default; @@ -57,17 +65,23 @@ template class NDView : public ArrayExpr shape) : buffer_(buffer), strides_(c_strides(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 &shape) - // : buffer_(buffer), strides_(c_strides(make_array(shape))), shape_(make_array(shape)), - // size_(std::accumulate(std::begin(shape), std::end(shape), 1, std::multiplies<>())) {} + // : buffer_(buffer), + // strides_(c_strides(make_array(shape))), + // shape_(make_array(shape)), + // size_(std::accumulate(std::begin(shape), std::end(shape), 1, + // std::multiplies<>())) {} - template std::enable_if_t operator()(Ix... index) { + template + std::enable_if_t operator()(Ix... index) { return buffer_[element_offset(strides_, index...)]; } - template std::enable_if_t operator()(Ix... index) const { + template + std::enable_if_t operator()(Ix... index) const { return buffer_[element_offset(strides_, index...)]; } @@ -92,18 +106,37 @@ template class NDView : public ArrayExpr &other, const T tolerance) const { + if (shape_ != other.shape_) + return false; + + using SignedT = typename make_signed::type; + + for (uint32_t i = 0; i != size_; ++i) + if (std::abs(static_cast(buffer_[i]) - + static_cast(other.buffer_[i])) > tolerance) + return false; + + return true; + } + NDView &operator+=(const T val) { return elemenwise(val, std::plus()); } NDView &operator-=(const T val) { return elemenwise(val, std::minus()); } - NDView &operator*=(const T val) { return elemenwise(val, std::multiplies()); } - NDView &operator/=(const T val) { return elemenwise(val, std::divides()); } + NDView &operator*=(const T val) { + return elemenwise(val, std::multiplies()); + } + NDView &operator/=(const T val) { + return elemenwise(val, std::divides()); + } - NDView &operator/=(const NDView &other) { return elemenwise(other, std::divides()); } + NDView &operator/=(const NDView &other) { + return elemenwise(other, std::divides()); + } - - template - NDView& operator=(const std::array &arr) { - if(size() != static_cast(arr.size())) - throw std::runtime_error(LOCATION + "Array and NDView size mismatch"); + template NDView &operator=(const std::array &arr) { + if (size() != static_cast(arr.size())) + throw std::runtime_error(LOCATION + + "Array and NDView size mismatch"); std::copy(arr.begin(), arr.end(), begin()); return *this; } @@ -147,13 +180,15 @@ template class NDView : public ArrayExpr shape_{}; uint64_t size_{}; - template NDView &elemenwise(T val, BinaryOperation op) { + template + NDView &elemenwise(T val, BinaryOperation op) { for (uint64_t i = 0; i != size_; ++i) { buffer_[i] = op(buffer_[i], val); } return *this; } - template NDView &elemenwise(const NDView &other, BinaryOperation op) { + template + 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 void NDView::print_all() const { } } - template -std::ostream& operator <<(std::ostream& os, const NDView& arr){ +std::ostream &operator<<(std::ostream &os, const NDView &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& arr){ return os; } +template NDView make_view(std::vector &vec) { + return NDView(vec.data(), {static_cast(vec.size())}); +} -template -NDView make_view(std::vector& vec){ - return NDView(vec.data(), {static_cast(vec.size())}); +template +void save(NDView img, const std::string &pathname) { + std::ofstream f; + f.open(pathname, std::ios::binary); + f.write(reinterpret_cast(img.data()), img.size() * sizeof(T)); + f.close(); } } // namespace aare \ No newline at end of file diff --git a/include/aare/defs.hpp b/include/aare/defs.hpp index 01d291b..ea02059 100644 --- a/include/aare/defs.hpp +++ b/include/aare/defs.hpp @@ -3,16 +3,16 @@ #include "aare/Dtype.hpp" #include -#include #include #include #include +#include #include #include +#include #include #include - /** * @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 struct t_xy { }; using xy = t_xy; - /** - * @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 module_pixel_0; }; -struct ROI{ +struct ROI { ssize_t xmin{}; ssize_t xmax{}; ssize_t ymin{}; ssize_t ymax{}; - + ssize_t height() const { return ymax - ymin; } ssize_t width() const { return xmax - xmin; } bool contains(ssize_t x, ssize_t y) const { return x >= xmin && x < xmax && y >= ymin && y < ymax; } - }; - +}; using dynamic_shape = std::vector; -//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; +template > struct make_signed { + using type = T; +}; + +template struct make_signed { + using type = std::make_signed_t; +}; + } // namespace aare \ No newline at end of file diff --git a/src/AngleCalibration.cpp b/src/AngleCalibration.cpp index 28a1fcd..632b8d5 100644 --- a/src/AngleCalibration.cpp +++ b/src/AngleCalibration.cpp @@ -2,6 +2,53 @@ namespace aare { +AngleCalibration::AngleCalibration( + std::shared_ptr mythen_detector_, + std::shared_ptr flat_field_, + std::shared_ptr 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 AngleCalibration::get_centers() const { return centers; } + +std::vector AngleCalibration::get_conversions() const { + return conversions; +} + +std::vector AngleCalibration::get_offsets() const { return offsets; } + +NDView AngleCalibration::get_new_photon_counts() const { + return new_photon_counts.view(); +} + +NDView 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 normal_distances(centers.size()); @@ -68,14 +115,14 @@ parameters AngleCalibration::convert_to_EE_parameters() { } std::tuple -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 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 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 = - global_to_local_strip_index_conversion(strip_index); + const size_t local_strip_index) const { 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 size_t local_strip_index = - global_to_local_strip_index_conversion(strip_index); + const double angle, const size_t local_strip_index) const { 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(std::array{num_bins}); + // TODO: maybe group these into a 2d array - better cache usage NDArray bin_counts(std::array{num_bins}, 0.0); NDArray new_statistical_weights(std::array{num_bins}, 0.0); - NDArray new_errors(std::array{num_bins}, 0.0); NDArray 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 bin_counts, NDView new_statistical_weights, NDView new_errors, - NDView inverse_normalized_flatfield) { + NDView 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 = - angular_strip_width_from_DG_parameters(strip_index); + 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 / @@ -324,8 +345,11 @@ void AngleCalibration::redistribute_photon_counts_to_fixed_angle_bins( } } -void AngleCalibration::write_to_file(const std::string &filename) { - std::ofstream output_file(filename); +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!" @@ -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::epsilon()) { + if (new_photon_counts[i] <= std::numeric_limits::epsilon() && + store_nonzero_bins) { continue; } diff --git a/src/AngleCalibration.test.cpp b/src/AngleCalibration.test.cpp index 162a69c..071dde0 100644 --- a/src/AngleCalibration.test.cpp +++ b/src/AngleCalibration.test.cpp @@ -18,101 +18,18 @@ using namespace aare; -template -NDArray read_into_array(const std::string &filename, - const std::array size) { - std::string word; - NDArray 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 > struct safe_make_signed { - using type = T; -}; - -template struct safe_make_signed { - using type = std::make_signed_t; -}; - -template -bool check_equality_of_arrays(NDView array1, NDView array2) { - bool equal = true; - - using SignedT = typename safe_make_signed::type; - - for (ssize_t i = 0; i < array1.size(); ++i) { - if (std::abs(static_cast(array1[i]) - - static_cast(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 mythen_detector_, - std::shared_ptr flat_field_, - std::shared_ptr mythen_file_reader_) - : AngleCalibration(mythen_detector_, flat_field_, mythen_file_reader_) { - } - ~AngleCalibrationTestClass() = default; - - std::vector get_centers() { return centers; } - - std::vector get_conversions() { return conversions; } - - std::vector 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 mythen_detector_ptr = std::make_shared(); - std::shared_ptr flat_field_ptr = - std::make_shared(mythen_detector_ptr); + AngleCalibration anglecalibration(mythen_detector_ptr, + std::shared_ptr{}, + std::shared_ptr{}); - std::shared_ptr mythen_file_reader_ptr = - std::make_shared(fpath, - "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"; + std::string filename = test_data_path() / "AngleCalibration_Test_Data" / + "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 mythen_detector_ptr = std::make_shared(); @@ -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]") { - auto fpath = test_data_path() / "AngleCalibration_Test_Data"; - - REQUIRE(std::filesystem::exists(fpath)); - - std::shared_ptr mythen_detector_ptr = - std::make_shared(); - - 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 expected_flatfield = read_into_array( - 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::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 mythen_detector_ptr = - std::make_shared(); - - 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 expected_inverseflatfield = read_into_array( - 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]") { +TEST_CASE("compare result with python code", "[.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( - "cpp_new_photon_counts.xye"); // TODO adjust output path + // 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( + expected_filename_errors, std::array{new_num_bins}); + + auto python_output_photons = load( + expected_filename_photons, std::array{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]") { - auto expected_filename = - test_data_path() / "AngleCalibration_Test_Data" / "out2.xye"; - - REQUIRE(std::filesystem::exists(expected_filename)); - - auto expected_array = read_into_array( - expected_filename.string(), std::array{61440, 3}); - - auto filename = - std::filesystem::current_path() / "../build/cpp_new_photon_counts.xye"; - - // auto array = load(filename, std::array{61440, 3}); - auto array = read_into_array(filename.string(), - std::array{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)); +TEST_CASE("check conversion from DG to EE parameters", "[.anglecalibration]") { std::shared_ptr mythen_detector_ptr = std::make_shared(); - std::shared_ptr flat_field_ptr = - std::make_shared(mythen_detector_ptr); + AngleCalibration anglecalibration(mythen_detector_ptr, + std::shared_ptr{}, + std::shared_ptr{}); - std::shared_ptr mythen_file_reader_ptr = - std::make_shared(fpath, - "ang1up_22keV_LaB60p3mm_48M_a_0"); - - AngleCalibration anglecalibration(mythen_detector_ptr, flat_field_ptr, - 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; + // 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, 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)); }