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