Jungfraujoch/common/BraggIntegrationSettings.cpp

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

#include <cmath>

#include "BraggIntegrationSettings.h"
#include "JFJochException.h"

#define check_max(param, val, max) if ((val) > (max)) throw JFJochException(JFJochExceptionCategory::InputParameterAboveMax, param)
#define check_min(param, val, min) if ((val) < (min)) throw JFJochException(JFJochExceptionCategory::InputParameterBelowMin, param)
#define check_finite(param, val) if (!std::isfinite(val)) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, param)

BraggIntegrationSettings &BraggIntegrationSettings::R1(float input) {
    check_finite("Integration radius R1", input);
    check_min("Integration radius R1", input, 0.1);
    check_max("Integration radius R1", input, 20.0);
    r_1 = input;
    return *this;
}

BraggIntegrationSettings &BraggIntegrationSettings::R2(float input) {
    check_finite("Background inner radius R2", input);
    check_min("Background inner radius R2", input, 0.1);
    check_max("Background inner radius R2", input, 30.0);

    if (input <= r_1)
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "Background inner radius (R2) must be larger than integration radius (R1)");

    r_2 = input;
    return *this;
}

BraggIntegrationSettings &BraggIntegrationSettings::R3(float input) {
    check_finite("Background outer radius R3", input);
    check_min("Background outer radius R3", input, 0.1);
    check_max("Background outer radius R3", input, 40.0);

    if (input <= r_2)
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "Background outer radius (R3) must be larger than background inner radius (R2)");

    r_3 = input;
    return *this;
}

BraggIntegrationSettings &BraggIntegrationSettings::StencilKSigma(float input) {
    check_finite("Integration stencil elongation", input);
    check_min("Integration stencil elongation", input, 0.0);
    check_max("Integration stencil elongation", input, 10.0);
    stencil_k_sigma = input;
    return *this;
}

BraggIntegrationSettings &BraggIntegrationSettings::DMinLimit_A(std::optional<float> input) {
    if (input) {
        check_finite("Minimum d-spacing", *input);
        check_min("Minimum d-spacing", *input, 0.5);
        check_max("Minimum d-spacing", *input, 100.0);
    }
    d_min_limit_A = input;
    return *this;
}

BraggIntegrationSettings & BraggIntegrationSettings::FixedProfileRadius_recipA(std::optional<float> input) {
    if (input) {
        check_finite("Profile radius", input.value());
        check_min("Profile radius [A^-1]", input.value(), 0.000001);
        check_max("Profile radius [A^-1]", input.value(), 0.01);
    }
    fixed_profile_radius = input;
    return *this;

}

std::optional<float> BraggIntegrationSettings::GetFixedProfileRadius_recipA() const {
    return fixed_profile_radius;
}

BraggIntegrationSettings & BraggIntegrationSettings::ForcedPredictionMosaicity_deg(std::optional<float> input) {
    if (input) {
        check_finite("Prediction mosaicity", input.value());
        check_min("Prediction mosaicity [deg]", input.value(), 0.001);
        check_max("Prediction mosaicity [deg]", input.value(), 10.0);
    }
    forced_prediction_mosaicity_deg = input;
    return *this;
}

std::optional<float> BraggIntegrationSettings::GetForcedPredictionMosaicity_deg() const {
    return forced_prediction_mosaicity_deg;
}

BraggIntegrationSettings &BraggIntegrationSettings::Integrator(IntegratorMode input) {
    integrator_mode = input;
    return *this;
}

IntegratorMode BraggIntegrationSettings::GetIntegrator() const {
    return integrator_mode;
}

float BraggIntegrationSettings::GetR1() const {
    return r_1;
}

float BraggIntegrationSettings::GetR2() const {
    return r_2;
}

float BraggIntegrationSettings::GetR3() const {
    return r_3;
}

float BraggIntegrationSettings::GetStencilKSigma() const {
    return stencil_k_sigma;
}

std::optional<float> BraggIntegrationSettings::GetDMinLimit_A() const {
    return d_min_limit_A;
}

float BraggIntegrationSettings::GetMinimumSigmaInRegardsToI() const {
    return minimum_sigma_in_regards_to_i;
}

BraggIntegrationSettings &BraggIntegrationSettings::BackgroundTrimFraction(float input) {
    check_finite("Background trim fraction", input);
    check_min("Background trim fraction", input, 0.0);
    check_max("Background trim fraction", input, 0.49);   // must leave a central majority after trimming
    bkg_trim_fraction = input;
    if (input > 0.0f)
        bkg_clip_nsigma = 0.0f;      // the two ring estimators are alternatives, not a stack
    return *this;
}

float BraggIntegrationSettings::GetBackgroundTrimFraction() const {
    return bkg_trim_fraction;
}

BraggIntegrationSettings &BraggIntegrationSettings::MaxHKL(std::optional<int> input) {
    if (input) {
        check_min("Maximum hkl index", *input, 1);
        // The GPU predictor launches one thread per candidate, so the cost is (2n+1)^3: 511 is 1.1e9
        // candidates per frame, already far past the point where prediction dominates a run.
        check_max("Maximum hkl index", *input, 511);
    }
    max_hkl = input;
    return *this;
}

std::optional<int> BraggIntegrationSettings::GetMaxHKL() const {
    return max_hkl;
}

BraggIntegrationSettings &BraggIntegrationSettings::BackgroundClipNSigma(float input) {
    check_finite("Background clip nsigma", input);
    check_min("Background clip nsigma", input, 0.0);
    bkg_clip_nsigma = input;
    if (input > 0.0f)
        bkg_trim_fraction = 0.0f;    // the two ring estimators are alternatives, not a stack
    return *this;
}

float BraggIntegrationSettings::GetBackgroundClipNSigma() const {
    return bkg_clip_nsigma;
}

BraggIntegrationSettings &BraggIntegrationSettings::BackgroundRadialCorrection(std::optional<bool> input) {
    bkg_radial_correction = input;
    return *this;
}

std::optional<bool> BraggIntegrationSettings::GetBackgroundRadialCorrection() const {
    return bkg_radial_correction;
}

BraggIntegrationSettings &BraggIntegrationSettings::Overlap(OverlapMode input) {
    overlap_mode = input;
    return *this;
}

OverlapMode BraggIntegrationSettings::GetOverlap() const {
    return overlap_mode;
}

BraggIntegrationSettings &BraggIntegrationSettings::OverlapMinPeak(float input) {
    check_finite("Overlap minimum peak fraction", input);
    check_min("Overlap minimum peak fraction", input, 0.0);
    check_max("Overlap minimum peak fraction", input, 1.0);
    overlap_min_peak = input;
    return *this;
}

float BraggIntegrationSettings::GetOverlapMinPeak() const {
    return overlap_min_peak;
}

BraggIntegrationSettings &BraggIntegrationSettings::FlightPath(FlightPathMedium input) {
    flight_path = input;
    return *this;
}

FlightPathMedium BraggIntegrationSettings::GetFlightPath() const {
    return flight_path;
}

BraggIntegrationSettings &BraggIntegrationSettings::Footprint(const SpotFootprint &input) {
    if (input.sigma_rad.size() != input.sigma_tan.size()
        || input.sigma_rad.size() > static_cast<size_t>(SpotFootprint::MAX_BINS)
        || (!input.empty() && !(input.bin_px > 0.0f)))
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Invalid spot footprint table");
    footprint = input;
    return *this;
}

const SpotFootprint &BraggIntegrationSettings::GetFootprint() const {
    return footprint;
}