From 923a21014270ec43404fe87d647e4ce889f30d70 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 20 Oct 2026 12:04:39 +0200 Subject: [PATCH] WIP: SACLA MPCCD (Cheetah-SACLA HDF5) stills reader for rugnux - not built, not tested Reads the per-hit /tag-*/data int16 8192x512 stacks (eight 1024x512 sensors), places the sensors on one flat image from /metadata posx/posy/angle (or a refined CrystFEL .geom via --crystfel-geometry), masks gaps, 5 px sensor edges and the amplifier-shield rows, converts ADU to photons at 0.00125 ADU/eV, and decodes chunks (shuffle+deflate) outside the HDF5 lock. Co-Authored-By: Claude Opus 5.5 --- reader/CMakeLists.txt | 2 + reader/JFJochMPCCDReader.cpp | 350 +++++++++++++++++++++++++++++++++++ reader/JFJochMPCCDReader.h | 78 ++++++++ rugnux/rugnux_cli.cpp | 23 ++- 4 files changed, 451 insertions(+), 2 deletions(-) create mode 100644 reader/JFJochMPCCDReader.cpp create mode 100644 reader/JFJochMPCCDReader.h diff --git a/reader/CMakeLists.txt b/reader/CMakeLists.txt index 0d6350514..16f137359 100644 --- a/reader/CMakeLists.txt +++ b/reader/CMakeLists.txt @@ -12,6 +12,8 @@ ADD_LIBRARY(JFJochReader STATIC MarCCD.h JFJochMarCCDReader.cpp JFJochMarCCDReader.h + JFJochMPCCDReader.cpp + JFJochMPCCDReader.h SMV.cpp SMV.h JFJochSMVReader.cpp diff --git a/reader/JFJochMPCCDReader.cpp b/reader/JFJochMPCCDReader.cpp new file mode 100644 index 000000000..a783fab24 --- /dev/null +++ b/reader/JFJochMPCCDReader.cpp @@ -0,0 +1,350 @@ +// 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 new file mode 100644 index 000000000..bce7419c2 --- /dev/null +++ b/reader/JFJochMPCCDReader.h @@ -0,0 +1,78 @@ +// 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 2a1ff9f88..dc46c49f6 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -28,6 +28,7 @@ #include "../reader/JFJochCBFReader.h" #include "../reader/JFJochMarCCDReader.h" #include "../reader/JFJochSMVReader.h" +#include "../reader/JFJochMPCCDReader.h" #include "../reader/JFJochStartMessageReader.h" #include "BrokerFeed.h" #include "../writer/HDF5NXmx.h" @@ -147,6 +148,10 @@ void print_usage() { std::cout << " --no-refine-tilt Do not refine the detector tilt: hold rot1/rot2 at the header value and fit only the beam centre and the distance. For a calibration handed to a program that cannot express a tilted detector, e.g. XDS" << std::endl; 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 << std::endl; + std::cout << " Detector mask" << std::endl; std::cout << " --detect-beam-stop[=N|off] Find the beam stop and its holder in a projection of N images and add them to the pixel mask (bit 9), so nothing shadowed by them is integrated. ON by default (60 images); =off disables. Reflections behind the stop are attenuated but not flagged, so they are integrated low with a plausible sigma and no existing rejection catches them" << std::endl; std::cout << " --beam-center-check[=off] Measure the beam centre from the isotropy of the scattered background on EVERY run, report how far the file's value is from it, and index a second first pass with it to see whether the two centres give the same lattice. The fit reads the projection --detect-beam-stop already builds, so it costs no extra frames. On a run that indexes, the measured centre is adopted where the file's centre indexes nothing and the measured one indexes a majority, and where the two centres return cells related by an integer volume factor, in either direction, and the measured centre's lattice carries materially more of the pooled spots against its own chance level - the axis harmonic a centre error along the spindle produces, which nothing downstream repairs. Where the indexing cannot tell the two centres apart - the same cell with a different metric symmetry, or the same lattice found at both only by reading less of each frame - the first pass is run at both and the run keeps the one that merges better (on the same lattice: the one that measures more reflections at I/sigma >= 2). ON by default; =off disables" << std::endl; @@ -303,6 +308,7 @@ enum { OPT_REJECT_OUTLIERS, OPT_REFERENCE_COLUMN, OPT_MODEL, + OPT_CRYSTFEL_GEOMETRY, OPT_MR, OPT_EXPORT_UNMERGED, OPT_NO_EXPORT_UNMERGED, @@ -369,6 +375,7 @@ static option long_options[] = { {"reference-mtz", required_argument, nullptr, 'z'}, {"reference-column", required_argument, nullptr, OPT_REFERENCE_COLUMN}, {"model", required_argument, nullptr, OPT_MODEL}, + {"crystfel-geometry", required_argument, nullptr, OPT_CRYSTFEL_GEOMETRY}, {"mr", required_argument, nullptr, OPT_MR}, {"export-unmerged", no_argument, nullptr, OPT_EXPORT_UNMERGED}, {"no-export-unmerged", no_argument, nullptr, OPT_NO_EXPORT_UNMERGED}, @@ -987,6 +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 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 @@ -1215,6 +1223,9 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { case OPT_MODEL: model_pdb = optarg; break; + case OPT_CRYSTFEL_GEOMETRY: + crystfel_geometry = optarg; + break; case OPT_MR: mr_model = optarg; break; @@ -1718,17 +1729,19 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { SimulatedLive cbf_reader; SimulatedLive marccd_reader; SimulatedLive smv_reader; + SimulatedLive mpccd_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; + bool input_is_cbf = false, input_is_marccd = false, input_is_smv = false, input_is_mpccd = false; try { input_is_cbf = 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) { cbf_reader.ReadFiles(input_file); reader_ptr = &cbf_reader; @@ -1744,6 +1757,12 @@ 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 @@ -1918,7 +1937,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) { + if (input_is_cbf || input_is_marccd || input_is_smv || input_is_mpccd) { 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);