diff --git a/docs/release_notes.rst b/docs/release_notes.rst index fa73677..4423fab 100644 --- a/docs/release_notes.rst +++ b/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) diff --git a/scripts/eta.py b/scripts/eta.py index 74b28c7..e013c42 100644 --- a/scripts/eta.py +++ b/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)): diff --git a/scripts/resolution.py b/scripts/resolution.py new file mode 100644 index 0000000..4c84804 --- /dev/null +++ b/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() diff --git a/src/hexsample/cli.py b/src/hexsample/cli.py index 947ec5b..05fd707 100644 --- a/src/hexsample/cli.py +++ b/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", diff --git a/src/hexsample/clustering.py b/src/hexsample/clustering.py index 7947680..2474cab 100644 --- a/src/hexsample/clustering.py +++ b/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 diff --git a/src/hexsample/pipeline.py b/src/hexsample/pipeline.py index c82448d..ac96a9e 100644 --- a/src/hexsample/pipeline.py +++ b/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) diff --git a/src/hexsample/resolution.py b/src/hexsample/resolution.py new file mode 100644 index 0000000..b5d1561 --- /dev/null +++ b/src/hexsample/resolution.py @@ -0,0 +1,252 @@ +# Copyright (C) 2023--2025 the hexsample team. +# +# For the license terms see the file LICENSE, distributed along with this +# software. +# +# This program is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 2 of the License, or (at your +# option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# General Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +"""Resolution analysis facilities. +""" + +from typing import Tuple + +import numpy as np +from aptapy.hist import Histogram1d, Histogram2d +from aptapy.plotting import plt + +from .fileio import ReconInputFile +from .hexagon import HexagonalGrid + + +def dist_residual(input_file: ReconInputFile) -> np.ndarray: + """Calculate the distance residuals between reconstructed and Monte Carlo true positions. + + Arguments + --------- + input_file : ReconInputFile + The input file to analyze. + + Returns + ------- + dr : np.ndarray + The distance residuals. + """ + # Access the Monte Carlo true positions + x_mc = input_file.mc_column("absx") + y_mc = input_file.mc_column("absy") + # Access the reconstructed positions + x = input_file.column("posx") + y = input_file.column("posy") + # Calculate the distance residuals + dr = np.sqrt((x - x_mc) ** 2 + (y - y_mc) ** 2) + return dr + + +def dist_from_pixel_center(input_file: ReconInputFile) -> np.ndarray: + """Calculate the reconstructed distance from the pixel center. + + Arguments + --------- + input_file : ReconInputFile + The input file to analyze. + + Returns + ------- + dr0 : np.ndarray + The reconstructed distance from the pixel center. + """ + # Create the hexagonal grid to get pixel centers + layout = input_file.digi_header["layout"] + num_cols = input_file.digi_header["num_cols"] + num_rows = input_file.digi_header["num_rows"] + grid = HexagonalGrid(layout=layout, num_cols=num_cols, num_rows=num_rows) + # Access Monte Carlo true positions + x_mc = input_file.mc_column("absx") + y_mc = input_file.mc_column("absy") + # Calculate Monte Carlo true pixel centers + x0, y0 = grid.pixel_to_world(*grid.world_to_pixel(x_mc, y_mc)) + # Access the reconstructed positions + x = input_file.column("posx") + y = input_file.column("posy") + # Calculate the reconstructed distance from the pixel center + dr0 = np.sqrt((x - x0) ** 2 + (y - y0) ** 2) + return dr0 + + +def hist_distance_residuals(input_file: ReconInputFile, num_neighbors: int = 0, + max_neighbors: int = -1) -> Histogram1d: + """Create the histogram of distance residuals for the given number of neighbors. + + Arguments + --------- + input_file : ReconInputFile + The input file to analyze. + num_neighbors : int, optional + The number of neighbors to be considered. Default is 0. + max_neighbors : int, optional + The maximum number of neighbors to be considered. If max_neighbors is specified, it has + priority over num_neighbors. Default is -1 (not used). + + Returns + ------- + hist : Histogram1d + The histogram of pitch normalized distance residuals. + """ + # Select the cluster sizes we want and create the mask + size = input_file.column("cluster_size") + mask = size <= max_neighbors + 1 if max_neighbors >= 0 else size == num_neighbors + 1 + pitch = input_file.digi_header["pitch"] + # Calculate the distance residuals normalized to pitch + dr = dist_residual(input_file) / pitch + # Create the histogram to calculate the EEF. The binning is taken in a way that + # spans all the pitch. + dr_binning = np.linspace(0., 1., 101) + hist = Histogram1d(dr_binning) + hist.fill(dr[mask]) + return hist + + +def eef(x: np.ndarray, input_file: ReconInputFile, num_neighbors: int = 0, + max_neighbors: int = -1) -> np.ndarray: + """Calculate the Encircled Energy Function (EEF) for a given cluster size. + + Arguments + --------- + x : np.ndarray + The normalized distance values where the EEF will be evaluated. + input_file : ReconInputFile + The input file to analyze. + num_neighbors : int, optional + The number of neighbors to be considered. Default is 0. + max_neighbors : int, optional + The maximum number of neighbors to be considered. If max_neighbors is specified, it has + priority over num_neighbors. Default is -1 (not used). + + Returns + ------- + eef : np.ndarray + The Encircled Energy Function evaluated at the given normalized distance values. + """ + hist = hist_distance_residuals(input_file, num_neighbors, max_neighbors) + return hist.cdf(x) + + +def hew(input_file: ReconInputFile, num_neighbors: int = 0, + max_neighbors: int = -1) -> float: + """Calculate the Half Energy Width (HEW) for a given cluster size. + + Arguments + --------- + input_file : ReconInputFile + The input file to analyze. + num_neighbors : int, optional + The number of neighbors to be considered. Default is 0. + max_neighbors : int, optional + The maximum number of neighbors to be considered. If max_neighbors is specified, it has + priority over num_neighbors. Default is -1 (not used). + + Returns + ------- + hew : float + The Half Energy Width (HEW) in pitch normalized units. + """ + hist = hist_distance_residuals(input_file, num_neighbors, max_neighbors) + return hist.ppf(0.5) + + +def eef_size_scan(x: np.ndarray, input_file: ReconInputFile) -> None: + """Plot the Encircled Energy Function for one, two, three pixel events and for all the + reconstructed events on the same figure. + + Arguments + --------- + x : np.ndarray + The pitch normalized distance values where the EEF will be evaluated. + input_file : ReconInputFile + The input file to analyze. + """ + xlabel = r"$r/p$" + ylabel = "Encircled Energy Fraction" + color = "black" + # Plot the EEFs for the different cluster sizes + plt.plot(x, eef(x, input_file, 0), + label=f"1 pix (HEW={hew(input_file, 0):.2f})", + linestyle=":", color=color) + plt.plot(x, eef(x, input_file, 1), + label=f"2 pix (HEW={hew(input_file, 1):.2f})", + linestyle="-", color=color) + plt.plot(x, eef(x, input_file, 2), + label=f"3 pix (HEW={hew(input_file, 2):.2f})", + linestyle="--", color=color) + plt.plot(x, eef(x, input_file, max_neighbors=6), + label=f"All events (HEW={hew(input_file, max_neighbors=6):.2f})", + linestyle="-.", color=color) + # Plot the line corresponding to 50% of the energy + plt.hlines(0.5, x[0], x[-1], colors="0.4", linestyles="--") + # Finalize the plot + plt.xlim(x[0], x[-1]) + plt.ylim(0, 1) + plt.xlabel(xlabel) + plt.ylabel(ylabel) + plt.legend() + + +def resolution_spatial_dependence(input_file: ReconInputFile, num_neighbors: int = 1, + max_neighbors: int = -1) -> Tuple[np.ndarray, np.ndarray]: + """Calculate the spatial dependence of the resolution. This is estimated by calculating the + Half Energy Width (HEW) as a function of the reconstructed distance from the true pixel center. + + Arguments + --------- + input_file : ReconInputFile + The input file to analyze. + num_neighbors : int, optional + The number of neighbors to be considered. Default is 1. + max_neighbors : int, optional + The maximum number of neighbors to be considered. If max_neighbors is specified, it has + priority over num_neighbors. Default is -1 (not used). + + Returns + ------- + bin_centers : np.ndarray + The bin centers of the reconstructed distance from the pixel center. + hews : np.ndarray + The HEW calculated for each bin of the reconstructed distance from the pixel center. + """ + # Select the cluster sizes we want and create the mask + size = input_file.column("cluster_size") + mask = size <= max_neighbors + 1 if max_neighbors >= 0 else size == num_neighbors + 1 + # Calculate the pitch normalized distance from pixel center and distance residuals + pitch = input_file.digi_header["pitch"] + dr0 = dist_from_pixel_center(input_file)[mask] / pitch + dr = dist_residual(input_file)[mask] / pitch + # Create the 2D histogram + xedges = np.linspace(0., 1., 101) + if dr.size == 0: + # If there are no events, raise an error + raise ValueError("No events found for the given cluster size.") + yedges = np.linspace(min(dr), max(dr), 101) + hist = Histogram2d(xedges, yedges) + hist.fill(dr0, dr) + # Calculate the hew for each slice of the histogram + bin_centers = hist.bin_centers() + hews = np.zeros(len(bin_centers)) + for i in range(len(bin_centers)): + xslice = hist.slice1d(i) + if xslice.content.sum() == 0: + hews[i] = np.nan + else: + hews[i] = xslice.ppf(0.5) + return bin_centers, hews diff --git a/src/hexsample/tasks.py b/src/hexsample/tasks.py index 578ef88..77e9a06 100644 --- a/src/hexsample/tasks.py +++ b/src/hexsample/tasks.py @@ -163,6 +163,7 @@ class ReconstructionDefaults: suffix: str = "recon" zero_sup_threshold: int = 0 num_neighbors: int = 2 + max_neighbors: int = -1 pos_recon_algorithm: str = "centroid" eta_2pix_rad: float = 0.127 eta_3pix_rad0: float = 0.513 @@ -175,6 +176,7 @@ def reconstruct( suffix: str = ReconstructionDefaults.suffix, zero_sup_threshold: int = ReconstructionDefaults.zero_sup_threshold, num_neighbors: int = ReconstructionDefaults.num_neighbors, + max_neighbors: int = ReconstructionDefaults.max_neighbors, pos_recon_algorithm: str = ReconstructionDefaults.pos_recon_algorithm, eta_2pix_rad: float = ReconstructionDefaults.eta_2pix_rad, eta_3pix_rad0: float = ReconstructionDefaults.eta_3pix_rad0, @@ -203,6 +205,10 @@ def reconstruct( num_neighbors : int The number of neighbor pixels to be used for the clustering. + + max_neighbors : int + The maximum number of neighbor pixels to be used for the clustering. If max_neighbors is + specified (i.e. different from -1), it has priority over num_neighbors. pos_recon_algorithm : str The position reconstruction algorithm to use. @@ -235,16 +241,23 @@ def reconstruct( raise RuntimeError(f"Unsupported readout mode: {readout_mode}") logger.info(f"Readout chip: {readout}") + # Define the effective number of neighbors to be used for the clustering. If max_neighbors is + # specified (i.e. different from -1), it has priority over num_neighbors. It is necessary to + # define it here because rectangular readout doesn't have a fixed number of neighbors, contrary + # to the circular. + effective_neighbors = max_neighbors if max_neighbors >= 0 else num_neighbors # Run the actual reconstruction. - clustering = ClusteringNN(readout, zero_sup_threshold, num_neighbors) + clustering = ClusteringNN(readout, zero_sup_threshold, effective_neighbors) output_file_path = input_file_path.replace(".h5", f"_{suffix}.h5") output_file = ReconOutputFile(output_file_path) if header_kwargs is not None: output_file.update_header(**header_kwargs) output_file.update_digi_header(**input_file.header) + # Create a list of acceptable cluster sizes. + size = list(range(1, max_neighbors + 2)) if max_neighbors >= 0 else [num_neighbors + 1] for i, event in tqdm(enumerate(input_file)): cluster = clustering.run(event) - if num_neighbors == 0 or cluster.size() == num_neighbors: + if cluster.size() in size: # Need to pass the recon method and other stuff as argument to ReconEvent args = event.trigger_id, event.timestamp(), event.livetime, cluster recon_event = ReconEvent(*args, pos_recon_algorithm, readout.pitch, @@ -378,6 +391,5 @@ def quicklook(input_file_path: str) -> None: binning = np.linspace((y-y_mc).min(), (y-y_mc).max(), 100) histy = Histogram1d(binning, xlabel=r"$y - y_{MC}$ [cm]").fill(y-y_mc) histy.plot() - input_file.close() plt.show()