From 5b78a4de4e2df9809a4480d811b24aaab8036566 Mon Sep 17 00:00:00 2001
From: Augusto Cattafesta <augustocattafesta@gmail.com>
Date: Fri, 6 Feb 2026 14:51:55 +0100
Subject: [PATCH 01/13] Just a draft.

---
 src/hexsample/clustering.py | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/src/hexsample/clustering.py b/src/hexsample/clustering.py
index 7947680..ed5c8af 100644
--- a/src/hexsample/clustering.py
+++ b/src/hexsample/clustering.py
@@ -168,6 +168,22 @@ def zero_suppress(self, array: np.ndarray) -> np.ndarray:
         out[out <= self.zero_sup_threshold] = 0
         return out
 
+    def topology_suppress(self, pha, col, row):
+        # Check if third pixel is neighbor of the second
+        out = pha.copy()
+        if (col[2], row[2]) not in self.grid.neighbors(col[1], row[1]):
+            out[2] = 0
+        # The same for the fourth pixel
+        if (col[3], row[3]) not in self.grid.neighbors(col[1], row[1]):
+            out[3] = 0
+        # Not sure if we want to do this check for 5th and 6th pixels
+        if (col[4], row[4]) not in self.grid.neighbors(col[2], row[2]):
+            if (col[4], row[4]) not in self.grid.neighbors(col[3], row[3]):
+                out[4] = 0
+        
+        return out
+
+
     def run(self, event: DigiEventRectangular) -> Cluster:
         """Workhorse method to be reimplemented by derived classes.
         """

From 939b5f79923898fd62a52322e6ab5a7911036e27 Mon Sep 17 00:00:00 2001
From: Augusto Cattafesta <augustocattafesta@gmail.com>
Date: Fri, 6 Feb 2026 16:54:08 +0100
Subject: [PATCH 02/13] Draft.

---
 src/hexsample/clustering.py | 37 ++++++++++++++++++++++++++-----------
 1 file changed, 26 insertions(+), 11 deletions(-)

diff --git a/src/hexsample/clustering.py b/src/hexsample/clustering.py
index ed5c8af..199cd79 100644
--- a/src/hexsample/clustering.py
+++ b/src/hexsample/clustering.py
@@ -168,20 +168,35 @@ def zero_suppress(self, array: np.ndarray) -> np.ndarray:
         out[out <= self.zero_sup_threshold] = 0
         return out
 
+    def are_neighbors(self, pix0, pix1) -> bool:
+        """Check if two pixels are neighbors.
+        """
+        return pix1 in self.grid.neighbors(*pix0)
+
+
+
     def topology_suppress(self, pha, col, row):
         # Check if third pixel is neighbor of the second
-        out = pha.copy()
-        if (col[2], row[2]) not in self.grid.neighbors(col[1], row[1]):
-            out[2] = 0
-        # The same for the fourth pixel
-        if (col[3], row[3]) not in self.grid.neighbors(col[1], row[1]):
-            out[3] = 0
-        # Not sure if we want to do this check for 5th and 6th pixels
-        if (col[4], row[4]) not in self.grid.neighbors(col[2], row[2]):
-            if (col[4], row[4]) not in self.grid.neighbors(col[3], row[3]):
-                out[4] = 0
+        ind = [0, 1]
+        pix = list(zip(col, row))
+        # Check whether the third pixel is neighbor of the second pixel
+        if self.are_neighbors(pix[ind[1]], pix[2]):
+            ind.append(2)
+        # Check whether the fourth pixel is neighbor of the second pixel
+        if self.are_neighbors(pix[ind[1]], pix[3]):
+            ind.append(3)
         
-        return out
+        if self.are_neighbors(pix[2], pix[4]) or self.are_neighbors(pix[3], pix[4]):
+            ind.append(4)
+        if self.are_neighbors(pix[2], pix[5]) or self.are_neighbors(pix[3], pix[5]):
+            ind.append(5)
+        ind = np.concatenate(ind, np.setdiff1d(np.arange(7), ind))
+
+
+
+        
+
+
 
 
     def run(self, event: DigiEventRectangular) -> Cluster:

From 6b78b2b2c48ea70e082a1f7b435b06f5f20bdaf8 Mon Sep 17 00:00:00 2001
From: Augusto Cattafesta <augustocattafesta@gmail.com>
Date: Fri, 6 Feb 2026 17:24:30 +0100
Subject: [PATCH 03/13] Draft.

---
 src/hexsample/clustering.py | 35 +++++++++++++++++++++--------------
 1 file changed, 21 insertions(+), 14 deletions(-)

