[Point Detector] Add distribution get_pdf functionality

include/openmc/angle_energy.h

@@ -16,6 +16,8 @@ class AngleEnergy {
 public:
   virtual void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const = 0;
+  virtual double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const = 0;
   virtual ~AngleEnergy() = default;
 };
 

include/openmc/chain.h

@@ -70,6 +70,8 @@ class DecayPhotonAngleEnergy : public AngleEnergy {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   const Distribution* photon_energy_;

include/openmc/distribution.h

@@ -9,6 +9,7 @@
 #include "pugixml.hpp"
 
 #include "openmc/constants.h"
+#include "openmc/error.h"
 #include "openmc/memory.h" // for unique_ptr
 #include "openmc/span.h"
 #include "openmc/vector.h" // for vector
@@ -23,6 +24,10 @@ class Distribution {
 public:
   virtual ~Distribution() = default;
   virtual double sample(uint64_t* seed) const = 0;
+  virtual double pdf(double x) const
+  {
+    fatal_error("pdf not available for this Distribution type");
+  }
 
   //! Return integral of distribution
   //! \return Integral of distribution
@@ -111,6 +116,7 @@ class Uniform : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double pdf(double x) const override;
 
   double a() const { return a_; }
   double b() const { return b_; }
@@ -135,6 +141,7 @@ class PowerLaw : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double pdf(double x) const override;
 
   double a() const { return std::pow(offset_, ninv_); }
   double b() const { return std::pow(offset_ + span_, ninv_); }
@@ -204,6 +211,7 @@ class Normal : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double pdf(double x) const override;
 
   double mean_value() const { return mean_value_; }
   double std_dev() const { return std_dev_; }
@@ -227,6 +235,7 @@ class Tabular : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double pdf(double x) const override;
 
   // properties
   vector<double>& x() { return x_; }

include/openmc/distribution_angle.h

@@ -25,6 +25,7 @@ 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 pdf(double E, double mu) const;
 
   //! Determine whether angle distribution is empty
   //! \return Whether distribution is empty

include/openmc/reaction_product.h

@@ -49,6 +49,9 @@ class ReactionProduct {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(double E_in, double& E_out, double& mu, uint64_t* seed) const;
 
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const;
+
   ParticleType particle_;      //!< Particle type
   EmissionMode emission_mode_; //!< Emission mode
   double decay_rate_; //!< Decay rate (for delayed neutron precursors) in [1/s]

include/openmc/secondary_correlated.h

@@ -40,6 +40,8 @@ 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;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
   // energy property
   vector<double>& energy() { return energy_; }

include/openmc/secondary_kalbach.h

@@ -31,6 +31,8 @@ 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;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   //! Outgoing energy/angle at a single incoming energy

include/openmc/secondary_nbody.h

@@ -27,6 +27,8 @@ 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;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   int n_bodies_;      //!< Number of particles distributed

include/openmc/secondary_thermal.h

@@ -32,6 +32,8 @@ 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 sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   const CoherentElasticXS& xs_; //!< Coherent elastic scattering cross section
@@ -55,6 +57,8 @@ 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 sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   double debye_waller_;
@@ -80,6 +84,8 @@ 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 sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   const vector<double>& energy_;  //!< Energies at which cosines are tabulated
@@ -106,6 +112,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 sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   const vector<double>& energy_; //!< Incident energies
@@ -134,6 +142,8 @@ 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 sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   //! Secondary energy/angle distribution
@@ -169,6 +179,8 @@ class MixedElasticAE : 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 sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   CoherentElasticAE coherent_dist_;         //!< Coherent distribution

include/openmc/secondary_uncorrelated.h

@@ -31,6 +31,8 @@ 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;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
   // Accessors
   AngleDistribution& angle() { return angle_; }

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 sample_energy_and_pdf(const NuclideMicroXS& micro_xs, double E_in,
+    double mu, double& E_out, uint64_t* seed) const;
 
 private:
   struct Reaction {

src/chain.cpp

@@ -74,6 +74,13 @@ void DecayPhotonAngleEnergy::sample(
   mu = Uniform(-1., 1.).sample(seed);
 }
 
+double DecayPhotonAngleEnergy::sample_energy_and_pdf(
+  double E_in, double mu, double& E_out, uint64_t* seed) const
+{
+  E_out = photon_energy_->sample(seed);
+  return 0.5;
+}
+
 //==============================================================================
 // Global variables
 //==============================================================================

src/distribution.cpp

@@ -12,6 +12,7 @@
 #include "openmc/math_functions.h"
 #include "openmc/random_dist.h"
 #include "openmc/random_lcg.h"
+#include "openmc/search.h"
 #include "openmc/xml_interface.h"
 
 namespace openmc {
@@ -156,6 +157,14 @@ double Uniform::sample(uint64_t* seed) const
   return a_ + prn(seed) * (b_ - a_);
 }
 
+double Uniform::pdf(double x) const
+{
+  if (x <= b_ && x >= a_)
+    return 1 / (b_ - a_);
+  else
+    return 0;
+}
+
 //==============================================================================
 // PowerLaw implementation
 //==============================================================================
@@ -182,6 +191,21 @@ double PowerLaw::sample(uint64_t* seed) const
   return std::pow(offset_ + prn(seed) * span_, ninv_);
 }
 
+double PowerLaw::pdf(double x) const
+{
+  // Use accessors
+  double a_val = this->a();
+  double b_val = this->b();
+  double n_val = this->n();
+
+  if (x < a_val || x > b_val) {
+    return 0.0; // outside support
+  }
+
+  return (n_val + 1.0) * std::pow(x, n_val) /
+         (std::pow(b_val, n_val + 1.0) - std::pow(a_val, n_val + 1.0));
+}
+
 //==============================================================================
 // Maxwell implementation
 //==============================================================================
@@ -236,6 +260,13 @@ double Normal::sample(uint64_t* seed) const
   return normal_variate(mean_value_, std_dev_, seed);
 }
 
+double Normal::pdf(double x) const
+{
+  double exponent = -0.5 * std::pow((x - mean_value_) / std_dev_, 2);
+  double coefficient = 1 / (std_dev_ * std::sqrt(2 * PI));
+  return coefficient * std::exp(exponent);
+}
+
 //==============================================================================
 // Tabular implementation
 //==============================================================================
@@ -356,6 +387,43 @@ double Tabular::sample(uint64_t* seed) const
   }
 }
 
