Pulse height tally for photons (#2452)

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

docs/source/usersguide/tallies.rst

@@ -317,6 +317,11 @@ The following tables show all valid scores:
     |                      |particle. This corresponds to MT=444 produced by   |
     |                      |NJOY's HEATR module.                               |
     +----------------------+---------------------------------------------------+
+    |pulse-height          |The energy deposited by an entire photon's history |
+    |                      |(including its progeny). Units are eV per source   | 
+    |                      |particle. Note that this score can only be combined|
+    |                      |with a cell filter and an energy filter.           |   
+    +----------------------+---------------------------------------------------+
 
 .. _usersguide_tally_normalization:
 

include/openmc/constants.h

@@ -279,7 +279,7 @@ enum class MgxsType {
 
 enum class TallyResult { VALUE, SUM, SUM_SQ, SIZE };
 
-enum class TallyType { VOLUME, MESH_SURFACE, SURFACE };
+enum class TallyType { VOLUME, MESH_SURFACE, SURFACE, PULSE_HEIGHT };
 
 enum class TallyEstimator { ANALOG, TRACKLENGTH, COLLISION };
 
@@ -307,7 +307,8 @@ enum TallyScore {
   SCORE_INVERSE_VELOCITY = -13,   // flux-weighted inverse velocity
   SCORE_FISS_Q_PROMPT = -14,      // prompt fission Q-value
   SCORE_FISS_Q_RECOV = -15,       // recoverable fission Q-value
-  SCORE_DECAY_RATE = -16          // delayed neutron precursor decay rate
+  SCORE_DECAY_RATE = -16,         // delayed neutron precursor decay rate
+  SCORE_PULSE_HEIGHT = -17        // pulse-height
 };
 
 // Global tally parameters

include/openmc/particle.h

@@ -66,6 +66,10 @@ public:
   void event_revive_from_secondary();
   void event_death();
 
+  //! pulse-height recording
+  void pht_collision_energy();
+  void pht_secondary_particles();
+
   //! Cross a surface and handle boundary conditions
   void cross_surface();
 

include/openmc/particle_data.h

@@ -313,6 +313,8 @@ private:
 
   vector<NuBank> nu_bank_; // bank of most recently fissioned particles
 
+  vector<double> pht_storage_; // interim pulse-height results
+
   // Global tally accumulators
   double keff_tally_absorption_ {0.0};
   double keff_tally_collision_ {0.0};
@@ -444,6 +446,7 @@ public:
   decltype(tracks_)& tracks() { return tracks_; }
   decltype(nu_bank_)& nu_bank() { return nu_bank_; }
   NuBank& nu_bank(int i) { return nu_bank_[i]; }
+  vector<double>& pht_storage() { return pht_storage_; }
 
   double& keff_tally_absorption() { return keff_tally_absorption_; }
   double& keff_tally_collision() { return keff_tally_collision_; }

include/openmc/tallies/tally.h

@@ -169,8 +169,10 @@ public:
   // We need to have quick access to some filters.  The following gives indices
   // for various filters that could be in the tally or C_NONE if they are not
   // present.
+  int energy_filter_ {C_NONE};
   int energyout_filter_ {C_NONE};
   int delayedgroup_filter_ {C_NONE};
+  int cell_filter_ {C_NONE};
 
   vector<Trigger> triggers_;
 
@@ -206,6 +208,8 @@ extern vector<int> active_tracklength_tallies;
 extern vector<int> active_collision_tallies;
 extern vector<int> active_meshsurf_tallies;
 extern vector<int> active_surface_tallies;
+extern vector<int> active_pulse_height_tallies;
+extern vector<int> pulse_height_cells;
 } // namespace model
 
 namespace simulation {

include/openmc/tallies/tally_scoring.h

@@ -94,9 +94,16 @@ void score_tracklength_tally(Particle& p, double distance);
 //! Score surface or mesh-surface tallies for particle currents.
 //
 //! \param p The particle being tracked
-//! \param tallies A vector of tallies to score to
+//! \param tallies A vector of the indices of the tallies to score to
 void score_surface_tally(Particle& p, const vector<int>& tallies);
 
+//! Score the pulse-height tally
+//! This is triggered at the end of every particle history
+//
+//! \param p The particle being tracked
+//! \param tallies A vector of the indices of the tallies to score to
+void score_pulse_height_tally(Particle& p, const vector<int>& tallies);
+
 } // namespace openmc
 
 #endif // OPENMC_TALLIES_TALLY_SCORING_H

openmc/lib/tally.py

@@ -104,13 +104,13 @@ _SCORES = {
     -5: 'absorption', -6: 'fission', -7: 'nu-fission', -8: 'kappa-fission',
     -9: 'current', -10: 'events', -11: 'delayed-nu-fission',
     -12: 'prompt-nu-fission', -13: 'inverse-velocity', -14: 'fission-q-prompt',
-    -15: 'fission-q-recoverable', -16: 'decay-rate'
+    -15: 'fission-q-recoverable', -16: 'decay-rate', -17: 'pulse-height'
 }
 _ESTIMATORS = {
     0: 'analog', 1: 'tracklength', 2: 'collision'
 }
 _TALLY_TYPES = {
-    0: 'volume', 1: 'mesh-surface', 2: 'surface'
+    0: 'volume', 1: 'mesh-surface', 2: 'surface', 3: 'pulse-height'
 }
 
 

src/output.cpp

@@ -594,6 +594,7 @@ const std::unordered_map<int, const char*> score_names = {
   {SCORE_FISS_Q_PROMPT, "Prompt fission power"},
   {SCORE_FISS_Q_RECOV, "Recoverable fission power"},
   {SCORE_CURRENT, "Current"},
+  {SCORE_PULSE_HEIGHT, "pulse-height"},
 };
 
 //! Create an ASCII output file showing all tally results.

src/particle.cpp

@@ -17,6 +17,7 @@
 #include "openmc/message_passing.h"
 #include "openmc/mgxs_interface.h"
 #include "openmc/nuclide.h"
+#include "openmc/particle_data.h"
 #include "openmc/photon.h"
 #include "openmc/physics.h"
 #include "openmc/physics_mg.h"
@@ -290,6 +291,11 @@ void Particle::event_collide()
     }
   }
 
