Add muon pair production

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);
     }

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
 ------------------
 

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)

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<Ownership W, MemSpace M>
+struct MuPairProductionTableData
+{
+    //// TYPES ////
+
+    template<class T>
+    using Items = Collection<T, W, M>;
+
+    //// MEMBER DATA ////
+
+    ItemRange<real_type> logz_grid;
+    Items<TwodGridData> grids;
+
+    // Backend data
+    Items<real_type> 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<Ownership W2, MemSpace M2>
+    MuPairProductionTableData&
+    operator=(MuPairProductionTableData<W2, M2> 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<Ownership W, MemSpace M>
+struct MuPairProductionData
+{
+    //// MEMBER DATA ////
+
+    //! Particle IDs
+    MuPairProductionIds ids;
+
+    //! Electron mass [MeV / c^2]
+    units::MevMass electron_mass;
+
+    // Sampling table storage
+    MuPairProductionTableData<W, M> 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<Ownership W2, MemSpace M2>
+    MuPairProductionData& operator=(MuPairProductionData<W2, M2> const& other)
+    {
+        CELER_EXPECT(other);
+        ids = other.ids;
+        electron_mass = other.electron_mass;
+        table = other.table;
+        return *this;
+    }
+};
+
+//---------------------------------------------------------------------------//
+}  // namespace celeritas

src/celeritas/em/distribution/MuBremsPPAngularDistribution.hh → 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<class Engine>
@@ -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<Energy>(inc_energy) / value_as<Mass>(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<class Engine>
-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_);