Add "minimal" and "safety plus" Urban MSC step limit algorithms (#1146)

* Remove redundant comparison in Urban MSC step limit
* Add minimal step limit algorithm for Urban MSC
* Add safety plus step limit algorithm for Urban MSC
* Add Urban MSC step limit test
* Import MSC step limit algorithm from Geant4
* Address review feedback

src/celeritas/Types.cc

@@ -97,6 +97,21 @@ char const* to_cstring(TrackOrder value)
     return to_cstring_impl(value);
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Get a string corresponding to the MSC step limit algorithm.
+ */
+char const* to_cstring(MscStepLimitAlgorithm value)
+{
+    static EnumStringMapper<MscStepLimitAlgorithm> const to_cstring_impl{
+        "minimal",
+        "safety",
+        "safety_plus",
+        "distance_to_boundary",
+    };
+    return to_cstring_impl(value);
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Checks that the TrackOrder will sort tracks by actions applied at the given

src/celeritas/Types.hh

@@ -160,6 +160,17 @@ enum class TrackOrder
     size_
 };
 
+//---------------------------------------------------------------------------//
+//! Algorithm used to calculate the multiple scattering step limit
+enum class MscStepLimitAlgorithm
+{
+    minimal,
+    safety,
+    safety_plus,
+    distance_to_boundary,
+    size_,
+};
+
 //---------------------------------------------------------------------------//
 // HELPER STRUCTS
 //---------------------------------------------------------------------------//
@@ -196,6 +207,9 @@ char const* to_cstring(ActionOrder);
 // Get a string corresponding to a track ordering policy
 char const* to_cstring(TrackOrder);
 
+// Get a string corresponding to the MSC step limit algorithm
+char const* to_cstring(MscStepLimitAlgorithm value);
+
 // Checks that the TrackOrder will sort tracks by actions applied at the given
 // ActionOrder
 bool is_action_sorted(ActionOrder action, TrackOrder track);

src/celeritas/em/UrbanMscParams.cc

@@ -61,6 +61,7 @@ UrbanMscParams::from_import(ParticleParams const& particles,
     opts.lambda_limit = import.em_params.msc_lambda_limit;
     opts.safety_fact = import.em_params.msc_safety_factor;
     opts.range_fact = import.em_params.msc_range_factor;
+    opts.step_limit_algorithm = import.em_params.msc_step_algorithm;
 
     return std::make_shared<UrbanMscParams>(
         particles, materials, import.msc_models, opts);
@@ -111,6 +112,16 @@ UrbanMscParams::UrbanMscParams(ParticleParams const& particles,
     host_data.params.lambda_limit = options.lambda_limit;
     host_data.params.range_fact = options.range_fact;
     host_data.params.safety_fact = options.safety_fact;
+    host_data.params.step_limit_algorithm = options.step_limit_algorithm;
+    if (host_data.params.step_limit_algorithm
+        == MscStepLimitAlgorithm::distance_to_boundary)
+    {
+        CELER_LOG(warning) << "Unsupported step limit algorithm '"
+                           << to_cstring(host_data.params.step_limit_algorithm)
+                           << "': defaulting to '"
+                           << to_cstring(MscStepLimitAlgorithm::safety) << "'";
+        host_data.params.step_limit_algorithm = MscStepLimitAlgorithm::safety;
+    }
 
     // Filter MSC data by model and particle type
     std::vector<ImportMscModel const*> urban_data(particles.size(), nullptr);

src/celeritas/em/UrbanMscParams.hh

@@ -48,6 +48,8 @@ class UrbanMscParams final : public ParamsDataInterface<UrbanMscData>
         real_type lambda_limit{1 * units::millimeter};  //!< Lambda limit
         real_type range_fact{0.04};  //!< Range factor for e-/e+
         real_type safety_fact{0.6};  //!< Safety factor
+        MscStepLimitAlgorithm step_limit_algorithm{
+            MscStepLimitAlgorithm::safety};
     };
 
   public:

src/celeritas/em/data/UrbanMscData.hh

@@ -39,6 +39,7 @@ struct UrbanMscParameters
     real_type geom_limit{5e-8 * units::millimeter};  //!< minimum step
     Energy low_energy_limit{0};
     Energy high_energy_limit{0};
+    MscStepLimitAlgorithm step_limit_algorithm{MscStepLimitAlgorithm::safety};
 
     //! A scale factor for the range
     static CELER_CONSTEXPR_FUNCTION real_type dtrl() { return 5e-2; }

src/celeritas/em/msc/UrbanMsc.hh

@@ -18,8 +18,9 @@
 #include "detail/MscStepFromGeo.hh"  // IWYU pragma: associated
 #include "detail/MscStepToGeo.hh"  // IWYU pragma: associated
 #include "detail/UrbanMscHelper.hh"  // IWYU pragma: associated
+#include "detail/UrbanMscMinimalStepLimit.hh"  // IWYU pragma: associated
+#include "detail/UrbanMscSafetyStepLimit.hh"  // IWYU pragma: associated
 #include "detail/UrbanMscScatter.hh"  // IWYU pragma: associated
-#include "detail/UrbanMscStepLimit.hh"  // IWYU pragma: associated
 
 namespace celeritas
 {
@@ -130,17 +131,33 @@ CELER_FUNCTION void UrbanMsc::limit_step(CoreTrackView const& track)
             }
         }
 
-        displaced = true;
-        detail::UrbanMscStepLimit calc_limit(msc_params_,
-                                             msc_helper,
-                                             par.energy(),
-                                             &phys,
-                                             phys.material_id(),
-                                             geo.is_on_boundary(),
-                                             safety,
-                                             sim.step_length());
         auto rng = track.make_rng_engine();
+        displaced = true;
+
+        if (msc_params_.params.step_limit_algorithm
+            == MscStepLimitAlgorithm::minimal)
+        {
+            // Calculate step limit using "minimal" algorithm
+            detail::UrbanMscMinimalStepLimit calc_limit(msc_params_,
+                                                        msc_helper,
+                                                        &phys,
+                                                        geo.is_on_boundary(),
+                                                        sim.step_length());
+            return calc_limit(rng);
+        }
+
+        // Calculate step limit using "safety" or "safety plus" algorithm
+        detail::UrbanMscSafetyStepLimit calc_limit(msc_params_,
+                                                   msc_helper,
+                                                   par.energy(),
+                                                   &phys,
+                                                   phys.material_id(),
+                                                   geo.is_on_boundary(),
+                                                   safety,
+                                                   sim.step_length());
         return calc_limit(rng);
+
+        // TODO: "distance to boundary" step limit algorithm
     }();
     CELER_ASSERT(true_path <= sim.step_length());
 

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

