diff --git a/examples/1_Microphysics_Simulation.ipynb b/examples/1_Microphysics_Simulation.ipynb
index d4a8bc14..69be5b50 100644
--- a/examples/1_Microphysics_Simulation.ipynb
+++ b/examples/1_Microphysics_Simulation.ipynb
@@ -87,15 +87,14 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "### Clustering and Volume Cuts\n",
+    "## Clustering, Physics, and Cuts\n",
     "\n",
-    "In the next step, all interactions with the same `cluster_ids` are merged. The energy of the interactions is summed up and the position and time is calculated as the weighted average of the positions of the individual interactions. The interaction type of the cluster is determined either by the interaction with the highest energy or by the first interaction in the cluster. The interaction type is later used to choose the correct emmision model. \n",
+    "### Clustering\n",
     "\n",
-    "Following the clustering, the `VolumeProperties` plugin is used to assign physical properties to different detector regions, like the xenon density and the ability to create S2s. The plugin also assigns a `vol_id` to each interaction based on its volume. The volume definitions are stored in a dictionary in fuse.common. \n",
+    "In the first step, all interactions with the same `cluster_id` are merged.  \n",
+    "The cluster energy is computed as the sum of the individual interaction energies, while the cluster position and time are calculated as energy-weighted averages of the contributing interactions.\n",
     "\n",
-    "After the volume properties have been assigned, the `VolumeSelection` plugin is used to select only interactions in the detector regions of interest. By default, these are the TPC and the region below the cathode. \n",
-    "\n",
-    "The next step is the `SelectionMerger` plugin. This plugin merges all interactions that pass the selection into a new data set, the `interactions_in_roi`. By default, the only selection applied is the volume selection, as defined in the `DefaultSimulation` plugin. However, more complex selections can be defined, for example using the `EnergyCut` plugin to select only events with a certain energy range. The selection logic can be used by registring - as an example - the `LowEnergySimulation` plugin, which inherits from `SelectionMerger` and applies both the volume selection and an energy cut defined in the `EnergyCut` plugin."
+    "The interaction type of a cluster is derived from the constituent interactions (for example from the interaction with the highest energy or from the first interaction in the cluster). This interaction type is later used to select the appropriate emission model.\n"
    ]
   },
   {
@@ -104,18 +103,18 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "st.make(run_number, \"clustered_interactions\")\n",
-    "st.make(run_number, \"volume_properties\")\n",
-    "st.make(run_number, \"interactions_in_roi\")"
+    "st.make(run_number, \"clustered_interactions\")"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "### Electric Field and Emission Model\n",
+    "### Volume Properties\n",
+    "\n",
+    "After clustering, the VolumeProperties plugin assigns physical properties to each interaction based on its location in the detector. These properties include, for example, the local xenon density and whether charge extraction and S2 production are possible.\n",
     "\n",
-    "The aim of this simulation part is to model the scintillation and ionization processes at the interaction site. First we need to estimate the electric field strength at the interaction position. This is done in the `ElectricField` plugin using a simulated field map. The field values can be accessed using the target `electric_field_values`. Next we can estimate the number of produced electrons and photons using an emission model. The default implementation of fuse uses the `NestYields` plugin where `nestpy` is used. fuse also provides alternative plugins where the yields are calculated using BBF or a beta response model. These plugins should only be used if you know what you are doing. The result of the emission model can be accessed using the target `quanta`. "
+    "In addition, a vol_id is assigned to each interaction, identifying the detector volume it belongs to. The volume definitions are stored in a dictionary in fuse.common."
    ]
   },
   {
@@ -124,15 +123,22 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "st.make(run_number, \"electric_field_values\")\n",
-    "st.make(run_number, \"quanta\")"
+    "st.make(run_number, \"volume_properties\")"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Finally we can collect the simulation results of the last few steps using the `MicroPhysicsSummary` plugin. This plugin is a `strax.MergeOnlyPlugin` and does not do any calculations. It just merges the results of the previous plugins and can be accessed using the target `microphysics_summary`."
+    "### Electric Field and Yields\n",
+    "\n",
+    "The next step is to model scintillation and ionization at the interaction site.\n",
+    "\n",
+    "First, the local electric field strength is estimated using a simulated field map in the ElectricField plugin. The resulting field values can be accessed via the electric_field_values target.\n",
+    "\n",
+    "Next, an emission model is used to estimate the number of produced photons and electrons. By default, fuse uses the NestYields plugin, which relies on nestpy. Alternative emission models (e.g. BBF or beta-response models) are also available, but should only be used with care.\n",
+    "\n",
+    "The output of the emission model is stored in the quanta target."
    ]
   },
   {
@@ -141,9 +147,31 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "st.make(run_number, \"microphysics_summary\")\n",
+    "st.make(run_number, \"electric_field_values\")\n",
+    "st.make(run_number, \"quanta\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Cuts and Selections\n",
+    "\n",
+    "After all physical quantities have been computed, selection criteria are evaluated.\n",
+    "The most basic selection is the VolumeSelection, which flags interactions located in detector regions of interest (by default the TPC and the region below the cathode).\n",
     "\n",
-    "microphysics_summary = st.get_df(run_number, [\"microphysics_summary\"])"
+    "Additional cuts can optionally be applied, such as an EnergyCut or an NRCut, depending on the simulation configuration. These cut plugins do not modify the data themselves, but instead produce boolean masks indicating which interactions pass each selection."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### MicroPhysicsSummary\n",
+    "\n",
+    "Finally, the MicroPhysicsSummary plugin collects the results of the previous steps. This plugin merges all relevant inputs (clustered interactions, volume properties, electric field values, quanta, and cut results) into a single data set.\n",
+    "\n",
+    "At the same time, it applies the configured selections and cuts, filtering out interactions that do not pass all required criteria. The resulting filtered and merged output can be accessed via the microphysics_summary target."
    ]
   },
   {
@@ -152,7 +180,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "microphysics_summary.head()"
+    "st.make(run_number, \"microphysics_summary\")\n",
+    "microphysics_summary = st.get_df(run_number, \"microphysics_summary\")"
    ]
   }
  ],
