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