Gaussian Energy Broadening

docs/source/methods/tallies.rst

@@ -207,6 +207,42 @@ the change-in-angle), we must use an analog estimator.
 
 .. TODO: Add description of surface current tallies
 
+.. _methods_geb:
+
+--------------------------
+Gaussian Energy Broadening
+--------------------------
+
+Radiation detectors are only able to resolve energies with a finite accuracy,
+which is inherent to each individual detector and depends on the incident energy
+of the particle. A detector's energy resolution :math:`R(E)` is the ratio of the
+Full Width at Half Maximum (FWHM) to the energy of the incident particle,
+:math:`E`:
+
+.. math::
+  R(E) = \frac{\text{FWHM}(E)}{E} = \frac{a + b\sqrt{E + cE^2}}{E}.
+
+The parameters :math:`a`, :math:`b`, and :math:`c` have units of eV,
+eV\ :sup:`1/2`, and eV\ :sup:`-1`, respectively. In general, these parameters
+are detector dependent, and are determined for a particular detector by fitting
+the peaks of known radioactive sources.
+
+To simulate the effect of finite detector resolution in a Monte Carlo
+simulation, we can use a technique known as Gaussian Energy Broadening. With
+this method, the particle's incident energy, :math:`E`, is used to sample a
+new energy, :math:`E'`, which is sampled from a `normal distribution`_ with
+an average of :math:`\mu = E` and a standard deviation of
+
+.. math::
+    \sigma(E) = \frac{\text{FWHM}(E)}{2\sqrt{2\ln{2}}} =
+    \frac{a + b\sqrt{E + cE^2}}{2\sqrt{2\ln{2}}}.
+
+:math:`E'` is only used for tallying, and will therefore only affect tally
+results which have an energy filter or an energy function filter. At each
+collision where a particle contributes to a tally with Gaussian Energy
+Broadening, a new effective energy is sampled for the purposes of tallying. The
+particle is never assigned this sampled energy.
+
 .. _tallies_statistics:
 
 ----------

docs/source/usersguide/tallies.rst

@@ -325,6 +325,21 @@ The following tables show all valid scores:
 
 .. _usersguide_tally_normalization:
 
+--------------------------
+Gaussian Energy Broadening
+--------------------------
+If simulating a radiation detector, :ref:`methods_geb` can be applied to the
+results in the simulation. If the detector parameters :math:`a`, :math:`b`, and
+:math:`c` are known, then Gaussian Energy Broadening can be applied to a tally
+in the following manner::
+
+  tally.gaussian_broadening = (a,b,c)
+
+All scores associated with the tally will then be broadened with these
+parameters. If you would like to simulate multiple detectors, each with
+different broadening parameters, a different :attr:`Tally` instance will be
+needed for each one.
+
 ------------------------------
 Normalization of Tally Results
 ------------------------------

include/openmc/random_dist.h

@@ -50,8 +50,9 @@ extern "C" double watt_spectrum(double a, double b, uint64_t* seed);
 //==============================================================================
 //! Samples an energy from the Gaussian distribution.
 //!
-//! Samples from a normal distribution with a given mean and standard deviation
-//! The PDF is defined as s(x) = (1/2*sigma*sqrt(2) * e-((mu-x)/2*sigma)^2
+//! Samples from a Normal distribution with a given mean and standard deviation
+//! The PDF is defined as
+//! s(x) = (1 / (sigma * sqrt(2 * pi))) * exp(-(1/2) * ((mu-x) / sigma)^2)
 //! Its sampled according to
 //! http://www-pdg.lbl.gov/2009/reviews/rpp2009-rev-monte-carlo-techniques.pdf
 //! section 33.4.4

include/openmc/random_lcg.h