diff --git a/fuse/context.py b/fuse/context.py
index 0767922b..f6fd86b6 100644
--- a/fuse/context.py
+++ b/fuse/context.py
@@ -43,11 +43,10 @@
 # Plugins to simulate microphysics (remaining)
 microphysics_plugins_remaining = [
     fuse.micro_physics.cuts_and_selections.VolumeSelection,
-    fuse.micro_physics.cuts_and_selections.DefaultSimulation,
     fuse.micro_physics.ElectricField,
     fuse.micro_physics.NestYields,
-    fuse.micro_physics.MicroPhysicsSummary,
     fuse.micro_physics.VolumeProperties,
+    fuse.micro_physics.DefaultSimulation,
 ]
 
 # Plugins to simulate S1 signals
diff --git a/fuse/plugins/micro_physics/cuts_and_selections/__init__.py b/fuse/plugins/micro_physics/cuts_and_selections/__init__.py
index 0eff7eb5..a4e7593f 100644
--- a/fuse/plugins/micro_physics/cuts_and_selections/__init__.py
+++ b/fuse/plugins/micro_physics/cuts_and_selections/__init__.py
@@ -1,6 +1,3 @@
-from . import apply_selections
-from .apply_selections import *
-
 from . import detector_volumes
 from .detector_volumes import *
 
