Add Seltzer-Berger interactor and kernel

src/CMakeLists.txt

@@ -93,6 +93,7 @@ if(CELERITAS_USE_CUDA)
     physics/em/detail/MollerBhabha.cu
     physics/em/detail/LivermorePE.cu
     physics/em/detail/Rayleigh.cu
+    physics/em/detail/SeltzerBerger.cu
     random/detail/RngStateInit.cu
     sim/detail/SimStateInit.cu
   )

src/physics/em/SeltzerBergerModel.cc

@@ -44,8 +44,7 @@ SeltzerBergerModel::SeltzerBergerModel(ModelId               id,
                    << this->label() << ")");
 
     // Save particle properties
-    host_data.electron_mass
-        = particles.get(host_data.ids.electron).mass().value();
+    host_data.electron_mass = particles.get(host_data.ids.electron).mass();
 
     // Load differential cross sections
     make_builder(&host_data.differential_xs.elements)

src/physics/em/detail/BetheHeitlerInteractor.hh

@@ -36,6 +36,12 @@ namespace detail
  */
 class BetheHeitlerInteractor
 {
+  public:
+    //!@{
+    //! Type aliases
+    using MevMass = units::MevMass;
+    //!@}
+
   public:
     //! Construct sampler from shared and state data
     inline CELER_FUNCTION
@@ -76,13 +82,6 @@ class BetheHeitlerInteractor
     // Gamma energy divided by electron mass * csquared
     const BetheHeitlerPointers& shared_;
 
-    // Sample outgoing particles directions.
-    // Based on the G4ModifiedTsai sampler, a simplified sampler that does not
-    // require exact momentum conservation (due to neglecting nucleus recoil).
-    template<class Engine>
-    inline CELER_FUNCTION real_type sample_cos_theta(real_type kinetic_energy,
-                                                     Engine&   rng);
-
     // Incident gamma energy
     const units::MevEnergy inc_energy_;
 

src/physics/em/detail/BetheHeitlerInteractor.i.hh

@@ -8,6 +8,7 @@
 
 #include "base/ArrayUtils.hh"
 #include "base/Constants.hh"
+#include "TsaiUrbanDistribution.hh"
 #include "random/distributions/BernoulliDistribution.hh"
 #include "random/distributions/GenerateCanonical.hh"
 #include "random/distributions/UniformRealDistribution.hh"
@@ -149,34 +150,20 @@ CELER_FUNCTION Interaction BetheHeitlerInteractor::operator()(Engine& rng)
     real_type                          phi = sample_phi(rng);
 
     // Electron
-    real_type cost = this->sample_cos_theta(secondaries[0].energy.value(), rng);
+    TsaiUrbanDistribution sample_electron_angle(
+        secondaries[0].energy, MevMass{1 / shared_.inv_electron_mass});
+    real_type cost = sample_electron_angle(rng);
     secondaries[0].direction
         = rotate(from_spherical(cost, phi), inc_direction_);
     // Positron
-    cost = this->sample_cos_theta(secondaries[1].energy.value(), rng);
+    TsaiUrbanDistribution sample_positron_angle(
+        secondaries[1].energy, MevMass{1 / shared_.inv_electron_mass});
+    cost = sample_positron_angle(rng);
     secondaries[1].direction
         = rotate(from_spherical(cost, phi + constants::pi), inc_direction_);
     return result;
 }
 
