|
|
|||
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
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated by the 2.3.7 LXR engine. The LXR team |
|