diff --git a/fuse/plugins/micro_physics/cuts_and_selections/apply_selections.py b/fuse/plugins/micro_physics/cuts_and_selections/apply_selections.py
deleted file mode 100644
index b497ed1c..00000000
--- a/fuse/plugins/micro_physics/cuts_and_selections/apply_selections.py
+++ /dev/null
@@ -1,91 +0,0 @@
-import strax
-import numpy as np
-
-from ....plugin import FuseBasePlugin
-
-from ....dtypes import (
-    primary_positions_fields,
-    cluster_positions_fields,
-    cluster_id_fields,
-    cluster_misc_fields,
-    volume_properties_fields,
-)
-
-export, __all__ = strax.exporter()
-
-
-@export
-class SelectionMerger(FuseBasePlugin):
-    """Merge cuts/selections and stamp per-volume constants.
-
-    The selection logic is given as a string expression over boolean fields
-    in the `clustered_interactions` data. The expression may use '&', '|', '~',
-    and parentheses. For example, to select interactions in the fiducial volume
-    and with energy between 1 and 10 keV, use:
-        "volume_selection & energy_range_cut"
-    """
-
-    __version__ = "1.0.0"
-    save_when = strax.SaveWhen.TARGET
-    provides = "interactions_in_roi"
-    data_kind = "interactions_in_roi"
-
-    selection_logic = "volume_selection"
-
-    dtype = (
-        cluster_positions_fields
-        + cluster_id_fields
-        + cluster_misc_fields
-        + primary_positions_fields
-        + volume_properties_fields
-        + strax.time_fields
-    )
-
-    @staticmethod
-    def _eval_logic(arr, expr: str) -> np.ndarray:
-        """Tiny evaluator for '&', '|', '~', and parentheses over boolean
-        fields."""
-        # Build an environment from available boolean-likes.
-        # (We expose ALL fields; the expression must reference valid names.)
-        env = {name: arr[name].astype(bool, copy=False) for name in arr.dtype.names}
-        # Safe eval: strip builtins.
-        out = eval(expr, {"__builtins__": None}, env)
-        # Normalize to a 1D boolean numpy array
-        return np.asarray(out, dtype=np.bool_)
-
-    def compute(self, clustered_interactions):
-        if len(clustered_interactions) == 0:
-            return np.zeros(0, dtype=self.dtype)
-
-        mask = self._eval_logic(clustered_interactions, self.selection_logic)
-        if not mask.any():
-            return np.zeros(0, dtype=self.dtype)
-
-        # Filter and project to final dtype (drops predicate fields automatically)
-        reduced = clustered_interactions[mask]
-        out = np.zeros((len(reduced),), dtype=self.dtype)
-        strax.copy_to_buffer(reduced, out, "_copy_selected")
-        return out
-
-
-@export
-class DefaultSimulation(SelectionMerger):
-    depends_on = (
-        "clustered_interactions",
-        "volume_properties",
-        "volume_selection",
-    )
-    __version__ = "1.0.2"
-    selection_logic = "volume_selection"
-
-
-@export
-class LowEnergySimulation(SelectionMerger):
-    depends_on = (
-        "clustered_interactions",
-        "volume_properties",
-        "volume_selection",
-        "energy_range_cut",
-    )
-    __version__ = "1.0.1"
-    selection_logic = "volume_selection & energy_range_cut"
diff --git a/fuse/plugins/micro_physics/cuts_and_selections/detector_volumes.py b/fuse/plugins/micro_physics/cuts_and_selections/detector_volumes.py
index f1ce3cec..4fa7798f 100644
--- a/fuse/plugins/micro_physics/cuts_and_selections/detector_volumes.py
+++ b/fuse/plugins/micro_physics/cuts_and_selections/detector_volumes.py
@@ -7,7 +7,10 @@
 logging.basicConfig(handlers=[logging.StreamHandler()])
 log = logging.getLogger("fuse.micro_physics.detector_volumes")
 
+export, __all__ = strax.exporter()
 
+
+@export
 class VolumeSelection(strax.CutPlugin):
     """Plugin that evaluates if interactions are in a defined detector
     volume."""
diff --git a/fuse/plugins/micro_physics/cuts_and_selections/physics_cases.py b/fuse/plugins/micro_physics/cuts_and_selections/physics_cases.py
index 025576a3..88e68380 100644
--- a/fuse/plugins/micro_physics/cuts_and_selections/physics_cases.py
+++ b/fuse/plugins/micro_physics/cuts_and_selections/physics_cases.py
@@ -1,8 +1,12 @@
 import straxen
 import strax
 import numpy as np
