msc/detail/UrbanMscScatter.hh

0001 // Copyright 2021-2024 UT-Battelle, LLC, and other Celeritas developers.
0002 // See the top-level COPYRIGHT file for details.
0003 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0004 //---------------------------------------------------------------------------//
0005 //! \file celeritas/em/msc/detail/UrbanMscScatter.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <cmath>
0010
0011 #include "corecel/Macros.hh"
0012 #include "corecel/Types.hh"
0013 #include "corecel/math/Algorithms.hh"
0014 #include "corecel/math/ArrayOperators.hh"
0015 #include "corecel/math/ArrayUtils.hh"
0016 #include "celeritas/Constants.hh"
0017 #include "celeritas/Quantities.hh"
0018 #include "celeritas/Types.hh"
0019 #include "celeritas/Units.hh"
0020 #include "celeritas/em/data/UrbanMscData.hh"
0021 #include "celeritas/geo/GeoTrackView.hh"
0022 #include "celeritas/grid/PolyEvaluator.hh"
0023 #include "celeritas/mat/MaterialView.hh"
0024 #include "celeritas/phys/Interaction.hh"
0025 #include "celeritas/phys/ParticleTrackView.hh"
0026 #include "celeritas/phys/PhysicsTrackView.hh"
0027 #include "celeritas/random/distribution/BernoulliDistribution.hh"
0028 #include "celeritas/random/distribution/GenerateCanonical.hh"
0029 #include "celeritas/random/distribution/UniformRealDistribution.hh"
0030
0031 #include "MscStepFromGeo.hh"
0032 #include "UrbanMscHelper.hh"
0033 #include "UrbanPositronCorrector.hh"
0034
0035 namespace celeritas
0036 {
0037 namespace detail
0038 {
0039 //---------------------------------------------------------------------------//
0040 /*!
0041  * Sample angular change and lateral displacement with the Urban multiple
0042  * scattering model.
0043
0044  * \note This code performs the same method as in
0045  * G4VMultipleScattering::AlongStepDoIt and G4UrbanMscModel::SampleScattering
0046  * of the Geant4 10.7 release.
0047  */
0048 class UrbanMscScatter
0049 {
0050   public:
0051     //!@{
0052     //! \name Type aliases
0053     using Energy = units::MevEnergy;
0054     using Mass = units::MevMass;
0055     using UrbanMscRef = NativeCRef<UrbanMscData>;
0056     //!@}
0057
0058   public:
0059     // TODO: improve this
0060     static inline CELER_FUNCTION real_type
0061     calc_displacement(real_type true_path, real_type geom_path);
0062
0063     // Construct with shared and state data
0064     inline CELER_FUNCTION UrbanMscScatter(UrbanMscRef const& shared,
0065                                           UrbanMscHelper const& helper,
0066                                           ParticleTrackView const& particle,
0067                                           PhysicsTrackView const& physics,
0068                                           MaterialView const& material,
0069                                           Real3 const& dir,
0070                                           real_type safety,
0071                                           MscStep const& input);
0072
0073     // Sample the final true step length, position and direction by msc
0074     template<class Engine>
0075     inline CELER_FUNCTION MscInteraction operator()(Engine& rng);
0076
0077   private:
0078     //// DATA ////
0079
0080     // Shared constant data
0081     UrbanMscRef const& shared_;
0082     // Urban MSC material data
0083     UrbanMscMaterialData const& msc_;
0084     // Urban MSC helper class
0085     UrbanMscHelper const& helper_;
0086     // Material data
0087     MaterialView const& material_;
0088
0089     real_type inc_energy_;
0090     bool is_positron_;
0091     Real3 const& inc_direction_;
0092     real_type safety_;
0093
0094     // Results from UrbanMSCStepLimit
0095     bool is_displaced_;
0096     real_type geom_path_;
0097     real_type true_path_;
0098     real_type limit_min_;
0099
0100     // Calculated values for sampling
0101     bool skip_sampling_;
0102     real_type end_energy_;
0103     real_type tau_{0};
0104     real_type xmean_{0};
0105     real_type x2mean_{0};
0106     real_type theta0_{-1};
0107
0108     //// COMMON PROPERTIES ////
0109
0110     //! The minimum step length for geometry 0.05 nm
0111     static CELER_CONSTEXPR_FUNCTION real_type geom_min()
0112     {
0113         return 5e-9 * units::centimeter;
0114     }
0115
0116     //// HELPER FUNCTIONS ////
0117
0118     // Sample the angle, cos(theta), of the multiple scattering
0119     template<class Engine>
0120     inline CELER_FUNCTION real_type sample_cos_theta(Engine& rng) const;
0121
0122     // Sample consine(theta) with a large angle scattering
0123     template<class Engine>
0124     inline CELER_FUNCTION real_type simple_scattering(Engine& rng) const;
0125
0126     // Calculate the theta0 of the Highland formula
0127     inline CELER_FUNCTION real_type compute_theta0() const;
0128
0129     // Update direction and position after the multiple scattering
0130     template<class Engine>
0131     inline CELER_FUNCTION Real3 sample_displacement_dir(Engine& rng,
0132                                                         real_type phi) const;
0133 };
0134
0135 //---------------------------------------------------------------------------//
0136 // INLINE DEFINITIONS
0137 //---------------------------------------------------------------------------//
0138 CELER_FUNCTION real_type UrbanMscScatter::calc_displacement(real_type geom_path,
0139                                                             real_type true_path)
0140 {
0141     CELER_EXPECT(true_path >= geom_path);
0142
0143     // true^2 - geo^2
0144     real_type rmax2 = diffsq(true_path, geom_path);
0145
0146     // 0.73 is (roughly) the expected value of a distribution of the mean
0147     // radius given rmax "based on single scattering results"
0148     // https://github.com/Geant4/geant4/blame/28a70706e0edf519b16e864ebf1d2f02a00ba596/source/processes/electromagnetic/standard/src/G4UrbanMscModel.cc#L1142
0149     constexpr real_type mean_radius_frac{0.73};
0150
0151     return mean_radius_frac * std::sqrt(rmax2);
0152 }
0153
0154 //---------------------------------------------------------------------------//
0155 /*!
0156  * Construct with shared and state data.
0157  *
0158  * This function also precalculates distribution-independent quantities, e.g.
0159  * converting the geometrical path length to the true path.
0160  */
0161 CELER_FUNCTION
0162 UrbanMscScatter::UrbanMscScatter(UrbanMscRef const& shared,
0163                                  UrbanMscHelper const& helper,
0164                                  ParticleTrackView const& particle,
0165                                  PhysicsTrackView const& physics,
0166                                  MaterialView const& material,
0167                                  Real3 const& dir,
0168                                  real_type safety,
0169                                  MscStep const& input)
0170     : shared_(shared)
0171     , msc_(shared.material_data[material.material_id()])
0172     , helper_(helper)
0173     , material_(material)
0174     , inc_energy_(value_as<Energy>(particle.energy()))
0175     , is_positron_(particle.particle_id() == shared.ids.positron)
0176     , inc_direction_(dir)
0177     , safety_(safety)
0178     , is_displaced_(input.is_displaced)
0179     , geom_path_(input.geom_path)
0180     , true_path_(input.true_path)
0181     , limit_min_(physics.msc_range().limit_min)
0182 {
0183     CELER_EXPECT(particle.particle_id() == shared.ids.electron
0184                  || particle.particle_id() == shared.ids.positron);
0185     CELER_EXPECT(safety_ >= 0);
0186     CELER_EXPECT(geom_path_ > 0);
0187     CELER_EXPECT(true_path_ >= geom_path_);
0188     CELER_EXPECT(limit_min_ >= UrbanMscParameters::limit_min_fix()
0189                  || !is_displaced_);
0190     CELER_EXPECT(!is_displaced_ || safety > 0);
0191
0192     skip_sampling_ = [this, &helper, &physics] {
0193         if (true_path_ == physics.dedx_range())
0194         {
0195             // Range-limited step (particle stops)
0196             // TODO: probably redundant with low 'end energy'
0197             return true;
0198         }
0199         if (true_path_ < shared_.params.geom_limit)
0200         {
0201             // Very small step (NOTE: with the default values in UrbanMscData,
0202             // this is redundant witih the tau_small comparison below if MFP >=
0203             // 0.005 cm)
0204             return true;
0205         }
0206
0207         // Lazy calculation of end energy
0208         end_energy_ = value_as<Energy>(helper.calc_end_energy(true_path_));
0209
0210         if (Energy{end_energy_} < shared_.params.min_sampling_energy())
0211         {
0212             // Ending energy is very low
0213             return true;
0214         }
0215         if (true_path_ <= helper_.msc_mfp() * shared_.params.tau_small)
0216         {
0217             // Very small MFP travelled
0218             return true;
0219         }
0220         return false;
0221     }();
0222
0223     // TODO: there are several different sampling strategies for angle change:
0224     // - very small step/very low energy endpoint: no scattering
0225     // - very small mfp: (probably impossible because of condition above):
0226     //   forward scatter
0227     // - very large mfp: exiting angle is isotropic
0228     // - large energy loss: "simple_scattering"
0229
0230     if (!skip_sampling_)
0231     {
0232         // Calculate number of mean free paths traveled
0233         tau_ = true_path_ / [this, &helper] {
0234             // Calculate the average MFP assuming the cross section varies
0235             // linearly over the step
0236             real_type lambda = helper_.msc_mfp();
0237             real_type lambda_end = helper.calc_msc_mfp(Energy{end_energy_});
0238             if (std::fabs(lambda - lambda_end) < lambda * real_type(0.01))
0239             {
0240                 // Cross section is almost constant over the step: avoid
0241                 // numerical explosion
0242                 return helper_.msc_mfp();
0243             }
0244             return (lambda - lambda_end) / std::log(lambda / lambda_end);
0245         }();
0246
0247         if (tau_ < shared_.params.tau_big)
0248         {
0249             // Eq. 8.2 and \f$ \cos^2\theta \f$ term in Eq. 8.3 in PRM
0250             xmean_ = std::exp(-tau_);
0251             x2mean_ = (1 + 2 * std::exp(real_type(-2.5) * tau_)) / 3;
0252
0253             // MSC "true path" step limit
0254             if (CELER_UNLIKELY(limit_min_ == 0))
0255             {
0256                 // Unlikely: MSC range cache wasn't initialized by
0257                 // UrbanMscStepLimit
0258                 CELER_ASSERT(!is_displaced_);
0259                 limit_min_ = UrbanMscParameters::limit_min();
0260             }
0261             limit_min_ = min(limit_min_, physics.scalars().lambda_limit);
0262
0263             // TODO: theta0_ calculation could be done externally, eliminating
0264             // many of the class member data
0265             theta0_ = this->compute_theta0();
0266
0267             if (theta0_ < real_type(1e-8))
0268             {
0269                 // Arbitrarily (?) small angle change (theta_0^2 < 1e-16):
0270                 // skip sampling angular distribution if width of direction
0271                 // distribution is too narrow
0272                 if (!is_displaced_)
0273                 {
0274                     // No angular sampling and no displacement => no change
0275                     skip_sampling_ = true;
0276                 }
0277                 else
0278                 {
0279                     theta0_ = 0;
0280                 }
0281             }
0282         }
0283     }
0284 }
0285
0286 //---------------------------------------------------------------------------//
0287 /*!
0288  * Sample the angular distribution and the lateral displacement by multiple
0289  * scattering.
0290  */
0291 template<class Engine>
0292 CELER_FUNCTION auto UrbanMscScatter::operator()(Engine& rng) -> MscInteraction
0293 {
0294     if (skip_sampling_)
0295     {
0296         // Do not sample scattering at the last or at a small step
0297         return {inc_direction_, {0, 0, 0}, MscInteraction::Action::unchanged};
0298     }
0299
0300     // Sample polar angle cosine
0301     real_type costheta = [this, &rng] {
0302         if (theta0_ <= 0)
0303         {
0304             // Very small outgoing angular distribution
0305             return real_type{1};
0306         }
0307         if (tau_ >= shared_.params.tau_big)
0308         {
0309             // Long mean free path: exiting direction is isotropic
0310             UniformRealDistribution<real_type> sample_isotropic(-1, 1);
0311             return sample_isotropic(rng);
0312         }
0313         if (2 * end_energy_ < inc_energy_ || theta0_ > constants::pi / 6)
0314         {
0315             // Large energy loss over the step or large angle distribution
0316             // width
0317             return this->simple_scattering(rng);
0318         }
0319         // No special cases match:
0320         return this->sample_cos_theta(rng);
0321     }();
0322     CELER_ASSERT(std::fabs(costheta) <= 1);
0323
0324     // Sample azimuthal angle, used for displacement and exiting angle
0325     real_type phi
0326         = UniformRealDistribution<real_type>(0, 2 * constants::pi)(rng);
0327
0328     MscInteraction result;
0329     result.action = MscInteraction::Action::scattered;
0330     // This should only be needed to silence compiler warning, since the
0331     // displacement should be ignored since our action result is
0332     // 'scattered'
0333     result.displacement = {0, 0, 0};
0334
0335     if (is_displaced_)
0336     {
0337         // Calculate displacement length
0338         real_type length = this->calc_displacement(geom_path_, true_path_);
0339         // Don't displace further than safety (minus a tolerance)
0340         length = min(length, (1 - shared_.params.safety_tol) * safety_);
0341
0342         if (length >= shared_.params.geom_limit)
0343         {
0344             // Displacement distance is large enough to worry about
0345             result.displacement = this->sample_displacement_dir(rng, phi);
0346             result.displacement *= length;
0347             result.action = MscInteraction::Action::displaced;
0348         }
0349     }
0350
0351     // Calculate direction and return
0352     result.direction = rotate(from_spherical(costheta, phi), inc_direction_);
0353     return result;
0354 }
0355
0356 //---------------------------------------------------------------------------//
0357 /*!
0358  * Sample the scattering angle at the end of the true step length.
0359  *
0360  * The scattering angle \f$\theta\f$ and true step length, \f$t\f$ are
0361  * described in G4UrbanMscModel::SampleCosineTheta of the
0362  * Geant4 10.7 release. See also, CERN-OPEN-2006-077 by L. Urban.
0363  *
0364  * The mean value of \f$u = \cos\theta\f$ follows \f$\exp(-t/\lambda_{1})\f$
0365  * and the variance is written as \f$\frac{1+2e^{-\kappa r}}{3} - e^{-2r}\f$
0366  * where \f$r = t/\lambda_{1}\f$ and \f$\kappa = \lambda_{1}/\lambda_{2}\f$.
0367  * The \f$\cos\theta\f$ is sampled according to a model function of \f$u\f$,
0368  * \f[
0369  *   g(u) = q [ p g_1(u) + (1-p) g_2(u) ] - (1 - q) g_3(u)
0370  * \f]
0371  * where \f$p,q \in [0,1]\f$ and the functions \f$g_i\f$ have been chosen as
0372  * \f[
0373  *   g_1(u) = c_1 e^{-a(1-u)},
0374  *   g_2(u) = \frac{c_2}{(b-u)^d},
0375  *   g_3(u) = c_3
0376  * \f]
0377  * with normalization constants, \f$d\f$. For small angles, \f$g_1\f$ is
0378  * nearly Gaussian, \f$ \exp(-\frac{\theta^{2}}{2\theta_{0}^{2}}), \f$
0379  * if \f$\theta_0 \approx 1/1\f$, while \f$g_2\f$ has a Rutherford-like tail
0380  * for large \f$\theta\f$, if \f$b \approx 1\f$ and \f$d\f$ is not far from 2.
0381  */
0382 template<class Engine>
0383 CELER_FUNCTION real_type UrbanMscScatter::sample_cos_theta(Engine& rng) const
0384 {
0385     CELER_EXPECT(xmean_ > 0 && x2mean_ > 0);
0386
0387     // Evaluate parameters for the tail distribution
0388     real_type xsi = [this] {
0389         using PolyQuad = PolyEvaluator<real_type, 2>;
0390
0391         real_type maxtau
0392             = true_path_ < limit_min_ ? limit_min_ / helper_.msc_mfp() : tau_;
0393         // note: 0 < u <= sqrt(2) when shared_.params.tau_big == 8
0394         real_type u = fastpow(maxtau, 1 / real_type(6));
0395         // Number of radiation lengths traveled by the average MFP over this
0396         // step
0397         real_type radlen_mfp = true_path_
0398                                / (tau_ * material_.radiation_length());
0399         real_type result = PolyQuad(msc_.tail_coeff)(u)
0400                            + msc_.tail_corr * std::log(radlen_mfp);
0401         // The tail should not be too big
0402         return max(result, real_type(1.9));
0403     }();
0404
0405     real_type ea = std::exp(-xsi);
0406     // Mean of cos\theta computed from the distribution g_1(cos\theta)
0407     // small theta => x = theta0^2
0408     // large xsi => xmean_1 = 1 - x
0409     // small tau => xmean = 1
0410     real_type x = ipow<2>(2 * std::sin(real_type(0.5) * theta0_));
0411     real_type xmean_1 = 1 - x * (1 + (xsi * ea) / (1 - ea));
0412
0413     if (xmean_1 <= real_type(0.999) * xmean_)
0414     {
0415         return this->simple_scattering(rng);
0416     }
0417
0418     // From continuity of derivatives
0419     real_type c = [xsi] {
0420         if (std::fabs(xsi - 3) < real_type(0.001))
0421         {
0422             return real_type(3.001);
0423         }
0424         else if (std::fabs(xsi - 2) < real_type(0.001))
0425         {
0426             return real_type(2.001);
0427         }
0428         return xsi;
0429     }();
0430     real_type b1 = 2 + (c - xsi) * x;
0431     real_type d = fastpow(c * x / b1, c - 1);
0432     real_type x0 = 1 - xsi * x;
0433
0434     // Mean of cos\theta computed from the distribution g_2(cos\theta)
0435     real_type xmean_2 = (x0 + d - (c * x - b1 * d) / (c - 2)) / (1 - d);
0436
0437     real_type f2x0 = (c - 1) / (c * (1 - d));
0438     real_type prob = f2x0 / (ea / (1 - ea) + f2x0);
0439
0440     // Eq. 8.14 in the PRM: note that can be greater than 1
0441     real_type qprob = xmean_ / (prob * xmean_1 + (1 - prob) * xmean_2);
0442     // Sampling of cos(theta)
0443     if (generate_canonical(rng) >= qprob)
0444     {
0445         // Sample \f$ \cos\theta \f$ from \f$ g_3(\cos\theta) \f$
0446         return UniformRealDistribution<real_type>(-1, 1)(rng);
0447     }
0448
0449     // Note: prob is sometime a little greater than one
0450     if (generate_canonical(rng) < prob)
0451     {
0452         // Sample \f$ \cos\theta \f$ from \f$ g_1(\cos\theta) \f$
0453         UniformRealDistribution<real_type> sample_inner(ea, 1);
0454         return 1 + std::log(sample_inner(rng)) * x;
0455     }
0456     else
0457     {
0458         // Sample \f$ \cos\theta \f$ from \f$ g_2(\cos\theta) \f$
0459         real_type var = (1 - d) * generate_canonical(rng);
0460         if (var < real_type(0.01) * d)
0461         {
0462             var /= (d * (c - 1));
0463             return -1
0464                    + var * (1 - real_type(0.5) * var * c) * (2 + (c - xsi) * x);
0465         }
0466         else
0467         {
0468             return x * (c - xsi - c * fastpow(var + d, -1 / (c - 1))) + 1;
0469         }
0470     }
0471 }
0472
0473 //---------------------------------------------------------------------------//
0474 /*!
0475  * Sample the large angle scattering using 2 model functions.
0476  *
0477  * \param rng Random number generator
0478  */
0479 template<class Engine>
0480 CELER_FUNCTION real_type UrbanMscScatter::simple_scattering(Engine& rng) const
0481 {
0482     real_type a = (2 * xmean_ + 9 * x2mean_ - 3)
0483                   / (2 * xmean_ - 3 * x2mean_ + 1);
0484     BernoulliDistribution sample_pow{(a + 2) * xmean_ / a};
0485
0486     // Sample cos(theta)
0487     real_type result{};
0488     do
0489     {
0490         real_type rdm = generate_canonical(rng);
0491         result = 2 * (sample_pow(rng) ? fastpow(rdm, 1 / (a + 1)) : rdm) - 1;
0492     } while (std::fabs(result) > 1);
0493
0494     return result;
0495 }
0496
0497 //---------------------------------------------------------------------------//
0498 /*!
0499  * Calculate the width of an approximate Gaussian projected angle distribution
0500  * using a modified Highland-Lynch-Dahl formula.
0501  *
0502  * All particles take the width
0503  * of the central part from a parameterization similar to the original Highland
0504  * formula, Particle Physics Booklet, July 2002, eq. 26.10.
0505  * \f[
0506  *   \theta_0 = \frac{13.6\rm{MeV}}{\beta c p} z_{ch} \sqrt{\ln(t/X_o)} c
0507  * \f]
0508  * where \f$p, \beta\c, z_{ch}\f$, \f$t/X_0\f$ and \f$c\f$ are the momentum,
0509  * velocity, charge number of the incident particle, the true path length in
0510  * radiation length unit and the correction term, respectively. For details,
0511  * see the section 8.1.5 of the Geant4 10.7 Physics Reference Manual.
0512  */
0513 CELER_FUNCTION
0514 real_type UrbanMscScatter::compute_theta0() const
0515 {
0516     real_type const mass = value_as<Mass>(shared_.electron_mass);
0517     real_type true_path = max(limit_min_, true_path_);
0518     real_type y = true_path / material_.radiation_length();
0519
0520     // Correction for the positron
0521     if (is_positron_)
0522     {
0523         detail::UrbanPositronCorrector calc_correction{material_.zeff()};
0524         y *= calc_correction(std::sqrt(inc_energy_ * end_energy_) / mass);
0525     }
0526     CELER_ASSERT(y > 0);
0527
0528     // TODO for hadrons: multiply abs(charge)
0529     real_type invbetacp
0530         = std::sqrt((inc_energy_ + mass) * (end_energy_ + mass)
0531                     / (inc_energy_ * (inc_energy_ + 2 * mass) * end_energy_
0532                        * (end_energy_ + 2 * mass)));
0533     constexpr units::MevEnergy c_highland{13.6};
0534     real_type theta0 = c_highland.value() * std::sqrt(y) * invbetacp;
0535
0536     // Correction factor from e- scattering data
0537     theta0 *= PolyEvaluator<real_type, 1>(msc_.theta_coeff)(std::log(y));
0538
0539     if (true_path_ < limit_min_)
0540     {
0541         // Apply correction if true path is very small
0542         theta0 *= std::sqrt(true_path_ / limit_min_);
0543     }
0544
0545     // Very small path lengths can result in a negative e- scattering
0546     // correction: clamp to zero so that too-small paths result in no change
0547     // in angle
0548     return max<real_type>(theta0, 0);
0549 }
0550
0551 //---------------------------------------------------------------------------//
0552 /*!
0553  * Sample the displacement direction using G4UrbanMscModel::SampleDisplacement
0554  * (simple and fast sampling based on single scattering results) and update
0555  * direction and position of the particle.
0556  *
0557  * A simple distribution for the unit direction on the lateral (x-y) plane,
0558  * \f$ Phi = \phi \pm \psi \f$ where \f$ psi \sim \exp(-\beta*v) \f$ and
0559  * \f$\beta\f$ is determined from the requirement that the distribution should
0560  * give the same mean value that is obtained from the single scattering
0561  * simulation.
0562  *
0563  * \param rng Random number generator
0564  * \param phi the azimuthal angle of the multiple scattering.
0565  */
0566 template<class Engine>
0567 CELER_FUNCTION Real3
0568 UrbanMscScatter::sample_displacement_dir(Engine& rng, real_type phi) const
0569 {
0570     // Sample a unit direction of the displacement
0571     constexpr real_type cbeta = 2.160;
0572     // cbeta1 = 1 - std::exp(-cbeta * constants::pi);
0573     constexpr real_type cbeta1 = 0.9988703417569197;
0574
0575     real_type psi = -std::log(1 - generate_canonical(rng) * cbeta1) / cbeta;
0576     phi += BernoulliDistribution(0.5)(rng) ? psi : -psi;
0577
0578     Real3 displacement{std::cos(phi), std::sin(phi), 0};
0579
0580     // Rotate along the incident particle direction
0581     displacement = rotate(displacement, inc_direction_);
0582     return displacement;
0583 }
0584
0585 //---------------------------------------------------------------------------//
0586 }  // namespace detail
0587 }  // namespace celeritas