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