+import numba
 
+export, __all__ = strax.exporter()
 
+
+@export
 class EnergyCut(strax.CutPlugin):
     """Plugin evaluates if the sum of the events energy is below a
     threshold."""
@@ -67,3 +71,122 @@ def group_interaction_energies_by_event_id(energy, event_id):
 
     unique_event_id, split_position = np.unique(event_id_sorted, return_index=True)
     return np.split(energy_sorted, split_position[1:]), unique_event_id
+
+
+@export
+class NRCut(strax.CutPlugin):
+    """Plugin which filters the microphysics summary for valid NR events."""
+
+    depends_on = ["clustered_interactions", "quanta"]
+
+    __version__ = "0.0.2"
+
+    provides = "nr_cut"
+    cut_name = "nr_cut"
+    cut_description = "Selects valid NR events in microphysics summary"
+
+    g1_photon_yield = straxen.URLConfig(
+        default=0.1,
+        type=(int, float),
+        help="Scaled g1 x 0.8 to account for corrections [pe/ph]",
+    )
+    g2_electron_yield = straxen.URLConfig(
+        default=13.4,
+        type=(int, float),
+        help="Scaled g2 x 0.8 to account for corrections [pe/e]",
+    )
+
+    max_s1_area = straxen.URLConfig(
+        default=700,
+        type=(int, float),
+        help="Max S1 area [pe] for NR roi",
+    )
+    max_s2_area = straxen.URLConfig(
+        default=3 * 10**4,
+        type=(int, float),
+        help="Max S2 area [pe] for NR roi",
+    )
+
+    def cut_by(self, clustered_interactions):
+
+        vertex_to_keep = filter_events(
+            clustered_interactions,
+            self.g1_photon_yield,
+            self.g2_electron_yield,
+            self.max_s1_area,
+            self.max_s2_area,
+        )
+
+        mask = vertex_to_keep.astype(bool)
+
+        self.log.info(f"Keeping {np.sum(mask)} out of {len(mask)} interactions after NR cut")
+
+        return mask
+
+
+@numba.njit
+def filter_events(mps, g1, g2, max_s1, max_s2):
+    """Small function to filter microphysics for valid NR events.
+
+    We cut all events which are overshadowed by other events outside of
+    our ROI excluding delayed deexcitations. To account for missing
+    detector corrections one should scale g1/g2 down.
+    """
+
+    event_ids = np.unique(mps["eventid"])
+
+    vertex_to_keep = np.ones(len(mps))
+
+    vertex_i = 0
+
+    for event_i in event_ids:
+        start_index = vertex_i
+        max_photons = 0
+        max_electrons = 0
+        prompt_photons = 0
+        number_of_nr_interactions = 0
+        start_time = mps["time"][vertex_i]
+        loop_broke = False
+
+        for vertex_i in range(start_index, len(mps)):
+            vertex = mps[vertex_i]
+
+            _is_a_new_event = event_i < vertex["eventid"]
+            if _is_a_new_event:
+                # Next event starts break for loop and check next event
+                loop_broke = True
+                break
+
+            # Is prompt vertex:
+            _is_prompt = (vertex["time"] - start_time) < 200  # ns
+            if _is_prompt:
+                prompt_photons += vertex["photons"]
+
+            # Ignore and drop vertex if too much delayed within event:
+            _vertex_is_delayed = (vertex["time"] - start_time) > 3_000_000  # ns (3 ms)
+            if _vertex_is_delayed:
+                vertex_to_keep[vertex_i] = 0
+                continue
+
+            max_photons = max(max_photons, vertex["photons"])
+            max_electrons = max(max_electrons, vertex["electrons"])
+            _is_nr = vertex["nestid"] == 0
+            if _is_nr:
+                number_of_nr_interactions += 1
+
+        # Determine the end index for the current event
+        # If we broke, vertex_i points to the next event, so end is vertex_i
+        # If we didn't break, vertex_i is the last vertex, so end is vertex_i + 1
+        end_index = vertex_i if loop_broke else vertex_i + 1
+
+        # Check if the largest interaction is still within ROI:
+        _is_in_nr_roi = (
+            (max_photons * g1 < max_s1)
+            & (max_electrons * g2 < max_s2)
+            & (number_of_nr_interactions > 0)
+            & (prompt_photons * g1 < max_s1)
+        )
+
+        if not _is_in_nr_roi:
+            vertex_to_keep[start_index:end_index] = 0
+    return vertex_to_keep
diff --git a/fuse/plugins/micro_physics/electric_field.py b/fuse/plugins/micro_physics/electric_field.py
index a447719e..3df0e23e 100644
--- a/fuse/plugins/micro_physics/electric_field.py
+++ b/fuse/plugins/micro_physics/electric_field.py
@@ -15,9 +15,9 @@ class ElectricField(FuseBasePlugin):
 
     __version__ = "0.2.2"
 