-template<class Engine>
-CELER_FUNCTION real_type
-BetheHeitlerInteractor::sample_cos_theta(real_type kinetic_energy, Engine& rng)
-{
-    real_type umax = 2 * (1 + kinetic_energy * shared_.inv_electron_mass);
-    real_type u;
-    do
-    {
-        real_type uu
-            = -std::log(generate_canonical(rng) * generate_canonical(rng));
-        u = uu
-            * (BernoulliDistribution(0.25)(rng) ? real_type(1.6)
-                                                : real_type(1.6 / 3));
-    } while (u > umax);
-
-    return 1 - 2 * ipow<2>(u / umax);
-}
-
 CELER_FUNCTION real_type
 BetheHeitlerInteractor::impact_parameter(real_type eps) const
 {

src/physics/em/detail/SeltzerBerger.cu

@@ -0,0 +1,106 @@
+//---------------------------------*-CUDA-*----------------------------------//
+// Copyright 2020 UT-Battelle, LLC, and other Celeritas developers.
+// See the top-level COPYRIGHT file for details.
+// SPDX-License-Identifier: (Apache-2.0 OR MIT)
+//---------------------------------------------------------------------------//
+//! \file SeltzerBerger.cu
+//---------------------------------------------------------------------------//
+#include "SeltzerBerger.hh"
+
+#include "base/Assert.hh"
+#include "base/KernelParamCalculator.cuda.hh"
+#include "random/RngEngine.hh"
+#include "physics/base/ModelInterface.hh"
+#include "physics/base/ParticleTrackView.hh"
+#include "physics/base/PhysicsTrackView.hh"
+#include "base/StackAllocator.hh"
+#include "physics/material/MaterialTrackView.hh"
+#include "SeltzerBergerInteractor.hh"
+
+namespace celeritas
+{
+namespace detail
+{
+namespace
+{
+//---------------------------------------------------------------------------//
+// KERNELS
+//---------------------------------------------------------------------------//
+/*!
+ * Interact using the Seltzer-Berger model on applicable tracks.
+ */
+__global__ void seltzer_berger_interact_kernel(
+    const SeltzerBergerDeviceRef&              device_pointers,
+    const ModelInteractRefs<MemSpace::device>& interaction)
+{
+    auto tid = celeritas::KernelParamCalculator::thread_id();
+    if (!(tid < interaction.states.size()))
+        return;
+
+    StackAllocator<Secondary> allocate_secondaries(
+        interaction.states.secondaries);
+    ParticleTrackView particle(
+        interaction.params.particle, interaction.states.particle, tid);
+
+    // Setup for ElementView access
+    MaterialTrackView material(
+        interaction.params.material, interaction.states.material, tid);
+
+    PhysicsTrackView physics(interaction.params.physics,
+                             interaction.states.physics,
+                             particle.particle_id(),
+                             material.material_id(),
+                             tid);
+
+    // This interaction only applies if the Seltzer-Berger model was selected
+    if (physics.model_id() != device_pointers.ids.model)
+        return;
+
+    CutoffView cutoffs(interaction.params.cutoffs, material.material_id());
+
+    // Cache the associated MaterialView as function calls to MaterialTrackView
+    // are expensive
+    MaterialView    material_view = material.material_view();
+    const ElementId element_id{0};
+
+    // Assume only a single element in the material, for now
+    CELER_ASSERT(material_view.num_elements() == 1);
+    SeltzerBergerInteractor interact(device_pointers,
+                                     particle,
+                                     interaction.states.direction[tid],
+                                     cutoffs,
+                                     allocate_secondaries,
+                                     material_view,
+                                     element_id);
+
+    RngEngine rng(interaction.states.rng, tid);
+    interaction.states.interactions[tid] = interact(rng);
+    CELER_ENSURE(interaction.states.interactions[tid]);
+}
+
+} // namespace
+
+//---------------------------------------------------------------------------//
+// LAUNCHERS
+//---------------------------------------------------------------------------//
+/*!
+ * Launch the Seltzer-Berger interaction.
+ */
+void seltzer_berger_interact(
+    const SeltzerBergerDeviceRef&              device_pointers,
+    const ModelInteractRefs<MemSpace::device>& interaction)
+{
+    CELER_EXPECT(device_pointers);
+    CELER_EXPECT(interaction);
+
+    static const KernelParamCalculator calc_kernel_params(
+        seltzer_berger_interact_kernel, "seltzer_berger_interact");
+    auto params = calc_kernel_params(interaction.states.size());
+    seltzer_berger_interact_kernel<<<params.grid_size, params.block_size>>>(
+        device_pointers, interaction);
+    CELER_CUDA_CHECK_ERROR();
+}
+
+//---------------------------------------------------------------------------//
+} // namespace detail
+} // namespace celeritas

src/physics/em/detail/SeltzerBerger.hh

@@ -121,13 +121,14 @@ struct SeltzerBergerIds
 template<Ownership W, MemSpace M>
 struct SeltzerBergerData
 {
+    using MevMass = units::MevMass;
     //// MEMBER DATA ////
 
     //! IDs in a separate struct for readability/easier copying
     SeltzerBergerIds ids;
 
-    //! Electron mass [MevMass]
-    real_type electron_mass;
+    //! Electron mass [MeV / c^2]
+    MevMass electron_mass;
 
     // Differential cross section storage
     SeltzerBergerTableData<W, M> differential_xs;
@@ -137,7 +138,7 @@ struct SeltzerBergerData
     //! Whether the data is assigned
     explicit inline CELER_FUNCTION operator bool() const
     {
-        return ids && electron_mass > 0 && differential_xs;
+        return ids && electron_mass.value() > 0 && differential_xs;
     }
 
     //! Assign from another set of data
@@ -156,14 +157,16 @@ using SeltzerBergerDeviceRef
     = SeltzerBergerData<Ownership::const_reference, MemSpace::device>;
 using SeltzerBergerHostRef
     = SeltzerBergerData<Ownership::const_reference, MemSpace::host>;
+using SeltzerBergerNativeRef
+    = SeltzerBergerData<Ownership::const_reference, MemSpace::native>;
 
 //---------------------------------------------------------------------------//
 // KERNEL LAUNCHERS
 //---------------------------------------------------------------------------//
 
 // Launch the Seltzer-Berger interaction
 void seltzer_berger_interact(
-    const SeltzerBergerDeviceRef&              device_pointers,
+    const SeltzerBergerNativeRef&              shared,
     const ModelInteractRefs<MemSpace::device>& interaction);
 
 //---------------------------------------------------------------------------//

src/physics/em/detail/SeltzerBergerInteractor.hh

@@ -0,0 +1,95 @@
+//----------------------------------*-C++-*----------------------------------//
+// Copyright 2020 UT-Battelle, LLC, and other Celeritas developers.
+// See the top-level COPYRIGHT file for details.
+// SPDX-License-Identifier: (Apache-2.0 OR MIT)
+//---------------------------------------------------------------------------//
+//! \file SeltzerBergerInteractor.hh
+//---------------------------------------------------------------------------//
+#pragma once
+
+#include "base/Macros.hh"
+#include "base/Types.hh"
+#include "physics/base/CutoffView.hh"
+#include "physics/base/Interaction.hh"
+#include "physics/base/ParticleTrackView.hh"
+#include "physics/base/Secondary.hh"
+#include "base/StackAllocator.hh"
+#include "physics/base/Units.hh"
+#include "physics/material/ElementView.hh"
+#include "physics/material/MaterialView.hh"
+#include "SeltzerBerger.hh"
+#include "SBEnergyDistribution.hh"
+#include "SBPositronXsCorrector.hh"
+
+namespace celeritas
+{
+namespace detail
+{
+//---------------------------------------------------------------------------//
+/*!
+ * Seltzer-Berger model for electron and positron bremsstrahlung processes.
+ *
+ * Given an incoming electron or positron of sufficient energy (as per
+ * CutOffView), this class provides the energy loss of these particles due to
+ * radiation of photons in the field of a nucleus. This model improves accuracy
+ * using cross sections based on interpolation of published tables from Seltzer
+ * and Berger given in Nucl. Instr. and Meth. in Phys. Research B, 12(1):95–134
+ * (1985) and Atomic Data and Nuclear Data Tables, 35():345 (1986). The cross
+ * sections are obtained from SBEnergyDistribution and are appropriately scaled
+ * in the case of positrons via SBPositronXsCorrector (to be done).
+ *
+ * \note This interactor performs an analogous sampling as in Geant4's
+ * G4SeltzerBergerModel, documented in 10.2.1 of the Geant Physics Reference
+ * (release 10.6). The implementation is based on Geant4 10.4.3.
+ */
+class SeltzerBergerInteractor
+{
+  public:
+    //!@{
+    //! Type aliases
+    using Energy   = units::MevEnergy;
+    using EnergySq = SBEnergyDistribution::EnergySq;
+    using Momentum = units::MevMomentum;
+    //!@}
+
+  public:
+    //! Construct sampler from device/shared and state data
+    inline CELER_FUNCTION
+    SeltzerBergerInteractor(const SeltzerBergerNativeRef& shared,
+                            const ParticleTrackView&      particle,
+                            const Real3&                  inc_direction,
+                            const CutoffView&             cutoffs,
+                            StackAllocator<Secondary>&    allocate,
+                            const MaterialView&           material,
+                            const ElementId&              element_id);
+
+    // Sample an interaction with the given RNG
+    template<class Engine>
+    inline CELER_FUNCTION Interaction operator()(Engine& rng);
+
+  private:
+    // Device (host CPU or GPU device) references
+    const SeltzerBergerNativeRef& shared_;
+    // Type of particle
+    const ParticleId particle_id_;
+    // Incident particle energy
+    const Energy inc_energy_;
+    // Incident particle direction
+    const Momentum inc_momentum_;
+    // Incident particle direction
+    const Real3& inc_direction_;
+    // Interactor thresholds
+    const CutoffView& cutoffs_;
+    // Allocate space for a secondary particle
+    StackAllocator<Secondary>& allocate_;
+    // Material in which interaction occurs
+    const MaterialView& material_;
+    // Element in which interaction occurs
+    const ElementId& element_id_;
+};
+
+//---------------------------------------------------------------------------//
+} // namespace detail
+} // namespace celeritas
+
+#include "SeltzerBergerInteractor.i.hh"