em/xs/EPlusGGMacroXsCalculator.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/EPlusGGMacroXsCalculator.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <cmath>
0010
0011 #include "corecel/Assert.hh"
0012 #include "corecel/Macros.hh"
0013 #include "corecel/Types.hh"
0014 #include "corecel/cont/Span.hh"
0015 #include "corecel/math/Algorithms.hh"
0016 #include "corecel/math/PolyEvaluator.hh"
0017 #include "celeritas/Constants.hh"
0018 #include "celeritas/UnitTypes.hh"
0019 #include "celeritas/em/data/EPlusGGData.hh"
0020 #include "celeritas/mat/MaterialView.hh"
0021
0022 namespace celeritas
0023 {
0024 //---------------------------------------------------------------------------//
0025 /*!
0026  * Calculates the macroscopic cross section of positron annihilation.
0027  *
0028  * The Heitler formula (section 10.3.2 of the Geant4 Physics Reference Manual,
0029  * Release 10.6) is used to compute the macroscopic cross section for positron
0030  * annihilation on the fly at the given energy.
0031  */
0032 class EPlusGGMacroXsCalculator
0033 {
0034   public:
0035     //!@{
0036     //! \name Type aliases
0037     using Energy = units::MevEnergy;
0038     using XsUnits = units::Native;  // [1/len]
0039     //!@}
0040
0041   public:
0042     // Construct with shared data and material
0043     inline CELER_FUNCTION
0044     EPlusGGMacroXsCalculator(EPlusGGData const& shared,
0045                              MaterialView const& material);
0046
0047     // Compute cross section on the fly at the given energy
0048     inline CELER_FUNCTION real_type operator()(Energy energy) const;
0049
0050     //! Minimum energy for Heitler formula validity
0051     static CELER_CONSTEXPR_FUNCTION Energy min_energy()
0052     {
0053         return units::MevEnergy{1e-6};
0054     }
0055
0056   private:
0057     real_type const electron_mass_;
0058     real_type const electron_density_;
0059 };
0060
0061 //---------------------------------------------------------------------------//
0062 // INLINE DEFINITIONS
0063 //---------------------------------------------------------------------------//
0064 /*!
0065  * Construct with material.
0066  */
0067 CELER_FUNCTION
0068 EPlusGGMacroXsCalculator::EPlusGGMacroXsCalculator(EPlusGGData const& shared,
0069                                                    MaterialView const& material)
0070     : electron_mass_(value_as<units::MevMass>(shared.electron_mass))
0071     , electron_density_(material.electron_density())
0072 {
0073 }
0074
0075 //---------------------------------------------------------------------------//
0076 /*!
0077  * Compute macroscopic cross section in native units.
0078  */
0079 CELER_FUNCTION real_type EPlusGGMacroXsCalculator::operator()(Energy energy) const
0080 {
0081     using constants::pi;
0082     using constants::r_electron;
0083     using PolyQuad = PolyEvaluator<real_type, 2>;
0084
0085     real_type const gamma
0086         = celeritas::max(energy.value(), value_as<Energy>(this->min_energy()))
0087           / electron_mass_;
0088     real_type const sqrt_gg2 = std::sqrt(gamma * (gamma + 2));
0089
0090     real_type result
0091         = pi * ipow<2>(r_electron) * electron_density_
0092           * (PolyQuad{1, 4, 1}(gamma + 1) * std::log(gamma + 1 + sqrt_gg2)
0093              - (gamma + 4) * sqrt_gg2)
0094           / (gamma * ipow<2>(gamma + 2));
0095
0096     CELER_ENSURE(result >= 0);
0097     return result;
0098 }
0099
0100 //---------------------------------------------------------------------------//
0101 }  // namespace celeritas