Merge 7ee8cc267b into 61c8a59cff

CMakeLists.txt

@@ -355,6 +355,7 @@ endif()
 #===============================================================================
 
 list(APPEND libopenmc_SOURCES
+  src/angle_energy.cpp
   src/atomic_mass.cpp
   src/bank.cpp
   src/boundary_condition.cpp

include/openmc/angle_energy.h

@@ -27,12 +27,20 @@ public:
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  virtual double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const = 0;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  virtual double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const = 0;
   virtual ~AngleEnergy() = default;
 };
 
+double get_jac_and_transform(
+  double E_in, double& mu, double& E_out, uint64_t* seed, double awr);
+
+double get_jac_and_transform_impl(
+  double E_com, double& mu, double& E_out, uint64_t* seed, double awr);
+
 } // namespace openmc
 
 #endif // OPENMC_ANGLE_ENERGY_H

include/openmc/chain.h

@@ -76,9 +76,11 @@ public:
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const override;
 
 private:
   const Distribution* photon_energy_;

include/openmc/reaction_product.h

@@ -60,9 +60,11 @@ public:
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const;
 
   ParticleType particle_;      //!< Particle type
   EmissionMode emission_mode_; //!< Emission mode

include/openmc/secondary_correlated.h

@@ -53,9 +53,11 @@ public:
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const override;
 
   // energy property
   vector<double>& energy() { return energy_; }

include/openmc/secondary_kalbach.h

@@ -46,9 +46,11 @@ public:
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const override;
 
 private:
   //! Outgoing energy/angle at a single incoming energy

include/openmc/secondary_nbody.h

@@ -39,9 +39,11 @@ public:
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const override;
 
 private:
   int n_bodies_;      //!< Number of particles distributed

include/openmc/secondary_thermal.h

@@ -39,9 +39,11 @@ public:
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const override;
 
 private:
   const CoherentElasticXS& xs_; //!< Coherent elastic scattering cross section
@@ -74,9 +76,11 @@ public:
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const override;
 
 private:
   double debye_waller_;
@@ -108,9 +112,11 @@ public:
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const override;
 
 private:
   const vector<double>& energy_;  //!< Energies at which cosines are tabulated
@@ -149,9 +155,11 @@ public:
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const override;
 
 private:
   const vector<double>& energy_; //!< Incident energies
@@ -196,9 +204,11 @@ public:
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const override;
 
 private:
   //! Secondary energy/angle distribution
@@ -246,9 +256,11 @@ public:
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const override;
 
 private:
   CoherentElasticAE coherent_dist_;         //!< Coherent distribution

include/openmc/secondary_uncorrelated.h

@@ -37,9 +37,11 @@ public:
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const override;
 
   // Accessors
   AngleDistribution& angle() { return angle_; }

include/openmc/thermal.h

