alongstep/detail/FluctELoss.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/global/alongstep/detail/FluctELoss.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "celeritas/em/data/FluctuationData.hh"
0011 #include "celeritas/em/distribution/EnergyLossHelper.hh"
0012 #include "celeritas/em/distribution/EnergyLossTraits.hh"
0013 #include "celeritas/global/CoreTrackView.hh"
0014 #include "celeritas/phys/PhysicsStepUtils.hh"
0015
0016 namespace celeritas
0017 {
0018 namespace detail
0019 {
0020 //---------------------------------------------------------------------------//
0021 /*!
0022  * Apply energy loss (with fluctuations) to a track.
0023  */
0024 class FluctELoss
0025 {
0026   public:
0027     //!@{
0028     //! \name Type aliases
0029     using Energy = EnergyLossHelper::Energy;
0030     using ParamsRef = NativeCRef<FluctuationData>;
0031     //!@}
0032
0033   public:
0034     // Construct with fluctuation data
0035     inline explicit CELER_FUNCTION FluctELoss(ParamsRef const& params);
0036
0037     // Whether energy loss can be used for this track
0038     inline CELER_FUNCTION bool is_applicable(CoreTrackView const&) const;
0039
0040     // Apply to the track
0041     inline CELER_FUNCTION Energy calc_eloss(CoreTrackView const& track,
0042                                             real_type step,
0043                                             bool apply_cut);
0044
0045     //! Indicate that we can lose all energy before hitting the dE/dx range
0046     static CELER_CONSTEXPR_FUNCTION bool imprecise_range() { return true; }
0047
0048   private:
0049     //// DATA ////
0050
0051     //! Reference to fluctuation data
0052     ParamsRef const fluct_params_;
0053
0054     //// HELPER FUNCTIONS ////
0055
0056     template<EnergyLossFluctuationModel M>
0057     inline CELER_FUNCTION Energy
0058     sample_energy_loss(EnergyLossHelper const& helper, RngEngine& rng);
0059 };
0060
0061 //---------------------------------------------------------------------------//
0062 // INLINE DEFINITIONS
0063 //---------------------------------------------------------------------------//
0064 /*!
0065  * Construct with reference to fluctuation data.
0066  */
0067 CELER_FUNCTION FluctELoss::FluctELoss(ParamsRef const& params)
0068     : fluct_params_{params}
0069 {
0070     CELER_EXPECT(fluct_params_);
0071 }
0072
0073 //---------------------------------------------------------------------------//
0074 /*!
0075  * Whether energy loss is used for this track.
0076  */
0077 CELER_FUNCTION bool FluctELoss::is_applicable(CoreTrackView const& track) const
0078 {
0079     // The track can be marked as `errored` *within* the along-step kernel,
0080     // during propagation
0081     if (track.make_sim_view().status() == TrackStatus::errored)
0082         return false;
0083
0084     // Energy loss grid ID will be 'false' if inapplicable
0085     auto ppid = track.make_physics_view().eloss_ppid();
0086     return static_cast<bool>(ppid);
0087 }
0088
0089 //---------------------------------------------------------------------------//
0090 /*!
0091  * Apply energy loss to the given track.
0092  *
0093  * - Before and after slowing down we apply a tracking cut to cull low-energy
0094  *   charged particles.
0095  * - If energy loss fluctuations are enabled, we apply those based on the mean
0096  *   energy loss.
0097  * - If the sampled energy loss is greater than or equal to the particle's
0098  *   energy, we reduce it to the particle energy (if energy cuts are to be
0099  *   applied) or to the mean energy loss (if cuts are prohibited due to this
0100  *   being a non-physics-based step).
0101  */
0102 CELER_FUNCTION auto FluctELoss::calc_eloss(CoreTrackView const& track,
0103                                            real_type step,
0104                                            bool apply_cut) -> Energy
0105 {
0106     CELER_EXPECT(step > 0);
0107
0108     auto particle = track.make_particle_view();
0109     auto phys = track.make_physics_view();
0110
0111     if (apply_cut && particle.energy() < phys.scalars().lowest_electron_energy)
0112     {
0113         // Deposit all energy immediately when we start below the tracking cut
0114         return particle.energy();
0115     }
0116
0117     // Calculate mean energy loss
0118     auto eloss = calc_mean_energy_loss(particle, phys, step);
0119     CELER_EXPECT(eloss > zero_quantity());
0120
0121     if (eloss < particle.energy())
0122     {
0123         // Apply energy loss fluctuations
0124         auto cutoffs = track.make_cutoff_view();
0125         auto material = track.make_material_view();
0126
0127         EnergyLossHelper loss_helper(
0128             fluct_params_, cutoffs, material, particle, eloss, step);
0129
0130         auto rng = track.make_rng_engine();
0131         switch (loss_helper.model())
0132         {
0133 #define ASU_SAMPLE_ELOSS(MODEL)                                              \
0134     case EnergyLossFluctuationModel::MODEL:                                  \
0135         eloss = this->sample_energy_loss<EnergyLossFluctuationModel::MODEL>( \
0136             loss_helper, rng);                                               \
0137         break
0138             ASU_SAMPLE_ELOSS(none);
0139             ASU_SAMPLE_ELOSS(gamma);
0140             ASU_SAMPLE_ELOSS(gaussian);
0141             ASU_SAMPLE_ELOSS(urban);
0142 #undef ASU_SAMPLE_ELOSS
0143         }
0144
0145         if (eloss >= particle.energy())
0146         {
0147             // Sampled energy loss can be greater than actual remaining energy
0148             // because the range calculation is based on the *mean* energy
0149             // loss.
0150             if (apply_cut)
0151             {
0152                 // Clamp to actual particle energy so that it stops
0153                 eloss = particle.energy();
0154             }
0155             else
0156             {
0157                 // Don't go to zero energy at geometry boundaries: just use the
0158                 // mean loss which should be positive because this isn't a
0159                 // range-limited step.
0160                 eloss = loss_helper.mean_loss();
0161                 CELER_ASSERT(eloss < particle.energy());
0162             }
0163         }
0164     }
0165
0166     if (apply_cut
0167         && (particle.energy() - eloss <= phys.scalars().lowest_electron_energy))
0168     {
0169         // Deposit all energy when we end below the tracking cut
0170         return particle.energy();
0171     }
0172
0173     CELER_ASSERT(eloss <= particle.energy());
0174     CELER_ENSURE(eloss != particle.energy() || apply_cut
0175                  || track.make_sim_view().post_step_action()
0176                         == phys.scalars().range_action());
0177     return eloss;
0178 }
0179
0180 //---------------------------------------------------------------------------//
0181 template<EnergyLossFluctuationModel M>
0182 CELER_FUNCTION auto
0183 FluctELoss::sample_energy_loss(EnergyLossHelper const& helper,
0184                                RngEngine& rng) -> Energy
0185 {
0186     CELER_EXPECT(helper.model() == M);
0187
0188     using Distribution = typename EnergyLossTraits<M>::type;
0189
0190     Distribution sample_eloss{helper};
0191     return sample_eloss(rng);
0192 }
0193
0194 //---------------------------------------------------------------------------//
0195 }  // namespace detail
0196 }  // namespace celeritas