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