Fix not reseting filter_matches corrupt pulse height results (#4019)

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

src/tallies/tally_scoring.cpp

@@ -2730,6 +2730,9 @@ void score_pulse_height_tally(Particle& p, const vector<int>& tallies)
   int orig_cell = p.coord(0).cell();
   double orig_E_last = p.E_last();
 
+  // Set particle in top level
+  p.n_coord() = 1;
+
   for (auto i_tally : tallies) {
     auto& tally {*model::tallies[i_tally]};
 
@@ -2740,7 +2743,6 @@ void score_pulse_height_tally(Particle& p, const vector<int>& tallies)
 
     for (auto cell_id : cells) {
       // Temporarily change cell of particle
-      p.n_coord() = 1;
       p.coord(0).cell() = cell_id;
 
       // Determine index of cell in model::pulse_height_cells
@@ -2756,20 +2758,20 @@ void score_pulse_height_tally(Particle& p, const vector<int>& tallies)
       // we skip the assume_separate break below.
       auto filter_iter = FilterBinIter(tally, p);
       auto end = FilterBinIter(tally, true, &p.filter_matches());
-      if (filter_iter == end)
+      if (filter_iter != end) {
-        continue;
 
-      // Loop over filter bins.
+        // Loop over filter bins.
-      for (; filter_iter != end; ++filter_iter) {
+        for (; filter_iter != end; ++filter_iter) {
-        auto filter_index = filter_iter.index_;
+          auto filter_index = filter_iter.index_;
-        auto filter_weight = filter_iter.weight_;
+          auto filter_weight = filter_iter.weight_;
 
-        // Loop over scores.
+          // Loop over scores.
-        for (auto score_index = 0; score_index < tally.scores_.size();
+          for (auto score_index = 0; score_index < tally.scores_.size();
-             ++score_index) {
+               ++score_index) {
 #pragma omp atomic
-          tally.results_(filter_index, score_index, TallyResult::VALUE) +=
+            tally.results_(filter_index, score_index, TallyResult::VALUE) +=
-            filter_weight;
+              filter_weight;
+          }
         }
       }
 
@@ -2777,10 +2779,10 @@ void score_pulse_height_tally(Particle& p, const vector<int>& tallies)
       for (auto& match : p.filter_matches())
         match.bins_present_ = false;
     }
-    // Restore cell/energy
-    p.n_coord() = orig_n_coord;
-    p.coord(0).cell() = orig_cell;
-    p.E_last() = orig_E_last;
   }
+  // Restore cell/energy
+  p.n_coord() = orig_n_coord;
+  p.coord(0).cell() = orig_cell;
+  p.E_last() = orig_E_last;
 }
 } // namespace openmc

tests/unit_tests/test_pulse_height.py

@@ -0,0 +1,59 @@
+import numpy as np
+import pytest
+import openmc
+
+
+@pytest.fixture
+def model():
+    openmc.reset_auto_ids()
+    model = openmc.Model()
+
+    # Define materials
+    NaI = openmc.Material()
+    NaI.set_density('g/cm3', 3.7)
+    NaI.add_element('Na', 1.0)
+    NaI.add_element('I', 1.0)
+
+    # Define geometry: two spheres in each other
+    s1 = openmc.Sphere(r=1)
+    s2 = openmc.Sphere(r=2, boundary_type='vacuum')
+    inner_sphere = openmc.Cell(name='inner sphere', fill=NaI, region=-s1)
+    outer_sphere = openmc.Cell(name='outer sphere', fill=NaI, region=+s1 & -s2)
+    model.geometry = openmc.Geometry([inner_sphere, outer_sphere])
+
+    # Define settings
+    model.settings.run_mode = 'fixed source'
+    model.settings.batches = 1
+    model.settings.particles = 10000
+    model.settings.photon_transport = True
+    model.settings.source = openmc.IndependentSource(
+        energy=openmc.stats.delta_function(1e6),
+        particle='photon'
+    )
+
+    # Define tallies
+    energy_filter = openmc.EnergyFilter([1e3, 1e7])    
+
+    tally1 = openmc.Tally()
+    tally1.scores = ['pulse-height']
+    cell_filter1 = openmc.CellFilter([inner_sphere, outer_sphere])
+    tally1.filters = [cell_filter1, energy_filter]
+
+    tally2 = openmc.Tally()
+    tally2.scores = ['pulse-height']
+    cell_filter2 = openmc.CellFilter([outer_sphere, inner_sphere])
+    tally2.filters = [cell_filter2, energy_filter]
+
+    model.tallies = [tally1, tally2]
+    return model
+
+
+def test_pulse_height(model, run_in_tmpdir):
+    sp_path = model.run()
+    sp = openmc.StatePoint(sp_path)
+    t1 = sp.tallies[1].mean.squeeze()
+    t2 = sp.tallies[2].mean.squeeze()
+    
+    np.testing.assert_array_equal(t1, t2[::-1])
+
+