+  if (!model::active_pulse_height_tallies.empty() &&
+      type() == ParticleType::photon) {
+    pht_collision_energy();
+  }
+
   // Reset banked weight during collision
   n_bank() = 0;
   n_bank_second() = 0;
@@ -354,6 +360,25 @@ void Particle::event_revive_from_secondary()
     secondary_bank().pop_back();
     n_event() = 0;
 
+    // Subtract secondary particle energy from interim pulse-height results
+    if (!model::active_pulse_height_tallies.empty() &&
+        this->type() == ParticleType::photon) {
+      // Since the birth cell of the particle has not been set we
+      // have to determine it before the energy of the secondary particle can be
+      // removed from the pulse-height of this cell.
+      if (coord(n_coord() - 1).cell == C_NONE) {
+        if (!exhaustive_find_cell(*this)) {
+          mark_as_lost("Could not find the cell containing particle " +
+                       std::to_string(id()));
+          return;
+        }
+        // Set birth cell attribute
+        if (cell_born() == C_NONE)
+          cell_born() = coord(n_coord() - 1).cell;
+      }
+      pht_secondary_particles();
+    }
+
     // Enter new particle in particle track file
     if (write_track())
       add_particle_track(*this);
@@ -387,6 +412,10 @@ void Particle::event_death()
   keff_tally_tracklength() = 0.0;
   keff_tally_leakage() = 0.0;
 
