interactor/detail/CascadeCollider.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/neutron/interactor/detail/CascadeCollider.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/Macros.hh"
0010 #include "corecel/Types.hh"
0011 #include "corecel/grid/NonuniformGrid.hh"
0012 #include "corecel/grid/TwodGridCalculator.hh"
0013 #include "corecel/grid/TwodGridData.hh"
0014 #include "corecel/math/Algorithms.hh"
0015 #include "corecel/math/ArrayOperators.hh"
0016 #include "corecel/math/ArrayUtils.hh"
0017 #include "corecel/random/distribution/GenerateCanonical.hh"
0018 #include "corecel/random/distribution/UniformRealDistribution.hh"
0019 #include "celeritas/Quantities.hh"
0020 #include "celeritas/Types.hh"
0021 #include "celeritas/grid/InverseCdfFinder.hh"
0022 #include "celeritas/phys/FourVector.hh"
0023
0024 #include "CascadeParticle.hh"
0025
0026 namespace celeritas
0027 {
0028 namespace detail
0029 {
0030 //---------------------------------------------------------------------------//
0031 /*!
0032  * Sample final state for a nucleon-nucleon intra-nucleus cascade collision.
0033  *
0034  * This samples the final state of outgoing particles from the two-body
0035  * intra-nucleus nucleon-nucleon collision in the center of mass (c.m.) frame
0036  * and returns them after converting momentum to the lab frame. It performs the
0037  * same sampling routine as in Geant4's \c G4ElementaryParticleCollider, mainly
0038  * implemented in collide and generateSCMfinalState methods. The
0039  * \f$ \cos\theta \f$ distribution in c.m. is inversely sampled using the
0040  * tabulated cumulative distribution function (c.d.f) data in the kinetic
0041  * energy and the cosine bins which are implemented in \c
0042  * G4CascadeFinalStateAlgorithm::GenerateTwoBody and
0043  * \c  G4NumIntTwoBodyAngDst::GetCosTheta methods.
0044  */
0045 class CascadeCollider
0046 {
0047   public:
0048     //!@{
0049     //! \name Type aliases
0050     using FinalState = Array<CascadeParticle, 2>;
0051     //!@}
0052
0053   public:
0054     // Construct with shared data and colliding particles
0055     inline CELER_FUNCTION CascadeCollider(NeutronInelasticRef const& shared,
0056                                           CascadeParticle const& bullet,
0057                                           CascadeParticle const& target);
0058
0059     // Sample the final state of the two-body intra-nucleus cascade collision
0060     template<class Engine>
0061     inline CELER_FUNCTION FinalState operator()(Engine& rng);
0062
0063   private:
0064     //// TYPES ////
0065
0066     using Mass = units::MevMass;
0067     using Momentum = units::MevMomentum;
0068     using Grid = NonuniformGrid<real_type>;
0069     using UniformRealDist = UniformRealDistribution<real_type>;
0070
0071     //// DATA ////
0072
0073     // Shared constant data
0074     NeutronInelasticRef const& shared_;
0075     // Participating cascade particles
0076     CascadeParticle const& bullet_;
0077     CascadeParticle const& target_;
0078
0079     // Id for intra-nuclear channels
0080     ChannelId ch_id_;
0081     // Boost vector in the center of mass frame [1/c]
0082     Real3 cm_velocity_;
0083     // Momentum magnitude in the center of mass frame [MeV/c]
0084     real_type cm_p_;
0085     // Kinetic energy in the target rest frame [MeV]
0086     real_type kin_energy_;
0087
0088     // Sampler
0089     UniformRealDist sample_phi_;
0090
0091     //// CONSTANTS ////
0092
0093     //! A criteria [1/c] for coplanarity in the Lorentz transformation
0094     static CELER_CONSTEXPR_FUNCTION real_type epsilon()
0095     {
0096         return real_type{1e-10};
0097     }
0098
0099     //// HELPER FUNCTIONS ////
0100
0101     //! Calculate the momentum magnitude in the center of mass frame
0102     inline CELER_FUNCTION real_type calc_cm_p(FourVector const& v) const;
0103 };
0104
0105 //---------------------------------------------------------------------------//
0106 // Inline DEFINITIONS
0107 //---------------------------------------------------------------------------//
0108 /*!
0109  * Construct from share data and colliding nucleons in the lab frame.
0110  */
0111 CELER_FUNCTION
0112 CascadeCollider::CascadeCollider(NeutronInelasticRef const& shared,
0113                                  CascadeParticle const& bullet,
0114                                  CascadeParticle const& target)
0115     : shared_(shared)
0116     , bullet_(bullet)
0117     , target_(target)
0118     , ch_id_(ChannelId{
0119           static_cast<size_type>((bullet_.type == target_.type) ? 0 : 1)})
0120     , sample_phi_(0, static_cast<real_type>(2 * constants::pi))
0121 {
0122     // Initialize the boost velocity and momentum in the center of mass frame
0123     FourVector sum_four_vec = bullet_.four_vec + target_.four_vec;
0124     cm_velocity_ = boost_vector(sum_four_vec);
0125     cm_p_ = this->calc_cm_p(sum_four_vec);
0126
0127     // Calculate the kinetic energy in the target rest frame
0128     FourVector bullet_p = bullet_.four_vec;
0129     boost(-boost_vector(target_.four_vec), &bullet_p);
0130     kin_energy_ = bullet_p.energy - norm(bullet_p);
0131 }
0132
0133 //---------------------------------------------------------------------------//
0134 /*!
0135  * Sample using the given RNG.
0136  */
0137 template<class Engine>
0138 CELER_FUNCTION auto CascadeCollider::operator()(Engine& rng) -> FinalState
0139 {
0140     // Sample cos\theta of outgoing particles in the center of mass frame
0141     real_type cdf = generate_canonical(rng);
0142     TwodGridData const& cdf_grid = shared_.angular_cdf[ch_id_];
0143     Grid energy_grid(cdf_grid.x, shared_.reals);
0144     real_type cos_theta{0};
0145
0146     if (kin_energy_ < energy_grid.back())
0147     {
0148         // Find cos\theta from tabulated angular data for a given CDF
0149         cos_theta = InverseCdfFinder(
0150             Grid(cdf_grid.y, shared_.reals),
0151             TwodGridCalculator(cdf_grid, shared_.reals)(kin_energy_))(cdf);
0152     }
0153     else
0154     {
0155         // Sample the angle outside tabulated data (unlikely)
0156         real_type slope = 2 * energy_grid.back() * ipow<2>(cm_p_);
0157         cos_theta = std::log(1 + cdf * std::expm1(2 * slope)) / slope - 1;
0158     }
0159
0160     // Sample the momentum of outgoing particles in the center of mass frame
0161     // Rotate the momentum along the reference z-axis
0162     auto fv = FourVector::from_mass_momentum(
0163         bullet_.mass,
0164         Momentum{cm_p_},
0165         from_spherical(cos_theta, sample_phi_(rng)));
0166
0167     // Find the final state of outgoing particles
0168     FinalState result = {bullet_, target_};
0169
0170     FourVector cm_momentum = target_.four_vec;
0171     boost(-cm_velocity_, &cm_momentum);
0172
0173     Real3 cm_dir = make_unit_vector(-cm_momentum.mom);
0174     real_type vel_parallel = dot_product(cm_velocity_, cm_dir);
0175
0176     Real3 vscm = cm_velocity_;
0177     axpy(-vel_parallel, cm_dir, &vscm);
0178
0179     if (norm(vscm) > this->epsilon())
0180     {
0181         vscm = make_unit_vector(vscm);
0182         Real3 vxcm = make_unit_vector(cross_product(cm_dir, cm_velocity_));
0183         if (norm(vxcm) > this->epsilon())
0184         {
0185             for (int i = 0; i < 3; ++i)
0186             {
0187                 result[0].four_vec.mom[i] = fv.mom[0] * vscm[i]
0188                                             + fv.mom[1] * vxcm[i]
0189                                             + fv.mom[2] * cm_dir[i];
0190             }
0191             result[0].four_vec.energy = fv.energy;
0192         }
0193     }
0194     else
0195     {
0196         // Degenerate if velocity perpendicular to the c.m. momentum is small
0197         result[0].four_vec = fv;
0198     }
0199
0200     result[1].four_vec = {{-result[0].four_vec.mom},
0201                           hypot(cm_p_, value_as<Mass>(target_.mass))};
0202
0203     // Convert the final state to the lab frame
0204     for (auto i : range(2))
0205     {
0206         boost(cm_velocity_, &result[i].four_vec);
0207     }
0208
0209     return result;
0210 }
0211
0212 //---------------------------------------------------------------------------//
0213 /*!
0214  * Calculate the momentum magnitude of outgoing particles in the center of mass
0215  * (c.m.) frame. See Geant4's G4VHadDecayAlgorithm::TwoBodyMomentum method.
0216  */
0217 CELER_FUNCTION real_type CascadeCollider::calc_cm_p(FourVector const& v) const
0218 {
0219     // The total energy in c.m.
0220     real_type m0 = norm(v);
0221
0222     real_type m1 = value_as<Mass>(bullet_.mass);
0223     real_type m2 = value_as<Mass>(target_.mass);
0224
0225     real_type pc_sq = diffsq(m0, m1 - m2) * diffsq(m0, m1 + m2);
0226
0227     return std::sqrt(clamp_to_nonneg(pc_sq)) / (2 * m0);
0228 }
0229
0230 //---------------------------------------------------------------------------//
0231 }  // namespace detail
0232 }  // namespace celeritas