MLE final model

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)...

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()
 
 

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])

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.

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)

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")