From 783e620a840416449a9845bff01ac96e71a37a01 Mon Sep 17 00:00:00 2001 From: Augusto Cattafesta Date: Mon, 19 Jan 2026 13:55:17 +0100 Subject: [PATCH 01/18] Script to analyze noise added. --- src/hexsample/pipeline.py | 7 +++++ src/hexsample/tasks.py | 55 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 62 insertions(+) diff --git a/src/hexsample/pipeline.py b/src/hexsample/pipeline.py index d7c094fb..d42e7f42 100644 --- a/src/hexsample/pipeline.py +++ b/src/hexsample/pipeline.py @@ -67,3 +67,10 @@ def quicklook(**kwargs) -> None: """ input_file_path = kwargs["input_file"] return tasks.quicklook(input_file_path) + + +def resolution(**kwargs) -> None: + """Analyze the resolution from a recon file. + """ + input_file_path = kwargs["input_file"] + return tasks.resolution(input_file_path) diff --git a/src/hexsample/tasks.py b/src/hexsample/tasks.py index 195a0489..91e0534d 100644 --- a/src/hexsample/tasks.py +++ b/src/hexsample/tasks.py @@ -339,3 +339,58 @@ def quicklook(input_file_path: str) -> None: input_file.close() plt.show() + + +class ResolutionDefaults: + """Default parameters for the resolution 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. + """ + + +def resolution(input_file_path: str) -> None: + """Calculate the distribution of the distance between reconstructed and true position + and calculate its mean, standard deviation and FWHM. + + Arguments + --------- + input_file_path : str + The path to the input recon file. + + Returns + ------- + mean : float + The mean of the distance distribution. + stddev : float + The standard deviation of the distance distribution. + fwhm : float + The FWHM of the distance distribution. + """ + name, args = current_call() + logger.info(f"Running {__name__}.{name} with arguments {args}...") + input_file = ReconInputFile(input_file_path) + header = input_file.digi_header + # Get the reconstructed and true positions + x_mc = input_file.mc_column("absx") + y_mc = input_file.mc_column("absy") + x = input_file.column("posx") + y = input_file.column("posy") + # Calculate the distance between reconstructed and true position + dr = np.sqrt((x - x_mc)**2 + (y - y_mc)**2) / header["pitch"] + # Create the histogram of the distance + edges = np.linspace(0., 1., 101) + hist = Histogram1d(edges, xlabel=r"$|\vec{r} - \vec{r}_{mc}|$ / pitch") + hist.fill(dr) + mean, stddev = hist.binned_statistics() + # Calculate the FWHM if possible + try: + fwhm = hist.fwhm() + except ValueError: + fwhm = np.nan + plt.figure("Distance resolution") + hist.plot() + + input_file.close() + + return mean, stddev, fwhm From a05e3afd88ebdc5c025c0a1f2274c12f6967800d Mon Sep 17 00:00:00 2001 From: Augusto Cattafesta Date: Mon, 19 Jan 2026 13:55:40 +0100 Subject: [PATCH 02/18] Script to analyze noise added. --- scripts/analyze_noise.py | 101 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 101 insertions(+) create mode 100644 scripts/analyze_noise.py diff --git a/scripts/analyze_noise.py b/scripts/analyze_noise.py new file mode 100644 index 00000000..d31c110e --- /dev/null +++ b/scripts/analyze_noise.py @@ -0,0 +1,101 @@ +# 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. + +"""Analyze the noise to find the best the zero suppression threshold. +""" + +import argparse +import pathlib + +import numpy as np +from aptapy.plotting import plt + +from hexsample.fileio import digi_input_file_class, peek_readout_type +from hexsample.pipeline import reconstruct, resolution +from hexsample.readout import HexagonalReadoutMode + +ARGPARSER = argparse.ArgumentParser() +ARGPARSER.add_argument("input_file", type=str, help="path to the input file") +ARGPARSER.add_argument("--num_neighbors", type=int, default=2, + help="number of neighbors for the clustering (default: 2)") +ARGPARSER.add_argument("--pos_recon_algorithm", type=str, default="eta", + help="position reconstruction algorithm (default: eta)") + + +def analyze_noise(**kwargs): + """ Analyze the effect of the noise and zero suppression threshold on the reconstucted position + resolution, by calculating the bias and the standard deviation of the position residuals + distribution for different values of the threshold. + """ + # Open the file to get the header information + input_file_path = str(kwargs["input_file"]) + if not input_file_path.endswith(".h5"): + raise RuntimeError(f"Input file {input_file_path} does not look like a HDF5 file") + readout_mode = peek_readout_type(input_file_path) + if readout_mode is not HexagonalReadoutMode.CIRCULAR: + raise RuntimeError("Only CIRCULAR readout is supported.") + file_type = digi_input_file_class(readout_mode) + input_file = file_type(input_file_path) + header = input_file.header + input_file.close() + + enc = header["enc"] + # Scan different zero-suppression thresholds + n_files = 10 + zero_sup_threshold_grid = np.linspace(0, 5 * enc, n_files, dtype=int) + results = np.zeros((n_files, 3)) + default_folder = pathlib.Path.home() / "hexsampledata" + for i, zero_sup_threshold in enumerate(zero_sup_threshold_grid): + # Define the arguments for the reconstruction + recon_kwargs = dict( + input_file=input_file_path, + suffix=(f"recon_enc{enc}_thr{zero_sup_threshold}_" + f"{kwargs['pos_recon_algorithm']}_nn{kwargs['num_neighbors']}"), + zero_sup_threshold=zero_sup_threshold, + num_neighbors=kwargs["num_neighbors"], + pos_recon_algorithm=kwargs["pos_recon_algorithm"], + ) + # Check if reconstructed files already exist + file_name = default_folder / ( + f"{pathlib.Path(recon_kwargs['input_file']).stem}_" + f"{recon_kwargs['suffix']}.h5") + if not file_name.with_suffix(".h5").exists(): + output_file_path = reconstruct(**recon_kwargs) + else: + output_file_path = str(file_name.with_suffix(".h5")) + # Analyze the resolution of the reconstructed file + res_kwargs = dict( + input_file=output_file_path, + ) + results[i] = resolution(**res_kwargs) + + plt.figure("Position resolution vs zero-suppression threshold") + plt.plot(zero_sup_threshold_grid / enc, results[:, 0], ".k", label="Mean") + plt.plot(zero_sup_threshold_grid / enc, results[:, 1], ".r", label="Stddev") + # plt.plot(zero_sup_threshold_grid / enc, results[:, 2], ".b", label="FWHM") + plt.xlabel("Zero-suppression threshold / enc") + plt.ylabel("Position resolution [pitch]") + plt.legend() + + plt.show() + + +if __name__ == "__main__": + analyze_noise(**vars(ARGPARSER.parse_args())) + From d4bbf93160fb0d044e8d8993460eaf71ead8dc71 Mon Sep 17 00:00:00 2001 From: Augusto Cattafesta Date: Mon, 2 Feb 2026 11:12:13 +0100 Subject: [PATCH 03/18] Minor --- src/hexsample/cli.py | 2 ++ src/hexsample/tasks.py | 8 +++++++- 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/src/hexsample/cli.py b/src/hexsample/cli.py index 88cfd988..4a2165f8 100644 --- a/src/hexsample/cli.py +++ b/src/hexsample/cli.py @@ -259,6 +259,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=0, + 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/tasks.py b/src/hexsample/tasks.py index c8bfab29..46a92511 100644 --- a/src/hexsample/tasks.py +++ b/src/hexsample/tasks.py @@ -162,6 +162,7 @@ class ReconstructionDefaults: suffix: str = "recon" zero_sup_threshold: int = 0 num_neighbors: int = 2 + max_neighbors: int = 0 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, @@ -245,7 +247,11 @@ def reconstruct( output_file.update_digi_header(**input_file.header) for i, event in tqdm(enumerate(input_file)): cluster = clustering.run(event) - if num_neighbors == 0 or cluster.size() == num_neighbors: + # Create a list of acceptable cluster sizes. If max_neighbors is 0, take num_neighbors + # only. Otherwise take all sizes from 1 to max_neighbors. + # I'd like to use cluster size instead of number of neighbors + neighbors = list(range(1, max_neighbors + 2)) if max_neighbors > 0 else [num_neighbors + 1] + if cluster.size() in neighbors: # 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, From c11f6c4b2cb81dbe8697a8e334dce1b15dd1c961 Mon Sep 17 00:00:00 2001 From: Augusto Cattafesta Date: Mon, 2 Feb 2026 14:13:40 +0100 Subject: [PATCH 04/18] Max neighbors argument added. --- src/hexsample/cli.py | 2 +- src/hexsample/clustering.py | 50 ++++++++++++++++++------------------- src/hexsample/pipeline.py | 3 ++- src/hexsample/tasks.py | 24 +++++++++++++----- 4 files changed, 46 insertions(+), 33 deletions(-) diff --git a/src/hexsample/cli.py b/src/hexsample/cli.py index 4a2165f8..17a0d6e2 100644 --- a/src/hexsample/cli.py +++ b/src/hexsample/cli.py @@ -259,7 +259,7 @@ 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=0, + 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") diff --git a/src/hexsample/clustering.py b/src/hexsample/clustering.py index ecfaf63c..78b69ce3 100644 --- a/src/hexsample/clustering.py +++ b/src/hexsample/clustering.py @@ -94,7 +94,7 @@ def n_versor(self) -> np.ndarray: def eta(self, eta_2pix_rad: float, eta_3pix_rad0: float, eta_3pix_rad1: float, eta_3pix_theta0: float, eta_3pix_theta1: 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 1, return centroid. Arguments --------- @@ -112,31 +112,31 @@ def eta(self, eta_2pix_rad: float, eta_3pix_rad0: float, eta_3pix_rad1: float, pitch : float The pitch of the pixels. """ - n = self.n_versor() - _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. - dr = Probit().evaluate(_eta[0], 0.5, eta_2pix_rad) * pitch - x_recon = self.x[0] + dr * n[0] - y_recon = self.y[0] + dr * n[1] - elif self.size() == 3: - # For 3-pixel events we estimate both dr and theta using the eta function - # calibrations. - eta_sum = _eta[0] + _eta[1] - eta_diff = (_eta[0] - _eta[1]) / eta_sum - dr = Probit().evaluate(eta_sum, eta_3pix_rad0, eta_3pix_rad1) * pitch - theta = Exponential().evaluate(eta_diff, eta_3pix_theta0, eta_3pix_theta1) - # We need to determine the sign of theta depending on the cluster orientation. - theta = theta * np.sign(np.cross(n, np.array([self.x[1] - self.x[0], - self.y[1] - self.y[0]]))) - x_recon = self.x[0] + dr * (np.cos(theta) * n[0] - np.sin(theta) * n[1]) - y_recon = self.y[0] + dr * (np.sin(theta) * n[0] + np.cos(theta) * n[1]) + if self.size() not in [2, 3]: + return self.centroid() else: - raise RuntimeError("Cluster must contain 2 or 3 pixels to reconstruct position using " - "eta function") - return x_recon, y_recon + n = self.n_versor() + _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. + dr = Probit().evaluate(_eta[0], 0.5, eta_2pix_rad) * pitch + x_recon = self.x[0] + dr * n[0] + y_recon = self.y[0] + dr * n[1] + elif self.size() == 3: + # For 3-pixel events we estimate both dr and theta using the eta function + # calibrations. + eta_sum = _eta[0] + _eta[1] + eta_diff = (_eta[0] - _eta[1]) / eta_sum + dr = Probit().evaluate(eta_sum, eta_3pix_rad0, eta_3pix_rad1) * pitch + theta = Exponential().evaluate(eta_diff, eta_3pix_theta0, eta_3pix_theta1) + # We need to determine the sign of theta depending on the cluster orientation. + theta = theta * np.sign(np.cross(n, np.array([self.x[1] - self.x[0], + self.y[1] - self.y[0]]))) + x_recon = self.x[0] + dr * (np.cos(theta) * n[0] - np.sin(theta) * n[1]) + y_recon = self.y[0] + dr * (np.sin(theta) * n[0] + np.cos(theta) * n[1]) + return x_recon, y_recon @dataclass diff --git a/src/hexsample/pipeline.py b/src/hexsample/pipeline.py index 5bd8b633..48b007b3 100644 --- a/src/hexsample/pipeline.py +++ b/src/hexsample/pipeline.py @@ -48,13 +48,14 @@ 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) eta_3pix_theta1 = kwargs.get("eta_3pix_theta1", defaults.eta_3pix_theta1) - 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, \ eta_3pix_theta1 return tasks.reconstruct(*args, kwargs) diff --git a/src/hexsample/tasks.py b/src/hexsample/tasks.py index 46a92511..42b21b6e 100644 --- a/src/hexsample/tasks.py +++ b/src/hexsample/tasks.py @@ -162,7 +162,7 @@ class ReconstructionDefaults: suffix: str = "recon" zero_sup_threshold: int = 0 num_neighbors: int = 2 - max_neighbors: int = 0 + max_neighbors: int = -1 pos_recon_algorithm: str = "centroid" eta_2pix_rad: float = 0.127 eta_3pix_rad0: float = 0.513 @@ -245,13 +245,13 @@ def reconstruct( 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. If max_neighbors is -1, take num_neighbors + # only. Otherwise take all sizes from 1 to max_neighbors + 1. When specified, max_neighbors + # gets priority over num_neighbors. + 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) - # Create a list of acceptable cluster sizes. If max_neighbors is 0, take num_neighbors - # only. Otherwise take all sizes from 1 to max_neighbors. - # I'd like to use cluster size instead of number of neighbors - neighbors = list(range(1, max_neighbors + 2)) if max_neighbors > 0 else [num_neighbors + 1] - if cluster.size() in 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, @@ -328,6 +328,13 @@ def quicklook(input_file_path: str) -> None: plt.xlabel("Energy [eV]") plt.legend() + size_hist = create_histogram(input_file, "cluster_size", mc=False) + plt.figure("Cluster size distribution") + size_hist.plot() + plt.xlabel("Cluster size") + plt.ylabel("Counts") + + # Plotting the reconstructed x and y position and the true position. plt.figure("Reconstructed photons position") binning = np.linspace(-5. * 0.2, 5. * 0.2, 100) @@ -351,6 +358,11 @@ 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() + plt.figure("distance resolution") + dr = np.sqrt((x - x_mc)**2 + (y - y_mc)**2) + binning = np.linspace(dr.min(), dr.max(), 100) + histdr = Histogram1d(binning, xlabel=r"$\sqrt{(x - x_{MC})^2 + (y - y_{MC})^2}$ [cm]").fill(dr) + histdr.plot() input_file.close() plt.show() From 4c192d4d4810bf8745022f22719023c2316b55f7 Mon Sep 17 00:00:00 2001 From: Augusto Cattafesta Date: Mon, 2 Feb 2026 16:29:35 +0100 Subject: [PATCH 05/18] Noise study. --- scripts/eta.py | 17 +++++++++++++++-- src/hexsample/clustering.py | 3 ++- src/hexsample/tasks.py | 21 ++++++++++++++++----- 3 files changed, 33 insertions(+), 8 deletions(-) diff --git a/scripts/eta.py b/scripts/eta.py index 4a082834..3d22ae5e 100644 --- a/scripts/eta.py +++ b/scripts/eta.py @@ -2,7 +2,7 @@ """ import numpy as np -from aptapy.hist import Histogram2d +from aptapy.hist import Histogram1d, Histogram2d from aptapy.modeling import AbstractFitModel from aptapy.models import Exponential, Probit from aptapy.plotting import plt @@ -146,6 +146,20 @@ def calibrate_2pix(eta: np.ndarray, photon_pos: np.ndarray, versor: np.ndarray plt.ylabel("dr / p") plt.legend() + eta_min = min(eta)[0] + print(eta_min) + eta_hist_binning = np.linspace(eta_min, 0.5, 101) + eta_hist = Histogram1d(eta_hist_binning, xlabel="eta", ylabel="Counts") + eta_hist.fill(eta) + plt.figure("eta_2pix_distribution") + eta_hist.plot() + + from scipy.special import ndtri + + sigma_x = model.sigma.value * np.sqrt(2 / np.pi) * np.exp(0.5 * abs(ndtri(eta_binning))**2) + plt.figure("dr_vs_eta_2pix_derivative") + plt.plot(eta_binning, sigma_x, label="d(dr/p)/d(eta)") + return model @@ -257,7 +271,6 @@ 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 - header["zero_sup_threshold"] = 0 args = HexagonalLayout(header["layout"]), header["num_cols"], header["num_rows"],\ header["pitch"], header["enc"], header["gain"], header["zero_sup_threshold"] readout = HexagonalReadoutCircular(*args) diff --git a/src/hexsample/clustering.py b/src/hexsample/clustering.py index 78b69ce3..6aca6156 100644 --- a/src/hexsample/clustering.py +++ b/src/hexsample/clustering.py @@ -94,7 +94,8 @@ def n_versor(self) -> np.ndarray: def eta(self, eta_2pix_rad: float, eta_3pix_rad0: float, eta_3pix_rad1: float, eta_3pix_theta0: float, eta_3pix_theta1: float, pitch: float) -> Tuple[float, float]: """Return the cluster reconstructed position using the eta function calibrated for 2 - and 3 pixel clusters. If cluster size is 1, return centroid. + and 3 pixel clusters. If cluster size is not 2 or 3, reconstruct the position with the + centroid. Arguments --------- diff --git a/src/hexsample/tasks.py b/src/hexsample/tasks.py index 42b21b6e..98ebd4d5 100644 --- a/src/hexsample/tasks.py +++ b/src/hexsample/tasks.py @@ -328,11 +328,11 @@ def quicklook(input_file_path: str) -> None: plt.xlabel("Energy [eV]") plt.legend() - size_hist = create_histogram(input_file, "cluster_size", mc=False) - plt.figure("Cluster size distribution") - size_hist.plot() - plt.xlabel("Cluster size") - plt.ylabel("Counts") + # size_hist = create_histogram(input_file, "cluster_size", mc=False) + # plt.figure("Cluster size distribution") + # size_hist.plot() + # plt.xlabel("Cluster size") + # plt.ylabel("Counts") # Plotting the reconstructed x and y position and the true position. @@ -364,6 +364,17 @@ def quicklook(input_file_path: str) -> None: histdr = Histogram1d(binning, xlabel=r"$\sqrt{(x - x_{MC})^2 + (y - y_{MC})^2}$ [cm]").fill(dr) histdr.plot() + from .hexagon import HexagonalGrid + grid = HexagonalGrid() + x0, y0 = grid.pixel_to_world(*grid.world_to_pixel(x_mc, y_mc)) + dr_abs = np.sqrt((x - x0)**2 + (y - y0)**2) / grid.pitch + bins = np.linspace(0, 1, 100) + hist = Histogram2d(bins, bins) + # I need the recon distance from the hit pixel center + hist.fill(dr_abs, dr / grid.pitch) + hist_mean, hist_sigma = hist.project_statistics() + plt.figure("Reconstructed vs true distance from pixel center") + plt.plot(hist_mean.bin_centers(), hist_mean.content.flatten(), 'o', label='Mean') input_file.close() plt.show() From 7cc76b81a44e9d1f3e83a7bd8b9adc0db09342d8 Mon Sep 17 00:00:00 2001 From: Augusto Cattafesta Date: Tue, 3 Feb 2026 16:34:46 +0100 Subject: [PATCH 06/18] Script to study the resolution. --- scripts/eta.py | 9 +- scripts/resolution.py | 357 +++++++++++++++++++++++++++++++++++++++++ src/hexsample/tasks.py | 20 ++- 3 files changed, 378 insertions(+), 8 deletions(-) create mode 100644 scripts/resolution.py diff --git a/scripts/eta.py b/scripts/eta.py index 3d22ae5e..a5a06c80 100644 --- a/scripts/eta.py +++ b/scripts/eta.py @@ -147,7 +147,6 @@ def calibrate_2pix(eta: np.ndarray, photon_pos: np.ndarray, versor: np.ndarray plt.legend() eta_min = min(eta)[0] - print(eta_min) eta_hist_binning = np.linspace(eta_min, 0.5, 101) eta_hist = Histogram1d(eta_hist_binning, xlabel="eta", ylabel="Counts") eta_hist.fill(eta) @@ -156,9 +155,13 @@ def calibrate_2pix(eta: np.ndarray, photon_pos: np.ndarray, versor: np.ndarray from scipy.special import ndtri - sigma_x = model.sigma.value * np.sqrt(2 / np.pi) * np.exp(0.5 * abs(ndtri(eta_binning))**2) + sigma_x = model.sigma.value * np.sqrt(2 * np.pi) * np.exp(0.5 * abs(ndtri(eta_binning))**2) * np.sqrt(eta_binning**2 + (1 - eta_binning)**2) * 30 / 1600 plt.figure("dr_vs_eta_2pix_derivative") - plt.plot(eta_binning, sigma_x, label="d(dr/p)/d(eta)") + plt.plot(model(eta_binning), sigma_x, label="d(dr/p)/d(eta)") + + plt.figure("dr_distr") + plt.hist(model(eta.astype(float)), bins=100, label="dr / p") + plt.legend() return model diff --git a/scripts/resolution.py b/scripts/resolution.py new file mode 100644 index 00000000..ed34847e --- /dev/null +++ b/scripts/resolution.py @@ -0,0 +1,357 @@ +import argparse +from pathlib import Path + +import numpy as np +from aptapy.hist import Histogram1d, Histogram2d +from aptapy.plotting import plt + +from hexsample.fileio import ReconInputFile +from hexsample.hexagon import HexagonalGrid, HexagonalLayout +from hexsample.pipeline import reconstruct, simulate + +__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") + +RESOLUTION_DIR = Path.home() / "hexsampledata" / "resolution" +if not RESOLUTION_DIR.exists(): + RESOLUTION_DIR.mkdir(parents=True, exist_ok=True) + + +def distance(input_file_path: str, grid): + pitch = grid.pitch + input_file = ReconInputFile(input_file_path) + # Monte Carlo true positions + x_mc = input_file.mc_column("absx") + y_mc = input_file.mc_column("absy") + # Monte Carlo true pixel centers + x0, y0 = grid.pixel_to_world(*grid.world_to_pixel(x_mc, y_mc)) + x = input_file.column("posx") + y = input_file.column("posy") + # Calulate the reconstructed position from the center of the pixel + dr0 = np.sqrt((x - x0) ** 2 + (y - y0) ** 2) / pitch + # Calculate the residual from the Monte Carlo true position + dr = np.sqrt((x - x_mc) ** 2 + (y - y_mc) ** 2) / pitch + input_file.close() + + return dr, dr0 + +def _estimate_loc(hist: Histogram2d) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Estimate the location of the distribution in each slice of the histogram. + + Arguments + --------- + hist : Histogram2d + The 2D histogram to analyze. + + Returns + ------- + centers : np.ndarray + The bin centers of the first axis. + loc : np.ndarray + The estimated location of the distribution in each slice. + err : np.ndarray + The estimated uncertainty on the location in each slice. + """ + centers = hist.bin_centers(axis=0) + loc = np.zeros(centers.shape) + err = np.zeros(centers.shape) + for i in range(len(centers)): + slice_ = hist.slice1d(i) + n = slice_.content.sum() + if n == 0: + loc[i] = np.nan + err[i] = np.nan + continue + # Estimate the location as the median of the distribution and its error + x = np.repeat(slice_.bin_centers(), slice_.content.astype(int)) + loc[i] = np.mean(x) + err[i] = np.std(x) / np.sqrt(n) # Approximate std error of median + # Remove NaN entries + mask = ~np.isnan(loc) + centers = centers[mask] + loc = loc[mask] + err = err[mask] + return centers, loc, err + + +def resolution(**kwargs): + # Define the simulation file path + enc = kwargs["enc"] + zero_sup_threshold = kwargs["zero_sup_threshold"] + if enc is None or zero_sup_threshold is None: + raise RuntimeError("Both enc and zero_sup_threshold arguments are required.") + file_prefix = f"simulation_resolution_{enc}enc_zsup{zero_sup_threshold}_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, + # This is the cause of the resolution degradation. If in simulate it is set to 30, + # then some events are reconstructed very badly. Furthermore, when setting it to 30, + # if the reconstruction is done with a zero suppression threshold of 0, the results + # are diffrent from the ones obtained with a threshold of 30, all without noise. + zero_sup_threshold=0, + readout_mode="circular", + pitch=0.005, + layout=HexagonalLayout.ODD_R, + num_cols=304, + num_rows=352, + gain=1. + ) + # Reconstruct the simulated file with different algorithms + # We start with centroid for all pixels + all_pixel_centroid_prefix = f"recon_allpix_centroid" + all_pixel_centroid = RESOLUTION_DIR / f"{file_prefix}_{all_pixel_centroid_prefix}.h5" + if not all_pixel_centroid.exists(): + reconstruct( + input_file=str(simulation_path), + suffix=all_pixel_centroid_prefix, + zero_sup_threshold=zero_sup_threshold, + max_neighbors=6, + pos_recon_algorithm="centroid" + ) + # All pixels reconstructed with the best algorithm (eta for 2 and 3, centroid otherwise) + all_pixel_best_prefix = f"recon_allpix_best" + all_pixel_best = RESOLUTION_DIR / f"{file_prefix}_{all_pixel_best_prefix}.h5" + if not all_pixel_best.exists(): + reconstruct( + input_file=str(simulation_path), + suffix=all_pixel_best_prefix, + zero_sup_threshold=zero_sup_threshold, + max_neighbors=6, + pos_recon_algorithm="eta" + ) + # Only one-pixel events with centroid + one_pixel_centroid_prefix = f"recon_1pix_centroid" + one_pixel_centroid = RESOLUTION_DIR / f"{file_prefix}_{one_pixel_centroid_prefix}.h5" + if not one_pixel_centroid.exists(): + reconstruct( + input_file=str(simulation_path), + suffix=one_pixel_centroid_prefix, + zero_sup_threshold=zero_sup_threshold, + max_neighbors=0, + pos_recon_algorithm="centroid" + ) + # Only two pixel events with centroid + two_pixel_centroid_prefix = f"recon_2pix_centroid" + two_pixel_centroid = RESOLUTION_DIR / f"{file_prefix}_{two_pixel_centroid_prefix}.h5" + if not two_pixel_centroid.exists(): + reconstruct( + input_file=str(simulation_path), + suffix=two_pixel_centroid_prefix, + zero_sup_threshold=zero_sup_threshold, + num_neighbors=1, + pos_recon_algorithm="centroid" + ) + # Only two-pixel events with eta + two_pixel_eta_prefix = f"recon_2pix_eta" + two_pixel_eta = RESOLUTION_DIR / f"{file_prefix}_{two_pixel_eta_prefix}.h5" + if not two_pixel_eta.exists(): + reconstruct( + input_file=str(simulation_path), + suffix=two_pixel_eta_prefix, + zero_sup_threshold=zero_sup_threshold, + num_neighbors=1, + pos_recon_algorithm="eta" + ) + # Only three-pixel events with centroid + three_pixel_centroid_prefix = f"recon_3pix_centroid" + three_pixel_centroid = RESOLUTION_DIR / f"{file_prefix}_{three_pixel_centroid_prefix}.h5" + if not three_pixel_centroid.exists(): + reconstruct( + input_file=str(simulation_path), + suffix=three_pixel_centroid_prefix, + zero_sup_threshold=zero_sup_threshold, + num_neighbors=2, + pos_recon_algorithm="centroid" + ) + # Only three-pixel events with eta + three_pixel_eta_prefix = f"recon_3pix_eta" + three_pixel_eta = RESOLUTION_DIR / f"{file_prefix}_{three_pixel_eta_prefix}.h5" + if not three_pixel_eta.exists(): + reconstruct( + input_file=str(simulation_path), + suffix=three_pixel_eta_prefix, + zero_sup_threshold=zero_sup_threshold, + num_neighbors=2, + pos_recon_algorithm="eta" + ) + # Define the hexagonal grid for pitch calculation and pixel center determination + grid = HexagonalGrid() + dr_binning = np.linspace(0., 0.4, 101) + + xlabel = "Distance residual [pitch normalized]" + ylabel = "EEF" + # 1 pixel events resolution + dr_1pix, _ = distance(str(one_pixel_centroid), grid) + hist_1pix = Histogram1d(dr_binning, xlabel=xlabel) + hist_1pix.fill(dr_1pix) + plt.figure("1pix_residuals") + hew = hist_1pix.ppf(0.5) + plt.plot(dr_binning, hist_1pix.cdf(dr_binning), "-k", label=f"Centroid (HEW: {hew:.2f})") + plt.hlines(0.5, dr_binning[0], dr_binning[-1], colors="0.4", linestyles="-.") + plt.xlabel(xlabel) + plt.ylabel(ylabel) + plt.xlim(dr_binning[0], dr_binning[-1]) + plt.ylim(0, 1) + plt.legend() + + # 2 pixel events resolution with eta and centroid + dr_2pix_eta, dr0_2pix_eta = distance(str(two_pixel_eta), grid) + dr_2pix_centroid, dr0_2pix_centroid = distance(str(two_pixel_centroid), grid) + hist_2pix_eta = Histogram1d(dr_binning, xlabel=xlabel) + hist_2pix_eta.fill(dr_2pix_eta) + hist_2pix_centroid = Histogram1d(dr_binning, xlabel=xlabel) + hist_2pix_centroid.fill(dr_2pix_centroid) + + + + + # # MOMENTANEOUS PLOT FOR TWO PIXEL EVENTS TAILS + # mask = dr_2pix_eta > 0.5 + # pitch = grid.pitch + # input_file = ReconInputFile(two_pixel_eta) + # # Monte Carlo true positions + # x_mc = input_file.mc_column("absx") + # y_mc = input_file.mc_column("absy") + # # Monte Carlo true pixel centers + # x0, y0 = grid.pixel_to_world(*grid.world_to_pixel(x_mc, y_mc)) + # x = (input_file.column("posx") - x0) / pitch + # y = (input_file.column("posy") - y0) / pitch + # id = input_file.column("trigger_id") + # mask = id == id[mask][0] + # plt.figure("2pix_residuals_scatter") + # hist = Histogram2d(np.linspace(-0.6, 0.6, 100), np.linspace(-0.6, 0.6, 100)) + # xx = (x_mc - x0) / pitch + # yy = (y_mc - y0) / pitch + # hist.fill(xx, yy) + # hist.plot() + # plt.scatter(x[mask], y[mask], s=10, alpha=1, color="blue", marker="*") + # plt.scatter(xx[mask], yy[mask], s=10, alpha=1, color="green", marker="*") + # plt.xlabel("Reconstructed x from pixel center [pitch normalized]") + # plt.ylabel("Reconstructed y from pixel center [pitch normalized]") + # print(id[mask], dr_2pix_eta[mask]) + # input_file.close() + + + + plt.figure("2pix_residuals") + hew_eta = hist_2pix_eta.ppf(0.5) + hew_centroid = hist_2pix_centroid.ppf(0.5) + plt.plot(dr_binning, hist_2pix_eta.cdf(dr_binning), "-k", label=f"$\\eta$ (HEW: {hew_eta:.2f})") + plt.plot(dr_binning, hist_2pix_centroid.cdf(dr_binning), "--k", label=f"Centroid (HEW: {hew_centroid:.2f})") + plt.hlines(0.5, dr_binning[0], dr_binning[-1], colors="0.4", linestyles="-.") + plt.xlabel(xlabel) + plt.ylabel(ylabel) + plt.xlim(dr_binning[0], dr_binning[-1]) + plt.ylim(0, 1) + plt.legend() + + # 3 pixel events resolution with eta and centroid + dr_3pix_eta, dr0_3pix_eta = distance(str(three_pixel_eta), grid) + dr_3pix_centroid, dr0_3pix_centroid = distance(str(three_pixel_centroid), grid) + hist_3pix_eta = Histogram1d(dr_binning, xlabel=xlabel) + hist_3pix_eta.fill(dr_3pix_eta) + hist_3pix_centroid = Histogram1d(dr_binning, xlabel=xlabel) + hist_3pix_centroid.fill(dr_3pix_centroid) + plt.figure("3pix_residuals") + hew_3pix_eta = hist_3pix_eta.ppf(0.5) + hew_3pix_centroid = hist_3pix_centroid.ppf(0.5) + plt.plot(dr_binning, hist_3pix_eta.cdf(dr_binning), "-k", label=f"$\\eta$ (HEW: {hew_3pix_eta:.2f})") + plt.plot(dr_binning, hist_3pix_centroid.cdf(dr_binning), "--k", label=f"Centroid (HEW: {hew_3pix_centroid:.2f})") + plt.hlines(0.5, dr_binning[0], dr_binning[-1], colors="0.4", linestyles="-.") + plt.xlabel(xlabel) + plt.ylabel(ylabel) + plt.xlim(dr_binning[0], dr_binning[-1]) + plt.ylim(0, 1) + plt.legend() + + # All pixel events resolution with best algorithm and centroid + dr_allpix_best, _ = distance(str(all_pixel_best), grid) + dr_allpix_centroid, _ = distance(str(all_pixel_centroid), grid) + hist_allpix_best = Histogram1d(dr_binning, xlabel=xlabel) + hist_allpix_best.fill(dr_allpix_best) + hist_allpix_centroid = Histogram1d(dr_binning, xlabel=xlabel) + hist_allpix_centroid.fill(dr_allpix_centroid) + plt.figure("allpix_residuals") + hew_allpix_best = hist_allpix_best.ppf(0.5) + hew_allpix_centroid = hist_allpix_centroid.ppf(0.5) + plt.plot(dr_binning, hist_allpix_best.cdf(dr_binning), "-k", label=f"Best (HEW: {hew_allpix_best:.2f})") + plt.plot(dr_binning, hist_allpix_centroid.cdf(dr_binning), "--k", label=f"Centroid (HEW: {hew_allpix_centroid:.2f})") + plt.hlines(0.5, dr_binning[0], dr_binning[-1], colors="0.4", linestyles="-.") + plt.xlabel(xlabel) + plt.ylabel(ylabel) + plt.xlim(dr_binning[0], dr_binning[-1]) + plt.ylim(0, 1) + plt.legend() + + # Now we study the dependence of the resolution on the recon distance from the pixel center + # For two pixel events eta and centroid + hist_2d_2pix_eta = Histogram2d( + dr_binning, dr_binning, + xlabel="Reconstructed distance from pixel center [pitch normalized]", + ylabel="Distance residual [pitch normalized]") + hist_2d_2pix_eta.fill(dr0_2pix_eta, dr_2pix_eta) + plt.figure("2pix_eta_residuals_vs_dr0") + hist_2d_2pix_eta.plot() + hist_2d_2pix_centroid = Histogram2d( + dr_binning, dr_binning, + xlabel="Reconstructed distance from pixel center [pitch normalized]", + ylabel="Distance residual [pitch normalized]") + hist_2d_2pix_centroid.fill(dr0_2pix_centroid, dr_2pix_centroid) + plt.figure("2pix_centroid_residuals_vs_dr0") + hist_2d_2pix_centroid.plot() + # Take the slices and calculate the mean and stddev + centers_2pix_eta, loc_2pix_eta, err_2pix_eta = _estimate_loc(hist_2d_2pix_eta) + centers_2pix_centroid, loc_2pix_centroid, err_2pix_centroid = _estimate_loc(hist_2d_2pix_centroid) + plt.figure("2pix_resolution_distance_dependency") + plt.errorbar( + centers_2pix_eta, loc_2pix_eta, yerr=err_2pix_eta, + fmt=".r", label="2pix ETA") + plt.errorbar( + centers_2pix_centroid, loc_2pix_centroid, yerr=err_2pix_centroid, + fmt=".b", label="2pix CENTROID") + plt.xlabel("Reconstructed distance from pixel center [pitch normalized]") + plt.ylabel("Mean distance residual [pitch normalized]") + plt.legend() + + # For three pixel events eta and centroid + hist_2d_3pix_eta = Histogram2d( + dr_binning, dr_binning, + xlabel="Reconstructed distance from pixel center [pitch normalized]", + ylabel="Distance residual [pitch normalized]") + hist_2d_3pix_eta.fill(dr0_3pix_eta, dr_3pix_eta) + plt.figure("3pix_eta_residuals_vs_dr0") + hist_2d_3pix_eta.plot() + hist_2d_3pix_centroid = Histogram2d( + dr_binning, dr_binning, + xlabel="Reconstructed distance from pixel center [pitch normalized]", + ylabel="Distance residual [pitch normalized]") + hist_2d_3pix_centroid.fill(dr0_3pix_centroid, dr_3pix_centroid) + plt.figure("3pix_centroid_residuals_vs_dr0") + hist_2d_3pix_centroid.plot() + # Take the slices and calculate the mean and stddev + centers_3pix_eta, loc_3pix_eta, err_3pix_eta = _estimate_loc(hist_2d_3pix_eta) + centers_3pix_centroid, loc_3pix_centroid, err_3pix_centroid = _estimate_loc(hist_2d_3pix_centroid) + plt.figure("3pix_resolution_distance_dependency") + plt.errorbar( + centers_3pix_eta, loc_3pix_eta, yerr=err_3pix_eta, + fmt=".r", label="3pix ETA") + plt.errorbar( + centers_3pix_centroid, loc_3pix_centroid, yerr=err_3pix_centroid, + fmt=".b", label="3pix CENTROID") + plt.xlabel("Reconstructed distance from pixel center [pitch normalized]") + plt.ylabel("Mean distance residual [pitch normalized]") + plt.legend() + + +if __name__ == "__main__": + resolution(**vars(HXETA_ARGPARSER.parse_args())) + plt.show() diff --git a/src/hexsample/tasks.py b/src/hexsample/tasks.py index 98ebd4d5..13bc364f 100644 --- a/src/hexsample/tasks.py +++ b/src/hexsample/tasks.py @@ -359,22 +359,32 @@ def quicklook(input_file_path: str) -> None: histy = Histogram1d(binning, xlabel=r"$y - y_{MC}$ [cm]").fill(y-y_mc) histy.plot() plt.figure("distance resolution") - dr = np.sqrt((x - x_mc)**2 + (y - y_mc)**2) + dr = np.sqrt((x - x_mc)**2 + (y - y_mc)**2) / 0.005 binning = np.linspace(dr.min(), dr.max(), 100) - histdr = Histogram1d(binning, xlabel=r"$\sqrt{(x - x_{MC})^2 + (y - y_{MC})^2}$ [cm]").fill(dr) + histdr = Histogram1d(binning, xlabel=r"$\sqrt{(x - x_{MC})^2 + (y - y_{MC})^2}$ [cm]") + histdr.fill(dr) + print(f"dr histogram FWHM: {histdr.fwhm()}") + histdr.plot() from .hexagon import HexagonalGrid grid = HexagonalGrid() x0, y0 = grid.pixel_to_world(*grid.world_to_pixel(x_mc, y_mc)) dr_abs = np.sqrt((x - x0)**2 + (y - y0)**2) / grid.pitch - bins = np.linspace(0, 1, 100) + bins = np.linspace(0, 0.7, 200) hist = Histogram2d(bins, bins) # I need the recon distance from the hit pixel center - hist.fill(dr_abs, dr / grid.pitch) + hist.fill(dr_abs, dr) + plt.figure() + hist.plot() hist_mean, hist_sigma = hist.project_statistics() plt.figure("Reconstructed vs true distance from pixel center") - plt.plot(hist_mean.bin_centers(), hist_mean.content.flatten(), 'o', label='Mean') + yy = hist_mean.content.flatten() + mask = yy > 0 + yy = yy[mask] + xx = hist_mean.bin_centers()[mask] + yy_err = hist_mean.errors.flatten()[mask] + plt.errorbar(xx, yy, fmt='.k', label='Mean') input_file.close() plt.show() From 045bdcca273b4ab8b7bf5fdd04f8360210e1e6e4 Mon Sep 17 00:00:00 2001 From: Augusto Cattafesta Date: Wed, 4 Feb 2026 11:47:37 +0100 Subject: [PATCH 07/18] Minor. --- scripts/resolution.py | 79 +++++++++++++++++-------------------------- 1 file changed, 31 insertions(+), 48 deletions(-) diff --git a/scripts/resolution.py b/scripts/resolution.py index ed34847e..11330bcd 100644 --- a/scripts/resolution.py +++ b/scripts/resolution.py @@ -15,11 +15,16 @@ 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 distance(input_file_path: str, grid): pitch = grid.pitch @@ -39,6 +44,7 @@ def distance(input_file_path: str, grid): return dr, dr0 + def _estimate_loc(hist: Histogram2d) -> tuple[np.ndarray, np.ndarray, np.ndarray]: """Estimate the location of the distribution in each slice of the histogram. @@ -84,7 +90,8 @@ def resolution(**kwargs): zero_sup_threshold = kwargs["zero_sup_threshold"] if enc is None or zero_sup_threshold is None: raise RuntimeError("Both enc and zero_sup_threshold arguments are required.") - file_prefix = f"simulation_resolution_{enc}enc_zsup{zero_sup_threshold}_hexagonal" + zero_sup_thr_simulation = 0 + file_prefix = f"simulation_resolution_{enc}enc_zsupsim{zero_sup_thr_simulation}_hexagonal" simulation_path = RESOLUTION_DIR / f"{file_prefix}.h5" # Run simulation if the output file does not already exist if not simulation_path.exists(): @@ -93,21 +100,22 @@ def resolution(**kwargs): output_file=str(simulation_path), beam="hexagonal", enc=enc, - # This is the cause of the resolution degradation. If in simulate it is set to 30, + # This may be the cause of the resolution degradation. If in simulate it is set to 30, # then some events are reconstructed very badly. Furthermore, when setting it to 30, # if the reconstruction is done with a zero suppression threshold of 0, the results # are diffrent from the ones obtained with a threshold of 30, all without noise. - zero_sup_threshold=0, + zero_sup_threshold=zero_sup_thr_simulation, readout_mode="circular", pitch=0.005, layout=HexagonalLayout.ODD_R, num_cols=304, num_rows=352, - gain=1. + gain=1., + random_seed=0 ) # Reconstruct the simulated file with different algorithms # We start with centroid for all pixels - all_pixel_centroid_prefix = f"recon_allpix_centroid" + all_pixel_centroid_prefix = f"recon_zsuprec{zero_sup_threshold}_allpix_centroid" all_pixel_centroid = RESOLUTION_DIR / f"{file_prefix}_{all_pixel_centroid_prefix}.h5" if not all_pixel_centroid.exists(): reconstruct( @@ -118,7 +126,7 @@ def resolution(**kwargs): pos_recon_algorithm="centroid" ) # All pixels reconstructed with the best algorithm (eta for 2 and 3, centroid otherwise) - all_pixel_best_prefix = f"recon_allpix_best" + all_pixel_best_prefix = f"recon_zsuprec{zero_sup_threshold}_allpix_best" all_pixel_best = RESOLUTION_DIR / f"{file_prefix}_{all_pixel_best_prefix}.h5" if not all_pixel_best.exists(): reconstruct( @@ -129,7 +137,7 @@ def resolution(**kwargs): pos_recon_algorithm="eta" ) # Only one-pixel events with centroid - one_pixel_centroid_prefix = f"recon_1pix_centroid" + one_pixel_centroid_prefix = f"recon_zsuprec{zero_sup_threshold}_1pix_centroid" one_pixel_centroid = RESOLUTION_DIR / f"{file_prefix}_{one_pixel_centroid_prefix}.h5" if not one_pixel_centroid.exists(): reconstruct( @@ -140,7 +148,7 @@ def resolution(**kwargs): pos_recon_algorithm="centroid" ) # Only two pixel events with centroid - two_pixel_centroid_prefix = f"recon_2pix_centroid" + two_pixel_centroid_prefix = f"recon_zsuprec{zero_sup_threshold}_2pix_centroid" two_pixel_centroid = RESOLUTION_DIR / f"{file_prefix}_{two_pixel_centroid_prefix}.h5" if not two_pixel_centroid.exists(): reconstruct( @@ -151,7 +159,7 @@ def resolution(**kwargs): pos_recon_algorithm="centroid" ) # Only two-pixel events with eta - two_pixel_eta_prefix = f"recon_2pix_eta" + two_pixel_eta_prefix = f"recon_zsuprec{zero_sup_threshold}_2pix_eta" two_pixel_eta = RESOLUTION_DIR / f"{file_prefix}_{two_pixel_eta_prefix}.h5" if not two_pixel_eta.exists(): reconstruct( @@ -162,7 +170,7 @@ def resolution(**kwargs): pos_recon_algorithm="eta" ) # Only three-pixel events with centroid - three_pixel_centroid_prefix = f"recon_3pix_centroid" + three_pixel_centroid_prefix = f"recon_zsuprec{zero_sup_threshold}_3pix_centroid" three_pixel_centroid = RESOLUTION_DIR / f"{file_prefix}_{three_pixel_centroid_prefix}.h5" if not three_pixel_centroid.exists(): reconstruct( @@ -173,7 +181,7 @@ def resolution(**kwargs): pos_recon_algorithm="centroid" ) # Only three-pixel events with eta - three_pixel_eta_prefix = f"recon_3pix_eta" + three_pixel_eta_prefix = f"recon_zsuprec{zero_sup_threshold}_3pix_eta" three_pixel_eta = RESOLUTION_DIR / f"{file_prefix}_{three_pixel_eta_prefix}.h5" if not three_pixel_eta.exists(): reconstruct( @@ -185,7 +193,7 @@ def resolution(**kwargs): ) # Define the hexagonal grid for pitch calculation and pixel center determination grid = HexagonalGrid() - dr_binning = np.linspace(0., 0.4, 101) + dr_binning = np.linspace(0., 0.6, 101) xlabel = "Distance residual [pitch normalized]" ylabel = "EEF" @@ -193,7 +201,7 @@ def resolution(**kwargs): dr_1pix, _ = distance(str(one_pixel_centroid), grid) hist_1pix = Histogram1d(dr_binning, xlabel=xlabel) hist_1pix.fill(dr_1pix) - plt.figure("1pix_residuals") + fig_1pix = plt.figure("1pix_residuals") hew = hist_1pix.ppf(0.5) plt.plot(dr_binning, hist_1pix.cdf(dr_binning), "-k", label=f"Centroid (HEW: {hew:.2f})") plt.hlines(0.5, dr_binning[0], dr_binning[-1], colors="0.4", linestyles="-.") @@ -210,39 +218,7 @@ def resolution(**kwargs): hist_2pix_eta.fill(dr_2pix_eta) hist_2pix_centroid = Histogram1d(dr_binning, xlabel=xlabel) hist_2pix_centroid.fill(dr_2pix_centroid) - - - - - # # MOMENTANEOUS PLOT FOR TWO PIXEL EVENTS TAILS - # mask = dr_2pix_eta > 0.5 - # pitch = grid.pitch - # input_file = ReconInputFile(two_pixel_eta) - # # Monte Carlo true positions - # x_mc = input_file.mc_column("absx") - # y_mc = input_file.mc_column("absy") - # # Monte Carlo true pixel centers - # x0, y0 = grid.pixel_to_world(*grid.world_to_pixel(x_mc, y_mc)) - # x = (input_file.column("posx") - x0) / pitch - # y = (input_file.column("posy") - y0) / pitch - # id = input_file.column("trigger_id") - # mask = id == id[mask][0] - # plt.figure("2pix_residuals_scatter") - # hist = Histogram2d(np.linspace(-0.6, 0.6, 100), np.linspace(-0.6, 0.6, 100)) - # xx = (x_mc - x0) / pitch - # yy = (y_mc - y0) / pitch - # hist.fill(xx, yy) - # hist.plot() - # plt.scatter(x[mask], y[mask], s=10, alpha=1, color="blue", marker="*") - # plt.scatter(xx[mask], yy[mask], s=10, alpha=1, color="green", marker="*") - # plt.xlabel("Reconstructed x from pixel center [pitch normalized]") - # plt.ylabel("Reconstructed y from pixel center [pitch normalized]") - # print(id[mask], dr_2pix_eta[mask]) - # input_file.close() - - - - plt.figure("2pix_residuals") + fig_2pix = plt.figure("2pix_residuals") hew_eta = hist_2pix_eta.ppf(0.5) hew_centroid = hist_2pix_centroid.ppf(0.5) plt.plot(dr_binning, hist_2pix_eta.cdf(dr_binning), "-k", label=f"$\\eta$ (HEW: {hew_eta:.2f})") @@ -261,7 +237,7 @@ def resolution(**kwargs): hist_3pix_eta.fill(dr_3pix_eta) hist_3pix_centroid = Histogram1d(dr_binning, xlabel=xlabel) hist_3pix_centroid.fill(dr_3pix_centroid) - plt.figure("3pix_residuals") + fig_3pix = plt.figure("3pix_residuals") hew_3pix_eta = hist_3pix_eta.ppf(0.5) hew_3pix_centroid = hist_3pix_centroid.ppf(0.5) plt.plot(dr_binning, hist_3pix_eta.cdf(dr_binning), "-k", label=f"$\\eta$ (HEW: {hew_3pix_eta:.2f})") @@ -280,7 +256,7 @@ def resolution(**kwargs): hist_allpix_best.fill(dr_allpix_best) hist_allpix_centroid = Histogram1d(dr_binning, xlabel=xlabel) hist_allpix_centroid.fill(dr_allpix_centroid) - plt.figure("allpix_residuals") + fig_allpix = plt.figure("allpix_residuals") hew_allpix_best = hist_allpix_best.ppf(0.5) hew_allpix_centroid = hist_allpix_centroid.ppf(0.5) plt.plot(dr_binning, hist_allpix_best.cdf(dr_binning), "-k", label=f"Best (HEW: {hew_allpix_best:.2f})") @@ -351,6 +327,13 @@ def resolution(**kwargs): plt.ylabel("Mean distance residual [pitch normalized]") plt.legend() + if kwargs["save"]: + fig_format = "png" + fig_1pix.savefig(FIGURES_DIR / f"1pix_resolution_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) + fig_2pix.savefig(FIGURES_DIR / f"2pix_resolution_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) + fig_3pix.savefig(FIGURES_DIR / f"3pix_resolution_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) + fig_allpix.savefig(FIGURES_DIR / f"allpix_resolution_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) + if __name__ == "__main__": resolution(**vars(HXETA_ARGPARSER.parse_args())) From d393ff7cca4058916fda42cc898aeb1494f46b15 Mon Sep 17 00:00:00 2001 From: Augusto Cattafesta Date: Wed, 4 Feb 2026 14:04:19 +0100 Subject: [PATCH 08/18] Minor. --- src/hexsample/pipeline.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/hexsample/pipeline.py b/src/hexsample/pipeline.py index efa0f3f0..f062249b 100644 --- a/src/hexsample/pipeline.py +++ b/src/hexsample/pipeline.py @@ -54,9 +54,8 @@ def reconstruct(**kwargs) -> str: 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) - eta_3pix_theta1 = kwargs.get("eta_3pix_theta1", defaults.eta_3pix_theta1) 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, \ + pos_recon_algorithm, eta_2pix_rad, eta_3pix_rad0, eta_3pix_rad1, eta_3pix_theta0 return tasks.reconstruct(*args, kwargs) From cebf311ecacdef139ca92ee1918fc46f31170959 Mon Sep 17 00:00:00 2001 From: Augusto Cattafesta Date: Thu, 5 Feb 2026 16:35:30 +0100 Subject: [PATCH 09/18] Resolution module and script refactoring. --- scripts/eta.py | 7 +- scripts/resolution.py | 339 ++++++++++++------------------------ src/hexsample/resolution.py | 166 ++++++++++++++++++ 3 files changed, 284 insertions(+), 228 deletions(-) create mode 100644 src/hexsample/resolution.py diff --git a/scripts/eta.py b/scripts/eta.py index 8b881728..6a89ba2b 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") @@ -299,12 +301,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 index 11330bcd..79932b6e 100644 --- a/scripts/resolution.py +++ b/scripts/resolution.py @@ -8,13 +8,14 @@ from hexsample.fileio import ReconInputFile from hexsample.hexagon import HexagonalGrid, HexagonalLayout from hexsample.pipeline import reconstruct, simulate +from hexsample.resolution import eef, hew __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("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" @@ -26,25 +27,6 @@ FIGURES_DIR.mkdir(parents=True, exist_ok=True) -def distance(input_file_path: str, grid): - pitch = grid.pitch - input_file = ReconInputFile(input_file_path) - # Monte Carlo true positions - x_mc = input_file.mc_column("absx") - y_mc = input_file.mc_column("absy") - # Monte Carlo true pixel centers - x0, y0 = grid.pixel_to_world(*grid.world_to_pixel(x_mc, y_mc)) - x = input_file.column("posx") - y = input_file.column("posy") - # Calulate the reconstructed position from the center of the pixel - dr0 = np.sqrt((x - x0) ** 2 + (y - y0) ** 2) / pitch - # Calculate the residual from the Monte Carlo true position - dr = np.sqrt((x - x_mc) ** 2 + (y - y_mc) ** 2) / pitch - input_file.close() - - return dr, dr0 - - def _estimate_loc(hist: Histogram2d) -> tuple[np.ndarray, np.ndarray, np.ndarray]: """Estimate the location of the distribution in each slice of the histogram. @@ -85,254 +67,159 @@ def _estimate_loc(hist: Histogram2d) -> tuple[np.ndarray, np.ndarray, np.ndarray def resolution(**kwargs): - # Define the simulation file path enc = kwargs["enc"] zero_sup_threshold = kwargs["zero_sup_threshold"] - if enc is None or zero_sup_threshold is None: - raise RuntimeError("Both enc and zero_sup_threshold arguments are required.") - zero_sup_thr_simulation = 0 - file_prefix = f"simulation_resolution_{enc}enc_zsupsim{zero_sup_thr_simulation}_hexagonal" + 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, + num_events=10000, output_file=str(simulation_path), beam="hexagonal", enc=enc, - # This may be the cause of the resolution degradation. If in simulate it is set to 30, - # then some events are reconstructed very badly. Furthermore, when setting it to 30, - # if the reconstruction is done with a zero suppression threshold of 0, the results - # are diffrent from the ones obtained with a threshold of 30, all without noise. - zero_sup_threshold=zero_sup_thr_simulation, + zero_sup_threshold=0, readout_mode="circular", pitch=0.005, layout=HexagonalLayout.ODD_R, num_cols=304, num_rows=352, gain=1., - random_seed=0 ) # Reconstruct the simulated file with different algorithms - # We start with centroid for all pixels - all_pixel_centroid_prefix = f"recon_zsuprec{zero_sup_threshold}_allpix_centroid" - all_pixel_centroid = RESOLUTION_DIR / f"{file_prefix}_{all_pixel_centroid_prefix}.h5" - if not all_pixel_centroid.exists(): + 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( input_file=str(simulation_path), - suffix=all_pixel_centroid_prefix, + suffix=centroid_prefix, zero_sup_threshold=zero_sup_threshold, max_neighbors=6, pos_recon_algorithm="centroid" ) # All pixels reconstructed with the best algorithm (eta for 2 and 3, centroid otherwise) - all_pixel_best_prefix = f"recon_zsuprec{zero_sup_threshold}_allpix_best" - all_pixel_best = RESOLUTION_DIR / f"{file_prefix}_{all_pixel_best_prefix}.h5" - if not all_pixel_best.exists(): + 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( input_file=str(simulation_path), - suffix=all_pixel_best_prefix, + suffix=best_prefix, zero_sup_threshold=zero_sup_threshold, max_neighbors=6, pos_recon_algorithm="eta" ) - # Only one-pixel events with centroid - one_pixel_centroid_prefix = f"recon_zsuprec{zero_sup_threshold}_1pix_centroid" - one_pixel_centroid = RESOLUTION_DIR / f"{file_prefix}_{one_pixel_centroid_prefix}.h5" - if not one_pixel_centroid.exists(): - reconstruct( - input_file=str(simulation_path), - suffix=one_pixel_centroid_prefix, - zero_sup_threshold=zero_sup_threshold, - max_neighbors=0, - pos_recon_algorithm="centroid" - ) - # Only two pixel events with centroid - two_pixel_centroid_prefix = f"recon_zsuprec{zero_sup_threshold}_2pix_centroid" - two_pixel_centroid = RESOLUTION_DIR / f"{file_prefix}_{two_pixel_centroid_prefix}.h5" - if not two_pixel_centroid.exists(): - reconstruct( - input_file=str(simulation_path), - suffix=two_pixel_centroid_prefix, - zero_sup_threshold=zero_sup_threshold, - num_neighbors=1, - pos_recon_algorithm="centroid" - ) - # Only two-pixel events with eta - two_pixel_eta_prefix = f"recon_zsuprec{zero_sup_threshold}_2pix_eta" - two_pixel_eta = RESOLUTION_DIR / f"{file_prefix}_{two_pixel_eta_prefix}.h5" - if not two_pixel_eta.exists(): - reconstruct( - input_file=str(simulation_path), - suffix=two_pixel_eta_prefix, - zero_sup_threshold=zero_sup_threshold, - num_neighbors=1, - pos_recon_algorithm="eta" - ) - # Only three-pixel events with centroid - three_pixel_centroid_prefix = f"recon_zsuprec{zero_sup_threshold}_3pix_centroid" - three_pixel_centroid = RESOLUTION_DIR / f"{file_prefix}_{three_pixel_centroid_prefix}.h5" - if not three_pixel_centroid.exists(): - reconstruct( - input_file=str(simulation_path), - suffix=three_pixel_centroid_prefix, - zero_sup_threshold=zero_sup_threshold, - num_neighbors=2, - pos_recon_algorithm="centroid" - ) - # Only three-pixel events with eta - three_pixel_eta_prefix = f"recon_zsuprec{zero_sup_threshold}_3pix_eta" - three_pixel_eta = RESOLUTION_DIR / f"{file_prefix}_{three_pixel_eta_prefix}.h5" - if not three_pixel_eta.exists(): - reconstruct( - input_file=str(simulation_path), - suffix=three_pixel_eta_prefix, - zero_sup_threshold=zero_sup_threshold, - num_neighbors=2, - pos_recon_algorithm="eta" - ) - # Define the hexagonal grid for pitch calculation and pixel center determination - grid = HexagonalGrid() - dr_binning = np.linspace(0., 0.6, 101) - + # Study the EEF + x = np.linspace(0, 0.6, 101) xlabel = "Distance residual [pitch normalized]" - ylabel = "EEF" - # 1 pixel events resolution - dr_1pix, _ = distance(str(one_pixel_centroid), grid) - hist_1pix = Histogram1d(dr_binning, xlabel=xlabel) - hist_1pix.fill(dr_1pix) - fig_1pix = plt.figure("1pix_residuals") - hew = hist_1pix.ppf(0.5) - plt.plot(dr_binning, hist_1pix.cdf(dr_binning), "-k", label=f"Centroid (HEW: {hew:.2f})") - plt.hlines(0.5, dr_binning[0], dr_binning[-1], colors="0.4", linestyles="-.") - plt.xlabel(xlabel) - plt.ylabel(ylabel) - plt.xlim(dr_binning[0], dr_binning[-1]) + ylabel = "Encircled Energy Fraction" + # Plot the centroid eef for all cluster sizes + centroid_recon_file = ReconInputFile(str(centroid_path)) + plt.figure(f"centroid_eef_{enc}enc_{zero_sup_threshold}zsup") + plt.plot(x, eef(x, centroid_recon_file, 0), + label=f"1 pix (HEW={hew(centroid_recon_file, 0):.2f})", + linestyle="-", color="black") + plt.plot(x, eef(x, centroid_recon_file, 1), + label=f"2 pix (HEW={hew(centroid_recon_file, 1):.2f})", + linestyle="--", color="black") + plt.plot(x, eef(x, centroid_recon_file, 2), + label=f"3 pix (HEW={hew(centroid_recon_file, 2):.2f})", + linestyle=":", color="black") + plt.plot(x, eef(x, centroid_recon_file, max_neighbors=6), + label=f"All pix (HEW={hew(centroid_recon_file, max_neighbors=6):.2f})", + linestyle="-.", color="black") + plt.xlim(x[0], x[-1]) plt.ylim(0, 1) - plt.legend() - - # 2 pixel events resolution with eta and centroid - dr_2pix_eta, dr0_2pix_eta = distance(str(two_pixel_eta), grid) - dr_2pix_centroid, dr0_2pix_centroid = distance(str(two_pixel_centroid), grid) - hist_2pix_eta = Histogram1d(dr_binning, xlabel=xlabel) - hist_2pix_eta.fill(dr_2pix_eta) - hist_2pix_centroid = Histogram1d(dr_binning, xlabel=xlabel) - hist_2pix_centroid.fill(dr_2pix_centroid) - fig_2pix = plt.figure("2pix_residuals") - hew_eta = hist_2pix_eta.ppf(0.5) - hew_centroid = hist_2pix_centroid.ppf(0.5) - plt.plot(dr_binning, hist_2pix_eta.cdf(dr_binning), "-k", label=f"$\\eta$ (HEW: {hew_eta:.2f})") - plt.plot(dr_binning, hist_2pix_centroid.cdf(dr_binning), "--k", label=f"Centroid (HEW: {hew_centroid:.2f})") - plt.hlines(0.5, dr_binning[0], dr_binning[-1], colors="0.4", linestyles="-.") plt.xlabel(xlabel) plt.ylabel(ylabel) - plt.xlim(dr_binning[0], dr_binning[-1]) - plt.ylim(0, 1) plt.legend() - - # 3 pixel events resolution with eta and centroid - dr_3pix_eta, dr0_3pix_eta = distance(str(three_pixel_eta), grid) - dr_3pix_centroid, dr0_3pix_centroid = distance(str(three_pixel_centroid), grid) - hist_3pix_eta = Histogram1d(dr_binning, xlabel=xlabel) - hist_3pix_eta.fill(dr_3pix_eta) - hist_3pix_centroid = Histogram1d(dr_binning, xlabel=xlabel) - hist_3pix_centroid.fill(dr_3pix_centroid) - fig_3pix = plt.figure("3pix_residuals") - hew_3pix_eta = hist_3pix_eta.ppf(0.5) - hew_3pix_centroid = hist_3pix_centroid.ppf(0.5) - plt.plot(dr_binning, hist_3pix_eta.cdf(dr_binning), "-k", label=f"$\\eta$ (HEW: {hew_3pix_eta:.2f})") - plt.plot(dr_binning, hist_3pix_centroid.cdf(dr_binning), "--k", label=f"Centroid (HEW: {hew_3pix_centroid:.2f})") - plt.hlines(0.5, dr_binning[0], dr_binning[-1], colors="0.4", linestyles="-.") - plt.xlabel(xlabel) - plt.ylabel(ylabel) - plt.xlim(dr_binning[0], dr_binning[-1]) + # Plot the eta + centroid eef for all cluster sizes + best_recon_file = ReconInputFile(str(best_path)) + plt.figure(f"best_eef_{enc}enc_{zero_sup_threshold}zsup") + plt.plot(x, eef(x, best_recon_file, 0), + label=f"1 pix (HEW={hew(best_recon_file, 0):.2f})", + linestyle="-", color="black") + plt.plot(x, eef(x, best_recon_file, 1), + label=f"2 pix (HEW={hew(best_recon_file, 1):.2f})", + linestyle="--", color="black") + plt.plot(x, eef(x, best_recon_file, 2), + label=f"3 pix (HEW={hew(best_recon_file, 2):.2f})", + linestyle=":", color="black") + plt.plot(x, eef(x, best_recon_file, max_neighbors=6), + label=f"All pix (HEW={hew(best_recon_file, max_neighbors=6):.2f})", + linestyle="-.", color="black") + plt.xlim(x[0], x[-1]) plt.ylim(0, 1) - plt.legend() - - # All pixel events resolution with best algorithm and centroid - dr_allpix_best, _ = distance(str(all_pixel_best), grid) - dr_allpix_centroid, _ = distance(str(all_pixel_centroid), grid) - hist_allpix_best = Histogram1d(dr_binning, xlabel=xlabel) - hist_allpix_best.fill(dr_allpix_best) - hist_allpix_centroid = Histogram1d(dr_binning, xlabel=xlabel) - hist_allpix_centroid.fill(dr_allpix_centroid) - fig_allpix = plt.figure("allpix_residuals") - hew_allpix_best = hist_allpix_best.ppf(0.5) - hew_allpix_centroid = hist_allpix_centroid.ppf(0.5) - plt.plot(dr_binning, hist_allpix_best.cdf(dr_binning), "-k", label=f"Best (HEW: {hew_allpix_best:.2f})") - plt.plot(dr_binning, hist_allpix_centroid.cdf(dr_binning), "--k", label=f"Centroid (HEW: {hew_allpix_centroid:.2f})") - plt.hlines(0.5, dr_binning[0], dr_binning[-1], colors="0.4", linestyles="-.") plt.xlabel(xlabel) plt.ylabel(ylabel) - plt.xlim(dr_binning[0], dr_binning[-1]) - plt.ylim(0, 1) plt.legend() + centroid_recon_file.close() + best_recon_file.close() # Now we study the dependence of the resolution on the recon distance from the pixel center # For two pixel events eta and centroid - hist_2d_2pix_eta = Histogram2d( - dr_binning, dr_binning, - xlabel="Reconstructed distance from pixel center [pitch normalized]", - ylabel="Distance residual [pitch normalized]") - hist_2d_2pix_eta.fill(dr0_2pix_eta, dr_2pix_eta) - plt.figure("2pix_eta_residuals_vs_dr0") - hist_2d_2pix_eta.plot() - hist_2d_2pix_centroid = Histogram2d( - dr_binning, dr_binning, - xlabel="Reconstructed distance from pixel center [pitch normalized]", - ylabel="Distance residual [pitch normalized]") - hist_2d_2pix_centroid.fill(dr0_2pix_centroid, dr_2pix_centroid) - plt.figure("2pix_centroid_residuals_vs_dr0") - hist_2d_2pix_centroid.plot() - # Take the slices and calculate the mean and stddev - centers_2pix_eta, loc_2pix_eta, err_2pix_eta = _estimate_loc(hist_2d_2pix_eta) - centers_2pix_centroid, loc_2pix_centroid, err_2pix_centroid = _estimate_loc(hist_2d_2pix_centroid) - plt.figure("2pix_resolution_distance_dependency") - plt.errorbar( - centers_2pix_eta, loc_2pix_eta, yerr=err_2pix_eta, - fmt=".r", label="2pix ETA") - plt.errorbar( - centers_2pix_centroid, loc_2pix_centroid, yerr=err_2pix_centroid, - fmt=".b", label="2pix CENTROID") - plt.xlabel("Reconstructed distance from pixel center [pitch normalized]") - plt.ylabel("Mean distance residual [pitch normalized]") - plt.legend() + # hist_2d_2pix_eta = Histogram2d( + # dr_binning, dr_binning, + # xlabel="Reconstructed distance from pixel center [pitch normalized]", + # ylabel="Distance residual [pitch normalized]") + # hist_2d_2pix_eta.fill(dr0_2pix_eta, dr_2pix_eta) + # plt.figure("2pix_eta_residuals_vs_dr0") + # hist_2d_2pix_eta.plot() + # hist_2d_2pix_centroid = Histogram2d( + # dr_binning, dr_binning, + # xlabel="Reconstructed distance from pixel center [pitch normalized]", + # ylabel="Distance residual [pitch normalized]") + # hist_2d_2pix_centroid.fill(dr0_2pix_centroid, dr_2pix_centroid) + # plt.figure("2pix_centroid_residuals_vs_dr0") + # hist_2d_2pix_centroid.plot() + # # Take the slices and calculate the mean and stddev + # centers_2pix_eta, loc_2pix_eta, err_2pix_eta = _estimate_loc(hist_2d_2pix_eta) + # centers_2pix_centroid, loc_2pix_centroid, err_2pix_centroid = _estimate_loc(hist_2d_2pix_centroid) + # plt.figure("2pix_resolution_distance_dependency") + # plt.errorbar( + # centers_2pix_eta, loc_2pix_eta, yerr=err_2pix_eta, + # fmt=".r", label="2pix ETA") + # plt.errorbar( + # centers_2pix_centroid, loc_2pix_centroid, yerr=err_2pix_centroid, + # fmt=".b", label="2pix CENTROID") + # plt.xlabel("Reconstructed distance from pixel center [pitch normalized]") + # plt.ylabel("Mean distance residual [pitch normalized]") + # plt.legend() - # For three pixel events eta and centroid - hist_2d_3pix_eta = Histogram2d( - dr_binning, dr_binning, - xlabel="Reconstructed distance from pixel center [pitch normalized]", - ylabel="Distance residual [pitch normalized]") - hist_2d_3pix_eta.fill(dr0_3pix_eta, dr_3pix_eta) - plt.figure("3pix_eta_residuals_vs_dr0") - hist_2d_3pix_eta.plot() - hist_2d_3pix_centroid = Histogram2d( - dr_binning, dr_binning, - xlabel="Reconstructed distance from pixel center [pitch normalized]", - ylabel="Distance residual [pitch normalized]") - hist_2d_3pix_centroid.fill(dr0_3pix_centroid, dr_3pix_centroid) - plt.figure("3pix_centroid_residuals_vs_dr0") - hist_2d_3pix_centroid.plot() - # Take the slices and calculate the mean and stddev - centers_3pix_eta, loc_3pix_eta, err_3pix_eta = _estimate_loc(hist_2d_3pix_eta) - centers_3pix_centroid, loc_3pix_centroid, err_3pix_centroid = _estimate_loc(hist_2d_3pix_centroid) - plt.figure("3pix_resolution_distance_dependency") - plt.errorbar( - centers_3pix_eta, loc_3pix_eta, yerr=err_3pix_eta, - fmt=".r", label="3pix ETA") - plt.errorbar( - centers_3pix_centroid, loc_3pix_centroid, yerr=err_3pix_centroid, - fmt=".b", label="3pix CENTROID") - plt.xlabel("Reconstructed distance from pixel center [pitch normalized]") - plt.ylabel("Mean distance residual [pitch normalized]") - plt.legend() - - if kwargs["save"]: - fig_format = "png" - fig_1pix.savefig(FIGURES_DIR / f"1pix_resolution_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) - fig_2pix.savefig(FIGURES_DIR / f"2pix_resolution_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) - fig_3pix.savefig(FIGURES_DIR / f"3pix_resolution_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) - fig_allpix.savefig(FIGURES_DIR / f"allpix_resolution_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) + # # For three pixel events eta and centroid + # hist_2d_3pix_eta = Histogram2d( + # dr_binning, dr_binning, + # xlabel="Reconstructed distance from pixel center [pitch normalized]", + # ylabel="Distance residual [pitch normalized]") + # hist_2d_3pix_eta.fill(dr0_3pix_eta, dr_3pix_eta) + # plt.figure("3pix_eta_residuals_vs_dr0") + # hist_2d_3pix_eta.plot() + # hist_2d_3pix_centroid = Histogram2d( + # dr_binning, dr_binning, + # xlabel="Reconstructed distance from pixel center [pitch normalized]", + # ylabel="Distance residual [pitch normalized]") + # hist_2d_3pix_centroid.fill(dr0_3pix_centroid, dr_3pix_centroid) + # plt.figure("3pix_centroid_residuals_vs_dr0") + # hist_2d_3pix_centroid.plot() + # # Take the slices and calculate the mean and stddev + # centers_3pix_eta, loc_3pix_eta, err_3pix_eta = _estimate_loc(hist_2d_3pix_eta) + # centers_3pix_centroid, loc_3pix_centroid, err_3pix_centroid = _estimate_loc(hist_2d_3pix_centroid) + # plt.figure("3pix_resolution_distance_dependency") + # plt.errorbar( + # centers_3pix_eta, loc_3pix_eta, yerr=err_3pix_eta, + # fmt=".r", label="3pix ETA") + # plt.errorbar( + # centers_3pix_centroid, loc_3pix_centroid, yerr=err_3pix_centroid, + # fmt=".b", label="3pix CENTROID") + # plt.xlabel("Reconstructed distance from pixel center [pitch normalized]") + # plt.ylabel("Mean distance residual [pitch normalized]") + # plt.legend() + + # if kwargs["save"]: + # fig_format = "png" + # fig_1pix.savefig(FIGURES_DIR / f"1pix_resolution_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) + # fig_2pix.savefig(FIGURES_DIR / f"2pix_resolution_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) + # fig_3pix.savefig(FIGURES_DIR / f"3pix_resolution_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) + # fig_allpix.savefig(FIGURES_DIR / f"allpix_resolution_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) if __name__ == "__main__": diff --git a/src/hexsample/resolution.py b/src/hexsample/resolution.py new file mode 100644 index 00000000..e8824da9 --- /dev/null +++ b/src/hexsample/resolution.py @@ -0,0 +1,166 @@ +# 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. +""" + +import numpy as np +from aptapy.hist import Histogram1d + +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") + # Calulate the reconstructed distance from the pixel center + dr0 = np.sqrt((x - x0) ** 2 + (y - y0) ** 2) + return dr0 + + +def hist_distance_res(input_file: ReconInputFile, num_neighbors: int = 0, + max_neighbors: int = -1) -> Histogram1d: + """Create the histogram of distance residuals for the given numqber 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") + if max_neighbors >= 0: + mask = size <= max_neighbors + 1 + else: + mask = 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 + # it covers the full range of dr. + 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_res(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_res(input_file, num_neighbors, max_neighbors) + return hist.ppf(0.5) From f67c63c252b957bff0f7e6c1cc052cb61ca713a6 Mon Sep 17 00:00:00 2001 From: Augusto Cattafesta Date: Thu, 5 Feb 2026 17:05:15 +0100 Subject: [PATCH 10/18] Minor. --- scripts/analyze_noise.py | 101 ------------------------------------ scripts/resolution.py | 91 +++++++++++--------------------- src/hexsample/resolution.py | 42 +++++++++++++-- src/hexsample/tasks.py | 55 -------------------- 4 files changed, 66 insertions(+), 223 deletions(-) delete mode 100644 scripts/analyze_noise.py diff --git a/scripts/analyze_noise.py b/scripts/analyze_noise.py deleted file mode 100644 index d31c110e..00000000 --- a/scripts/analyze_noise.py +++ /dev/null @@ -1,101 +0,0 @@ -# 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. - -"""Analyze the noise to find the best the zero suppression threshold. -""" - -import argparse -import pathlib - -import numpy as np -from aptapy.plotting import plt - -from hexsample.fileio import digi_input_file_class, peek_readout_type -from hexsample.pipeline import reconstruct, resolution -from hexsample.readout import HexagonalReadoutMode - -ARGPARSER = argparse.ArgumentParser() -ARGPARSER.add_argument("input_file", type=str, help="path to the input file") -ARGPARSER.add_argument("--num_neighbors", type=int, default=2, - help="number of neighbors for the clustering (default: 2)") -ARGPARSER.add_argument("--pos_recon_algorithm", type=str, default="eta", - help="position reconstruction algorithm (default: eta)") - - -def analyze_noise(**kwargs): - """ Analyze the effect of the noise and zero suppression threshold on the reconstucted position - resolution, by calculating the bias and the standard deviation of the position residuals - distribution for different values of the threshold. - """ - # Open the file to get the header information - input_file_path = str(kwargs["input_file"]) - if not input_file_path.endswith(".h5"): - raise RuntimeError(f"Input file {input_file_path} does not look like a HDF5 file") - readout_mode = peek_readout_type(input_file_path) - if readout_mode is not HexagonalReadoutMode.CIRCULAR: - raise RuntimeError("Only CIRCULAR readout is supported.") - file_type = digi_input_file_class(readout_mode) - input_file = file_type(input_file_path) - header = input_file.header - input_file.close() - - enc = header["enc"] - # Scan different zero-suppression thresholds - n_files = 10 - zero_sup_threshold_grid = np.linspace(0, 5 * enc, n_files, dtype=int) - results = np.zeros((n_files, 3)) - default_folder = pathlib.Path.home() / "hexsampledata" - for i, zero_sup_threshold in enumerate(zero_sup_threshold_grid): - # Define the arguments for the reconstruction - recon_kwargs = dict( - input_file=input_file_path, - suffix=(f"recon_enc{enc}_thr{zero_sup_threshold}_" - f"{kwargs['pos_recon_algorithm']}_nn{kwargs['num_neighbors']}"), - zero_sup_threshold=zero_sup_threshold, - num_neighbors=kwargs["num_neighbors"], - pos_recon_algorithm=kwargs["pos_recon_algorithm"], - ) - # Check if reconstructed files already exist - file_name = default_folder / ( - f"{pathlib.Path(recon_kwargs['input_file']).stem}_" - f"{recon_kwargs['suffix']}.h5") - if not file_name.with_suffix(".h5").exists(): - output_file_path = reconstruct(**recon_kwargs) - else: - output_file_path = str(file_name.with_suffix(".h5")) - # Analyze the resolution of the reconstructed file - res_kwargs = dict( - input_file=output_file_path, - ) - results[i] = resolution(**res_kwargs) - - plt.figure("Position resolution vs zero-suppression threshold") - plt.plot(zero_sup_threshold_grid / enc, results[:, 0], ".k", label="Mean") - plt.plot(zero_sup_threshold_grid / enc, results[:, 1], ".r", label="Stddev") - # plt.plot(zero_sup_threshold_grid / enc, results[:, 2], ".b", label="FWHM") - plt.xlabel("Zero-suppression threshold / enc") - plt.ylabel("Position resolution [pitch]") - plt.legend() - - plt.show() - - -if __name__ == "__main__": - analyze_noise(**vars(ARGPARSER.parse_args())) - diff --git a/scripts/resolution.py b/scripts/resolution.py index 79932b6e..51bf186f 100644 --- a/scripts/resolution.py +++ b/scripts/resolution.py @@ -2,13 +2,13 @@ from pathlib import Path import numpy as np -from aptapy.hist import Histogram1d, Histogram2d +from aptapy.hist import Histogram2d from aptapy.plotting import plt from hexsample.fileio import ReconInputFile -from hexsample.hexagon import HexagonalGrid, HexagonalLayout +from hexsample.hexagon import HexagonalLayout from hexsample.pipeline import reconstruct, simulate -from hexsample.resolution import eef, hew +from hexsample.resolution import eef_size_scan __description__ = "" @@ -86,72 +86,36 @@ def resolution(**kwargs): num_rows=352, gain=1., ) - # Reconstruct the simulated file with different algorithms + 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( - input_file=str(simulation_path), - suffix=centroid_prefix, - zero_sup_threshold=zero_sup_threshold, - max_neighbors=6, - pos_recon_algorithm="centroid" - ) - # All pixels reconstructed with the best algorithm (eta for 2 and 3, centroid otherwise) + 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( - input_file=str(simulation_path), - suffix=best_prefix, - zero_sup_threshold=zero_sup_threshold, - max_neighbors=6, - pos_recon_algorithm="eta" - ) - # Study the EEF + 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) - xlabel = "Distance residual [pitch normalized]" - ylabel = "Encircled Energy Fraction" # Plot the centroid eef for all cluster sizes - centroid_recon_file = ReconInputFile(str(centroid_path)) - plt.figure(f"centroid_eef_{enc}enc_{zero_sup_threshold}zsup") - plt.plot(x, eef(x, centroid_recon_file, 0), - label=f"1 pix (HEW={hew(centroid_recon_file, 0):.2f})", - linestyle="-", color="black") - plt.plot(x, eef(x, centroid_recon_file, 1), - label=f"2 pix (HEW={hew(centroid_recon_file, 1):.2f})", - linestyle="--", color="black") - plt.plot(x, eef(x, centroid_recon_file, 2), - label=f"3 pix (HEW={hew(centroid_recon_file, 2):.2f})", - linestyle=":", color="black") - plt.plot(x, eef(x, centroid_recon_file, max_neighbors=6), - label=f"All pix (HEW={hew(centroid_recon_file, max_neighbors=6):.2f})", - linestyle="-.", color="black") - plt.xlim(x[0], x[-1]) - plt.ylim(0, 1) - plt.xlabel(xlabel) - plt.ylabel(ylabel) - plt.legend() + 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_recon_file = ReconInputFile(str(best_path)) - plt.figure(f"best_eef_{enc}enc_{zero_sup_threshold}zsup") - plt.plot(x, eef(x, best_recon_file, 0), - label=f"1 pix (HEW={hew(best_recon_file, 0):.2f})", - linestyle="-", color="black") - plt.plot(x, eef(x, best_recon_file, 1), - label=f"2 pix (HEW={hew(best_recon_file, 1):.2f})", - linestyle="--", color="black") - plt.plot(x, eef(x, best_recon_file, 2), - label=f"3 pix (HEW={hew(best_recon_file, 2):.2f})", - linestyle=":", color="black") - plt.plot(x, eef(x, best_recon_file, max_neighbors=6), - label=f"All pix (HEW={hew(best_recon_file, max_neighbors=6):.2f})", - linestyle="-.", color="black") - plt.xlim(x[0], x[-1]) - plt.ylim(0, 1) - plt.xlabel(xlabel) - plt.ylabel(ylabel) - plt.legend() + best_fig = plt.figure(f"best_eef_{enc}enc_{zero_sup_threshold}zsup") + eef_size_scan(x, best_recon_file) + + + + + centroid_recon_file.close() best_recon_file.close() @@ -214,8 +178,11 @@ def resolution(**kwargs): # plt.ylabel("Mean distance residual [pitch normalized]") # plt.legend() - # if kwargs["save"]: - # fig_format = "png" + if kwargs["save"]: + fig_format = "png" + centroid_fig.savefig(FIGURES_DIR / f"centroid_eef_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) + best_fig.savefig(FIGURES_DIR / f"best_eef_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) + # # fig # fig_1pix.savefig(FIGURES_DIR / f"1pix_resolution_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) # fig_2pix.savefig(FIGURES_DIR / f"2pix_resolution_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) # fig_3pix.savefig(FIGURES_DIR / f"3pix_resolution_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) diff --git a/src/hexsample/resolution.py b/src/hexsample/resolution.py index e8824da9..5233ab39 100644 --- a/src/hexsample/resolution.py +++ b/src/hexsample/resolution.py @@ -22,6 +22,7 @@ import numpy as np from aptapy.hist import Histogram1d +from aptapy.plotting import plt from .fileio import ReconInputFile from .hexagon import HexagonalGrid @@ -82,7 +83,7 @@ def dist_from_pixel_center(input_file: ReconInputFile) -> np.ndarray: return dr0 -def hist_distance_res(input_file: ReconInputFile, num_neighbors: int = 0, +def hist_distance_res(input_file: ReconInputFile, num_neighbors: int = 0, max_neighbors: int = -1) -> Histogram1d: """Create the histogram of distance residuals for the given numqber of neighbors. @@ -103,10 +104,7 @@ def hist_distance_res(input_file: ReconInputFile, num_neighbors: int = 0, """ # Select the cluster sizes we want and create the mask size = input_file.column("cluster_size") - if max_neighbors >= 0: - mask = size <= max_neighbors + 1 - else: - mask = size == num_neighbors + 1 + 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 @@ -164,3 +162,37 @@ def hew(input_file: ReconInputFile, num_neighbors: int = 0, """ hist = hist_distance_res(input_file, num_neighbors, max_neighbors) return hist.ppf(0.5) + + +def eef_size_scan(x: np.ndarray, input_file: ReconInputFile) -> None: + """Plot the EEF 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/p" + ylabel = "Encircled Energy Fraction" + # 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="black") + plt.plot(x, eef(x, input_file, 1), + label=f"2 pix (HEW={hew(input_file, 1):.2f})", + linestyle="-", color="black") + plt.plot(x, eef(x, input_file, 2), + label=f"3 pix (HEW={hew(input_file, 2):.2f})", + linestyle="--", color="black") + plt.plot(x, eef(x, input_file, max_neighbors=6), + label=f"All pix (HEW={hew(input_file, max_neighbors=6):.2f})", + linestyle="-.", color="black") + # Finalize the plot + plt.xlim(x[0], x[-1]) + plt.ylim(0, 1) + plt.xlabel(xlabel) + plt.ylabel(ylabel) + plt.legend() diff --git a/src/hexsample/tasks.py b/src/hexsample/tasks.py index aa6f5bd6..e814a683 100644 --- a/src/hexsample/tasks.py +++ b/src/hexsample/tasks.py @@ -385,58 +385,3 @@ def quicklook(input_file_path: str) -> None: plt.errorbar(xx, yy, fmt='.k', label='Mean') input_file.close() plt.show() - - -class ResolutionDefaults: - """Default parameters for the resolution 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. - """ - - -def resolution(input_file_path: str) -> None: - """Calculate the distribution of the distance between reconstructed and true position - and calculate its mean, standard deviation and FWHM. - - Arguments - --------- - input_file_path : str - The path to the input recon file. - - Returns - ------- - mean : float - The mean of the distance distribution. - stddev : float - The standard deviation of the distance distribution. - fwhm : float - The FWHM of the distance distribution. - """ - name, args = current_call() - logger.info(f"Running {__name__}.{name} with arguments {args}...") - input_file = ReconInputFile(input_file_path) - header = input_file.digi_header - # Get the reconstructed and true positions - x_mc = input_file.mc_column("absx") - y_mc = input_file.mc_column("absy") - x = input_file.column("posx") - y = input_file.column("posy") - # Calculate the distance between reconstructed and true position - dr = np.sqrt((x - x_mc)**2 + (y - y_mc)**2) / header["pitch"] - # Create the histogram of the distance - edges = np.linspace(0., 1., 101) - hist = Histogram1d(edges, xlabel=r"$|\vec{r} - \vec{r}_{mc}|$ / pitch") - hist.fill(dr) - mean, stddev = hist.binned_statistics() - # Calculate the FWHM if possible - try: - fwhm = hist.fwhm() - except ValueError: - fwhm = np.nan - plt.figure("Distance resolution") - hist.plot() - - input_file.close() - - return mean, stddev, fwhm From f152bdca0bc70ac0b5ac10115e1e937d2e79e7c6 Mon Sep 17 00:00:00 2001 From: Augusto Cattafesta Date: Thu, 5 Feb 2026 18:28:28 +0100 Subject: [PATCH 11/18] Spatial resolution dependence. --- scripts/resolution.py | 85 ++++++------------------------------- src/hexsample/resolution.py | 70 +++++++++++++++++++++++++----- 2 files changed, 74 insertions(+), 81 deletions(-) diff --git a/scripts/resolution.py b/scripts/resolution.py index 51bf186f..7eb2e5f9 100644 --- a/scripts/resolution.py +++ b/scripts/resolution.py @@ -8,7 +8,7 @@ from hexsample.fileio import ReconInputFile from hexsample.hexagon import HexagonalLayout from hexsample.pipeline import reconstruct, simulate -from hexsample.resolution import eef_size_scan +from hexsample.resolution import eef_size_scan, resolution_spatial_dependence __description__ = "" @@ -74,7 +74,7 @@ def resolution(**kwargs): # Run simulation if the output file does not already exist if not simulation_path.exists(): simulate( - num_events=10000, + num_events=100000, output_file=str(simulation_path), beam="hexagonal", enc=enc, @@ -112,82 +112,25 @@ def resolution(**kwargs): 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) + spatial_dependence_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() - # Now we study the dependence of the resolution on the recon distance from the pixel center - # For two pixel events eta and centroid - # hist_2d_2pix_eta = Histogram2d( - # dr_binning, dr_binning, - # xlabel="Reconstructed distance from pixel center [pitch normalized]", - # ylabel="Distance residual [pitch normalized]") - # hist_2d_2pix_eta.fill(dr0_2pix_eta, dr_2pix_eta) - # plt.figure("2pix_eta_residuals_vs_dr0") - # hist_2d_2pix_eta.plot() - # hist_2d_2pix_centroid = Histogram2d( - # dr_binning, dr_binning, - # xlabel="Reconstructed distance from pixel center [pitch normalized]", - # ylabel="Distance residual [pitch normalized]") - # hist_2d_2pix_centroid.fill(dr0_2pix_centroid, dr_2pix_centroid) - # plt.figure("2pix_centroid_residuals_vs_dr0") - # hist_2d_2pix_centroid.plot() - # # Take the slices and calculate the mean and stddev - # centers_2pix_eta, loc_2pix_eta, err_2pix_eta = _estimate_loc(hist_2d_2pix_eta) - # centers_2pix_centroid, loc_2pix_centroid, err_2pix_centroid = _estimate_loc(hist_2d_2pix_centroid) - # plt.figure("2pix_resolution_distance_dependency") - # plt.errorbar( - # centers_2pix_eta, loc_2pix_eta, yerr=err_2pix_eta, - # fmt=".r", label="2pix ETA") - # plt.errorbar( - # centers_2pix_centroid, loc_2pix_centroid, yerr=err_2pix_centroid, - # fmt=".b", label="2pix CENTROID") - # plt.xlabel("Reconstructed distance from pixel center [pitch normalized]") - # plt.ylabel("Mean distance residual [pitch normalized]") - # plt.legend() - - # # For three pixel events eta and centroid - # hist_2d_3pix_eta = Histogram2d( - # dr_binning, dr_binning, - # xlabel="Reconstructed distance from pixel center [pitch normalized]", - # ylabel="Distance residual [pitch normalized]") - # hist_2d_3pix_eta.fill(dr0_3pix_eta, dr_3pix_eta) - # plt.figure("3pix_eta_residuals_vs_dr0") - # hist_2d_3pix_eta.plot() - # hist_2d_3pix_centroid = Histogram2d( - # dr_binning, dr_binning, - # xlabel="Reconstructed distance from pixel center [pitch normalized]", - # ylabel="Distance residual [pitch normalized]") - # hist_2d_3pix_centroid.fill(dr0_3pix_centroid, dr_3pix_centroid) - # plt.figure("3pix_centroid_residuals_vs_dr0") - # hist_2d_3pix_centroid.plot() - # # Take the slices and calculate the mean and stddev - # centers_3pix_eta, loc_3pix_eta, err_3pix_eta = _estimate_loc(hist_2d_3pix_eta) - # centers_3pix_centroid, loc_3pix_centroid, err_3pix_centroid = _estimate_loc(hist_2d_3pix_centroid) - # plt.figure("3pix_resolution_distance_dependency") - # plt.errorbar( - # centers_3pix_eta, loc_3pix_eta, yerr=err_3pix_eta, - # fmt=".r", label="3pix ETA") - # plt.errorbar( - # centers_3pix_centroid, loc_3pix_centroid, yerr=err_3pix_centroid, - # fmt=".b", label="3pix CENTROID") - # plt.xlabel("Reconstructed distance from pixel center [pitch normalized]") - # plt.ylabel("Mean distance residual [pitch normalized]") - # plt.legend() - if kwargs["save"]: fig_format = "png" centroid_fig.savefig(FIGURES_DIR / f"centroid_eef_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) best_fig.savefig(FIGURES_DIR / f"best_eef_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) - # # fig - # fig_1pix.savefig(FIGURES_DIR / f"1pix_resolution_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) - # fig_2pix.savefig(FIGURES_DIR / f"2pix_resolution_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) - # fig_3pix.savefig(FIGURES_DIR / f"3pix_resolution_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) - # fig_allpix.savefig(FIGURES_DIR / f"allpix_resolution_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) - + spatial_dependence_fig.savefig(FIGURES_DIR / f"resolution_spatial_dependence_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) if __name__ == "__main__": resolution(**vars(HXETA_ARGPARSER.parse_args())) diff --git a/src/hexsample/resolution.py b/src/hexsample/resolution.py index 5233ab39..12557437 100644 --- a/src/hexsample/resolution.py +++ b/src/hexsample/resolution.py @@ -20,8 +20,9 @@ """Resolution analysis facilities. """ +from typing import Tuple import numpy as np -from aptapy.hist import Histogram1d +from aptapy.hist import Histogram1d, Histogram2d from aptapy.plotting import plt from .fileio import ReconInputFile @@ -83,7 +84,7 @@ def dist_from_pixel_center(input_file: ReconInputFile) -> np.ndarray: return dr0 -def hist_distance_res(input_file: ReconInputFile, num_neighbors: int = 0, +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 numqber of neighbors. @@ -137,7 +138,7 @@ def eef(x: np.ndarray, input_file: ReconInputFile, num_neighbors: int = 0, eef : np.ndarray The Encircled Energy Function evaluated at the given normalized distance values. """ - hist = hist_distance_res(input_file, num_neighbors, max_neighbors) + hist = hist_distance_residuals(input_file, num_neighbors, max_neighbors) return hist.cdf(x) @@ -160,7 +161,7 @@ def hew(input_file: ReconInputFile, num_neighbors: int = 0, hew : float The Half Energy Width (HEW) in pitch normalized units. """ - hist = hist_distance_res(input_file, num_neighbors, max_neighbors) + hist = hist_distance_residuals(input_file, num_neighbors, max_neighbors) return hist.ppf(0.5) @@ -175,24 +176,73 @@ def eef_size_scan(x: np.ndarray, input_file: ReconInputFile) -> None: input_file : ReconInputFile The input file to analyze. """ - xlabel = "r/p" + 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="black") + linestyle=":", color=color) plt.plot(x, eef(x, input_file, 1), label=f"2 pix (HEW={hew(input_file, 1):.2f})", - linestyle="-", color="black") + linestyle="-", color=color) plt.plot(x, eef(x, input_file, 2), label=f"3 pix (HEW={hew(input_file, 2):.2f})", - linestyle="--", color="black") + linestyle="--", color=color) plt.plot(x, eef(x, input_file, max_neighbors=6), - label=f"All pix (HEW={hew(input_file, max_neighbors=6):.2f})", - linestyle="-.", color="black") + 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) + 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 From 64a2c4895ebc30987946498df008171294875387 Mon Sep 17 00:00:00 2001 From: Augusto Cattafesta Date: Fri, 6 Feb 2026 09:23:46 +0100 Subject: [PATCH 12/18] Minor. --- src/hexsample/resolution.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/hexsample/resolution.py b/src/hexsample/resolution.py index 12557437..b6e18e6d 100644 --- a/src/hexsample/resolution.py +++ b/src/hexsample/resolution.py @@ -166,8 +166,8 @@ def hew(input_file: ReconInputFile, num_neighbors: int = 0, def eef_size_scan(x: np.ndarray, input_file: ReconInputFile) -> None: - """Plot the EEF for one, two, three pixel events and for all the reconstructed events on - the same figure. + """Plot the Encircled Energy Function for one, two, three pixel events and for all the + reconstructed events on the same figure. Arguments --------- From b52fcc2a874ac1451dc94bde1a7ada7ba31471b6 Mon Sep 17 00:00:00 2001 From: Augusto Cattafesta Date: Fri, 6 Feb 2026 11:17:41 +0100 Subject: [PATCH 13/18] Study resolution with different sup_threshold. --- scripts/resolution.py | 44 ++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 43 insertions(+), 1 deletion(-) diff --git a/scripts/resolution.py b/scripts/resolution.py index 7eb2e5f9..78ed886c 100644 --- a/scripts/resolution.py +++ b/scripts/resolution.py @@ -8,7 +8,7 @@ from hexsample.fileio import ReconInputFile from hexsample.hexagon import HexagonalLayout from hexsample.pipeline import reconstruct, simulate -from hexsample.resolution import eef_size_scan, resolution_spatial_dependence +from hexsample.resolution import eef, eef_size_scan, resolution_spatial_dependence __description__ = "" @@ -126,11 +126,53 @@ def resolution(**kwargs): 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 i, zero_sup_ratio in enumerate(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() + + eef_zsup_best_fig = plt.figure(f"eef_vs_zsup_best_enc{enc}") + print("Reconstructing files for eta algorithm...") + for i, zero_sup_ratio in enumerate(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() + if kwargs["save"]: fig_format = "png" centroid_fig.savefig(FIGURES_DIR / f"centroid_eef_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) best_fig.savefig(FIGURES_DIR / f"best_eef_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) spatial_dependence_fig.savefig(FIGURES_DIR / f"resolution_spatial_dependence_{enc}enc_{zero_sup_threshold}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())) From 713e3d4db951ad5f5f6311366eb17857c6d46e35 Mon Sep 17 00:00:00 2001 From: Augusto Cattafesta Date: Fri, 6 Feb 2026 12:04:35 +0100 Subject: [PATCH 14/18] Minor. --- scripts/resolution.py | 53 +++++++++++-------------------------- src/hexsample/resolution.py | 2 +- 2 files changed, 16 insertions(+), 39 deletions(-) diff --git a/scripts/resolution.py b/scripts/resolution.py index 78ed886c..74a36be8 100644 --- a/scripts/resolution.py +++ b/scripts/resolution.py @@ -2,7 +2,6 @@ from pathlib import Path import numpy as np -from aptapy.hist import Histogram2d from aptapy.plotting import plt from hexsample.fileio import ReconInputFile @@ -27,46 +26,19 @@ FIGURES_DIR.mkdir(parents=True, exist_ok=True) -def _estimate_loc(hist: Histogram2d) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """Estimate the location of the distribution in each slice of the histogram. +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 diffrent + cluster sizes and reconstruction algorithm for a given ENC and zero suppression threshold. - Arguments - --------- - hist : Histogram2d - The 2D histogram to analyze. + 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. - Returns - ------- - centers : np.ndarray - The bin centers of the first axis. - loc : np.ndarray - The estimated location of the distribution in each slice. - err : np.ndarray - The estimated uncertainty on the location in each slice. + Finally, for the given ENC, the EEF is plotted for different zero suppression thresholds and + for both the centroid and eta algorithms. """ - centers = hist.bin_centers(axis=0) - loc = np.zeros(centers.shape) - err = np.zeros(centers.shape) - for i in range(len(centers)): - slice_ = hist.slice1d(i) - n = slice_.content.sum() - if n == 0: - loc[i] = np.nan - err[i] = np.nan - continue - # Estimate the location as the median of the distribution and its error - x = np.repeat(slice_.bin_centers(), slice_.content.astype(int)) - loc[i] = np.mean(x) - err[i] = np.std(x) / np.sqrt(n) # Approximate std error of median - # Remove NaN entries - mask = ~np.isnan(loc) - centers = centers[mask] - loc = loc[mask] - err = err[mask] - return centers, loc, err - - -def resolution(**kwargs): enc = kwargs["enc"] zero_sup_threshold = kwargs["zero_sup_threshold"] file_prefix = f"simulation_resolution_{enc}enc_hexagonal" @@ -143,6 +115,11 @@ def resolution(**kwargs): 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...") diff --git a/src/hexsample/resolution.py b/src/hexsample/resolution.py index b6e18e6d..5d5e2eef 100644 --- a/src/hexsample/resolution.py +++ b/src/hexsample/resolution.py @@ -85,7 +85,7 @@ def dist_from_pixel_center(input_file: ReconInputFile) -> np.ndarray: def hist_distance_residuals(input_file: ReconInputFile, num_neighbors: int = 0, - max_neighbors: int = -1) -> Histogram1d: + max_neighbors: int = -1) -> Histogram1d: """Create the histogram of distance residuals for the given numqber of neighbors. Arguments From 0134b3ed57fd070955aa3fe95b9451d4d2e8e8bf Mon Sep 17 00:00:00 2001 From: Augusto Cattafesta Date: Fri, 6 Feb 2026 12:08:08 +0100 Subject: [PATCH 15/18] Minor. --- src/hexsample/pipeline.py | 7 ------- src/hexsample/tasks.py | 36 +----------------------------------- 2 files changed, 1 insertion(+), 42 deletions(-) diff --git a/src/hexsample/pipeline.py b/src/hexsample/pipeline.py index f062249b..f7d6f48f 100644 --- a/src/hexsample/pipeline.py +++ b/src/hexsample/pipeline.py @@ -71,10 +71,3 @@ def quicklook(**kwargs) -> None: """ input_file_path = kwargs["input_file"] return tasks.quicklook(input_file_path) - - -def resolution(**kwargs) -> None: - """Analyze the resolution from a recon file. - """ - input_file_path = kwargs["input_file"] - return tasks.resolution(input_file_path) diff --git a/src/hexsample/tasks.py b/src/hexsample/tasks.py index e814a683..9c37af77 100644 --- a/src/hexsample/tasks.py +++ b/src/hexsample/tasks.py @@ -326,13 +326,6 @@ def quicklook(input_file_path: str) -> None: plt.xlabel("Energy [eV]") plt.legend() - # size_hist = create_histogram(input_file, "cluster_size", mc=False) - # plt.figure("Cluster size distribution") - # size_hist.plot() - # plt.xlabel("Cluster size") - # plt.ylabel("Counts") - - # Plotting the reconstructed x and y position and the true position. plt.figure("Reconstructed photons position") binning = np.linspace(-5. * 0.2, 5. * 0.2, 100) @@ -356,32 +349,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() - plt.figure("distance resolution") - dr = np.sqrt((x - x_mc)**2 + (y - y_mc)**2) / 0.005 - binning = np.linspace(dr.min(), dr.max(), 100) - histdr = Histogram1d(binning, xlabel=r"$\sqrt{(x - x_{MC})^2 + (y - y_{MC})^2}$ [cm]") - histdr.fill(dr) - print(f"dr histogram FWHM: {histdr.fwhm()}") - - histdr.plot() - - from .hexagon import HexagonalGrid - grid = HexagonalGrid() - x0, y0 = grid.pixel_to_world(*grid.world_to_pixel(x_mc, y_mc)) - dr_abs = np.sqrt((x - x0)**2 + (y - y0)**2) / grid.pitch - bins = np.linspace(0, 0.7, 200) - hist = Histogram2d(bins, bins) - # I need the recon distance from the hit pixel center - hist.fill(dr_abs, dr) - plt.figure() - hist.plot() - hist_mean, hist_sigma = hist.project_statistics() - plt.figure("Reconstructed vs true distance from pixel center") - yy = hist_mean.content.flatten() - mask = yy > 0 - yy = yy[mask] - xx = hist_mean.bin_centers()[mask] - yy_err = hist_mean.errors.flatten()[mask] - plt.errorbar(xx, yy, fmt='.k', label='Mean') input_file.close() - plt.show() + plt.show() \ No newline at end of file From a9e2c9a638fbb03a1c0a2a355e6eaa740d0f9191 Mon Sep 17 00:00:00 2001 From: Augusto Cattafesta Date: Fri, 6 Feb 2026 14:02:17 +0100 Subject: [PATCH 16/18] Linting. --- scripts/eta.py | 19 +--------------- scripts/resolution.py | 32 +++++++++++++++++--------- src/hexsample/clustering.py | 45 ++++++++++++++++++++----------------- src/hexsample/resolution.py | 1 + src/hexsample/tasks.py | 2 +- 5 files changed, 49 insertions(+), 50 deletions(-) diff --git a/scripts/eta.py b/scripts/eta.py index 6a89ba2b..4ee914cf 100644 --- a/scripts/eta.py +++ b/scripts/eta.py @@ -4,7 +4,7 @@ from pathlib import Path import numpy as np -from aptapy.hist import Histogram1d, Histogram2d +from aptapy.hist import Histogram2d from aptapy.modeling import AbstractFitModel from aptapy.models import Probit from aptapy.plotting import last_line_color, plt @@ -159,23 +159,6 @@ def calibrate_2pix(eta: np.ndarray, photon_pos: np.ndarray, versors: np.ndarray, if kwargs.get("save", False): fig.savefig(RESULTS_DIR / "2pix_cal.pdf", format="pdf") - eta_min = min(eta)[0] - eta_hist_binning = np.linspace(eta_min, 0.5, 101) - eta_hist = Histogram1d(eta_hist_binning, xlabel="eta", ylabel="Counts") - eta_hist.fill(eta) - plt.figure("eta_2pix_distribution") - eta_hist.plot() - - from scipy.special import ndtri - - sigma_x = model.sigma.value * np.sqrt(2 * np.pi) * np.exp(0.5 * abs(ndtri(eta_binning))**2) * np.sqrt(eta_binning**2 + (1 - eta_binning)**2) * 30 / 1600 - plt.figure("dr_vs_eta_2pix_derivative") - plt.plot(model(eta_binning), sigma_x, label="d(dr/p)/d(eta)") - - plt.figure("dr_distr") - plt.hist(model(eta.astype(float)), bins=100, label="dr / p") - plt.legend() - return model diff --git a/scripts/resolution.py b/scripts/resolution.py index 74a36be8..bffa395e 100644 --- a/scripts/resolution.py +++ b/scripts/resolution.py @@ -13,9 +13,12 @@ # 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") +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(): @@ -85,7 +88,7 @@ def resolution(**kwargs): eef_size_scan(x, best_recon_file) # Plot the spatial dependence of the resolution for both algorithms (all events) - spatial_dependence_fig = plt.figure("resolution_spatial_dependence") + 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), @@ -104,7 +107,7 @@ def resolution(**kwargs): max_neighbors=6) eef_zsup_centroid_fig = plt.figure(f"eef_vs_zsup_centroid_enc{enc}") print("Reconstructing files for centroid algorithm...") - for i, zero_sup_ratio in enumerate(zero_sup_ratios): + 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" @@ -123,7 +126,7 @@ def resolution(**kwargs): eef_zsup_best_fig = plt.figure(f"eef_vs_zsup_best_enc{enc}") print("Reconstructing files for eta algorithm...") - for i, zero_sup_ratio in enumerate(zero_sup_ratios): + 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" @@ -141,13 +144,20 @@ def resolution(**kwargs): plt.ylim(0, 1) plt.legend() + # Save figures, if requested if kwargs["save"]: fig_format = "png" - centroid_fig.savefig(FIGURES_DIR / f"centroid_eef_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) - best_fig.savefig(FIGURES_DIR / f"best_eef_{enc}enc_{zero_sup_threshold}zsup.{fig_format}", format=fig_format) - spatial_dependence_fig.savefig(FIGURES_DIR / f"resolution_spatial_dependence_{enc}enc_{zero_sup_threshold}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) + 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) diff --git a/src/hexsample/clustering.py b/src/hexsample/clustering.py index 526618e2..2474cab2 100644 --- a/src/hexsample/clustering.py +++ b/src/hexsample/clustering.py @@ -128,29 +128,34 @@ 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 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 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) + theta = Probit().evaluate((eta_diff + 1)/2, 0, eta_3pix_theta0) / r + # Reconstructing the position using r and theta + 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: - 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. - 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. - 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) - theta = Probit().evaluate((eta_diff + 1)/2, 0, eta_3pix_theta0) / r - # Reconstructing the position using r and theta - 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]) - return x_recon, y_recon + # 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 @dataclass diff --git a/src/hexsample/resolution.py b/src/hexsample/resolution.py index 5d5e2eef..dcc1b9fa 100644 --- a/src/hexsample/resolution.py +++ b/src/hexsample/resolution.py @@ -21,6 +21,7 @@ """ from typing import Tuple + import numpy as np from aptapy.hist import Histogram1d, Histogram2d from aptapy.plotting import plt diff --git a/src/hexsample/tasks.py b/src/hexsample/tasks.py index 9c37af77..3a0a03f3 100644 --- a/src/hexsample/tasks.py +++ b/src/hexsample/tasks.py @@ -350,4 +350,4 @@ def quicklook(input_file_path: str) -> None: histy = Histogram1d(binning, xlabel=r"$y - y_{MC}$ [cm]").fill(y-y_mc) histy.plot() input_file.close() - plt.show() \ No newline at end of file + plt.show() From b10c8e26d4ebcbdf0832e56be84a2a34349ff1f5 Mon Sep 17 00:00:00 2001 From: Augusto Cattafesta Date: Fri, 6 Feb 2026 14:27:05 +0100 Subject: [PATCH 17/18] Minor. --- scripts/eta.py | 2 +- scripts/resolution.py | 4 ++-- src/hexsample/resolution.py | 9 ++++++--- src/hexsample/tasks.py | 15 +++++++++++---- 4 files changed, 20 insertions(+), 10 deletions(-) diff --git a/scripts/eta.py b/scripts/eta.py index 4ee914cf..e013c424 100644 --- a/scripts/eta.py +++ b/scripts/eta.py @@ -24,7 +24,7 @@ # 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, +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") diff --git a/scripts/resolution.py b/scripts/resolution.py index bffa395e..4c848048 100644 --- a/scripts/resolution.py +++ b/scripts/resolution.py @@ -33,8 +33,8 @@ 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 diffrent - cluster sizes and reconstruction algorithm for a given ENC and zero suppression threshold. + 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. diff --git a/src/hexsample/resolution.py b/src/hexsample/resolution.py index dcc1b9fa..b5d15610 100644 --- a/src/hexsample/resolution.py +++ b/src/hexsample/resolution.py @@ -80,14 +80,14 @@ def dist_from_pixel_center(input_file: ReconInputFile) -> np.ndarray: # Access the reconstructed positions x = input_file.column("posx") y = input_file.column("posy") - # Calulate the reconstructed distance from the pixel center + # 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 numqber of neighbors. + """Create the histogram of distance residuals for the given number of neighbors. Arguments --------- @@ -111,7 +111,7 @@ def hist_distance_residuals(input_file: ReconInputFile, num_neighbors: int = 0, # 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 - # it covers the full range of dr. + # spans all the pitch. dr_binning = np.linspace(0., 1., 101) hist = Histogram1d(dr_binning) hist.fill(dr[mask]) @@ -234,6 +234,9 @@ def resolution_spatial_dependence(input_file: ReconInputFile, num_neighbors: int 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) diff --git a/src/hexsample/tasks.py b/src/hexsample/tasks.py index 3a0a03f3..e8703c40 100644 --- a/src/hexsample/tasks.py +++ b/src/hexsample/tasks.py @@ -204,6 +204,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. @@ -236,16 +240,19 @@ 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. If max_neighbors is -1, take num_neighbors - # only. Otherwise take all sizes from 1 to max_neighbors + 1. When specified, max_neighbors - # gets priority over num_neighbors. + # 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) From e641bab07f505b1cbf849c0d8ed75cc0f3cf0ac1 Mon Sep 17 00:00:00 2001 From: Luca Baldini Date: Mon, 9 Feb 2026 10:36:03 +0100 Subject: [PATCH 18/18] Updated release notes. --- docs/release_notes.rst | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/docs/release_notes.rst b/docs/release_notes.rst index fa736772..4423fabc 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)