Implement filter for secondary particle production binned by energy (#3453)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

docs/source/pythonapi/base.rst

@@ -132,6 +132,7 @@ Constructing Tallies
    openmc.MeshSurfaceFilter
    openmc.EnergyFilter
    openmc.EnergyoutFilter
+   openmc.ParticleProductionFilter
    openmc.MuFilter
    openmc.MuSurfaceFilter
    openmc.PolarFilter

include/openmc/constants.h

@@ -301,7 +301,7 @@ enum class TallyEstimator { ANALOG, TRACKLENGTH, COLLISION };
 enum class TallyEvent { SURFACE, LATTICE, KILL, SCATTER, ABSORB };
 
 // Tally score type -- if you change these, make sure you also update the
-// _SCORES dictionary in openmc/capi/tally.py
+// _SCORES dictionary in openmc/lib/tally.py
 //
 // These are kept as a normal enum and made negative, since variables which
 // store one of these enum values usually also may be responsible for storing

include/openmc/particle_data.h

@@ -535,6 +535,11 @@ private:
 
   vector<SourceSite> secondary_bank_;
 
+  // Keep track of how many secondary particles were created in the collision
+  // and what the starting index is in the secondary bank for this particle
+  int n_secondaries_ {0};
+  int secondary_bank_index_ {0};
+
   int64_t current_work_;
 
   vector<double> flux_derivs_;
@@ -689,7 +694,20 @@ public:
 
   // secondary particle bank
   SourceSite& secondary_bank(int i) { return secondary_bank_[i]; }
+  const SourceSite& secondary_bank(int i) const { return secondary_bank_[i]; }
   decltype(secondary_bank_)& secondary_bank() { return secondary_bank_; }
+  decltype(secondary_bank_) const& secondary_bank() const
+  {
+    return secondary_bank_;
+  }
+
+  // Number of secondaries created in a collision
+  int& n_secondaries() { return n_secondaries_; }
+  const int& n_secondaries() const { return n_secondaries_; }
+
+  // Starting index in secondary bank for this collision
+  int& secondary_bank_index() { return secondary_bank_index_; }
+  const int& secondary_bank_index() const { return secondary_bank_index_; }
 
   // Current simulation work index
   int64_t& current_work() { return current_work_; }

include/openmc/tallies/filter.h

@@ -39,6 +39,7 @@ enum class FilterType {
   MUSURFACE,
   PARENT_NUCLIDE,
   PARTICLE,
+  PARTICLE_PRODUCTION,
   POLAR,
   SPHERICAL_HARMONICS,
   SPATIAL_LEGENDRE,

include/openmc/tallies/filter_energy.h

@@ -1,6 +1,7 @@
 #ifndef OPENMC_TALLIES_FILTER_ENERGY_H
 #define OPENMC_TALLIES_FILTER_ENERGY_H
 
+#include "openmc/particle.h"
 #include "openmc/span.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/vector.h"
@@ -72,5 +73,36 @@ public:
   std::string text_label(int bin) const override;
 };
 
+//==============================================================================
+//! Bins the outgoing energy of secondary particles
+//!
+//! This filter can be used to get the photon production matrix for multigroup
+//! photon transport, the energy distribution of secondary neutrons, etc. Unlike
+//! other energy filters, the weight that is applied is equal to the weight of
+//! the secondary particle. Thus, to get secondary production it should be used
+//! in conjunction with the "events" score.
+//==============================================================================
+
+class ParticleProductionFilter : public EnergyFilter {
+public:
+  //----------------------------------------------------------------------------
+  // Methods
+
+  std::string type_str() const override { return "particleproduction"; }
+  FilterType type() const override { return FilterType::PARTICLE_PRODUCTION; }
+
+  void get_all_bins(const Particle& p, TallyEstimator estimator,
+    FilterMatch& match) const override;
+
+  std::string text_label(int bin) const override;
+
+  void from_xml(pugi::xml_node node) override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+protected:
+  ParticleType secondary_type_; //!< Type of secondary particle to filter
+};
+
 } // namespace openmc
 #endif // OPENMC_TALLIES_FILTER_ENERGY_H

openmc/filter.py

