interactor/detail/RBEnergySampler.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/RBEnergySampler.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/RelativisticBremData.hh"
0016 #include "celeritas/em/xs/RBDiffXsCalculator.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 #include "celeritas/random/distribution/ReciprocalDistribution.hh"
0022 #include "celeritas/random/distribution/RejectionSampler.hh"
0023
0024 #include "PhysicsConstants.hh"
0025
0026 namespace celeritas
0027 {
0028 namespace detail
0029 {
0030 //---------------------------------------------------------------------------//
0031 /*!
0032  * Sample the bremsstrahlung photon energy from the relativistic model.
0033  *
0034  * Based on \c G4eBremsstrahlungRelModel of the Geant4 10.7 release.
0035  */
0036 class RBEnergySampler
0037 {
0038   public:
0039     //!@{
0040     //! \name Type aliases
0041     using Energy = units::MevEnergy;
0042     //!@}
0043
0044   public:
0045     // Construct with shared and state data
0046     inline CELER_FUNCTION RBEnergySampler(RelativisticBremRef const& shared,
0047                                           ParticleTrackView const& particle,
0048                                           CutoffView const& cutoffs,
0049                                           MaterialView const& material,
0050                                           ElementComponentId elcomp_id);
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
0059     // Differential cross section calcuator
0060     RBDiffXsCalculator calc_dxsec_;
0061     // Square of minimum of incident particle energy and cutoff
0062     real_type tmin_sq_;
0063     // Square of production cutoff for gammas
0064     real_type tmax_sq_;
0065 };
0066
0067 //---------------------------------------------------------------------------//
0068 // INLINE DEFINITIONS
0069 //---------------------------------------------------------------------------//
0070 /*!
0071  * Construct from incident particle and energy.
0072  */
0073 CELER_FUNCTION
0074 RBEnergySampler::RBEnergySampler(RelativisticBremRef const& shared,
0075                                  ParticleTrackView const& particle,
0076                                  CutoffView const& cutoffs,
0077                                  MaterialView const& material,
0078                                  ElementComponentId elcomp_id)
0079     : calc_dxsec_(shared, particle, material, elcomp_id)
0080 {
0081     // Min and max kinetic energy limits for sampling the secondary photon
0082     tmin_sq_ = ipow<2>(value_as<Energy>(
0083         min(cutoffs.energy(shared.ids.gamma), particle.energy())));
0084     tmax_sq_ = ipow<2>(
0085         value_as<Energy>(min(detail::high_energy_limit(), particle.energy())));
0086
0087     CELER_ENSURE(tmax_sq_ >= tmin_sq_);
0088 }
0089
0090 //---------------------------------------------------------------------------//
0091 /*!
0092  * Sample the bremsstrahlung photon energy.
0093  */
0094 template<class Engine>
0095 CELER_FUNCTION auto RBEnergySampler::operator()(Engine& rng) -> Energy
0096 {
0097     real_type density_corr = calc_dxsec_.density_correction();
0098     ReciprocalDistribution<real_type> sample_exit_esq(tmin_sq_ + density_corr,
0099                                                       tmax_sq_ + density_corr);
0100
0101     // Sampled energy and corresponding cross section for rejection
0102     real_type gamma_energy{0};
0103     real_type dsigma{0};
0104
0105     do
0106     {
0107         gamma_energy = std::sqrt(sample_exit_esq(rng) - density_corr);
0108         dsigma = calc_dxsec_(Energy{gamma_energy});
0109     } while (RejectionSampler(dsigma, calc_dxsec_.maximum_value())(rng));
0110
0111     return Energy{gamma_energy};
0112 }
0113
0114 //---------------------------------------------------------------------------//
0115 }  // namespace detail
0116 }  // namespace celeritas