msc/detail/MscStepToGeo.hh

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