diff --git a/src/accel/SharedParams.cc b/src/accel/SharedParams.cc index 101ed459c3..e8d9214508 100644 --- a/src/accel/SharedParams.cc +++ b/src/accel/SharedParams.cc @@ -68,7 +68,6 @@ #include "celeritas/phys/ParticleParams.hh" #include "celeritas/phys/PhysicsParams.hh" #include "celeritas/phys/Process.hh" -#include "celeritas/phys/ProcessBuilder.hh" #include "celeritas/setup/FrameworkInput.hh" #include "celeritas/track/SimParams.hh" #include "celeritas/track/TrackInitParams.hh" diff --git a/src/celeritas/CMakeLists.txt b/src/celeritas/CMakeLists.txt index da32d79623..11eeadfb54 100644 --- a/src/celeritas/CMakeLists.txt +++ b/src/celeritas/CMakeLists.txt @@ -36,6 +36,7 @@ list(APPEND PUBLIC_DEPS ${_core_geo_deps}) list(APPEND SOURCES Types.cc TypesIO.json.cc + decay/DecayProcess.cc em/detail/Utils.cc em/params/AtomicRelaxationParams.cc em/params/FluctuationParams.cc @@ -118,6 +119,7 @@ list(APPEND SOURCES optical/surface/SurfacePhysicsParams.cc optical/surface/SurfaceSteppingAction.cc phys/CutoffParams.cc + phys/detail/DecayTableInserter.cc phys/detail/EnergyMaxXsCalculator.cc phys/GeneratorInterface.cc phys/GeneratorRegistry.cc @@ -328,6 +330,7 @@ celeritas_polysource(alongstep/AlongStepNeutralAction) celeritas_polysource(alongstep/AlongStepUniformMscAction) celeritas_polysource(alongstep/AlongStepRZMapFieldMscAction) celeritas_polysource(alongstep/AlongStepCylMapFieldMscAction) +celeritas_polysource(decay/channel/MuDecayChannel) celeritas_polysource(em/model/BetheHeitlerModel) celeritas_polysource(em/model/BetheBlochModel) celeritas_polysource(em/model/BraggModel) diff --git a/src/celeritas/Types.cc b/src/celeritas/Types.cc index 42419553b8..3353e8d50c 100644 --- a/src/celeritas/Types.cc +++ b/src/celeritas/Types.cc @@ -90,6 +90,18 @@ char const* to_cstring(NuclearFormFactorType value) return to_cstring_impl(value); } +//---------------------------------------------------------------------------// +/*! + * Get a string corresponding to the decay channel. + */ +char const* to_cstring(DecayChannelType value) +{ + static EnumStringMapper const to_cstring_impl{ + "muon", + }; + return to_cstring_impl(value); +} + //---------------------------------------------------------------------------// /*! * Get a string corresponding to the wavelength shifting time model selection. diff --git a/src/celeritas/Types.hh b/src/celeritas/Types.hh index 9232c0a776..345860ad42 100644 --- a/src/celeritas/Types.hh +++ b/src/celeritas/Types.hh @@ -68,6 +68,9 @@ using TrackId = OpaqueId; //! Opaque index of particle-nucleon cascade channel using ChannelId = OpaqueId; +//! Opaque index of decay channel +using DecayChannelId = OpaqueId; + //! Opaque index to one elemental component datum in a particular material using ElementComponentId = OpaqueId; @@ -202,6 +205,14 @@ enum class NuclearFormFactorType size_ }; +//---------------------------------------------------------------------------// +//! Decay channel type +enum class DecayChannelType +{ + muon, + size_ +}; + //---------------------------------------------------------------------------// //! Optical photon wavelength shifting time model enum class WlsTimeProfile @@ -277,6 +288,9 @@ char const* to_cstring(MscStepLimitAlgorithm value); // Get a string corresponding to the nuclear form factor model char const* to_cstring(NuclearFormFactorType value); +// Get a string corresponding to the decay channel +char const* to_cstring(DecayChannelType value); + // Get a string corresponding to the interpolation method char const* to_cstring(InterpolationType value); diff --git a/src/celeritas/decay/DecayData.hh b/src/celeritas/decay/DecayData.hh new file mode 100644 index 0000000000..29b17ba121 --- /dev/null +++ b/src/celeritas/decay/DecayData.hh @@ -0,0 +1,54 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/decay/DecayData.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/Macros.hh" +#include "corecel/Types.hh" +#include "corecel/data/Collection.hh" +#include "celeritas/Types.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Data for a decay interactor. + * + * This stores the particle IDs of the daughters for a specific decay channel. + */ +struct DecayChannelData +{ + //! Daughter particle IDs + ItemRange daughters; + + //! Whether the data is assigned + explicit CELER_FUNCTION operator bool() const + { + return !daughters.empty(); + } +}; + +//---------------------------------------------------------------------------// +/*! + * Decay channels for a particle type. + */ +struct DecayTableData +{ + //! Decay channels + ItemRange channel_ids; + //! Branching ratio of each decay channel + ItemRange branching_ratios; + + //! Whether the data is assigned + explicit CELER_FUNCTION operator bool() const + { + return !channel_ids.empty() + && branching_ratios.size() == channel_ids.size(); + } +}; + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/decay/DecayProcess.cc b/src/celeritas/decay/DecayProcess.cc new file mode 100644 index 0000000000..6114c712fa --- /dev/null +++ b/src/celeritas/decay/DecayProcess.cc @@ -0,0 +1,90 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/decay/DecayProcess.cc +//---------------------------------------------------------------------------// +#include "DecayProcess.hh" + +#include + +#include "corecel/Assert.hh" +#include "celeritas/decay/DecayData.hh" +#include "celeritas/decay/channel/MuDecayChannel.hh" +#include "celeritas/phys/PDGNumber.hh" +#include "celeritas/phys/ParticleParams.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Construct from particles and input data. + */ +DecayProcess::DecayProcess(SPConstParticles particles, + inp::DecayPhysics const& input) + : particles_(std::move(particles)), input_(input) +{ + CELER_EXPECT(particles_); + CELER_VALIDATE(input_, << "no decay tables are present"); +} + +//---------------------------------------------------------------------------// +/*! + * Construct the decay channels. + */ +auto DecayProcess::build_channels(ActionIdIter start_id) const -> VecChannel +{ + VecChannel result; + std::set types; + for (auto const& [pdg, table] : input_.tables) + { + for (auto const& channel : table) + { + // Identify the unique channel types + auto [iter, inserted] = types.insert(channel.type); + if (inserted) + { + // Construct an action for each channel + switch (*iter) + { + case DecayChannelType::muon: + result.push_back(std::make_shared( + *start_id++, input_)); + break; + default: + CELER_NOT_IMPLEMENTED("Decay channel type"); + } + } + } + } + CELER_VALIDATE(!result.empty(), + << "no supported channels for decay process"); + return result; +} + +//---------------------------------------------------------------------------// +/*! + * Get the decay table for a particle. + */ +auto DecayProcess::decay_table(ParticleId pid) const -> DecayTable +{ + auto iter = input_.tables.find(particles_->id_to_pdg(pid)); + if (iter == input_.tables.end()) + { + return {}; + } + return iter->second; +} + +//---------------------------------------------------------------------------// +/*! + * Whether the decay process applies to the particle. + */ +bool DecayProcess::is_applicable(ParticleId pid) const +{ + return input_.tables.find(particles_->id_to_pdg(pid)) + != input_.tables.end(); +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/decay/DecayProcess.hh b/src/celeritas/decay/DecayProcess.hh new file mode 100644 index 0000000000..5259c29360 --- /dev/null +++ b/src/celeritas/decay/DecayProcess.hh @@ -0,0 +1,64 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/decay/DecayProcess.hh +//---------------------------------------------------------------------------// +#pragma once + +#include + +#include "celeritas/Types.hh" +#include "celeritas/decay/channel/DecayChannel.hh" +#include "celeritas/inp/Physics.hh" +#include "celeritas/phys/ImportedProcessAdapter.hh" +#include "celeritas/phys/Process.hh" + +namespace celeritas +{ +class ParticleParams; + +//---------------------------------------------------------------------------// +/*! + * Process for decay. + */ +class DecayProcess final : public Process +{ + public: + //!@{ + //! \name Type aliases + using SPConstChannel = std::shared_ptr; + using SPConstParticles = std::shared_ptr; + using VecChannel = std::vector; + using DecayTable = inp::DecayPhysics::DecayTable; + //!@} + + public: + // Construct from particles and input data + DecayProcess(SPConstParticles, inp::DecayPhysics const&); + + // Construct the decay channels + VecChannel build_channels(ActionIdIter start_id) const; + + // Get the decay table for a particle + DecayTable decay_table(ParticleId) const; + + // Whether the decay process applies to the particle + bool is_applicable(ParticleId) const; + + //!@{ + //! \name Process interface + + //! Whether the process applies when the particle is stopped + bool applies_at_rest() const final { return true; } + //! Name of the process + std::string_view label() const final { return "Decay"; } + //!@} + + private: + SPConstParticles particles_; + inp::DecayPhysics const& input_; +}; + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/decay/channel/DecayChannel.hh b/src/celeritas/decay/channel/DecayChannel.hh new file mode 100644 index 0000000000..9135406b7b --- /dev/null +++ b/src/celeritas/decay/channel/DecayChannel.hh @@ -0,0 +1,30 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/decay/channel/DecayChannel.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/Types.hh" +#include "celeritas/Types.hh" +#include "celeritas/global/ActionInterface.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Abstract base class representing a decay channel. + * + * This class is responsible for validating the input decay table data and + * launching a kernel for the decay interactor. + */ +class DecayChannel : public CoreStepActionInterface +{ + public: + //! Dependency ordering of the action + StepActionOrder order() const final { return StepActionOrder::post; } +}; + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/decay/channel/MuDecayChannel.cc b/src/celeritas/decay/channel/MuDecayChannel.cc new file mode 100644 index 0000000000..fafb873718 --- /dev/null +++ b/src/celeritas/decay/channel/MuDecayChannel.cc @@ -0,0 +1,78 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/decay/channel/MuDecayChannel.cc +//---------------------------------------------------------------------------// +#include "MuDecayChannel.hh" + +#include "celeritas/Quantities.hh" +#include "celeritas/decay/executor/MuDecayExecutor.hh" +#include "celeritas/em/interactor/detail/PhysicsConstants.hh" +#include "celeritas/global/ActionLauncher.hh" +#include "celeritas/global/CoreParams.hh" +#include "celeritas/global/CoreState.hh" +#include "celeritas/global/TrackExecutor.hh" +#include "celeritas/phys/InteractionApplier.hh" +#include "celeritas/phys/PDGNumber.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Construct from action ID and input decay table data. + */ +MuDecayChannel::MuDecayChannel(ActionId id, inp::DecayPhysics const& input) + : StaticConcreteAction(id, "muon-decay-channel", "muon decay") +{ + CELER_EXPECT(id); + CELER_EXPECT(input); + + // Validate the input data + for (auto const& [pdg, table] : input.tables) + { + for (auto const& channel : table) + { + // Loop over decay tables and find all channels of this type + if (channel.type == DecayChannelType::muon) + { + // Check that the input data matches the data expected by the + // muon decay interactor + auto const& daughters = channel.daughters; + CELER_VALIDATE(daughters.size() == 1, + << "muon decay only supports one daughter: " + "neutrinos are currently neglected"); + CELER_VALIDATE( + (pdg == pdg::mu_minus() && daughters[0] == pdg::electron()) + || (pdg == pdg::mu_plus() + && daughters[0] == pdg::positron()), + << "expected decay channel mu- -> e- or mu+ -> e+"); + } + } + } +} + +//---------------------------------------------------------------------------// +/*! + * Interact with host data. + */ +void MuDecayChannel::step(CoreParams const& params, CoreStateHost& state) const +{ + auto execute + = make_action_track_executor(params.ptr(), + state.ptr(), + this->action_id(), + InteractionApplier{MuDecayExecutor{}}); + return launch_action(*this, params, state, execute); +} + +//---------------------------------------------------------------------------// +#if !CELER_USE_DEVICE +void MuDecayChannel::step(CoreParams const&, CoreStateDevice&) const +{ + CELER_NOT_CONFIGURED("CUDA OR HIP"); +} +#endif + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/decay/channel/MuDecayChannel.cu b/src/celeritas/decay/channel/MuDecayChannel.cu new file mode 100644 index 0000000000..ac85d12f97 --- /dev/null +++ b/src/celeritas/decay/channel/MuDecayChannel.cu @@ -0,0 +1,34 @@ +//------------------------------ -*- cuda -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/decay/channel/MuDecayChannel.cu +//---------------------------------------------------------------------------// +#include "MuDecayChannel.hh" + +#include "celeritas/decay/executor/MuDecayExecutor.hh" +#include "celeritas/global/ActionLauncher.device.hh" +#include "celeritas/global/CoreParams.hh" +#include "celeritas/global/CoreState.hh" +#include "celeritas/global/TrackExecutor.hh" +#include "celeritas/phys/InteractionApplier.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Interact with device data. + */ +void MuDecayChannel::step(CoreParams const& params, CoreStateDevice& state) const +{ + auto execute + = make_action_track_executor(params.ptr(), + state.ptr(), + this->action_id(), + InteractionApplier{MuDecayExecutor{}}); + static ActionLauncher const launch_kernel(*this); + launch_kernel(*this, params, state, execute); +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/decay/channel/MuDecayChannel.hh b/src/celeritas/decay/channel/MuDecayChannel.hh new file mode 100644 index 0000000000..e8508b8ca8 --- /dev/null +++ b/src/celeritas/decay/channel/MuDecayChannel.hh @@ -0,0 +1,36 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/decay/channel/MuDecayChannel.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "celeritas/inp/Physics.hh" + +#include "DecayChannel.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Set up and launch an action to model muon decay kinematics. + */ +class MuDecayChannel : public DecayChannel, public StaticConcreteAction +{ + public: + // Construct from action ID and imported data + MuDecayChannel(ActionId, inp::DecayPhysics const&); + + //!@{ + //! \name StepAction interface + + // Apply the interaction kernel on host + void step(CoreParams const&, CoreStateHost&) const final; + // Apply the interaction kernel on device + void step(CoreParams const&, CoreStateDevice&) const final; + //!@} +}; + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/decay/data/MuDecayData.hh b/src/celeritas/decay/data/MuDecayData.hh deleted file mode 100644 index 3495e37caf..0000000000 --- a/src/celeritas/decay/data/MuDecayData.hh +++ /dev/null @@ -1,49 +0,0 @@ -//------------------------------- -*- C++ -*- -------------------------------// -// Copyright Celeritas contributors: see top-level COPYRIGHT file for details -// SPDX-License-Identifier: (Apache-2.0 OR MIT) -//---------------------------------------------------------------------------// -//! \file celeritas/decay/data/MuDecayData.hh -//---------------------------------------------------------------------------// -#pragma once - -#include "corecel/Macros.hh" -#include "corecel/Types.hh" -#include "celeritas/Quantities.hh" -#include "celeritas/Types.hh" - -namespace celeritas -{ -//---------------------------------------------------------------------------// -/*! - * Device data for creating a muon decay interactor. - */ -struct MuDecayData -{ - //!@{ - //! \name Type aliases - using Mass = units::MevMass; - //!@} - - //! Particle identifiers - ParticleId mu_minus_id; - ParticleId mu_plus_id; - ParticleId electron_id; - ParticleId positron_id; - - //! Muon/anti-muon mass [MeV] - Mass muon_mass; - - //! Electron/positron mass [MeV] - Mass electron_mass; - - //! Check whether the data is assigned - explicit CELER_FUNCTION operator bool() const - { - return mu_minus_id && mu_plus_id && electron_id && positron_id - && muon_mass > zero_quantity() - && electron_mass > zero_quantity(); - } -}; - -//---------------------------------------------------------------------------// -} // namespace celeritas diff --git a/src/celeritas/decay/executor/MuDecayExecutor.hh b/src/celeritas/decay/executor/MuDecayExecutor.hh index c7c188c938..1eec125b52 100644 --- a/src/celeritas/decay/executor/MuDecayExecutor.hh +++ b/src/celeritas/decay/executor/MuDecayExecutor.hh @@ -6,7 +6,6 @@ //---------------------------------------------------------------------------// #pragma once -#include "celeritas/decay/data/MuDecayData.hh" #include "celeritas/decay/interactor/MuDecayInteractor.hh" #include "celeritas/global/CoreTrackView.hh" #include "celeritas/phys/Interaction.hh" @@ -18,8 +17,6 @@ struct MuDecayExecutor { inline CELER_FUNCTION Interaction operator()(celeritas::CoreTrackView const& track); - - MuDecayData data; }; //---------------------------------------------------------------------------// @@ -31,8 +28,10 @@ CELER_FUNCTION Interaction MuDecayExecutor::operator()(CoreTrackView const& trac auto allocate_secondaries = track.physics_step().make_secondary_allocator(); auto particle = track.particle(); auto const& dir = track.geometry().dir(); + auto daughters + = track.physics().daughters(track.physics_step().decay_channel()); - MuDecayInteractor interact(data, particle, dir, allocate_secondaries); + MuDecayInteractor interact(daughters, particle, dir, allocate_secondaries); auto rng = track.rng(); return interact(rng); } diff --git a/src/celeritas/decay/interactor/MuDecayInteractor.hh b/src/celeritas/decay/interactor/MuDecayInteractor.hh index a51b000acd..e13f635159 100644 --- a/src/celeritas/decay/interactor/MuDecayInteractor.hh +++ b/src/celeritas/decay/interactor/MuDecayInteractor.hh @@ -14,7 +14,6 @@ #include "corecel/random/distribution/RejectionSampler.hh" #include "geocel/random/IsotropicDistribution.hh" #include "celeritas/Quantities.hh" -#include "celeritas/decay/data/MuDecayData.hh" #include "celeritas/phys/FourVector.hh" #include "celeritas/phys/Interaction.hh" #include "celeritas/phys/ParticleTrackView.hh" @@ -64,7 +63,7 @@ class MuDecayInteractor public: // Construct with shared and state data inline CELER_FUNCTION - MuDecayInteractor(MuDecayData const& shared, + MuDecayInteractor(Span daughters, ParticleTrackView const& particle, Real3 const& inc_direction, StackAllocator& allocate); @@ -82,14 +81,12 @@ class MuDecayInteractor //// DATA //// - // Constant data - MuDecayData const& shared_; + // Daughter particle properties + ParticleView daughter_; // Incident muon energy Energy const inc_energy_; // Allocate space for secondary particles (electron only) StackAllocator& allocate_; - // Define decay channel based on muon or anti-muon primary - ParticleId sec_id_; // Incident muon four vector FourVector inc_fourvec_; // Maximum electron energy [MeV] @@ -121,23 +118,18 @@ class MuDecayInteractor * using the physics manual definition. */ CELER_FUNCTION -MuDecayInteractor::MuDecayInteractor(MuDecayData const& shared, +MuDecayInteractor::MuDecayInteractor(Span daughters, ParticleTrackView const& particle, Real3 const& inc_direction, StackAllocator& allocate) - : shared_(shared) + : daughter_(particle.particle_view(daughters.front())) , inc_energy_(particle.energy()) , allocate_(allocate) - , sec_id_((particle.particle_id() == shared_.mu_minus_id) - ? shared_.electron_id - : shared_.positron_id) , inc_fourvec_{FourVector::from_particle(particle, inc_direction)} - , max_energy_(real_type{0.5} * shared_.muon_mass.value() - - shared_.electron_mass.value()) + , max_energy_(real_type{0.5} * value_as(particle.mass()) + - value_as(daughter_.mass())) { - CELER_EXPECT(shared_); - CELER_EXPECT(particle.particle_id() == shared_.mu_minus_id - || particle.particle_id() == shared_.mu_plus_id); + CELER_EXPECT(daughters.size() == 1); } //---------------------------------------------------------------------------// @@ -163,26 +155,26 @@ CELER_FUNCTION Interaction MuDecayInteractor::operator()(Engine& rng) { electron_nu_energy_frac = generate_canonical(rng); } while (RejectionSampler( - electron_nu_energy_frac * (real_type{1} - electron_nu_energy_frac), + electron_nu_energy_frac * (1 - electron_nu_energy_frac), real_type{0.25})(rng)); electron_energy_frac = generate_canonical(rng); - } while (electron_nu_energy_frac + electron_energy_frac < real_type{1}); + } while (electron_nu_energy_frac + electron_energy_frac < 1); // Decay isotropically in rest frame and boost secondaries to the lab frame auto charged_lep_fv = this->to_lab_frame( IsotropicDistribution{}(rng), - this->calc_momentum(electron_energy_frac, shared_.electron_mass), - shared_.electron_mass); + this->calc_momentum(electron_energy_frac, daughter_.mass()), + daughter_.mass()); // Return charged lepton only Interaction result = Interaction::from_absorption(); - result.secondaries = {secondaries, 1}; - result.secondaries[0].particle_id = sec_id_; + secondaries[0].particle_id = daughter_.particle_id(); // Interaction stores kinetic energy; FourVector stores total energy - result.secondaries[0].energy - = Energy{charged_lep_fv.energy - shared_.electron_mass.value()}; - result.secondaries[0].direction = make_unit_vector(charged_lep_fv.mom); + secondaries[0].energy + = Energy{charged_lep_fv.energy - value_as(daughter_.mass())}; + secondaries[0].direction = make_unit_vector(charged_lep_fv.mom); + result.secondaries = {secondaries, 1}; return result; } diff --git a/src/celeritas/em/process/BremsstrahlungProcess.hh b/src/celeritas/em/process/BremsstrahlungProcess.hh index bfbc39b503..b0c396b6f9 100644 --- a/src/celeritas/em/process/BremsstrahlungProcess.hh +++ b/src/celeritas/em/process/BremsstrahlungProcess.hh @@ -21,7 +21,7 @@ namespace celeritas /*! * Bremsstrahlung process for electrons and positrons. */ -class BremsstrahlungProcess : public Process +class BremsstrahlungProcess : public InteractionProcess { public: //!@{ diff --git a/src/celeritas/em/process/ComptonProcess.hh b/src/celeritas/em/process/ComptonProcess.hh index dacf7fd361..77344d1c0f 100644 --- a/src/celeritas/em/process/ComptonProcess.hh +++ b/src/celeritas/em/process/ComptonProcess.hh @@ -19,7 +19,7 @@ namespace celeritas /*! * Compton scattering process for gammas. */ -class ComptonProcess : public Process +class ComptonProcess : public InteractionProcess { public: //!@{ diff --git a/src/celeritas/em/process/CoulombScatteringProcess.hh b/src/celeritas/em/process/CoulombScatteringProcess.hh index 0705fbce65..fb6d758cbf 100644 --- a/src/celeritas/em/process/CoulombScatteringProcess.hh +++ b/src/celeritas/em/process/CoulombScatteringProcess.hh @@ -22,7 +22,7 @@ namespace celeritas /*! * Coulomb scattering process for electrons off of atoms. */ -class CoulombScatteringProcess : public Process +class CoulombScatteringProcess : public InteractionProcess { public: //!@{ diff --git a/src/celeritas/em/process/EIonizationProcess.hh b/src/celeritas/em/process/EIonizationProcess.hh index 8233292d4c..4f438d054e 100644 --- a/src/celeritas/em/process/EIonizationProcess.hh +++ b/src/celeritas/em/process/EIonizationProcess.hh @@ -19,7 +19,7 @@ namespace celeritas /*! * Ionization process for electrons and positrons. */ -class EIonizationProcess : public Process +class EIonizationProcess : public InteractionProcess { public: //!@{ diff --git a/src/celeritas/em/process/EPlusAnnihilationProcess.hh b/src/celeritas/em/process/EPlusAnnihilationProcess.hh index 53f1f3d1df..fabc0e38c3 100644 --- a/src/celeritas/em/process/EPlusAnnihilationProcess.hh +++ b/src/celeritas/em/process/EPlusAnnihilationProcess.hh @@ -19,7 +19,7 @@ namespace celeritas /*! * Annihiliation process for positrons. */ -class EPlusAnnihilationProcess final : public Process +class EPlusAnnihilationProcess final : public InteractionProcess { public: //!@{ diff --git a/src/celeritas/em/process/GammaConversionProcess.hh b/src/celeritas/em/process/GammaConversionProcess.hh index bdef82129c..ccb1012880 100644 --- a/src/celeritas/em/process/GammaConversionProcess.hh +++ b/src/celeritas/em/process/GammaConversionProcess.hh @@ -19,7 +19,7 @@ namespace celeritas /*! * Conversion of gammas to electrons and positrons. */ -class GammaConversionProcess : public Process +class GammaConversionProcess : public InteractionProcess { public: //!@{ diff --git a/src/celeritas/em/process/MuBremsstrahlungProcess.hh b/src/celeritas/em/process/MuBremsstrahlungProcess.hh index 2444b2ab71..4353d9bd83 100644 --- a/src/celeritas/em/process/MuBremsstrahlungProcess.hh +++ b/src/celeritas/em/process/MuBremsstrahlungProcess.hh @@ -20,7 +20,7 @@ namespace celeritas /*! * Bremsstrahlung process for muons. */ -class MuBremsstrahlungProcess : public Process +class MuBremsstrahlungProcess : public InteractionProcess { public: //!@{ diff --git a/src/celeritas/em/process/MuIonizationProcess.hh b/src/celeritas/em/process/MuIonizationProcess.hh index 1e81aa4152..33981bfbe5 100644 --- a/src/celeritas/em/process/MuIonizationProcess.hh +++ b/src/celeritas/em/process/MuIonizationProcess.hh @@ -19,7 +19,7 @@ namespace celeritas /*! * Ionization process for muons. */ -class MuIonizationProcess : public Process +class MuIonizationProcess : public InteractionProcess { public: //!@{ diff --git a/src/celeritas/em/process/MuPairProductionProcess.hh b/src/celeritas/em/process/MuPairProductionProcess.hh index e8cd7e9a10..dbd5e3fec7 100644 --- a/src/celeritas/em/process/MuPairProductionProcess.hh +++ b/src/celeritas/em/process/MuPairProductionProcess.hh @@ -21,7 +21,7 @@ namespace celeritas /*! * Electron-positron pair production process for muons. */ -class MuPairProductionProcess : public Process +class MuPairProductionProcess : public InteractionProcess { public: //!@{ diff --git a/src/celeritas/em/process/PhotoelectricProcess.hh b/src/celeritas/em/process/PhotoelectricProcess.hh index e416cc7106..d5ddff0ff9 100644 --- a/src/celeritas/em/process/PhotoelectricProcess.hh +++ b/src/celeritas/em/process/PhotoelectricProcess.hh @@ -22,7 +22,7 @@ namespace celeritas /*! * Photoelectric effect process for gammas. */ -class PhotoelectricProcess : public Process +class PhotoelectricProcess : public InteractionProcess { public: //!@{ diff --git a/src/celeritas/em/process/RayleighProcess.hh b/src/celeritas/em/process/RayleighProcess.hh index 05df11c087..989615d9a5 100644 --- a/src/celeritas/em/process/RayleighProcess.hh +++ b/src/celeritas/em/process/RayleighProcess.hh @@ -20,7 +20,7 @@ namespace celeritas /*! * Rayleigh scattering process for gammas. */ -class RayleighProcess : public Process +class RayleighProcess : public InteractionProcess { public: //!@{ diff --git a/src/celeritas/ext/GeantImporter.cc b/src/celeritas/ext/GeantImporter.cc index c06fcb9c6a..ecd5b0698d 100644 --- a/src/celeritas/ext/GeantImporter.cc +++ b/src/celeritas/ext/GeantImporter.cc @@ -18,6 +18,8 @@ #include #include #include +#include +#include #include #include #include @@ -29,6 +31,7 @@ #include #include #include +#include #include #include #include @@ -54,6 +57,7 @@ #include #include #include +#include #include #include #include @@ -171,6 +175,8 @@ struct ProcessFilter return (which & DataSelection::em); case G4ProcessType::fOptical: return (which & DataSelection::optical); + case G4ProcessType::fDecay: + return (which & (DataSelection::em | DataSelection::hadron)); case G4ProcessType::fHadronic: return (which & DataSelection::hadron); default: @@ -388,6 +394,25 @@ to_form_factor_type(G4NuclearFormfactorType const& form_factor_type) CELER_ASSERT_UNREACHABLE(); } +//---------------------------------------------------------------------------// +/*! + * Get a \c DecayChannelType from a \c G4VDecayChannel subclass. + */ +DecayChannelType to_decay_channel_type(G4VDecayChannel const* channel) +{ + CELER_EXPECT(channel); + + DecayChannelType result{DecayChannelType::size_}; + if (dynamic_cast(channel)) + { + result = DecayChannelType::muon; + } + CELER_VALIDATE(result != DecayChannelType::size_, + << "unsupported decay channel type '" + << channel->GetKinematicsName() << "'"); + return result; +} + //---------------------------------------------------------------------------// /*! * Return a populated \c inp::Particle vector. @@ -929,6 +954,50 @@ std::vector import_regions() return result; } +//---------------------------------------------------------------------------// +/*! + * Import a decay table from a Geant4 particle definition. + */ +inp::DecayPhysics::DecayTable import_decay_table(G4ParticleDefinition const& p) +{ + inp::DecayPhysics::DecayTable result; + + auto* g4_table = p.GetDecayTable(); + CELER_ASSERT(g4_table); + + for (auto channel_idx : range(g4_table->entries())) + { + auto* g4_channel = (*g4_table)[channel_idx]; + CELER_ASSERT(g4_channel); + + inp::DecayChannel channel; + channel.type = to_decay_channel_type(g4_channel); + channel.branching_ratio = g4_channel->GetBR(); + + // Exclude the neutrino daughters for muon decay + //! \todo Remove once neutrinos are included in \c MuDecayInteractor + size_type num_daughters = channel.type == DecayChannelType::muon + ? 1 + : g4_channel->GetNumberOfDaughters(); + + for (auto daughter_idx : range(num_daughters)) + { + auto const* daughter = g4_channel->GetDaughter(daughter_idx); + CELER_VALIDATE(daughter, + << "decay process is enabled for particle " + << p.GetParticleName() << ", but daughter '" + << g4_channel->GetDaughterName(daughter_idx) + << "' is not defined"); + channel.daughters.push_back(PDGNumber(daughter->GetPDGEncoding())); + } + result.push_back(std::move(channel)); + } + + CELER_LOG(debug) << "Imported decay table for particle " + << p.GetParticleName(); + return result; +} + //---------------------------------------------------------------------------// /*! * Return a populated \c ImportProcess vector. @@ -947,6 +1016,7 @@ auto import_processes(GeantImporter::DataSelection selected, auto& processes = imported.processes; auto& msc_models = imported.msc_models; auto& optical_models = imported.optical_models; + auto& decay = imported.decay; static celeritas::TypeDemangler const demangle_process; std::unordered_map visited; @@ -956,6 +1026,17 @@ auto import_processes(GeantImporter::DataSelection selected, auto append_process = [&](G4ParticleDefinition const& particle, G4VProcess const& process) -> void { + if (dynamic_cast(&process)) + { + // Store decay table if decay process is enabled for this particle + CELER_LOG(debug) << "Saving process '" << process.GetProcessName() + << "' for particle " << particle.GetParticleName() + << " (" << particle.GetPDGEncoding() << ')'; + decay.tables.insert({PDGNumber(particle.GetPDGEncoding()), + import_decay_table(particle)}); + return; + } + // Check for duplicate processes auto [prev, inserted] = visited.insert({&process, &particle}); diff --git a/src/celeritas/ext/RootInterfaceLinkDef.h b/src/celeritas/ext/RootInterfaceLinkDef.h index f3b24b7cea..c5b2f9fdc9 100644 --- a/src/celeritas/ext/RootInterfaceLinkDef.h +++ b/src/celeritas/ext/RootInterfaceLinkDef.h @@ -45,6 +45,8 @@ // Input data +#pragma link C++ class celeritas::inp::DecayChannel+; +#pragma link C++ class celeritas::inp::DecayPhysics+; #pragma link C++ class celeritas::inp::FresnelReflection+; #pragma link C++ class celeritas::inp::GaussianRoughness+; #pragma link C++ class celeritas::inp::Grid+; diff --git a/src/celeritas/inp/Physics.hh b/src/celeritas/inp/Physics.hh index 89659979ec..7b692bf729 100644 --- a/src/celeritas/inp/Physics.hh +++ b/src/celeritas/inp/Physics.hh @@ -14,6 +14,7 @@ #include "corecel/io/Label.hh" #include "celeritas/Types.hh" #include "celeritas/phys/AtomicNumber.hh" +#include "celeritas/phys/PDGNumber.hh" #include "PhysicsProcess.hh" #include "ProcessBuilder.hh" @@ -95,6 +96,32 @@ struct OpticalPhysics } }; +//---------------------------------------------------------------------------// +/*! + * Branching ratio and daughters for a decay channel. + */ +struct DecayChannel +{ + DecayChannelType type{DecayChannelType::size_}; + double branching_ratio{}; + std::vector daughters; +}; + +//---------------------------------------------------------------------------// +/*! + * Decay processes and options. + */ +struct DecayPhysics +{ + using DecayTable = std::vector; + + //! Decay channels for particles for which decay is enabled + std::unordered_map tables; + + //! Whether the data are assigned + explicit operator bool() const { return !tables.empty(); } +}; + //---------------------------------------------------------------------------// /*! * Set up physics options. diff --git a/src/celeritas/inp/ProcessBuilder.hh b/src/celeritas/inp/ProcessBuilder.hh index 6a8c73bd57..2656cf4e84 100644 --- a/src/celeritas/inp/ProcessBuilder.hh +++ b/src/celeritas/inp/ProcessBuilder.hh @@ -18,7 +18,7 @@ namespace celeritas class ImportedProcesses; class MaterialParams; class ParticleParams; -class Process; +class InteractionProcess; namespace inp { @@ -38,7 +38,8 @@ struct ProcessBuilderInput //!@{ //! \name User builder type aliases using ProcessBuilderFunction - = std::function(ProcessBuilderInput const&)>; + = std::function( + ProcessBuilderInput const&)>; using ProcessBuilderMap = std::unordered_map; //!@} diff --git a/src/celeritas/io/ImportData.hh b/src/celeritas/io/ImportData.hh index 93bba504dd..f3d1de0136 100644 --- a/src/celeritas/io/ImportData.hh +++ b/src/celeritas/io/ImportData.hh @@ -113,7 +113,7 @@ struct ImportData // Physics groups inp::OpticalPhysics optical_physics; - + inp::DecayPhysics decay; //!@} }; diff --git a/src/celeritas/neutron/process/NeutronElasticProcess.hh b/src/celeritas/neutron/process/NeutronElasticProcess.hh index 6181244275..1186765e95 100644 --- a/src/celeritas/neutron/process/NeutronElasticProcess.hh +++ b/src/celeritas/neutron/process/NeutronElasticProcess.hh @@ -22,7 +22,7 @@ namespace celeritas /*! * Elastic scattering process for neutrons. */ -class NeutronElasticProcess : public Process +class NeutronElasticProcess : public InteractionProcess { public: //!@{ diff --git a/src/celeritas/neutron/process/NeutronInelasticProcess.hh b/src/celeritas/neutron/process/NeutronInelasticProcess.hh index 01e8e991d9..fa73435895 100644 --- a/src/celeritas/neutron/process/NeutronInelasticProcess.hh +++ b/src/celeritas/neutron/process/NeutronInelasticProcess.hh @@ -22,7 +22,7 @@ namespace celeritas /*! * Inelastic interaction process for neutrons. */ -class NeutronInelasticProcess : public Process +class NeutronInelasticProcess : public InteractionProcess { public: //!@{ diff --git a/src/celeritas/phys/ImportedProcessAdapter.hh b/src/celeritas/phys/ImportedProcessAdapter.hh index b47f582794..6797576cdb 100644 --- a/src/celeritas/phys/ImportedProcessAdapter.hh +++ b/src/celeritas/phys/ImportedProcessAdapter.hh @@ -77,8 +77,8 @@ class ImportedProcessAdapter //! \name Type aliases using SPConstImported = std::shared_ptr; using SPConstParticles = std::shared_ptr; - using XsGrid = Process::XsGrid; - using EnergyLossGrid = Process::EnergyLossGrid; + using XsGrid = InteractionProcess::XsGrid; + using EnergyLossGrid = InteractionProcess::EnergyLossGrid; using SpanConstPDG = Span; //!@} diff --git a/src/celeritas/phys/ParticleTrackView.hh b/src/celeritas/phys/ParticleTrackView.hh index c360562fc0..7a0da1ebf5 100644 --- a/src/celeritas/phys/ParticleTrackView.hh +++ b/src/celeritas/phys/ParticleTrackView.hh @@ -77,8 +77,12 @@ class ParticleTrackView //// STATIC PROPERTIES (requires an indirection) //// + // Get static particle properties for the current state CELER_FORCEINLINE_FUNCTION ParticleView particle_view() const; + // Get static particle properties for another particle type + CELER_FORCEINLINE_FUNCTION ParticleView particle_view(ParticleId) const; + // Rest mass [MeV / c^2] CELER_FORCEINLINE_FUNCTION Mass mass() const; @@ -117,6 +121,9 @@ class ParticleTrackView // Relativistic momentum squared [MeV^2 / c^2] inline CELER_FUNCTION MomentumSq momentum_sq() const; + // Decay length [len] + inline CELER_FUNCTION real_type decay_length() const; + private: ParamsRef const& params_; StateRef const& states_; @@ -219,6 +226,16 @@ CELER_FUNCTION ParticleView ParticleTrackView::particle_view() const return ParticleView(params_, states_.particle_id[track_slot_]); } +//---------------------------------------------------------------------------// +/*! + * Get static particle properties for another particle type. + */ +CELER_FUNCTION ParticleView ParticleTrackView::particle_view(ParticleId id) const +{ + CELER_EXPECT(id < params_.size()); + return ParticleView(params_, id); +} + //---------------------------------------------------------------------------// /*! * Rest mass [MeV / c^2]. @@ -405,5 +422,18 @@ CELER_FUNCTION auto ParticleTrackView::momentum() const -> Momentum return Momentum{std::sqrt(this->momentum_sq().value())}; } +//---------------------------------------------------------------------------// +/*! + * Decay length \f$ d = \beta \gamma c \tau \f$ [len]. + */ +CELER_FUNCTION real_type ParticleTrackView::decay_length() const +{ + CELER_EXPECT(this->mass() > zero_quantity()); + CELER_EXPECT(this->decay_constant() != constants::stable_decay_constant); + + return constants::c_light * value_as(this->momentum()) + / (value_as(this->mass()) * this->decay_constant()); +} + //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/celeritas/phys/PhysicsData.hh b/src/celeritas/phys/PhysicsData.hh index c0e27c7fa3..082e760828 100644 --- a/src/celeritas/phys/PhysicsData.hh +++ b/src/celeritas/phys/PhysicsData.hh @@ -12,6 +12,7 @@ #include "corecel/data/StackAllocatorData.hh" #include "celeritas/Quantities.hh" #include "celeritas/Types.hh" +#include "celeritas/decay/DecayData.hh" #include "celeritas/em/data/AtomicRelaxationData.hh" #include "celeritas/em/data/EPlusGGData.hh" #include "celeritas/em/data/LivermorePEData.hh" @@ -132,6 +133,7 @@ struct ProcessGroup UniformTable energy_loss; //!< Process-integrated energy loss UniformTable range; //!< Process-integrated range UniformTable inverse_range; //!< Inverse process-integrated range + DecayTableData decay; //!< Decay channels and branching ratios ParticleProcessId at_rest; //!< ID of the particle's at-rest process //! True if assigned and valid @@ -243,6 +245,10 @@ struct PhysicsParamsScalars ActionId::size_type model_to_action{}; //! Number of physics models ModelId::size_type num_models{}; + //! Number of decay channels + ModelId::size_type num_channels{}; + //! ID of the decay process + ProcessId decay; // User-configurable constants real_type min_eprime_over_e{}; //!< xi [unitless] @@ -300,7 +306,7 @@ struct PhysicsParamsScalars //! Indicate an interaction failed to allocate memory CELER_FORCEINLINE_FUNCTION ActionId failure_action() const { - return ActionId{model_to_action + num_models}; + return ActionId{model_to_action + num_models + num_channels}; } }; @@ -323,6 +329,8 @@ struct PhysicsParamsData using ParticleItems = Collection; template using ParticleModelItems = Collection; + template + using DecayChannelItems = Collection; //// DATA //// @@ -349,6 +357,12 @@ struct PhysicsParamsData Items pmodel_ids; Items process_ids; + // Decay table storage + Items daughters; + Items channel_ids; + DecayChannelItems channels; + DecayChannelItems actions; + // Backend storage Items reals; @@ -385,6 +399,11 @@ struct PhysicsParamsData pmodel_ids = other.pmodel_ids; process_ids = other.process_ids; + daughters = other.daughters; + channel_ids = other.channel_ids; + channels = other.channels; + actions = other.actions; + reals = other.reals; return *this; @@ -454,6 +473,7 @@ struct PhysicsStateData StateItems state; //!< Track state [track] StateItems msc_step; //!< Internal MSC data [track] + StateItems decay_channel; //!< Sampled channel [track] Items per_process_xs; //!< XS [track][particle process] @@ -478,6 +498,7 @@ struct PhysicsStateData CELER_EXPECT(other); state = other.state; msc_step = other.msc_step; + decay_channel = other.decay_channel; per_process_xs = other.per_process_xs; @@ -501,6 +522,10 @@ inline void resize(PhysicsStateData* state, CELER_EXPECT(params.scalars.max_particle_processes > 0); resize(&state->state, size); resize(&state->msc_step, size); + if (params.scalars.decay) + { + resize(&state->decay_channel, size); + } resize(&state->per_process_xs, size * params.scalars.max_particle_processes); resize(&state->relaxation, params.hardwired.relaxation, size); diff --git a/src/celeritas/phys/PhysicsParams.cc b/src/celeritas/phys/PhysicsParams.cc index b498400d18..e384d74ea1 100644 --- a/src/celeritas/phys/PhysicsParams.cc +++ b/src/celeritas/phys/PhysicsParams.cc @@ -52,6 +52,7 @@ #include "PhysicsData.hh" #include "Process.hh" +#include "detail/DecayTableInserter.hh" #include "detail/DiscreteSelectAction.hh" #include "detail/EnergyMaxXsCalculator.hh" #include "detail/PreStepAction.hh" @@ -87,6 +88,7 @@ bool ignore_particle(PDGNumber pdg) */ PhysicsParams::PhysicsParams(Input inp) : processes_(std::move(inp.processes)) + , decay_(std::move(inp.decay)) , relaxation_(std::move(inp.relaxation)) { CELER_EXPECT(!processes_.empty()); @@ -136,6 +138,9 @@ PhysicsParams::PhysicsParams(Input inp) // Emit models for associated processes models_ = this->build_models(inp.action_registry); + // Emit channels for the decay process + channels_ = this->build_channels(inp.action_registry); + // Place "failure" *after* all the model IDs auto failure_action = make_shared( action_reg.next_id(), @@ -149,7 +154,7 @@ PhysicsParams::PhysicsParams(Input inp) HostValue host_data; this->build_options(inp.options, &host_data); this->build_ids(*inp.particles, &host_data); - this->build_tables(inp.options, *inp.materials, &host_data); + this->build_tables(inp, &host_data); this->build_model_tables(*inp.materials, &host_data); // Add step limiter if being used (TODO: remove this hack from physics) @@ -205,6 +210,9 @@ auto PhysicsParams::processes(ParticleId id) const -> SpanConstProcessId //---------------------------------------------------------------------------// // HELPER FUNCTIONS //---------------------------------------------------------------------------// +/*! + * Construct models and add to the action registry. + */ auto PhysicsParams::build_models(ActionRegistry* mgr) const -> VecModel { VecModel models; @@ -231,6 +239,33 @@ auto PhysicsParams::build_models(ActionRegistry* mgr) const -> VecModel return models; } +//---------------------------------------------------------------------------// +/*! + * Construct decay channels and add to the action registry. + */ +auto PhysicsParams::build_channels(ActionRegistry* reg) const -> VecChannel +{ + if (!decay_) + { + // No decay process + return {}; + } + + // Construct channels, assigning each channel ID + auto id_iter = Process::ActionIdIter{reg->next_id()}; + auto result = decay_->build_channels(id_iter); + CELER_ASSERT(!result.empty()); + for (auto const& channel : result) + { + CELER_ASSERT(channel); + CELER_ASSERT(channel->action_id() == *id_iter++); + + // Add channel to action registry + reg->insert(channel); + } + return result; +} + //---------------------------------------------------------------------------// /*! * Construct on-device particle-dependent physics options. @@ -315,7 +350,7 @@ void PhysicsParams::build_ids(ParticleParams const& particles, data->scalars.model_to_action = this->model(ModelId{0})->action_id().get(); // Note: use map to keep ProcessId sorted - std::vector particle_models(particles.size()); + std::vector particle_to_proc(particles.size()); std::vector temp_model_ids; ParticleModelId::size_type pm_idx{0}; @@ -339,13 +374,26 @@ void PhysicsParams::build_ids(ParticleParams const& particles, << " MeV) to be less than upper energy limit (" << value_as(applic.upper) << " MeV) for model " << m.label()); - particle_models[applic.particle.get()][process_id].push_back( + particle_to_proc[applic.particle.get()][process_id].push_back( {value_as(applic.lower), value_as(applic.upper), ParticleModelId{pm_idx++}}); temp_model_ids.push_back(mid); } } + if (decay_) + { + // Store decay process ID + data->scalars.decay = ProcessId(processes_.size()); + for (auto pid : range(ParticleId(particles.size()))) + { + if (decay_->is_applicable(pid)) + { + // Add null models for decay process + particle_to_proc[pid.get()][data->scalars.decay].push_back({}); + } + } + } make_builder(&data->model_ids) .insert_back(temp_model_ids.begin(), temp_model_ids.end()); @@ -355,13 +403,14 @@ void PhysicsParams::build_ids(ParticleParams const& particles, CollectionBuilder pmodel_ids(&data->pmodel_ids); DedupeCollectionBuilder reals(&data->reals); - process_groups.reserve(particle_models.size()); + process_groups.reserve(particle_to_proc.size()); // Loop over particle IDs, set ProcessGroup ProcessId::size_type max_particle_processes = 0; for (auto par_id : range(ParticleId{particles.size()})) { - auto& process_to_models = particle_models[par_id.get()]; + // Get map of process ID to model range for this particle + auto& process_to_models = particle_to_proc[par_id.get()]; if (process_to_models.empty() && !ignore_particle(particles.id_to_pdg(par_id))) { @@ -380,8 +429,15 @@ void PhysicsParams::build_ids(ParticleParams const& particles, // Add process ID temp_processes.push_back(pid_models.first); + if (pid_models.second.empty()) + { + // Skip constructing model data for decay, which has no models + CELER_ASSERT(pid_models.first == data->scalars.decay); + temp_model_groups.push_back({}); + continue; + } + std::vector& models = pid_models.second; - CELER_ASSERT(!models.empty()); // Construct model data std::vector temp_energy_grid; @@ -431,6 +487,7 @@ void PhysicsParams::build_ids(ParticleParams const& particles, } data->scalars.max_particle_processes = max_particle_processes; data->scalars.num_models = this->num_models(); + data->scalars.num_channels = this->channels_.size(); // Assign hardwired models that do on-the-fly xs calculation for (auto model_id : range(ModelId(this->num_models()))) @@ -493,9 +550,7 @@ void PhysicsParams::build_hardwired() /*! * Construct cross section data. */ -void PhysicsParams::build_tables(Options const& opts, - MaterialParams const& mats, - HostValue* data) const +void PhysicsParams::build_tables(Input const& inp, HostValue* data) const { CELER_EXPECT(*data); @@ -510,6 +565,7 @@ void PhysicsParams::build_tables(Options const& opts, XsGridInserter insert_xs(&data->reals, &data->xs_grids); UniformGridInserter insert_uniform(&data->reals, &data->uniform_grids); + detail::DecayTableInserter insert_decay(inp.particles, channels_, *data); Applicability applic; for (auto particle_id : range(ParticleId(data->process_groups.size()))) @@ -531,93 +587,120 @@ void PhysicsParams::build_tables(Options const& opts, // Loop over per-particle processes for (auto pp_idx : range(process_ids.size())) { - // Get energy bounds for this process - auto energy_grid = data->reals[model_groups[pp_idx].energy]; - applic.lower = Energy{energy_grid.front()}; - applic.upper = Energy{energy_grid.back()}; - CELER_ASSERT(applic.lower < applic.upper); - - Process const& proc = *this->process(process_ids[pp_idx]); - // Grid IDs for each grid type, each material - std::vector macro_xs_ids(mats.size()); - std::vector energy_loss_ids(mats.size()); - std::vector range_ids(mats.size()); + std::vector macro_xs_ids(inp.materials->size()); + std::vector energy_loss_ids(inp.materials->size()); + std::vector range_ids(inp.materials->size()); std::vector inverse_range_ids; - - // Energy of maximum cross section for each material - detail::EnergyMaxXsCalculator calc_integral_xs(opts, proc); std::vector energy_max_xs; - if (calc_integral_xs) - { - energy_max_xs.resize(mats.size()); - } - if (proc.applies_at_rest()) + auto set_at_rest = [&](Process const& p) { + if (p.applies_at_rest()) + { + /* \todo For now assume only one process per particle + * applies at rest. If a particle has multiple + * at-rest processes, we will need to select the + * process with the shortest time to interaction in + * \c select_discrete_interaction. + */ + CELER_VALIDATE(!process_group.at_rest, + << "particle " + << inp.particles->id_to_label(particle_id) + << " has multiple at-rest processes"); + + // Discrete interaction can occur at rest + process_group.at_rest = ParticleProcessId(pp_idx); + } + }; + + ProcessId process_id = process_ids[pp_idx]; + if (process_id == data->scalars.decay) { - /* \todo For now assume only one process per particle applies - * at rest. If a particle has multiple at-rest processes, we - * will need to check which process has the shortest time to - * interaction and choose that process in \c - * select_discrete_interaction. - */ - CELER_VALIDATE(!process_group.at_rest, - << "particle ID " << particle_id.get() - << " has multiple at-rest processes"); - - // Discrete interaction can occur at rest - process_group.at_rest = ParticleProcessId(pp_idx); - } + // Check whether discrete interaction can occur at rest + CELER_ASSERT(decay_); + set_at_rest(*decay_); - // Loop over materials - for (auto mat_idx : range(mats.size())) + // Build the decay table for this particle + auto table = decay_->decay_table(particle_id); + process_group.decay = insert_decay(table); + } + else { - applic.material = PhysMatId(mat_idx); + // Build the cross section and energy loss grids + InteractionProcess const& proc = *this->process(process_id); + + // Check whether discrete interaction can occur at rest + set_at_rest(proc); + + // Get energy bounds for this process + auto energy_grid = data->reals[model_groups[pp_idx].energy]; + applic.lower = Energy{energy_grid.front()}; + applic.upper = Energy{energy_grid.back()}; + CELER_ASSERT(applic.lower < applic.upper); - // Construct macroscopic cross section grid - auto macro_xs = proc.macro_xs(applic); - if (macro_xs) + // Energy of maximum cross section for each material + detail::EnergyMaxXsCalculator calc_integral_xs(inp.options, + proc); + if (calc_integral_xs) { - macro_xs_ids[mat_idx] = insert_xs(macro_xs); + energy_max_xs.resize(inp.materials->size()); } - // Construct energy loss grid - auto energy_loss = proc.energy_loss(applic); - if (energy_loss) + // Loop over materials + for (auto mat_idx : range(inp.materials->size())) { - energy_loss_ids[mat_idx] = insert_uniform(energy_loss); + applic.material = PhysMatId(mat_idx); - // Construct range grid from energy loss - auto range = RangeGridCalculator(BC::geant)(energy_loss); - range_ids[mat_idx] = insert_uniform(range); + // Construct macroscopic cross section grid + auto macro_xs = proc.macro_xs(applic); + if (macro_xs) + { + macro_xs_ids[mat_idx] = insert_xs(macro_xs); + } - if (range.interpolation.type - == InterpolationType::cubic_spline) + // Construct energy loss grid + auto energy_loss = proc.energy_loss(applic); + if (energy_loss) { - // Build the inverse range grid if cubic spline - // interpolation is used - inverse_range_ids.resize(mats.size(), UniformGridId{}); - - // The range and energy values are not inverted on the - // grid, but the derivatives are calculated using the - // inverted grid. - auto inverse_range - = data->uniform_grids[range_ids[mat_idx]]; - auto derivative - = SplineDerivCalculator(BC::geant).calc_from_inverse( - inverse_range, make_const_ref(*data).reals); - inverse_range.derivative = reals.insert_back( - derivative.begin(), derivative.end()); - inverse_range_ids[mat_idx] - = uniform_grids.push_back(inverse_range); + energy_loss_ids[mat_idx] = insert_uniform(energy_loss); + + // Construct range grid from energy loss + auto range + = RangeGridCalculator(BC::geant)(energy_loss); + range_ids[mat_idx] = insert_uniform(range); + + if (range.interpolation.type + == InterpolationType::cubic_spline) + { + // Build the inverse range grid if cubic spline + // interpolation is used + inverse_range_ids.resize(inp.materials->size(), + UniformGridId{}); + + // The range and energy values are not inverted on + // the grid, but the derivatives are calculated + // using the inverted grid. + auto inverse_range + = data->uniform_grids[range_ids[mat_idx]]; + auto derivative + = SplineDerivCalculator(BC::geant) + .calc_from_inverse( + inverse_range, + make_const_ref(*data).reals); + inverse_range.derivative = reals.insert_back( + derivative.begin(), derivative.end()); + inverse_range_ids[mat_idx] + = uniform_grids.push_back(inverse_range); + } } - } - if (calc_integral_xs) - { - // Find and store the energy of the largest cross section - // for this material if the integral approach is used - energy_max_xs[mat_idx] = calc_integral_xs(macro_xs); + if (calc_integral_xs) + { + // Find and store the energy of the largest cross + // section for this material if the integral approach + // is used + energy_max_xs[mat_idx] = calc_integral_xs(macro_xs); + } } } diff --git a/src/celeritas/phys/PhysicsParams.hh b/src/celeritas/phys/PhysicsParams.hh index 976d3b48c6..b986187b99 100644 --- a/src/celeritas/phys/PhysicsParams.hh +++ b/src/celeritas/phys/PhysicsParams.hh @@ -20,6 +20,7 @@ #include "celeritas/Quantities.hh" #include "celeritas/Types.hh" #include "celeritas/Units.hh" +#include "celeritas/decay/DecayProcess.hh" #include "celeritas/global/ActionInterface.hh" #include "Model.hh" @@ -62,10 +63,12 @@ class PhysicsParams final : public ParamsDataInterface //! \name Type aliases using SPConstParticles = std::shared_ptr; using SPConstMaterials = std::shared_ptr; - using SPConstProcess = std::shared_ptr; + using SPConstProcess = std::shared_ptr; + using SPConstDecay = std::shared_ptr; using SPConstModel = std::shared_ptr; using SPConstRelaxation = std::shared_ptr; + using VecChannel = DecayProcess::VecChannel; using VecProcess = std::vector; using SpanConstProcessId = Span; using ActionIdRange = Range; @@ -78,6 +81,7 @@ class PhysicsParams final : public ParamsDataInterface SPConstParticles particles; SPConstMaterials materials; VecProcess processes; + SPConstDecay decay; //!< Optional decay process SPConstRelaxation relaxation; //!< Optional atomic relaxation ActionRegistry* action_registry = nullptr; @@ -111,6 +115,12 @@ class PhysicsParams final : public ParamsDataInterface // Get the process for the given model inline ProcessId process_id(ModelId id) const; + //! Get decay process + SPConstDecay decay_process() const { return decay_; } + + //! Get decay channels + VecChannel const& decay_channels() const { return channels_; } + // Get the action IDs for all models inline ActionIdRange model_actions() const; @@ -142,6 +152,8 @@ class PhysicsParams final : public ParamsDataInterface // Host metadata/access VecProcess processes_; VecModel models_; + SPConstDecay decay_; + VecChannel channels_; SPConstRelaxation relaxation_; // Host/device storage and reference @@ -152,10 +164,11 @@ class PhysicsParams final : public ParamsDataInterface private: VecModel build_models(ActionRegistry*) const; + VecChannel build_channels(ActionRegistry*) const; void build_options(Options const&, HostValue*) const; void build_particle_options(ParticleOptions const&, ParticleScalars*) const; void build_ids(ParticleParams const&, HostValue*) const; - void build_tables(Options const&, MaterialParams const&, HostValue*) const; + void build_tables(Input const& inp, HostValue*) const; void build_model_tables(MaterialParams const&, HostValue*) const; void build_hardwired(); }; diff --git a/src/celeritas/phys/PhysicsParamsOutput.cc b/src/celeritas/phys/PhysicsParamsOutput.cc index fb0daf83fb..9b2084a4fe 100644 --- a/src/celeritas/phys/PhysicsParamsOutput.cc +++ b/src/celeritas/phys/PhysicsParamsOutput.cc @@ -69,9 +69,24 @@ void PhysicsParamsOutput::output(JsonPimpl* j) const Process const& p = *physics_->process(id); label.push_back(p.label()); } + if (auto p = physics_->decay_process()) + { + label.push_back(p->label()); + } obj["processes"] = {{"label", std::move(label)}}; } + // Save decay channels + { + auto label = json::array(); + + for (auto const& channel : physics_->decay_channels()) + { + label.push_back(channel->label()); + } + obj["decay_channels"] = {{"label", std::move(label)}}; + } + // Save options { auto const& scalars = physics_->host_ref().scalars; diff --git a/src/celeritas/phys/PhysicsStepUtils.hh b/src/celeritas/phys/PhysicsStepUtils.hh index 916f776236..502d931a2e 100644 --- a/src/celeritas/phys/PhysicsStepUtils.hh +++ b/src/celeritas/phys/PhysicsStepUtils.hh @@ -116,7 +116,13 @@ calc_physics_step_limit(MaterialTrackView const& material, for (auto ppid : range(ParticleProcessId{physics.num_particle_processes()})) { real_type process_xs = 0; - if (auto const& process = physics.integral_xs_process(ppid)) + if (physics.process(ppid) == physics.scalars().decay) + { + // Calculate the macroscopic cross section as the inverse of the + // particle's decay length + process_xs = 1 / particle.decay_length(); + } + else if (auto const& process = physics.integral_xs_process(ppid)) { // If the integral approach is used and this particle has an energy // loss process, estimate the maximum cross section over the step @@ -314,6 +320,16 @@ select_discrete_interaction(MaterialView const& material, return physics.scalars().integral_rejection_action(); } + if (physics.process(ppid) == physics.scalars().decay) + { + // The decay process was selected: sample a decay channel according to + // the branching ratios and return the action corresponding to the + // channel type + auto channel_id = physics.select_decay_channel(rng); + pstep.decay_channel(channel_id); + return physics.channel_to_action(channel_id); + } + // Find the model that applies at the particle energy auto find_model = physics.make_model_finder(ppid); auto pmid = find_model(particle.energy()); diff --git a/src/celeritas/phys/PhysicsStepView.hh b/src/celeritas/phys/PhysicsStepView.hh index c8a05dd502..3e3adff02f 100644 --- a/src/celeritas/phys/PhysicsStepView.hh +++ b/src/celeritas/phys/PhysicsStepView.hh @@ -56,6 +56,9 @@ class PhysicsStepView // Set the sampled element inline CELER_FUNCTION void element(ElementComponentId); + // Set the sampled decay channel + inline CELER_FUNCTION void decay_channel(DecayChannelId); + // Save MSC step data inline CELER_FUNCTION void msc_step(MscStep const&); @@ -77,6 +80,9 @@ class PhysicsStepView // Sampled element for discrete interaction CELER_FORCEINLINE_FUNCTION ElementComponentId element() const; + // Sampled decay channel + CELER_FORCEINLINE_FUNCTION DecayChannelId decay_channel() const; + // Mutable access to MSC step data (TODO: hack) inline CELER_FUNCTION MscStep& msc_step(); @@ -148,6 +154,16 @@ CELER_FUNCTION void PhysicsStepView::element(ElementComponentId elcomp_id) this->state().element = elcomp_id; } +//---------------------------------------------------------------------------// +/*! + * Set the sampled decay channel. + */ +CELER_FUNCTION void PhysicsStepView::decay_channel(DecayChannelId channel_id) +{ + CELER_EXPECT(track_slot_ < states_.decay_channel.size()); + states_.decay_channel[track_slot_] = channel_id; +} + //---------------------------------------------------------------------------// /*! * Save MSC step limit data. @@ -222,6 +238,16 @@ CELER_FUNCTION ElementComponentId PhysicsStepView::element() const return this->state().element; } +//---------------------------------------------------------------------------// +/*! + * Sampled decay channel. + */ +CELER_FUNCTION DecayChannelId PhysicsStepView::decay_channel() const +{ + CELER_EXPECT(track_slot_ < states_.decay_channel.size()); + return states_.decay_channel[track_slot_]; +} + //---------------------------------------------------------------------------// /*! * Mutable access to MSC step data (TODO: hack) diff --git a/src/celeritas/phys/PhysicsTrackView.hh b/src/celeritas/phys/PhysicsTrackView.hh index 49019625f6..62c75b0c40 100644 --- a/src/celeritas/phys/PhysicsTrackView.hh +++ b/src/celeritas/phys/PhysicsTrackView.hh @@ -135,9 +135,17 @@ class PhysicsTrackView TabulatedElementSelector make_element_selector(UniformTableId, Energy) const; + // Sample a decay channel according to the branching ratios + template + inline CELER_FUNCTION DecayChannelId select_decay_channel(Engine&) const; + // ID of the particle's at-rest process inline CELER_FUNCTION ParticleProcessId at_rest_process() const; + // Get the daughters for a decay channel + inline CELER_FUNCTION Span + daughters(DecayChannelId) const; + //// PARAMETER DATA //// // Convert an action to a model ID for diagnostics, empty if not a model @@ -146,6 +154,9 @@ class PhysicsTrackView // Convert a selected model ID into a simulation action ID inline CELER_FUNCTION ActionId model_to_action(ModelId) const; + // Convert a selected decay channel ID into a simulation action ID + inline CELER_FUNCTION ActionId channel_to_action(DecayChannelId) const; + // Get the model ID corresponding to the given ParticleModelId inline CELER_FUNCTION ModelId model_id(ParticleModelId) const; @@ -603,6 +614,23 @@ PhysicsTrackView::make_element_selector(UniformTableId table_id, energy}; } +//---------------------------------------------------------------------------// +/*! + * Sample a decay channel according to the branching ratios. + */ +template +CELER_FUNCTION DecayChannelId +PhysicsTrackView::select_decay_channel(Engine& rng) const +{ + auto const& table = this->process_group().decay; + CELER_ASSERT(table); + auto branching_ratios = params_.reals[table.branching_ratios]; + auto idx = make_selector( + [&branching_ratios](size_type i) { return branching_ratios[i]; }, + branching_ratios.size())(rng); + return params_.channel_ids[table.channel_ids[idx]]; +} + //---------------------------------------------------------------------------// /*! * ID of the particle's at-rest process. @@ -615,6 +643,21 @@ CELER_FUNCTION ParticleProcessId PhysicsTrackView::at_rest_process() const return this->process_group().at_rest; } +//---------------------------------------------------------------------------// +/*! + * Get the daughters for a decay channel. + */ +CELER_FUNCTION Span +PhysicsTrackView::daughters(DecayChannelId channel_id) const +{ + CELER_EXPECT(channel_id < params_.channel_ids.size()); + auto const& table = this->process_group().decay; + CELER_ASSERT(table); + auto const& channel = params_.channels[channel_id]; + CELER_ASSERT(channel); + return params_.daughters[channel.daughters]; +} + //---------------------------------------------------------------------------// /*! * Convert an action to a model ID for diagnostics, false if not a model. @@ -643,6 +686,21 @@ CELER_FUNCTION ActionId PhysicsTrackView::model_to_action(ModelId model) const return ActionId{model.unchecked_get() + params_.scalars.model_to_action}; } +//---------------------------------------------------------------------------// +/*! + * Convert a sampled decay channel ID into a simulation action ID. + * + * Unlike the model-to-action ID conversion, there isn't a one-to-one + * correspondence between the channel ID and the action ID. Multiple channel + * IDs can map to a single channel action. + */ +CELER_FUNCTION ActionId +PhysicsTrackView::channel_to_action(DecayChannelId channel) const +{ + CELER_ASSERT(channel < params_.actions.size()); + return params_.actions[channel]; +} + //---------------------------------------------------------------------------// /*! * Get the model ID corresponding to the given ParticleModelId. diff --git a/src/celeritas/phys/Process.hh b/src/celeritas/phys/Process.hh index 96477854e9..d9962d0d42 100644 --- a/src/celeritas/phys/Process.hh +++ b/src/celeritas/phys/Process.hh @@ -19,6 +19,36 @@ namespace celeritas struct Applicability; class Model; +//---------------------------------------------------------------------------// +/*! + * An interface for physics processes. + */ +class Process +{ + public: + //!@{ + //! \name Type aliases + using ActionIdIter = RangeIter; + //!@} + + public: + // Virtual destructor for polymorphic deletion + virtual ~Process(); + + //! Whether the process applies when the particle is stopped + virtual bool applies_at_rest() const = 0; + + //! Name of the process + virtual std::string_view label() const = 0; + + protected: + //!@{ + //! Allow construction and assignment only through daughter classes + Process() = default; + CELER_DEFAULT_COPY_MOVE(Process); + //!@} +}; + //---------------------------------------------------------------------------// /*! * An interface/factory method for creating models. @@ -33,22 +63,18 @@ class Model; * Each process has an interaction ("post step doit") and may have both energy * loss and range limiters. */ -class Process +class InteractionProcess : virtual public Process { public: //!@{ //! \name Type aliases using SPConstModel = std::shared_ptr; using VecModel = std::vector; - using ActionIdIter = RangeIter; using XsGrid = inp::XsGrid; using EnergyLossGrid = inp::UniformGrid; //!@} public: - // Virtual destructor for polymorphic deletion - virtual ~Process(); - //! Construct the models associated with this process virtual VecModel build_models(ActionIdIter start_id) const = 0; @@ -60,19 +86,6 @@ class Process //! Whether the integral method can be used to sample interaction length virtual bool supports_integral_xs() const = 0; - - //! Whether the process applies when the particle is stopped - virtual bool applies_at_rest() const = 0; - - //! Name of the process - virtual std::string_view label() const = 0; - - protected: - //!@{ - //! Allow construction and assignment only through daughter classes - Process() = default; - CELER_DEFAULT_COPY_MOVE(Process); - //!@} }; //---------------------------------------------------------------------------// diff --git a/src/celeritas/phys/ProcessBuilder.hh b/src/celeritas/phys/ProcessBuilder.hh index cd882b47a6..94b128b327 100644 --- a/src/celeritas/phys/ProcessBuilder.hh +++ b/src/celeritas/phys/ProcessBuilder.hh @@ -58,7 +58,7 @@ class ProcessBuilder //!@{ //! \name Type aliases using IPC = ImportProcessClass; - using SPProcess = std::shared_ptr; + using SPProcess = std::shared_ptr; using SPConstParticle = std::shared_ptr; using SPConstMaterial = std::shared_ptr; using SPConstImported = std::shared_ptr; diff --git a/src/celeritas/phys/detail/DecayTableInserter.cc b/src/celeritas/phys/detail/DecayTableInserter.cc new file mode 100644 index 0000000000..585339e876 --- /dev/null +++ b/src/celeritas/phys/detail/DecayTableInserter.cc @@ -0,0 +1,130 @@ +//------------------------------ -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/phys/detail/DecayTableInserter.cc +//---------------------------------------------------------------------------// +#include "DecayTableInserter.hh" + +#include "corecel/io/Logger.hh" +#include "celeritas/decay/channel/MuDecayChannel.hh" +#include "celeritas/phys/ParticleParams.hh" + +namespace celeritas +{ +namespace detail +{ +//---------------------------------------------------------------------------// +/*! + * Construct with particles, decay channels and pointer to host data. + */ +DecayTableInserter::DecayTableInserter(SPConstParticles particles, + VecChannel const& channels, + Data& data) + : particles_(particles) + , reals_{&data.reals} + , daughters_{&data.daughters} + , channel_ids_{&data.channel_ids} + , channels_{&data.channels} + , actions_{&data.actions} +{ + CELER_EXPECT(particles_); + + // Build a mapping of channel type to action ID + for (auto const& channel : channels) + { + if (dynamic_cast(channel.get())) + { + channel_to_action_.insert({DCT::muon, channel->action_id()}); + } + else + { + CELER_NOT_IMPLEMENTED("Decay channels other than muon"); + } + } +} + +//---------------------------------------------------------------------------// +/*! + * Construct decay table for a single particle. + */ +DecayTableData DecayTableInserter::operator()(DecayTable const& inp) +{ + if (inp.empty()) + { + // No decay process for this particle + return {}; + } + + DecayTableData result; + + double accum_br{0}; + std::vector branching_ratios; + std::vector channel_ids; + std::vector channels; + std::vector actions; + + branching_ratios.reserve(inp.size()); + channel_ids.reserve(inp.size()); + channels.reserve(inp.size()); + actions.reserve(inp.size()); + + for (auto const& ch_inp : inp) + { + CELER_VALIDATE(ch_inp.type != DCT::size_, + << "invalid decay channel type"); + CELER_VALIDATE(ch_inp.branching_ratio > 0, + << "invalid branching_ratio=" << ch_inp.branching_ratio + << " (should be positive)"); + CELER_VALIDATE(!ch_inp.daughters.empty(), + << "decay channel must have daughters"); + + // Get the particle ID from the PDG and store the daughters + DecayChannelData channel; + std::vector daughters; + daughters.reserve(ch_inp.daughters.size()); + for (auto pdg : ch_inp.daughters) + { + daughters.push_back(particles_->find(pdg)); + } + channel.daughters + = daughters_.insert_back(daughters.begin(), daughters.end()); + channel_ids.push_back(channels_.push_back(channel)); + + // Store the branching ratio and find the channel's action ID + branching_ratios.push_back(ch_inp.branching_ratio); + actions.push_back(channel_to_action_[ch_inp.type]); + + // Calculate the sum of the branching ratios + accum_br += ch_inp.branching_ratio; + } + result.channel_ids + = channel_ids_.insert_back(channel_ids.begin(), channel_ids.end()); + + if (!(soft_equal(1, accum_br))) + { + CELER_LOG(warning) + << "branching ratios for decay channels should sum to 1 " + "but instead sum tp " + << accum_br; + } + double norm = 1 / accum_br; + for (auto& br : branching_ratios) + { + // Renormalize the branching ratios + br *= norm; + } + result.branching_ratios + = reals_.insert_back(branching_ratios.begin(), branching_ratios.end()); + + // Append to mapping from channel ID to action ID + actions_.insert_back(actions.begin(), actions.end()); + + CELER_ENSURE(channels_.size() == actions_.size()); + CELER_ENSURE(result); + return result; +} + +//---------------------------------------------------------------------------// +} // namespace detail +} // namespace celeritas diff --git a/src/celeritas/phys/detail/DecayTableInserter.hh b/src/celeritas/phys/detail/DecayTableInserter.hh new file mode 100644 index 0000000000..50950bfd22 --- /dev/null +++ b/src/celeritas/phys/detail/DecayTableInserter.hh @@ -0,0 +1,63 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/phys/detail/DecayTableInserter.hh +//---------------------------------------------------------------------------// +#pragma once + +#include +#include + +#include "corecel/data/CollectionBuilder.hh" +#include "celeritas/Types.hh" +#include "celeritas/decay/DecayData.hh" +#include "celeritas/decay/DecayProcess.hh" +#include "celeritas/decay/channel/DecayChannel.hh" +#include "celeritas/inp/Physics.hh" +#include "celeritas/phys/PhysicsData.hh" + +namespace celeritas +{ +class ParticleParams; + +namespace detail +{ +//---------------------------------------------------------------------------// +/*! + * Construct a decay table from input data. + */ +class DecayTableInserter +{ + public: + //!@{ + //! \name Type aliases + using Data = HostVal; + using DecayTable = inp::DecayPhysics::DecayTable; + using SPConstParticles = std::shared_ptr; + using VecChannel = DecayProcess::VecChannel; + //!@} + + public: + // Construct with particles, decay channels and pointer to host data + DecayTableInserter(SPConstParticles, VecChannel const&, Data&); + + // Construct decay table for a single particle + DecayTableData operator()(DecayTable const& inp); + + private: + using DCT = DecayChannelType; + using MapChannelAction = std::unordered_map; + + SPConstParticles particles_; + CollectionBuilder reals_; + CollectionBuilder daughters_; + CollectionBuilder channel_ids_; + CollectionBuilder channels_; + CollectionBuilder actions_; + MapChannelAction channel_to_action_; +}; + +//---------------------------------------------------------------------------// +} // namespace detail +} // namespace celeritas diff --git a/src/celeritas/phys/detail/EnergyMaxXsCalculator.cc b/src/celeritas/phys/detail/EnergyMaxXsCalculator.cc index e43c28d716..f05e148987 100644 --- a/src/celeritas/phys/detail/EnergyMaxXsCalculator.cc +++ b/src/celeritas/phys/detail/EnergyMaxXsCalculator.cc @@ -22,7 +22,7 @@ namespace detail * Construct with physics options and process. */ EnergyMaxXsCalculator::EnergyMaxXsCalculator(PhysicsOptions const& opts, - Process const& proc) + InteractionProcess const& proc) : use_integral_xs_(!opts.disable_integral_xs && proc.supports_integral_xs()) , is_annihilation_(dynamic_cast(&proc)) { diff --git a/src/celeritas/phys/detail/EnergyMaxXsCalculator.hh b/src/celeritas/phys/detail/EnergyMaxXsCalculator.hh index 847e1171fb..fffa0e3d16 100644 --- a/src/celeritas/phys/detail/EnergyMaxXsCalculator.hh +++ b/src/celeritas/phys/detail/EnergyMaxXsCalculator.hh @@ -25,7 +25,7 @@ class EnergyMaxXsCalculator { public: // Construct with physics options and process - EnergyMaxXsCalculator(PhysicsOptions const&, Process const&); + EnergyMaxXsCalculator(PhysicsOptions const&, InteractionProcess const&); // Calculate the energy of the maximum cross section in the grid real_type operator()(inp::XsGrid const&) const; diff --git a/src/celeritas/setup/Problem.cc b/src/celeritas/setup/Problem.cc index 3d319743f1..b107a60654 100644 --- a/src/celeritas/setup/Problem.cc +++ b/src/celeritas/setup/Problem.cc @@ -31,6 +31,7 @@ #include "celeritas/alongstep/AlongStepGeneralLinearAction.hh" #include "celeritas/alongstep/AlongStepRZMapFieldMscAction.hh" #include "celeritas/alongstep/AlongStepUniformMscAction.hh" +#include "celeritas/decay/DecayProcess.hh" #include "celeritas/em/params/UrbanMscParams.hh" #include "celeritas/em/params/WentzelOKVIParams.hh" #include "celeritas/ext/GeantPhysicsOptions.hh" @@ -99,7 +100,7 @@ auto build_physics_processes(inp::EmPhysics const& em, { // TODO: process builder should be deleted; instead it should get // p.physics.em or whatever - std::vector> result; + std::vector> result; ProcessBuilder build_process( imported, params.particle, params.material, em.user_processes); for (auto pc : ProcessBuilder::get_all_process_classes(imported.processes)) @@ -166,6 +167,11 @@ auto build_physics(inp::Problem const& p, // Build processes input.processes = build_physics_processes(p.physics.em, params, imported); + if (imported.decay) + { + input.decay + = std::make_shared(params.particle, imported.decay); + } return std::make_shared(std::move(input)); } diff --git a/test/celeritas/data/four-steel-slabs.root b/test/celeritas/data/four-steel-slabs.root index 239c31a907..b1c0e2c32c 100644 Binary files a/test/celeritas/data/four-steel-slabs.root and b/test/celeritas/data/four-steel-slabs.root differ diff --git a/test/celeritas/data/four-steel-slabs.root-dump.json b/test/celeritas/data/four-steel-slabs.root-dump.json index 9fb7735d2f..958b5e8082 100644 --- a/test/celeritas/data/four-steel-slabs.root-dump.json +++ b/test/celeritas/data/four-steel-slabs.root-dump.json @@ -53,7 +53,7 @@ "name" : "G4_STAINLESS-STEEL", "state" : 1, "temperature" : 293.15, - "number_density" : 86993489258991547580416, + "number_density" : 86993489258991530803200, "elements" : [{ "_typename" : "celeritas::ImportMatElemComponent", "element_id" : 0, @@ -356,7 +356,7 @@ "_typename" : "celeritas::EnumArray", "data_" : [-9.21034037197618, 4.60517018598809] }, - "y" : [0.0919755519795959, 128.588033594672], + "y" : [0.0919755519795958, 128.588033594672], "interpolation" : { "_typename" : "celeritas::inp::Interpolation", "type" : 0, @@ -574,5 +574,9 @@ "dielectric_metal" : [] } } + }, + "decay" : { + "_typename" : "celeritas::inp::DecayPhysics", + "tables" : [] } } diff --git a/test/celeritas/data/simple-cms.root b/test/celeritas/data/simple-cms.root index 92a67975e2..95fcc6e68a 100644 Binary files a/test/celeritas/data/simple-cms.root and b/test/celeritas/data/simple-cms.root differ diff --git a/test/celeritas/decay/MuDecay.test.cc b/test/celeritas/decay/MuDecay.test.cc index df8e50d92f..7771ec0210 100644 --- a/test/celeritas/decay/MuDecay.test.cc +++ b/test/celeritas/decay/MuDecay.test.cc @@ -23,18 +23,14 @@ class MuDecayInteractorTest : public InteractorHostTestBase void SetUp() override { auto const& params = *this->particle_params(); - data_.electron_id = params.find(pdg::electron()); - data_.positron_id = params.find(pdg::positron()); - data_.mu_minus_id = params.find(pdg::mu_minus()); - data_.mu_plus_id = params.find(pdg::mu_plus()); - data_.electron_mass = params.get(data_.electron_id).mass(); - data_.muon_mass = params.get(data_.mu_minus_id).mass(); - + mu_minus_daughters_.push_back(params.find(pdg::electron())); + mu_plus_daughters_.push_back(params.find(pdg::positron())); this->set_inc_direction({0, 0, 1}); } protected: - MuDecayData data_; + std::vector mu_minus_daughters_; + std::vector mu_plus_daughters_; struct TestResult { @@ -46,7 +42,7 @@ class MuDecayInteractorTest : public InteractorHostTestBase { this->resize_secondaries(num_samples); this->set_inc_particle(pdg::mu_minus(), energy); - MuDecayInteractor interact(data_, + MuDecayInteractor interact(make_span(mu_minus_daughters_), this->particle_track(), this->direction(), this->secondary_allocator()); @@ -75,13 +71,15 @@ TEST_F(MuDecayInteractorTest, basic) { auto const& params = *this->particle_params(); auto const at_rest = MevEnergy{0}; - auto const max_lepton_energy = real_type{0.5} * data_.muon_mass.value() - - data_.electron_mass.value(); + auto const max_lepton_energy + = real_type{0.5} + * params.get(params.find(pdg::mu_minus())).mass().value() + - params.get(params.find(pdg::electron())).mass().value(); // Anti-muon decay { this->set_inc_particle(pdg::mu_plus(), at_rest); - MuDecayInteractor interact(data_, + MuDecayInteractor interact(make_span(mu_plus_daughters_), this->particle_track(), this->direction(), this->secondary_allocator()); @@ -97,7 +95,7 @@ TEST_F(MuDecayInteractorTest, basic) // Muon decay { this->set_inc_particle(pdg::mu_minus(), at_rest); - MuDecayInteractor interact(data_, + MuDecayInteractor interact(make_span(mu_minus_daughters_), this->particle_track(), this->direction(), this->secondary_allocator()); diff --git a/test/celeritas/phys/MockProcess.hh b/test/celeritas/phys/MockProcess.hh index d5fb821b61..d86678524d 100644 --- a/test/celeritas/phys/MockProcess.hh +++ b/test/celeritas/phys/MockProcess.hh @@ -57,7 +57,7 @@ using MevPerCmLoss = RealQuantity; * The given applicability vector has one element per model that it will * create. Each model can have a different particle type and/or energy range. */ -class MockProcess : public Process +class MockProcess : public InteractionProcess { public: //!@{ diff --git a/test/celeritas/phys/Physics.test.cc b/test/celeritas/phys/Physics.test.cc index 9e19179128..618d8bc17d 100644 --- a/test/celeritas/phys/Physics.test.cc +++ b/test/celeritas/phys/Physics.test.cc @@ -149,7 +149,7 @@ TEST_F(PhysicsParamsTest, output) auto j = nlohmann::json::parse(to_string(out)); j["sizes"].erase("reals"); EXPECT_JSON_EQ( - R"json({"_category":"internal","_label":"physics","models":{"label":["mock-model-1","mock-model-2","mock-model-3","mock-model-4","mock-model-5","mock-model-6","mock-model-7","mock-model-8","mock-model-9","mock-model-10","mock-model-11"],"process_id":[0,0,1,2,2,2,3,3,4,4,5]},"options":{"fixed_step_limiter":0.0,"heavy.lowest_energy":[0.001,"MeV"],"heavy.max_step_over_range":0.2,"heavy.min_range":0.010000000000000002,"light.lowest_energy":[0.001,"MeV"],"light.max_step_over_range":0.2,"light.min_range":0.1,"linear_loss_limit":0.01,"min_eprime_over_e":0.8},"processes":{"label":["scattering","absorption","purrs","hisses","meows","barks"]},"sizes":{"integral_xs":8,"model_groups":8,"model_ids":11,"process_groups":5,"process_ids":8,"uniform_grid_ids":57,"uniform_grids":57,"uniform_tables":44,"xs_grid_ids":32,"xs_grids":32,"xs_tables":8}})json", + R"json({"_category":"internal","_label":"physics","decay_channels":{"label":[]},"models":{"label":["mock-model-1","mock-model-2","mock-model-3","mock-model-4","mock-model-5","mock-model-6","mock-model-7","mock-model-8","mock-model-9","mock-model-10","mock-model-11"],"process_id":[0,0,1,2,2,2,3,3,4,4,5]},"options":{"fixed_step_limiter":0.0,"heavy.lowest_energy":[0.001,"MeV"],"heavy.max_step_over_range":0.2,"heavy.min_range":0.010000000000000002,"light.lowest_energy":[0.001,"MeV"],"light.max_step_over_range":0.2,"light.min_range":0.1,"linear_loss_limit":0.01,"min_eprime_over_e":0.8},"processes":{"label":["scattering","absorption","purrs","hisses","meows","barks"]},"sizes":{"integral_xs":8,"model_groups":8,"model_ids":11,"process_groups":5,"process_ids":8,"uniform_grid_ids":57,"uniform_grids":57,"uniform_tables":44,"xs_grid_ids":32,"xs_grids":32,"xs_tables":8}})json", j.dump()); }