Extend ParticleProductionFilter to support multiple particle types (#3780)

CMakeLists.txt

@@ -427,6 +427,7 @@ list(APPEND libopenmc_SOURCES
   src/tallies/filter_musurface.cpp
   src/tallies/filter_parent_nuclide.cpp
   src/tallies/filter_particle.cpp
+  src/tallies/filter_particle_production.cpp
   src/tallies/filter_polar.cpp
   src/tallies/filter_sph_harm.cpp
   src/tallies/filter_sptl_legendre.cpp

docs/source/io_formats/tallies.rst

@@ -142,9 +142,9 @@ attributes/sub-elements:
 
   :type:
     The type of the filter. Accepted options are "cell", "cellfrom",
-    "cellborn", "surface", "material", "universe", "energy", "energyout", "mu",
-    "polar", "azimuthal", "mesh", "distribcell", "delayedgroup",
-    "energyfunction", and "particle".
+    "cellborn", "surface", "material", "universe", "energy", "energyout",
+    "mu", "polar", "azimuthal", "mesh", "distribcell", "delayedgroup",
+    "energyfunction", "particle", and "particleproduction".
 
   :bins:
      A description of the bins for each type of filter can be found in
@@ -321,6 +321,32 @@ should be set to:
   A list of particle identifiers to tally, specified as strings (e.g.,
   ``neutron``, ``photon``, ``He4``) or as integer PDG numbers.
 
+:particleproduction:
+  This filter tallies secondary particles produced in reactions, binned by
+  particle type and, optionally, by energy. Unlike other energy filters, the
+  weight applied is the weight of the secondary particle. To obtain secondary
+  particle production rates, use this filter with the ``events`` score.
+
+  The filter uses the following sub-elements instead of ``bins``:
+
+  :particles:
+    A space-separated list of secondary particle types to tally (e.g.,
+    ``photon``, ``neutron``, ``electron``).
+
+  :energies:
+    An optional monotonically increasing list of energy boundaries in [eV]
+    for binning the secondary particle energies. If omitted, total production
+    is tallied without energy binning.
+
+  For example, to tally photon and neutron production in three energy groups:
+
+  .. code-block:: xml
+
+      <filter id="1" type="particleproduction">
+        <particles>photon neutron</particles>
+        <energies>0.0 1.0e5 1.0e6 20.0e6</energies>
+      </filter>
+
 ------------------
 ``<mesh>`` Element
 ------------------

include/openmc/tallies/filter_energy.h

@@ -73,36 +73,5 @@ 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

include/openmc/tallies/filter_particle_production.h

@@ -0,0 +1,53 @@
+#ifndef OPENMC_TALLIES_FILTER_PARTICLE_PRODUCTION_H
+#define OPENMC_TALLIES_FILTER_PARTICLE_PRODUCTION_H
+
+#include <unordered_map>
+
+#include "openmc/particle.h"
+#include "openmc/tallies/filter.h"
+#include "openmc/vector.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Bins the outgoing energy of secondary particles
+//!
+//! This filter bins secondary particles by type and, optionally, by energy. It
+//! can be used to get the photon production matrix for multigroup photon
+//! transport, the energy distribution of secondary neutrons, etc. 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 Filter {
+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;
+
+private:
+  //----------------------------------------------------------------------------
+  // Data members
+
+  vector<ParticleType> secondary_types_; //!< Types of secondary particles
+
+  //! Map from PDG number to index in secondary_types_ for O(1) lookup
+  std::unordered_map<int32_t, int> type_to_index_;
+
+  vector<double> energy_bins_; //!< Energy bin boundaries (optional)
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_PARTICLE_PRODUCTION_H

openmc/filter.py

@@ -24,9 +24,9 @@ _FILTER_TYPES = (
     'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy',
     'energyout', 'mu', 'musurface', 'polar', 'azimuthal', 'distribcell', 'delayedgroup',
     'energyfunction', 'cellfrom', 'materialfrom', 'legendre', 'spatiallegendre',
-    'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'cellinstance',
-    'collision', 'time', 'parentnuclide', 'weight', 'meshborn', 'meshsurface',
-    'meshmaterial',
+    'sphericalharmonics', 'zernike', 'zernikeradial', 'particle',
+    'particleproduction', 'cellinstance', 'collision', 'time', 'parentnuclide',
+    'weight', 'meshborn', 'meshsurface', 'meshmaterial',
 )
 
 _CURRENT_NAMES = (
@@ -1748,15 +1748,16 @@ class EnergyoutFilter(EnergyFilter):
     """
 
 
-class ParticleProductionFilter(EnergyFilter):
-    """Bins tally events based on energy of secondary particles.
+class ParticleProductionFilter(Filter):
+    """Bins tally events based on secondary particle type and energy.
 
-    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.
+    This filter bins secondary particles (e.g., photons, electrons, or recoils)
+    produced in a reaction by particle type and, optionally, by energy. 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`.
 
@@ -1764,63 +1765,206 @@ class ParticleProductionFilter(EnergyFilter):
 
     Parameters
     ----------
-    particle : str, int, openmc.ParticleType
-        Type of secondary particle to tally ('photon', 'neutron', etc.)
-    values : Iterable of Real
+    particles : str, int, openmc.ParticleType, or iterable thereof
+        Type(s) of secondary particle(s) to tally ('photon', 'neutron', etc.)
+    energies : Iterable of Real or str, optional
         A list of energy boundaries in [eV]; each successive pair defines a bin.
+        Alternatively, the name of the group structure can be given as a string
+        (must be a key in :data:`openmc.mgxs.GROUP_STRUCTURES`). If not
+        provided, the filter tallies total secondary particle production without
+        energy binning.
     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
+    particles : list of openmc.ParticleType
+        The secondary particle types this filter applies to
+    energies : numpy.ndarray or None
+        Energy boundaries in [eV], or None if no energy binning
+    bins : list
+        A list of bins; each element fully describes one bin. When energies are
+        specified, each element is a tuple ``(particle, energy_low,
+        energy_high)``. When no energies are specified, each element is a
+        particle name string.
     num_bins : int
-        Number of filter bins
-    particle : str
-        The secondary particle type this filter applies to
+        Total number of filter bins
+    num_energy_bins : int
+        Number of energy bins (1 if no energies specified)
+    shape : tuple of int
+        Shape of the filter as (n_particles, n_energy_bins)
     """
 
-    def __init__(self, particle, values, filter_id=None):
-        super().__init__(values, filter_id)
-        self.particle = particle
+    def __init__(self, particles, energies=None, filter_id=None):
+        self.particles = particles
+        self.energies = energies
+        self.id = filter_id
 
     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('\tParticles',
+            [str(p) for p in self.particles])
+        if self.energies is not None:
+            string += '{: <16}=\t{}\n'.format('\tEnergies', self.energies)
         string += '{: <16}=\t{}\n'.format('\tID', self.id)
         return string
 
     @property
-    def particle(self) -> openmc.ParticleType:
-        return self._particle
+    def particles(self):
+        return self._particles
 
-    @particle.setter
-    def particle(self, particle):
-        self._particle = openmc.ParticleType(particle)
+    @particles.setter
+    def particles(self, particles):
+        if isinstance(particles, (str, int, openmc.ParticleType)):
+            self._particles = [openmc.ParticleType(particles)]
+        else:
+            self._particles = [openmc.ParticleType(p) for p in particles]
+
+    @property
+    def energies(self):
+        return self._energies
+
+    @energies.setter
+    def energies(self, energies):
+        if energies is None:
+            self._energies = None
+        elif isinstance(energies, str):
+            cv.check_value('energies', energies,
+                           openmc.mgxs.GROUP_STRUCTURES.keys())
+            self._energies = np.array(
+                openmc.mgxs.GROUP_STRUCTURES[energies.upper()])
+        else:
+            energies = np.asarray(energies, dtype=float)
+            cv.check_length('energies', energies, 2)
+            for i in range(len(energies) - 1):
+                if energies[i + 1] <= energies[i]:
+                    raise ValueError("Energy bins must be monotonically "
+                                     "increasing.")
+            self._energies = energies
+
+    @property
+    def bins(self):
+        if self.energies is None:
+            return [str(p) for p in self.particles]
+        else:
+            result = []
+            energy_pairs = np.vstack(
+                (self.energies[:-1], self.energies[1:])).T
+            for particle in self.particles:
+                for e_low, e_high in energy_pairs:
+                    result.append((str(particle), e_low, e_high))
+            return result
+
+    @bins.setter
+    def bins(self, bins):
+        # bins is set indirectly through particles/energies
+        pass
+
+    def check_bins(self, bins):
+        pass
+
+    @property
+    def num_energy_bins(self):
+        if self.energies is None:
+            return 1
+        else:
+            return len(self.energies) - 1
+
+    @property
+    def num_bins(self):
+        return len(self.particles) * self.num_energy_bins
+
+    @property
+    def shape(self):
+        return (len(self.particles), self.num_energy_bins)
 
     def to_xml_element(self):
-        element = super().to_xml_element()
-        subelement = ET.SubElement(element, 'particle')
-        subelement.text = str(self.particle)
+        element = ET.Element('filter')
+        element.set('id', str(self.id))
+        element.set('type', self.short_name.lower())
+
+        subelement = ET.SubElement(element, 'particles')
+        subelement.text = ' '.join(str(p) for p in self.particles)
+
+        if self.energies is not None:
+            subelement = ET.SubElement(element, 'energies')
+            subelement.text = ' '.join(str(e) for e in self.energies)
+
         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)
+        particles = get_text(elem, 'particles').split()
+
+        bins_elem = elem.find('energies')
+        if bins_elem is not None:
+            energies = [float(x) for x in bins_elem.text.split()]
+        else:
+            energies = None
+
+        return cls(particles, energies=energies, 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)
+
+        # Read particle types
+        particles = [b.decode() if isinstance(b, bytes) else b
+                     for b in group['particles'][()]]
+
+        # Read energy bins if present
+        if 'energies' in group:
+            energies = group['energies'][()]
+        else:
+            energies = None
+
+        return cls(particles, energies=energies, filter_id=filter_id)
+
+    def get_pandas_dataframe(self, data_size, stride, **kwargs):
+        """Builds a Pandas DataFrame for the Filter's bins.
+
+        This method constructs a Pandas DataFrame object for the filter with
+        columns annotated by filter bin information. This is a helper method for
+        :meth:`Tally.get_pandas_dataframe`.
+
+        Parameters
+        ----------
+        data_size : int
+            The total number of bins in the tally corresponding to this filter
+        stride : int
+            Stride in memory for the filter
+
+        Returns
+        -------
+        pandas.DataFrame
+            A Pandas DataFrame with columns for particle type and, if energy
+            bins are specified, energy bin boundaries.
+
+        See also
+        --------
+        Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe()
+
+        """
+        filter_dict = {}
+        key = self.short_name.lower()
+        n_ebins = self.num_energy_bins
+
+        # Particle column — outer dimension (changes slowest)
+        particle_names = [str(p) for p in self.particles]
+        filter_dict[key, 'particle'] = _repeat_and_tile(
+            np.array(particle_names), n_ebins * stride, data_size)
+
+        # Energy columns only if energies were specified
+        if self.energies is not None:
+            energy_pairs = np.vstack(
+                (self.energies[:-1], self.energies[1:])).T
+            filter_dict[key, 'energy low [eV]'] = _repeat_and_tile(
+                energy_pairs[:, 0], stride, data_size)
+            filter_dict[key, 'energy high [eV]'] = _repeat_and_tile(
+                energy_pairs[:, 1], stride, data_size)
+
+        return pd.DataFrame(filter_dict)
 
 
 class TimeFilter(RealFilter):

src/tallies/filter.cpp

@@ -31,6 +31,7 @@
 #include "openmc/tallies/filter_musurface.h"
 #include "openmc/tallies/filter_parent_nuclide.h"
 #include "openmc/tallies/filter_particle.h"
+#include "openmc/tallies/filter_particle_production.h"
 #include "openmc/tallies/filter_polar.h"
 #include "openmc/tallies/filter_sph_harm.h"
 #include "openmc/tallies/filter_sptl_legendre.h"

src/tallies/filter_energy.cpp

@@ -117,50 +117,6 @@ 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/filter_particle_production.cpp

@@ -0,0 +1,108 @@
+#include "openmc/tallies/filter_particle_production.h"
+
+#include <fmt/core.h>
+
+#include "openmc/search.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+//==============================================================================
+// 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++) {
+    const auto& site = p.secondary_bank(bank_idx);
+
+    // Find which particle-type slot this secondary belongs to
+    auto it = type_to_index_.find(site.particle.pdg_number());
+    if (it == type_to_index_.end())
+      continue;
+
+    int particle_idx = it->second;
+    if (energy_bins_.empty()) {
+      // No energy binning, just particle type
+      match.bins_.push_back(particle_idx);
+      match.weights_.push_back(site.wgt);
+    } else {
+      // Bin the energy
+      if (site.E >= energy_bins_.front() && site.E <= energy_bins_.back()) {
+        int n_energies = static_cast<int>(energy_bins_.size()) - 1;
+        auto energy_idx =
+          lower_bound_index(energy_bins_.begin(), energy_bins_.end(), site.E);
+        match.bins_.push_back(particle_idx * n_energies + energy_idx);
+        match.weights_.push_back(site.wgt);
+      }
+    }
+  }
+}
+
+std::string ParticleProductionFilter::text_label(int bin) const
+{
+  if (energy_bins_.empty()) {
+    return fmt::format("Secondary {}", secondary_types_.at(bin).str());
+  } else {
+    int n_energies = static_cast<int>(energy_bins_.size()) - 1;
+    int particle_idx = bin / n_energies;
+    int energy_idx = bin % n_energies;
+    return fmt::format("Secondary {}, Energy [{}, {})",
+      secondary_types_.at(particle_idx).str(), energy_bins_.at(energy_idx),
+      energy_bins_.at(energy_idx + 1));
+  }
+}
+
+void ParticleProductionFilter::from_xml(pugi::xml_node node)
+{
+  // Read energy bins if present (optional)
+  if (check_for_node(node, "energies")) {
+    auto bins = get_node_array<double>(node, "energies");
+    for (int64_t i = 1; i < bins.size(); ++i) {
+      if (bins[i] <= bins[i - 1]) {
+        throw std::runtime_error {
+          "Energy bins must be monotonically increasing."};
+      }
+    }
+    energy_bins_.assign(bins.begin(), bins.end());
+  }
+
+  // Read particle types (required)
+  auto names = get_node_array<std::string>(node, "particles");
+  for (const auto& name : names) {
+    int idx = secondary_types_.size();
+    secondary_types_.emplace_back(name);
+    type_to_index_[secondary_types_.back().pdg_number()] = idx;
+  }
+
+  // Compute total bins
+  if (energy_bins_.empty()) {
+    n_bins_ = secondary_types_.size();
+  } else {
+    n_bins_ = secondary_types_.size() * (energy_bins_.size() - 1);
+  }
+}
+
+void ParticleProductionFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+
+  // Write energy bins if present
+  if (!energy_bins_.empty()) {
+    write_dataset(filter_group, "energies", energy_bins_);
+  }
+
+  // Write particle types
+  vector<std::string> names;
+  for (const auto& pt : secondary_types_) {
+    names.push_back(pt.str());
+  }
+  write_dataset(filter_group, "particles", names);
+}
+
+} // namespace openmc

