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