celeritas/random/IsotopeSelector.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/random/IsotopeSelector.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/Assert.hh"
0010 #include "corecel/Macros.hh"
0011 #include "corecel/Types.hh"
0012 #include "corecel/random/distribution/GenerateCanonical.hh"
0013 #include "celeritas/mat/ElementView.hh"
0014
0015 namespace celeritas
0016 {
0017 //---------------------------------------------------------------------------//
0018 /*!
0019  * Make a weighted random selection of an isotope.
0020  */
0021 class IsotopeSelector
0022 {
0023   public:
0024     // Construct with element information
0025     inline CELER_FUNCTION IsotopeSelector(ElementView const& element);
0026
0027     // Sample with the given RNG
0028     template<class Engine>
0029     inline CELER_FUNCTION IsotopeComponentId operator()(Engine& rng) const;
0030
0031   private:
0032     ElementView element_;
0033 };
0034
0035 //---------------------------------------------------------------------------//
0036 // INLINE DEFINITIONS
0037 //---------------------------------------------------------------------------//
0038 /*!
0039  * Construct with element.
0040  */
0041 CELER_FUNCTION IsotopeSelector::IsotopeSelector(ElementView const& element)
0042     : element_(element)
0043 {
0044     CELER_EXPECT(element_.num_isotopes());
0045 }
0046
0047 //---------------------------------------------------------------------------//
0048 /*!
0049  * Sample the isotope with the given RNG.
0050  */
0051 template<class Engine>
0052 CELER_FUNCTION IsotopeComponentId IsotopeSelector::operator()(Engine& rng) const
0053 {
0054     auto const& isotopes = element_.isotopes();
0055     real_type cumulative = -generate_canonical(rng);
0056     size_type imax = isotopes.size() - 1;
0057     size_type i = 0;
0058     for (; i < imax; ++i)
0059     {
0060         cumulative += isotopes[i].fraction;
0061         if (cumulative > 0)
0062             break;
0063     }
0064     return IsotopeComponentId{i};
0065 }
0066
0067 //---------------------------------------------------------------------------//
0068 }  // namespace celeritas