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