diff --git a/src/hexsample/clustering.py b/src/hexsample/clustering.py
index 199cd79..4af1352 100644
--- a/src/hexsample/clustering.py
+++ b/src/hexsample/clustering.py
@@ -173,31 +173,33 @@ def are_neighbors(self, pix0, pix1) -> bool:
         """
         return pix1 in self.grid.neighbors(*pix0)
 
-
-
     def topology_suppress(self, pha, col, row):
-        # Check if third pixel is neighbor of the second
+        print()
+        print(pha)
+        n = np.count_nonzero(pha)
+        if n <= 2:
+            print("<= 2 PIXEL")
+            return pha
         ind = [0, 1]
         pix = list(zip(col, row))
         # Check whether the third pixel is neighbor of the second pixel
-        if self.are_neighbors(pix[ind[1]], pix[2]):
+        if self.are_neighbors(pix[1], pix[2]):
             ind.append(2)
         # Check whether the fourth pixel is neighbor of the second pixel
-        if self.are_neighbors(pix[ind[1]], pix[3]):
+        if self.are_neighbors(pix[1], pix[3]):
             ind.append(3)
-        
+
+
         if self.are_neighbors(pix[2], pix[4]) or self.are_neighbors(pix[3], pix[4]):
             ind.append(4)
         if self.are_neighbors(pix[2], pix[5]) or self.are_neighbors(pix[3], pix[5]):
             ind.append(5)
-        ind = np.concatenate(ind, np.setdiff1d(np.arange(7), ind))
-
-
-
-        
-
-
-
+        diff = np.setdiff1d(np.arange(0, len(pha)), ind)
+        out = pha.copy()
+        out[diff] = 0
+        if not np.array_equal(out, pha):
+            print("OUT: ", out)
+        return out
 
     def run(self, event: DigiEventRectangular) -> Cluster:
         """Workhorse method to be reimplemented by derived classes.
@@ -269,6 +271,11 @@ def run(self, event) -> Cluster:
         # This is useless for the circular readout because in that case all
         # neighbors are used for track reconstruction.
         mask = idx[:self.num_neighbors + 1]
+
+
+        pha = self.topology_suppress(pha[mask], col[mask], row[mask])
+
+
         # If there's any zero left in the target pixels, get rid of it.
         mask = mask[pha[mask] > 0]
         # Trim the relevant arrays.

From 79faf996ef3d0ee538cefcbe9a08e750e1a5890b Mon Sep 17 00:00:00 2001
From: Augusto Cattafesta <augustocattafesta@gmail.com>
Date: Fri, 6 Feb 2026 17:33:42 +0100
Subject: [PATCH 04/13] Draft.

---
 src/hexsample/clustering.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/hexsample/clustering.py b/src/hexsample/clustering.py
index 4af1352..c39f6ec 100644
--- a/src/hexsample/clustering.py
+++ b/src/hexsample/clustering.py
@@ -190,9 +190,9 @@ def topology_suppress(self, pha, col, row):
             ind.append(3)
 
 
-        if self.are_neighbors(pix[2], pix[4]) or self.are_neighbors(pix[3], pix[4]):
+        if self.are_neighbors(pix[ind[-1]], pix[4]) or (self.are_neighbors(pix[ind[-2]], pix[4]) and len(ind) > 2):
             ind.append(4)
-        if self.are_neighbors(pix[2], pix[5]) or self.are_neighbors(pix[3], pix[5]):
+        if self.are_neighbors(pix[ind[-1]], pix[5]) or (self.are_neighbors(pix[ind[-2]], pix[5]) and len(ind) > 2):
             ind.append(5)
         diff = np.setdiff1d(np.arange(0, len(pha)), ind)
         out = pha.copy()

From 403c7501b06ac3f1f96f4ac3811f64abf58b7d25 Mon Sep 17 00:00:00 2001
From: Augusto Cattafesta <augustocattafesta@gmail.com>
Date: Sat, 7 Feb 2026 09:44:22 +0100
Subject: [PATCH 05/13] Topology suppression.

---
 src/hexsample/clustering.py | 51 ++++++++++++++++++++++++++++++++++++-
 1 file changed, 50 insertions(+), 1 deletion(-)

diff --git a/src/hexsample/clustering.py b/src/hexsample/clustering.py
index c39f6ec..e13042a 100644
--- a/src/hexsample/clustering.py
+++ b/src/hexsample/clustering.py
@@ -171,9 +171,57 @@ def zero_suppress(self, array: np.ndarray) -> np.ndarray:
     def are_neighbors(self, pix0, pix1) -> bool:
         """Check if two pixels are neighbors.
         """
