neutron/interactor/ChipsNeutronElasticInteractor.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/ChipsNeutronElasticInteractor.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/Macros.hh"
0011 #include "corecel/Types.hh"
0012 #include "corecel/math/ArrayOperators.hh"
0013 #include "celeritas/Constants.hh"
0014 #include "celeritas/Quantities.hh"
0015 #include "celeritas/mat/IsotopeView.hh"
0016 #include "celeritas/neutron/data/NeutronElasticData.hh"
0017 #include "celeritas/phys/FourVector.hh"
0018 #include "celeritas/phys/Interaction.hh"
0019 #include "celeritas/phys/ParticleTrackView.hh"
0020 #include "celeritas/random/distribution/UniformRealDistribution.hh"
0021
0022 #include "detail/MomentumTransferSampler.hh"
0023
0024 namespace celeritas
0025 {
0026 //---------------------------------------------------------------------------//
0027 /*!
0028  * Perform neutron elastic scattering based on the CHIPS (Chiral invariant
0029  * phase space) model.
0030  *
0031  * \note This performs the sampling procedure as in G4HadronElastic,
0032  * G4ChipsElasticModel and G4ChipsNeutronElasticXS, as partly documented
0033  * in section 21.1.3 of the Geant4 Physics Reference (release 11.2).
0034  */
0035 class ChipsNeutronElasticInteractor
0036 {
0037   public:
0038     //!@{
0039     //! \name Type aliases
0040     using Energy = units::MevEnergy;
0041     using Mass = units::MevMass;
0042     using Momentum = units::MevMomentum;
0043     //!@}
0044
0045   public:
0046     // Construct from shared and state data
0047     inline CELER_FUNCTION
0048     ChipsNeutronElasticInteractor(NeutronElasticRef const& shared,
0049                                   ParticleTrackView const& particle,
0050                                   Real3 const& inc_direction,
0051                                   IsotopeView const& target);
0052
0053     // Sample an interaction with the given RNG
0054     template<class Engine>
0055     inline CELER_FUNCTION Interaction operator()(Engine& rng);
0056
0057   private:
0058     //// TYPES ////
0059
0060     using UniformRealDist = UniformRealDistribution<real_type>;
0061
0062     //// DATA ////
0063
0064     // Constant shared data
0065     NeutronElasticRef const& shared_;
0066     // Incident neutron direction
0067     Real3 const& inc_direction_;
0068     // Target nucleus
0069     IsotopeView const& target_;
0070
0071     // Values of neutron mass (MevMass) and energy (MevEnergy)
0072     real_type neutron_mass_;
0073     real_type neutron_energy_;
0074     Momentum neutron_p_;
0075
0076     // Sampler
0077     UniformRealDist sample_phi_;
0078     detail::MomentumTransferSampler sample_momentum_square_;
0079 };
0080
0081 //---------------------------------------------------------------------------//
0082 // INLINE DEFINITIONS
0083 //---------------------------------------------------------------------------//
0084 /*!
0085  * Construct with shared and state data, and a target nucleus.
0086  */
0087 CELER_FUNCTION ChipsNeutronElasticInteractor::ChipsNeutronElasticInteractor(
0088     NeutronElasticRef const& shared,
0089     ParticleTrackView const& particle,
0090     Real3 const& inc_direction,
0091     IsotopeView const& target)
0092     : shared_(shared)
0093     , inc_direction_(inc_direction)
0094     , target_(target)
0095     , neutron_mass_(value_as<Mass>(shared_.neutron_mass))
0096     , neutron_energy_(neutron_mass_ + value_as<Energy>(particle.energy()))
0097     , neutron_p_(particle.momentum())
0098     , sample_phi_(0, 2 * constants::pi)
0099     , sample_momentum_square_(shared_, target_, neutron_p_)
0100 {
0101     CELER_EXPECT(particle.particle_id() == shared_.neutron);
0102 }
0103 //---------------------------------------------------------------------------//
0104 /*!
0105  * Sample the final state of the neutron-nucleus elastic scattering.
0106  *
0107  * The scattering angle (\f$ cos \theta \f$) in the elastic neutron-nucleus
0108  * scattering is expressed in terms of the momentum transfer (\f$ Q^{2} \f$),
0109  * \f[
0110  *  cos \theta = 1 - \frac{Q^{2}}{2 |\vec{k}_i|^{2}}
0111  * \f]
0112  * where \f$ \vec{k}_i \f$ is the momentum of the incident neutron in the
0113  * center of mass frame and the momentum transfer (\f$ Q^{2} \f$) is calculated
0114  * according to the CHIPS (Chiral Invariant Phase Space) model (see references
0115  * in model/ChipsNeutronElasticModel.cc and detail/MomentumTransferSampler.hh).
0116  * The final direction of the scattered neutron in the laboratory frame is
0117  * then transformed by the Lorentz boost of the initial four vector of the
0118  * neutron-nucleus system.
0119  */
0120 template<class Engine>
0121 CELER_FUNCTION Interaction ChipsNeutronElasticInteractor::operator()(Engine& rng)
0122 {
0123     // Scattered neutron with respect to the axis of incident direction
0124     Interaction result;
0125
0126     // The momentum magnitude in the c.m. frame
0127     real_type target_mass = value_as<Mass>(target_.nuclear_mass());
0128
0129     real_type cm_p = value_as<units::MevMomentum>(neutron_p_)
0130                      / std::sqrt(1 + ipow<2>(neutron_mass_ / target_mass)
0131                                  + 2 * neutron_energy_ / target_mass);
0132
0133     // Sample the scattered direction from the invariant momentum transfer
0134     // squared (\f$ -t = Q^{2} \f$) in the c.m. frame
0135     real_type cos_theta
0136         = 1 - real_type(0.5) * sample_momentum_square_(rng) / ipow<2>(cm_p);
0137     CELER_ASSERT(std::fabs(cos_theta) <= 1);
0138
0139     // Boost to the center of mass (c.m.) frame
0140     Real3 cm_mom = cm_p * from_spherical(cos_theta, sample_phi_(rng));
0141     FourVector nlv1(
0142         {cm_mom, std::sqrt(ipow<2>(cm_p) + ipow<2>(neutron_mass_))});
0143
0144     FourVector lv({{0, 0, value_as<units::MevMomentum>(neutron_p_)},
0145                    neutron_energy_ + target_mass});
0146     boost(boost_vector(lv), &nlv1);
0147
0148     result.direction = rotate(make_unit_vector(nlv1.mom), inc_direction_);
0149
0150     // Kinetic energy of the scattered neutron and the recoiled nucleus
0151     lv.energy -= nlv1.energy;
0152     result.energy = Energy(nlv1.energy - neutron_mass_);
0153     real_type recoil_energy = clamp_to_nonneg(lv.energy - target_mass);
0154     result.energy_deposition = Energy{recoil_energy};
0155
0156     /*!
0157      * \todo Create a secondary ion(Z, N) if recoil_energy > recoil_threshold
0158      * with energy = recoil_energy and direction = lv.mom - nlv1.mom
0159      *
0160      * Note: the tracking of the secondary ion is only needed when there is
0161      * a detail simulation of radiative decay for the recoiled nucleus.
0162      */
0163
0164     CELER_ENSURE(result.action == Interaction::Action::scattered);
0165
0166     return result;
0167 }
0168
0169 //---------------------------------------------------------------------------//
0170 }  // namespace celeritas