diff --git a/broker/OpenAPIConvert.cpp b/broker/OpenAPIConvert.cpp index 2845c9d0a..8ef834086 100644 --- a/broker/OpenAPIConvert.cpp +++ b/broker/OpenAPIConvert.cpp @@ -2,6 +2,7 @@ // SPDX-License-Identifier: GPL-3.0-only #include "../common/JFJochMath.h" +#include "../common/GitInfo.h" #include "OpenAPIConvert.h" // From https://en.cppreference.com/w/cpp/string/byte/tolower @@ -1183,6 +1184,89 @@ org::openapitools::server::model::Analysis_settings Convert(const AnalysisSettin return ret; } +org::openapitools::server::model::Powder_calibration_output Convert(const CalibrationResult &input) { + org::openapitools::server::model::Powder_calibration_output ret; + const auto &g = input.geometry; + + // The geometry under this API's own property names, so it can be PUT back without translating a + // field. Same content as the .json rugnux writes (WriteCalibrationJson). + org::openapitools::server::model::Dataset_settings settings; + settings.setBeamXPxl(g.GetBeamX_pxl()); + settings.setBeamYPxl(g.GetBeamY_pxl()); + settings.setDetectorDistanceMm(g.GetDetectorDistance_mm()); + settings.setIncidentEnergyKeV(WVL_1A_IN_KEV / g.GetWavelength_A()); + // The rotations go in together or not at all: omitting them would state a FLAT detector, since the + // schema's default is 0.0, rather than leave the tilt unstated. + if (g.GetPoniRot1_rad() != 0.0f || g.GetPoniRot2_rad() != 0.0f || g.GetPoniRot3_rad() != 0.0f) { + settings.setPoniRot1Rad(g.GetPoniRot1_rad()); + settings.setPoniRot2Rad(g.GetPoniRot2_rad()); + settings.setPoniRot3Rad(g.GetPoniRot3_rad()); + } + ret.setDatasetSettings(settings); + + org::openapitools::server::model::Powder_calibration_quality quality; + // The verdict first: a caller that reads nothing else must still not mistake a non-fit for a fit. + quality.setConverged(input.converged); + if (!input.converged) + quality.setNotConvergedReason(input.reason); + quality.setRingPoints(input.ring_points); + quality.setRmsRadialPxl(input.rms_radial_pxl); + quality.setBeamSigmaPxl(input.beam_sigma_pxl); + // Where the beam actually lands, which is not the PONI above once the detector is tilted. + const auto [direct_x, direct_y] = g.GetDirectBeam_pxl(); + quality.setDirectBeamXPxl(direct_x); + quality.setDirectBeamYPxl(direct_y); + quality.setHeaderDistanceMm(input.header_distance_mm); + quality.setTiltRefined(input.tilt_refined); + quality.setTiltSignificance(input.tilt_significance); + if (input.seed_distance_mm > 0.0f) + quality.setRingSeedDistanceMm(input.seed_distance_mm); + + if (input.uncertainty.valid) { + org::openapitools::server::model::Powder_calibration_fit_sigma sigma; + sigma.setBeamXPxl(input.uncertainty.sigma_beam_x_pxl); + sigma.setBeamYPxl(input.uncertainty.sigma_beam_y_pxl); + sigma.setDetectorDistanceMm(input.uncertainty.sigma_distance_mm); + if (input.uncertainty.sigma_rot1_rad > 0.0 || input.uncertainty.sigma_rot2_rad > 0.0) { + sigma.setPoniRot1Rad(input.uncertainty.sigma_rot1_rad); + sigma.setPoniRot2Rad(input.uncertainty.sigma_rot2_rad); + sigma.setCorrelationBeamXRot1(input.uncertainty.corr_beam_x_rot1); + sigma.setCorrelationBeamYRot2(input.uncertainty.corr_beam_y_rot2); + } + quality.setFitSigma(sigma); + } + + if (input.spots_geometry) { + org::openapitools::server::model::Powder_calibration_spot_check cross; + cross.setBeamXPxl(input.spots_geometry->GetBeamX_pxl()); + cross.setBeamYPxl(input.spots_geometry->GetBeamY_pxl()); + cross.setDetectorDistanceMm(input.spots_geometry->GetDetectorDistance_mm()); + cross.setDisagreementPxl(input.spots_disagreement_pxl); + quality.setSpotCrossCheck(cross); + } + + ret.setCalibration(quality); + ret.setJfjochVersion(jfjoch_version()); + return ret; +} + +org::openapitools::server::model::Grid_scan_crystal Convert(const GridScanCrystal &input) { + org::openapitools::server::model::Grid_scan_crystal ret; + ret.setNx(input.nx); + ret.setNy(input.ny); + ret.setXUm(input.x_um); + ret.setYUm(input.y_um); + ret.setImageNumber(input.image_number); + ret.setMajorUm(input.major_um); + ret.setMinorUm(input.minor_um); + ret.setAngleDeg(input.angle_deg); + ret.setScore(input.score); + ret.setIceScore(input.ice_score); + ret.setResA(input.res_A); + ret.setNImages(input.n_images); + return ret; +} + org::openapitools::server::model::Scan_result Convert(const ScanResult& input) { org::openapitools::server::model::Scan_result ret; ret.setFilePrefix(input.file_prefix); @@ -1260,6 +1344,21 @@ org::openapitools::server::model::Scan_result Convert(const ScanResult& input) { if (input.rotation_crystal_system && input.rotation_centering) ret.setRotationBravais(BravaisSymbol(*input.rotation_crystal_system, *input.rotation_centering)); + // The two mode-specific payloads. Optional and independent, exactly as the rotation_* members + // above are: present when that mode ran, absent otherwise. + if (input.grid) { + org::openapitools::server::model::Grid_scan_result grid; + std::vector crystals; + for (const auto &c : input.grid->crystals) + crystals.emplace_back(Convert(c)); + grid.setCrystals(crystals); + grid.setBeamSizeXUm(input.grid->beam_size_x_um); + grid.setBeamSizeYUm(input.grid->beam_size_y_um); + ret.setGrid(grid); + } + if (input.calibration) + ret.setCalibration(Convert(input.calibration.value())); + return ret; } diff --git a/broker/OpenAPIConvert.h b/broker/OpenAPIConvert.h index 9c0cbbf8b..7d0df3d57 100644 --- a/broker/OpenAPIConvert.h +++ b/broker/OpenAPIConvert.h @@ -30,6 +30,8 @@ #include "gen/model/Indexing_settings.h" #include "gen/model/Bragg_integration_settings.h" #include "gen/model/Analysis_settings.h" +#include "gen/model/Grid_scan_result.h" +#include "gen/model/Powder_calibration_output.h" #include "gen/model/Scan_result.h" #include "../common/JFJochMessages.h" @@ -56,6 +58,9 @@ org::openapitools::server::model::Bragg_integration_settings Convert(const Bragg AnalysisSettings Convert(const org::openapitools::server::model::Analysis_settings &input); org::openapitools::server::model::Analysis_settings Convert(const AnalysisSettings &input); +org::openapitools::server::model::Grid_scan_crystal Convert(const GridScanCrystal &input); +org::openapitools::server::model::Powder_calibration_output Convert(const CalibrationResult &input); + org::openapitools::server::model::Measurement_statistics Convert(const MeasurementStatistics &input); DetectorSettings Convert(const org::openapitools::server::model::Detector_settings &input); diff --git a/broker/jfjoch_api.yaml b/broker/jfjoch_api.yaml index f8e883b4c..52b5a5b6b 100644 --- a/broker/jfjoch_api.yaml +++ b/broker/jfjoch_api.yaml @@ -1624,6 +1624,10 @@ components: description: | Bravais lattice type of the global rotation-indexing solution: crystal-system letter + centering, e.g. "tP", "oC", "cF", "hR", "mP". Empty if not determined. + grid: + $ref: '#/components/schemas/grid_scan_result' + calibration: + $ref: '#/components/schemas/powder_calibration_output' images: type: array items: diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h index fb98b1ec5..2eec2f8d6 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -23,6 +23,7 @@ #include "CrystalLattice.h" #include "IndexingSettings.h" #include "AnalysisSettings.h" +#include "GridScanResult.h" #include "XrayFluorescenceSpectrum.h" #include "../gemmi_gph/gemmi/symmetry.hpp" @@ -404,6 +405,9 @@ struct EndMessage { std::optional rotation_lattice_type; std::optional rotation_lattice; std::vector rotation_extra_lattices; + // The crystals a grid scan found, best score first. Empty when the run was not a raster, or when + // the raster found nothing; a raster can find any number, so nothing here caps it at one. + std::vector grid_crystals; std::optional unit_cell; // Space group determined by the offline analysis (overrides the start message when writing the // master, since it is only known after merging). Both spellings are carried and a reader prefers diff --git a/common/ScanResult.h b/common/ScanResult.h index a7a96b49a..4bc550856 100644 --- a/common/ScanResult.h +++ b/common/ScanResult.h @@ -11,6 +11,8 @@ #include "CrystalLattice.h" #include "UnitCell.h" +#include "GridScanResult.h" +#include "../image_analysis/geom_refinement/PowderCalibration.h" // CalibrationResult struct ScanResultElem { int64_t number = -1; @@ -55,6 +57,13 @@ struct ScanResult { // Bravais lattice type of the global rotation-indexing solution (from RotationIndexer). std::optional rotation_crystal_system; std::optional rotation_centering; + + // What the run produced beyond the per-image list, one member per analysis mode that has an + // answer of its own. Same idiom as the three rotation_* members above: the payload of a mode is a + // top-level optional, present when that mode ran and absent otherwise, so images stays required + // and every existing /result/scan client is untouched. + std::optional grid; // AnalysisMode::Grid + std::optional calibration; // AnalysisMode::PowderCalibration }; // Two-letter Bravais lattice symbol: crystal-system letter + centering, e.g. "tP", "oC", "cF", "hR". diff --git a/docs/CBOR.md b/docs/CBOR.md index 8e26d963f..1b98d6f63 100644 --- a/docs/CBOR.md +++ b/docs/CBOR.md @@ -315,6 +315,16 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s | - system | string | Crystal system: triclinic, monoclinic, orthorhombic, tetragonal, trigonal, hexagonal, cubic | | | rotation_lattice | Array(9 * float) | Real-space lattice basis, flattened 3x3 in row-major order | | | rotation_extra_lattices | Array(Array(9*float)) | Additional indexed lattices (orientation variants); present only if found | | +| grid_crystals | Array(Map) | Crystals found by a grid scan, best score first; present only for a raster that found any. Any number of entries - a consumer must not assume at most one | | +| - nx, ny | float | Crystal centre as a fractional grid coordinate, along the fast and the slow grid axis | | +| - x_um, y_um | float | The same centre as a signed offset along those two axes \[um\] | | +| - image_number | int64 | Image ordinal of the grid point at the centre | | +| - major_um, minor_um | float | Extent along the crystal's major and minor principal axes \[um\] | | +| - angle_deg | float | Major axis from the +x grid axis, counter-clockwise \[deg\]. An AXIS, not a direction: in \[0, 180) and wrapping, so 179 and 0 are adjacent and these do not average arithmetically | | +| - score | float | Diffraction score the crystal was ranked on | | +| - ice_score | float | Ice-ring score over the crystal's grid points | | +| - res_A | float | Resolution estimate over the crystal's grid points \[A\] | | +| - n_images | int64 | Grid points that fell in this crystal | | | data_collection_efficiency_image | Array(float) | Per-image data collection efficiency. Missing values are encoded as 0 or 1 depending on producer context | | | spot_count | Array(int32) | Per-image spot count | | | spot_count_ice_ring | Array(int32) | Per-image number of spots within identified ice-ring resolution ranges | | diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 5011548b5..0246d8e42 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -8,6 +8,7 @@ * Rotation MX analysis is not offered online - it is absent from the REST API and refused by the broker - and a rotation sweep collected under an MX mode is logged as being analysed frame by frame as stills. * `rugnux --mode` accepts `mx_rotation` and `mx_stills` beside `mx`, which continues to choose between them from the goniometer. * The stream and the written file record which analysis produced them, as `analysis_mode` in the CBOR start message and `/entry/MX/analysis_mode` in the HDF5 master. +* The crystals a grid scan found travel through the CBOR end message (`grid_crystals`), the HDF5 master (`/entry/MX/crystals`) and `GET /result/scan`, which also carries a powder-calibration result. * `rugnux --model` reports CC(model, data) - the correlation of the merged intensities with the placed, scaled model - by resolution shell, on the same shells as CC1/2, with the reflection count and a significance for each. * `rugnux --model` fits the model's scale, anisotropic B and bulk-solvent parameters on the working reflections only, so the R-free it reports is measured against a model no free reflection helped scale. * The bulk-solvent parameters of `rugnux --model` are searched over their physically meaningful range instead of being fitted without bounds, so a model is never scaled with a solvent term that has silently switched itself off. diff --git a/docs/HDF5.md b/docs/HDF5.md index f45bf27d7..42a3b0552 100644 --- a/docs/HDF5.md +++ b/docs/HDF5.md @@ -375,6 +375,23 @@ variants. | `indexedLatticeCount` | | per-image lattice count summary (master). *Note: data files use `indexingLatticeCount`; readers accept either.* | | `reindexMatrix` | | change of basis from the setting the per-image data are in to the setting of `/entry/sample/unit_cell` (`[9]`, `int32`, flattened 3×3, row major) — see below | +**Grid-scan crystals** (master, `/entry/MX/crystals`, written when a raster found any). Parallel 1-D +datasets, all of length *N*, one per field, sorted by `score` descending — crystals are numbered +rather than named, so a group per crystal would put the index in the path. A raster can find any +number of them: nothing may assume at most one. + +| Dataset | Units | Meaning | +|---------|-------|---------| +| `nx`, `ny` | | crystal centre as a fractional grid coordinate, along the fast and the slow grid axis | +| `xUm`, `yUm` | µm | the same centre as a signed offset along those two axes | +| `imageNumber` | | image ordinal of the grid point at the centre | +| `majorUm`, `minorUm` | µm | extent along the crystal's major and minor principal axes | +| `angleDeg` | ° | major axis from the +x grid axis, counter-clockwise. An **axis, not a direction**: it lies in [0, 180) and wraps, so 179° and 0° are adjacent. Two of these do not average arithmetically — combine them on the doubled angle (mean of 2θ, halved), or a pair of nearly parallel needles comes out as a right angle | +| `score` | | diffraction score the crystal was ranked on | +| `iceScore` | ratio | ice-ring score over the crystal's grid points | +| `resA` | Å | resolution estimate over the crystal's grid points | +| `nImages` | | grid points that fell in this crystal | + **Reindex matrix.** The per-image `h`, `k`, `l` and `latticeIndexed` are written as each image is processed, in the setting that image was *indexed* in. The space group is only settled afterwards, by the merge, and settling it can re-seat the lattice into the group's conventional setting — so diff --git a/frame_serialize/CBORStream2Deserializer.cpp b/frame_serialize/CBORStream2Deserializer.cpp index ac27948f3..cd92d1aa1 100644 --- a/frame_serialize/CBORStream2Deserializer.cpp +++ b/frame_serialize/CBORStream2Deserializer.cpp @@ -595,6 +595,33 @@ namespace { return lm; } + GridScanCrystal GetCBORGridScanCrystal(CborValue &value) { + GridScanCrystal c{}; + CborValue map_value; + cborErr(cbor_value_enter_container(&value, &map_value)); + while (!cbor_value_at_end(&map_value)) { + auto key = GetCBORString(map_value); + if (key == "nx") c.nx = GetCBORFloat(map_value); + else if (key == "ny") c.ny = GetCBORFloat(map_value); + else if (key == "x_um") c.x_um = GetCBORFloat(map_value); + else if (key == "y_um") c.y_um = GetCBORFloat(map_value); + else if (key == "image_number") c.image_number = GetCBORInt(map_value); + else if (key == "major_um") c.major_um = GetCBORFloat(map_value); + else if (key == "minor_um") c.minor_um = GetCBORFloat(map_value); + else if (key == "angle_deg") c.angle_deg = GetCBORFloat(map_value); + else if (key == "score") c.score = GetCBORFloat(map_value); + else if (key == "ice_score") c.ice_score = GetCBORFloat(map_value); + else if (key == "res_A") c.res_A = GetCBORFloat(map_value); + else if (key == "n_images") c.n_images = GetCBORInt(map_value); + else + // A field this build does not know. Step over its value, or the loop would read the + // payload as the next key and never reach the end of the map. + cbor_value_advance(&map_value); + } + cborErr(cbor_value_leave_container(&value, &map_value)); + return c; + } + XrayFluorescenceSpectrum GetCBORFluorescenceSpectrum(CborValue &value) { std::vector data, energy; CborValue map_value; @@ -1553,7 +1580,14 @@ namespace { message.rotation_lattice = CrystalLattice(tmp); } else if (key == "rotation_lattice_type") message.rotation_lattice_type = GetCBORLatticeMessage(value); - else if (key == "rotation_extra_lattices") { + else if (key == "grid_crystals") { + size_t array_len = GetCBORArrayLen(value); + CborValue array_value; + cborErr(cbor_value_enter_container(&value, &array_value)); + for (size_t i = 0; i < array_len; i++) + message.grid_crystals.push_back(GetCBORGridScanCrystal(array_value)); + cborErr(cbor_value_leave_container(&value, &array_value)); + } else if (key == "rotation_extra_lattices") { size_t array_len = GetCBORArrayLen(value); CborValue array_value; cborErr(cbor_value_enter_container(&value, &array_value)); diff --git a/frame_serialize/CBORStream2Serializer.cpp b/frame_serialize/CBORStream2Serializer.cpp index 98bc4c003..55235b14c 100644 --- a/frame_serialize/CBORStream2Serializer.cpp +++ b/frame_serialize/CBORStream2Serializer.cpp @@ -302,6 +302,39 @@ inline void CBOR_ENC(CborEncoder &encoder, const char* key, const std::vector& crystals) { + CborEncoder arrayEncoder; + + cborErr(cbor_encode_text_stringz(&encoder, key)); + cborErr(cbor_encoder_create_array(&encoder, &arrayEncoder, crystals.size())); + + for (const auto &c: crystals) + CBOR_ENC(arrayEncoder, c); + + cborErr(cbor_encoder_close_container(&encoder, &arrayEncoder)); +} + inline void CBOR_ENC(CborEncoder &encoder, const char* key, const std::vector& refs) { CborEncoder arrayEncoder, mapEncoder; @@ -822,6 +855,8 @@ void CBORStream2Serializer::SerializeSequenceEnd(const EndMessage& message) { CBOR_ENC_FLOAT_ARRAY_NOKEY(arrayEncoder, el.GetVector()); cborErr(cbor_encoder_close_container(&mapEncoder, &arrayEncoder)); } + if (!message.grid_crystals.empty()) + CBOR_ENC(mapEncoder, "grid_crystals", message.grid_crystals); CBOR_ENC(mapEncoder, "data_collection_efficiency_image", message.data_collection_efficiency); CBOR_ENC(mapEncoder, "spot_count", message.spot_count); CBOR_ENC(mapEncoder, "spot_count_ice_ring", message.spot_count_ice_ring); diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index e56c7275c..f3d4a90b8 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -784,6 +784,47 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen dataset->file_analysis_mode = AnalysisModeFromName(master_file->GetString("/entry/MX/analysis_mode", "")); + // Grid-scan crystals: parallel 1-D datasets, one per field, all of length N. nx sets N and + // the rest are read into that length, so a file written by a build that had one field + // fewer still re-opens with the crystals it does have. + if (master_file->Exists("/entry/MX/crystals/nx")) { + auto nx = master_file->ReadOptVector("/entry/MX/crystals/nx"); + auto ny = master_file->ReadOptVector("/entry/MX/crystals/ny"); + auto x_um = master_file->ReadOptVector("/entry/MX/crystals/xUm"); + auto y_um = master_file->ReadOptVector("/entry/MX/crystals/yUm"); + auto image_number = master_file->ReadOptVector("/entry/MX/crystals/imageNumber"); + auto major_um = master_file->ReadOptVector("/entry/MX/crystals/majorUm"); + auto minor_um = master_file->ReadOptVector("/entry/MX/crystals/minorUm"); + auto angle_deg = master_file->ReadOptVector("/entry/MX/crystals/angleDeg"); + auto score = master_file->ReadOptVector("/entry/MX/crystals/score"); + auto ice_score = master_file->ReadOptVector("/entry/MX/crystals/iceScore"); + auto res_A = master_file->ReadOptVector("/entry/MX/crystals/resA"); + auto n_images = master_file->ReadOptVector("/entry/MX/crystals/nImages"); + + const size_t n = nx.size(); + auto f = [n](const std::vector &v, size_t i) { return v.size() == n ? v[i] : 0.0f; }; + auto i64 = [n](const std::vector &v, size_t i) { + return v.size() == n ? v[i] : int64_t(0); + }; + + dataset->grid_crystals.resize(n); + for (size_t i = 0; i < n; i++) { + auto &c = dataset->grid_crystals[i]; + c.nx = nx[i]; + c.ny = f(ny, i); + c.x_um = f(x_um, i); + c.y_um = f(y_um, i); + c.image_number = i64(image_number, i); + c.major_um = f(major_um, i); + c.minor_um = f(minor_um, i); + c.angle_deg = f(angle_deg, i); + c.score = f(score, i); + c.ice_score = f(ice_score, i); + c.res_A = f(res_A, i); + c.n_images = i64(n_images, i); + } + } + auto indexing = master_file->GetString("/entry/MX/indexing_algorithm", "none"); if (indexing == "fft" || indexing == "FFT (CUDA)" || indexing == "FFT (FFTW)") dataset->experiment.IndexingAlgorithm(IndexingAlgorithmEnum::FFT); diff --git a/reader/JFJochReaderDataset.h b/reader/JFJochReaderDataset.h index f458cc08d..1db1532e5 100644 --- a/reader/JFJochReaderDataset.h +++ b/reader/JFJochReaderDataset.h @@ -36,6 +36,10 @@ struct JFJochReaderDataset { // written before the mode existed. std::optional file_analysis_mode; + // The crystals a grid scan found (/entry/MX/crystals), best score first. Empty for a file that is + // not a raster, or whose raster found nothing. + std::vector grid_crystals; + std::string jfjoch_release; // Change of basis (3x3 integers, row major) from the setting the per-image reflections and diff --git a/tests/AnalysisSettingsTest.cpp b/tests/AnalysisSettingsTest.cpp index 1996c1eba..0142b1095 100644 --- a/tests/AnalysisSettingsTest.cpp +++ b/tests/AnalysisSettingsTest.cpp @@ -135,9 +135,42 @@ TEST_CASE("AnalysisMode_CBORStartAbsentMode", "[AnalysisMode][CBOR]") { CHECK_FALSE(deserialized->start_message->analysis_mode.has_value()); } +TEST_CASE("AnalysisMode_CBORGridCrystals", "[AnalysisMode][CBOR]") { + std::vector buffer(1024 * 1024); + CBORStream2Serializer serializer(buffer.data(), buffer.size()); + + EndMessage message{}; + message.max_image_number = 0; + // Two crystals, not one: a raster can find any number and nothing may assume at most one. + message.grid_crystals.push_back(GridScanCrystal{ + .nx = 3.5f, .ny = 7.25f, .x_um = -12.5f, .y_um = 40.0f, .image_number = 143, + .major_um = 22.5f, .minor_um = 8.0f, .angle_deg = 179.5f, .score = 0.92f, + .ice_score = 0.03f, .res_A = 1.85f, .n_images = 17}); + message.grid_crystals.push_back(GridScanCrystal{ + .nx = 11.0f, .ny = 2.0f, .x_um = 60.0f, .y_um = -30.0f, .image_number = 44, + .major_um = 9.0f, .minor_um = 7.5f, .angle_deg = 0.5f, .score = 0.41f, + .ice_score = 0.30f, .res_A = 3.2f, .n_images = 4}); + REQUIRE_NOTHROW(serializer.SerializeSequenceEnd(message)); + + auto deserialized = CBORStream2Deserialize(buffer.data(), serializer.GetBufferSize()); + REQUIRE(deserialized); + REQUIRE(deserialized->end_message); + const auto &out = deserialized->end_message->grid_crystals; + REQUIRE(out.size() == 2); + CHECK(out[0].nx == Catch::Approx(3.5)); + CHECK(out[0].image_number == 143); + CHECK(out[0].major_um == Catch::Approx(22.5)); + CHECK(out[0].angle_deg == Catch::Approx(179.5)); + CHECK(out[0].res_A == Catch::Approx(1.85)); + CHECK(out[0].n_images == 17); + CHECK(out[1].score == Catch::Approx(0.41)); + CHECK(out[1].ice_score == Catch::Approx(0.30)); + CHECK(out[1].y_um == Catch::Approx(-30.0)); +} + TEST_CASE("AnalysisMode_HDF5MasterRoundTrip", "[AnalysisMode][HDF5][Full]") { - // The mode is dataset-wide metadata: written to the master and read back from it, so a stored - // file re-opens knowing what produced it. + // The mode and the crystal list are dataset-wide metadata: written to the master, read back from + // it, so a stored raster re-opens knowing what produced it and what it found. DiffractionExperiment x(DetJF(1)); x.FilePrefix("test_analysis_mode").ImagesPerTrigger(1).OverwriteExistingFiles(true); x.BeamX_pxl(100).BeamY_pxl(200).DetectorDistance_mm(150) @@ -160,6 +193,14 @@ TEST_CASE("AnalysisMode_HDF5MasterRoundTrip", "[AnalysisMode][HDF5][Full]") { EndMessage end_message; end_message.max_image_number = 1; + end_message.grid_crystals.push_back(GridScanCrystal{ + .nx = 3.5f, .ny = 7.25f, .x_um = -12.5f, .y_um = 40.0f, .image_number = 143, + .major_um = 22.5f, .minor_um = 8.0f, .angle_deg = 179.5f, .score = 0.92f, + .ice_score = 0.03f, .res_A = 1.85f, .n_images = 17}); + end_message.grid_crystals.push_back(GridScanCrystal{ + .nx = 11.0f, .ny = 2.0f, .x_um = 60.0f, .y_um = -30.0f, .image_number = 44, + .major_um = 9.0f, .minor_um = 7.5f, .angle_deg = 0.5f, .score = 0.41f, + .ice_score = 0.30f, .res_A = 3.2f, .n_images = 4}); REQUIRE_NOTHROW(file_set.WriteHDF5(end_message)); file_set.Finalize(); } @@ -169,6 +210,13 @@ TEST_CASE("AnalysisMode_HDF5MasterRoundTrip", "[AnalysisMode][HDF5][Full]") { auto dataset = reader.GetDataset(); REQUIRE(dataset->file_analysis_mode.has_value()); CHECK(*dataset->file_analysis_mode == AnalysisMode::Grid); + REQUIRE(dataset->grid_crystals.size() == 2); + CHECK(dataset->grid_crystals[0].nx == Catch::Approx(3.5)); + CHECK(dataset->grid_crystals[0].image_number == 143); + CHECK(dataset->grid_crystals[0].angle_deg == Catch::Approx(179.5)); + CHECK(dataset->grid_crystals[0].n_images == 17); + CHECK(dataset->grid_crystals[1].score == Catch::Approx(0.41)); + CHECK(dataset->grid_crystals[1].res_A == Catch::Approx(3.2)); } remove("test_analysis_mode_master.h5"); remove("test_analysis_mode_data_000001.h5"); diff --git a/writer/HDF5NXmx.cpp b/writer/HDF5NXmx.cpp index cbcfa2b34..38ae883a8 100644 --- a/writer/HDF5NXmx.cpp +++ b/writer/HDF5NXmx.cpp @@ -1171,6 +1171,47 @@ void NXmx::Finalize(const EndMessage &end) { if (end.rotation_lattice_type) SaveScalar(*hdf5_file, "/entry/MX/rotationLatticeNiggliClass", end.rotation_lattice_type->niggli_class); + // The crystals a grid scan found: dataset-wide metadata, so it goes in the master. Written as + // parallel 1-D datasets of length N under /entry/MX/crystals, one per field, because crystals + // are numbered rather than named - a group per crystal would put the index in the path and make + // the list unreadable without walking it. Sorted by score, best first, as it arrives. + if (!end.grid_crystals.empty()) { + const size_t n = end.grid_crystals.size(); + std::vector nx(n), ny(n), x_um(n), y_um(n), major_um(n), minor_um(n), + angle_deg(n), score(n), ice_score(n), res_A(n); + std::vector image_number(n), n_images(n); + for (size_t i = 0; i < n; i++) { + const auto &c = end.grid_crystals[i]; + nx[i] = c.nx; + ny[i] = c.ny; + x_um[i] = c.x_um; + y_um[i] = c.y_um; + image_number[i] = c.image_number; + major_um[i] = c.major_um; + minor_um[i] = c.minor_um; + angle_deg[i] = c.angle_deg; + score[i] = c.score; + ice_score[i] = c.ice_score; + res_A[i] = c.res_A; + n_images[i] = c.n_images; + } + HDF5Group(*hdf5_file, "/entry/MX/crystals").NXClass("NXcollection"); + SaveVector(*hdf5_file, "/entry/MX/crystals/nx", nx); + SaveVector(*hdf5_file, "/entry/MX/crystals/ny", ny); + SaveVector(*hdf5_file, "/entry/MX/crystals/xUm", x_um)->Units("um"); + SaveVector(*hdf5_file, "/entry/MX/crystals/yUm", y_um)->Units("um"); + SaveVector(*hdf5_file, "/entry/MX/crystals/imageNumber", image_number); + SaveVector(*hdf5_file, "/entry/MX/crystals/majorUm", major_um)->Units("um"); + SaveVector(*hdf5_file, "/entry/MX/crystals/minorUm", minor_um)->Units("um"); + // An axis, not a direction: in [0, 180), so 179 deg and 0 deg are adjacent and these do + // not average arithmetically. + SaveVector(*hdf5_file, "/entry/MX/crystals/angleDeg", angle_deg)->Units("deg"); + SaveVector(*hdf5_file, "/entry/MX/crystals/score", score); + SaveVector(*hdf5_file, "/entry/MX/crystals/iceScore", ice_score); + SaveVector(*hdf5_file, "/entry/MX/crystals/resA", res_A)->Units("Angstrom"); + SaveVector(*hdf5_file, "/entry/MX/crystals/nImages", n_images); + } + // The setting the per-image reflections and lattices were written in is the one they were // indexed in, which is not always the setting of the cell above: the space group is chosen // after the images have gone to file, and choosing it can re-seat the lattice. Write the