executor/detail/MuonicAtomSelector.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/mucf/executor/detail/MuonicAtomSelector.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/cont/EnumArray.hh"
0010 #include "corecel/random/distribution/BernoulliDistribution.hh"
0011 #include "celeritas/Types.hh"
0012 #include "celeritas/mucf/Types.hh"
0013
0014 namespace celeritas
0015 {
0016 namespace detail
0017 {
0018 //---------------------------------------------------------------------------//
0019 /*!
0020  * Select a muonic atom given the mixture of dt in the material.
0021  *
0022  * This class assumes that the material is hydrogen and that the capture \em
0023  * happened to a deuterium or tritium via a simple isotopic fraction selection.
0024  *
0025  * It is needed to \em correct the probability of a deuterium or tritium
0026  * capture, since the isotopic fraction sampling is not sufficient: tritium has
0027  * a higher mass and thus has a biased capture rate.
0028  *
0029  * This effect is calculated using the \f$ q_\text{1S} \f$ formula
0030  * \citet{bom-experimentaldt-2005, https://doi.org/10.1134/1.1926428}
0031  * \f[
0032  * q_\text{1s} = \frac{1}{1 + 2.9 C_t},
0033  * \f]
0034  * where \f$ C_t \f$  is the relative tritium isotope concentration and
0035  * \f$ q_\text{1s} \f$ is the fraction muonic deuterium atoms in the ground
0036  * state. This expression allows calculating the probability of forming a
0037  * muonic deuterium atom via
0038  * \f[
0039  * P_\text{d} = C_d \times q\text{1s}.
0040  * \f]
0041  *
0042  * If a selected uniform random number is \f$ x \leq P_\text{d} \f$, a muonic
0043  * deuterium is formed. Otherwise, a muonic tritium is selected.
0044  */
0045 class MuonicAtomSelector
0046 {
0047   public:
0048     //! Construct with deuterium fraction in the material
0049     inline CELER_FUNCTION MuonicAtomSelector(real_type deuterium_fraction);
0050
0051     // Select muonic atom
0052     template<class Engine>
0053     inline CELER_FUNCTION MucfMuonicAtom operator()(Engine& rng);
0054
0055   private:
0056     real_type deuterium_probability_;
0057 };
0058
0059 //---------------------------------------------------------------------------//
0060 // INLINE DEFINITIONS
0061 //---------------------------------------------------------------------------//
0062 /*!
0063  * Construct with deuterium fraction in the material.
0064  */
0065 CELER_FUNCTION MuonicAtomSelector::MuonicAtomSelector(
0066     real_type deuterium_fraction)
0067 {
0068     using namespace celeritas::literals;
0069
0070     CELER_EXPECT(deuterium_fraction >= 0 && deuterium_fraction <= 1);
0071
0072     real_type tritium_fraction = 1_r - deuterium_fraction;
0073     real_type const q1s = 1_r / (1_r + 2.9_r * tritium_fraction);
0074     deuterium_probability_ = deuterium_fraction * q1s;
0075
0076     CELER_ENSURE(deuterium_probability_ >= 0 && deuterium_probability_ <= 1);
0077 }
0078
0079 //---------------------------------------------------------------------------//
0080 /*!
0081  * Select a muonic atom.
0082  */
0083 template<class Engine>
0084 CELER_FUNCTION MucfMuonicAtom MuonicAtomSelector::operator()(Engine& rng)
0085 {
0086     return BernoulliDistribution(deuterium_probability_)(rng)
0087                ? MucfMuonicAtom::deuterium
0088                : MucfMuonicAtom::tritium;
0089 }
0090
0091 //---------------------------------------------------------------------------//
0092 }  // namespace detail
0093 }  // namespace celeritas