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