interactor/detail/CascadeCollider.hh

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