leonarski_fandClaude Opus 5.5
e4eaeba874
tools/battery/still/agipd_to_int16.py converts an AGIPD 1M Cheetah CXI file into the only form rugnux
supports for AGIPD: each module's eight 64-row ASICs placed 66 rows apart (the 2-pixel ASIC gaps
inserted), one 526 x 128 panel per module; the per-frame mask and mask_shared (bad bits = the
geometry's mask_bad) and the geometry's bad regions folded into the data as -32768; values
round(value / ADU per photon) as int16 (clipped to +-32767; no pixel clipped in either set);
bitshuffle/LZ4, one chunk per frame at /entry_1/data_1/data. It writes the matching CrystFEL geometry
(16 panels, each the rigid fit of its ASICs, module coffset kept) and refuses a module whose ASICs are
not on one straight line at a whole-pixel pitch within 0.5 px (worst here: 0.16 px for CXIDB 83,
0.002 px for 152).
Reader: a stored value equal to the smallest value of its signed integer type (INT16_MIN, INT32_MIN)
is no measurement. The reader decodes to float and divides by the gain, so -32768 used to come out as a
count; the CrystFEL parser takes only module geometry, so a flag_lessthan in the .geom could not say it.
Still battery: euxfel83 / euxfel152 read conv/*_int16.h5 with conv/agipd_int16.geom and --gain 1.
Full sets, -N 16, --refine-geometry=off, --model (other jobs on the machine, load 7-23):
indexed d_min CC1/2 CC_HALF_SPLIT R_split R_FREE main pass wall
CXIDB 83 original 90.9% 1.84 0.641 0.617 0.482 0.280 89.3 s (31/s) 105 s
converted 94.1% 1.84 0.686 0.644 0.470 0.277 6.5 s (429/s) 20 s
CXIDB 152 original 17.5% 5.37 0.521 0.479 0.590 0.308 131.7 s (29/s) 153 s
converted 62.8% 2.92 0.876 0.891 0.354 0.284 13.2 s (287/s) 27 s
Files 4.7 -> 1.3 GB and 5.7 -> 1.1 GB. The 152 gain is the mask: the reader took only Cheetah's
generic bad bits, while the geometry's mask_bad (0xffff) also marks three broken modules (bit 32) and
its bad region a fourth; the converter run with the generic bits instead gives 17.5% / 3.41 A again.
With geometry refinement on (16 rigid modules): 83 94.2%, 1.75 A, CC1/2 0.665, R_FREE 0.284, 49 s;
152 65.1%, 2.94 A, CC1/2 0.843, R_FREE 0.278, 27 s - mixed, the manifest keeps it off.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu
189 lines
9.1 KiB
C++
189 lines
9.1 KiB
C++
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include <catch2/catch_all.hpp>
|
|
#include <cstdint>
|
|
#include <filesystem>
|
|
#include <hdf5.h>
|
|
#include <string>
|
|
#include <vector>
|
|
|
|
#include "../reader/JFJochXFELReader.h"
|
|
|
|
// Synthetic SACLA Cheetah and Cheetah CXI files, written with the HDF5 C API the way those pipelines
|
|
// write them: whole frames as single deflated (and byte-shuffled) chunks, so the direct-chunk path is
|
|
// what is read.
|
|
namespace {
|
|
|
|
void WriteFrames(hid_t loc, const std::string &name, hid_t type, const void *data, std::vector<hsize_t> dims) {
|
|
const hid_t space = H5Screate_simple(static_cast<int>(dims.size()), dims.data(), nullptr);
|
|
std::vector<hsize_t> chunk = dims;
|
|
if (chunk.size() == 3)
|
|
chunk[0] = 1;
|
|
const hid_t dcpl = H5Pcreate(H5P_DATASET_CREATE);
|
|
H5Pset_chunk(dcpl, static_cast<int>(chunk.size()), chunk.data());
|
|
H5Pset_shuffle(dcpl);
|
|
H5Pset_deflate(dcpl, 4);
|
|
const hid_t ds = H5Dcreate2(loc, name.c_str(), type, space, H5P_DEFAULT, dcpl, H5P_DEFAULT);
|
|
H5Dwrite(ds, type, H5S_ALL, H5S_ALL, H5P_DEFAULT, data);
|
|
H5Dclose(ds);
|
|
H5Pclose(dcpl);
|
|
H5Sclose(space);
|
|
}
|
|
|
|
void WriteVector(hid_t loc, const std::string &name, const std::vector<double> &v) {
|
|
const hsize_t n = v.size();
|
|
const hid_t space = H5Screate_simple(1, &n, nullptr);
|
|
const hid_t ds = H5Dcreate2(loc, name.c_str(), H5T_NATIVE_DOUBLE, space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
|
|
H5Dwrite(ds, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, v.data());
|
|
H5Dclose(ds);
|
|
H5Sclose(space);
|
|
}
|
|
|
|
void WriteScalar(hid_t loc, const std::string &name, double v) {
|
|
const hid_t space = H5Screate(H5S_SCALAR);
|
|
const hid_t ds = H5Dcreate2(loc, name.c_str(), H5T_NATIVE_DOUBLE, space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
|
|
H5Dwrite(ds, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &v);
|
|
H5Dclose(ds);
|
|
H5Sclose(space);
|
|
}
|
|
|
|
int32_t PixelAt(JFJochXFELReader &reader, int64_t image, size_t col, size_t row) {
|
|
auto raw = reader.GetRawImage(image);
|
|
REQUIRE(raw);
|
|
return reinterpret_cast<const int32_t *>(raw->image_buffer.data())[row * raw->image.GetWidth() + col];
|
|
}
|
|
|
|
} // namespace
|
|
|
|
TEST_CASE("JFJochXFELReader_SACLA", "[portable]") {
|
|
// Two sensors of 8 rows x 4 columns stacked in a 16 x 4 frame, 50 um pixels. The first lies square
|
|
// with its first corner on the beam; the second is turned by 90 degrees with its corner 20 pixels
|
|
// along x (posx = -1000 um), so its fast axis runs along y and its slow axis along -x.
|
|
const std::string name = "xfel_reader_sacla.h5";
|
|
{
|
|
const hid_t f = H5Fcreate(name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
|
|
const hid_t meta = H5Gcreate2(f, "metadata", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
|
|
WriteVector(meta, "posx_in_um", {0, -1000});
|
|
WriteVector(meta, "posy_in_um", {0, 0});
|
|
WriteVector(meta, "posz_in_um", {0, 0});
|
|
WriteVector(meta, "angle_in_rad", {0, 90}); // degrees, whatever the name says
|
|
WriteVector(meta, "pixelsizex_in_um", {50, 50});
|
|
WriteVector(meta, "pixelsizey_in_um", {50, 50});
|
|
WriteScalar(meta, "distance_in_mm", 60.0);
|
|
std::vector<uint16_t> pixelmask(16 * 4, 0);
|
|
pixelmask[1] = 1; // first sensor, row 0, column 1
|
|
WriteFrames(meta, "pixelmask", H5T_NATIVE_UINT16, pixelmask.data(), {16, 4});
|
|
H5Gclose(meta);
|
|
for (int shot = 0; shot < 2; shot++) {
|
|
const std::string group = "tag-" + std::to_string(100 + shot);
|
|
const hid_t g = H5Gcreate2(f, group.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
|
|
std::vector<int16_t> frame(16 * 4, 10);
|
|
frame[0] = static_cast<int16_t>(20 + shot); // first sensor (fs 0, ss 0)
|
|
frame[8 * 4] = -5; // second sensor (fs 0, ss 0): below zero
|
|
frame[8 * 4 + 1] = 33; // second sensor (fs 1, ss 0)
|
|
frame[9 * 4] = 44; // second sensor (fs 0, ss 1)
|
|
frame[2 * 4] = INT16_MIN; // first sensor (fs 0, ss 2): no measurement
|
|
WriteFrames(g, "data", H5T_NATIVE_INT16, frame.data(), {16, 4});
|
|
WriteScalar(g, "photon_energy_ev", 10000.0 + 10.0 * shot);
|
|
H5Gclose(g);
|
|
}
|
|
H5Fclose(f);
|
|
}
|
|
|
|
REQUIRE(JFJochXFELReader::CanReadPanelTable(name));
|
|
CHECK_FALSE(JFJochXFELReader::CanReadCXI(name));
|
|
JFJochXFELReader reader;
|
|
REQUIRE_NOTHROW(reader.ReadPanelTable({name}));
|
|
CHECK(reader.GetNumberOfImages() == 2);
|
|
CHECK_FALSE(reader.CountsInPhotons());
|
|
const auto x = reader.GetDataset()->experiment;
|
|
CHECK(x.GetIncidentEnergy_keV() == Catch::Approx(10.005));
|
|
CHECK(x.GetDetectorDistance_mm() == Catch::Approx(60.0));
|
|
CHECK(x.GetPixelSize_mm() == Catch::Approx(0.05));
|
|
// The first sensor's first pixel centre is at (0.5, 0.5) from the beam, which the flat image puts on
|
|
// a whole column: the beam at 1.5, the pixel at column and row 2.
|
|
CHECK(x.GetBeamX_pxl() == Catch::Approx(1.5));
|
|
CHECK(x.GetBeamY_pxl() == Catch::Approx(1.5));
|
|
CHECK(PixelAt(reader, 0, 2, 2) == 20);
|
|
CHECK(PixelAt(reader, 1, 2, 2) == 21);
|
|
// Second sensor: (fs 0, ss 0) centre at x = 20 - 0.5, y = 0.5; fs runs along +y, ss along -x.
|
|
CHECK(PixelAt(reader, 0, 21, 2) == -5);
|
|
CHECK(PixelAt(reader, 0, 21, 3) == 33);
|
|
CHECK(PixelAt(reader, 0, 20, 2) == 44);
|
|
// The smallest int16 is no measurement, not a count.
|
|
CHECK(PixelAt(reader, 0, 2, 4) == INT32_MIN);
|
|
// The file's mask holds the first sensor's (fs 1, ss 0), at column 3, row 2.
|
|
const auto &mask = reader.GetDataset()->pixel_mask->GetMask();
|
|
CHECK(mask[2 * x.GetXPixelsNum() + 3] != 0);
|
|
CHECK(mask[2 * x.GetXPixelsNum() + 2] == 0);
|
|
|
|
reader.ApplyGain(2.0f);
|
|
CHECK(reader.CountsInPhotons());
|
|
CHECK(PixelAt(reader, 0, 2, 2) == 10);
|
|
reader.Close();
|
|
std::filesystem::remove(name);
|
|
}
|
|
|
|
TEST_CASE("JFJochXFELReader_CXI", "[portable]") {
|
|
// A single-panel Cheetah CXI file: three 6 x 6 frames of signed 32-bit counts, 75 um pixels,
|
|
// 9 keV stated in joules, a mask_shared with one bad pixel.
|
|
const std::string name = "xfel_reader_test.cxi";
|
|
{
|
|
const hid_t f = H5Fcreate(name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
|
|
for (const char *g: {"entry_1", "entry_1/data_1", "entry_1/instrument_1", "entry_1/instrument_1/detector_1",
|
|
"entry_1/instrument_1/source_1"})
|
|
H5Gclose(H5Gcreate2(f, g, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT));
|
|
std::vector<int32_t> frames(3 * 36, 5);
|
|
frames[0] = -7; // frame 0, pixel (0, 0): dark-subtracted below zero
|
|
frames[36 + 7] = 70000; // frame 1, pixel (1, 1)
|
|
WriteFrames(f, "entry_1/data_1/data", H5T_NATIVE_INT32, frames.data(), {3, 6, 6});
|
|
const std::string det = "entry_1/instrument_1/detector_1/";
|
|
WriteVector(f, det + "x_pixel_size", {75e-6, 75e-6, 75e-6});
|
|
WriteVector(f, det + "distance", {0.1, 0.1, 0.1});
|
|
WriteVector(f, "entry_1/instrument_1/source_1/energy", std::vector<double>(3, 9000.0 * 1.602176634e-19));
|
|
std::vector<uint16_t> mask(36, 0);
|
|
mask[35] = 8; // dead
|
|
WriteFrames(f, det + "mask_shared", H5T_NATIVE_UINT16, mask.data(), {6, 6});
|
|
H5Fclose(f);
|
|
}
|
|
|
|
REQUIRE(JFJochXFELReader::CanReadCXI(name));
|
|
CHECK_FALSE(JFJochXFELReader::CanReadPanelTable(name));
|
|
CHECK_FALSE(JFJochXFELReader::CXIRefusal(name).has_value());
|
|
JFJochXFELReader reader;
|
|
REQUIRE_NOTHROW(reader.ReadCXI({name}));
|
|
CHECK(reader.GetNumberOfImages() == 3);
|
|
const auto x = reader.GetDataset()->experiment;
|
|
CHECK(x.GetIncidentEnergy_keV() == Catch::Approx(9.0));
|
|
CHECK(x.GetPixelSize_mm() == Catch::Approx(0.075));
|
|
CHECK(x.GetDetectorDistance_mm() == Catch::Approx(100.0));
|
|
// Corner at (-3, -3): the pixel (fs, ss) lands at column and row 2 + fs, 2 + ss.
|
|
CHECK(PixelAt(reader, 0, 2, 2) == -7);
|
|
CHECK(PixelAt(reader, 1, 3, 3) == 70000);
|
|
CHECK(PixelAt(reader, 2, 4, 4) == 5);
|
|
const auto &m = reader.GetDataset()->pixel_mask->GetMask();
|
|
CHECK(m[7 * x.GetXPixelsNum() + 7] != 0);
|
|
CHECK(m[6 * x.GetXPixelsNum() + 6] == 0);
|
|
reader.Close();
|
|
std::filesystem::remove(name);
|
|
}
|
|
|
|
TEST_CASE("JFJochXFELReader_CXI_TiledRefused", "[portable]") {
|
|
// A Cheetah CXI slab of a tiled detector is not square, and is not read as one flat panel.
|
|
const std::string name = "xfel_reader_tiled.cxi";
|
|
{
|
|
const hid_t f = H5Fcreate(name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
|
|
for (const char *g: {"entry_1", "entry_1/data_1"})
|
|
H5Gclose(H5Gcreate2(f, g, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT));
|
|
std::vector<int16_t> frames(2 * 8 * 4, 1);
|
|
WriteFrames(f, "entry_1/data_1/data", H5T_NATIVE_INT16, frames.data(), {2, 8, 4});
|
|
H5Fclose(f);
|
|
}
|
|
CHECK_FALSE(JFJochXFELReader::CanReadCXI(name));
|
|
const auto refusal = JFJochXFELReader::CXIRefusal(name);
|
|
REQUIRE(refusal.has_value());
|
|
CHECK(refusal->find("tiled detector") != std::string::npos);
|
|
std::filesystem::remove(name);
|
|
}
|