msc/detail/UrbanMscSafetyStepLimit.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2021-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/UrbanMscSafetyStepLimit.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/grid/PolyEvaluator.hh"
0019 #include "celeritas/phys/Interaction.hh"
0020 #include "celeritas/phys/ParticleTrackView.hh"
0021 #include "celeritas/phys/PhysicsTrackView.hh"
0022 #include "celeritas/random/distribution/NormalDistribution.hh"
0023
0024 #include "UrbanMscHelper.hh"
0025
0026 namespace celeritas
0027 {
0028 namespace detail
0029 {
0030 //---------------------------------------------------------------------------//
0031 /*!
0032  * Sample a step limit for the Urban MSC model using the "safety" algorithm.
0033  *
0034  * This distribution is to be used for tracks that have non-negligible steps
0035  * and are near the boundary. Otherwise, no displacement or step limiting is
0036  * needed.
0037  *
0038  * \note This code performs the same method as in ComputeTruePathLengthLimit
0039  * of G4UrbanMscModel, as documented in section 8.1.6 of the Geant4 10.7
0040  * Physics Reference Manual or CERN-OPEN-2006-077 by L. Urban.
0041  */
0042 class UrbanMscSafetyStepLimit
0043 {
0044   public:
0045     //!@{
0046     //! \name Type aliases
0047     using Energy = units::MevEnergy;
0048     using UrbanMscRef = NativeCRef<UrbanMscData>;
0049     //!@}
0050
0051   public:
0052     // Construct with shared and state data
0053     inline CELER_FUNCTION UrbanMscSafetyStepLimit(UrbanMscRef const& shared,
0054                                                   UrbanMscHelper const& helper,
0055                                                   Energy inc_energy,
0056                                                   PhysicsTrackView* physics,
0057                                                   MaterialId matid,
0058                                                   bool on_boundary,
0059                                                   real_type safety,
0060                                                   real_type phys_step);
0061
0062     // Apply the step limitation algorithm for the e-/e+ MSC with the RNG
0063     template<class Engine>
0064     inline CELER_FUNCTION real_type operator()(Engine& rng);
0065
0066   private:
0067     //// DATA ////
0068
0069     // Shared constant data
0070     UrbanMscRef const& shared_;
0071     // Urban MSC helper class
0072     UrbanMscHelper const& helper_;
0073
0074     // Physical step limitation up to this point
0075     real_type max_step_{};
0076     // Cached approximation for the minimum step length
0077     real_type limit_min_{};
0078     // Limit based on the range and safety
0079     real_type limit_{};
0080
0081     //// COMMON PROPERTIES ////
0082
0083     //! Minimum range for an empirical step-function approach
0084     static CELER_CONSTEXPR_FUNCTION real_type min_range()
0085     {
0086         return 1e-3 * units::centimeter;
0087     }
0088
0089     //! Maximum step over the range
0090     static CELER_CONSTEXPR_FUNCTION real_type max_step_over_range()
0091     {
0092         return 0.35;
0093     }
0094
0095     //// HELPER FUNCTIONS ////
0096
0097     // Calculate the minimum of the true path length limit
0098     inline CELER_FUNCTION real_type calc_limit_min(UrbanMscMaterialData const&,
0099                                                    Energy const) const;
0100 };
0101
0102 //---------------------------------------------------------------------------//
0103 // INLINE DEFINITIONS
0104 //---------------------------------------------------------------------------//
0105 /*!
0106  * Construct with shared and state data.
0107  */
0108 CELER_FUNCTION
0109 UrbanMscSafetyStepLimit::UrbanMscSafetyStepLimit(UrbanMscRef const& shared,
0110                                                  UrbanMscHelper const& helper,
0111                                                  Energy inc_energy,
0112                                                  PhysicsTrackView* physics,
0113                                                  MaterialId matid,
0114                                                  bool on_boundary,
0115                                                  real_type safety,
0116                                                  real_type phys_step)
0117     : shared_(shared), helper_(helper), max_step_(phys_step)
0118 {
0119     CELER_EXPECT(safety >= 0);
0120     CELER_EXPECT(safety < helper_.max_step());
0121     CELER_EXPECT(max_step_ > shared_.params.limit_min_fix());
0122     CELER_EXPECT(max_step_ <= physics->dedx_range());
0123
0124     bool use_safety_plus = physics->scalars().step_limit_algorithm
0125                            == MscStepLimitAlgorithm::safety_plus;
0126     real_type const range = physics->dedx_range();
0127     auto const& msc_range = physics->msc_range();
0128
0129     if (!msc_range || on_boundary)
0130     {
0131         MscRange new_range;
0132         // Initialize MSC range cache on the first step in a volume
0133         // TODO for hadrons/muons: this value is hard-coded for electrons
0134         new_range.range_factor = physics->scalars().range_factor;
0135         // XXX the 1 MFP limitation is applied to the *geo* step, not the true
0136         // step, so this isn't quite right (See UrbanMsc.hh)
0137         new_range.range_init = use_safety_plus
0138                                    ? range
0139                                    : max<real_type>(range, helper_.msc_mfp());
0140         if (helper_.msc_mfp() > physics->scalars().lambda_limit)
0141         {
0142             real_type c = use_safety_plus ? 0.84 : 0.75;
0143             new_range.range_factor *= c
0144                                       + (1 - c) * helper_.msc_mfp()
0145                                             / physics->scalars().lambda_limit;
0146         }
0147         new_range.limit_min
0148             = this->calc_limit_min(shared_.material_data[matid], inc_energy);
0149
0150         // Store persistent range properties within this tracking volume
0151         physics->msc_range(new_range);
0152         // Range is a reference so should be updated
0153         CELER_ASSERT(msc_range);
0154     }
0155     limit_min_ = msc_range.limit_min;
0156
0157     limit_ = range;
0158     if (safety < range)
0159     {
0160         limit_ = max<real_type>(msc_range.range_factor * msc_range.range_init,
0161                                 physics->scalars().safety_factor * safety);
0162     }
0163     limit_ = max<real_type>(limit_, limit_min_);
0164
0165     if (use_safety_plus)
0166     {
0167         real_type rho = UrbanMscSafetyStepLimit::min_range();
0168         if (range > rho)
0169         {
0170             // Calculate the scaled step range \f$ s = \alpha r + \rho (1 -
0171             // \alpha) (2 - \frac{\rho}{r}) \f$, where \f$ \alpha \f$ is the
0172             // maximum step over the range and \f$ \rho \f$ is the minimum
0173             // range
0174             real_type alpha = UrbanMscSafetyStepLimit::max_step_over_range();
0175             real_type limit_step = alpha * range
0176                                    + rho * (1 - alpha) * (2 - rho / range);
0177             max_step_ = min(max_step_, limit_step);
0178         }
0179     }
0180 }
0181
0182 //---------------------------------------------------------------------------//
0183 /*!
0184  * Sample the true path length using the Urban multiple scattering model.
0185  */
0186 template<class Engine>
0187 CELER_FUNCTION real_type UrbanMscSafetyStepLimit::operator()(Engine& rng)
0188 {
0189     if (max_step_ <= limit_)
0190     {
0191         // Skip sampling if the physics step is limiting
0192         return max_step_;
0193     }
0194     if (limit_ == limit_min_)
0195     {
0196         // Skip sampling below the minimum step limit
0197         return limit_min_;
0198     }
0199
0200     // Randomize the limit if this step should be determined by msc
0201     NormalDistribution<real_type> sample_gauss(
0202         limit_, real_type(0.1) * (limit_ - limit_min_));
0203     real_type sampled_limit = sample_gauss(rng);
0204
0205     // Keep sampled limit between the minimum value and maximum step
0206     return clamp(sampled_limit, limit_min_, max_step_);
0207 }
0208
0209 //---------------------------------------------------------------------------//
0210 /*!
0211  * Calculate the minimum of the true path length limit.
0212  *
0213  * See \c G4UrbanMscModel::ComputeTlimitmin .
0214  */
0215 CELER_FUNCTION real_type UrbanMscSafetyStepLimit::calc_limit_min(
0216     UrbanMscMaterialData const& msc, Energy const inc_energy) const
0217 {
0218     using PolyQuad = PolyEvaluator<real_type, 2>;
0219
0220     // Calculate minimum step
0221     PolyQuad calc_min_mfp(2, msc.stepmin_coeff[0], msc.stepmin_coeff[1]);
0222     real_type xm = helper_.msc_mfp() / calc_min_mfp(inc_energy.value());
0223
0224     // Scale based on particle type and effective atomic number
0225     xm *= helper_.scaled_zeff();
0226
0227     if (inc_energy < shared_.params.min_scaling_energy())
0228     {
0229         // Energy is below a pre-defined limit
0230         xm *= (real_type(0.5)
0231                + real_type(0.5) * value_as<Energy>(inc_energy)
0232                      / value_as<Energy>(shared_.params.min_scaling_energy()));
0233     }
0234
0235     return max<real_type>(xm, shared_.params.limit_min_fix());
0236 }
0237
0238 //---------------------------------------------------------------------------//
0239 }  // namespace detail
0240 }  // namespace celeritas