+        # Eliminala non serve, non abbrevia
         return pix1 in self.grid.neighbors(*pix0)
 
+
     def topology_suppress(self, pha, col, row):
+        print()
+        print(pha)
+        # If the cluster size is 1 or 2, we don't need to check the topology
+        if np.count_nonzero(pha) <= 2:
+            return pha
+        ind = [0, 1]
+        pix = list(zip(col, row))
+        # Otherwise, we start checking if pixel 3 and 4 are neighbors of pixel 2
+        for i in range(2, 4):
+            if pix[i] in self.grid.neighbors(*pix[1]):
+                ind.append(i)
+        # If both pixel 3 and 4 are neighbors of pixel 2, we check if pixel 5 and 6 are neighbors
+        # of pixel 3 or 4.
+        if len(ind) == 4:
+            for i in range(4, 6):
+                if pix[i] in self.grid.neighbors(*pix[2]) or pix[i] in self.grid.neighbors(*pix[3]):
+                    ind.append(i)
+        # If only one between pixel 3 and 4 is a neighbor of pixel 2, we should check if pixel 5
+        # is a neighbor of pixel 2.
+        elif len(ind) == 3:
+            if pix[4] in self.grid.neighbors(*pix[1]):
+                ind.append(4)
+                # If this condition is true, then we have to check if pixel 6 is a neighbor of the
+                # third good pixel (3 or 4) and pixel 5 (the fourth good pixel).
+                if pix[5] in self.grid.neighbors(*pix[ind[-1]]) or pix[5] in self.grid.neighbors(*pix[ind[-2]]):
+                    ind.append(5)
+            else:
+                # If pixel 5 is not a neighbor of pixel 2, we check if pixel 6 is a neighbor of pixel 2.
+                if pix[5] in self.grid.neighbors(*pix[1]):
+                    ind.append(5)
+        # If none of pixel 3 and 4 is a neighbor of pixel 2, we have to check if pixel 5 and 6 are
+        # neighbors of pixel 2.
+        elif len(ind) == 2:
+            for i in range(4, 6):
+                if pix[i] in self.grid.neighbors(*pix[1]):
+                    ind.append(i)
+        # Finally, we set to zero the pixels that are not in the good topology.
+        diff = np.setdiff1d(np.arange(0, len(pha)), ind)
+        out = pha.copy()
+        out[diff] = 0
+        if not np.array_equal(out, pha):
+            print("OUT: ", out)
+        return out
+                    
+
+    def old_topology_suppress(self, pha, col, row):
         print()
         print(pha)
         n = np.count_nonzero(pha)
@@ -183,6 +231,7 @@ def topology_suppress(self, pha, col, row):
         ind = [0, 1]
         pix = list(zip(col, row))
         # Check whether the third pixel is neighbor of the second pixel
+        # if pix[0] in grid.neighbors(*pix[1]):
         if self.are_neighbors(pix[1], pix[2]):
             ind.append(2)
         # Check whether the fourth pixel is neighbor of the second pixel
@@ -192,7 +241,7 @@ def topology_suppress(self, pha, col, row):
 
         if self.are_neighbors(pix[ind[-1]], pix[4]) or (self.are_neighbors(pix[ind[-2]], pix[4]) and len(ind) > 2):
             ind.append(4)
-        if self.are_neighbors(pix[ind[-1]], pix[5]) or (self.are_neighbors(pix[ind[-2]], pix[5]) and len(ind) > 2):
+        if self.are_neighbors(pix[ind[-2]], pix[5]) or (self.are_neighbors(pix[ind[-1]], pix[5]) and len(ind) > 2):
             ind.append(5)
         diff = np.setdiff1d(np.arange(0, len(pha)), ind)
         out = pha.copy()

From fec22bfc675ae55a043e5376bb37dac8ed570091 Mon Sep 17 00:00:00 2001
From: Augusto Cattafesta <augustocattafesta@gmail.com>
Date: Sat, 7 Feb 2026 09:54:17 +0100
Subject: [PATCH 06/13] Minor.

---
 src/hexsample/clustering.py | 60 +++++++------------------------------
 1 file changed, 11 insertions(+), 49 deletions(-)

diff --git a/src/hexsample/clustering.py b/src/hexsample/clustering.py
index e13042a..b96f7ce 100644
--- a/src/hexsample/clustering.py
+++ b/src/hexsample/clustering.py
@@ -168,16 +168,7 @@ def zero_suppress(self, array: np.ndarray) -> np.ndarray:
         out[out <= self.zero_sup_threshold] = 0
         return out
 
