Resolution analysis

docs/release_notes.rst

@@ -4,6 +4,17 @@ Release notes
 =============
 
 * Added rough support for DigiEventCircular objects in the event display.
+* New `resolution.py` module incorporating some common functions that were
+  previously living in the script area.
+* New script to calculate the EEF for events reconstructed with the
+  different algorithms, disaggregated by cluster size.
+* `--max-neighbors` option added to the cli interface.
+* eta reconstructions modified to use the centroid for events with more than
+  three pixels.
+* Pull requests merged and issues closed:
+
+  - https://github.com/lucabaldini/hexsample/pull/94
+  - https://github.com/lucabaldini/hexsample/pull/92
 
 
 Version 0.13.2 (2026-02-04)

scripts/eta.py

@@ -24,6 +24,8 @@
 # Parser object.
 HXETA_ARGPARSER = argparse.ArgumentParser(description=__description__)
 HXETA_ARGPARSER.add_argument("input_file", type=str, help="path to the input file")
+HXETA_ARGPARSER.add_argument("zero_sup_threshold", type=int,
+                             help="zero suppression threshold in electrons")
 HXETA_ARGPARSER.add_argument("--save", action="store_true",
                             help="save the calibration plots to the results directory")
 
@@ -282,12 +284,13 @@ def hxeta(**kwargs) -> tuple[AbstractFitModel, AbstractFitModel, AbstractFitMode
     file_type = digi_input_file_class(readout_mode)
     input_file = file_type(input_file_path)
     header = input_file.header
+    zero_sup_threshold = kwargs["zero_sup_threshold"]
     args = HexagonalLayout(header["layout"]), header["num_cols"], header["num_rows"],\
-        header["pitch"], header["enc"], header["gain"], header["zero_sup_threshold"]
+        header["pitch"], header["enc"], header["gain"], zero_sup_threshold
     readout = HexagonalReadoutCircular(*args)
     nneighbors = 6
     logger.info(f"Readout chip: {readout}")
-    clustering = ClusteringNN(readout, header["zero_sup_threshold"], nneighbors)
+    clustering = ClusteringNN(readout, zero_sup_threshold, nneighbors)
     # Create all the lists we need to fill
     size, x0, y0, absx, absy, eta, versors = [[] for _ in range(7)]
     for i, event in tqdm(enumerate(input_file)):

scripts/resolution.py

@@ -0,0 +1,166 @@
+import argparse
+from pathlib import Path
+
+import numpy as np
+from aptapy.plotting import plt
+
+from hexsample.fileio import ReconInputFile
+from hexsample.hexagon import HexagonalLayout
+from hexsample.pipeline import reconstruct, simulate
+from hexsample.resolution import eef, eef_size_scan, resolution_spatial_dependence
+
+__description__ = ""
+
+# Parser object.
+HXETA_ARGPARSER = argparse.ArgumentParser(description=__description__)
+HXETA_ARGPARSER.add_argument("enc", type=int,
+                             help="equivalent noise charge in electrons")
+HXETA_ARGPARSER.add_argument("zero_sup_threshold", type=int,
+                             help="zero suppression threshold in electrons")
+HXETA_ARGPARSER.add_argument("--save", action="store_true",
+                             help="save the figures")
+
+RESOLUTION_DIR = Path.home() / "hexsampledata" / "resolution"
+if not RESOLUTION_DIR.exists():
+    RESOLUTION_DIR.mkdir(parents=True, exist_ok=True)
+
+FIGURES_DIR = Path.home() / "hexsample_figures" / "resolution"
+if not FIGURES_DIR.exists():
+    FIGURES_DIR.mkdir(parents=True, exist_ok=True)
+
+
+def resolution(**kwargs):
+    """Run the resolution analysis. This analysis consists of multiple steps to study different
+    aspects of the resolution.
+
+    First the Encircled Energy Function (EEF) is created for different cluster sizes and
+    reconstruction algorithms for a given ENC and zero suppression threshold.
+
+    Then the spatial dependence of the resolution is studied by plotting the HEW as a function of
+    the reconstructed distance from the true pixel center.
+
+    Finally, for the given ENC, the EEF is plotted for different zero suppression thresholds and
+    for both the centroid and eta algorithms.
+    """
+    enc = kwargs["enc"]
+    zero_sup_threshold = kwargs["zero_sup_threshold"]
+    file_prefix = f"simulation_resolution_{enc}enc_hexagonal"
+    simulation_path = RESOLUTION_DIR / f"{file_prefix}.h5"
+    # Run simulation if the output file does not already exist
+    if not simulation_path.exists():
+        simulate(
+                num_events=100000,
+                output_file=str(simulation_path),
+                beam="hexagonal",
+                enc=enc,
+                zero_sup_threshold=0,
+                readout_mode="circular",
+                pitch=0.005,
+                layout=HexagonalLayout.ODD_R,
+                num_cols=304,
+                num_rows=352,
+                gain=1.,
+        )
+    recon_kwargs = dict(input_file=str(simulation_path),
+                        zero_sup_threshold=zero_sup_threshold,
+                        max_neighbors=6)
+    # Reconstruct the simulated file with different algorithms, first with centroid
+    centroid_prefix = f"recon_zsuprec{zero_sup_threshold}_centroid"
+    centroid_path = RESOLUTION_DIR / f"{file_prefix}_{centroid_prefix}.h5"
+    if not centroid_path.exists():
+        reconstruct(suffix=centroid_prefix, pos_recon_algorithm="centroid", **recon_kwargs)
+    # Reconstruct with the best algorithm (eta for 2 and 3, centroid otherwise)
+    best_prefix = f"recon_zsuprec{zero_sup_threshold}_best"
+    best_path = RESOLUTION_DIR / f"{file_prefix}_{best_prefix}.h5"
+    if not best_path.exists():
+        reconstruct(suffix=best_prefix, pos_recon_algorithm="eta",**recon_kwargs)
+
+    # Open the recon files
+    centroid_recon_file = ReconInputFile(str(centroid_path))
+    best_recon_file = ReconInputFile(str(best_path))
+    # Define the plotting range
+    x = np.linspace(0, 0.6, 101)
+    # Plot the centroid eef for all cluster sizes
+    centroid_fig = plt.figure(f"centroid_eef_{enc}enc_{zero_sup_threshold}zsup")
+    eef_size_scan(x, centroid_recon_file)
+    # Plot the eta + centroid eef for all cluster sizes
+    best_fig = plt.figure(f"best_eef_{enc}enc_{zero_sup_threshold}zsup")
+    eef_size_scan(x, best_recon_file)
+
+    # Plot the spatial dependence of the resolution for both algorithms (all events)
+    sp_dep_fig = plt.figure("resolution_spatial_dependence")
+    plt.plot(*resolution_spatial_dependence(best_recon_file, max_neighbors=6),
+             ".k", label=r"$\eta$ + centroid")
+    plt.plot(*resolution_spatial_dependence(centroid_recon_file, max_neighbors=6),
+             "vk", label="centroid", markersize=4.)
+    plt.xlabel(r"$r_0 / p$")
+    plt.ylabel("Half Energy Width")
+    plt.legend()
+
+    # Close the recon files
+    centroid_recon_file.close()
+    best_recon_file.close()
+
+    # Study the resolution as a function of zero sup threshold for both algorithms
+    zero_sup_ratios = np.linspace(0, 3, 4)
+    recon_kwargs = dict(input_file=str(simulation_path),
+                        max_neighbors=6)
+    eef_zsup_centroid_fig = plt.figure(f"eef_vs_zsup_centroid_enc{enc}")
+    print("Reconstructing files for centroid algorithm...")
+    for zero_sup_ratio in zero_sup_ratios:
+        zsup = int(zero_sup_ratio * enc)
+        suffix = f"recon_zsuprec{zsup}_centroid"
+        file_path = RESOLUTION_DIR / f"{file_prefix}_{suffix}.h5"
+        if not file_path.exists():
+            reconstruct(suffix=suffix, pos_recon_algorithm="centroid", zero_sup_threshold=zsup,
+                        **recon_kwargs)
+        # Open file and plot EEF
+        recon_file = ReconInputFile(str(file_path))
+        plt.plot(x, eef(x, recon_file, max_neighbors=6), label=f"zsup/enc {zero_sup_ratio}")
+        recon_file.close()
+    plt.xlabel(xlabel = r"$r/p$")
+    plt.ylabel("Encircled Energy Fraction")
+    plt.xlim(x[0], x[-1])
+    plt.ylim(0, 1)
+    plt.legend()
+
+    eef_zsup_best_fig = plt.figure(f"eef_vs_zsup_best_enc{enc}")
+    print("Reconstructing files for eta algorithm...")
+    for zero_sup_ratio in zero_sup_ratios:
+        zsup = int(zero_sup_ratio * enc)
+        suffix = f"recon_zsuprec{zsup}_best"
+        file_path = RESOLUTION_DIR / f"{file_prefix}_{suffix}.h5"
+        if not file_path.exists():
+            reconstruct(suffix=suffix, pos_recon_algorithm="eta", zero_sup_threshold=zsup,
+                        **recon_kwargs)
+        # Open file and plot EEF
+        recon_file = ReconInputFile(str(file_path))
+        plt.plot(x, eef(x, recon_file, max_neighbors=6), label=f"zsup/enc {zero_sup_ratio}")
+        recon_file.close()
+
+    plt.xlabel(xlabel = r"$r/p$")
+    plt.ylabel("Encircled Energy Fraction")
+    plt.xlim(x[0], x[-1])
+    plt.ylim(0, 1)
+    plt.legend()
+
+    # Save figures, if requested
+    if kwargs["save"]:
+        fig_format = "png"
+        zsup_th = zero_sup_threshold
+        centroid_fig.savefig(FIGURES_DIR / f"centroid_eef_{enc}enc_{zsup_th}zsup.{fig_format}",
+                             format=fig_format)
+        best_fig.savefig(FIGURES_DIR / f"best_eef_{enc}enc_{zsup_th}zsup.{fig_format}",
+                         format=fig_format)
+        sp_dep_fig.savefig(FIGURES_DIR / f"res_spatial_depend_{enc}enc_{zsup_th}zsup.{fig_format}",
+                           format=fig_format)
+        eef_zsup_centroid_fig.savefig(FIGURES_DIR / f"eef_vs_zsup_centroid_{enc}enc.{fig_format}",
+                                      format=fig_format)
+        eef_zsup_best_fig.savefig(FIGURES_DIR / f"eef_vs_zsup_best_{enc}enc.{fig_format}",
+                                  format=fig_format)
+
+
+
+if __name__ == "__main__":
+    resolution(**vars(HXETA_ARGPARSER.parse_args()))
+    plt.show()

src/hexsample/cli.py

@@ -265,6 +265,8 @@ def add_recon_options(self, parser: argparse.ArgumentParser) -> None:
                            help="zero-suppression threshold in ADC counts")
         group.add_argument("--num_neighbors", type=int, default=2,
                            help="number of neighbors to be considered (0--6)")
