celeritas/mat/IsotopeSelector.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2023-2024 UT-Battelle, LLC, and other Celeritas developers.
0003 // See the top-level COPYRIGHT file for details.
0004 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0005 //---------------------------------------------------------------------------//
0006 //! \file celeritas/mat/IsotopeSelector.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/Assert.hh"
0011 #include "corecel/Macros.hh"
0012 #include "corecel/Types.hh"
0013 #include "celeritas/random/distribution/GenerateCanonical.hh"
0014
0015 #include "ElementView.hh"
0016
0017 namespace celeritas
0018 {
0019 //---------------------------------------------------------------------------//
0020 /*!
0021  * Make a weighted random selection of an isotope.
0022  */
0023 class IsotopeSelector
0024 {
0025   public:
0026     // Construct with element information
0027     inline CELER_FUNCTION IsotopeSelector(ElementView const& element);
0028
0029     // Sample with the given RNG
0030     template<class Engine>
0031     inline CELER_FUNCTION IsotopeComponentId operator()(Engine& rng) const;
0032
0033   private:
0034     ElementView element_;
0035 };
0036
0037 //---------------------------------------------------------------------------//
0038 // INLINE DEFINITIONS
0039 //---------------------------------------------------------------------------//
0040 /*!
0041  * Construct with element.
0042  */
0043 CELER_FUNCTION IsotopeSelector::IsotopeSelector(ElementView const& element)
0044     : element_(element)
0045 {
0046     CELER_EXPECT(element_.num_isotopes());
0047 }
0048
0049 //---------------------------------------------------------------------------//
0050 /*!
0051  * Sample the isotope with the given RNG.
0052  */
0053 template<class Engine>
0054 CELER_FUNCTION IsotopeComponentId IsotopeSelector::operator()(Engine& rng) const
0055 {
0056     auto const& isotopes = element_.isotopes();
0057     real_type cumulative = -generate_canonical(rng);
0058     size_type imax = isotopes.size() - 1;
0059     size_type i = 0;
0060     for (; i < imax; ++i)
0061     {
0062         cumulative += isotopes[i].fraction;
0063         if (cumulative > 0)
0064             break;
0065     }
0066     return IsotopeComponentId{i};
0067 }
0068
0069 //---------------------------------------------------------------------------//
0070 }  // namespace celeritas