@@ -126,9 +126,9 @@ class Filter(IDManagerMixin, metaclass=FilterMeta):
     def __gt__(self, other):
         if type(self) is not type(other):
             if self.short_name in _FILTER_TYPES and \
-                other.short_name in _FILTER_TYPES:
+                    other.short_name in _FILTER_TYPES:
                 delta = _FILTER_TYPES.index(self.short_name) - \
-                        _FILTER_TYPES.index(other.short_name)
+                    _FILTER_TYPES.index(other.short_name)
                 return delta > 0
             else:
                 return False
@@ -275,7 +275,6 @@ class Filter(IDManagerMixin, metaclass=FilterMeta):
             if filter_type == subclass.short_name.lower():
                 return subclass.from_xml_element(elem, **kwargs)
 
-
     def can_merge(self, other):
         """Determine if filter can be merged with another.
 
@@ -426,6 +425,7 @@ class Filter(IDManagerMixin, metaclass=FilterMeta):
 
 class WithIDFilter(Filter):
     """Abstract parent for filters of types with IDs (Cell, Material, etc.)."""
+
     def __init__(self, bins, filter_id=None):
         bins = np.atleast_1d(bins)
 
@@ -629,6 +629,7 @@ class CellInstanceFilter(Filter):
     DistribcellFilter
 
     """
+
     def __init__(self, bins, filter_id=None):
         self.bins = bins
         self.id = filter_id
@@ -749,6 +750,7 @@ class ParticleFilter(Filter):
         The number of filter bins
 
     """
+
     def __eq__(self, other):
         if type(self) is not type(other):
             return False
@@ -1103,6 +1105,7 @@ class MeshMaterialFilter(MeshFilter):
         Unique identifier for the filter
 
     """
+
     def __init__(self, mesh: openmc.MeshBase, bins, filter_id=None):
         self.mesh = mesh
         self.bins = bins
@@ -1437,6 +1440,7 @@ class RealFilter(Filter):
         The number of filter bins
 
     """
+
     def __init__(self, values, filter_id=None):
         self.values = np.asarray(values)
         self.bins = np.vstack((self.values[:-1], self.values[1:])).T
@@ -1712,7 +1716,8 @@ class EnergyFilter(RealFilter):
 
         """
 
-        cv.check_value('group_structure', group_structure, openmc.mgxs.GROUP_STRUCTURES.keys())
+        cv.check_value('group_structure', group_structure,
+                       openmc.mgxs.GROUP_STRUCTURES.keys())
         return cls(openmc.mgxs.GROUP_STRUCTURES[group_structure.upper()])
 
 
@@ -1742,6 +1747,82 @@ class EnergyoutFilter(EnergyFilter):
 
     """
 
+
+class ParticleProductionFilter(EnergyFilter):
+    """Bins tally events based on energy of secondary particles.
+
+    This filter bins the energies of secondary particles (e.g., photons,
+    electrons, or recoils) produced in a reaction. This is useful for
+    constructing production matrices or analyzing secondary particle spectra.
+    Note that unlike other energy filters, the weight that is applied is equal
+    to the weight of the secondary particle. Thus, to obtain secondary particle
+    production, it should be used in conjunction with the "events" score.
+
+    The incident particle type can be filtered using :class:`ParticleFilter`.
+
+    .. versionadded:: 0.15.4
+
+    Parameters
+    ----------
+    particle : str, int, openmc.ParticleType
+        Type of secondary particle to tally ('photon', 'neutron', etc.)
+    values : Iterable of Real
+        A list of energy boundaries in [eV]; each successive pair defines a bin.
+    filter_id : int, optional
+        Unique identifier for the filter
+
+    Attributes
+    ----------
+    values : numpy.ndarray
+        Energy boundaries in [eV]
+    bins : numpy.ndarray
+        Array of (low, high) energy bin pairs
+    num_bins : int
+        Number of filter bins
+    particle : str
+        The secondary particle type this filter applies to
+    """
+
+    def __init__(self, particle, values, filter_id=None):
+        super().__init__(values, filter_id)
+        self.particle = particle
+
+    def __repr__(self):
+        string = type(self).__name__ + '\n'
+        string += '{: <16}=\t{}\n'.format('\tParticle', self.particle)
+        string += '{: <16}=\t{}\n'.format('\tValues', self.values)
+        string += '{: <16}=\t{}\n'.format('\tID', self.id)
+        return string
+
+    @property
+    def particle(self) -> openmc.ParticleType:
+        return self._particle
+
+    @particle.setter
+    def particle(self, particle):
+        self._particle = openmc.ParticleType(particle)
+
+    def to_xml_element(self):
+        element = super().to_xml_element()
+        subelement = ET.SubElement(element, 'particle')
+        subelement.text = str(self.particle)
+        return element
+
+    @classmethod
+    def from_xml_element(cls, elem, **kwargs):
+        filter_id = int(elem.get('id'))
+        values = [float(x) for x in get_text(elem, 'bins').split()]
+        particle = get_text(elem, 'particle')
+        return cls(particle, values, filter_id=filter_id)
+
+    @classmethod
+    def from_hdf5(cls, group, **kwargs):
+        filter_id = int(group.name.split('/')[-1].lstrip('filter '))
+        bins = group['bins'][()]
+        particle = group['particle'][()].decode()
+        return cls(particle, bins, filter_id=filter_id)
+
+
 class TimeFilter(RealFilter):
     """Bins tally events based on the particle's time.
 
