Add helper functions for sampling/calculating direction in interactors (

#1301)

* Add helper function for calculating exiting direction in interactors
* Add helper to sample azimuthal angle, calculate Cartesian vector, and rotate direction
* Address review feedback
* Add PhysicsUtils test
* Rename to InteractionUtils and add assertions

src/celeritas/em/interactor/BetheHeitlerInteractor.hh

@@ -302,8 +302,8 @@ CELER_FUNCTION Interaction BetheHeitlerInteractor::operator()(Engine& rng)
 
     // Sample secondary directions.
     // Note that momentum is not exactly conserved.
-    UniformRealDistribution<real_type> sample_phi(0, 2 * constants::pi);
-    real_type phi = sample_phi(rng);
+    real_type phi
+        = UniformRealDistribution<real_type>(0, 2 * constants::pi)(rng);
 
     // Electron
     TsaiUrbanDistribution sample_electron_angle(secondaries[0].energy,

src/celeritas/em/interactor/CoulombScatteringInteractor.hh

@@ -19,6 +19,7 @@
 #include "celeritas/mat/MaterialView.hh"
 #include "celeritas/phys/CutoffView.hh"
 #include "celeritas/phys/Interaction.hh"
+#include "celeritas/phys/InteractionUtils.hh"
 #include "celeritas/phys/ParticleTrackView.hh"
 
 #include "detail/PhysicsConstants.hh"
@@ -143,9 +144,7 @@ CELER_FUNCTION Interaction CoulombScatteringInteractor::operator()(Engine& rng)
 
     // Sample the new direction
     real_type cos_theta = sample_angle_(rng);
-    UniformRealDistribution<real_type> sample_phi(0, 2 * constants::pi);
-    result.direction
-        = rotate(inc_direction_, from_spherical(cos_theta, sample_phi(rng)));
+    result.direction = ExitingDirectionSampler{cos_theta, inc_direction_}(rng);
 
     // Recoil energy is kinetic energy transfered to the atom
     real_type inc_energy = value_as<Energy>(particle_.energy());

src/celeritas/em/interactor/EPlusGGInteractor.hh

@@ -15,6 +15,7 @@
 #include "celeritas/Quantities.hh"
 #include "celeritas/em/data/EPlusGGData.hh"
 #include "celeritas/phys/Interaction.hh"
+#include "celeritas/phys/InteractionUtils.hh"
 #include "celeritas/phys/ParticleTrackView.hh"
 #include "celeritas/phys/Secondary.hh"
 #include "celeritas/random/distribution/BernoulliDistribution.hh"
@@ -154,19 +155,12 @@ CELER_FUNCTION Interaction EPlusGGInteractor::operator()(Engine& rng)
         real_type const eplus_moment = std::sqrt(inc_energy_ * total_energy);
 
         // Sample and save outgoing secondary data
-        UniformRealDistribution<real_type> sample_phi(0, 2 * constants::pi);
-
         secondaries[0].energy = Energy{gamma_energy};
         secondaries[0].direction
-            = rotate(from_spherical(cost, sample_phi(rng)), inc_direction_);
-
+            = ExitingDirectionSampler{cost, inc_direction_}(rng);
         secondaries[1].energy = Energy{total_energy - gamma_energy};
-        for (int i = 0; i < 3; ++i)
-        {
-            secondaries[1].direction[i] = eplus_moment * inc_direction_[i]
-                                          - inc_energy_ * inc_direction_[i];
-        }
-        secondaries[1].direction = make_unit_vector(secondaries[1].direction);
+        secondaries[1].direction = calc_exiting_direction(
+            {eplus_moment, inc_direction_}, {inc_energy_, inc_direction_});
     }
 
     return result;

src/celeritas/em/interactor/KleinNishinaInteractor.hh

@@ -15,6 +15,7 @@
 #include "celeritas/Quantities.hh"
 #include "celeritas/em/data/KleinNishinaData.hh"
 #include "celeritas/phys/Interaction.hh"
+#include "celeritas/phys/InteractionUtils.hh"
 #include "celeritas/phys/ParticleTrackView.hh"
 #include "celeritas/phys/Secondary.hh"
 #include "celeritas/random/distribution/BernoulliDistribution.hh"
@@ -163,10 +164,8 @@ CELER_FUNCTION Interaction KleinNishinaInteractor::operator()(Engine& rng)
     result.secondaries = {electron_secondary, 1};
 
     // Sample azimuthal direction and rotate the outgoing direction
-    UniformRealDistribution<real_type> sample_phi(0, 2 * constants::pi);
-    result.direction
-        = rotate(from_spherical(1 - one_minus_costheta, sample_phi(rng)),
-                 result.direction);
+    result.direction = ExitingDirectionSampler{1 - one_minus_costheta,
+                                               result.direction}(rng);
 
     // Construct secondary energy by neglecting electron binding energy
     electron_secondary->energy
@@ -183,14 +182,9 @@ CELER_FUNCTION Interaction KleinNishinaInteractor::operator()(Engine& rng)
     // Outgoing secondary is an electron
     electron_secondary->particle_id = shared_.ids.electron;
     // Calculate exiting electron direction via conservation of momentum
-    for (int i = 0; i < 3; ++i)
-    {
-        electron_secondary->direction[i]
-            = inc_direction_[i] * inc_energy_.value()
-              - result.direction[i] * result.energy.value();
-    }
     electron_secondary->direction
-        = make_unit_vector(electron_secondary->direction);
+        = calc_exiting_direction({inc_energy_.value(), inc_direction_},
+                                 {result.energy.value(), result.direction});
 
     return result;
 }