+  if (!model::active_pulse_height_tallies.empty()) {
+    score_pulse_height_tally(*this, model::active_pulse_height_tallies);
+  }
+
   // Record the number of progeny created by this particle.
   // This data will be used to efficiently sort the fission bank.
   if (settings::run_mode == RunMode::EIGENVALUE) {
@@ -395,6 +424,41 @@ void Particle::event_death()
   }
 }
 
+void Particle::pht_collision_energy()
+{
+  // Adds the energy particles lose in a collision to the pulse-height
+
+  // determine index of cell in pulse_height_cells
+  auto it = std::find(model::pulse_height_cells.begin(),
+    model::pulse_height_cells.end(), coord(n_coord() - 1).cell);
+
+  if (it != model::pulse_height_cells.end()) {
+    int index = std::distance(model::pulse_height_cells.begin(), it);
+    pht_storage()[index] += E_last() - E();
+
+    // If the energy of the particle is below the cutoff, it will not be sampled
+    // so its energy is added to the pulse-height in the cell
+    int photon = static_cast<int>(ParticleType::photon);
+    if (E() < settings::energy_cutoff[photon]) {
+      pht_storage()[index] += E();
+    }
+  }
+}
+
+void Particle::pht_secondary_particles()
+{
+  // Removes the energy of secondary produced particles from the pulse-height
+
+  // determine index of cell in pulse_height_cells
+  auto it = std::find(model::pulse_height_cells.begin(),
+    model::pulse_height_cells.end(), cell_born());
+
+  if (it != model::pulse_height_cells.end()) {
+    int index = std::distance(model::pulse_height_cells.begin(), it);
+    pht_storage()[index] -= E();
+  }
+}
+
 void Particle::cross_surface()
 {
   int i_surface = std::abs(surface());

src/particle_data.cpp

@@ -53,6 +53,11 @@ ParticleData::ParticleData()
   // Create microscopic cross section caches
   neutron_xs_.resize(data::nuclides.size());
   photon_xs_.resize(data::elements.size());
+
+  // Creates the pulse-height storage for the particle
+  if (!model::pulse_height_cells.empty()) {
+    pht_storage_.resize(model::pulse_height_cells.size(), 0.0);
+  }
 }
 
 TrackState ParticleData::get_track_state() const

src/reaction.cpp

@@ -171,6 +171,7 @@ std::unordered_map<int, std::string> REACTION_NAME_MAP {
   {SCORE_INVERSE_VELOCITY, "inverse-velocity"},
   {SCORE_FISS_Q_PROMPT, "fission-q-prompt"},
   {SCORE_FISS_Q_RECOV, "fission-q-recoverable"},
+  {SCORE_PULSE_HEIGHT, "pulse-height"},
   // Normal ENDF-based reactions
   {TOTAL_XS, "(n,total)"},
   {ELASTIC, "(n,elastic)"},

src/simulation.cpp

@@ -550,6 +550,9 @@ void initialize_history(Particle& p, int64_t index_source)
   // Reset weight window ratio
   p.ww_factor() = 0.0;
 
+  // Reset pulse_height_storage
+  std::fill(p.pht_storage().begin(), p.pht_storage().end(), 0);
+
   // set random number seed
   int64_t particle_seed =
     (simulation::total_gen + overall_generation() - 1) * settings::n_particles +

src/tallies/tally.cpp

@@ -3,6 +3,7 @@
 #include "openmc/array.h"
 #include "openmc/capi.h"
 #include "openmc/constants.h"
+#include "openmc/container_util.h"
 #include "openmc/error.h"
 #include "openmc/file_utils.h"
 #include "openmc/mesh.h"
@@ -18,6 +19,7 @@
 #include "openmc/tallies/derivative.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/tallies/filter_cell.h"
+#include "openmc/tallies/filter_cellborn.h"
 #include "openmc/tallies/filter_cellfrom.h"
 #include "openmc/tallies/filter_collision.h"
 #include "openmc/tallies/filter_delayedgroup.h"
@@ -55,6 +57,8 @@ vector<int> active_tracklength_tallies;
 vector<int> active_collision_tallies;
 vector<int> active_meshsurf_tallies;
 vector<int> active_surface_tallies;
+vector<int> active_pulse_height_tallies;
+vector<int> pulse_height_cells;
 } // namespace model
 
 namespace simulation {
@@ -175,6 +179,16 @@ Tally::Tally(pugi::xml_node node)
   }
 
   // Check if tally is compatible with particle type
