Back to home page

EIC code displayed by LXR

 
 

    


Warning, file /include/celeritas/em/msc/detail/MscStepFromGeo.hh was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

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/MscStepFromGeo.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008 
0009 #include <cmath>
0010 
0011 #include "corecel/math/Algorithms.hh"
0012 #include "celeritas/Types.hh"
0013 #include "celeritas/em/data/UrbanMscData.hh"
0014 #include "celeritas/phys/Interaction.hh"
0015 
0016 namespace celeritas
0017 {
0018 namespace detail
0019 {
0020 //---------------------------------------------------------------------------//
0021 /*!
0022  * Convert the "geometrical" path traveled to the true path with MSC.
0023  *
0024  * The "true" step length is the physical path length taken along the
0025  * geometrical path length (which is straight without magnetic fields or curved
0026  * with them), accounting for the extra distance taken between
0027  * along-step elastic collisions.
0028  *
0029  * The transformation can be written as
0030  * \f[
0031  *     t(z) = \langle t \rangle = -\lambda_{1} \log(1 - \frac{z}{\lambda_{1}})
0032  * \f]
0033  * or \f$ t(z) = \frac{1}{\alpha} [ 1 - (1-\alpha w z)^{1/w}] \f$ if the
0034  * geom path is small, where \f$ w = 1 + \frac{1}{\alpha \lambda_{10}}\f$.
0035  *
0036  * \param true_path the proposed step before transportation.
0037  * \param gstep the proposed step after transportation.
0038  * \param alpha variable from UrbanMscStepLimit.
0039  */
0040 class MscStepFromGeo
0041 {
0042   public:
0043     //!@{
0044     //! \name Type aliases
0045     using MscParameters = UrbanMscParameters;
0046     //!@}
0047 
0048   public:
0049     // Construct with path length data
0050     inline CELER_FUNCTION MscStepFromGeo(MscParameters const& params_,
0051                                          MscStep const& step,
0052                                          real_type range,
0053                                          real_type lambda);
0054 
0055     // Calculate the true path
0056     inline CELER_FUNCTION real_type operator()(real_type gstep) const;
0057 
0058   private:
0059     // Urban MSC parameters
0060     MscParameters const& params_;
0061     //! Physical path length limited and transformed by MSC
0062     real_type true_step_;
0063     //! Scaled slope of the MFP
0064     real_type alpha_;
0065     //! Range for checking edge cases
0066     real_type range_;
0067     //! MFP at the start of step
0068     real_type lambda_;
0069 };
0070 
0071 //---------------------------------------------------------------------------//
0072 // INLINE DEFINITIONS
0073 //---------------------------------------------------------------------------//
0074 /*!
0075  * Construct from MSC data.
0076  */
0077 CELER_FUNCTION MscStepFromGeo::MscStepFromGeo(MscParameters const& params,
0078                                               MscStep const& step,
0079                                               real_type range,
0080                                               real_type lambda)
0081     : params_(params)
0082     , true_step_(step.true_path)
0083     , alpha_(step.alpha)
0084     , range_(range)
0085     , lambda_(lambda)
0086 {
0087     CELER_EXPECT(lambda_ > 0);
0088     CELER_EXPECT(true_step_ <= range_);
0089     // TODO: expect that the "skip sampling" conditions are false.
0090 }
0091 
0092 //---------------------------------------------------------------------------//
0093 /*!
0094  * Calculate and return the true path from the input "geometrical" step.
0095  *
0096  * This should only be applied if the step was geometry-limited. Otherwise, the
0097  * original interaction length is correct.
0098  */
0099 CELER_FUNCTION real_type MscStepFromGeo::operator()(real_type gstep) const
0100 {
0101     CELER_EXPECT(gstep >= 0 && gstep <= true_step_);
0102 
0103     if (gstep < params_.min_step_transform)
0104     {
0105         // Geometrical path length is true path length for a very small step
0106         return gstep;
0107     }
0108 
0109     real_type tstep = [&] {
0110         if (alpha_ == MscStep::small_step_alpha())
0111         {
0112             // Cross section was assumed to be constant over the step:
0113             // z = lambda * (1 - exp(-g / lambda))
0114             // => g = -lambda * log(1 - g / lambda)
0115             real_type tstep = -lambda_ * std::log1p(-gstep / lambda_);
0116             if (tstep < params_.min_step_transform)
0117             {
0118                 // Very small step: did not transform path length
0119                 return gstep;
0120             }
0121             return tstep;
0122         }
0123 
0124         real_type w = 1 + 1 / (alpha_ * lambda_);
0125         real_type x = alpha_ * w * gstep;  // = (1 - (1 - alpha * true)^w)
0126         // Range-limited step results in x = 1, which gives the correct
0127         // result in the inverted equation below for alpha = 1/range.
0128         // x >= 1 corresponds to the stopping MFP being <= 0: x=1 meaning
0129         // ending MFP of zero is correct when the step is the range. Precision
0130         // loss means this conversion also may result in x > 1, which we guard
0131         // against.
0132         x = min(x, real_type(1));
0133         // Near x=1, (1 - (1-x)^(1/w)) suffers from numerical precision loss;
0134         // the maximum value of the expression should be range_.
0135         // TODO: we should use the action ID to avoid applying this
0136         // transformation if not range-limited.
0137         real_type temp = 1 - fastpow(1 - x, 1 / w);
0138         real_type result = temp / alpha_;
0139         return min(result, range_);
0140     }();
0141 
0142     // The result can be no less than the geometry path (effectively no
0143     // scattering took place) and no more than the original calculated true
0144     // path (if the step is path-limited).
0145     // The result can be slightly outside the bounds due to numerical
0146     // imprecision and edge cases:
0147     // - a few ULP due to inexactness of log1p/expm1
0148     // - small travel distance results in roundoff error inside (1 -
0149     //   alpha/true)^w
0150     // - gstep is just above min_step_transform so that the updated
0151     // recalculation of
0152     //   tstep is just below min_step_transform
0153     return clamp(tstep, gstep, true_step_);
0154 }
0155 
0156 //---------------------------------------------------------------------------//
0157 }  // namespace detail
0158 }  // namespace celeritas