mirror of
https://github.com/slsdetectorgroup/aare.git
synced 2025-06-18 18:27:13 +02:00
some more tests and refactoring
This commit is contained in:
@ -69,7 +69,7 @@ class MythenDetectorSpecifications {
|
||||
} 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 <= std::stoi(line.substr(pos + 1, line_size - pos));
|
||||
++i)
|
||||
bad_channels(i) = true;
|
||||
}
|
||||
@ -200,17 +200,21 @@ class FlatField {
|
||||
flat_field.begin(), flat_field.end(),
|
||||
std::make_pair<double, ssize_t>(0.0, 0),
|
||||
[&tolerance](std::pair<double, ssize_t> acc, const auto &element) {
|
||||
return element < tolerance ? acc
|
||||
: std::make_pair(acc.first + element,
|
||||
acc.second + 1);
|
||||
return element == 0 ? acc
|
||||
: std::make_pair(acc.first + element,
|
||||
acc.second + 1);
|
||||
});
|
||||
|
||||
std::cout << "sum: " << sum << std::endl;
|
||||
std::cout << "count: " << count << std::endl;
|
||||
return sum / count;
|
||||
}
|
||||
|
||||
NDArray<double, 1> inverse_normalized_flatfield(double tolerance = 0.001) {
|
||||
double mean = calculate_mean(tolerance);
|
||||
|
||||
std::cout << "mean: " << mean << std::endl;
|
||||
|
||||
NDArray<double, 1> inverse_normalized_flatfield(flat_field.shape());
|
||||
|
||||
/*
|
||||
@ -223,7 +227,7 @@ class FlatField {
|
||||
|
||||
for (ssize_t i = 0; i < flat_field.size(); ++i) {
|
||||
inverse_normalized_flatfield[i] =
|
||||
(flat_field[i] == 0 ? 0.0 : mean / flat_field[i]);
|
||||
(flat_field[i] <= tolerance ? 0.0 : mean / flat_field[i]);
|
||||
if (inverse_normalized_flatfield[i] < tolerance)
|
||||
mythen_detector->get_bad_channels()[i] = true;
|
||||
}
|
||||
@ -274,20 +278,22 @@ class AngleCalibration {
|
||||
|
||||
exposure_rate = 1. / mythen_detector->exposure_time();
|
||||
|
||||
num_bins = mythen_detector->max_angle() / histogram_bin_width -
|
||||
mythen_detector->min_angle() /
|
||||
histogram_bin_width; // TODO only works if negative
|
||||
// and positive angle
|
||||
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
|
||||
}
|
||||
|
||||
/** set the histogram bin width [degrees] */
|
||||
void set_histogram_bin_width(double bin_width) {
|
||||
histogram_bin_width = bin_width;
|
||||
|
||||
num_bins = mythen_detector->max_angle() / histogram_bin_width -
|
||||
mythen_detector->min_angle() /
|
||||
histogram_bin_width; // TODO only works if negative
|
||||
// and positive angle
|
||||
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() { return histogram_bin_width; }
|
||||
@ -295,6 +301,12 @@ class AngleCalibration {
|
||||
/** reads the historical Detector Group (DG) parameters from file **/
|
||||
void read_initial_calibration_from_file(const std::string &filename);
|
||||
|
||||
std::vector<double> get_centers() { return centers; }
|
||||
|
||||
std::vector<double> get_conversions() { return conversions; }
|
||||
|
||||
std::vector<double> get_offsets() { return offsets; }
|
||||
|
||||
/** converts DG parameters to easy EE parameters e.g.geometric
|
||||
* parameters */
|
||||
parameters convert_to_EE_parameters();
|
||||
@ -302,6 +314,10 @@ class AngleCalibration {
|
||||
std::tuple<double, double, double>
|
||||
convert_to_EE_parameters(const size_t module_index);
|
||||
|
||||
std::tuple<double, double, double>
|
||||
convert_to_EE_parameters(const double center, const double conversion,
|
||||
const double offset);
|
||||
|
||||
/** converts DG parameters to BC parameters e.g. best computing
|
||||
* parameters */
|
||||
parameters convert_to_BC_parameters();
|
||||
@ -311,22 +327,38 @@ class AngleCalibration {
|
||||
* @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);
|
||||
const double module_center_distance, const double normal_distance,
|
||||
const double angle, const size_t strip_index,
|
||||
const double distance_to_strip = 0);
|
||||
|
||||
/** 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);
|
||||
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);
|
||||
|
||||
/** calculated the strip width expressed as angle [degrees]
|
||||
* @param strip_index gloabl strip index of detector
|
||||
*/
|
||||
double angular_strip_width(const size_t strip_index);
|
||||
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_EE_parameters(
|
||||
const double module_center_distance, const double normal_distance,
|
||||
const double angle, const size_t strip_index);
|
||||
|
||||
/** converts global strip index to local strip index of that module */
|
||||
size_t
|
||||
@ -352,8 +384,8 @@ class AngleCalibration {
|
||||
*/
|
||||
void redistribute_photon_counts_to_fixed_angle_bins(
|
||||
const MythenFrame &frame, NDView<double, 1> bin_counts,
|
||||
NDView<double, 1> new_statistical_weights,
|
||||
NDView<double, 1> new_errors);
|
||||
NDView<double, 1> new_statistical_weights, NDView<double, 1> new_errors,
|
||||
NDView<double, 1> inverse_nromalized_flatfield);
|
||||
|
||||
void write_to_file(const std::string &filename);
|
||||
|
||||
@ -385,7 +417,7 @@ class AngleCalibration {
|
||||
|
||||
ssize_t num_bins{};
|
||||
|
||||
double exposure_rate;
|
||||
double exposure_rate{};
|
||||
|
||||
std::shared_ptr<MythenFileReader>
|
||||
mythen_file_reader; // TODO replace by FileInterface ptr
|
||||
|
@ -423,7 +423,7 @@ template <typename T, ssize_t Ndim>
|
||||
void save(NDArray<T, Ndim> &img, std::string &pathname) {
|
||||
std::ofstream f;
|
||||
f.open(pathname, std::ios::binary);
|
||||
f.write(img.buffer(), img.size() * sizeof(T));
|
||||
f.write(reinterpret_cast<char *>(img.buffer()), img.size() * sizeof(T));
|
||||
f.close();
|
||||
}
|
||||
|
||||
@ -433,7 +433,7 @@ NDArray<T, Ndim> load(const std::string &pathname,
|
||||
NDArray<T, Ndim> img{shape};
|
||||
std::ifstream f;
|
||||
f.open(pathname, std::ios::binary);
|
||||
f.read(img.buffer(), img.size() * sizeof(T));
|
||||
f.read(reinterpret_cast<char *>(img.buffer()), img.size() * sizeof(T));
|
||||
f.close();
|
||||
return img;
|
||||
}
|
||||
|
@ -57,28 +57,32 @@ parameters AngleCalibration::convert_to_EE_parameters() {
|
||||
std::vector<double> angles(centers.size());
|
||||
|
||||
for (size_t i = 0; i < centers.size(); ++i) {
|
||||
auto [normal_distance, module_center_distance, angle] =
|
||||
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(normal_distances, module_center_distances, angles);
|
||||
return std::make_tuple(module_center_distances, normal_distances, angles);
|
||||
}
|
||||
|
||||
std::tuple<double, double, double>
|
||||
AngleCalibration::convert_to_EE_parameters(const size_t module_index) {
|
||||
const double normal_distance =
|
||||
centers[module_index] * MythenDetectorSpecifications::pitch();
|
||||
const double module_center_distance =
|
||||
MythenDetectorSpecifications::pitch() /
|
||||
std::abs(conversions[module_index]);
|
||||
const double angle =
|
||||
offsets[module_index] + 180.0 / M_PI * centers[module_index] *
|
||||
std::abs(conversions[module_index]);
|
||||
return convert_to_EE_parameters(centers[module_index],
|
||||
conversions[module_index],
|
||||
offsets[module_index]);
|
||||
}
|
||||
|
||||
return std::make_tuple(normal_distance, module_center_distance, angle);
|
||||
std::tuple<double, double, double> AngleCalibration::convert_to_EE_parameters(
|
||||
const double center, const double conversion, const double offset) {
|
||||
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(
|
||||
@ -92,7 +96,7 @@ size_t AngleCalibration::global_to_local_strip_index_conversion(
|
||||
// switch if indexing is in clock-wise direction
|
||||
local_strip_index =
|
||||
std::signbit(conversions[module_index])
|
||||
? MythenDetectorSpecifications::strips_per_module() -
|
||||
? MythenDetectorSpecifications::strips_per_module() - 1 -
|
||||
local_strip_index
|
||||
: local_strip_index;
|
||||
|
||||
@ -106,45 +110,81 @@ 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 size_t strip_index, const double distance_to_strip) {
|
||||
|
||||
return offset + 180.0 / M_PI *
|
||||
(center * conversion -
|
||||
atan((center - strip_index) * conversion));
|
||||
(center * std::abs(conversion) -
|
||||
atan((center - (strip_index + distance_to_strip)) *
|
||||
std::abs(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) {
|
||||
const double module_center_distance, const double normal_distance,
|
||||
const double angle, const size_t strip_index,
|
||||
const double distance_to_strip) {
|
||||
|
||||
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
|
||||
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(const size_t strip_index) {
|
||||
double AngleCalibration::angular_strip_width_from_DG_parameters(
|
||||
const size_t strip_index) {
|
||||
|
||||
const size_t module_index =
|
||||
strip_index / MythenDetectorSpecifications::strips_per_module();
|
||||
|
||||
const auto [normal_distance, module_center_distance, angle] =
|
||||
convert_to_EE_parameters(module_index);
|
||||
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);
|
||||
|
||||
return 180.0 / M_PI *
|
||||
std::abs(diffraction_angle_from_EE_parameters(
|
||||
normal_distance, module_center_distance, angle,
|
||||
local_strip_index - 0.5) -
|
||||
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 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);
|
||||
|
||||
return std::abs(diffraction_angle_from_EE_parameters(
|
||||
module_center_distance, normal_distance, angle,
|
||||
local_strip_index, -0.5) -
|
||||
diffraction_angle_from_EE_parameters(
|
||||
normal_distance, module_center_distance, angle,
|
||||
local_strip_index + 0.5));
|
||||
module_center_distance, normal_distance, angle,
|
||||
local_strip_index, 0.5));
|
||||
|
||||
// TODO: again not sure about division order - taking abs anyway
|
||||
}
|
||||
|
||||
@ -158,27 +198,37 @@ void AngleCalibration::calculate_fixed_bin_angle_width_histogram(
|
||||
|
||||
NDArray<double, 1> bin_counts(std::array<ssize_t, 1>{num_bins}, 0.0);
|
||||
NDArray<double, 1> new_statistical_weights(std::array<ssize_t, 1>{num_bins},
|
||||
1.0);
|
||||
0.0);
|
||||
|
||||
NDArray<double, 1> new_errors(std::array<ssize_t, 1>{num_bins}, 0.0);
|
||||
|
||||
NDArray<double, 1> 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, new_photon_counts.view(), new_statistical_weights.view(),
|
||||
new_errors.view());
|
||||
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] = bin_counts[i] / new_statistical_weights[i];
|
||||
new_photon_count_errors[i] = 1.0 / std::sqrt(bin_counts[i]);
|
||||
new_photon_counts[i] = (new_statistical_weights[i] <=
|
||||
std::numeric_limits<double>::epsilon())
|
||||
? 0.
|
||||
: bin_counts[i] / new_statistical_weights[i];
|
||||
new_photon_count_errors[i] =
|
||||
(bin_counts[i] <= std::numeric_limits<double>::epsilon())
|
||||
? 0.
|
||||
: 1.0 / std::sqrt(bin_counts[i]);
|
||||
}
|
||||
}
|
||||
|
||||
void AngleCalibration::redistribute_photon_counts_to_fixed_angle_bins(
|
||||
const MythenFrame &frame, NDView<double, 1> bin_counts,
|
||||
NDView<double, 1> new_statistical_weights, NDView<double, 1> new_errors) {
|
||||
NDView<double, 1> new_statistical_weights, NDView<double, 1> new_errors,
|
||||
NDView<double, 1> inverse_normalized_flatfield) {
|
||||
|
||||
ssize_t channel = 0; // TODO handle mask - FlatField still 1d
|
||||
|
||||
@ -186,12 +236,17 @@ void AngleCalibration::redistribute_photon_counts_to_fixed_angle_bins(
|
||||
throw std::runtime_error("wrong number of strips read");
|
||||
}
|
||||
|
||||
NDArray<double, 1> inverse_normalized_flatfield =
|
||||
flat_field->inverse_normalized_flatfield();
|
||||
// NDArray<double, 1> diffraction_angles(
|
||||
// std::array<ssize_t, 1>{mythen_detector->num_strips()}, 0.0);
|
||||
|
||||
// NDArray<double, 1> angle_widths(
|
||||
// std::array<ssize_t, 1>{mythen_detector->num_strips()}, 0.0);
|
||||
|
||||
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;
|
||||
|
||||
for (ssize_t strip_index = 0; strip_index < mythen_detector->num_strips();
|
||||
++strip_index) {
|
||||
|
||||
@ -219,11 +274,16 @@ void AngleCalibration::redistribute_photon_counts_to_fixed_angle_bins(
|
||||
diffraction_angle += (frame.detector_angle + mythen_detector->dtt0() +
|
||||
mythen_detector->bloffset());
|
||||
|
||||
// diffraction_angles[strip_index] = diffraction_angle;
|
||||
|
||||
if (diffraction_angle < mythen_detector->min_angle() ||
|
||||
diffraction_angle > mythen_detector->max_angle())
|
||||
continue;
|
||||
|
||||
double angle_covered_by_strip = angular_strip_width(strip_index);
|
||||
double angle_covered_by_strip =
|
||||
angular_strip_width_from_DG_parameters(strip_index);
|
||||
|
||||
// angle_widths[strip_index] = angle_covered_by_strip;
|
||||
|
||||
double photon_count_per_bin = histogram_bin_width *
|
||||
corrected_photon_count /
|
||||
@ -258,10 +318,10 @@ void AngleCalibration::redistribute_photon_counts_to_fixed_angle_bins(
|
||||
double bin_coverage_factor = bin_coverage / histogram_bin_width;
|
||||
|
||||
ssize_t bin_index = bin - num_bins1;
|
||||
if (bin_coverage > 0.0001) {
|
||||
// TODO: maybe have this threshold configurable
|
||||
if (bin_coverage >= 0.0001) {
|
||||
new_statistical_weights(bin_index) +=
|
||||
statistical_weights * bin_coverage_factor -
|
||||
1.0; //- 1 to avoid division by zero - initiallized with 1.
|
||||
statistical_weights * bin_coverage_factor;
|
||||
bin_counts(bin_index) += statistical_weights *
|
||||
bin_coverage_factor *
|
||||
photon_count_per_bin;
|
||||
@ -271,6 +331,9 @@ void AngleCalibration::redistribute_photon_counts_to_fixed_angle_bins(
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// std::string filename = "angle_widths.txt";
|
||||
// save(angle_widths, filename);
|
||||
}
|
||||
|
||||
void AngleCalibration::write_to_file(const std::string &filename) {
|
||||
@ -281,11 +344,17 @@ void AngleCalibration::write_to_file(const std::string &filename) {
|
||||
<< std::endl; // TODO: replace with log
|
||||
}
|
||||
|
||||
output_file << std::fixed << std::setprecision(6);
|
||||
output_file << std::fixed << std::setprecision(15);
|
||||
|
||||
for (ssize_t i = 0; i < num_bins; ++i) {
|
||||
if (new_photon_counts[i] == 0)
|
||||
if (new_photon_counts[i] <= std::numeric_limits<double>::epsilon()) {
|
||||
continue;
|
||||
output_file << mythen_detector->min_angle() + i * histogram_bin_width
|
||||
}
|
||||
|
||||
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;
|
||||
}
|
||||
|
@ -9,12 +9,58 @@
|
||||
|
||||
#include "test_config.hpp"
|
||||
|
||||
#include <iomanip>
|
||||
|
||||
#include <catch2/catch_all.hpp>
|
||||
#include <catch2/catch_test_macros.hpp>
|
||||
#include <catch2/matchers/catch_matchers_floating_point.hpp>
|
||||
|
||||
using namespace aare;
|
||||
|
||||
template <typename T, ssize_t Ndim = 1>
|
||||
NDArray<T, Ndim> read_into_array(const std::string &filename,
|
||||
const std::array<ssize_t, 1> size) {
|
||||
std::string word;
|
||||
NDArray<T, Ndim> 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 <typename T, ssize_t Ndim>
|
||||
bool check_equality_of_arrays(NDView<T, Ndim> array1, NDView<T, Ndim> array2) {
|
||||
bool equal = true;
|
||||
for (ssize_t i = 0; i < array1.size(); ++i) {
|
||||
if (std::abs(array1[i] - array2[i]) > 1e-10) {
|
||||
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:
|
||||
@ -144,6 +190,110 @@ 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<MythenDetectorSpecifications> mythen_detector_ptr =
|
||||
std::make_shared<MythenDetectorSpecifications>();
|
||||
|
||||
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<double, 1> expected_flatfield = read_into_array<double, 1>(
|
||||
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<double>::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<MythenDetectorSpecifications> mythen_detector_ptr =
|
||||
std::make_shared<MythenDetectorSpecifications>();
|
||||
|
||||
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<double, 1> expected_inverseflatfield = read_into_array<double, 1>(
|
||||
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]") {
|
||||
|
||||
@ -197,3 +347,93 @@ TEST_CASE("calculate new fixed angle width bins histogram",
|
||||
anglecalibration.write_to_file(
|
||||
"cpp_new_photon_counts.xye"); // TODO adjust output path
|
||||
}
|
||||
|
||||
TEST_CASE("check diffraction angles") {
|
||||
auto expected_diffraction_angle_filename = test_data_path() /
|
||||
"AngleCalibration_Test_Data" /
|
||||
"diffraction_angle.txt";
|
||||
|
||||
REQUIRE(std::filesystem::exists(expected_diffraction_angle_filename));
|
||||
|
||||
auto expected_angles =
|
||||
read_into_array<double, 1>(expected_diffraction_angle_filename.string(),
|
||||
std::array<ssize_t, 1>{61440});
|
||||
|
||||
auto diffraction_angle_filename =
|
||||
std::filesystem::current_path() / "../build/diffraction_angles.txt";
|
||||
|
||||
auto angles = load<double, 1>(diffraction_angle_filename,
|
||||
std::array<ssize_t, 1>{61440});
|
||||
|
||||
CHECK(check_equality_of_arrays(angles.view(), expected_angles.view()));
|
||||
}
|
||||
|
||||
TEST_CASE("check angle widths") {
|
||||
auto expected_filename =
|
||||
test_data_path() / "AngleCalibration_Test_Data" / "angle_width.txt";
|
||||
|
||||
REQUIRE(std::filesystem::exists(expected_filename));
|
||||
|
||||
auto expected_array = read_into_array<double, 1>(
|
||||
expected_filename.string(), std::array<ssize_t, 1>{61440});
|
||||
|
||||
auto filename =
|
||||
std::filesystem::current_path() / "../build/angle_widths.txt";
|
||||
|
||||
auto array = load<double, 1>(filename, std::array<ssize_t, 1>{61440});
|
||||
|
||||
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));
|
||||
|
||||
std::shared_ptr<MythenDetectorSpecifications> mythen_detector_ptr =
|
||||
std::make_shared<MythenDetectorSpecifications>();
|
||||
|
||||
std::shared_ptr<FlatField> flat_field_ptr =
|
||||
std::make_shared<FlatField>(mythen_detector_ptr);
|
||||
|
||||
std::shared_ptr<MythenFileReader> mythen_file_reader_ptr =
|
||||
std::make_shared<MythenFileReader>(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;
|
||||
|
||||
double diffraction_angle_DG_param =
|
||||
anglecalibration.diffraction_angle_from_DG_parameters(
|
||||
center, conversion, offset, 1);
|
||||
|
||||
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);
|
||||
|
||||
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);
|
||||
|
||||
double strip_width_EE_param =
|
||||
anglecalibration.angular_strip_width_from_EE_parameters(
|
||||
distance_center, normal_distance, angle, global_strip_index);
|
||||
|
||||
CHECK(strip_width_DG_param == Catch::Approx(strip_width_EE_param));
|
||||
}
|
||||
|
Reference in New Issue
Block a user