diff --git a/SOAP/particle_selection/halo_properties.py b/SOAP/particle_selection/halo_properties.py
index 3198ae4a..15e8220e 100644
--- a/SOAP/particle_selection/halo_properties.py
+++ b/SOAP/particle_selection/halo_properties.py
@@ -25,6 +25,7 @@ def __init__(self, cellgrid):
         self.softening_of_parttype = {
             "PartType0": cellgrid.baryon_softening,
             "PartType1": cellgrid.dark_matter_softening,
+            "PartType3": cellgrid.baryon_softening,
             "PartType4": cellgrid.baryon_softening,
             "PartType5": cellgrid.baryon_softening,
             "PartType6": cellgrid.nu_softening,
diff --git a/SOAP/particle_selection/subhalo_properties.py b/SOAP/particle_selection/subhalo_properties.py
index cbb68649..dce0d82b 100644
--- a/SOAP/particle_selection/subhalo_properties.py
+++ b/SOAP/particle_selection/subhalo_properties.py
@@ -179,6 +179,14 @@ def dm_mask_sh(self) -> NDArray[bool]:
         """
         return self.types == 1
 
+    @lazy_property
+    def sink_mask_sh(self) -> NDArray[bool]:
+        """
+        Mask used to mask out sink particles that belong to this subhalo in
+        arrays containing all particles, e.g. self.mass.
+        """
+        return self.types == 3
+
     @lazy_property
     def star_mask_sh(self) -> NDArray[bool]:
         """
@@ -217,6 +225,13 @@ def Ndm(self) -> int:
         """
         return self.dm_mask_sh.sum()
 
+    @lazy_property
+    def Nsink(self) -> int:
+        """
+        Number of sink particles in the subhalo.
+        """
+        return self.sink_mask_sh.sum()
+
     @lazy_property
     def Nstar(self) -> int:
         """
@@ -2540,6 +2555,7 @@ def __init__(
         self.particle_properties = {
             "PartType0": ["Coordinates", "Masses", "Velocities", "GroupNr_bound"],
             "PartType1": ["Coordinates", "Masses", "Velocities", "GroupNr_bound"],
+            "PartType3": ["Coordinates", "Masses", "Velocities", "GroupNr_bound"],
             "PartType4": ["Coordinates", "Masses", "Velocities", "GroupNr_bound"],
             "PartType5": [
                 "Coordinates",
@@ -2603,6 +2619,7 @@ def calculate(self, input_halo, search_radius, data, halo_result):
             {
                 "BoundSubhalo/NumberOfDarkMatterParticles": part_props.Ndm,
                 "BoundSubhalo/NumberOfGasParticles": part_props.Ngas,
+                "BoundSubhalo/NumberOfSinkParticles": part_props.Nsink,
                 "BoundSubhalo/NumberOfStarParticles": part_props.Nstar,
                 "BoundSubhalo/NumberOfBlackHoleParticles": part_props.Nbh,
             },
@@ -2612,7 +2629,7 @@ def calculate(self, input_halo, search_radius, data, halo_result):
         # If not, we need to try again with a larger search radius.
         # For HBT this should not happen since we use the radius of the most distant
         # bound particle.
-        Ntot = part_props.Ngas + part_props.Ndm + part_props.Nstar + part_props.Nbh
+        Ntot = part_props.Ngas + part_props.Ndm + part_props.Nsink + part_props.Nstar + part_props.Nbh
         Nexpected = input_halo["nr_bound_part"]
         if Ntot < Nexpected:
             # Try again with a larger search radius