@@ -67,9 +67,11 @@ public:
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
   double sample_energy_and_pdf(const NuclideMicroXS& micro_xs, double E_in,
-    double mu, double& E_out, uint64_t* seed) const;
+    double mu, double& E_out, uint64_t* seed, bool is_com, double awr) const;
 
 private:
   struct Reaction {

src/angle_energy.cpp

@@ -0,0 +1,50 @@
+#include "openmc/angle_energy.h"
+
+#include <algorithm> // for clamp
+#include <cmath>     // for sqrt
+
+#include "openmc/constants.h"
+#include "openmc/random_lcg.h"
+
+namespace openmc {
+
+double get_jac_and_transform(
+  double E_in, double& mu, double& E_out, uint64_t* seed, double awr)
+{
+  double E_com = E_in / ((awr + 1.0) * (awr + 1.0));
+  return get_jac_and_transform_impl(E_com, mu, E_out, seed, awr);
+}
+
+double get_jac_and_transform_impl(
+  double E_com, double& mu, double& E_out, uint64_t* seed, double awr)
+{
+  double E_cm = E_out;
+  double mu_lab = mu;
+  double D = mu_lab * mu_lab - 1.0 + E_cm / E_com;
+  if (D <= 0.0)
+    return 0.0;
+  D = std::sqrt(D);
+
+  if ((mu_lab <= 0.0) && (E_cm <= E_com))
+    return 0.0;
+  double E_out1 = E_com * (mu_lab + D) * (mu_lab + D);
+  double mult;
+  if (E_cm > E_com) {
+    mult = 1.0;
+    E_out = E_out1;
+  } else {
+    mult = 2.0;
+    if (prn(seed) < 0.5) {
+      E_out = E_out1;
+    } else {
+      E_out = E_com * (mu_lab - D) * (mu_lab - D);
+    }
+  }
+  mu = mu_lab * std::sqrt(E_out / E_cm) - std::sqrt(E_com / E_cm);
+
+  if ((std::abs(mu) > 1.0) && (std::abs(mu) < 1.0 + FP_PRECISION))
+    mu = std::clamp(mu, -1.0, 1.0);
+  return mult * E_out / (D * std::sqrt(E_cm * E_com));
+}
+
+} // namespace openmc

src/chain.cpp

@@ -74,9 +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
+double DecayPhotonAngleEnergy::sample_energy_and_pdf(double E_in, double mu,
+  double& E_out, uint64_t* seed, bool is_com, double awr) const
 {
+  if (is_com)
+    fatal_error(
+      "DecayPhotonAngleEnergy distribution must be in lab coordinates");
+
   E_out = photon_energy_->sample(seed).first;
   return 0.5;
 }

src/reaction_product.cpp

@@ -134,10 +134,11 @@ void ReactionProduct::sample(
   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
+double ReactionProduct::sample_energy_and_pdf(double E_in, double mu,
+  double& E_out, uint64_t* seed, bool is_com, double awr) const
 {
-  return sample_dist(E_in, seed).sample_energy_and_pdf(E_in, mu, E_out, seed);
+  return sample_dist(E_in, seed)
+    .sample_energy_and_pdf(E_in, mu, E_out, seed, is_com, awr);
 }
 
 } // namespace openmc

src/secondary_correlated.cpp

@@ -260,10 +260,15 @@ void CorrelatedAngleEnergy::sample(
   mu = sample_dist(E_in, E_out, seed).sample(seed).first;
 }
 
-double CorrelatedAngleEnergy::sample_energy_and_pdf(
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+double CorrelatedAngleEnergy::sample_energy_and_pdf(double E_in, double mu,
+  double& E_out, uint64_t* seed, bool is_com, double awr) const
 {
-  return sample_dist(E_in, E_out, seed).evaluate(mu);
+  auto& dist = sample_dist(E_in, E_out, seed);
+  double jac = 1.0;
+  if (is_com)
+    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
+
+  return jac * dist.evaluate(mu);
 }
 
 } // namespace openmc

src/secondary_kalbach.cpp

@@ -232,14 +232,18 @@ void KalbachMann::sample(
     mu = std::log(r1 * std::exp(km_a) + (1.0 - r1) * std::exp(-km_a)) / km_a;
   }
 }
