leonarski_fandjungfrau
4dc2534dbf
Build Packages / build:rpm (rocky9_sls9) (push) Successful in 18m57s
Build Packages / Unit tests (push) Skipped
Build Packages / build:windows:nocuda (push) Successful in 16m55s
Build Packages / build:windows:cuda (push) Successful in 18m48s
Build Packages / build:viewer-tgz:cpu (push) Successful in 13m10s
Build Packages / build:viewer-tgz:cuda (push) Successful in 14m45s
Build Packages / build:rpm (rocky8_nocuda) (push) Successful in 22m23s
Build Packages / build:rpm (rocky9_nocuda) (push) Successful in 20m12s
Build Packages / build:rpm (ubuntu2204_nocuda) (push) Successful in 23m7s
Build Packages / build:rpm (ubuntu2404_nocuda) (push) Successful in 20m43s
Build Packages / build:rpm (rocky8_sls9) (push) Successful in 23m9s
Build Packages / XDS test (durin plugin) (push) Successful in 12m26s
Build Packages / build:rpm (rocky9) (push) Successful in 24m58s
Build Packages / Generate python client (push) Successful in 50s
Build Packages / build:rpm (ubuntu2404) (push) Successful in 23m20s
Build Packages / Create release (push) Skipped
Build Packages / XDS test (JFJoch plugin) (push) Successful in 12m37s
Build Packages / build:rpm (rocky8) (push) Successful in 27m58s
Build Packages / build:rpm (ubuntu2204) (push) Successful in 25m38s
Build Packages / Build documentation (push) Successful in 59s
Build Packages / DIALS test (push) Successful in 23m16s
Build Packages / XDS test (neggia plugin) (push) Successful in 6m38s
**Files written by Jungfraujoch now import correctly in DIALS, XDS and pyFAI.** A tilted detector, a grid scan, a still recorded at a goniometer position, and saturated or unreadable pixels were each described in a way that a third-party program acted on wrongly. If you process Jungfraujoch data outside Jungfraujoch, prefer this release to any earlier one. * HDF5: the detector tilt (`rot1`/`rot2`/`rot3`) is exported correctly in the NXmx transformation chain; untilted geometries are unaffected. * HDF5: a still recorded at a goniometer position is no longer read back as a single image, and a grid scan records a stationary spindle so a program that requires a rotation axis can open it. * HDF5: the sample transformation chain is written in mounting order, with a Smargon head position told apart from the spindle, one entry per image, `module_offset` as a float unit vector, and `offset_units` on every offset. * HDF5: saturated, underloaded and unreadable pixels are described so a downstream program masks them - `saturation_value`, `underload_value`, `error_value` and `bit_depth_readout` are written correctly, and a data file missing next to a VDS master reads as the error marker rather than as zero counts. * HDF5: the rotation axis is read back under whatever name it carries, and `mirror_y` records whether the assembled image is mirrored in Y relative to the detector's raw readout. * A grid scan and a goniometer axis can both be set; they are no longer alternatives. * `images_per_file` is chosen from the acquisition when it is not given: a rotation sweep of at most 20000 images goes into a single data file, a grid scan splits on whole fast-axis rows, and stills and serial keep 1000. * The writer refuses a stream whose start message declares a different pixel format than its images carry, and a DECTRIS detector sending signed images is no longer declared unsigned. * The image stream can carry the sample transformation chain (`transformations`, in the END message); a producer that does not send it gets the same chain built by the writer. * rugnux: fixing the space group with `-S` no longer prevents the lattice from being found - a lattice indexed in a different setting is reindexed into that group's own setting, and a run whose crystal does not have that group's lattice stops and names the cell it indexed as, rather than reporting statistics that cannot describe it. * rugnux: the per-image resolution estimate now predicts the resolution the merged data reach rather than the highest-resolution spot found, and is reported as `SPOT_RESOLUTION_ESTIMATE`. * rugnux: two runs of the same command on the same images produce the same merged intensities; the azimuthal profile written alongside them is not yet reproducible in the same way. * rugnux: the offline lattice refinement is bounded by iterations rather than by a wall clock, so a loaded machine can no longer refine to a different lattice; a live acquisition keeps its real-time bound. * rugnux: the detector-frame modulation correction is fitted on a grid spanning the detector, so whether it is applied no longer depends on how far integration reached. * rugnux: the geometry pre-pass no longer writes `<prefix>_01.mtz`, `_01.cif`, `_01.hkl` and `_01_image.dat`; the refined second pass writes those files under `<prefix>`, and that is the result to use. * rugnux: `_process.h5` describes the pixel format of the images it links to, and is written on a thread of its own. * rugnux: the detector geometry is also logged in XDS's convention (`ORGX`/`ORGY`, detector axis vectors, rotation axis), so it can be compared with an XDS refinement. * rugnux: an image integrated in pyFAI through the `.poni` file written by `--mode calibration` comes out with the correct azimuth, and the file declares pyFAI's `orientation`, which needs pyFAI 2024.01 or newer. Radial integration is unchanged. * rugnux: a rotation run is substantially faster throughout - beam-stop detection, first-pass indexing, geometry refinement, integration, scaling and merging - and observations outside the scaling resolution range are dropped as they are ingested. The refined geometry, the space group chosen and the merged statistics are unchanged. * Faster spot finding and indexing, on the broker as well as in rugnux; the spots found and the lattices indexed are unchanged. * A run reserves substantially less GPU memory: nothing is allocated for buffers that are never read, and a worker builds only the engines it uses. * rugnux: with `-N` left at its default the per-image loop of `--mode mx` uses at most 16 workers per GPU, rather than one per hardware thread; an explicit `-N` is obeyed as given. * CUDA 12 builds now contain device code for Volta, so the RHEL 8 packages and the portable Linux `.tgz` run on a V100; the CUDA 13 artefacts (RHEL 9, Ubuntu, Windows) remain Turing and newer. * The build resolves a single Eigen for the whole project, and refuses to configure if Ceres picks up a different one; a build that mixed two Eigen versions was undefined behaviour and crashed at -O2. * Documentation: a security page, and the supported GPU generations and minimum NVIDIA driver version of every released artefact. **Breaking change to OpenAPI** - regenerate the client (`jfjoch-client` 1.0.0-rc.162, `frontend/src/client`): * `dataset_settings.images_per_file` is no longer `default: 1000` and no longer accepts `0`; it is optional, and its minimum is 1. A client sending `0` (previously "one file for the whole run") is now rejected - omit the field instead, which for a rotation sweep gives the same single file. * `file_writer_format` now defaults to `NXmxVDS`, matching the server's own default and the layout recommended for DIALS, XDS and CrystFEL. A generated client that fills in schema defaults and does not set the format explicitly will write VDS masters where it previously wrote legacy ones; set `NXmxLegacy` explicitly to keep them. --------- Co-authored-by: jungfrau <jungfrau@mx-aare-test.psi.ch> Reviewed-on: #72 Co-authored-by: Filip Leonarski <filip.leonarski@psi.ch>
674 lines
27 KiB
C++
674 lines
27 KiB
C++
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include <catch2/catch_all.hpp>
|
|
#include <iostream>
|
|
|
|
#include "../image_analysis/geom_refinement/XtalOptimizer.h"
|
|
#include "../image_analysis/bragg_prediction/BraggPrediction.h"
|
|
|
|
TEST_CASE("XtalOptimizer") {
|
|
DiffractionExperiment exp_i;
|
|
exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV)
|
|
.BeamX_pxl(1000)
|
|
.BeamY_pxl(1000)
|
|
.PoniRot1_rad(0.01)
|
|
.PoniRot2_rad(0.02)
|
|
.DetectorDistance_mm(200);
|
|
|
|
CrystalLattice latt_i(40,40,80,90,90,90);
|
|
|
|
BraggPredictionSettings prediction_settings{
|
|
.high_res_A = 1.5,
|
|
.ewald_dist_cutoff = 0.001};
|
|
BraggPrediction prediction;
|
|
auto count = prediction.Calc(exp_i, latt_i, prediction_settings);
|
|
|
|
std::vector<SpotToSave> spots;
|
|
|
|
for (int i = 0; i < count; ++i) {
|
|
auto refl = prediction.GetReflections().at(i);
|
|
spots.push_back(SpotToSave{refl.predicted_x, refl.predicted_y});
|
|
}
|
|
|
|
XtalOptimizerData xtal_opt{};
|
|
xtal_opt.max_time = 30.0;
|
|
xtal_opt.latt = CrystalLattice(40.2,39.4,80.2, 90,91, 89);
|
|
xtal_opt.geom.BeamX_pxl(1010).BeamY_pxl(995).DetectorDistance_mm(200)
|
|
.PoniRot1_rad(0.01).PoniRot2_rad(0.02);
|
|
xtal_opt.crystal_system = gemmi::CrystalSystem::Triclinic;
|
|
|
|
auto start = std::chrono::high_resolution_clock::now();
|
|
REQUIRE(XtalOptimizer(xtal_opt, {spots}));
|
|
auto end = std::chrono::high_resolution_clock::now();
|
|
std::cout << "XtalOptimizer took " << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
|
|
<< " microseconds" << std::endl;
|
|
|
|
auto uc_i = latt_i.GetUnitCell();
|
|
auto uc_o = xtal_opt.latt.GetUnitCell();
|
|
|
|
std::cout << "Beam center: " << xtal_opt.geom.GetBeamX_pxl() << " " << xtal_opt.geom.GetBeamY_pxl() << std::endl;
|
|
std::cout << "Unit cell: " << uc_o.a << " " << uc_o.b << " " << uc_o.c << std::endl;
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamX_pxl() - exp_i.GetBeamX_pxl()) < 0.05);
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamY_pxl() - exp_i.GetBeamY_pxl()) < 0.05);
|
|
|
|
CHECK(fabsf(uc_i.a - uc_o.a) < 0.1);
|
|
CHECK(fabsf(uc_i.b - uc_o.b) < 0.1);
|
|
CHECK(fabsf(uc_i.c - uc_o.c) < 0.2);
|
|
CHECK(fabsf(uc_i.alpha - uc_o.alpha) < 0.1);
|
|
CHECK(fabsf(uc_i.beta - uc_o.beta) < 0.1);
|
|
CHECK(fabsf(uc_i.gamma - uc_o.gamma) < 0.1);
|
|
}
|
|
|
|
|
|
TEST_CASE("XtalOptimizer_NoBeamCenter") {
|
|
DiffractionExperiment exp_i;
|
|
exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV)
|
|
.BeamX_pxl(1000)
|
|
.BeamY_pxl(1000)
|
|
.PoniRot1_rad(0.01)
|
|
.PoniRot2_rad(0.02)
|
|
.DetectorDistance_mm(200);
|
|
|
|
CrystalLattice latt_i(40,50,80,90,95,90);
|
|
|
|
BraggPredictionSettings prediction_settings{
|
|
.high_res_A = 1.5,
|
|
.ewald_dist_cutoff = 0.001};
|
|
BraggPrediction prediction;
|
|
auto count = prediction.Calc(exp_i, latt_i, prediction_settings);
|
|
|
|
std::vector<SpotToSave> spots;
|
|
|
|
for (int i = 0; i < count; ++i) {
|
|
auto refl = prediction.GetReflections().at(i);
|
|
spots.push_back(SpotToSave{refl.predicted_x, refl.predicted_y});
|
|
}
|
|
|
|
XtalOptimizerData xtal_opt{};
|
|
xtal_opt.latt = CrystalLattice(40.2,49.4,80.2, 90,94, 89);
|
|
xtal_opt.geom.BeamX_pxl(999.8).BeamY_pxl(1000.2).DetectorDistance_mm(200)
|
|
.PoniRot1_rad(0.01).PoniRot2_rad(0.02);
|
|
xtal_opt.crystal_system = gemmi::CrystalSystem::Monoclinic;
|
|
xtal_opt.refine_beam_center = false;
|
|
xtal_opt.max_time = 30.0;
|
|
|
|
auto start = std::chrono::high_resolution_clock::now();
|
|
REQUIRE(XtalOptimizer(xtal_opt, {spots}));
|
|
auto end = std::chrono::high_resolution_clock::now();
|
|
std::cout << "XtalOptimizer took " << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
|
|
<< " microseconds" << std::endl;
|
|
|
|
auto uc_i = latt_i.GetUnitCell();
|
|
auto uc_o = xtal_opt.latt.GetUnitCell();
|
|
|
|
std::cout << "Beam center: " << xtal_opt.geom.GetBeamX_pxl() << " " << xtal_opt.geom.GetBeamY_pxl() << std::endl;
|
|
std::cout << "Unit cell: " << uc_o.a << " " << uc_o.b << " " << uc_o.c << std::endl;
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamX_pxl() - 999.8) < 0.01);
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamY_pxl() - 1000.2) < 0.01);
|
|
|
|
CHECK(fabsf(uc_i.a - uc_o.a) < 0.1);
|
|
CHECK(fabsf(uc_i.b - uc_o.b) < 0.1);
|
|
CHECK(fabsf(uc_i.c - uc_o.c) < 0.2);
|
|
CHECK(fabsf(uc_i.alpha - uc_o.alpha) < 0.1);
|
|
CHECK(fabsf(uc_i.beta - uc_o.beta) < 0.1);
|
|
CHECK(fabsf(uc_i.gamma - uc_o.gamma) < 0.1);
|
|
}
|
|
|
|
TEST_CASE("XtalOptimizer_orthorombic") {
|
|
DiffractionExperiment exp_i;
|
|
exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV)
|
|
.BeamX_pxl(1000)
|
|
.BeamY_pxl(1000)
|
|
.PoniRot1_rad(0.01)
|
|
.PoniRot2_rad(0.02)
|
|
.DetectorDistance_mm(200);
|
|
|
|
CrystalLattice latt_i(40,50,80,90,90,90);
|
|
|
|
BraggPredictionSettings prediction_settings{
|
|
.high_res_A = 1.5,
|
|
.ewald_dist_cutoff = 0.001};
|
|
BraggPrediction prediction;
|
|
auto count = prediction.Calc(exp_i, latt_i, prediction_settings);
|
|
|
|
std::vector<SpotToSave> spots;
|
|
|
|
for (int i = 0; i < count; ++i) {
|
|
auto refl = prediction.GetReflections().at(i);
|
|
spots.push_back(SpotToSave{refl.predicted_x, refl.predicted_y});
|
|
}
|
|
XtalOptimizerData xtal_opt{};
|
|
xtal_opt.latt = CrystalLattice(40.2,49.6,80.3, 90,91, 89);
|
|
xtal_opt.geom.BeamX_pxl(1005).BeamY_pxl(997).DetectorDistance_mm(200)
|
|
.PoniRot1_rad(0.01).PoniRot2_rad(0.02);
|
|
xtal_opt.max_time = 30.0;
|
|
xtal_opt.crystal_system = gemmi::CrystalSystem::Orthorhombic;
|
|
auto start = std::chrono::high_resolution_clock::now();
|
|
REQUIRE(XtalOptimizer(xtal_opt, {spots}));
|
|
auto end = std::chrono::high_resolution_clock::now();
|
|
std::cout << "XtalOptimizer took " << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
|
|
<< " microseconds" << std::endl;
|
|
|
|
auto uc_i = latt_i.GetUnitCell();
|
|
auto uc_o = xtal_opt.latt.GetUnitCell();
|
|
|
|
std::cout << "Beam center: " << xtal_opt.geom.GetBeamX_pxl() << " " << xtal_opt.geom.GetBeamY_pxl() << std::endl;
|
|
std::cout << "Unit cell: " << uc_o.a << " " << uc_o.b << " " << uc_o.c << std::endl;
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamX_pxl() - exp_i.GetBeamX_pxl()) < 0.1);
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamY_pxl() - exp_i.GetBeamY_pxl()) < 0.1);
|
|
|
|
CHECK(fabsf(uc_i.a - uc_o.a) < 0.1);
|
|
CHECK(fabsf(uc_i.b - uc_o.b) < 0.1);
|
|
CHECK(fabsf(uc_i.c - uc_o.c) < 0.2);
|
|
CHECK(fabs(uc_o.alpha - 90) < 0.02);
|
|
CHECK(fabs(uc_o.beta - 90) < 0.02);
|
|
CHECK(fabs(uc_o.gamma - 90) < 0.02);
|
|
}
|
|
|
|
TEST_CASE("XtalOptimizer_triclinic") {
|
|
DiffractionExperiment exp_i;
|
|
exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV)
|
|
.BeamX_pxl(1000)
|
|
.BeamY_pxl(1000)
|
|
.PoniRot1_rad(0.01)
|
|
.PoniRot2_rad(0.02)
|
|
.DetectorDistance_mm(200);
|
|
|
|
CrystalLattice latt_i(40,55,120,95,97,100);
|
|
|
|
BraggPredictionSettings prediction_settings{
|
|
.high_res_A = 1.5,
|
|
.ewald_dist_cutoff = 0.001,
|
|
};
|
|
BraggPrediction prediction;
|
|
auto count = prediction.Calc(exp_i, latt_i, prediction_settings);
|
|
|
|
std::vector<SpotToSave> spots;
|
|
|
|
for (int i = 0; i < count; ++i) {
|
|
auto refl = prediction.GetReflections().at(i);
|
|
spots.push_back(SpotToSave{refl.predicted_x, refl.predicted_y});
|
|
}
|
|
XtalOptimizerData xtal_opt{};
|
|
xtal_opt.latt = CrystalLattice(40.1,54.9,121, 95,97, 99.5);
|
|
xtal_opt.geom.BeamX_pxl(997).BeamY_pxl(1005).DetectorDistance_mm(200)
|
|
.PoniRot1_rad(0.01).PoniRot2_rad(0.02);
|
|
xtal_opt.crystal_system = gemmi::CrystalSystem::Triclinic;
|
|
xtal_opt.max_time = 36.0;
|
|
auto start = std::chrono::high_resolution_clock::now();
|
|
REQUIRE(XtalOptimizer(xtal_opt, {spots}));
|
|
auto end = std::chrono::high_resolution_clock::now();
|
|
std::cout << "XtalOptimizer took " << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
|
|
<< " microseconds" << std::endl;
|
|
|
|
auto uc_i = latt_i.GetUnitCell();
|
|
auto uc_o = xtal_opt.latt.GetUnitCell();
|
|
|
|
std::cout << "Beam center: " << xtal_opt.geom.GetBeamX_pxl() << " " << xtal_opt.geom.GetBeamY_pxl() << std::endl;
|
|
std::cout << "Unit cell: " << uc_o.a << " " << uc_o.b << " " << uc_o.c << std::endl;
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamX_pxl() - exp_i.GetBeamX_pxl()) < 0.2);
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamY_pxl() - exp_i.GetBeamY_pxl()) < 0.2);
|
|
|
|
CHECK(fabsf(uc_i.a - uc_o.a) < 0.1);
|
|
CHECK(fabsf(uc_i.b - uc_o.b) < 0.1);
|
|
CHECK(fabsf(uc_i.c - uc_o.c) < 0.5);
|
|
CHECK(fabsf(uc_i.alpha - uc_o.alpha) < 0.1);
|
|
CHECK(fabsf(uc_i.beta - uc_o.beta) < 0.1);
|
|
CHECK(fabsf(uc_i.gamma - uc_o.gamma) < 0.1);
|
|
}
|
|
|
|
|
|
TEST_CASE("XtalOptimizer_tetragonal") {
|
|
DiffractionExperiment exp_i;
|
|
exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV)
|
|
.BeamX_pxl(1000)
|
|
.BeamY_pxl(1000)
|
|
.PoniRot1_rad(0.01)
|
|
.PoniRot2_rad(0.02)
|
|
.DetectorDistance_mm(200);
|
|
|
|
CrystalLattice latt_i(40,40,80,90,90,90);
|
|
|
|
BraggPredictionSettings prediction_settings{
|
|
.high_res_A = 1.5,
|
|
.ewald_dist_cutoff = 0.001};
|
|
BraggPrediction prediction;
|
|
auto count = prediction.Calc(exp_i, latt_i, prediction_settings);
|
|
|
|
std::vector<SpotToSave> spots;
|
|
|
|
for (int i = 0; i < count; ++i) {
|
|
auto refl = prediction.GetReflections().at(i);
|
|
spots.push_back(SpotToSave{refl.predicted_x, refl.predicted_y});
|
|
}
|
|
|
|
XtalOptimizerData xtal_opt{};
|
|
xtal_opt.latt = CrystalLattice(40.6,39.3,80.5, 90,91, 89);
|
|
xtal_opt.geom.BeamX_pxl(1010).BeamY_pxl(995).DetectorDistance_mm(200)
|
|
.PoniRot1_rad(0.01).PoniRot2_rad(0.02);
|
|
xtal_opt.crystal_system = gemmi::CrystalSystem::Tetragonal;
|
|
xtal_opt.max_time = 30.0;
|
|
auto start = std::chrono::high_resolution_clock::now();
|
|
REQUIRE(XtalOptimizer(xtal_opt, {spots}));
|
|
auto end = std::chrono::high_resolution_clock::now();
|
|
std::cout << "XtalOptimizer took " << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
|
|
<< " microseconds" << std::endl;
|
|
|
|
auto uc_i = latt_i.GetUnitCell();
|
|
auto uc_o = xtal_opt.latt.GetUnitCell();
|
|
|
|
std::cout << "Beam center: " << xtal_opt.geom.GetBeamX_pxl() << " " << xtal_opt.geom.GetBeamY_pxl() << std::endl;
|
|
std::cout << "Unit cell: " << uc_o.a << " " << uc_o.b << " " << uc_o.c << std::endl;
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamX_pxl() - exp_i.GetBeamX_pxl()) < 0.1);
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamY_pxl() - exp_i.GetBeamY_pxl()) < 0.1);
|
|
|
|
CHECK(fabsf(uc_i.a - uc_o.a) < 0.1);
|
|
CHECK(fabsf(uc_i.b - uc_o.b) < 0.1);
|
|
CHECK(fabsf(uc_i.c - uc_o.c) < 0.5);
|
|
CHECK(fabs(uc_o.alpha - 90) < 0.02);
|
|
CHECK(fabs(uc_o.beta - 90) < 0.02);
|
|
CHECK(fabs(uc_o.gamma - 90) < 0.02);
|
|
}
|
|
|
|
TEST_CASE("XtalOptimizer_hexagonal") {
|
|
DiffractionExperiment exp_i;
|
|
exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV)
|
|
.BeamX_pxl(1000)
|
|
.BeamY_pxl(1000)
|
|
.PoniRot1_rad(0.01)
|
|
.PoniRot2_rad(0.02)
|
|
.DetectorDistance_mm(200);
|
|
|
|
CrystalLattice latt_i(40,40,70,90,90,120);
|
|
|
|
BraggPredictionSettings prediction_settings{
|
|
.high_res_A = 1.5,
|
|
.ewald_dist_cutoff = 0.001};
|
|
BraggPrediction prediction;
|
|
auto count = prediction.Calc(exp_i, latt_i, prediction_settings);
|
|
|
|
std::vector<SpotToSave> spots;
|
|
|
|
for (int i = 0; i < count; ++i) {
|
|
auto refl = prediction.GetReflections().at(i);
|
|
spots.push_back(SpotToSave{refl.predicted_x, refl.predicted_y});
|
|
}
|
|
|
|
XtalOptimizerData xtal_opt{};
|
|
xtal_opt.latt = CrystalLattice(39.5,39.8,70.1, 90,90, 119.5);
|
|
xtal_opt.geom.BeamX_pxl(1007).BeamY_pxl(990).DetectorDistance_mm(200)
|
|
.PoniRot1_rad(0.01).PoniRot2_rad(0.02);
|
|
xtal_opt.crystal_system = gemmi::CrystalSystem::Hexagonal;
|
|
xtal_opt.max_time = 60.0;
|
|
auto start = std::chrono::high_resolution_clock::now();
|
|
bool ret = XtalOptimizer(xtal_opt, {spots});
|
|
auto end = std::chrono::high_resolution_clock::now();
|
|
std::cout << "XtalOptimizer took " << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
|
|
<< " microseconds" << std::endl;
|
|
REQUIRE(ret);
|
|
auto uc_i = latt_i.GetUnitCell();
|
|
auto uc_o = xtal_opt.latt.GetUnitCell();
|
|
|
|
std::cout << "Beam center: " << xtal_opt.geom.GetBeamX_pxl() << " " << xtal_opt.geom.GetBeamY_pxl() << std::endl;
|
|
std::cout << "Unit cell: " << uc_o.a << " " << uc_o.b << " " << uc_o.c << " " << uc_o.alpha << " " << uc_o.beta
|
|
<< " " << uc_o.gamma << std::endl;
|
|
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamX_pxl() - exp_i.GetBeamX_pxl()) < 0.1);
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamY_pxl() - exp_i.GetBeamY_pxl()) < 0.1);
|
|
|
|
CHECK(fabsf(uc_i.a - uc_o.a) < 0.1);
|
|
CHECK(fabsf(uc_i.b - uc_o.b) < 0.1);
|
|
CHECK(fabsf(uc_i.c - uc_o.c) < 0.2);
|
|
CHECK(fabs(uc_o.alpha - 90) < 0.02);
|
|
CHECK(fabs(uc_o.beta - 90) < 0.01);
|
|
CHECK(fabs(uc_o.gamma - 120) < 0.01);
|
|
}
|
|
|
|
TEST_CASE("XtalOptimizer_hexagonal_unconstrained") {
|
|
DiffractionExperiment exp_i;
|
|
exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV)
|
|
.BeamX_pxl(1000)
|
|
.BeamY_pxl(1000)
|
|
.PoniRot1_rad(0.01)
|
|
.PoniRot2_rad(0.02)
|
|
.DetectorDistance_mm(200);
|
|
|
|
CrystalLattice latt_i(40,40,70,90,90,120);
|
|
|
|
BraggPredictionSettings prediction_settings{
|
|
.high_res_A = 1.5,
|
|
.ewald_dist_cutoff = 0.001};
|
|
BraggPrediction prediction;
|
|
auto count = prediction.Calc(exp_i, latt_i, prediction_settings);
|
|
|
|
std::vector<SpotToSave> spots;
|
|
|
|
for (int i = 0; i < count; ++i) {
|
|
auto refl = prediction.GetReflections().at(i);
|
|
spots.push_back(SpotToSave{refl.predicted_x, refl.predicted_y});
|
|
}
|
|
|
|
XtalOptimizerData xtal_opt{};
|
|
xtal_opt.latt = CrystalLattice(39.9,39.8,70.1, 90,90, 120);
|
|
xtal_opt.geom.BeamX_pxl(1002).BeamY_pxl(998).DetectorDistance_mm(200)
|
|
.PoniRot1_rad(0.01).PoniRot2_rad(0.02);
|
|
xtal_opt.crystal_system = gemmi::CrystalSystem::Triclinic;
|
|
xtal_opt.max_time = 30.0;
|
|
auto start = std::chrono::high_resolution_clock::now();
|
|
REQUIRE(XtalOptimizer(xtal_opt, {spots}));
|
|
auto end = std::chrono::high_resolution_clock::now();
|
|
std::cout << "XtalOptimizer took " << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
|
|
<< " microseconds" << std::endl;
|
|
|
|
auto uc_i = latt_i.GetUnitCell();
|
|
auto uc_o = xtal_opt.latt.GetUnitCell();
|
|
|
|
std::cout << "Beam center: " << xtal_opt.geom.GetBeamX_pxl() << " " << xtal_opt.geom.GetBeamY_pxl() << std::endl;
|
|
std::cout << "Unit cell: " << uc_o.a << " " << uc_o.b << " " << uc_o.c << " " << uc_o.alpha << " " << uc_o.beta
|
|
<< " " << uc_o.gamma << std::endl;
|
|
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamX_pxl() - exp_i.GetBeamX_pxl()) < 0.3);
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamY_pxl() - exp_i.GetBeamY_pxl()) < 0.3);
|
|
|
|
CHECK(fabsf(uc_i.a - uc_o.a) < 0.1);
|
|
CHECK(fabsf(uc_i.b - uc_o.b) < 0.1);
|
|
CHECK(fabsf(uc_i.c - uc_o.c) < 0.2);
|
|
CHECK(fabs(uc_o.alpha - 90) < 0.1);
|
|
CHECK(fabs(uc_o.beta - 90) < 0.1);
|
|
CHECK(fabs(uc_o.gamma - 120) < 0.1);
|
|
}
|
|
|
|
TEST_CASE("XtalOptimizer_cubic") {
|
|
DiffractionExperiment exp_i;
|
|
exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV)
|
|
.BeamX_pxl(1000)
|
|
.BeamY_pxl(1000)
|
|
.PoniRot1_rad(0.01)
|
|
.PoniRot2_rad(0.02)
|
|
.DetectorDistance_mm(200);
|
|
|
|
CrystalLattice latt_i(Coord(40,0,0),
|
|
Coord(0, 40 / sqrt(2), -40 / sqrt(2)),
|
|
Coord(0, 40 / sqrt(2), 40 / sqrt(2)));
|
|
auto uc_i = latt_i.GetUnitCell();
|
|
|
|
BraggPredictionSettings prediction_settings{
|
|
.high_res_A = 1.5,
|
|
.ewald_dist_cutoff = 0.001};
|
|
BraggPrediction prediction;
|
|
auto count = prediction.Calc(exp_i, latt_i, prediction_settings);
|
|
|
|
std::vector<SpotToSave> spots;
|
|
|
|
for (int i = 0; i < count; ++i) {
|
|
auto refl = prediction.GetReflections().at(i);
|
|
spots.push_back(SpotToSave{refl.predicted_x, refl.predicted_y});
|
|
}
|
|
|
|
XtalOptimizerData xtal_opt{};
|
|
xtal_opt.latt = CrystalLattice(Coord(39,0,0),
|
|
Coord(0, 39.5 / sqrt(2), -40.5 / sqrt(2)),
|
|
Coord(0, 39.2 / sqrt(2), 39.7 / sqrt(2)));
|
|
xtal_opt.geom.BeamX_pxl(1007).BeamY_pxl(990).DetectorDistance_mm(200)
|
|
.PoniRot1_rad(0.01).PoniRot2_rad(0.02);
|
|
xtal_opt.crystal_system = gemmi::CrystalSystem::Cubic;
|
|
xtal_opt.max_time = 30.0;
|
|
auto start = std::chrono::high_resolution_clock::now();
|
|
REQUIRE(XtalOptimizer(xtal_opt, {spots}));
|
|
auto end = std::chrono::high_resolution_clock::now();
|
|
std::cout << "XtalOptimizer took " << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
|
|
<< " microseconds" << std::endl;
|
|
|
|
auto uc_o = xtal_opt.latt.GetUnitCell();
|
|
|
|
std::cout << "Beam center: " << xtal_opt.geom.GetBeamX_pxl() << " " << xtal_opt.geom.GetBeamY_pxl() << std::endl;
|
|
std::cout << "Unit cell: " << uc_o.a << " " << uc_o.b << " " << uc_o.c << std::endl;
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamX_pxl() - exp_i.GetBeamX_pxl()) < 0.1);
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamY_pxl() - exp_i.GetBeamY_pxl()) < 0.1);
|
|
|
|
CHECK(fabsf(uc_i.a - uc_o.a) < 0.1);
|
|
CHECK(fabsf(uc_i.b - uc_o.b) < 0.1);
|
|
CHECK(fabsf(uc_i.c - uc_o.c) < 0.2);
|
|
CHECK(fabs(uc_o.alpha - 90) < 0.02);
|
|
CHECK(fabs(uc_o.beta - 90) < 0.02);
|
|
CHECK(fabs(uc_o.gamma - 90) < 0.02);
|
|
}
|
|
|
|
TEST_CASE("XtalOptimizer_monoclinic") {
|
|
DiffractionExperiment exp_i;
|
|
exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV)
|
|
.BeamX_pxl(1000)
|
|
.BeamY_pxl(1000)
|
|
.PoniRot1_rad(0.01)
|
|
.PoniRot2_rad(0.02)
|
|
.DetectorDistance_mm(200);
|
|
|
|
CrystalLattice latt_i(50,60,70,90,115,90);
|
|
|
|
auto uc_i = latt_i.GetUnitCell();
|
|
|
|
BraggPredictionSettings prediction_settings{
|
|
.high_res_A = 1.5,
|
|
.ewald_dist_cutoff = 0.001
|
|
};
|
|
BraggPrediction prediction;
|
|
auto count = prediction.Calc(exp_i, latt_i, prediction_settings);
|
|
|
|
std::vector<SpotToSave> spots;
|
|
|
|
for (int i = 0; i < count; ++i) {
|
|
auto refl = prediction.GetReflections().at(i);
|
|
spots.push_back(SpotToSave{refl.predicted_x, refl.predicted_y});
|
|
}
|
|
|
|
XtalOptimizerData xtal_opt{};
|
|
xtal_opt.latt = CrystalLattice(49.5, 60.5, 69.8, 90, 113.5, 90);
|
|
xtal_opt.geom.BeamX_pxl(1007).BeamY_pxl(990).DetectorDistance_mm(200)
|
|
.PoniRot1_rad(0.01).PoniRot2_rad(0.02);
|
|
xtal_opt.crystal_system = gemmi::CrystalSystem::Monoclinic;
|
|
xtal_opt.max_time = 30.0;
|
|
auto start = std::chrono::high_resolution_clock::now();
|
|
REQUIRE(XtalOptimizer(xtal_opt, {spots}));
|
|
auto end = std::chrono::high_resolution_clock::now();
|
|
std::cout << "XtalOptimizer took " << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
|
|
<< " microseconds" << std::endl;
|
|
|
|
auto uc_o = xtal_opt.latt.GetUnitCell();
|
|
|
|
std::cout << "Beam center: " << xtal_opt.geom.GetBeamX_pxl() << " " << xtal_opt.geom.GetBeamY_pxl() << std::endl;
|
|
std::cout << "Unit cell: " << uc_o.a << " " << uc_o.b << " " << uc_o.c << " " << uc_o.alpha << " " << uc_o.beta << " " << uc_o.gamma << std::endl;
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamX_pxl() - exp_i.GetBeamX_pxl()) < 0.2);
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamY_pxl() - exp_i.GetBeamY_pxl()) < 0.2);
|
|
|
|
CHECK(fabsf(uc_i.a - uc_o.a) < 0.1);
|
|
CHECK(fabsf(uc_i.b - uc_o.b) < 0.1);
|
|
CHECK(fabsf(uc_i.c - uc_o.c) < 0.2);
|
|
CHECK(fabs(uc_o.alpha - 90) < 0.05);
|
|
CHECK(fabs(uc_o.beta - uc_i.beta) < 0.05);
|
|
CHECK(fabs(uc_o.gamma - 90) < 0.05);
|
|
}
|
|
|
|
TEST_CASE("XtalOptimizer_rotation") {
|
|
// Geometry
|
|
DiffractionExperiment exp_i;
|
|
exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV)
|
|
.BeamX_pxl(1000)
|
|
.BeamY_pxl(1000)
|
|
.PoniRot1_rad(0.01)
|
|
.PoniRot2_rad(0.02)
|
|
.DetectorDistance_mm(200);
|
|
|
|
// Base lattice (non-pathological)
|
|
CrystalLattice latt_base(40, 50, 80, 90, 95, 90);
|
|
auto uc_ref = latt_base.GetUnitCell();
|
|
|
|
// Rotation axis: around X with 3 deg per image
|
|
GoniometerAxis axis("omega", 0.0f, 3.0f, Coord(1,0,0), std::nullopt);
|
|
|
|
BraggPredictionSettings prediction_settings{
|
|
.high_res_A = 1.5,
|
|
.ewald_dist_cutoff = 0.002
|
|
};
|
|
|
|
size_t nimages = 10;
|
|
|
|
std::vector<std::vector<SpotToSave>> spots(nimages);
|
|
BraggPrediction prediction;
|
|
|
|
// Predict reflections for images at 0-30 deg.
|
|
for (int img = 0; img < nimages; ++img) {
|
|
// For a rotated image, per-image lattice is obtained as Multiply(rot.transpose())
|
|
const float angle_deg = axis.GetAngle_deg(img) + axis.GetWedge_deg() / 2.0f;
|
|
const RotMatrix rot = axis.GetTransformationAngle(angle_deg);
|
|
const CrystalLattice latt_img = latt_base.Multiply(rot.transpose());
|
|
|
|
const auto n = prediction.Calc(exp_i, latt_img, prediction_settings);
|
|
for (int i = 0; i < n; ++i) {
|
|
const auto& r = prediction.GetReflections().at(i);
|
|
SpotToSave s{};
|
|
s.x = r.predicted_x;
|
|
s.y = r.predicted_y;
|
|
s.image = img; // provide image index for rotation-aware refinement
|
|
s.phi = angle_deg;
|
|
s.intensity = 1.0f; // minimal positive value
|
|
s.ice_ring = false;
|
|
s.indexed = true;
|
|
spots[img].push_back(s);
|
|
}
|
|
}
|
|
|
|
// Seed slightly perturbed geometry and lattice; provide rotation axis for refinement.
|
|
// The beam-centre component PARALLEL to the spindle (here X, the rotation axis) is a gauge-weak
|
|
// direction that XtalOptimizer now deliberately restrains toward the header rather than refining
|
|
// freely (see the BeamComponentPrior in XtalOptimizer.cpp). So only the perpendicular component
|
|
// (Y) is seeded off-truth to exercise beam-centre recovery; X is seeded at its true value.
|
|
// The distance is seeded at its true value too: XtalOptimizer does not refine it - the rotation
|
|
// post-refinement fits it globally instead, where the cell is held fixed and the two are no longer
|
|
// degenerate - so a distance seeded off-truth here would simply be absorbed by the cell.
|
|
XtalOptimizerData xtal_opt{};
|
|
xtal_opt.latt = CrystalLattice(39.7f, 50.6f, 79.6f, 90.0f, 94.5f, 90.5f);
|
|
xtal_opt.geom.BeamX_pxl(1000).BeamY_pxl(997).DetectorDistance_mm(200.0)
|
|
.PoniRot1_rad(0.01).PoniRot2_rad(0.02);
|
|
xtal_opt.crystal_system = gemmi::CrystalSystem::Monoclinic;
|
|
xtal_opt.axis = axis;
|
|
xtal_opt.min_spots = 200;
|
|
xtal_opt.refine_beam_center = true;
|
|
xtal_opt.refine_detector_angles = false;
|
|
xtal_opt.max_time = 30.0;
|
|
|
|
auto t0 = std::chrono::high_resolution_clock::now();
|
|
REQUIRE(XtalOptimizer(xtal_opt, spots));
|
|
auto t1 = std::chrono::high_resolution_clock::now();
|
|
std::cout << "XtalOptimizer (rotation 4 images) took "
|
|
<< std::chrono::duration_cast<std::chrono::microseconds>(t1 - t0).count()
|
|
<< " microseconds" << std::endl;
|
|
|
|
const auto uc_out = xtal_opt.latt.GetUnitCell();
|
|
|
|
// Geometry checks
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamX_pxl() - exp_i.GetBeamX_pxl()) < 0.2f);
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamY_pxl() - exp_i.GetBeamY_pxl()) < 0.2f);
|
|
|
|
// Lattice checks
|
|
CHECK(fabsf(uc_ref.a - uc_out.a) < 0.2f);
|
|
CHECK(fabsf(uc_ref.b - uc_out.b) < 0.2f);
|
|
CHECK(fabsf(uc_ref.c - uc_out.c) < 0.4f);
|
|
CHECK(fabsf(uc_ref.alpha - uc_out.alpha) < 0.2f);
|
|
CHECK(fabsf(uc_ref.beta - uc_out.beta) < 0.2f);
|
|
CHECK(fabsf(uc_ref.gamma - uc_out.gamma) < 0.2f);
|
|
}
|
|
|
|
TEST_CASE("XtalOptimizer_refine_rotation_axis") {
|
|
// Geometry
|
|
DiffractionExperiment exp_i;
|
|
exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV)
|
|
.BeamX_pxl(1000)
|
|
.BeamY_pxl(1000)
|
|
.PoniRot1_rad(0.01)
|
|
.PoniRot2_rad(0.02)
|
|
.DetectorDistance_mm(200);
|
|
|
|
// Base lattice (non-pathological)
|
|
CrystalLattice latt_base(40, 50, 80, 90, 95, 90);
|
|
auto uc_ref = latt_base.GetUnitCell();
|
|
|
|
// Rotation axis: around X with 3 deg per image
|
|
GoniometerAxis axis("omega", 0.0f, 3.0f, Coord(1,0,0), std::nullopt);
|
|
|
|
BraggPredictionSettings prediction_settings{
|
|
.high_res_A = 1.5,
|
|
.ewald_dist_cutoff = 0.002
|
|
};
|
|
|
|
|
|
BraggPrediction prediction;
|
|
|
|
const size_t nimages = 10;
|
|
std::vector<std::vector<SpotToSave>> spots(nimages);
|
|
|
|
// Predict reflections for images at 0-30 deg.
|
|
for (int img = 0; img < nimages; ++img) {
|
|
// For a rotated image, per-image lattice is obtained as Multiply(rot.transpose())
|
|
const float angle_deg = axis.GetAngle_deg(img) + axis.GetWedge_deg() / 2.0f;
|
|
const RotMatrix rot = axis.GetTransformationAngle(angle_deg);
|
|
const CrystalLattice latt_img = latt_base.Multiply(rot.transpose());
|
|
|
|
const auto n = prediction.Calc(exp_i, latt_img, prediction_settings);
|
|
for (int i = 0; i < n; ++i) {
|
|
const auto& r = prediction.GetReflections().at(i);
|
|
SpotToSave s{};
|
|
s.x = r.predicted_x;
|
|
s.y = r.predicted_y;
|
|
s.image = img; // provide image index for rotation-aware refinement
|
|
s.intensity = 1.0f; // minimal positive value
|
|
s.phi = angle_deg;
|
|
s.ice_ring = false;
|
|
s.indexed = true;
|
|
spots.at(img).push_back(s);
|
|
}
|
|
}
|
|
|
|
// Seed slightly perturbed geometry and lattice; provide rotation axis for refinement
|
|
XtalOptimizerData xtal_opt{};
|
|
xtal_opt.latt = CrystalLattice(39.7f, 50.6f, 79.6f, 90.0f, 94.5f, 90.5f);
|
|
xtal_opt.geom.BeamX_pxl(1003).BeamY_pxl(997).DetectorDistance_mm(200.0)
|
|
.PoniRot1_rad(0.01).PoniRot2_rad(0.02);
|
|
xtal_opt.crystal_system = gemmi::CrystalSystem::Monoclinic;
|
|
xtal_opt.axis = GoniometerAxis("omega", 0.0f, 3.0f,
|
|
Coord(0.8, 0.05, 0.05).Normalize(),
|
|
std::nullopt);
|
|
xtal_opt.min_spots = 200;
|
|
xtal_opt.refine_beam_center = true;
|
|
xtal_opt.refine_detector_angles = false;
|
|
xtal_opt.refine_rotation_axis = true;
|
|
xtal_opt.max_time = 30.0;
|
|
|
|
auto t0 = std::chrono::high_resolution_clock::now();
|
|
REQUIRE(XtalOptimizer(xtal_opt, spots));
|
|
auto t1 = std::chrono::high_resolution_clock::now();
|
|
std::cout << "XtalOptimizer (rotation 4 images) took "
|
|
<< std::chrono::duration_cast<std::chrono::microseconds>(t1 - t0).count()
|
|
<< " microseconds" << std::endl;
|
|
|
|
const auto uc_out = xtal_opt.latt.GetUnitCell();
|
|
|
|
// Geometry checks
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamX_pxl() - exp_i.GetBeamX_pxl()) < 0.2f);
|
|
CHECK(fabsf(xtal_opt.geom.GetBeamY_pxl() - exp_i.GetBeamY_pxl()) < 0.2f);
|
|
|
|
// Lattice checks
|
|
CHECK(fabsf(uc_ref.a - uc_out.a) < 0.2f);
|
|
CHECK(fabsf(uc_ref.b - uc_out.b) < 0.2f);
|
|
CHECK(fabsf(uc_ref.c - uc_out.c) < 0.3f);
|
|
CHECK(fabsf(uc_ref.alpha - uc_out.alpha) < 0.2f);
|
|
CHECK(fabsf(uc_ref.beta - uc_out.beta) < 0.2f);
|
|
CHECK(fabsf(uc_ref.gamma - uc_out.gamma) < 0.2f);
|
|
CHECK(fabsf(xtal_opt.axis->GetAxis().x - 1.0) < 0.01f);
|
|
CHECK(fabsf(xtal_opt.axis->GetAxis().y) < 0.01f);
|
|
CHECK(fabsf(xtal_opt.axis->GetAxis().z) < 0.01f);
|
|
}
|
|
|
|
// --- helpers for lattice sanity tests ---
|
|
#include <Eigen/Dense>
|