src/tallies/tally_scoring.cpp

@@ -903,10 +903,9 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_EVENTS:
-// Simply count the number of scoring events
-#pragma omp atomic
-      tally.results_(filter_index, score_index, TallyResult::VALUE) += 1.0;
-      continue;
+      // Simply count the number of scoring events
+      score = 1.0;
+      break;
 
     case ELASTIC:
       if (!p.type().is_neutron())
@@ -1513,10 +1512,9 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_EVENTS:
-// Simply count the number of scoring events
-#pragma omp atomic
-      tally.results_(filter_index, score_index, TallyResult::VALUE) += 1.0;
-      continue;
+      // Simply count the number of scoring events
+      score = 1.0;
+      break;
 
     case ELASTIC:
       if (!p.type().is_neutron())
@@ -2291,10 +2289,9 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_EVENTS:
-// Simply count the number of scoring events
-#pragma omp atomic
-      tally.results_(filter_index, score_index, TallyResult::VALUE) += 1.0;
-      continue;
+      // Simply count the number of scoring events
+      score = 1.0;
+      break;
 
     default:
       continue;

tests/regression_tests/mg_tallies/results_true.dat

@@ -238,8 +238,8 @@ tally 2:
 3.994127E-08
 4.815195E+06
 6.359497E+12