-double KalbachMann::sample_energy_and_pdf(
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+double KalbachMann::sample_energy_and_pdf(double E_in, double mu, double& E_out,
+  uint64_t* seed, bool is_com, double awr) const
 {
   double km_r, km_a;
   sample_params(E_in, E_out, km_a, km_r, seed);
 
+  double jac = 1.0;
+  if (is_com)
+    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
+
   // https://docs.openmc.org/en/latest/methods/neutron_physics.html#equation-KM-pdf-angle
-  return km_a / (2 * std::sinh(km_a)) *
+  return jac * km_a / (2 * std::sinh(km_a)) *
          (std::cosh(km_a * mu) + km_r * std::sinh(km_a * mu));
 }
 

src/secondary_nbody.cpp

@@ -72,11 +72,14 @@ void NBodyPhaseSpace::sample(
   E_out = sample_energy(E_in, seed);
 }
 
-double NBodyPhaseSpace::sample_energy_and_pdf(
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+double NBodyPhaseSpace::sample_energy_and_pdf(double E_in, double mu,
+  double& E_out, uint64_t* seed, bool is_com, double awr) const
 {
   E_out = sample_energy(E_in, seed);
-  return 0.5;
+  double jac = 1.0;
+  if (is_com)
+    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
+  return jac * 0.5;
 }
 
 } // namespace openmc

src/secondary_thermal.cpp

@@ -53,11 +53,16 @@ void CoherentElasticAE::sample(
   mu = 1.0 - 2.0 * energies[k] / E_in;
 }
 
-double CoherentElasticAE::sample_energy_and_pdf(
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+double CoherentElasticAE::sample_energy_and_pdf(double E_in, double mu,
+  double& E_out, uint64_t* seed, bool is_com, double awr) const
 {
   // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7-1)
   E_out = E_in;
+
+  double jac = 1.0;
+  if (is_com)
+    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
+
   const auto& factors = xs_.factors();
 
   if (E_in < bragg_edges_.front())
@@ -68,8 +73,9 @@ double CoherentElasticAE::sample_energy_and_pdf(
   double E = 0.5 * (1 - mu) * E_in;
   double C = 0.5 * E_in / factors[i];
 
-  return C * get_pdf_discrete(bragg_edges_.slice(tensor::range(i + 1)),
-               factors_diff_.slice(tensor::range(i + 1)), E, 0.0, E_in);
+  return jac * C *
+         get_pdf_discrete(bragg_edges_.slice(tensor::range(i + 1)),
+           factors_diff_.slice(tensor::range(i + 1)), E, 0.0, E_in);
 }
 
 //==============================================================================
@@ -90,15 +96,20 @@ void IncoherentElasticAE::sample(
   double c = 2 * E_in * debye_waller_;
   mu = std::log(1.0 + prn(seed) * (std::exp(2.0 * c) - 1)) / c - 1.0;
 }
-double IncoherentElasticAE::sample_energy_and_pdf(
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+
+double IncoherentElasticAE::sample_energy_and_pdf(double E_in, double mu,
+  double& E_out, uint64_t* seed, bool is_com, double awr) const
 {
   E_out = E_in;
 
+  double jac = 1.0;
+  if (is_com)
+    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
+
   // Sample angle by inverting the distribution in ENDF-102, Eq. 7.4
   double c = 2 * E_in * debye_waller_;
   double A = c / (1 - std::exp(-2.0 * c)); // normalization factor
-  return A * std::exp(-c * (1 - mu));
+  return jac * A * std::exp(-c * (1 - mu));
 }
 
 //==============================================================================
@@ -155,8 +166,8 @@ void IncoherentElasticAEDiscrete::sample(
   E_out = E_in;
 }
 
-double IncoherentElasticAEDiscrete::sample_energy_and_pdf(
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+double IncoherentElasticAEDiscrete::sample_energy_and_pdf(double E_in,
+  double mu, double& E_out, uint64_t* seed, bool is_com, double awr) const
 {
   // Get index and interpolation factor for elastic grid
   int i;
@@ -165,8 +176,12 @@ double IncoherentElasticAEDiscrete::sample_energy_and_pdf(
   // Energy doesn't change in elastic scattering
   E_out = E_in;
 
-  return get_pdf_discrete_interpolated(
-    mu_out_.slice(i, tensor::all), mu_out_.slice(i + 1, tensor::all), f, mu);
+  double jac = 1.0;
+  if (is_com)
+    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
+
+  return jac * get_pdf_discrete_interpolated(mu_out_.slice(i, tensor::all),
+                 mu_out_.slice(i + 1, tensor::all), f, mu);
 }
 
 //==============================================================================
@@ -253,8 +268,8 @@ void IncoherentInelasticAEDiscrete::sample(
   mu = (1 - f) * mu_ijk + f * mu_i1jk;
 }
 
-double IncoherentInelasticAEDiscrete::sample_energy_and_pdf(
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+double IncoherentInelasticAEDiscrete::sample_energy_and_pdf(double E_in,
+  double mu, double& E_out, uint64_t* seed, bool is_com, double awr) const
 {
   // Get index and interpolation factor for inelastic grid
   int i;
@@ -263,8 +278,12 @@ double IncoherentInelasticAEDiscrete::sample_energy_and_pdf(
   int j;
   sample_params(E_in, E_out, j, seed);
 
-  return get_pdf_discrete_interpolated(mu_out_.slice(i, j, tensor::all),
-    mu_out_.slice(i + 1, j, tensor::all), f, mu);
+  double jac = 1.0;
+  if (is_com)
+    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
+
+  return jac * get_pdf_discrete_interpolated(mu_out_.slice(i, j, tensor::all),
+                 mu_out_.slice(i + 1, j, tensor::all), f, mu);
 }
 
 //==============================================================================
@@ -403,8 +422,8 @@ void IncoherentInelasticAE::sample(
   mu += std::min(mu - mu_left, mu_right - mu) * (prn(seed) - 0.5);
 }
 
-double IncoherentInelasticAE::sample_energy_and_pdf(
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+double IncoherentInelasticAE::sample_energy_and_pdf(double E_in, double mu,
+  double& E_out, uint64_t* seed, bool is_com, double awr) const
 {
   double f;
   int l, j;
@@ -412,8 +431,12 @@ double IncoherentInelasticAE::sample_energy_and_pdf(
 
   const auto& mu_l = distribution_[l].mu;
 
-  return get_pdf_discrete_interpolated(
-    mu_l.slice(j, tensor::all), mu_l.slice(j + 1, tensor::all), f, mu);
+  double jac = 1.0;
+  if (is_com)
+    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
+
+  return jac * get_pdf_discrete_interpolated(mu_l.slice(j, tensor::all),
+                 mu_l.slice(j + 1, tensor::all), f, mu);
 }
 
 //==============================================================================
@@ -458,10 +481,11 @@ void MixedElasticAE::sample(
   sample_dist(E_in, seed).sample(E_in, E_out, mu, seed);
 }
 
-double MixedElasticAE::sample_energy_and_pdf(
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+double MixedElasticAE::sample_energy_and_pdf(double E_in, double mu,
+  double& E_out, uint64_t* seed, bool is_com, double awr) const
 {
-  return sample_dist(E_in, seed).sample_energy_and_pdf(E_in, mu, E_out, seed);
+  return sample_dist(E_in, seed)
+    .sample_energy_and_pdf(E_in, mu, E_out, seed, is_com, awr);
 }
 
 } // namespace openmc

src/secondary_uncorrelated.cpp

@@ -65,8 +65,8 @@ void UncorrelatedAngleEnergy::sample(
   E_out = energy_->sample(E_in, seed);
 }
 
-double UncorrelatedAngleEnergy::sample_energy_and_pdf(
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+double UncorrelatedAngleEnergy::sample_energy_and_pdf(double E_in, double mu,
+  double& E_out, uint64_t* seed, bool is_com, double awr) const
 {
   // Sample outgoing energy
   if (energy_ != nullptr) {
@@ -75,11 +75,15 @@ double UncorrelatedAngleEnergy::sample_energy_and_pdf(
     E_out = E_in;
   }
 
+  double jac = 1.0;
+  if (is_com)
+    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
+
   if (!angle_.empty()) {
-    return angle_.evaluate(E_in, mu);
+    return jac * angle_.evaluate(E_in, mu);
   } else {
     // no angle distribution given => assume isotropic for all energies
-    return 0.5;
+    return jac * 0.5;
   }
 }
 

src/thermal.cpp

@@ -314,10 +314,11 @@ void ThermalData::sample(const NuclideMicroXS& micro_xs, double E,
 }
 
 double ThermalData::sample_energy_and_pdf(const NuclideMicroXS& micro_xs,
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+  double E_in, double mu, double& E_out, uint64_t* seed, bool is_com,
+  double awr) const
 {
   return sample_dist(micro_xs, E_in, seed)
-    .sample_energy_and_pdf(E_in, mu, E_out, seed);
+    .sample_energy_and_pdf(E_in, mu, E_out, seed, is_com, awr);
 }
 
 void free_memory_thermal()