Calculate fof extent

SOAP/core/combine_chunks.py

@@ -175,7 +175,10 @@ def combine_chunks(
     # Add metadata for FOF properties
     fof_metadata = []
     if (args.fof_group_filename != "") and (args.halo_format == "HBTplus"):
-        for fofkey in ["Centres", "Masses", "Sizes"]:
+        fof_keys = ["Centres", "Masses", "Sizes"]
+        if (args.fof_radius_filename != ""):
+            fof_keys.append('Radii')
+        for fofkey in fof_keys:
             prop = PropertyTable.full_property_list[f"FOF/{fofkey}"]
             name = f"InputHalos/{prop.name}"
             size = prop.shape
@@ -488,6 +491,43 @@ def combine_chunks(
             comm=comm_world,
         )
 
+        # Handle radius differently since SWIFT did not always output radius
+        # Assumes the FOF radii files are the same order as the other FOFs
+        read_radii = 'InputHalos/FOF/Radii' in [m[0] for m in fof_metadata]
+        if read_radii:
+            # Open file in parallel
+            pf = PartialFormatter()
+            fof_filename = pf.format(
+                args.fof_radius_filename, snap_nr=args.snapshot_nr, file_nr=None
+            )
+            fof_file = phdf5.MultiFile(
+                fof_filename,
+                file_nr_attr=("Header", "NumFilesPerSnapshot"),
+                comm=comm_world,
+            )
+
+            fof_radii = np.zeros(keep.shape[0], dtype=np.float32)
+            fof_radii[keep] = psort.fetch_elements(
+                fof_file.read("Groups/Radii"), indices, comm=comm_world
+            )
+            prop = PropertyTable.full_property_list[f"FOF/Radii"]
+            soap_radii_unit = cellgrid.get_unit(prop.unit)
+            physical = prop.output_physical
+            a_exponent = prop.a_scale_exponent
+            if not physical:
+                soap_radii_unit = soap_radii_unit * cellgrid.get_unit("a") ** a_exponent
+            fof_radii = (fof_radii * fof_com_unit).to(soap_radii_unit)
+            phdf5.collective_write(
+                outfile,
+                "InputHalos/FOF/Radii",
+                fof_radii,
+                create_dataset=False,
+                comm=comm_world,
+            )
+        else:
+            if comm_world.Get_rank() == 0:
+                print("Groups/Radii not found in FOF files")
+
     # Calculate the index in the SOAP output of the host field halo (VR) or the central subhalo of the host FOF group (HBTplus)
     if "SOAP/HostHaloIndex" in soap_props:
         with MPITimer("Calculate and write host index of each satellite", comm_world):

SOAP/core/soap_args.py

@@ -122,6 +122,7 @@ def get_soap_args(comm):
     args.halo_basename = all_args["HaloFinder"]["filename"]
     args.halo_format = all_args["HaloFinder"]["type"]
     args.fof_group_filename = all_args["HaloFinder"].get("fof_filename", "")
+    args.fof_radius_filename = all_args["HaloFinder"].get("fof_radius_filename", "")
     args.output_file = all_args["HaloProperties"]["filename"]
     args.snapshot_nr = all_args["Parameters"]["snap_nr"]
     args.chunks = all_args["Parameters"]["chunks"]
@@ -179,5 +180,13 @@ def get_soap_args(comm):
             if not os.path.exists(fof_filename):
                 print(f"Could not find FOF group catalogue: {fof_filename}")
                 comm.Abort(1)
+        if args.fof_radius_filename != "":
+            assert args.fof_group_filename != ""
+            fof_filename = args.fof_radius_filename.format(
+                snap_nr=args.snapshot_nr, file_nr=0
+            )
+            if not os.path.exists(fof_filename):
+                print(f"Could not find FOF radius catalogue: {fof_filename}")
+                comm.Abort(1)
 
     return args

SOAP/property_table.py

@@ -4232,6 +4232,18 @@ class PropertyTable:
             output_physical=True,
             a_scale_exponent=None,
         ),
+        "FOF/Radii": Property(
+            name="FOF/Radii",
+            shape=1,
+            dtype=np.float32,
+            unit="snap_length",
+            description="Radius of the particle furthest from the FOF centre of mass. Zero for satellite and hostless subhalos. Missing for older runs.",
+            lossy_compression_filter="FMantissa9",
+            dmo_property=True,
+            particle_properties=[],
+            output_physical=False,
+            a_scale_exponent=1,
+        ),
         # SOAP properties
         "SOAP/SubhaloRankByBoundMass": Property(
             name="SOAP/SubhaloRankByBoundMass",
@@ -4484,8 +4496,8 @@ def generate_tex_files(self, output_dir: str):
             prop = self.properties[prop_name]
             footnotes = self.get_footnotes(prop_name)
             prop_outputname = f"{prop['name']}{footnotes}"
-            if prop_outputname.split('/')[0] in ['HBTplus', 'VR', 'FOF']:
-                prop_outputname = 'InputHalos/' + prop_outputname
+            if prop_outputname.split("/")[0] in ["HBTplus", "VR", "FOF"]:
+                prop_outputname = "InputHalos/" + prop_outputname
             prop_outputname = word_wrap_name(prop_outputname)
             prop_shape = f'{prop["shape"]}'
             prop_dtype = prop["dtype"]