-1.373000E+00
-4.448330E-01
+1.340228E+00
+4.256704E-01
 1.788012E-04
 8.768717E-09
 3.941968E+00
@@ -260,8 +260,8 @@ tally 2:
 1.493588E-07
 1.064740E+07
 2.378109E+13
-2.753000E+00
-1.542179E+00
+2.698926E+00
+1.484990E+00
 3.953669E-04
 3.279028E-08
 8.155605E+00
@@ -282,8 +282,8 @@ tally 2:
 1.319523E-07
 1.005866E+07
 2.100960E+13
-2.829000E+00
-1.702143E+00
+2.795871E+00
+1.662994E+00
 3.735055E-04
 2.896884E-08
 8.365907E+00
@@ -304,8 +304,8 @@ tally 2:
 2.314228E-07
 1.328300E+07
 3.684741E+13
-3.246000E+00
-2.132054E+00
+3.162533E+00
+2.024101E+00
 4.932336E-04
 5.080661E-08
 9.080077E+00
@@ -326,8 +326,8 @@ tally 2:
 1.984988E-07
 1.206905E+07
 3.160522E+13
-3.759000E+00
-2.944479E+00
+3.670108E+00
+2.815133E+00
 4.481564E-04
 4.357848E-08
 1.076370E+01
@@ -348,8 +348,8 @@ tally 2:
 1.474939E-07
 1.003129E+07
 2.348415E+13