@@ -1909,7 +1990,7 @@ class DistribcellFilter(Filter):
         # Make sure there is only 1 bin.
         if not len(bins) == 1:
             msg = (f'Unable to add bins "{bins}" to a DistribcellFilter since '
-                  'only a single distribcell can be used per tally')
+                   'only a single distribcell can be used per tally')
             raise ValueError(msg)
 
         # Check the type and extract the id, if necessary.
@@ -2062,7 +2143,7 @@ class DistribcellFilter(Filter):
         # requests Summary geometric information
         filter_bins = _repeat_and_tile(
             np.arange(self.num_bins), stride, data_size)
-        df = pd.DataFrame({self.short_name.lower() : filter_bins})
+        df = pd.DataFrame({self.short_name.lower(): filter_bins})
 
         # Concatenate with DataFrame of distribcell instance IDs
         if level_df is not None:
@@ -2101,6 +2182,7 @@ class MuFilter(RealFilter):
         The number of filter bins
 
     """
+
     def __init__(self, values, filter_id=None):
         if isinstance(values, Integral):
             values = np.linspace(-1., 1., values + 1)
@@ -2266,6 +2348,7 @@ class DelayedGroupFilter(Filter):
         The number of filter bins
 
     """
+
     def check_bins(self, bins):
         # Check the bin values.
         for g in bins:
@@ -2334,9 +2417,9 @@ class EnergyFunctionFilter(Filter):
     def __gt__(self, other):
         if type(self) is not type(other):
             if self.short_name in _FILTER_TYPES and \
-                other.short_name in _FILTER_TYPES:
+                    other.short_name in _FILTER_TYPES:
                 delta = _FILTER_TYPES.index(self.short_name) - \
-                        _FILTER_TYPES.index(other.short_name)
+                    _FILTER_TYPES.index(other.short_name)
                 return delta > 0
             else:
                 return False
@@ -2346,9 +2429,9 @@ class EnergyFunctionFilter(Filter):
     def __lt__(self, other):
         if type(self) is not type(other):
             if self.short_name in _FILTER_TYPES and \
-                other.short_name in _FILTER_TYPES:
+                    other.short_name in _FILTER_TYPES:
                 delta = _FILTER_TYPES.index(self.short_name) - \
-                        _FILTER_TYPES.index(other.short_name)
+                    _FILTER_TYPES.index(other.short_name)
                 return delta < 0
             else:
                 return False
@@ -2359,14 +2442,16 @@ class EnergyFunctionFilter(Filter):
         string = type(self).__name__ + '\n'
         string += '{: <16}=\t{}\n'.format('\tEnergy', self.energy)
         string += '{: <16}=\t{}\n'.format('\tInterpolant', self.y)
-        string += '{: <16}=\t{}\n'.format('\tInterpolation', self.interpolation)
+        string += '{: <16}=\t{}\n'.format('\tInterpolation',
+                                          self.interpolation)
         return hash(string)
 
     def __repr__(self):
         string = type(self).__name__ + '\n'
         string += '{: <16}=\t{}\n'.format('\tEnergy', self.energy)
         string += '{: <16}=\t{}\n'.format('\tInterpolant', self.y)