+        group.add_argument("--max_neighbors", type=int, default=-1,
+                           help="maximum number of neighbors to be considered")
         group.add_argument("--pos_recon_algorithm", choices=["centroid", "eta", "dnn", "gnn"],
                            type=str, default="centroid", help="How to reconstruct position")
         group.add_argument("--eta_2pix_rad",

src/hexsample/clustering.py

@@ -112,7 +112,8 @@ def versors(self) -> Tuple[np.ndarray, np.ndarray]:
     def eta(self, eta_2pix_rad: float, eta_3pix_rad0: float, eta_3pix_rad1: float,
             eta_3pix_theta0: float, pitch: float) -> Tuple[float, float]:
         """Return the cluster reconstructed position using the eta function calibrated for 2
-        and 3 pixel clusters.
+        and 3 pixel clusters. If cluster size is not 2 or 3, reconstruct the position with the
+        centroid.
 
         Arguments
         ---------
@@ -127,18 +128,21 @@ def eta(self, eta_2pix_rad: float, eta_3pix_rad0: float, eta_3pix_rad1: float,
         pitch : float
             The pitch of the pixels.
         """
+        # Return the centroid position if it's not possible to use the eta function
+        if self.size() not in (2, 3):
+            return self.centroid()
+        # Calculate versors and eta.
         u, v = self.versors()
         _eta = self.calculate_eta()
 
         if self.size() == 2:
             # For 2-pixel events we estimate the position along the line that connects the
-            # two pixels using the eta function calibration.
+            # two pixels using the probit function.
             r = Probit().evaluate(_eta[0], 0.5, eta_2pix_rad)
             x_recon = self.x[0] + r * pitch * u[0]
             y_recon = self.y[0] + r * pitch * u[1]
         elif self.size() == 3:
-            # For 3-pixel events we estimate both r and theta using the eta function
-            # calibrations.
+            # For 3-pixel events we estimate both r and theta using the probit function.
             eta_sum = _eta[0] + _eta[1]
             eta_diff = (_eta[0] - _eta[1]) / eta_sum
             r = Probit().evaluate(eta_sum, eta_3pix_rad0, eta_3pix_rad1)
@@ -147,8 +151,10 @@ def eta(self, eta_2pix_rad: float, eta_3pix_rad0: float, eta_3pix_rad1: float,
             x_recon = self.x[0] + r * pitch * (np.cos(theta) * u[0] + np.sin(theta) * v[0])
             y_recon = self.y[0] + r * pitch * (np.cos(theta) * u[1] + np.sin(theta) * v[1])
         else:
-            raise RuntimeError("Cluster must contain 2 or 3 pixels to reconstruct position using "
-                               "eta function")
+            # This condition should never be reached because of the check at the beginning of the
+            # method, but it's here for safety.
+            raise RuntimeError("Cluster must contain 2 or 3 pixels to reconstruct position with" \
+                               " eta function")
         return x_recon, y_recon
 
 

src/hexsample/pipeline.py

@@ -48,12 +48,13 @@ def reconstruct(**kwargs) -> str:
     suffix = kwargs.get("suffix", defaults.suffix)
     zero_sup_threshold = kwargs.get("zero_sup_threshold", defaults.zero_sup_threshold)
     num_neighbors = kwargs.get("num_neighbors", defaults.num_neighbors)
+    max_neighbors = kwargs.get("max_neighbors", defaults.max_neighbors)
     pos_recon_algorithm = kwargs.get("pos_recon_algorithm", defaults.pos_recon_algorithm)
     eta_2pix_rad = kwargs.get("eta_2pix_rad", defaults.eta_2pix_rad)
     eta_3pix_rad0 = kwargs.get("eta_3pix_rad0", defaults.eta_3pix_rad0)
     eta_3pix_rad1 = kwargs.get("eta_3pix_rad1", defaults.eta_3pix_rad1)
     eta_3pix_theta0 = kwargs.get("eta_3pix_theta0", defaults.eta_3pix_theta0)
-    args = input_file_path, suffix, zero_sup_threshold, num_neighbors, \
+    args = input_file_path, suffix, zero_sup_threshold, num_neighbors, max_neighbors, \
            pos_recon_algorithm, eta_2pix_rad, eta_3pix_rad0, eta_3pix_rad1, eta_3pix_theta0
     return tasks.reconstruct(*args, kwargs)