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
0039     ChipsNeutronElasticInteractor(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()(Engine& rng)
0116 {
0117     // Scattered neutron with respect to the axis of incident direction
0118     Interaction result;
0119
0120     // The momentum magnitude in the c.m. frame
0121     real_type target_mass = value_as<Mass>(target_.nuclear_mass());
0122
0123     real_type cm_p = value_as<Momentum>(neutron_p_)
0124                      / std::sqrt(1 + ipow<2>(neutron_mass_ / target_mass)
0125                                  + 2 * neutron_energy_ / target_mass);
0126
0127     // Sample the scattered direction from the invariant momentum transfer
0128     // squared (\f$ -t = Q^{2} \f$) in the c.m. frame
0129     real_type cos_theta
0130         = 1 - real_type(0.5) * sample_momentum_square_(rng) / ipow<2>(cm_p);
0131     CELER_ASSERT(std::fabs(cos_theta) <= 1);
0132
0133     // Boost to the center of mass (c.m.) frame
0134     auto nlv1 = FourVector::from_mass_momentum(
0135         Mass{neutron_mass_},
0136         Momentum{cm_p},
0137         from_spherical(cos_theta, sample_phi_(rng)));
0138
0139     FourVector lv({{0, 0, value_as<Momentum>(neutron_p_)},
0140                    neutron_energy_ + target_mass});
0141     boost(boost_vector(lv), &nlv1);
0142
0143     result.direction = rotate(make_unit_vector(nlv1.mom), inc_direction_);
0144
0145     // Kinetic energy of the scattered neutron and the recoiled nucleus
0146     lv.energy -= nlv1.energy;
0147     result.energy = Energy(nlv1.energy - neutron_mass_);
0148     real_type recoil_energy = clamp_to_nonneg(lv.energy - target_mass);
0149     result.energy_deposition = Energy{recoil_energy};
0150
0151     /*!
0152      * \todo Create a secondary ion(Z, N) if recoil_energy > recoil_threshold
0153      * with energy = recoil_energy and direction = lv.mom - nlv1.mom
0154      *
0155      * Note: the tracking of the secondary ion is only needed when there is
0156      * a detail simulation of radiative decay for the recoiled nucleus.
0157      */
0158
0159     CELER_ENSURE(result.action == Interaction::Action::scattered);
0160
0161     return result;
0162 }
0163
0164 //---------------------------------------------------------------------------//
0165 }  // namespace celeritas