em/xs/EPlusGGMacroXsCalculator.hh

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