From 023e08e80df961c2014d03c8975fdc610f0b4252 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sun, 17 May 2026 16:08:48 +0200 Subject: [PATCH] Integrate consensus unit cell handling --- common/CMakeLists.txt | 2 + common/IndexingSettings.cpp | 4 + common/IndexingSettings.h | 2 + common/JFJochMessages.h | 1 + common/UnitCell.cpp | 71 ++++++++++++++++ common/UnitCell.h | 17 ++-- docs/CBOR.md | 1 + frame_serialize/CBORStream2Deserializer.cpp | 2 + frame_serialize/CBORStream2Serializer.cpp | 1 + image_analysis/IndexAndRefine.cpp | 74 ++++++++++++++++- image_analysis/IndexAndRefine.h | 4 + image_analysis/scale_merge/Merge.cpp | 28 +++++++ image_analysis/scale_merge/Merge.h | 11 ++- receiver/JFJochReceiver.cpp | 1 + tests/CMakeLists.txt | 1 + tests/UnitCellTest.cpp | 89 +++++++++++++++++++++ tools/jfjoch_process.cpp | 33 ++++++-- writer/HDF5NXmx.cpp | 12 +-- 18 files changed, 328 insertions(+), 26 deletions(-) create mode 100644 common/UnitCell.cpp create mode 100644 tests/UnitCellTest.cpp diff --git a/common/CMakeLists.txt b/common/CMakeLists.txt index a15a652b..0f5b5574 100644 --- a/common/CMakeLists.txt +++ b/common/CMakeLists.txt @@ -128,6 +128,8 @@ ADD_LIBRARY(JFJochCommon STATIC BrokerStatus.h ScalingSettings.cpp ScalingSettings.h + UnitCell.cpp + UnitCell.h ) TARGET_LINK_LIBRARIES(JFJochCommon JFJochLogger Compression JFCalibration gemmi Threads::Threads -lrt ) diff --git a/common/IndexingSettings.cpp b/common/IndexingSettings.cpp index af1bb46b..7cd61a83 100644 --- a/common/IndexingSettings.cpp +++ b/common/IndexingSettings.cpp @@ -152,6 +152,10 @@ float IndexingSettings::GetUnitCellDistTolerance() const { return unit_cell_dist_tolerance_vs_reference; } +float IndexingSettings::GetUnitCellAngleTolerance_deg() const { + return unit_cell_angle_tolerance_deg; +} + GeomRefinementAlgorithmEnum IndexingSettings::GetGeomRefinementAlgorithm() const { return refinement; } diff --git a/common/IndexingSettings.h b/common/IndexingSettings.h index 3bcee607..35eda787 100644 --- a/common/IndexingSettings.h +++ b/common/IndexingSettings.h @@ -20,6 +20,7 @@ class IndexingSettings { float indexing_tolerance = 0.1; float max_angle_from_ewald_deg = 2.0; float unit_cell_dist_tolerance_vs_reference = 0.05; // relative + static constexpr float unit_cell_angle_tolerance_deg = 5.0; // degree int64_t indexing_threads = 4; int64_t viable_cell_min_spots = 9; @@ -64,6 +65,7 @@ public: [[nodiscard]] float GetFFT_MaxAngle_deg() const; [[nodiscard]] int64_t GetIndexingThreads() const; [[nodiscard]] float GetUnitCellDistTolerance() const; + [[nodiscard]] float GetUnitCellAngleTolerance_deg() const; [[nodiscard]] bool GetIndexIceRings() const; [[nodiscard]] bool GetRotationIndexing() const; [[nodiscard]] float GetRotationIndexingMinAngularRange_deg() const; diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h index 5463829e..44a617d3 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -322,6 +322,7 @@ struct EndMessage { std::optional rotation_lattice_type; std::optional rotation_lattice; + std::optional unit_cell; // Vectors with end result: std::vector data_collection_efficiency; diff --git a/common/UnitCell.cpp b/common/UnitCell.cpp new file mode 100644 index 00000000..30da2878 --- /dev/null +++ b/common/UnitCell.cpp @@ -0,0 +1,71 @@ +// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "UnitCell.h" + +#include +#include + +bool UnitCell::is_finite() const { + return std::isfinite(a) + && std::isfinite(b) + && std::isfinite(c) + && std::isfinite(alpha) + && std::isfinite(beta) + && std::isfinite(gamma); +} + +bool UnitCell::is_close(const UnitCell &ref, float dist_tolerance, float angle_tolerance_deg) const { + if (!is_finite() || !ref.is_finite()) + return false; + + constexpr float min_length = 1e-6f; + + if (std::fabs(a - ref.a) / std::max(std::fabs(ref.a), min_length) > dist_tolerance) + return false; + if (std::fabs(b - ref.b) / std::max(std::fabs(ref.b), min_length) > dist_tolerance) + return false; + if (std::fabs(c - ref.c) / std::max(std::fabs(ref.c), min_length) > dist_tolerance) + return false; + + if (std::fabs(alpha - ref.alpha) > angle_tolerance_deg) + return false; + if (std::fabs(beta - ref.beta) > angle_tolerance_deg) + return false; + if (std::fabs(gamma - ref.gamma) > angle_tolerance_deg) + return false; + + return true; +} + +std::ostream &operator<<(std::ostream &output, const UnitCell &in) { + output << in.a << " " << in.b << " " << in.c; + output << " " << in.alpha << " " << in.beta << " " << in.gamma; + return output; +} + +std::optional MeanUnitCell(const std::vector &cells) { + if (cells.empty()) + return {}; + + UnitCell ret{}; + for (const auto &cell: cells) { + ret.a += cell.a; + ret.b += cell.b; + ret.c += cell.c; + ret.alpha += cell.alpha; + ret.beta += cell.beta; + ret.gamma += cell.gamma; + } + + const auto scale = 1.0f / static_cast(cells.size()); + ret.a *= scale; + ret.b *= scale; + ret.c *= scale; + ret.alpha *= scale; + ret.beta *= scale; + ret.gamma *= scale; + + return ret; +} + diff --git a/common/UnitCell.h b/common/UnitCell.h index 0216b491..3aa196dc 100644 --- a/common/UnitCell.h +++ b/common/UnitCell.h @@ -1,10 +1,11 @@ // SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only -#ifndef JUNGFRAUJOCH_UNITCELL_H -#define JUNGFRAUJOCH_UNITCELL_H +#pragma once #include +#include +#include struct UnitCell { float a; @@ -13,12 +14,10 @@ struct UnitCell { float alpha; float beta; float gamma; + + bool is_finite() const; + bool is_close(const UnitCell &ref, float dist_tolerance, float angle_tolerance_deg) const; }; -inline std::ostream &operator<<( std::ostream &output, const UnitCell &in ) { - output << in.a << " " << in.b << " " << in.c; - output << " " << in.alpha << " " << in.beta << " " << in.gamma; - return output; -} - -#endif //JUNGFRAUJOCH_UNITCELL_H +std::ostream &operator<<(std::ostream &output, const UnitCell &in); +std::optional MeanUnitCell(const std::vector &cells); diff --git a/docs/CBOR.md b/docs/CBOR.md index e9abbd27..b14d9d83 100644 --- a/docs/CBOR.md +++ b/docs/CBOR.md @@ -273,6 +273,7 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s | max_receiver_delay | uint64 | Internal performance of Jungfraujoch | | | bkg_estimate | float | Mean background estimate for the whole run | | | indexing_rate | float | Mean indexing rate for the whole run | | +| unit_cell | object (optional) | Unit cell of the system, based on the actual experiment: a, b, c \[angstrom\] and alpha, beta, gamma \[degree\] | | | rotation_lattice_type | object | Bravais lattice classification of the total rotation solution over the run, if available; same schema as `lattice_type` | | | - centering | string | One-letter centering code: P, A, B, C, I, F, or R | | | - niggli_class | int64 | Integer identifier for the Niggli-reduced Bravais class | | diff --git a/frame_serialize/CBORStream2Deserializer.cpp b/frame_serialize/CBORStream2Deserializer.cpp index fe8c1e87..db1c1026 100644 --- a/frame_serialize/CBORStream2Deserializer.cpp +++ b/frame_serialize/CBORStream2Deserializer.cpp @@ -1376,6 +1376,8 @@ namespace { GetCBORInt32Array(value, message.integrated_reflections); else if (key == "niggli_class") GetCBORUInt8Array(value, message.niggli_class); + else if (key == "unit_cell") + message.unit_cell = ProcessUnitCellElement(value); else if (key == "pixel_sum") GetCBORInt64Array(value, message.pixel_sum); else if (key == "rotation_lattice") { diff --git a/frame_serialize/CBORStream2Serializer.cpp b/frame_serialize/CBORStream2Serializer.cpp index 776969ad..8834a467 100644 --- a/frame_serialize/CBORStream2Serializer.cpp +++ b/frame_serialize/CBORStream2Serializer.cpp @@ -719,6 +719,7 @@ void CBORStream2Serializer::SerializeSequenceEnd(const EndMessage& message) { CBOR_ENC(mapEncoder, "integrated_reflections", message.integrated_reflections); CBOR_ENC(mapEncoder, "niggli_class", message.niggli_class); CBOR_ENC(mapEncoder, "pixel_sum", message.pixel_sum); + CBOR_ENC(mapEncoder, "unit_cell", message.unit_cell); cborErr(cbor_encoder_close_container(&encoder, &mapEncoder)); diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 94abb0eb..0f05f180 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -263,7 +263,10 @@ void IndexAndRefine::ProcessImage(DataMessage &msg, if (!AnalyzeIndexing(msg, outcome.experiment, *outcome.lattice_candidate)) return; - unit_cells[msg.number] = outcome.lattice_candidate->GetUnitCell(); + { + std::unique_lock ul(reflections_mutex); + unit_cells[msg.number] = outcome.lattice_candidate->GetUnitCell(); + } msg.lattice_type = outcome.symmetry; if (spot_finding_settings.quick_integration) @@ -312,3 +315,72 @@ const std::vector &IndexAndRefine::GetImageCC() const { return scale_cc; } +const std::vector > & IndexAndRefine::GetUnitCells() const { + return unit_cells; +} + +std::optional IndexAndRefine::GetConsensusUnitCell() const { + const auto dist_tolerance = experiment.GetIndexingSettings().GetUnitCellDistTolerance(); + const auto angle_tolerance = experiment.GetIndexingSettings().GetUnitCellAngleTolerance_deg(); + + if (rotation_indexer) { + auto result = rotation_indexer->GetLattice(); + if (!result) + return {}; + return result->lattice.GetUnitCell(); + } + + std::vector cells; + { + std::unique_lock ul(reflections_mutex); + cells.reserve(unit_cells.size()); + for (const auto &cell: unit_cells) { + if (cell && cell->is_finite()) + cells.emplace_back(*cell); + } + } + + if (cells.empty()) + return {}; + + if (experiment.GetUnitCell()) { + std::vector accepted; + accepted.reserve(cells.size()); + + for (const auto &cell: cells) { + if (cell.is_close(*experiment.GetUnitCell(), dist_tolerance, angle_tolerance)) + accepted.emplace_back(cell); + } + + return MeanUnitCell(accepted); + } + + size_t best_count = 0; + UnitCell best_reference{}; + + for (const auto &ref: cells) { + size_t count = 0; + for (const auto &cell: cells) { + if (cell.is_close(ref, dist_tolerance, angle_tolerance)) + ++count; + } + + if (count > best_count) { + best_count = count; + best_reference = ref; + } + } + + if (best_count == 0) + return {}; + + std::vector accepted; + accepted.reserve(best_count); + + for (const auto &cell: cells) { + if (cell.is_close(best_reference, dist_tolerance, angle_tolerance)) + accepted.emplace_back(cell); + } + + return MeanUnitCell(accepted); +} \ No newline at end of file diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index 0e5553e1..216f25cb 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -70,6 +70,10 @@ public: ScalingResult ScaleAllImages(const std::vector &reference, size_t nthreads = 0); std::optional Finalize(); + + // Not thread safe, need to be run after processing is all done const std::vector> &GetReflections() const; const std::vector &GetImageCC() const; + std::optional GetConsensusUnitCell() const; + const std::vector > &GetUnitCells() const; }; diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 9cf01704..650a0e99 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -36,6 +36,33 @@ size_t CalcMergeMask(const DiffractionExperiment &x, return n_rejected; } +size_t CalcMergeMaskUnitCell(const DiffractionExperiment &x, + const UnitCell &ref_cell, + const std::vector > &unit_cells, + std::vector &mask) { + if (mask.size() != unit_cells.size()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Mismatch in vector size"); + + size_t n_rejected = 0; + + const auto dtol = x.GetIndexingSettings().GetUnitCellDistTolerance(); + const auto atol = x.GetIndexingSettings().GetUnitCellAngleTolerance_deg(); + for (size_t i = 0; i < unit_cells.size(); ++i) { + if (mask[i]) { + // Counter should not include trivial case (non-indexed image or weird things) + // It should only count where indexing did happen, but cell was too far from reference + if (!unit_cells[i] || !unit_cells[i]->is_finite()) + mask[i] = 0; + else if (!unit_cells[i]->is_close(ref_cell, dtol, atol)) { + mask[i] = 0; + ++n_rejected; + } + } + } + + return n_rejected; +} + std::vector MergeAll(const DiffractionExperiment &x, const std::vector > &reflections, const std::vector &merge_mask) { @@ -346,3 +373,4 @@ void MergeStatistics::Print(Logger &logger) const { } logger.Info(""); } + diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index b235aa46..f9bfe0b3 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -31,9 +31,14 @@ struct MergeResult { MergeStatistics statistics; }; -size_t CalcMergeMask(const DiffractionExperiment &x, - const std::vector &scale_cc, - std::vector &result_mask); +size_t CalcMergeMaskCC(const DiffractionExperiment &x, + const std::vector &scale_cc, + std::vector &result_mask); + +size_t CalcMergeMaskUnitCell(const DiffractionExperiment &x, + const UnitCell &reference_cell, + const std::vector> &unit_cells, + std::vector &mask); std::vector MergeAll(const DiffractionExperiment &x, const std::vector > &reflections, diff --git a/receiver/JFJochReceiver.cpp b/receiver/JFJochReceiver.cpp index 46bbd53b..8997ef9e 100644 --- a/receiver/JFJochReceiver.cpp +++ b/receiver/JFJochReceiver.cpp @@ -171,6 +171,7 @@ void JFJochReceiver::SendEndMessage() { }; rotation_indexing_lattice = rotation_indexer_ret->lattice; } + message.unit_cell = indexer.GetConsensusUnitCell(); for (int i = 0; i < adu_histogram_module.size(); i++) message.adu_histogram["module" + std::to_string(i)] = adu_histogram_module[i]->GetHistogram(); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 1761d5ac..b38c8e97 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -68,6 +68,7 @@ ADD_EXECUTABLE(jfjoch_test SearchSpaceGroupTest.cpp XDSPluginTest.cpp MergeScaleTest.cpp + UnitCellTest.cpp ) target_link_libraries(jfjoch_test Catch2WithMain JFJochBroker JFJochReceiver JFJochReader JFJochWriter diff --git a/tests/UnitCellTest.cpp b/tests/UnitCellTest.cpp new file mode 100644 index 00000000..b89cd3d0 --- /dev/null +++ b/tests/UnitCellTest.cpp @@ -0,0 +1,89 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include "../common/UnitCell.h" + +#include +#include + +TEST_CASE("UnitCell_is_finite") { + UnitCell cell{50, 60, 70, 90, 90, 120}; + CHECK(cell.is_finite()); + + cell.a = std::numeric_limits::quiet_NaN(); + CHECK_FALSE(cell.is_finite()); + + cell = UnitCell{50, 60, 70, 90, 90, 120}; + cell.beta = std::numeric_limits::infinity(); + CHECK_FALSE(cell.is_finite()); + + cell = UnitCell{50, 60, 70, 90, -std::numeric_limits::infinity(), 120}; + CHECK_FALSE(cell.is_finite()); +} + +TEST_CASE("UnitCell_is_close_accepts_cells_within_tolerance") { + const UnitCell ref{100, 120, 150, 90, 95, 120}; + + CHECK(UnitCell{100, 120, 150, 90, 95, 120}.is_close(ref, 0.05f, 5.0f)); + + CHECK(UnitCell{105, 126, 157.5, 95, 100, 125}.is_close(ref, 0.05f, 5.0f)); + CHECK(UnitCell{95, 114, 142.5, 85, 90, 115}.is_close(ref, 0.05f, 5.0f)); +} + +TEST_CASE("UnitCell_is_close_rejects_length_outliers") { + const UnitCell ref{100, 120, 150, 90, 95, 120}; + + CHECK_FALSE(UnitCell{105.1f, 120, 150, 90, 95, 120}.is_close(ref, 0.05f, 5.0f)); + CHECK_FALSE(UnitCell{100, 126.1f, 150, 90, 95, 120}.is_close(ref, 0.05f, 5.0f)); + CHECK_FALSE(UnitCell{100, 120, 157.6f, 90, 95, 120}.is_close(ref, 0.05f, 5.0f)); + + CHECK_FALSE(UnitCell{94.9f, 120, 150, 90, 95, 120}.is_close(ref, 0.05f, 5.0f)); + CHECK_FALSE(UnitCell{100, 113.9f, 150, 90, 95, 120}.is_close(ref, 0.05f, 5.0f)); + CHECK_FALSE(UnitCell{100, 120, 142.4f, 90, 95, 120}.is_close(ref, 0.05f, 5.0f)); +} + +TEST_CASE("UnitCell_is_close_rejects_angle_outliers") { + const UnitCell ref{100, 120, 150, 90, 95, 120}; + + CHECK_FALSE(UnitCell{100, 120, 150, 95.1f, 95, 120}.is_close(ref, 0.05f, 5.0f)); + CHECK_FALSE(UnitCell{100, 120, 150, 90, 100.1f, 120}.is_close(ref, 0.05f, 5.0f)); + CHECK_FALSE(UnitCell{100, 120, 150, 90, 95, 125.1f}.is_close(ref, 0.05f, 5.0f)); + + CHECK_FALSE(UnitCell{100, 120, 150, 84.9f, 95, 120}.is_close(ref, 0.05f, 5.0f)); + CHECK_FALSE(UnitCell{100, 120, 150, 90, 89.9f, 120}.is_close(ref, 0.05f, 5.0f)); + CHECK_FALSE(UnitCell{100, 120, 150, 90, 95, 114.9f}.is_close(ref, 0.05f, 5.0f)); +} + +TEST_CASE("UnitCell_is_close_rejects_non_finite_cells") { + const UnitCell ref{100, 120, 150, 90, 95, 120}; + + UnitCell cell{100, 120, 150, 90, 95, 120}; + cell.a = std::numeric_limits::quiet_NaN(); + CHECK_FALSE(cell.is_close(ref, 0.05f, 5.0f)); + + UnitCell bad_ref{100, 120, 150, 90, 95, 120}; + bad_ref.gamma = std::numeric_limits::infinity(); + CHECK_FALSE(ref.is_close(bad_ref, 0.05f, 5.0f)); +} + +TEST_CASE("UnitCell_MeanUnitCell") { + const std::vector cells{ + {100, 120, 150, 90, 95, 120}, + {102, 122, 152, 92, 97, 122}, + {98, 118, 148, 88, 93, 118} + }; + + const auto mean = MeanUnitCell(cells); + REQUIRE(mean.has_value()); + + CHECK(mean->a == Catch::Approx(100)); + CHECK(mean->b == Catch::Approx(120)); + CHECK(mean->c == Catch::Approx(150)); + CHECK(mean->alpha == Catch::Approx(90)); + CHECK(mean->beta == Catch::Approx(95)); + CHECK(mean->gamma == Catch::Approx(120)); + + CHECK_FALSE(MeanUnitCell({}).has_value()); +} \ No newline at end of file diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index e16e6ef2..4551aecf 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -600,7 +600,25 @@ int main(int argc, char **argv) { logger.Info("Rotation Indexing found lattice"); } - // --- Optional: run scaling (mosaicity refinement) on accumulated reflections --- + std::vector merging_mask_uc(images_to_process, 1); + + auto consensus_start_time = std::chrono::steady_clock::now(); + const auto consensus_cell = indexer.GetConsensusUnitCell(); + auto consensus_end_time = std::chrono::steady_clock::now(); + auto consensus_duration = std::chrono::duration(consensus_end_time - consensus_start_time).count(); + if (consensus_cell) { + logger.Info("Consensus unit cell found in {:.2f} ms", consensus_duration * 1e3); + logger.Info("UC: a={:.2f} b={:.2f} c={:.2f} alpha={:.2f} beta={:.2f} gamma={:.2f}", + consensus_cell->a, consensus_cell->b, consensus_cell->c, + consensus_cell->alpha, consensus_cell->beta, consensus_cell->gamma); + + auto rejected_uc = CalcMergeMaskUnitCell(experiment, *consensus_cell, indexer.GetUnitCells(), merging_mask_uc); + if (rejected_uc > 0) + logger.Info("Rejected {} images for merging due to unit cell being too far from consensus", rejected_uc); + } else + logger.Info("Consensus unit cell not found - calculation tool {:.2f} ms", consensus_duration * 1e3); + end_msg.unit_cell = consensus_cell; + if (run_scaling || !reference_data.empty()) { logger.Info("Running scaling (mosaicity refinement) ..."); @@ -613,7 +631,7 @@ int main(int argc, char **argv) { auto scale_start = std::chrono::steady_clock::now(); for (int i = 0; i < scaling_iter; i++) { auto iter_start = std::chrono::steady_clock::now(); - auto merge_result = MergeAll(experiment, indexer.GetReflections()); + auto merge_result = MergeAll(experiment, indexer.GetReflections(), merging_mask_uc); scale_result = indexer.ScaleAllImages(merge_result); scale_result.SaveToFile(output_prefix + "_iter" + std::to_string(i) + "_scale.dat"); auto iter_end = std::chrono::steady_clock::now(); @@ -635,13 +653,12 @@ int main(int argc, char **argv) { auto merge_start = std::chrono::steady_clock::now(); - std::vector merge_mask(images_to_process, 1); - auto rejected = CalcMergeMask(experiment, scale_cc, merge_mask); - if (rejected > 0) - logger.Info("Rejected {} images due to low CC with reference", rejected); + auto rejected_cc = CalcMergeMaskCC(experiment, scale_cc, merging_mask_uc); + if (rejected_cc > 0) + logger.Info("Rejected {} images for merging due to low CC with reference", rejected_cc); - auto merged_reflections = MergeAll(experiment, indexer.GetReflections(), merge_mask); - auto merged_statistics = MergeStats(experiment, merged_reflections, indexer.GetReflections(), merge_mask); + auto merged_reflections = MergeAll(experiment, indexer.GetReflections(), merging_mask_uc); + auto merged_statistics = MergeStats(experiment, merged_reflections, indexer.GetReflections(), merging_mask_uc); auto merge_end = std::chrono::steady_clock::now(); double merge_time = std::chrono::duration(merge_end - merge_start).count(); diff --git a/writer/HDF5NXmx.cpp b/writer/HDF5NXmx.cpp index 8728dd16..13c03e39 100644 --- a/writer/HDF5NXmx.cpp +++ b/writer/HDF5NXmx.cpp @@ -562,19 +562,21 @@ void NXmx::Sample(const StartMessage &start, const EndMessage &end) { } std::optional unit_cell; - std::optional set_unit_cell; - if (end.rotation_lattice) + std::optional input_unit_cell; + if (end.unit_cell) + unit_cell = end.unit_cell; + else if (end.rotation_lattice) unit_cell = end.rotation_lattice->GetUnitCell(); else if (start.unit_cell) { unit_cell = start.unit_cell; - set_unit_cell = start.unit_cell; + input_unit_cell = start.unit_cell; } if (unit_cell) SaveUnitCell(group, "unit_cell", unit_cell.value()); - if (set_unit_cell) - SaveUnitCell(group, "input_unit_cell", unit_cell.value()); + if (input_unit_cell) + SaveUnitCell(group, "input_unit_cell", input_unit_cell.value()); if (end.rotation_lattice) { group.SaveVector("ub_matrix",