+double Tabular::pdf(double x) const
+{
+  // get PDF value at x
+
+  int i;
+  std::size_t n = x_.size();
+  if (x < x_[0]) {
+    return 0;
+  } else if (x > x_[n - 1]) {
+    return 0;
+  } else {
+    i = lower_bound_index(x_.begin(), x_.end(), x);
+  }
+
+  // Determine bounding PDF values
+  double x_i = x_[i];
+  double p_i = p_[i];
+
+  switch (interp_) {
+  case Interpolation::histogram:
+    // Histogram interpolation
+    return p_i;
+
+  case Interpolation::lin_lin: {
+    // Linear-linear interpolation
+    double x_i1 = x_[i + 1];
+    double p_i1 = p_[i + 1];
+
+    double m = (p_i1 - p_i) / (x_i1 - x_i);
+    return (m == 0.0) ? p_i : p_i + (x - x_i) * m;
+  }
+
+  default:
+    fatal_error("Unsupported interpolation type in PDF evaluation.");
+  }
+}
+
 //==============================================================================
 // Equiprobable implementation
 //==============================================================================

src/distribution_angle.cpp

@@ -95,4 +95,32 @@ double AngleDistribution::sample(double E, uint64_t* seed) const
   return mu;
 }
 
+double AngleDistribution::pdf(double E, double mu) const
+{
+  // Determine number of incoming energies
+  auto n = energy_.size();
+
+  // Find energy bin and calculate interpolation factor -- if the energy is
+  // outside the range of the tabulated energies, choose the first or last bins
+  int i;
+  double r;
+  if (E < energy_[0]) {
+    i = 0;
+    r = 0.0;
+  } else if (E > energy_[n - 1]) {
+    i = n - 2;
+    r = 1.0;
+  } else {
+    i = lower_bound_index(energy_.begin(), energy_.end(), E);
+    r = (E - energy_[i]) / (energy_[i + 1] - energy_[i]);
+  }
+
+  double pdf = 0.0;
+  if (r > 0.0)
+    pdf += r * distribution_[i + 1]->pdf(mu);
+  if (r < 1.0)
+    pdf += (1.0 - r) * distribution_[i]->pdf(mu);
+  return pdf;
+}
+
 } // namespace openmc

src/reaction_product.cpp

@@ -129,4 +129,33 @@ void ReactionProduct::sample(
   }
 }
 
+double ReactionProduct::sample_energy_and_pdf(
+  double E_in, double mu, double& E_out, uint64_t* seed) const
+{
+
+  int distribution_index;
+  auto n = applicability_.size();
+  if (n > 1) {
+    double prob = 0.0;
+    double c = prn(seed);
+    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, seed);
+        distribution_index = i;
+        break;
+      }
+    }
+  } else {
+    distribution_index = 0;
+  }
+  // now extract pdf
+
+  AngleEnergy* angleEnergyPtr = distribution_[distribution_index].get();
+  return angleEnergyPtr->sample_energy_and_pdf(E_in, mu, E_out, seed);
+}
+
 } // namespace openmc