em/distribution/MuPPEnergyDistribution.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/distribution/MuPPEnergyDistribution.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/Constants.hh"
0010 #include "corecel/Macros.hh"
0011 #include "corecel/Types.hh"
0012 #include "corecel/grid/FindInterp.hh"
0013 #include "corecel/grid/Interpolator.hh"
0014 #include "corecel/grid/NonuniformGrid.hh"
0015 #include "corecel/grid/TwodGridCalculator.hh"
0016 #include "corecel/math/Algorithms.hh"
0017 #include "corecel/random/distribution/UniformRealDistribution.hh"
0018 #include "celeritas/Quantities.hh"
0019 #include "celeritas/em/data/MuPairProductionData.hh"
0020 #include "celeritas/grid/InverseCdfFinder.hh"
0021 #include "celeritas/mat/ElementView.hh"
0022 #include "celeritas/phys/CutoffView.hh"
0023 #include "celeritas/phys/ParticleTrackView.hh"
0024
0025 namespace celeritas
0026 {
0027 //---------------------------------------------------------------------------//
0028 /*!
0029  * Sample the electron and positron energies for muon pair production.
0030  *
0031  * The energy transfer to the electron-positron pair is sampled using inverse
0032  * transform sampling on a tabulated CDF. The CDF is calculated on a 2D grid,
0033  * where the x-axis is the log of the incident muon energy and the y-axis is
0034  * the log of the ratio of the energy transfer to the incident particle energy.
0035  * Because the shape of the distribution depends only weakly on the atomic
0036  * number, the CDF is calculated for a hardcoded set of points equally spaced
0037  * in \f$ \log Z \f$ and linearly interpolated.
0038  *
0039  * The formula used for the differential cross section is valid when the
0040  * maximum energy transfer to the electron-positron pair lies between \f$
0041  * \epsilon_{\text{min}} = 4 m \f$, where \f$ m \f$ is the electron mass, and
0042  * \f[
0043    \epsilon_{\text{max}} = E + \frac{3 \sqrt{e}}{4} \mu Z^{1/3}),
0044  * \f]
0045  * where \f$ E = T + \mu \f$ is the total muon energy, \f$ \mu \f$ is the muon
0046  * mass, \f$ e \f$ is Euler's number, and \f$ Z \f$ is the atomic number.
0047  *
0048  * The maximum energy partition between the electron and positron is calculated
0049  * as
0050  * \f[
0051    r_{\text{max}} = \left[1 - 6 \frac{\mu^2}{E (E - \epsilon)} \right] \sqrt{1
0052    - \epsilon_{\text{min}} / \epsilon}.
0053  * \f]
0054  * The partition \f$ r \f$ is then sampled uniformly in \f$ [-r_{\text{max}},
0055  * r_{\text{max}}) \f$.
0056  */
0057 class MuPPEnergyDistribution
0058 {
0059   public:
0060     //!@{
0061     //! \name Type aliases
0062     using Mass = units::MevMass;
0063     using Energy = units::MevEnergy;
0064     //!@}
0065
0066     //! Sampled secondary energies
0067     struct PairEnergy
0068     {
0069         Energy electron;
0070         Energy positron;
0071     };
0072
0073   public:
0074     // Construct from shared and incident particle data
0075     inline CELER_FUNCTION
0076     MuPPEnergyDistribution(NativeCRef<MuPairProductionData> const& shared,
0077                            ParticleTrackView const& particle,
0078                            CutoffView const& cutoffs,
0079                            ElementView const& element);
0080
0081     template<class Engine>
0082     inline CELER_FUNCTION PairEnergy operator()(Engine& rng);
0083
0084     //! Minimum energy of the electron-positron pair [MeV].
0085     CELER_FUNCTION Energy min_pair_energy() const
0086     {
0087         return Energy(min_pair_energy_);
0088     }
0089
0090     //! Maximum energy of the electron-positron pair [MeV].
0091     CELER_FUNCTION Energy max_pair_energy() const
0092     {
0093         return Energy(max_pair_energy_);
0094     }
0095
0096   private:
0097     //// DATA ////
0098
0099     // CDF table for sampling the pair energy
0100     NativeCRef<MuPairProductionTableData> const& table_;
0101     // Incident particle energy [MeV]
0102     real_type inc_energy_;
0103     // Incident particle total energy [MeV]
0104     real_type total_energy_;
0105     // Square of the muon mass
0106     real_type inc_mass_sq_;
0107     // Secondary mass
0108     real_type electron_mass_;
0109     // Minimum energy transfer to electron/positron pair [MeV]
0110     real_type min_pair_energy_;
0111     // Maximum energy transfer to electron/positron pair [MeV]
0112     real_type max_pair_energy_;
0113     // Minimum incident particle kinetic energy [MeV]
0114     real_type min_energy_;
0115     // Log Z grid interpolation for the target element
0116     FindInterp<real_type> logz_interp_;
0117     // Coefficient for calculating the pair energy
0118     real_type coeff_;
0119     // Lower bound on the ratio of the pair energy to the incident energy
0120     real_type y_min_;
0121     // Upper bound on the ratio of the pair energy to the incident energy
0122     real_type y_max_;
0123
0124     //// HELPER FUNCTIONS ////
0125
0126     // Sample the scaled energy and interpolate in log Z
0127     template<class Engine>
0128     inline CELER_FUNCTION real_type sample_scaled_energy(Engine& rng) const;
0129
0130     // Calculate the scaled energy for a given Z grid and sampled CDF
0131     inline CELER_FUNCTION real_type calc_scaled_energy(size_type z_idx,
0132                                                        real_type u) const;
0133 };
0134
0135 //---------------------------------------------------------------------------//
0136 // INLINE DEFINITIONS
0137 //---------------------------------------------------------------------------//
0138 /*!
0139  * Construct from shared and incident particle data.
0140  *
0141  * The incident energy *must* be within the bounds of the sampling table data.
0142  */
0143 CELER_FUNCTION
0144 MuPPEnergyDistribution::MuPPEnergyDistribution(
0145     NativeCRef<MuPairProductionData> const& shared,
0146     ParticleTrackView const& particle,
0147     CutoffView const& cutoffs,
0148     ElementView const& element)
0149     : table_(shared.table)
0150     , inc_energy_(value_as<Energy>(particle.energy()))
0151     , total_energy_(value_as<Energy>(particle.total_energy()))
0152     , inc_mass_sq_(ipow<2>(value_as<Mass>(particle.mass())))
0153     , electron_mass_(value_as<Mass>(shared.electron_mass))
0154     , min_pair_energy_(4 * value_as<Mass>(shared.electron_mass))
0155     , max_pair_energy_(inc_energy_
0156                        + value_as<Mass>(particle.mass())
0157                              * (1
0158                                 - real_type(0.75) * constants::sqrt_euler
0159                                       * element.cbrt_z()))
0160     , min_energy_(max(value_as<Energy>(cutoffs.energy(shared.ids.positron)),
0161                       min_pair_energy_))
0162 {
0163     CELER_EXPECT(max_pair_energy_ > min_energy_);
0164
0165     NonuniformGrid logz_grid(table_.logz_grid, table_.reals);
0166     logz_interp_ = find_interp(logz_grid, element.log_z());
0167
0168     NonuniformGrid y_grid(
0169         table_.grids[ItemId<TwodGridData>(logz_interp_.index)].y, table_.reals);
0170     coeff_ = std::log(min_pair_energy_ / inc_energy_) / y_grid.front();
0171
0172     // Compute the bounds on the ratio of the pair energy to incident energy
0173     y_min_ = std::log(min_energy_ / inc_energy_) / coeff_;
0174     y_max_ = std::log(max_pair_energy_ / inc_energy_) / coeff_;
0175
0176     // Check that the bounds are within the grid bounds
0177     CELER_ASSERT(y_min_ >= y_grid.front()
0178                  || soft_equal(y_grid.front(), y_min_));
0179     CELER_ASSERT(y_max_ <= y_grid.back() || soft_equal(y_grid.back(), y_max_));
0180     y_min_ = max(y_min_, y_grid.front());
0181     y_max_ = min(y_max_, y_grid.back());
0182 }
0183
0184 //---------------------------------------------------------------------------//
0185 /*!
0186  * Sample the exiting pair energy.
0187  */
0188 template<class Engine>
0189 CELER_FUNCTION auto MuPPEnergyDistribution::operator()(Engine& rng)
0190     -> PairEnergy
0191 {
0192     // Sample the energy transfer
0193     real_type pair_energy
0194         = inc_energy_ * std::exp(coeff_ * this->sample_scaled_energy(rng));
0195     CELER_ASSERT(pair_energy >= min_energy_ && pair_energy <= max_pair_energy_);
0196
0197     // Sample the energy partition between the electron and positron
0198     real_type r_max = (1
0199                        - 6 * inc_mass_sq_
0200                              / (total_energy_ * (total_energy_ - pair_energy)))
0201                       * std::sqrt(1 - min_pair_energy_ / pair_energy);
0202     real_type r = UniformRealDistribution(-r_max, r_max)(rng);
0203
0204     // Calculate the electron and positron energies
0205     PairEnergy result;
0206     real_type half_energy = pair_energy * real_type(0.5);
0207     result.electron = Energy((1 - r) * half_energy - electron_mass_);
0208     result.positron = Energy((1 + r) * half_energy - electron_mass_);
0209
0210     CELER_ENSURE(result.electron > zero_quantity());
0211     CELER_ENSURE(result.positron > zero_quantity());
0212     return result;
0213 }
0214
0215 //---------------------------------------------------------------------------//
0216 /*!
0217  * Sample the scaled energy and interpolate in log Z.
0218  */
0219 template<class Engine>
0220 CELER_FUNCTION real_type
0221 MuPPEnergyDistribution::sample_scaled_energy(Engine& rng) const
0222 {
0223     real_type u = generate_canonical(rng);
0224     LinearInterpolator<real_type> interp_energy{
0225         {0, this->calc_scaled_energy(logz_interp_.index, u)},
0226         {1, this->calc_scaled_energy(logz_interp_.index + 1, u)}};
0227     return interp_energy(logz_interp_.fraction);
0228 }
0229
0230 //---------------------------------------------------------------------------//
0231 /*!
0232  * Calculate the scaled energy for a given Z grid and sampled CDF value.
0233  */
0234 CELER_FUNCTION real_type
0235 MuPPEnergyDistribution::calc_scaled_energy(size_type z_idx, real_type u) const
0236 {
0237     CELER_EXPECT(z_idx < table_.grids.size());
0238     CELER_EXPECT(u >= 0 && u < 1);
0239
0240     TwodGridData const& cdf_grid = table_.grids[ItemId<TwodGridData>(z_idx)];
0241     auto calc_cdf
0242         = TwodGridCalculator(cdf_grid, table_.reals)(std::log(inc_energy_));
0243
0244     // Get the sampled CDF value between the y bounds
0245     real_type cdf = LinearInterpolator<real_type>{{0, calc_cdf(y_min_)},
0246                                                   {1, calc_cdf(y_max_)}}(u);
0247
0248     // Find the grid index of the sampled CDF value
0249     return InverseCdfFinder(NonuniformGrid(cdf_grid.y, table_.reals),
0250                             std::move(calc_cdf))(cdf);
0251 }
0252
0253 //---------------------------------------------------------------------------//
0254 }  // namespace celeritas