src/celeritas/em/interactor/LivermorePEInteractor.hh

@@ -19,9 +19,9 @@
 #include "celeritas/grid/PolyEvaluator.hh"
 #include "celeritas/phys/CutoffView.hh"
 #include "celeritas/phys/Interaction.hh"
+#include "celeritas/phys/InteractionUtils.hh"
 #include "celeritas/phys/ParticleTrackView.hh"
 #include "celeritas/phys/Secondary.hh"
-#include "celeritas/random/distribution/UniformRealDistribution.hh"
 
 #include "AtomicRelaxationHelper.hh"
 
@@ -351,8 +351,7 @@ CELER_FUNCTION Real3 LivermorePEInteractor::sample_direction(Engine& rng) const
 
     // Sample the azimuthal angle and calculate the direction of the
     // photoelectron
-    UniformRealDistribution<real_type> sample_phi(0, 2 * constants::pi);
-    return rotate(from_spherical(1 - nu, sample_phi(rng)), inc_direction_);
+    return ExitingDirectionSampler{1 - nu, inc_direction_}(rng);
 }
 
 //---------------------------------------------------------------------------//

src/celeritas/em/interactor/MollerBhabhaInteractor.hh

@@ -21,6 +21,7 @@
 #include "celeritas/em/distribution/MollerEnergyDistribution.hh"
 #include "celeritas/phys/CutoffView.hh"
 #include "celeritas/phys/Interaction.hh"
+#include "celeritas/phys/InteractionUtils.hh"
 #include "celeritas/phys/ParticleTrackView.hh"
 #include "celeritas/phys/Secondary.hh"
 #include "celeritas/random/distribution/BernoulliDistribution.hh"
@@ -164,39 +165,21 @@ CELER_FUNCTION Interaction MollerBhabhaInteractor::operator()(Engine& rng)
         = secondary_energy
           * (total_energy + value_as<Mass>(shared_.electron_mass))
           / (secondary_momentum * inc_momentum_);
-
-    secondary_cos_theta = celeritas::min<real_type>(secondary_cos_theta, 1);
     CELER_ASSERT(secondary_cos_theta >= -1 && secondary_cos_theta <= 1);
 
-    // Sample phi isotropically
-    UniformRealDistribution<real_type> sample_phi(0, 2 * constants::pi);
-
-    // Create cartesian direction from the sampled theta and phi
-    Real3 secondary_direction = rotate(
-        from_spherical(secondary_cos_theta, sample_phi(rng)), inc_direction_);
-
-    // Calculate incident particle final direction
-    Real3 inc_exiting_direction;
-    for (int i = 0; i < 3; ++i)
-    {
-        real_type inc_momentum_ijk = inc_momentum_ * inc_direction_[i];
-        real_type secondary_momentum_ijk = secondary_momentum
-                                           * secondary_direction[i];
-        inc_exiting_direction[i] = inc_momentum_ijk - secondary_momentum_ijk;
-    }
-    inc_exiting_direction = make_unit_vector(inc_exiting_direction);
+    // Assign values to the secondary particle
+    electron_secondary->particle_id = shared_.ids.electron;
+    electron_secondary->energy = Energy{secondary_energy};
+    electron_secondary->direction
+        = ExitingDirectionSampler{secondary_cos_theta, inc_direction_}(rng);
 
     // Construct interaction for change to primary (incident) particle
