interactor/detail/MomentumTransferSampler.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/MomentumTransferSampler.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <cmath>
0010
0011 #include "corecel/Macros.hh"
0012 #include "corecel/Types.hh"
0013 #include "corecel/math/Algorithms.hh"
0014 #include "corecel/random/distribution/BernoulliDistribution.hh"
0015 #include "corecel/random/distribution/UniformRealDistribution.hh"
0016 #include "celeritas/Quantities.hh"
0017 #include "celeritas/Types.hh"
0018 #include "celeritas/mat/IsotopeView.hh"
0019
0020 namespace celeritas
0021 {
0022 namespace detail
0023 {
0024 //---------------------------------------------------------------------------//
0025 /*!
0026  * Sample the momentum transfer (\f$ Q^{2} = -t \f$) of the neutron-nucleus
0027  * elastic scattering.
0028  *
0029  * \note This performs the same sampling routine as in Geant4's
0030  *  G4ChipsElasticModel and G4ChipsNeutronElasticXS where the differential
0031  *  cross section for quark-exchange and the amplitude of the scattering are
0032  *  parameterized in terms of the neutron momentum (GeV/c) in the lab frame.
0033  */
0034 class MomentumTransferSampler
0035 {
0036   public:
0037     //!@{
0038     //! \name Type aliases
0039     using Mass = units::MevMass;
0040     using Momentum = units::MevMomentum;
0041     using AtomicMassNumber = AtomicNumber;
0042     //!@}
0043
0044   public:
0045     // Construct with shared and target data, and the neutron momentum
0046     inline CELER_FUNCTION
0047     MomentumTransferSampler(NeutronElasticRef const& shared,
0048                             IsotopeView const& target,
0049                             Momentum neutron_p);
0050
0051     // Sample the momentum transfer
0052     template<class Engine>
0053     inline CELER_FUNCTION real_type operator()(Engine& rng);
0054
0055   private:
0056     //// TYPES ////
0057
0058     using UniformRealDist = UniformRealDistribution<real_type>;
0059
0060     //// DATA ////
0061
0062     ChipsDiffXsCoefficients::ChipsArray const& par_;
0063
0064     // Mass of neutron and target
0065     Mass neutron_mass_;
0066     Mass target_mass_;
0067
0068     // Atomic mass number (A) of the target
0069     AtomicMassNumber amass_;
0070     bool heavy_target_{false};
0071     // Momentum magnitude of the incident neutron [MeV/c]
0072     Momentum neutron_p_;
0073     // Maximum momentum transfer for the elastic scattering [GeV^2/c^2]
0074     real_type max_q_sq_;
0075     // Parameters for the CHIPS differential cross section
0076     ExchangeParameters par_q_sq_;
0077
0078     //// HELPER FUNCTIONS ////
0079
0080     // Sample the slope
0081     template<class Engine>
0082     inline CELER_FUNCTION real_type sample_q_sq(real_type radius, Engine& rng)
0083     {
0084         return -std::log(1 - radius * generate_canonical(rng));
0085     }
0086
0087     // Calculate the maximum momentum transfer
0088     inline CELER_FUNCTION real_type calc_max_q_sq(Momentum p) const;
0089
0090     // Calculate parameters used in the quark-exchange process
0091     inline CELER_FUNCTION ExchangeParameters calc_par_q_sq(Momentum p) const;
0092
0093     //// COMMON PROPERTIES ////
0094
0095     // Covert from clhep::MeV value to clhep::GeV value
0096     static CELER_CONSTEXPR_FUNCTION real_type to_gev() { return 1e-3; }
0097
0098     // S-wave limit for neutron, log(p) < -4.3 (GeV/c) (kinetic energy < 0.1
0099     // MeV)
0100     static CELER_CONSTEXPR_FUNCTION Momentum s_wave_limit()
0101     {
0102         return Momentum{13.568559012200934};  // = exp(-4.3) * 1000
0103     }
0104
0105     // Limit of the slope square
0106     static CELER_CONSTEXPR_FUNCTION real_type tolerance_slope_sq()
0107     {
0108         return 1e-7;
0109     }
0110 };
0111
0112 //---------------------------------------------------------------------------//
0113 // Inline DEFINITIONS
0114 //---------------------------------------------------------------------------//
0115 /*!
0116  * Construct from the incident momentum in the lab frame and the target(A).
0117  *
0118  * \note The incident neutron momentum, and neutron and nucleus masses are
0119  *  converted to the GeV value .
0120  */
0121 CELER_FUNCTION
0122 MomentumTransferSampler::MomentumTransferSampler(NeutronElasticRef const& shared,
0123                                                  IsotopeView const& target,
0124                                                  Momentum neutron_p)
0125     : par_(shared.coeffs[target.isotope_id()].par)
0126     , neutron_mass_(shared.neutron_mass)
0127     , target_mass_(target.nuclear_mass())
0128     , amass_(target.atomic_mass_number())
0129     , heavy_target_(amass_.get() > 6)
0130     , neutron_p_(neutron_p)
0131     , max_q_sq_(calc_max_q_sq(neutron_p_))
0132     , par_q_sq_(calc_par_q_sq(neutron_p_))
0133 {
0134 }
0135
0136 //---------------------------------------------------------------------------//
0137 /*!
0138  * Sample the momentum transfer of the neutron-nucleus elastic scattering
0139  * based on G4ChipsElasticModel and G4ChipsNeutronElasticXS of the Geant4
0140  * 11.2 release.
0141  *
0142  */
0143 template<class Engine>
0144 CELER_FUNCTION auto
0145 MomentumTransferSampler::operator()(Engine& rng) -> real_type
0146 {
0147     // Sample \f$ Q^{2} \f$ below S-wave limit
0148     if (neutron_p_ < this->s_wave_limit())
0149     {
0150         return max_q_sq_ * generate_canonical(rng) / ipow<2>(this->to_gev());
0151     }
0152
0153     // Sample \f$ Q^{2} \f$
0154     real_type q_sq = 0;
0155
0156     if (amass_ == AtomicMassNumber{1})
0157     {
0158         // Special case for the \f$ n + p \rightarrow n + p \f$ channel
0159         real_type const r[2] = {-std::expm1(-max_q_sq_ * par_q_sq_.slope[0]),
0160                                 -std::expm1(-max_q_sq_ * par_q_sq_.slope[1])};
0161
0162         real_type const mi[2]
0163             = {r[0] * par_q_sq_.expnt[0],
0164                r[1] * par_q_sq_.expnt[1] / par_q_sq_.slope[1]};
0165
0166         // Sample by t-channel and u-channel (charge exchange)
0167         q_sq = (BernoulliDistribution(mi[0], mi[1])(rng))
0168                    ? sample_q_sq(r[0], rng) / par_q_sq_.slope[0]
0169                    : max_q_sq_ - sample_q_sq(r[1], rng) / par_q_sq_.slope[1];
0170     }
0171     else
0172     {
0173         // Sample \f$ Q^{2} \f$ for \f$ n + A \rightarrow n + A \f$
0174         constexpr real_type one_third = 1 / real_type(3);
0175         constexpr real_type one_fifth{0.2};
0176         constexpr real_type one_seventh = 1 / real_type(7);
0177
0178         real_type const r[4]
0179             = {-std::expm1(-max_q_sq_
0180                            * (par_q_sq_.slope[0] + max_q_sq_ * par_q_sq_.ss)),
0181                -std::expm1(
0182                    -(heavy_target_ ? ipow<5>(max_q_sq_) : ipow<3>(max_q_sq_))
0183                    * par_q_sq_.slope[1]),
0184                -std::expm1(-(heavy_target_ ? ipow<7>(max_q_sq_) : max_q_sq_)
0185                            * par_q_sq_.slope[2]),
0186                -std::expm1(-max_q_sq_ * par_q_sq_.slope[3])};
0187
0188         real_type mi[6] = {};
0189         for (auto i : range(4))
0190         {
0191             mi[i] = r[i] * par_q_sq_.expnt[i];
0192         }
0193         mi[4] = mi[0] + mi[1];
0194         mi[5] = mi[2] + mi[4];
0195
0196         real_type rand = (mi[3] + mi[5]) * generate_canonical(rng);
0197         if (rand < mi[0])
0198         {
0199             real_type tss = 2 * par_q_sq_.ss;
0200             q_sq = this->sample_q_sq(r[0], rng) / par_q_sq_.slope[0];
0201             if (std::fabs(tss) > this->tolerance_slope_sq())
0202             {
0203                 q_sq = (std::sqrt(par_q_sq_.slope[0]
0204                                   * (par_q_sq_.slope[0] + 2 * tss * q_sq))
0205                         - par_q_sq_.slope[0])
0206                        / tss;
0207             }
0208         }
0209         else if (rand < mi[4])
0210         {
0211             q_sq = clamp_to_nonneg(this->sample_q_sq(r[1], rng)
0212                                    / par_q_sq_.slope[1]);
0213             q_sq = std::pow(q_sq, heavy_target_ ? one_fifth : one_third);
0214         }
0215         else if (rand < mi[5])
0216         {
0217             q_sq = clamp_to_nonneg(this->sample_q_sq(r[2], rng)
0218                                    / par_q_sq_.slope[2]);
0219             if (heavy_target_)
0220             {
0221                 q_sq = std::pow(q_sq, one_seventh);
0222             }
0223         }
0224         else
0225         {
0226             q_sq = this->sample_q_sq(r[3], rng) / par_q_sq_.slope[3];
0227             // u reduced for light-A (starts from 0)
0228             if (!heavy_target_)
0229             {
0230                 q_sq = max_q_sq_ - q_sq;
0231             }
0232         }
0233     }
0234     return clamp(q_sq, real_type{0}, max_q_sq_) / ipow<2>(this->to_gev());
0235 }
0236
0237 //---------------------------------------------------------------------------//
0238 /*!
0239  * Returns the maximum \f$ -t = Q^{2} \f$ (value in \f$ clhep::GeV^{2} \f$) of
0240  * the elastic scattering between the incident neutron of which momentum is
0241  * \c p (value in clhep::GeV) and the target nucleus (A = Z + N) where Z > 0.
0242  *
0243  * For the neutron-nucleus scattering, the maximum momentum transfer is
0244  * \f$ Q_{max}^{2} = 4 M_{A}^{2} * p^{2}/(mds) \f$ where the Mandelstam mds
0245  * is \f$ mds = 2 M_{A} E_{neutron} + M_{neutron}^{2} + M_{A}^{2} \f$ (value
0246  * in clhep::GeV square. The the neutron-neutron channel, the maximum momentum
0247  * is \f$ Q^{2} = E_{neutron} * M_{neutron} - M_{neutron}^{2} \f$ when the
0248  * collision angle is 90 degree in the center of mass system, is currently
0249  * excluded, but may be supported if there is a user case.
0250  */
0251 CELER_FUNCTION real_type MomentumTransferSampler::calc_max_q_sq(Momentum p) const
0252 {
0253     // Momentum and mass square of the incident neutron
0254     real_type target_mass = value_as<Mass>(target_mass_);
0255
0256     real_type p_sq = ipow<2>(value_as<units::MevMomentum>(p));
0257     real_type m_sq = ipow<2>(value_as<Mass>(neutron_mass_));
0258     real_type a_sq = ipow<2>(target_mass);
0259
0260     // Return the maximum momentum transfer in the value of clhep::GeV square
0261     return 4 * a_sq * p_sq * ipow<2>(this->to_gev())
0262            / (2 * target_mass * std::sqrt(p_sq + m_sq) + m_sq + a_sq);
0263 }
0264
0265 //---------------------------------------------------------------------------//
0266 /*!
0267  * Calculate parameters of the neutron-nucleus elastic scattering.
0268  *
0269  * \param neutron_p the neutron momentum in the lab frame (value in clhep::GeV
0270  * unit).
0271  */
0272 CELER_FUNCTION
0273 auto MomentumTransferSampler::calc_par_q_sq(Momentum neutron_p) const
0274     -> ExchangeParameters
0275 {
0276     // ExchangeParameters
0277     real_type p = value_as<units::MevMomentum>(neutron_p) * this->to_gev();
0278     real_type lp = std::log(p);
0279     real_type sp = std::sqrt(p);
0280     real_type p2 = ipow<2>(p);
0281     real_type p3 = p2 * p;
0282     real_type p4 = p3 * p;
0283
0284     ExchangeParameters result;
0285
0286     if (amass_ == AtomicMassNumber{1})
0287     {
0288         // Special case for the \f$ n + p \rightarrow n + p \f$ channel
0289         real_type dl1 = lp - par_[3];
0290
0291         constexpr real_type np_el[15] = {
0292             6.75,
0293             0.14,
0294             13,
0295             0.14,
0296             0.6,
0297             0.00013,
0298             75,
0299             0.001,
0300             7.2,
0301             4.32,
0302             0.012,
0303             2.5,
0304             12,
0305             0.34,
0306             18,
0307         };
0308
0309         result.expnt[0] = (np_el[0] + np_el[1] * ipow<2>(dl1) + np_el[2] / p)
0310                               / (1 + np_el[3] / p4)
0311                           + np_el[4] / (p4 + np_el[5]);
0312         result.slope[0] = (np_el[8] + np_el[9] / (ipow<2>(p4) + np_el[10] * p3))
0313                           / (1 + np_el[11] / p4);
0314         result.expnt[1] = (np_el[6] + np_el[7] / p4 / p) / p3;
0315         result.slope[1] = np_el[12] / (p * sp + np_el[13]);
0316         result.ss = np_el[14];
0317     }
0318     else
0319     {
0320         real_type p5 = p4 * p;
0321         real_type p6 = p5 * p;
0322         real_type p8 = p6 * p2;
0323         real_type p10 = p8 * p2;
0324         real_type p12 = p10 * p2;
0325         real_type p16 = ipow<2>(p8);
0326         real_type dl = lp - real_type(5);
0327
0328         if (!heavy_target_)
0329         {
0330             real_type pah = std::pow(p, real_type(0.5) * amass_.get());
0331             real_type pa = ipow<2>(pah);
0332             real_type pa2 = ipow<2>(pa);
0333             result.expnt[0] = par_[0] / (1 + par_[1] * p4 * pa)
0334                               + par_[2] / (p4 + par_[3] * p4 / pa2)
0335                               + (par_[4] * ipow<2>(dl) + par_[5])
0336                                     / (1 + par_[6] / p2);
0337             result.slope[0] = (par_[7] + par_[8] * p2) / (p4 + par_[9] / pah)
0338                               + par_[10];
0339             result.ss = par_[11] / (1 + par_[12] / p2)
0340                         + par_[13] / (p6 / pa + par_[14] / p16);
0341             result.expnt[1] = par_[15] / (pa / p2 + par_[16] / p4) + par_[17];
0342             result.slope[1] = par_[18] * std::pow(p, par_[19])
0343                               + par_[20] / (p8 + par_[21] / p16);
0344             result.expnt[2] = par_[22] / (pa * p + par_[23] / pa) + par_[24];
0345             result.slope[2] = par_[25] / (p3 + par_[26] / p6)
0346                               + par_[27] / (1 + par_[28] / p2);
0347             result.expnt[3]
0348                 = p2
0349                   * (pah * par_[29] * std::exp(-pah * par_[30])
0350                      + par_[31] / (1 + par_[32] * std::pow(p, par_[33])));
0351             result.slope[3] = par_[34] * pa / p2 / (1 + pa * par_[35]);
0352         }
0353         else
0354         {
0355             result.expnt[0] = par_[0] / (1 + par_[1] / p4)
0356                               + par_[2] / (p4 + par_[3] / p2)
0357                               + par_[4] / (p5 + par_[5] / p16);
0358             result.slope[0] = (par_[6] / p8 + par_[10])
0359                                   / (p + par_[7] / std::pow(p, par_[11]))
0360                               + par_[8] / (1 + par_[9] / p4);
0361             result.ss = par_[12] / (p4 / std::pow(p, par_[14]) + par_[13] / p4);
0362             result.expnt[1] = par_[15] / p4
0363                                   / (std::pow(p, par_[16]) + par_[17] / p12)
0364                               + par_[18];
0365             result.slope[1] = par_[19] / std::pow(p, par_[20])
0366                               + par_[21] / std::pow(p, par_[22]);
0367             result.expnt[2] = par_[23] / std::pow(p, par_[26])
0368                                   / (1 + par_[27] / p12)
0369                               + par_[24] / (1. + par_[25] / p6);
0370             result.slope[2] = par_[28] / p8 + par_[29] / p2
0371                               + par_[30] / (1 + par_[31] / p8);
0372             result.expnt[3]
0373                 = (par_[32] / p4 + par_[37] / p) / (1 + par_[33] / p10)
0374                   + (par_[34] + par_[35] * dl * dl) / (1 + par_[36] / p12);
0375             result.slope[3] = par_[38] / (1 + par_[39] / p)
0376                               + par_[40] * p4 / (1 + par_[41] * p5);
0377         }
0378     }
0379
0380     return result;
0381 }
0382
0383 //---------------------------------------------------------------------------//
0384 }  // namespace detail
0385 }  // namespace celeritas