interactor/detail/MucfInteractorUtils.hh

0001 //------------------------------- -*- C++ -*- -------------------------------//
0002 // Copyright Celeritas contributors: see top-level COPYRIGHT file for details
0003 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0004 //---------------------------------------------------------------------------//
0005 //! \file celeritas/mucf/interactor/detail/MucfInteractorUtils.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/Macros.hh"
0010 #include "corecel/Types.hh"
0011 #include "corecel/math/Algorithms.hh"
0012 #include "corecel/math/ArrayUtils.hh"
0013 #include "corecel/random/distribution/GenerateCanonical.hh"
0014 #include "corecel/random/distribution/IsotropicDistribution.hh"
0015 #include "celeritas/Quantities.hh"
0016 #include "celeritas/Units.hh"
0017 #include "celeritas/grid/NonuniformGridCalculator.hh"
0018 #include "celeritas/phys/Secondary.hh"
0019
0020 namespace celeritas
0021 {
0022 namespace detail
0023 {
0024 //---------------------------------------------------------------------------//
0025 /*!
0026  * Calculate momentum magnitude from particle energy and mass via
0027  * \f$ p = \sqrt{K^2 + 2mK} \f$ .
0028  */
0029 inline CELER_FUNCTION real_type calc_momentum(units::MevEnergy const energy,
0030                                               units::MevMass const mass)
0031 {
0032     return std::sqrt(ipow<2>(value_as<units::MevEnergy>(energy))
0033                      + 2 * value_as<units::MevMass>(mass)
0034                            * value_as<units::MevEnergy>(energy));
0035 }
0036
0037 //---------------------------------------------------------------------------//
0038 /*!
0039  * Calculate kinetic energy given a particle's momentum and mass via
0040  * \f$ K = \sqrt{p^2 - m^2} - m\f$ .
0041  */
0042 inline CELER_FUNCTION units::MevEnergy calc_kinetic_energy(
0043     Real3 const& momentum_vec, units::MevMass const mass)
0044 {
0045     real_type momentum_mag = norm(momentum_vec);
0046     return units::MevEnergy{std::sqrt(ipow<2>(momentum_mag)
0047                                       + ipow<2>(value_as<units::MevMass>(mass)))
0048                             - value_as<units::MevMass>(mass)};
0049 }
0050
0051 //---------------------------------------------------------------------------//
0052 /*!
0053  * Sample muCF secondary with isotropic random direction and known outgoing
0054  * energy.
0055  */
0056 template<class Engine>
0057 inline CELER_FUNCTION Secondary sample_mucf_secondary(
0058     ParticleId pid, units::MevEnergy energy, Engine& rng)
0059 {
0060     Secondary result;
0061     result.particle_id = pid;
0062     result.energy = energy;
0063     result.direction = IsotropicDistribution<real_type>()(rng);
0064     return result;
0065 }
0066
0067 //---------------------------------------------------------------------------//
0068 /*!
0069  * Sample a muon secondary with isotropic random direction and energy from the
0070  * provided CDF.
0071  *
0072  * \note The muon grid range is [0,1) and its domain is the energy in MeV.
0073  */
0074 template<class Engine>
0075 inline CELER_FUNCTION Secondary sample_mucf_muon(
0076     ParticleId pid, NonuniformGridCalculator sample_energy, Engine& rng)
0077 {
0078     Secondary result;
0079     result.particle_id = pid;
0080     result.energy = units::MevEnergy{sample_energy(generate_canonical(rng))};
0081     result.direction = IsotropicDistribution<real_type>()(rng);
0082     return result;
0083 }
0084
0085 //---------------------------------------------------------------------------//
0086 /*!
0087  * Return opposite direction.
0088  */
0089 inline CELER_FUNCTION Real3 opposite(Real3 const& vec)
0090 {
0091     Real3 result;
0092     for (auto i : range(3))
0093     {
0094         result[i] = -vec[i];
0095     }
0096     return result;
0097 }
0098
0099 //---------------------------------------------------------------------------//
0100 /*!
0101  * Return the momentum of the third secondary from an at-rest three-body
0102  * sampling, once the other two are known.
0103  *
0104  * \todo This may be expanded to do a full three-body energy + momentum
0105  * conservation.
0106  */
0107 inline CELER_FUNCTION Secondary calc_third_secondary(
0108     Secondary sec_a,
0109     units::MevMass const mass_a,
0110     Secondary sec_b,
0111     units::MevMass const mass_b,
0112     ParticleId pid_c,
0113     units::MevMass const mass_c)
0114 {
0115     auto const momentum_a_mag = calc_momentum(sec_a.energy, mass_a);
0116     auto const momentum_b_mag = calc_momentum(sec_b.energy, mass_b);
0117
0118     // Calculate direction from momentum conservation: p3 = - (p1 + p2)
0119     Real3 momentum_vec;
0120     for (auto i : range(3))
0121     {
0122         momentum_vec[i] = -(sec_a.direction[i] * momentum_a_mag
0123                             + sec_b.direction[i] * momentum_b_mag);
0124     }
0125
0126     Secondary result;
0127     result.particle_id = pid_c;
0128     result.direction = make_unit_vector(momentum_vec);
0129     result.energy = calc_kinetic_energy(momentum_vec, mass_c);
0130     return result;
0131 }
0132
0133 //---------------------------------------------------------------------------//
0134 }  // namespace detail
0135 }  // namespace celeritas