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