Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-02-22 10:31:29

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/PhysicsParams.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009 
0010 #include <memory>
0011 #include <utility>
0012 #include <vector>
0013 
0014 #include "corecel/Assert.hh"
0015 #include "corecel/Types.hh"
0016 #include "corecel/cont/Range.hh"
0017 #include "corecel/cont/Span.hh"
0018 #include "corecel/data/CollectionMirror.hh"
0019 #include "corecel/data/ParamsDataInterface.hh"
0020 #include "celeritas/Quantities.hh"
0021 #include "celeritas/Types.hh"
0022 #include "celeritas/Units.hh"
0023 #include "celeritas/global/ActionInterface.hh"
0024 
0025 #include "Model.hh"
0026 #include "PhysicsData.hh"
0027 #include "Process.hh"
0028 
0029 namespace celeritas
0030 {
0031 class ActionRegistry;
0032 class AtomicRelaxationParams;
0033 class MaterialParams;
0034 class ParticleParams;
0035 
0036 //---------------------------------------------------------------------------//
0037 /*!
0038  * Physics configuration options.
0039  *
0040  * Input options are:
0041  * - \c min_range: below this value, there is no extra transformation from
0042  *   particle range to step length.
0043  * - \c max_step_over_range: at higher energy (longer range), gradually
0044  *   decrease the maximum step length until it's this fraction of the tabulated
0045  *   range.
0046  * - \c fixed_step_limiter: if nonzero, prevent any tracks from taking a step
0047  *   longer than this length.
0048  * - \c min_eprime_over_e: energy scaling fraction used to estimate the maximum
0049  *   cross section over the step in the integral approach for energy loss
0050  *   processes.
0051  * - \c linear_loss_limit: if the mean energy loss along a step is greater than
0052  *   this fractional value of the pre-step kinetic energy, recalculate the
0053  *   energy loss.
0054  * - \c lowest_electron_energy: lowest kinetic energy for electrons/positrons
0055  * - \c lambda_limit: limit on the MSC mean free path.
0056  * - \c range_factor: used in the MSC step limitation algorithm to restrict the
0057  *   step size to \f$ f_r \cdot max(r, \lambda) \f$ at the start of a track or
0058  *   after entering a volume, where \f$ f_r \f$ is the range factor, \f$ r \f$
0059  *   is the range, and \f$ \lambda \f$ is the mean free path.
0060  * - \c safety_factor: used in the MSC step limitation algorithm to restrict
0061  *   the step size to \f$ f_s s \f$, where \f$ f_s \f$ is the safety factor and
0062  *   \f$  s \f$ is the safety distance.
0063  * - \c step_limit_algorithm: algorithm used to determine the MSC step limit.
0064  * - \c secondary_stack_factor: the number of secondary slots per track slot
0065  *   allocated.
0066  * - \c disable_integral_xs: for particles with energy loss processes, the
0067  *   particle energy changes over the step, so the assumption that the cross
0068  *   section is constant is no longer valid. By default, many charged particle
0069  *   processes use MC integration to sample the discrete interaction length
0070  *   with the correct probability. Disable this integral approach for all
0071  *   processes.
0072  *
0073  * NOTE: min_range/max_step_over_range are not accessible through Geant4, and
0074  * they can also be set to be different for electrons, mu/hadrons, and ions
0075  * (they are set in Geant4 with \c G4EmParameters::SetStepFunction()).
0076  */
0077 struct PhysicsParamsOptions
0078 {
0079     using Energy = units::MevEnergy;
0080 
0081     //!@{
0082     //! \name Range calculation
0083     real_type min_range = 1 * units::millimeter;
0084     real_type max_step_over_range = 0.2;
0085     real_type fixed_step_limiter = 0;
0086     //!@}
0087 
0088     //!@{
0089     //! \name Energy loss
0090     real_type min_eprime_over_e = 0.8;
0091     real_type linear_loss_limit = 0.01;
0092     Energy lowest_electron_energy = Energy{0.001};
0093     //!@}
0094 
0095     //!@{
0096     //! \name Multiple scattering
0097     real_type lambda_limit = 1 * units::millimeter;
0098     real_type range_factor = 0.04;
0099     real_type safety_factor = 0.6;
0100     MscStepLimitAlgorithm step_limit_algorithm{MscStepLimitAlgorithm::safety};
0101     //!@}
0102 
0103     real_type secondary_stack_factor = 3;
0104     bool disable_integral_xs = false;
0105 };
0106 
0107 //---------------------------------------------------------------------------//
0108 /*!
0109  * Manage physics processes and models.
0110  *
0111  * The physics params takes a vector of processes and sets up the processes and
0112  * models. It constructs data and mappings of data:
0113  * - particle type and process to tabulated values of cross sections etc,
0114  * - particle type to applicable processes
0115  *
0116  * During construction it constructs models and their corresponding list of
0117  * \c ActionId values, as well as the tables of cross section data. Besides the
0118  * individual interaction kernels, the physics parameters manage additional
0119  * actions:
0120  * - "pre-step": calculate physics step limits
0121  * - "along-step": propagate, apply energy loss, multiple scatter
0122  * - "range": limit step by energy loss
0123  * - "discrete-select": sample a process for a discrete interaction, or reject
0124  *   due to integral cross sectionl
0125  * - "integral-rejected": do not apply a discrete interaction
0126  * - "failure": model failed to allocate secondaries
0127  */
0128 class PhysicsParams final : public ParamsDataInterface<PhysicsParamsData>
0129 {
0130   public:
0131     //!@{
0132     //! \name Type aliases
0133     using SPConstParticles = std::shared_ptr<ParticleParams const>;
0134     using SPConstMaterials = std::shared_ptr<MaterialParams const>;
0135     using SPConstProcess = std::shared_ptr<Process const>;
0136     using SPConstModel = std::shared_ptr<Model const>;
0137     using SPConstRelaxation = std::shared_ptr<AtomicRelaxationParams const>;
0138 
0139     using VecProcess = std::vector<SPConstProcess>;
0140     using SpanConstProcessId = Span<ProcessId const>;
0141     using ActionIdRange = Range<ActionId>;
0142     using Options = PhysicsParamsOptions;
0143     //!@}
0144 
0145     //! Physics parameter construction arguments
0146     struct Input
0147     {
0148         SPConstParticles particles;
0149         SPConstMaterials materials;
0150         VecProcess processes;
0151         SPConstRelaxation relaxation;  //!< Optional atomic relaxation
0152         ActionRegistry* action_registry = nullptr;
0153 
0154         Options options;
0155     };
0156 
0157   public:
0158     // Construct with processes and helper classes
0159     explicit PhysicsParams(Input);
0160 
0161     //// HOST ACCESSORS ////
0162 
0163     //! Number of models
0164     ModelId::size_type num_models() const { return models_.size(); }
0165 
0166     //! Number of processes
0167     ProcessId::size_type num_processes() const { return processes_.size(); }
0168 
0169     // Number of particle types
0170     inline ParticleId::size_type num_particles() const;
0171 
0172     // Maximum number of processes that apply to any one particle
0173     inline ProcessId::size_type max_particle_processes() const;
0174 
0175     // Get a model
0176     inline SPConstModel const& model(ModelId) const;
0177 
0178     // Get a process
0179     inline SPConstProcess const& process(ProcessId) const;
0180 
0181     // Get the process for the given model
0182     inline ProcessId process_id(ModelId id) const;
0183 
0184     // Get the action IDs for all models
0185     inline ActionIdRange model_actions() const;
0186 
0187     // Get the processes that apply to a particular particle
0188     SpanConstProcessId processes(ParticleId) const;
0189 
0190     //! Access physics properties on the host
0191     HostRef const& host_ref() const final { return data_.host_ref(); }
0192 
0193     //! Access physics properties on the device
0194     DeviceRef const& device_ref() const final { return data_.device_ref(); }
0195 
0196   private:
0197     using SPAction = std::shared_ptr<StaticConcreteAction>;
0198     using VecModel = std::vector<std::pair<SPConstModel, ProcessId>>;
0199     using HostValue = celeritas::HostVal<PhysicsParamsData>;
0200 
0201     // Kernels/actions
0202     SPAction pre_step_action_;
0203     SPAction msc_action_;
0204     SPAction range_action_;
0205     SPAction discrete_action_;
0206     SPAction integral_rejection_action_;
0207     SPAction failure_action_;
0208     SPAction fixed_step_action_;
0209 
0210     // Host metadata/access
0211     VecProcess processes_;
0212     VecModel models_;
0213     SPConstRelaxation relaxation_;
0214 
0215     // Host/device storage and reference
0216     CollectionMirror<PhysicsParamsData> data_;
0217 
0218   private:
0219     VecModel build_models(ActionRegistry*) const;
0220     void build_options(Options const& opts, HostValue* data) const;
0221     void build_ids(ParticleParams const& particles, HostValue* data) const;
0222     void build_xs(Options const& opts,
0223                   MaterialParams const& mats,
0224                   HostValue* data) const;
0225     void build_model_xs(MaterialParams const& mats, HostValue* data) const;
0226 };
0227 
0228 //---------------------------------------------------------------------------//
0229 // INLINE DEFINITIONS
0230 //---------------------------------------------------------------------------//
0231 /*!
0232  * Number of particle types.
0233  */
0234 auto PhysicsParams::num_particles() const -> ParticleId::size_type
0235 {
0236     return this->host_ref().process_ids.size();
0237 }
0238 
0239 //---------------------------------------------------------------------------//
0240 /*!
0241  * Number of particle types.
0242  */
0243 auto PhysicsParams::max_particle_processes() const -> ProcessId::size_type
0244 {
0245     return this->host_ref().scalars.max_particle_processes;
0246 }
0247 
0248 //---------------------------------------------------------------------------//
0249 /*!
0250  * Get a model.
0251  */
0252 auto PhysicsParams::model(ModelId id) const -> SPConstModel const&
0253 {
0254     CELER_EXPECT(id < this->num_models());
0255     return models_[id.get()].first;
0256 }
0257 
0258 //---------------------------------------------------------------------------//
0259 /*!
0260  * Get a process.
0261  */
0262 auto PhysicsParams::process(ProcessId id) const -> SPConstProcess const&
0263 {
0264     CELER_EXPECT(id < this->num_processes());
0265     return processes_[id.get()];
0266 }
0267 
0268 //---------------------------------------------------------------------------//
0269 /*!
0270  * Get the process ID of the given model.
0271  */
0272 ProcessId PhysicsParams::process_id(ModelId id) const
0273 {
0274     CELER_EXPECT(id < this->num_models());
0275     return models_[id.get()].second;
0276 }
0277 
0278 //---------------------------------------------------------------------------//
0279 /*!
0280  * Get the action kernel IDs for all models.
0281  */
0282 auto PhysicsParams::model_actions() const -> ActionIdRange
0283 {
0284     auto offset = host_ref().scalars.model_to_action;
0285     return {ActionId{offset}, ActionId{offset + this->num_models()}};
0286 }
0287 
0288 //---------------------------------------------------------------------------//
0289 }  // namespace celeritas