Generalize Poisson distribution to allow zeros

doc/CMakeLists.txt

@@ -288,6 +288,8 @@ file(GLOB _DOXYGEN_SOURCE
   "${PROJECT_SOURCE_DIR}/src/celeritas/optical/surface/*.hh"
   "${PROJECT_SOURCE_DIR}/src/celeritas/optical/surface/model/*.hh"
   "${PROJECT_SOURCE_DIR}/test/*.hh"
+  "${PROJECT_SOURCE_DIR}/test/corecel/*.hh"
+  "${PROJECT_SOURCE_DIR}/test/corecel/random/*.hh"
 )
 
 # IMPORTANT: the target name and USE_STAMP_FILE generate a dependency used by

doc/development/testing.rst

@@ -20,48 +20,72 @@ functions to better isolate conditionals from each other.
 
 .. _cyclomatic complexity: https://en.wikipedia.org/wiki/Cyclomatic_complexity
 
-Running CTest
--------------
-
-When configured with ``CELERITAS_BUILD_TESTS`` (see :ref:`configuration`),
-CTest_ will be automatically configured. Running through CTest sets special
-environment variables for data, disabling GPUs, or testing the code itself.
-CTest can be run either through ``ninja test`` or by manually invoking
-``ctest``.  Two useful ways to run are ``ctest -V -R <test regex>``, which will
-run one or more tests that match the given regular expression (such as
-``corecel/math/``), and ``ctest -j --output-on-failure`` which runs in parallel
-and prints only test failures.
-
-.. _CTest: https://cmake.org/cmake/help/latest/manual/ctest.1.html
-
-Using GoogleTest
-----------------
+Writing unit tests
+------------------
 
 `Google test`_ is very well documented, and because so much testing code exists
 in AI training data, tools like ChatGPT and Copilot are very good at writing a
 first pass at test code.
+
 Celeritas defines a base class test harness with some utility functions:
 
 .. doxygenclass:: celeritas::test::Test
 
-as well as several macros that simplify testing with floating-point data (and
-vectors thereof):
+Testing real numbers and JSON data
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Celeritas defines special test macros that simplify testing with floating-point data (and vectors thereof):
 
 .. doxygendefine:: EXPECT_VEC_EQ
 .. doxygendefine:: EXPECT_REAL_EQ
 .. doxygendefine:: EXPECT_SOFT_EQ
 .. doxygendefine:: EXPECT_SOFT_NEAR
 .. doxygendefine:: EXPECT_VEC_SOFT_EQ
 .. doxygendefine:: EXPECT_VEC_NEAR
+.. doxygendefine:: EXPECT_REF_EQ
 .. doxygendefine:: EXPECT_JSON_EQ
+
+Expected data can be printed with another macro:
+
 .. doxygendefine:: PRINT_EXPECTED
 
 For more details on the test harnesses, especially the hierarchy used for
 setting up physics problems for testing, see the ``celeritas::test`` namespace
 in the Doxygen developer documentation.
 
-You can run most tests manually from the build directory and filter so that
-only a subset of tests run::
+Testing random number distributions
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. doxygenclass:: celeritas::test::DiagnosticRngEngine
+
+.. doxygenclass:: celeritas::test::Histogram
+
+.. doxygenclass:: celeritas::test::HistogramSampler
+
+Testing strings and logger output
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. doxygenclass:: celeritas::test::ScopedLogStorer
+
+.. doxygenclass:: celeritas::test::StringSimplifier
+
+
+Running unit tests
+------------------
+
+When configured with ``CELERITAS_BUILD_TESTS`` (see :ref:`configuration`),
+CTest_ will be automatically configured. Running through CTest sets special
+environment variables for data, disabling GPUs, or testing the code itself.
+CTest can be run either through ``ninja test`` or by manually invoking
+``ctest``.  Two useful ways to run are ``ctest -V -R <test regex>``, which will
+run one or more tests that match the given regular expression (such as
+``corecel/math/``), and ``ctest -j --output-on-failure`` which runs in parallel
+and prints only test failures.
+
+.. _CTest: https://cmake.org/cmake/help/latest/manual/ctest.1.html
+
+You can *also* run most tests manually from the build directory and filter so
+that only a subset of tests run::
 
    $ ./test/celeritas/global_Stepper --gtest_filter=SimpleComptonTest.host
 