-    depends_on = "interactions_in_roi"
+    depends_on = "clustered_interactions"
     provides = "electric_field_values"
-    data_kind = "interactions_in_roi"
+    data_kind = "clustered_interactions"
 
     save_when = strax.SaveWhen.TARGET
 
@@ -34,16 +34,16 @@ class ElectricField(FuseBasePlugin):
         help="Map of the electric field in the detector",
     )
 
-    def compute(self, interactions_in_roi):
-        if len(interactions_in_roi) == 0:
+    def compute(self, clustered_interactions):
+        if len(clustered_interactions) == 0:
             return np.zeros(0, dtype=self.dtype)
 
-        electric_field_array = np.zeros(len(interactions_in_roi), dtype=self.dtype)
-        electric_field_array["time"] = interactions_in_roi["time"]
-        electric_field_array["endtime"] = interactions_in_roi["endtime"]
+        electric_field_array = np.zeros(len(clustered_interactions), dtype=self.dtype)
+        electric_field_array["time"] = clustered_interactions["time"]
+        electric_field_array["endtime"] = clustered_interactions["endtime"]
 
-        r = np.sqrt(interactions_in_roi["x"] ** 2 + interactions_in_roi["y"] ** 2)
-        positions = np.stack((r, interactions_in_roi["z"]), axis=1)
+        r = np.sqrt(clustered_interactions["x"] ** 2 + clustered_interactions["y"] ** 2)
+        positions = np.stack((r, clustered_interactions["z"]), axis=1)
         electric_field_array["e_field"] = self.efield_map(positions)
 
         # Clip negative values to 0
diff --git a/fuse/plugins/micro_physics/microphysics_summary.py b/fuse/plugins/micro_physics/microphysics_summary.py
index 10ae617d..d850b0b6 100644
--- a/fuse/plugins/micro_physics/microphysics_summary.py
+++ b/fuse/plugins/micro_physics/microphysics_summary.py
@@ -1,19 +1,113 @@
 import strax
+import numpy as np
+
+from ...plugin import FuseBasePlugin
+
+from ...dtypes import (
+    primary_positions_fields,
+    cluster_positions_fields,
+    cluster_id_fields,
+    cluster_misc_fields,
+    volume_properties_fields,
+    quanta_fields,
+    electric_fields,
+)
 
 export, __all__ = strax.exporter()
 
 
 @export
-class MicroPhysicsSummary(strax.MergeOnlyPlugin):
-    """MergeOnlyPlugin that summarizes the fuse microphysics simulation results
-    into a single output."""
+class SelectionMerger(FuseBasePlugin):
+    """Merge cuts/selections and stamp per-volume constants.
+
+    The selection logic is given as a string expression over boolean fields
+    in the `clustered_interactions` data. The expression may use '&', '|', '~',
+    and parentheses. For example, to select interactions in the fiducial volume
+    and with energy between 1 and 10 keV, use:
+        "volume_selection & energy_range_cut"
+    """
+
+    __version__ = "1.0.0"
+    save_when = strax.SaveWhen.TARGET
+    provides = "microphysics_summary"
+    data_kind = "interactions_in_roi"
+
+    selection_logic = "volume_selection"
+
+    dtype = (
+        cluster_positions_fields
+        + cluster_id_fields
+        + cluster_misc_fields
+        + primary_positions_fields
+        + volume_properties_fields
+        + quanta_fields
+        + electric_fields
+        + strax.time_fields
+    )
 
