detectortemporalflux2d.cpp

Go to the documentation of this file.

/**
 * @file detectortemporalflux2d.cpp
 * @brief Implements the detectortemporalflux2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/detectortemporalflux2d.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 "detectortemporalflux2d.h"
 
#include "core/math/complexfield2d.h"
 
#include <algorithm>
#include <cmath>
#include <complex>
#include <random>
#include <stdexcept>
 
namespace qpe {
namespace {
 
std::vector<std::size_t> detectorColumns(const ComplexField2D &field,
                                         float detectorX,
                                         float thickness)
{
    std::vector<std::size_t> columns;
    if (field.empty()) {
        return columns;
    }
 
    const double cx = 0.5 * static_cast<double>(field.width() - 1u);
    const double dx = static_cast<double>(field.dx());
    const double halfThickness = 0.5 * std::max<double>(static_cast<double>(thickness), dx);
 
    for (std::size_t x = 1u; x + 1u < field.width(); ++x) {
        const double physicalX = (static_cast<double>(x) - cx) * dx;
        if (std::fabs(physicalX - static_cast<double>(detectorX)) <= halfThickness) {
            columns.push_back(x);
        }
    }
 
    if (columns.empty() && field.width() >= 3u) {
        long nearest = static_cast<long>(std::lround(static_cast<double>(detectorX) / dx + cx));
        nearest = std::max<long>(1, std::min<long>(static_cast<long>(field.width()) - 2, nearest));
        columns.push_back(static_cast<std::size_t>(nearest));
    }
    return columns;
}
 
} // namespace
 
void DetectorTemporalFlux2D::clear()
{
    profile_.clear();
    height_ = 0u;
    dy_ = 0.0f;
    detectorX_ = 0.0f;
    totalPositiveCrossingProbability_ = 0.0;
}
 
void DetectorTemporalFlux2D::accumulate(const ComplexField2D &field,
                                         float hbar,
                                         float mass,
                                         float dt,
                                         const DetectorTemporalFlux2DParams &params)
{
    if (field.empty()) {
        return;
    }
    if (!(hbar > 0.0f) || !std::isfinite(static_cast<double>(hbar))) {
        throw std::invalid_argument("DetectorTemporalFlux2D: hbar must be finite and positive");
    }
    if (!(mass > 0.0f) || !std::isfinite(static_cast<double>(mass))) {
        throw std::invalid_argument("DetectorTemporalFlux2D: mass must be finite and positive");
    }
    if (!(dt > 0.0f) || !std::isfinite(static_cast<double>(dt))) {
        throw std::invalid_argument("DetectorTemporalFlux2D: dt must be finite and positive");
    }
 
    if (height_ != field.height() || dy_ != field.dy() ||
        std::fabs(static_cast<double>(detectorX_ - params.detectorX)) > 1.0e-6) {
        profile_.assign(field.height(), 0.0);
        height_ = field.height();
        dy_ = field.dy();
        detectorX_ = params.detectorX;
        totalPositiveCrossingProbability_ = 0.0;
    }
 
    const std::vector<std::size_t> columns =
        detectorColumns(field, params.detectorX, params.thickness);
    if (columns.empty() || field.height() < 3u) {
        return;
    }
 
    const double inv2Dx = 1.0 / (2.0 * static_cast<double>(field.dx()));
    const double prefactor = static_cast<double>(hbar) / static_cast<double>(mass);
    const double stepDt = static_cast<double>(dt);
    const double invColumns = 1.0 / static_cast<double>(columns.size());
    double stepProbability = 0.0;
 
    for (std::size_t y = 1u; y + 1u < field.height(); ++y) {
        double currentAverage = 0.0;
        for (std::size_t i = 0u; i < columns.size(); ++i) {
            const std::size_t x = columns[i];
            const std::complex<double> psi(
                static_cast<double>(field(x, y).real()),
                static_cast<double>(field(x, y).imag()));
            const std::complex<double> psiXp(
                static_cast<double>(field(x + 1u, y).real()),
                static_cast<double>(field(x + 1u, y).imag()));
            const std::complex<double> psiXm(
                static_cast<double>(field(x - 1u, y).real()),
                static_cast<double>(field(x - 1u, y).imag()));
            const std::complex<double> dPsiDx = (psiXp - psiXm) * inv2Dx;
            currentAverage += prefactor * std::imag(std::conj(psi) * dPsiDx);
        }
        currentAverage *= invColumns;
        const double positiveCurrent = std::max(0.0, currentAverage);
        if (positiveCurrent > 0.0 && std::isfinite(positiveCurrent)) {
            const double weightIncrement = positiveCurrent * stepDt;
            profile_[y] += weightIncrement;
            stepProbability += positiveCurrent * static_cast<double>(field.dy()) * stepDt;
        }
    }
 
    totalPositiveCrossingProbability_ += stepProbability;
}
 
std::vector<PositionMeasurement2DPoint> DetectorTemporalFlux2D::sample(std::size_t shots) const
{
    return sampleWithSeed(shots, 0u, false);
}
 
std::vector<PositionMeasurement2DPoint> DetectorTemporalFlux2D::sampleSeeded(std::size_t shots,
                                                                              std::uint32_t seed) const
{
    return sampleWithSeed(shots, seed, true);
}
 
std::vector<PositionMeasurement2DPoint> DetectorTemporalFlux2D::sampleWithSeed(
    std::size_t shots,
    std::uint32_t seed,
    bool deterministic) const
{
    std::vector<PositionMeasurement2DPoint> result;
    if (shots == 0u || empty()) {
        return result;
    }
 
    double sum = 0.0;
    for (std::size_t i = 0u; i < profile_.size(); ++i) {
        if (std::isfinite(profile_[i]) && profile_[i] > 0.0) {
            sum += profile_[i];
        }
    }
    if (!(sum > 0.0)) {
        return result;
    }
 
    std::mt19937 engine;
    if (deterministic) {
        engine.seed(seed);
    } else {
        std::random_device rd;
        std::seed_seq sequence{rd(), rd(), rd(), rd(),
                               static_cast<unsigned int>(shots),
                               static_cast<unsigned int>(profile_.size())};
        engine.seed(sequence);
    }
 
    std::discrete_distribution<std::size_t> distribution(profile_.begin(), profile_.end());
    result.reserve(shots);
    const double cy = 0.5 * static_cast<double>(height_ - 1u);
    for (std::size_t i = 0u; i < shots; ++i) {
        const std::size_t yIndex = distribution(engine);
        PositionMeasurement2DPoint point;
        point.x = detectorX_;
        point.y = static_cast<float>((static_cast<double>(yIndex) - cy) * static_cast<double>(dy_));
        result.push_back(point);
    }
    return result;
}
 
bool DetectorTemporalFlux2D::empty() const noexcept
{
    return profile_.empty() || height_ == 0u || !(dy_ > 0.0f) ||
           !(totalPositiveCrossingProbability_ > 0.0);
}
 
double DetectorTemporalFlux2D::totalPositiveCrossingProbability() const noexcept
{
    return totalPositiveCrossingProbability_;
}
 
float DetectorTemporalFlux2D::detectorX() const noexcept
{
    return detectorX_;
}
 
float DetectorTemporalFlux2D::dy() const noexcept
{
    return dy_;
}
 
std::size_t DetectorTemporalFlux2D::sampleCountY() const noexcept
{
    return height_;
}
 
const std::vector<double> &DetectorTemporalFlux2D::profile() const noexcept
{
    return profile_;
}
 
} // namespace qpe