Jungfraujoch/reader/JFJochMarCCDReader.cpp

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#include "JFJochMarCCDReader.h"

#include <cmath>
#include <future>
#include <thread>

#include "../common/JFJochException.h"
#include "../common/JFJochMath.h"
#include "../common/Logger.h"

namespace {

// The base rotation axis in the internal frame (x along increasing detector column, y along
// increasing row, z along the beam). A marCCD header names the circle that turned but never states
// a direction for it, so this is the convention an NXmx master writes for the same instruments, and
// a file that needs the other sign is settled from the data by the run's axis-sign rescue - the
// same arrangement JFJochCBFReader makes for a miniCBF that states no axis table.
const Coord ASSUMED_BASE_AXIS(-1.0f, 0.0f, 0.0f);

}  // namespace

bool JFJochMarCCDReader::CanRead(const std::string &path) {
    return marccd::CanRead(path);
}

void JFJochMarCCDReader::ReadFiles(const std::string &path) {
    files_ = marccd::CollectSweep(path);
    if (files_.empty())
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "No marCCD images found for " + path);

    header0_ = marccd::ReadHeader(files_[0]);
    // A pixel size is what makes the file an IMAGE: every resolution, every scattering vector and
    // the beam centre in millimetres scale by it, and a default of 0 collapses all of them without
    // a word.
    if (!(header0_.pixel_x_m > 0.0) || !(header0_.pixel_y_m > 0.0))
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              files_[0] + " states no pixel size in its marCCD header");
    if (!(header0_.wavelength_A > 0.0))
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              files_[0] + " states no wavelength in its marCCD header");
    if (!(header0_.distance_m > 0.0))
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              files_[0] + " states no detector distance in its marCCD header");

    dataset_ = std::make_shared<JFJochReaderDataset>();
    dataset_->experiment = default_experiment;

    DetectorSetup detector = DetDECTRIS(header0_.nx, header0_.ny,
                                        header0_.detector.empty() ? "marCCD" : header0_.detector, {});
    // Not rounded to whole micrometres, as the miniCBF path can afford to be: a PILATUS pixel is
    // exactly 172 um, but these are 73.242 um, and rounding that to 73 is a 0.33% scale error on
    // every cell edge the run reports.
    detector.PixelSize_um(static_cast<float>(header0_.pixel_x_m * 1e6));
    // A CCD has no sensor thickness worth correcting for: the phosphor converts at the surface and
    // the fibre optic carries light, not X-rays, so the parallax correction a silicon sensor needs
    // does not apply. Left at zero, which is what the geometry means by "no depth".
    detector.SensorThickness_um(0);
    if (header0_.saturated_value > 0)
        detector.SaturationLimit(SaturationLimitFromValue(header0_.saturated_value));
    else
        Logger("MarCCDReader").Warning("{} states no saturated value, so no pixel will be called "
                                       "saturated; if this detector overloads, its strongest "
                                       "reflections will be integrated as if they were valid.",
                                       files_[0]);
    // Images are handed out as signed 32-bit whatever the file stored, so that is the depth the
    // rest of the code must see; the real overflow is the header's saturated value, set above.
    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<float>(header0_.beam_x_px));
    dataset_->experiment.BeamY_pxl(static_cast<float>(header0_.beam_y_px));
    dataset_->experiment.DetectorDistance_mm(static_cast<float>(header0_.distance_m * 1000.0));

    // A detector swung out on a 2theta arm. The arm turns the detector about the sample and so
    // carries the square-on geometry with it: the header's distance stays the distance along the
    // detector normal and the beam centre stays the point of normal incidence, which is exactly
    // what the PONI convention wants, so the swing is a PONI rotation and nothing else changes.
    // The arm turns about the same axis as the spindle on the geometries these headers describe.
    if (header0_.two_theta_deg != 0.0) {
        float rot1 = 0, rot2 = 0, rot3 = 0;
        PoniAnglesFromMatrix(RotMatrix(static_cast<float>(header0_.two_theta_deg * PI / 180.0),
                                       ASSUMED_BASE_AXIS), rot1, rot2, rot3);
        dataset_->experiment.PoniRot1_rad(rot1).PoniRot2_rad(rot2).PoniRot3_rad(rot3);
    }

    dataset_->experiment.IncidentEnergy_keV(WVL_1A_IN_KEV / static_cast<float>(header0_.wavelength_A));
    dataset_->experiment.FrameTime(
        std::chrono::duration_cast<std::chrono::nanoseconds>(
            std::chrono::duration<double>(header0_.exposure_s)),
        std::chrono::duration_cast<std::chrono::nanoseconds>(
            std::chrono::duration<double>(header0_.exposure_s)));

    // The rotation angle of every image, from its own header. Reading one costs a 4 kB read, so on
    // a sweep of several thousand frames this is worth spreading over the cores, as the CBF path
    // does for the same reason.
    std::vector<double> angles(files_.size());
    {
        const size_t nthreads = std::min<size_t>(std::max(1u, std::thread::hardware_concurrency()), 8);
        std::vector<std::future<void>> futures;
        for (size_t t = 0; t < nthreads; t++)
            futures.push_back(std::async(std::launch::async, [&, t] {
                for (size_t i = t; i < files_.size(); i += nthreads)
                    angles[i] = marccd::ReadHeader(files_[i]).start_angle_deg;
            }));
        for (auto &f : futures)
            f.get();
    }

    double increment = header0_.angle_increment_deg;
    if (files_.size() > 1) {
        // Prefer the measured step over the header's nominal one, and unwrap a sweep that passes 360.
        double d = angles[1] - angles[0];
        if (d < -180.0) d += 360.0;
        if (std::abs(d) > 1e-6) increment = d;
    }
    dataset_->experiment.Goniometer(GoniometerAxis(header0_.axis_name,
                                                   static_cast<float>(angles.front()),
                                                   static_cast<float>(increment),
                                                   ASSUMED_BASE_AXIS, {}));

    dataset_->error_value = -1;
    dataset_->experiment.ImagesPerTrigger(static_cast<int64_t>(files_.size()));
    // A CCD frame stores no untrusted-pixel marker - every value is a real reading, and the
    // detector has no module gaps - so the sweep starts with nothing masked.
    dataset_->pixel_mask = std::make_shared<const PixelMask>(static_cast<size_t>(header0_.nx),
                                                            static_cast<size_t>(header0_.ny));

    SetStartMessage(dataset_);
}

