msc/detail/UrbanPositronCorrector.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2023-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/msc/detail/UrbanPositronCorrector.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/Assert.hh"
0011 #include "corecel/Macros.hh"
0012 #include "corecel/math/Algorithms.hh"
0013 #include "celeritas/Types.hh"
0014
0015 namespace celeritas
0016 {
0017 namespace detail
0018 {
0019 //---------------------------------------------------------------------------//
0020 /*!
0021  * Calculate a correction on theta0 for positrons.
0022  */
0023 class UrbanPositronCorrector
0024 {
0025   public:
0026     // Construct with effective Z
0027     explicit inline CELER_FUNCTION UrbanPositronCorrector(real_type zeff);
0028
0029     // Calculate correction with the ratio of avg energy to electron mass
0030     inline CELER_FUNCTION real_type operator()(real_type y) const;
0031
0032   private:
0033     real_type a_, b_, c_, d_, mult_;
0034 };
0035
0036 //---------------------------------------------------------------------------//
0037 // INLINE DEFINITIONS
0038 //---------------------------------------------------------------------------//
0039 /*!
0040  * Construct with effective Z.
0041  */
0042 CELER_FUNCTION UrbanPositronCorrector::UrbanPositronCorrector(real_type zeff)
0043 {
0044     CELER_EXPECT(zeff >= 1);
0045     using PolyLin = PolyEvaluator<real_type, 1>;
0046     using PolyQuad = PolyEvaluator<real_type, 2>;
0047
0048     a_ = PolyLin(0.994, -4.08e-3)(zeff);
0049     b_ = PolyQuad(7.16, 52.6, 365)(1 / zeff);
0050     c_ = PolyLin(1, -4.47e-3)(zeff);
0051     d_ = real_type(1.21e-3) * zeff;
0052     mult_ = PolyQuad(1.41125, -1.86427e-2, 1.84035e-4)(zeff);
0053 }
0054
0055 //---------------------------------------------------------------------------//
0056 /*!
0057  * Calculate the correction.
0058  *
0059  * The input parameter is a unitless ratio of the geometric mean energy to the
0060  * electron mass : \f[
0061    y = \frac{\sqrt{E_i E_f}}{m_e}
0062   \f]
0063  */
0064 CELER_FUNCTION real_type UrbanPositronCorrector::operator()(real_type y) const
0065 {
0066     CELER_EXPECT(y > 0);
0067     real_type corr = [this, y] {
0068         constexpr real_type xl = 0.6;
0069         constexpr real_type xh = 0.9;
0070         constexpr real_type e = 113;
0071
0072         // x = sqrt((y^2 + 2y) / (y^2 + 2y + 1))
0073         real_type x = std::sqrt(y * (y + 2) / ipow<2>(y + 1));
0074         if (x < xl)
0075         {
0076             return a_ * (1 - std::exp(-b_ * x));
0077         }
0078         if (x > xh)
0079         {
0080             return c_ + d_ * std::exp(e * (x - 1));
0081         }
0082         real_type yl = a_ * (1 - std::exp(-b_ * xl));
0083         real_type yh = c_ + d_ * std::exp(e * (xh - 1));
0084         real_type y0 = (yh - yl) / (xh - xl);
0085         real_type y1 = yl - y0 * xl;
0086         return y0 * x + y1;
0087     }();
0088
0089     return corr * mult_;
0090 }
0091
0092 //---------------------------------------------------------------------------//
0093 }  // namespace detail
0094 }  // namespace celeritas