-2.577000E+00
-1.509201E+00
+2.523640E+00
+1.451617E+00
 3.724888E-04
 3.238084E-08
 7.654439E+00
@@ -370,8 +370,8 @@ tally 2:
 1.993067E-07
 1.242911E+07
 3.173385E+13
-3.266000E+00
-2.172442E+00
+3.198836E+00
+2.085466E+00
 4.615266E-04
 4.375584E-08
 9.265396E+00
@@ -392,8 +392,8 @@ tally 2:
 2.277785E-07
 1.329858E+07
 3.626717E+13
-3.485000E+00
-2.510947E+00
+3.467276E+00
+2.484377E+00
 4.938123E-04
 5.000655E-08
 9.871966E+00
@@ -414,8 +414,8 @@ tally 2:
 1.293793E-09
 7.073280E+05
 2.059993E+11
-3.000000E-01
-4.514200E-02
+2.978788E-01
+4.498561E-02
 2.626500E-05
 2.840396E-10
 9.197485E-01
@@ -436,8 +436,8 @@ tally 2:
 2.205650E-10
 2.942995E+05
 3.511863E+10
-1.300000E-01
-5.822000E-03
+1.270281E-01
+5.517180E-03
 1.092814E-05
 4.842290E-11
 3.746117E-01
@@ -908,8 +908,8 @@ tally 12:
 3.994127E-08
 4.815195E+06
 6.359497E+12