src/celeritas/em/distribution/EnergyLossGaussianDistribution.hh

@@ -11,7 +11,6 @@
 #include "corecel/Assert.hh"
 #include "corecel/Macros.hh"
 #include "corecel/Types.hh"
-#include "corecel/math/Algorithms.hh"
 #include "corecel/random/distribution/NormalDistribution.hh"
 #include "corecel/random/distribution/TruncatedDistribution.hh"
 

src/celeritas/em/distribution/EnergyLossUrbanDistribution.hh

@@ -90,7 +90,10 @@ class EnergyLossUrbanDistribution
     static CELER_CONSTEXPR_FUNCTION real_type rate() { return 0.56; }
 
     //! Number of collisions above which to use faster sampling from Gaussian
-    static CELER_CONSTEXPR_FUNCTION size_type max_collisions() { return 8; }
+    static CELER_CONSTEXPR_FUNCTION size_type collision_poisson_threshold()
+    {
+        return 8;
+    }
 
     //! Threshold number of excitations used in width correction
     static CELER_CONSTEXPR_FUNCTION real_type exc_thresh() { return 42; }
@@ -263,7 +266,7 @@ EnergyLossUrbanDistribution::sample_excitation_loss(Engine& rng)
 
     for (int i : range(2))
     {
-        if (xs_exc_[i] > this->max_collisions())
+        if (xs_exc_[i] > this->collision_poisson_threshold())
         {
             // When the number of collisions is large, use faster approach
             // of sampling from a Gaussian
@@ -275,12 +278,12 @@ EnergyLossUrbanDistribution::sample_excitation_loss(Engine& rng)
             // The loss due to excitation is \f$ \Delta E_{exc} = n_1 E_1 + n_2
             // E_2 \f$, where the number of collisions \f$ n_i \f$ is sampled
             // from a Poisson distribution with mean \f$ \Sigma_i \f$
-            unsigned int n = PoissonDistribution<real_type>(xs_exc_[i])(rng);
+            auto n = PoissonDistributionKnuth<real_type>(xs_exc_[i])(rng);
             if (n > 0)
             {
-                UniformRealDistribution<real_type> sample_fraction(n - 1,
-                                                                   n + 1);
-                result += sample_fraction(rng) * binding_energy_[i];
+                UniformRealDistribution<real_type> sample_smooth(-1, 1);
+                result += (static_cast<real_type>(n) + sample_smooth(rng))
+                          * binding_energy_[i];
             }
         }
     }
@@ -312,15 +315,16 @@ EnergyLossUrbanDistribution::sample_ionization_loss(Engine& rng)
     // Mean number of collisions in the fast simulation interval
     real_type mean_num_coll = 0;
 
