From b935f361bc010afe783d5d141fd1e27407cf813b Mon Sep 17 00:00:00 2001 From: Amanda Lund Date: Mon, 1 Jul 2024 16:11:08 -0500 Subject: [PATCH] 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 --- .../em/interactor/BetheHeitlerInteractor.hh | 4 +- .../interactor/CoulombScatteringInteractor.hh | 5 +- .../em/interactor/EPlusGGInteractor.hh | 14 ++-- .../em/interactor/KleinNishinaInteractor.hh | 16 ++--- .../em/interactor/LivermorePEInteractor.hh | 5 +- .../em/interactor/MollerBhabhaInteractor.hh | 37 +++------- .../em/interactor/MuBetheBlochInteractor.hh | 13 ++-- .../interactor/MuBremsstrahlungInteractor.hh | 36 ++++------ .../em/interactor/RayleighInteractor.hh | 6 +- .../interactor/detail/BremFinalStateHelper.hh | 29 +++----- src/celeritas/phys/InteractionUtils.hh | 67 ++++++++++++++++++ test/celeritas/CMakeLists.txt | 1 + test/celeritas/phys/InteractionUtils.test.cc | 69 +++++++++++++++++++ 13 files changed, 191 insertions(+), 111 deletions(-) create mode 100644 src/celeritas/phys/InteractionUtils.hh create mode 100644 test/celeritas/phys/InteractionUtils.test.cc diff --git a/src/celeritas/em/interactor/BetheHeitlerInteractor.hh b/src/celeritas/em/interactor/BetheHeitlerInteractor.hh index fd6e5c6cb0..1e707392be 100644 --- a/src/celeritas/em/interactor/BetheHeitlerInteractor.hh +++ b/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 sample_phi(0, 2 * constants::pi); - real_type phi = sample_phi(rng); + real_type phi + = UniformRealDistribution(0, 2 * constants::pi)(rng); // Electron TsaiUrbanDistribution sample_electron_angle(secondaries[0].energy, diff --git a/src/celeritas/em/interactor/CoulombScatteringInteractor.hh b/src/celeritas/em/interactor/CoulombScatteringInteractor.hh index 89929a9ea9..bea6e8d178 100644 --- a/src/celeritas/em/interactor/CoulombScatteringInteractor.hh +++ b/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 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(particle_.energy()); diff --git a/src/celeritas/em/interactor/EPlusGGInteractor.hh b/src/celeritas/em/interactor/EPlusGGInteractor.hh index a03168a7eb..4e67d2e8b8 100644 --- a/src/celeritas/em/interactor/EPlusGGInteractor.hh +++ b/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 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; diff --git a/src/celeritas/em/interactor/KleinNishinaInteractor.hh b/src/celeritas/em/interactor/KleinNishinaInteractor.hh index 6ca76fb12b..8c31ceb89a 100644 --- a/src/celeritas/em/interactor/KleinNishinaInteractor.hh +++ b/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 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; } diff --git a/src/celeritas/em/interactor/LivermorePEInteractor.hh b/src/celeritas/em/interactor/LivermorePEInteractor.hh index eaa753f6a1..ba3cce14c8 100644 --- a/src/celeritas/em/interactor/LivermorePEInteractor.hh +++ b/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 sample_phi(0, 2 * constants::pi); - return rotate(from_spherical(1 - nu, sample_phi(rng)), inc_direction_); + return ExitingDirectionSampler{1 - nu, inc_direction_}(rng); } //---------------------------------------------------------------------------// diff --git a/src/celeritas/em/interactor/MollerBhabhaInteractor.hh b/src/celeritas/em/interactor/MollerBhabhaInteractor.hh index fa7e9feafe..af21ef3e80 100644 --- a/src/celeritas/em/interactor/MollerBhabhaInteractor.hh +++ b/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(shared_.electron_mass)) / (secondary_momentum * inc_momentum_); - - secondary_cos_theta = celeritas::min(secondary_cos_theta, 1); CELER_ASSERT(secondary_cos_theta >= -1 && secondary_cos_theta <= 1); - // Sample phi isotropically - UniformRealDistribution 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; } diff --git a/src/celeritas/em/interactor/MuBetheBlochInteractor.hh b/src/celeritas/em/interactor/MuBetheBlochInteractor.hh index 6ec2b1577d..4bd000a0cb 100644 --- a/src/celeritas/em/interactor/MuBetheBlochInteractor.hh +++ b/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 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; diff --git a/src/celeritas/em/interactor/MuBremsstrahlungInteractor.hh b/src/celeritas/em/interactor/MuBremsstrahlungInteractor.hh index 5db43b4124..8ae9baeb10 100644 --- a/src/celeritas/em/interactor/MuBremsstrahlungInteractor.hh +++ b/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 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 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(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(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(particle_.momentum()), inc_direction_}, + {epsilon, secondary->direction}); + result.secondaries = {secondary, 1}; return result; } diff --git a/src/celeritas/em/interactor/RayleighInteractor.hh b/src/celeritas/em/interactor/RayleighInteractor.hh index 36899f84e7..2991952932 100644 --- a/src/celeritas/em/interactor/RayleighInteractor.hh +++ b/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 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; diff --git a/src/celeritas/em/interactor/detail/BremFinalStateHelper.hh b/src/celeritas/em/interactor/detail/BremFinalStateHelper.hh index 823f9a5408..92c12b40a9 100644 --- a/src/celeritas/em/interactor/detail/BremFinalStateHelper.hh +++ b/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 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 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; diff --git a/src/celeritas/phys/InteractionUtils.hh b/src/celeritas/phys/InteractionUtils.hh new file mode 100644 index 0000000000..b24975122f --- /dev/null +++ b/src/celeritas/phys/InteractionUtils.hh @@ -0,0 +1,67 @@ +//----------------------------------*-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/phys/InteractionUtils.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/Types.hh" +#include "corecel/math/ArrayUtils.hh" +#include "geocel/Types.hh" +#include "celeritas/Constants.hh" +#include "celeritas/random/distribution/UniformRealDistribution.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +//! Particle momenntum +struct Momentum +{ + real_type magnitude; + Real3 const& direction; +}; + +//---------------------------------------------------------------------------// +/*! + * Calculate exiting direction via conservation of momentum. + */ +inline CELER_FUNCTION Real3 calc_exiting_direction(Momentum inc_momentum, + Momentum out_momentum) +{ + CELER_EXPECT(inc_momentum.magnitude > 0); + CELER_EXPECT(out_momentum.magnitude > 0); + + Real3 result; + for (int i = 0; i < 3; ++i) + { + result[i] = inc_momentum.direction[i] * inc_momentum.magnitude + - out_momentum.direction[i] * out_momentum.magnitude; + } + return make_unit_vector(result); +} + +//---------------------------------------------------------------------------// +/*! + * Sample an exiting direction from a polar cosine and incident direction. + * + * Combine an already-sampled change in polar cosine (dot product of incident + * and exiting) with a sampled uniform azimuthal direction, and apply that + * rotation to the original track's incident direction. + */ +struct ExitingDirectionSampler +{ + real_type costheta; + Real3 const& direction; + + template + inline CELER_FUNCTION Real3 operator()(Engine& rng) + { + UniformRealDistribution sample_phi(0, 2 * constants::pi); + return rotate(from_spherical(costheta, sample_phi(rng)), direction); + } +}; + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/test/celeritas/CMakeLists.txt b/test/celeritas/CMakeLists.txt index b403f34752..61ed2c0f9d 100644 --- a/test/celeritas/CMakeLists.txt +++ b/test/celeritas/CMakeLists.txt @@ -298,6 +298,7 @@ celeritas_add_test(neutron/NeutronInelastic.test.cc ${_needs_double}) celeritas_add_test(phys/CutoffParams.test.cc) celeritas_add_device_test(phys/Particle) celeritas_add_device_test(phys/Physics) +celeritas_add_test(phys/InteractionUtils.test.cc) celeritas_add_test(phys/PhysicsStepUtils.test.cc) celeritas_add_test(phys/PrimaryGenerator.test.cc LINK_LIBRARIES nlohmann_json::nlohmann_json) diff --git a/test/celeritas/phys/InteractionUtils.test.cc b/test/celeritas/phys/InteractionUtils.test.cc new file mode 100644 index 0000000000..97fa13f495 --- /dev/null +++ b/test/celeritas/phys/InteractionUtils.test.cc @@ -0,0 +1,69 @@ +//----------------------------------*-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/phys/InteractionUtils.test.cc +//---------------------------------------------------------------------------// +#include "celeritas/phys/InteractionUtils.hh" + +#include +#include +#include + +#include "TestMacros.hh" +#include "celeritas_test.hh" + +namespace celeritas +{ +namespace test +{ +//---------------------------------------------------------------------------// +// TESTS +//---------------------------------------------------------------------------// + +TEST(InteractionUtilsTest, calc_exiting_direction) +{ + Real3 inc_dir = {1, 0, 0}; + Real3 out_dir = {0, 1, 0}; + + EXPECT_VEC_SOFT_EQ(Real3({real_type(1.0 / std::sqrt(2)), + real_type(-1.0 / std::sqrt(2)), + 0}), + calc_exiting_direction({10, inc_dir}, {10, out_dir})); + + EXPECT_VEC_SOFT_EQ(Real3({real_type(1.0 / std::sqrt(101)), + real_type(-10.0 / std::sqrt(101)), + 0}), + calc_exiting_direction({1, inc_dir}, {10, out_dir})); +} + +TEST(InteractionUtilsTest, sample_exiting_direction) +{ + std::mt19937 rng; + Real3 inc_dir = make_unit_vector(Real3{1, 2, 3}); + + std::vector out_dir; + for (real_type costheta : {-1.0, 0.9, 0.1, 0.0, 0.1, 0.9, 1.0}) + { + Real3 result = ExitingDirectionSampler{costheta, inc_dir}(rng); + EXPECT_SOFT_EQ(costheta, dot_product(result, inc_dir)); + out_dir.insert(out_dir.end(), result.begin(), result.end()); + } + if (CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE) + { + static double const expected_out_dir[] + = {-0.26726124191242, -0.53452248382485, -0.80178372573727, + 0.65567203926594, 0.47242330622799, 0.58899099879154, + 0.54966700236953, 0.66690027070903, -0.5031006017034, + -0.81475551489392, 0.56962591956018, -0.10816544140881, + -0.93194807728444, 0.21402106250404, 0.29269056365125, + 0.20505130592048, 0.12514344585105, 0.97071781682466, + 0.26726124191242, 0.53452248382485, 0.80178372573727}; + EXPECT_VEC_SOFT_EQ(expected_out_dir, out_dir); + } +} + +//---------------------------------------------------------------------------// +} // namespace test +} // namespace celeritas