Fixed the X/Y convention and clipping bugs

Co-authored-by: Copilot <copilot@github.com>
This commit is contained in:
2026-06-10 16:11:44 +02:00
co-authored by Copilot
parent 8716f175dd
commit 73b709a563
+45 -13
View File
@@ -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: