em/xs/WentzelTransportXsCalculator.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/em/xs/WentzelTransportXsCalculator.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/Assert.hh"
0010 #include "corecel/Macros.hh"
0011 #include "corecel/Types.hh"
0012 #include "corecel/cont/Span.hh"
0013 #include "celeritas/mat/MaterialView.hh"
0014
0015 #include "WentzelHelper.hh"
0016
0017 namespace celeritas
0018 {
0019 //---------------------------------------------------------------------------//
0020 /*!
0021  * Calculate the transport cross section for the Wentzel OK and VI model.
0022  *
0023  * \note This performs the same calculation as the Geant4 method
0024  * G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom.
0025  */
0026 class WentzelTransportXsCalculator
0027 {
0028   public:
0029     //!@{
0030     //! \name Type aliases
0031     using XsUnits = units::Native;  // [len^2]
0032     using Mass = units::MevMass;
0033     using MomentumSq = units::MevMomentumSq;
0034     //!@}
0035
0036   public:
0037     // Construct with particle and precalculatad Wentzel data
0038     inline CELER_FUNCTION
0039     WentzelTransportXsCalculator(ParticleTrackView const& particle,
0040                                  WentzelHelper const& helper);
0041
0042     // Calculate the transport cross section for the given angle [len^2]
0043     inline CELER_FUNCTION real_type operator()(real_type cos_thetamax) const;
0044
0045   private:
0046     //// DATA ////
0047
0048     AtomicNumber z_;
0049     real_type screening_coeff_;
0050     real_type cos_thetamax_elec_;
0051     real_type beta_sq_;
0052     real_type kin_factor_;
0053
0054     //// HELPER FUNCTIONS ////
0055
0056     // Calculate xs contribution from scattering off electrons or nucleus
0057     real_type calc_xs_contribution(real_type cos_thetamax) const;
0058
0059     //! Limit on (1 - \c cos_thetamax) / \c screening_coeff
0060     static CELER_CONSTEXPR_FUNCTION real_type limit() { return 0.1; }
0061 };
0062
0063 //---------------------------------------------------------------------------//
0064 // INLINE DEFINITIONS
0065 //---------------------------------------------------------------------------//
0066 /*!
0067  * Construct with particle and precalculatad Wentzel data.
0068  *
0069  * \c beta_sq should be calculated from the incident particle energy and mass.
0070  * \c screening_coeff and \c cos_thetamax_elec are calculated using the Wentzel
0071  * OK and VI model in \c WentzelHelper and depend on properties of the incident
0072  * particle, the energy cutoff in the current material, and the target element.
0073  */
0074 CELER_FUNCTION
0075 WentzelTransportXsCalculator::WentzelTransportXsCalculator(
0076     ParticleTrackView const& particle, WentzelHelper const& helper)
0077     : z_(helper.atomic_number())
0078     , screening_coeff_(2 * helper.screening_coefficient())
0079     , cos_thetamax_elec_(helper.cos_thetamax_electron())
0080     , beta_sq_(particle.beta_sq())
0081     , kin_factor_(helper.kin_factor())
0082 {
0083 }
0084
0085 //---------------------------------------------------------------------------//
0086 /*!
0087  * Calculate the transport cross section for the given angle [len^2].
0088  */
0089 CELER_FUNCTION real_type
0090 WentzelTransportXsCalculator::operator()(real_type cos_thetamax) const
0091 {
0092     CELER_EXPECT(cos_thetamax <= 1);
0093
0094     // Sum xs contributions from scattering off electrons and nucleus
0095     real_type xs_nuc = this->calc_xs_contribution(cos_thetamax);
0096     real_type xs_elec = cos_thetamax_elec_ > cos_thetamax
0097                             ? this->calc_xs_contribution(cos_thetamax_elec_)
0098                             : xs_nuc;
0099     real_type result = kin_factor_ * (xs_elec + z_.get() * xs_nuc);
0100
0101     CELER_ENSURE(result >= 0);
0102     return result;
0103 }
0104
0105 //---------------------------------------------------------------------------//
0106 /*!
0107  * Calculate contribution to xs from scattering off electrons or nucleus.
0108  */
0109 CELER_FUNCTION real_type WentzelTransportXsCalculator::calc_xs_contribution(
0110     real_type cos_thetamax) const
0111 {
0112     real_type result;
0113     real_type const spin = real_type(0.5);
0114     real_type x = (1 - cos_thetamax) / screening_coeff_;
0115     if (x < WentzelTransportXsCalculator::limit())
0116     {
0117         real_type x_sq = ipow<2>(x);
0118         result = real_type(0.5) * x_sq
0119                  * ((1 - real_type(4) / 3 * x + real_type(1.5) * x_sq)
0120                     - screening_coeff_ * spin * beta_sq_ * x
0121                           * (real_type(2) / 3 - x));
0122     }
0123     else
0124     {
0125         real_type x_1 = x / (1 + x);
0126         real_type log_x = std::log(1 + x);
0127         result = log_x - x_1
0128                  - screening_coeff_ * spin * beta_sq_ * (x + x_1 - 2 * log_x);
0129     }
0130     return clamp_to_nonneg(result);
0131 }
0132
0133 //---------------------------------------------------------------------------//
0134 }  // namespace celeritas