+    @staticmethod
+    def _eval_logic(arr, expr: str) -> np.ndarray:
+        """Tiny evaluator for '&', '|', '~', and parentheses over boolean
+        fields."""
+        # Build an environment from available boolean-likes.
+        # (We expose ALL fields; the expression must reference valid names.)
+        env = {name: arr[name].astype(bool, copy=False) for name in arr.dtype.names}
+        # Safe eval: strip builtins.
+        out = eval(expr, {"__builtins__": None}, env)
+        # Normalize to a 1D boolean numpy array
+        return np.asarray(out, dtype=np.bool_)
+
+    def compute(self, clustered_interactions):
+        if len(clustered_interactions) == 0:
+            return np.zeros(0, dtype=self.dtype)
+
+        mask = self._eval_logic(clustered_interactions, self.selection_logic)
+        if not mask.any():
+            return np.zeros(0, dtype=self.dtype)
+
+        # Filter and project to final dtype (drops predicate fields automatically)
+        reduced = clustered_interactions[mask]
+        out = np.zeros((len(reduced),), dtype=self.dtype)
+        strax.copy_to_buffer(reduced, out, "_copy_selected")
+        return out
+
+
+@export
+class DefaultSimulation(SelectionMerger):
     depends_on = (
-        "interactions_in_roi",
+        "clustered_interactions",
+        "volume_properties",
+        "electric_field_values",
         "quanta",
+        "volume_selection",
+    )
+    __version__ = "1.0.2"
+    selection_logic = "volume_selection"
+
+
+@export
+class LowEnergySimulation(SelectionMerger):
+    depends_on = (
+        "clustered_interactions",
+        "volume_properties",
         "electric_field_values",
+        "quanta",
+        "volume_selection",
+        "energy_range_cut",
     )
-    rechunk_on_save = False
-    save_when = strax.SaveWhen.ALWAYS
-    provides = "microphysics_summary"
-    __version__ = "0.1.0"
+    __version__ = "1.0.1"
+    selection_logic = "volume_selection & energy_range_cut"
+
+
+@export
+class NRSimulation(SelectionMerger):
+    depends_on = (
+        "clustered_interactions",
+        "volume_properties",
+        "electric_field_values",
+        "quanta",
+        "volume_selection",
+        "nr_cut",
+    )
+    __version__ = "1.0.1"
+    selection_logic = "volume_selection & nr_cut"
diff --git a/fuse/plugins/micro_physics/yields.py b/fuse/plugins/micro_physics/yields.py
index ec2f0412..60461b13 100644
--- a/fuse/plugins/micro_physics/yields.py
+++ b/fuse/plugins/micro_physics/yields.py
@@ -20,9 +20,13 @@ class NestYields(FuseBasePlugin):
 
     __version__ = "0.3.0"
 
-    depends_on = ("interactions_in_roi", "electric_field_values")
+    depends_on = (
+        "clustered_interactions",
+        "electric_field_values",
+        "volume_properties",
+    )
     provides = "quanta"
-    data_kind = "interactions_in_roi"
+    data_kind = "clustered_interactions"
 
     dtype = quanta_fields + strax.time_fields
 
@@ -140,9 +144,9 @@ def update_nest_width_parameters(self):
 
         return free_parameters
 
-    def compute(self, interactions_in_roi):
+    def compute(self, clustered_interactions):
 
-        if len(interactions_in_roi) == 0:
+        if len(clustered_interactions) == 0:
             return np.zeros(0, dtype=self.dtype)
 
         # set the global nest random generator with self.short_seed
@@ -152,21 +156,21 @@ def compute(self, interactions_in_roi):
         # increment the seed. Next chunk we will use the modified seed to generate random numbers
         self.short_seed += 1
 
