em/xs/NuclearFormFactors.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/xs/NuclearFormFactors.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <cmath>
0010
0011 #include "corecel/math/Algorithms.hh"
0012 #include "celeritas/Quantities.hh"
0013 #include "celeritas/phys/AtomicNumber.hh"
0014
0015 namespace celeritas
0016 {
0017 //---------------------------------------------------------------------------//
0018 //! Helper traits used in nuclear form factors
0019 struct NuclearFormFactorTraits
0020 {
0021     //!@{
0022     //! \name Type aliases
0023     using AtomicMassNumber = AtomicNumber;
0024     using Momentum = units::MevMomentum;
0025     using MomentumSq = units::MevMomentumSq;
0026     using InvMomentum = RealQuantity<UnitInverse<Momentum::unit_type>>;
0027     using InvMomentumSq = RealQuantity<UnitInverse<MomentumSq::unit_type>>;
0028     using FFType = NuclearFormFactorType;
0029     //!@}
0030
0031     //! Momentum transfer prefactor: 1 fm / hbar
0032     static CELER_CONSTEXPR_FUNCTION InvMomentum fm_par_hbar()
0033     {
0034         return native_value_to<InvMomentum>(units::femtometer
0035                                             / constants::hbar_planck);
0036     }
0037 };
0038
0039 //---------------------------------------------------------------------------//
0040 /*!
0041  * Exponential nuclear form factor.
0042  *
0043  * This nuclear form factor corresponds \c NuclearFormFactorType::exponential
0044  * and assumes the nuclear charge decays exponentially from its center. This
0045  * assumes a parameterization of the atomic nucleus valid for light and medium
0046  * atomic nuclei from Eq. 7 of \citet{butkevich-msc-2002,
0047  * https://doi.org/10.1016/S0168-9002(02)00478-3} : \f[
0048    R_N = 1.27A^{0.27} \,\mathrm{fm}
0049  * \f] with a special case for the proton radius, \f$ R_p = 0.85 \f$ fm.
0050  *
0051  * See \citet{leroy-2016, https://doi-org.ezproxy.cern.ch/10.1142/9167} Eq.
0052  * 2.262.
0053  *
0054  * \todo Instead of using this coarse parameterization, we should add nuclear
0055  * radius to the isotope properties for a more accurate treatment, and
0056  * construct these classes directly from the atomic radius.
0057  */
0058 class ExpNuclearFormFactor : public NuclearFormFactorTraits
0059 {
0060   public:
0061     //! Form factor type corresponding to this distribution
0062     static CELER_CONSTEXPR_FUNCTION FFType ff_type()
0063     {
0064         return FFType::exponential;
0065     }
0066
0067     // Construct with atomic mass number
0068     explicit inline CELER_FUNCTION
0069     ExpNuclearFormFactor(AtomicMassNumber a_mass);
0070
0071     // Construct with precalculated form factor
0072     explicit inline CELER_FUNCTION
0073     ExpNuclearFormFactor(InvMomentumSq prefactor);
0074
0075     // Calculate from square of target momentum
0076     inline CELER_FUNCTION real_type operator()(MomentumSq target_momsq) const;
0077
0078     // Calculate from target momentum
0079     inline CELER_FUNCTION real_type operator()(Momentum target_mom) const;
0080
0081     //! Nuclear form prefactor for the selected isotope
0082     CELER_FORCEINLINE_FUNCTION InvMomentumSq prefactor() const
0083     {
0084         return InvMomentumSq{prefactor_};
0085     }
0086
0087   private:
0088     real_type prefactor_;  // Function of nuclear radius [(MeV/c)^{-2}]
0089 };
0090
0091 //---------------------------------------------------------------------------//
0092 /*!
0093  * Gaussian nuclear form factor.
0094  *
0095  * This nuclear form factor corresponds \c NuclearFormFactorType::gaussian and
0096  * assumes a Gaussian distribution of nuclear charge: see \citet{leroy-2016,
0097  * https://doi-org.ezproxy.cern.ch/10.1142/9167} Eq. 2.264. Its
0098  * prefactor has the same value as the exponential.
0099  */
0100 class GaussianNuclearFormFactor : public ExpNuclearFormFactor
0101 {
0102   public:
0103     //! Form factor type corresponding to this distribution
0104     static CELER_CONSTEXPR_FUNCTION FFType ff_type()
0105     {
0106         return FFType::gaussian;
0107     }
0108
0109     using ExpNuclearFormFactor::ExpNuclearFormFactor;
0110
0111     // Calculate from square of target momentum
0112     inline CELER_FUNCTION real_type operator()(MomentumSq target_momsq) const;
0113
0114     // Calculate from target momentum
0115     inline CELER_FUNCTION real_type operator()(Momentum target_mom) const;
0116 };
0117
0118 //---------------------------------------------------------------------------//
0119 /*!
0120  * Uniform-uniform folded nuclear form factor.
0121  *
0122  * This nuclear form factor corresponds \c NuclearFormFactorType::flat and
0123  * assumes a uniform nuclear charge at the center with a smoothly decreasing
0124  * charge at the surface. From \citet{leroy-2016,
0125  * https://doi-org.ezproxy.cern.ch/10.1142/9167} Eq. 2.265, this leads to a
0126  * form factor: \f[
0127  * F(q) = F'(x(R_0, q)) F'(x(R_1, q))
0128  * \f]
0129  * where \f$ x \equiv q R / \hbar \f$ uses the effective nuclear radius
0130  * \f$ R_0 = 1.2 A^{1/3} \,\mathrm{fm} \f$ and nuclear surface skin
0131  * \f$ R_1 = 2.0 \,\mathrm{fm} \f$,
0132  * and
0133  * \f[
0134  * F'(x) = \frac{3}{x^3} ( \sin x - x \cos x)
0135  * \f]
0136  * is the form factor for a uniformly charged sphere.
0137  *
0138  * \warning This form factor suffers from catastrophic numerical cancellation
0139  * for small radii and momenta so should only be used for large nuclei or large
0140  * momentum transfers.
0141  *
0142  * See also
0143  * \citet{helm-scattering-1956, https://doi.org/10.1103/PhysRev.104.1466},
0144  * \citet{fernandez-msc-1993, https://doi.org/10.1016/0168-583X(93)95827-R}
0145  */
0146 class UUNuclearFormFactor : public NuclearFormFactorTraits
0147 {
0148   public:
0149     //! Form factor type corresponding to this distribution
0150     static CELER_CONSTEXPR_FUNCTION FFType ff_type() { return FFType::flat; }
0151
0152     // Construct with atomic mass number
0153     explicit inline CELER_FUNCTION UUNuclearFormFactor(AtomicMassNumber a_mass);
0154
0155     // Calculate from square of target momentum
0156     inline CELER_FUNCTION real_type operator()(MomentumSq target_momsq) const;
0157
0158     // Calculate from target momentum
0159     inline CELER_FUNCTION real_type operator()(Momentum target_mom) const;
0160
0161   private:
0162     real_type nucl_radius_fm_;
0163
0164     // Effective nuclear skin radius: 2 fm
0165     static CELER_CONSTEXPR_FUNCTION real_type skin_radius_fm() { return 2; }
0166 };
0167
0168 //---------------------------------------------------------------------------//
0169 // INLINE DEFINITIONS
0170 //---------------------------------------------------------------------------//
0171 /*!
0172  * Construct from atomic mass number.
0173  */
0174 CELER_FUNCTION
0175 ExpNuclearFormFactor::ExpNuclearFormFactor(AtomicMassNumber a_mass)
0176 {
0177     CELER_EXPECT(a_mass);
0178     real_type nucl_radius_fm = [a_mass] {
0179         if (CELER_UNLIKELY(a_mass == AtomicMassNumber{1}))
0180         {
0181             // Special case for proton radius
0182             return real_type{0.85};
0183         }
0184         return real_type{1.27}
0185                * fastpow(real_type(a_mass.get()), real_type{0.27});
0186     }();
0187     prefactor_ = ipow<2>(nucl_radius_fm * value_as<InvMomentum>(fm_par_hbar()))
0188                  * (real_type{1} / 12);
0189     CELER_ENSURE(prefactor_ > 0);
0190 }
0191
0192 //---------------------------------------------------------------------------//
0193 /*!
0194  * Construct with precalculated form factor.
0195  */
0196 CELER_FUNCTION
0197 ExpNuclearFormFactor::ExpNuclearFormFactor(InvMomentumSq prefactor)
0198     : prefactor_{prefactor.value()}
0199 {
0200     CELER_EXPECT(prefactor_ > 0);
0201 }
0202
0203 //---------------------------------------------------------------------------//
0204 /*!
0205  * Calculate the exponential folded form factor from the square momentum.
0206  */
0207 CELER_FUNCTION real_type
0208 ExpNuclearFormFactor::operator()(MomentumSq target_momsq) const
0209 {
0210     CELER_EXPECT(target_momsq >= zero_quantity());
0211     return 1 / ipow<2>(1 + prefactor_ * target_momsq.value());
0212 }
0213
0214 //---------------------------------------------------------------------------//
0215 /*!
0216  * Calculate the exponential folded form factor.
0217  */
0218 CELER_FUNCTION real_type ExpNuclearFormFactor::operator()(Momentum target_mom) const
0219 {
0220     return (*this)(MomentumSq{ipow<2>(target_mom.value())});
0221 }
0222
0223 //---------------------------------------------------------------------------//
0224 /*!
0225  * Calculate the gaussian folded form factor.
0226  */
0227 CELER_FUNCTION real_type
0228 GaussianNuclearFormFactor::operator()(MomentumSq target_momsq) const
0229 {
0230     CELER_EXPECT(target_momsq >= zero_quantity());
0231     return std::exp(-2 * value_as<InvMomentumSq>(this->prefactor())
0232                     * target_momsq.value());
0233 }
0234
0235 //---------------------------------------------------------------------------//
0236 /*!
0237  * Calculate the gaussian folded form factor by squaring the momentum.
0238  */
0239 CELER_FUNCTION real_type
0240 GaussianNuclearFormFactor::operator()(Momentum target_mom) const
0241 {
0242     return (*this)(MomentumSq{ipow<2>(target_mom.value())});
0243 }
0244
0245 //---------------------------------------------------------------------------//
0246 /*!
0247  * Construct from atomic mass number.
0248  */
0249 CELER_FUNCTION UUNuclearFormFactor::UUNuclearFormFactor(AtomicMassNumber a_mass)
0250 {
0251     CELER_EXPECT(a_mass);
0252     nucl_radius_fm_ = real_type{1.2}
0253                       * fastpow(real_type(a_mass.get()), real_type{1} / 3);
0254 }
0255
0256 //---------------------------------------------------------------------------//
0257 /*!
0258  * Calculate the uniform-uniform folded form factor by calculating momentum.
0259  */
0260 CELER_FUNCTION real_type
0261 UUNuclearFormFactor::operator()(MomentumSq target_momsq) const
0262 {
0263     CELER_EXPECT(target_momsq >= zero_quantity());
0264     return (*this)(Momentum{std::sqrt(target_momsq.value())});
0265 }
0266
0267 //---------------------------------------------------------------------------//
0268 /*!
0269  * Calculate the uniform-uniform folded form factor.
0270  */
0271 CELER_FUNCTION real_type UUNuclearFormFactor::operator()(Momentum target_mom) const
0272 {
0273     auto sphere_ff = [&target_mom](real_type r) {
0274         // x = q R / hbar
0275         // r units: fm
0276         real_type x = value_as<Momentum>(target_mom)
0277                       * (r * value_as<InvMomentum>(fm_par_hbar()));
0278         return (3 / ipow<3>(x)) * std::fma(-x, std::cos(x), std::sin(x));
0279     };
0280
0281     // Due to catastrophic error for small x, clamp the result to 1
0282     return min(sphere_ff(nucl_radius_fm_)
0283                    * sphere_ff(UUNuclearFormFactor::skin_radius_fm()),
0284                real_type{1});
0285 }
0286
0287 //---------------------------------------------------------------------------//
0288 }  // namespace celeritas