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

@@ -75,9 +75,9 @@ constexpr double MASS_NEUTRON_EV {
   939.56542052e6};                             // mass of a neutron in eV/c^2
 constexpr double MASS_PROTON {1.007276466621}; // mass of a proton in amu
 constexpr double MASS_ELECTRON_EV {
-  0.51099895000e6}; // electron mass energy equivalent in eV/c^2
+  0.51099895000e6};      // electron mass energy equivalent in eV/c^2
 constexpr double FINE_STRUCTURE {
-  137.035999084}; // inverse fine structure constant
+  137.035999084};        // inverse fine structure constant
 constexpr double PLANCK_C {
   1.2398419839593942e4}; // Planck's constant times c in eV-Angstroms
 constexpr double AMU {1.66053906660e-27};      // 1 amu in kg
@@ -278,9 +278,9 @@ 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 };
+enum class TallyEstimator { ANALOG, TRACKLENGTH, COLLISION, PULSE_HEIGHT };
 
 enum class TallyEvent { SURFACE, LATTICE, KILL, SCATTER, ABSORB };
 
@@ -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/mcpl_interface.h

@@ -36,7 +36,7 @@ vector<SourceSite> mcpl_source_sites(std::string path);
 //!                         calculate_parallel_index_vector on
 //!                         source_bank.size().
 void write_mcpl_source_point(const char* filename,
-  gsl::span<SourceSite> source_bank, vector<int64_t> const& bank_index);
+  gsl::span<SourceSite> source_bank, const vector<int64_t>& bank_index);
 } // namespace openmc
 
 #endif // OPENMC_MCPL_INTERFACE_H

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

