msc/detail/UrbanMscMinimalStepLimit.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2024 UT-Battelle, LLC, and other Celeritas developers.
0003 // See the top-level COPYRIGHT file for details.
0004 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0005 //---------------------------------------------------------------------------//
0006 //! \file celeritas/em/msc/detail/UrbanMscMinimalStepLimit.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include <cmath>
0011
0012 #include "corecel/Macros.hh"
0013 #include "corecel/Types.hh"
0014 #include "corecel/math/Algorithms.hh"
0015 #include "celeritas/Quantities.hh"
0016 #include "celeritas/Types.hh"
0017 #include "celeritas/em/data/UrbanMscData.hh"
0018 #include "celeritas/phys/PhysicsTrackView.hh"
0019 #include "celeritas/random/distribution/NormalDistribution.hh"
0020
0021 #include "UrbanMscHelper.hh"
0022
0023 namespace celeritas
0024 {
0025 namespace detail
0026 {
0027 //---------------------------------------------------------------------------//
0028 /*!
0029  * Sample a step limit for the Urban MSC model using the "minimal" algorithm.
0030  *
0031  * \note This code performs the same method as in ComputeTruePathLengthLimit
0032  * of G4UrbanMscModel, as documented in section 8.1.6 of the Geant4 10.7
0033  * Physics Reference Manual or CERN-OPEN-2006-077 by L. Urban.
0034  */
0035 class UrbanMscMinimalStepLimit
0036 {
0037   public:
0038     // Construct with shared and state data
0039     inline CELER_FUNCTION
0040     UrbanMscMinimalStepLimit(NativeCRef<UrbanMscData> const& shared,
0041                              UrbanMscHelper const& helper,
0042                              PhysicsTrackView* physics,
0043                              bool on_boundary,
0044                              real_type phys_step);
0045
0046     // Apply the step limitation algorithm for e-/e+ MSC
0047     template<class Engine>
0048     inline CELER_FUNCTION real_type operator()(Engine& rng);
0049
0050   private:
0051     //// DATA ////
0052
0053     // Physical step limitation up to this point
0054     real_type max_step_{};
0055     // Cached approximation for the minimum step length
0056     real_type limit_min_{};
0057     // Limit based on the range
0058     real_type limit_{};
0059 };
0060
0061 //---------------------------------------------------------------------------//
0062 // INLINE DEFINITIONS
0063 //---------------------------------------------------------------------------//
0064 /*!
0065  * Construct with shared and state data.
0066  */
0067 CELER_FUNCTION
0068 UrbanMscMinimalStepLimit::UrbanMscMinimalStepLimit(
0069     NativeCRef<UrbanMscData> const& shared,
0070     UrbanMscHelper const& helper,
0071     PhysicsTrackView* physics,
0072     bool on_boundary,
0073     real_type phys_step)
0074     : max_step_(phys_step)
0075 {
0076     CELER_EXPECT(max_step_ > shared.params.limit_min_fix());
0077     CELER_EXPECT(max_step_ <= physics->dedx_range());
0078
0079     auto const& msc_range = physics->msc_range();
0080
0081     if (!msc_range)
0082     {
0083         // Store initial range properties if this is the track's first step
0084         MscRange new_range;
0085         new_range.range_init = numeric_limits<real_type>::infinity();
0086         new_range.range_factor = physics->scalars().range_factor;
0087         new_range.limit_min = 10 * shared.params.limit_min_fix();
0088         physics->msc_range(new_range);
0089         CELER_ASSERT(msc_range);
0090     }
0091     limit_min_ = msc_range.limit_min;
0092
0093     if (on_boundary)
0094     {
0095         // Update the MSC range for the new volume
0096         MscRange new_range = msc_range;
0097         new_range.range_init = msc_range.range_factor
0098                                * max(physics->dedx_range(), helper.msc_mfp());
0099         new_range.range_init = max(new_range.range_init, limit_min_);
0100         physics->msc_range(new_range);
0101         CELER_ASSERT(msc_range);
0102     }
0103     limit_ = msc_range.range_init;
0104 }
0105
0106 //---------------------------------------------------------------------------//
0107 /*!
0108  * Sample the true path length using the Urban multiple scattering model.
0109  */
0110 template<class Engine>
0111 CELER_FUNCTION real_type UrbanMscMinimalStepLimit::operator()(Engine& rng)
0112 {
0113     if (max_step_ <= limit_)
0114     {
0115         // Skip sampling if the physics step is limiting
0116         return max_step_;
0117     }
0118     if (limit_ == limit_min_)
0119     {
0120         // Skip sampling below the minimum step limit
0121         return limit_min_;
0122     }
0123
0124     // Randomize the limit if this step should be determined by MSC
0125     NormalDistribution<real_type> sample_gauss(
0126         limit_, real_type(0.1) * (limit_ - limit_min_));
0127     real_type sampled_limit = sample_gauss(rng);
0128
0129     // Keep sampled limit between the minimum value and maximum step
0130     return clamp(sampled_limit, limit_min_, max_step_);
0131 }
0132
0133 //---------------------------------------------------------------------------//
0134 }  // namespace detail
0135 }  // namespace celeritas