From c7b01f3062fe4b86aa78c83e3fcdeee6300de1d3 Mon Sep 17 00:00:00 2001 From: Lukas Winkler Date: Thu, 15 Dec 2022 16:44:14 +0100 Subject: [PATCH] fix HMF --- halo_mass_functions.py | 52 +++++++++++++++++++++++++++++++++--------- hmf_test.py | 40 ++++++++++++++++++-------------- readfiles.py | 8 +++++++ 3 files changed, 72 insertions(+), 28 deletions(-) diff --git a/halo_mass_functions.py b/halo_mass_functions.py index 7e5cf27..319e06e 100644 --- a/halo_mass_functions.py +++ b/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) - - halo_masses: np.ndarray = halos["Mvir"].to_numpy() * 1e10 * 0.67742 + if "gadget" in waveform: + halos = read_gadget_halos(dir) + 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_ylabel("dndlnM") + ax.set_xlabel("Mass ($h^{-1} M_{\odot}$)") + 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) diff --git a/hmf_test.py b/hmf_test.py index c69dbb9..f5dc6b7 100644 --- a/hmf_test.py +++ b/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") - print(data) - masses_200crit = data.masses.mass_200crit - masses_200crit.convert_to_units("msun") - lowest_halo_mass = 1e9 * unyt.msun - highest_halo_mass = 8e15 * unyt.msun - bin_centers, mass_function, error = tools.create_mass_function( - masses_200crit, lowest_halo_mass, highest_halo_mass, (100*Mpc)**3 - ) - plt.loglog(bin_centers,mass_function) - # halos = read_velo_halos(dir) + # data = load(str(dir) + "/vroutput.properties.0") + # print(data) + # masses_200crit = data.masses.mass_200crit + # masses_200crit.convert_to_units("msun") + # lowest_halo_mass = 1e9 * unyt.msun + # highest_halo_mass = 8e15 * unyt.msun + # bin_centers, mass_function, error = tools.create_mass_function( + # masses_200crit, lowest_halo_mass, highest_halo_mass, (100 * Mpc) ** 3 + # ) # - # 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() diff --git a/readfiles.py b/readfiles.py index 9fd31cd..4d743cd 100644 --- a/readfiles.py +++ b/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 / ""