Tally spectrum of secondary particles

include/openmc/particle_data.h

@@ -469,6 +469,7 @@ class ParticleData : public GeometryState {
   CacheDataMG mg_xs_cache_;
 
   ParticleType type_ {ParticleType::neutron};
+  ParticleType type_last_ {ParticleType::neutron};
 
   double E_;
   double E_last_;
@@ -566,6 +567,8 @@ class ParticleData : public GeometryState {
   // Particle type (n, p, e, gamma, etc)
   ParticleType& type() { return type_; }
   const ParticleType& type() const { return type_; }
+  ParticleType& type_last() { return type_last_; }
+  const ParticleType& type_last() const { return type_last_; }
 
   // Current particle energy, energy before collision,
   // and corresponding multigroup group indices. Energy

include/openmc/tallies/filter.h

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

include/openmc/tallies/filter_particle.h

@@ -48,5 +48,45 @@ class ParticleFilter : public Filter {
   vector<ParticleType> particles_;
 };
 
+//==============================================================================
+//! Bins by type of outgoing particle (e.g. neutron, photon).
+//==============================================================================
+
+class ParticleOutFilter : public Filter {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors, destructors
+
+  ~ParticleOutFilter() = default;
+
+  //----------------------------------------------------------------------------
+  // Methods
+
+  std::string type_str() const override { return "particleout"; }
+  FilterType type() const override { return FilterType::PARTICLE_OUT; }
+
+  void from_xml(pugi::xml_node node) override;
+
+  void get_all_bins(const Particle& p, TallyEstimator estimator,
+    FilterMatch& match) const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  //----------------------------------------------------------------------------
+  // Accessors
+
+  const vector<ParticleType>& particles() const { return particles_; }
+
+  void set_particles(span<ParticleType> particles);
+
+private:
+  //----------------------------------------------------------------------------
+  // Data members
+
+  vector<ParticleType> particles_;
+};
+
 } // namespace openmc
 #endif // OPENMC_TALLIES_FILTER_PARTICLE_H

include/openmc/tallies/tally.h

@@ -202,6 +202,7 @@ extern std::unordered_map<int, int> tally_map;
 extern vector<unique_ptr<Tally>> tallies;
 extern vector<int> active_tallies;
 extern vector<int> active_analog_tallies;
+extern vector<int> active_particleout_analog_tallies;
 extern vector<int> active_tracklength_tallies;
 extern vector<int> active_timed_tracklength_tallies;
 extern vector<int> active_collision_tallies;

include/openmc/tallies/tally_scoring.h

@@ -72,14 +72,16 @@ void score_collision_tally(Particle& p);
 //! Analog tallies are triggered at every collision, not every event.
 //
 //! \param p The particle being tracked
-void score_analog_tally_ce(Particle& p);
+//! \param tallies A vector of the indices of the tallies to score to
+void score_analog_tally_ce(Particle& p, const vector<int>& tallies);
 
 //! Score tallies based on a simple count of events (for multigroup).
 //
 //! Analog tallies are triggered at every collision, not every event.
 //
 //! \param p The particle being tracked
-void score_analog_tally_mg(Particle& p);
+//! \param tallies A vector of the indices of the tallies to score to
+void score_analog_tally_mg(Particle& p, const vector<int>& tallies);
 
 //! Score tallies using a tracklength estimate of the flux.
 //

openmc/filter.py

@@ -24,9 +24,9 @@
     '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', 'particleout', 
+    'cellinstance', 'collision', 'time', 'parentnuclide', 'weight', 'meshborn',
+    'meshsurface', 'meshmaterial',
 )
 
 _CURRENT_NAMES = (
@@ -785,6 +785,29 @@ def from_xml_element(cls, elem, **kwargs):
         filter_id = int(get_text(elem, "id"))
         bins = get_elem_list(elem, "bins", str) or []
         return cls(bins, filter_id=filter_id)
+
+
+class ParticleoutFilter(ParticleFilter):
+    """Bins tally events based on the outgoing particle type.
+
+    Parameters
+    ----------
+    bins : str, or sequence of str
+        The particles to tally represented as strings ('neutron', 'photon',
+        'electron', 'positron').
+    filter_id : int
+        Unique identifier for the filter
+
+    Attributes
+    ----------
+    bins : sequence of str
+        The particles to tally
+    id : int
+        Unique identifier for the filter
+    num_bins : Integral
+        The number of filter bins
+
+    """        
 
 
 class ParentNuclideFilter(ParticleFilter):

openmc/lib/filter.py

@@ -22,9 +22,9 @@
     'EnergyFilter', 'EnergyoutFilter', 'EnergyFunctionFilter', 'LegendreFilter',
     'MaterialFilter', 'MaterialFromFilter', 'MeshFilter', 'MeshBornFilter',
     'MeshMaterialFilter', 'MeshSurfaceFilter', 'MuFilter', 'MuSurfaceFilter',
-    'ParentNuclideFilter', 'ParticleFilter', 'PolarFilter', 'SphericalHarmonicsFilter',
-    'SpatialLegendreFilter', 'SurfaceFilter', 'TimeFilter', 'UniverseFilter',
-    'WeightFilter', 'ZernikeFilter', 'ZernikeRadialFilter', 'filters'
+    'ParentNuclideFilter', 'ParticleFilter', 'ParticleoutFilter', 'PolarFilter',
+    'SphericalHarmonicsFilter', 'SpatialLegendreFilter', 'SurfaceFilter', 'TimeFilter',
+    'UniverseFilter', 'WeightFilter', 'ZernikeFilter', 'ZernikeRadialFilter', 'filters'
 ]
 
 # Tally functions
@@ -564,6 +564,10 @@ def bins(self):
         return [ParticleType(i) for i in particle_i]
 
 