-        string += '{: <16}=\t{}\n'.format('\tInterpolation', self.interpolation)
+        string += '{: <16}=\t{}\n'.format('\tInterpolation',
+                                          self.interpolation)
         string += '{: <16}=\t{}\n'.format('\tID', self.id)
         return string
 
@@ -2452,10 +2537,12 @@ class EnergyFunctionFilter(Filter):
     @interpolation.setter
     def interpolation(self, val):
         cv.check_type('interpolation', val, str)
-        cv.check_value('interpolation', val, self.INTERPOLATION_SCHEMES.values())
+        cv.check_value('interpolation', val,
+                       self.INTERPOLATION_SCHEMES.values())
 
         if val == 'quadratic' and len(self.energy) < 3:
-            raise ValueError('Quadratic interpolation requires 3 or more values.')
+            raise ValueError(
+                'Quadratic interpolation requires 3 or more values.')
 
         if val == 'cubic' and len(self.energy) < 4:
             raise ValueError('Cubic interpolation requires 3 or more values.')

src/particle.cpp

@@ -87,6 +87,9 @@ bool Particle::create_secondary(
     return false;
   }
 
+  // Increment number of secondaries created (for ParticleProductionFilter)
+  n_secondaries()++;
+
   auto& bank = secondary_bank().emplace_back();
   bank.particle = type;
   bank.wgt = wgt;
@@ -336,6 +339,7 @@ void Particle::event_cross_surface()
 
 void Particle::event_collide()
 {
+
   // Score collision estimate of keff
   if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) {
     keff_tally_collision() += wgt() * macro_xs().nu_fission / macro_xs().total;
@@ -383,6 +387,11 @@ void Particle::event_collide()
   n_bank() = 0;
   bank_second_E() = 0.0;
   wgt_bank() = 0.0;
+
+  // Clear number of secondaries in this collision. This is
+  // distinct from the number of created neutrons n_bank() above!
+  n_secondaries() = 0;
+
   zero_delayed_bank();
 
   // Reset fission logical

src/physics.cpp

@@ -44,6 +44,7 @@ void collision(Particle& p)
 {
   // Add to collision counter for particle
   ++(p.n_collision());
+  p.secondary_bank_index() = p.secondary_bank().size();
 
   // Sample reaction for the material the particle is in
   switch (p.type().pdg_number()) {
@@ -253,6 +254,7 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
       }
     } else {
       p.secondary_bank().push_back(site);
+      p.n_secondaries()++;
     }
 
     // Increment the number of neutrons born delayed

src/physics_mg.cpp

@@ -27,6 +27,7 @@ void collision_mg(Particle& p)
 {
   // Add to the collision counter for the particle
   p.n_collision()++;
+  p.secondary_bank_index() = p.secondary_bank().size();
 
   // Sample the reaction type
   sample_reaction(p);
@@ -200,6 +201,7 @@ void create_fission_sites(Particle& p)
       }
     } else {
       p.secondary_bank().push_back(site);
+      p.n_secondaries()++;
     }
 
     // Set the delayed group on the particle as well

src/tallies/filter.cpp

@@ -146,6 +146,8 @@ Filter* Filter::create(const std::string& type, int32_t id)
     return Filter::create<ParentNuclideFilter>(id);
   } else if (type == "particle") {
     return Filter::create<ParticleFilter>(id);
+  } else if (type == "particleproduction") {
+    return Filter::create<ParticleProductionFilter>(id);
   } else if (type == "polar") {
     return Filter::create<PolarFilter>(id);
   } else if (type == "surface") {

src/tallies/filter_delayedgroup.cpp

@@ -27,7 +27,7 @@ void DelayedGroupFilter::set_groups(span<int> groups)
     } else if (group > MAX_DELAYED_GROUPS) {
       throw std::invalid_argument {
         "Encountered delayedgroup bin with index " + std::to_string(group) +
-        " which is greater than MAX_DELATED_GROUPS (" +
+        " which is greater than MAX_DELAYED_GROUPS (" +
         std::to_string(MAX_DELAYED_GROUPS) + ")"};
     }
     groups_.push_back(group);
