NR filter and edit cuts and selections logic

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\")"
    ]
   }
  ],

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

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 *
 

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."""

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