em/xs/WentzelMacroXsCalculator.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/WentzelMacroXsCalculator.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/Assert.hh"
0011 #include "corecel/Macros.hh"
0012 #include "corecel/Types.hh"
0013 #include "celeritas/em/data/WentzelOKVIData.hh"
0014 #include "celeritas/em/data/WentzelVIMscData.hh"
0015 #include "celeritas/mat/MaterialView.hh"
0016 #include "celeritas/phys/ParticleTrackView.hh"
0017
0018 #include "WentzelHelper.hh"
0019
0020 namespace celeritas
0021 {
0022 //---------------------------------------------------------------------------//
0023 /*!
0024  * Calculate the total cross section for the Wentzel VI MSC model.
0025  *
0026  * \note This performs the same calculation of the total cross section (\c
0027  * xtsec) as the Geant4 method
0028  * G4WentzelVIModel::ComputeTransportXSectionPerVolume.
0029  */
0030 class WentzelMacroXsCalculator
0031 {
0032   public:
0033     //!@{
0034     //! \name Type aliases
0035     using Energy = units::MevEnergy;
0036     using XsUnits = units::Native;  // [1/len]
0037     //!@}
0038
0039   public:
0040     // Construct with particle, material, and precalculatad Wentzel data
0041     inline CELER_FUNCTION
0042     WentzelMacroXsCalculator(ParticleTrackView const& particle,
0043                              MaterialView const& material,
0044                              NativeCRef<WentzelVIMscData> const& data,
0045                              NativeCRef<WentzelOKVIData> const& wentzel,
0046                              Energy cutoff);
0047
0048     // Compute the total cross section for the given angle
0049     inline CELER_FUNCTION real_type operator()(real_type cos_theta) const;
0050
0051   private:
0052     ParticleTrackView const& particle_;
0053     MaterialView const& material_;
0054     NativeCRef<WentzelOKVIData> const& wentzel_;
0055     CoulombIds const& ids_;
0056     Energy cutoff_;
0057 };
0058
0059 //---------------------------------------------------------------------------//
0060 // INLINE DEFINITIONS
0061 //---------------------------------------------------------------------------//
0062 /*!
0063  * Construct with shared model and material data.
0064  */
0065 CELER_FUNCTION
0066 WentzelMacroXsCalculator::WentzelMacroXsCalculator(
0067     ParticleTrackView const& particle,
0068     MaterialView const& material,
0069     NativeCRef<WentzelVIMscData> const& data,
0070     NativeCRef<WentzelOKVIData> const& wentzel,
0071     Energy cutoff)
0072     : particle_(particle)
0073     , material_(material)
0074     , wentzel_(wentzel)
0075     , ids_(data.ids)
0076     , cutoff_(cutoff)
0077 {
0078 }
0079
0080 //---------------------------------------------------------------------------//
0081 /*!
0082  * Compute the total cross section for the given angle.
0083  */
0084 CELER_FUNCTION real_type
0085 WentzelMacroXsCalculator::operator()(real_type cos_theta) const
0086 {
0087     real_type result = 0;
0088
0089     for (auto elcomp_id : range(ElementComponentId(material_.num_elements())))
0090     {
0091         AtomicNumber z = material_.make_element_view(elcomp_id).atomic_number();
0092         WentzelHelper helper(particle_, material_, z, wentzel_, ids_, cutoff_);
0093
0094         real_type cos_thetamax = helper.cos_thetamax_nuclear();
0095         if (cos_thetamax < cos_theta)
0096         {
0097             result += material_.elements()[elcomp_id.get()].fraction
0098                       * (helper.calc_xs_nuclear(cos_theta, cos_thetamax)
0099                          + helper.calc_xs_electron(cos_theta, cos_thetamax));
0100         }
0101     }
0102     result *= material_.number_density();
0103
0104     CELER_ENSURE(result >= 0);
0105     return result;
0106 }
0107
0108 //---------------------------------------------------------------------------//
0109 }  // namespace celeritas