detectorfluxanalysis2d.cpp

Go to the documentation of this file.

/**
 * @file detectorfluxanalysis2d.cpp
 * @brief Implements the detectorfluxanalysis2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/detectorfluxanalysis2d.cpp. This file is part of Quantum Phenomena Explorer (QPE),
 * maintained by Officina Turini. See LICENSE_IT.txt / LICENSE.txt for the
 * QPE freeware source-available terms and THIRD_PARTY_NOTICES.txt for
 * third-party components.
 */
 
#include "detectorfluxanalysis2d.h"
 
#include <algorithm>
#include <cmath>
 
namespace qpe {
namespace {
 
/**
 * @brief Returns normalize.
 * @param values The values value used by the operation.
 */
bool normalize(std::vector<double> &values)
{
    double sum = 0.0;
    for (std::size_t i = 0u; i < values.size(); ++i) {
        if (!std::isfinite(values[i]) || values[i] < 0.0) {
            values[i] = 0.0;
        }
        sum += values[i];
    }
    if (!(sum > 0.0)) {
        return false;
    }
    const double inv = 1.0 / sum;
    for (std::size_t i = 0u; i < values.size(); ++i) {
        values[i] *= inv;
    }
    return true;
}
 
std::size_t profileIndexToBin(std::size_t index,
                              std::size_t sampleCount,
                              std::size_t bins)
{
    if (sampleCount <= 1u || bins <= 1u) {
        return 0u;
    }
    const double u = static_cast<double>(index) /
                     static_cast<double>(sampleCount - 1u);
    std::size_t bin = static_cast<std::size_t>(u * static_cast<double>(bins));
    if (bin >= bins) {
        bin = bins - 1u;
    }
    return bin;
}
 
} // namespace
 
std::vector<double> DetectorFluxAnalysis2D::rebinNormalizedProfile(
    const std::vector<double> &fluxProfile,
    std::size_t bins)
{
    if (bins == 0u || fluxProfile.empty()) {
        return std::vector<double>();
    }
 
    std::vector<double> result(bins, 0.0);
    for (std::size_t i = 0u; i < fluxProfile.size(); ++i) {
        const double v = fluxProfile[i];
        if (std::isfinite(v) && v > 0.0) {
            result[profileIndexToBin(i, fluxProfile.size(), bins)] += v;
        }
    }
    if (!normalize(result)) {
        result.clear();
    }
    return result;
}
 
DetectorFluxAnalysis2DResult DetectorFluxAnalysis2D::analyze(
    const std::vector<double> &fluxProfile,
    float dy,
    const std::vector<PositionMeasurement2DPoint> &hits,
    std::size_t bins)
{
    DetectorFluxAnalysis2DResult result;
    if (bins == 0u || fluxProfile.empty() || !(dy > 0.0f)) {
        return result;
    }
 
    result.expectedProbability = rebinNormalizedProfile(fluxProfile, bins);
    result.expectedValid = !result.expectedProbability.empty();
 
    result.observedProbability.assign(bins, 0.0);
    if (!hits.empty()) {
        const double yMin = -0.5 * static_cast<double>(fluxProfile.size() - 1u) *
                            static_cast<double>(dy);
        const double yMax =  0.5 * static_cast<double>(fluxProfile.size() - 1u) *
                            static_cast<double>(dy);
        const double span = yMax - yMin;
        if (span > 0.0) {
            for (std::size_t i = 0u; i < hits.size(); ++i) {
                const double y = static_cast<double>(hits[i].y);
                if (y < yMin || y > yMax) {
                    continue;
                }
                double u = (y - yMin) / span;
                u = std::max(0.0, std::min(1.0, u));
                std::size_t bin = static_cast<std::size_t>(u * static_cast<double>(bins));
                if (bin >= bins) {
                    bin = bins - 1u;
                }
                result.observedProbability[bin] += 1.0;
            }
        }
    }
    result.observedValid = normalize(result.observedProbability);
    if (!result.observedValid) {
        result.observedProbability.clear();
    }
 
    if (result.expectedValid && result.observedValid) {
        double distance = 0.0;
        for (std::size_t i = 0u; i < bins; ++i) {
            distance += std::fabs(result.expectedProbability[i] -
                                  result.observedProbability[i]);
        }
        result.totalVariationDistance = 0.5 * distance;
    }
    return result;
}
 
} // namespace qpe