celeritas/phys/PhysicsTrackView.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/PhysicsTrackView.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/Config.hh"
0010
0011 #include "corecel/Assert.hh"
0012 #include "corecel/Macros.hh"
0013 #include "corecel/Types.hh"
0014 #include "celeritas/Quantities.hh"
0015 #include "celeritas/Types.hh"
0016 #include "celeritas/em/xs/EPlusGGMacroXsCalculator.hh"
0017 #include "celeritas/em/xs/ElectroNuclearMicroXsCalculator.hh"
0018 #include "celeritas/em/xs/GammaNuclearMicroXsCalculator.hh"
0019 #include "celeritas/em/xs/LivermorePEMicroXsCalculator.hh"
0020 #include "celeritas/grid/GridIdFinder.hh"
0021 #include "celeritas/grid/XsCalculator.hh"
0022 #include "celeritas/mat/MaterialView.hh"
0023 #include "celeritas/neutron/xs/NeutronElasticMicroXsCalculator.hh"
0024 #include "celeritas/random/TabulatedElementSelector.hh"
0025
0026 #include "MacroXsCalculator.hh"
0027 #include "ParticleTrackView.hh"
0028 #include "PhysicsData.hh"
0029
0030 namespace celeritas
0031 {
0032 //---------------------------------------------------------------------------//
0033 /*!
0034  * Physics data for a track.
0035  *
0036  * The physics track view provides an interface for data and operations
0037  * common to most processes and models.
0038  */
0039 class PhysicsTrackView
0040 {
0041   public:
0042     //!@{
0043     //! \name Type aliases
0044     using Initializer_t = PhysicsTrackInitializer;
0045     using PhysicsParamsRef = NativeCRef<PhysicsParamsData>;
0046     using PhysicsStateRef = NativeRef<PhysicsStateData>;
0047     using Energy = units::MevEnergy;
0048     using ModelFinder = GridIdFinder<Energy, ParticleModelId>;
0049     using UniformTableId = OpaqueId<UniformTable>;
0050     //!@}
0051
0052   public:
0053     // Construct from params, states, and per-state IDs
0054     inline CELER_FUNCTION PhysicsTrackView(PhysicsParamsRef const& params,
0055                                            PhysicsStateRef const& states,
0056                                            ParticleTrackView const& particle,
0057                                            PhysMatId material,
0058                                            TrackSlotId tid);
0059
0060     // Initialize the track view
0061     inline CELER_FUNCTION PhysicsTrackView& operator=(Initializer_t const&);
0062
0063     // Set the remaining MFP to interaction
0064     inline CELER_FUNCTION void interaction_mfp(real_type);
0065
0066     // Reset the remaining MFP to interaction
0067     inline CELER_FUNCTION void reset_interaction_mfp();
0068
0069     // Set the energy loss range for the current material and particle energy
0070     inline CELER_FUNCTION void dedx_range(real_type);
0071
0072     // Set the range properties for multiple scattering
0073     inline CELER_FUNCTION void msc_range(MscRange const&);
0074
0075     //// DYNAMIC PROPERTIES (pure accessors, free) ////
0076
0077     // Whether the remaining MFP has been calculated
0078     CELER_FORCEINLINE_FUNCTION bool has_interaction_mfp() const;
0079
0080     // Remaining MFP to interaction [1]
0081     CELER_FORCEINLINE_FUNCTION real_type interaction_mfp() const;
0082
0083     // Energy loss range for the current material and particle energy
0084     CELER_FORCEINLINE_FUNCTION real_type dedx_range() const;
0085
0086     // Range properties for multiple scattering
0087     CELER_FORCEINLINE_FUNCTION MscRange const& msc_range() const;
0088
0089     // Current material identifier
0090     CELER_FORCEINLINE_FUNCTION PhysMatId material_id() const;
0091
0092     //// PROCESSES (depend on particle type and possibly material) ////
0093
0094     // Number of processes that apply to this track
0095     inline CELER_FUNCTION ParticleProcessId::size_type
0096     num_particle_processes() const;
0097
0098     // Process ID for the given within-particle process index
0099     inline CELER_FUNCTION ProcessId process(ParticleProcessId) const;
0100
0101     // Get macro xs table, null if not present for this particle/material
0102     inline CELER_FUNCTION XsGridId macro_xs_grid(ParticleProcessId) const;
0103
0104     // Get energy loss table, null if not present for this particle/material
0105     inline CELER_FUNCTION UniformGridId energy_loss_grid() const;
0106
0107     // Get range table, null if not present for this particle/material
0108     inline CELER_FUNCTION UniformGridId range_grid() const;
0109
0110     // Get inverse range table, null if not present
0111     inline CELER_FUNCTION UniformGridId inverse_range_grid() const;
0112
0113     // Get data for processes that use the integral approach
0114     inline CELER_FUNCTION IntegralXsProcess const& integral_xs_process(
0115         ParticleProcessId ppid) const;
0116
0117     // Calculate macroscopic cross section for the process
0118     inline CELER_FUNCTION real_type calc_xs(ParticleProcessId ppid,
0119                                             MaterialView const& material,
0120                                             Energy energy) const;
0121
0122     // Estimate maximum macroscopic cross section for the process over the step
0123     inline CELER_FUNCTION real_type calc_max_xs(
0124         IntegralXsProcess const& process,
0125         ParticleProcessId ppid,
0126         MaterialView const& material,
0127         Energy energy) const;
0128
0129     // Models that apply to the given process ID
0130     inline
0131         CELER_FUNCTION ModelFinder make_model_finder(ParticleProcessId) const;
0132
0133     // Return CDF table data for the given particle/model/material
0134     inline CELER_FUNCTION UniformTableId cdf_table(ParticleModelId) const;
0135
0136     // Construct an element selector
0137     inline CELER_FUNCTION TabulatedElementSelector make_element_selector(
0138         UniformTableId, Energy) const;
0139
0140     // ID of the particle's at-rest process
0141     inline CELER_FUNCTION ParticleProcessId at_rest_process() const;
0142
0143     //// PARAMETER DATA ////
0144
0145     // Convert an action to a model ID for diagnostics, empty if not a model
0146     inline CELER_FUNCTION ModelId action_to_model(ActionId) const;
0147
0148     // Convert a selected model ID into a simulation action ID
0149     inline CELER_FUNCTION ActionId model_to_action(ModelId) const;
0150
0151     // Get the model ID corresponding to the given ParticleModelId
0152     inline CELER_FUNCTION ModelId model_id(ParticleModelId) const;
0153
0154     // Calculate scaled step range
0155     inline CELER_FUNCTION real_type range_to_step(real_type range) const;
0156
0157     // Access scalar properties
0158     CELER_FORCEINLINE_FUNCTION PhysicsParamsScalars const& scalars() const;
0159
0160     // Access particle-dependent scalar properties
0161     CELER_FORCEINLINE_FUNCTION ParticleScalars const& particle_scalars() const;
0162
0163     // Number of particle types
0164     inline CELER_FUNCTION size_type num_particles() const;
0165
0166     // Construct a grid calculator from a physics table
0167     template<class T>
0168     inline CELER_FUNCTION T make_calculator(UniformGridId) const;
0169
0170     //// HACKS ////
0171
0172     // Get hardwired model, null if not present
0173     inline CELER_FUNCTION ModelId hardwired_model(ParticleProcessId ppid,
0174                                                   Energy energy) const;
0175
0176   private:
0177     PhysicsParamsRef const& params_;
0178     PhysicsStateRef const& states_;
0179     ParticleId const particle_;
0180     PhysMatId const material_;
0181     TrackSlotId const track_slot_;
0182     bool is_heavy_;
0183
0184     //// IMPLEMENTATION HELPER FUNCTIONS ////
0185
0186     CELER_FORCEINLINE_FUNCTION PhysicsTrackState& state();
0187     CELER_FORCEINLINE_FUNCTION PhysicsTrackState const& state() const;
0188     CELER_FORCEINLINE_FUNCTION ProcessGroup const& process_group() const;
0189     inline CELER_FUNCTION UniformGridId uniform_grid(UniformTable const&) const;
0190 };
0191
0192 //---------------------------------------------------------------------------//
0193 // INLINE DEFINITIONS
0194 //---------------------------------------------------------------------------//
0195 /*!
0196  * Construct from shared and state data.
0197  *
0198  * Particle and material IDs are derived from other class states.
0199  */
0200 CELER_FUNCTION PhysicsTrackView::PhysicsTrackView(
0201     PhysicsParamsRef const& params,
0202     PhysicsStateRef const& states,
0203     ParticleTrackView const& particle,
0204     PhysMatId mid,
0205     TrackSlotId tid)
0206     : params_(params)
0207     , states_(states)
0208     , particle_(particle.particle_id())
0209     , material_(mid)
0210     , track_slot_(tid)
0211     , is_heavy_(particle.is_heavy())
0212 {
0213     CELER_EXPECT(track_slot_);
0214 }
0215
0216 //---------------------------------------------------------------------------//
0217 /*!
0218  * Initialize the track view.
0219  */
0220 CELER_FUNCTION PhysicsTrackView& PhysicsTrackView::operator=(
0221     Initializer_t const&)
0222 {
0223     this->state().interaction_mfp = 0;
0224     this->state().msc_range = {};
0225     return *this;
0226 }
0227
0228 //---------------------------------------------------------------------------//
0229 /*!
0230  * Set the distance to the next interaction, in mean free paths.
0231  *
0232  * This value will be decremented at each step.
0233  */
0234 CELER_FUNCTION void PhysicsTrackView::interaction_mfp(real_type mfp)
0235 {
0236     CELER_EXPECT(mfp > 0);
0237     this->state().interaction_mfp = mfp;
0238 }
0239
0240 //---------------------------------------------------------------------------//
0241 /*!
0242  * Set the distance to the next interaction, in mean free paths.
0243  *
0244  * This value will be decremented at each step.
0245  */
0246 CELER_FUNCTION void PhysicsTrackView::reset_interaction_mfp()
0247 {
0248     this->state().interaction_mfp = 0;
0249 }
0250
0251 //---------------------------------------------------------------------------//
0252 /*!
0253  * Set the energy loss range for the current material and particle energy.
0254  *
0255  * This value will be calculated once at the beginning of each step.
0256  */
0257 CELER_FUNCTION void PhysicsTrackView::dedx_range(real_type range)
0258 {
0259     CELER_EXPECT(range > 0);
0260     this->state().dedx_range = range;
0261 }
0262
0263 //---------------------------------------------------------------------------//
0264 /*!
0265  * Set the range properties for multiple scattering.
0266  *
0267  * These values will be calculated at the first step in every tracking volume.
0268  */
0269 CELER_FUNCTION void PhysicsTrackView::msc_range(MscRange const& msc_range)
0270 {
0271     this->state().msc_range = msc_range;
0272 }
0273
0274 //---------------------------------------------------------------------------//
0275 /*!
0276  * Current material identifier.
0277  */
0278 CELER_FUNCTION PhysMatId PhysicsTrackView::material_id() const
0279 {
0280     return material_;
0281 }
0282
0283 //---------------------------------------------------------------------------//
0284 /*!
0285  * Whether the remaining MFP has been calculated.
0286  */
0287 CELER_FUNCTION bool PhysicsTrackView::has_interaction_mfp() const
0288 {
0289     return this->state().interaction_mfp > 0;
0290 }
0291
0292 //---------------------------------------------------------------------------//
0293 /*!
0294  * Remaining MFP to interaction.
0295  */
0296 CELER_FUNCTION real_type PhysicsTrackView::interaction_mfp() const
0297 {
0298     real_type mfp = this->state().interaction_mfp;
0299     CELER_ENSURE(mfp >= 0);
0300     return mfp;
0301 }
0302
0303 //---------------------------------------------------------------------------//
0304 /*!
0305  * Energy loss range.
0306  */
0307 CELER_FUNCTION real_type PhysicsTrackView::dedx_range() const
0308 {
0309     real_type range = this->state().dedx_range;
0310     CELER_ENSURE(range > 0);
0311     return range;
0312 }
0313
0314 //---------------------------------------------------------------------------//
0315 /*!
0316  * Persistent range properties for multiple scattering within a same volume.
0317  */
0318 CELER_FUNCTION MscRange const& PhysicsTrackView::msc_range() const
0319 {
0320     return this->state().msc_range;
0321 }
0322
0323 //---------------------------------------------------------------------------//
0324 /*!
0325  * Number of processes that apply to this track.
0326  */
0327 CELER_FUNCTION ParticleProcessId::size_type
0328 PhysicsTrackView::num_particle_processes() const
0329 {
0330     return this->process_group().size();
0331 }
0332
0333 //---------------------------------------------------------------------------//
0334 /*!
0335  * Process ID for the given within-particle process index.
0336  */
0337 CELER_FUNCTION ProcessId PhysicsTrackView::process(ParticleProcessId ppid) const
0338 {
0339     CELER_EXPECT(ppid < this->num_particle_processes());
0340     return params_.process_ids[this->process_group().processes[ppid.get()]];
0341 }
0342
0343 //---------------------------------------------------------------------------//
0344 /*!
0345  * Return macro xs value grid data for the given process if available.
0346  */
0347 CELER_FUNCTION XsGridId PhysicsTrackView::macro_xs_grid(
0348     ParticleProcessId ppid) const
0349 {
0350     CELER_EXPECT(ppid < this->num_particle_processes());
0351     auto table_id = this->process_group().macro_xs[ppid.get()];
0352     CELER_ASSERT(table_id);
0353     auto const& table = params_.xs_tables[table_id];
0354     if (!table)
0355     {
0356         // No table for this process
0357         return {};
0358     }
0359     CELER_ASSERT(material_ < table.grids.size());
0360     auto grid_id = table.grids[material_.get()];
0361     if (!grid_id)
0362     {
0363         // No table for this particular material
0364         return {};
0365     }
0366     return params_.xs_grid_ids[grid_id];
0367 }
0368
0369 //---------------------------------------------------------------------------//
0370 /*!
0371  * Return the energy loss grid data if available.
0372  */
0373 CELER_FUNCTION UniformGridId PhysicsTrackView::energy_loss_grid() const
0374 {
0375     return this->uniform_grid(this->process_group().energy_loss);
0376 }
0377
0378 //---------------------------------------------------------------------------//
0379 /*!
0380  * Return the range grid data if available.
0381  */
0382 CELER_FUNCTION UniformGridId PhysicsTrackView::range_grid() const
0383 {
0384     return this->uniform_grid(this->process_group().range);
0385 }
0386
0387 //---------------------------------------------------------------------------//
0388 /*!
0389  * Return the inverse range grid data if available.
0390  *
0391  * If spline interpolation is used, the inverse grid is explicitly stored with
0392  * the derivatives calculated using the range as the x values and the energy as
0393  * the y values.
0394  *
0395  * The grid and values are identical to the range grid (i.e., not inverted)
0396  * even if the inverse grid is explicitly stored: the inversion is done in the
0397  * \c InverseRangeCalculator .
0398  */
0399 CELER_FUNCTION UniformGridId PhysicsTrackView::inverse_range_grid() const
0400 {
0401     if (auto const& grid
0402         = this->uniform_grid(this->process_group().inverse_range))
0403     {
0404         return grid;
0405     }
0406     return this->range_grid();
0407 }
0408
0409 //---------------------------------------------------------------------------//
0410 /*!
0411  * Get data for processes that use the integral approach.
0412  *
0413  * Particles that have energy loss processes will have a different energy at
0414  * the pre- and post-step points. This means the assumption that the cross
0415  * section is constant along the step is no longer valid. Instead, Monte Carlo
0416  * integration can be used to sample the interaction length for the discrete
0417  * process with the correct probability from the exact distribution,
0418  * \f[
0419      p = 1 - \exp \left( -\int_{E_0}^{E_1} n \sigma(E) \dif s \right),
0420  * \f]
0421  * where \f$ E_0 \f$ is the pre-step energy, \f$ E_1 \f$ is the post-step
0422  * energy, \em n is the atom density, and \em s is the interaction length.
0423  *
0424  * At the start of the step, the maximum value of the cross section over the
0425  * step \f$ \sigma_{max} \f$ is estimated and used as the macroscopic cross
0426  * section for the process rather than \f$ \sigma_{E_0} \f$. After the step,
0427  * the new value of the cross section \f$ \sigma(E_1) \f$ is calculated, and
0428  * the discrete interaction for the process occurs with probability
0429  * \f[
0430      p = \frac{\sigma(E_1)}{\sigma_{\max}}.
0431  * \f]
0432  *
0433  * See section 7.4 of the Geant4 Physics Reference (release 10.6) for details.
0434  */
0435 CELER_FUNCTION auto PhysicsTrackView::integral_xs_process(
0436     ParticleProcessId ppid) const -> IntegralXsProcess const&
0437 {
0438     CELER_EXPECT(ppid < this->num_particle_processes());
0439     return params_.integral_xs[this->process_group().integral_xs[ppid.get()]];
0440 }
0441
0442 //---------------------------------------------------------------------------//
0443 /*!
0444  * Calculate macroscopic cross section for the process.
0445  */
0446 CELER_FUNCTION real_type PhysicsTrackView::calc_xs(
0447     ParticleProcessId ppid, MaterialView const& material, Energy energy) const
0448 {
0449     real_type result = 0;
0450
0451     if (auto model_id = this->hardwired_model(ppid, energy))
0452     {
0453         // Calculate macroscopic cross section on the fly for special
0454         // hardwired processes.
0455         if (model_id == params_.hardwired.ids.livermore_pe)
0456         {
0457             result = MacroXsCalculator<LivermorePEMicroXsCalculator>(
0458                 params_.hardwired.livermore_pe, material)(energy);
0459         }
0460
0461         else if (model_id == params_.hardwired.ids.eplusgg)
0462         {
0463             result = EPlusGGMacroXsCalculator(params_.hardwired.eplusgg,
0464                                               material)(energy);
0465         }
0466         else if (model_id == params_.hardwired.ids.electro_vd)
0467         {
0468             result = MacroXsCalculator<ElectroNuclearMicroXsCalculator>(
0469                 params_.hardwired.electro_vd, material)(energy);
0470         }
0471         else if (model_id == params_.hardwired.ids.bertini_qgs)
0472         {
0473             result = MacroXsCalculator<GammaNuclearMicroXsCalculator>(
0474                 params_.hardwired.bertini_qgs, material)(energy);
0475         }
0476         else if (model_id == params_.hardwired.ids.chips)
0477         {
0478             result = MacroXsCalculator<NeutronElasticMicroXsCalculator>(
0479                 params_.hardwired.chips, material)(energy);
0480         }
0481     }
0482     else if (auto grid_id = this->macro_xs_grid(ppid))
0483     {
0484         // Calculate cross section from the tabulated data
0485         CELER_ASSERT(grid_id < params_.xs_grids.size());
0486         result = XsCalculator(params_.xs_grids[grid_id], params_.reals)(energy);
0487     }
0488
0489     CELER_ENSURE(result >= 0);
0490     return result;
0491 }
0492
0493 //---------------------------------------------------------------------------//
0494 /*!
0495  * Estimate maximum macroscopic cross section for the process over the step.
0496  *
0497  * If this is a particle with an energy loss process, this returns the
0498  * estimate of the maximum cross section over the step. If the energy of the
0499  * global maximum of the cross section (calculated at initialization) is in the
0500  * interval \f$ [\xi E_0, E_0) \f$, where \f$ E_0 \f$ is the pre-step energy
0501  * and \f$ \xi \f$ is \c min_eprime_over_e (defined by default as \f$ \xi = 1 -
0502  * \alpha \f$, where \f$ \alpha \f$ is \c max_step_over_range), \f$
0503  * \sigma_{\max} \f$ is set to the global maximum.  Otherwise, \f$
0504  * \sigma_{\max} = \max( \sigma(E_0), \sigma(\xi E_0) ) \f$. If the cross
0505  * section is not monotonic in the interval \f$ [\xi E_0, E_0) \f$ and the
0506  * interval does not contain the global maximum, the post-step cross section
0507  * \f$ \sigma(E_1) \f$ may be larger than \f$ \sigma_{\max} \f$.
0508  */
0509 CELER_FUNCTION real_type PhysicsTrackView::calc_max_xs(
0510     IntegralXsProcess const& process,
0511     ParticleProcessId ppid,
0512     MaterialView const& material,
0513     Energy energy) const
0514 {
0515     CELER_EXPECT(process);
0516     CELER_EXPECT(material_ < process.energy_max_xs.size());
0517
0518     real_type energy_max_xs
0519         = params_.reals[process.energy_max_xs[material_.get()]];
0520     real_type energy_xi = energy.value() * params_.scalars.min_eprime_over_e;
0521     if (energy_max_xs >= energy_xi && energy_max_xs < energy.value())
0522     {
0523         return this->calc_xs(ppid, material, Energy{energy_max_xs});
0524     }
0525     return max(this->calc_xs(ppid, material, energy),
0526                this->calc_xs(ppid, material, Energy{energy_xi}));
0527 }
0528
0529 //---------------------------------------------------------------------------//
0530 /*!
0531  * Get a hardwired model for on-the-fly cross section calculation.
0532  *
0533  * This returns the model ID that applies to the given process ID and energy if
0534  * the process is hardwired to calculate macroscopic cross sections on the fly.
0535  * If the result is null, tables should be used for this process/energy.
0536  */
0537 CELER_FUNCTION ModelId PhysicsTrackView::hardwired_model(
0538     ParticleProcessId ppid, Energy energy) const
0539 {
0540     ProcessId process = this->process(ppid);
0541     if ((process == params_.hardwired.ids.photoelectric
0542          && energy < LivermoreElement::tabulated_threshold())
0543         || (process == params_.hardwired.ids.annihilation)
0544         || (process == params_.hardwired.ids.gamma_nuclear)
0545         || (process == params_.hardwired.ids.neutron_elastic))
0546     {
0547         auto find_model = this->make_model_finder(ppid);
0548         return this->model_id(find_model(energy));
0549     }
0550     // Not a hardwired process
0551     return {};
0552 }
0553
0554 //---------------------------------------------------------------------------//
0555 /*!
0556  * Models that apply to the given process ID.
0557  */
0558 CELER_FUNCTION auto PhysicsTrackView::make_model_finder(
0559     ParticleProcessId ppid) const -> ModelFinder
0560 {
0561     CELER_EXPECT(ppid < this->num_particle_processes());
0562     ModelGroup const& md
0563         = params_.model_groups[this->process_group().models[ppid.get()]];
0564     return ModelFinder(params_.reals[md.energy], params_.pmodel_ids[md.model]);
0565 }
0566
0567 //---------------------------------------------------------------------------//
0568 /*!
0569  * Return value table data for the given particle/model/material.
0570  *
0571  * A null result means either the model is material independent or the material
0572  * only has one element, so no cross section CDF tables are stored.
0573  */
0574 CELER_FUNCTION auto PhysicsTrackView::cdf_table(ParticleModelId pmid) const
0575     -> UniformTableId
0576 {
0577     CELER_EXPECT(pmid < params_.model_cdf.size());
0578
0579     // Get the CDF table for the given particle and model
0580     ModelCdfTable const& model_cdf = params_.model_cdf[pmid];
0581     if (!model_cdf)
0582     {
0583         // No tables stored for this model
0584         return {};
0585     }
0586
0587     // Get the value table ID for the current material
0588     CELER_ASSERT(material_ < model_cdf.tables.size());
0589     auto table_id = model_cdf.tables[material_.get()];
0590     CELER_ASSERT(table_id < params_.uniform_tables.size());
0591     if (!params_.uniform_tables[table_id])
0592     {
0593         // No tables stored for this material
0594         return {};
0595     }
0596     return table_id;
0597 }
0598
0599 //---------------------------------------------------------------------------//
0600 /*!
0601  * Construct an element selector to sample an element from tabulated xs data.
0602  */
0603 CELER_FUNCTION TabulatedElementSelector PhysicsTrackView::make_element_selector(
0604     UniformTableId table_id, Energy energy) const
0605 {
0606     CELER_EXPECT(table_id < params_.uniform_tables.size());
0607     auto const& table = params_.uniform_tables[table_id];
0608     return TabulatedElementSelector{table,
0609                                     params_.uniform_grids,
0610                                     params_.uniform_grid_ids,
0611                                     params_.reals,
0612                                     energy};
0613 }
0614
0615 //---------------------------------------------------------------------------//
0616 /*!
0617  * ID of the particle's at-rest process.
0618  *
0619  * If the particle can have a discrete interaction at rest, this returns the \c
0620  * ParticleProcessId of that process. Otherwise, it returns a null ID.
0621  */
0622 CELER_FUNCTION ParticleProcessId PhysicsTrackView::at_rest_process() const
0623 {
0624     return this->process_group().at_rest;
0625 }
0626
0627 //---------------------------------------------------------------------------//
0628 /*!
0629  * Convert an action to a model ID for diagnostics, false if not a model.
0630  */
0631 CELER_FUNCTION ModelId PhysicsTrackView::action_to_model(ActionId action) const
0632 {
0633     if (!action)
0634         return ModelId{};
0635
0636     // Rely on unsigned rollover if action ID is less than the first model
0637     ModelId::size_type result = action - params_.scalars.first_model_action;
0638     if (result >= params_.scalars.num_models)
0639         return ModelId{};
0640
0641     return ModelId{result};
0642 }
0643
0644 //---------------------------------------------------------------------------//
0645 /*!
0646  * Convert a selected model ID into a simulation action ID.
0647  */
0648 CELER_FUNCTION ActionId PhysicsTrackView::model_to_action(ModelId model) const
0649 {
0650     CELER_ASSERT(model < params_.scalars.num_models);
0651     return params_.scalars.first_model_action + model.unchecked_get();
0652 }
0653
0654 //---------------------------------------------------------------------------//
0655 /*!
0656  * Get the model ID corresponding to the given ParticleModelId.
0657  */
0658 CELER_FUNCTION ModelId PhysicsTrackView::model_id(ParticleModelId pmid) const
0659 {
0660     CELER_EXPECT(pmid < params_.model_ids.size());
0661     return params_.model_ids[pmid];
0662 }
0663
0664 //---------------------------------------------------------------------------//
0665 /*!
0666  * Calculate scaled step range.
0667  *
0668  * This is the updated step function given by Eq. 7.4 of Geant4 Physics
0669  * Reference Manual, Release 10.6: \f[
0670    s = \alpha r + \rho (1 - \alpha) (2 - \frac{\rho}{r})
0671  \f]
0672  * where alpha is \c max_step_over_range and rho is \c min_range .
0673  *
0674  * Below \c min_range, no step scaling is applied, but the step can still
0675  * be arbitrarily small.
0676  *
0677  * \todo Rename \c calc_eloss_step_limit . This step limiter allows tuning the
0678  * accuracy loss from approximating a constant cross section along the step. We
0679  * should also split this into limiting the \em actual range (where the energy
0680  * goes to zero or the minimum allowable tracking range) versus a dE/dx
0681  * limiter.
0682  */
0683 CELER_FUNCTION real_type PhysicsTrackView::range_to_step(real_type range) const
0684 {
0685     CELER_ASSERT(range >= 0);
0686     auto const& scalars = this->particle_scalars();
0687     real_type const rho = scalars.min_range;
0688     if (range < rho * (1 + celeritas::sqrt_tol<real_type>))
0689     {
0690         // Small range returns the step. The fudge factor avoids floating point
0691         // error in the interpolation below while preserving the near-linear
0692         // behavior for range = rho + epsilon.
0693         return range;
0694     }
0695
0696     real_type const alpha = scalars.max_step_over_range;
0697     real_type step = alpha * range + rho * (1 - alpha) * (2 - rho / range);
0698     CELER_ENSURE(step > 0 && step <= range);
0699     return step;
0700 }
0701
0702 //---------------------------------------------------------------------------//
0703 /*!
0704  * Access scalar properties (options, IDs).
0705  */
0706 CELER_FORCEINLINE_FUNCTION PhysicsParamsScalars const&
0707 PhysicsTrackView::scalars() const
0708 {
0709     return params_.scalars;
0710 }
0711
0712 //---------------------------------------------------------------------------//
0713 /*!
0714  * Access particle-dependent scalar properties.
0715  *
0716  * These properties are different for light particles (electrons/positrons) and
0717  * heavy particles (muons/hadrons).
0718  */
0719 CELER_FORCEINLINE_FUNCTION ParticleScalars const&
0720 PhysicsTrackView::particle_scalars() const
0721 {
0722     if (is_heavy_)
0723     {
0724         return params_.scalars.heavy;
0725     }
0726     return params_.scalars.light;
0727 }
0728
0729 //---------------------------------------------------------------------------//
0730 /*!
0731  * Number of particle types.
0732  */
0733 CELER_FUNCTION size_type PhysicsTrackView::num_particles() const
0734 {
0735     return params_.process_groups.size();
0736 }
0737
0738 //---------------------------------------------------------------------------//
0739 /*!
0740  * Construct a grid calculator of the given type.
0741  *
0742  * The calculator must take two arguments: a reference to \c UniformGridRecord
0743  * and a reference to the backend reals storage.
0744  */
0745 template<class T>
0746 CELER_FUNCTION T PhysicsTrackView::make_calculator(UniformGridId id) const
0747 {
0748     CELER_EXPECT(id < params_.uniform_grids.size());
0749     return T{params_.uniform_grids[id], params_.reals};
0750 }
0751
0752 //---------------------------------------------------------------------------//
0753 // IMPLEMENTATION HELPER FUNCTIONS
0754 //---------------------------------------------------------------------------//
0755 /*!
0756  * Return value grid data for the given table ID if available.
0757  *
0758  * If the result is not null, it can be used to instantiate a
0759  * grid Calculator.
0760  *
0761  * If the result is null, it's likely because the process doesn't have the
0762  * associated value (e.g. if the table type is "energy_loss" and the process is
0763  * not a slowing-down process).
0764  */
0765 CELER_FUNCTION UniformGridId PhysicsTrackView::uniform_grid(
0766     UniformTable const& table) const
0767 {
0768     if (!table)
0769     {
0770         // No table for this process
0771         return {};
0772     }
0773     CELER_ASSERT(material_ < table.grids.size());
0774     auto grid_id = table.grids[material_.get()];
0775     if (!grid_id)
0776     {
0777         // No table for this particular material
0778         return {};
0779     }
0780     return params_.uniform_grid_ids[grid_id];
0781 }
0782
0783 //! Get the thread-local state (mutable)
0784 CELER_FUNCTION PhysicsTrackState& PhysicsTrackView::state()
0785 {
0786     return states_.state[track_slot_];
0787 }
0788
0789 //! Get the thread-local state (const)
0790 CELER_FUNCTION PhysicsTrackState const& PhysicsTrackView::state() const
0791 {
0792     return states_.state[track_slot_];
0793 }
0794
0795 //! Get the group of processes that apply to the particle
0796 CELER_FUNCTION ProcessGroup const& PhysicsTrackView::process_group() const
0797 {
0798     CELER_EXPECT(particle_ < params_.process_groups.size());
0799     return params_.process_groups[particle_];
0800 }
0801
0802 //---------------------------------------------------------------------------//
0803 }  // namespace celeritas