+  if (!settings::photon_transport) {
+    for (int score : scores_) {
+      switch (score) {
+      case SCORE_PULSE_HEIGHT:
+        fatal_error(
+          "For pulse-height tallies, photon transport needs to be activated.");
+        break;
+      }
+    }
+  }
   if (settings::photon_transport) {
     if (particle_filter_index == C_NONE) {
       for (int score : scores_) {
@@ -409,10 +423,14 @@ void Tally::add_filter(Filter* filter)
     return;
 
   // Keep track of indices for special filters
-  if (dynamic_cast<const EnergyoutFilter*>(filter)) {
+  if (filter->type() == FilterType::ENERGY_OUT) {
     energyout_filter_ = filters_.size();
-  } else if (dynamic_cast<const DelayedGroupFilter*>(filter)) {
+  } else if (filter->type() == FilterType::DELAYED_GROUP) {
     delayedgroup_filter_ = filters_.size();
+  } else if (filter->type() == FilterType::CELL) {
+    cell_filter_ = filters_.size();
+  } else if (filter->type() == FilterType::ENERGY) {
+    energy_filter_ = filters_.size();
   }
   filters_.push_back(filter_idx);
 }
@@ -454,17 +472,23 @@ void Tally::set_scores(const vector<std::string>& scores)
   bool cellfrom_present = false;
   bool surface_present = false;
   bool meshsurface_present = false;
+  bool non_cell_energy_present = false;
   for (auto i_filt : filters_) {
     const auto* filt {model::tally_filters[i_filt].get()};
-    if (dynamic_cast<const LegendreFilter*>(filt)) {
+    // Checking for only cell and energy filters for pulse-height tally
+    if (!(filt->type() == FilterType::CELL ||
+          filt->type() == FilterType::ENERGY)) {
+      non_cell_energy_present = true;
+    }
+    if (filt->type() == FilterType::LEGENDRE) {
       legendre_present = true;
-    } else if (dynamic_cast<const CellFromFilter*>(filt)) {
+    } else if (filt->type() == FilterType::CELLFROM) {
       cellfrom_present = true;
-    } else if (dynamic_cast<const CellFilter*>(filt)) {
+    } else if (filt->type() == FilterType::CELL) {
       cell_present = true;
-    } else if (dynamic_cast<const SurfaceFilter*>(filt)) {
+    } else if (filt->type() == FilterType::SURFACE) {
       surface_present = true;
-    } else if (dynamic_cast<const MeshSurfaceFilter*>(filt)) {
+    } else if (filt->type() == FilterType::MESH_SURFACE) {
       meshsurface_present = true;
     }
   }
@@ -536,6 +560,31 @@ void Tally::set_scores(const vector<std::string>& scores)
       if (settings::photon_transport)
         estimator_ = TallyEstimator::COLLISION;
       break;
+
+    case SCORE_PULSE_HEIGHT:
+      if (non_cell_energy_present) {
+        fatal_error("Pulse-height tallies are not compatible with filters "
+                    "other than CellFilter and EnergyFilter");
+      }
+      type_ = TallyType::PULSE_HEIGHT;
+
+      // Collecting indices of all cells covered by the filters in the pulse
+      // height tally in global variable pulse_height_cells
+      for (const auto& i_filt : filters_) {
+        auto cell_filter =
+          dynamic_cast<CellFilter*>(model::tally_filters[i_filt].get());
+        if (cell_filter) {
+          const auto& cells = cell_filter->cells();
+          for (int i = 0; i < cell_filter->n_bins(); i++) {
+            int cell_index = cells[i];
+            if (!contains(model::pulse_height_cells, cell_index)) {
+              model::pulse_height_cells.push_back(cell_index);
+            }
+          }
+        }
+      }
+
+      break;
     }
 
     scores_.push_back(score);
@@ -953,6 +1002,7 @@ void setup_active_tallies()
   model::active_collision_tallies.clear();
   model::active_meshsurf_tallies.clear();
   model::active_surface_tallies.clear();
+  model::active_pulse_height_tallies.clear();
 
   for (auto i = 0; i < model::tallies.size(); ++i) {
     const auto& tally {*model::tallies[i]};
@@ -980,6 +1030,11 @@ void setup_active_tallies()
 
       case TallyType::SURFACE:
         model::active_surface_tallies.push_back(i);
+        break;
+
+      case TallyType::PULSE_HEIGHT:
+        model::active_pulse_height_tallies.push_back(i);
+        break;
       }
     }
   }
@@ -1001,6 +1056,7 @@ void free_memory_tally()
   model::active_collision_tallies.clear();
   model::active_meshsurf_tallies.clear();
   model::active_surface_tallies.clear();
+  model::active_pulse_height_tallies.clear();
 
   model::tally_map.clear();
 }