-    if (xs_ion_ > this->max_collisions())
+    if (xs_ion_ > this->collision_poisson_threshold())
     {
         // When the number of collisions is large, fast sampling from a
         // Gaussian is used in the lower portion of the energy loss interval.
         // See PHYS332 section 2.4: Fast simulation for \f$ n_3 \ge 16 \f$
 
         // Calculate the maximum value of \f$ \alpha \f$ (Eq. 25)
-        alpha = (xs_ion_ + this->max_collisions()) * energy_ratio
-                / (this->max_collisions() * energy_ratio + xs_ion_);
+        alpha
+            = (xs_ion_ + this->collision_poisson_threshold()) * energy_ratio
+              / (this->collision_poisson_threshold() * energy_ratio + xs_ion_);
 
         // Mean energy loss for a single collision of this type (Eq. 14)
         real_type const mean_loss_coll = alpha * std::log(alpha) / (alpha - 1);

src/celeritas/em/msc/detail/UrbanMscMinimalStepLimit.hh

@@ -6,13 +6,10 @@
 //---------------------------------------------------------------------------//
 #pragma once
 
-#include <cmath>
-
 #include "corecel/Macros.hh"
 #include "corecel/Types.hh"
 #include "corecel/math/Algorithms.hh"
 #include "corecel/random/distribution/NormalDistribution.hh"
-#include "celeritas/Quantities.hh"
 #include "celeritas/Types.hh"
 #include "celeritas/em/data/UrbanMscData.hh"
 #include "celeritas/phys/PhysicsTrackView.hh"

src/celeritas/em/msc/detail/UrbanMscScatter.hh

@@ -526,7 +526,7 @@ UrbanMscScatter::compute_theta0(ParticleTrackView const& particle) const
     constexpr units::MevEnergy c_highland{13.6};
     real_type theta0
         = c_highland.value()
-          * std::abs(value_as<units::ElementaryCharge>(particle.charge()))
+          * std::fabs(value_as<units::ElementaryCharge>(particle.charge()))
           * std::sqrt(y) * invbetacp;
 
     // Correction factor from e- scattering data

src/celeritas/ext/detail/GeantProcessImporter.cc

@@ -387,14 +387,15 @@ inp::UniformGrid import_physics_log_vector(G4PhysicsVector const& pv,
     double const x_scaling = native_value_from_clhep(units[0]);
     double const y_scaling = native_value_from_clhep(units[1]);
     auto size = pv.GetVectorLength();
+    CELER_ASSERT(size > 1);
 
     inp::UniformGrid grid;
     grid.x = {std::log(pv.Energy(0) * x_scaling),
               std::log(pv.Energy(size - 1) * x_scaling)};
     grid.y.resize(size);
 
-    double delta
-        = fastpow(pv.Energy(size - 1) / pv.Energy(0), 1.0 / (size - 1));
+    double delta = fastpow(pv.Energy(size - 1) / pv.Energy(0),
+                           1.0 / static_cast<double>(size - 1));
     for (auto i : range(size))
     {
         // Check that the grid has log spacing

src/celeritas/optical/gen/CherenkovOffload.hh

@@ -69,7 +69,7 @@ class CherenkovOffload
     real_type step_length_;
     OffloadPreStepData const& pre_step_;
     optical::GeneratorStepData post_step_;
-    real_type num_photons_per_len_;
+    PoissonDistribution<real_type> sample_num_photons_{0};
 };
 
 //---------------------------------------------------------------------------//
@@ -98,11 +98,12 @@ CherenkovOffload::CherenkovOffload(units::ElementaryCharge charge,
 
     using namespace celeritas::literals;
 
-    units::LightSpeed beta(
+    // NOTE: this uses the average beta, not the average number of photons
+    units::LightSpeed avg_beta(
         0.5_r * (pre_step_.speed.value() + post_step_.speed.value()));
-
     CherenkovDndxCalculator calculate_dndx(pre_mat, shared, charge_);
-    num_photons_per_len_ = calculate_dndx(beta);
+    sample_num_photons_
+        = PoissonDistribution{calculate_dndx(avg_beta) * step_length_};
 }
 
 //---------------------------------------------------------------------------//
@@ -137,14 +138,8 @@ template<class Generator>
 CELER_FUNCTION optical::GeneratorDistributionData
 CherenkovOffload::operator()(Generator& rng)
 {
-    if (num_photons_per_len_ == 0)
-    {
-        return {};
-    }
-
     optical::GeneratorDistributionData data;
-    data.num_photons = PoissonDistribution<real_type>(num_photons_per_len_
-                                                      * step_length_)(rng);
+    data.num_photons = sample_num_photons_(rng);
     if (data.num_photons > 0)
     {
         data.type = GeneratorType::cherenkov;

src/celeritas/optical/gen/ScintillationData.hh

@@ -48,9 +48,10 @@ struct ScintDistributionRecord
 /*!
  * Unnormalized scintillation spectrum as a sum of independent components.
  *
- * \todo The yield and resolution scale should live together (used for sampling
- * the number of photons) and be separated from the normalized spectrum (used
- * during generation, represented as a sum of components).
+ * \todo The yield and resolution scale should live together (used by \c
+ * ScintillationOffload for sampling the number of photons) and be separated
+ * from the normalized spectrum (used during generation by \c
+ * ScintillationGenerator, represented as a sum of components).
  *
  * - \c yield_per_energy is the average number of photons released by a unit of
  *   locally deposited energy.
@@ -64,11 +65,10 @@ struct ScintSpectrumRecord
     ItemRange<real_type> yield_pdf;
     ItemRange<ScintDistributionRecord> components;
 
-    //! Whether all data are assigned and valid
+    //! Whether scintillation is present for this material
     explicit CELER_FUNCTION operator bool() const
     {
-        return yield_per_energy > 0 && !yield_pdf.empty()
-               && yield_pdf.size() == components.size();
+        return yield_per_energy > 0;
     }
 };
 

src/celeritas/optical/gen/ScintillationOffload.hh

@@ -10,6 +10,7 @@
 #include "corecel/Macros.hh"
 #include "corecel/random/distribution/NormalDistribution.hh"
 #include "corecel/random/distribution/PoissonDistribution.hh"
+#include "corecel/random/distribution/RoundedNonnegDistribution.hh"
 #include "celeritas/Quantities.hh"
 #include "celeritas/Types.hh"
 #include "celeritas/phys/ParticleTrackView.hh"
@@ -63,17 +64,15 @@ class ScintillationOffload
 
   private:
     units::ElementaryCharge charge_;
-    real_type step_length_;
+    real_type step_length_{};
     OffloadPreStepData const& pre_step_;
     optical::GeneratorStepData post_step_;
     NativeCRef<ScintillationData> const& shared_;
-    real_type continuous_edep_fraction_;
+    real_type continuous_edep_fraction_{};
     real_type mean_num_photons_{0};
 
-    static CELER_CONSTEXPR_FUNCTION real_type poisson_threshold()
-    {
-        return 10;
-    }
+    // Use scaled gaussian above this average lambda
+    static constexpr real_type poisson_threshold = 10;
 };
 
 //---------------------------------------------------------------------------//
@@ -111,6 +110,7 @@ CELER_FUNCTION ScintillationOffload::ScintillationOffload(
     //! \todo Use visible energy deposition when Birks law is implemented
     if (spectrum)
     {
+        // This material is a scintillator
         mean_num_photons_ = spectrum.yield_per_energy
                             * energy_deposition.value();
     }
@@ -128,20 +128,24 @@ ScintillationOffload::operator()(Generator& rng)
 {
     // Material-only sampling
     optical::GeneratorDistributionData result;
-    if (mean_num_photons_ > poisson_threshold())
+    if (mean_num_photons_ > poisson_threshold)
     {
         using namespace celeritas::literals;
 
         real_type sigma = shared_.resolution_scale[pre_step_.material]
                           * std::sqrt(mean_num_photons_);
-        result.num_photons = static_cast<size_type>(clamp_to_nonneg(
-            NormalDistribution<real_type>(mean_num_photons_, sigma)(rng)
-            + 0.5_r));
+        result.num_photons
+            = RoundedNonnegDistribution<NormalDistribution<real_type>>(
+                mean_num_photons_, sigma)(rng);
     }
     else if (mean_num_photons_ > 0)
     {
-        result.num_photons = static_cast<size_type>(
-            PoissonDistribution<real_type>(mean_num_photons_)(rng));
+        result.num_photons
+            = PoissonDistributionKnuth<real_type>(mean_num_photons_)(rng);
+    }
+    else
+    {
+        result.num_photons = 0;
     }
 
     if (result.num_photons > 0)

src/celeritas/optical/gen/detail/ScintSpectrumInserter.hh

@@ -139,6 +139,9 @@ auto ScintSpectrumInserter::operator()(inp::ScintillationMaterial const& mat)
     spectrum.yield_per_energy = total_yield;
     spectrum.components = {begin_components, scint_records_.size_id()};
     spectrum.yield_pdf = reals_.insert_back(yield_pdf.begin(), yield_pdf.end());
+    CELER_ASSERT(total_yield > 0);
+    CELER_ASSERT(!spectrum.components.empty());
+    CELER_ASSERT(spectrum.components.size() == spectrum.yield_pdf.size());
 
     CELER_ENSURE(spectrum.components.size() == mat.components.size());
     return spectra_.push_back(std::move(spectrum));