uint64_t JFJochMarCCDReader::GetNumberOfImages() const {
    return files_.size();
}

void JFJochMarCCDReader::Close() {
    files_.clear();
    dataset_.reset();
}

template <class Buffer>
CompressedImage JFJochMarCCDReader::DecodeInto(int64_t image_number, Buffer &buffer,
                                               std::vector<uint8_t> &scratch) const {
    if (image_number < 0 || static_cast<size_t>(image_number) >= files_.size())
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "Image number out of range");

    const size_t npixel = static_cast<size_t>(header0_.nx) * static_cast<size_t>(header0_.ny);
    buffer.resize(npixel * sizeof(int32_t));

    const auto h = marccd::ReadInto(files_[image_number],
                                    reinterpret_cast<int32_t *>(buffer.data()), npixel, scratch);
    if (h.nx != header0_.nx || h.ny != header0_.ny)
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "marCCD image size differs from the first image of the sweep");

    return CompressedImage(buffer.data(), buffer.size(),
                           static_cast<size_t>(header0_.nx), static_cast<size_t>(header0_.ny),
                           CompressedImageMode::Int32, CompressionAlgorithm::NO_COMPRESSION);
}

bool JFJochMarCCDReader::LoadImage_i(std::shared_ptr<JFJochReaderDataset> &dataset,
                                     DataMessage &message,
                                     std::vector<uint8_t> &buffer,
                                     int64_t image_number,
                                     bool update_dataset) {
    (void) update_dataset;
    if (!dataset)
        return false;

    std::vector<uint8_t> scratch;
    message.image = DecodeInto(image_number, buffer, scratch);
    message.number = image_number;
    return true;
}

bool JFJochMarCCDReader::ReadRawImage(int64_t image_number, JFJochReaderRawImage &image) {
    image.image = DecodeInto(image_number, image.image_buffer, image.read_buffer);
    return true;
}

std::vector<SpotToSave> JFJochMarCCDReader::ReadSpots(int64_t) const {
    return {};   // a raw marCCD file stores no analysis results
}