celeritas/phys/InteractionApplier.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/phys/InteractionApplier.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/Macros.hh"
0010 #include "corecel/cont/Span.hh"
0011 #include "corecel/sys/KernelTraits.hh"
0012 #include "celeritas/Quantities.hh"
0013 #include "celeritas/Types.hh"
0014 #include "celeritas/geo/GeoFwd.hh"
0015 #include "celeritas/global/CoreTrackView.hh"
0016 #include "celeritas/track/SimTrackView.hh"
0017
0018 #include "CutoffView.hh"
0019 #include "Interaction.hh"
0020 #include "ParticleTrackView.hh"
0021 #include "ParticleView.hh"
0022 #include "PhysicsData.hh"
0023 #include "PhysicsStepView.hh"
0024 #include "PhysicsTrackView.hh"
0025 #include "Secondary.hh"
0026
0027 namespace celeritas
0028 {
0029 //---------------------------------------------------------------------------//
0030 /*!
0031  * Wrap an Interaction executor to apply it to a track.
0032  *
0033  * The function F must take a \c CoreTrackView and return a \c Interaction
0034  */
0035 template<class F>
0036 struct InteractionApplierBaseImpl
0037 {
0038     //// DATA ////
0039
0040     F sample_interaction;
0041
0042     //// METHODS ////
0043
0044     CELER_FUNCTION void operator()(celeritas::CoreTrackView const&);
0045 };
0046
0047 //---------------------------------------------------------------------------//
0048 /*!
0049  *
0050  * This class is partially specialized with a second template argument to
0051  * extract any launch bounds from the functor class. TODO: we could probably
0052  * inherit from a helper class to pull in those constants (if available).
0053  */
0054 template<class F, typename = void>
0055 struct InteractionApplier : public InteractionApplierBaseImpl<F>
0056 {
0057     CELER_FUNCTION InteractionApplier(F&& f)
0058         : InteractionApplierBaseImpl<F>{celeritas::forward<F>(f)}
0059     {
0060     }
0061 };
0062
0063 template<class F>
0064 struct InteractionApplier<F, std::enable_if_t<kernel_max_blocks_min_warps<F>>>
0065     : public InteractionApplierBaseImpl<F>
0066 {
0067     static constexpr int max_block_size = F::max_block_size;
0068     static constexpr int min_warps_per_eu = F::min_warps_per_eu;
0069
0070     CELER_FUNCTION InteractionApplier(F&& f)
0071         : InteractionApplierBaseImpl<F>{celeritas::forward<F>(f)}
0072     {
0073     }
0074 };
0075
0076 template<class F>
0077 struct InteractionApplier<F, std::enable_if_t<kernel_max_blocks<F>>>
0078     : public InteractionApplierBaseImpl<F>
0079 {
0080     static constexpr int max_block_size = F::max_block_size;
0081
0082     CELER_FUNCTION InteractionApplier(F&& f)
0083         : InteractionApplierBaseImpl<F>{celeritas::forward<F>(f)}
0084     {
0085     }
0086 };
0087
0088 //---------------------------------------------------------------------------//
0089 // DEDUCTION GUIDES
0090 //---------------------------------------------------------------------------//
0091 template<class F>
0092 CELER_CTAD_FUNCTION InteractionApplier(F&&) -> InteractionApplier<F>;
0093
0094 //---------------------------------------------------------------------------//
0095 // INLINE DEFINITIONS
0096 //---------------------------------------------------------------------------//
0097 /*!
0098  * Sample an interaction and apply to the track view.
0099  *
0100  * The given track *must* be an active track with the correct step limit action
0101  * ID.
0102  */
0103 template<class F>
0104 CELER_FUNCTION void InteractionApplierBaseImpl<F>::operator()(
0105     celeritas::CoreTrackView const& track)
0106 {
0107     Interaction result = this->sample_interaction(track);
0108
0109     auto sim = track.sim();
0110     if (CELER_UNLIKELY(result.action == Interaction::Action::failed))
0111     {
0112         auto phys = track.physics();
0113         // Particle already moved to the collision site, but an out-of-memory
0114         // (allocation failure) occurred. Someday we can add error handling,
0115         // but for now use the "failure" action in the physics and set the step
0116         // limit to zero since it needs to interact again at this location.
0117         sim.step_limit({0, phys.scalars().failure_action()});
0118         return;
0119     }
0120     else if (!result.changed())
0121     {
0122         return;
0123     }
0124
0125     // Scattered or absorbed
0126     {
0127         // Update post-step energy
0128         auto particle = track.particle();
0129         particle.energy(result.energy);
0130     }
0131
0132     if (result.action != Interaction::Action::absorbed)
0133     {
0134         // Update direction
0135         track.geometry().set_dir(result.direction);
0136     }
0137     else
0138     {
0139         // Mark particle as dead
0140         sim.status(TrackStatus::killed);
0141     }
0142
0143     using Energy = units::MevEnergy;
0144     using MassCSq = units::MevMass;
0145
0146     real_type deposition = value_as<Energy>(result.energy_deposition);
0147     auto cutoff = track.cutoff();
0148     if (cutoff.apply_post_interaction())
0149     {
0150         for (auto& secondary : result.secondaries)
0151         {
0152             if (cutoff.apply(secondary))
0153             {
0154                 // Secondary is an electron, positron or gamma with energy
0155                 // below the production cut -- deposit the energy locally
0156                 // and clear the secondary
0157                 deposition += value_as<Energy>(secondary.energy) * sim.weight();
0158                 auto sec_par = track.particle_record(secondary.particle_id);
0159                 if (sec_par.is_antiparticle())
0160                 {
0161                     // Conservation of energy for positrons
0162                     // TODO: This is unphysical (but is consistent with
0163                     // Geant4): we should ONLY do this if the gamma tracking
0164                     // cut is also below 511 keV. We should probably just stop
0165                     // the positron and deposit its energy locally, but leave
0166                     // the secondary.
0167                     deposition += 2 * value_as<MassCSq>(sec_par.mass());
0168                 }
0169                 secondary = {};
0170
0171                 // Mark as a cut track
0172                 using namespace celeritas::literals;
0173                 atomic_add(
0174                     &const_cast<CoreTrackView&>(track).counters().num_cut,
0175                     1_sz);
0176             }
0177         }
0178     }
0179
0180     // Deposit energy and save secondaries
0181     auto phys = track.physics_step();
0182     if (deposition > 0)
0183     {
0184         phys.deposit_energy(Energy{deposition});
0185     }
0186     if (!result.secondaries.empty())
0187     {
0188         phys.secondaries(result.secondaries);
0189     }
0190 }
0191
0192 //---------------------------------------------------------------------------//
0193 }  // namespace celeritas