[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_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 evaluate(double E, double mu) const;
 
   //! Determine whether angle distribution is empty
   //! \return Whether distribution is empty

include/openmc/distribution_multi.h

@@ -6,6 +6,7 @@
 #include "pugixml.hpp"
 
 #include "openmc/distribution.h"
+#include "openmc/error.h"
 #include "openmc/position.h"
 
 namespace openmc {
@@ -29,6 +30,11 @@ class UnitSphereDistribution {
   //! \return (sampled Direction, sample weight)
   virtual std::pair<Direction, double> sample(uint64_t* seed) const = 0;
 
+  virtual double evaluate(Direction u) const
+  {
+    fatal_error("evaluate not available for this UnitSphereDistribution type");
+  }
+
   Direction u_ref_ {0.0, 0.0, 1.0}; //!< reference direction
 };
 
@@ -52,6 +58,8 @@ class PolarAzimuthal : public UnitSphereDistribution {
   //! \return (sampled Direction, value of the PDF at this Direction)
   std::pair<Direction, double> sample_as_bias(uint64_t* seed) const;
 
+  double evaluate(Direction u) const override;
+
   // Observing pointers
   Distribution* mu() const { return mu_.get(); }
   Distribution* phi() const { return phi_.get(); }
@@ -87,6 +95,8 @@ class Isotropic : public UnitSphereDistribution {
   //! \return (sampled direction, sample weight)
   std::pair<Direction, double> sample(uint64_t* seed) const override;
 
+  double evaluate(Direction u) const override;
+
   // Set or get bias distribution
   void set_bias(std::unique_ptr<PolarAzimuthal> bias)
   {

include/openmc/reaction_product.h

@@ -49,6 +49,11 @@ class ReactionProduct {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(double E_in, double& E_out, double& mu, uint64_t* seed) const;
 
+  AngleEnergy& sample_dist(double E_in, 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,9 @@ 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;
+  Distribution& sample_dist(double E_in, double& E_out, uint64_t* seed) const;
+  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,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 sample_params(double E_in, double& E_out, double& km_a, double& km_r,
+    uint64_t* seed) const;
+  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,9 @@ 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(double E_in, uint64_t* seed) const;
+  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

@@ -6,6 +6,7 @@
 
 #include "openmc/angle_energy.h"
 #include "openmc/endf.h"
+#include "openmc/search.h"
 #include "openmc/secondary_correlated.h"
 #include "openmc/vector.h"
 
@@ -32,6 +33,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 +58,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 +85,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 +113,9 @@ 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;
+  void sample_params(double E_in, double& E_out, int& j, uint64_t* seed) const;
+  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 +144,10 @@ 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;
+  void sample_params(double E_in, double& E_out, double& f, int& l, int& j,
+    uint64_t* seed) const;
+  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 +183,9 @@ 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;
+  const AngleEnergy& sample_dist(double E_in, uint64_t* seed) const;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   CoherentElasticAE coherent_dist_;         //!< Coherent distribution
@@ -178,6 +195,52 @@ class MixedElasticAE : public AngleEnergy {
   const Function1D& incoherent_xs_;      //!< Polymorphic ref. to incoherent XS
 };
 
+struct DoubleVector {
+  double data;
+  const double& operator[](size_t index) const { return data; }
+};
+
+template<typename T>
+double get_pdf_discrete(const vector<double>& mu, const T& w, double mu_0)
+{
+  // Make sure mu is in range [-1,1]
+  if (std::abs(mu_0) > 1.0)
+    mu_0 = std::copysign(1.0, mu_0);
+  double a0;
+  double a1;
+  double b0;
+  double b1;
+  int32_t ai = -1;
+  int32_t bi = -1;
+  if (mu_0 > mu[0]) {
+    ai = lower_bound_index(mu.begin(), mu.end(), mu_0);
+    a0 = mu[ai];
+    a1 = (ai > 1) ? mu[ai - 1] : -1.0;
+  } else {
+    a0 = -1.0;
+    a1 = -1.0;
+  }
+  if (mu_0 < mu[mu.size() - 1]) {
+    bi = upper_bound_index(mu.begin(), mu.end(), mu_0);
+    b0 = mu[bi];
+    b1 = (bi < mu.size() - 1) ? mu[bi + 1] : 1.0;
+  } else {
+    b0 = 1.0;
+    b1 = 1.0;
+  }
+
+  //  Calculate Delta_a and Delta_b
+  double delta_a = 0.5 * std::min(b0 - a0, a0 - a1);
+  double delta_b = 0.5 * std::min(b1 - b0, b0 - a0);
+
+  if (mu_0 < a0 + delta_a)
+    return w[ai] / (2.0 * delta_a);
+  else if (mu_0 + delta_b < b0)
+    return w[bi] / (2.0 * delta_b);
+  else
+    return 0.0;
+}
+
 } // namespace openmc
 
 #endif // OPENMC_SECONDARY_THERMAL_H

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

@@ -51,7 +51,11 @@ class ThermalData {
   //! \param[out] mu Outgoing scattering angle cosine
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(const NuclideMicroXS& micro_xs, double E_in, double* E_out,
-    double* mu, uint64_t* seed);
+    double* mu, uint64_t* seed) const;
+  AngleEnergy& sample_dist(
+    const NuclideMicroXS& micro_xs, double E, uint64_t* seed) const;
+  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).first;
 }
 
+double DecayPhotonAngleEnergy::sample_energy_and_pdf(
+  double E_in, double mu, double& E_out, uint64_t* seed) const
+{
+  E_out = photon_energy_->sample(seed).first;
+  return 0.5;
+}
+
 //==============================================================================
 // Global variables
 //==============================================================================

src/distribution_angle.cpp

@@ -83,4 +83,19 @@ double AngleDistribution::sample(double E, uint64_t* seed) const
   return mu;
 }
 
+double AngleDistribution::evaluate(double E, double mu) const
+{
+  // Find energy bin and calculate interpolation factor
+  int i;
+  double r;
+  get_energy_index(energy_, E, i, r);
+
+  double pdf = 0.0;
+  if (r > 0.0)
+    pdf += r * distribution_[i + 1]->evaluate(mu);
+  if (r < 1.0)
+    pdf += (1.0 - r) * distribution_[i]->evaluate(mu);
+  return pdf;
+}
+
 } // namespace openmc

src/distribution_multi.cpp

@@ -44,6 +44,7 @@ UnitSphereDistribution::UnitSphereDistribution(pugi::xml_node node)
       fatal_error("Angular distribution reference direction must have "
                   "three parameters specified.");
     u_ref_ = Direction(u_ref.data());
+    u_ref_ /= u_ref_.norm();
   }
 }
 
