Add support for circular readout events in the event display

docs/release_notes.rst

@@ -3,18 +3,18 @@
 Release notes
 =============
 
+* Added rough support for DigiEventCircular objects in the event display.
+
 
 Version 0.13.2 (2026-02-04)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-
 * Modified the model to fit the angular coordinate for three-pixel events eta reconstruction.
 * Pull requests merged and issues closed:
 
   - https://github.com/lucabaldini/hexsample/pull/91
 
 
-
 Version 0.13.1 (2026-01-30)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 

src/hexsample/cli.py

@@ -103,6 +103,7 @@ def __init__(self) -> None:
             formatter_class=self._FORMATTER_CLASS)
         self.add_input_file(display)
         self.add_logging_level(display)
+        self.add_display_options(display)
         display.set_defaults(runner=pipeline.display)
 
         # Run the quicklook?
@@ -162,6 +163,13 @@ def add_suffix(parser: argparse.ArgumentParser, default: str) -> None:
         parser.add_argument("--suffix", type=str, default=default,
                             help="suffix for the output file")
 
+    @staticmethod
+    def add_zero_sup_threshold(parser: argparse.ArgumentParser, default: int) -> None:
+        """Add an option for the zero-suppression threshold.
+        """
+        parser.add_argument("--zero_sup_threshold", type=int, default=default,
+                            help="zero-suppression threshold in ADC counts")
+
     @staticmethod
     def add_source_options(parser: argparse.ArgumentParser) -> None:
         """Add an option group for to a given (sub-)parser to define the basic
@@ -246,10 +254,8 @@ def add_readout_options(parser: argparse.ArgumentParser) -> None:
         group.add_argument("--trg_threshold", type=float,
                            default=readout.HexagonalReadoutBase.trg_threshold,
                            help="trigger threshold in electron equivalent")
-        group.add_argument("--zero_sup_threshold", type=int,
-                           default=readout.HexagonalReadoutBase.zero_sup_threshold,
-                           help="zero suppression threshold in ADC counts")
-
+        CliArgumentParser.add_zero_sup_threshold(group,
+                           default=readout.HexagonalReadoutBase.zero_sup_threshold)
 
     def add_recon_options(self, parser: argparse.ArgumentParser) -> None:
         """Add an option group for the reconstruction properties.
@@ -287,6 +293,16 @@ def add_recon_options(self, parser: argparse.ArgumentParser) -> None:
         group.add_argument("--model_path", type=str,
                            help="path of the model to use, in case of custom model")
 
+    def add_display_options(self, parser: argparse.ArgumentParser) -> None:
+        """Add an option group for the event display.
+        """
+        group = parser.add_argument_group("display", "Event display configuration")
+        CliArgumentParser.add_zero_sup_threshold(group,
+                           default=tasks.DisplayDefaults.zero_sup_threshold)
+        group.add_argument("--event_id", type=int,
+                           default=tasks.DisplayDefaults.event_id,
+                           help="ID of the event to display")
+
     def run(self) -> None:
         """Run the actual command tied to the specific options.
         """

src/hexsample/digi.py

@@ -238,6 +238,7 @@ def from_digi(cls, file_row: np.ndarray):
 
     def ascii(self, pha_width: int = 5) -> str:
         """Ascii representation.
+
         In the specific case of this class, the ascii representation is simply a px
         (that is the highest PHA pixel), because the neighbor position is not accessible
         by the DigiEvent.

src/hexsample/display.py

@@ -28,8 +28,11 @@
 from matplotlib.collections import PatchCollection
 from matplotlib.patches import RegularPolygon
 
-from .digi import DigiEventRectangular
+from .clustering import ClusteringNN
+from .digi import DigiEventBase, DigiEventCircular, DigiEventRectangular
 from .hexagon import HexagonalGrid
+from .mc import MonteCarloEvent
+from .readout import HexagonalReadoutBase
 from .roi import RegionOfInterest
 
 
@@ -98,29 +101,6 @@ def show():
         HexagonalGridDisplay.setup_gca()
         plt.show()
 
-    def pha_to_colors(self, pha: np.array, zero_sup_threshold: float = None) -> np.array:
-        """Convert the pha values to colors for display purposes.
-        """
-        values = pha.flatten()
-        values += self.color_map_offset
-        if zero_sup_threshold is not None:
-            values[values <= zero_sup_threshold + self.color_map_offset] = -1.
-        values = values / float(values.max())
-        return self.color_map(values)
-
-    @staticmethod
-    def brightness(color: np.array) -> np.array:
-        """Quick and dirty proxy for the brighness of a given array of colors.
-
-        See https://stackoverflow.com/questions/9733288
-        and also
-        https://stackoverflow.com/questions/30820962
-        for how to split in columns the array of colors.
-        """
-        # pylint: disable = invalid-name
-        r, g, b, _ = color.T
-        return (299 * r + 587 * g + 114 * b) / 1000
-
     def draw(self, offset: Tuple[float, float] = (0., 0.), pixel_labels: bool = False,
         **kwargs) -> HexagonCollection:
         """Draw the full grid display.