@@ -0,0 +1,135 @@
+//----------------------------------*-C++-*----------------------------------//
+// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers.
+// See the top-level COPYRIGHT file for details.
+// SPDX-License-Identifier: (Apache-2.0 OR MIT)
+//---------------------------------------------------------------------------//
+//! \file celeritas/em/msc/detail/UrbanMscMinimalStepLimit.hh
+//---------------------------------------------------------------------------//
+#pragma once
+
+#include <cmath>
+
+#include "corecel/Macros.hh"
+#include "corecel/Types.hh"
+#include "corecel/math/Algorithms.hh"
+#include "celeritas/Quantities.hh"
+#include "celeritas/Types.hh"
+#include "celeritas/em/data/UrbanMscData.hh"
+#include "celeritas/phys/PhysicsTrackView.hh"
+#include "celeritas/random/distribution/NormalDistribution.hh"
+
+#include "UrbanMscHelper.hh"
+
+namespace celeritas
+{
+namespace detail
+{
+//---------------------------------------------------------------------------//
+/*!
+ * Sample a step limit for the Urban MSC model using the "minimal" algorithm.
+ *
+ * \note This code performs the same method as in ComputeTruePathLengthLimit
+ * of G4UrbanMscModel, as documented in section 8.1.6 of the Geant4 10.7
+ * Physics Reference Manual or CERN-OPEN-2006-077 by L. Urban.
+ */
+class UrbanMscMinimalStepLimit
+{
+  public:
+    // Construct with shared and state data
+    inline CELER_FUNCTION
+    UrbanMscMinimalStepLimit(NativeCRef<UrbanMscData> const& shared,
+                             UrbanMscHelper const& helper,
+                             PhysicsTrackView* physics,
+                             bool on_boundary,
+                             real_type phys_step);
+
+    // Apply the step limitation algorithm for e-/e+ MSC
+    template<class Engine>
+    inline CELER_FUNCTION real_type operator()(Engine& rng);
+
+  private:
+    //// DATA ////
+
+    // Physical step limitation up to this point
+    real_type max_step_{};
+    // Cached approximation for the minimum step length
+    real_type limit_min_{};
+    // Limit based on the range
+    real_type limit_{};
+};
+
+//---------------------------------------------------------------------------//
+// INLINE DEFINITIONS
+//---------------------------------------------------------------------------//
+/*!
+ * Construct with shared and state data.
+ */
+CELER_FUNCTION
+UrbanMscMinimalStepLimit::UrbanMscMinimalStepLimit(
+    NativeCRef<UrbanMscData> const& shared,
+    UrbanMscHelper const& helper,
+    PhysicsTrackView* physics,
+    bool on_boundary,
+    real_type phys_step)
+    : max_step_(phys_step)
+{
+    CELER_EXPECT(max_step_ > shared.params.limit_min_fix());
+    CELER_EXPECT(max_step_ <= physics->dedx_range());
+
+    auto const& msc_range = physics->msc_range();
+
+    if (!msc_range)
+    {
+        // Store initial range properties if this is the track's first step
+        MscRange new_range;
+        new_range.range_init = numeric_limits<real_type>::infinity();
+        new_range.range_fact = shared.params.range_fact;
+        new_range.limit_min = 10 * shared.params.limit_min_fix();
+        physics->msc_range(new_range);
+        CELER_ASSERT(msc_range);
+    }
+    limit_min_ = msc_range.limit_min;
+
+    if (on_boundary)
+    {
+        // Update the MSC range for the new volume
+        MscRange new_range = msc_range;
+        new_range.range_init = msc_range.range_fact
+                               * max(physics->dedx_range(), helper.msc_mfp());
+        new_range.range_init = max(new_range.range_init, limit_min_);
+        physics->msc_range(new_range);
+        CELER_ASSERT(msc_range);
+    }
+    limit_ = msc_range.range_init;
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Sample the true path length using the Urban multiple scattering model.
+ */
+template<class Engine>
+CELER_FUNCTION real_type UrbanMscMinimalStepLimit::operator()(Engine& rng)
+{
+    if (max_step_ <= limit_)
+    {
+        // Skip sampling if the physics step is limiting
+        return max_step_;
+    }
+    if (limit_ == limit_min_)
+    {
+        // Skip sampling below the minimum step limit
+        return limit_min_;
+    }
+
+    // Randomize the limit if this step should be determined by MSC
+    NormalDistribution<real_type> sample_gauss(
+        limit_, real_type(0.1) * (limit_ - limit_min_));
+    real_type sampled_limit = sample_gauss(rng);
+
+    // Keep sampled limit between the minimum value and maximum step
+    return clamp(sampled_limit, limit_min_, max_step_);
+}
+
+//---------------------------------------------------------------------------//
+}  // namespace detail
+}  // namespace celeritas