From f0e1ed3e2a35d3063b629ffbdec5c782ccf10655 Mon Sep 17 00:00:00 2001 From: Alice Date: Fri, 20 Jun 2025 21:28:35 +0200 Subject: [PATCH] moved anglecalibration stuff to own repository still keep some changes --- CMakeLists.txt | 33 -- acquire_flatfield_data.py | 64 --- include/aare/AngleCalibration.hpp | 164 -------- include/aare/CustomFiles.hpp | 59 --- include/aare/FlatField.hpp | 215 ---------- include/aare/Hdf5FileReader.hpp | 212 ---------- include/aare/MythenDetectorSpecifications.hpp | 157 -------- include/aare/MythenFileReader.hpp | 82 ---- src/AngleCalibration.cpp | 377 ------------------ src/AngleCalibration.test.cpp | 265 ------------ src/CustomFiles.cpp | 89 ----- src/Hdf5FileReader.test.cpp | 110 ----- src/MythenFileReader.test.cpp | 33 -- 13 files changed, 1860 deletions(-) delete mode 100644 acquire_flatfield_data.py delete mode 100644 include/aare/AngleCalibration.hpp delete mode 100644 include/aare/CustomFiles.hpp delete mode 100644 include/aare/FlatField.hpp delete mode 100644 include/aare/Hdf5FileReader.hpp delete mode 100644 include/aare/MythenDetectorSpecifications.hpp delete mode 100644 include/aare/MythenFileReader.hpp delete mode 100644 src/AngleCalibration.cpp delete mode 100644 src/AngleCalibration.test.cpp delete mode 100644 src/CustomFiles.cpp delete mode 100644 src/Hdf5FileReader.test.cpp delete mode 100644 src/MythenFileReader.test.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 33c033b..64c08ad 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -62,7 +62,6 @@ option(AARE_FETCH_CATCH "Use FetchContent to download catch2" ON) option(AARE_FETCH_JSON "Use FetchContent to download nlohmann::json" ON) option(AARE_FETCH_ZMQ "Use FetchContent to download libzmq" ON) option(AARE_FETCH_LMFIT "Use FetchContent to download lmfit" ON) -option(AARE_FETCH_HDF5 "Use FetchContent to download hdf5-devel" OFF) #Convenience option to use system libraries only (no FetchContent) @@ -74,7 +73,6 @@ if(AARE_SYSTEM_LIBRARIES) set(AARE_FETCH_CATCH OFF CACHE BOOL "Disabled FetchContent for catch2" FORCE) set(AARE_FETCH_JSON OFF CACHE BOOL "Disabled FetchContent for nlohmann::json" FORCE) set(AARE_FETCH_ZMQ OFF CACHE BOOL "Disabled FetchContent for libzmq" FORCE) - set(AARE_FETCH_HDF5 OFF CACHE BOOL "Disabled FetchContent for hdf5" FORCE) # Still fetch lmfit when setting AARE_SYSTEM_LIBRARIES since this is not available # on conda-forge endif() @@ -221,23 +219,6 @@ else() find_package(nlohmann_json 3.11.3 REQUIRED) endif() -if(AARE_FETCH_HDF5) - message(FATAL_ERROR "Fetching HDF5 via FetchContent is not supported here. Please install it via your system. - For Ubuntu: sudo apt install libhdf5-dev - For Red Hat: sudo dnf install hdf5-devel - For MacOS: brew install hdf5") -else() - find_package(HDF5 QUIET COMPONENTS CXX) - if (HDF5_FOUND) - message(STATUS "Found HDF5: ${HDF5_INCLUDE_DIRS}") - else() - message(FATAL_ERROR "HDF5 was NOT found! Please install it via your system - For Ubuntu: sudo apt install libhdf5-dev - For Red Hat: sudo dnf install hdf5-devel - For MacOS: brew install hdf5") - endif() -endif() - include(GNUInstallDirs) # If conda build, always set lib dir to 'lib' @@ -359,7 +340,6 @@ endif() ###------------------------------------------------------------------------------------------ set(PUBLICHEADERS - include/aare/AngleCalibration.hpp include/aare/ArrayExpr.hpp include/aare/CalculateEta.hpp include/aare/Cluster.hpp @@ -367,7 +347,6 @@ set(PUBLICHEADERS include/aare/ClusterFile.hpp include/aare/CtbRawFile.hpp include/aare/ClusterVector.hpp - include/aare/CustomFiles.hpp include/aare/decode.hpp include/aare/defs.hpp include/aare/Dtype.hpp @@ -375,14 +354,10 @@ set(PUBLICHEADERS include/aare/Fit.hpp include/aare/FileInterface.hpp include/aare/FilePtr.hpp - include/aare/FlatField.hpp include/aare/Frame.hpp include/aare/GainMap.hpp include/aare/geo_helpers.hpp - include/aare/Hdf5FileReader.hpp include/aare/JungfrauDataFile.hpp - include/aare/MythenDetectorSpecifications.hpp - include/aare/MythenFileReader.hpp include/aare/NDArray.hpp include/aare/NDView.hpp include/aare/NumpyFile.hpp @@ -398,9 +373,7 @@ set(PUBLICHEADERS set(SourceFiles - ${CMAKE_CURRENT_SOURCE_DIR}/src/AngleCalibration.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/CtbRawFile.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/src/CustomFiles.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/defs.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/decode.cpp @@ -417,7 +390,6 @@ set(SourceFiles ${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/RawSubFile.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/RawMasterFile.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/src/utils/task.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/utils/ifstream_helpers.cpp ) @@ -436,13 +408,11 @@ target_link_libraries( PUBLIC fmt::fmt nlohmann_json::nlohmann_json - HDF5::HDF5 ${STD_FS_LIB} # from helpers.cmake PRIVATE aare_compiler_flags Threads::Threads $ - ) set_target_properties(aare_core PROPERTIES @@ -457,14 +427,12 @@ endif() if(AARE_TESTS) set(TestSources ${CMAKE_CURRENT_SOURCE_DIR}/src/algorithm.test.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/src/AngleCalibration.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/defs.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/decode.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/geo_helpers.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/RawMasterFile.test.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/src/Hdf5FileReader.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/NDArray.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/NDView.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFinder.test.cpp @@ -475,7 +443,6 @@ if(AARE_TESTS) ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFinderMT.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Pedestal.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/JungfrauDataFile.test.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/src/MythenFileReader.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.test.cpp diff --git a/acquire_flatfield_data.py b/acquire_flatfield_data.py deleted file mode 100644 index 63af684..0000000 --- a/acquire_flatfield_data.py +++ /dev/null @@ -1,64 +0,0 @@ -from slsdet import Mythen3,timingMode,detectorSettings,runStatus,dacIndex,scanParameters -#from patterntools.zmqreceiver import ZmqReceiver -#from detConf_module import * -import numpy as np -from slsdet.lookup import view, find -from epics import caput, caget -import time - - -d = Mythen3() -#d.detsize = [10, 1] - -d.hostname = 'localhost:2000+localhost:2002' - -d.rx_hostname = 'localhost' - -#d.rx_tcpport = 2012 - -d.udp_dstport = 5006 - -d.udp_dstip = 'auto' - -d.udp_srcip = 'auto' - -#rx=makeReceiver(d) - -nmod=len(d.hostname) -print("num modules: ", nmod) - -d.frames = 3 -d.exptime = 5 - -d.fwrite=1 -#d.counters = [0, 1] -print(d.counters) - -d.fpath = "~/Documents/tmp/Flatfieldacquisition" -d.fwrite = 1 - -#where do i set the strips? -#d.detsize = 2 do I have to set this - -d.startReceiver() -d.startDetector() -for angle in [87,2]: - d.fname = 'run_'+str(angle) - caput('BL11I-MO-DIFF-01:DELTA.VAL',angle,wait=False) - print("moving detector to ",angle) - d.acquire() - time.sleep(d.exptime) - while d.status != runStatus.IDLE: - time.sleep(0.01) - - nf=np.min(d.rx_framescaught) - ang=caget('BL11I-MO-DIFF-01:DELTA.RBV') - print("angle is: ", ang) - print("caught frames: ", nf) - -d.stopReceiver() - -#dd, hh = rx[imod].receive_one_frame() what is this - does this exist in cpp? - - - diff --git a/include/aare/AngleCalibration.hpp b/include/aare/AngleCalibration.hpp deleted file mode 100644 index a96b4b9..0000000 --- a/include/aare/AngleCalibration.hpp +++ /dev/null @@ -1,164 +0,0 @@ -#pragma once -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include - -#include "FlatField.hpp" -#include "MythenDetectorSpecifications.hpp" -#include "MythenFileReader.hpp" -#include "NDArray.hpp" - -namespace aare { - -using parameters = - std::tuple, std::vector, std::vector>; - -class AngleCalibration { - - public: - AngleCalibration( - std::shared_ptr mythen_detector_, - std::shared_ptr flat_field_, - std::shared_ptr mythen_file_reader_); - - /** set the histogram bin width [degrees] */ - void set_histogram_bin_width(double bin_width); - - double get_histogram_bin_width() const; - - ssize_t get_new_num_bins() const; - - /** reads the historical Detector Group (DG) parameters from file **/ - void read_initial_calibration_from_file(const std::string &filename); - - std::vector get_centers() const; - std::vector get_conversions() const; - - std::vector get_offsets() const; - - 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() const; - - std::tuple - 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; - - /** converts DG parameters to BC parameters e.g. best computing - * parameters */ - parameters convert_to_BC_parameters() const; - - /** - * calculates new histogram with fixed sized angle bins - * for several acquisitions at different detector angles for given frame - * indices - * @param start_frame_index, end_frame_index gives range of frames - */ - void - calculate_fixed_bin_angle_width_histogram(const size_t start_frame_index, - const size_t end_frame_index); - - void write_to_file(const std::string &filename, - const bool store_nonzero_bins = false, - const std::filesystem::path &filepath = - std::filesystem::current_path()) const; - - /** calculates diffraction angle from EE module parameters (used in - * Beer's Law) - * @param strip_index local strip index of module - */ - double diffraction_angle_from_EE_parameters( - const double module_center_distance, const double normal_distance, - const double angle, const size_t strip_index, - const double distance_to_strip = 0) const; - - /** calculates diffraction angle from EE module parameters (used in - * Beer's Law) - * @param center module center - * @param conversion module conversion - * @param offset module offset - * @param strip_index local strip index of module - * @param distance_to_strip distance to strip given by strip_index and - * module -> note needs to be small enough to be in the respective module - */ - double diffraction_angle_from_DG_parameters( - const double center, const double conversion, const double offset, - const size_t strip_index, const double distance_to_strip = 0) const; - - /** calculated the strip width expressed as angle [degrees] - * @param strip_index local strip index of module - */ - double angular_strip_width_from_DG_parameters( - const double center, const double conversion, const double offset, - const size_t local_strip_index) const; - - double angular_strip_width_from_EE_parameters( - const double module_center_distance, const double normal_distance, - const double angle, const size_t local_strip_index) const; - - protected: - /** converts global strip index to local strip index of that module */ - size_t global_to_local_strip_index_conversion( - const size_t global_strip_index) const; - - /** - * redistributes photon counts with of histogram using one bin per strip - * to histogram with fixed size angle bins - * @param frame MythenFrame storing data from image - * @param bin_counts accumulate new photon counts - * @param new_statistical_weights accumulate new statistical weights - * @param new_errors accumulate new_errors - */ - void redistribute_photon_counts_to_fixed_angle_bins( - const MythenFrame &frame, NDView bin_counts, - NDView new_statistical_weights, NDView new_errors, - NDView inverse_nromalized_flatfield) const; - - private: - // TODO: Design maybe have a struct with three vectors, store all three - // sets of parameters as member variables - - // TODO: check if interpretation and units are correct - // historical DG parameters - // TODO change to NDArray - std::vector centers; // orthogonal projection of sample onto - // detector (given in strip number) [mm] - // D/pitch - std::vector conversions; // pitch/(normal distance from sample - // to detector (R)) [mm] - // //used for easy conversion - std::vector - offsets; // position of strip zero relative to sample [degrees] phi - // - 180/pi*D/R TODO: expected an arcsin(D/R)? - - std::shared_ptr mythen_detector; - - std::shared_ptr flat_field; - - NDArray new_photon_counts; - NDArray new_photon_count_errors; - - double histogram_bin_width = 0.0036; // [degrees] - - ssize_t num_bins{}; - - std::shared_ptr - mythen_file_reader; // TODO replace by FileInterface ptr -}; - -} // namespace aare diff --git a/include/aare/CustomFiles.hpp b/include/aare/CustomFiles.hpp deleted file mode 100644 index aefb9a0..0000000 --- a/include/aare/CustomFiles.hpp +++ /dev/null @@ -1,59 +0,0 @@ -#include "Dtype.hpp" -#include "FileInterface.hpp" -#include -#include -#include -#include -#include - -namespace aare { - -// TODO a lot to overload -class CustomMythenFile : public FileInterface { - - public: - CustomMythenFile(const std::string &filename, const ssize_t rows, - const ssize_t cols = 1, const std::string &mode = "r"); - - ~CustomMythenFile(); - - Frame read_frame() override; - - Frame read_frame(size_t frame_number) override; - - std::vector read_n(size_t n_frames) override; - - void read_into(std::byte *image_buf) override; - - void read_into(std::byte *image_buf, size_t n_frames) override; - - size_t frame_number(size_t frame_index) override; - - size_t bytes_per_frame() override; - - size_t pixels_per_frame() override; - - void seek(size_t frame_number) override; - - size_t tell() override; - - size_t total_frames() const override; - - size_t rows() const override; - size_t cols() const override; - - size_t bitdepth() const override; - - DetectorType detector_type() const override; - - private: - std::string m_filename{}; - std::ifstream m_file{}; - ssize_t m_num_strips{}; - // uint8_t m_num_counts{}; TODO extend - ssize_t m_rows{}; - ssize_t m_cols{}; - static const Dtype m_dtype; - static const DetectorType det_type = DetectorType::Mythen3; -}; -} // namespace aare \ No newline at end of file diff --git a/include/aare/FlatField.hpp b/include/aare/FlatField.hpp deleted file mode 100644 index 367cf7c..0000000 --- a/include/aare/FlatField.hpp +++ /dev/null @@ -1,215 +0,0 @@ - -/** - * stores flatfield for angle calibration - */ - -#pragma once -#include -#include - -#include -#include -#include -#include -#include -#include - -#include "File.hpp" -#include "MythenDetectorSpecifications.hpp" -#include "NDArray.hpp" - -namespace aare { - -template struct custom_file_compatibility { - static constexpr bool value = - std::is_base_of::value && - std::is_constructible::value; -}; - -template -constexpr bool custom_file_compatibility_v = - custom_file_compatibility::value; - -class FlatField { - - public: - FlatField(std::shared_ptr mythen_detector_) - : mythen_detector(mythen_detector_) { - - flat_field = NDArray( - std::array{mythen_detector->num_strips()}, 0); - } - - /** - * @brief sums up the photon counts for multiple acquisitions - * @param file_path: path to filesystem - the filesystem should contain - * multiple acqisitions for different detector angles acquired using - * slsReceiver and mythen3Detector //TODO: constructor needs to be the same - * - ugly - */ - // TODO unsure about design - maybe scientist has to give file with paths - // one wants to accumulate - what to do with strange file format? - // TODO: adjust for different counters!!!! - void create_flatfield_from_rawfilesystem( - const std::filesystem::path &file_path) { - - try { - for (const auto &file_in_path : - std::filesystem::recursive_directory_iterator(file_path)) { - if (file_in_path.is_regular_file()) { - std::string filename = - file_in_path.path().filename().string(); - if (filename.find("master") != std::string::npos) { - File f(file_in_path); - auto frames = f.read_n(f.total_frames()); - for (const auto &frame : frames) { - if (frame.rows() * frame.cols() != - mythen_detector->num_strips() * - mythen_detector->num_counters()) { - throw std::runtime_error(fmt::format( - "sizes mismatch. Expect a size of " - "{} - frame has a size of {}", - mythen_detector->num_strips() * - mythen_detector->num_counters(), - frame.rows() * frame.cols())); - } - for (ssize_t row = 0; row < frame.rows(); ++row) - for (ssize_t col = 0; col < frame.cols(); - ++col) { - flat_field(row * frame.cols() + col) += - *reinterpret_cast( - frame.pixel_ptr( - row, - col)); // TODO inefficient as - // one has to copy twice - // into frame and into - // flat_field - } - } - } - } - } - } catch (const std::filesystem::filesystem_error &e) { - std::cerr << "Filesystem error: " << e.what() - << '\n'; // TODO replace with log - } catch (const std::exception &e) { - std::cerr << "Runtime error: " << e.what() << '\n'; - } - } - - /** - * @brief sums up the photon counts for multiple acquisitions - * @param filelist: path to file that stores the file paths to the aquired - * data the list should contain multiple acquisitions for different detector - * angles - * @tparam CustomFile: Fileclass that inherits from aare::FileInterface - * class needs to overload read_frame and constructor needs to have the - * signature CustomFile(std::string filename, ssize_t rows, ssize_t cols, - * std::string reading_mode) rows, and cols define the intended image size - * to read - */ - template , void>> - void create_flatfield_from_filelist(const std::filesystem::path &filelist) { - std::ifstream file_filelist(filelist); - - try { - std::string filename; - while (std::getline(file_filelist, filename)) { - CustomFile file(filename, mythen_detector->num_strips(), 1); - Frame frame = file.read_frame(); - if (frame.rows() * frame.cols() != - mythen_detector->num_strips()) { - throw std::runtime_error( - fmt::format("sizes mismatch. Expect a size of " - "{} - frame has a size of {}", - mythen_detector->num_strips(), - frame.rows() * frame.cols())); - } - for (ssize_t row = 0; row < frame.rows(); ++row) - for (ssize_t col = 0; col < frame.cols(); ++col) { - flat_field(row * frame.cols() + col) += - *reinterpret_cast(frame.pixel_ptr( - row, - col)); // TODO inefficient as one has to copy - // twice into frame and into flat_field - } - } - file_filelist.close(); - } catch (const std::exception &e) { - std::cerr << "Error: " << e.what() << '\n'; - } - } - - /** - * @tparam CustomFile: Fileclass that inherits from aare::FileInterface - * class needs to overload read_frame and constructor needs to have the - * signature CustomFile(std::string filename, ssize_t rows, ssize_t cols, - * std::string reading_mode) rows, and cols define the intended image size - * to read - */ - - template , void>> - void read_flatfield_from_file(const std::string &filename) { - CustomFile file(filename, mythen_detector->num_strips(), 1); - file.read_into(reinterpret_cast(flat_field.data())); - } - - NDView get_flatfield() const { return flat_field.view(); } - - // TODO: remove tolerance - double calculate_mean(double tolerance) 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.0001) 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.0001) 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 deleted file mode 100644 index 3b2fe5f..0000000 --- a/include/aare/Hdf5FileReader.hpp +++ /dev/null @@ -1,212 +0,0 @@ -/************************************************ - * @file HD5FFileReader.hpp - * @short HDF5FileReader based on H5File object - ***********************************************/ - -#pragma once - -#include "Frame.hpp" -#include "NDArray.hpp" -#include -#include -#include -#include - -namespace aare { - -// return std::type_info -inline const std::type_info &deduce_cpp_type(const H5::DataType datatype) { - if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_UINT8.getId())) { - return typeid(uint8_t); - } else if (H5Tequal(datatype.getId(), - H5::PredType::NATIVE_UINT16.getId())) { - return typeid(uint16_t); - } else if (H5Tequal(datatype.getId(), - H5::PredType::NATIVE_UINT32.getId())) { - return typeid(uint32_t); - } else if (H5Tequal(datatype.getId(), - H5::PredType::NATIVE_UINT64.getId())) { - return typeid(uint64_t); - } else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT8.getId())) { - return typeid(int8_t); - } else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT16.getId())) { - return typeid(int16_t); - } else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT32.getId())) { - return typeid(int32_t); - } else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT64.getId())) { - return typeid(int64_t); - } else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_INT.getId())) { - return typeid(int); - } else if (H5Tequal(datatype.getId(), H5::PredType::IEEE_F64LE.getId())) { - return typeid(double); - } else if (H5Tequal(datatype.getId(), H5::PredType::IEEE_F32LE.getId())) { - return typeid(float); - } else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_FLOAT.getId())) { - return typeid(float); - } else if (H5Tequal(datatype.getId(), - H5::PredType::NATIVE_DOUBLE.getId())) { - return typeid(float); - } else if (H5Tequal(datatype.getId(), H5::PredType::NATIVE_CHAR.getId()) && - datatype.getId() == H5::PredType::NATIVE_CHAR.getId()) { - return typeid(char); - } else { - throw std::runtime_error("c++ type cannot be deduced"); - } -} - -struct Subset { - std::vector shape; - std::vector offset; // index where data subset should start -}; - -class HDF5Dataset { - - public: - HDF5Dataset(const std::string &datasetname_, const H5::DataSet dataset_) - : datasetname(datasetname_), dataset(dataset_) { - datatype = dataset.getDataType(); - - cpp_type = &deduce_cpp_type(datatype); - - dataspace = dataset.getSpace(); - rank = dataspace.getSimpleExtentNdims(); // number of dimensions - - shape.resize(rank); - dataspace.getSimpleExtentDims(shape.data(), nullptr); - } - - hsize_t get_shape(ssize_t index) const { return shape[index]; } - - std::vector get_shape() const { return shape; } - - H5::DataType get_datatype() const { return datatype; } - - const std::type_info *get_cpp_type() const { return cpp_type; } - - /** - * Reads subset of dataset into the buffer - * e.g. to read one 2d frame pass Subset({shape[1], shape[2]}, {frame_index, - * 0,0}) - */ - void - read_into_buffer(std::byte *buffer, - std::optional subset = std::nullopt) const { - - if (subset) { - // TODO treat scalar cases - if (static_cast(subset->offset.size()) != rank) { - throw std::runtime_error("provide an offset for" + - std::to_string(rank) + "dimensions"); - } - for (ssize_t i = 0; i < rank; ++i) { - hsize_t size = - i < rank - static_cast(subset->shape.size()) - ? 0 - : subset->shape[i - (rank - subset->shape.size())]; - if ((size + subset->offset[i]) > shape[i]) { - throw std::runtime_error( - "subset is too large or offset is out of bounds"); - } - } - - H5::DataSpace memspace(static_cast(subset->shape.size()), - subset->shape.data()); - - dataspace.selectHyperslab(H5S_SELECT_SET, subset->shape.data(), - subset->offset.data()); - dataset.read(buffer, datatype, memspace, dataspace); - } else { - dataset.read(buffer, datatype); - } - } - - Frame store_as_frame() const { - uint32_t rows{}, cols{}; - if (rank == 1) { - rows = 1; - // TODO overflow - cols = static_cast(shape[0]); - } else if (rank == 2) { - rows = static_cast(shape[0]); - cols = static_cast(shape[1]); - } else { - throw std::runtime_error("Frame only supports 2d images"); - } - - Frame frame(rows, cols, Dtype(*cpp_type)); - - read_into_buffer(frame.data()); - - return frame; - } - - template - NDArray store_as_ndarray() const { - if (NDim != rank) { - std::cout - << "Warning: dataset dimension and array dimension mismatch" - << std::endl; // TODO: replace with log - still valid if we - // want subset - } - if (typeid(T) != *cpp_type) { - std::cout << "Warning: dataset and array type mismatch" - << std::endl; - } - std::array array_shape{}; - std::transform( - shape.begin(), shape.end(), array_shape.begin(), - [](const auto dim) { return static_cast(dim); }); - - aare::NDArray dataset_array(array_shape); - - read_into_buffer(reinterpret_cast(dataset_array.data())); - - return dataset_array; - } - - // getMemDataSize() - - private: - std::string datasetname{}; - H5::DataSet dataset; - H5::DataSpace dataspace; - H5::DataType datatype; - const std::type_info *cpp_type; - ssize_t rank{}; - std::vector shape{}; -}; - -class HDF5FileReader { - - public: - HDF5FileReader() = default; - - void open_file(const std::string &filename_) { - filename = filename_; - try { - file = H5::H5File(filename, H5F_ACC_RDONLY); - } catch (H5::Exception &e) { - std::cerr << "Error: " << e.getDetailMsg() << std::endl; - } - } - - void close_file() { file.close(); } - - HDF5Dataset get_dataset(const std::string &dataset_name) const { - H5::DataSet dataset; - try { - dataset = file.openDataSet(dataset_name); - } catch (H5::Exception &e) { - std::cerr << "Error: " << e.getDetailMsg() << std::endl; - } - - // TODO use optional to handle error - return HDF5Dataset(dataset_name, dataset); - } - - private: - std::string filename{}; - H5::H5File file; -}; - -} // namespace aare \ No newline at end of file diff --git a/include/aare/MythenDetectorSpecifications.hpp b/include/aare/MythenDetectorSpecifications.hpp deleted file mode 100644 index 0466833..0000000 --- a/include/aare/MythenDetectorSpecifications.hpp +++ /dev/null @@ -1,157 +0,0 @@ -#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 num_counters = 1, - double bloffset = 1.532) - : max_modules_(max_modules), num_counters_(num_counters), - 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); - } - - // TODO templated on filereader - 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; - } - } - - // TODO template on filereader - 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_; } - - size_t num_counters() const { return num_counters_; } - - 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; - - size_t num_counters_ = 1; - - 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 deleted file mode 100644 index 6d78da3..0000000 --- a/include/aare/MythenFileReader.hpp +++ /dev/null @@ -1,82 +0,0 @@ -/************************************************ - * @file MythenFileReader.hpp - * @short minimal file reader to read mythen files - ***********************************************/ - -#include -#include -#include - -#include "Hdf5FileReader.hpp" -#include "NDArray.hpp" - -namespace aare { - -struct MythenFrame { - NDArray photon_counts; - double detector_angle{}; - // double reference_intensity{}; not needed - std::array channel_mask{}; -}; - -/** minimal version for a mythen file reader */ -class MythenFileReader : public HDF5FileReader { - - public: - MythenFileReader(const std::filesystem::path &file_path_, - const std::string &file_prefix_) - : m_base_path(file_path_), file_prefix(file_prefix_) {}; - - MythenFrame read_frame(ssize_t frame_index) { - // TODO not a good design fixed number of digits in file name for frame - // number -> pad with zeros - // not even sure if files have the same name - std::string current_file_name = - m_base_path / (file_prefix + std::to_string(frame_index) + ".h5"); - - MythenFrame frame; - open_file(current_file_name); - - auto dataset_photon_count = - get_dataset("/entry/instrument/detector/data"); - - frame.photon_counts = - dataset_photon_count.store_as_ndarray(); - - ++frame.photon_counts; // Why though? - - auto dataset_detector_angle = - get_dataset("/entry/instrument/NDAttributes/DetectorAngle"); - - dataset_detector_angle.read_into_buffer( - reinterpret_cast(&frame.detector_angle)); - - auto dataset_channel_number = - get_dataset("/entry/instrument/NDAttributes/CounterMask"); - - uint8_t channel_number; - - dataset_channel_number.read_into_buffer( - reinterpret_cast(&channel_number)); - - std::bitset<3> binary_channel_numbers(channel_number); // 1 0 0 - - // binary_channel_numbers.flip(); // TODO not sure where most - // significant - // bit is ask Anna again - - frame.channel_mask = std::array{binary_channel_numbers[0], - binary_channel_numbers[1], - binary_channel_numbers[2]}; - - close_file(); - - return frame; - } - - private: - std::filesystem::path m_base_path{}; - std::string file_prefix{}; -}; - -} // namespace aare \ No newline at end of file diff --git a/src/AngleCalibration.cpp b/src/AngleCalibration.cpp deleted file mode 100644 index 632b8d5..0000000 --- a/src/AngleCalibration.cpp +++ /dev/null @@ -1,377 +0,0 @@ -#include "aare/AngleCalibration.hpp" - -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) { - - std::string line; - uint32_t module_number{}; - - try { - std::ifstream file(filename, std::ios_base::in); - if (!file.good()) { - throw std::logic_error("file does not exist"); - } - - std::stringstream file_buffer; - file_buffer << file.rdbuf(); - - while (file_buffer >> line) { - if (line == "module") { - file_buffer >> line; - module_number = std::stoi(line); - } - if (line == "center") { - file_buffer >> line; - centers.insert(centers.begin() + module_number, - std::stod(line)); - } - if (line == "conversion") { - file_buffer >> line; - conversions.insert(conversions.begin() + module_number, - std::stod(line)); - } - if (line == "offset") { - file_buffer >> line; - offsets.insert(offsets.begin() + module_number, - std::stod(line)); - } - } - - file.close(); - } catch (const std::exception &e) { - std::cerr << "Error: " << e.what() << std::endl; - } -} - -parameters AngleCalibration::convert_to_EE_parameters() const { - - // normal distance between sample and detector (R) - std::vector normal_distances(centers.size()); - // distances between intersection point of sample normal and module origin - // (D) - std::vector module_center_distances(centers.size()); - // angles between undiffracted beam and orthogonal sample projection on - // detector (phi) - std::vector angles(centers.size()); - - for (size_t i = 0; i < centers.size(); ++i) { - auto [module_center_distance, normal_distance, angle] = - convert_to_EE_parameters(i); - normal_distances[i] = normal_distance; - module_center_distances[i] = module_center_distance; - angles[i] = angle; - } - - return std::make_tuple(module_center_distances, normal_distances, angles); -} - -std::tuple -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 { - const double module_center_distance = - center * MythenDetectorSpecifications::pitch(); - const double normal_distance = - MythenDetectorSpecifications::pitch() / std::abs(conversion); - const double angle = offset + 180.0 / M_PI * center * std::abs(conversion); - - return std::make_tuple(module_center_distance, normal_distance, angle); -} - -size_t AngleCalibration::global_to_local_strip_index_conversion( - const size_t global_strip_index) const { - const size_t module_index = - global_strip_index / MythenDetectorSpecifications::strips_per_module(); - // local strip index in module - size_t local_strip_index = - global_strip_index - - module_index * MythenDetectorSpecifications::strips_per_module(); - // switch if indexing is in clock-wise direction - local_strip_index = - std::signbit(conversions[module_index]) - ? MythenDetectorSpecifications::strips_per_module() - 1 - - local_strip_index - : local_strip_index; - - return local_strip_index; -} - -/* -parameters -AngleCalibration::convert_to_BC_parameters() {} -*/ - -double AngleCalibration::diffraction_angle_from_DG_parameters( - const double center, const double conversion, const double offset, - const size_t strip_index, const double distance_to_strip) const { - - return offset + 180.0 / M_PI * - (center * std::abs(conversion) - - atan((center - (strip_index + distance_to_strip)) * - std::abs(conversion))); -} - -double AngleCalibration::diffraction_angle_from_EE_parameters( - const double module_center_distance, const double normal_distance, - const double angle, const size_t strip_index, - const double distance_to_strip) const { - - return angle - 180.0 / M_PI * - atan((module_center_distance - - MythenDetectorSpecifications::pitch() * - (strip_index + distance_to_strip)) / - normal_distance); // TODO: why is it minus - // is it defined counter - // clockwise? thought - // should have a flipped - // sign -} - -double AngleCalibration::angular_strip_width_from_DG_parameters( - const double center, const double conversion, const double offset, - const size_t local_strip_index) const { - - return std::abs(diffraction_angle_from_DG_parameters( - center, conversion, offset, local_strip_index, -0.5) - - diffraction_angle_from_DG_parameters( - center, conversion, offset, local_strip_index, 0.5)); -} - -double AngleCalibration::angular_strip_width_from_EE_parameters( - const double module_center_distance, const double normal_distance, - const double angle, const size_t local_strip_index) const { - - return std::abs(diffraction_angle_from_EE_parameters( - module_center_distance, normal_distance, angle, - local_strip_index, -0.5) - - diffraction_angle_from_EE_parameters( - module_center_distance, normal_distance, angle, - local_strip_index, 0.5)); - - // TODO: again not sure about division order - taking abs anyway -} - -void AngleCalibration::calculate_fixed_bin_angle_width_histogram( - const size_t start_frame_index, const size_t end_frame_index) { - - new_photon_counts = NDArray(std::array{num_bins}); - - 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 = - flat_field->inverse_normalized_flatfield(); - - for (size_t frame_index = start_frame_index; frame_index < end_frame_index; - ++frame_index) { - MythenFrame frame = mythen_file_reader->read_frame(frame_index); - redistribute_photon_counts_to_fixed_angle_bins( - frame, bin_counts.view(), new_statistical_weights.view(), - new_errors.view(), inverse_normalized_flatfield.view()); - } - - for (ssize_t i = 0; i < new_photon_counts.size(); ++i) { - new_photon_counts[i] = (new_statistical_weights[i] <= - std::numeric_limits::epsilon()) - ? 0. - : bin_counts[i] / new_statistical_weights[i]; - new_photon_count_errors[i] = - (bin_counts[i] <= std::numeric_limits::epsilon()) - ? 0. - : 1.0 / std::sqrt(bin_counts[i]); - } -} - -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) const { - - ssize_t channel = 0; // TODO handle mask - FlatField still 1d - - if (frame.photon_counts.shape()[0] != mythen_detector->num_strips()) { - throw std::runtime_error("wrong number of strips read"); - } - - ssize_t num_bins1 = mythen_detector->min_angle() / histogram_bin_width; - ssize_t num_bins2 = mythen_detector->max_angle() / histogram_bin_width; - - std::cout << "position: " << frame.detector_angle - << std::endl; // replace with log - - double exposure_rate = 1. / mythen_detector->exposure_time(); - - for (ssize_t strip_index = 0; strip_index < mythen_detector->num_strips(); - ++strip_index) { - - size_t module_index = - strip_index / MythenDetectorSpecifications::strips_per_module(); - - if (mythen_detector->get_bad_channels()[strip_index] || - !mythen_detector->get_connected_modules()[module_index]) - continue; - - double poisson_error = std::sqrt(frame.photon_counts(strip_index)) * - inverse_normalized_flatfield(strip_index) * - exposure_rate; - - double corrected_photon_count = - frame.photon_counts(strip_index) * - inverse_normalized_flatfield(strip_index) * exposure_rate; - - size_t local_strip_index = - global_to_local_strip_index_conversion(strip_index); - - double diffraction_angle = diffraction_angle_from_DG_parameters( - centers[module_index], conversions[module_index], - offsets[module_index], local_strip_index); - - diffraction_angle += (frame.detector_angle + mythen_detector->dtt0() + - mythen_detector->bloffset()); - - if (diffraction_angle < mythen_detector->min_angle() || - diffraction_angle > mythen_detector->max_angle()) - continue; - - double angle_covered_by_strip = angular_strip_width_from_DG_parameters( - centers[module_index], conversions[module_index], - offsets[module_index], local_strip_index); - - double photon_count_per_bin = histogram_bin_width * - corrected_photon_count / - angle_covered_by_strip; - double error_photon_count_per_bin = - histogram_bin_width * poisson_error / angle_covered_by_strip; - - double statistical_weights = - 1.0 / std::pow(error_photon_count_per_bin, 2); // 1./sigma² - - double strip_boundary_left = - diffraction_angle - 0.5 * angle_covered_by_strip; - double strip_boundary_right = - diffraction_angle + 0.5 * angle_covered_by_strip; - - ssize_t left_bin_index = std::max( - num_bins1, - static_cast( - std::floor(strip_boundary_left / histogram_bin_width) - 1)); - ssize_t right_bin_index = std::min( - num_bins2, - static_cast( - std::ceil(strip_boundary_right / histogram_bin_width) + 1)); - - // TODO should it be < or <= - for (ssize_t bin = left_bin_index; bin <= right_bin_index; ++bin) { - double bin_coverage = std::min(strip_boundary_right, - (bin + 0.5) * histogram_bin_width) - - std::max(strip_boundary_left, - (bin - 0.5) * histogram_bin_width); - - double bin_coverage_factor = bin_coverage / histogram_bin_width; - - ssize_t bin_index = bin - num_bins1; - // TODO: maybe have this threshold configurable - if (bin_coverage >= 0.0001) { - new_statistical_weights(bin_index) += - statistical_weights * bin_coverage_factor; - bin_counts(bin_index) += statistical_weights * - bin_coverage_factor * - photon_count_per_bin; - new_errors(bin_index) += statistical_weights * - bin_coverage_factor * - std::pow(photon_count_per_bin, 2); - } - } - } -} - -void AngleCalibration::write_to_file( - const std::string &filename, const bool store_nonzero_bins, - const std::filesystem::path &filepath) const { - - std::ofstream output_file(filepath / filename); - - if (!output_file) { - std::cerr << "Error opening file!" - << std::endl; // TODO: replace with log - } - - output_file << std::fixed << std::setprecision(15); - - for (ssize_t i = 0; i < num_bins; ++i) { - if (new_photon_counts[i] <= std::numeric_limits::epsilon() && - store_nonzero_bins) { - continue; - } - - output_file << std::floor(mythen_detector->min_angle() / - histogram_bin_width) * - histogram_bin_width + - i * histogram_bin_width - << " " << new_photon_counts[i] << " " - << new_photon_count_errors[i] << std::endl; - } - output_file.close(); -} - -} // namespace aare diff --git a/src/AngleCalibration.test.cpp b/src/AngleCalibration.test.cpp deleted file mode 100644 index dc455e3..0000000 --- a/src/AngleCalibration.test.cpp +++ /dev/null @@ -1,265 +0,0 @@ -/************************************************ - * @file AngleCalibration.test.cpp - * @short test case for angle calibration class - ***********************************************/ - -#include "aare/AngleCalibration.hpp" -#include "aare/CustomFiles.hpp" -#include "aare/FlatField.hpp" - -#include - -#include "test_config.hpp" - -#include -#include - -#include -#include -#include - -using namespace aare; - -TEST_CASE("read initial angle calibration file", - "[anglecalibration] [.files]") { - - std::shared_ptr mythen_detector_ptr = - std::make_shared(); - - AngleCalibration anglecalibration(mythen_detector_ptr, - std::shared_ptr{}, - std::shared_ptr{}); - - std::string filename = test_data_path() / "AngleCalibration_Test_Data" / - "Angcal_2E_Feb2023_P29.off"; - - REQUIRE(std::filesystem::exists(filename)); - - anglecalibration.read_initial_calibration_from_file(filename); - - auto centers = anglecalibration.get_centers(); - auto conversions = anglecalibration.get_conversions(); - auto offsets = anglecalibration.get_offsets(); - - std::cout.precision(17); - - CHECK(centers.size() == 48); - CHECK(conversions.size() == 48); - CHECK(offsets.size() == 48); - - CHECK(centers[9] == Catch::Approx(660.342326)); - CHECK(offsets[47] == Catch::Approx(5.8053312)); - CHECK(conversions[27] == Catch::Approx(-0.6581179125e-4)); -} - -TEST_CASE("read bad channels", - "[anglecalibration][mythenspecifications][.files]") { - - MythenDetectorSpecifications mythen_detector; - - std::string bad_channels_filename = test_data_path() / - "AngleCalibration_Test_Data" / - "bc2023_003_RING.chans"; - - REQUIRE(std::filesystem::exists(bad_channels_filename)); - - mythen_detector.read_bad_channels_from_file(bad_channels_filename); - - CHECK(mythen_detector.get_bad_channels().size() == 61440); - - CHECK(mythen_detector.get_bad_channels()[61437] == true); - CHECK(std::all_of(mythen_detector.get_bad_channels().begin() + 30720, - mythen_detector.get_bad_channels().begin() + 61439, - [](const bool element) { return element; })); -} - -TEST_CASE("read unconnected modules", - "[anglecalibration][mythenspecifications][.files]") { - - MythenDetectorSpecifications mythen_detector; - - std::string unconnected_modules_filename = - test_data_path() / "AngleCalibration_Test_Data" / "ModOut.txt"; - - REQUIRE(std::filesystem::exists(unconnected_modules_filename)); - - mythen_detector.read_unconnected_modules_from_file( - unconnected_modules_filename); - - CHECK(mythen_detector.get_connected_modules().size() == 48); - - CHECK(std::all_of(mythen_detector.get_connected_modules().begin(), - mythen_detector.get_connected_modules().end(), - [](const bool element) { return element; })); -} - -TEST_CASE("read flatfield", "[anglecalibration][flatfield][.files]") { - - std::shared_ptr mythen_detector_ptr = - std::make_shared(); - - FlatField flatfield(mythen_detector_ptr); - - std::string flatfield_filename = - test_data_path() / "AngleCalibration_Test_Data" / - "Flatfield_E22p0keV_T11000eV_up_48M_a_LONG_Feb2023_open_WS_SUMC.raw"; - - REQUIRE(std::filesystem::exists(flatfield_filename)); - - flatfield.read_flatfield_from_file(flatfield_filename); - - auto flatfield_data = flatfield.get_flatfield(); - - CHECK(flatfield_data.size() == 61440); - - CHECK(flatfield_data[0] == 0); - CHECK(flatfield_data[21] == 4234186); -} - -TEST_CASE("create flatfield", "[anglecalibration], [flatfield], [.files]") { - - ssize_t n_modules = 2; - double exposure_time = 5.0; - size_t n_counters = 1; - std::shared_ptr mythen_detector_ptr = - std::make_shared(n_modules, exposure_time, - n_counters); - - FlatField flatfield(mythen_detector_ptr); - - auto data_path = test_data_path() / "AngleCalibration_Test_Data" / - "Flatfieldacquisition"; - - REQUIRE(std::filesystem::exists(data_path)); - - CHECK_NOTHROW(flatfield.create_flatfield_from_rawfilesystem(data_path)); - - auto flatfield_data = flatfield.get_flatfield(); - - CHECK(flatfield_data.size() == n_modules * 1280); - - CHECK(flatfield_data[0] == 0); - CHECK(flatfield_data[21] == - 2 * 3 * - 21); // virtual data 2 angles, 3 frames - 21 as increasing numbers -} - -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 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); - - std::string unconnected_modules_filename = fpath / "ModOut.txt"; - - REQUIRE(std::filesystem::exists(unconnected_modules_filename)); - - mythen_detector_ptr->read_unconnected_modules_from_file( - unconnected_modules_filename); - - std::shared_ptr flat_field_ptr = - std::make_shared(mythen_detector_ptr); - - std::string flatfield_filename = - fpath / - "Flatfield_E22p0keV_T11000eV_up_48M_a_LONG_Feb2023_open_WS_SUMC.raw"; - - REQUIRE(std::filesystem::exists(flatfield_filename)); - - flat_field_ptr->read_flatfield_from_file( - flatfield_filename); - - std::shared_ptr 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); - - std::string initial_angles_filename = fpath / "Angcal_2E_Feb2023_P29.off"; - - REQUIRE(std::filesystem::exists(initial_angles_filename)); - - anglecalibration.read_initial_calibration_from_file( - initial_angles_filename); - - anglecalibration.calculate_fixed_bin_angle_width_histogram(320, 340); - - // anglecalibration.write_to_file("cpp_new_photon_counts.xye"); - - auto expected_filename_photons = - test_data_path() / "AngleCalibration_Test_Data" / "new_photons.bin"; - - REQUIRE(std::filesystem::exists(expected_filename_photons)); - - auto expected_filename_errors = - test_data_path() / "AngleCalibration_Test_Data" / "new_errors.bin"; - - REQUIRE(std::filesystem::exists(expected_filename_errors)); - - ssize_t new_num_bins = anglecalibration.get_new_num_bins(); - - auto python_output_errors = load( - 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("check conversion from DG to EE parameters", "[anglecalibration]") { - - std::shared_ptr mythen_detector_ptr = - std::make_shared(); - - AngleCalibration anglecalibration(mythen_detector_ptr, - std::shared_ptr{}, - std::shared_ptr{}); - - // DG test parameters - const double center = 642.197591224993; - const double conversion = 0.657694036246975e-4; - const double offset = 5.004892881251670; - const ssize_t local_strip_index = 1; - - double diffraction_angle_DG_param = - anglecalibration.diffraction_angle_from_DG_parameters( - center, conversion, offset, local_strip_index); - - auto [distance_center, normal_distance, angle] = - anglecalibration.convert_to_EE_parameters(center, conversion, offset); - - double diffraction_angle_EE_param = - anglecalibration.diffraction_angle_from_EE_parameters( - distance_center, normal_distance, angle, local_strip_index); - - CHECK(diffraction_angle_EE_param == - Catch::Approx(diffraction_angle_DG_param)); - - double strip_width_DG_param = - anglecalibration.angular_strip_width_from_DG_parameters( - center, conversion, offset, local_strip_index); - - double strip_width_EE_param = - anglecalibration.angular_strip_width_from_EE_parameters( - distance_center, normal_distance, angle, local_strip_index); - - CHECK(strip_width_DG_param == Catch::Approx(strip_width_EE_param)); -} diff --git a/src/CustomFiles.cpp b/src/CustomFiles.cpp deleted file mode 100644 index a14b987..0000000 --- a/src/CustomFiles.cpp +++ /dev/null @@ -1,89 +0,0 @@ -#include "aare/CustomFiles.hpp" - -namespace aare { - -CustomMythenFile::CustomMythenFile(const std::string &filename, - const ssize_t rows, const ssize_t cols, - const std::string &mode) - : m_filename(filename), m_rows(rows), m_cols(cols) { - - m_mode = mode; - if (m_mode == "r") { - try { - m_file.open(m_filename, std::ios_base::in); - if (!m_file.good()) { - throw std::logic_error("file does not exist"); - } - } catch (std::exception &e) { - - std::cerr << "Error: " << e.what() - << std::endl; // TODO replace with log - } - - } else { - throw std::runtime_error(LOCATION + - "Unsupported mode. Can only read RawFiles."); - } -} - -CustomMythenFile::~CustomMythenFile() { m_file.close(); } - -Frame CustomMythenFile::read_frame() { - auto f = Frame(m_rows, m_cols, m_dtype); - uint32_t *frame_buffer = reinterpret_cast(f.data()); - uint32_t strip_index, photon_count; - while (m_file >> strip_index >> photon_count) { - *frame_buffer = photon_count; - ++frame_buffer; - } - return f; -} - -void CustomMythenFile::read_into(std::byte *image_buf) { - uint32_t strip_index, photon_count; - while (m_file >> strip_index >> photon_count) { - std::memcpy(image_buf, &photon_count, sizeof(photon_count)); - image_buf += sizeof(photon_count); - } -} - -size_t CustomMythenFile::bytes_per_frame() { - return m_num_strips * m_dtype.bytes(); // TODO do i want m_counts? -} - -Frame CustomMythenFile::read_frame(size_t frame_number) { - return read_frame(); // maybe give count as frame_number -} - -std::vector CustomMythenFile::read_n(size_t n_frames) { - std::vector vec; - vec.reserve(1); - vec.push_back(read_frame()); - return vec; // std::vector{read_frame()}; -} - -void CustomMythenFile::read_into(std::byte *image_buf, size_t n_frames) { - read_into(image_buf); -} - -size_t CustomMythenFile::frame_number(size_t frame_index) { return 1; } - -size_t CustomMythenFile::pixels_per_frame() { return m_rows * m_cols; } - -void CustomMythenFile::seek(size_t frame_number) {} - -size_t CustomMythenFile::tell() { return 1; } - -size_t CustomMythenFile::total_frames() const { return 1; } - -size_t CustomMythenFile::rows() const { return m_rows; } - -size_t CustomMythenFile::cols() const { return m_cols; } - -size_t CustomMythenFile::bitdepth() const { return m_dtype.bitdepth(); } - -DetectorType CustomMythenFile::detector_type() const { return det_type; } - -const Dtype CustomMythenFile::m_dtype(Dtype::TypeIndex::UINT32); - -} // namespace aare \ No newline at end of file diff --git a/src/Hdf5FileReader.test.cpp b/src/Hdf5FileReader.test.cpp deleted file mode 100644 index d9659db..0000000 --- a/src/Hdf5FileReader.test.cpp +++ /dev/null @@ -1,110 +0,0 @@ -/************************************************ - * @file Hdf5FileReader.test.cpp - * @short test case for reading hdf5 files - ***********************************************/ - -#include - -#include "test_config.hpp" - -#include - -#include "aare/Hdf5FileReader.hpp" -#include "aare/NDArray.hpp" - -#include -#include -#include - -using namespace aare; - -TEST_CASE("read hdf5 file", "[.hdf5file][.files]") { - - // TODO generalize datasetpath - std::string filename = test_data_path() / "AngleCalibration_Test_Data" / - "ang1up_22keV_LaB60p3mm_48M_a_0320.h5"; - - REQUIRE(std::filesystem::exists(filename)); - - HDF5FileReader file_reader; - - file_reader.open_file(filename); - - auto dataset = file_reader.get_dataset("/entry/data/data"); - - auto shape = dataset.get_shape(); - - CHECK(shape[0] == 61440); - - auto type = dataset.get_datatype(); - - const std::type_info *type_info = dataset.get_cpp_type(); - - CHECK(*type_info == typeid(uint32_t)); - - SECTION("read dataset into NDArray") { - - NDArray dataset_array = - dataset.store_as_ndarray(); - - CHECK(dataset_array(0) == 866); - CHECK(dataset_array(61439) == 1436); - } - - SECTION("read dataset into Frame") { - Frame frame = dataset.store_as_frame(); - CHECK(*(reinterpret_cast(frame.pixel_ptr(0, 0))) == 866); - CHECK(*(reinterpret_cast(frame.pixel_ptr(0, 61439))) == - 1436); - } - SECTION("read subset of dataset") { - Frame frame(1, 10, Dtype(typeid(uint32_t))); - - Subset subset{std::vector{10}, std::vector{10}}; - - dataset.read_into_buffer(frame.data(), subset); - - CHECK(*(reinterpret_cast(frame.pixel_ptr(0, 0))) == 664); - CHECK(*(reinterpret_cast(frame.pixel_ptr(0, 9))) == 654); - } - /* - SECTION("read scalar") { - } - */ -} - -TEST_CASE("test datatypes", "[.hdf5file]") { - - auto [dtype, expected_type_info] = GENERATE( - std::make_tuple(H5::DataType(H5::PredType::NATIVE_INT), &typeid(int)), - std::make_tuple(H5::DataType(H5::PredType::NATIVE_INT8), - &typeid(int8_t)), - std::make_tuple(H5::DataType(H5::PredType::NATIVE_UINT16), - &typeid(uint16_t)), - std::make_tuple(H5::DataType(H5::PredType::NATIVE_INT16), - &typeid(int16_t)), - std::make_tuple(H5::DataType(H5::PredType::STD_U32LE), - &typeid(uint32_t)), - std::make_tuple(H5::DataType(H5::PredType::STD_I32LE), - &typeid(int32_t)), - std::make_tuple(H5::DataType(H5::PredType::NATIVE_INT32), - &typeid(int32_t)), - std::make_tuple(H5::DataType(H5::PredType::IEEE_F64LE), - &typeid(double)), - std::make_tuple(H5::DataType(H5::PredType::IEEE_F32LE), &typeid(float)), - std::make_tuple(H5::DataType(H5::PredType::NATIVE_FLOAT), - &typeid(float)), - std::make_tuple(H5::DataType(H5::PredType::NATIVE_DOUBLE), - &typeid(double)), - std::make_tuple(H5::DataType(H5::PredType::NATIVE_CHAR), - &typeid(int8_t))); - - const std::type_info &type_info = deduce_cpp_type(dtype); - - CHECK(type_info == *expected_type_info); - - // TODO: handle bit swapping - REQUIRE_THROWS(deduce_cpp_type( - H5::DataType(H5::PredType::IEEE_F32BE))); // does not convert from big - // to little endian -} diff --git a/src/MythenFileReader.test.cpp b/src/MythenFileReader.test.cpp deleted file mode 100644 index a5d0a46..0000000 --- a/src/MythenFileReader.test.cpp +++ /dev/null @@ -1,33 +0,0 @@ -/************************************************ - * @file MythenFileReader.test.cpp - * @short test case for angle calibration class - ***********************************************/ - -#include "aare/MythenFileReader.hpp" - -#include - -#include "test_config.hpp" - -#include -#include -#include - -using namespace aare; - -TEST_CASE("test mythenfile_reader", "[.mythenfilereader][.files]") { - - auto fpath = test_data_path() / "AngleCalibration_Test_Data"; - - REQUIRE(std::filesystem::exists(fpath)); - - MythenFileReader file_reader(fpath, "ang1up_22keV_LaB60p3mm_48M_a_0"); - - auto frame = file_reader.read_frame(320); - - CHECK(frame.detector_angle == 0.99955); - - CHECK(frame.channel_mask == std::array{0, 0, 1}); - - CHECK(frame.photon_counts.size() == 61440); -}