@@ -39,6 +39,10 @@ void DelayedGroupFilter::set_groups(span<int> groups)
 void DelayedGroupFilter::get_all_bins(
   const Particle& p, TallyEstimator estimator, FilterMatch& match) const
 {
+  // Note that the bin is set to zero here, but bins outside zero are
+  // tallied to regardless. This is because that logic has to be handled
+  // in the scoring code instead where looping over the delayed
+  // group takes place (tally_scoring.cpp).
   match.bins_.push_back(0);
   match.weights_.push_back(1.0);
 }

src/tallies/filter_energy.cpp

@@ -117,6 +117,50 @@ std::string EnergyoutFilter::text_label(int bin) const
     "Outgoing Energy [{}, {})", bins_.at(bin), bins_.at(bin + 1));
 }
 
+//==============================================================================
+// ParticleProductionFilter implementation
+//==============================================================================
+
+void ParticleProductionFilter::get_all_bins(
+  const Particle& p, TallyEstimator estimator, FilterMatch& match) const
+{
+  int start_idx = p.secondary_bank_index();
+  int end_idx = start_idx + p.n_secondaries();
+
+  // Loop over secondary bank entries
+  for (int bank_idx = start_idx; bank_idx < end_idx; bank_idx++) {
+    // Check if this is the correct type of secondary, then match its energy if
+    // it's the right type
+    const auto& site = p.secondary_bank(bank_idx);
+    if (site.particle == secondary_type_) {
+      if (site.E >= bins_.front() && site.E <= bins_.back()) {
+        auto bin = lower_bound_index(bins_.begin(), bins_.end(), site.E);
+        match.bins_.push_back(bin);
+        match.weights_.push_back(site.wgt);
+      }
+    }
+  }
+}
+
+std::string ParticleProductionFilter::text_label(int bin) const
+{
+  return fmt::format("Secondary {}, Energy [{}, {})", secondary_type_.str(),
+    bins_.at(bin), bins_.at(bin + 1));
+}
+
+void ParticleProductionFilter::from_xml(pugi::xml_node node)
+{
+  EnergyFilter::from_xml(node);
+  std::string p = get_node_value(node, "particle");
+  secondary_type_ = ParticleType {p};
+}
+
+void ParticleProductionFilter::to_statepoint(hid_t filter_group) const
+{
+  EnergyFilter::to_statepoint(filter_group);
+  write_dataset(filter_group, "particle", secondary_type_.str());
+}
+
 //==============================================================================
 // C-API functions
 //==============================================================================

src/tallies/tally.cpp

@@ -171,6 +171,8 @@ Tally::Tally(pugi::xml_node node)
                filt_type == FilterType::ZERNIKE ||
                filt_type == FilterType::ZERNIKE_RADIAL) {
       estimator_ = TallyEstimator::COLLISION;
+    } else if (filt_type == FilterType::PARTICLE_PRODUCTION) {
+      estimator_ = TallyEstimator::ANALOG;
     }
   }
 

tests/regression_tests/model_xml/photon_production_inputs_true.dat

@@ -41,6 +41,16 @@
     <filter id="2" type="particle">
       <bins>neutron photon electron positron</bins>
     </filter>
+    <filter id="5" type="particle">
+      <bins>neutron</bins>
+    </filter>
+    <filter id="3" type="energy">
+      <bins>0.0 20000000.0</bins>
+    </filter>
+    <filter id="4" type="particleproduction">
+      <bins>0.0 100000.0 300000.0 500000.0 2000000.0 20000000.0</bins>
+      <particle>photon</particle>
+    </filter>
     <tally id="1">
       <filters>1 2</filters>
       <scores>current</scores>
@@ -63,5 +73,10 @@
       <scores>total heating (n,gamma)</scores>
       <estimator>analog</estimator>
     </tally>
+    <tally id="5">
+      <filters>5 3 4</filters>
+      <scores>events</scores>
+      <estimator>analog</estimator>
+    </tally>
   </tallies>
 </model>

tests/regression_tests/photon_production/inputs_true.dat

