diff --git a/app/celer-export-geant.cc b/app/celer-export-geant.cc index 567af14440..a25bf78484 100644 --- a/app/celer-export-geant.cc +++ b/app/celer-export-geant.cc @@ -104,6 +104,7 @@ void run(std::string const& gdml_filename, { ImportDataTrimmer::Input options; options.mupp = true; + options.max_size = 16; ImportDataTrimmer trim(options); trim(imported); } diff --git a/doc/implementation/em-physics.rst b/doc/implementation/em-physics.rst index 4a9ad9b36a..3286b50585 100644 --- a/doc/implementation/em-physics.rst +++ b/doc/implementation/em-physics.rst @@ -61,6 +61,8 @@ The following table summarizes the EM processes and models in Celeritas. | | | Mu Bethe--Bloch | :cpp:class:`celeritas::MuHadIonizationInteractor` | 200 keV--100 TeV | | +---------------------+-----------------------------+-----------------------------------------------------+--------------------------+ | | Bremsstrahlung | Mu bremsstrahlung | :cpp:class:`celeritas::MuBremsstrahlungInteractor` | 0--100 TeV | + | +---------------------+-----------------------------+-----------------------------------------------------+--------------------------+ + | | Pair production | Mu pair production | :cpp:class:`celeritas::MuPairProductionInteractor` | 0.85 GeV--100 TeV | +----------------+---------------------+-----------------------------+-----------------------------------------------------+--------------------------+ | :math:`\mu^+` | Ionization | Bragg | :cpp:class:`celeritas::MuHadIonizationInteractor` | 0--200 keV | | + +-----------------------------+-----------------------------------------------------+--------------------------+ @@ -69,6 +71,8 @@ The following table summarizes the EM processes and models in Celeritas. | | | Mu Bethe--Bloch | :cpp:class:`celeritas::MuHadIonizationInteractor` | 200 keV--100 TeV | | +---------------------+-----------------------------+-----------------------------------------------------+--------------------------+ | | Bremsstrahlung | Mu bremsstrahlung | :cpp:class:`celeritas::MuBremsstrahlungInteractor` | 0--100 TeV | + | +---------------------+-----------------------------+-----------------------------------------------------+--------------------------+ + | | Pair production | Mu pair production | :cpp:class:`celeritas::MuPairProductionInteractor` | 0.85 GeV--100 TeV | +----------------+---------------------+-----------------------------+-----------------------------------------------------+--------------------------+ .. only:: latex @@ -120,6 +124,8 @@ The following table summarizes the EM processes and models in Celeritas. & & Mu Bethe--Bloch & \texttt{\scriptsize celeritas::MuHadIonizationInteractor} & 200 keV -- 100 TeV \\ \cline{2-5} & Bremsstrahlung & Mu bremsstrahlung & \texttt{\scriptsize celeritas::MuBremsstrahlungInteractor} & 0--100 TeV \\ + \cline{2-5} + & Pair production & Mu pair production & \texttt{\scriptsize celeritas::MuPairProductionInteractor} & 0.85 GeV--100 TeV \\ \hline \multirow{3}{*}{$\mu^+$} & \multirow{2}{*}{Ionization} & Bragg & \texttt{\scriptsize celeritas::MuHadIonizationInteractor} & 0--200 keV \\ \cline{3-5} @@ -128,6 +134,8 @@ The following table summarizes the EM processes and models in Celeritas. & & Mu Bethe--Bloch & \texttt{\scriptsize celeritas::MuHadIonizationInteractor} & 200 keV -- 100 TeV \\ \cline{2-5} & Bremsstrahlung & Mu bremsstrahlung & \texttt{\scriptsize celeritas::MuBremsstrahlungInteractor} & 0--100 TeV \\ + \cline{2-5} + & Pair production & Mu pair production & \texttt{\scriptsize celeritas::MuPairProductionInteractor} & 0.85 GeV--100 TeV \\ \bottomrule \end{tabular} \end{threeparttable} @@ -207,7 +215,7 @@ rejection sampling. Muon bremsstrahlung and pair production use a simple distribution to sample the exiting polar angles. -.. doxygenclass:: celeritas::MuBremsPPAngularDistribution +.. doxygenclass:: celeritas::MuAngularDistribution Photon scattering ----------------- @@ -221,6 +229,7 @@ Conversion/annihilation/photoelectric .. doxygenclass:: celeritas::BetheHeitlerInteractor .. doxygenclass:: celeritas::EPlusGGInteractor .. doxygenclass:: celeritas::LivermorePEInteractor +.. doxygenclass:: celeritas::MuPairProductionInteractor .. doxygenclass:: celeritas::AtomicRelaxation @@ -230,6 +239,11 @@ on the fly (as opposed to pre-tabulated cross sections). .. doxygenclass:: celeritas::EPlusGGMacroXsCalculator .. doxygenclass:: celeritas::LivermorePEMicroXsCalculator +The energy transfer for muon pair production is sampled using the inverse +transform method with tabulated CDFs. + +.. doxygenclass:: celeritas::MuPPEnergyDistribution + Coulomb scattering ------------------ diff --git a/src/celeritas/CMakeLists.txt b/src/celeritas/CMakeLists.txt index c7d653f0fb..6b8a7f9f56 100644 --- a/src/celeritas/CMakeLists.txt +++ b/src/celeritas/CMakeLists.txt @@ -30,6 +30,7 @@ list(APPEND SOURCES em/process/GammaConversionProcess.cc em/process/MuIonizationProcess.cc em/process/MuBremsstrahlungProcess.cc + em/process/MuPairProductionProcess.cc em/process/PhotoelectricProcess.cc em/process/RayleighProcess.cc ext/GeantOpticalPhysicsOptions.cc @@ -279,6 +280,7 @@ celeritas_polysource(em/model/LivermorePEModel) celeritas_polysource(em/model/MollerBhabhaModel) celeritas_polysource(em/model/MuBetheBlochModel) celeritas_polysource(em/model/MuBremsstrahlungModel) +celeritas_polysource(em/model/MuPairProductionModel) celeritas_polysource(em/model/RayleighModel) celeritas_polysource(em/model/RelativisticBremModel) celeritas_polysource(em/model/SeltzerBergerModel) diff --git a/src/celeritas/em/data/MuPairProductionData.hh b/src/celeritas/em/data/MuPairProductionData.hh new file mode 100644 index 0000000000..1c5d6511cf --- /dev/null +++ b/src/celeritas/em/data/MuPairProductionData.hh @@ -0,0 +1,125 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/em/data/MuPairProductionData.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/Macros.hh" +#include "corecel/Types.hh" +#include "corecel/grid/TwodGridData.hh" +#include "celeritas/Quantities.hh" +#include "celeritas/Types.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * IDs used by muon pair production. + */ +struct MuPairProductionIds +{ + ParticleId mu_minus; + ParticleId mu_plus; + ParticleId electron; + ParticleId positron; + + //! Whether the IDs are assigned + explicit CELER_FUNCTION operator bool() const + { + return mu_minus && mu_plus && electron && positron; + } +}; + +//---------------------------------------------------------------------------// +/*! + * Sampling table for electron-positron pair production by muons. + * + * The value grids are organized by atomic number, where the Z grid is a + * hardcoded, problem-independent set of atomic numbers equally spaced in \f$ + * \log Z \f$ that is linearly interpolated on. Each 2D grid is: + * - x: logarithm of the energy [MeV] of the incident charged particle + * - y: logarithm of the ratio of the energy transfer to the incident particle + * energy + * - value: CDF calculated from the differential cross section + */ +template +struct MuPairProductionTableData +{ + //// TYPES //// + + template + using Items = Collection; + + //// MEMBER DATA //// + + ItemRange logz_grid; + Items grids; + + // Backend data + Items reals; + + //// MEMBER FUNCTIONS //// + + //! Whether the data is assigned + explicit CELER_FUNCTION operator bool() const + { + return !reals.empty() && !logz_grid.empty() + && logz_grid.size() == grids.size(); + } + + //! Assign from another set of data + template + MuPairProductionTableData& + operator=(MuPairProductionTableData const& other) + { + CELER_EXPECT(other); + reals = other.reals; + logz_grid = other.logz_grid; + grids = other.grids; + return *this; + } +}; + +//---------------------------------------------------------------------------// +/*! + * Constant data for the muon pair production interactor. + */ +template +struct MuPairProductionData +{ + //// MEMBER DATA //// + + //! Particle IDs + MuPairProductionIds ids; + + //! Electron mass [MeV / c^2] + units::MevMass electron_mass; + + // Sampling table storage + MuPairProductionTableData table; + + //// MEMBER FUNCTIONS //// + + //! Whether all data are assigned and valid + explicit CELER_FUNCTION operator bool() const + { + return ids && electron_mass > zero_quantity() && table; + } + + //! Assign from another set of data + template + MuPairProductionData& operator=(MuPairProductionData const& other) + { + CELER_EXPECT(other); + ids = other.ids; + electron_mass = other.electron_mass; + table = other.table; + return *this; + } +}; + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/em/distribution/MuBremsPPAngularDistribution.hh b/src/celeritas/em/distribution/MuAngularDistribution.hh similarity index 82% rename from src/celeritas/em/distribution/MuBremsPPAngularDistribution.hh rename to src/celeritas/em/distribution/MuAngularDistribution.hh index 68ac8f87dc..e0819830b3 100644 --- a/src/celeritas/em/distribution/MuBremsPPAngularDistribution.hh +++ b/src/celeritas/em/distribution/MuAngularDistribution.hh @@ -3,7 +3,7 @@ // See the top-level COPYRIGHT file for details. // SPDX-License-Identifier: (Apache-2.0 OR MIT) //---------------------------------------------------------------------------// -//! \file celeritas/em/distribution/MuBremsPPAngularDistribution.hh +//! \file celeritas/em/distribution/MuAngularDistribution.hh //---------------------------------------------------------------------------// #pragma once @@ -42,7 +42,7 @@ namespace celeritas * G4ModifiedMephi class and documented in section 11.2.4 of the Geant4 Physics * Reference (release 11.2). */ -class MuBremsPPAngularDistribution +class MuAngularDistribution { public: //!@{ @@ -53,9 +53,8 @@ class MuBremsPPAngularDistribution public: // Construct with incident and secondary particle quantities - inline CELER_FUNCTION MuBremsPPAngularDistribution(Energy inc_energy, - Mass inc_mass, - Energy energy); + inline CELER_FUNCTION + MuAngularDistribution(Energy inc_energy, Mass inc_mass, Energy energy); // Sample the cosine of the polar angle of the secondary template @@ -75,9 +74,9 @@ class MuBremsPPAngularDistribution * Construct with incident and secondary particle. */ CELER_FUNCTION -MuBremsPPAngularDistribution::MuBremsPPAngularDistribution(Energy inc_energy, - Mass inc_mass, - Energy energy) +MuAngularDistribution::MuAngularDistribution(Energy inc_energy, + Mass inc_mass, + Energy energy) : gamma_(1 + value_as(inc_energy) / value_as(inc_mass)) { real_type r_max_sq = ipow<2>( @@ -93,7 +92,7 @@ MuBremsPPAngularDistribution::MuBremsPPAngularDistribution(Energy inc_energy, * Sample the cosine of the polar angle of the secondary. */ template -CELER_FUNCTION real_type MuBremsPPAngularDistribution::operator()(Engine& rng) +CELER_FUNCTION real_type MuAngularDistribution::operator()(Engine& rng) { real_type a = generate_canonical(rng) * a_over_xi_; return std::cos(std::sqrt(a / (1 - a)) / gamma_); diff --git a/src/celeritas/em/distribution/MuPPEnergyDistribution.hh b/src/celeritas/em/distribution/MuPPEnergyDistribution.hh new file mode 100644 index 0000000000..d738e9a1de --- /dev/null +++ b/src/celeritas/em/distribution/MuPPEnergyDistribution.hh @@ -0,0 +1,254 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/em/distribution/MuPPEnergyDistribution.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/Macros.hh" +#include "corecel/Types.hh" +#include "corecel/grid/FindInterp.hh" +#include "corecel/grid/Interpolator.hh" +#include "corecel/grid/NonuniformGrid.hh" +#include "corecel/grid/TwodGridCalculator.hh" +#include "corecel/math/Algorithms.hh" +#include "celeritas/Quantities.hh" +#include "celeritas/em/data/MuPairProductionData.hh" +#include "celeritas/grid/InverseCdfFinder.hh" +#include "celeritas/mat/ElementView.hh" +#include "celeritas/phys/CutoffView.hh" +#include "celeritas/phys/ParticleTrackView.hh" +#include "celeritas/random/distribution/UniformRealDistribution.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Sample the electron and positron energies for muon pair production. + * + * The energy transfer to the electron-positron pair is sampled using inverse + * transform sampling on a tabulated CDF. The CDF is calculated on a 2D grid, + * where the x-axis is the log of the incident muon energy and the y-axis is + * the log of the ratio of the energy transfer to the incident particle energy. + * Because the shape of the distribution depends only weakly on the atomic + * number, the CDF is calculated for a hardcoded set of points equally spaced + * in \f$ \log Z \f$ and linearly interpolated. + * + * The formula used for the differential cross section is valid when the + * maximum energy transfer to the electron-positron pair lies between \f$ + * \epsilon_{\text{min}} = 4 m \f$, where \f$ m \f$ is the electron mass, and + * \f[ + \epsilon_{\text{max}} = E + \frac{3 \sqrt{e}}{4} \mu Z^{1/3}), + * \f] + * where \f$ E = T + \mu \f$ is the total muon energy, \f$ \mu \f$ is the muon + * mass, \f$ e \f$ is Euler's number, and \f$ Z \f$ is the atomic number. + * + * The maximum energy partition between the electron and positron is calculated + * as + * \f[ + r_{\text{max}} = \left[1 - 6 \frac{\mu^2}{E (E - \epsilon)} \right] \sqrt{1 + - \epsilon_{\text{min}} / \epsilon}. + * \f] + * The partition \f$ r \f$ is then sampled uniformly in \f$ [-r_{\text{max}}, + * r_{\text{max}}) \f$. + */ +class MuPPEnergyDistribution +{ + public: + //!@{ + //! \name Type aliases + using Mass = units::MevMass; + using Energy = units::MevEnergy; + //!@} + + //! Sampled secondary energies + struct PairEnergy + { + Energy electron; + Energy positron; + }; + + public: + // Construct from shared and incident particle data + inline CELER_FUNCTION + MuPPEnergyDistribution(NativeCRef const& shared, + ParticleTrackView const& particle, + CutoffView const& cutoffs, + ElementView const& element); + + template + inline CELER_FUNCTION PairEnergy operator()(Engine& rng); + + //! Minimum energy of the electron-positron pair [MeV]. + CELER_FUNCTION Energy min_pair_energy() const + { + return Energy(min_pair_energy_); + } + + //! Maximum energy of the electron-positron pair [MeV]. + CELER_FUNCTION Energy max_pair_energy() const + { + return Energy(max_pair_energy_); + } + + private: + //// DATA //// + + // CDF table for sampling the pair energy + NativeCRef const& table_; + // Incident particle energy [MeV] + real_type inc_energy_; + // Incident particle total energy [MeV] + real_type total_energy_; + // Square of the muon mass + real_type inc_mass_sq_; + // Secondary mass + real_type electron_mass_; + // Minimum energy transfer to electron/positron pair [MeV] + real_type min_pair_energy_; + // Maximum energy transfer to electron/positron pair [MeV] + real_type max_pair_energy_; + // Minimum incident particle kinetic energy [MeV] + real_type min_energy_; + // Log Z grid interpolation for the target element + FindInterp logz_interp_; + // Coefficient for calculating the pair energy + real_type coeff_; + // Lower bound on the ratio of the pair energy to the incident energy + real_type y_min_; + // Upper bound on the ratio of the pair energy to the incident energy + real_type y_max_; + + //// HELPER FUNCTIONS //// + + // Sample the scaled energy and interpolate in log Z + template + inline CELER_FUNCTION real_type sample_scaled_energy(Engine& rng) const; + + // Calculate the scaled energy for a given Z grid and sampled CDF + inline CELER_FUNCTION real_type calc_scaled_energy(size_type z_idx, + real_type u) const; +}; + +//---------------------------------------------------------------------------// +// INLINE DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Construct from shared and incident particle data. + * + * The incident energy *must* be within the bounds of the sampling table data. + */ +CELER_FUNCTION +MuPPEnergyDistribution::MuPPEnergyDistribution( + NativeCRef const& shared, + ParticleTrackView const& particle, + CutoffView const& cutoffs, + ElementView const& element) + : table_(shared.table) + , inc_energy_(value_as(particle.energy())) + , total_energy_(value_as(particle.total_energy())) + , inc_mass_sq_(ipow<2>(value_as(particle.mass()))) + , electron_mass_(value_as(shared.electron_mass)) + , min_pair_energy_(4 * value_as(shared.electron_mass)) + , max_pair_energy_(inc_energy_ + + value_as(particle.mass()) + * (1 + - real_type(0.75) * std::sqrt(constants::euler) + * element.cbrt_z())) + , min_energy_(max(value_as(cutoffs.energy(shared.ids.positron)), + min_pair_energy_)) +{ + CELER_EXPECT(max_pair_energy_ > min_energy_); + + NonuniformGrid logz_grid(table_.logz_grid, table_.reals); + logz_interp_ = find_interp(logz_grid, element.log_z()); + + NonuniformGrid y_grid( + table_.grids[ItemId(logz_interp_.index)].y, table_.reals); + coeff_ = std::log(min_pair_energy_ / inc_energy_) / y_grid.front(); + + // Compute the bounds on the ratio of the pair energy to incident energy + y_min_ = std::log(min_energy_ / inc_energy_) / coeff_; + y_max_ = std::log(max_pair_energy_ / inc_energy_) / coeff_; + + // Check that the bounds are within the grid bounds + CELER_ASSERT(y_min_ >= y_grid.front() + || soft_equal(y_grid.front(), y_min_)); + CELER_ASSERT(y_max_ <= y_grid.back() || soft_equal(y_grid.back(), y_max_)); + y_min_ = max(y_min_, y_grid.front()); + y_max_ = min(y_max_, y_grid.back()); +} + +//---------------------------------------------------------------------------// +/*! + * Sample the exiting pair energy. + */ +template +CELER_FUNCTION auto MuPPEnergyDistribution::operator()(Engine& rng) + -> PairEnergy +{ + // Sample the energy transfer + real_type pair_energy + = inc_energy_ * std::exp(coeff_ * this->sample_scaled_energy(rng)); + CELER_ASSERT(pair_energy >= min_energy_ && pair_energy <= max_pair_energy_); + + // Sample the energy partition between the electron and positron + real_type r_max = (1 + - 6 * inc_mass_sq_ + / (total_energy_ * (total_energy_ - pair_energy))) + * std::sqrt(1 - min_pair_energy_ / pair_energy); + real_type r = UniformRealDistribution(-r_max, r_max)(rng); + + // Calculate the electron and positron energies + PairEnergy result; + real_type half_energy = pair_energy * real_type(0.5); + result.electron = Energy((1 - r) * half_energy - electron_mass_); + result.positron = Energy((1 + r) * half_energy - electron_mass_); + + CELER_ENSURE(result.electron > zero_quantity()); + CELER_ENSURE(result.positron > zero_quantity()); + return result; +} + +//---------------------------------------------------------------------------// +/*! + * Sample the scaled energy and interpolate in log Z. + */ +template +CELER_FUNCTION real_type +MuPPEnergyDistribution::sample_scaled_energy(Engine& rng) const +{ + real_type u = generate_canonical(rng); + LinearInterpolator interp_energy{ + {0, this->calc_scaled_energy(logz_interp_.index, u)}, + {1, this->calc_scaled_energy(logz_interp_.index + 1, u)}}; + return interp_energy(logz_interp_.fraction); +} + +//---------------------------------------------------------------------------// +/*! + * Calculate the scaled energy for a given Z grid and sampled CDF value. + */ +CELER_FUNCTION real_type +MuPPEnergyDistribution::calc_scaled_energy(size_type z_idx, real_type u) const +{ + CELER_EXPECT(z_idx < table_.grids.size()); + CELER_EXPECT(u >= 0 && u < 1); + + TwodGridData const& cdf_grid = table_.grids[ItemId(z_idx)]; + auto calc_cdf + = TwodGridCalculator(cdf_grid, table_.reals)(std::log(inc_energy_)); + + // Get the sampled CDF value between the y bounds + real_type cdf = LinearInterpolator{{0, calc_cdf(y_min_)}, + {1, calc_cdf(y_max_)}}(u); + + // Find the grid index of the sampled CDF value + return InverseCdfFinder(NonuniformGrid(cdf_grid.y, table_.reals), + std::move(calc_cdf))(cdf); +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/em/executor/MuPairProductionExecutor.hh b/src/celeritas/em/executor/MuPairProductionExecutor.hh new file mode 100644 index 0000000000..efe928dd0b --- /dev/null +++ b/src/celeritas/em/executor/MuPairProductionExecutor.hh @@ -0,0 +1,58 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/em/executor/MuPairProductionExecutor.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/Assert.hh" +#include "corecel/Macros.hh" +#include "celeritas/em/data/MuPairProductionData.hh" +#include "celeritas/em/interactor/MuPairProductionInteractor.hh" +#include "celeritas/geo/GeoTrackView.hh" +#include "celeritas/global/CoreTrackView.hh" +#include "celeritas/mat/MaterialTrackView.hh" +#include "celeritas/phys/Interaction.hh" +#include "celeritas/phys/PhysicsStepView.hh" +#include "celeritas/random/RngEngine.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +struct MuPairProductionExecutor +{ + inline CELER_FUNCTION Interaction + operator()(celeritas::CoreTrackView const& track); + + NativeCRef params; +}; + +//---------------------------------------------------------------------------// +/*! + * Sample muon pair production from the current track. + */ +CELER_FUNCTION Interaction +MuPairProductionExecutor::operator()(CoreTrackView const& track) +{ + auto cutoff = track.make_cutoff_view(); + auto particle = track.make_particle_view(); + auto elcomp_id = track.make_physics_step_view().element(); + CELER_ASSERT(elcomp_id); + auto element + = track.make_material_view().make_material_view().make_element_view( + elcomp_id); + auto allocate_secondaries + = track.make_physics_step_view().make_secondary_allocator(); + auto const& dir = track.make_geo_view().dir(); + + MuPairProductionInteractor interact( + params, particle, cutoff, element, dir, allocate_secondaries); + + auto rng = track.make_rng_engine(); + return interact(rng); +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/em/interactor/MuBremsstrahlungInteractor.hh b/src/celeritas/em/interactor/MuBremsstrahlungInteractor.hh index 053ab575a3..b081cab047 100644 --- a/src/celeritas/em/interactor/MuBremsstrahlungInteractor.hh +++ b/src/celeritas/em/interactor/MuBremsstrahlungInteractor.hh @@ -14,7 +14,7 @@ #include "celeritas/Constants.hh" #include "celeritas/Quantities.hh" #include "celeritas/em/data/MuBremsstrahlungData.hh" -#include "celeritas/em/distribution/MuBremsPPAngularDistribution.hh" +#include "celeritas/em/distribution/MuAngularDistribution.hh" #include "celeritas/em/xs/MuBremsDiffXsCalculator.hh" #include "celeritas/mat/ElementView.hh" #include "celeritas/mat/MaterialView.hh" @@ -144,7 +144,7 @@ CELER_FUNCTION Interaction MuBremsstrahlungInteractor::operator()(Engine& rng) } while (RejectionSampler{gamma_energy * calc_dcs_(Energy{gamma_energy}), envelope_}(rng)); - MuBremsPPAngularDistribution sample_costheta( + MuAngularDistribution sample_costheta( particle_.energy(), particle_.mass(), Energy{gamma_energy}); // Update kinematics of the final state and return this interaction diff --git a/src/celeritas/em/interactor/MuPairProductionInteractor.hh b/src/celeritas/em/interactor/MuPairProductionInteractor.hh new file mode 100644 index 0000000000..469c3c62f8 --- /dev/null +++ b/src/celeritas/em/interactor/MuPairProductionInteractor.hh @@ -0,0 +1,176 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/em/interactor/MuPairProductionInteractor.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/Macros.hh" +#include "corecel/Types.hh" +#include "corecel/data/StackAllocator.hh" +#include "corecel/math/Algorithms.hh" +#include "corecel/math/ArrayOperators.hh" +#include "corecel/math/ArrayUtils.hh" +#include "celeritas/Constants.hh" +#include "celeritas/Quantities.hh" +#include "celeritas/em/data/MuPairProductionData.hh" +#include "celeritas/em/distribution/MuAngularDistribution.hh" +#include "celeritas/em/distribution/MuPPEnergyDistribution.hh" +#include "celeritas/mat/ElementView.hh" +#include "celeritas/phys/CutoffView.hh" +#include "celeritas/phys/Interaction.hh" +#include "celeritas/phys/ParticleTrackView.hh" +#include "celeritas/phys/Secondary.hh" +#include "celeritas/random/distribution/UniformRealDistribution.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Perform electron-positron pair production by muons. + * + * \note This performs the same sampling routine as in Geant4's + * G4MuPairProductionModel and as documented in the Geant4 Physics Reference + * Manual (Release 11.1) section 11.3. + */ +class MuPairProductionInteractor +{ + public: + // Construct with shared and state data + inline CELER_FUNCTION + MuPairProductionInteractor(NativeCRef const& shared, + ParticleTrackView const& particle, + CutoffView const& cutoffs, + ElementView const& element, + Real3 const& inc_direction, + StackAllocator& allocate); + + // Sample an interaction with the given RNG + template + inline CELER_FUNCTION Interaction operator()(Engine& rng); + + private: + //// TYPES //// + + using Energy = units::MevEnergy; + using Mass = units::MevMass; + using Momentum = units::MevMomentum; + using UniformRealDist = UniformRealDistribution; + + //// DATA //// + + // Shared model data + NativeCRef const& shared_; + // Allocate space for the secondary particle + StackAllocator& allocate_; + // Incident direction + Real3 const& inc_direction_; + // Incident particle energy [MeV] + Energy inc_energy_; + // Incident particle mass + Mass inc_mass_; + // Incident particle momentum [MeV / c] + real_type inc_momentum_; + // Sample the azimuthal angle + UniformRealDist sample_phi_; + // Sampler for the electron-positron pair energy + MuPPEnergyDistribution sample_energy_; + + //// HELPER FUNCTIONS //// + + // Calculate the secondary particle momentum from the sampled energy + inline CELER_FUNCTION real_type calc_momentum(Energy) const; +}; + +//---------------------------------------------------------------------------// +// INLINE DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Construct with shared and state data. + */ +CELER_FUNCTION MuPairProductionInteractor::MuPairProductionInteractor( + NativeCRef const& shared, + ParticleTrackView const& particle, + CutoffView const& cutoffs, + ElementView const& element, + Real3 const& inc_direction, + StackAllocator& allocate) + : shared_(shared) + , allocate_(allocate) + , inc_direction_(inc_direction) + , inc_energy_(particle.energy()) + , inc_mass_(particle.mass()) + , inc_momentum_(value_as(particle.momentum())) + , sample_phi_(0, 2 * constants::pi) + , sample_energy_(shared, particle, cutoffs, element) +{ + CELER_EXPECT(particle.particle_id() == shared.ids.mu_minus + || particle.particle_id() == shared.ids.mu_plus); +} + +//---------------------------------------------------------------------------// +/*! + * Simulate electron-posiitron pair production by muons. + */ +template +CELER_FUNCTION Interaction MuPairProductionInteractor::operator()(Engine& rng) +{ + // Allocate secondary electron and positron + Secondary* secondaries = allocate_(2); + if (secondaries == nullptr) + { + // Failed to allocate space for a secondary + return Interaction::from_failure(); + } + + // Sample the electron and positron energies + auto energy = sample_energy_(rng); + Energy pair_energy = energy.electron + energy.positron; + + // Sample the secondary directions + MuAngularDistribution sample_costheta(inc_energy_, inc_mass_, pair_energy); + real_type phi = sample_phi_(rng); + + // Create the secondary electron + Secondary& electron = secondaries[0]; + electron.particle_id = shared_.ids.electron; + electron.energy = energy.electron; + electron.direction + = rotate(from_spherical(sample_costheta(rng), phi), inc_direction_); + + // Create the secondary electron + Secondary& positron = secondaries[1]; + positron.particle_id = shared_.ids.positron; + positron.energy = energy.positron; + positron.direction + = rotate(from_spherical(sample_costheta(rng), phi + constants::pi), + inc_direction_); + + // Construct interaction for change to the incident muon + Interaction result; + result.secondaries = {secondaries, 2}; + result.energy = inc_energy_ - pair_energy; + result.direction = make_unit_vector( + inc_momentum_ * inc_direction_ + - this->calc_momentum(electron.energy) * electron.direction + - this->calc_momentum(positron.energy) * positron.direction); + + return result; +} + +//---------------------------------------------------------------------------// +/*! + * Calculate the secondary particle momentum from the sampled energy. + */ +CELER_FUNCTION real_type +MuPairProductionInteractor::calc_momentum(Energy energy) const +{ + return std::sqrt(ipow<2>(value_as(energy)) + + 2 * value_as(shared_.electron_mass) + * value_as(energy)); +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/em/model/MuPairProductionModel.cc b/src/celeritas/em/model/MuPairProductionModel.cc new file mode 100644 index 0000000000..83921570db --- /dev/null +++ b/src/celeritas/em/model/MuPairProductionModel.cc @@ -0,0 +1,181 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/em/model/MuPairProductionModel.cc +//---------------------------------------------------------------------------// +#include "MuPairProductionModel.hh" + +#include +#include + +#include "corecel/Config.hh" + +#include "corecel/cont/Range.hh" +#include "corecel/data/Collection.hh" +#include "corecel/data/CollectionBuilder.hh" +#include "corecel/data/HyperslabIndexer.hh" +#include "corecel/sys/ScopedMem.hh" +#include "celeritas/em/executor/MuPairProductionExecutor.hh" +#include "celeritas/em/interactor/detail/PhysicsConstants.hh" +#include "celeritas/global/ActionLauncher.hh" +#include "celeritas/global/CoreParams.hh" +#include "celeritas/global/TrackExecutor.hh" +#include "celeritas/grid/TwodGridBuilder.hh" +#include "celeritas/io/ImportProcess.hh" +#include "celeritas/phys/InteractionApplier.hh" +#include "celeritas/phys/PDGNumber.hh" +#include "celeritas/phys/ParticleParams.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Construct from model ID and shared data. + */ +MuPairProductionModel::MuPairProductionModel( + ActionId id, + ParticleParams const& particles, + SPConstImported data, + ImportMuPairProductionTable const& imported) + : StaticConcreteAction( + id, "pair-prod-muon", "interact by e-/e+ pair production by muons") + , imported_(data, + particles, + ImportProcessClass::mu_pair_prod, + ImportModelClass::mu_pair_prod, + {pdg::mu_minus(), pdg::mu_plus()}) +{ + CELER_EXPECT(id); + + ScopedMem record_mem("MuPairProductionModel.construct"); + + HostVal host_data; + + // Save IDs + host_data.ids.mu_minus = particles.find(pdg::mu_minus()); + host_data.ids.mu_plus = particles.find(pdg::mu_plus()); + host_data.ids.electron = particles.find(pdg::electron()); + host_data.ids.positron = particles.find(pdg::positron()); + CELER_VALIDATE(host_data.ids, + << "missing particles (required for '" + << this->description() << "')"); + + // Save particle properties + host_data.electron_mass = particles.get(host_data.ids.electron).mass(); + + // Build sampling table + CELER_VALIDATE(imported, + << "sampling table (required for '" << this->description() + << "') is empty"); + this->build_table(imported, &host_data.table); + + // Move to mirrored data, copying to device + data_ = CollectionMirror{std::move(host_data)}; + + CELER_ENSURE(data_); +} + +//---------------------------------------------------------------------------// +/*! + * Particle types and energy ranges that this model applies to. + */ +auto MuPairProductionModel::applicability() const -> SetApplicability +{ + Applicability mu_minus_applic; + mu_minus_applic.particle = this->host_ref().ids.mu_minus; + mu_minus_applic.lower = zero_quantity(); + mu_minus_applic.upper = detail::high_energy_limit(); + + Applicability mu_plus_applic = mu_minus_applic; + mu_plus_applic.particle = this->host_ref().ids.mu_plus; + + return {mu_minus_applic, mu_plus_applic}; +} + +//---------------------------------------------------------------------------// +/*! + * Get the microscopic cross sections for the given particle and material. + */ +auto MuPairProductionModel::micro_xs(Applicability applic) const + -> MicroXsBuilders +{ + return imported_.micro_xs(std::move(applic)); +} + +//---------------------------------------------------------------------------// +/*! + * Interact with host data. + */ +void MuPairProductionModel::step(CoreParams const& params, + CoreStateHost& state) const +{ + auto execute = make_action_track_executor( + params.ptr(), + state.ptr(), + this->action_id(), + InteractionApplier{MuPairProductionExecutor{this->host_ref()}}); + return launch_action(*this, params, state, execute); +} + +//---------------------------------------------------------------------------// +#if !CELER_USE_DEVICE +void MuPairProductionModel::step(CoreParams const&, CoreStateDevice&) const +{ + CELER_NOT_CONFIGURED("CUDA OR HIP"); +} +#endif + +//---------------------------------------------------------------------------// +/*! + * Construct sampling table. + */ +void MuPairProductionModel::build_table( + ImportMuPairProductionTable const& imported, HostTable* table) const +{ + CELER_EXPECT(imported); + CELER_EXPECT(table); + + // Build 2D sampling table + TwodGridBuilder build_grid{&table->reals}; + CollectionBuilder grids{&table->grids}; + for (auto const& pvec : imported.physics_vectors) + { + CELER_VALIDATE(pvec, + << "invalid grid in sampling table for '" + << this->description() << "'"); + + Array dims{static_cast(pvec.x.size()), + static_cast(pvec.y.size())}; + HyperslabIndexer index(dims); + + // Normalize the CDF + std::vector cdf(pvec.value.size()); + for (size_type i : range(dims[0])) + { + double norm = 1 / pvec.value[index(i, dims[1] - 1)]; + for (size_type j : range(dims[1])) + { + cdf[index(i, j)] = pvec.value[index(i, j)] * norm; + } + } + grids.push_back( + build_grid(make_span(pvec.x), make_span(pvec.y), make_span(cdf))); + } + + // Build log Z grid + std::vector log_z; + log_z.reserve(imported.atomic_number.size()); + for (auto z : imported.atomic_number) + { + log_z.push_back(std::log(z)); + } + table->logz_grid + = make_builder(&table->reals).insert_back(log_z.begin(), log_z.end()); + + CELER_ENSURE(*table); +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/em/model/MuPairProductionModel.cu b/src/celeritas/em/model/MuPairProductionModel.cu new file mode 100644 index 0000000000..9d642dcffd --- /dev/null +++ b/src/celeritas/em/model/MuPairProductionModel.cu @@ -0,0 +1,36 @@ +//---------------------------------*-CUDA-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/em/model/MuPairProductionModel.cu +//---------------------------------------------------------------------------// +#include "MuPairProductionModel.hh" + +#include "celeritas/em/executor/MuPairProductionExecutor.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 MuPairProductionModel::step(CoreParams const& params, + CoreStateDevice& state) const +{ + auto execute = make_action_track_executor( + params.ptr(), + state.ptr(), + this->action_id(), + InteractionApplier{MuPairProductionExecutor{this->device_ref()}}); + static ActionLauncher const launch_kernel(*this); + launch_kernel(*this, params, state, execute); +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/em/model/MuPairProductionModel.hh b/src/celeritas/em/model/MuPairProductionModel.hh new file mode 100644 index 0000000000..ae4bf45dbe --- /dev/null +++ b/src/celeritas/em/model/MuPairProductionModel.hh @@ -0,0 +1,73 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/em/model/MuPairProductionModel.hh +//---------------------------------------------------------------------------// +#pragma once + +#include +#include + +#include "corecel/data/CollectionMirror.hh" +#include "celeritas/Quantities.hh" +#include "celeritas/em/data/MuPairProductionData.hh" +#include "celeritas/io/ImportMuPairProductionTable.hh" +#include "celeritas/phys/ImportedModelAdapter.hh" +#include "celeritas/phys/ImportedProcessAdapter.hh" +#include "celeritas/phys/Model.hh" + +namespace celeritas +{ +class ParticleParams; + +//---------------------------------------------------------------------------// +/*! + * Set up and launch the muon pair production model. + */ +class MuPairProductionModel final : public Model, public StaticConcreteAction +{ + public: + //!@{ + //! \name Type aliases + using HostRef = HostCRef; + using DeviceRef = DeviceCRef; + using SPConstImported = std::shared_ptr; + //!@} + + public: + // Construct from model ID and other necessary data + MuPairProductionModel(ActionId, + ParticleParams const&, + SPConstImported, + ImportMuPairProductionTable const&); + + // Particle types and energy ranges that this model applies to + SetApplicability applicability() const final; + + // Get the microscopic cross sections for the given particle and material + MicroXsBuilders micro_xs(Applicability) const final; + + // Apply the interaction kernel on device + void step(CoreParams const&, CoreStateHost&) const final; + + // Apply the interaction kernel + void step(CoreParams const&, CoreStateDevice&) const final; + + //!@{ + //! Access model data + HostRef const& host_ref() const { return data_.host_ref(); } + DeviceRef const& device_ref() const { return data_.device_ref(); } + //!@} + + private: + CollectionMirror data_; + ImportedModelAdapter imported_; + + using HostTable = HostVal; + void build_table(ImportMuPairProductionTable const&, HostTable*) const; +}; + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/em/process/MuPairProductionProcess.cc b/src/celeritas/em/process/MuPairProductionProcess.cc new file mode 100644 index 0000000000..98f8c781d8 --- /dev/null +++ b/src/celeritas/em/process/MuPairProductionProcess.cc @@ -0,0 +1,71 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/em/process/MuPairProductionProcess.cc +//---------------------------------------------------------------------------// +#include "MuPairProductionProcess.hh" + +#include + +#include "corecel/Assert.hh" +#include "corecel/cont/Range.hh" +#include "celeritas/em/model/MuPairProductionModel.hh" +#include "celeritas/io/ImportProcess.hh" +#include "celeritas/phys/PDGNumber.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Construct from host data. + */ +MuPairProductionProcess::MuPairProductionProcess(SPConstParticles particles, + SPConstImported process_data, + Options options, + SPConstImportTable table) + : particles_(std::move(particles)) + , imported_(process_data, + particles_, + ImportProcessClass::mu_pair_prod, + {pdg::mu_minus(), pdg::mu_plus()}) + , options_(options) + , table_(std::move(table)) +{ + CELER_EXPECT(particles_); + CELER_EXPECT(table_); +} + +//---------------------------------------------------------------------------// +/*! + * Construct the models associated with this process. + */ +auto MuPairProductionProcess::build_models(ActionIdIter start_id) const + -> VecModel +{ + return {std::make_shared( + *start_id++, *particles_, imported_.processes(), *table_)}; +} + +//---------------------------------------------------------------------------// +/*! + * Get the interaction cross sections for the given energy range. + */ +auto MuPairProductionProcess::step_limits(Applicability applic) const + -> StepLimitBuilders +{ + return imported_.step_limits(std::move(applic)); +} + +//---------------------------------------------------------------------------// +/*! + * Name of the process. + */ +std::string_view MuPairProductionProcess::label() const +{ + return "Muon electron-positron pair production"; +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/em/process/MuPairProductionProcess.hh b/src/celeritas/em/process/MuPairProductionProcess.hh new file mode 100644 index 0000000000..a889b35c01 --- /dev/null +++ b/src/celeritas/em/process/MuPairProductionProcess.hh @@ -0,0 +1,70 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/em/process/MuPairProductionProcess.hh +//---------------------------------------------------------------------------// +#pragma once + +#include + +#include "celeritas/io/ImportMuPairProductionTable.hh" +#include "celeritas/phys/Applicability.hh" +#include "celeritas/phys/AtomicNumber.hh" +#include "celeritas/phys/ImportedProcessAdapter.hh" +#include "celeritas/phys/ParticleParams.hh" +#include "celeritas/phys/Process.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Electron-positron pair production process for muons. + */ +class MuPairProductionProcess : public Process +{ + public: + //!@{ + //! \name Type aliases + using SPConstParticles = std::shared_ptr; + using SPConstImported = std::shared_ptr; + using SPConstImportTable + = std::shared_ptr; + //!@} + + // Options for the pair production process + struct Options + { + bool use_integral_xs{true}; //!> Use integral method for sampling + //! discrete interaction length + }; + + public: + // Construct from pair production data + MuPairProductionProcess(SPConstParticles particles, + SPConstImported process_data, + Options options, + SPConstImportTable table); + + // Construct the models associated with this process + VecModel build_models(ActionIdIter start_id) const final; + + // Get the interaction cross sections for the given energy range + StepLimitBuilders step_limits(Applicability range) const final; + + //! Whether to use the integral method to sample interaction length + bool use_integral_xs() const final { return options_.use_integral_xs; } + + // Name of the process + std::string_view label() const final; + + private: + SPConstParticles particles_; + ImportedProcessAdapter imported_; + Options options_; + SPConstImportTable table_; +}; + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/ext/detail/GeantProcessImporter.cc b/src/celeritas/ext/detail/GeantProcessImporter.cc index 7442bbd3a8..85ec6855ec 100644 --- a/src/celeritas/ext/detail/GeantProcessImporter.cc +++ b/src/celeritas/ext/detail/GeantProcessImporter.cc @@ -33,6 +33,7 @@ #include "corecel/Assert.hh" #include "corecel/cont/Range.hh" +#include "corecel/data/HyperslabIndexer.hh" #include "corecel/io/Logger.hh" #include "celeritas/UnitTypes.hh" #include "celeritas/io/ImportUnits.hh" @@ -428,6 +429,11 @@ import_physics_vector(G4PhysicsVector const& g4v, Array units) //---------------------------------------------------------------------------// /*! * Import a 2D physics vector. + * + * \note In Geant4 the values are stored as a vector of vectors indexed as + * [y][x]. Because the Celeritas \c TwodGridCalculator and \c + * TwodSubgridCalculator expect the y grid values to be on the inner dimension, + * the table is inverted during import so that the x and y grids are swapped. */ ImportPhysics2DVector import_physics_2dvector(G4Physics2DVector const& g4pv, Array units) @@ -437,18 +443,22 @@ ImportPhysics2DVector import_physics_2dvector(G4Physics2DVector const& g4pv, double const y_scaling = native_value_from_clhep(units[1]); double const v_scaling = native_value_from_clhep(units[2]); + Array dims{static_cast(g4pv.GetLengthY()), + static_cast(g4pv.GetLengthX())}; + HyperslabIndexer<2> index(dims); + ImportPhysics2DVector pv; - pv.x.resize(g4pv.GetLengthX()); - pv.y.resize(g4pv.GetLengthY()); - pv.value.resize(pv.x.size() * pv.y.size()); + pv.x.resize(dims[0]); + pv.y.resize(dims[1]); + pv.value.resize(dims[0] * dims[1]); - for (auto i : range(pv.x.size())) + for (auto i : range(dims[0])) { - pv.x[i] = g4pv.GetX(i) * x_scaling; - for (auto j : range(pv.y.size())) + pv.x[i] = g4pv.GetY(i) * y_scaling; + for (auto j : range(dims[1])) { - pv.y[j] = g4pv.GetY(j) * y_scaling; - pv.value[pv.y.size() * i + j] = g4pv.GetValue(i, j) * v_scaling; + pv.y[j] = g4pv.GetX(j) * x_scaling; + pv.value[index(i, j)] = g4pv.GetValue(j, i) * v_scaling; } } CELER_ENSURE(pv); diff --git a/src/celeritas/grid/InverseCdfFinder.hh b/src/celeritas/grid/InverseCdfFinder.hh index a0623edcdb..4ad3002acf 100644 --- a/src/celeritas/grid/InverseCdfFinder.hh +++ b/src/celeritas/grid/InverseCdfFinder.hh @@ -12,6 +12,7 @@ #include "corecel/cont/Range.hh" #include "corecel/grid/Interpolator.hh" #include "corecel/math/Algorithms.hh" +#include "corecel/math/SoftEqual.hh" namespace celeritas { @@ -52,7 +53,8 @@ CELER_FUNCTION InverseCdfFinder::InverseCdfFinder(G&& grid, C&& calc_cdf) , calc_cdf_(celeritas::forward(calc_cdf)) { CELER_EXPECT(grid_.size() >= 2); - CELER_EXPECT(calc_cdf_[0] == 0 && calc_cdf_[grid_.size() - 1] == 1); + CELER_EXPECT(calc_cdf_[0] == 0 + && soft_equal(calc_cdf_[grid_.size() - 1], real_type(1))); } //---------------------------------------------------------------------------// diff --git a/src/celeritas/io/ImportDataTrimmer.cc b/src/celeritas/io/ImportDataTrimmer.cc index c77c857a29..66362b6185 100644 --- a/src/celeritas/io/ImportDataTrimmer.cc +++ b/src/celeritas/io/ImportDataTrimmer.cc @@ -23,6 +23,9 @@ struct ImportDataTrimmer::GridFilterer // Whether to keep the data at this index inline bool operator()(size_type i) const; + + // Whether to trim the data + explicit operator bool() const { return stride > 0; }; }; //---------------------------------------------------------------------------// @@ -54,6 +57,10 @@ void ImportDataTrimmer::operator()(ImportData& data) (*this)(data.atomic_relaxation_data); (*this)(data.optical_materials); + + this->for_each(data.elements); + this->for_each(data.geo_materials); + this->for_each(data.phys_materials); } if (options_.physics) @@ -62,6 +69,16 @@ void ImportDataTrimmer::operator()(ImportData& data) (*this)(data.particles); (*this)(data.processes); (*this)(data.msc_models); + + this->for_each(data.processes); + this->for_each(data.msc_models); + this->for_each(data.sb_data); + this->for_each(data.livermore_pe_data); + this->for_each(data.neutron_elastic_data); + this->for_each(data.atomic_relaxation_data); + + this->for_each(data.optical_models); + this->for_each(data.optical_materials); } if (options_.mupp) @@ -69,20 +86,6 @@ void ImportDataTrimmer::operator()(ImportData& data) // Reduce the resolution of the muon pair production table (*this)(data.mu_pair_production_data); } - - this->for_each(data.elements); - this->for_each(data.geo_materials); - this->for_each(data.phys_materials); - - this->for_each(data.processes); - this->for_each(data.msc_models); - this->for_each(data.sb_data); - this->for_each(data.livermore_pe_data); - this->for_each(data.neutron_elastic_data); - this->for_each(data.atomic_relaxation_data); - - this->for_each(data.optical_models); - this->for_each(data.optical_materials); } //---------------------------------------------------------------------------// @@ -264,7 +267,7 @@ void ImportDataTrimmer::operator()(ImportPhysics2DVector& data) { for (auto j : range(y_filter.orig_size)) { - if (x_filter(i) && y_filter(j)) + if ((!x_filter || x_filter(i)) && (!y_filter || y_filter(j))) { new_value.push_back(*src); } @@ -272,6 +275,7 @@ void ImportDataTrimmer::operator()(ImportPhysics2DVector& data) } } CELER_ASSERT(src == data.value.cend()); + CELER_ASSERT(new_value.size() == data.x.size() * data.y.size()); data.value = std::move(new_value); @@ -285,7 +289,8 @@ void ImportDataTrimmer::operator()(ImportPhysics2DVector& data) template void ImportDataTrimmer::operator()(std::vector& data) { - if (options_.max_size == numeric_limits::max()) + auto filter = this->make_filterer(data.size()); + if (!filter) { // Don't trim return; @@ -294,7 +299,6 @@ void ImportDataTrimmer::operator()(std::vector& data) std::vector result; result.reserve(std::min(options_.max_size + 1, data.size())); - auto filter = this->make_filterer(data.size()); for (auto i : range(data.size())) { if (filter(i)) @@ -312,8 +316,14 @@ void ImportDataTrimmer::operator()(std::vector& data) template void ImportDataTrimmer::operator()(std::map& data) { - std::map result; auto filter = this->make_filterer(data.size()); + if (!filter) + { + // Don't trim + return; + } + + std::map result; auto iter = data.begin(); for (auto i : range(filter.orig_size)) { @@ -369,6 +379,7 @@ auto ImportDataTrimmer::make_filterer(size_type vec_size) const -> GridFilterer */ bool ImportDataTrimmer::GridFilterer::operator()(size_type i) const { + CELER_EXPECT(stride > 0); return i % stride == 0 || i + 1 == orig_size; } diff --git a/src/celeritas/io/ImportMuPairProductionTable.hh b/src/celeritas/io/ImportMuPairProductionTable.hh index e04ebe17ac..3490edfff4 100644 --- a/src/celeritas/io/ImportMuPairProductionTable.hh +++ b/src/celeritas/io/ImportMuPairProductionTable.hh @@ -19,10 +19,10 @@ namespace celeritas * * This 3-dimensional table is used to sample the energy transfer to the * electron-positron pair, \f$ \epsilon_p \f$. The outer grid stores the atomic - * number using 5 equally spaced points in \f$ \log Z \f$; the x grid tabulates - * the ratio \f$ \log \epsilon_p / T \f$, where \f$ T \f$ is the incident muon - * energy; the y grid stores the incident muon energy using equal spacing in - * \f$ \log T \f$. + * number using 5 equally spaced points in \f$ \log Z \f$; the x grid stores + * the logarithm of the incident muon energy \f$ T \f$ using equal spacing + * in\f$ \log T \f$; the y grid stores the ratio \f$ \log \epsilon_p / T \f$. + * The values are the unnormalized CDF. */ struct ImportMuPairProductionTable { diff --git a/src/celeritas/phys/ProcessBuilder.cc b/src/celeritas/phys/ProcessBuilder.cc index ce7857a6bf..d8b1e05e77 100644 --- a/src/celeritas/phys/ProcessBuilder.cc +++ b/src/celeritas/phys/ProcessBuilder.cc @@ -21,6 +21,7 @@ #include "celeritas/em/process/GammaConversionProcess.hh" #include "celeritas/em/process/MuBremsstrahlungProcess.hh" #include "celeritas/em/process/MuIonizationProcess.hh" +#include "celeritas/em/process/MuPairProductionProcess.hh" #include "celeritas/em/process/PhotoelectricProcess.hh" #include "celeritas/em/process/RayleighProcess.hh" #include "celeritas/io/ImportData.hh" @@ -90,6 +91,8 @@ ProcessBuilder::ProcessBuilder(ImportData const& data, read_neutron_elastic_ = make_imported_element_loader(data.neutron_elastic_data); } + mu_pairprod_table_ = std::make_shared( + data.mu_pair_production_data); } //---------------------------------------------------------------------------// @@ -140,6 +143,7 @@ auto ProcessBuilder::operator()(IPC ipc) -> SPProcess {IPC::e_ioni, &ProcessBuilder::build_eioni}, {IPC::mu_brems, &ProcessBuilder::build_mubrems}, {IPC::mu_ioni, &ProcessBuilder::build_muioni}, + {IPC::mu_pair_prod, &ProcessBuilder::build_mupairprod}, {IPC::neutron_elastic, &ProcessBuilder::build_neutron_elastic}, {IPC::photoelectric, &ProcessBuilder::build_photoelectric}, {IPC::rayleigh, &ProcessBuilder::build_rayleigh}, @@ -274,6 +278,16 @@ auto ProcessBuilder::build_muioni() -> SPProcess this->particle(), this->imported(), options); } +//---------------------------------------------------------------------------// +auto ProcessBuilder::build_mupairprod() -> SPProcess +{ + MuPairProductionProcess::Options options; + options.use_integral_xs = use_integral_xs_; + + return std::make_shared( + this->particle(), this->imported(), options, mu_pairprod_table_); +} + //---------------------------------------------------------------------------// /*! * Warn and return a null process. diff --git a/src/celeritas/phys/ProcessBuilder.hh b/src/celeritas/phys/ProcessBuilder.hh index 10c187ef61..6bf2318e41 100644 --- a/src/celeritas/phys/ProcessBuilder.hh +++ b/src/celeritas/phys/ProcessBuilder.hh @@ -28,6 +28,7 @@ class MaterialParams; class ParticleParams; struct ImportData; struct ImportLivermorePE; +struct ImportMuPairProductionTable; //---------------------------------------------------------------------------// //! Options used for constructing built-in Celeritas processes @@ -115,6 +116,7 @@ class ProcessBuilder std::function read_sb_; std::function read_livermore_; std::function read_neutron_elastic_; + std::shared_ptr mu_pairprod_table_; BremsModelSelection selection_; bool brem_combined_; @@ -135,6 +137,7 @@ class ProcessBuilder auto build_eioni() -> SPProcess; auto build_mubrems() -> SPProcess; auto build_muioni() -> SPProcess; + auto build_mupairprod() -> SPProcess; auto build_msc() -> SPProcess; auto build_neutron_elastic() -> SPProcess; auto build_photoelectric() -> SPProcess; diff --git a/src/corecel/data/HyperslabIndexer.hh b/src/corecel/data/HyperslabIndexer.hh index 57bf7e52ed..b4f97d5579 100644 --- a/src/corecel/data/HyperslabIndexer.hh +++ b/src/corecel/data/HyperslabIndexer.hh @@ -22,6 +22,9 @@ namespace celeritas * Indexing is in standard C iteration order, such that final dimension * "changes fastest". For example, when indexing into a 3D grid (N=3) with * coords (i=0, j=0, k=1) the resulting index will be 1. + * + * \todo + * https://github.com/celeritas-project/celeritas/pull/1518#discussion_r1856614365 */ template class HyperslabIndexer @@ -42,6 +45,10 @@ class HyperslabIndexer //! Convert N-dimensional coordinates to an index inline CELER_FUNCTION size_type operator()(Coords const& coords) const; + //! Convert N-dimensional coordinates to an index + template + inline CELER_FUNCTION size_type operator()(Args... args) const; + private: //// DATA //// @@ -115,6 +122,20 @@ CELER_FUNCTION size_type HyperslabIndexer::operator()(Coords const& coords) c return result; } +//---------------------------------------------------------------------------// +/*! + * Convert N-dimensional coordinates to an index. + */ +template +template +CELER_FUNCTION size_type HyperslabIndexer::operator()(Args... args) const +{ + static_assert(sizeof...(Args) == N, + "Number of coordinates must match number of dimensions"); + + return this->operator()(Coords{static_cast(args)...}); +} + //---------------------------------------------------------------------------// /*! * Construct from array denoting the sizes of each dimension. diff --git a/test/celeritas/CMakeLists.txt b/test/celeritas/CMakeLists.txt index af009b4047..2d25a7e7e3 100644 --- a/test/celeritas/CMakeLists.txt +++ b/test/celeritas/CMakeLists.txt @@ -130,7 +130,7 @@ celeritas_add_test(decay/MuDecay.test.cc ${_needs_double}) # Distributions celeritas_add_test(em/distribution/EnergyLossHelper.test.cc ${_fixme_single}) -celeritas_add_test(em/distribution/MuBremsPPAngularDistribution.test.cc ${_needs_double}) +celeritas_add_test(em/distribution/MuAngularDistribution.test.cc ${_needs_double}) celeritas_add_test(em/distribution/TsaiUrbanDistribution.test.cc ${_needs_double}) # Cross section calculation @@ -147,6 +147,7 @@ celeritas_add_test(em/LivermorePE.test.cc ${_needs_double}) celeritas_add_test(em/MollerBhabha.test.cc ${_needs_double}) celeritas_add_test(em/MuBetheBloch.test.cc ${_needs_double}) celeritas_add_test(em/MuBremsstrahlung.test.cc ${_needs_double}) +celeritas_add_test(em/MuPairProduction.test.cc ${_needs_double} ${_needs_root}) celeritas_add_test(em/Rayleigh.test.cc ${_needs_double}) celeritas_add_test(em/RelativisticBrem.test.cc ${_needs_double}) celeritas_add_test(em/SeltzerBerger.test.cc ${_needs_double}) diff --git a/test/celeritas/data/four-steel-slabs.geant.json b/test/celeritas/data/four-steel-slabs.geant.json index 697e6e8746..c6c5813fa1 100644 --- a/test/celeritas/data/four-steel-slabs.geant.json +++ b/test/celeritas/data/four-steel-slabs.geant.json @@ -3,7 +3,8 @@ "coulomb_scattering": true, "muon": { "ionization": true, - "bremsstrahlung": true + "bremsstrahlung": true, + "pair_production": true }, "eloss_fluctuation": true, "em_bins_per_decade": 7, diff --git a/test/celeritas/data/four-steel-slabs.root b/test/celeritas/data/four-steel-slabs.root index d2913600cc..327ed8ba57 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/lar-sphere.root b/test/celeritas/data/lar-sphere.root index 5e1c200c6d..180236d03b 100644 Binary files a/test/celeritas/data/lar-sphere.root and b/test/celeritas/data/lar-sphere.root differ diff --git a/test/celeritas/data/simple-cms.root b/test/celeritas/data/simple-cms.root index aa078f7959..5bfa4d4f99 100644 Binary files a/test/celeritas/data/simple-cms.root and b/test/celeritas/data/simple-cms.root differ diff --git a/test/celeritas/em/MuPairProduction.test.cc b/test/celeritas/em/MuPairProduction.test.cc new file mode 100644 index 0000000000..104bf41311 --- /dev/null +++ b/test/celeritas/em/MuPairProduction.test.cc @@ -0,0 +1,383 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/em/MuPairProduction.test.cc +//---------------------------------------------------------------------------// +#include "corecel/cont/Range.hh" +#include "corecel/math/ArrayUtils.hh" +#include "celeritas/RootTestBase.hh" +#include "celeritas/Quantities.hh" +#include "celeritas/em/distribution/MuPPEnergyDistribution.hh" +#include "celeritas/em/interactor/MuPairProductionInteractor.hh" +#include "celeritas/em/model/MuPairProductionModel.hh" +#include "celeritas/io/ImportData.hh" +#include "celeritas/mat/MaterialTrackView.hh" +#include "celeritas/phys/CutoffView.hh" +#include "celeritas/phys/InteractionIO.hh" +#include "celeritas/phys/InteractorHostTestBase.hh" + +#include "celeritas_test.hh" + +namespace celeritas +{ +namespace test +{ +//---------------------------------------------------------------------------// +// TEST HARNESS +//---------------------------------------------------------------------------// + +class MuPairProductionTest : public InteractorHostBase, public RootTestBase +{ + protected: + void SetUp() override + { + using namespace units; + + // Set up shared material data + MaterialParams::Input mat_inp; + mat_inp.elements = {{AtomicNumber{29}, AmuMass{63.546}, {}, "Cu"}}; + mat_inp.materials = { + {native_value_from(MolCcDensity{0.141}), + 293.0, + MatterState::solid, + {{ElementId{0}, 1.0}}, + "Cu"}, + }; + this->set_material_params(mat_inp); + + // Set 1 keV cutoff + CutoffParams::Input cut_inp; + cut_inp.materials = this->material_params(); + cut_inp.particles = this->particle_params(); + cut_inp.cutoffs = {{pdg::positron(), {{MevEnergy{0.001}, 0.1234}}}}; + this->set_cutoff_params(cut_inp); + + // Construct model + auto imported = std::make_shared( + this->imported_data().processes); + model_ = std::make_shared( + ActionId{0}, + *this->particle_params(), + imported, + this->imported_data().mu_pair_production_data); + + // Set default particle to 10 GeV muon + this->set_inc_particle(pdg::mu_minus(), MevEnergy{1e4}); + this->set_inc_direction({0, 0, 1}); + this->set_material("Cu"); + } + + void sanity_check(Interaction const& interaction) const + { + // Check change to parent track + EXPECT_GT(this->particle_track().energy().value(), + interaction.energy.value()); + EXPECT_LT(0, interaction.energy.value()); + EXPECT_SOFT_EQ(1.0, norm(interaction.direction)); + EXPECT_EQ(Action::scattered, interaction.action); + + // Check secondaries + ASSERT_EQ(2, interaction.secondaries.size()); + auto const& electron = interaction.secondaries[0]; + EXPECT_TRUE(electron); + EXPECT_GT(this->particle_track().energy(), electron.energy); + EXPECT_LT(zero_quantity(), electron.energy); + EXPECT_SOFT_EQ(1.0, norm(electron.direction)); + EXPECT_EQ(model_->host_ref().ids.electron, electron.particle_id); + + auto const& positron = interaction.secondaries[1]; + EXPECT_TRUE(positron); + EXPECT_GT(this->particle_track().energy(), positron.energy); + EXPECT_LT(zero_quantity(), positron.energy); + EXPECT_SOFT_EQ(1.0, norm(positron.direction)); + EXPECT_EQ(model_->host_ref().ids.positron, positron.particle_id); + + // Check conservation between primary and secondaries + // this->check_conservation(interaction); + this->check_energy_conservation(interaction); + } + + //! \note These tests use a trimmed element table + std::string_view geometry_basename() const final + { + return "four-steel-slabs"; + } + + SPConstTrackInit build_init() override { CELER_ASSERT_UNREACHABLE(); } + SPConstAction build_along_step() override { CELER_ASSERT_UNREACHABLE(); } + + protected: + std::shared_ptr model_; +}; + +//---------------------------------------------------------------------------// +// TESTS +//---------------------------------------------------------------------------// + +TEST_F(MuPairProductionTest, distribution) +{ + int num_samples = 10000; + int num_bins = 12; + + real_type two_me + = 2 * value_as(model_->host_ref().electron_mass); + + // Get view to the current element + auto element + = this->material_track().make_material_view().make_element_view( + ElementComponentId{0}); + + // Get the production cuts + auto cutoff = this->cutoff_params()->get(MaterialId{0}); + + RandomEngine& rng = InteractorHostBase::rng(); + + std::vector counters; + std::vector min_energy; + std::vector max_energy; + std::vector avg_energy; + std::vector avg_energy_fraction; + for (real_type energy : {1e3, 1e4, 1e5, 1e6, 1e7}) + { + this->set_inc_particle(pdg::mu_minus(), MevEnergy(energy)); + + MuPPEnergyDistribution sample( + model_->host_ref(), this->particle_track(), cutoff, element); + real_type min = value_as(sample.min_pair_energy()) - two_me; + real_type max = value_as(sample.max_pair_energy()) - two_me; + + real_type sum_energy = 0; + real_type energy_fraction = 0; + std::vector count(num_bins); + for ([[maybe_unused]] int i : range(num_samples)) + { + // TODO: test energy partition + auto energy = sample(rng); + auto r = value_as(energy.electron + energy.positron); + ASSERT_GE(r, min); + ASSERT_LE(r, max); + int bin = int(std::log(r / min) / std::log(max / min) * num_bins); + CELER_ASSERT(bin >= 0 && bin < num_bins); + ++count[bin]; + sum_energy += r; + energy_fraction += value_as(energy.electron) / r; + } + counters.insert(counters.end(), count.begin(), count.end()); + min_energy.push_back(min); + max_energy.push_back(max); + avg_energy.push_back(sum_energy / num_samples); + avg_energy_fraction.push_back(energy_fraction / num_samples); + } + + static int const expected_counters[] = { + 162, 1066, 2154, 2691, 2020, 1115, 518, 196, 60, 13, 5, 0, + 270, 931, 1869, 2496, 2103, 1514, 534, 212, 51, 15, 5, 0, + 208, 811, 1608, 2146, 2099, 1668, 947, 387, 101, 20, 3, 2, + 203, 782, 1564, 1987, 2015, 1717, 1058, 489, 161, 16, 8, 0, + 197, 767, 1548, 1948, 1965, 1607, 1116, 605, 200, 43, 4, 0, + }; + static double const expected_min_energy[] = { + 1.0219978922, + 1.0219978922, + 1.0219978922, + 1.0219978922, + 1.0219978922, + }; + static double const expected_max_energy[] = { + 703.23539643546, + 9703.2353964355, + 99703.235396435, + 999703.23539644, + 9999703.2353964, + }; + static double const expected_avg_energy[] = { + 11.551429194638, + 42.5605253686, + 216.32003549047, + 1093.2367878798, + 6041.1513546243, + }; + static double const expected_avg_energy_fraction[] = { + 0.50427974126835, + 0.50112476955009, + 0.49759901188728, + 0.50543113944062, + 0.50102592483879, + }; + EXPECT_VEC_EQ(expected_counters, counters); + EXPECT_VEC_SOFT_EQ(expected_min_energy, min_energy); + EXPECT_VEC_SOFT_EQ(expected_max_energy, max_energy); + EXPECT_VEC_SOFT_EQ(expected_avg_energy, avg_energy); + EXPECT_VEC_SOFT_EQ(expected_avg_energy_fraction, avg_energy_fraction); +} + +TEST_F(MuPairProductionTest, basic) +{ + // Reserve 8 secondaries, two for each sample + int const num_samples = 4; + this->resize_secondaries(2 * num_samples); + + // Get view to the current element + auto element + = this->material_track().make_material_view().make_element_view( + ElementComponentId{0}); + + // Get the production cuts + auto cutoff = this->cutoff_params()->get(MaterialId{0}); + + // Create the interactor + MuPairProductionInteractor interact(model_->host_ref(), + this->particle_track(), + cutoff, + element, + this->direction(), + this->secondary_allocator()); + RandomEngine& rng = InteractorHostBase::rng(); + + std::vector pair_energy; + std::vector costheta; + + // Produce four samples from the original incident energy + for (int i : range(num_samples)) + { + Interaction result = interact(rng); + SCOPED_TRACE(result); + this->sanity_check(result); + + EXPECT_EQ(result.secondaries.data(), + this->secondary_allocator().get().data() + + result.secondaries.size() * i); + + pair_energy.push_back(value_as( + result.secondaries[0].energy + result.secondaries[1].energy)); + costheta.push_back(dot_product(result.secondaries[0].direction, + result.secondaries[1].direction)); + } + + EXPECT_EQ(2 * num_samples, this->secondary_allocator().get().size()); + + // Note: these are "gold" values based on the host RNG. + static double const expected_pair_energy[] = { + 5.1981351222035, + 21.411122079708, + 39.340205211007, + 1.2067098240449, + }; + static double const expected_costheta[] = { + 0.99992128683238, + 0.97331314773255, + 0.9996196536095, + 0.99925389709579, + }; + EXPECT_VEC_SOFT_EQ(expected_pair_energy, pair_energy); + EXPECT_VEC_SOFT_EQ(expected_costheta, costheta); + + // Next sample should fail because we're out of secondary buffer space + { + Interaction result = interact(rng); + EXPECT_EQ(0, result.secondaries.size()); + EXPECT_EQ(Action::failed, result.action); + } +} + +TEST_F(MuPairProductionTest, stress_test) +{ + unsigned int const num_samples = 10000; + std::vector avg_engine_samples; + std::vector avg_electron_energy; + std::vector avg_positron_energy; + std::vector avg_costheta; + + // Get view to the current element + auto element + = this->material_track().make_material_view().make_element_view( + ElementComponentId{0}); + + // Get the production cuts + auto cutoff = this->cutoff_params()->get(MaterialId{0}); + + for (real_type inc_e : {1e3, 1e4, 1e5, 1e6, 1e7}) + { + SCOPED_TRACE("Incident energy: " + std::to_string(inc_e)); + this->set_inc_particle(pdg::mu_minus(), MevEnergy{inc_e}); + + RandomEngine& rng = InteractorHostBase::rng(); + RandomEngine::size_type num_particles_sampled = 0; + double electron_energy = 0; + double positron_energy = 0; + double costheta = 0; + + // Loop over several incident directions + for (Real3 const& inc_dir : + {Real3{0, 0, 1}, Real3{1, 0, 0}, Real3{1e-9, 0, 1}, Real3{1, 1, 1}}) + { + SCOPED_TRACE("Incident direction: " + to_string(inc_dir)); + this->set_inc_direction(inc_dir); + this->resize_secondaries(2 * num_samples); + + // Create the interactor + MuPairProductionInteractor interact(model_->host_ref(), + this->particle_track(), + cutoff, + element, + this->direction(), + this->secondary_allocator()); + + // Loop over many particles + for (unsigned int i = 0; i < num_samples; ++i) + { + Interaction result = interact(rng); + this->sanity_check(result); + + electron_energy + += value_as(result.secondaries[0].energy); + positron_energy + += value_as(result.secondaries[1].energy); + costheta += dot_product(result.secondaries[0].direction, + result.secondaries[1].direction); + } + EXPECT_EQ(2 * num_samples, + this->secondary_allocator().get().size()); + num_particles_sampled += num_samples; + } + avg_engine_samples.push_back(real_type(rng.count()) + / num_particles_sampled); + avg_electron_energy.push_back(electron_energy / num_particles_sampled); + avg_positron_energy.push_back(positron_energy / num_particles_sampled); + avg_costheta.push_back(costheta / num_particles_sampled); + } + + // Gold values for average number of calls to RNG + static double const expected_avg_engine_samples[] = {10, 10, 10, 10, 10}; + static double const expected_avg_electron_energy[] = { + 5.9452433822303, + 20.776282536509, + 98.201429477115, + 555.92710681765, + 2855.9810205079, + }; + static double const expected_avg_positron_energy[] = { + 5.8651456808897, + 21.483133310816, + 100.92564951414, + 546.95048450797, + 2824.7431627774, + }; + static double const expected_avg_costheta[] = { + 0.94178280008365, + 0.99880165151034, + 0.99998776687485, + 0.99999983141391, + 0.99999999832285, + }; + EXPECT_VEC_SOFT_EQ(expected_avg_engine_samples, avg_engine_samples); + EXPECT_VEC_SOFT_EQ(expected_avg_electron_energy, avg_electron_energy); + EXPECT_VEC_SOFT_EQ(expected_avg_positron_energy, avg_positron_energy); + EXPECT_VEC_SOFT_EQ(expected_avg_costheta, avg_costheta); +} + +//---------------------------------------------------------------------------// +} // namespace test +} // namespace celeritas diff --git a/test/celeritas/em/distribution/MuBremsPPAngularDistribution.test.cc b/test/celeritas/em/distribution/MuAngularDistribution.test.cc similarity index 88% rename from test/celeritas/em/distribution/MuBremsPPAngularDistribution.test.cc rename to test/celeritas/em/distribution/MuAngularDistribution.test.cc index 28b4894f4b..fdab5e59ef 100644 --- a/test/celeritas/em/distribution/MuBremsPPAngularDistribution.test.cc +++ b/test/celeritas/em/distribution/MuAngularDistribution.test.cc @@ -3,9 +3,9 @@ // See the top-level COPYRIGHT file for details. // SPDX-License-Identifier: (Apache-2.0 OR MIT) //---------------------------------------------------------------------------// -//! \file celeritas/em/distribution/MuBremsPPAngularDistribution.test.cc +//! \file celeritas/em/distribution/MuAngularDistribution.test.cc //---------------------------------------------------------------------------// -#include "celeritas/em/distribution/MuBremsPPAngularDistribution.hh" +#include "celeritas/em/distribution/MuAngularDistribution.hh" #include @@ -19,7 +19,7 @@ namespace celeritas namespace test { //---------------------------------------------------------------------------// -TEST(MuBremsPPAngularDistributionTest, costheta_dist) +TEST(MuAngularDistributionTest, costheta_dist) { using Energy = units::MevEnergy; using Mass = units::MevMass; @@ -33,7 +33,7 @@ TEST(MuBremsPPAngularDistributionTest, costheta_dist) { for (real_type eps : {0.001, 0.01, 0.1}) { - MuBremsPPAngularDistribution sample_costheta( + MuAngularDistribution sample_costheta( Energy{inc_e}, muon_mass, Energy{eps * inc_e}); real_type costheta_sum = 0; diff --git a/test/celeritas/ext/GeantImporter.test.cc b/test/celeritas/ext/GeantImporter.test.cc index 48bd047f6c..94474ad0fa 100644 --- a/test/celeritas/ext/GeantImporter.test.cc +++ b/test/celeritas/ext/GeantImporter.test.cc @@ -1123,18 +1123,6 @@ TEST_F(FourSteelSlabsEmStandard, mu_pair_production_data) real_type const tol = geant4_version < Version(11, 1, 0) ? 1e-12 : 0.03; static double const expected_table_x[] = { - -6.1928487397154, - 0, - -6.1928487397154, - 0, - -6.1928487397154, - 0, - -6.1928487397154, - 0, - -6.1928487397154, - 0, - }; - static double const expected_table_y[] = { 6.9077552789821, 18.420680743952, 6.9077552789821, @@ -1146,6 +1134,18 @@ TEST_F(FourSteelSlabsEmStandard, mu_pair_production_data) 6.9077552789821, 18.420680743952, }; + static double const expected_table_y[] = { + -6.1928487397154, + 0, + -6.1928487397154, + 0, + -6.1928487397154, + 0, + -6.1928487397154, + 0, + -6.1928487397154, + 0, + }; static double const expected_table_value[] = { 0, 0.24363843626056, diff --git a/test/celeritas/ext/RootImporter.test.cc b/test/celeritas/ext/RootImporter.test.cc index 7b2e890086..20ad50f58f 100644 --- a/test/celeritas/ext/RootImporter.test.cc +++ b/test/celeritas/ext/RootImporter.test.cc @@ -151,7 +151,7 @@ TEST_F(RootImporterTest, phys_materials) TEST_F(RootImporterTest, processes) { auto const& processes = this->imported_data().processes; - EXPECT_EQ(15, processes.size()); + EXPECT_EQ(17, processes.size()); auto find_process = [&processes](PDGNumber pdg, ImportProcessClass ipc) { return std::find_if(processes.begin(), diff --git a/test/celeritas/ext/RootJsonDumper.test.cc b/test/celeritas/ext/RootJsonDumper.test.cc index 3f0ad37c4a..15fc3cfcfe 100644 --- a/test/celeritas/ext/RootJsonDumper.test.cc +++ b/test/celeritas/ext/RootJsonDumper.test.cc @@ -42,6 +42,7 @@ TEST_F(RootJsonDumperTest, all) ImportDataTrimmer::Input inp; inp.materials = true; inp.physics = true; + inp.mupp = true; inp.max_size = 2; ImportDataTrimmer trim{inp}; trim(imported); @@ -170,8 +171,6 @@ TEST_F(RootJsonDumperTest, all) "range" : 0.1 }}] }], -"optical_models" : [], -"optical_materials" : [], "regions" : [{ "_typename" : "celeritas::ImportRegion", "name" : "DefaultRegionForTheWorld", @@ -325,8 +324,18 @@ TEST_F(RootJsonDumperTest, all) "atomic_relaxation_data" : [], "mu_pair_production_data" : { "_typename" : "celeritas::ImportMuPairProductionTable", - "atomic_number" : [], - "physics_vectors" : [] + "atomic_number" : [1, 92], + "physics_vectors" : [{ + "_typename" : "celeritas::ImportPhysics2DVector", + "x" : [6.90775527898214, 18.4206807439524], + "y" : [-6.19284873971536, 0], + "value" : [0, 4.0853712905423e-28, 0, 2.43638436260562e-25] + }, { + "_typename" : "celeritas::ImportPhysics2DVector", + "x" : [6.90775527898214, 18.4206807439524], + "y" : [-6.19284873971536, 0], + "value" : [0, 3.30246663127583e-24, 0, 7.93413967608228e-22] + }] }, "em_params" : { "_typename" : "celeritas::ImportEmParameters", @@ -375,6 +384,8 @@ TEST_F(RootJsonDumperTest, all) "_typename" : "celeritas::ImportOpticalParameters", "scintillation_by_particle" : false }, +"optical_models" : [], +"optical_materials" : [], "units" : "cgs" })json", str); diff --git a/test/celeritas/phys/InteractorHostTestBase.cc b/test/celeritas/phys/InteractorHostTestBase.cc index 263dc158df..e0e7d70991 100644 --- a/test/celeritas/phys/InteractorHostTestBase.cc +++ b/test/celeritas/phys/InteractorHostTestBase.cc @@ -24,7 +24,7 @@ namespace test /*! * Initialize secondary allocation on construction. */ -InteractorHostTestBase::InteractorHostTestBase() +InteractorHostBase::InteractorHostBase() { using namespace constants; using namespace units; @@ -117,13 +117,13 @@ InteractorHostTestBase::InteractorHostTestBase() /*! * Default destructor. */ -InteractorHostTestBase::~InteractorHostTestBase() = default; +InteractorHostBase::~InteractorHostBase() = default; //---------------------------------------------------------------------------// /*! * Helper to make dummy ImportProcess data . */ -ImportProcess InteractorHostTestBase::make_import_process( +ImportProcess InteractorHostBase::make_import_process( PDGNumber particle, PDGNumber secondary, ImportProcessClass ipc, @@ -159,7 +159,7 @@ ImportProcess InteractorHostTestBase::make_import_process( /*! * Set particle parameters. */ -void InteractorHostTestBase::set_material_params(MaterialParams::Input inp) +void InteractorHostBase::set_material_params(MaterialParams::Input inp) { CELER_EXPECT(!inp.materials.empty()); @@ -172,7 +172,7 @@ void InteractorHostTestBase::set_material_params(MaterialParams::Input inp) /*! * Initialize the incident track's material */ -void InteractorHostTestBase::set_material(std::string const& name) +void InteractorHostBase::set_material(std::string const& name) { CELER_EXPECT(material_params_); @@ -190,7 +190,7 @@ void InteractorHostTestBase::set_material(std::string const& name) /*! * Set particle parameters. */ -void InteractorHostTestBase::set_particle_params(ParticleParams::Input inp) +void InteractorHostBase::set_particle_params(ParticleParams::Input inp) { CELER_EXPECT(!inp.empty()); particle_params_ = std::make_shared(std::move(inp)); @@ -202,7 +202,7 @@ void InteractorHostTestBase::set_particle_params(ParticleParams::Input inp) /*! * Set cutoff parameters. */ -void InteractorHostTestBase::set_cutoff_params(CutoffParams::Input inp) +void InteractorHostBase::set_cutoff_params(CutoffParams::Input inp) { CELER_EXPECT(inp.materials && inp.particles); cutoff_params_ = std::make_shared(std::move(inp)); @@ -212,7 +212,7 @@ void InteractorHostTestBase::set_cutoff_params(CutoffParams::Input inp) /*! * Set imported processes. */ -void InteractorHostTestBase::set_imported_processes( +void InteractorHostBase::set_imported_processes( std::vector inp) { CELER_EXPECT(!inp.empty()); @@ -223,7 +223,7 @@ void InteractorHostTestBase::set_imported_processes( /*! * Initialize the incident particle data */ -void InteractorHostTestBase::set_inc_particle(PDGNumber pdg, MevEnergy energy) +void InteractorHostBase::set_inc_particle(PDGNumber pdg, MevEnergy energy) { CELER_EXPECT(particle_params_); CELER_EXPECT(pdg); @@ -244,7 +244,7 @@ void InteractorHostTestBase::set_inc_particle(PDGNumber pdg, MevEnergy energy) /*! * Set an incident direction (and normalize it). */ -void InteractorHostTestBase::set_inc_direction(Real3 const& dir) +void InteractorHostBase::set_inc_direction(Real3 const& dir) { CELER_EXPECT(norm(dir) > 0); @@ -255,7 +255,7 @@ void InteractorHostTestBase::set_inc_direction(Real3 const& dir) /*! * Resize secondaries. */ -void InteractorHostTestBase::resize_secondaries(int count) +void InteractorHostBase::resize_secondaries(int count) { CELER_EXPECT(count > 0); secondaries_ = StateStore(count); @@ -266,7 +266,7 @@ void InteractorHostTestBase::resize_secondaries(int count) /*! * Check for energy and momentum conservation in the interaction. */ -void InteractorHostTestBase::check_conservation(Interaction const& interaction) const +void InteractorHostBase::check_conservation(Interaction const& interaction) const { ASSERT_NE(interaction.action, Action::failed); @@ -278,7 +278,7 @@ void InteractorHostTestBase::check_conservation(Interaction const& interaction) /*! * Check for energy conservation in the interaction. */ -void InteractorHostTestBase::check_energy_conservation( +void InteractorHostBase::check_energy_conservation( Interaction const& interaction) const { // Sum of exiting kinetic energy @@ -296,7 +296,8 @@ void InteractorHostTestBase::check_energy_conservation( exit_energy += s.energy.value(); // Account for positron production - if (s && s.particle_id == particle_params_->find(pdg::positron())) + if (s && s.particle_id == particle_params_->find(pdg::positron()) + && interaction.action == Action::absorbed) { exit_energy += 2 * particle_params_->get(s.particle_id).mass().value(); @@ -311,7 +312,7 @@ void InteractorHostTestBase::check_energy_conservation( /*! * Check for momentum conservation in the interaction. */ -void InteractorHostTestBase::check_momentum_conservation( +void InteractorHostBase::check_momentum_conservation( Interaction const& interaction) const { CollectionStateStore temp_store( diff --git a/test/celeritas/phys/InteractorHostTestBase.hh b/test/celeritas/phys/InteractorHostTestBase.hh index 591837d414..76ea1e8228 100644 --- a/test/celeritas/phys/InteractorHostTestBase.hh +++ b/test/celeritas/phys/InteractorHostTestBase.hh @@ -49,7 +49,7 @@ namespace test * \todo Since this now uses Collection objects it's generally safe to use this * to test Models as well as device code -- think about renaming it. */ -class InteractorHostTestBase : public Test +class InteractorHostBase { public: //!@{ @@ -64,8 +64,8 @@ class InteractorHostTestBase : public Test public: //!@{ //! Initialize and destroy - InteractorHostTestBase(); - ~InteractorHostTestBase(); + InteractorHostBase(); + ~InteractorHostBase(); //!@} // Helper to make dummy ImportProcess @@ -189,6 +189,10 @@ class InteractorHostTestBase : public Test std::shared_ptr sa_view_; }; +class InteractorHostTestBase : public InteractorHostBase, public Test +{ +}; + //---------------------------------------------------------------------------// } // namespace test } // namespace celeritas diff --git a/test/corecel/data/HyperslabIndexer.test.cc b/test/corecel/data/HyperslabIndexer.test.cc index 6a2589f819..b607c58b0e 100644 --- a/test/corecel/data/HyperslabIndexer.test.cc +++ b/test/corecel/data/HyperslabIndexer.test.cc @@ -30,6 +30,7 @@ TEST(HyperslabIndexerTest, 2D) { Array coords{a, b}; EXPECT_EQ(index, to_index(coords)); + EXPECT_EQ(index, to_index(a, b)); EXPECT_EQ(coords, to_coords(index)); index++; } @@ -51,6 +52,7 @@ TEST(HyperslabIndexerTest, 3D) { Array coords{a, b, c}; EXPECT_EQ(index, to_index(coords)); + EXPECT_EQ(index, to_index(a, b, c)); EXPECT_EQ(coords, to_coords(index)); index++; } @@ -75,6 +77,7 @@ TEST(HyperslabIndexerTest, 4D) { Array coords{a, b, c, d}; EXPECT_EQ(index, to_index(coords)); + EXPECT_EQ(index, to_index(a, b, c, d)); EXPECT_EQ(coords, to_coords(index)); index++; } @@ -98,6 +101,7 @@ TEST(HyperslabIndexerTest, 5D_with_ones) { Array coords{a, 0, b, 0, c}; EXPECT_EQ(index, to_index(coords)); + EXPECT_EQ(index, to_index(a, 0, b, 0, c)); EXPECT_EQ(coords, to_coords(index)); index++; }