openmc-dev · itay-space · Aug 27, 2025 · Aug 28, 2025 · Aug 28, 2025 · Aug 28, 2025
diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h
@@ -23,7 +23,7 @@ class Distribution {
 public:
   virtual ~Distribution() = default;
   virtual double sample(uint64_t* seed) const = 0;
-
+  virtual double get_pdf(double x) const = 0;
   //! Return integral of distribution
   //! \return Integral of distribution
   virtual double integral() const { return 1.0; };
@@ -84,6 +84,10 @@ class Discrete : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  //! Calculate the probability density function (PDF) at a given value
+  //! \param x The value at which to evaluate the PDF
+  //! \return The value of the PDF at the given point
+  double get_pdf(double x) const;
 
   double integral() const override { return di_.integral(); };
 
@@ -111,6 +115,7 @@ class Uniform : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double get_pdf(double x) const;
 
   double a() const { return a_; }
   double b() const { return b_; }
@@ -135,6 +140,7 @@ class PowerLaw : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double get_pdf(double x) const;
 
   double a() const { return std::pow(offset_, ninv_); }
   double b() const { return std::pow(offset_ + span_, ninv_); }
@@ -160,6 +166,7 @@ class Maxwell : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double get_pdf(double x) const;
 
   double theta() const { return theta_; }
 
@@ -180,6 +187,7 @@ class Watt : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double get_pdf(double x) const;
 
   double a() const { return a_; }
   double b() const { return b_; }
@@ -204,6 +212,7 @@ class Normal : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double get_pdf(double x) const;
 
   double mean_value() const { return mean_value_; }
   double std_dev() const { return std_dev_; }
@@ -227,8 +236,9 @@ class Tabular : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double get_pdf(double x) const;
 
-  // properties
+  //  properties
   vector<double>& x() { return x_; }
   const vector<double>& x() const { return x_; }
   const vector<double>& p() const { return p_; }
@@ -263,6 +273,7 @@ class Equiprobable : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double get_pdf(double x) const;
 
   const vector<double>& x() const { return x_; }
 
@@ -282,6 +293,7 @@ class Mixture : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double get_pdf(double x) const;
 
   double integral() const override { return integral_; }
 

diff --git a/include/openmc/distribution_angle.h b/include/openmc/distribution_angle.h
@@ -25,10 +25,11 @@ class AngleDistribution {
   //! \param[inout] seed pseudorandom number seed pointer
   //! \return Cosine of the angle in the range [-1,1]
   double sample(double E, uint64_t* seed) const;
-
+  double get_pdf(double E, double mu, uint64_t* seed) const;
   //! Determine whether angle distribution is empty
   //! \return Whether distribution is empty
   bool empty() const { return energy_.empty(); }
+  double get_energy(int num) { return energy_[num]; }
 
 private:
   vector<double> energy_;

diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h
@@ -488,6 +488,7 @@ class ParticleData : public GeometryState {
   int event_nuclide_;
   int event_mt_;
   int delayed_group_ {0};
+  int event_index_mt_;
   int parent_nuclide_ {-1};
 
   int n_bank_ {0};
@@ -622,6 +623,7 @@ class ParticleData : public GeometryState {
   const int& event_nuclide() const { return event_nuclide_; }
   int& event_mt() { return event_mt_; } // MT number of collision
   const int& event_mt() const { return event_mt_; }
+  int& event_index_mt() { return event_index_mt_; }
   int& delayed_group() { return delayed_group_; } // delayed group
   const int& parent_nuclide() const { return parent_nuclide_; }
   int& parent_nuclide() { return parent_nuclide_; } // Parent nuclide

diff --git a/include/openmc/reaction_product.h b/include/openmc/reaction_product.h
@@ -11,8 +11,9 @@
 #include "openmc/endf.h"
 #include "openmc/memory.h" // for unique_ptr
 #include "openmc/particle.h"
+#include "openmc/tallies/filter.h"
+#include "openmc/tallies/tally.h"
 #include "openmc/vector.h" // for vector
-
 namespace openmc {
 
 //==============================================================================
@@ -48,7 +49,10 @@ class ReactionProduct {
   //! \param[out] mu Outgoing cosine with respect to current direction
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(double E_in, double& E_out, double& mu, uint64_t* seed) const;
-
+  void get_pdf(int i_tally, double E_in, double& E_out, uint64_t* seed,
+    Particle& p, std::vector<double>& mu_cm, std::vector<double>& Js,
+    std::vector<Particle>& ghost_particles,
+    std::vector<double>& pdfs_lab) const;
   ParticleType particle_;      //!< Particle type
   EmissionMode emission_mode_; //!< Emission mode
   double decay_rate_; //!< Decay rate (for delayed neutron precursors) in [1/s]

diff --git a/include/openmc/secondary_correlated.h b/include/openmc/secondary_correlated.h
@@ -10,6 +10,7 @@
 #include "openmc/angle_energy.h"
 #include "openmc/distribution.h"
 #include "openmc/endf.h"
+#include "openmc/particle.h"
 #include "openmc/vector.h"
 
 namespace openmc {
@@ -40,6 +41,10 @@ class CorrelatedAngleEnergy : public AngleEnergy {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  void get_pdf(double det_pos[4], double E_in, double& E_out, uint64_t* seed,
+    Particle& p, std::vector<double>& mu_cm, std::vector<double>& Js,
+    std::vector<Particle>& ghost_particles,
+    std::vector<double>& pdfs_lab) const;
 
   // energy property
   vector<double>& energy() { return energy_; }

diff --git a/include/openmc/secondary_kalbach.h b/include/openmc/secondary_kalbach.h
@@ -10,6 +10,7 @@
 #include "openmc/angle_energy.h"
 #include "openmc/constants.h"
 #include "openmc/endf.h"
+#include "openmc/particle.h"
 #include "openmc/vector.h"
 
 namespace openmc {
@@ -31,6 +32,10 @@ class KalbachMann : public AngleEnergy {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  void get_pdf(double det_pos[4], double E_in, double& E_out, uint64_t* seed,
+    Particle& p, std::vector<double>& mu_cm, std::vector<double>& Js,
+    std::vector<Particle>& ghost_particles,
+    std::vector<double>& pdfs_lab) const;
 
 private:
   //! Outgoing energy/angle at a single incoming energy

diff --git a/include/openmc/secondary_nbody.h b/include/openmc/secondary_nbody.h
@@ -7,6 +7,7 @@
 #include "hdf5.h"
 
 #include "openmc/angle_energy.h"
+#include "openmc/particle.h"
 
 namespace openmc {
 
@@ -27,6 +28,10 @@ class NBodyPhaseSpace : public AngleEnergy {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  void get_pdf(double det_pos[4], double E_in, double& E_out, uint64_t* seed,
+    Particle& p, std::vector<double>& mu_cm, std::vector<double>& Js,
+    std::vector<Particle>& ghost_particles,
+    std::vector<double>& pdfs_lab) const;
 
 private:
   int n_bodies_;      //!< Number of particles distributed

diff --git a/include/openmc/secondary_thermal.h b/include/openmc/secondary_thermal.h
@@ -32,6 +32,7 @@ class CoherentElasticAE : public AngleEnergy {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double get_pdf(double E_in, double& E_out, double& mu, uint64_t* seed) const;
 
 private:
   const CoherentElasticXS& xs_; //!< Coherent elastic scattering cross section
@@ -55,6 +56,7 @@ class IncoherentElasticAE : public AngleEnergy {
   //! \param[inout] seed Pseudorandom number seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double get_pdf(double E_in, double& E_out, double& mu, uint64_t* seed) const;
 
 private:
   double debye_waller_;
@@ -80,6 +82,7 @@ class IncoherentElasticAEDiscrete : public AngleEnergy {
   //! \param[inout] seed Pseudorandom number seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double get_pdf(double E_in, double& E_out, double& mu, uint64_t* seed) const;
 
 private:
   const vector<double>& energy_;  //!< Energies at which cosines are tabulated
@@ -106,6 +109,8 @@ class IncoherentInelasticAEDiscrete : public AngleEnergy {
   //! \param[inout] seed Pseudorandom number seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double get_pdf(
+    double E_in, double& E_out, double& mu, uint64_t* seed, int l = -1) const;
 
 private:
   const vector<double>& energy_; //!< Incident energies
@@ -134,6 +139,7 @@ class IncoherentInelasticAE : public AngleEnergy {
   //! \param[inout] seed Pseudorandom number seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double get_pdf(double E_in, double& E_out, double& mu, uint64_t* seed) const;
 
 private:
   //! Secondary energy/angle distribution

diff --git a/include/openmc/secondary_uncorrelated.h b/include/openmc/secondary_uncorrelated.h
@@ -10,6 +10,7 @@
 #include "openmc/distribution_angle.h"
 #include "openmc/distribution_energy.h"
 #include "openmc/memory.h"
+#include "openmc/particle.h"
 #include "openmc/vector.h"
 
 namespace openmc {
@@ -31,6 +32,10 @@ class UncorrelatedAngleEnergy : public AngleEnergy {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  void get_pdf(double det_pos[4], double E_in, double& E_out, uint64_t* seed,
+    Particle& p, std::vector<double>& mu_cm, std::vector<double>& Js,
+    std::vector<Particle>& ghost_particles,
+    std::vector<double>& pdfs_lab) const;
 
   // Accessors
   AngleDistribution& angle() { return angle_; }

diff --git a/include/openmc/thermal.h b/include/openmc/thermal.h
@@ -52,6 +52,8 @@ class ThermalData {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(const NuclideMicroXS& micro_xs, double E_in, double* E_out,
     double* mu, uint64_t* seed);
+  double get_pdf(const NuclideMicroXS& micro_xs, double E, double& E_out,
+    double mu, uint64_t* seed);
 
 private:
   struct Reaction {