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
0046     RBEnergySampler(NativeCRef<RelativisticBremData> 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 calculator
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(NativeCRef<RelativisticBremData> 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     , tmin_sq_(ipow<2>(value_as<Energy>(
0081           min(cutoffs.energy(shared.ids.gamma), particle.energy()))))
0082     , tmax_sq_(ipow<2>(value_as<Energy>(
0083           min(detail::high_energy_limit(), particle.energy()))))
0084 {
0085     CELER_ENSURE(tmax_sq_ >= tmin_sq_);
0086 }
0087
0088 //---------------------------------------------------------------------------//
0089 /*!
0090  * Sample the bremsstrahlung photon energy.
0091  */
0092 template<class Engine>
0093 CELER_FUNCTION auto RBEnergySampler::operator()(Engine& rng) -> Energy
0094 {
0095     real_type density_corr = calc_dxsec_.density_correction();
0096     ReciprocalDistribution<real_type> sample_exit_esq(tmin_sq_ + density_corr,
0097                                                       tmax_sq_ + density_corr);
0098
0099     // Sampled energy and corresponding cross section for rejection
0100     real_type gamma_energy{0};
0101     real_type dsigma{0};
0102
0103     do
0104     {
0105         gamma_energy = std::sqrt(sample_exit_esq(rng) - density_corr);
0106         dsigma = calc_dxsec_(Energy{gamma_energy});
0107     } while (RejectionSampler(dsigma, calc_dxsec_.maximum_value())(rng));
0108
0109     return Energy{gamma_energy};
0110 }
0111
0112 //---------------------------------------------------------------------------//
0113 }  // namespace detail
0114 }  // namespace celeritas