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