@@ -9,11 +9,12 @@ namespace openmc {
 // Module constants.
 //==============================================================================
 
-constexpr int N_STREAMS {4};
+constexpr int N_STREAMS {5};
 constexpr int STREAM_TRACKING {0};
 constexpr int STREAM_SOURCE {1};
 constexpr int STREAM_URR_PTABLE {2};
 constexpr int STREAM_VOLUME {3};
+constexpr int STREAM_TALLYING {4};
 constexpr int64_t DEFAULT_SEED {1};
 
 //==============================================================================

include/openmc/tallies/tally.h

@@ -178,6 +178,19 @@ class Tally {
 
   int deriv_ {C_NONE}; //!< Index of a TallyDerivative object for diff tallies.
 
+  // This struct contains the information for applying Gaussian Energy
+  // Broadening to a tally. The parameters a, b, and c are used to calculate
+  // the FWHM of the Gaussian distribtuion from which the incident energy for
+  // scoring purposes will be sampled.
+  struct GaussianEnergyBroadening {
+    bool active = false;
+    double a = 0., b = 0., c = 0.;
+
+    void apply(Particle& p) const;
+  };
+
+  GaussianEnergyBroadening gaussian_energy_broadening_;
+
 private:
   //----------------------------------------------------------------------------
   // Private data.

openmc/tallies.py

@@ -56,8 +56,13 @@ class Tally(IDManagerMixin):
         Name of the tally
     multiply_density : bool
         Whether reaction rates should be multiplied by atom density
-
+        
         .. versionadded:: 0.14.0
+    gaussian_broadening : 3-tuple of float
+        Parameters for applying Gaussian energy broadening to the tally. The
+        parameters have units of eV, eV^(1/2), and eV^(-1), respectively.
+
+        .. versionadded:: 0.14.1
     filters : list of openmc.Filter
         List of specified filters for the tally
     nuclides : list of str
@@ -116,6 +121,7 @@ def __init__(self, tally_id=None, name=''):
         self._triggers = cv.CheckedList(openmc.Trigger, 'tally triggers')
         self._derivative = None
         self._multiply_density = True
+        self._gaussian_broadening = None
 
         self._num_realizations = 0
         self._with_summary = False
@@ -164,6 +170,24 @@ def multiply_density(self, value):
         cv.check_type('multiply density', value, bool)
         self._multiply_density = value
 
+    @property
+    def gaussian_broadening(self):
+        return self._gaussian_broadening
+
+    @gaussian_broadening.setter
+    def gaussian_broadening(self, geb):
+        if not isinstance(geb, tuple):
+            raise TypeError(f'Gaussian broadening on Tally must be a '
+                            f'3-tuple of floats.')
+        elif len(geb) != 3:
+            raise TypeError(f'Gaussian broadening on Tally must be a '
+                            f'3-tuple of floats.')
+        for p in geb:
+            if p < 0.:
+                raise ValueError(f'All Gaussian broadening parameters '
+                                 f'on a Tally must be >= 0.')
+        self._gaussian_broadening = geb
+
     @property
     def filters(self):
         return self._filters
@@ -848,6 +872,11 @@ def to_xml_element(self):
         if not self.multiply_density:
             element.set("multiply_density", str(self.multiply_density).lower())
 
+        # Gaussian Energy Broadenning
+        if self.gaussian_broadening is not None:
+            subelement = ET.SubElement(element, "gaussian-energy-broadening")
+            subelement.text = ' '.join(str(p) for p in self.gaussian_broadening)
+
         # Optional Tally filters
         if len(self.filters) > 0:
             subelement = ET.SubElement(element, "filters")
@@ -941,6 +970,11 @@ def from_xml_element(cls, elem, **kwargs):
             deriv_id = int(deriv_elem.text)
             tally.derivative = kwargs['derivatives'][deriv_id]
 
+        # Read Gaussian Energy Broadening
+        geb_elem = elem.find('gaussian-energy-broadening')
+        if geb_elem is not None:
+            tally.gaussian_broadening = (float(p) for p in geb_elem.text.split())
+
         return tally
 
     def contains_filter(self, filter_type):

src/tallies/tally.cpp

@@ -11,6 +11,7 @@
 #include "openmc/mgxs_interface.h"
 #include "openmc/nuclide.h"
 #include "openmc/particle.h"
+#include "openmc/random_dist.h"
 #include "openmc/reaction.h"
 #include "openmc/reaction_product.h"
 #include "openmc/settings.h"
@@ -317,6 +318,47 @@ Tally::Tally(pugi::xml_node node)
     }
   }
 
+  // =======================================================================
+  // SET GAUSSIAN ENERGY BROADENING
+  if (check_for_node(node, "gaussian-energy-broadening")) {
+    if (settings::run_CE == false) {
+      // Cannot use Gaussian broadening in a MG problem
+      throw std::runtime_error {
+        fmt::format("Cannot use Gaussian energy broadening in a multi-group "
+                    "problem. Found on tally {}",
+          id_)};
+    }
+
+    auto params = get_node_array<double>(node, "gaussian-energy-broadening");
+
+    if (params.size() != 3) {
+      throw std::runtime_error {
+        fmt::format("Gaussian energy broadening was given {} parameters "
+                    "instead of 3 on tally {}",
+          params.size(), id_)};
+    }
+
+    gaussian_energy_broadening_.active = true;
+    gaussian_energy_broadening_.a = params[0];
+    gaussian_energy_broadening_.b = params[1];
+    gaussian_energy_broadening_.c = params[2];
+
+    if (gaussian_energy_broadening_.a < 0.) {
+      throw std::runtime_error {fmt::format(
+        "Gaussian energy broadening parameter a is < 0 on tally {}", id_)};
+    }
+
+    if (gaussian_energy_broadening_.b < 0.) {
+      throw std::runtime_error {fmt::format(
+        "Gaussian energy broadening parameter b is < 0 on tally {}", id_)};
+    }
+
+    if (gaussian_energy_broadening_.c < 0.) {
+      throw std::runtime_error {fmt::format(
+        "Gaussian energy broadening parameter c is < 0 on tally {}", id_)};
+    }
+  }
+
 #ifdef LIBMESH
   // ensure a tracklength tally isn't used with a libMesh filter
   for (auto i : this->filters_) {
@@ -834,6 +876,46 @@ std::string Tally::nuclide_name(int nuclide_idx) const
   return data::nuclides.at(nuclide)->name_;
 }
 
+//==============================================================================
+// GaussianEnergyBroadening object implementation
+//==============================================================================
+
+void Tally::GaussianEnergyBroadening::apply(Particle& p) const
+{
+  // This method applies GEB to a particle, changing E_last
+
+  if (active == false)
+    return;
+
+  // If the particle has a zero energy or (more likely) its contribution to a
+  // pulse height tally was 0, we don't perform GEB
+  if (p.E_last() <= 0.)
+    return;
+
+  // Calculate the FWHM
+  const double FWHM =
+    a + b * std::sqrt(p.E_last() + c * p.E_last() * p.E_last());
+
+  // Calculate sigma of the gaussian
+  constexpr double sigma_coeff = 1. / (2. * std::sqrt(2. * std::log(2.)));
+  const double sigma = sigma_coeff * FWHM;
+
+  // Save a copy of the original RNG stream, and set the stream for GEB
+  const int orig_stream = p.stream();
+  p.stream() = STREAM_TALLYING;
+
+  // Sample new energy
+  const double E = normal_variate(p.E_last(), sigma, p.current_seed());
+
+  // Reset the original RNG stream of the particle
+  p.stream() = orig_stream;
+
+  // Set E_last to the sampled energy. We only do this however if we
+  // sampled a valid energy.
+  if (E > 0.)
+    p.E_last() = E;
+}
+
 //==============================================================================
 // Non-member functions
 //==============================================================================