msc/detail/MscStepToGeo.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/msc/detail/MscStepToGeo.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <cmath>
0010
0011 #include "corecel/math/Algorithms.hh"
0012 #include "corecel/math/Quantity.hh"
0013 #include "corecel/math/SoftEqual.hh"
0014 #include "celeritas/Quantities.hh"
0015 #include "celeritas/Types.hh"
0016 #include "celeritas/em/data/UrbanMscData.hh"
0017 #include "celeritas/phys/Interaction.hh"
0018
0019 #include "UrbanMscHelper.hh"
0020
0021 namespace celeritas
0022 {
0023 namespace detail
0024 {
0025 //---------------------------------------------------------------------------//
0026 /*!
0027  * Convert the "true" path traveled to a geometrical approximation.
0028  *
0029  * This takes the physical step limit---whether limited by range, physics
0030  * interaction, or an MSC step limiter---and converts it to a "geometrical"
0031  * path length, which is a smooth curve (straight if in the absence of a
0032  * magnetic field).
0033  *
0034  * The mean value of the geometrical path length \f$ z \f$ (the first moment)
0035  * corresponding to a given true path length \f$ t \f$ is given by
0036  * \f[
0037  *     \langle z \rangle = \lambda_{1} [ 1 - \exp({-\frac{t}{\lambda_{1}}})]
0038  * \f]
0039  * where \f$\lambda_{1}\f$ is the first transport mean free path. Due to the
0040  * fact that \f$\lambda_{1}\f$ depends on the kinetic energy of the path and
0041  * decreases along the step, the path length correction is approximated as
0042  * \f[
0043  *     \lambda_{1} (t') = \lambda_{1,i} (1 - \alpha t') \,, 0 \le t' \le t
0044  * \f]
0045  * where the transport MFP at the start of the step is \f$ \lambda_{1,i} \def
0046  * \lambda_{1}(0)\f$ , and a linearity coefficient is defined as \f$ \alpha
0047  * \def \frac{\lambda_{1,i} - \lambda_{1,f}}{t \lambda_{1,i}} \f$, where  \f$
0048  * \lambda_{1,f}  \def \lambda_{1}(t) \f$. For low-energy or range-limited
0049  * steps, \f$ \alpha = 1/r_0 \f$ because the cross section at zero energy is
0050  * effectively infinite.
0051  *
0052  * Since the MSC cross section generally decreases as the energy increases,
0053  * \f$ \lambda_{1,i} \f$ will be larger than \f$ \lambda_{1,f} \f$ and
0054  * \f$ \alpha \f$ will be positive.
0055  * However, in the Geant4 Urban MSC model, different methods
0056  * are used to calculate the cross section above and below 10 MeV. In the
0057  * higher energy region the cross sections are identical for electrons and
0058  * positrons, resulting in a discontinuity in the positron cross section at 10
0059  * MeV. This means on fine energy grids it's possible for the cross section to
0060  * be *increasing* with energy just above the 10 MeV threshold and therefore
0061  * for \f$ \alpha \f$ is negative.
0062  *
0063  * The resulting geometrical step length *can* be greater than 1 MFP: it's the
0064  * MSC step limiter's job to add additional restrictions.
0065  *
0066  * \note This performs the same method as in ComputeGeomPathLength of
0067  * G4UrbanMscModel of the Geant4 10.7 release.
0068  */
0069 class MscStepToGeo
0070 {
0071   public:
0072     //!@{
0073     //! \name Type aliases
0074     using Energy = units::MevEnergy;
0075     using Mass = units::MevMass;
0076     using UrbanMscRef = NativeCRef<UrbanMscData>;
0077     //!@}
0078
0079     struct result_type
0080     {
0081         real_type step{};  //!< Geometrical step length [len]
0082         real_type alpha{0};  //!< Scaled MFP slope [1/len]
0083     };
0084
0085   public:
0086     // Construct from MSC data
0087     inline CELER_FUNCTION MscStepToGeo(UrbanMscRef const& shared,
0088                                        UrbanMscHelper const& helper,
0089                                        Energy energy,
0090                                        real_type lambda,
0091                                        real_type range);
0092
0093     // Calculate the geometrical step
0094     inline CELER_FUNCTION result_type operator()(real_type tstep) const;
0095
0096   private:
0097     UrbanMscRef const& shared_;
0098     UrbanMscHelper const& helper_;
0099     real_type energy_;
0100     real_type lambda_;
0101     real_type range_;
0102 };
0103
0104 //---------------------------------------------------------------------------//
0105 // INLINE DEFINITIONS
0106 //---------------------------------------------------------------------------//
0107 /*!
0108  * Construct from MSC data.
0109  */
0110 CELER_FUNCTION MscStepToGeo::MscStepToGeo(UrbanMscRef const& shared,
0111                                           UrbanMscHelper const& helper,
0112                                           Energy energy,
0113                                           real_type lambda,
0114                                           real_type range)
0115     : shared_(shared)
0116     , helper_(helper)
0117     , energy_(energy.value())
0118     , lambda_(lambda)
0119     , range_(range)
0120 {
0121     CELER_EXPECT(energy_ > 0);
0122     CELER_EXPECT(lambda_ > 0);
0123     CELER_EXPECT(range_ > 0);
0124 }
0125
0126 //---------------------------------------------------------------------------//
0127 /*!
0128  * Calculate the geometry step length for a given true step length.
0129  */
0130 CELER_FUNCTION auto
0131 MscStepToGeo::operator()(real_type tstep) const -> result_type
0132 {
0133     CELER_EXPECT(tstep >= 0 && tstep <= range_);
0134
0135     result_type result;
0136     result.alpha = MscStep::small_step_alpha();
0137     if (tstep < shared_.params.min_step_transform)
0138     {
0139         // Very small step: do not transform step length
0140         result.step = tstep;
0141     }
0142     else if (tstep < range_ * shared_.params.small_range_frac)
0143     {
0144         // Small enough distance to assume cross section is constant
0145         // over the step: z = lambda * (1 - exp(-tau))
0146         result.step = -lambda_ * std::expm1(-tstep / lambda_);
0147     }
0148     else
0149     {
0150         // Eq 8.8: mfp_slope = (lambda_f / lambda_i) = (1 - alpha * tstep)
0151         real_type mfp_slope;
0152         if (energy_ < value_as<Mass>(shared_.electron_mass) || tstep == range_)
0153         {
0154             // Low energy or range-limited step
0155             // (For range-limited step, the final cross section and thus MFP
0156             // slope are zero).
0157             result.alpha = 1 / range_;
0158             // Use max to avoid slightly negative slope due to roundoff for
0159             // range-limited steps
0160             mfp_slope = max<real_type>(1 - result.alpha * tstep, 0);
0161         }
0162         else
0163         {
0164             // Calculate the energy at the end of a physics-limited step
0165             real_type rfinal = range_ - tstep;
0166             Energy endpoint_energy = helper_.calc_inverse_range(rfinal);
0167             real_type lambda1 = helper_.calc_msc_mfp(endpoint_energy);
0168
0169             // Calculate the geometric path assuming the cross section is
0170             // linear between the start and end energy. Eq 8.10+1
0171             result.alpha = (lambda_ - lambda1) / (lambda_ * tstep);
0172             mfp_slope = lambda1 / lambda_;
0173         }
0174
0175         // Eq 8.10 with simplifications
0176         real_type w = 1 + 1 / (result.alpha * lambda_);
0177         result.step = (1 - fastpow(mfp_slope, w)) / (result.alpha * w);
0178     }
0179
0180     // For extremely large lambda we can slightly exceed the step
0181     CELER_ENSURE(result.step <= tstep || soft_equal(result.step, tstep));
0182     result.step = min(result.step, tstep);
0183     return result;
0184 }
0185
0186 //---------------------------------------------------------------------------//
0187 }  // namespace detail
0188 }  // namespace celeritas