Addition of a collision tracking feature (#3417)

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

CMakeLists.txt

@@ -338,6 +338,7 @@ list(APPEND libopenmc_SOURCES
   src/cell.cpp
   src/chain.cpp
   src/cmfd_solver.cpp
+  src/collision_track.cpp
   src/cross_sections.cpp
   src/dagmc.cpp
   src/distribution.cpp

docs/source/io_formats/collision_track.rst

@@ -0,0 +1,46 @@
+.. _io_collision_track:
+
+===========================
+Collision Track File Format
+===========================
+
+When collision tracking is enabled with ``mcpl=false`` (the default), OpenMC
+writes binary data to an HDF5 file named ``collision_track.h5``. The same data
+may also be written after each batch when multiple files are requested
+(``collision_track.N.h5``) or when the run is performed in parallel. The file
+contains the information needed to reconstruct each recorded collision.
+
+The current revision of the collision track file format is 1.0.
+
+**/**
+
+:Attributes:
+            - **filetype** (*char[]*) -- String indicating the type of file.
+              For collision-track files the value is ``"collision_track"``.
+
+:Datasets:
+
+          - **collision_track_bank** (Compound type) -- Collision information
+            for each stored event. Each entry in the dataset corresponds to one
+            collision and contains the following fields:
+
+              - ``r`` (*double[3]*) -- Position of the collision in [cm].
+              - ``u`` (*double[3]*) -- Direction unit vector immediately after the collision.
+              - ``E`` (*double*) -- Incident particle energy before the collision in [eV].
+              - ``dE`` (*double*) -- Energy loss over the collision (:math:`E_\text{before} - E_\text{after}`) in [eV].
+              - ``time`` (*double*) -- Time of the collision in [s].
+              - ``wgt`` (*double*) -- Particle weight at the collision.
+              - ``event_mt`` (*int*) -- ENDF MT number identifying the reaction.
+              - ``delayed_group`` (*int*) -- Delayed neutron group index (non-zero for delayed events).
+              - ``cell_id`` (*int*) -- ID of the cell in which the collision occurred.
+              - ``nuclide_id`` (*int*) -- ZA identifier of the nuclide (ZZZAAAM format).
+              - ``material_id`` (*int*) -- ID of the material containing the collision site.
+              - ``universe_id`` (*int*) -- ID of the universe containing the collision site.
+              - ``n_collision`` (*int*) -- Collision counter for the particle history.
+              - ``particle`` (*int*) -- Particle type (0=neutron, 1=photon, 2=electron, 3=positron).
+              - ``parent_id`` (*int64*) -- Unique ID of the parent particle.
+              - ``progeny_id`` (*int64*) -- Progeny ID of the particle.
+
+In an MPI run, OpenMC writes the combined dataset by gathering collision-track
+entries from all ranks before flushing them to disk, so the final file appears
+as though it were produced serially.

docs/source/io_formats/index.rst

@@ -44,6 +44,7 @@ Output Files
 
    statepoint
    source
+   collision_track
    summary
    properties
    depletion_results

docs/source/io_formats/settings.rst

@@ -20,6 +20,85 @@ source neutrons.
 
   *Default*: None
 
+-----------------------------
+``<collision_track>`` Element
+-----------------------------
+
+The ``<collision_track>`` element indicates to track information about particle
+collisions based on a set of criteria and store these events in a file named
+``collision_track.h5``. This file records details such as the position of the
+interaction, direction of the incoming particle, incident energy and deposited
+energy, weight, time of the interaction, and the delayed neutron group (0 for
+prompt neutrons). Additional information such as the cell ID, material ID,
+universe ID, nuclide ZAID, particle type, and event MT number are also stored.
+Users can specify one or more criterion to filter collisions. If no criteria are
+specified, it defaults to tracking all collisions across the model.
+
+.. warning::
+    Storing all collisions can be very memory intensive. For more targeted
+    tracking, users can employ a variety of parameters such as ``cell_ids``,
+    ``reactions``, ``universe_ids``, ``material_ids``, ``nuclides``, and
+    ``deposited_E_threshold`` to refine the selection of particle interactions
+    to be banked.
+
+This element can contain one or more of the following attributes or
+sub-elements:
+
+  :max_collisions:
+    An integer indicating the maximum number of collisions to be banked per file.
+
+    *Default*: 1000
+
+  :max_collision_track_files:
+    An integer indicating the number of collision_track files to be used.
+
+    *Default*: 1
+
+  :mcpl:
+    An optional boolean to enable MCPL_-format instead of the native HDF5-based
+    format. If activated, the output file name and type is changed to
+    ``collision_track.mcpl``.
+
+    *Default*: false
+
+    .. _MCPL: https://mctools.github.io/mcpl/mcpl.pdf
+
+  :cell_ids:
+    A list of integers representing cell IDs to define specific cells in which
+    collisions are to be banked.
+
+    *Default*: None
+
+  :universe_ids:
+    A list of integers representing the universe IDs to define specific
+    universes in which collisions are to be banked.
+
+    *Default*: None
+
+  :material_ids:
+    A list of integers representing the material IDs to define specific
+    materials in which collisions are to be banked.
+
+    *Default*: None
+
+  :nuclides:
+    A list of strings representing the nuclide, to define specific
+    define specific target nuclide collisions to be banked.
+
+    *Default*: None
+
+  :reactions:
+    A list of integers representing the ENDF-6 format MT numbers or strings
+    (e.g. (n,fission)) to define specific reaction types to be banked.
+
+    *Default*: None
+
+  :deposited_E_threshold:
+    A float defining the minimum deposited energy per collision (in eV) to
+    trigger banking.
+
+    *Default*: 0.0
+
 ----------------------------------
 ``<confidence_intervals>`` Element
 ----------------------------------

docs/source/pythonapi/base.rst

@@ -216,6 +216,9 @@ Post-processing
    :nosignatures:
    :template: myfunction.rst
 
+   openmc.read_collision_track_file
+   openmc.read_collision_track_hdf5
+   openmc.read_collision_track_mcpl
    openmc.voxel_to_vtk
 
 The following classes and functions are used for functional expansion reconstruction.

docs/source/usersguide/settings.rst

@@ -756,6 +756,62 @@ instance, whereas the :meth:`openmc.Track.filter` method returns a new
             track_files = [f"tracks_p{rank}.h5" for rank in range(32)]
             openmc.Tracks.combine(track_files, "tracks.h5")
 
+Collision Track File
+---------------------
+
+OpenMC can generate a collision track file that contains detailed collision
+information (position, direction, energy, deposited energy, time, weight, cell
+ID, material ID, universe ID, nuclide ZAID, particle type, particle delayed
+group and particle ID) for each particle collision depending on user-defined
+parameters. To invoke this feature, set the
+:attr:`~openmc.Settings.collision_track` attribute as shown in this example::
+
+  settings.collision_track = {
+      "max_collisions": 300,
+      "reactions": ["(n,fission)", "(n,2n)"],
+      "material_ids": [1,2],
+      "nuclides": ["U238", "O16"],
+      "cell_ids": [5, 12]
+  }
+
+In this example, collision track information is written to the
+collision_track.h5 file at the end of the simulation. The file contains
+300 recorded collisions that occurred in materials with IDs 1 or 2, involving
+fission or (n,2n) reactions on the nuclides U-238 or O-16, within cells
+with IDs 5 and 12.
+The file can be read using :func:`openmc.read_collision_track_file`.
+The example below shows how to extract the data from the collision_track
+feature and displays the fields stored in the file:
+
+>>> data = openmc.read_collision_track_file('collision_track.h5')
+>>> data.dtype
+    dtype([('r', [('x', '<f8'), ('y', '<f8'), ('z', '<f8')]),
+    ('u', [('x', '<f8'), ('y', '<f8'), ('z', '<f8')]), ('E', '<f8'),
+    ('dE', '<f8'), ('time', '<f8'), ('wgt', '<f8'), ('event_mt', '<i4'),
+    ('delayed_group', '<i4'), ('cell_id', '<i4'), ('nuclide_id', '<i4'),
+    ('material_id', '<i4'), ('universe_id', '<i4'), ('n_collision', '<i4'),
+    ('particle', '<i4'), ('parent_id', '<i8'), ('progeny_id', '<i8')])
+
+
+The full list of fields is as follows:
+
+  :r: Position (each direction in [cm])
+  :u: Direction
+  :E: Energy in [eV]
+  :dE: Energy deposited during collision in [eV]
+  :time: Time in [s]
+  :wgt: Weight of the particle
+  :event_mt: Reaction MT number
+  :delayed_group: Delayed group of the particle
+  :cell_id: Cell ID
+  :nuclide_id: Nuclide ID (10000×Z + 10×A + M)
+  :material_id: Material ID
+  :universe_id: Universe ID
+  :n_collision: Number of collision suffered by the particle
+  :particle: Particle type
+  :parent_id: Source particle ID
+  :progeny_id: Progeny ID
+
 -----------------------
 Restarting a Simulation
 -----------------------

include/openmc/bank.h

@@ -20,6 +20,8 @@ extern vector<SourceSite> source_bank;
 
 extern SharedArray<SourceSite> surf_source_bank;
 
+extern SharedArray<CollisionTrackSite> collision_track_bank;
+
 extern SharedArray<SourceSite> fission_bank;
 
 extern vector<vector<int>> ifp_source_delayed_group_bank;

include/openmc/bank_io.h

@@ -0,0 +1,103 @@
+#ifndef OPENMC_BANK_IO_H
+#define OPENMC_BANK_IO_H
+
+#include "hdf5.h"
+
+#include "openmc/message_passing.h"
+#include "openmc/span.h"
+#include "openmc/vector.h"
+
+#include <algorithm>
+
+#ifdef OPENMC_MPI
+#include <mpi.h>
+#endif
+
+namespace openmc {
+
+template<typename SiteType>
+void write_bank_dataset(const char* dataset_name, hid_t group_id,
+  span<SiteType> bank, const vector<int64_t>& bank_index, hid_t banktype
+#ifdef OPENMC_MPI
+  ,
+  MPI_Datatype mpi_dtype
+#endif
+)
+{
+  int64_t dims_size = bank_index.back();
+  int64_t count_size = bank_index[mpi::rank + 1] - bank_index[mpi::rank];
+
+#ifdef PHDF5
+  hsize_t dims[] {static_cast<hsize_t>(dims_size)};
+  hid_t dspace = H5Screate_simple(1, dims, nullptr);
+  hid_t dset = H5Dcreate(group_id, dataset_name, banktype, dspace, H5P_DEFAULT,
+    H5P_DEFAULT, H5P_DEFAULT);
+
+  hsize_t count[] {static_cast<hsize_t>(count_size)};
+  hid_t memspace = H5Screate_simple(1, count, nullptr);
+
+  hsize_t start[] {static_cast<hsize_t>(bank_index[mpi::rank])};
+  H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
+
+  hid_t plist = H5Pcreate(H5P_DATASET_XFER);
+  H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
+
+  H5Dwrite(dset, banktype, memspace, dspace, plist, bank.data());
+
+  H5Sclose(dspace);
+  H5Sclose(memspace);
+  H5Dclose(dset);
+  H5Pclose(plist);
+#else
+  if (mpi::master) {
+    hsize_t dims[] {static_cast<hsize_t>(dims_size)};
+    hid_t dspace = H5Screate_simple(1, dims, nullptr);
+    hid_t dset = H5Dcreate(group_id, dataset_name, banktype, dspace,
+      H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+
+#ifdef OPENMC_MPI
+    vector<SiteType> temp_bank {bank.begin(), bank.end()};
+#endif
+
+    for (int i = 0; i < mpi::n_procs; ++i) {
+      hsize_t count[] {static_cast<hsize_t>(bank_index[i + 1] - bank_index[i])};
+      hid_t memspace = H5Screate_simple(1, count, nullptr);
+
+#ifdef OPENMC_MPI
+      if (i > 0) {
+        MPI_Recv(bank.data(), count[0], mpi_dtype, i, i, mpi::intracomm,
+          MPI_STATUS_IGNORE);
+      }
+#endif
+
+      hid_t dspace_rank = H5Dget_space(dset);
+      hsize_t start[] {static_cast<hsize_t>(bank_index[i])};
+      H5Sselect_hyperslab(
+        dspace_rank, H5S_SELECT_SET, start, nullptr, count, nullptr);
+
+      H5Dwrite(dset, banktype, memspace, dspace_rank, H5P_DEFAULT, bank.data());
+
+      H5Sclose(memspace);
+      H5Sclose(dspace_rank);
+    }
+
+    H5Dclose(dset);
+
+#ifdef OPENMC_MPI
+    std::copy(temp_bank.begin(), temp_bank.end(), bank.begin());
+#endif
+  }
+#ifdef OPENMC_MPI
+  else {
+    if (!bank.empty()) {
+      MPI_Send(
+        bank.data(), bank.size(), mpi_dtype, 0, mpi::rank, mpi::intracomm);
+    }
+  }
+#endif
+#endif
+}
+
+} // namespace openmc
+
+#endif // OPENMC_BANK_IO_H

include/openmc/collision_track.h

@@ -0,0 +1,23 @@
+#ifndef OPENMC_COLLISION_TRACK_H
+#define OPENMC_COLLISION_TRACK_H
+
+#include <string>
+
+namespace openmc {
+
+class Particle;
+
+//! Reserve space in the collision track bank according to user settings.
+void collision_track_reserve_bank();
+
+//! Write collision track data to disk when the bank is full or the batch ends.
+void collision_track_flush_bank();
+
+//! Record the current particle as a collision-track entry when applicable.
+//!
+//! \param particle Particle whose collision should be recorded if eligible
+void collision_track_record(Particle& particle);
+
+} // namespace openmc
+
+#endif // OPENMC_COLLISION_TRACK_H

include/openmc/constants.h

@@ -34,6 +34,7 @@ constexpr array<int, 2> VERSION_VOXEL {2, 0};
 constexpr array<int, 2> VERSION_MGXS_LIBRARY {1, 0};
 constexpr array<int, 2> VERSION_PROPERTIES {1, 1};
 constexpr array<int, 2> VERSION_WEIGHT_WINDOWS {1, 0};
+constexpr array<int, 2> VERSION_COLLISION_TRACK {1, 0};
 
 // ============================================================================
 // ADJUSTABLE PARAMETERS

include/openmc/mcpl_interface.h

@@ -38,6 +38,21 @@ vector<SourceSite> mcpl_source_sites(std::string path);
 void write_mcpl_source_point(const char* filename, span<SourceSite> source_bank,
   const vector<int64_t>& bank_index);
 
+//! Write an MCPL collision track file
+//!
+//! This function writes collision track data to an MCPL file. Additional
+//! collision-specific metadata (such as energy deposition, material info, etc.)
+//! is stored in the file header as blob data.
+//!
+//! \param[in] filename     Path to MCPL file
+//! \param[in] collision_track_bank  Vector of CollisionTrackSites to write to
+//!                         file for this MPI rank.
+//! \param[in] bank_index   Pointer to vector of site index ranges over all
+//!                         MPI ranks.
+void write_mcpl_collision_track(const char* filename,
+  span<CollisionTrackSite> collision_track_bank,
+  const vector<int64_t>& bank_index);
+
 //! Check if MCPL functionality is available
 bool is_mcpl_interface_available();
 

include/openmc/message_passing.h

@@ -18,6 +18,7 @@ extern bool master;
 
 #ifdef OPENMC_MPI
 extern MPI_Datatype source_site;
+extern MPI_Datatype collision_track_site;
 extern MPI_Comm intracomm;
 #endif
 

include/openmc/particle_data.h

@@ -56,6 +56,25 @@ struct SourceSite {
   int64_t progeny_id;
 };
 
+struct CollisionTrackSite {
+  Position r;
+  Direction u;
+  double E;
+  double dE;
+  double time {0.0};
+  double wgt {1.0};
+  int event_mt {0};
+  int delayed_group {0};
+  int cell_id {0};
+  int nuclide_id;
+  int material_id {0};
+  int universe_id {0};
+  int n_collision {0};
+  ParticleType particle;
+  int64_t parent_id;
+  int64_t progeny_id;
+};
+
 //! State of a particle used for particle track files
 struct TrackState {
   Position r;           //!< Position in [cm]

include/openmc/settings.h

@@ -32,6 +32,24 @@ enum class IFPParameter {
   GenerationTime,
 };
 
+struct CollisionTrackConfig {
+  bool mcpl_write {false}; //!< Write collision tracks using MCPL?
+  std::unordered_set<int>
+    cell_ids; //!< Cell ids where collisions will be written
+  std::unordered_set<int>
+    mt_numbers; //!< MT Numbers where collisions will be written
+  std::unordered_set<int>
+    universe_ids; //!< Universe IDs where collisions will be written
+  std::unordered_set<int>
+    material_ids; //!< Material IDs where collisions will be written
+  std::unordered_set<std::string>
+    nuclides; //!< Nuclides where collisions will be written
+  double deposited_energy_threshold {0.0}; //!< Minimum deposited energy [eV]
+  int64_t max_collisions {
+    1000};               //!< Maximum events recorded per collision track file
+  int64_t max_files {1}; //!< Maximum number of collision track files
+};
+
 //==============================================================================
 // Global variable declarations
 //==============================================================================
@@ -41,6 +59,7 @@ namespace settings {
 // Boolean flags
 extern bool assume_separate;      //!< assume tallies are spatially separate?
 extern bool check_overlaps;       //!< check overlaps in geometry?
+extern bool collision_track;      //!< flag to use collision track feature?
 extern bool confidence_intervals; //!< use confidence intervals for results?
 extern bool
   create_fission_neutrons; //!< create fission neutrons (fixed source)?
@@ -145,6 +164,7 @@ extern std::unordered_set<int>
   statepoint_batch; //!< Batches when state should be written
 extern std::unordered_set<int>
   source_write_surf_id; //!< Surface ids where sources will be written
+extern CollisionTrackConfig collision_track_config;
 extern double source_rejection_fraction; //!< Minimum fraction of source sites
                                          //!< that must be accepted
 extern double free_gas_threshold;        //!< Threshold multiplier for free gas

include/openmc/simulation.h

@@ -22,6 +22,7 @@ constexpr int STATUS_EXIT_ON_TRIGGER {2};
 
 namespace simulation {
 
+extern int ct_current_file;   //!< current collision track file index
 extern "C" int current_batch; //!< current batch
 extern "C" int current_gen;   //!< current fission generation
 extern "C" bool initialized;  //!< has simulation been initialized?

include/openmc/urr.h

@@ -27,10 +27,10 @@ public:
     double heating;
   };
 
-  Interpolation interp_;          //!< interpolation type
-  int inelastic_flag_;            //!< inelastic competition flag
-  int absorption_flag_;           //!< other absorption flag
-  bool multiply_smooth_;          //!< multiply by smooth cross section?
+  Interpolation interp_; //!< interpolation type
+  int inelastic_flag_;   //!< inelastic competition flag
+  int absorption_flag_;  //!< other absorption flag
+  bool multiply_smooth_; //!< multiply by smooth cross section?
 
   vector<double> energy_; //!< incident energies
   auto n_energy() const { return energy_.size(); }

openmc/settings.py

@@ -6,7 +6,7 @@ from numbers import Integral, Real
 from pathlib import Path
 
 import lxml.etree as ET
-
+import warnings
 import openmc
 import openmc.checkvalue as cv
 from openmc.checkvalue import PathLike
@@ -47,6 +47,21 @@ class Settings:
         half-width of the 95% two-sided confidence interval. If False,
         uncertainties on tally results will be reported as the sample standard
         deviation.
+    collision_track : dict
+        Options for writing collision information. Acceptable keys are:
+
+        :max_collisions: Maximum number of collisions to be banked per file. (int)
+        :max_collision_track_files: Maximum number of collision_track files. (int)
+        :mcpl: Output in the form of an MCPL-file. (bool)
+        :cell_ids: List of cell IDs to define cells in which collisions should be banked. (list of int)
+        :universe_ids: List of universe IDs to define universes in which collisions should be banked. (list of int)
+        :material_ids: List of material IDs to define materials in which collisions should be banked. (list of int)
+        :nuclides: List of nuclides to define nuclides in which collisions should be banked.
+                    (ex: ["I135m", "U233"] ). (list of str)
+        :reactions: List of reaction to define specific reactions that should be banked
+                    (ex: ["(n,fission)", 2, "(n,2n)"] ). (list of str or int)
+        :deposited_E_threshold: Number to define the minimum deposited energy during
+                     per collision to trigger banking. (float)
     create_fission_neutrons : bool
         Indicate whether fission neutrons should be created or not.
     cutoff : dict
@@ -395,6 +410,9 @@ class Settings:
         # Iterated Fission Probability
         self._ifp_n_generation = None
 
+        # Collision track feature
+        self._collision_track = {}
+
         # Output options
         self._statepoint = {}
         self._sourcepoint = {}
@@ -434,7 +452,8 @@ class Settings:
         self._max_particle_events = None
         self._write_initial_source = None
         self._weight_windows = WeightWindowsList()
-        self._weight_window_generators = cv.CheckedList(WeightWindowGenerator, 'weight window generators')
+        self._weight_window_generators = cv.CheckedList(
+            WeightWindowGenerator, 'weight window generators')
         self._weight_windows_on = None
         self._weight_windows_file = None
         self._weight_window_checkpoints = {}
@@ -475,8 +494,9 @@ class Settings:
 
     @generations_per_batch.setter
     def generations_per_batch(self, generations_per_batch: int):
-        cv.check_type('generations per patch', generations_per_batch, Integral)
-        cv.check_greater_than('generations per batch', generations_per_batch, 0)
+        cv.check_type('generations per batch', generations_per_batch, Integral)
+        cv.check_greater_than('generations per batch',
+                              generations_per_batch, 0)
         self._generations_per_batch = generations_per_batch
 
     @property
@@ -506,7 +526,8 @@ class Settings:
     @rel_max_lost_particles.setter
     def rel_max_lost_particles(self, rel_max_lost_particles: float):
         cv.check_type('rel_max_lost_particles', rel_max_lost_particles, Real)
-        cv.check_greater_than('rel_max_lost_particles', rel_max_lost_particles, 0)
+        cv.check_greater_than('rel_max_lost_particles',
+                              rel_max_lost_particles, 0)
         cv.check_less_than('rel_max_lost_particles', rel_max_lost_particles, 1)
         self._rel_max_lost_particles = rel_max_lost_particles
 
@@ -516,8 +537,10 @@ class Settings:
 
     @max_write_lost_particles.setter
     def max_write_lost_particles(self, max_write_lost_particles: int):
-        cv.check_type('max_write_lost_particles', max_write_lost_particles, Integral)
-        cv.check_greater_than('max_write_lost_particles', max_write_lost_particles, 0)
+        cv.check_type('max_write_lost_particles',
+                      max_write_lost_particles, Integral)
+        cv.check_greater_than('max_write_lost_particles',
+                              max_write_lost_particles, 0)
         self._max_write_lost_particles = max_write_lost_particles
 
     @property
@@ -538,12 +561,12 @@ class Settings:
     def keff_trigger(self, keff_trigger: dict):
         if not isinstance(keff_trigger, dict):
             msg = f'Unable to set a trigger on keff from "{keff_trigger}" ' \
-                  'which is not a Python dictionary'
+                'which is not a Python dictionary'
             raise ValueError(msg)
 
         elif 'type' not in keff_trigger:
             msg = f'Unable to set a trigger on keff from "{keff_trigger}" ' \
-                  'which does not have a "type" key'
+                'which does not have a "type" key'
             raise ValueError(msg)
 
         elif keff_trigger['type'] not in ['variance', 'std_dev', 'rel_err']:
@@ -553,7 +576,7 @@ class Settings:
 
         elif 'threshold' not in keff_trigger:
             msg = f'Unable to set a trigger on keff from "{keff_trigger}" ' \
-                  'which does not have a "threshold" key'
+                'which does not have a "threshold" key'
             raise ValueError(msg)
 
         elif not isinstance(keff_trigger['threshold'], Real):
@@ -570,7 +593,7 @@ class Settings:
     @energy_mode.setter
     def energy_mode(self, energy_mode: str):
         cv.check_value('energy mode', energy_mode,
-                    ['continuous-energy', 'multi-group'])
+                       ['continuous-energy', 'multi-group'])
         self._energy_mode = energy_mode
 
     @property
@@ -593,7 +616,8 @@ class Settings:
     def source(self, source: SourceBase | Iterable[SourceBase]):
         if not isinstance(source, MutableSequence):
             source = [source]
-        self._source = cv.CheckedList(SourceBase, 'source distributions', source)
+        self._source = cv.CheckedList(
+            SourceBase, 'source distributions', source)
 
     @property
     def confidence_intervals(self) -> bool:
@@ -610,7 +634,8 @@ class Settings:
 
     @electron_treatment.setter
     def electron_treatment(self, electron_treatment: str):
-        cv.check_value('electron treatment', electron_treatment, ['led', 'ttb'])
+        cv.check_value('electron treatment',
+                       electron_treatment, ['led', 'ttb'])
         self._electron_treatment = electron_treatment
 
     @property
@@ -704,7 +729,8 @@ class Settings:
     @trigger_max_batches.setter
     def trigger_max_batches(self, trigger_max_batches: int):
         cv.check_type('trigger maximum batches', trigger_max_batches, Integral)
-        cv.check_greater_than('trigger maximum batches', trigger_max_batches, 0)
+        cv.check_greater_than('trigger maximum batches',
+                              trigger_max_batches, 0)
         self._trigger_max_batches = trigger_max_batches
 
     @property
@@ -713,8 +739,10 @@ class Settings:
 
     @trigger_batch_interval.setter
     def trigger_batch_interval(self, trigger_batch_interval: int):
-        cv.check_type('trigger batch interval', trigger_batch_interval, Integral)
-        cv.check_greater_than('trigger batch interval', trigger_batch_interval, 0)
+        cv.check_type('trigger batch interval',
+                      trigger_batch_interval, Integral)
+        cv.check_greater_than('trigger batch interval',
+                              trigger_batch_interval, 0)
         self._trigger_batch_interval = trigger_batch_interval
 
     @property
@@ -798,19 +826,22 @@ class Settings:
 
     @surf_source_write.setter
     def surf_source_write(self, surf_source_write: dict):
-        cv.check_type("surface source writing options", surf_source_write, Mapping)
+        cv.check_type("surface source writing options",
+                      surf_source_write, Mapping)
         for key, value in surf_source_write.items():
             cv.check_value(
                 "surface source writing key",
                 key,
-                ("surface_ids", "max_particles", "max_source_files", "mcpl", "cell", "cellfrom", "cellto"),
+                ("surface_ids", "max_particles", "max_source_files",
+                 "mcpl", "cell", "cellfrom", "cellto"),
             )
             if key == "surface_ids":
                 cv.check_type(
                     "surface ids for source banking", value, Iterable, Integral
                 )
                 for surf_id in value:
-                    cv.check_greater_than("surface id for source banking", surf_id, 0)
+                    cv.check_greater_than(
+                        "surface id for source banking", surf_id, 0)
 
             elif key == "mcpl":
                 cv.check_type("write to an MCPL-format file", value, bool)
@@ -827,6 +858,79 @@ class Settings:
 
         self._surf_source_write = surf_source_write
 
+    @property
+    def collision_track(self) -> dict:
+        return self._collision_track
+
+    @collision_track.setter
+    def collision_track(self, collision_track: dict):
+        cv.check_type('Collision tracking options', collision_track, Mapping)
+        for key, value in collision_track.items():
+            cv.check_value('collision_track key', key,
+                           ('cell_ids', 'reactions', 'universe_ids', 'material_ids', 'nuclides',
+                            'deposited_E_threshold', 'max_collisions', 'max_collision_track_files', 'mcpl'))
+            if key == 'cell_ids':
+                cv.check_type('cell ids for collision tracking data banking', value,
+                              Iterable, Integral)
+                for cell_id in value:
+                    cv.check_greater_than('cell id for collision tracking data banking',
+                                          cell_id, 0)
+            elif key == 'reactions':
+                cv.check_type('MT numbers for collision tracking data banking', value,
+                              Iterable)
+                for reaction in value:
+                    if isinstance(reaction, int):
+                        cv.check_greater_than(
+                            'MT number for collision tracking data banking', reaction, 0
+                        )
+                    elif isinstance(reaction, str):
+                        # check against allowed strings? so far let C++ code handle it
+                        pass
+                    else:
+                        raise TypeError(
+                            f"MT number for collision tracking data banking must be a positive int or string, "
+                            f"got {type(reaction).__name__}")
+            elif key == 'universe_ids':
+                cv.check_type('universe ids for collision tracking data banking', value,
+                              Iterable, Integral)
+                for universe_id in value:
+                    cv.check_greater_than('universe id for collision tracking data banking',
+                                          universe_id, 0)
+            elif key == 'material_ids':
+                cv.check_type('material ids for collision tracking data banking', value,
+                              Iterable, Integral)
+                for material_id in value:
+                    cv.check_greater_than('material id for collision tracking data banking',
+                                          material_id, 0)
+            elif key == 'nuclides':
+                cv.check_type('nuclides for collision tracking data banking', value,
+                              Iterable, str)
+                for nuclide in value:
+                    # If nuclide name doesn't look valid, give a warning
+                    try:
+                        openmc.data.zam(nuclide)
+                    except ValueError:
+                        warnings.warn(f"Nuclide {nuclide} is not valid")
+            elif key == 'deposited_E_threshold':
+                cv.check_type('Deposited Energy Threshold for collision tracking data banking',
+                              value, Real)
+                cv.check_greater_than('Deposited Energy Threshold for collision tracking data banking',
+                                      value, 0)
+            elif key == 'max_collisions':
+                cv.check_type('maximum collisions banks per file',
+                              value, Integral)
+                cv.check_greater_than('maximum collisions banks in collision tracking',
+                                      value, 0)
+            elif key == 'max_collision_track_files':
+                cv.check_type('maximum collisions banks',
+                              value, Integral)
+                cv.check_greater_than('maximum number of collision_track files ',
+                                      value, 0)
+            elif key == 'mcpl':
+                cv.check_type('write to an MCPL-format file', value, bool)
+
+        self._collision_track = collision_track
+
     @property
     def no_reduce(self) -> bool:
         return self._no_reduce
@@ -943,7 +1047,7 @@ class Settings:
     def cutoff(self, cutoff: dict):
         if not isinstance(cutoff, Mapping):
             msg = f'Unable to set cutoff from "{cutoff}" which is not a '\
-                  'Python dictionary'
+                'Python dictionary'
             raise ValueError(msg)
         for key in cutoff:
             if key == 'weight':
@@ -961,7 +1065,7 @@ class Settings:
                 cv.check_greater_than('energy cutoff', cutoff[key], 0.0)
             else:
                 msg = f'Unable to set cutoff to "{key}" which is unsupported ' \
-                      'by OpenMC'
+                    'by OpenMC'
 
         self._cutoff = cutoff
 
@@ -1130,12 +1234,14 @@ class Settings:
     @weight_window_checkpoints.setter
     def weight_window_checkpoints(self, weight_window_checkpoints: dict):
         for key in weight_window_checkpoints.keys():
-            cv.check_value('weight_window_checkpoints', key, ('collision', 'surface'))
+            cv.check_value('weight_window_checkpoints',
+                           key, ('collision', 'surface'))
         self._weight_window_checkpoints = weight_window_checkpoints
 
     @property
     def max_splits(self):
-        raise AttributeError('max_splits has been deprecated. Please use max_history_splits instead')
+        raise AttributeError(
+            'max_splits has been deprecated. Please use max_history_splits instead')
 
     @property
     def max_history_splits(self) -> int:
@@ -1184,7 +1290,8 @@ class Settings:
     def weight_window_generators(self, wwgs):
         if not isinstance(wwgs, MutableSequence):
             wwgs = [wwgs]
-        self._weight_window_generators = cv.CheckedList(WeightWindowGenerator, 'weight window generators', wwgs)
+        self._weight_window_generators = cv.CheckedList(
+            WeightWindowGenerator, 'weight window generators', wwgs)
 
     @property
     def random_ray(self) -> dict:
@@ -1221,7 +1328,8 @@ class Settings:
                 for mesh, domains in value:
                     cv.check_type('mesh', mesh, MeshBase)
                     cv.check_type('domains', domains, Iterable)
-                    valid_types = (openmc.Material, openmc.Cell, openmc.Universe)
+                    valid_types = (openmc.Material,
+                                   openmc.Cell, openmc.Universe)
                     for domain in domains:
                         if not isinstance(domain, valid_types):
                             raise ValueError(
@@ -1255,9 +1363,12 @@ class Settings:
 
     @source_rejection_fraction.setter
     def source_rejection_fraction(self, source_rejection_fraction: float):
-        cv.check_type('source_rejection_fraction', source_rejection_fraction, Real)
-        cv.check_greater_than('source_rejection_fraction', source_rejection_fraction, 0)
-        cv.check_less_than('source_rejection_fraction', source_rejection_fraction, 1)
+        cv.check_type('source_rejection_fraction',
+                      source_rejection_fraction, Real)
+        cv.check_greater_than('source_rejection_fraction',
+                              source_rejection_fraction, 0)
+        cv.check_less_than('source_rejection_fraction',
+                           source_rejection_fraction, 1)
         self._source_rejection_fraction = source_rejection_fraction
 
     @property
@@ -1422,6 +1533,45 @@ class Settings:
                     subelement = ET.SubElement(element, key)
                     subelement.text = str(self._surf_source_write[key])
 
+    def _create_collision_track_subelement(self, root):
+        if self._collision_track:
+            element = ET.SubElement(root, "collision_track")
+            if 'cell_ids' in self._collision_track:
+                subelement = ET.SubElement(element, "cell_ids")
+                subelement.text = ' '.join(
+                    str(x) for x in self._collision_track['cell_ids'])
+            if 'reactions' in self._collision_track:
+                subelement = ET.SubElement(element, "reactions")
+                subelement.text = ' '.join(
+                    str(x) for x in self._collision_track['reactions'])
+            if 'universe_ids' in self._collision_track:
+                subelement = ET.SubElement(element, "universe_ids")
+                subelement.text = ' '.join(
+                    str(x) for x in self._collision_track['universe_ids'])
+            if 'material_ids' in self._collision_track:
+                subelement = ET.SubElement(element, "material_ids")
+                subelement.text = ' '.join(
+                    str(x) for x in self._collision_track['material_ids'])
+            if 'nuclides' in self._collision_track:
+                subelement = ET.SubElement(element, "nuclides")
+                subelement.text = ' '.join(
+                    str(x) for x in self._collision_track['nuclides'])
+            if 'deposited_E_threshold' in self._collision_track:
+                subelement = ET.SubElement(element, "deposited_E_threshold")
+                subelement.text = str(
+                    self._collision_track['deposited_E_threshold'])
+            if 'max_collisions' in self._collision_track:
+                subelement = ET.SubElement(element, "max_collisions")
+                subelement.text = str(self._collision_track['max_collisions'])
+            if 'max_collision_track_files' in self._collision_track:
+                subelement = ET.SubElement(
+                    element, "max_collision_track_files")
+                subelement.text = str(
+                    self._collision_track['max_collision_track_files'])
+            if 'mcpl' in self._collision_track:
+                subelement = ET.SubElement(element, "mcpl")
+                subelement.text = str(self._collision_track['mcpl']).lower()
+
     def _create_confidence_intervals(self, root):
         if self._confidence_intervals is not None:
             element = ET.SubElement(root, "confidence_intervals")
@@ -1477,8 +1627,8 @@ class Settings:
         # use default heuristic for entropy mesh if not set by user
         if self.entropy_mesh.dimension is None:
             if self.particles is None:
-                raise RuntimeError("Number of particles must be set in order to " \
-                    "use entropy mesh dimension heuristic")
+                raise RuntimeError("Number of particles must be set in order to "
+                                   "use entropy mesh dimension heuristic")
             else:
                 n = ceil((self.particles / 20.0)**(1.0 / 3.0))
                 d = len(self.entropy_mesh.lower_left)
@@ -1568,7 +1718,8 @@ class Settings:
         path = f"./mesh[@id='{self.ufs_mesh.id}']"
         if root.find(path) is None:
             root.append(self.ufs_mesh.to_xml_element())
-            if mesh_memo is not None: mesh_memo.add(self.ufs_mesh.id)
+            if mesh_memo is not None:
+                mesh_memo.add(self.ufs_mesh.id)
 
     def _create_use_decay_photons_subelement(self, root):
         if self._use_decay_photons is not None:
@@ -1694,11 +1845,13 @@ class Settings:
 
         if 'collision' in self._weight_window_checkpoints:
             subelement = ET.SubElement(element, "collision")
-            subelement.text = str(self._weight_window_checkpoints['collision']).lower()
+            subelement.text = str(
+                self._weight_window_checkpoints['collision']).lower()
 
         if 'surface' in self._weight_window_checkpoints:
             subelement = ET.SubElement(element, "surface")
-            subelement.text = str(self._weight_window_checkpoints['surface']).lower()
+            subelement.text = str(
+                self._weight_window_checkpoints['surface']).lower()
 
     def _create_max_history_splits_subelement(self, root):
         if self._max_history_splits is not None:
@@ -1730,6 +1883,9 @@ class Settings:
                         for domain in domains:
                             domain_elem = ET.SubElement(mesh_elem, 'domain')
                             domain_elem.set('id', str(domain.id))
+                            domain_elem.set(
+                                'type', domain.__class__.__name__.lower())
+                        if mesh_memo is not None and mesh.id not in mesh_memo:
                             domain_elem.set('type', domain.__class__.__name__.lower())
                         # See if a <mesh> element already exists -- if not, add it
                         path = f"./mesh[@id='{mesh.id}']"
@@ -1881,6 +2037,25 @@ class Settings:
                     value = int(value)
                 self.surf_source_write[key] = value
 
+    def _collision_track_from_xml_element(self, root):
+        elem = root.find('collision_track')
+        if elem is not None:
+            for key in ('cell_ids', 'reactions', 'universe_ids', 'material_ids', 'nuclides',
+                        'deposited_E_threshold', 'max_collisions', "max_collision_track_files", 'mcpl'):
+                value = get_text(elem, key)
+                if value is not None:
+                    if key in ('cell_ids', 'universe_ids', 'material_ids'):
+                        value = [int(x) for x in value.split()]
+                    elif key in ('reactions', 'nuclides'):
+                        value = value.split()
+                    elif key in ('max_collisions', 'max_collision_track_files'):
+                        value = int(value)
+                    elif key == 'deposited_E_threshold':
+                        value = float(value)
+                    elif key == 'mcpl':
+                        value = value in ('true', '1')
+                    self.collision_track[key] = value
+
     def _confidence_intervals_from_xml_element(self, root):
         text = get_text(root, 'confidence_intervals')
         if text is not None:
@@ -2180,7 +2355,8 @@ class Settings:
                             elif domain_type == 'universe':
                                 domain = openmc.Universe(domain_id)
                             domains.append(domain)
-                        self.random_ray['source_region_meshes'].append((mesh, domains))
+                        self.random_ray['source_region_meshes'].append(
+                            (mesh, domains))
 
     def _use_decay_photons_from_xml_element(self, root):
         text = get_text(root, 'use_decay_photons')
@@ -2223,6 +2399,7 @@ class Settings:
         self._create_sourcepoint_subelement(element)
         self._create_surf_source_read_subelement(element)
         self._create_surf_source_write_subelement(element)
+        self._create_collision_track_subelement(element)
         self._create_confidence_intervals(element)
         self._create_electron_treatment_subelement(element)
         self._create_energy_mode_subelement(element)
@@ -2336,6 +2513,7 @@ class Settings:
         settings._sourcepoint_from_xml_element(elem)
         settings._surf_source_read_from_xml_element(elem)
         settings._surf_source_write_from_xml_element(elem)
+        settings._collision_track_from_xml_element(elem)
         settings._confidence_intervals_from_xml_element(elem)
         settings._electron_treatment_from_xml_element(elem)
         settings._energy_mode_from_xml_element(elem)

openmc/source.py

@@ -6,7 +6,6 @@ from numbers import Real
 from pathlib import Path
 import warnings
 from typing import Any
-from pathlib import Path
 
 import lxml.etree as ET
 import numpy as np
@@ -107,10 +106,12 @@ class SourceBase(ABC):
                 cv.check_type('fissionable', value, bool)
                 self._constraints['fissionable'] = value
             elif key == 'rejection_strategy':
-                cv.check_value('rejection strategy', value, ('resample', 'kill'))
+                cv.check_value('rejection strategy',
+                               value, ('resample', 'kill'))
                 self._constraints['rejection_strategy'] = value
             else:
-                raise ValueError(f'Unknown key in constraints dictionary: {key}')
+                raise ValueError(
+                    f'Unknown key in constraints dictionary: {key}')
 
     @abstractmethod
     def populate_xml_element(self, element):
@@ -144,13 +145,16 @@ class SourceBase(ABC):
                 dt_elem = ET.SubElement(constraints_elem, "domain_type")
                 dt_elem.text = constraints["domain_type"]
                 id_elem = ET.SubElement(constraints_elem, "domain_ids")
-                id_elem.text = ' '.join(str(uid) for uid in constraints["domain_ids"])
+                id_elem.text = ' '.join(str(uid)
+                                        for uid in constraints["domain_ids"])
             if "time_bounds" in constraints:
                 dt_elem = ET.SubElement(constraints_elem, "time_bounds")
-                dt_elem.text = ' '.join(str(t) for t in constraints["time_bounds"])
+                dt_elem.text = ' '.join(str(t)
+                                        for t in constraints["time_bounds"])
             if "energy_bounds" in constraints:
                 dt_elem = ET.SubElement(constraints_elem, "energy_bounds")
-                dt_elem.text = ' '.join(str(E) for E in constraints["energy_bounds"])
+                dt_elem.text = ' '.join(str(E)
+                                        for E in constraints["energy_bounds"])
             if "fissionable" in constraints:
                 dt_elem = ET.SubElement(constraints_elem, "fissionable")
                 dt_elem.text = str(constraints["fissionable"]).lower()
@@ -199,7 +203,8 @@ class SourceBase(ABC):
             elif source_type == 'mesh':
                 return MeshSource.from_xml_element(elem, meshes)
             else:
-                raise ValueError(f'Source type {source_type} is not recognized')
+                raise ValueError(
+                    f'Source type {source_type} is not recognized')
 
     @staticmethod
     def _get_constraints(elem: ET.Element) -> dict[str, Any]:
@@ -316,7 +321,8 @@ class IndependentSource(SourceBase):
         time: openmc.stats.Univariate | None = None,
         strength: float = 1.0,
         particle: str = 'neutron',
-        domains: Sequence[openmc.Cell | openmc.Material | openmc.Universe] | None = None,
+        domains: Sequence[openmc.Cell | openmc.Material |
+                          openmc.Universe] | None = None,
         constraints: dict[str, Any] | None = None
     ):
         if domains is not None:
@@ -527,11 +533,12 @@ class MeshSource(SourceBase):
         'fissionable', and 'rejection_strategy'.
 
     """
+
     def __init__(
             self,
             mesh: MeshBase,
             sources: Sequence[SourceBase],
-            constraints: dict[str, Any] | None  = None,
+            constraints: dict[str, Any] | None = None,
     ):
         super().__init__(strength=None, constraints=constraints)
         self.mesh = mesh
@@ -577,7 +584,8 @@ class MeshSource(SourceBase):
 
         elif isinstance(self.mesh, UnstructuredMesh):
             if s.ndim > 1:
-                raise ValueError('Sources must be a 1-D array for unstructured mesh')
+                raise ValueError(
+                    'Sources must be a 1-D array for unstructured mesh')
 
         self._sources = s
         for src in self._sources:
@@ -646,7 +654,8 @@ class MeshSource(SourceBase):
         mesh_id = int(get_text(elem, 'mesh'))
         mesh = meshes[mesh_id]
 
-        sources = [SourceBase.from_xml_element(e) for e in elem.iterchildren('source')]
+        sources = [SourceBase.from_xml_element(
+            e) for e in elem.iterchildren('source')]
         constraints = cls._get_constraints(elem)
         return cls(mesh, sources, constraints=constraints)
 
@@ -656,7 +665,8 @@ def Source(*args, **kwargs):
     A function for backward compatibility of sources. Will be removed in the
     future. Please update to IndependentSource.
     """
-    warnings.warn("This class is deprecated in favor of 'IndependentSource'", FutureWarning)
+    warnings.warn(
+        "This class is deprecated in favor of 'IndependentSource'", FutureWarning)
     return openmc.IndependentSource(*args, **kwargs)
 
 
@@ -703,6 +713,7 @@ class CompiledSource(SourceBase):
         'fissionable', and 'rejection_strategy'.
 
     """
+
     def __init__(
         self,
         library: PathLike,
@@ -914,7 +925,8 @@ class ParticleType(IntEnum):
         try:
             return cls[value.upper()]
         except KeyError:
-            raise ValueError(f"Invalid string for creation of {cls.__name__}: {value}")
+            raise ValueError(
+                f"Invalid string for creation of {cls.__name__}: {value}")
 
     @classmethod
     def from_pdg_number(cls, pdg_number: int) -> ParticleType:
@@ -984,6 +996,7 @@ class SourceParticle:
         Type of the particle
 
     """
+
     def __init__(
         self,
         r: Iterable[float] = (0., 0., 0.),
@@ -1042,7 +1055,8 @@ def write_source_file(
     openmc.SourceParticle
 
     """
-    cv.check_iterable_type("source particles", source_particles, SourceParticle)
+    cv.check_iterable_type(
+        "source particles", source_particles, SourceParticle)
     pl = ParticleList(source_particles)
     pl.export_to_hdf5(filename, **kwargs)
 
@@ -1104,7 +1118,8 @@ class ParticleList(list):
             for particle in f.particles:
                 # Determine particle type based on the PDG number
                 try:
-                    particle_type = ParticleType.from_pdg_number(particle.pdgcode)
+                    particle_type = ParticleType.from_pdg_number(
+                        particle.pdgcode)
                 except ValueError:
                     particle_type = "UNKNOWN"
 
@@ -1241,3 +1256,133 @@ def read_source_file(filename: PathLike) -> ParticleList:
         return ParticleList.from_hdf5(filename)
     else:
         return ParticleList.from_mcpl(filename)
+
+
+def read_collision_track_hdf5(filename):
+    """Read a collision track file in HDF5 format.
+
+    Parameters
+    ----------
+    filename : str or path-like
+        Path to the HDF5 collision track file.
+
+    Returns
+    -------
+    numpy.ndarray
+        Structured array containing collision track data.
+
+    See Also
+    --------
+    read_collision_track_mcpl
+    read_collision_track_file
+    """
+
+    with h5py.File(filename, 'r') as file:
+        data = file['collision_track_bank'][:]
+
+    return data
+
+
+def read_collision_track_mcpl(file_path):
+    """Read a collision track file in MCPL format.
+
+    Parameters
+    ----------
+    file_path : str or path-like
+        Path to the MCPL collision track file.
+
+    Returns
+    -------
+    numpy.ndarray
+        Structured array of particle collision track information, including
+        position, direction, energy, weight, reaction data, and identifiers.
+
+    See Also
+    --------
+    read_collision_track_hdf5
+    read_collision_track_file
+    """
+    import mcpl
+    myfile = mcpl.MCPLFile(file_path)
+    data = {
+        'r': [],  # for position (x, y, z)
+        'u': [],  # for direction (ux, uy, uz)
+        'E': [], 'dE': [], 'time': [],
+        'wgt': [], 'event_mt': [], 'delayed_group': [],
+        'cell_id': [], 'nuclide_id': [], 'material_id': [],
+        'universe_id': [], 'n_collision': [], 'particle': [],
+        'parent_id': [], 'progeny_id': []
+    }
+
+    # Read and collect data from the MCPL file
+    for i, p in enumerate(myfile.particles):
+        if f'blob_{i}' in myfile.blobs:
+            blob_data = myfile.blobs[f'blob_{i}']
+            decoded_str = blob_data.decode('utf-8')
+            pairs = decoded_str.split(';')
+            values_dict = {k.strip(): v.strip()
+                           for k, v in (pair.split(':') for pair in pairs if pair.strip())}
+
+            data['r'].append((p.x, p.y, p.z))  # Append as tuple
+            data['u'].append((p.ux, p.uy, p.uz))  # Append as tuple
+            data['E'].append(p.ekin * 1e6)
+            data['dE'].append(float(values_dict.get('dE', 0)))
+            data['time'].append(p.time * 1e-3)
+            data['wgt'].append(p.weight)
+            data['event_mt'].append(int(values_dict.get('event_mt', 0)))
+            data['delayed_group'].append(
+                int(values_dict.get('delayed_group', 0)))
+            data['cell_id'].append(int(values_dict.get('cell_id', 0)))
+            data['nuclide_id'].append(int(values_dict.get('nuclide_id', 0)))
+            data['material_id'].append(int(values_dict.get('material_id', 0)))
+            data['universe_id'].append(int(values_dict.get('universe_id', 0)))
+            data['n_collision'].append(int(values_dict.get('n_collision', 0)))
+            data['particle'].append(ParticleType.from_pdg_number(p.pdgcode))
+            data['parent_id'].append(int(values_dict.get('parent_id', 0)))
+            data['progeny_id'].append(int(values_dict.get('progeny_id', 0)))
+
+    dtypes = [
+        ('r', [('x', 'f8'), ('y', 'f8'), ('z', 'f8')]),
+        ('u', [('x', 'f8'), ('y', 'f8'), ('z', 'f8')]),
+        ('E', 'f8'), ('dE', 'f8'), ('time', 'f8'), ('wgt', 'f8'),
+        ('event_mt', 'f8'), ('delayed_group', 'i4'), ('cell_id', 'i4'),
+        ('nuclide_id', 'i4'), ('material_id', 'i4'), ('universe_id', 'i4'),
+        ('n_collision', 'i4'), ('particle', 'i4'),
+        ('parent_id', 'i8'), ('progeny_id', 'i8')
+    ]
+
+    structured_array = np.zeros(len(data['r']), dtype=dtypes)
+    for key in data:
+        structured_array[key] = data[key]  # Assign data
+
+    return structured_array
+
+
+def read_collision_track_file(filename):
+    """Read a collision track file (HDF5 or MCPL) and return its data.
+
+    Parameters
+    ----------
+    filename : str or path-like
+        Path to the collision track file to read. Must end with
+        ``.h5`` or ``.mcpl``.
+
+    Returns
+    -------
+    numpy.ndarray
+        Structured array containing collision track data.
+
+    See Also
+    --------
+    read_collision_track_hdf5
+    read_collision_track_mcpl
+    """
+
+    filename = Path(filename)
+    if filename.suffix not in ('.h5', '.mcpl'):
+        raise ValueError('Collision track file must have a .h5 or .mcpl extension.')
+
+    if filename.suffix == '.h5':
+        return read_collision_track_hdf5(filename)
+    else:
+        return read_collision_track_mcpl(filename)

src/bank.cpp

@@ -20,6 +20,8 @@ vector<SourceSite> source_bank;
 
 SharedArray<SourceSite> surf_source_bank;
 
+SharedArray<CollisionTrackSite> collision_track_bank;
+
 // The fission bank is allocated as a SharedArray, rather than a vector, as it
 // will be shared by all threads in the simulation. It will be allocated to a
 // fixed maximum capacity in the init_fission_bank() function. Then, Elements
@@ -50,6 +52,7 @@ void free_memory_bank()
 {
   simulation::source_bank.clear();
   simulation::surf_source_bank.clear();
+  simulation::collision_track_bank.clear();
   simulation::fission_bank.clear();
   simulation::progeny_per_particle.clear();
   simulation::ifp_source_delayed_group_bank.clear();

src/collision_track.cpp

@@ -0,0 +1,238 @@
+#include "openmc/collision_track.h"
+
+#include <algorithm>
+#include <string>
+
+#include <fmt/format.h>
+
+#include "openmc/bank.h"
+#include "openmc/bank_io.h"
+#include "openmc/cell.h"
+#include "openmc/constants.h"
+#include "openmc/error.h"
+#include "openmc/file_utils.h"
+#include "openmc/hdf5_interface.h"
+#include "openmc/material.h"
+#include "openmc/mcpl_interface.h"
+#include "openmc/message_passing.h"
+#include "openmc/nuclide.h"
+#include "openmc/output.h"
+#include "openmc/particle.h"
+#include "openmc/settings.h"
+#include "openmc/simulation.h"
+#include "openmc/universe.h"
+
+#ifdef OPENMC_MPI
+#include <mpi.h>
+#endif
+
+namespace openmc {
+
+namespace {
+
+hid_t h5_collision_track_banktype()
+{
+  hid_t postype = H5Tcreate(H5T_COMPOUND, sizeof(Position));
+  H5Tinsert(postype, "x", HOFFSET(Position, x), H5T_NATIVE_DOUBLE);
+  H5Tinsert(postype, "y", HOFFSET(Position, y), H5T_NATIVE_DOUBLE);
+  H5Tinsert(postype, "z", HOFFSET(Position, z), H5T_NATIVE_DOUBLE);
+
+  hid_t banktype = H5Tcreate(H5T_COMPOUND, sizeof(CollisionTrackSite));
+
+  H5Tinsert(banktype, "r", HOFFSET(CollisionTrackSite, r), postype);
+  H5Tinsert(banktype, "u", HOFFSET(CollisionTrackSite, u), postype);
+  H5Tinsert(banktype, "E", HOFFSET(CollisionTrackSite, E), H5T_NATIVE_DOUBLE);
+  H5Tinsert(banktype, "dE", HOFFSET(CollisionTrackSite, dE), H5T_NATIVE_DOUBLE);
+  H5Tinsert(
+    banktype, "time", HOFFSET(CollisionTrackSite, time), H5T_NATIVE_DOUBLE);
+  H5Tinsert(
+    banktype, "wgt", HOFFSET(CollisionTrackSite, wgt), H5T_NATIVE_DOUBLE);
+  H5Tinsert(banktype, "event_mt", HOFFSET(CollisionTrackSite, event_mt),
+    H5T_NATIVE_INT);
+  H5Tinsert(banktype, "delayed_group",
+    HOFFSET(CollisionTrackSite, delayed_group), H5T_NATIVE_INT);
+  H5Tinsert(
+    banktype, "cell_id", HOFFSET(CollisionTrackSite, cell_id), H5T_NATIVE_INT);
+  H5Tinsert(banktype, "nuclide_id", HOFFSET(CollisionTrackSite, nuclide_id),
+    H5T_NATIVE_INT);
+  H5Tinsert(banktype, "material_id", HOFFSET(CollisionTrackSite, material_id),
+    H5T_NATIVE_INT);
+  H5Tinsert(banktype, "universe_id", HOFFSET(CollisionTrackSite, universe_id),
+    H5T_NATIVE_INT);
+  H5Tinsert(banktype, "n_collision", HOFFSET(CollisionTrackSite, n_collision),
+    H5T_NATIVE_INT);
+  H5Tinsert(banktype, "particle", HOFFSET(CollisionTrackSite, particle),
+    H5T_NATIVE_INT);
+  H5Tinsert(banktype, "parent_id", HOFFSET(CollisionTrackSite, parent_id),
+    H5T_NATIVE_INT64);
+  H5Tinsert(banktype, "progeny_id", HOFFSET(CollisionTrackSite, progeny_id),
+    H5T_NATIVE_INT64);
+  H5Tclose(postype);
+  return banktype;
+}
+
+void write_collision_track_bank(hid_t group_id,
+  openmc::span<CollisionTrackSite> collision_track_bank,
+  const openmc::vector<int64_t>& bank_index)
+{
+  hid_t banktype = h5_collision_track_banktype();
+#ifdef OPENMC_MPI
+  write_bank_dataset("collision_track_bank", group_id, collision_track_bank,
+    bank_index, banktype, mpi::collision_track_site);
+#else
+  write_bank_dataset("collision_track_bank", group_id, collision_track_bank,
+    bank_index, banktype);
+#endif
+
+  H5Tclose(banktype);
+}
+
+void write_h5_collision_track(const char* filename,
+  openmc::span<CollisionTrackSite> collision_track_bank,
+  const openmc::vector<int64_t>& bank_index)
+{
+#ifdef PHDF5
+  bool parallel = true;
+#else
+  bool parallel = false;
+#endif
+
+  if (!filename)
+    fatal_error("write_h5_collision_track filename needs a nonempty name.");
+
+  std::string filename_(filename);
+  const auto extension = get_file_extension(filename_);
+  if (extension.empty()) {
+    filename_.append(".h5");
+  } else if (extension != "h5") {
+    warning("write_h5_collision_track was passed a file extension differing "
+            "from .h5, but an hdf5 file will be written.");
+  }
+
+  hid_t file_id;
+  if (mpi::master || parallel) {
+    file_id = file_open(filename_.c_str(), 'w', true);
+
+    // Write filetype and version info
+    write_attribute(file_id, "filetype", "collision_track");
+    write_attribute(file_id, "version", VERSION_COLLISION_TRACK);
+  }
+
+  write_collision_track_bank(file_id, collision_track_bank, bank_index);
+
+  if (mpi::master || parallel)
+    file_close(file_id);
+}
+
+} // namespace
+
+bool should_record_event(int id_cell, int mt_event, const std::string& nuclide,
+  int id_universe, int id_material, double energy_loss)
+{
+  auto matches_filter = [](const auto& filter_set, const auto& value) {
+    return filter_set.empty() || filter_set.count(value) > 0;
+  };
+
+  const auto& cfg = settings::collision_track_config;
+  return simulation::current_batch > settings::n_inactive &&
+         !simulation::collision_track_bank.full() &&
+         matches_filter(cfg.cell_ids, id_cell) &&
+         matches_filter(cfg.mt_numbers, mt_event) &&
+         matches_filter(cfg.universe_ids, id_universe) &&
+         matches_filter(cfg.material_ids, id_material) &&
+         matches_filter(cfg.nuclides, nuclide) &&
+         (cfg.deposited_energy_threshold == 0 ||
+           cfg.deposited_energy_threshold < energy_loss);
+}
+
+void collision_track_reserve_bank()
+{
+  simulation::collision_track_bank.reserve(
+    settings::collision_track_config.max_collisions);
+}
+
+void collision_track_flush_bank()
+{
+  const auto& cfg = settings::collision_track_config;
+  if (simulation::ct_current_file > cfg.max_files)
+    return;
+
+  bool last_batch = (simulation::current_batch == settings::n_batches);
+  if (!simulation::collision_track_bank.full() && !last_batch)
+    return;
+
+  auto size = simulation::collision_track_bank.size();
+  if (size == 0 && !last_batch)
+    return;
+
+  auto collision_track_work_index = mpi::calculate_parallel_index_vector(size);
+  openmc::span<CollisionTrackSite> collisiontrackbankspan(
+    simulation::collision_track_bank.begin(), size);
+
+  std::string ext = cfg.mcpl_write ? "mcpl" : "h5";
+  auto filename = fmt::format("{}collision_track.{}.{}", settings::path_output,
+    simulation::ct_current_file, ext);
+
+  if (cfg.max_files == 1 || (simulation::ct_current_file == 1 && last_batch)) {
+    filename = settings::path_output + "collision_track." + ext;
+  }
+  write_message("Creating {}...", filename, 4);
+
+  if (cfg.mcpl_write) {
+    write_mcpl_collision_track(
+      filename.c_str(), collisiontrackbankspan, collision_track_work_index);
+  } else {
+    write_h5_collision_track(
+      filename.c_str(), collisiontrackbankspan, collision_track_work_index);
+  }
+
+  simulation::collision_track_bank.clear();
+  if (!last_batch && cfg.max_files >= 1) {
+    collision_track_reserve_bank();
+  }
+  ++simulation::ct_current_file;
+}
+
+void collision_track_record(Particle& particle)
+{
+  int cell_index = particle.lowest_coord().cell();
+  if (cell_index == C_NONE)
+    return;
+
+  int cell_id = model::cells[cell_index]->id_;
+  const auto* nuclide_ptr = data::nuclides[particle.event_nuclide()].get();
+  std::string nuclide = nuclide_ptr->name_;
+  int universe_id = model::universes[particle.lowest_coord().universe()]->id_;
+  double delta_E = particle.E_last() - particle.E();
+  int material_index = particle.material();
+  if (material_index == C_NONE)
+    return;
+
+  int material_id = model::materials[material_index]->id_;
+
+  if (!should_record_event(cell_id, particle.event_mt(), nuclide, universe_id,
+        material_id, delta_E))
+    return;
+
+  CollisionTrackSite site;
+  site.r = particle.r();
+  site.u = particle.u();
+  site.E = particle.E_last();
+  site.dE = delta_E;
+  site.time = particle.time();
+  site.wgt = particle.wgt();
+  site.event_mt = particle.event_mt();
+  site.delayed_group = particle.delayed_group();
+  site.cell_id = cell_id;
+  site.nuclide_id =
+    10000 * nuclide_ptr->Z_ + 10 * nuclide_ptr->A_ + nuclide_ptr->metastable_;
+  site.material_id = material_id;
+  site.universe_id = universe_id;
+  site.n_collision = particle.n_collision();
+  site.particle = particle.type();
+  site.parent_id = particle.id();
+  site.progeny_id = particle.n_progeny();
+  simulation::collision_track_bank.thread_safe_append(site);
+}
+
+} // namespace openmc

src/finalize.cpp

@@ -3,6 +3,7 @@
 #include "openmc/bank.h"
 #include "openmc/capi.h"
 #include "openmc/cmfd_solver.h"
+#include "openmc/collision_track.h"
 #include "openmc/constants.h"
 #include "openmc/cross_sections.h"
 #include "openmc/dagmc.h"
@@ -76,6 +77,7 @@ int openmc_finalize()
   // Reset global variables
   settings::assume_separate = false;
   settings::check_overlaps = false;
+  settings::collision_track_config = CollisionTrackConfig {};
   settings::confidence_intervals = false;
   settings::create_fission_neutrons = true;
   settings::create_delayed_neutrons = true;
@@ -177,6 +179,9 @@ int openmc_finalize()
   if (mpi::source_site != MPI_DATATYPE_NULL) {
     MPI_Type_free(&mpi::source_site);
   }
+  if (mpi::collision_track_site != MPI_DATATYPE_NULL) {
+    MPI_Type_free(&mpi::collision_track_site);
+  }
 #endif
 
   openmc_reset_random_ray();

src/initialize.cpp

@@ -178,6 +178,37 @@ void initialize_mpi(MPI_Comm intracomm)
     MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_LONG, MPI_LONG};
   MPI_Type_create_struct(11, blocks, disp, types, &mpi::source_site);
   MPI_Type_commit(&mpi::source_site);
+
+  CollisionTrackSite bc;
+  MPI_Aint dispc[16];
+  MPI_Get_address(&bc.r, &dispc[0]);             // double
+  MPI_Get_address(&bc.u, &dispc[1]);             // double
+  MPI_Get_address(&bc.E, &dispc[2]);             // double
+  MPI_Get_address(&bc.dE, &dispc[3]);            // double
+  MPI_Get_address(&bc.time, &dispc[4]);          // double
+  MPI_Get_address(&bc.wgt, &dispc[5]);           // double
+  MPI_Get_address(&bc.event_mt, &dispc[6]);      // int
+  MPI_Get_address(&bc.delayed_group, &dispc[7]); // int
+  MPI_Get_address(&bc.cell_id, &dispc[8]);       // int
+  MPI_Get_address(&bc.nuclide_id, &dispc[9]);    // int
+  MPI_Get_address(&bc.material_id, &dispc[10]);  // int
+  MPI_Get_address(&bc.universe_id, &dispc[11]);  // int
+  MPI_Get_address(&bc.n_collision, &dispc[12]);  // int
+  MPI_Get_address(&bc.particle, &dispc[13]);     // int
+  MPI_Get_address(&bc.parent_id, &dispc[14]);    // int64_t
+  MPI_Get_address(&bc.progeny_id, &dispc[15]);   // int64_t
+  for (int i = 15; i >= 0; --i) {
+    dispc[i] -= dispc[0];
+  }
+
+  int blocksc[] = {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
+  MPI_Datatype typesc[] = {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE,
+    MPI_DOUBLE, MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_INT,
+    MPI_INT, MPI_INT, MPI_INT, MPI_INT64_T, MPI_INT64_T};
+
+  MPI_Type_create_struct(
+    16, blocksc, dispc, typesc, &mpi::collision_track_site);
+  MPI_Type_commit(&mpi::collision_track_site);
 }
 #endif // OPENMC_MPI
 

src/mcpl_interface.cpp

@@ -18,6 +18,7 @@
 #include <cstring>
 #include <memory>
 #include <mutex>
+#include <sstream>
 #include <stdexcept>
 #include <string>
 
@@ -61,6 +62,8 @@ using mcpl_hdr_nparticles_fpt = uint64_t (*)(mcpl_file_t* file_handle);
 using mcpl_read_fpt = const mcpl_particle_repr_t* (*)(mcpl_file_t* file_handle);
 using mcpl_close_file_fpt = void (*)(mcpl_file_t* file_handle);
 
+using mcpl_hdr_add_data_fpt = void (*)(mcpl_outfile_t* file_handle,
+  const char* key, int32_t ldata, const char* data);
 using mcpl_create_outfile_fpt = mcpl_outfile_t* (*)(const char* filename);
 using mcpl_hdr_set_srcname_fpt = void (*)(
   mcpl_outfile_t* outfile_handle, const char* srcname);
@@ -110,6 +113,7 @@ struct McplApi {
   mcpl_close_file_fpt close_file;
   mcpl_create_outfile_fpt create_outfile;
   mcpl_hdr_set_srcname_fpt hdr_set_srcname;
+  mcpl_hdr_add_data_fpt hdr_add_data;
   mcpl_add_particle_fpt add_particle;
   mcpl_close_outfile_fpt close_outfile;
   mcpl_hdr_add_stat_sum_fpt hdr_add_stat_sum;
@@ -146,6 +150,8 @@ struct McplApi {
       load_symbol_platform("mcpl_create_outfile"));
     hdr_set_srcname = reinterpret_cast<mcpl_hdr_set_srcname_fpt>(
       load_symbol_platform("mcpl_hdr_set_srcname"));
+    hdr_add_data = reinterpret_cast<mcpl_hdr_add_data_fpt>(
+      load_symbol_platform("mcpl_hdr_add_data"));
     add_particle = reinterpret_cast<mcpl_add_particle_fpt>(
       load_symbol_platform("mcpl_add_particle"));
     close_outfile = reinterpret_cast<mcpl_close_outfile_fpt>(
@@ -545,4 +551,153 @@ void write_mcpl_source_point(const char* filename, span<SourceSite> source_bank,
   }
 }
 
+// Collision track feature with MCPL
+void write_mcpl_collision_track_internal(mcpl_outfile_t* file_id,
+  span<CollisionTrackSite> collision_track_bank,
+  const vector<int64_t>& bank_index_all_ranks)
+{
+  if (mpi::master) {
+    if (!file_id) {
+      fatal_error("MCPL: Internal error - master rank called "
+                  "write_mcpl_source_bank_internal with null file_id.");
+    }
+    vector<CollisionTrackSite> receive_buffer;
+    vector<CollisionTrackSite> all_sites;
+    all_sites.reserve(static_cast<size_t>(bank_index_all_ranks.back()));
+    vector<std::string> all_blobs;
+    all_blobs.reserve(static_cast<size_t>(bank_index_all_ranks.back()));
+
+    for (int rank_idx = 0; rank_idx < mpi::n_procs; ++rank_idx) {
+      size_t num_sites_on_rank = static_cast<size_t>(
+        bank_index_all_ranks[rank_idx + 1] - bank_index_all_ranks[rank_idx]);
+      if (num_sites_on_rank == 0)
+        continue;
+
+      span<const CollisionTrackSite> sites_to_process;
+#ifdef OPENMC_MPI
+      if (rank_idx == mpi::rank) {
+        sites_to_process = openmc::span<const CollisionTrackSite>(
+          collision_track_bank.data(), num_sites_on_rank);
+      } else {
+        receive_buffer.resize(num_sites_on_rank);
+        MPI_Recv(receive_buffer.data(), num_sites_on_rank,
+          mpi::collision_track_site, rank_idx, rank_idx, mpi::intracomm,
+          MPI_STATUS_IGNORE);
+        sites_to_process = openmc::span<const CollisionTrackSite>(
+          receive_buffer.data(), num_sites_on_rank);
+      }
+#else
+      sites_to_process = openmc::span<const CollisionTrackSite>(
+        collision_track_bank.data(), num_sites_on_rank);
+#endif
+
+      for (const auto& site : sites_to_process) {
+        std::ostringstream custom_data_stream;
+        custom_data_stream << " dE : " << site.dE
+                           << " ; event_mt : " << site.event_mt
+                           << " ; delayed_group : " << site.delayed_group
+                           << " ; cell_id : " << site.cell_id
+                           << " ; nuclide_id : " << site.nuclide_id
+                           << " ; material_id : " << site.material_id
+                           << " ; universe_id : " << site.universe_id
+                           << " ; n_collision : " << site.n_collision
+                           << " ; parent_id : " << site.parent_id
+                           << " ; progeny_id : " << site.progeny_id;
+
+        all_blobs.push_back(custom_data_stream.str());
+        all_sites.push_back(site);
+      }
+    }
+
+    for (size_t idx = 0; idx < all_blobs.size(); ++idx) {
+      const auto& blob = all_blobs[idx];
+      std::string key = "blob_" + std::to_string(idx);
+      g_mcpl_api->hdr_add_data(file_id, key.c_str(), blob.size(), blob.c_str());
+    }
+
+    for (const auto& site : all_sites) {
+      mcpl_particle_repr_t p_repr {};
+      p_repr.position[0] = site.r.x;
+      p_repr.position[1] = site.r.y;
+      p_repr.position[2] = site.r.z;
+      p_repr.direction[0] = site.u.x;
+      p_repr.direction[1] = site.u.y;
+      p_repr.direction[2] = site.u.z;
+      p_repr.ekin = site.E * 1e-6;
+      p_repr.time = site.time * 1e3;
+      p_repr.weight = site.wgt;
+      switch (site.particle) {
+      case ParticleType::neutron:
+        p_repr.pdgcode = 2112;
+        break;
+      case ParticleType::photon:
+        p_repr.pdgcode = 22;
+        break;
+      case ParticleType::electron:
+        p_repr.pdgcode = 11;
+        break;
+      case ParticleType::positron:
+        p_repr.pdgcode = -11;
+        break;
+      default:
+        continue;
+      }
+      g_mcpl_api->add_particle(file_id, &p_repr);
+    }
+  } else {
+#ifdef OPENMC_MPI
+    if (!collision_track_bank.empty()) {
+      MPI_Send(collision_track_bank.data(), collision_track_bank.size(),
+        mpi::collision_track_site, 0, mpi::rank, mpi::intracomm);
+    }
+#endif
+  }
+}
+
+void write_mcpl_collision_track(const char* filename,
+  span<CollisionTrackSite> collision_track_bank,
+  const vector<int64_t>& bank_index)
+{
+  ensure_mcpl_ready_or_fatal();
+
+  std::string filename_(filename);
+  const auto extension = get_file_extension(filename_);
+  if (extension.empty()) {
+    filename_.append(".mcpl");
+  } else if (extension != "mcpl") {
+    warning(fmt::format("Specified filename '{}' has an extension '.{}', but "
+                        "an MCPL file (.mcpl) will be written using this name.",
+      filename, extension));
+  }
+
+  mcpl_outfile_t* file_id = nullptr;
+
+  if (mpi::master) {
+    file_id = g_mcpl_api->create_outfile(filename_.c_str());
+    if (!file_id) {
+      fatal_error(fmt::format(
+        "MCPL: Failed to create output file '{}'. Check permissions and path.",
+        filename_));
+    }
+    std::string src_line;
+    if (VERSION_DEV) {
+      src_line = fmt::format("OpenMC {}.{}.{}-dev{}", VERSION_MAJOR,
+        VERSION_MINOR, VERSION_RELEASE, VERSION_COMMIT_COUNT);
+    } else {
+      src_line = fmt::format(
+        "OpenMC {}.{}.{}", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE);
+    }
+
+    g_mcpl_api->hdr_set_srcname(file_id, src_line.c_str());
+  }
+  write_mcpl_collision_track_internal(
+    file_id, collision_track_bank, bank_index);
+
+  if (mpi::master) {
+    if (file_id) {
+      g_mcpl_api->close_outfile(file_id);
+    }
+  }
+}
+
 } // namespace openmc

src/message_passing.cpp

@@ -10,6 +10,7 @@ bool master {true};
 #ifdef OPENMC_MPI
 MPI_Comm intracomm {MPI_COMM_NULL};
 MPI_Datatype source_site {MPI_DATATYPE_NULL};
+MPI_Datatype collision_track_site {MPI_DATATYPE_NULL};
 #endif
 
 extern "C" bool openmc_master()

src/particle.cpp

@@ -8,6 +8,7 @@
 #include "openmc/bank.h"
 #include "openmc/capi.h"
 #include "openmc/cell.h"
+#include "openmc/collision_track.h"
 #include "openmc/constants.h"
 #include "openmc/dagmc.h"
 #include "openmc/error.h"
@@ -351,6 +352,11 @@ void Particle::event_collide()
     collision_mg(*this);
   }
 
+  // Collision track feature to recording particle interaction
+  if (settings::collision_track) {
+    collision_track_record(*this);
+  }
+
   // Score collision estimator tallies -- this is done after a collision
   // has occurred rather than before because we need information on the
   // outgoing energy for any tallies with an outgoing energy filter

src/physics.cpp

@@ -122,6 +122,7 @@ void sample_neutron_reaction(Particle& p)
           "with k-effective close to or greater than one.");
       }
     }
+    p.event_mt() = rx.mt_;
   }
 
   // Create secondary photons
@@ -663,7 +664,9 @@ void absorption(Particle& p, int i_nuclide)
 
       p.wgt() = 0.0;
       p.event() = TallyEvent::ABSORB;
-      p.event_mt() = N_DISAPPEAR;
+      if (!p.fission()) {
+        p.event_mt() = N_DISAPPEAR;
+      }
     }
   }
 }

src/settings.cpp

@@ -11,6 +11,7 @@
 #endif
 
 #include "openmc/capi.h"
+#include "openmc/collision_track.h"
 #include "openmc/constants.h"
 #include "openmc/container_util.h"
 #include "openmc/distribution.h"
@@ -26,6 +27,7 @@
 #include "openmc/plot.h"
 #include "openmc/random_lcg.h"
 #include "openmc/random_ray/random_ray.h"
+#include "openmc/reaction.h"
 #include "openmc/simulation.h"
 #include "openmc/source.h"
 #include "openmc/string_utils.h"
@@ -45,6 +47,7 @@ namespace settings {
 // Default values for boolean flags
 bool assume_separate {false};
 bool check_overlaps {false};
+bool collision_track {false};
 bool cmfd_run {false};
 bool confidence_intervals {false};
 bool create_delayed_neutrons {true};
@@ -128,6 +131,7 @@ std::unordered_set<int> statepoint_batch;
 double source_rejection_fraction {0.05};
 double free_gas_threshold {400.0};
 std::unordered_set<int> source_write_surf_id;
+CollisionTrackConfig collision_track_config {};
 int64_t ssw_max_particles;
 int64_t ssw_max_files;
 int64_t ssw_cell_id {C_NONE};
@@ -925,8 +929,72 @@ void read_settings_xml(pugi::xml_node root)
     }
   }
 
-  // If source is not separate and is to be written out in the statepoint file,
-  // make sure that the sourcepoint batch numbers are contained in the
+  // Check if the user has specified to write specific collisions
+  if (check_for_node(root, "collision_track")) {
+    settings::collision_track = true;
+    // Get collision track node
+    xml_node node_ct = root.child("collision_track");
+    collision_track_config = CollisionTrackConfig {};
+
+    // Determine cell ids at which crossing particles are to be banked
+    if (check_for_node(node_ct, "cell_ids")) {
+      auto temp = get_node_array<int>(node_ct, "cell_ids");
+      for (const auto& b : temp) {
+        collision_track_config.cell_ids.insert(b);
+      }
+    }
+    if (check_for_node(node_ct, "reactions")) {
+      auto temp = get_node_array<std::string>(node_ct, "reactions");
+      for (const auto& b : temp) {
+        int reaction_int = reaction_type(b);
+        if (reaction_int > 0) {
+          collision_track_config.mt_numbers.insert(reaction_int);
+        }
+      }
+    }
+    if (check_for_node(node_ct, "universe_ids")) {
+      auto temp = get_node_array<int>(node_ct, "universe_ids");
+      for (const auto& b : temp) {
+        collision_track_config.universe_ids.insert(b);
+      }
+    }
+    if (check_for_node(node_ct, "material_ids")) {
+      auto temp = get_node_array<int>(node_ct, "material_ids");
+      for (const auto& b : temp) {
+        collision_track_config.material_ids.insert(b);
+      }
+    }
+    if (check_for_node(node_ct, "nuclides")) {
+      auto temp = get_node_array<std::string>(node_ct, "nuclides");
+      for (const auto& b : temp) {
+        collision_track_config.nuclides.insert(b);
+      }
+    }
+    if (check_for_node(node_ct, "deposited_E_threshold")) {
+      collision_track_config.deposited_energy_threshold =
+        std::stod(get_node_value(node_ct, "deposited_E_threshold"));
+    }
+    // Get maximum number of particles to be banked per collision
+    if (check_for_node(node_ct, "max_collisions")) {
+      collision_track_config.max_collisions =
+        std::stoll(get_node_value(node_ct, "max_collisions"));
+    } else {
+      warning("A maximum number of collisions needs to be specified. "
+              "By default the code sets 'max_collisions' parameter equals to "
+              "1000.");
+    }
+    // Get maximum number of collision_track files to be created
+    if (check_for_node(node_ct, "max_collision_track_files")) {
+      collision_track_config.max_files =
+        std::stoll(get_node_value(node_ct, "max_collision_track_files"));
+    }
+    if (check_for_node(node_ct, "mcpl")) {
+      collision_track_config.mcpl_write = get_node_value_bool(node_ct, "mcpl");
+    }
+  }
+
+  // If source is not separate and is to be written out in the statepoint
+  // file, make sure that the sourcepoint batch numbers are contained in the
   // statepoint list
   if (!source_separate) {
     for (const auto& b : sourcepoint_batch) {
@@ -1171,8 +1239,8 @@ void read_settings_xml(pugi::xml_node root)
       variance_reduction::weight_windows_generators.emplace_back(
         std::make_unique<WeightWindowsGenerator>(node_wwg));
     }
-    // if any of the weight windows are intended to be generated otf, make sure
-    // they're applied
+    // if any of the weight windows are intended to be generated otf, make
+    // sure they're applied
     for (const auto& wwg : variance_reduction::weight_windows_generators) {
       if (wwg->on_the_fly_) {
         settings::weight_windows_on = true;

src/simulation.cpp

@@ -2,6 +2,7 @@
 
 #include "openmc/bank.h"
 #include "openmc/capi.h"
+#include "openmc/collision_track.h"
 #include "openmc/container_util.h"
 #include "openmc/eigenvalue.h"
 #include "openmc/error.h"
@@ -9,7 +10,6 @@
 #include "openmc/geometry_aux.h"
 #include "openmc/ifp.h"
 #include "openmc/material.h"
-#include "openmc/mcpl_interface.h"
 #include "openmc/message_passing.h"
 #include "openmc/nuclide.h"
 #include "openmc/output.h"
@@ -118,6 +118,7 @@ int openmc_simulation_init()
   // Reset global variables -- this is done before loading state point (as that
   // will potentially populate k_generation and entropy)
   simulation::current_batch = 0;
+  simulation::ct_current_file = 1;
   simulation::ssw_current_file = 1;
   simulation::k_generation.clear();
   simulation::entropy.clear();
@@ -297,6 +298,7 @@ namespace openmc {
 
 namespace simulation {
 
+int ct_current_file;
 int current_batch;
 int current_gen;
 bool initialized {false};
@@ -347,6 +349,11 @@ void allocate_banks()
     // Allocate surface source bank
     simulation::surf_source_bank.reserve(settings::ssw_max_particles);
   }
+
+  if (settings::collision_track) {
+    // Allocate collision track bank
+    collision_track_reserve_bank();
+  }
 }
 
 void initialize_batch()
@@ -490,6 +497,10 @@ void finalize_batch()
       ++simulation::ssw_current_file;
     }
   }
+  // Write collision track file if requested
+  if (settings::collision_track) {
+    collision_track_flush_bank();
+  }
 }
 
 void initialize_generation()

src/state_point.cpp

@@ -9,6 +9,7 @@
 #include <fmt/core.h>
 
 #include "openmc/bank.h"
+#include "openmc/bank_io.h"
 #include "openmc/capi.h"
 #include "openmc/constants.h"
 #include "openmc/eigenvalue.h"
@@ -642,91 +643,12 @@ void write_source_bank(hid_t group_id, span<SourceSite> source_bank,
 {
   hid_t banktype = h5banktype();
 
-  // Set total and individual process dataspace sizes for source bank
-  int64_t dims_size = bank_index.back();
-  int64_t count_size = bank_index[mpi::rank + 1] - bank_index[mpi::rank];
-
-#ifdef PHDF5
-  // Set size of total dataspace for all procs and rank
-  hsize_t dims[] {static_cast<hsize_t>(dims_size)};
-  hid_t dspace = H5Screate_simple(1, dims, nullptr);
-  hid_t dset = H5Dcreate(group_id, "source_bank", banktype, dspace, H5P_DEFAULT,
-    H5P_DEFAULT, H5P_DEFAULT);
-
-  // Create another data space but for each proc individually
-  hsize_t count[] {static_cast<hsize_t>(count_size)};
-  hid_t memspace = H5Screate_simple(1, count, nullptr);
-
-  // Select hyperslab for this dataspace
-  hsize_t start[] {static_cast<hsize_t>(bank_index[mpi::rank])};
-  H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
-
-  // Set up the property list for parallel writing
-  hid_t plist = H5Pcreate(H5P_DATASET_XFER);
-  H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
-
-  // Write data to file in parallel
-  H5Dwrite(dset, banktype, memspace, dspace, plist, source_bank.data());
-
-  // Free resources
-  H5Sclose(dspace);
-  H5Sclose(memspace);
-  H5Dclose(dset);
-  H5Pclose(plist);
-
+#ifdef OPENMC_MPI
+  write_bank_dataset("source_bank", group_id, source_bank, bank_index, banktype,
+    mpi::source_site);
 #else
-
-  if (mpi::master) {
-    // Create dataset big enough to hold all source sites
-    hsize_t dims[] {static_cast<hsize_t>(dims_size)};
-    hid_t dspace = H5Screate_simple(1, dims, nullptr);
-    hid_t dset = H5Dcreate(group_id, "source_bank", banktype, dspace,
-      H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-
-    // Save source bank sites since the array is overwritten below
-#ifdef OPENMC_MPI
-    vector<SourceSite> temp_source {source_bank.begin(), source_bank.end()};
-#endif
-
-    for (int i = 0; i < mpi::n_procs; ++i) {
-      // Create memory space
-      hsize_t count[] {static_cast<hsize_t>(bank_index[i + 1] - bank_index[i])};
-      hid_t memspace = H5Screate_simple(1, count, nullptr);
-
-#ifdef OPENMC_MPI
-      // Receive source sites from other processes
-      if (i > 0)
-        MPI_Recv(source_bank.data(), count[0], mpi::source_site, i, i,
-          mpi::intracomm, MPI_STATUS_IGNORE);
-#endif
-
-      // Select hyperslab for this dataspace
-      dspace = H5Dget_space(dset);
-      hsize_t start[] {static_cast<hsize_t>(bank_index[i])};
-      H5Sselect_hyperslab(
-        dspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
-
-      // Write data to hyperslab
-      H5Dwrite(
-        dset, banktype, memspace, dspace, H5P_DEFAULT, source_bank.data());
-
-      H5Sclose(memspace);
-      H5Sclose(dspace);
-    }
-
-    // Close all ids
-    H5Dclose(dset);
-
-#ifdef OPENMC_MPI
-    // Restore state of source bank
-    std::copy(temp_source.begin(), temp_source.end(), source_bank.begin());
-#endif
-  } else {
-#ifdef OPENMC_MPI
-    MPI_Send(source_bank.data(), count_size, mpi::source_site, 0, mpi::rank,
-      mpi::intracomm);
-#endif
-  }
+  write_bank_dataset(
+    "source_bank", group_id, source_bank, bank_index, banktype);
 #endif
 
   H5Tclose(banktype);

tests/regression_tests/collision_track/case_1_Reactions/inputs_true.dat

@@ -0,0 +1,58 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" depletable="true">
+      <density value="11.0" units="g/cm3"/>
+      <nuclide name="U234" ao="0.0004524"/>
+      <nuclide name="U235" ao="0.0506068"/>
+      <nuclide name="U238" ao="0.948709"/>
+      <nuclide name="U236" ao="0.0002318"/>
+      <nuclide name="O16" ao="2.0"/>
+    </material>
+    <material id="11">
+      <density value="1.0" units="g/cm3"/>
+      <nuclide name="H1" ao="2.0"/>
+      <nuclide name="O16" ao="1.0"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="5" fill="77" region="-5 4 -7 6 -9 8" universe="1"/>
+    <cell id="8" material="11" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="1"/>
+    <cell id="22" material="1" region="-1 2 -3" universe="77"/>
+    <cell id="33" material="11" region="1 | -2 | 3" universe="77"/>
+    <surface id="1" type="z-cylinder" coeffs="0.0 0.0 2.0"/>
+    <surface id="2" type="z-plane" coeffs="-2.0"/>
+    <surface id="3" type="z-plane" coeffs="2.0"/>
+    <surface id="4" type="x-plane" coeffs="-3.0"/>
+    <surface id="5" type="x-plane" coeffs="3.0"/>
+    <surface id="6" type="y-plane" coeffs="-3.0"/>
+    <surface id="7" type="y-plane" coeffs="3.0"/>
+    <surface id="8" type="z-plane" coeffs="-3.0"/>
+    <surface id="9" type="z-plane" coeffs="3.0"/>
+    <surface id="10" type="x-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="11" type="x-plane" boundary="vacuum" coeffs="4.0"/>
+    <surface id="12" type="y-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="13" type="y-plane" boundary="vacuum" coeffs="4.0"/>
+    <surface id="14" type="z-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="15" type="z-plane" boundary="vacuum" coeffs="4.0"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source type="independent" strength="1.0" particle="neutron">
+      <space type="box">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+      <constraints>
+        <fissionable>true</fissionable>
+      </constraints>
+    </source>
+    <collision_track>
+      <reactions>(n,fission) 101</reactions>
+      <max_collisions>300</max_collisions>
+    </collision_track>
+    <seed>1</seed>
+  </settings>
+  </model>

tests/regression_tests/collision_track/case_1_Reactions/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+5.642735E-02 1.494035E-02

tests/regression_tests/collision_track/case_2_Cell_ID/inputs_true.dat

@@ -0,0 +1,58 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" depletable="true">
+      <density value="11.0" units="g/cm3"/>
+      <nuclide name="U234" ao="0.0004524"/>
+      <nuclide name="U235" ao="0.0506068"/>
+      <nuclide name="U238" ao="0.948709"/>
+      <nuclide name="U236" ao="0.0002318"/>
+      <nuclide name="O16" ao="2.0"/>
+    </material>
+    <material id="11">
+      <density value="1.0" units="g/cm3"/>
+      <nuclide name="H1" ao="2.0"/>
+      <nuclide name="O16" ao="1.0"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="5" fill="77" region="-5 4 -7 6 -9 8" universe="1"/>
+    <cell id="8" material="11" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="1"/>
+    <cell id="22" material="1" region="-1 2 -3" universe="77"/>
+    <cell id="33" material="11" region="1 | -2 | 3" universe="77"/>
+    <surface id="1" type="z-cylinder" coeffs="0.0 0.0 2.0"/>
+    <surface id="2" type="z-plane" coeffs="-2.0"/>
+    <surface id="3" type="z-plane" coeffs="2.0"/>
+    <surface id="4" type="x-plane" coeffs="-3.0"/>
+    <surface id="5" type="x-plane" coeffs="3.0"/>
+    <surface id="6" type="y-plane" coeffs="-3.0"/>
+    <surface id="7" type="y-plane" coeffs="3.0"/>
+    <surface id="8" type="z-plane" coeffs="-3.0"/>
+    <surface id="9" type="z-plane" coeffs="3.0"/>
+    <surface id="10" type="x-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="11" type="x-plane" boundary="vacuum" coeffs="4.0"/>
+    <surface id="12" type="y-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="13" type="y-plane" boundary="vacuum" coeffs="4.0"/>
+    <surface id="14" type="z-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="15" type="z-plane" boundary="vacuum" coeffs="4.0"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source type="independent" strength="1.0" particle="neutron">
+      <space type="box">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+      <constraints>
+        <fissionable>true</fissionable>
+      </constraints>
+    </source>
+    <collision_track>
+      <cell_ids>22</cell_ids>
+      <max_collisions>300</max_collisions>
+    </collision_track>
+    <seed>1</seed>
+  </settings>
+  </model>

tests/regression_tests/collision_track/case_2_Cell_ID/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+5.642735E-02 1.494035E-02

tests/regression_tests/collision_track/case_3_Material_ID/inputs_true.dat

@@ -0,0 +1,58 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" depletable="true">
+      <density value="11.0" units="g/cm3"/>
+      <nuclide name="U234" ao="0.0004524"/>
+      <nuclide name="U235" ao="0.0506068"/>
+      <nuclide name="U238" ao="0.948709"/>
+      <nuclide name="U236" ao="0.0002318"/>
+      <nuclide name="O16" ao="2.0"/>
+    </material>
+    <material id="11">
+      <density value="1.0" units="g/cm3"/>
+      <nuclide name="H1" ao="2.0"/>
+      <nuclide name="O16" ao="1.0"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="5" fill="77" region="-5 4 -7 6 -9 8" universe="1"/>
+    <cell id="8" material="11" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="1"/>
+    <cell id="22" material="1" region="-1 2 -3" universe="77"/>
+    <cell id="33" material="11" region="1 | -2 | 3" universe="77"/>
+    <surface id="1" type="z-cylinder" coeffs="0.0 0.0 2.0"/>
+    <surface id="2" type="z-plane" coeffs="-2.0"/>
+    <surface id="3" type="z-plane" coeffs="2.0"/>
+    <surface id="4" type="x-plane" coeffs="-3.0"/>
+    <surface id="5" type="x-plane" coeffs="3.0"/>
+    <surface id="6" type="y-plane" coeffs="-3.0"/>
+    <surface id="7" type="y-plane" coeffs="3.0"/>
+    <surface id="8" type="z-plane" coeffs="-3.0"/>
+    <surface id="9" type="z-plane" coeffs="3.0"/>
+    <surface id="10" type="x-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="11" type="x-plane" boundary="vacuum" coeffs="4.0"/>
+    <surface id="12" type="y-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="13" type="y-plane" boundary="vacuum" coeffs="4.0"/>
+    <surface id="14" type="z-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="15" type="z-plane" boundary="vacuum" coeffs="4.0"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source type="independent" strength="1.0" particle="neutron">
+      <space type="box">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+      <constraints>
+        <fissionable>true</fissionable>
+      </constraints>
+    </source>
+    <collision_track>
+      <material_ids>1</material_ids>
+      <max_collisions>300</max_collisions>
+    </collision_track>
+    <seed>1</seed>
+  </settings>
+  </model>

tests/regression_tests/collision_track/case_3_Material_ID/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+5.642735E-02 1.494035E-02

tests/regression_tests/collision_track/case_4_Nuclide_ID/inputs_true.dat

@@ -0,0 +1,58 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" depletable="true">
+      <density value="11.0" units="g/cm3"/>
+      <nuclide name="U234" ao="0.0004524"/>
+      <nuclide name="U235" ao="0.0506068"/>
+      <nuclide name="U238" ao="0.948709"/>
+      <nuclide name="U236" ao="0.0002318"/>
+      <nuclide name="O16" ao="2.0"/>
+    </material>
+    <material id="11">
+      <density value="1.0" units="g/cm3"/>
+      <nuclide name="H1" ao="2.0"/>
+      <nuclide name="O16" ao="1.0"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="5" fill="77" region="-5 4 -7 6 -9 8" universe="1"/>
+    <cell id="8" material="11" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="1"/>
+    <cell id="22" material="1" region="-1 2 -3" universe="77"/>
+    <cell id="33" material="11" region="1 | -2 | 3" universe="77"/>
+    <surface id="1" type="z-cylinder" coeffs="0.0 0.0 2.0"/>
+    <surface id="2" type="z-plane" coeffs="-2.0"/>
+    <surface id="3" type="z-plane" coeffs="2.0"/>
+    <surface id="4" type="x-plane" coeffs="-3.0"/>
+    <surface id="5" type="x-plane" coeffs="3.0"/>
+    <surface id="6" type="y-plane" coeffs="-3.0"/>
+    <surface id="7" type="y-plane" coeffs="3.0"/>
+    <surface id="8" type="z-plane" coeffs="-3.0"/>
+    <surface id="9" type="z-plane" coeffs="3.0"/>
+    <surface id="10" type="x-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="11" type="x-plane" boundary="vacuum" coeffs="4.0"/>
+    <surface id="12" type="y-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="13" type="y-plane" boundary="vacuum" coeffs="4.0"/>
+    <surface id="14" type="z-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="15" type="z-plane" boundary="vacuum" coeffs="4.0"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source type="independent" strength="1.0" particle="neutron">
+      <space type="box">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+      <constraints>
+        <fissionable>true</fissionable>
+      </constraints>
+    </source>
+    <collision_track>
+      <nuclides>O16 U235</nuclides>
+      <max_collisions>300</max_collisions>
+    </collision_track>
+    <seed>1</seed>
+  </settings>
+  </model>

tests/regression_tests/collision_track/case_4_Nuclide_ID/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+5.642735E-02 1.494035E-02

tests/regression_tests/collision_track/case_5_Universe_ID/inputs_true.dat

@@ -0,0 +1,59 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" depletable="true">
+      <density value="11.0" units="g/cm3"/>
+      <nuclide name="U234" ao="0.0004524"/>
+      <nuclide name="U235" ao="0.0506068"/>
+      <nuclide name="U238" ao="0.948709"/>
+      <nuclide name="U236" ao="0.0002318"/>
+      <nuclide name="O16" ao="2.0"/>
+    </material>
+    <material id="11">
+      <density value="1.0" units="g/cm3"/>
+      <nuclide name="H1" ao="2.0"/>
+      <nuclide name="O16" ao="1.0"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="5" fill="77" region="-5 4 -7 6 -9 8" universe="1"/>
+    <cell id="8" material="11" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="1"/>
+    <cell id="22" material="1" region="-1 2 -3" universe="77"/>
+    <cell id="33" material="11" region="1 | -2 | 3" universe="77"/>
+    <surface id="1" type="z-cylinder" coeffs="0.0 0.0 2.0"/>
+    <surface id="2" type="z-plane" coeffs="-2.0"/>
+    <surface id="3" type="z-plane" coeffs="2.0"/>
+    <surface id="4" type="x-plane" coeffs="-3.0"/>
+    <surface id="5" type="x-plane" coeffs="3.0"/>
+    <surface id="6" type="y-plane" coeffs="-3.0"/>
+    <surface id="7" type="y-plane" coeffs="3.0"/>
+    <surface id="8" type="z-plane" coeffs="-3.0"/>
+    <surface id="9" type="z-plane" coeffs="3.0"/>
+    <surface id="10" type="x-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="11" type="x-plane" boundary="vacuum" coeffs="4.0"/>
+    <surface id="12" type="y-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="13" type="y-plane" boundary="vacuum" coeffs="4.0"/>
+    <surface id="14" type="z-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="15" type="z-plane" boundary="vacuum" coeffs="4.0"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source type="independent" strength="1.0" particle="neutron">
+      <space type="box">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+      <constraints>
+        <fissionable>true</fissionable>
+      </constraints>
+    </source>
+    <collision_track>
+      <cell_ids>22</cell_ids>
+      <universe_ids>77</universe_ids>
+      <max_collisions>300</max_collisions>
+    </collision_track>
+    <seed>1</seed>
+  </settings>
+  </model>

tests/regression_tests/collision_track/case_5_Universe_ID/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+5.642735E-02 1.494035E-02

tests/regression_tests/collision_track/case_6_deposited_energy_threshold/inputs_true.dat

@@ -0,0 +1,58 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" depletable="true">
+      <density value="11.0" units="g/cm3"/>
+      <nuclide name="U234" ao="0.0004524"/>
+      <nuclide name="U235" ao="0.0506068"/>
+      <nuclide name="U238" ao="0.948709"/>
+      <nuclide name="U236" ao="0.0002318"/>
+      <nuclide name="O16" ao="2.0"/>
+    </material>
+    <material id="11">
+      <density value="1.0" units="g/cm3"/>
+      <nuclide name="H1" ao="2.0"/>
+      <nuclide name="O16" ao="1.0"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="5" fill="77" region="-5 4 -7 6 -9 8" universe="1"/>
+    <cell id="8" material="11" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="1"/>
+    <cell id="22" material="1" region="-1 2 -3" universe="77"/>
+    <cell id="33" material="11" region="1 | -2 | 3" universe="77"/>
+    <surface id="1" type="z-cylinder" coeffs="0.0 0.0 2.0"/>
+    <surface id="2" type="z-plane" coeffs="-2.0"/>
+    <surface id="3" type="z-plane" coeffs="2.0"/>
+    <surface id="4" type="x-plane" coeffs="-3.0"/>
+    <surface id="5" type="x-plane" coeffs="3.0"/>
+    <surface id="6" type="y-plane" coeffs="-3.0"/>
+    <surface id="7" type="y-plane" coeffs="3.0"/>
+    <surface id="8" type="z-plane" coeffs="-3.0"/>
+    <surface id="9" type="z-plane" coeffs="3.0"/>
+    <surface id="10" type="x-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="11" type="x-plane" boundary="vacuum" coeffs="4.0"/>
+    <surface id="12" type="y-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="13" type="y-plane" boundary="vacuum" coeffs="4.0"/>
+    <surface id="14" type="z-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="15" type="z-plane" boundary="vacuum" coeffs="4.0"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source type="independent" strength="1.0" particle="neutron">
+      <space type="box">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+      <constraints>
+        <fissionable>true</fissionable>
+      </constraints>
+    </source>
+    <collision_track>
+      <deposited_E_threshold>550000.0</deposited_E_threshold>
+      <max_collisions>300</max_collisions>
+    </collision_track>
+    <seed>1</seed>
+  </settings>
+  </model>

tests/regression_tests/collision_track/case_6_deposited_energy_threshold/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+5.642735E-02 1.494035E-02

tests/regression_tests/collision_track/case_7_all_parameters_used_together/inputs_true.dat

@@ -0,0 +1,63 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" depletable="true">
+      <density value="11.0" units="g/cm3"/>
+      <nuclide name="U234" ao="0.0004524"/>
+      <nuclide name="U235" ao="0.0506068"/>
+      <nuclide name="U238" ao="0.948709"/>
+      <nuclide name="U236" ao="0.0002318"/>
+      <nuclide name="O16" ao="2.0"/>
+    </material>
+    <material id="11">
+      <density value="1.0" units="g/cm3"/>
+      <nuclide name="H1" ao="2.0"/>
+      <nuclide name="O16" ao="1.0"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="5" fill="77" region="-5 4 -7 6 -9 8" universe="1"/>
+    <cell id="8" material="11" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="1"/>
+    <cell id="22" material="1" region="-1 2 -3" universe="77"/>
+    <cell id="33" material="11" region="1 | -2 | 3" universe="77"/>
+    <surface id="1" type="z-cylinder" coeffs="0.0 0.0 2.0"/>
+    <surface id="2" type="z-plane" coeffs="-2.0"/>
+    <surface id="3" type="z-plane" coeffs="2.0"/>
+    <surface id="4" type="x-plane" coeffs="-3.0"/>
+    <surface id="5" type="x-plane" coeffs="3.0"/>
+    <surface id="6" type="y-plane" coeffs="-3.0"/>
+    <surface id="7" type="y-plane" coeffs="3.0"/>
+    <surface id="8" type="z-plane" coeffs="-3.0"/>
+    <surface id="9" type="z-plane" coeffs="3.0"/>
+    <surface id="10" type="x-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="11" type="x-plane" boundary="vacuum" coeffs="4.0"/>
+    <surface id="12" type="y-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="13" type="y-plane" boundary="vacuum" coeffs="4.0"/>
+    <surface id="14" type="z-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="15" type="z-plane" boundary="vacuum" coeffs="4.0"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source type="independent" strength="1.0" particle="neutron">
+      <space type="box">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+      <constraints>
+        <fissionable>true</fissionable>
+      </constraints>
+    </source>
+    <collision_track>
+      <cell_ids>22 33</cell_ids>
+      <reactions>elastic 18 (n,disappear)</reactions>
+      <universe_ids>77</universe_ids>
+      <material_ids>1 11</material_ids>
+      <nuclides>U238 U235 H1 U234</nuclides>
+      <deposited_E_threshold>100000.0</deposited_E_threshold>
+      <max_collisions>300</max_collisions>
+    </collision_track>
+    <seed>1</seed>
+  </settings>
+  </model>

tests/regression_tests/collision_track/case_7_all_parameters_used_together/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+5.642735E-02 1.494035E-02

tests/regression_tests/collision_track/case_8_2threads/inputs_true.dat

@@ -0,0 +1,57 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" depletable="true">
+      <density value="11.0" units="g/cm3"/>
+      <nuclide name="U234" ao="0.0004524"/>
+      <nuclide name="U235" ao="0.0506068"/>
+      <nuclide name="U238" ao="0.948709"/>
+      <nuclide name="U236" ao="0.0002318"/>
+      <nuclide name="O16" ao="2.0"/>
+    </material>
+    <material id="11">
+      <density value="1.0" units="g/cm3"/>
+      <nuclide name="H1" ao="2.0"/>
+      <nuclide name="O16" ao="1.0"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="5" fill="77" region="-5 4 -7 6 -9 8" universe="1"/>
+    <cell id="8" material="11" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="1"/>
+    <cell id="22" material="1" region="-1 2 -3" universe="77"/>
+    <cell id="33" material="11" region="1 | -2 | 3" universe="77"/>
+    <surface id="1" type="z-cylinder" coeffs="0.0 0.0 2.0"/>
+    <surface id="2" type="z-plane" coeffs="-2.0"/>
+    <surface id="3" type="z-plane" coeffs="2.0"/>
+    <surface id="4" type="x-plane" coeffs="-3.0"/>
+    <surface id="5" type="x-plane" coeffs="3.0"/>
+    <surface id="6" type="y-plane" coeffs="-3.0"/>
+    <surface id="7" type="y-plane" coeffs="3.0"/>
+    <surface id="8" type="z-plane" coeffs="-3.0"/>
+    <surface id="9" type="z-plane" coeffs="3.0"/>
+    <surface id="10" type="x-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="11" type="x-plane" boundary="vacuum" coeffs="4.0"/>
+    <surface id="12" type="y-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="13" type="y-plane" boundary="vacuum" coeffs="4.0"/>
+    <surface id="14" type="z-plane" boundary="vacuum" coeffs="-4.0"/>
+    <surface id="15" type="z-plane" boundary="vacuum" coeffs="4.0"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source type="independent" strength="1.0" particle="neutron">
+      <space type="box">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+      <constraints>
+        <fissionable>true</fissionable>
+      </constraints>
+    </source>
+    <collision_track>
+      <max_collisions>200</max_collisions>
+    </collision_track>
+    <seed>1</seed>
+  </settings>
+  </model>

tests/regression_tests/collision_track/case_8_2threads/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+5.642735E-02 1.494035E-02

tests/regression_tests/collision_track/test.py

@@ -0,0 +1,255 @@
+"""Test the 'collision_track' setting.
+
+Results
+-------
+
+All results are generated using only 1 MPI process.
+
+All results are generated using 1 thread except for "test_consistency_low_realization_number".
+This specific test verifies that when the number of realization (i.e., point being candidate
+to be stored) is lower than the capacity, results are reproducible even with multiple
+threads (i.e., there is no potential thread competition that would produce different
+results in that case).
+
+All results are generated using the history-based mode except for cases e01 to e03.
+
+All results are visually verified using the '_visualize.py' script in the regression test folder.
+
+OpenMC models
+-------------
+
+Four OpenMC models with CSG-only geometries are used to cover the transmission, vacuum,
+reflective and periodic Boundary Conditions (BC):
+
+- model_1: cylindrical core in 2 boxes (vacuum and transmission BC),
+
+# Test cases for simulation parameters using CSG-only geometries
+# ============================================================
+# Each test case is defined by a combination of folder name, model name, and specific parameters.
+# Below is a summary of the parameters used in the test cases:
+#
+# - max_collisions: Maximum number of particles to track in the simulation.
+# - reactions: List of MT numbers (reaction types- 2 for scattering, 18 for fission, 101 for absorbtion).
+# - cell_ids: IDs of specific cells in the model.
+# - mat_ids: Material IDs for filtering particles.
+# - nuclides: Nuclides for filtering particles.
+# - univ_ids: Universe IDs for filtering particles.
+# - E_threshold: Energy threshold for filtering particles (optional).
+#
+# The test cases are designed to validate the behavior of the simulation under various configurations.
+
+*: BC stands for Boundary Conditions, T for Transmission, R for Reflective, and V for Vacuum.
+
+An additional case, called 'case-a01', is used to check that the results are comparable when
+the number of threads is set to 2 if the number of realization is lower than the capacity.
+
+
+*: BC stands for Boundary Conditions, T for Transmission, and V for Vacuum.
+
+Notes:
+
+- The test cases list is non-exhaustive compared to the number of possible combinations.
+  Test cases have been selected based on use and internal code logic.
+
+
+
+TODO:
+
+- Test with a lattice.
+
+"""
+
+import os
+
+import openmc
+import openmc.lib
+import pytest
+
+from tests.testing_harness import CollisionTrackTestHarness
+from tests.regression_tests import config
+
+
+@pytest.fixture(scope="function")
+def two_threads(monkeypatch):
+    """Set the number of OMP threads to 2 for the test."""
+    monkeypatch.setenv("OMP_NUM_THREADS", "2")
+
+
+@pytest.fixture(scope="function")
+def single_process(monkeypatch):
+    """Set the number of MPI process to 1 for the test."""
+    monkeypatch.setitem(config, "mpi_np", "1")
+
+
+@pytest.fixture(scope="module")
+def model_1():
+    """Cylindrical core contained in a first box which is contained in a larger box.
+    A lower universe is used to describe the interior of the first box which
+    contains the core and its surrounding space.
+
+    """
+    openmc.reset_auto_ids()
+    model = openmc.Model()
+
+    # =============================================================================
+    # Materials
+    # =============================================================================
+
+    fuel = openmc.Material(material_id=1)
+    fuel.add_nuclide("U234", 0.0004524)
+    fuel.add_nuclide("U235", 0.0506068)
+    fuel.add_nuclide("U238", 0.9487090)
+    fuel.add_nuclide("U236", 0.0002318)
+    fuel.add_nuclide("O16", 2.0)
+    fuel.set_density("g/cm3", 11.0)
+
+    water = openmc.Material(material_id=11)
+    water.add_nuclide("H1", 2.0)
+    water.add_nuclide("O16", 1.0)
+    water.set_density("g/cm3", 1.0)
+
+    # =============================================================================
+    # Geometry
+    # =============================================================================
+
+    # -----------------------------------------------------------------------------
+    # Cylindrical core
+    # -----------------------------------------------------------------------------
+
+    # Parameters
+    core_radius = 2.0
+    core_height = 4.0
+
+    # Surfaces
+    core_cylinder = openmc.ZCylinder(r=core_radius)
+    core_lower_plane = openmc.ZPlane(-core_height / 2.0)
+    core_upper_plane = openmc.ZPlane(core_height / 2.0)
+
+    # Region
+    core_region = -core_cylinder & +core_lower_plane & -core_upper_plane
+
+    # Cells
+    core = openmc.Cell(fill=fuel, region=core_region, cell_id=22)
+    outside_core_region = +core_cylinder | -core_lower_plane | +core_upper_plane
+    outside_core = openmc.Cell(
+        fill=water, region=outside_core_region, cell_id=33)
+
+    # Universe
+    inside_box1_universe = openmc.Universe(
+        cells=[core, outside_core], universe_id=77)
+
+    # -----------------------------------------------------------------------------
+    # Box 1
+    # -----------------------------------------------------------------------------
+
+    # Parameters
+    box1_size = 6.0
+
+    # Surfaces
+    box1_rpp = openmc.model.RectangularParallelepiped(
+        -box1_size / 2.0, box1_size / 2.0,
+        -box1_size / 2.0, box1_size / 2.0,
+        -box1_size / 2.0, box1_size / 2.0,
+    )
+
+    # Cell
+    box1 = openmc.Cell(fill=inside_box1_universe, region=-box1_rpp, cell_id=5)
+
+    # -----------------------------------------------------------------------------
+    # Box 2
+    # -----------------------------------------------------------------------------
+
+    # Parameters
+    box2_size = 8
+
+    # Surfaces
+    box2_rpp = openmc.model.RectangularParallelepiped(
+        -box2_size / 2.0, box2_size / 2.0,
+        -box2_size / 2.0, box2_size / 2.0,
+        -box2_size / 2.0, box2_size / 2.0,
+        boundary_type="vacuum"
+    )
+
+    # Cell
+    box2 = openmc.Cell(fill=water, region=-box2_rpp & +box1_rpp, cell_id=8)
+
+    # Register geometry
+    model.geometry = openmc.Geometry([box1, box2])
+
+    # =============================================================================
+    # Settings
+    # =============================================================================
+
+    model.settings = openmc.Settings()
+    model.settings.particles = 100
+    model.settings.batches = 5
+    model.settings.inactive = 1
+    model.settings.seed = 1
+
+    bounds = [
+        -core_radius,
+        -core_radius,
+        -core_height / 2.0,
+        core_radius,
+        core_radius,
+        core_height / 2.0,
+    ]
+    distribution = openmc.stats.Box(bounds[:3], bounds[3:])
+    model.settings.source = openmc.IndependentSource(
+        space=distribution, constraints={'fissionable': True})
+
+    return model
+
+
+@pytest.mark.parametrize(
+    "folder, model_name, parameter",
+    [("case_1_Reactions", "model_1", {"max_collisions": 300, "reactions": ["(n,fission)", 101]}),
+        ("case_2_Cell_ID", "model_1", {
+         "max_collisions": 300, "cell_ids": [22]}),
+        ("case_3_Material_ID", "model_1", {
+         "max_collisions": 300, "material_ids": [1]}),
+        ("case_4_Nuclide_ID", "model_1", {
+         "max_collisions": 300, "nuclides": ["O16", "U235"]}),
+        ("case_5_Universe_ID", "model_1", {
+         "max_collisions": 300, "cell_ids": [22], "universe_ids": [77]}),
+        ("case_6_deposited_energy_threshold", "model_1", {
+         "max_collisions": 300, "deposited_E_threshold": 5.5e5}),
+        ("case_7_all_parameters_used_together", "model_1", {
+            "max_collisions": 300,
+            "reactions": ["elastic", 18, "(n,disappear)"],
+            "material_ids": [1, 11],
+            "universe_ids": [77],
+            "nuclides": ["U238", "U235", "H1", "U234"],
+            "cell_ids": [22, 33],
+            "deposited_E_threshold": 1e5})
+     ],
+)
+def test_collision_track_several_cases(
+    folder, model_name, parameter, request
+):
+    # Since for these tests the actual number of collisions recorded is < max_collisions,
+    # we can run them with 1 or 2 threads, and in history or event mode.
+    model = request.getfixturevalue(model_name)
+    model.settings.collision_track = parameter
+    harness = CollisionTrackTestHarness(
+        "statepoint.5.h5", model=model, workdir=folder
+    )
+    harness.main()
+
+
+@pytest.mark.skipif(config["event"], reason="Results from history-based mode.")
+def test_collision_track_2threads(model_1, two_threads, single_process):
+    # This test checks that the `max_collisions` setting is honored:
+    # no collisions beyond the specified limit should be recorded.
+    #
+    # For the result to be reproducible, the number of threads and
+    # the transport mode (history vs. event) must remain fixed.
+    assert os.environ["OMP_NUM_THREADS"] == "2"
+    assert config["mpi_np"] == "1"
+    model_1.settings.collision_track = {
+        "max_collisions": 200
+    }
+    harness = CollisionTrackTestHarness(
+        "statepoint.5.h5", model=model_1, workdir="case_8_2threads"
+    )
+    harness.main()

tests/testing_harness.py

@@ -68,7 +68,7 @@ class TestHarness:
         if config['mpi']:
             mpi_args = [config['mpiexec'], '-n', config['mpi_np']]
             openmc.run(openmc_exec=config['exe'], mpi_args=mpi_args,
-              event_based=config['event'])
+                       event_based=config['event'])
         else:
             openmc.run(openmc_exec=config['exe'], event_based=config['event'])
 
@@ -305,9 +305,12 @@ class PyAPITestHarness(TestHarness):
         else:
             self.execute_test()
 
-    def execute_test(self):
+    def execute_test(self, change_dir=False):
         """Build input XMLs, run OpenMC, and verify correct results."""
+        base_dir = os.getcwd() if change_dir else None
         try:
+            if change_dir:
+                os.chdir(self.workdir)
             self._build_inputs()
             inputs = self._get_inputs()
             self._write_inputs(inputs)
@@ -319,10 +322,15 @@ class PyAPITestHarness(TestHarness):
             self._compare_results()
         finally:
             self._cleanup()
+            if base_dir:
+                os.chdir(base_dir)
 
-    def update_results(self):
+    def update_results(self, change_dir=False):
         """Update results_true.dat and inputs_true.dat"""
+        base_dir = os.getcwd() if change_dir else None
         try:
+            if change_dir:
+                os.chdir(self.workdir)
             self._build_inputs()
             inputs = self._get_inputs()
             self._write_inputs(inputs)
@@ -334,6 +342,8 @@ class PyAPITestHarness(TestHarness):
             self._overwrite_results()
         finally:
             self._cleanup()
+            if base_dir:
+                os.chdir(base_dir)
 
     def _build_inputs(self):
         """Write input XML files."""
@@ -371,7 +381,8 @@ class PyAPITestHarness(TestHarness):
         """Delete XMLs, statepoints, tally, and test files."""
         super()._cleanup()
         output = ['materials.xml', 'geometry.xml', 'settings.xml',
-                  'tallies.xml', 'plots.xml', 'inputs_test.dat', 'model.xml']
+                  'tallies.xml', 'plots.xml', 'inputs_test.dat', 'model.xml',
+                  'collision_track.h5', 'collision_track.mcpl']
         for f in output:
             if os.path.exists(f):
                 os.remove(f)
@@ -389,6 +400,7 @@ class TolerantPyAPITestHarness(PyAPITestHarness):
     due to single precision usage (e.g., as in the random ray solver).
 
     """
+
     def _are_files_equal(self, actual_path, expected_path, tolerance):
         def isfloat(value):
             try:
@@ -428,7 +440,8 @@ class TolerantPyAPITestHarness(PyAPITestHarness):
 
     def _compare_results(self):
         """Make sure the current results agree with the reference."""
-        compare = self._are_files_equal('results_test.dat', 'results_true.dat', 1e-6)
+        compare = self._are_files_equal(
+            'results_test.dat', 'results_true.dat', 1e-6)
         if not compare:
             expected = open('results_true.dat').readlines()
             actual = open('results_test.dat').readlines()
@@ -476,6 +489,7 @@ class WeightWindowPyAPITestHarness(PyAPITestHarness):
 
 class PlotTestHarness(TestHarness):
     """Specialized TestHarness for running OpenMC plotting tests."""
+
     def __init__(self, plot_names, voxel_convert_checks=[]):
         super().__init__(None)
         self._plot_names = plot_names
@@ -523,3 +537,105 @@ class PlotTestHarness(TestHarness):
         outstr = sha512.hexdigest()
 
         return outstr
+
+
+class CollisionTrackTestHarness(PyAPITestHarness):
+    def __init__(self, statepoint_name, model=None, inputs_true=None, workdir=None):
+        super().__init__(statepoint_name, model, inputs_true)
+        self.workdir = workdir
+
+    def _test_output_created(self):
+        """Make sure collision_track.h5 has also been created."""
+        super()._test_output_created()
+        if self._model.settings.collision_track:
+            assert os.path.exists(
+                "collision_track.h5"
+            ), "collision_track file has not been created."
+
+    def _compare_output(self):
+        """Compare collision_track.h5 files."""
+        if self._model.settings.collision_track:
+            collision_track_true = self._return_collision_track_data(
+                "collision_track_true.h5")
+            collision_track_test = self._return_collision_track_data(
+                "collision_track.h5")
+            np.testing.assert_allclose(
+                collision_track_true, collision_track_test, rtol=1e-07)
+
+    def main(self):
+        """Accept commandline arguments and either run or update tests."""
+        if config["build_inputs"]:
+            self.build_inputs()
+        elif config["update"]:
+            self.update_results(change_dir=True)
+        else:
+            self.execute_test(change_dir=True)
+
+    def build_inputs(self):
+        """Build inputs."""
+        base_dir = os.getcwd()
+        try:
+            os.chdir(self.workdir)
+            self._build_inputs()
+        finally:
+            os.chdir(base_dir)
+
+    def _overwrite_results(self):
+        """Also add the 'collision_track.h5' file during overwriting."""
+        super()._overwrite_results()
+        if os.path.exists("collision_track.h5"):
+            shutil.copyfile("collision_track.h5", "collision_track_true.h5")
+
+    @staticmethod
+    def _return_collision_track_data(filepath):
+        """
+        Read a collision_track file and return a sorted array composed
+        of flatten arrays of collision information.
+
+        Parameters
+        ----------
+        filepath : str
+            Path to the collision_track file
+
+        Returns
+        -------
+        data : np.array
+            Sorted array composed of flatten arrays of collision_track data for
+            each collision information
+        """
+        data = []
+        keys = []
+
+        # Read source file
+        source = openmc.read_collision_track_file(filepath)
+        for src in source:
+            r = src['r']
+            u = src['u']
+            e = src['E']
+            de = src['dE']
+            time = src['time']
+            wgt = src['wgt']
+            delayed_group = src['delayed_group']
+            cell_id = src['cell_id']
+            nuclide_id = src['nuclide_id']
+            material_id = src['material_id']
+            universe_id = src['universe_id']
+            n_collision = src['n_collision']
+            event_mt = src['event_mt']
+            key = (
+                f"{r[0]:.10e} {r[1]:.10e} {r[2]:.10e} {u[0]:.10e} {u[1]:.10e} {u[2]:.10e}"
+                f"{e:.10e} {de:.10e}  {time:.10e} {wgt:.10e} {event_mt} {delayed_group} {cell_id}"
+                f"{nuclide_id} {material_id} {universe_id} {n_collision} "
+            )
+            keys.append(key)
+            values = [*r, *u, e, de, time, wgt, event_mt,
+                      delayed_group, cell_id, nuclide_id, material_id,
+                      universe_id, n_collision]
+            assert len(values) == 17
+            data.append(values)
+
+        data = np.array(data)
+        keys = np.array(keys)
+        sorted_idx = np.argsort(keys, kind='stable')
+
+        return data[sorted_idx]

tests/unit_tests/test_collision_track.py

@@ -0,0 +1,127 @@
+"""Test the 'collision_track' setting used to store particle information
+during specified collision conditions in a file for a given simulation."""
+
+import openmc
+import pytest
+import h5py
+import numpy as np
+
+from tests.testing_harness import CollisionTrackTestHarness as ctt
+
+
+@pytest.fixture(scope="module")
+def geometry():
+    """Simple hydrogen sphere geometry"""
+    openmc.reset_auto_ids()
+    material = openmc.Material(name="H1")
+    material.add_element("H", 1.0)
+    sphere = openmc.Sphere(r=1.0, boundary_type="vacuum")
+    cell = openmc.Cell(region=-sphere, fill=material)
+    return openmc.Geometry([cell])
+
+
+@pytest.mark.parametrize(
+    "parameter",
+    [
+        {"max_collisions": 200},
+        {"max_collisions": 200, "reactions": ["(n,disappear)"]},
+        {"max_collisions": 200, "cell_ids": [1]},
+        {"max_collisions": 200, "material_ids": [1]},
+        {"max_collisions": 200, "universe_ids": [1]},
+        {"max_collisions": 200, "nuclides": ["H1"]},
+        {"max_collisions": 200, "deposited_E_threshold": 200000.0},
+        {"max_collisions": 200, "mcpl": True}
+
+    ],
+)
+def test_xml_serialization(parameter, run_in_tmpdir):
+    """Check that the different use cases can be written and read in XML."""
+    settings = openmc.Settings()
+    settings.collision_track = parameter
+    settings.export_to_xml()
+
+    read_settings = openmc.Settings.from_xml()
+    assert read_settings.collision_track == parameter
+
+
+@pytest.fixture(scope="module")
+def model():
+    """Simple hydrogen sphere divided in two hemispheres
+    by a z-plane to form 2 cells."""
+    openmc.reset_auto_ids()
+    model = openmc.Model()
+
+    # Material
+    material = openmc.Material(name="H1")
+    material.add_element("H", 1.0)
+
+    # Geometry
+    radius = 1.0
+    sphere = openmc.Sphere(r=radius, boundary_type="reflective")
+    plane = openmc.ZPlane(0.0)
+    cell_1 = openmc.Cell(region=-sphere & -plane, fill=material, cell_id=1)
+    cell_2 = openmc.Cell(region=-sphere & +plane, fill=material, cell_id=2)
+    root = openmc.Universe(cells=[cell_1, cell_2])
+    model.geometry = openmc.Geometry(root)
+
+    # Settings
+    model.settings = openmc.Settings()
+    model.settings.run_mode = "fixed source"
+    model.settings.particles = 1
+    model.settings.batches = 1
+    model.settings.seed = 2
+
+    bounds = [-radius, -radius, -radius, radius, radius, radius]
+    distribution = openmc.stats.Box(bounds[:3], bounds[3:])
+    model.settings.source = openmc.IndependentSource(space=distribution)
+
+    return model
+
+
+def test_particle_location(run_in_tmpdir, model):
+    """Test the location of particles with respected to the "cell_ids"
+    and the location x, y, z of the particle itself. the upper sphere will
+    have positive z component and the bottom sphere a negative z compnent.
+
+    """
+    model.settings.collision_track = {
+        "max_collisions": 200,
+        "reactions": ["elastic"],
+        "cell_ids": [1, 2]
+    }
+    model.run()
+
+    with h5py.File("collision_track.h5", "r") as f:
+        source = f["collision_track_bank"]
+
+        assert len(source) == 60
+
+        # We want to verify that the collisions happenening are in the right cells
+        # and the position of the particle is either positive or negative relative
+        # to the z plane. In this case, we track the position of the particle
+        # relative to the cell_id already set.
+        for point in source:
+            if point['cell_id'] == 1:
+                assert point['r'][2] < 0.0  # z component negative
+            elif point['cell_id'] == 2:
+                assert point["r"][2] > 0.0  # z component positive
+            else:
+                assert False
+
+
+def test_format_similarity(run_in_tmpdir, model):
+    model.settings.collision_track = {"max_collisions": 200, "reactions": ['elastic'],
+                                      "cell_ids": [1, 2], "mcpl": False}
+    model.run()
+    data_h5 = ctt._return_collision_track_data('collision_track.h5')
+
+    model.settings.collision_track["mcpl"] = True
+    model.run()
+    data_mcpl = ctt._return_collision_track_data('collision_track.mcpl')
+
+    assert len(data_h5) == 60
+    assert len(data_mcpl) == 60
+
+    np.testing.assert_allclose(data_h5, data_mcpl, rtol=1e-05)
+    # tolerance not that low due to the strings that is saved in MCPL,
+    # not enough precision!