celeritas/grid/InverseRangeCalculator.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/grid/InverseRangeCalculator.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <cmath>
0010
0011 #include "corecel/Assert.hh"
0012 #include "corecel/cont/LdgSpan.hh"
0013 #include "corecel/data/Collection.hh"
0014 #include "corecel/grid/Interpolator.hh"
0015 #include "corecel/grid/NonuniformGrid.hh"
0016 #include "corecel/grid/SplineInterpolator.hh"
0017 #include "corecel/grid/UniformGrid.hh"
0018 #include "corecel/grid/UniformGridData.hh"
0019 #include "corecel/math/Algorithms.hh"
0020 #include "celeritas/Quantities.hh"
0021
0022 namespace celeritas
0023 {
0024 //---------------------------------------------------------------------------//
0025 /*!
0026  * Calculate the energy that would limit a particle to a particular range.
0027  *
0028  * This should provide the inverse of the result of \c RangeCalculator. The
0029  * given \c range is not allowed to be greater than the maximum range in the
0030  * physics data.
0031  *
0032  * The range must be monotonically increasing in energy, since it's defined as
0033  * the integral of the inverse of the stopping power (which is always
0034  * positive). For ranges shorter than the minimum energy in the table, the
0035  * resulting energy is scaled:
0036  * \f[
0037     E = E_\mathrm{min} \left( \frac{r}{r_\mathrm{min}} \right)^2
0038  * \f]
0039  * This scaling is the inverse of the off-the-end energy scaling in the
0040  * RangeCalculator.
0041  *
0042  * \todo Construct with \c UniformGridRecord
0043  */
0044 class InverseRangeCalculator
0045 {
0046   public:
0047     //!@{
0048     //! \name Type aliases
0049     using Energy = units::MevEnergy;
0050     using Values
0051         = Collection<real_type, Ownership::const_reference, MemSpace::native>;
0052     //!@}
0053
0054   public:
0055     // Construct from state-independent data
0056     inline CELER_FUNCTION InverseRangeCalculator(UniformGridRecord const& grid,
0057                                                  Values const& values);
0058
0059     // Find and interpolate from the energy
0060     inline CELER_FUNCTION Energy operator()(real_type range) const;
0061
0062   private:
0063     UniformGrid log_energy_;
0064     NonuniformGrid<real_type> range_;
0065     LdgSpan<real_type const> deriv_;
0066 };
0067
0068 //---------------------------------------------------------------------------//
0069 // INLINE DEFINITIONS
0070 //---------------------------------------------------------------------------//
0071 /*!
0072  * Construct from range data.
0073  *
0074  * The range is expected to be monotonically increasing with energy.
0075  * Lower-energy particles have shorter ranges.
0076  */
0077 CELER_FUNCTION InverseRangeCalculator::InverseRangeCalculator(
0078     UniformGridRecord const& grid, Values const& values)
0079     : log_energy_(grid.grid)
0080     , range_(grid.value, values)
0081     , deriv_(values[grid.derivative])
0082 {
0083     CELER_EXPECT(range_.size() == log_energy_.size());
0084 }
0085
0086 //---------------------------------------------------------------------------//
0087 /*!
0088  * Calculate the energy of a particle that has the given range.
0089  */
0090 CELER_FUNCTION auto InverseRangeCalculator::operator()(real_type range) const
0091     -> Energy
0092 {
0093     CELER_EXPECT(range >= 0 && range <= range_.back());
0094
0095     if (range < range_.front())
0096     {
0097         // Very short range:  this corresponds to "energy < emin" for range
0098         // calculation: range = r[0] * sqrt(E / E[0])
0099         return Energy{std::exp(log_energy_.front())
0100                       * ipow<2>(range / range_.front())};
0101     }
0102     // Range should *never* exceed the longest range (highest energy) since
0103     // that should have limited the step
0104     if (CELER_UNLIKELY(range >= range_.back()))
0105     {
0106         CELER_ASSERT(range == range_.back());
0107         return Energy{std::exp(log_energy_.back())};
0108     }
0109
0110     // Search for lower bin index
0111     auto idx = range_.find(range);
0112     CELER_ASSERT(idx + 1 < log_energy_.size());
0113
0114     real_type result;
0115     if (deriv_.empty())
0116     {
0117         // Interpolate: 'x' = range, y = log energy
0118         result = LinearInterpolator<real_type>(
0119             {range_[idx], std::exp(log_energy_[idx])},
0120             {range_[idx + 1], std::exp(log_energy_[idx + 1])})(range);
0121     }
0122     else
0123     {
0124         // Use cubic spline interpolation
0125         result = SplineInterpolator<real_type>(
0126             {range_[idx], std::exp(log_energy_[idx]), deriv_[idx]},
0127             {range_[idx + 1], std::exp(log_energy_[idx + 1]), deriv_[idx + 1]})(
0128             range);
0129     }
0130     return Energy{result};
0131 }
0132
0133 //---------------------------------------------------------------------------//
0134 }  // namespace celeritas