Merge pull request #39 from zxkjack123/point-detector-contrib

Fix point detector scoring: activate tallies, set ParticleRay state, fix Python reader

Author: GuySten

include/openmc/tallies/tally_scoring.h

@@ -158,6 +158,15 @@ void score_point_tally_impl(
     p.wgt_last() = p.wgt();
     p.wgt() *= attenuation;
 
+    // Set position state so PointFilter::get_all_bins can match:
+    // r() = detector position (for bin matching),
+    // r_last() = source/scatter position (for 1/r^2 weight)
+    p.r_last() = r;
+    p.r() = det;
+
+    // Set E_last so EnergyFilter::get_all_bins bins correctly
+    p.E_last() = E;
+
     double flux = p.wgt_last() * pdf;
     score_tracklength_tally_general(p, flux, model::active_point_tallies);
   }

openmc/filter.py

@@ -834,6 +834,28 @@ def bins(self, bins):
             cv.check_length(f'bins[{i}][0]', item[0], 3, 3)
             cv.check_type(f'bins[{i}][1]', item[1], Real)
         self._bins = bins
+
+    @classmethod
+    def from_hdf5(cls, group, **kwargs):
+        filter_id = int(group.name.split('/')[-1].lstrip('filter '))
+        flat = group['bins'][()]
+        # Reconstruct tuple structure: every 4 values = (x, y, z, r0)
+        bins = []
+        for i in range(0, len(flat), 4):
+            pos = (float(flat[i]), float(flat[i+1]), float(flat[i+2]))
+            r0 = float(flat[i+3])
+            bins.append((pos, r0))
+        out = cls(bins, filter_id=filter_id)
+        out._num_bins = group['n_bins'][()]
+        return out
+
+    def get_pandas_dataframe(self, data_size, stride, **kwargs):
+        import pandas as pd
+        labels = [f"({p[0]}, {p[1]}, {p[2]}) R0={r}" for (p, r) in self.bins]
+        filter_bins = np.repeat(labels, stride)
+        tile_factor = data_size // len(filter_bins)
+        filter_bins = np.tile(filter_bins, tile_factor)
+        return pd.DataFrame({self.short_name.lower(): filter_bins})
 
     def to_xml_element(self):
         """Return XML Element representing the Filter.

openmc/tallies.py

@@ -53,7 +53,7 @@
 _FILTER_CLASSES = (Filter, CrossFilter, AggregateFilter)
 
 # Valid types of estimators
-ESTIMATOR_TYPES = {'tracklength', 'collision', 'analog'}
+ESTIMATOR_TYPES = {'tracklength', 'collision', 'analog', 'next-event'}
 
 
 class Tally(IDManagerMixin):

src/state_point.cpp

@@ -214,6 +214,8 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source)
           write_dataset(tally_group, "estimator", "tracklength");
         } else if (tally->estimator_ == TallyEstimator::COLLISION) {
           write_dataset(tally_group, "estimator", "collision");
+        } else if (tally->estimator_ == TallyEstimator::NEXT_EVENT) {
+          write_dataset(tally_group, "estimator", "next-event");
         }
 
         write_dataset(tally_group, "n_realizations", tally->n_realizations_);

src/tallies/tally.cpp

@@ -31,6 +31,7 @@
 #include "openmc/tallies/filter_meshmaterial.h"
 #include "openmc/tallies/filter_meshsurface.h"
 #include "openmc/tallies/filter_particle.h"
+#include "openmc/tallies/filter_point.h"
 #include "openmc/tallies/filter_sph_harm.h"
 #include "openmc/tallies/filter_surface.h"
 #include "openmc/tallies/filter_time.h"
@@ -1199,6 +1200,8 @@ void setup_active_tallies()
   model::active_meshsurf_tallies.clear();
   model::active_surface_tallies.clear();
   model::active_pulse_height_tallies.clear();
+  model::active_point_tallies.clear();
+  model::active_point_detectors.clear();
   model::time_grid.clear();
 
   for (auto i = 0; i < model::tallies.size(); ++i) {
@@ -1240,6 +1243,16 @@ void setup_active_tallies()
       case TallyType::PULSE_HEIGHT:
         model::active_pulse_height_tallies.push_back(i);
         break;
+
+      case TallyType::POINT:
+        model::active_point_tallies.push_back(i);
+        // Populate the set of unique detector positions from PointFilter
+        if (auto pf = tally.get_filter<PointFilter>()) {
+          for (const auto& [pos, r0] : pf->detectors()) {
+            model::active_point_detectors.insert(pos);
+          }
+        }
+        break;
       }
     }
   }
@@ -1263,6 +1276,8 @@ void free_memory_tally()
   model::active_meshsurf_tallies.clear();
   model::active_surface_tallies.clear();
   model::active_pulse_height_tallies.clear();
+  model::active_point_tallies.clear();
+  model::active_point_detectors.clear();
   model::time_grid.clear();
 
   model::tally_map.clear();