Jungfraujoch/common/AzimuthalIntegrationProfile.cpp

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

#include "AzimuthalIntegrationProfile.h"
#include "JFJochException.h"

inline float sum_to_count(float sum, uint64_t count) {
    if (count == 0)
        return NAN;
    else
        return sum / (static_cast<float>(count));
}

AzimuthalIntegrationProfile::AzimuthalIntegrationProfile(const AzimuthalIntegration &mapping)
    : sum(mapping.GetBinNumber(), 0),
      sum2(mapping.GetBinNumber(), 0),
      count(mapping.GetBinNumber(), 0),
      bin_to_q(mapping.GetBinToQ()),
      bin_to_d(mapping.GetBinToD()),
      bin_to_2theta(mapping.GetBinToTwoTheta()),
      bin_to_phi(mapping.GetBinToPhi()),
      q_bins(mapping.GetQBinCount()),
      azim_bins(mapping.GetAzimuthalBinCount()) {
}

void AzimuthalIntegrationProfile::Clear(const AzimuthalIntegration &mapping) {
    std::unique_lock ul(m);

    bin_to_d = mapping.GetBinToD();
    bin_to_q = mapping.GetBinToQ();
    bin_to_2theta = mapping.GetBinToTwoTheta();
    bin_to_phi = mapping.GetBinToPhi();
    q_bins = mapping.GetQBinCount();
    azim_bins = mapping.GetAzimuthalBinCount();

    sum = std::vector<float>(mapping.GetBinNumber(), 0);
    count = std::vector<uint64_t>(mapping.GetBinNumber(), 0);
}

void AzimuthalIntegrationProfile::Add(int64_t bin, int64_t value) {
    if (bin < 0 || bin >= sum.size())
        return;
    std::unique_lock ul(m);
    sum[bin] += static_cast<float>(value);
    sum2[bin] += static_cast<float>(value * value);
    count[bin]++;
}

void AzimuthalIntegrationProfile::Add(const std::vector<float> &in_sum, const std::vector<uint32_t> &in_count) {
    std::unique_lock ul(m);
    if ((in_sum.size() == sum.size()) && (in_count.size() == count.size())) {
        for (int i = 0; i < sum.size(); i++) {
            sum[i] += in_sum[i];
            count[i] += in_count[i];
        }
    } else if (!in_sum.empty() && !in_count.empty())
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Mismatch in size of sum/count datasets");
}

std::vector<float> AzimuthalIntegrationProfile::GetResult() const {
    std::unique_lock ul(m);

    std::vector<float> rad_int_profile(sum.size(), 0);
    for (int i = 0; i < sum.size(); i++)
        rad_int_profile[i] = sum_to_count(sum[i], count[i]);
    return rad_int_profile;
}

std::vector<float> AzimuthalIntegrationProfile::GetResult1D() const {
    std::unique_lock ul(m);

    std::vector<float> rad_int_profile(q_bins, 0);
    for (int i = 0; i < sum.size(); i++)
        rad_int_profile[i % q_bins] = sum_to_count(sum[i], count[i]);
    return rad_int_profile;
}

void AzimuthalIntegrationProfile::SetTitle(const std::string &input) {
    title = input;
}

const std::vector<float> &AzimuthalIntegrationProfile::GetXAxis(PlotAzintUnit unit) const {
    switch (unit) {
        case PlotAzintUnit::TwoTheta_deg:
            return bin_to_2theta;
        case PlotAzintUnit::D_A:
            return bin_to_d;
        default:
        case PlotAzintUnit::Q_recipA:
            return bin_to_q;
    }
}

MultiLinePlot AzimuthalIntegrationProfile::GetPlot(bool force_1d, PlotAzintUnit plot_unit) const {
    MultiLinePlot ret;

    const std::vector<float> &x_coord = GetXAxis(plot_unit);

    if (azim_bins == 1)
        ret.AddPlot(MultiLinePlotStruct{.title = title, .x = x_coord, .y = GetResult()});
    else {
        if (force_1d) {
            std::vector<float> x_shortened(q_bins);
            for (int i = 0; i < q_bins; i++)
                x_shortened[i] = x_coord[i];
            ret.AddPlot(MultiLinePlotStruct{.title = title, .x = x_shortened, .y = GetResult1D()});
        } else {
            ret.AddPlot(MultiLinePlotStruct{.title = title, .x = x_coord, .y= bin_to_phi, .z = GetResult()});
        }
    }
    return ret;
}

float AzimuthalIntegrationProfile::GetMeanValueOfBins(uint16_t min_bin, uint16_t max_bin) const {
    std::unique_lock ul(m);

    float ret_sum = 0;
    uint64_t ret_count = 0;

    for (int i = 0; i < sum.size(); i++) {
        uint16_t q_bin = i % q_bins;
        if (q_bin >= min_bin && q_bin <= max_bin) {
            ret_sum += sum[i];
            ret_count += count[i];
        }
    }

    return sum_to_count(ret_sum, ret_count);
}

float AzimuthalIntegrationProfile::GetBkgEstimate(const AzimuthalIntegrationSettings &settings) const {
    auto min_bin = settings.QToBin(settings.GetBkgEstimateLowQ_recipA());
    auto max_bin = settings.QToBin(settings.GetBkgEstimateHighQ_recipA());
    return GetMeanValueOfBins(min_bin, max_bin);
}

AzimuthalIntegrationProfile &AzimuthalIntegrationProfile::operator+=(const AzimuthalIntegrationProfile &other) {
    if ((other.bin_to_q != bin_to_q) || (sum.size() != other.sum.size())) {
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "Error combining two radial integration profiles");
    }

    for (int i = 0; i < sum.size(); i++) {
        sum[i] += other.sum[i];
        count[i] += other.count[i];
    }

    return *this;
}

void AzimuthalIntegrationProfile::Add(const DeviceOutput &result) {
    std::unique_lock ul(m);

    if (sum.size() > FPGA_INTEGRATION_BIN_COUNT )
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "Error in getting result from FPGA");

    for (int i = 0; i < sum.size(); i++) {
        sum[i] += result.integration_result[i].sum;
        count[i] += result.integration_result[i].count;
    }
}