@@ -65,6 +66,7 @@ PolarAzimuthal::PolarAzimuthal(pugi::xml_node node)
       fatal_error("Angular distribution reference v direction must have "
                   "three parameters specified.");
     v_ref_ = Direction(v_ref.data());
+    v_ref_ /= v_ref_.norm();
   }
   w_ref_ = u_ref_.cross(v_ref_);
   if (check_for_node(node, "mu")) {
@@ -116,6 +118,13 @@ std::pair<Direction, double> PolarAzimuthal::sample_impl(
     weight};
 }
 
+double PolarAzimuthal::evaluate(Direction u) const
+{
+  double mu = u.dot(u_ref_);
+  double phi = std::acos(u.dot(v_ref_) / std::sqrt(1 - mu * mu));
+  return mu_->evaluate(mu) * phi_->evaluate(phi);
+}
+
 //==============================================================================
 // Isotropic implementation
 //==============================================================================
@@ -157,6 +166,11 @@ std::pair<Direction, double> Isotropic::sample(uint64_t* seed) const
   }
 }
 
+double Isotropic::evaluate(Direction u) const
+{
+  return 1.0 / (4.0 * PI);
+}
+
 //==============================================================================
 // Monodirectional implementation
 //==============================================================================

src/reaction_product.cpp

@@ -106,8 +106,7 @@ ReactionProduct::ReactionProduct(const ChainNuclide::Product& product)
     make_unique<DecayPhotonAngleEnergy>(chain_nuc->photon_energy()));
 }
 
-void ReactionProduct::sample(
-  double E_in, double& E_out, double& mu, uint64_t* seed) const
+AngleEnergy& ReactionProduct::sample_dist(double E_in, uint64_t* seed) const
 {
   auto n = applicability_.size();
   if (n > 1) {
@@ -118,15 +117,25 @@ void ReactionProduct::sample(
       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);
-        break;
-      }
+      if (c <= prob)
+        return *distribution_[i];
     }
   } else {
     // If only one distribution is present, go ahead and sample it
-    distribution_[0]->sample(E_in, E_out, mu, seed);
+    return *distribution_[0];
   }
 }
 
+void ReactionProduct::sample(
+  double E_in, double& E_out, double& mu, uint64_t* seed) const
+{
+  sample_dist(E_in, seed).sample(E_in, E_out, mu, seed);
+}
+
+double ReactionProduct::sample_energy_and_pdf(
+  double E_in, double mu, double& E_out, uint64_t* seed) const
+{
+  return sample_dist(E_in, seed).sample_energy_and_pdf(E_in, mu, E_out, seed);
+}
+
 } // namespace openmc