-1.373000E+00
-4.448330E-01
+1.340228E+00
+4.256704E-01
 1.788012E-04
 8.768717E-09
 2.806910E+00
@@ -928,8 +928,8 @@ tally 12:
 1.493588E-07
 1.064740E+07
 2.378109E+13
-2.753000E+00
-1.542179E+00
+2.698926E+00
+1.484990E+00
 3.953669E-04
 3.279028E-08
 2.869647E+00
@@ -948,8 +948,8 @@ tally 12:
 1.319523E-07
 1.005866E+07
 2.100960E+13
-2.829000E+00
-1.702143E+00
+2.795871E+00
+1.662994E+00
 3.735055E-04
 2.896884E-08
 3.141043E+00
@@ -968,8 +968,8 @@ tally 12:
 2.314228E-07
 1.328300E+07
 3.684741E+13
-3.246000E+00
-2.132054E+00
+3.162533E+00
+2.024101E+00
 4.932336E-04
 5.080661E-08
 3.682383E+00
@@ -988,8 +988,8 @@ tally 12:
 1.984988E-07
 1.206905E+07
 3.160522E+13
-3.759000E+00
-2.944479E+00
+3.670108E+00
+2.815133E+00
 4.481564E-04
 4.357848E-08
 2.634850E+00
@@ -1008,8 +1008,8 @@ tally 12:
 1.474939E-07
 1.003129E+07
 2.348415E+13