@@ -41,6 +41,16 @@
     <filter id="2" type="particle">
       <bins>neutron photon electron positron</bins>
     </filter>
+    <filter id="5" type="particle">
+      <bins>neutron</bins>
+    </filter>
+    <filter id="3" type="energy">
+      <bins>0.0 20000000.0</bins>
+    </filter>
+    <filter id="4" type="particleproduction">
+      <bins>0.0 100000.0 300000.0 500000.0 2000000.0 20000000.0</bins>
+      <particle>photon</particle>
+    </filter>
     <tally id="1">
       <filters>1 2</filters>
       <scores>current</scores>
@@ -63,5 +73,10 @@
       <scores>total heating (n,gamma)</scores>
       <estimator>analog</estimator>
     </tally>
+    <tally id="5">
+      <filters>5 3 4</filters>
+      <scores>events</scores>
+      <estimator>analog</estimator>
+    </tally>
   </tallies>
 </model>

tests/regression_tests/photon_production/results_true.dat

@@ -138,3 +138,14 @@ tally 4:
 2.173936E+08
 0.000000E+00
 0.000000E+00
+tally 5:
+2.000000E-02
+4.000000E-04
+8.200000E-03
+6.724000E-05
+5.280000E-02
+2.787840E-03
+3.546000E-01
+1.257412E-01
+5.063000E-01
+2.563397E-01

tests/regression_tests/photon_production/test.py

@@ -25,7 +25,8 @@ def model():
     inner_cyl_right.region = -cyl & +x_plane_center & -x_plane_right
     outer_cyl.region = ~(-cyl & +x_plane_left & -x_plane_right)
     inner_cyl_right.fill = mat
-    model.geometry = openmc.Geometry([inner_cyl_left, inner_cyl_right, outer_cyl])
+    model.geometry = openmc.Geometry(
+        [inner_cyl_left, inner_cyl_right, outer_cyl])
 
     source = openmc.IndependentSource()
     source.space = openmc.stats.Point((0, 0, 0))
@@ -38,17 +39,19 @@ def model():
     model.settings.batches = 1
     model.settings.photon_transport = True
     model.settings.electron_treatment = 'ttb'
-    model.settings.cutoff = {'energy_photon' : 1000.0}
+    model.settings.cutoff = {'energy_photon': 1000.0}
     model.settings.source = source
 
     surface_filter = openmc.SurfaceFilter(cyl)
-    particle_filter = openmc.ParticleFilter(['neutron', 'photon', 'electron', 'positron'])
+    particle_filter = openmc.ParticleFilter(
+        ['neutron', 'photon', 'electron', 'positron'])
     current_tally = openmc.Tally()
     current_tally.filters = [surface_filter, particle_filter]
     current_tally.scores = ['current']
     tally_tracklength = openmc.Tally()
     tally_tracklength.filters = [particle_filter]
-    tally_tracklength.scores = ['total', '(n,gamma)']  # heating doesn't work with tracklength
+    # heating doesn't work with tracklength
+    tally_tracklength.scores = ['total', '(n,gamma)']
     tally_tracklength.nuclides = ['Al27', 'total']
     tally_tracklength.estimator = 'tracklength'
     tally_collision = openmc.Tally()
@@ -61,8 +64,25 @@ def model():
     tally_analog.scores = ['total', 'heating', '(n,gamma)']
     tally_analog.nuclides = ['Al27', 'total']
     tally_analog.estimator = 'analog'
+
+    # This is an analog tally tracking the energy distribution of photons
+    # generated by neutrons. The sum of the tally should give the total
+    # number of photons generated per source neutron.
+    ene_filter = openmc.EnergyFilter([0.0, 20e6])  # incident neutron energy
+
+    # Track source energies of secondary gammas
+    ene2_filter = openmc.ParticleProductionFilter(
+        'photon', [0.0, 100e3, 300e3, 500e3, 2e6, 20e6])
+
+    neutron_only = openmc.ParticleFilter(['neutron'])
+    tally_gam_ene = openmc.Tally()
+    tally_gam_ene.filters = [neutron_only, ene_filter, ene2_filter]
+    tally_gam_ene.scores = ['events']
+    tally_gam_ene.estimator = 'analog'
+
     model.tallies.extend([current_tally, tally_tracklength,
-                          tally_collision, tally_analog])
+                          tally_collision, tally_analog,
+                          tally_gam_ene])
 
     return model
 