@@ -169,7 +149,8 @@ def draw_roi(self, roi: RegionOfInterest, offset: Tuple[float, float] = (0., 0.)
                 plt.text(x + dx - self._grid.pitch, y + dy, f"{row}", **fmt)
         return collection
 
-    def draw_digi_event(self, event: DigiEventRectangular, offset: Tuple[float, float] = (0., 0.),
+    def draw_digi_event_rectangular(self, event: DigiEventRectangular,
+        offset: Tuple[float, float] = (0., 0.),
         indices: bool = True, padding: bool = True, zero_sup_threshold: float = 0,
         values: bool = True, **kwargs) -> HexagonCollection:
         """Draw an actual event int the parent hexagonal grid.
@@ -179,18 +160,84 @@ def draw_digi_event(self, event: DigiEventRectangular, offset: Tuple[float, floa
         """
         # pylint: disable = invalid-name, too-many-arguments, too-many-locals
         collection = self.draw_roi(event.roi, offset, indices, padding, **kwargs)
-        face_color = self.pha_to_colors(event.pha, zero_sup_threshold)
-        collection.set_facecolor(face_color)
         if values:
-            # Draw the pixel values---note that we use black or white for the text
-            # color depending on the brightness of the pixel.
-            black = np.array([0., 0., 0., 1.])
-            white = np.array([1., 1., 1., 1.])
-            text_color = np.tile(black, len(face_color)).reshape(face_color.shape)
-            text_color[self.brightness(face_color) < 0.5] = white
-            fmt = dict(ha="center", va="center", fontsize="xx-small")
-            for x, y, value, color in zip(collection.x, collection.y,\
-                event.pha.flatten(), text_color):
+            # Draw the pixel values
+            fmt = dict(ha="center", va="center", fontsize="small")
+            for x, y, value in zip(collection.x, collection.y, event.pha.flatten()):
                 if value > zero_sup_threshold:
-                    plt.text(x, y, f"{value}", color=color, **fmt)
+                    plt.text(x, y, f"{value}", color="black", **fmt)
         return collection
+
+    def draw_digi_event_circular(self, event: DigiEventCircular,
+        offset: Tuple[float, float] = (0., 0.), zero_sup_threshold: float = 0,
+        values: bool = True, **kwargs) -> HexagonCollection:
+        """Display a digi event with circular readout.
+        """
+        dx, dy = offset
+        # This is shamelessly copied from clustering.py, and we should really
+        # have a function in event that is returning the physical coordinates
+        # and the pha values of all the pixels involved in the event.
+        col = [event.column]
+        row = [event.row]
+        adc_channel_order = [self._grid.adc_channel(event.column, event.row)]
+        # Taking the NN in logical coordinates ...
+        for _col, _row in self._grid.neighbors(event.column, event.row):
+            col.append(_col)
+            row.append(_row)
+            # ... transforming the coordinates of the NN in its corresponding ADC channel ...
+            adc_channel_order.append(self. _grid.adc_channel(_col, _row))
+        # ... reordering the pha array for the correspondence (col[i], row[i]) with pha[i].
+        pha = event.pha[adc_channel_order]
+        # Converting lists into numpy arrays
+        cols = np.array(col)
+        rows = np.array(row)
+        pha = np.array(pha)
+        x, y = self._grid.pixel_to_world(cols, rows)
+        args = x + dx, y + dy, 0.5 * self._grid.pitch, self._grid.hexagon_orientation()
+        collection = HexagonCollection(*args, **kwargs)
+        if values:
+            # Draw the pixel values
+            fmt = dict(ha="center", va="center", fontsize="small")
+            for x, y, value in zip(collection.x, collection.y, pha.flatten()):
+                if value > zero_sup_threshold:
+                    plt.text(x, y, f"{value}", color="black", **fmt)
+        plt.gca().add_collection(collection)
+        return collection
+
+    def draw_digi_event(self, event, zero_sup_threshold) -> HexagonCollection:
+        """Draw a digi event.
+
+        This is just dispatching the call to the proper method depending
+        on the event type.
+        """
+        if isinstance(event, DigiEventRectangular):
+            return self.draw_digi_event_rectangular(event, zero_sup_threshold=zero_sup_threshold)
+        if isinstance(event, DigiEventCircular):
+            return self.draw_digi_event_circular(event, zero_sup_threshold=zero_sup_threshold)
+        raise NotImplementedError(f"Cannot draw event of type {type(event)}.")
+
+    def draw_positions(self, mc_event: MonteCarloEvent, digi_event: DigiEventBase,
+                       readout: HexagonalReadoutBase, recon_defaults: object,
+                       zero_sup_threshold: int) -> None:
+        """Draw the Monte Carlo truth position and the reconstructed positions on top of the digi
+        event.
+        """
+        # Plot the Monte Carlo truth position.
+        plt.scatter(mc_event.absx, mc_event.absy, marker=".", s=100, label="Monte Carlo")
+        # Calculate the cluster from the digi event.
+        cluster = ClusteringNN(readout, zero_sup_threshold,
+                               num_neighbors=6).run(digi_event)
+        # Calculate and plot centroid position.
+        centroid_position = cluster.centroid()
+        plt.scatter(*centroid_position, marker="x", s=100, label="Centroid")
+        # Calculate and plot eta reconstructed position.
+        eta_recon_args = (recon_defaults.eta_2pix_rad, recon_defaults.eta_3pix_rad0,
+                          recon_defaults.eta_3pix_rad1, recon_defaults.eta_3pix_theta0)
+        try:
+            eta_position = cluster.eta(*eta_recon_args, pitch=readout.pitch)
+            # If cluster size is not 2 or 3, eta returns the centroid position, so we only
+            # plot it if it's different from the centroid.
+            plt.scatter(*eta_position, marker="+", s=100, label=r"$\eta$")
+        except RuntimeError:
+            pass
+        plt.legend()

src/hexsample/pipeline.py

@@ -62,7 +62,9 @@ def display(**kwargs) -> None:
     """Display events from a digi or recon file.
     """
     input_file_path = kwargs["input_file"]
-    return tasks.display(input_file_path)
+    zero_sup_threshold = kwargs.get("zero_sup_threshold", tasks.DisplayDefaults.zero_sup_threshold)
+    event_id = kwargs.get("event_id", tasks.DisplayDefaults.event_id)
+    return tasks.display(input_file_path, zero_sup_threshold, event_id)
 
 
 def quicklook(**kwargs) -> None:

src/hexsample/tasks.py

@@ -35,7 +35,6 @@
 from .clustering import ClusteringNN
 from .display import HexagonalGridDisplay
 from .fileio import (
-    DigiInputFileRectangular,
     ReconInputFile,
     ReconOutputFile,
     digi_input_file_class,
@@ -154,11 +153,13 @@ def simulate(
 
 @dataclass(frozen=True)
 class ReconstructionDefaults:
+
     """Default parameters for the reconstruction task.
 
     This is a small helper dataclass to help ensure consistency between the main task
     definition in this Python module and the command-line interface.
     """
+
     suffix: str = "recon"
     zero_sup_threshold: int = 0
     num_neighbors: int = 2
@@ -258,33 +259,67 @@ def reconstruct(
 
 
 class DisplayDefaults:
+
     """Default parameters for the display task.
 
     This is a small helper dataclass to help ensure consistency between the main task
     definition in this Python module and the command-line interface.
     """
 
+    zero_sup_threshold: int = 30
+    num_neighbors: int = 6
+    event_id: int = None
+
 
-def display(input_file_path: str) -> None:
+def display(
+        input_file_path: str,
+        zero_sup_threshold: int = DisplayDefaults.zero_sup_threshold,
+        event_id: int = DisplayDefaults.event_id
+        ) -> None:
     """Display events from a digi file.
 
     Arguments
     ---------
     file_path : str
         The path to the digi file.
+    zero_sup_threshold : int
+        The zero-suppression threshold to use when displaying the digi event.
+    event_id : int
+        The ID of the event to display. If None, display all events.
     """
     name, args = current_call()
     logger.info(f"Running {__name__}.{name} with arguments {args}...")
-    input_file = DigiInputFileRectangular(input_file_path)
+
+    # Note we cast the input file to string, in case it happens to be a pathlib.Path object.
+    input_file_path = str(input_file_path)
+    if not input_file_path.endswith(".h5"):
+        raise RuntimeError("Input file {input_file_path} does not look like a HDF5 file")
+
+    # It is necessary to extract the reaodut type because every readout type
+    # corresponds to a different DigiEvent type.
+    readout_mode = peek_readout_type(input_file_path)
+    # And we should get rid of all this crap when we store the readout type and all the
+    # relevant metadata in the hdf5 file in a sensible way.
+    file_type = digi_input_file_class(readout_mode)
+    input_file = file_type(input_file_path)
     header = input_file.header
     args = HexagonalLayout(header["layout"]), header["num_cols"], header["num_rows"],\
         header["pitch"], header["enc"], header["gain"]
-    readout = HexagonalReadoutRectangular(*args)
+    if readout_mode is HexagonalReadoutMode.RECTANGULAR:
+        readout = HexagonalReadoutRectangular(*args, padding=header["padding"])
+    elif readout_mode is HexagonalReadoutMode.CIRCULAR:
+        readout = HexagonalReadoutCircular(*args)
+    else:
+        raise RuntimeError(f"Unsupported readout mode: {readout_mode}")
     logger.info(f"Readout chip: {readout}")
     grid_display = HexagonalGridDisplay(readout)
-    for event in input_file:
-        print(event.ascii())
-        grid_display.draw_digi_event(event, zero_sup_threshold=0)
+    for i, event in enumerate(input_file):
+        if event_id is not None and i != event_id:
+            continue
+        recon_defaults = ReconstructionDefaults
+        mc_event = input_file.mc_event(i)
+        grid_display.draw_digi_event(event, zero_sup_threshold=zero_sup_threshold)
+        grid_display.draw_positions(mc_event, event, readout, recon_defaults, zero_sup_threshold)
         grid_display.show()
     input_file.close()