fix HMF

halo_mass_functions.py

@@ -10,6 +10,7 @@ from matplotlib.figure import Figure
 
 from paths import base_dir, has_1024_simulations
 from read_vr_files import read_velo_halos
+from readfiles import read_gadget_halos
 from utils import print_progress, figsize_from_page_fraction
 
 
@@ -46,17 +47,26 @@ def monofonic_tests():
     print(cosmology.getCurrent())
     x = None
     for i, waveform in enumerate(
-            ["DB2", "shannon", "shannon_rehalo", "resim_master", "resim_newrandom", "resim_newerrandom",
+            [
-             "resim_newswift"]):
+                "DB2", "shannon", "shannon_rehalo",
+                # "resim_master", "resim_newrandom", "resim_newerrandom",
+                # "resim_newswift",
+                "resim_mastergadget"
+            ]):
         for j, resolution in enumerate(resolutions):
             if (waveform == "shannon_rehalo" or "resim" in waveform) and resolution != 128:
                 continue
             print(waveform, resolution)
             dir = base_dir / f"{waveform}_{resolution}_100"
-            halos = read_velo_halos(dir)
+            if "gadget" in waveform:
-
+                halos = read_gadget_halos(dir)
-            halo_masses: np.ndarray = halos["Mvir"].to_numpy() * 1e10 * 0.67742
+                halo_masses: np.ndarray = halos["Masses"].to_numpy() * 1e10
+            else:
+                halos = read_velo_halos(dir)
+                halo_masses: np.ndarray = halos["Mvir"].to_numpy() * 1e10 * 0.67742
             # halo_masses = halo_masses[halo_masses > 0]
+            halo_masses.sort()
+            print(halo_masses[-4:])
 
             # halos = read_halo_file(dir/"fof_output_0004.hdf5")
             # halo_masses: np.ndarray = halos["Masses"].to_numpy() * 1e10 * 0.67742
@@ -68,7 +78,7 @@ def monofonic_tests():
                 number_densities,
                 lower_error_limit,
                 upper_error_limit,
-            ) = halo_mass_function(halo_masses)
+            ) = halo_mass_function(halo_masses, sim_volume=(100 * 0.67742) ** 3)
 
             ax.set_xscale("log")
             ax.set_yscale("log")
@@ -79,11 +89,12 @@ def monofonic_tests():
             name = f"{waveform} {resolution}"
             ax.step(
                 left_edges,
-                number_densities / (0.67742 ** 3),
+                number_densities,
                 where="post",
                 color=f"C{i + 1}",
                 linestyle=linestyles[j],
                 label=name,
+                zorder=10
             )
 
             # ax.fill_between(
@@ -99,12 +110,31 @@ def monofonic_tests():
         # break
 
     mfunc = mass_function.massFunction(
-        x, 1, mdef="vir", model="tinker08", q_out="dndlnM"
+        x, z=0, mdef="vir", model="tinker08", q_out="dndlnM"
     )
 
     ax.plot(x, mfunc, label="tinker08 (vir)", color="C0")
-    ax.set_xlabel("M")
+    ax.set_xlabel("Mass ($h^{-1} M_{\odot}$)")
-    ax.set_ylabel("dndlnM")
+    ax.set_ylabel(r"$\frac{dN}{d\log M}$ ($h^{3}Mpc^{-3}$)")
+
+    # s: GadgetHDFSnap = pynbody.load(str(base_dir / "resim_mastergadget_128_100" / "output" / "snapshot_019.hdf5"))
+    # s.physical_units()
+    # h3 = s.properties["h"] ** 3
+    # bin_center, bin_counts, err = pynbody.analysis.hmf.simulation_halo_mass_function(
+    #     s, log_M_min=10, log_M_max=15, delta_log_M=0.1
+    # )
+    # ax.errorbar(
+    #     bin_center, bin_counts * h3,
+    #     yerr=err,
+    #     fmt='o',
+    #     capthick=2,
+    #     elinewidth=2,
+    #     color='darkgoldenrod'
+    # )
+    # m, sig, dn_dlogm = pynbody.analysis.hmf.halo_mass_function(
+    #     s, log_M_min=10, log_M_max=15, delta_log_M=0.1, kern="ST"
+    # )
+    # ax.plot(m, dn_dlogm, color='darkmagenta', linewidth=2)
 
     ax.legend()
     fig.tight_layout()
