msc/detail/UrbanMscMinimalStepLimit.hh

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