em/interactor/RayleighInteractor.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/RayleighInteractor.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/Macros.hh"
0011 #include "corecel/Types.hh"
0012 #include "corecel/data/Collection.hh"
0013 #include "corecel/math/Algorithms.hh"
0014 #include "corecel/math/ArrayUtils.hh"
0015 #include "celeritas/Quantities.hh"
0016 #include "celeritas/Types.hh"
0017 #include "celeritas/em/data/RayleighData.hh"
0018 #include "celeritas/phys/Interaction.hh"
0019 #include "celeritas/phys/InteractionUtils.hh"
0020 #include "celeritas/phys/ParticleTrackView.hh"
0021 #include "celeritas/random/Selector.hh"
0022 #include "celeritas/random/distribution/GenerateCanonical.hh"
0023 #include "celeritas/random/distribution/IsotropicDistribution.hh"
0024
0025 namespace celeritas
0026 {
0027 //---------------------------------------------------------------------------//
0028 /*!
0029  * Apply the Livermore model of Rayleigh scattering to photons.
0030  *
0031  * \note This performs the same sampling routine as in Geant4's
0032  * G4LivermoreRayleighModel class, as documented in section 6.2.2 of the
0033  * Geant4 Physics Reference (release 10.6).
0034  */
0035 class RayleighInteractor
0036 {
0037   public:
0038     // Construct with shared and state data
0039     inline CELER_FUNCTION RayleighInteractor(RayleighRef const& shared,
0040                                              ParticleTrackView const& particle,
0041                                              Real3 const& inc_direction,
0042                                              ElementId element_id);
0043
0044     // Sample an interaction with the given RNG
0045     template<class Engine>
0046     inline CELER_FUNCTION Interaction operator()(Engine& rng);
0047
0048   private:
0049     //// DATA ////
0050
0051     // Shared constant physics properties
0052     RayleighRef const& shared_;
0053     // Incident gamma energy
0054     units::MevEnergy const inc_energy_;
0055     // Incident direction
0056     Real3 const& inc_direction_;
0057     // Id of element
0058     ElementId element_id_;
0059
0060     //// CONSTANTS ////
0061
0062     //! A point where the functional form of the form factor fit changes
0063     static CELER_CONSTEXPR_FUNCTION real_type fit_slice() { return 0.02; }
0064
0065     //// HELPER TYPES ////
0066
0067     //! Intermediate data for sampling input
0068     struct SampleInput
0069     {
0070         real_type factor{0};
0071         Real3 weight{0, 0, 0};
0072         Real3 prob{0, 0, 0};
0073     };
0074
0075     //// HELPER FUNCTIONS ////
0076
0077     //! Evaluate weights and probabilities for the angular sampling algorithm
0078     inline CELER_FUNCTION auto evaluate_weight_and_prob() const -> SampleInput;
0079 };
0080
0081 //---------------------------------------------------------------------------//
0082 // INLINE DEFINITIONS
0083 //---------------------------------------------------------------------------//
0084 /*!
0085  * Construct with shared and state data.
0086  */
0087 CELER_FUNCTION
0088 RayleighInteractor::RayleighInteractor(RayleighRef const& shared,
0089                                        ParticleTrackView const& particle,
0090                                        Real3 const& direction,
0091                                        ElementId el_id)
0092
0093     : shared_(shared)
0094     , inc_energy_(particle.energy())
0095     , inc_direction_(direction)
0096     , element_id_(el_id)
0097 {
0098     CELER_EXPECT(particle.particle_id() == shared_.gamma);
0099     CELER_EXPECT(element_id_ < shared_.params.size());
0100 }
0101
0102 //---------------------------------------------------------------------------//
0103 /*!
0104  * Sample the Rayleigh scattering angle using the G4LivermoreRayleighModel
0105  * and G4RayleighAngularGenerator of Geant4 6.10.
0106  */
0107 template<class Engine>
0108 CELER_FUNCTION Interaction RayleighInteractor::operator()(Engine& rng)
0109 {
0110     // Construct interaction for change to primary (incident) particle
0111     Interaction result;
0112     result.energy = inc_energy_;
0113
0114     SampleInput input = this->evaluate_weight_and_prob();
0115
0116     Real3 const& pb = shared_.params[element_id_].b;
0117     Real3 const& pn = shared_.params[element_id_].n;
0118
0119     constexpr real_type half = 0.5;
0120     real_type cost;
0121
0122     do
0123     {
0124         // Sample index from input.prob
0125         auto const index = celeritas::make_selector(
0126             [&input](size_type i) { return input.prob[i]; },
0127             input.prob.size())(rng);
0128
0129         real_type const w = input.weight[index];
0130         real_type const ninv = 1 / pn[index];
0131         real_type const b = pb[index];
0132
0133         // Sampling of scattering angle
0134         real_type x;
0135         real_type y = w * generate_canonical(rng);
0136
0137         if (y < fit_slice())
0138         {
0139             x = y * ninv
0140                 * (1 + half * (ninv + 1) * y * (1 - (ninv + 2) * y / 3));
0141         }
0142         else
0143         {
0144             x = fastpow(1 - y, -ninv) - 1;
0145         }
0146
0147         cost = 1 - 2 * x / (b * input.factor);
0148
0149     } while (2 * generate_canonical(rng) > 1 + ipow<2>(cost) || cost < -1);
0150
0151     // Scattered direction
0152     result.direction = ExitingDirectionSampler{cost, inc_direction_}(rng);
0153
0154     CELER_ENSURE(result.action == Interaction::Action::scattered);
0155     return result;
0156 }
0157
0158 CELER_FUNCTION
0159 auto RayleighInteractor::evaluate_weight_and_prob() const -> SampleInput
0160 {
0161     Real3 const& a = shared_.params[element_id_].a;
0162     Real3 const& b = shared_.params[element_id_].b;
0163     Real3 const& n = shared_.params[element_id_].n;
0164
0165     SampleInput input;
0166     input.factor = ipow<2>(units::centimeter
0167                            / (constants::c_light * constants::h_planck)
0168                            * native_value_from(inc_energy_));
0169
0170     Real3 x = b;
0171     axpy(input.factor, b, &x);
0172
0173     Real3 prob;
0174     for (int i = 0; i < 3; ++i)
0175     {
0176         input.weight[i] = (x[i] > this->fit_slice())
0177                               ? 1 - fastpow(1 + x[i], -n[i])
0178                               : n[i] * x[i]
0179                                     * (1
0180                                        - (n[i] - 1) / 2 * x[i]
0181                                              * (1 - (n[i] - 2) / 3 * x[i]));
0182
0183         prob[i] = input.weight[i] * a[i] / (b[i] * n[i]);
0184     }
0185
0186     real_type inv_sum = 1 / (prob[0] + prob[1] + prob[2]);
0187     axpy(inv_sum, prob, &input.prob);
0188
0189     return input;
0190 }
0191
0192 //---------------------------------------------------------------------------//
0193 }  // namespace celeritas