@@ -202,7 +232,7 @@ def hmf_from_rockstar_tree(file: Path):
     cosmology.setCosmology("aurigaCosmo")
     print(cosmology.getCurrent())
     mfunc = mass_function.massFunction(
-        left_edges, 1, mdef="vir", model="tinker08", q_out="dndlnM"
+        left_edges, z=0, mdef="vir", model="tinker08", q_out="dndlnM"
     )
 
     ax.plot(left_edges, mfunc)

hmf_test.py

@@ -1,28 +1,34 @@
 import matplotlib.pyplot as plt
-from velociraptor import load, tools
 
+import read_vr_files
 from paths import base_dir
-from read_vr_files import read_velo_halos
-from unyt import Mpc
-import unyt
 
-for i, waveform in enumerate(["DB2", "DB4", "DB8", "shannon"]):
+for i, waveform in enumerate(
+        ["DB2", "DB8", "shannon"]):
     for j, resolution in enumerate([128, 256, 512]):
+        if (waveform == "shannon_rehalo" or "resim" in waveform) and resolution != 128:
+            continue
         print(waveform, resolution)
         dir = base_dir / f"{waveform}_{resolution}_100"
-        data = load(str(dir)+"/vroutput.properties.0")
+        # data = load(str(dir) + "/vroutput.properties.0")
-        print(data)
+        # print(data)
-        masses_200crit = data.masses.mass_200crit
+        # masses_200crit = data.masses.mass_200crit
-        masses_200crit.convert_to_units("msun")
+        # masses_200crit.convert_to_units("msun")
-        lowest_halo_mass = 1e9 * unyt.msun
+        # lowest_halo_mass = 1e9 * unyt.msun
-        highest_halo_mass = 8e15 * unyt.msun
+        # highest_halo_mass = 8e15 * unyt.msun
-        bin_centers, mass_function, error = tools.create_mass_function(
+        # bin_centers, mass_function, error = tools.create_mass_function(
-            masses_200crit, lowest_halo_mass, highest_halo_mass, (100*Mpc)**3
+        #     masses_200crit, lowest_halo_mass, highest_halo_mass, (100 * Mpc) ** 3
-        )
+        # )
-        plt.loglog(bin_centers,mass_function)
-        # halos = read_velo_halos(dir)
         #
-        # print(halos["Rvir"].sort_values(ascending=False)[:4].to_numpy())
+        # name = f"{waveform} {resolution}"
 
+        # plt.loglog(bin_centers, mass_function, label=name)
+        # halos = read_halo_file(dir/"fof_output_0004.hdf5")
+        #
+        # print(halos["Masses"].sort_values(ascending=False)[:4].to_numpy())
 
+        halos = read_vr_files.read_velo_halos(dir)
+        print(halos["R_200mean"].sort_values(ascending=False)[:4].to_numpy())
+
+plt.legend()
 plt.show()

readfiles.py

@@ -75,6 +75,14 @@ def read_halo_file(file: Path) -> DataFrame:
     return df
 
 
+def read_gadget_halos(directory: Path):
+    reference_file = h5py.File(directory / "output" / "fof_subhalo_tab_019.hdf5")
+    df1 = pd.DataFrame(reference_file["Subhalo"]["SubhaloPos"], columns=["X", "Y", "Z"])
+    df2 = pd.DataFrame(reference_file["Subhalo"]["SubhaloMass"], columns=["Masses"])
+    df = df1.merge(df2, "outer", left_index=True, right_index=True)
+    return df
+
+
 def read_fof_file(path: Path):
     file = path / ""