+class ParticleoutFilter(ParticleFilter):
+    filter_type = 'particleout'
+
+
 class PolarFilter(Filter):
     filter_type = 'polar'
 

src/particle.cpp

@@ -38,6 +38,10 @@
 
 namespace openmc {
 
+namespace simulation {
+thread_local Particle tmp_particle;
+}
+
 //==============================================================================
 // Particle implementation
 //==============================================================================
@@ -88,6 +92,24 @@ bool Particle::create_secondary(
   bank.E = settings::run_CE ? E : g();
   bank.time = time();
   bank_second_E() += bank.E;
+
+  // Score tallies affected by secondary particles
+  if (!model::active_particleout_analog_tallies.empty()) {
+    // Create secondary particle for tallying purposes only
+    simulation::tmp_particle.from_source(&bank);
+    simulation::tmp_particle.u_last() = this->u();
+    simulation::tmp_particle.r_last() = this->r();
+    simulation::tmp_particle.E_last() = this->E();
+    simulation::tmp_particle.type_last() = this->type();
+
+    if (settings::run_CE) {
+      score_analog_tally_ce(
+        simulation::tmp_particle, model::active_particleout_analog_tallies);
+    } else {
+      score_analog_tally_mg(
+        simulation::tmp_particle, model::active_particleout_analog_tallies);
+    }
+  }
   return true;
 }
 
@@ -124,6 +146,7 @@ void Particle::from_source(const SourceSite* src)
 
   // Copy attributes from source bank site
   type() = src->particle;
+  type_last() = src->particle;
   wgt() = src->wgt;
   wgt_last() = src->wgt;
   r() = src->r;
@@ -160,6 +183,7 @@ void Particle::event_calculate_xs()
   // Store pre-collision particle properties
   wgt_last() = wgt();
   E_last() = E();
+  type_last() = type();
   u_last() = u();
   r_last() = r();
   time_last() = time();
@@ -353,9 +377,9 @@ void Particle::event_collide()
     score_collision_tally(*this);
   if (!model::active_analog_tallies.empty()) {
     if (settings::run_CE) {
-      score_analog_tally_ce(*this);
+      score_analog_tally_ce(*this, model::active_analog_tallies);
     } else {
-      score_analog_tally_mg(*this);
+      score_analog_tally_mg(*this, model::active_analog_tallies);
     }
   }
 

src/particle_restart.cpp

@@ -64,6 +64,7 @@ void read_particle_restart(Particle& p, RunMode& previous_run_mode)
   p.r_last() = p.r();
   p.u_last() = p.u();
   p.E_last() = p.E();
+  p.type_last() = p.type();
   p.g_last() = p.g();
   p.time_last() = p.time();
 

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 == "particleout") {
+    return Filter::create<ParticleOutFilter>(id);
   } else if (type == "polar") {
     return Filter::create<PolarFilter>(id);
   } else if (type == "surface") {

src/tallies/filter_particle.cpp

@@ -6,6 +6,10 @@
 
 namespace openmc {
 
+//==============================================================================
+// ParticleFilter implementation
+//==============================================================================
+
 void ParticleFilter::from_xml(pugi::xml_node node)
 {
   auto particles = get_node_array<std::string>(node, "bins");
@@ -34,8 +38,11 @@ void ParticleFilter::set_particles(span<ParticleType> particles)
 void ParticleFilter::get_all_bins(
   const Particle& p, TallyEstimator estimator, FilterMatch& match) const
 {
+  // Get the pre-collision type of the particle.
+  auto type = p.type_last();
+
   for (auto i = 0; i < particles_.size(); i++) {
-    if (particles_[i] == p.type()) {
+    if (particles_[i] == type) {
       match.bins_.push_back(i);
       match.weights_.push_back(1.0);
     }
@@ -58,6 +65,66 @@ std::string ParticleFilter::text_label(int bin) const
   return fmt::format("Particle: {}", particle_type_to_str(p));
 }
 
+//==============================================================================
+// ParticleOutFilter implementation
+//==============================================================================
+
+void ParticleOutFilter::from_xml(pugi::xml_node node)
+{
+  auto particles = get_node_array<std::string>(node, "bins");
+
+  // Convert to vector of ParticleType
+  vector<ParticleType> types;
+  for (auto& p : particles) {
+    types.push_back(str_to_particle_type(p));
+  }
+  this->set_particles(types);
+}
+
+void ParticleOutFilter::set_particles(span<ParticleType> particles)
+{
+  // Clear existing particles
+  particles_.clear();
+  particles_.reserve(particles.size());
+
+  // Set particles and number of bins
+  for (auto p : particles) {
+    particles_.push_back(p);
+  }
+  n_bins_ = particles_.size();
+}
+
+void ParticleOutFilter::get_all_bins(
+  const Particle& p, TallyEstimator estimator, FilterMatch& match) const
+{
+  for (auto i = 0; i < particles_.size(); i++) {
+    if (particles_[i] == p.type()) {
+      match.bins_.push_back(i);
+      match.weights_.push_back(1.0);
+    }
+  }
+}
+
+void ParticleOutFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  vector<std::string> particles;
+  for (auto p : particles_) {
+    particles.push_back(particle_type_to_str(p));
+  }
+  write_dataset(filter_group, "bins", particles);
+}
+
+std::string ParticleOutFilter::text_label(int bin) const
+{
+  const auto& p = particles_.at(bin);
+  return fmt::format("ParticleOut: {}", particle_type_to_str(p));
+}
+
+//==============================================================================
+// C-API functions
+//==============================================================================
+
 extern "C" int openmc_particle_filter_get_bins(int32_t idx, int bins[])
 {
   if (int err = verify_filter(idx))