-    real_type const inc_exiting_energy = inc_energy_ - secondary_energy;
     Interaction result;
-    result.energy = Energy{inc_exiting_energy};
+    result.energy = Energy{inc_energy_ - secondary_energy};
     result.secondaries = {electron_secondary, 1};
-    result.direction = inc_exiting_direction;
-
-    // Assign values to the secondary particle
-    electron_secondary[0].particle_id = shared_.ids.electron;
-    electron_secondary[0].energy = Energy{secondary_energy};
-    electron_secondary[0].direction = secondary_direction;
+    result.direction = calc_exiting_direction(
+        {inc_momentum_, inc_direction_},
+        {secondary_momentum, electron_secondary->direction});
 
     return result;
 }

src/celeritas/em/interactor/MuBetheBlochInteractor.hh

@@ -19,6 +19,7 @@
 #include "celeritas/em/data/MuBetheBlochData.hh"
 #include "celeritas/phys/CutoffView.hh"
 #include "celeritas/phys/Interaction.hh"
+#include "celeritas/phys/InteractionUtils.hh"
 #include "celeritas/phys/ParticleTrackView.hh"
 #include "celeritas/phys/Secondary.hh"
 #include "celeritas/random/distribution/BernoulliDistribution.hh"
@@ -209,20 +210,16 @@ CELER_FUNCTION Interaction MuBetheBlochInteractor::operator()(Engine& rng)
     CELER_ASSERT(costheta <= 1);
 
     // Sample and save outgoing secondary data
-    UniformRealDistribution<real_type> sample_phi(0, 2 * constants::pi);
     secondary->direction
-        = rotate(from_spherical(costheta, sample_phi(rng)), inc_direction_);
+        = ExitingDirectionSampler{costheta, inc_direction_}(rng);
     secondary->energy = Energy{secondary_energy};
     secondary->particle_id = shared_.electron;
 
     Interaction result;
     result.energy = Energy{inc_energy_ - secondary_energy};
-    for (int i = 0; i < 3; ++i)
-    {
-        result.direction[i] = inc_momentum_ * inc_direction_[i]
-                              - secondary_momentum * secondary->direction[i];
-    }
-    result.direction = make_unit_vector(result.direction);
+    result.direction
+        = calc_exiting_direction({inc_momentum_, inc_direction_},
+                                 {secondary_momentum, secondary->direction});
     result.secondaries = {secondary, 1};
 
     return result;

src/celeritas/em/interactor/MuBremsstrahlungInteractor.hh

@@ -17,11 +17,11 @@
 #include "celeritas/mat/ElementView.hh"
 #include "celeritas/mat/MaterialView.hh"
 #include "celeritas/phys/Interaction.hh"
+#include "celeritas/phys/InteractionUtils.hh"
 #include "celeritas/phys/ParticleTrackView.hh"
 #include "celeritas/phys/Secondary.hh"
 #include "celeritas/random/distribution/BernoulliDistribution.hh"
 #include "celeritas/random/distribution/ReciprocalDistribution.hh"
