diff --git a/CMakeLists.txt b/CMakeLists.txt
index ee5d2114e8..c2dff6e100 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -412,6 +412,7 @@ add_library(libopenmc SHARED
   src/position.cpp
   src/pugixml/pugixml_c.cpp
   src/random_lcg.cpp
+  src/reaction_product.cpp
   src/secondary_correlated.cpp
   src/secondary_kalbach.cpp
   src/secondary_nbody.cpp
diff --git a/src/particle.h b/src/particle.h
index 18244098e0..7409a270e8 100644
--- a/src/particle.h
+++ b/src/particle.h
@@ -21,6 +21,11 @@ constexpr int NEUTRON {1};
 constexpr int MAX_LOST_PARTICLES {10};
 constexpr double REL_MAX_LOST_PARTICLES {1.0e-6};
 
+//! Particle types
+enum class ParticleType {
+  neutron, photon, electron, positron
+};
+
 extern "C" {
 
   struct LocalCoord {
diff --git a/src/reaction_product.cpp b/src/reaction_product.cpp
new file mode 100644
index 0000000000..200d9b6bfc
--- /dev/null
+++ b/src/reaction_product.cpp
@@ -0,0 +1,103 @@
+#include "reaction_product.h"
+
+#include <memory> // for unique_ptr
+#include <string> // for string
+
+#include "hdf5_interface.h"
+#include "random_lcg.h"
+#include "secondary_correlated.h"
+#include "secondary_kalbach.h"
+#include "secondary_nbody.h"
+#include "secondary_uncorrelated.h"
+
+namespace openmc {
+
+ReactionProduct::ReactionProduct(hid_t group)
+{
+  // Read particle type
+  std::string temp;
+  read_attribute(group, "particle", temp);
+  if (temp == "neutron") {
+    particle_ = ParticleType::neutron;
+  } else if (temp == "photon") {
+    particle_ = ParticleType::photon;
+  }
+
+  // Read emission mode and decay rate
+  read_attribute(group, "emission_mode", temp);
+  if (temp == "prompt") {
+    emission_mode_ = EmissionMode::prompt;
+  } else if (temp == "delayed") {
+    emission_mode_ = EmissionMode::delayed;
+  } else if (temp == "total") {
+    emission_mode_ = EmissionMode::total;
+  }
+
+  // Read decay rate for delayed emission
+  if (emission_mode_ == EmissionMode::delayed)
+    read_attribute(group, "decay_rate", decay_rate_);
+
+  // Read secondary particle yield
+  hid_t yield = open_dataset(group, "yield");
+  read_attribute(yield, "type", temp);
+  if (temp == "Tabulated1D") {
+    yield_ = std::unique_ptr<Function1D>{new Tabulated1D{yield}};
+  } else if (temp == "Polynomial") {
+    yield_ = std::unique_ptr<Function1D>{new Polynomial{yield}};
+  }
+  close_dataset(yield);
+
+  int n;
+  read_attribute(group, "n_distribution", n);
+
+  for (int i = 0; i < n; ++i) {
+    std::string s {"distribution_"};
+    s.append(std::to_string(i));
+    hid_t dgroup = open_group(group, s.c_str());
+
+    // Read applicability
+    if (n > 1) {
+      hid_t app = open_dataset(dgroup, "applicability");
+      applicability_.emplace_back(app);
+      close_dataset(app);
+    }
+
+    // Determine distribution type and read data
+    read_attribute(dgroup, "type", temp);
+    if (temp == "uncorrelated") {
+      distribution_.emplace_back(new UncorrelatedAngleEnergy{dgroup});
+    } else if (temp == "correlated") {
+      distribution_.emplace_back(new CorrelatedAngleEnergy{dgroup});
+    } else if (temp == "nbody") {
+      distribution_.emplace_back(new NBodyPhaseSpace{dgroup});
+    } else if (temp == "kalbach-mann") {
+      distribution_.emplace_back(new KalbachMann{dgroup});
+    }
+
+    close_group(dgroup);
+  }
+}
+
+void ReactionProduct::sample(double E_in, double& E_out, double& mu) const
+{
+  auto n = applicability_.size();
+  if (n > 1) {
+    double prob = 0.0;
+    double c = prn();
+    for (int i = 0; i < n; ++i) {
+      // Determine probability that i-th energy distribution is sampled
+      prob += applicability_[i](E_in);
+
+      // If i-th distribution is sampled, sample energy from the distribution
+      if (c <= prob) {
+        distribution_[i]->sample(E_in, E_out, mu);
+        break;
+      }
+    }
+  } else {
+    // If only one distribution is present, go ahead and sample it
+    distribution_[0]->sample(E_in, E_out, mu);
+  }
+}
+
+}
\ No newline at end of file
diff --git a/src/reaction_product.h b/src/reaction_product.h
new file mode 100644
index 0000000000..c72657f002
--- /dev/null
+++ b/src/reaction_product.h
@@ -0,0 +1,37 @@
+#ifndef OPENMC_REACTION_PRODUCT_H
+#define OPENMC_REACTION_PRODUCT_H
+
+#include <memory> // for unique_ptr
+#include <vector> // for vector
+
+#include "hdf5.h"
+#include "angle_energy.h"
+#include "endf.h"
+#include "particle.h"
+
+namespace openmc {
+
+class ReactionProduct {
+public:
+  explicit ReactionProduct(hid_t group);
+  void sample(double E_in, double& E_out, double& mu) const;
+private:
+  enum class EmissionMode {
+    prompt,  // Prompt emission of secondary particle
+    total,   // Delayed emission of secondary particle
+    delayed  // Yield represents total emission (prompt + delayed)
+  };
+
+  using Secondary = std::unique_ptr<AngleEnergy>;
+
+  ParticleType particle_;
+  EmissionMode emission_mode_;
+  double decay_rate_;
+  std::unique_ptr<Function1D> yield_;
+  std::vector<Tabulated1D> applicability_;
+  std::vector<Secondary> distribution_;
+};
+
+}
+
+#endif // OPENMC_REACTION_PRODUCT_H