From 550f7c7d7b118de8c2b3088342dc9f73abb5d40d Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 8 Sep 2026 00:17:10 +0200 Subject: [PATCH] scan result: a grid scan's crystals and a calibration's geometry reach the stream, the file and the API A grid scan and a powder calibration each produce an answer that is not a per-image number, and neither had anywhere to go. ScanResult gains two top-level optionals, grid and calibration, in the same idiom the three rotation_* members already use: a mode's payload is present when that mode ran and absent otherwise. images stays required, so every existing /result/scan client - the python client and the beamline scripts on it - is untouched. No oneOf or discriminator. The spec contains none today, so polymorphism is unexercised across all three generators at once, and the only thing it would buy is making "exactly one payload" structurally true rather than true by convention. powder_calibration_output has been fully specified and referenced by nothing since it was written; the calibration member reuses it. The crystal list travels the same three hops a per-image quantity does. In the CBOR end message grid_crystals is an array of maps keyed by field name - the shape a spot already uses - so a reader that does not know a field skips it and the record survives a field being added; the deserializer steps over an unknown key rather than reading its payload as the next one. In the HDF5 master it is parallel 1-D datasets under /entry/MX/crystals, one per field, because crystals are numbered rather than named and a group per crystal would put the index in the path. The reader reads them back so a stored raster re-opens with what it found. Nothing caps the list at one. A raster can find several crystals and the ordering carries the answer: sorted by score, best first. angle_deg is an axis, not a direction. It lies in [0, 180) and wraps, so 179 and 0 are adjacent and two of these do not average arithmetically - a pair of nearly parallel needles averaged across the wrap comes out as a right angle. That is stated at every layer it passes through, since it is the kind of field a consumer will reach for without reading the producer. The producer side is not wired yet: nothing fills ScanResult::grid or EndMessage::grid_crystals here. The seam is receiver/JFJochReceiver.cpp, where rotation_lattice is filled in EndDataCollection and GetFinalStatistics. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N --- broker/OpenAPIConvert.cpp | 99 +++++++++++++++++++++ broker/OpenAPIConvert.h | 5 ++ broker/jfjoch_api.yaml | 4 + common/JFJochMessages.h | 4 + common/ScanResult.h | 9 ++ docs/CBOR.md | 10 +++ docs/CHANGELOG.md | 1 + docs/HDF5.md | 17 ++++ frame_serialize/CBORStream2Deserializer.cpp | 36 +++++++- frame_serialize/CBORStream2Serializer.cpp | 35 ++++++++ reader/HDF5MetadataSource.cpp | 41 +++++++++ reader/JFJochReaderDataset.h | 4 + tests/AnalysisSettingsTest.cpp | 52 ++++++++++- writer/HDF5NXmx.cpp | 41 +++++++++ 14 files changed, 355 insertions(+), 3 deletions(-) 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