tests/unit_tests/test_filters.py

@@ -17,14 +17,16 @@ def box_model():
     model.settings.particles = 100
     model.settings.batches = 10
     model.settings.inactive = 0
-    model.settings.source = openmc.IndependentSource(space=openmc.stats.Point())
+    model.settings.source = openmc.IndependentSource(
+        space=openmc.stats.Point())
     return model
 
 
 def test_cell_instance():
     c1 = openmc.Cell()
     c2 = openmc.Cell()
-    f = openmc.CellInstanceFilter([(c1, 0), (c1, 1), (c1, 2), (c2, 0), (c2, 1)])
+    f = openmc.CellInstanceFilter(
+        [(c1, 0), (c1, 1), (c1, 2), (c2, 0), (c2, 1)])
 
     # Make sure __repr__ works
     repr(f)
@@ -234,7 +236,8 @@ def test_first_moment(run_in_tmpdir, box_model):
         flux, scatter = sp.tallies[plain_tally.id].mean.ravel()
 
         # Check that first moment matches
-        first_score = lambda t: sp.tallies[t.id].mean.ravel()[0]
+        def first_score(t):
+            return sp.tallies[t.id].mean.ravel()[0]
         assert first_score(leg_tally) == scatter
         assert first_score(leg_sptl_tally) == scatter
         assert first_score(sph_scat_tally) == scatter
@@ -258,7 +261,8 @@ def test_lethargy_bin_width():
     assert len(f.lethargy_bin_width) == 175
     energy_bins = openmc.mgxs.GROUP_STRUCTURES['VITAMIN-J-175']
     assert f.lethargy_bin_width[0] == np.log10(energy_bins[1]/energy_bins[0])
-    assert f.lethargy_bin_width[-1] == np.log10(energy_bins[-1]/energy_bins[-2])
+    assert f.lethargy_bin_width[-1] == np.log10(
+        energy_bins[-1]/energy_bins[-2])
 
 
 def test_energyfunc():
@@ -292,7 +296,8 @@ def test_tabular_from_energyfilter():
     # 'histogram' is the default
     assert tab.interpolation == 'histogram'
 
-    tab = efilter.get_tabular(values=np.array([10, 10, 5]), interpolation='linear-linear')
+    tab = efilter.get_tabular(values=np.array(
+        [10, 10, 5]), interpolation='linear-linear')
     assert tab.interpolation == 'linear-linear'
 
 
@@ -314,6 +319,38 @@ def test_energy_filter():
         openmc.EnergyFilter([-1.2, 0.25, 0.5])
 
 
+def test_particle_production_filter():
+    energy_bins = [1e3, 1e4, 1e5, 1e6]
+    f = openmc.ParticleProductionFilter('photon', energy_bins)
+
+    assert f.particle == openmc.ParticleType.PHOTON
+    assert f.num_bins == 3
+    assert f.bins.shape == (3, 2)
+    assert np.allclose(f.values, energy_bins)
+
+    # __repr__ check
+    repr(f)
+
+    # to_xml_element()
+    elem = f.to_xml_element()
+    assert elem.tag == 'filter'
+    assert elem.attrib['type'] == 'particleproduction'
+    assert elem.find('particle').text == 'photon'
+    assert elem.find('bins').text.split()[0] == str(energy_bins[0])
+
+    # from_xml_element()
+    new_f = openmc.Filter.from_xml_element(elem)
+    assert new_f.id == f.id
+    assert new_f.particle == f.particle
+    assert np.allclose(new_f.bins, f.bins)
+
+    # pandas output
+    df = f.get_pandas_dataframe(data_size=3, stride=1)
+    assert df.shape[0] == 3
+    assert "particleproduction low [eV]" in df.columns
+    assert "particleproduction high [eV]" in df.columns
+
+
 def test_weight():
     f = openmc.WeightFilter([0.01, 0.1, 1.0, 10.0])
     expected_bins = [[0.01, 0.1], [0.1, 1.0], [1.0, 10.0]]