diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst
index 57242047a58..8a971c213c0 100644
--- a/Docs/source/usage/parameters.rst
+++ b/Docs/source/usage/parameters.rst
@@ -3064,6 +3064,23 @@ Details about the collision models can be found in the :ref:`theory section <mul
     With ``forward``, the scattering angle is set to zero, i.e. the products are emitted in the same direction as the reactant (in the center of mass frame).
     With ``backward``, the scattering angle is set to :math:`\pi`, i.e. the products are emitted in the opposite direction of the reactant (in the center of mass frame).
 
+.. pp:param:: <collision_name>.create_products
+    :type: ``bool``
+    :default: ``1``
+    :optional:
+
+    Only for ``nuclearfusion``. When true, the product particles are created, otherwise not.
+
+.. pp:param:: <collision_name>.save_particle_production
+    :type: ``bool``
+    :default: ``0``
+    :optional:
+
+    Only for ``nuclearfusion``.
+    When true, the integrated product particle density is saved in a MultiFab with the name ``<collision_name>_particle_production``.
+    The data can be written out by adding that name to the ``<diag_name>.fields_to_plot`` input parameter.
+    The option can be used in conjunction with ``<collision_name>.create_products`` to save only the product density and not create particles.
+
 .. pp:param:: <collision_name>.background_density
     :type: ``float``
 
diff --git a/Examples/Tests/nuclear_fusion/analysis_two_product_fusion.py b/Examples/Tests/nuclear_fusion/analysis_two_product_fusion.py
index 26f32c63a2e..ffc9b29e346 100755
--- a/Examples/Tests/nuclear_fusion/analysis_two_product_fusion.py
+++ b/Examples/Tests/nuclear_fusion/analysis_two_product_fusion.py
@@ -65,6 +65,7 @@
     reaction_type = "DT"
     reactant_species = ["deuterium", "tritium"]
     product_species = ["helium4", "neutron"]
+    collision_names = ["DTF1", "DTF2"]
     ntests = 2
     E_fusion = 17.58929696 * MeV_to_Joule
 else:
@@ -72,6 +73,7 @@
     reaction_type = "DD"
     reactant_species = ["deuterium", "hydrogen2"]
     product_species = ["helium3", "neutron"]
+    collision_names = ["DDNHeF1"]
     ntests = 1
     E_fusion = 3.26891111e6 * MeV_to_Joule
 
@@ -399,6 +401,19 @@ def check_macroparticle_number(
         atol=5.0 * std_macroparticle_number,
     )
 
+    if "particle_production" in data:
+        w_sum = data[product_species[0] + "_w_end"].sum()
+        n_sum = data["particle_production"].sum()
+        tolerance = 0.02
+        print(
+            f"Check particle production diagnostic for collision {data['collision_name']}:"
+        )
+        print(f"from particles  = {w_sum}")
+        print(f"from diagnostic = {n_sum}")
+        print(f"error = {np.abs(w_sum - n_sum) / w_sum}")
+        print(f"tolerance = {tolerance}")
+        assert is_close(w_sum, n_sum, rtol=tolerance)
+
     ## used in subsequent function
     return expected_fusion_number
 
@@ -545,6 +560,13 @@ def main():
         # General checks that are performed for all tests
         generic_check(data)
 
+        product_production_name = f"{collision_names[i - 1]}_particle_production"
+        if ("boxlib", product_production_name) in ds_end.field_list:
+            data["collision_name"] = collision_names[i - 1]
+            data["particle_production"] = field_data_end[
+                "boxlib", product_production_name
+            ].to_ndarray()
+
         # Checks that are specific to test number i
         eval("specific_check" + str(i) + "(data, dt)")
 
diff --git a/Examples/Tests/nuclear_fusion/inputs_test_3d_deuterium_tritium_fusion b/Examples/Tests/nuclear_fusion/inputs_test_3d_deuterium_tritium_fusion
index 759f136c845..1a1d5488a40 100644
--- a/Examples/Tests/nuclear_fusion/inputs_test_3d_deuterium_tritium_fusion
+++ b/Examples/Tests/nuclear_fusion/inputs_test_3d_deuterium_tritium_fusion
@@ -2,7 +2,9 @@
 ####### GENERAL PARAMETERS ######
 #################################
 ## With these parameters, each cell has a size of exactly 1 by 1 by 1
-max_step = 1
+# Run two steps. The first step to check the fusion collisions, the
+# second step to provide a comparison for the rerun.
+max_step = 2
 amr.n_cell = 8 8 16
 amr.max_grid_size = 8
 amr.blocking_factor = 8
@@ -121,15 +123,24 @@ DTF1.species = deuterium_1 tritium_1
 DTF1.product_species = helium4_1 neutron_1
 DTF1.type = nuclearfusion
 DTF1.event_multiplier = 1.e50
+DTF1.create_products = 1
+DTF1.save_particle_production = 1
 
 DTF2.species = deuterium_2 tritium_2
 DTF2.product_species = helium4_2 neutron_2
 DTF2.type = nuclearfusion
 DTF2.event_multiplier = 1.e15
 DTF2.probability_target_value = 0.02
+DTF2.create_products = 1
+DTF2.save_particle_production = 1
 
 # Diagnostics
-diagnostics.diags_names = diag1
+diagnostics.diags_names = diag1 checkpoint
 diag1.intervals = 1
 diag1.diag_type = Full
-diag1.fields_to_plot = rho
+diag1.fields_to_plot = rho DTF1_particle_production DTF2_particle_production
+
+checkpoint.format = checkpoint
+checkpoint.diag_type = Full
+checkpoint.intervals = 1:1:1
+checkpoint.dump_last_timestep = 0
diff --git a/Regression/Checksum/benchmarks_json/test_3d_deuterium_tritium_fusion.json b/Regression/Checksum/benchmarks_json/test_3d_deuterium_tritium_fusion.json
index a0000ee6f59..17e0893cc67 100644
--- a/Regression/Checksum/benchmarks_json/test_3d_deuterium_tritium_fusion.json
+++ b/Regression/Checksum/benchmarks_json/test_3d_deuterium_tritium_fusion.json
@@ -1,24 +1,12 @@
 {
-  "lev=0": {
-    "rho": 0.0
-  },
-  "neutron_2": {
-    "particle_momentum_x": 1.5369063360838545e-15,
-    "particle_momentum_y": 1.5327717119671177e-15,
-    "particle_momentum_z": 1.5632203763888702e-15,
-    "particle_position_x": 136756.17264787608,
-    "particle_position_y": 136453.48037878488,
-    "particle_position_z": 290503.22456411575,
-    "particle_weight": 6.347081228434342e+18
-  },
-  "neutron_1": {
-    "particle_momentum_x": 1.7270063957926637e-15,
-    "particle_momentum_y": 1.7295255445271788e-15,
-    "particle_momentum_z": 1.7442619148907942e-15,
-    "particle_position_x": 151546.7150677448,
-    "particle_position_y": 151695.5086129642,
-    "particle_position_z": 323004.4593236664,
-    "particle_weight": 4.337788155202713e-28
+  "deuterium_1": {
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 2.8872569136407634e-13,
+    "particle_position_x": 40958427.50992301,
+    "particle_position_y": 40959476.34450768,
+    "particle_position_z": 81921930.27522022,
+    "particle_weight": 1024.0000000000002
   },
   "deuterium_2": {
     "particle_momentum_x": 0.0,
@@ -29,23 +17,14 @@
     "particle_position_z": 8192362.405430986,
     "particle_weight": 1.0240001137714307e+30
   },
-  "tritium_1": {
-    "particle_momentum_x": 0.0,
-    "particle_momentum_y": 0.0,
-    "particle_momentum_z": 2.887256913640763e-13,
-    "particle_position_x": 40959200.11081588,
-    "particle_position_y": 40960650.407891415,
-    "particle_position_z": 81920772.7986121,
-    "particle_weight": 1024.0000000000002
-  },
-  "deuterium_1": {
-    "particle_momentum_x": 0.0,
-    "particle_momentum_y": 0.0,
-    "particle_momentum_z": 2.8872569136407634e-13,
-    "particle_position_x": 40958427.50992301,
-    "particle_position_y": 40959476.34450768,
-    "particle_position_z": 81921930.27522022,
-    "particle_weight": 1024.0000000000002
+  "helium4_1": {
+    "particle_momentum_x": 1.7270063957926637e-15,
+    "particle_momentum_y": 1.7295255445271788e-15,
+    "particle_momentum_z": 1.7442619148907942e-15,
+    "particle_position_x": 151546.7150677448,
+    "particle_position_y": 151695.5086129642,
+    "particle_position_z": 323004.4593236664,
+    "particle_weight": 4.337788155202713e-28
   },
   "helium4_2": {
     "particle_momentum_x": 1.5369063360838545e-15,
@@ -56,7 +35,12 @@
     "particle_position_z": 290503.22456411575,
     "particle_weight": 6.347081228434342e+18
   },
-  "helium4_1": {
+  "lev=0": {
+    "DTF1_particle_production": 4.415241048343040e-28,
+    "DTF2_particle_production": 6.295151973703336e+18,
+    "rho": 0.0
+  },
+  "neutron_1": {
     "particle_momentum_x": 1.7270063957926637e-15,
     "particle_momentum_y": 1.7295255445271788e-15,
     "particle_momentum_z": 1.7442619148907942e-15,
@@ -65,6 +49,24 @@
     "particle_position_z": 323004.4593236664,
     "particle_weight": 4.337788155202713e-28
   },
+  "neutron_2": {
+    "particle_momentum_x": 1.5369063360838545e-15,
+    "particle_momentum_y": 1.5327717119671177e-15,
+    "particle_momentum_z": 1.5632203763888702e-15,
+    "particle_position_x": 136756.17264787608,
+    "particle_position_y": 136453.48037878488,
+    "particle_position_z": 290503.22456411575,
+    "particle_weight": 6.347081228434342e+18
+  },
+  "tritium_1": {
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 2.887256913640763e-13,
+    "particle_position_x": 40959200.11081588,
+    "particle_position_y": 40960650.407891415,
+    "particle_position_z": 81920772.7986121,
+    "particle_weight": 1024.0000000000002
+  },
   "tritium_2": {
     "particle_momentum_x": 0.0,
     "particle_momentum_y": 0.0,
@@ -74,4 +76,4 @@
     "particle_position_z": 819126.8984535292,
     "particle_weight": 1.0239999999365294e+29
   }
-}
\ No newline at end of file
+}
diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H
index 604fa45369f..044935fb0ab 100644
--- a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H
+++ b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H
@@ -177,6 +177,11 @@ public:
     BinaryCollision ( BinaryCollision&& )                  = delete;
     BinaryCollision& operator= ( BinaryCollision&& )       = delete;
 
+    void AllocData () override
+    {
+        m_binary_collision_functor.AllocData();
+    }
+
     /** Perform the collisions
      *
      * @param cur_time Current time
@@ -296,7 +301,7 @@ public:
 
         ABLASTR_PROFILE("BinaryCollision::doCollisionsWithinTile");
 
-        const auto& binary_collision_functor = m_binary_collision_functor.executor();
+        const auto& binary_collision_functor = m_binary_collision_functor.executor(lev, mfi);
         const bool have_product_species = m_have_product_species;
 
         // Store product species data in vectors
@@ -326,6 +331,8 @@ public:
             global_debye_length_data = global_debye_length_fab.dataPtr();
         }
 
+        const amrex::Box tilebox = mfi.tilebox();
+
         amrex::Geometry const& geom_lev = WarpX::GetInstance().Geom(lev);
         // dV is level-specific: cell volume at this refinement level (smaller on fine levels).
         amrex::ParticleReal const dV = AMREX_D_TERM(geom_lev.CellSize(0), *geom_lev.CellSize(1), *geom_lev.CellSize(2));
@@ -652,6 +659,7 @@ public:
 
                     // Use a bisection algorithm to find the index of the cell in which this pair is located
                     const int i_cell = amrex::bisect( p_coll_offsets, 0, n_cells, ui_coll );
+                    const amrex::IntVect global_index = tilebox.atOffset(i_cell);
 
                     // The particles from species1 that are in the cell `i_cell` are
                     // given by the `indices_1[cell_start_1:cell_stop_1]`
@@ -691,7 +699,7 @@ public:
                         soa_1, soa_1, get_position_1, get_position_1,
                         n1, n1, T1, T1, global_lamdb,
                         q1, q1, m1, m1, dt, dV*volume_factor(i_cell), coll_idx,
-                        cell_start_pair, p_mask, p_pair_indices_1, p_pair_indices_2,
+                        cell_start_pair, global_index, p_mask, p_pair_indices_1, p_pair_indices_2,
                         p_pair_reaction_weight, p_product_data, engine);
                 }
             );
@@ -1238,6 +1246,7 @@ public:
 
                     // Use a bisection algorithm to find the index of the cell in which this pair is located
                     const int i_cell = amrex::bisect( p_coll_offsets, 0, n_cells, ui_coll );
+                    const amrex::IntVect global_index = tilebox.atOffset(i_cell);
 
                     // The particles from species1 that are in the cell `i_cell` are
                     // given by the `indices_1[cell_start_1:cell_stop_1]`
@@ -1284,7 +1293,7 @@ public:
                         soa_1, soa_2, get_position_1, get_position_2,
                         n1, n2, T1, T2, global_lamdb,
                         q1, q2, m1, m2, dt, dV*volume_factor(i_cell), coll_idx,
-                        cell_start_pair, p_mask, p_pair_indices_1, p_pair_indices_2,
+                        cell_start_pair, global_index, p_mask, p_pair_indices_1, p_pair_indices_2,
                         p_pair_reaction_weight, p_product_data, engine);
                 }
             );
diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H
index caefda2a027..c0e5c9b8a6d 100644
--- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H
@@ -8,6 +8,7 @@
 #ifndef WARPX_BREMSSTRAHLUNG_FUNC_H_
 #define WARPX_BREMSSTRAHLUNG_FUNC_H_
 
+#include "Particles/Collision/CollisionFuncBase.H"
 #include "Particles/Collision/BinaryCollision/BinaryCollisionUtils.H"
 #include "Particles/Pusher/GetAndSetPosition.H"
 #include "Particles/MultiParticleContainer.H"
@@ -31,6 +32,7 @@
  *  effectively created in the particle creation functor.
  */
 class BremsstrahlungFunc
+    : public CollisionFuncBase
 {
     // Define shortcuts for frequently-used type names
     using ParticleType = WarpXParticleContainer::ParticleType;
@@ -92,7 +94,8 @@ public:
             amrex::ParticleReal const  /*q1*/, amrex::ParticleReal const  /*q2*/,
             amrex::ParticleReal const m1, amrex::ParticleReal const m2,
             amrex::Real const dt, amrex::Real const /*dV*/, index_type coll_idx,
-            index_type const cell_start_pair, index_type* AMREX_RESTRICT p_mask,
+            index_type const cell_start_pair, amrex::IntVect const /*global_index*/,
+            index_type* AMREX_RESTRICT p_mask,
             index_type* AMREX_RESTRICT p_pair_indices_1, index_type* AMREX_RESTRICT p_pair_indices_2,
             amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight,
             amrex::ParticleReal* AMREX_RESTRICT p_product_data,
@@ -406,7 +409,7 @@ public:
 
     };
 
-    [[nodiscard]] Executor const& executor () const { return m_exe; }
+    [[nodiscard]] Executor const& executor (int const /*lev*/, amrex::MFIter const& /*mfi*/) const { return m_exe; }
 
     [[nodiscard]] bool use_global_debye_length() const { return m_use_global_debye_length; }
 
diff --git a/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H b/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H
index 3f516d8969e..382ce1a3c36 100644
--- a/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H
@@ -8,6 +8,7 @@
 #ifndef WARPX_PAIRWISE_COULOMB_COLLISION_FUNC_H_
 #define WARPX_PAIRWISE_COULOMB_COLLISION_FUNC_H_
 
+#include "Particles/Collision/CollisionFuncBase.H"
 #include "ElasticCollisionPerez.H"
 #include "Particles/Pusher/GetAndSetPosition.H"
 #include "Particles/WarpXParticleContainer.H"
@@ -24,6 +25,7 @@
  *        ElasticCollisionPerez. It also reads and contains the Coulomb logarithm.
  */
 class PairWiseCoulombCollisionFunc
+    : public CollisionFuncBase
 {
     // Define shortcuts for frequently-used type names
     using ParticleType = WarpXParticleContainer::ParticleType;
@@ -120,7 +122,8 @@ public:
             amrex::ParticleReal const  q1, amrex::ParticleReal const  q2,
             amrex::ParticleReal const  m1, amrex::ParticleReal const  m2,
             amrex::Real const  dt, amrex::Real const dV, index_type coll_idx,
-            index_type const /*cell_start_pair*/, index_type* /*p_mask*/,
+            index_type const /*cell_start_pair*/, amrex::IntVect const /*global_index*/,
+            index_type* /*p_mask*/,
             index_type* /*p_pair_indices_1*/, index_type* /*p_pair_indices_2*/,
             amrex::ParticleReal* /*p_pair_reaction_weight*/,
             amrex::ParticleReal* /*p_product_data*/,
@@ -144,7 +147,7 @@ public:
         bool m_isSameSpecies;
     };
 
-    [[nodiscard]] Executor const& executor () const { return m_exe; }
+    [[nodiscard]] Executor const& executor (int const /*lev*/, amrex::MFIter const& /*mfi*/) const { return m_exe; }
 
     [[nodiscard]] bool use_global_debye_length() const { return m_use_global_debye_length; }
 
diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H
index 546e4ce8a23..83a513842f3 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H
@@ -12,6 +12,7 @@
 
 #include "CollisionFilterFunc.H"
 
+#include "Particles/Collision/CollisionFuncBase.H"
 #include "Particles/Collision/BinaryCollision/BinaryCollisionUtils.H"
 #include "Particles/Collision/BinaryCollision/ParticleShufflers.H"
 #include "Particles/Collision/CollisionBase.H"
@@ -35,6 +36,7 @@
  * used for binary Coulomb collisions and the nuclear fusion module.
  */
 class DSMCFunc
+    : public CollisionFuncBase
 {
     // Define shortcuts for frequently-used type names
     using ParticleType = WarpXParticleContainer::ParticleType;
@@ -103,7 +105,8 @@ public:
             amrex::ParticleReal const  /*q1*/, amrex::ParticleReal const  /*q2*/,
             amrex::ParticleReal const  m1, amrex::ParticleReal const  m2,
             amrex::Real const  dt, amrex::Real const dV, index_type coll_idx,
-            index_type const cell_start_pair, index_type* AMREX_RESTRICT p_mask,
+            index_type const cell_start_pair, amrex::IntVect const /*global_index*/,
+            index_type* AMREX_RESTRICT p_mask,
             index_type* AMREX_RESTRICT p_pair_indices_1, index_type* AMREX_RESTRICT p_pair_indices_2,
             amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight,
             amrex::ParticleReal* /*p_product_data*/,
@@ -183,7 +186,7 @@ public:
         ScatteringProcess::Executor* m_scattering_processes_data;
     };
 
-    [[nodiscard]] Executor const& executor () const { return m_exe; }
+    [[nodiscard]] Executor const& executor (int const /*lev*/, amrex::MFIter const& /*mfi*/) const { return m_exe; }
 
     bool use_global_debye_length() { return false; }
 
diff --git a/Source/Particles/Collision/BinaryCollision/LinearBreitWheeler/LinearBreitWheelerCollisionFunc.H b/Source/Particles/Collision/BinaryCollision/LinearBreitWheeler/LinearBreitWheelerCollisionFunc.H
index a7b006048eb..e5447429a6f 100644
--- a/Source/Particles/Collision/BinaryCollision/LinearBreitWheeler/LinearBreitWheelerCollisionFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/LinearBreitWheeler/LinearBreitWheelerCollisionFunc.H
@@ -10,6 +10,7 @@
 
 #include "SingleLinearBreitWheelerCollisionEvent.H"
 
+#include "Particles/Collision/CollisionFuncBase.H"
 #include "Particles/Collision/BinaryCollision/BinaryCollisionUtils.H"
 #include "Particles/Pusher/GetAndSetPosition.H"
 #include "Particles/MultiParticleContainer.H"
@@ -33,6 +34,7 @@
  *  This functor also reads and stores the event multiplier.
  */
 class LinearBreitWheelerCollisionFunc
+    : public CollisionFuncBase
 {
     // Define shortcuts for frequently-used type names
     using ParticleType = WarpXParticleContainer::ParticleType;
@@ -146,7 +148,8 @@ public:
             amrex::ParticleReal const  /*q1*/, amrex::ParticleReal const  /*q2*/,
             amrex::ParticleReal const  /*m1*/, amrex::ParticleReal const  /*m2*/,
             amrex::Real const  dt, amrex::Real const dV, index_type coll_idx,
-            index_type const cell_start_pair, index_type* AMREX_RESTRICT p_mask,
+            index_type const cell_start_pair, amrex::IntVect const /*global_index*/,
+            index_type* AMREX_RESTRICT p_mask,
             index_type* AMREX_RESTRICT p_pair_indices_1, index_type* AMREX_RESTRICT p_pair_indices_2,
             amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight,
             amrex::ParticleReal* /*p_product_data*/,
@@ -224,7 +227,7 @@ public:
         bool m_need_product_data = false;
     };
 
-    [[nodiscard]] Executor const& executor () const { return m_exe; }
+    [[nodiscard]] Executor const& executor (int const /*lev*/, amrex::MFIter const& /*mfi*/) const { return m_exe; }
 
     bool use_global_debye_length() { return false; }
 
diff --git a/Source/Particles/Collision/BinaryCollision/LinearCompton/LinearComptonCollisionFunc.H b/Source/Particles/Collision/BinaryCollision/LinearCompton/LinearComptonCollisionFunc.H
index 0f2e24231f9..9b87953da3d 100644
--- a/Source/Particles/Collision/BinaryCollision/LinearCompton/LinearComptonCollisionFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/LinearCompton/LinearComptonCollisionFunc.H
@@ -10,6 +10,7 @@
 
 #include "SingleLinearComptonCollisionEvent.H"
 
+#include "Particles/Collision/CollisionFuncBase.H"
 #include "Particles/Collision/BinaryCollision/BinaryCollisionUtils.H"
 #include "Particles/Pusher/GetAndSetPosition.H"
 #include "Particles/MultiParticleContainer.H"
@@ -33,6 +34,7 @@
  *  This functor also reads and stores the event multiplier.
  */
 class LinearComptonCollisionFunc
+    : public CollisionFuncBase
 {
     // Define shortcuts for frequently-used type names
     using ParticleType = WarpXParticleContainer::ParticleType;
@@ -142,7 +144,8 @@ public:
             amrex::ParticleReal const  /*q1*/, amrex::ParticleReal const  /*q2*/,
             amrex::ParticleReal const  /*m1*/, amrex::ParticleReal const  /*m2*/,
             amrex::Real const  dt, amrex::Real const dV, index_type coll_idx,
-            index_type const cell_start_pair, index_type* AMREX_RESTRICT p_mask,
+            index_type const cell_start_pair, amrex::IntVect const /*global_index*/,
+            index_type* AMREX_RESTRICT p_mask,
             index_type* AMREX_RESTRICT p_pair_indices_1, index_type* AMREX_RESTRICT p_pair_indices_2,
             amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight,
             amrex::ParticleReal* /*p_product_data*/,
@@ -220,7 +223,7 @@ public:
         bool m_need_product_data = false;
     };
 
-    [[nodiscard]] Executor const& executor () const { return m_exe; }
+    [[nodiscard]] Executor const& executor (int const /*lev*/, amrex::MFIter const& /*mfi*/) const { return m_exe; }
 
     bool use_global_debye_length() { return false; }
 
diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H
index 4e2dd84b885..21fec12ff54 100644
--- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H
@@ -10,6 +10,7 @@
 
 #include "SingleNuclearFusionEvent.H"
 
+#include "Particles/Collision/CollisionFuncBase.H"
 #include "Particles/Collision/BinaryCollision/BinaryCollisionUtils.H"
 #include "Particles/Pusher/GetAndSetPosition.H"
 #include "Particles/MultiParticleContainer.H"
@@ -33,6 +34,7 @@
  *  This functor also reads and contains the fusion multiplier.
  */
 class NuclearFusionFunc
+    : public CollisionFuncBase
 {
     // Define shortcuts for frequently-used type names
     using ParticleType = WarpXParticleContainer::ParticleType;
@@ -77,11 +79,41 @@ public:
             pp_collision_name, "probability_target_value",
             m_probability_target_value);
 
+        bool create_products = true;
+        utils::parser::queryWithParser(
+            pp_collision_name, "create_products", create_products);
+        utils::parser::queryWithParser(
+            pp_collision_name, "save_particle_production", m_save_particle_production);
+
         m_exe.m_fusion_multiplier = m_fusion_multiplier;
         m_exe.m_probability_threshold = m_probability_threshold;
         m_exe.m_probability_target_value = m_probability_target_value;
         m_exe.m_fusion_type = m_fusion_type;
         m_exe.m_isSameSpecies = m_isSameSpecies;
+        m_exe.m_create_products = create_products;
+
+        if (m_save_particle_production) {
+            m_particle_production_mf_name = collision_name + "_particle_production";
+        }
+    }
+
+    void AllocData () override {
+        if (m_save_particle_production) {
+            WarpX & warpx = WarpX::GetInstance();
+            for (int lev = 0; lev <= warpx.finestLevel(); ++lev) {
+
+                amrex::BoxArray const & ba = warpx.boxArray(lev);
+                amrex::DistributionMapping const & dmap = warpx.DistributionMap(lev);
+                int const ncomps = 1;
+                amrex::IntVect const ng = amrex::IntVect::TheZeroVector();
+                amrex::Real const initial_value = 0.;
+                bool const remake = true;
+                bool const redistribute_on_remake = true;
+                bool const checkpoint_restart = true;
+                warpx.m_fields.alloc_init(m_particle_production_mf_name, lev, ba, dmap, ncomps, ng,
+                                          initial_value, remake, redistribute_on_remake, checkpoint_restart);
+            }
+        }
     }
 
     struct Executor {
@@ -111,6 +143,7 @@ public:
          * @param[in] dV is the volume of the corresponding cell.
          * @param[in] coll_idx is the collision index offset.
          * @param[in] cell_start_pair is the start index of the pairs in that cell.
+         * @param[in] global_index grid cell where collision is taking place
          * @param[out] p_mask is a mask that will be set to true if a fusion event occurs for a given
          * pair. It is only needed here to store information that will be used later on when actually
          * creating the product particles.
@@ -136,7 +169,8 @@ public:
             amrex::ParticleReal const  /*q1*/, amrex::ParticleReal const  /*q2*/,
             amrex::ParticleReal const  m1, amrex::ParticleReal const  m2,
             amrex::Real const  dt, amrex::Real const dV, index_type coll_idx,
-            index_type const cell_start_pair, index_type* AMREX_RESTRICT p_mask,
+            index_type const cell_start_pair, amrex::IntVect const global_index,
+            index_type* AMREX_RESTRICT p_mask,
             index_type* AMREX_RESTRICT p_pair_indices_1, index_type* AMREX_RESTRICT p_pair_indices_2,
             amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight,
             amrex::ParticleReal* /*p_product_data*/,
@@ -189,7 +223,10 @@ public:
                         m_fusion_multiplier, multiplier_ratio,
                         m_probability_threshold,
                         m_probability_target_value,
-                        m_fusion_type, engine);
+                        m_fusion_type, engine,
+                        m_create_products,
+                        global_index,
+                        m_particle_production);
 
                     // Remove pair reaction weight from the colliding particles' weights
                     if (p_mask[pair_index]) {
@@ -217,9 +254,22 @@ public:
         bool m_computeSpeciesTemperatures = false;
         bool m_need_product_data = false;
         bool m_isSameSpecies;
+
+        bool m_create_products = true;
+        amrex::Array4<amrex::Real> m_particle_production;
     };
 
-    [[nodiscard]] Executor const& executor () const { return m_exe; }
+    [[nodiscard]] Executor executor (int const lev, amrex::MFIter const& mfi) const {
+        // Note that a copy is needed since with m_save_particle_production, it will be modified
+        // in a way that is thread dependent.
+        Executor exe = m_exe;
+        if (m_save_particle_production) {
+            WarpX & warpx = WarpX::GetInstance();
+            amrex::MultiFab * particle_production_mf = warpx.m_fields.get(m_particle_production_mf_name, lev);
+            exe.m_particle_production = particle_production_mf->array(mfi);
+        }
+        return exe;
+    }
 
     bool use_global_debye_length() { return false; }
 
@@ -238,6 +288,9 @@ private:
     NuclearFusionType m_fusion_type;
     bool m_isSameSpecies;
 
+    bool m_save_particle_production = false;
+    std::string m_particle_production_mf_name;
+
     Executor m_exe;
 };
 
diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H
index b0fe66da37e..1633e288160 100644
--- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H
+++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H
@@ -48,6 +48,9 @@
  * to determine by how much the fusion multiplier is reduced
  * @param[in] fusion_type the physical fusion process to model
  * @param[in] engine the random engine.
+ * @param[in] create_products flags whether or not to create particles
+ * @param[in] global_index the grid cell where the collision is taking place
+ * @param[in] particle_production the particle production result array
  */
 template <typename index_type>
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
@@ -64,7 +67,10 @@ void SingleNuclearFusionEvent (const amrex::ParticleReal& u1x, const amrex::Part
                                const amrex::ParticleReal& probability_threshold,
                                const amrex::ParticleReal& probability_target_value,
                                const NuclearFusionType& fusion_type,
-                               const amrex::RandomEngine& engine)
+                               const amrex::RandomEngine& engine,
+                               const bool create_products,
+                               const amrex::IntVect global_index,
+                               const amrex::Array4<amrex::Real> & particle_production)
 {
     amrex::ParticleReal E_coll, v_coll, lab_to_COM_factor;
 
@@ -114,14 +120,23 @@ void SingleNuclearFusionEvent (const amrex::ParticleReal& u1x, const amrex::Part
     // std::expm1 is used since it maintains correctness for small exponent.
     const amrex::ParticleReal probability = -std::expm1(-probability_estimate);
 
+    const amrex::Real w_new = w_min/fusion_multiplier_eff;
+
+    // Save the particle production density if requested
+    if (particle_production.ok()) {
+        const amrex::Real new_products = probability*w_new/dV;
+        amrex::Gpu::Atomic::AddNoRet(&particle_production(global_index), new_products);
+
+    }
+
     // Get a random number
     const amrex::ParticleReal random_number = amrex::Random(engine);
 
-    // If we have a fusion event, set the mask the true and fill the product weight array
-    if (random_number < probability)
+    // If we have a fusion event and are creating products, set the mask the true and fill the product weight array
+    if (random_number < probability && create_products)
     {
         p_mask[pair_index] = true;
-        p_pair_reaction_weight[pair_index] = w_min/fusion_multiplier_eff;
+        p_pair_reaction_weight[pair_index] = w_new;
     }
     else
     {
diff --git a/Source/Particles/Collision/CollisionBase.H b/Source/Particles/Collision/CollisionBase.H
index de674cb6cd1..f5d8998b5c7 100644
--- a/Source/Particles/Collision/CollisionBase.H
+++ b/Source/Particles/Collision/CollisionBase.H
@@ -22,6 +22,8 @@ public:
 
     explicit CollisionBase (const std::string& collision_name);
 
+    virtual void AllocData () {}
+
     virtual void doCollisions (amrex::Real /*cur_time*/, amrex::Real /*dt*/, MultiParticleContainer* /*mypc*/ ){}
 
     CollisionBase(CollisionBase const &) = delete;
diff --git a/Source/Particles/Collision/CollisionFuncBase.H b/Source/Particles/Collision/CollisionFuncBase.H
new file mode 100644
index 00000000000..8ce6704bac8
--- /dev/null
+++ b/Source/Particles/Collision/CollisionFuncBase.H
@@ -0,0 +1,29 @@
+/* Copyright 2026 The WarpX Community
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef WARPX_COLLISION_FUNC_BASE_H_
+#define WARPX_COLLISION_FUNC_BASE_H_
+
+class CollisionFuncBase
+{
+public:
+
+    CollisionFuncBase () = default;
+    virtual ~CollisionFuncBase () = default;
+
+    CollisionFuncBase (const CollisionFuncBase&) = default;
+    CollisionFuncBase& operator= (const CollisionFuncBase&) = default;
+
+    CollisionFuncBase (CollisionFuncBase&&) noexcept = default;
+    CollisionFuncBase& operator= (CollisionFuncBase&&) noexcept = default;
+
+    //! Optional hook to allocate additional data structures.
+    virtual void AllocData () {}
+
+};
+
+#endif /* WARPX_COLLISION_FUNC_BASE_H_ */
diff --git a/Source/Particles/Collision/CollisionHandler.H b/Source/Particles/Collision/CollisionHandler.H
index 03d68df99c8..25281508cfe 100644
--- a/Source/Particles/Collision/CollisionHandler.H
+++ b/Source/Particles/Collision/CollisionHandler.H
@@ -25,6 +25,9 @@ class CollisionHandler
 public:
     explicit CollisionHandler (const MultiParticleContainer*  mypc);
 
+    /* Allocate data needed for collision */
+    void AllocData ();
+
     /* Perform all of the collisions */
     void doCollisions (int step, amrex::Real cur_time, amrex::Real dt, MultiParticleContainer* mypc);
 
diff --git a/Source/Particles/Collision/CollisionHandler.cpp b/Source/Particles/Collision/CollisionHandler.cpp
index 1301b819f92..dedceccf4b0 100644
--- a/Source/Particles/Collision/CollisionHandler.cpp
+++ b/Source/Particles/Collision/CollisionHandler.cpp
@@ -111,6 +111,14 @@ CollisionHandler::CollisionHandler(MultiParticleContainer const * const mypc)
 
 }
 
+/* \brief Allocate any data needed for the collision */
+void CollisionHandler::AllocData ()
+{
+    for (auto& collision : allcollisions) {
+        collision->AllocData();
+    }
+}
+
 /** Perform all collisions
  *
  * @param step Current iteration
diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp
index e32fcc53c39..4100a69f8fe 100644
--- a/Source/Particles/MultiParticleContainer.cpp
+++ b/Source/Particles/MultiParticleContainer.cpp
@@ -435,6 +435,8 @@ MultiParticleContainer::AllocData ()
     for (auto& pc : allcontainers) {
         pc->AllocData();
     }
+
+    collisionhandler->AllocData();
 }
 
 void