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