-#include "celeritas/random/distribution/UniformRealDistribution.hh"
 
 namespace celeritas
 {
@@ -113,8 +113,8 @@ template<class Engine>
 CELER_FUNCTION Interaction MuBremsstrahlungInteractor::operator()(Engine& rng)
 {
     // Allocate space for gamma
-    Secondary* secondaries = allocate_(1);
-    if (secondaries == nullptr)
+    Secondary* secondary = allocate_(1);
+    if (secondary == nullptr)
     {
         // Failed to allocate space for a secondary
         return Interaction::from_failure();
@@ -137,32 +137,20 @@ CELER_FUNCTION Interaction MuBremsstrahlungInteractor::operator()(Engine& rng)
     } while (!BernoulliDistribution(
         epsilon * this->differential_cross_section(epsilon) / func_1)(rng));
 
-    // Sample secondary direction.
-    UniformRealDistribution<real_type> phi(0, 2 * constants::pi);
-
-    real_type cost = this->sample_cos_theta(epsilon, rng);
-    Real3 gamma_dir = rotate(from_spherical(cost, phi(rng)), inc_direction_);
-
-    Real3 inc_direction;
-    for (int i = 0; i < 3; ++i)
-    {
-        inc_direction[i] = value_as<Momentum>(particle_.momentum())
-                               * inc_direction_[i]
-                           - epsilon * gamma_dir[i];
-    }
-    inc_direction = make_unit_vector(inc_direction);
+    // Save outgoing secondary data
+    secondary->particle_id = shared_.ids.gamma;
+    secondary->energy = Energy{epsilon};
+    secondary->direction = ExitingDirectionSampler{
+        this->sample_cos_theta(epsilon, rng), inc_direction_}(rng);
 
     // Construct interaction for change to primary (incident) particle
     Interaction result;
     result.energy
         = units::MevEnergy{value_as<Energy>(particle_.energy()) - epsilon};
-    result.direction = inc_direction;
-    result.secondaries = {secondaries, 1};
-
-    // Save outgoing secondary data
-    secondaries[0].particle_id = shared_.ids.gamma;
-    secondaries[0].energy = units::MevEnergy{epsilon};
-    secondaries[0].direction = gamma_dir;
+    result.direction = calc_exiting_direction(
+        {value_as<Momentum>(particle_.momentum()), inc_direction_},
+        {epsilon, secondary->direction});
+    result.secondaries = {secondary, 1};
 
     return result;
 }

src/celeritas/em/interactor/RayleighInteractor.hh

@@ -16,6 +16,7 @@
 #include "celeritas/Types.hh"
 #include "celeritas/em/data/RayleighData.hh"
 #include "celeritas/phys/Interaction.hh"
+#include "celeritas/phys/InteractionUtils.hh"
 #include "celeritas/phys/ParticleTrackView.hh"
 #include "celeritas/random/Selector.hh"
 #include "celeritas/random/distribution/GenerateCanonical.hh"
@@ -147,11 +148,8 @@ CELER_FUNCTION Interaction RayleighInteractor::operator()(Engine& rng)
 
     } while (2 * generate_canonical(rng) > 1 + ipow<2>(cost) || cost < -1);
 
-    UniformRealDistribution<real_type> sample_phi(0, 2 * constants::pi);
-
     // Scattered direction
-    result.direction
-        = rotate(from_spherical(cost, sample_phi(rng)), inc_direction_);
+    result.direction = ExitingDirectionSampler{cost, inc_direction_}(rng);
 
     CELER_ENSURE(result.action == Interaction::Action::scattered);
     return result;

src/celeritas/em/interactor/detail/BremFinalStateHelper.hh

@@ -9,6 +9,7 @@
 
 #include "celeritas/Quantities.hh"
 #include "celeritas/em/distribution/TsaiUrbanDistribution.hh"
+#include "celeritas/phys/InteractionUtils.hh"
 #include "celeritas/phys/Secondary.hh"
 #include "celeritas/random/distribution/UniformRealDistribution.hh"
 
@@ -93,29 +94,19 @@ BremFinalStateHelper::BremFinalStateHelper(Energy const& inc_energy,
 template<class Engine>
 CELER_FUNCTION Interaction BremFinalStateHelper::operator()(Engine& rng)
 {
-    // Construct interaction for change to parent (incoming) particle
-    Interaction result;
-    result.energy = exit_energy_;
-
     // Generate exiting gamma direction from isotropic azimuthal angle and
     // TsaiUrbanDistribution for polar angle (based on G4ModifiedTsai)
-    UniformRealDistribution<real_type> sample_phi(0, 2 * constants::pi);
-    real_type cost = sample_polar_angle_(rng);
-    secondary_->direction
-        = rotate(from_spherical(cost, sample_phi(rng)), inc_direction_);
-
-    // Update parent particle direction
-    for (int i = 0; i < 3; ++i)
-    {
-        real_type inc_momentum_i = inc_momentum_.value() * inc_direction_[i];
-        real_type gamma_momentum_i = gamma_energy_.value()
-                                     * secondary_->direction[i];
-        result.direction[i] = inc_momentum_i - gamma_momentum_i;
-    }
-    result.direction = make_unit_vector(result.direction);
-
+    secondary_->direction = ExitingDirectionSampler{sample_polar_angle_(rng),
+                                                    inc_direction_}(rng);
     secondary_->particle_id = gamma_id_;
     secondary_->energy = gamma_energy_;
+
+    // Construct interaction for change to parent (incoming) particle
+    Interaction result;
+    result.energy = exit_energy_;
+    result.direction = calc_exiting_direction(
+        {inc_momentum_.value(), inc_direction_},
+        {gamma_energy_.value(), secondary_->direction});
     result.secondaries = {secondary_, 1};
 
     return result;