-    def are_neighbors(self, pix0, pix1) -> bool:
-        """Check if two pixels are neighbors.
-        """
-        # Eliminala non serve, non abbrevia
-        return pix1 in self.grid.neighbors(*pix0)
-
-
     def topology_suppress(self, pha, col, row):
-        print()
-        print(pha)
         # If the cluster size is 1 or 2, we don't need to check the topology
         if np.count_nonzero(pha) <= 2:
             return pha
@@ -216,39 +207,12 @@ def topology_suppress(self, pha, col, row):
         diff = np.setdiff1d(np.arange(0, len(pha)), ind)
         out = pha.copy()
         out[diff] = 0
-        if not np.array_equal(out, pha):
-            print("OUT: ", out)
-        return out
-                    
-
-    def old_topology_suppress(self, pha, col, row):
-        print()
-        print(pha)
-        n = np.count_nonzero(pha)
-        if n <= 2:
-            print("<= 2 PIXEL")
-            return pha
-        ind = [0, 1]
-        pix = list(zip(col, row))
-        # Check whether the third pixel is neighbor of the second pixel
-        # if pix[0] in grid.neighbors(*pix[1]):
-        if self.are_neighbors(pix[1], pix[2]):
-            ind.append(2)
-        # Check whether the fourth pixel is neighbor of the second pixel
-        if self.are_neighbors(pix[1], pix[3]):
-            ind.append(3)
-
-
-        if self.are_neighbors(pix[ind[-1]], pix[4]) or (self.are_neighbors(pix[ind[-2]], pix[4]) and len(ind) > 2):
-            ind.append(4)
-        if self.are_neighbors(pix[ind[-2]], pix[5]) or (self.are_neighbors(pix[ind[-1]], pix[5]) and len(ind) > 2):
-            ind.append(5)
-        diff = np.setdiff1d(np.arange(0, len(pha)), ind)
-        out = pha.copy()
-        out[diff] = 0
-        if not np.array_equal(out, pha):
-            print("OUT: ", out)
-        return out
+        # Sort the pha array in decreasing order and reorder the col and row arrays
+        idx = np.argsort(-out)
+        out = out[idx]
+        col = col[idx]
+        row = row[idx]
+        return out, col, row
 
     def run(self, event: DigiEventRectangular) -> Cluster:
         """Workhorse method to be reimplemented by derived classes.
@@ -256,7 +220,6 @@ def run(self, event: DigiEventRectangular) -> Cluster:
         raise NotImplementedError
 
 
-
 @dataclass
 class ClusteringNN(ClusteringBase):
 
@@ -320,13 +283,12 @@ def run(self, event) -> Cluster:
         # This is useless for the circular readout because in that case all
         # neighbors are used for track reconstruction.
         mask = idx[:self.num_neighbors + 1]
-
-
-        pha = self.topology_suppress(pha[mask], col[mask], row[mask])
-
-
+        # Sort the arrays in decreasing order before applying the topology suppression.
+        pha, col, row = self.topology_suppress(pha[mask], col[mask], row[mask])
+        # The returned arrays have alredy been re-sorted in decreasing order, so we just need
+        # to remove the pixels with zero pha. 
         # If there's any zero left in the target pixels, get rid of it.
-        mask = mask[pha[mask] > 0]
+        mask = mask[pha > 0]
         # Trim the relevant arrays.
         col = col[mask]
         row = row[mask]

From 688e208bccae125ab1fbf646a43dfce0501a5437 Mon Sep 17 00:00:00 2001
From: Augusto Cattafesta <augustocattafesta@gmail.com>
Date: Sat, 7 Feb 2026 10:54:44 +0100
Subject: [PATCH 07/13] Function optimization.

---
 src/hexsample/clustering.py | 84 ++++++++++++++++++++++++++-----------
 1 file changed, 60 insertions(+), 24 deletions(-)

diff --git a/src/hexsample/clustering.py b/src/hexsample/clustering.py
index b96f7ce..72c09dc 100644
--- a/src/hexsample/clustering.py
+++ b/src/hexsample/clustering.py
@@ -168,51 +168,87 @@ def zero_suppress(self, array: np.ndarray) -> np.ndarray:
         out[out <= self.zero_sup_threshold] = 0
         return out
 
-    def topology_suppress(self, pha, col, row):
-        # If the cluster size is 1 or 2, we don't need to check the topology
-        if np.count_nonzero(pha) <= 2:
-            return pha
+    def topology_suppress(self, pha: np.ndarray, col: np.ndarray, row: np.ndarray
+                          ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
+        """Suppress pixels that do not have the right topology.
+
+        In a well-formed topology, the highest pixel is located at the center of the cluster.
+        For clusters larger than two pixels, the charge distribution should be symmetric with
+        respect to the axis connecting the two highest pixels. Furthermore, the charge must
+        decrease monotonically as the distance from the second-highest pixel increases.
+
+        Arguments
+        ---------
+        pha : np.ndarray
+            The array of pulse heights of the pixels in the cluster, ordered in decreasing order.
+        col : np.ndarray
+            The array of column indexes of the pixels in the cluster.
+        row : np.ndarray
+            The array of row indexes of the pixels in the cluster.
+        
+        Returns
+        -------
+        pha : np.ndarray
+            The array of pulse heights of the pixels in the cluster after topology suppression,
+            ordered in decreasing order.
+        col : np.ndarray
+            The array of column indexes of the pixels in the cluster after topology suppression,
+            ordered in decreasing order of pulse height.
+        row : np.ndarray
+            The array of row indexes of the pixels in the cluster after topology suppression,
+            ordered in decreasing order of pulse height.
+        """
+        # If the number of pixels above the zero suppression threshold or the cluster size is less
+        # than or equal to 2, we don't need to check the topology. We need to check the cluster
+        # size because with rectangular readout we might have clusters with less than 7 pixels.
+        if np.count_nonzero(pha) <= 2 or pha.size <= 2:
+            return pha, col, row
         ind = [0, 1]
         pix = list(zip(col, row))
-        # Otherwise, we start checking if pixel 3 and 4 are neighbors of pixel 2
-        for i in range(2, 4):
-            if pix[i] in self.grid.neighbors(*pix[1]):
-                ind.append(i)
+        # Otherwise, we start checking if pixel 3 and 4 are neighbors of pixel 2. We also need to
+        # check the cluster dimension
+        ind.extend([
+            i for i in range(2, 4)
+            if pha.size > i and pix[i] in self.grid.neighbors(*pix[1])
+            ])
         # If both pixel 3 and 4 are neighbors of pixel 2, we check if pixel 5 and 6 are neighbors
         # of pixel 3 or 4.
         if len(ind) == 4:
-            for i in range(4, 6):
-                if pix[i] in self.grid.neighbors(*pix[2]) or pix[i] in self.grid.neighbors(*pix[3]):
-                    ind.append(i)
+            ind.extend([
+                i for i in range(4, 6)
+                if pha.size > i and (
+                    pix[i] in self.grid.neighbors(*pix[2]) or
+                    pix[i] in self.grid.neighbors(*pix[3]))
+            ])
         # If only one between pixel 3 and 4 is a neighbor of pixel 2, we should check if pixel 5
         # is a neighbor of pixel 2.
         elif len(ind) == 3:
-            if pix[4] in self.grid.neighbors(*pix[1]):
+            if pha.size > 4 and pix[4] in self.grid.neighbors(*pix[1]):
                 ind.append(4)
                 # If this condition is true, then we have to check if pixel 6 is a neighbor of the
                 # third good pixel (3 or 4) and pixel 5 (the fourth good pixel).
-                if pix[5] in self.grid.neighbors(*pix[ind[-1]]) or pix[5] in self.grid.neighbors(*pix[ind[-2]]):
+                if pha.size > 5 and ((pix[5] in self.grid.neighbors(*pix[ind[-1]])) or
+                   (pix[5] in self.grid.neighbors(*pix[ind[-2]]))):
                     ind.append(5)
             else:
-                # If pixel 5 is not a neighbor of pixel 2, we check if pixel 6 is a neighbor of pixel 2.
-                if pix[5] in self.grid.neighbors(*pix[1]):
+                # If pixel 5 is not a neighbor of pixel 2, we check if pixel 6 is a neighbor of
+                # pixel 2.
+                if pha.size > 5 and pix[5] in self.grid.neighbors(*pix[1]):
                     ind.append(5)
         # If none of pixel 3 and 4 is a neighbor of pixel 2, we have to check if pixel 5 and 6 are
         # neighbors of pixel 2.
         elif len(ind) == 2:
-            for i in range(4, 6):
-                if pix[i] in self.grid.neighbors(*pix[1]):
-                    ind.append(i)
+            ind.extend([
+                i for i in range(4, 6)
+                if pha.size > i and pix[i] in self.grid.neighbors(*pix[1])
+            ])
         # Finally, we set to zero the pixels that are not in the good topology.
-        diff = np.setdiff1d(np.arange(0, len(pha)), ind)
         out = pha.copy()
+        diff = np.setdiff1d(np.arange(0, len(out)), ind)
         out[diff] = 0
         # Sort the pha array in decreasing order and reorder the col and row arrays
         idx = np.argsort(-out)
-        out = out[idx]
-        col = col[idx]
-        row = row[idx]
-        return out, col, row
+        return out[idx], col[idx], row[idx]
 
     def run(self, event: DigiEventRectangular) -> Cluster:
         """Workhorse method to be reimplemented by derived classes.
