diff --git a/src/hexsample/clustering.py b/src/hexsample/clustering.py index 1ee896f..b01445f 100644 --- a/src/hexsample/clustering.py +++ b/src/hexsample/clustering.py @@ -52,6 +52,7 @@ class Cluster: adc_to_ev: float pos_recon_algorithm: str recon_pars: Optional[dict] = None + _pulse_height = None def __post_init__(self) -> None: """Small cross check on the dimensions of the arrays passed in the constructor. @@ -67,7 +68,9 @@ def size(self) -> int: def pulse_height(self) -> float: """Return the total pulse height of the cluster. """ - return self.pha.sum() + if self._pulse_height is None: + self._pulse_height = self.pha.sum() + return self._pulse_height def energy(self) -> float: """Return the energy of the cluster in eV. @@ -146,6 +149,7 @@ def mle(self, m = mle(self.pha, equal_noise, position_cal.values, position_cal.bin_size, position_cal.xlims, position_cal.ylims, p0=p0) # Calculate the absolute position of the photon from the fit results. + self._pulse_height = m.values["q"] return self.x[0] + m.values["x"] * pitch, self.y[0] + m.values["y"] * pitch def position(self) -> Tuple[float, float]: @@ -297,12 +301,17 @@ def run(self, event) -> Optional[Cluster]: # and applying pedestal and gain correction. pha = (event.pha[adc_channel_order] - pedestal(col, row)) / gain(col, row) elif isinstance(event, DigiEventRectangular): + rows, cols = np.mgrid[event.roi.min_row:event.roi.max_row + 1, + event.roi.min_col:event.roi.max_col + 1] + cols = cols.flatten() + rows = rows.flatten() + event.pha = event.pha - pedestal(cols.ravel(), rows.ravel()).reshape(event.roi.shape()) seed_coords = event.highest_pixel() if self.readout.is_at_border(*seed_coords): return None neigh_coords = self.readout.neighbors(*seed_coords) col, row = np.vstack((seed_coords, neigh_coords)).T - pha = (event(col, row) - pedestal(col, row)) / gain(col, row) + pha = event(col, row) / gain(col, row) else: raise RuntimeError(f"Unsupported event type {type(event)} for clustering") # Zero suppressing the event (whatever the readout type)... diff --git a/src/hexsample/fileio.py b/src/hexsample/fileio.py index dd01322..292cc89 100644 --- a/src/hexsample/fileio.py +++ b/src/hexsample/fileio.py @@ -245,9 +245,9 @@ def _fill_recon_row(row: tables.tableextension.Row, event: ReconEvent) -> None: row["livetime"] = event.livetime #row["roi_size"] = event.roi_size row["cluster_size"] = event.cluster.size() + row["posx"], row["posy"] = event.position() row["adc"] = event.adc() row["energy"] = event.energy() - row["posx"], row["posy"] = event.position() row.append() diff --git a/src/hexsample/likelihood.py b/src/hexsample/likelihood.py index fd70d29..524cab0 100644 --- a/src/hexsample/likelihood.py +++ b/src/hexsample/likelihood.py @@ -187,50 +187,104 @@ def weighted_pha(pha: np.ndarray, f_interp: np.ndarray, inv_sigma2: np.ndarray) @njit -def nll_numba(x: float, y: float, pha: np.ndarray, f: np.ndarray, xbin0: float, ybin0: float, - bin_size: float, noise: np.ndarray) -> float: +def sherman_morrison_inverse(res: np.ndarray, f_interp: np.ndarray, total_pha: float, + noise: np.ndarray): + """Compute the chi2 (res^T * V^-1 * res) and log-determinant of the covariance + matrix V using the Sherman-Morrison formula and the matrix determinant lemma. + + This avoids the explicit inversion of the 7x7 matrix and allows to save space and time. + """ + # First we compute the diagonal elements of D, given by the square sum of electronic + # noise and statistical fluctuations, and their inverses. + d_matrix = np.zeros(7) + inv_d_matrix = np.zeros(7) + for i in range(7): + d_matrix[i] = (noise[i] ** 2) + total_pha * f_interp[i] + inv_d_matrix[i] = 1.0 / d_matrix[i] + # Then we calculate the term at the denominator. + sum_fk2_d = 0.0 + for k in range(7): + sum_fk2_d += (f_interp[k] ** 2) * inv_d_matrix[k] + # The final denominator term. + g = 1.0 / (1.0 - total_pha * sum_fk2_d + 1e-12) + # Calculate the log determinant of V using the matrix determinant lemma. + prod_d = 0.0 + for k in range(7): + prod_d += math.log(d_matrix[k]) + log_det_v = prod_d + math.log(abs(1.0 - total_pha * sum_fk2_d) + 1e-12) + # Calculate the quadratic form res^T * V^-1 * res using Sherman-Morrison algebra. + diagional_term = 0.0 + cross_term = 0.0 + for i in range(7): + diagional_term += (res[i] ** 2) * inv_d_matrix[i] + cross_term += res[i] * f_interp[i] * inv_d_matrix[i] + # Final chi2 value. + chi2 = diagional_term + g * total_pha * (cross_term ** 2) + # Calculate the vector w = V^-1 * res for the gradient calculation. + alpha = g * total_pha * cross_term + vettore_w = np.zeros(7) + for i in range(7): + vettore_w[i] = inv_d_matrix[i] * (res[i] + alpha * f_interp[i]) + dnll_df = np.zeros(7) + for i in range(7): + dnll_df[i] = total_pha * ( + 0.5 * inv_d_matrix[i] + - vettore_w[i] + - 0.5 * (vettore_w[i] ** 2) + + g * vettore_w[i] * f_interp[i] * (vettore_w[i] - inv_d_matrix[i] * res[i]) + ) + return chi2, log_det_v, dnll_df + + +@njit +def nll_numba(x: float, y: float, total_pha: float, pha: np.ndarray, f: np.ndarray, + xbin0: float, ybin0: float, bin_size: float, noise: np.ndarray) -> float: """Compute the negative log-likelihood for a given position (x, y). The model is based on the Gaussian diffusion of the charge cloud, and uses the precomputed - charge fractions in each pixel from the f map. The summed pha is profiled out to reduce the - dimensionality of the optimization. + charge fractions in each pixel from the f map. """ # Calculate the bin indices and fractional coordinates for the interpolation ix0, iy0, wx, wy = coordinates(x, y, xbin0, ybin0, bin_size, f.shape[1:]) # Interpolate the charge fractions for the 7 pixels in the cluster f_interp = interpolation(f, ix0, iy0, wx, wy) - # Calculate the inverse of the noise variance for each pixel - inv_sigma2 = 1.0 / (noise**2) - # Profile out the summed pha by finding the value that minimizes the NLL for fixed (x, y) - total_pha = weighted_pha(pha, f_interp, inv_sigma2) - # Now compute the NLL using the optimal energy - nll = 0.0 + # We need to invert the covariance matrix V. We use the Sherman-Morrison formula + # to avoid inverting a full 7x7 matrix. This is possible because V is a diagonal + # matrix D plus a rank 1 matrix. + res = np.zeros(7) for i in range(7): - mu = f_interp[i] * total_pha - res = pha[i] - mu - nll += 0.5 * (res**2 * inv_sigma2[i] + LOG2PI) + res[i] = pha[i] - (total_pha * f_interp[i]) + chi2, log_det_v, _ = sherman_morrison_inverse(res, f_interp, total_pha, noise) + # Finally we compute the negative log-likelihood using the chi2 and the log-determinant of V. + nll = 0.5 * chi2 + 0.5 * log_det_v return nll @njit -def nll_grad_numba(x: float, y: float, pha: np.ndarray, f: np.ndarray, xbin0: float, ybin0: float, - bin_size: float, noise: np.ndarray) -> np.ndarray: - """Compute the gradient of the negative log-likelihood with respect to the free parameters. +def nll_grad_numba(x: float, y: float, total_pha: float, pha: np.ndarray, f: np.ndarray, + xbin0: float, ybin0: float, bin_size: float, noise: np.ndarray) -> np.ndarray: + """Compute the complete analytical gradient of the exact multivariate NLL with + respect to x, y, and Q simultaneously (Returns array of size 3). """ # Calculate the bin indices and fractional coordinates for the interpolation ix0, iy0, wx, wy = coordinates(x, y, xbin0, ybin0, bin_size, f.shape[1:]) # Interpolate the charge fractions and their derivatives for the 7 pixels in the cluster f_interp, df_dx, df_dy = interpolation_derivatives(f, ix0, iy0, wx, wy, bin_size) - # Calculate the inverse of the noise variance for each pixel - inv_sigma2 = 1.0 / (noise**2) - # Profile out the summed pha by finding the value that minimizes the NLL for fixed (x, y) - total_pha = weighted_pha(pha, f_interp, inv_sigma2) - # Now compute the gradient using the optimal energy + # Calculate the residuals for the 7 pixels in the cluster. + res = np.zeros(7) + for i in range(7): + res[i] = pha[i] - (total_pha * f_interp[i]) + # We need to invert the covariance matrix V. We use the Sherman-Morrison formula + # to avoid inverting a full 7x7 matrix. + _, _, dnll_df = sherman_morrison_inverse(res, f_interp, total_pha, noise) + # Finally we compute the gradient of the negative log-likelihood with respect to x, y, and Q. gnll_x = 0.0 gnll_y = 0.0 for i in range(7): - mu = f_interp[i] * total_pha - d_loss_dmu = -(pha[i] - mu) * inv_sigma2[i] - gnll_x += d_loss_dmu * total_pha * df_dx[i] - gnll_y += d_loss_dmu * total_pha * df_dy[i] - return np.array([gnll_x, gnll_y]) + gnll_x += dnll_df[i] * df_dx[i] + gnll_y += dnll_df[i] * df_dy[i] + gnll_q = 0.0 + for i in range(7): + gnll_q -= (res[i] * f_interp[i]) / (noise[i] ** 2) + # Return the gradient as a numpy array of size 3. + return np.array([gnll_x, gnll_y, gnll_q]) diff --git a/src/hexsample/position.py b/src/hexsample/position.py index 370dbb0..4a7a13a 100644 --- a/src/hexsample/position.py +++ b/src/hexsample/position.py @@ -329,24 +329,27 @@ def mle( ymin, ymax = ylims # Define the objective functions for the optimization, which are the # negative log-likelihood... - def nll(x: float, y: float) -> float: - return nll_numba(x, y, pha, f, xmin, ymin, bin_size, noise) + def nll(x: float, y: float, q: float) -> float: + return nll_numba(x, y, q, pha, f, xmin, ymin, bin_size, noise) # ... and its gradient. - def nll_grad(x: float, y: float) -> Tuple[float, float]: - return nll_grad_numba(x, y, pha, f, xmin, ymin, bin_size, noise) + def nll_grad(x: float, y: float, q: float) -> Tuple[float, float]: + return nll_grad_numba(x, y, q, pha, f, xmin, ymin, bin_size, noise) # Assign a name to the free parameters. - parnames = ["x", "y"] + parnames = ["x", "y", "q"] # Initialize the minimizer. + p0 = (*p0, np.sum(pha)) m = Minuit(nll, *p0, grad=nll_grad, name=parnames) # Set the limits for the free parameters. m.limits["x"] = (xmin, xmax) m.limits["y"] = (ymin, ymax) + m.limits["q"] = (np.sum(pha) * 0.5, np.sum(pha) * 1.5) # Set the initial step sizes for the minimizer. We choose half the bin # size as default value. This value is automatically adjusted by the # minimizer during the fit, so it is not critical to choose a very # precise value. m.errors["x"] = bin_size / 2 m.errors["y"] = bin_size / 2 + m.errors["q"] = np.sum(pha) * 0.1 # Define the strategy for the minimizer. We use the higher strategy to # avoid numerical problems with the hessian. For an explanation of the # strategy levels, please refer to the iminuit documentation. diff --git a/src/hexsample/tasks.py b/src/hexsample/tasks.py index 5f631c8..68cf1b0 100644 --- a/src/hexsample/tasks.py +++ b/src/hexsample/tasks.py @@ -300,7 +300,7 @@ def reconstruct( equalization_matrix=equalization_matrix ) if pos_recon_algorithm == "mle": - clustering = ClusteringHex(readout, 0, pos_recon_algorithm, recon_pars) + clustering = ClusteringHex(readout, -5, pos_recon_algorithm, recon_pars) num_neighbors = 6 else: clustering = ClusteringNN( @@ -416,7 +416,7 @@ def calibrate_position( readout = create_readout(readout_mode, header, *readout_args) # To correctly analyze every type of event, we need a zero suppression threshold # of 0, because the calibration should be performed on zero-noise simulations. - clustering = ClusteringHex(readout, zero_sup_threshold=0.) + clustering = ClusteringHex(readout, zero_sup_threshold=-5.) clustering_nn = ClusteringNN(readout, zero_sup_threshold=zero_sup_threshold, num_neighbors=6, pos_recon_algorithm="centroid") # Initialize the position calibrator and run the event loop. @@ -997,16 +997,23 @@ def calibview( f"lower_quantile={lower_quantile}, upper_quantile={upper_quantile}" ) logger.info(f"Number of calibrated pixels after quality cuts: {np.sum(mask)}") + mask = mask & (matrix.values >= lower_bound) & (matrix.values <= upper_bound) # Plot the values matrix. plt.figure(f"Calibrated matrix: {matrix.metadata['file_name']}") plt.imshow(matrix.values, origin="upper", vmin=lower_bound, vmax=upper_bound) plt.xlabel("Column") plt.ylabel("Row") plt.colorbar(label=unit) + # Plot the entries matrix. + plt.figure(f"Entries matrix: {matrix.metadata['file_name']}") + plt.imshow(matrix.entries, origin="upper") + plt.xlabel("Column") + plt.ylabel("Row") + plt.colorbar(label="Entries") # Plot the distribution of the calibrated values. vals = matrix.values.flatten()[mask.flatten()] edges = np.linspace(lower_bound, upper_bound, 100) - vals_hist = Histogram1d(edges, label="Distribution", xlabel=unit).fill(vals) + vals_hist = Histogram1d(edges, label="Values distribution", xlabel=unit).fill(vals) plt.figure("Distribution of calibrated values") vals_hist.plot(statistics=True) plt.legend() @@ -1051,10 +1058,10 @@ def calibview( plt.xlabel(f"Calibrated values [{unit}]") plt.ylabel(f"Monte Carlo truth values [{mc_unit}]") # Plot the residuals distribution. - residuals = (vals - mc_vals[mask.flatten()]) / mc_vals[mask.flatten()] + residuals = vals - mc_vals[mask.flatten()] residual_edges = np.linspace(np.nanmin(residuals), np.nanmax(residuals), 100) residual_hist = Histogram1d( - residual_edges, label="Residuals", xlabel="Relative Residual" + residual_edges, label="Residuals distribution", xlabel=f"Residuals [{unit}]" ).fill(residuals) plt.figure("Relative residuals distribution") residual_hist.plot(statistics=True) diff --git a/tests/test_likelihood.py b/tests/test_likelihood.py index c193b3f..6618f0b 100644 --- a/tests/test_likelihood.py +++ b/tests/test_likelihood.py @@ -44,7 +44,7 @@ def test_nll_numba(): nll = np.zeros((len(x), len(y))) for i_x, _x in enumerate(x): for i_y, _y in enumerate(y): - nll[i_x, i_y] = nll_numba(_x, _y, pha, f, xbin0, ybin0, bin_size, sigma) + nll[i_x, i_y] = nll_numba(_x, _y, np.sum(pha),pha, f, xbin0, ybin0, bin_size, sigma) plt.figure("test_negative_log_likelihood") plt.imshow(nll.T, extent=(x[0], x[-1], y[0], y[-1]), origin="lower") plt.colorbar(label="Negative log-likelihood")