-        result = np.zeros(len(interactions_in_roi), dtype=self.dtype)
-        result["time"] = interactions_in_roi["time"]
-        result["endtime"] = interactions_in_roi["endtime"]
+        result = np.zeros(len(clustered_interactions), dtype=self.dtype)
+        result["time"] = clustered_interactions["time"]
+        result["endtime"] = clustered_interactions["endtime"]
 
         # Generate quanta:
-        if len(interactions_in_roi) > 0:
+        if len(clustered_interactions) > 0:
 
             photons, electrons, excitons = self.vectorized_get_quanta(
-                interactions_in_roi["ed"],
-                interactions_in_roi["nestid"],
-                interactions_in_roi["e_field"],
-                interactions_in_roi["A"],
-                interactions_in_roi["Z"],
-                interactions_in_roi["create_S2"],
-                interactions_in_roi["xe_density"],
+                clustered_interactions["ed"],
+                clustered_interactions["nestid"],
+                clustered_interactions["e_field"],
+                clustered_interactions["A"],
+                clustered_interactions["Z"],
+                clustered_interactions["create_S2"],
+                clustered_interactions["xe_density"],
             )
             result["photons"] = photons
             result["electrons"] = electrons
@@ -264,8 +268,9 @@ def process_yields(self, yields_result, create_s2):
 class BBFYields(FuseBasePlugin):
     __version__ = "0.1.1"
 
-    depends_on = ("interactions_in_roi", "electric_field_values")
+    depends_on = ("clustered_interactions", "electric_field_values", "volume_properties")
     provides = "quanta"
+    data_kind = "clustered_interactions"
 
     dtype = quanta_fields + strax.time_fields
 
@@ -274,17 +279,17 @@ def setup(self):
 
         self.bbfyields = BBF_quanta_generator(self.rng)
 
-    def compute(self, interactions_in_roi):
-        result = np.zeros(len(interactions_in_roi), dtype=self.dtype)
-        result["time"] = interactions_in_roi["time"]
-        result["endtime"] = interactions_in_roi["endtime"]
+    def compute(self, clustered_interactions):
+        result = np.zeros(len(clustered_interactions), dtype=self.dtype)
+        result["time"] = clustered_interactions["time"]
+        result["endtime"] = clustered_interactions["endtime"]
 
         # Generate quanta:
-        if len(interactions_in_roi) > 0:
+        if len(clustered_interactions) > 0:
             photons, electrons, excitons = self.bbfyields.get_quanta_vectorized(
-                energy=interactions_in_roi["ed"],
-                interaction=interactions_in_roi["nestid"],
-                field=interactions_in_roi["e_field"],
+                energy=clustered_interactions["ed"],
+                interaction=clustered_interactions["nestid"],
+                field=clustered_interactions["e_field"],
             )
 
             result["photons"] = photons
@@ -434,19 +439,37 @@ def NR_recomb(self, energy, field, par_dict):
 
         return 1.0 - np.log(1.0 + Ni * xi) / (Ni * xi)
 
+    def safe_binomial(self, n, p, *, where=""):
+        import warnings
+
+        # Make p safe for numpy binomial
+        p_safe = np.nan_to_num(p, nan=0.0, posinf=1.0, neginf=0.0)
+        p_safe = np.clip(p_safe, 0.0, 1.0)
+
+        # Make n safe too (binomial needs integer n >= 0)
+        n_safe = np.asarray(n)
+        n_safe = np.where(np.isfinite(n_safe), n_safe, 0)
+        n_safe = np.clip(n_safe, 0, None).astype(np.int64)
+
+        # Warn if we had to change anything (compact, informative)
+        if np.any(p_safe != p) or np.any(n_safe != n):
+            warnings.warn(f"BBF safe_binomial sanitized inputs {where}", RuntimeWarning)
+
+        return self.rng.binomial(n_safe, p_safe)
+
     def get_ER_quanta(self, energy, field, par_dict):
         Nq_mean = energy / par_dict["W"]
         Nq = np.clip(
             np.round(self.rng.normal(Nq_mean, np.sqrt(Nq_mean * par_dict["fano"]))), 0, np.inf
         ).astype(np.int64)
 
-        Ni = self.rng.binomial(Nq, 1.0 / (1.0 + par_dict["Nex/Ni"]))
+        Ni = self.safe_binomial(Nq, 1.0 / (1.0 + par_dict["Nex/Ni"]), where="get_ER_quanta")
 
         recomb = self.ER_recomb(energy, field, par_dict)
         drecomb = self.ER_drecomb(energy, par_dict)
         true_recomb = np.clip(self.rng.normal(recomb, drecomb), 0.0, 1.0)
 
-        Ne = self.rng.binomial(Ni, 1.0 - true_recomb)
+        Ne = self.safe_binomial(Ni, 1.0 - true_recomb, where="get_ER_quanta")
         Nph = Nq - Ne
         Nex = Nq - Ni
         return Nph, Ne, Nex
@@ -458,18 +481,18 @@ def get_NR_quanta(self, energy, field, par_dict):
         )
 
         quenching = self.NR_quenching(energy, par_dict)
-        Nq = self.rng.binomial(Nq, quenching)
+        Nq = self.safe_binomial(Nq, quenching, where="get_NR_quanta")
 
         ExIonRatio = self.NR_ExIonRatio(energy, field, par_dict)
-        Ni = self.rng.binomial(Nq, ExIonRatio / (1.0 + ExIonRatio))
+        Ni = self.safe_binomial(Nq, ExIonRatio / (1.0 + ExIonRatio), where="get_NR_quanta")
 
         penning_quenching = self.NR_Penning_quenching(energy, par_dict)
-        Nex = self.rng.binomial(Nq - Ni, penning_quenching)
+        Nex = self.safe_binomial(Nq - Ni, penning_quenching, where="get_NR_quanta")
 
         recomb = self.NR_recomb(energy, field, par_dict)
         if recomb < 0 or recomb > 1:
             return None, None
 
-        Ne = self.rng.binomial(Ni, 1.0 - recomb)
+        Ne = self.safe_binomial(Ni, 1.0 - recomb, where="get_NR_quanta")
         Nph = Ni + Nex - Ne
         return Nph, Ne, Nex
diff --git a/tests/test_MicroPhysics.py b/tests/test_MicroPhysics.py
index 9dab449e..d877edec 100644
--- a/tests/test_MicroPhysics.py
+++ b/tests/test_MicroPhysics.py
@@ -60,10 +60,6 @@ def test_FindCluster(self):
     def test_MergeCluster(self):
         self.test_context.make(self.run_number, "clustered_interactions")
 
-    @timeout_decorator.timeout(TIMEOUT, exception_message="VolumesMerger timed out")
-    def test_VolumesMerger(self):
-        self.test_context.make(self.run_number, "interactions_in_roi")
-
     @timeout_decorator.timeout(TIMEOUT, exception_message="ElectricField timed out")
     def test_ElectricField(self):
         self.test_context.make(self.run_number, "electric_field_values")
diff --git a/tests/test_MicroPhysics_with_lineage_clustering.py b/tests/test_MicroPhysics_with_lineage_clustering.py
index 97393c13..87721d94 100644
--- a/tests/test_MicroPhysics_with_lineage_clustering.py
+++ b/tests/test_MicroPhysics_with_lineage_clustering.py
@@ -59,10 +59,6 @@ def test_LineageClustering(self):
     def test_MergeLineage(self):
         self.test_context.make(self.run_number, "clustered_interactions")
 
-    @timeout_decorator.timeout(TIMEOUT, exception_message="VolumesMerger timed out")
-    def test_VolumesMerger(self):
-        self.test_context.make(self.run_number, "interactions_in_roi")
-
     @timeout_decorator.timeout(TIMEOUT, exception_message="ElectricField timed out")
     def test_ElectricField(self):
         self.test_context.make(self.run_number, "electric_field_values")