@@ -104,10 +104,10 @@ class LocalCoord {
 
 struct NuclideMicroXS {
   // Microscopic cross sections in barns
-  double total;      //!< total cross section
-  double absorption; //!< absorption (disappearance)
-  double fission;    //!< fission
-  double nu_fission; //!< neutron production from fission
+  double total;           //!< total cross section
+  double absorption;      //!< absorption (disappearance)
+  double fission;         //!< fission
+  double nu_fission;      //!< neutron production from fission
 
   double elastic;         //!< If sab_frac is not 1 or 0, then this value is
                           //!<   averaged over bound and non-bound nuclei
@@ -231,7 +231,7 @@ class ParticleData {
   vector<ElementMicroXS> photon_xs_;  //!< Microscopic photon cross sections
   MacroXS macro_xs_;                  //!< Macroscopic cross sections
 
-  int64_t id_;                                //!< Unique ID
+  int64_t id_;                        //!< Unique ID
   ParticleType type_ {ParticleType::neutron}; //!< Particle type (n, p, e, etc.)
 
   int n_coord_ {1};          //!< number of current coordinate levels
@@ -302,23 +302,25 @@ class ParticleData {
   // Secondary particle bank
   vector<SourceSite> secondary_bank_;
 
-  int64_t current_work_; // current work index
+  int64_t current_work_;               // current work index
 
-  vector<double> flux_derivs_; // for derivatives for this particle
+  vector<double> flux_derivs_;         // for derivatives for this particle
 
   vector<FilterMatch> filter_matches_; // tally filter matches
 
-  vector<TrackStateHistory> tracks_; // tracks for outputting to file
+  vector<TrackStateHistory> tracks_;   // tracks for outputting to file
 
-  vector<NuBank> nu_bank_; // bank of most recently fissioned particles
+  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};
   double keff_tally_tracklength_ {0.0};
   double keff_tally_leakage_ {0.0};
 
-  bool trace_ {false}; //!< flag to show debug information
+  bool trace_ {false};        //!< flag to show debug information
 
   double collision_distance_; // distance to particle's next closest collision
 
@@ -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/plot.h

@@ -29,7 +29,7 @@ namespace model {
 
 extern std::unordered_map<int, int> plot_map; //!< map of plot ids to index
 extern vector<std::unique_ptr<PlottableInterface>>
-  plots; //!< Plot instance container
+  plots;                                      //!< Plot instance container
 
 extern uint64_t plotter_seed; // Stream index used by the plotter
 
@@ -321,11 +321,11 @@ class ProjectionPlot : public PlottableInterface {
   // loop:
   static const int MAX_INTERSECTIONS = 1000000;
 
-  std::array<int, 2> pixels_;       // pixel dimension of resulting image
+  std::array<int, 2> pixels_;              // pixel dimension of resulting image
   double horizontal_field_of_view_ {70.0}; // horiz. f.o.v. in degrees
-  Position camera_position_;        // where camera is
-  Position look_at_;                // point camera is centered looking at
-  Direction up_ {0.0, 0.0, 1.0};    // which way is up
+  Position camera_position_;               // where camera is
+  Position look_at_;             // point camera is centered looking at
+  Direction up_ {0.0, 0.0, 1.0}; // which way is up
 
   // which color IDs should be wireframed. If empty, all cells are wireframed.
   vector<int> wireframe_ids_;

include/openmc/tallies/tally.h

@@ -86,9 +86,9 @@ class Tally {
   //----------------------------------------------------------------------------
   // Major public data members.
 
-  int id_ {C_NONE}; //!< User-defined identifier
+  int id_ {C_NONE};                    //!< User-defined identifier
 
-  std::string name_; //!< User-defined name
+  std::string name_;                   //!< User-defined name
 
   TallyType type_ {TallyType::VOLUME}; //!< e.g. volume, surface current
 
@@ -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'
+    0: 'analog', 1: 'tracklength', 2: 'collision', 3: 'pulse-height'
 }
 _TALLY_TYPES = {
-    0: 'volume', 1: 'mesh-surface', 2: 'surface'
+    0: 'volume', 1: 'mesh-surface', 2: 'surface', 3: 'pulse-height'
 }
 
 
@@ -169,7 +169,7 @@ class Tally(_FortranObjectWithID):
     id : int
         ID of the tally
     estimator: str
-        Estimator type of tally (analog, tracklength, collision)
+        Estimator type of tally (analog, tracklength, collision, pulse-height)
     filters : list
         List of tally filters
     mean : numpy.ndarray

openmc/tallies.py

@@ -33,7 +33,7 @@
 _FILTER_CLASSES = (openmc.Filter, openmc.CrossFilter, openmc.AggregateFilter)
 
 # Valid types of estimators
-ESTIMATOR_TYPES = ['tracklength', 'collision', 'analog']
+ESTIMATOR_TYPES = ['tracklength', 'collision', 'analog', 'pulse-height']
 
 
 class Tally(IDManagerMixin):

src/dagmc.cpp

@@ -609,7 +609,8 @@ std::pair<double, int32_t> DAGCell::distance(
     // into the implicit complement on the other side where no intersection will
     // be found. Treating this as a lost particle is problematic when plotting.
     // Instead, the infinite distance and invalid surface index are returned.
-    if (settings::run_mode == RunMode::PLOTTING) return {INFTY, -1};
+    if (settings::run_mode == RunMode::PLOTTING)
+      return {INFTY, -1};
 
     // the particle should be marked as lost immediately if an intersection
     // isn't found in a volume that is not the implicit complement. In the case
@@ -739,7 +740,8 @@ void check_dagmc_root_univ()
   }
 }
 
-int32_t next_cell(int32_t surf, int32_t curr_cell, int32_t univ) {
+int32_t next_cell(int32_t surf, int32_t curr_cell, int32_t univ)
+{
   auto surfp = dynamic_cast<DAGSurface*>(model::surfaces[surf - 1].get());
   auto cellp = dynamic_cast<DAGCell*>(model::cells[curr_cell].get());
   auto univp = static_cast<DAGUniverse*>(model::universes[univ].get());
@@ -750,11 +752,12 @@ int32_t next_cell(int32_t surf, int32_t curr_cell, int32_t univ) {
     cellp->dagmc_ptr()->entity_by_index(3, cellp->dag_index());
 
   moab::EntityHandle new_vol;
-  moab::ErrorCode rval = cellp->dagmc_ptr()->next_vol(surf_handle, curr_vol, new_vol);
-  if (rval != moab::MB_SUCCESS) return -1;
+  moab::ErrorCode rval =
+    cellp->dagmc_ptr()->next_vol(surf_handle, curr_vol, new_vol);
+  if (rval != moab::MB_SUCCESS)
+    return -1;
 
-  return cellp->dagmc_ptr()->index_by_handle(new_vol) +
-         univp->cell_idx_offset_;
+  return cellp->dagmc_ptr()->index_by_handle(new_vol) + univp->cell_idx_offset_;
 }
 
 } // namespace openmc

src/output.cpp

@@ -5,8 +5,8 @@
 #include <cstring>   // for strlen
 #include <ctime>     // for time, localtime
 #include <fstream>
-#include <iomanip> // for setw, setprecision, put_time
-#include <ios>     // for fixed, scientific, left
+#include <iomanip>   // for setw, setprecision, put_time
+#include <ios>       // for fixed, scientific, left
 #include <iostream>
 #include <sstream>
 #include <unordered_map>
@@ -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,12 @@ 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) {
+      pht_secondary_particles();
+    }
+
     // Enter new particle in particle track file
     if (write_track())
       add_particle_track(*this);
@@ -383,6 +395,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 +407,47 @@ void Particle::event_death()
   }
 }
 
+void Particle::pht_collision_energy()
+{
+  // 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);
+  int index = std::distance(model::pulse_height_cells.begin(), it);
+  // Adds the energy particles lose in a collision to the pulse-height at the
+  // cell index
+  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 the birth cell of the particle
+  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;
+  }
+  // determine index of cell in pulse_height_cells
+  auto it = std::find(model::pulse_height_cells.begin(),
+    model::pulse_height_cells.end(), cell_born());
+  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 +490,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();