-2.577000E+00
-1.509201E+00
+2.523640E+00
+1.451617E+00
 3.724888E-04
 3.238084E-08
 3.192584E+00
@@ -1028,8 +1028,8 @@ tally 12:
 1.993067E-07
 1.242911E+07
 3.173385E+13
-3.266000E+00
-2.172442E+00
+3.198836E+00
+2.085466E+00
 4.615266E-04
 4.375584E-08
 3.402213E+00
@@ -1048,8 +1048,8 @@ tally 12:
 2.277785E-07
 1.329858E+07
 3.626717E+13
-3.485000E+00
-2.510947E+00
+3.467276E+00
+2.484377E+00
 4.938123E-04
 5.000655E-08
 3.110378E-01
@@ -1068,8 +1068,8 @@ tally 12:
 1.293793E-09
 7.073280E+05
 2.059993E+11
-3.000000E-01
-4.514200E-02
+2.978788E-01
+4.498561E-02
 2.626500E-05
 2.840396E-10
 1.267544E-01
@@ -1088,8 +1088,8 @@ tally 12:
 2.205650E-10
 2.942995E+05
 3.511863E+10
-1.300000E-01
-5.822000E-03
+1.270281E-01
+5.517180E-03
 1.092814E-05
 4.842290E-11
 tally 13:

tests/regression_tests/model_xml/photon_production_inputs_true.dat

@@ -48,8 +48,8 @@
       <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>
+      <particles>neutron photon</particles>
+      <energies>0.0 100000.0 300000.0 500000.0 2000000.0 20000000.0</energies>
     </filter>
     <tally id="1">
       <filters>1 2</filters>

tests/regression_tests/photon_production/inputs_true.dat

@@ -48,8 +48,8 @@
       <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>
+      <particles>neutron photon</particles>
+      <energies>0.0 100000.0 300000.0 500000.0 2000000.0 20000000.0</energies>
     </filter>
     <tally id="1">
       <filters>1 2</filters>

tests/regression_tests/photon_production/results_true.dat

@@ -139,6 +139,16 @@ tally 4:
 0.000000E+00
 0.000000E+00
 tally 5:
+5.000000E-04
+2.500000E-07
+1.300000E-03
+1.690000E-06
+8.000000E-04
+6.400000E-07
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
 2.000000E-02
 4.000000E-04
 8.200000E-03

tests/regression_tests/photon_production/test.py

@@ -72,7 +72,7 @@ def model():
 
     # Track source energies of secondary gammas
     ene2_filter = openmc.ParticleProductionFilter(
-        'photon', [0.0, 100e3, 300e3, 500e3, 2e6, 20e6])
+        ['neutron', 'photon'], [0.0, 100e3, 300e3, 500e3, 2e6, 20e6])
 
     neutron_only = openmc.ParticleFilter(['neutron'])
     tally_gam_ene = openmc.Tally()

tests/unit_tests/test_filters.py

