diff --git a/CMakeLists.txt b/CMakeLists.txt index 1149975..6a726fa 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -471,6 +471,7 @@ 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/include/aare/AngleCalibration.hpp b/include/aare/AngleCalibration.hpp index f97bd93..65dcabe 100644 --- a/include/aare/AngleCalibration.hpp +++ b/include/aare/AngleCalibration.hpp @@ -1,6 +1,5 @@ #pragma once #include -#include #include #include @@ -162,13 +161,13 @@ class FlatField { : mythen_detector(mythen_detector_) { flat_field = NDArray( - std::array{mythen_detector->num_strips()}); + std::array{mythen_detector->num_strips()}, 0); } void read_flatfield_from_file(const std::string &filename) { std::string word; - uint32_t module_number{}; + uint32_t strip_number{}; try { std::ifstream file(filename, std::ios_base::in); @@ -181,11 +180,11 @@ class FlatField { while (file_buffer >> word) { - module_number = std::stoi(word); + strip_number = std::stoi(word); file_buffer >> word; - if (!mythen_detector->get_bad_channels()[module_number]) - flat_field[module_number] = std::stod(word); + if (!mythen_detector->get_bad_channels()[strip_number]) + flat_field[strip_number] = std::stod(word); } file.close(); @@ -194,6 +193,8 @@ class FlatField { } } + 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(), @@ -212,9 +213,20 @@ class FlatField { 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; }); + [&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] == 0 ? 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 @@ -228,14 +240,20 @@ class FlatField { 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; } - // TODO: update is bad channels private: NDArray flat_field; // TODO: should be 2d std::shared_ptr mythen_detector; @@ -277,7 +295,8 @@ class AngleCalibration { /** reads the historical Detector Group (DG) parameters from file **/ void read_initial_calibration_from_file(const std::string &filename); - /** converts DG parameters to easy EE parameters e.g.geometric parameters */ + /** converts DG parameters to easy EE parameters e.g.geometric + * parameters */ parameters convert_to_EE_parameters(); std::tuple @@ -287,16 +306,16 @@ class AngleCalibration { * parameters */ parameters convert_to_BC_parameters(); - /** calculates diffraction angle from EE module parameters (used in Beer's - * Law) + /** 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 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) + /** 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, @@ -324,8 +343,8 @@ class AngleCalibration { const size_t end_frame_index); /** - * redistributes photon counts with of histogram using one bin per strip to - * histogram with fixed size angle bins + * 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 @@ -345,15 +364,15 @@ class AngleCalibration { // 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 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 - 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)? + 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; diff --git a/include/aare/MythenFileReader.hpp b/include/aare/MythenFileReader.hpp index 1859fcb..0811f76 100644 --- a/include/aare/MythenFileReader.hpp +++ b/include/aare/MythenFileReader.hpp @@ -3,6 +3,7 @@ * @short minimal file reader to read mythen files ***********************************************/ +#include #include #include @@ -15,7 +16,7 @@ struct MythenFrame { NDArray photon_counts; double detector_angle{}; // double reference_intensity{}; not needed - std::array channel_mask{}; + std::array channel_mask{}; }; /** minimal version for a mythen file reader */ @@ -30,8 +31,8 @@ class MythenFileReader : public HDF5FileReader { // 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.string() + file_prefix + - std::to_string(frame_index) + ".h5"; + std::string current_file_name = + m_base_path / (file_prefix + std::to_string(frame_index) + ".h5"); open_file(current_file_name); @@ -63,9 +64,9 @@ class MythenFileReader : public HDF5FileReader { // significant // bit is ask Anna again - std::array channel_mask{binary_channel_numbers[2], - binary_channel_numbers[1], - binary_channel_numbers[0]}; + std::array channel_mask{binary_channel_numbers[0], + binary_channel_numbers[1], + binary_channel_numbers[2]}; close_file(); diff --git a/src/AngleCalibration.cpp b/src/AngleCalibration.cpp index 3dcef71..c447294 100644 --- a/src/AngleCalibration.cpp +++ b/src/AngleCalibration.cpp @@ -283,8 +283,11 @@ void AngleCalibration::write_to_file(const std::string &filename) { output_file << std::fixed << std::setprecision(6); for (ssize_t i = 0; i < num_bins; ++i) { - output_file << i * histogram_bin_width << " " << new_photon_counts[i] - << " " << new_photon_count_errors[i] << std::endl; + if (new_photon_counts[i] == 0) + continue; + output_file << mythen_detector->min_angle() + i * histogram_bin_width + << " " << new_photon_counts[i] << " " + << new_photon_count_errors[i] << std::endl; } output_file.close(); } diff --git a/src/AngleCalibration.test.cpp b/src/AngleCalibration.test.cpp index 922541e..3ba1aa2 100644 --- a/src/AngleCalibration.test.cpp +++ b/src/AngleCalibration.test.cpp @@ -7,21 +7,23 @@ #include +#include "test_config.hpp" + #include #include #include using namespace aare; -class AngleCalibrationTestClass : public aare::AngleCalibration { +class AngleCalibrationTestClass : public AngleCalibration { public: AngleCalibrationTestClass( std::shared_ptr mythen_detector_, std::shared_ptr flat_field_, std::shared_ptr mythen_file_reader_) - : aare::AngleCalibration(mythen_detector_, flat_field_, - mythen_file_reader_) {} + : AngleCalibration(mythen_detector_, flat_field_, mythen_file_reader_) { + } ~AngleCalibrationTestClass() = default; std::vector get_centers() { return centers; } @@ -32,7 +34,11 @@ class AngleCalibrationTestClass : public aare::AngleCalibration { }; TEST_CASE("read initial angle calibration file", - "[.anglecalibration][.fileread]") { + "[.anglecalibration] [.files]") { + + auto fpath = test_data_path() / "AngleCalibration_Test_Data"; + + REQUIRE(std::filesystem::exists(fpath)); std::shared_ptr mythen_detector_ptr = std::make_shared(); @@ -41,17 +47,13 @@ TEST_CASE("read initial angle calibration file", std::make_shared(mythen_detector_ptr); std::shared_ptr mythen_file_reader_ptr = - std::make_shared( - std::filesystem::path{ - "/home/mazzola/Documents/mythen3tools/beamline/TDATA"}, - "ang1up_22keV_LaB60p3mm_48M_a_0"); + 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 = - "/home/mazzol_a/Documents/mythen3tools/beamline/" - "TDATA/Angcal_2E_Feb2023_P29.off"; // TODO change path upload data + std::string filename = fpath / "Angcal_2E_Feb2023_P29.off"; REQUIRE(std::filesystem::exists(filename)); @@ -71,3 +73,127 @@ TEST_CASE("read initial angle calibration file", CHECK(offsets[47] == Catch::Approx(5.8053312)); CHECK(conversions[27] == Catch::Approx(-0.6581179125e-4)); } + +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"; + + 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]") { + auto fpath = test_data_path() / "AngleCalibration_Test_Data"; + + REQUIRE(std::filesystem::exists(fpath)); + + MythenDetectorSpecifications mythen_detector; + + std::string unconnected_modules_filename = fpath / "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]") { + auto fpath = test_data_path() / "AngleCalibration_Test_Data"; + + REQUIRE(std::filesystem::exists(fpath)); + + std::shared_ptr mythen_detector_ptr = + std::make_shared(); + + 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(); + + CHECK(flatfield_data.size() == 61440); + + CHECK(flatfield_data[0] == 0); + CHECK(flatfield_data[21] == 4234186); +} + +TEST_CASE("calculate new fixed angle width bins histogram", + "[.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"); // TODO adjust output path +} diff --git a/src/Hdf5FileReader.test.cpp b/src/Hdf5FileReader.test.cpp index 26f27a8..d9659db 100644 --- a/src/Hdf5FileReader.test.cpp +++ b/src/Hdf5FileReader.test.cpp @@ -5,6 +5,8 @@ #include +#include "test_config.hpp" + #include #include "aare/Hdf5FileReader.hpp" @@ -16,11 +18,11 @@ using namespace aare; -TEST_CASE("read hdf5 file", "[.hdf5file]") { +TEST_CASE("read hdf5 file", "[.hdf5file][.files]") { // TODO generalize datasetpath - std::string filename = "/home/mazzol_a/Documents/mythen3tools/beamline/" - "TDATA/ang1up_22keV_LaB60p3mm_48M_a_0320.h5"; + std::string filename = test_data_path() / "AngleCalibration_Test_Data" / + "ang1up_22keV_LaB60p3mm_48M_a_0320.h5"; REQUIRE(std::filesystem::exists(filename)); diff --git a/src/MythenFileReader.test.cpp b/src/MythenFileReader.test.cpp new file mode 100644 index 0000000..a5d0a46 --- /dev/null +++ b/src/MythenFileReader.test.cpp @@ -0,0 +1,33 @@ +/************************************************ + * @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); +}