@@ -286,7 +322,7 @@ def run(self, event) -> Cluster:
         # Sort the arrays in decreasing order before applying the topology suppression.
         pha, col, row = self.topology_suppress(pha[mask], col[mask], row[mask])
         # The returned arrays have alredy been re-sorted in decreasing order, so we just need
-        # to remove the pixels with zero pha. 
+        # to remove the pixels with zero pha.
         # If there's any zero left in the target pixels, get rid of it.
         mask = mask[pha > 0]
         # Trim the relevant arrays.

From 500354a92620c815152fe06a9e2eea8c6264b51a Mon Sep 17 00:00:00 2001
From: Augusto Cattafesta <augustocattafesta@gmail.com>
Date: Sat, 7 Feb 2026 11:13:04 +0100
Subject: [PATCH 08/13] Minor.

---
 src/hexsample/clustering.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/hexsample/clustering.py b/src/hexsample/clustering.py
index 72c09dc..e020e5d 100644
--- a/src/hexsample/clustering.py
+++ b/src/hexsample/clustering.py
@@ -321,7 +321,7 @@ def run(self, event) -> Cluster:
         mask = idx[:self.num_neighbors + 1]
         # Sort the arrays in decreasing order before applying the topology suppression.
         pha, col, row = self.topology_suppress(pha[mask], col[mask], row[mask])