@@ -1122,6 +1178,8 @@ extern "C" int openmc_tally_set_type(int32_t index, const char* type)
     model::tallies[index]->type_ = TallyType::MESH_SURFACE;
   } else if (strcmp(type, "surface") == 0) {
     model::tallies[index]->type_ = TallyType::SURFACE;
+  } else if (strcmp(type, "pulse-height") == 0) {
+    model::tallies[index]->type_ = TallyType::PULSE_HEIGHT;
   } else {
     set_errmsg(fmt::format("Unknown tally type: {}", type));
     return OPENMC_E_INVALID_ARGUMENT;

src/tallies/tally_scoring.cpp

@@ -15,6 +15,7 @@
 #include "openmc/string_utils.h"
 #include "openmc/tallies/derivative.h"
 #include "openmc/tallies/filter.h"
+#include "openmc/tallies/filter_cell.h"
 #include "openmc/tallies/filter_delayedgroup.h"
 #include "openmc/tallies/filter_energy.h"
 
@@ -2498,4 +2499,79 @@ void score_surface_tally(Particle& p, const vector<int>& tallies)
     match.bins_present_ = false;
 }
 
+void score_pulse_height_tally(Particle& p, const vector<int>& tallies)
+{
+  // The pulse height tally in OpenMC hijacks the logic of CellFilter and
+  // EnergyFilter to score specific quantities related to particle pulse height.
+  // This is achieved by setting the pulse-height cell of the tally to the cell
+  // of the particle being scored, and the energy to the particle's last
+  // recorded energy (E_last()). After the tally is scored, the values are reset
+  // to ensure proper accounting and avoid interference with subsequent
+  // calculations or tallies.
+
+  // Save original cell/energy information
+  int orig_n_coord = p.n_coord();
+  int orig_cell = p.coord(0).cell;
+  double orig_E_last = p.E_last();
+
+  for (auto i_tally : tallies) {
+    auto& tally {*model::tallies[i_tally]};
+
+    // Determine all CellFilter in the tally
+    for (const auto& filter : tally.filters()) {
+      auto cell_filter =
+        dynamic_cast<CellFilter*>(model::tally_filters[filter].get());
+      if (cell_filter != nullptr) {
+
+        const auto& cells = cell_filter->cells();
+        // Loop over all cells in the CellFilter
+        for (auto cell_index = 0; cell_index < cells.size(); ++cell_index) {
+          int cell_id = cells[cell_index];
+
+          // Temporarily change cell of particle
+          p.n_coord() = 1;
+          p.coord(0).cell = cell_id;
+
+          // Determine index of cell in model::pulse_height_cells
+          auto it = std::find(model::pulse_height_cells.begin(),
+            model::pulse_height_cells.end(), cell_id);
+          int index = std::distance(model::pulse_height_cells.begin(), it);
+
+          // Temporarily change energy of particle to pulse-height value
+          p.E_last() = p.pht_storage()[index];
+
+          // Initialize an iterator over valid filter bin combinations. If
+          // there are no valid combinations, use a continue statement to ensure
+          // 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)
+            continue;
+
+          // Loop over filter bins.
+          for (; filter_iter != end; ++filter_iter) {
+            auto filter_index = filter_iter.index_;
+            auto filter_weight = filter_iter.weight_;
+
+            // Loop over scores.
+            for (auto score_index = 0; score_index < tally.scores_.size();
+                 ++score_index) {
+#pragma omp atomic
+              tally.results_(filter_index, score_index, TallyResult::VALUE) +=
+                filter_weight;
+            }
+          }
+
+          // Reset all the filter matches for the next tally event.
+          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;
+  }
+}
 } // namespace openmc

tests/regression_tests/pulse_height/inputs_true.dat

@@ -0,0 +1,42 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1">
+      <density units="g/cc" value="3.7"/>
+      <nuclide ao="1.0" name="Na23"/>
+      <nuclide ao="1.0" name="I127"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" name="inner sphere" region="-1" universe="1"/>
+    <cell id="2" material="void" name="outer sphere" region="1 -2" universe="1"/>
+    <surface coeffs="0.0 0.0 0.0 1" id="1" type="sphere"/>
+    <surface boundary="vacuum" coeffs="0.0 0.0 0.0 2" id="2" type="sphere"/>
+  </geometry>
+  <settings>
+    <run_mode>fixed source</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <source particle="photon" strength="1.0" type="independent">
+      <space type="point">
+        <parameters>0.0 0.0 0.0</parameters>
+      </space>
+      <energy type="discrete">
+        <parameters>1000000.0 1.0</parameters>
+      </energy>
+    </source>
+    <photon_transport>true</photon_transport>
+  </settings>
+  <tallies>
+    <filter id="1" type="cell">
+      <bins>1</bins>
+    </filter>
+    <filter id="2" type="energy">
+      <bins>0.0 10000.0 20000.0 30000.0 40000.0 50000.0 60000.0 70000.0 80000.0 90000.0 100000.0 110000.0 120000.0 130000.0 140000.0 150000.0 160000.0 170000.0 180000.0 190000.0 200000.0 210000.0 220000.0 230000.0 240000.0 250000.0 260000.0 270000.0 280000.0 290000.0 300000.0 310000.0 320000.0 330000.0 340000.0 350000.0 360000.0 370000.0 380000.0 390000.0 400000.0 410000.0 420000.0 430000.0 440000.0 450000.0 460000.0 470000.0 480000.0 490000.0 500000.0 510000.0 520000.0 530000.0 540000.0 550000.0 560000.0 570000.0 580000.0 590000.0 600000.0 610000.0 620000.0 630000.0 640000.0 650000.0 660000.0 670000.0 680000.0 690000.0 700000.0 710000.0 720000.0 730000.0 740000.0 750000.0 760000.0 770000.0 780000.0 790000.0 800000.0 810000.0 820000.0 830000.0 840000.0 850000.0 860000.0 870000.0 880000.0 890000.0 900000.0 910000.0 920000.0 930000.0 940000.0 950000.0 960000.0 970000.0 980000.0 990000.0 1000000.0</bins>
+    </filter>
+    <tally id="1" name="pht tally">
+      <filters>1 2</filters>
+      <scores>pulse-height</scores>
+    </tally>
+  </tallies>
+</model>

tests/regression_tests/pulse_height/results_true.dat

@@ -0,0 +1,201 @@
+tally 1:
+4.140000E+00
+3.443000E+00
+1.000000E-02
+1.000000E-04
+1.000000E-02
+1.000000E-04
+1.000000E-02
+1.000000E-04
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+1.000000E-02
+1.000000E-04
+1.000000E-02
+1.000000E-04
+1.000000E-02
+1.000000E-04
+0.000000E+00
+0.000000E+00
+1.000000E-02
+1.000000E-04
+3.000000E-02
+5.000000E-04
+2.000000E-02
+4.000000E-04
+2.000000E-02
+4.000000E-04
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+3.000000E-02
+3.000000E-04
+1.000000E-02
+1.000000E-04
+1.000000E-02
+1.000000E-04
+1.000000E-02
+1.000000E-04
+1.000000E-02
+1.000000E-04
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+1.000000E-02
+1.000000E-04
+1.000000E-02
+1.000000E-04
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+1.000000E-02
+1.000000E-04
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+4.000000E-02
+6.000000E-04
+1.000000E-02
+1.000000E-04
+0.000000E+00
+0.000000E+00
+2.000000E-02
+2.000000E-04
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+1.000000E-02
+1.000000E-04
+0.000000E+00
+0.000000E+00
+1.000000E-02
+1.000000E-04
+0.000000E+00
+0.000000E+00
+1.000000E-02
+1.000000E-04
+2.000000E-02
+4.000000E-04
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+2.000000E-02
+2.000000E-04
+1.000000E-02
+1.000000E-04
+1.000000E-02
+1.000000E-04
+1.000000E-02
+1.000000E-04
+1.000000E-02
+1.000000E-04
+3.000000E-02
+5.000000E-04
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+1.000000E-02
+1.000000E-04
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+2.000000E-02
+2.000000E-04
+0.000000E+00
+0.000000E+00
+1.000000E-02
+1.000000E-04
+2.000000E-02
+2.000000E-04
+1.000000E-02
+1.000000E-04
+1.000000E-02
+1.000000E-04
+2.000000E-02
+2.000000E-04
+1.000000E-02
+1.000000E-04
+1.000000E-02
+1.000000E-04
+3.000000E-02
+5.000000E-04
+0.000000E+00
+0.000000E+00
+1.000000E-02
+1.000000E-04
+0.000000E+00
+0.000000E+00
+1.000000E-02
+1.000000E-04
+1.000000E-02
+1.000000E-04
+1.000000E-02
+1.000000E-04
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+2.000000E-02
+2.000000E-04
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+2.000000E-01
+8.600000E-03

tests/regression_tests/pulse_height/test.py

@@ -0,0 +1,51 @@
+import numpy as np
+import openmc
+import pytest
+
+from tests.testing_harness import PyAPITestHarness
+
+
+@pytest.fixture
+def sphere_model():
+    
+    model = openmc.model.Model()
+
+    # Define materials
+    NaI = openmc.Material()
+    NaI.set_density('g/cc', 3.7)
+    NaI.add_element('Na', 1.0)
+    NaI.add_element('I', 1.0)
+
+    model.materials = openmc.Materials([NaI])
+
+    # 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', region=+s1 & -s2)
+    model.geometry = openmc.Geometry([inner_sphere, outer_sphere])
+
+    # Define settings
+    model.settings.run_mode = 'fixed source'
+    model.settings.batches = 5
+    model.settings.particles = 100
+    model.settings.photon_transport = True
+    model.settings.source = openmc.Source(space=openmc.stats.Point(),
+                                          energy=openmc.stats.Discrete([1e6],[1]),
+                                          particle='photon')
+
+    # Define tallies
+    tally = openmc.Tally(name="pht tally")
+    tally.scores = ['pulse-height']
+    cell_filter = openmc.CellFilter(inner_sphere)
+    energy_filter = openmc.EnergyFilter(np.linspace(0, 1_000_000, 101))
+    tally.filters = [cell_filter, energy_filter]
+    model.tallies = [tally]
+
+    return model
+
+
+
+def test_pulse_height(sphere_model):
+    harness = PyAPITestHarness('statepoint.5.h5', sphere_model)
+    harness.main()