Pulse height tally for photons

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

@@ -278,7 +278,7 @@ enum class MgxsType {
 
 enum class TallyResult { VALUE, SUM, SUM_SQ };
 
-enum class TallyType { VOLUME, MESH_SURFACE, SURFACE };
+enum class TallyType { VOLUME, MESH_SURFACE, SURFACE, PULSE_HEIGHT };
 
 enum class TallyEstimator { ANALOG, TRACKLENGTH, COLLISION };
 
@@ -306,7 +306,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

@@ -63,6 +63,10 @@ class Particle : public ParticleData {
   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

@@ -312,6 +312,8 @@ class ParticleData {
 
   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};
@@ -443,6 +445,7 @@ class ParticleData {
   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

@@ -121,8 +121,10 @@ class Tally {
   // 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_;
 
@@ -155,6 +157,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

@@ -98,13 +98,13 @@
     -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"
@@ -286,6 +287,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;
@@ -350,6 +356,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);
@@ -383,6 +408,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) {
@@ -391,6 +420,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());
@@ -433,7 +497,9 @@ void Particle::cross_surface()
 #ifdef DAGMC
   // in DAGMC, we know what the next cell should be
   if (surf->geom_type_ == GeometryType::DAG) {
-    int32_t i_cell = next_cell(i_surface, cell_last(n_coord() - 1), lowest_coord().universe) - 1;
+    int32_t i_cell =
+      next_cell(i_surface, cell_last(n_coord() - 1), lowest_coord().universe) -
+      1;
     // save material and temp
     material_last() = material();
     sqrtkT_last() = sqrtkT();

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

@@ -533,6 +533,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 +