-        # The returned arrays have alredy been re-sorted in decreasing order, so we just need
+        # The returned arrays have already been re-sorted in decreasing order, so we just need
         # to remove the pixels with zero pha.
         # If there's any zero left in the target pixels, get rid of it.
         mask = mask[pha > 0]

From 6f7dd69a01ca46a07b0221c51fd071bb4db97992 Mon Sep 17 00:00:00 2001
From: Augusto Cattafesta <augustocattafesta@gmail.com>
Date: Mon, 9 Feb 2026 14:38:21 +0100
Subject: [PATCH 09/13] Suppression algorithmh modified.

---
 src/hexsample/clustering.py | 96 +++++++++++--------------------------
 1 file changed, 28 insertions(+), 68 deletions(-)

diff --git a/src/hexsample/clustering.py b/src/hexsample/clustering.py
index 7843105..224458b 100644
--- a/src/hexsample/clustering.py
+++ b/src/hexsample/clustering.py
@@ -174,14 +174,13 @@ def zero_suppress(self, array: np.ndarray) -> np.ndarray:
         out[out <= self.zero_sup_threshold] = 0
         return out
 
-    def topology_suppress(self, pha: np.ndarray, col: np.ndarray, row: np.ndarray
+    def position_suppress(self, pha: np.ndarray, col: np.ndarray, row: np.ndarray
                           ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
-        """Suppress pixels that do not have the right topology.
+        """Suppress pixels in the cluster that do not satisfy the position requirements.
 
-        In a well-formed topology, the highest pixel is located at the center of the cluster.
-        For clusters larger than two pixels, the charge distribution should be symmetric with
-        respect to the axis connecting the two highest pixels. Furthermore, the charge must
-        decrease monotonically as the distance from the second-highest pixel increases.
+        If a cluster has 2 or less pixels above threshold, don't do anything. If it has more than
+        2 pixels, only keep the highest neighbor of the second pixel apart from the first two
+        pixels.
 
         Arguments
         ---------
@@ -195,66 +194,35 @@ def topology_suppress(self, pha: np.ndarray, col: np.ndarray, row: np.ndarray
         Returns
         -------
         pha : np.ndarray
-            The array of pulse heights of the pixels in the cluster after topology suppression,
+            The array of pulse heights of the pixels in the cluster after position suppression,
             ordered in decreasing order.
         col : np.ndarray
-            The array of column indexes of the pixels in the cluster after topology suppression,
+            The array of column indexes of the pixels in the cluster after position suppression,
             ordered in decreasing order of pulse height.
         row : np.ndarray
-            The array of row indexes of the pixels in the cluster after topology suppression,
+            The array of row indexes of the pixels in the cluster after position suppression,
             ordered in decreasing order of pulse height.
         """
-        # If the number of pixels above the zero suppression threshold or the cluster size is less
-        # than or equal to 2, we don't need to check the topology. We need to check the cluster
-        # size because with rectangular readout we might have clusters with less than 7 pixels.
+        # If we have 2 or less pixels above threshold, we can't do anything, just remove the zeros.
         if np.count_nonzero(pha) <= 2 or pha.size <= 2:
-            return pha, col, row
-        ind = [0, 1]
+            mask = pha > 0
+            return pha[mask], col[mask], row[mask]
+        # For events with more than 2 pixels above threshold, we keep only the highest neighbor of
+        # the second pixel.
         pix = list(zip(col, row))
-        # Otherwise, we start checking if pixel 3 and 4 are neighbors of pixel 2. We also need to
-        # check the cluster dimension
-        ind.extend([
-            i for i in range(2, 4)
-            if pha.size > i and pix[i] in self.grid.neighbors(*pix[1])
-            ])
-        # If both pixel 3 and 4 are neighbors of pixel 2, we check if pixel 5 and 6 are neighbors
-        # of pixel 3 or 4.
-        if len(ind) == 4:
-            ind.extend([
-                i for i in range(4, 6)
-                if pha.size > i and (
-                    pix[i] in self.grid.neighbors(*pix[2]) or
-                    pix[i] in self.grid.neighbors(*pix[3]))
-            ])
-        # If only one between pixel 3 and 4 is a neighbor of pixel 2, we should check if pixel 5
-        # is a neighbor of pixel 2.
-        elif len(ind) == 3:
-            if pha.size > 4 and pix[4] in self.grid.neighbors(*pix[1]):
-                ind.append(4)
-                # If this condition is true, then we have to check if pixel 6 is a neighbor of the
-                # third good pixel (3 or 4) and pixel 5 (the fourth good pixel).
-                if pha.size > 5 and ((pix[5] in self.grid.neighbors(*pix[ind[-1]])) or
-                   (pix[5] in self.grid.neighbors(*pix[ind[-2]]))):
-                    ind.append(5)
-            else:
-                # If pixel 5 is not a neighbor of pixel 2, we check if pixel 6 is a neighbor of
-                # pixel 2.
-                if pha.size > 5 and pix[5] in self.grid.neighbors(*pix[1]):
-                    ind.append(5)
-        # If none of pixel 3 and 4 is a neighbor of pixel 2, we have to check if pixel 5 and 6 are
-        # neighbors of pixel 2.
-        elif len(ind) == 2:
-            ind.extend([
-                i for i in range(4, 6)
-                if pha.size > i and pix[i] in self.grid.neighbors(*pix[1])
-            ])
-        # Finally, we set to zero the pixels that are not in the good topology.
-        out = pha.copy()
-        diff = np.setdiff1d(np.arange(0, len(out)), ind)
-        out[diff] = 0
-        # Sort the pha array in decreasing order and reorder the col and row arrays
-        idx = np.argsort(-out)
-        return out[idx], col[idx], row[idx]
+        # We find the neighbors of the second pixel.
+        neighbors = set(self.grid.neighbors(col[1], row[1]))
+        # We always keep the first two pixels, plus the two neighbors of the second pixel.
+        mask = [True, True] + [p in neighbors for p in pix[2:]]
+        mask = mask & (pha > 0)
+        # Throw away the pixels that are not neighbors of the second pixel and that are zero.
+        out_pha = pha[mask]
+        out_col = col[mask]
+        out_row = row[mask]
+        # Sort the arrays in decreasing order of pulse height.
+        idx = np.argsort(-out_pha)[:3]
+        return out_pha[idx], out_col[idx], out_row[idx]
+
 
     def run(self, event: DigiEventRectangular) -> Cluster:
         """Workhorse method to be reimplemented by derived classes.
@@ -325,15 +293,7 @@ def run(self, event) -> Cluster:
         # This is useless for the circular readout because in that case all
         # neighbors are used for track reconstruction.
         mask = idx[:self.num_neighbors + 1]
-        # Sort the arrays in decreasing order before applying the topology suppression.
-        pha, col, row = self.topology_suppress(pha[mask], col[mask], row[mask])
-        # The returned arrays have already been re-sorted in decreasing order, so we just need
-        # to remove the pixels with zero pha.
-        # If there's any zero left in the target pixels, get rid of it.
-        mask = mask[pha > 0]
-        # Trim the relevant arrays.
-        col = col[mask]
-        row = row[mask]
-        pha = pha[mask]
+        # Sort the arrays in decreasing order before applying the position suppression.
+        pha, col, row = self.position_suppress(pha[mask], col[mask], row[mask])
         x, y = self.grid.pixel_to_world(col, row)
         return Cluster(x, y, pha)

From ac5251ddbe561b23de57a967a841fb2f06e7660d Mon Sep 17 00:00:00 2001
From: Augusto Cattafesta <augustocattafesta@gmail.com>
Date: Mon, 9 Feb 2026 14:42:06 +0100
Subject: [PATCH 10/13] Release notes.

---
 docs/release_notes.rst | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/docs/release_notes.rst b/docs/release_notes.rst
index 4423fab..af208c2 100644
--- a/docs/release_notes.rst
+++ b/docs/release_notes.rst
@@ -11,6 +11,8 @@ Release notes
 * `--max-neighbors` option added to the cli interface.
 * eta reconstructions modified to use the centroid for events with more than
   three pixels.
+* New suppression step added to accept at most three pixels in a cluster, with the condition
+  that all the pixels are neighbors.
 * Pull requests merged and issues closed:
 
   - https://github.com/lucabaldini/hexsample/pull/94

From 9a667c0c245580129c38043682c35dbe72505c36 Mon Sep 17 00:00:00 2001
From: Augusto Cattafesta <augustocattafesta@gmail.com>
Date: Mon, 9 Feb 2026 14:42:22 +0100
Subject: [PATCH 11/13] Release notes.

---
 docs/release_notes.rst | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/docs/release_notes.rst b/docs/release_notes.rst
index af208c2..a49a7d2 100644
--- a/docs/release_notes.rst
+++ b/docs/release_notes.rst
@@ -14,11 +14,12 @@ Release notes
 * New suppression step added to accept at most three pixels in a cluster, with the condition
   that all the pixels are neighbors.
 * Pull requests merged and issues closed:
-
+  - https://github.com/lucabaldini/hexsample/pull/9
   - https://github.com/lucabaldini/hexsample/pull/94
   - https://github.com/lucabaldini/hexsample/pull/92
 
 
+
 Version 0.13.2 (2026-02-04)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 

From 7d8964f5d05c4bac81301cd0ffd98a87a9895089 Mon Sep 17 00:00:00 2001
From: Augusto Cattafesta <augustocattafesta@gmail.com>
Date: Tue, 10 Feb 2026 15:36:15 +0100
Subject: [PATCH 12/13] Minor.

---
 docs/release_notes.rst | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/release_notes.rst b/docs/release_notes.rst
index a49a7d2..1bacb0e 100644
--- a/docs/release_notes.rst
+++ b/docs/release_notes.rst
@@ -14,6 +14,7 @@ Release notes
 * New suppression step added to accept at most three pixels in a cluster, with the condition
   that all the pixels are neighbors.
 * Pull requests merged and issues closed:
+
   - https://github.com/lucabaldini/hexsample/pull/9
   - https://github.com/lucabaldini/hexsample/pull/94
   - https://github.com/lucabaldini/hexsample/pull/92

From 0a9d1b2e4318b5735af54cca3d486495a0608951 Mon Sep 17 00:00:00 2001
From: Augusto Cattafesta <augustocattafesta@gmail.com>
Date: Tue, 10 Feb 2026 15:39:18 +0100
Subject: [PATCH 13/13] Release notes.

---
 docs/release_notes.rst      | 1 -
 src/hexsample/clustering.py | 7 ++++---
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/docs/release_notes.rst b/docs/release_notes.rst
index 1bacb0e..109827a 100644
--- a/docs/release_notes.rst
+++ b/docs/release_notes.rst
@@ -20,7 +20,6 @@ Release notes
   - https://github.com/lucabaldini/hexsample/pull/92
 
 
-
 Version 0.13.2 (2026-02-04)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
diff --git a/src/hexsample/clustering.py b/src/hexsample/clustering.py
index 224458b..ec5942b 100644
--- a/src/hexsample/clustering.py
+++ b/src/hexsample/clustering.py
@@ -178,9 +178,10 @@ def position_suppress(self, pha: np.ndarray, col: np.ndarray, row: np.ndarray
                           ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
         """Suppress pixels in the cluster that do not satisfy the position requirements.
 
-        If a cluster has 2 or less pixels above threshold, don't do anything. If it has more than
-        2 pixels, only keep the highest neighbor of the second pixel apart from the first two
-        pixels.
+        If the cluster contains 2 or fewer pixels, no action is taken. For clusters with 
+        more than 2 pixels, the algorithm retains the two most charged pixels and 
+        only one additional neighbor (the one with the highest charge) of the second 
+        pixel, discarding all others.
 
         Arguments
         ---------