em/msc/UrbanMsc.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2022-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/UrbanMsc.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/Assert.hh"
0011 #include "corecel/Macros.hh"
0012 #include "corecel/Types.hh"
0013 #include "corecel/math/ArrayOperators.hh"
0014 #include "celeritas/Types.hh"
0015 #include "celeritas/em/data/UrbanMscData.hh"
0016 #include "celeritas/global/CoreTrackView.hh"
0017
0018 #include "detail/MscStepFromGeo.hh"  // IWYU pragma: associated
0019 #include "detail/MscStepToGeo.hh"  // IWYU pragma: associated
0020 #include "detail/UrbanMscHelper.hh"  // IWYU pragma: associated
0021 #include "detail/UrbanMscMinimalStepLimit.hh"  // IWYU pragma: associated
0022 #include "detail/UrbanMscSafetyStepLimit.hh"  // IWYU pragma: associated
0023 #include "detail/UrbanMscScatter.hh"  // IWYU pragma: associated
0024
0025 namespace celeritas
0026 {
0027 //---------------------------------------------------------------------------//
0028 /*!
0029  * Apply Urban multiple scattering to a track.
0030  */
0031 class UrbanMsc
0032 {
0033   public:
0034     //!@{
0035     //! \name Type aliases
0036     using ParamsRef = NativeCRef<UrbanMscData>;
0037     //!@}
0038
0039   public:
0040     // Construct from MSC params
0041     explicit inline CELER_FUNCTION UrbanMsc(ParamsRef const& shared);
0042
0043     // Whether MSC applies to the current track
0044     inline CELER_FUNCTION bool
0045     is_applicable(CoreTrackView const&, real_type step) const;
0046
0047     // Update the physical and geometric step lengths
0048     inline CELER_FUNCTION void limit_step(CoreTrackView const&);
0049
0050     // Apply MSC
0051     inline CELER_FUNCTION void apply_step(CoreTrackView const&);
0052
0053   private:
0054     ParamsRef const shared_;
0055
0056     // Whether the step was limited by geometry
0057     static inline CELER_FUNCTION bool is_geo_limited(CoreTrackView const&);
0058 };
0059
0060 //---------------------------------------------------------------------------//
0061 // INLINE DEFINITIONS
0062 //---------------------------------------------------------------------------//
0063 /*!
0064  * Construct with defaults.
0065  */
0066 CELER_FUNCTION UrbanMsc::UrbanMsc(ParamsRef const& shared) : shared_(shared)
0067 {
0068     CELER_EXPECT(shared_);
0069 }
0070
0071 //---------------------------------------------------------------------------//
0072 /*!
0073  * Whether MSC applies to the current track.
0074  */
0075 CELER_FUNCTION bool
0076 UrbanMsc::is_applicable(CoreTrackView const& track, real_type step) const
0077 {
0078     if (step <= shared_.params.geom_limit)
0079         return false;
0080
0081     if (track.make_sim_view().status() != TrackStatus::alive)
0082         return false;
0083
0084     auto par = track.make_particle_view();
0085     if (par.particle_id() != shared_.ids.electron
0086         && par.particle_id() != shared_.ids.positron)
0087         return false;
0088
0089     return par.energy() > shared_.params.low_energy_limit
0090            && par.energy() < shared_.params.high_energy_limit;
0091 }
0092
0093 //---------------------------------------------------------------------------//
0094 /*!
0095  * Update the physical and geometric step lengths.
0096  */
0097 CELER_FUNCTION void UrbanMsc::limit_step(CoreTrackView const& track)
0098 {
0099     auto phys = track.make_physics_view();
0100     auto par = track.make_particle_view();
0101     auto sim = track.make_sim_view();
0102     detail::UrbanMscHelper msc_helper(shared_, par, phys);
0103
0104     bool displaced = false;
0105
0106     // Sample multiple scattering step length
0107     real_type const true_path = [&] {
0108         if (sim.step_length() <= shared_.params.limit_min_fix())
0109         {
0110             // Very short step: don't displace or limit
0111             return sim.step_length();
0112         }
0113
0114         auto geo = track.make_geo_view();
0115
0116         real_type safety = 0;
0117         if (!geo.is_on_boundary())
0118         {
0119             // Because the MSC behavior changes based on whether the *total*
0120             // track range is close to a boundary, rather than whether the next
0121             // steps are closer, we need to find the safety distance up to the
0122             // potential travel radius of the particle at its current energy.
0123             real_type const max_step = msc_helper.max_step();
0124             safety = geo.find_safety(max_step);
0125             if (safety >= max_step)
0126             {
0127                 // The nearest boundary is further than the maximum expected
0128                 // travel distance of the particle: don't displace or limit
0129                 return sim.step_length();
0130             }
0131         }
0132
0133         auto rng = track.make_rng_engine();
0134         displaced = true;
0135
0136         if (phys.scalars().step_limit_algorithm
0137             == MscStepLimitAlgorithm::minimal)
0138         {
0139             // Calculate step limit using "minimal" algorithm
0140             detail::UrbanMscMinimalStepLimit calc_limit(shared_,
0141                                                         msc_helper,
0142                                                         &phys,
0143                                                         geo.is_on_boundary(),
0144                                                         sim.step_length());
0145             return calc_limit(rng);
0146         }
0147
0148         // Calculate step limit using "safety" or "safety plus" algorithm
0149         detail::UrbanMscSafetyStepLimit calc_limit(shared_,
0150                                                    msc_helper,
0151                                                    par.energy(),
0152                                                    &phys,
0153                                                    phys.material_id(),
0154                                                    geo.is_on_boundary(),
0155                                                    safety,
0156                                                    sim.step_length());
0157         return calc_limit(rng);
0158
0159         // TODO: "distance to boundary" step limit algorithm
0160     }();
0161     CELER_ASSERT(true_path <= sim.step_length());
0162
0163     bool limited = (true_path < sim.step_length());
0164
0165     // Always apply the step transformation, even if the physical step wasn't
0166     // necessarily limited. This transformation will be reversed in
0167     // `apply_step` below.
0168     auto gp = [&] {
0169         detail::MscStepToGeo calc_geom_path(shared_,
0170                                             msc_helper,
0171                                             par.energy(),
0172                                             msc_helper.msc_mfp(),
0173                                             phys.dedx_range());
0174         auto gp = calc_geom_path(true_path);
0175
0176         // Limit geometrical step to 1 MSC MFP
0177         if (gp.step > msc_helper.msc_mfp())
0178         {
0179             gp.step = msc_helper.msc_mfp();
0180             limited = true;
0181         }
0182
0183         return gp;
0184     }();
0185     CELER_ASSERT(0 < gp.step && gp.step <= true_path);
0186
0187     // Save MSC step for later
0188     track.make_physics_step_view().msc_step([&] {
0189         MscStep result;
0190         result.is_displaced = displaced;
0191         result.true_path = true_path;
0192         result.geom_path = gp.step;
0193         result.alpha = gp.alpha;
0194         return result;
0195     }());
0196
0197     sim.step_length(gp.step);
0198     if (limited)
0199     {
0200         // Physical step was further limited by MSC
0201         sim.post_step_action(phys.scalars().msc_action());
0202     }
0203 }
0204
0205 //---------------------------------------------------------------------------//
0206 /*!
0207  * Apply MSC.
0208  */
0209 CELER_FUNCTION void UrbanMsc::apply_step(CoreTrackView const& track)
0210 {
0211     auto par = track.make_particle_view();
0212     auto geo = track.make_geo_view();
0213     auto phys = track.make_physics_view();
0214     auto sim = track.make_sim_view();
0215
0216     // Replace step with actual geometry distance traveled
0217     detail::UrbanMscHelper msc_helper(shared_, par, phys);
0218     auto msc_step = track.make_physics_step_view().msc_step();
0219     if (this->is_geo_limited(track))
0220     {
0221         // Convert geometrical distance to equivalent physical distance, which
0222         // will be greater than (or in edge cases equal to) that distance and
0223         // less than the original physical step limit.
0224         msc_step.geom_path = sim.step_length();
0225         detail::MscStepFromGeo geo_to_true(
0226             shared_.params, msc_step, phys.dedx_range(), msc_helper.msc_mfp());
0227         msc_step.true_path = geo_to_true(msc_step.geom_path);
0228         CELER_ASSERT(msc_step.true_path >= msc_step.geom_path);
0229
0230         // Disable displacement on boundary
0231         msc_step.is_displaced = false;
0232     }
0233
0234     // Update full path length traveled along the step based on MSC to
0235     // correctly calculate energy loss, step time, etc.
0236     sim.step_length(msc_step.true_path);
0237
0238     auto msc_result = [&] {
0239         real_type safety = 0;
0240         if (msc_step.is_displaced)
0241         {
0242             CELER_ASSERT(!geo.is_on_boundary());
0243             // Calculate the safety up to the maximum needed by UrbanMscScatter
0244             real_type displ = detail::UrbanMscScatter::calc_displacement(
0245                 msc_step.geom_path, msc_step.true_path);
0246             // The extra factor here is because UrbanMscScatter compares the
0247             // displacement length against the safety * (1 - eps), which we
0248             // bound here using [1 + 2*eps > 1/(1 - eps)], and we want to check
0249             // to at least the minimum geometry limit.
0250             // TODO: this is hacky and relies on UrbanMscScatter internals...
0251             displ = max(displ * (1 + 2 * shared_.params.safety_tol),
0252                         shared_.params.geom_limit);
0253             safety = geo.find_safety(displ);
0254             if (CELER_UNLIKELY(safety == 0))
0255             {
0256                 // The track is effectively on a boundary without being
0257                 // "logically" on, which is possible after tricky VecGeom
0258                 // boundary crossings.
0259                 msc_step.is_displaced = false;
0260             }
0261         }
0262
0263         auto mat = track.make_material_view().make_material_view();
0264         detail::UrbanMscScatter sample_scatter(
0265             shared_, msc_helper, par, phys, mat, geo.dir(), safety, msc_step);
0266
0267         auto rng = track.make_rng_engine();
0268         return sample_scatter(rng);
0269     }();
0270
0271     // Update direction and position
0272     if (msc_result.action != MscInteraction::Action::unchanged)
0273     {
0274         // Changing direction during a boundary crossing is OK
0275         geo.set_dir(msc_result.direction);
0276     }
0277     if (msc_result.action == MscInteraction::Action::displaced)
0278     {
0279         // Displacment during a boundary crossing is *not* OK
0280         CELER_ASSERT(!geo.is_on_boundary());
0281         geo.move_internal(geo.pos() + msc_result.displacement);
0282     }
0283 }
0284
0285 //---------------------------------------------------------------------------//
0286 /*!
0287  * Whether the step was limited by geometry.
0288  *
0289  * Usually the track is limited only if it's on the boundary (in which case
0290  * it should be "boundary action") but in rare circumstances the propagation
0291  * has to pause before the end of the step is reached.
0292  */
0293 CELER_FUNCTION bool UrbanMsc::is_geo_limited(CoreTrackView const& track)
0294 {
0295     auto sim = track.make_sim_view();
0296     return (sim.post_step_action() == track.boundary_action()
0297             || sim.post_step_action() == track.propagation_limit_action());
0298 }
0299
0300 //---------------------------------------------------------------------------//
0301 }  // namespace celeritas