@@ -321,12 +321,22 @@ def test_energy_filter():
 
 def test_particle_production_filter():
     energy_bins = [1e3, 1e4, 1e5, 1e6]
+
+    # --- Single particle with energy bins ---
     f = openmc.ParticleProductionFilter('photon', energy_bins)
 
-    assert f.particle == openmc.ParticleType.PHOTON
+    # particles getter always returns a list
+    assert isinstance(f.particles, list)
+    assert len(f.particles) == 1
+    assert f.particles[0] == openmc.ParticleType.PHOTON
+
+    assert f.num_energy_bins == 3
     assert f.num_bins == 3
-    assert f.bins.shape == (3, 2)
-    assert np.allclose(f.values, energy_bins)
+    assert f.shape == (1, 3)
+    assert len(f.bins) == 3
+    # Each bin is (particle_name, e_low, e_high)
+    assert f.bins[0] == ('photon', 1e3, 1e4)
+    assert f.bins[2] == ('photon', 1e5, 1e6)
 
     # __repr__ check
     repr(f)
@@ -335,20 +345,72 @@ def test_particle_production_filter():
     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])
+    assert elem.find('particles').text == 'photon'
+    assert elem.find('energies').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)
+    assert new_f.particles == f.particles
+    assert np.allclose(new_f.energies, f.energies)
 
-    # pandas output
+    # pandas output (with energy bins -> 3 MultiIndex columns)
     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
+    assert ('particleproduction', 'particle') in df.columns
+    assert ('particleproduction', 'energy low [eV]') in df.columns
+    assert ('particleproduction', 'energy high [eV]') in df.columns
+
+    # --- Multiple particles with energy bins ---
+    f2 = openmc.ParticleProductionFilter(['photon', 'neutron'], energy_bins)
+    assert len(f2.particles) == 2
+    assert f2.num_bins == 6  # 2 particles * 3 energy bins
+    assert f2.shape == (2, 3)
+    assert len(f2.bins) == 6
+    # First 3 bins are photon, next 3 are neutron
+    assert f2.bins[0] == ('photon', 1e3, 1e4)
+    assert f2.bins[3] == ('neutron', 1e3, 1e4)
+
+    df2 = f2.get_pandas_dataframe(data_size=6, stride=1)
+    assert df2.shape[0] == 6
+    assert list(df2[('particleproduction', 'particle')]) == \
+        ['photon'] * 3 + ['neutron'] * 3
+
+    # XML round-trip
+    elem2 = f2.to_xml_element()
+    new_f2 = openmc.Filter.from_xml_element(elem2)
+    assert len(new_f2.particles) == 2
+    assert np.allclose(new_f2.energies, energy_bins)
+
+    # --- Multiple particles without energy bins ---
+    f3 = openmc.ParticleProductionFilter(['photon', 'neutron', 'electron'])
+    assert f3.energies is None
+    assert f3.num_bins == 3
+    assert f3.num_energy_bins == 1
+    assert f3.shape == (3, 1)
+    assert f3.bins == ['photon', 'neutron', 'electron']
+
+    repr(f3)
+
+    df3 = f3.get_pandas_dataframe(data_size=3, stride=1)
+    assert df3.shape[0] == 3
+    assert ('particleproduction', 'particle') in df3.columns
+    # Should not have energy columns
+    assert ('particleproduction', 'energy low [eV]') not in df3.columns
+
+    # XML round-trip without energies
+    elem3 = f3.to_xml_element()
+    assert elem3.find('energies') is None
+    new_f3 = openmc.Filter.from_xml_element(elem3)
+    assert new_f3.energies is None
+    assert len(new_f3.particles) == 3
+
+    # --- Energies from group structure name ---
+    f4 = openmc.ParticleProductionFilter('photon', energies='CCFE-709')
+    expected = openmc.mgxs.GROUP_STRUCTURES['CCFE-709']
+    assert np.allclose(f4.energies, expected)
+    assert f4.num_energy_bins == len(expected) - 1
+    assert f4.num_bins == len(expected) - 1
 
 
 def test_weight():