interactor/detail/SBEnergySampler.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2021-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/interactor/detail/SBEnergySampler.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include <cmath>
0011
0012 #include "corecel/math/Algorithms.hh"
0013 #include "celeritas/Quantities.hh"
0014 #include "celeritas/Types.hh"
0015 #include "celeritas/em/data/SeltzerBergerData.hh"
0016 #include "celeritas/em/distribution/SBEnergyDistHelper.hh"
0017 #include "celeritas/em/distribution/SBEnergyDistribution.hh"
0018 #include "celeritas/mat/ElementView.hh"
0019 #include "celeritas/mat/MaterialView.hh"
0020 #include "celeritas/phys/CutoffView.hh"
0021 #include "celeritas/phys/ParticleTrackView.hh"
0022
0023 #include "PhysicsConstants.hh"
0024 #include "SBPositronXsCorrector.hh"
0025
0026 namespace celeritas
0027 {
0028 namespace detail
0029 {
0030 //---------------------------------------------------------------------------//
0031 /*!
0032  * Sample the bremsstrahlung photon energy from the SeltzerBerger model.
0033  */
0034 class SBEnergySampler
0035 {
0036   public:
0037     //!@{
0038     //! \name Type aliases
0039     using Energy = units::MevEnergy;
0040     using Mass = units::MevMass;
0041     using SBTable = NativeCRef<SeltzerBergerTableData>;
0042     //!@}
0043
0044   public:
0045     // Construct with shared and state data
0046     inline CELER_FUNCTION SBEnergySampler(SBTable const& differential_xs,
0047                                           ParticleTrackView const& particle,
0048                                           Energy gamma_cutoff,
0049                                           MaterialView const& material,
0050                                           ElementComponentId elcomp_id,
0051                                           bool is_electron);
0052
0053     // Sample the bremsstrahlung photon energy with the given RNG
0054     template<class Engine>
0055     inline CELER_FUNCTION Energy operator()(Engine& rng);
0056
0057   private:
0058     //// DATA ////
0059     // Differential cross section table
0060     SBTable const& differential_xs_;
0061     // Incident particle energy
0062     Energy inc_energy_;
0063     // Production cutoff for gammas
0064     Energy gamma_cutoff_;
0065     // Material in which interaction occurs
0066     MaterialView const& material_;
0067     // Element in which interaction occurs
0068     ElementComponentId elcomp_id_;
0069     // Incident particle mass
0070     Mass inc_mass_;
0071     // Incident particle identification flag
0072     bool is_electron_;
0073     // Density correction
0074     real_type density_correction_;
0075 };
0076
0077 //---------------------------------------------------------------------------//
0078 // INLINE DEFINITIONS
0079 //---------------------------------------------------------------------------//
0080 /*!
0081  * Construct from incident particle and energy.
0082  */
0083 CELER_FUNCTION
0084 SBEnergySampler::SBEnergySampler(SBTable const& differential_xs,
0085                                  ParticleTrackView const& particle,
0086                                  Energy gamma_cutoff,
0087                                  MaterialView const& material,
0088                                  ElementComponentId elcomp_id,
0089                                  bool is_electron)
0090     : differential_xs_(differential_xs)
0091     , inc_energy_(value_as<Energy>(particle.energy()))
0092     , gamma_cutoff_(gamma_cutoff)
0093     , material_(material)
0094     , elcomp_id_(elcomp_id)
0095     , inc_mass_(value_as<Mass>(particle.mass()))
0096     , is_electron_(is_electron)
0097 {
0098     // Density correction
0099     real_type density_factor = material.electron_density() * migdal_constant();
0100     density_correction_ = density_factor
0101                           * ipow<2>(value_as<Energy>(particle.total_energy()));
0102 }
0103
0104 //---------------------------------------------------------------------------//
0105 /*!
0106  * Sample the exiting energy by doing a table lookup and rejection.
0107  */
0108 template<class Engine>
0109 CELER_FUNCTION auto SBEnergySampler::operator()(Engine& rng) -> Energy
0110 {
0111     // Outgoing photon secondary energy sampler
0112     Energy gamma_exit_energy;
0113
0114     // Helper class preprocesses cross section bounds and calculates
0115     // distribution
0116     SBEnergyDistHelper sb_helper(
0117         differential_xs_,
0118         inc_energy_,
0119         material_.element_id(elcomp_id_),
0120         SBEnergyDistHelper::EnergySq{density_correction_},
0121         gamma_cutoff_);
0122
0123     if (is_electron_)
0124     {
0125         // Rejection sample without modifying cross section
0126         SBEnergyDistribution<SBElectronXsCorrector> sample_gamma_energy(
0127             sb_helper, {});
0128         gamma_exit_energy = sample_gamma_energy(rng);
0129     }
0130     else
0131     {
0132         SBEnergyDistribution<SBPositronXsCorrector> sample_gamma_energy(
0133             sb_helper,
0134             {inc_mass_,
0135              material_.make_element_view(elcomp_id_),
0136              gamma_cutoff_,
0137              inc_energy_});
0138         gamma_exit_energy = sample_gamma_energy(rng);
0139     }
0140
0141     return gamma_exit_energy;
0142 }
0143
0144 //---------------------------------------------------------------------------//
0145 }  // namespace detail
0146 }  // namespace celeritas