From 73b709a56310ce743372945f5c84be5d70d98d5d Mon Sep 17 00:00:00 2001 From: "xiangyu.xie" Date: Wed, 10 Jun 2026 16:11:44 +0200 Subject: [PATCH] Fixed the X/Y convention and clipping bugs Co-authored-by: Copilot --- etaInterpolationFromClusters.py | 58 +++++++++++++++++++++++++-------- 1 file changed, 45 insertions(+), 13 deletions(-) diff --git a/etaInterpolationFromClusters.py b/etaInterpolationFromClusters.py index ed751dc..71154a6 100644 --- a/etaInterpolationFromClusters.py +++ b/etaInterpolationFromClusters.py @@ -1,4 +1,7 @@ from EtaInterpolationFunctions import * +import matplotlib +matplotlib.use('Agg') + import numpy as np from array import array import matplotlib.pyplot as plt @@ -7,7 +10,11 @@ from multiprocessing import Pool from argparse import ArgumentParser parser = ArgumentParser() -parser.add_argument('--runname', '-r', type=str, default='2603MaxIV_Edge3Filters_pos0_12keV', choices=['2603MaxIV_Edge3Filters_pos0_12keV', '2603MaxIV_FlatField3Filters_pos0_12keV'], help='Name of the run, used to locate the cluster files and save the results') +parser.add_argument('--runname', '-r', type=str, default='2603MaxIV_Edge3Filters_pos0_12keV', choices=[ + '2603MaxIV_Edge3Filters_pos0_12keV', '2603MaxIV_FlatField3Filters_pos0_12keV', + '2603MaxIV_Edge2Filters_12keV', '2603MaxIV_Flat2Filters_12keV', + '2603MaxIV_Edge1Filters_12keV', '2603MaxIV_Flat1Filters_12keV' + ], help='Name of the run, used to locate the cluster files and save the results') args = parser.parse_args() RUNNAME = args.runname @@ -19,16 +26,21 @@ NEtaBins = 201 ### generate eta map def get_eta_hist2D(clusterFile): with h5py.File(clusterFile, 'r') as f: - _eta_hist2D = np.zeros((NEtaBins, NEtaBins)) ### shape (201, 201), Y, X + _eta_hist2D = np.zeros((NEtaBins, NEtaBins)) ### shape (201, 201), X, Y, according build_xy_lut_Rosenblatt's convention clusters = f['clusters'][:] ### shape (N, 3, 3), N, Y, X sum_clusters = np.sum(clusters, axis=(1, 2)) ### shape (N,) - position_weights = np.array([-1, 0, 1]) - etaX = ((np.sum(clusters, axis=1) * position_weights).sum(axis=-1) / sum_clusters + 0.5) % 1 - etaY = ((np.sum(clusters, axis=2) * position_weights).sum(axis=-1) / sum_clusters + 0.5) % 1 + position_weights = np.array([-1, 0, 1]) + 0.5 - _eta_hist2D += np.histogram2d(etaY, etaX, bins=NEtaBins, range=[[0, 1], [0, 1]])[0] + etaX = (np.sum(clusters, axis=1) * position_weights).sum(axis=-1) / sum_clusters + etaY = (np.sum(clusters, axis=2) * position_weights).sum(axis=-1) / sum_clusters + + etaX = np.clip(etaX, 0, 1) + etaY = np.clip(etaY, 0, 1) + + _eta_hist2D += np.histogram2d(etaX, etaY, bins=NEtaBins, range=[[0, 1], [0, 1]])[0] ### X, Y convention according to build_xy_lut_Rosenblatt's convention return _eta_hist2D +print('Building eta histogram from clusters...') eta_hist2D = np.zeros((NEtaBins, NEtaBins)) with Pool(16) as pool: results = pool.map(get_eta_hist2D, clusterFiles) @@ -50,18 +62,38 @@ def reconstruct_position(clusterFile): _subpixel_positions = np.zeros((interpolationBins, interpolationBins)) with h5py.File(clusterFile, 'r') as f: clusters = f['clusters'][:] ### shape (N, 3, 3), N, Y, X - referencePoints = f['referencePoint'][:] ### shape (N, 2), N, (y, x) + referencePoints = f['referencePoint'][:] ### x, y + sum_clusters = np.sum(clusters, axis=(1, 2)) ### shape (N,) position_weights = np.array([-1, 0, 1]) + 0.5 - etaX = (np.sum(clusters, axis=1) * position_weights).sum(axis=-1) / sum_clusters - etaY = (np.sum(clusters, axis=2) * position_weights).sum(axis=-1) / sum_clusters + + etaX = (np.sum(clusters, axis=1) * position_weights).sum(axis=-1) / sum_clusters + etaY = (np.sum(clusters, axis=2) * position_weights).sum(axis=-1) / sum_clusters + + etaX = np.clip(etaX, 0, 1) + etaY = np.clip(etaY, 0, 1) + x_frac, y_frac = bilinear_xy_lookup(etaX, etaY, U_tab, V_tab) - x = referencePoints[:, 0] + x_frac + clusters.shape[1]//2 - y = referencePoints[:, 1] + y_frac + clusters.shape[2]//2 - _image += np.histogram2d(y, x, bins=(Roi[1]-Roi[0])*interpolationBins, range=[[Roi[0], Roi[1]], [Roi[2], Roi[3]]])[0] - _subpixel_positions += np.histogram2d(x%1, y%1, bins=interpolationBins, range=[[0, 1], [0, 1]])[0] + + x_ref = referencePoints[:, 0] + y_ref = referencePoints[:, 1] + + x = x_ref + x_frac + clusters.shape[1]//2 + y = y_ref + y_frac + clusters.shape[2]//2 + + _image += np.histogram2d( + y, x, + bins=(Roi[1]-Roi[0])*interpolationBins, + range=[[Roi[0], Roi[1]], [Roi[2], Roi[3]]] + )[0] + _subpixel_positions += np.histogram2d( + y%1, x%1, + bins=interpolationBins, + range=[[0, 1], [0, 1]] + )[0] return (_image, _subpixel_positions) +print('Reconstructing positions from clusters and accumulating hits...') reconstructed_image = np.zeros(((Roi[1]-Roi[0])*interpolationBins, (Roi[3]-Roi[2])*interpolationBins)) subpixel_positions = np.zeros((interpolationBins, interpolationBins)) with Pool(16) as pool: