From ce528f3a14981d0de1109de03c9e6f8df90886cc Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 20 Oct 2026 13:03:12 +0200 Subject: [PATCH] rugnux: read stills through a CrystFEL geometry file (--crystfel-geometry) A CrystFEL .geom (current and old syntax) is parsed for panels (pixel ranges, fs/ss, corners), clen/coffset (a number, or an HDF5 path read as mm), res, adu_per_eV/adu_per_photon, max_adu, no_index, bad regions (panel, whole-frame or lab-plane) and the HDF5 mask (mask/mask_file, mask_good/mask_bad). JFJochCrystFELReader reads the frames the geometry's data path names - a 3D stack (Cheetah CXI), one 2D frame per file, or one group per shot (/%/data) - from a file or a directory, puts the panels on one flat image at their rounded positions with gaps and bad regions masked, and converts to photons where a gain is stated (otherwise the pre-scan measures it). One wavelength and one camera length per run: values taken from the files are averaged over the first file. Panel positions are used as given. Frames stored as one deflate(+shuffle) chunk are read raw under the HDF5 lock and decoded outside it. Replaces the MPCCD-specific reader of the previous commit. Co-Authored-By: Claude Opus 5.5 --- reader/CMakeLists.txt | 6 +- reader/CrystFELGeometry.cpp | 138 ++++++++++++ reader/CrystFELGeometry.h | 57 +++++ reader/JFJochCrystFELReader.cpp | 358 ++++++++++++++++++++++++++++++++ reader/JFJochCrystFELReader.h | 81 ++++++++ reader/JFJochMPCCDReader.cpp | 350 ------------------------------- reader/JFJochMPCCDReader.h | 78 ------- rugnux/rugnux_cli.cpp | 31 ++- 8 files changed, 653 insertions(+), 446 deletions(-) create mode 100644 reader/CrystFELGeometry.cpp create mode 100644 reader/CrystFELGeometry.h create mode 100644 reader/JFJochCrystFELReader.cpp create mode 100644 reader/JFJochCrystFELReader.h delete mode 100644 reader/JFJochMPCCDReader.cpp delete mode 100644 reader/JFJochMPCCDReader.h diff --git a/reader/CMakeLists.txt b/reader/CMakeLists.txt index 16f137359..fe0f0d7ec 100644 --- a/reader/CMakeLists.txt +++ b/reader/CMakeLists.txt @@ -12,8 +12,10 @@ ADD_LIBRARY(JFJochReader STATIC MarCCD.h JFJochMarCCDReader.cpp JFJochMarCCDReader.h - JFJochMPCCDReader.cpp - JFJochMPCCDReader.h + CrystFELGeometry.cpp + CrystFELGeometry.h + JFJochCrystFELReader.cpp + JFJochCrystFELReader.h SMV.cpp SMV.h JFJochSMVReader.cpp diff --git a/reader/CrystFELGeometry.cpp b/reader/CrystFELGeometry.cpp new file mode 100644 index 000000000..f377fbd17 --- /dev/null +++ b/reader/CrystFELGeometry.cpp @@ -0,0 +1,138 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "CrystFELGeometry.h" + +#include +#include +#include + +#include "../common/JFJochException.h" + +namespace crystfel { + +namespace { + +std::string Trim(const std::string &s) { + const auto b = s.find_first_not_of(" \t\r"); + if (b == std::string::npos) + return ""; + return s.substr(b, s.find_last_not_of(" \t\r") - b + 1); +} + +// "+0.000390x -1.000000y", "-x", "0.5x+0.5y": the in-plane direction of a panel axis. +void ParseDirection(const std::string &txt, double &x, double &y) { + x = 0; + y = 0; + static const std::regex term(R"(([+-]?)\s*([0-9.]*(?:[eE][+-]?[0-9]+)?)\s*([xyz]))"); + for (auto it = std::sregex_iterator(txt.begin(), txt.end(), term); it != std::sregex_iterator(); ++it) { + double v = (*it)[2].str().empty() ? 1.0 : std::stod((*it)[2].str()); + if ((*it)[1] == "-") + v = -v; + if ((*it)[3] == "x") + x = v; + else if ((*it)[3] == "y") + y = v; + } +} + +// Sets one panel property; false if the key is not one. +bool SetPanel(Panel &p, const std::string &key, const std::string &value) { + if (key == "min_fs") p.min_fs = std::stoll(value); + else if (key == "max_fs") p.max_fs = std::stoll(value); + else if (key == "min_ss") p.min_ss = std::stoll(value); + else if (key == "max_ss") p.max_ss = std::stoll(value); + else if (key == "fs") ParseDirection(value, p.fs_x, p.fs_y); + else if (key == "ss") ParseDirection(value, p.ss_x, p.ss_y); + else if (key == "corner_x") p.corner_x = std::stod(value); + else if (key == "corner_y") p.corner_y = std::stod(value); + else if (key == "res") p.res = std::stod(value); + else if (key == "clen") p.clen = value; + else if (key == "coffset") p.coffset = std::stod(value); + else if (key == "adu_per_eV") p.adu_per_eV = std::stod(value); + else if (key == "adu_per_photon") p.adu_per_photon = std::stod(value); + else if (key == "max_adu") p.max_adu = std::stod(value); + else if (key == "no_index") p.no_index = std::stoi(value) != 0; + else + return false; + return true; +} + +} // namespace + +Geometry Read(const std::string &filename) { + std::ifstream in(filename); + if (!in) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Cannot open " + filename); + + Geometry geom; + // Global panel keys hold for every panel that appears after them, as in the old syntax; a file in + // the current syntax states them first, so this reads both the same. + Panel defaults; + std::map panel_index; + std::map bad_index; + + std::string line; + while (std::getline(in, line)) { + line = line.substr(0, line.find(';')); + const auto eq = line.find('='); + if (eq == std::string::npos) + continue; + const std::string key = Trim(line.substr(0, eq)); + const std::string value = Trim(line.substr(eq + 1)); + const auto slash = key.find('/'); + + if (slash == std::string::npos) { + if (key == "data") geom.data = value; + else if (key == "photon_energy") geom.photon_energy = value; + else if (key == "mask") geom.mask = value; + else if (key == "mask_file") geom.mask_file = value; + else if (key == "mask_good") geom.mask_good = std::stoul(value, nullptr, 0); + else if (key == "mask_bad") geom.mask_bad = std::stoul(value, nullptr, 0); + else SetPanel(defaults, key, value); // anything else (rigid groups, dim, ...) is not needed + continue; + } + + const std::string name = key.substr(0, slash); + const std::string prop = key.substr(slash + 1); + + if (name.starts_with("bad")) { + if (!bad_index.contains(name)) { + bad_index[name] = geom.bad.size(); + geom.bad.push_back({}); + } + auto &b = geom.bad[bad_index[name]]; + if (prop == "panel") b.panel = value; + else if (prop == "min_x") { b.in_lab = true; b.min_x = std::stod(value); } + else if (prop == "max_x") { b.in_lab = true; b.max_x = std::stod(value); } + else if (prop == "min_y") { b.in_lab = true; b.min_y = std::stod(value); } + else if (prop == "max_y") { b.in_lab = true; b.max_y = std::stod(value); } + else if (prop == "min_fs") b.min_fs = std::stoll(value); + else if (prop == "max_fs") b.max_fs = std::stoll(value); + else if (prop == "min_ss") b.min_ss = std::stoll(value); + else if (prop == "max_ss") b.max_ss = std::stoll(value); + continue; + } + + if (!panel_index.contains(name)) { + panel_index[name] = geom.panels.size(); + geom.panels.push_back(defaults); + geom.panels.back().name = name; + } + SetPanel(geom.panels[panel_index[name]], prop, value); + } + + if (geom.panels.empty()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, filename + " defines no panels"); + for (const auto &p: geom.panels) { + if (!(p.res > 0)) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + filename + ": panel " + p.name + " has no pixel size (res)"); + if (p.clen.empty()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + filename + ": panel " + p.name + " has no camera length (clen)"); + } + return geom; +} + +} // namespace crystfel diff --git a/reader/CrystFELGeometry.h b/reader/CrystFELGeometry.h new file mode 100644 index 000000000..9253f9825 --- /dev/null +++ b/reader/CrystFELGeometry.h @@ -0,0 +1,57 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include +#include + +// A CrystFEL detector geometry file (.geom), as far as it places pixels and says what they are worth: +// the panels with their pixel ranges in the stored frame, fast/slow-scan directions and corners, the +// camera length, the pixel size, the gain and the regions to leave out. Both the current syntax and +// the older one - global keys that only take effect for the panels that follow them, bad regions +// stated in whole-frame coordinates - are read. +// +// Every length is in pixels except clen and coffset, which are in metres. A value CrystFEL takes from +// the image file instead (clen = /path, photon_energy = /path) is kept as that path, for the reader +// to resolve. +namespace crystfel { + +struct Panel { + std::string name; + int64_t min_fs = 0, max_fs = 0, min_ss = 0, max_ss = 0; // inclusive, in the stored frame + double fs_x = 1, fs_y = 0, ss_x = 0, ss_y = 1; + double corner_x = 0, corner_y = 0; + double res = 0; // pixels per metre + std::string clen; // metres, or an HDF5 path + double coffset = 0; // metres + std::optional adu_per_eV; + std::optional adu_per_photon; + std::optional max_adu; + bool no_index = false; +}; + +// A region to leave out: either a pixel range on one panel (or, in the old syntax, on the whole +// frame), or a box in the laboratory plane, in pixels. +struct BadRegion { + std::string panel; // empty: fs/ss are whole-frame coordinates + bool in_lab = false; + double min_x = 0, max_x = 0, min_y = 0, max_y = 0; + int64_t min_fs = 0, max_fs = 0, min_ss = 0, max_ss = 0; +}; + +struct Geometry { + std::vector panels; + std::vector bad; + std::string data = "/data/data"; + std::string photon_energy; // eV, or an HDF5 path + std::string mask; // HDF5 path, or empty + std::string mask_file; // file holding it, or empty for the image file itself + uint32_t mask_good = 0, mask_bad = 0; +}; + +Geometry Read(const std::string &filename); + +} // namespace crystfel diff --git a/reader/JFJochCrystFELReader.cpp b/reader/JFJochCrystFELReader.cpp new file mode 100644 index 000000000..77fb47f26 --- /dev/null +++ b/reader/JFJochCrystFELReader.cpp @@ -0,0 +1,358 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "JFJochCrystFELReader.h" + +#include +#include +#include +#include + +#include + +#include "../common/JFJochException.h" +#include "../common/Logger.h" +#include "../writer/HDF5Objects.h" + +namespace { + +bool IsPath(const std::string &s) { + return !s.empty() && s[0] == '/'; +} + +std::string ReplaceEvent(std::string path, const std::string &group) { + const auto pos = path.find('%'); + if (pos != std::string::npos) + path.replace(pos, 1, group); + return path; +} + +} // namespace + +JFJochCrystFELReader::JFJochCrystFELReader() = default; + +JFJochCrystFELReader::~JFJochCrystFELReader() { + Close(); +} + +// One number for the whole run from a value the geometry takes from the image files: a scalar, a +// per-frame vector, or one scalar per shot group. +double JFJochCrystFELReader::MeanOverFirstFile(const std::string &path) { + auto &file = *open_files_[0]; + std::vector values; + if (path.find('%') != std::string::npos) { + for (const auto &f: frames_) + if (f.file == 0) + values.push_back(HDF5DataSet(file, ReplaceEvent(path, f.group)).ReadScalar()); + } else if (file.GetDimension(path).empty()) + values.push_back(HDF5DataSet(file, path).ReadScalar()); + else + values = file.ReadVector(path); + if (values.empty()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "No values in " + path); + return std::accumulate(values.begin(), values.end(), 0.0) / static_cast(values.size()); +} + +void JFJochCrystFELReader::ReadFiles(const std::string &path, const std::string &geometry_file) { + Close(); + geom_ = crystfel::Read(geometry_file); + + if (std::filesystem::is_directory(path)) { + for (const auto &e: std::filesystem::directory_iterator(path)) + if (e.path().extension() == ".h5" || e.path().extension() == ".cxi") + files_.push_back(e.path().string()); + std::sort(files_.begin(), files_.end()); + } else + files_.push_back(path); + if (files_.empty()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "No .h5 or .cxi files in " + path); + + std::unique_lock ul(hdf5_mutex); + + // Every frame of every file, in file order. + const auto pct = geom_.data.find('%'); + for (size_t f = 0; f < files_.size(); f++) { + open_files_.emplace_back(std::make_unique(files_[f])); + auto &file = *open_files_.back(); + if (pct != std::string::npos) { + const std::string parent = geom_.data.substr(0, pct); + for (const auto &group: file.FindLeafs(parent.empty() ? "/" : parent)) { + const auto dataset = ReplaceEvent(geom_.data, group); + if (file.Exists(dataset)) + frames_.push_back({f, dataset, -1, group}); + } + } else { + const auto dim = file.GetDimension(geom_.data); + if (dim.size() == 3) { + for (hsize_t i = 0; i < dim[0]; i++) + frames_.push_back({f, geom_.data, static_cast(i), ""}); + } else if (dim.size() == 2) + frames_.push_back({f, geom_.data, -1, ""}); + } + } + if (frames_.empty()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "No frames at " + geom_.data + " in " + path); + + const auto dim = open_files_[frames_[0].file]->GetDimension(frames_[0].dataset); + frame_height_ = dim[dim.size() - 2]; + frame_width_ = dim[dim.size() - 1]; + + const auto &p0 = geom_.panels[0]; + const double clen_m = IsPath(p0.clen) ? MeanOverFirstFile(p0.clen) * 1e-3 : std::stod(p0.clen); + const double distance_mm = (clen_m + p0.coffset) * 1000.0; + if (geom_.photon_energy.empty()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + geometry_file + " states no photon_energy"); + const double photon_energy_eV = IsPath(geom_.photon_energy) ? MeanOverFirstFile(geom_.photon_energy) + : std::stod(geom_.photon_energy); + for (const auto &p: geom_.panels) + if (p.res != p0.res || p.clen != p0.clen || p.coffset != p0.coffset) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + geometry_file + ": panels at different distances or with different " + "pixel sizes cannot be put on one flat detector"); + + std::vector file_mask; + if (!geom_.mask.empty()) { + if (!geom_.mask_file.empty()) { + HDF5ReadOnlyFile mask_file(geom_.mask_file); + file_mask = mask_file.ReadVector(geom_.mask, {0, 0}, {frame_height_, frame_width_}); + } else + file_mask = open_files_[0]->ReadVector(geom_.mask, {0, 0}, {frame_height_, frame_width_}); + } + CheckDirectChunks(); + ul.unlock(); + + // CrystFEL's laboratory frame and the flat image's (column, row, beam) are both right-handed with + // the beam along z, so the point (X, Y) lands at column beam_x + X, row beam_y + Y - and CrystFEL's + // horizontal x stays the image's horizontal, where the polarisation correction expects it. A pixel + // centre is at its index, so beam_x is chosen to put the first panel's pixel centres on whole columns. + auto position = [](const crystfel::Panel &p, double fs, double ss) { + return std::pair{p.corner_x + fs * p.fs_x + ss * p.ss_x, p.corner_y + fs * p.fs_y + ss * p.ss_y}; + }; + double max_x = -1e9, min_x = 1e9, max_y = -1e9, min_y = 1e9; + for (const auto &p: geom_.panels) + for (double fs: {0.0, double(p.max_fs - p.min_fs + 1)}) + for (double ss: {0.0, double(p.max_ss - p.min_ss + 1)}) { + const auto [x, y] = position(p, fs, ss); + max_x = std::max(max_x, x); + min_x = std::min(min_x, x); + max_y = std::max(max_y, y); + min_y = std::min(min_y, y); + } + const auto [x0, y0] = position(p0, 0.5, 0.5); + const double beam_x = std::ceil(x0 - min_x) + 1 - x0; + const double beam_y = std::ceil(y0 - min_y) + 1 - y0; + width_ = static_cast(std::ceil(beam_x + max_x)) + 2; + height_ = static_cast(std::ceil(beam_y + max_y)) + 2; + + // Everything starts as a gap; the panels' pixels are put on top, less the bad regions. + std::vector mask(width_ * height_, 1u << PixelMask::ModuleGapPixelBit); + destination_.assign(frame_width_ * frame_height_, -1); + adu_per_photon_.assign(frame_width_ * frame_height_, 1.0f); + counts_in_photons_ = true; + std::optional saturation; + + for (const auto &p: geom_.panels) { + if (p.max_fs >= static_cast(frame_width_) || p.max_ss >= static_cast(frame_height_)) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Panel " + p.name + " lies outside the stored frame"); + double gain = 1.0; + if (p.adu_per_photon.has_value()) + gain = p.adu_per_photon.value(); + else if (p.adu_per_eV.has_value()) + gain = p.adu_per_eV.value() * photon_energy_eV; + else + counts_in_photons_ = false; + if (p.max_adu.has_value()) + saturation = std::min(saturation.value_or(INT64_MAX), p.max_adu.value() / gain); + if (p.no_index) + continue; + + for (int64_t ss = p.min_ss; ss <= p.max_ss; ss++) { + for (int64_t fs = p.min_fs; fs <= p.max_fs; fs++) { + const auto [x, y] = position(p, fs - p.min_fs + 0.5, ss - p.min_ss + 0.5); + bool bad = false; + for (const auto &b: geom_.bad) { + if (b.in_lab) + bad |= x >= b.min_x && x <= b.max_x && y >= b.min_y && y <= b.max_y; + else if (b.panel.empty() || b.panel == p.name) + bad |= fs >= b.min_fs && fs <= b.max_fs && ss >= b.min_ss && ss <= b.max_ss; + } + if (bad) + continue; + + const size_t src = ss * frame_width_ + fs; + const auto col = std::lround(beam_x + x); + const auto row = std::lround(beam_y + y); + const size_t dst = static_cast(row) * width_ + static_cast(col); + destination_[src] = static_cast(dst); + adu_per_photon_[src] = static_cast(gain); + + const uint32_t m = src < file_mask.size() ? file_mask[src] : 0; + const bool masked = (geom_.mask_bad != 0 && (m & geom_.mask_bad) != 0) + || (geom_.mask_good != 0 && (m & geom_.mask_good) != geom_.mask_good); + mask[dst] = masked ? (1u << PixelMask::ErrorPixelBit) : 0; + } + } + } + + dataset_ = std::make_shared(); + dataset_->experiment = default_experiment; + DetectorSetup detector = DetDECTRIS(static_cast(width_), static_cast(height_), + "CrystFEL geometry", {}); + detector.PixelSize_um(static_cast(1e6 / p0.res)); + // Nothing a parallax or sensor-efficiency correction needs is stated in a geometry file. + detector.SensorThickness_um(0); + if (saturation.has_value()) + detector.SaturationLimit(static_cast(std::floor(saturation.value()))); + detector.BitDepthImage(32); + detector.MinFrameTime(std::chrono::microseconds(0)); + detector.MinCountTime(std::chrono::microseconds(0)); + detector.ReadOutTime(std::chrono::nanoseconds(0)); + dataset_->experiment.Detector(detector); + dataset_->experiment.BeamX_pxl(static_cast(beam_x)); + dataset_->experiment.BeamY_pxl(static_cast(beam_y)); + dataset_->experiment.DetectorDistance_mm(static_cast(distance_mm)); + dataset_->experiment.IncidentEnergy_keV(static_cast(photon_energy_eV / 1000.0)); + dataset_->experiment.ImagesPerTrigger(static_cast(frames_.size())); + dataset_->pixel_mask = std::make_shared(mask); + + Logger("CrystFELReader").Info("{} frames in {} files, {} panels; {} x {} image, beam at ({:.2f}, {:.2f}) px, " + "{:.3f} mm, {:.1f} eV, {}{}", frames_.size(), files_.size(), geom_.panels.size(), + width_, height_, beam_x, beam_y, distance_mm, photon_energy_eV, + counts_in_photons_ ? "converted to photons" : "in detector units (no gain stated)", + direct_chunks_ ? ", chunks decoded in parallel" : ""); + SetStartMessage(dataset_); +} + +void JFJochCrystFELReader::CheckDirectChunks() { + const auto &frame = frames_[0]; + HDF5DataSet dataset(*open_files_[frame.file], frame.dataset); + + const hid_t dcpl = H5Dget_create_plist(dataset.GetID()); + const int rank = frame.index >= 0 ? 3 : 2; + std::vector chunk(rank, 0); + bool ok = H5Pget_layout(dcpl) == H5D_CHUNKED && H5Pget_chunk(dcpl, rank, chunk.data()) == rank + && chunk[rank - 1] == frame_width_ && chunk[rank - 2] == frame_height_ && (rank == 2 || chunk[0] == 1); + bool deflate = false; + shuffled_ = false; + for (int i = 0; ok && i < H5Pget_nfilters(dcpl); i++) { + unsigned int flags; + size_t nelements = 0; + const H5Z_filter_t filter = H5Pget_filter2(dcpl, i, &flags, &nelements, nullptr, 0, nullptr, nullptr); + if (filter == H5Z_FILTER_DEFLATE) + deflate = true; + else if (filter == H5Z_FILTER_SHUFFLE) + shuffled_ = true; + else + ok = false; + } + H5Pclose(dcpl); + + const hid_t type = H5Dget_type(dataset.GetID()); + stored_size_ = H5Tget_size(type); + stored_float_ = H5Tget_class(type) == H5T_FLOAT; + stored_signed_ = H5Tget_sign(type) == H5T_SGN_2; + ok &= (H5Tget_class(type) == H5T_INTEGER && (stored_size_ == 2 || stored_size_ == 4)) + || (stored_float_ && stored_size_ == 4); + H5Tclose(type); + + direct_chunks_ = ok && deflate; +} + +bool JFJochCrystFELReader::CountsInPhotons() const { + return counts_in_photons_; +} + +uint64_t JFJochCrystFELReader::GetNumberOfImages() const { + return frames_.size(); +} + +void JFJochCrystFELReader::Close() { + std::unique_lock ul(hdf5_mutex); + open_files_.clear(); + files_.clear(); + frames_.clear(); + dataset_.reset(); +} + +template +void JFJochCrystFELReader::DecodeInto(int64_t image_number, Buffer &buffer) const { + const auto &frame = frames_.at(image_number); + std::vector stored(frame_width_ * frame_height_); + if (direct_chunks_) { + std::vector compressed; + { + std::unique_lock ul(hdf5_mutex); + HDF5DataSet dataset(*open_files_[frame.file], frame.dataset); + if (frame.index >= 0) + dataset.ReadDirectChunk(compressed, {static_cast(frame.index), 0, 0}); + else + dataset.ReadDirectChunk(compressed, {0, 0}); + } + const size_t n = stored.size(); + std::vector raw(n * stored_size_); + uLongf raw_size = raw.size(); + if (uncompress(raw.data(), &raw_size, compressed.data(), compressed.size()) != Z_OK || raw_size != raw.size()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Cannot inflate " + frame.dataset + " of " + files_[frame.file]); + if (shuffled_) { + std::vector tmp(raw.size()); + for (size_t i = 0; i < n; i++) + for (size_t b = 0; b < stored_size_; b++) + tmp[i * stored_size_ + b] = raw[b * n + i]; + raw.swap(tmp); + } + for (size_t i = 0; i < n; i++) { + const uint8_t *v = raw.data() + i * stored_size_; + if (stored_float_) + stored[i] = *reinterpret_cast(v); + else if (stored_size_ == 2) + stored[i] = stored_signed_ ? *reinterpret_cast(v) : *reinterpret_cast(v); + else + stored[i] = stored_signed_ ? *reinterpret_cast(v) : *reinterpret_cast(v); + } + } else { + std::unique_lock ul(hdf5_mutex); + HDF5DataSet dataset(*open_files_[frame.file], frame.dataset); + if (frame.index >= 0) + dataset.ReadVector(stored, {static_cast(frame.index), 0, 0}, {1, frame_height_, frame_width_}); + else + dataset.ReadVector(stored, {0, 0}, {frame_height_, frame_width_}); + } + + buffer.assign(width_ * height_ * sizeof(int32_t), 0); + auto out = reinterpret_cast(buffer.data()); + for (size_t i = 0; i < stored.size(); i++) + if (destination_[i] >= 0) + out[destination_[i]] = static_cast(std::lround(stored[i] / adu_per_photon_[i])); +} + +bool JFJochCrystFELReader::LoadImage_i(std::shared_ptr &dataset, + DataMessage &message, + std::vector &buffer, + int64_t image_number, + bool update_dataset) { + (void) update_dataset; + if (!dataset || image_number < 0 || static_cast(image_number) >= frames_.size()) + return false; + DecodeInto(image_number, buffer); + message.image = CompressedImage(buffer.data(), buffer.size(), width_, height_, + CompressedImageMode::Int32, CompressionAlgorithm::NO_COMPRESSION); + message.number = image_number; + return true; +} + +bool JFJochCrystFELReader::ReadRawImage(int64_t image_number, JFJochReaderRawImage &image) { + if (image_number < 0 || static_cast(image_number) >= frames_.size()) + return false; + DecodeInto(image_number, image.image_buffer); + image.image = CompressedImage(image.image_buffer.data(), image.image_buffer.size(), width_, height_, + CompressedImageMode::Int32, CompressionAlgorithm::NO_COMPRESSION); + return true; +} + +std::vector JFJochCrystFELReader::ReadSpots(int64_t) const { + return {}; +} diff --git a/reader/JFJochCrystFELReader.h b/reader/JFJochCrystFELReader.h new file mode 100644 index 000000000..4eb111382 --- /dev/null +++ b/reader/JFJochCrystFELReader.h @@ -0,0 +1,81 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include + +#include "JFJochReader.h" +#include "CrystFELGeometry.h" + +class HDF5ReadOnlyFile; + +// Reads serial (stills) frames from HDF5 files laid out the way CrystFEL reads them, placed on the +// detector by a CrystFEL geometry file: Cheetah CXI files with the frames stacked along the first +// dimension, files holding one frame each, or files with one group per shot (data = /%/data). The +// panels are put on one flat image at their geometry-file positions, rounded to whole pixels, with +// everything between them and every bad region masked. Pixel values are converted to photons with the +// geometry's gain where it states one, and are otherwise left in detector units. +// +// One wavelength and one camera length hold for the whole run: where the geometry takes either from +// the image files, it is the mean over the first file's frames. +// +// A path is one file, or a directory, every .h5/.cxi file of which is read in name order. +class JFJochCrystFELReader : public JFJochReader { + struct Frame { + size_t file; + std::string dataset; + int64_t index; // along the first dimension of a 3D dataset; -1 for a 2D one + std::string group; // the shot group the geometry's % stands for, if any + }; + + crystfel::Geometry geom_; + std::vector files_; + std::vector> open_files_; + std::vector frames_; + std::shared_ptr dataset_; + + size_t frame_width_ = 0, frame_height_ = 0; // the stored frame + size_t width_ = 0, height_ = 0; // the flat image + // For every stored pixel, where it lands on the flat image (-1: nowhere) and what one photon reads. + std::vector destination_; + std::vector adu_per_photon_; + bool counts_in_photons_ = false; + + // Where every frame is stored as one deflated (and optionally byte-shuffled) chunk of a plain + // number type - the Cheetah outputs - the chunk is read raw under the HDF5 lock and decoded + // outside it, so decoding runs on every reading thread instead of one at a time. + bool direct_chunks_ = false; + bool shuffled_ = false; + bool stored_float_ = false, stored_signed_ = false; + size_t stored_size_ = 0; + + void CheckDirectChunks(); + + double MeanOverFirstFile(const std::string &path); + + bool LoadImage_i(std::shared_ptr &dataset, + DataMessage &message, + std::vector &buffer, + int64_t image_number, + bool update_dataset) override; + + template + void DecodeInto(int64_t image_number, Buffer &buffer) const; + +public: + JFJochCrystFELReader(); + ~JFJochCrystFELReader() override; + + void ReadFiles(const std::string &path, const std::string &geometry_file); + + // False where the geometry states no gain, so the images are in detector units, not photons. + [[nodiscard]] bool CountsInPhotons() const; + + [[nodiscard]] uint64_t GetNumberOfImages() const override; + void Close() override; + bool ReadRawImage(int64_t image_number, JFJochReaderRawImage &image) override; + [[nodiscard]] std::vector ReadSpots(int64_t image) const override; +}; diff --git a/reader/JFJochMPCCDReader.cpp b/reader/JFJochMPCCDReader.cpp deleted file mode 100644 index a783fab24..000000000 --- a/reader/JFJochMPCCDReader.cpp +++ /dev/null @@ -1,350 +0,0 @@ -// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#include "JFJochMPCCDReader.h" - -#include -#include -#include -#include -#include -#include - -#include - -#include "../common/JFJochException.h" -#include "../common/Logger.h" -#include "../writer/HDF5Objects.h" - -namespace { - -constexpr size_t PANELS = 8; -constexpr size_t PANEL_FS = 512; // fast axis, the stored image's columns -constexpr size_t PANEL_SS = 1024; // slow axis: the eight sensors are stacked along it -constexpr size_t STACK_PIXELS = PANELS * PANEL_FS * PANEL_SS; - -// Pixels on each sensor that measure nothing usable, as the geometry files the SACLA Cheetah -// pipeline generates mask them: a 5-pixel frame on every side, and the last 23 rows of the slow -// axis, which sit under the amplifier shield. -constexpr size_t EDGE = 5; -constexpr size_t SHIELD_FIRST_SS = 1001; - -// The Cheetah-SACLA output is gain-corrected so that 1 eV deposited reads 0.00125 ADU - the -// adu_per_eV its geometry generator writes for every panel - about 10 ADU per 8 keV photon. -constexpr double ADU_PER_EV = 0.00125; - -bool IsMPCCDFile(const std::string &path) { - if (!std::filesystem::is_regular_file(path)) - return false; - bool ret = false; - std::unique_lock ul(hdf5_mutex); - H5E_BEGIN_TRY { - if (H5Fis_hdf5(path.c_str()) > 0) { - try { - HDF5ReadOnlyFile file(path); - ret = file.Exists("/metadata/sensor_id") && file.Exists("/metadata/posx_in_um") - && file.Exists("/metadata/angle_in_rad"); - } catch (...) { - ret = false; - } - } - } H5E_END_TRY; - return ret; -} - -std::vector CollectFiles(const std::string &path) { - std::vector ret; - if (std::filesystem::is_directory(path)) { - for (const auto &e: std::filesystem::directory_iterator(path)) { - const auto p = e.path().string(); - if (e.path().extension() == ".h5" && IsMPCCDFile(p)) - ret.push_back(p); - } - std::sort(ret.begin(), ret.end()); - } else if (IsMPCCDFile(path)) - ret.push_back(path); - return ret; -} - -// "+0.000390x -1.000000y" as CrystFEL writes a panel's fs and ss directions. -void ParseDirection(const std::string &txt, double &x, double &y) { - x = 0; - y = 0; - static const std::regex term(R"(([+-]?\s*[0-9.eE+-]*)\s*([xyz]))"); - for (auto it = std::sregex_iterator(txt.begin(), txt.end(), term); it != std::sregex_iterator(); ++it) { - std::string num = (*it)[1].str(); - num.erase(std::remove(num.begin(), num.end(), ' '), num.end()); - double v = 1.0; - if (num == "-") - v = -1.0; - else if (!num.empty() && num != "+") - v = std::stod(num); - if ((*it)[2] == "x") - x = v; - else if ((*it)[2] == "y") - y = v; - } -} - -} // namespace - -JFJochMPCCDReader::JFJochMPCCDReader() = default; -JFJochMPCCDReader::~JFJochMPCCDReader() { - Close(); -} - -bool JFJochMPCCDReader::CanRead(const std::string &path) { - if (std::filesystem::is_directory(path)) { - for (const auto &e: std::filesystem::directory_iterator(path)) - if (e.path().extension() == ".h5") - return IsMPCCDFile(e.path().string()); - return false; - } - return IsMPCCDFile(path); -} - -void JFJochMPCCDReader::ReadFiles(const std::string &path) { - Close(); - files_ = CollectFiles(path); - if (files_.empty()) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - "No SACLA MPCCD (Cheetah) files found for " + path); - - std::array panels; - double distance_mm = 0, photon_energy_eV = 0; - std::vector file_mask; - { - std::unique_lock ul(hdf5_mutex); - for (size_t f = 0; f < files_.size(); f++) { - open_files_.emplace_back(std::make_unique(files_[f])); - for (const auto &name: open_files_.back()->FindLeafs("/")) - if (name.starts_with("tag-")) - shots_.push_back({f, "/" + name}); - } - if (shots_.empty()) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - "No images (tag-* groups) in " + path); - - auto &file = *open_files_[0]; - const auto posx = file.ReadVector("/metadata/posx_in_um"); - const auto posy = file.ReadVector("/metadata/posy_in_um"); - // Degrees, whatever the name says: the sensors stand at +-90. - const auto angle = file.ReadVector("/metadata/angle_in_rad"); - const auto pixel = file.ReadVector("/metadata/pixelsizex_in_um"); - if (posx.size() != PANELS || posy.size() != PANELS || angle.size() != PANELS || pixel.size() != PANELS) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - files_[0] + " does not describe eight MPCCD sensors"); - // The same placement the SACLA pipeline's geometry generator writes for CrystFEL. - for (size_t i = 0; i < PANELS; i++) { - const double a = angle[i] * M_PI / 180.0; - panels[i].corner_x = -posx[i] / pixel[i]; - panels[i].corner_y = posy[i] / pixel[i]; - panels[i].fs_x = std::cos(a); - panels[i].fs_y = std::sin(a); - panels[i].ss_x = -std::sin(a); - panels[i].ss_y = std::cos(a); - } - distance_mm = HDF5DataSet(file, "/metadata/distance_in_mm").ReadScalar(); - photon_energy_eV = HDF5DataSet(file, shots_[0].group + "/photon_energy_ev").ReadScalar(); - if (file.Exists("/metadata/pixelmask")) - file_mask = file.ReadVector("/metadata/pixelmask"); - } - - if (!(photon_energy_eV > 0) || !(distance_mm > 0)) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - files_[0] + " states no photon energy or no camera length"); - Build(panels, distance_mm, photon_energy_eV, file_mask); -} - -void JFJochMPCCDReader::ApplyCrystFELGeometry(const std::string &geom_file) { - if (!dataset_) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "No MPCCD files read"); - - std::ifstream in(geom_file); - if (!in) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Cannot open " + geom_file); - - std::map keys; - std::string line; - while (std::getline(in, line)) { - line = line.substr(0, line.find(';')); - const auto eq = line.find('='); - if (eq == std::string::npos) - continue; - auto trim = [](std::string s) { - s.erase(0, s.find_first_not_of(" \t")); - s.erase(s.find_last_not_of(" \t\r") + 1); - return s; - }; - keys[trim(line.substr(0, eq))] = trim(line.substr(eq + 1)); - } - - std::array panels; - double coffset = 0; - for (size_t i = 0; i < PANELS; i++) { - const std::string q = "q" + std::to_string(i + 1) + "/"; - if (!keys.contains(q + "corner_x") || !keys.contains(q + "fs") || !keys.contains(q + "ss")) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - geom_file + " does not describe panels q1..q8"); - panels[i].corner_x = std::stod(keys[q + "corner_x"]); - panels[i].corner_y = std::stod(keys[q + "corner_y"]); - ParseDirection(keys[q + "fs"], panels[i].fs_x, panels[i].fs_y); - ParseDirection(keys[q + "ss"], panels[i].ss_x, panels[i].ss_y); - if (keys.contains(q + "coffset")) - coffset = std::stod(keys[q + "coffset"]); - } - if (!keys.contains("clen")) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, geom_file + " states no clen"); - const double distance_mm = (std::stod(keys["clen"]) + coffset) * 1000.0; - - std::vector file_mask; - { - std::unique_lock ul(hdf5_mutex); - if (open_files_[0]->Exists("/metadata/pixelmask")) - file_mask = open_files_[0]->ReadVector("/metadata/pixelmask"); - } - Build(panels, distance_mm, dataset_->experiment.GetIncidentEnergy_keV() * 1000.0, file_mask); -} - -// CrystFEL's frame has y up and x to the left looking downstream; the flat image has columns to the -// right and rows down. The two differ by a half turn about the beam, which keeps the handedness, so -// a position (X, Y) lands on column beam_x - X and row beam_y - Y. -void JFJochMPCCDReader::Build(const std::array &panels, double distance_mm, - double photon_energy_eV, const std::vector &file_mask) { - auto position = [&](size_t p, double fs, double ss) { - return std::pair{panels[p].corner_x + fs * panels[p].fs_x + ss * panels[p].ss_x, - panels[p].corner_y + fs * panels[p].fs_y + ss * panels[p].ss_y}; - }; - double max_x = -1e9, min_x = 1e9, max_y = -1e9, min_y = 1e9; - for (size_t p = 0; p < PANELS; p++) - for (double fs: {0.0, double(PANEL_FS)}) - for (double ss: {0.0, double(PANEL_SS)}) { - const auto [x, y] = position(p, fs, ss); - max_x = std::max(max_x, x); - min_x = std::min(min_x, x); - max_y = std::max(max_y, y); - min_y = std::min(min_y, y); - } - const double beam_x = std::ceil(max_x) + 1; - const double beam_y = std::ceil(max_y) + 1; - width_ = static_cast(beam_x - std::floor(min_x)) + 2; - height_ = static_cast(beam_y - std::floor(min_y)) + 2; - - // Everything starts as a gap; each sensor's usable pixels are then placed on top. - std::vector mask(width_ * height_, 1u << PixelMask::ModuleGapPixelBit); - destination_.assign(STACK_PIXELS, -1); - for (size_t p = 0; p < PANELS; p++) { - for (size_t ss = 0; ss < PANEL_SS; ss++) { - for (size_t fs = 0; fs < PANEL_FS; fs++) { - const size_t src = (p * PANEL_SS + ss) * PANEL_FS + fs; - const auto [x, y] = position(p, fs + 0.5, ss + 0.5); - const auto col = std::lround(beam_x - x); - const auto row = std::lround(beam_y - y); - const size_t dst = static_cast(row) * width_ + static_cast(col); - if (fs < EDGE || fs >= PANEL_FS - EDGE || ss < EDGE || ss >= SHIELD_FIRST_SS) - continue; - destination_[src] = static_cast(dst); - mask[dst] = (src < file_mask.size() && file_mask[src] != 0) ? (1u << PixelMask::ErrorPixelBit) : 0; - } - } - } - - adu_per_photon_ = ADU_PER_EV * photon_energy_eV; - - dataset_ = std::make_shared(); - dataset_->experiment = default_experiment; - DetectorSetup detector = DetDECTRIS(static_cast(width_), static_cast(height_), "SACLA MPCCD", {}); - detector.PixelSize_um(50.0f); - // A CCD: nothing for a parallax or sensor-efficiency correction to work with is stated in the file. - detector.SensorThickness_um(0); - // The pipeline stores signed 16-bit ADU, so that is where a pixel stops counting. - detector.SaturationLimit(static_cast(std::floor(INT16_MAX / adu_per_photon_))); - detector.BitDepthImage(32); - detector.MinFrameTime(std::chrono::microseconds(0)); - detector.MinCountTime(std::chrono::microseconds(0)); - detector.ReadOutTime(std::chrono::nanoseconds(0)); - dataset_->experiment.Detector(detector); - dataset_->experiment.BeamX_pxl(static_cast(beam_x)); - dataset_->experiment.BeamY_pxl(static_cast(beam_y)); - dataset_->experiment.DetectorDistance_mm(static_cast(distance_mm)); - dataset_->experiment.IncidentEnergy_keV(static_cast(photon_energy_eV / 1000.0)); - dataset_->experiment.ImagesPerTrigger(static_cast(shots_.size())); - dataset_->pixel_mask = std::make_shared(mask); - - Logger("MPCCDReader").Info("{} shots in {} files; {} x {} image, beam at ({}, {}) px, {:.2f} mm, " - "{:.1f} eV, {:.2f} ADU per photon", shots_.size(), files_.size(), - width_, height_, beam_x, beam_y, distance_mm, photon_energy_eV, adu_per_photon_); - SetStartMessage(dataset_); -} - -uint64_t JFJochMPCCDReader::GetNumberOfImages() const { - return shots_.size(); -} - -void JFJochMPCCDReader::Close() { - std::unique_lock ul(hdf5_mutex); - open_files_.clear(); - files_.clear(); - shots_.clear(); - dataset_.reset(); -} - -// The stored chunk is the whole image, byte-shuffled and deflated. Only reading it needs the HDF5 lock; -// it is inflated, unshuffled and placed outside it, so the decoding runs on as many threads as read. -template -void JFJochMPCCDReader::DecodeInto(int64_t image_number, Buffer &buffer) const { - if (image_number < 0 || static_cast(image_number) >= shots_.size()) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Image number out of range"); - const auto &shot = shots_[image_number]; - - std::vector compressed; - { - std::unique_lock ul(hdf5_mutex); - HDF5DataSet dataset(*open_files_[shot.file], shot.group + "/data"); - dataset.ReadDirectChunk(compressed, {0, 0}); - } - - std::vector shuffled(STACK_PIXELS * sizeof(int16_t)); - uLongf out_size = shuffled.size(); - if (uncompress(shuffled.data(), &out_size, compressed.data(), compressed.size()) != Z_OK - || out_size != shuffled.size()) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - "Cannot inflate image " + shot.group + " of " + files_[shot.file]); - - buffer.assign(width_ * height_ * sizeof(int32_t), 0); - auto out = reinterpret_cast(buffer.data()); - for (size_t i = 0; i < STACK_PIXELS; i++) { - if (destination_[i] < 0) - continue; - const auto adu = static_cast(shuffled[i] | (shuffled[STACK_PIXELS + i] << 8)); - out[destination_[i]] = static_cast(std::lround(adu / adu_per_photon_)); - } -} - -bool JFJochMPCCDReader::LoadImage_i(std::shared_ptr &dataset, - DataMessage &message, - std::vector &buffer, - int64_t image_number, - bool update_dataset) { - (void) update_dataset; - if (!dataset || image_number < 0 || static_cast(image_number) >= shots_.size()) - return false; - DecodeInto(image_number, buffer); - message.image = CompressedImage(buffer.data(), buffer.size(), width_, height_, - CompressedImageMode::Int32, CompressionAlgorithm::NO_COMPRESSION); - message.number = image_number; - return true; -} - -bool JFJochMPCCDReader::ReadRawImage(int64_t image_number, JFJochReaderRawImage &image) { - if (image_number < 0 || static_cast(image_number) >= shots_.size()) - return false; - DecodeInto(image_number, image.image_buffer); - image.image = CompressedImage(image.image_buffer.data(), image.image_buffer.size(), width_, height_, - CompressedImageMode::Int32, CompressionAlgorithm::NO_COMPRESSION); - return true; -} - -std::vector JFJochMPCCDReader::ReadSpots(int64_t) const { - return {}; -} diff --git a/reader/JFJochMPCCDReader.h b/reader/JFJochMPCCDReader.h deleted file mode 100644 index bce7419c2..000000000 --- a/reader/JFJochMPCCDReader.h +++ /dev/null @@ -1,78 +0,0 @@ -// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#pragma once - -#include -#include -#include -#include - -#include "JFJochReader.h" - -class HDF5ReadOnlyFile; - -// Reads serial (stills) data from the SACLA multi-port CCD (MPCCD) octal detector, as the facility's -// Cheetah pipeline writes it: one HDF5 file per run part, one /tag-/data group per hit, each -// an int16 8192 x 512 image of the eight 1024 x 512 sensors stacked along the slow axis. The pixels -// are dark-subtracted and gain-corrected ADU, so no pedestal is needed; they are converted to -// photons here. /metadata states each sensor's position and in-plane angle, from which the sensors -// are placed on one flat image with the gaps between them masked. -// -// A path is one such file, or a directory: every Cheetah-SACLA file in it is read as one run, in -// name order. -class JFJochMPCCDReader : public JFJochReader { -public: - // One sensor as CrystFEL describes a panel: the position of its first pixel's corner and the - // directions of its fast and slow axes, in pixels, in CrystFEL's laboratory frame. - struct Panel { - double corner_x = 0, corner_y = 0; - double fs_x = 0, fs_y = 0, ss_x = 0, ss_y = 0; - }; - -private: - struct Shot { - size_t file; - std::string group; - }; - - std::vector files_; - std::vector> open_files_; - std::vector shots_; - std::shared_ptr dataset_; - - size_t width_ = 0, height_ = 0; - // For every pixel of the 8192 x 512 stack, where it lands on the flat image (-1: not placed). - std::vector destination_; - double adu_per_photon_ = 1.0; - - void Build(const std::array &panels, double distance_mm, double photon_energy_eV, - const std::vector &file_mask); - - bool LoadImage_i(std::shared_ptr &dataset, - DataMessage &message, - std::vector &buffer, - int64_t image_number, - bool update_dataset) override; - - template - void DecodeInto(int64_t image_number, Buffer &buffer) const; - -public: - JFJochMPCCDReader(); - ~JFJochMPCCDReader() override; - - // True for a Cheetah-SACLA MPCCD file, or a directory holding at least one. - static bool CanRead(const std::string &path); - - void ReadFiles(const std::string &path); - - // Replaces the sensor placement and the camera length /metadata states with a refined CrystFEL - // geometry file for the same 512 x 8192 layout (panels q1..q8). - void ApplyCrystFELGeometry(const std::string &geom_file); - - [[nodiscard]] uint64_t GetNumberOfImages() const override; - void Close() override; - bool ReadRawImage(int64_t image_number, JFJochReaderRawImage &image) override; - [[nodiscard]] std::vector ReadSpots(int64_t image) const override; -}; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index dc46c49f6..956ff62a8 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -28,7 +28,7 @@ #include "../reader/JFJochCBFReader.h" #include "../reader/JFJochMarCCDReader.h" #include "../reader/JFJochSMVReader.h" -#include "../reader/JFJochMPCCDReader.h" +#include "../reader/JFJochCrystFELReader.h" #include "../reader/JFJochStartMessageReader.h" #include "BrokerFeed.h" #include "../writer/HDF5NXmx.h" @@ -149,7 +149,7 @@ void print_usage() { std::cout << std::endl; std::cout << " Detector geometry" << std::endl; - std::cout << " --crystfel-geometry SACLA MPCCD input only: place the eight sensors and set the camera length from this refined CrystFEL geometry file (panels q1..q8 of the 512x8192 layout) instead of the nominal placement the files' /metadata states" << std::endl; + std::cout << " --crystfel-geometry Read - an HDF5/CXI file, or a directory of them - as CrystFEL does with this geometry file: frames from its data path, panels placed at its corners, camera length, pixel size, gain, photon energy and bad regions from it. One wavelength and one camera length for the run; panel positions are taken as given" << std::endl; std::cout << std::endl; std::cout << " Detector mask" << std::endl; @@ -994,7 +994,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { float false_pixels_per_frame = SpotFindingSettings{}.false_pixels_per_frame; std::string ref_mtz; std::string ref_column; - std::string crystfel_geometry; // --crystfel-geometry: refined MPCCD sensor placement + std::string crystfel_geometry; // --crystfel-geometry: read the input as CrystFEL would with this geometry std::string model_pdb; // --model: PDB to validate merged intensities against (R-free + maps) std::string mr_model; // --mr: --model, placed by molecular replacement straight away bool export_unmerged = true; // on by default; --no-export-unmerged turns it off @@ -1729,20 +1729,25 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { SimulatedLive cbf_reader; SimulatedLive marccd_reader; SimulatedLive smv_reader; - SimulatedLive mpccd_reader; + SimulatedLive crystfel_reader; JFJochStartMessageReader live_reader; std::unique_ptr broker_feed; JFJochReader *reader_ptr = nullptr; // Which format this is. Asked INSIDE the try: these open and read the file - a directory whose // first plausible member cannot be read throws here - and an input the program cannot read has to // come out as the message below, not as an unhandled exception with nothing printed at all. - bool input_is_cbf = false, input_is_marccd = false, input_is_smv = false, input_is_mpccd = false; + bool input_is_cbf = false, input_is_marccd = false, input_is_smv = false, input_is_crystfel = false; try { - input_is_cbf = JFJochCBFReader::CanRead(input_file); + input_is_crystfel = !crystfel_geometry.empty(); + input_is_cbf = !input_is_crystfel && JFJochCBFReader::CanRead(input_file); input_is_marccd = !input_is_cbf && JFJochMarCCDReader::CanRead(input_file); input_is_smv = !input_is_cbf && !input_is_marccd && JFJochSMVReader::CanRead(input_file); - input_is_mpccd = !input_is_cbf && !input_is_marccd && !input_is_smv && JFJochMPCCDReader::CanRead(input_file); - if (input_is_cbf) { + if (input_is_crystfel) { + crystfel_reader.ReadFiles(input_file, crystfel_geometry); + reader_ptr = &crystfel_reader; + logger.Info("Read {} frames from {} with the CrystFEL geometry {}", + crystfel_reader.GetNumberOfImages(), input_file, crystfel_geometry); + } else if (input_is_cbf) { cbf_reader.ReadFiles(input_file); reader_ptr = &cbf_reader; logger.Info("Read {} CBF images of the sweep containing {}", @@ -1757,12 +1762,6 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { reader_ptr = &smv_reader; logger.Info("Read {} SMV images of the sweep containing {}", smv_reader.GetNumberOfImages(), input_file); - } else if (input_is_mpccd) { - mpccd_reader.ReadFiles(input_file); - if (!crystfel_geometry.empty()) - mpccd_reader.ApplyCrystFELGeometry(crystfel_geometry); - reader_ptr = &mpccd_reader; - logger.Info("Read {} SACLA MPCCD images from {}", mpccd_reader.GetNumberOfImages(), input_file); } else if (std::filesystem::is_directory(input_file)) { // A directory is never an HDF5 master, so handing it to the HDF5 reader only buries the // real answer - that no one-file-per-image sweep was found in it - under a page of HDF5 @@ -1937,7 +1936,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { "runs no space-group search. Use --mode mx, which makes the P1 merge the " "ledger is folded from."); // Re-scaling reads reflections a previous run integrated, which only a _process.h5 holds. - if (input_is_cbf || input_is_marccd || input_is_smv || input_is_mpccd) { + if (input_is_cbf || input_is_marccd || input_is_smv || input_is_crystfel) { logger.Error("--mode scale needs the integrated reflections in a _process.h5; " "a raw image sweep holds none. Run --mode mx on it first."); exit(EXIT_FAILURE); @@ -3168,7 +3167,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { config.detect_beam_stop = detect_beam_stop; config.beam_center_check = beam_center_check; config.adaptive_integration_radius = adaptive_integration_radius; - config.counts_in_adu = input_is_smv || input_is_marccd; + config.counts_in_adu = input_is_smv || input_is_marccd || (input_is_crystfel && !crystfel_reader.CountsInPhotons()); config.rotation_postrefine_geometry = rotation_postrefine_geometry; // Measure the spot budget from the data unless the user pinned it. config.measure_spot_budget = !max_spot_count_override.has_value();