detectorscreenmeasurement2d.cpp

Go to the documentation of this file.

/**
 * @file detectorscreenmeasurement2d.cpp
 * @brief Implements the detectorscreenmeasurement2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/detectorscreenmeasurement2d.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 "detectorscreenmeasurement2d.h"
 
#include "core/math/complexfield2d.h"
 
#include <algorithm>
#include <cmath>
#include <random>
 
namespace qpe {
namespace {
 
std::vector<PositionMeasurement2DPoint> sampleWithEngine(
    const ComplexField2D &field,
    const DetectorScreenMeasurement2DParams &params,
    std::size_t shots,
    std::mt19937 &engine)
{
    std::vector<PositionMeasurement2DPoint> result;
    if (field.empty() || shots == 0u || !(params.thickness > 0.0f)) {
        return result;
    }
 
    const double cx = 0.5 * static_cast<double>(field.width() - 1u);
    const double cy = 0.5 * static_cast<double>(field.height() - 1u);
    const double dx = static_cast<double>(field.dx());
    const double dy = static_cast<double>(field.dy());
    const double halfThickness = 0.5 * static_cast<double>(params.thickness);
 
    std::vector<double> weights(field.height(), 0.0);
    double total = 0.0;
    bool selectedAnyColumn = false;
 
    for (std::size_t x = 0; x < field.width(); ++x) {
        const double physicalX = (static_cast<double>(x) - cx) * dx;
        if (std::fabs(physicalX - static_cast<double>(params.detectorX)) > halfThickness) {
            continue;
        }
        selectedAnyColumn = true;
        for (std::size_t y = 0; y < field.height(); ++y) {
            const ComplexField2D::Value v = field(x, y);
            const double re = static_cast<double>(v.real());
            const double im = static_cast<double>(v.imag());
            const double p = re * re + im * im;
            if (std::isfinite(p) && p > 0.0) {
                weights[y] += p;
                total += p;
            }
        }
    }
 
    // For a strip thinner than one grid cell, use the nearest detector column.
    if (!selectedAnyColumn) {
        const double gridX = static_cast<double>(params.detectorX) / dx + cx;
        long nearest = static_cast<long>(std::lround(gridX));
        nearest = std::max<long>(0, std::min<long>(static_cast<long>(field.width()) - 1, nearest));
        const std::size_t x = static_cast<std::size_t>(nearest);
        for (std::size_t y = 0; y < field.height(); ++y) {
            const ComplexField2D::Value v = field(x, y);
            const double re = static_cast<double>(v.real());
            const double im = static_cast<double>(v.imag());
            const double p = re * re + im * im;
            if (std::isfinite(p) && p > 0.0) {
                weights[y] += p;
                total += p;
            }
        }
    }
 
    if (!(total > 0.0) || !std::isfinite(total)) {
        return result;
    }
 
    std::discrete_distribution<std::size_t> distribution(weights.begin(), weights.end());
    result.reserve(shots);
    for (std::size_t shot = 0; shot < shots; ++shot) {
        const std::size_t yIndex = distribution(engine);
        PositionMeasurement2DPoint point;
        point.x = params.detectorX;
        point.y = static_cast<float>((static_cast<double>(yIndex) - cy) * dy);
        result.push_back(point);
    }
 
    return result;
}
 
} // namespace
 
std::vector<PositionMeasurement2DPoint> DetectorScreenMeasurement2D::sample(
    const ComplexField2D &field,
    const DetectorScreenMeasurement2DParams &params,
    std::size_t shots)
{
    std::random_device rd;
    std::seed_seq seed{rd(), rd(), rd(), rd(),
                       static_cast<unsigned int>(field.width()),
                       static_cast<unsigned int>(field.height()),
                       static_cast<unsigned int>(shots)};
    std::mt19937 engine(seed);
    return sampleWithEngine(field, params, shots, engine);
}
 
std::vector<PositionMeasurement2DPoint> DetectorScreenMeasurement2D::sampleSeeded(
    const ComplexField2D &field,
    const DetectorScreenMeasurement2DParams &params,
    std::size_t shots,
    std::uint32_t seed)
{
    std::mt19937 engine(seed);
    return sampleWithEngine(field, params, shots, engine);
}
 
std::vector<std::size_t> DetectorScreenMeasurement2D::histogramY(
    const std::vector<PositionMeasurement2DPoint> &hits,
    float yMin,
    float yMax,
    std::size_t bins)
{
    std::vector<std::size_t> histogram;
    if (!(yMax > yMin) || bins == 0u) {
        return histogram;
    }
    histogram.assign(bins, 0u);
    const double span = static_cast<double>(yMax - yMin);
    for (std::size_t i = 0; i < hits.size(); ++i) {
        const double u = (static_cast<double>(hits[i].y) - static_cast<double>(yMin)) / span;
        if (u < 0.0 || u > 1.0) {
            continue;
        }
        std::size_t index = static_cast<std::size_t>(u * static_cast<double>(bins));
        if (index >= bins) {
            index = bins - 1u;
        }
        ++histogram[index];
    }
    return histogram;
}
 
} // namespace qpe