celeritas/phys/ParticleTrackView.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2020-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/phys/ParticleTrackView.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include <cmath>
0011
0012 #include "corecel/Assert.hh"
0013 #include "corecel/Macros.hh"
0014 #include "corecel/Types.hh"
0015 #include "corecel/math/Algorithms.hh"
0016 #include "corecel/sys/ThreadId.hh"
0017 #include "celeritas/Quantities.hh"
0018
0019 #include "ParticleData.hh"
0020 #include "ParticleView.hh"
0021
0022 namespace celeritas
0023 {
0024 //---------------------------------------------------------------------------//
0025 /*!
0026  * Read/write view to the physical properties of a single particle track.
0027  *
0028  * These functions should be used in each physics Process or Interactor or
0029  * anything else that needs to access particle properties. Assume that all
0030  * these functions are expensive: when using them as accessors, locally store
0031  * the results rather than calling the function repeatedly. If any of the
0032  * calculations prove to be hot spots we will experiment with cacheing some of
0033  * the variables.
0034  */
0035 class ParticleTrackView
0036 {
0037   public:
0038     //!@{
0039     //! \name Type aliases
0040     using ParticleParamsRef = celeritas::NativeCRef<ParticleParamsData>;
0041     using ParticleStateRef = celeritas::NativeRef<ParticleStateData>;
0042     using Energy = units::MevEnergy;
0043     using Initializer_t = ParticleTrackInitializer;
0044     //!@}
0045
0046   public:
0047     // Construct from "dynamic" state and "static" particle definitions
0048     inline CELER_FUNCTION ParticleTrackView(ParticleParamsRef const& params,
0049                                             ParticleStateRef const& states,
0050                                             TrackSlotId id);
0051
0052     // Initialize the particle
0053     inline CELER_FUNCTION ParticleTrackView&
0054     operator=(Initializer_t const& other);
0055
0056     // Change the particle's energy [MeV]
0057     inline CELER_FUNCTION void energy(Energy);
0058
0059     // Reduce the particle's energy [MeV]
0060     inline CELER_FUNCTION void subtract_energy(Energy);
0061
0062     //// DYNAMIC PROPERTIES (pure accessors, free) ////
0063
0064     // Unique particle type identifier
0065     CELER_FORCEINLINE_FUNCTION ParticleId particle_id() const;
0066
0067     // Kinetic energy [MeV]
0068     CELER_FORCEINLINE_FUNCTION Energy energy() const;
0069
0070     // Whether the particle is stopped (zero kinetic energy)
0071     CELER_FORCEINLINE_FUNCTION bool is_stopped() const;
0072
0073     //// STATIC PROPERTIES (requires an indirection) ////
0074
0075     CELER_FORCEINLINE_FUNCTION ParticleView particle_view() const;
0076
0077     // Rest mass [MeV / c^2]
0078     CELER_FORCEINLINE_FUNCTION units::MevMass mass() const;
0079
0080     // Charge [elemental charge e+]
0081     CELER_FORCEINLINE_FUNCTION units::ElementaryCharge charge() const;
0082
0083     // Decay constant [1 / time]
0084     CELER_FORCEINLINE_FUNCTION real_type decay_constant() const;
0085
0086     // Whether it is an antiparticle
0087     CELER_FORCEINLINE_FUNCTION bool is_antiparticle() const;
0088
0089     // Whether the particle is stable
0090     CELER_FORCEINLINE_FUNCTION bool is_stable() const;
0091
0092     //// DERIVED PROPERTIES (indirection plus calculation) ////
0093
0094     // Kinetic energy plus rest energy [MeV]
0095     inline CELER_FUNCTION Energy total_energy() const;
0096
0097     // Square of fraction of lightspeed [unitless]
0098     inline CELER_FUNCTION real_type beta_sq() const;
0099
0100     // Speed [1/c]
0101     inline CELER_FUNCTION units::LightSpeed speed() const;
0102
0103     // Lorentz factor [unitless]
0104     inline CELER_FUNCTION real_type lorentz_factor() const;
0105
0106     // Relativistic momentum [MeV / c]
0107     inline CELER_FUNCTION units::MevMomentum momentum() const;
0108
0109     // Relativistic momentum squared [MeV^2 / c^2]
0110     inline CELER_FUNCTION units::MevMomentumSq momentum_sq() const;
0111
0112   private:
0113     ParticleParamsRef const& params_;
0114     ParticleStateRef const& states_;
0115     TrackSlotId const track_slot_;
0116 };
0117
0118 //---------------------------------------------------------------------------//
0119 // INLINE DEFINITIONS
0120 //---------------------------------------------------------------------------//
0121 /*!
0122  * Construct from dynamic and static particle properties.
0123  */
0124 CELER_FUNCTION
0125 ParticleTrackView::ParticleTrackView(ParticleParamsRef const& params,
0126                                      ParticleStateRef const& states,
0127                                      TrackSlotId tid)
0128     : params_(params), states_(states), track_slot_(tid)
0129 {
0130     CELER_EXPECT(track_slot_ < states_.size());
0131 }
0132
0133 //---------------------------------------------------------------------------//
0134 /*!
0135  * Initialize the particle.
0136  */
0137 CELER_FUNCTION ParticleTrackView&
0138 ParticleTrackView::operator=(Initializer_t const& other)
0139 {
0140     CELER_EXPECT(other.particle_id < params_.size());
0141     CELER_EXPECT(other.energy >= zero_quantity());
0142     states_.particle_id[track_slot_] = other.particle_id;
0143     states_.particle_energy[track_slot_] = other.energy.value();
0144     return *this;
0145 }
0146
0147 //---------------------------------------------------------------------------//
0148 /*!
0149  * Change the particle's kinetic energy.
0150  *
0151  * This should only be used when the particle is in a valid state. For HEP
0152  * applications, the new energy should always be less than the starting energy.
0153  */
0154 CELER_FUNCTION
0155 void ParticleTrackView::energy(Energy quantity)
0156 {
0157     CELER_EXPECT(this->particle_id());
0158     CELER_EXPECT(quantity >= zero_quantity());
0159     states_.particle_energy[track_slot_] = quantity.value();
0160 }
0161
0162 //---------------------------------------------------------------------------//
0163 /*!
0164  * Reduce the particle's energy without undergoing a collision [MeV].
0165  */
0166 CELER_FUNCTION void ParticleTrackView::subtract_energy(Energy eloss)
0167 {
0168     CELER_EXPECT(eloss >= zero_quantity());
0169     CELER_EXPECT(eloss <= this->energy());
0170     // TODO: save a read/write by only saving if eloss is positive?
0171     states_.particle_energy[track_slot_] -= eloss.value();
0172 }
0173
0174 //---------------------------------------------------------------------------//
0175 // DYNAMIC PROPERTIES
0176 //---------------------------------------------------------------------------//
0177 /*!
0178  * Unique particle type identifier.
0179  */
0180 CELER_FUNCTION ParticleId ParticleTrackView::particle_id() const
0181 {
0182     return states_.particle_id[track_slot_];
0183 }
0184
0185 //---------------------------------------------------------------------------//
0186 /*!
0187  * Kinetic energy [MeV].
0188  */
0189 CELER_FUNCTION auto ParticleTrackView::energy() const -> Energy
0190 {
0191     return Energy{states_.particle_energy[track_slot_]};
0192 }
0193
0194 //---------------------------------------------------------------------------//
0195 /*!
0196  * Whether the track is stopped (zero kinetic energy).
0197  */
0198 CELER_FUNCTION bool ParticleTrackView::is_stopped() const
0199 {
0200     return this->energy() == zero_quantity();
0201 }
0202
0203 //---------------------------------------------------------------------------//
0204 // STATIC PROPERTIES
0205 //---------------------------------------------------------------------------//
0206 /*!
0207  * Get static particle properties for the current state.
0208  */
0209 CELER_FUNCTION ParticleView ParticleTrackView::particle_view() const
0210 {
0211     return ParticleView(params_, states_.particle_id[track_slot_]);
0212 }
0213
0214 //---------------------------------------------------------------------------//
0215 /*!
0216  * Rest mass [MeV / c^2].
0217  */
0218 CELER_FUNCTION units::MevMass ParticleTrackView::mass() const
0219 {
0220     return this->particle_view().mass();
0221 }
0222
0223 //---------------------------------------------------------------------------//
0224 /*!
0225  * Elementary charge.
0226  */
0227 CELER_FUNCTION units::ElementaryCharge ParticleTrackView::charge() const
0228 {
0229     return this->particle_view().charge();
0230 }
0231
0232 //---------------------------------------------------------------------------//
0233 /*!
0234  * Decay constant in native units.
0235  */
0236 CELER_FUNCTION real_type ParticleTrackView::decay_constant() const
0237 {
0238     return this->particle_view().decay_constant();
0239 }
0240
0241 //---------------------------------------------------------------------------//
0242 /*!
0243  * Whether it is an antiparticle.
0244  */
0245 CELER_FUNCTION bool ParticleTrackView::is_antiparticle() const
0246 {
0247     return this->particle_view().is_antiparticle();
0248 }
0249
0250 //---------------------------------------------------------------------------//
0251 /*!
0252  * Whether particle is stable.
0253  */
0254 CELER_FUNCTION bool ParticleTrackView::is_stable() const
0255 {
0256     return this->decay_constant() == constants::stable_decay_constant;
0257 }
0258
0259 //---------------------------------------------------------------------------//
0260 // COMBINED PROPERTIES
0261 //---------------------------------------------------------------------------//
0262 /*!
0263  * Kinetic energy plus rest energy [MeV].
0264  */
0265 CELER_FUNCTION auto ParticleTrackView::total_energy() const -> Energy
0266 {
0267     return Energy(value_as<Energy>(this->energy())
0268                   + value_as<units::MevMass>(this->mass()));
0269 }
0270
0271 //---------------------------------------------------------------------------//
0272 /*!
0273  * Square of \f$ \beta \f$, which is the fraction of lightspeed [unitless].
0274  *
0275  * Beta is calculated using the equality
0276  * \f[
0277  * pc/E = \beta ; \quad \beta^2 = \frac{p^2 c^2}{E^2}
0278  * \f].
0279  * Using
0280  * \f[
0281  * E^2 = p^2 c^2 + m^2 c^4
0282  * \f]
0283  * and
0284  * \f[
0285  * E = K + mc^2
0286  * \f]
0287  *
0288  * the square of the fraction of lightspeed speed is
0289  * \f[
0290  * \beta^2 = 1 - \left( \frac{mc^2}{K + mc^2} \right)^2
0291  *         = 1 - \gamma^{-2}
0292  * \f]
0293  *
0294  * where \f$ \gamma \f$ is the Lorentz factor (see below).
0295  *
0296  * By choosing not to divide out the mass, this expression will work for
0297  * massless particles.
0298  *
0299  * \todo For nonrelativistic particles (low kinetic energy) this expression
0300  * may be inaccurate due to rounding error. It is however guaranteed to be in
0301  * the exact range [0, 1].
0302  */
0303 CELER_FUNCTION real_type ParticleTrackView::beta_sq() const
0304 {
0305     // Rest mass as energy
0306     real_type const mcsq = this->mass().value();
0307     // Inverse of lorentz factor (safe for m=0)
0308     real_type inv_gamma = mcsq / (this->energy().value() + mcsq);
0309
0310     return 1 - ipow<2>(inv_gamma);
0311 }
0312
0313 //---------------------------------------------------------------------------//
0314 /*!
0315  * Speed [1/c].
0316  *
0317  * Speed is calculated using beta so that the expression works for massless
0318  * particles.
0319  */
0320 CELER_FUNCTION units::LightSpeed ParticleTrackView::speed() const
0321 {
0322     return units::LightSpeed{std::sqrt(this->beta_sq())};
0323 }
0324
0325 //---------------------------------------------------------------------------//
0326 /*!
0327  * Lorentz factor [unitless].
0328  *
0329  * The Lorentz factor can be viewed as a transformation from
0330  * classical quantities to relativistic quantities. It's defined as
0331  * \f[
0332   \gamma = \frac{1}{\sqrt{1 - v^2 / c^2}}
0333   \f]
0334  *
0335  * Its value is infinite for the massless photon, and greater than or equal to
0336  * 1 otherwise.
0337  *
0338  * Gamma can also be calculated from the total (rest + kinetic) energy
0339  * \f[
0340   E = \gamma mc^2 = K + mc^2
0341   \f]
0342  * which we ues here since \em K and \em m are the primary stored quantities of
0343  * the particles:
0344  * \f[
0345   \gamma = 1 + \frac{K}{mc^2}
0346  * \f]
0347  */
0348 CELER_FUNCTION real_type ParticleTrackView::lorentz_factor() const
0349 {
0350     CELER_EXPECT(this->mass() > zero_quantity());
0351
0352     real_type k_over_mc2 = this->energy().value() / this->mass().value();
0353     return 1 + k_over_mc2;
0354 }
0355
0356 //---------------------------------------------------------------------------//
0357 /*!
0358  * Square of relativistic momentum [MeV^2 / c^2].
0359  *
0360  * Total energy:
0361  * \f[
0362  * E = K + mc^2
0363  * \f]
0364  * Relation between energy and momentum:
0365  * \f[
0366  * E^2 = p^2 c^2 + m^2 c^4
0367  * \f]
0368  * therefore
0369  * \f[
0370  * p^2 = \frac{E^2}{c^2} - m^2 c^2
0371  * \f]
0372  * or
0373  * \f[
0374  * p^2 = \frac{K^2}{c^2} + 2 * m * K
0375  * \f]
0376  */
0377 CELER_FUNCTION units::MevMomentumSq ParticleTrackView::momentum_sq() const
0378 {
0379     real_type const energy = this->energy().value();
0380     real_type result = ipow<2>(energy) + 2 * this->mass().value() * energy;
0381     CELER_ENSURE(result >= 0);
0382     return units::MevMomentumSq{result};
0383 }
0384
0385 //---------------------------------------------------------------------------//
0386 /*!
0387  * Relativistic momentum [MeV / c].
0388  *
0389  * This is calculated by taking the root of the square of the momentum.
0390  */
0391 CELER_FUNCTION units::MevMomentum ParticleTrackView::momentum() const
0392 {
0393     return units::MevMomentum{std::sqrt(this->momentum_sq().value())};
0394 }
0395
0396 //---------------------------------------------------------------------------//
0397 }  // namespace celeritas