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