add schulze bhmf plot (f57ea368) · Commits · Akke Esa Tapio Viitanen / lsst_inaf_agile

src/scripts/plots/plot_black_hole_mass_function.py

+217 −147

Original line number	Diff line number	Diff line
		@@ -7,25 +7,24 @@
		Plots the black hole mass function
		"""

		import matplotlib as mpl
		import matplotlib.pyplot as plt
		import numpy as np
		from lsst_inaf_agile.catalog_combined import CatalogCombined

		mpl.style.use("etc/agile.mplstyle")
		plt.style.use("etc/agile.mplstyle")


		def plot_black_hole_mass_function_local(
		catalogs,
		bins=None,
		add_ctk=False,
		use_occupation_fraction=False,
		catalog, bins=None, add_ctk=False, use_occupation_fraction=False, reference="ananna"
		):
		"""
		bins = dict of key: (bin_lo, bin_hi) pairs to plot the BHMF in bins
		"""
		if reference == "ananna":
		zmin, zmax = 0.21, 0.30
		elif reference == "schulze":
		zmin, zmax = 1.00, 2.00
		else:
		raise ValueError

		from lsst_inaf_agile import ananna2022, ueda2014
		from lsst_inaf_agile import ananna2022, schulze2015, ueda2014

		def get_bins_cens(amin, amax, dbin=0.20):
		bins = np.arange(amin, amax, dbin)
		@@ -37,9 +36,25 @@ def plot_black_hole_mass_function_local(
		axes = np.atleast_2d(axes)

		# Define the selection
		def get_select(log_mbh, log_lambda_edd, fix_for_ctk=False, z=None, log_LX_2_10=None):
		def get_select(
		log_mbh, log_lambda_edd, fix_for_ctk=False, z=None, log_LX_2_10=None, log_lambda_edd_limit=-3
		):
		# Ananna-like sample selection
		select = (log_mbh >= 6.5) * (log_mbh < 10.5) * (log_lambda_edd >= -3.0) * (log_lambda_edd < 1.0)
		if reference == "ananna":
		select = (
		(log_mbh >= 6.5)
		& (log_mbh < 10.5)
		& (log_lambda_edd >= log_lambda_edd_limit)
		& (log_lambda_edd < 1.0)
		)
		elif reference == "schulze":
		select = (
		(log_mbh >= 7.0)
		& (log_mbh < 10.5)
		& (log_lambda_edd >= -2)
		& (log_lambda_edd < 1.0)
		& (catalog["is_optical_type2"] == 0)
		)

		# NOTE: fix for CTK
		if fix_for_ctk:
		@@ -54,13 +69,11 @@ def plot_black_hole_mass_function_local(

		return select

		for col, (zmin, zmax) in enumerate([(0.21, 0.30)]):
		# Plot the mock
		for catalog in catalogs:
		zvec = catalog["Z"][np.isfinite(catalog["Z"])]
		if zvec.min() < zmin:
		zmin = zvec.min()
		print(zmin, zmax, catalog)
		## Plot the mock
		# zvec = catalog["Z"][np.isfinite(catalog["Z"])]
		# if zvec.min() < zmin:
		# zmin = zvec.min()
		# print(zmin, zmax, catalog)

		def my_plot(
		catalog,
		@@ -70,10 +83,11 @@ def plot_black_hole_mass_function_local(
		my_mbh,
		my_log_lambda_edd,
		label,
		color,
		color=None,
		select_bin=None,
		fix_for_ctk=True,
		add_label="",
		log_lambda_edd_limit=-3,
		):
		if select_bin is None:
		select_bin = np.ones_like(my_mbh, dtype=bool)
		@@ -82,14 +96,19 @@ def plot_black_hole_mass_function_local(
		(0.5, "all AGN", "s", np.ones_like(catalog["is_optical_type2"], dtype=bool)),
		):
		select = get_select(
		my_mbh, my_log_lambda_edd, fix_for_ctk, catalog["log_LX_2_10"], catalog["Z"]
		my_mbh,
		my_log_lambda_edd,
		fix_for_ctk,
		catalog["log_LX_2_10"],
		catalog["Z"],
		log_lambda_edd_limit,
		)
		select *= select_type
		if select_bin is not None:
		select *= select_bin

		# plot bhmf
		x, dx, y, dy = catalog.get_luminosity_function(
		x, _, y, dy = catalog.get_luminosity_function(
		"",
		bins=bins_bh,
		zmin=zmin,
		@@ -100,20 +119,20 @@ def plot_black_hole_mass_function_local(
		)
		my_select = y > 0.0

		axes[0, col].errorbar(
		axes[0, 0].errorbar(
		x[my_select],
		y[my_select],
		dy[my_select],
		label=r"Mock catalog, $24\,\mathrm{deg}^2$" + add_label,
		marker=marker,
		color="k",
		color="k" if color is None else color,
		alpha=alpha,
		linestyle="none",
		)

		# plot bhmf
		if use_occupation_fraction:
		x, dx, y, dy = catalog.get_luminosity_function(
		x, _, y, dy = catalog.get_luminosity_function(
		"",
		bins=bins_bh,
		zmin=zmin,
		@@ -124,7 +143,7 @@ def plot_black_hole_mass_function_local(
		)
		my_select = y > 0.0

		axes[0, col].errorbar(
		axes[0, 0].errorbar(
		x[my_select],
		y[my_select],
		dy[my_select],
		@@ -135,7 +154,6 @@ def plot_black_hole_mass_function_local(
		linestyle="none",
		)

		# Handle case add_ctk
		my_mbh = catalog["MBH"]

		# NOTE: Dan's email on 2025/05: need to be consistent with the
		@@ -144,18 +162,35 @@ def plot_black_hole_mass_function_local(
		# p(lambda) with the extrapolation. To compare with Ananna, we
		# estimated Eddington ratio using the same kbol, and the MBH
		# from the continuity equation for each object individually
		my_log_lambda = np.log10(25) + catalog["log_LX_2_10"] - np.log10(1.26e38) - my_mbh
		my_log_lambda0 = np.log10(25) + np.abs(catalog["log_LX_2_10"]) - np.log10(1.26e38) - my_mbh
		my_log_lambda1 = catalog["log_lambda_SAR"] - 34
		my_log_lambda2 = catalog["log_lambda_SAR"] - 33
		my_log_lambda3 = catalog["log_lambda_SAR"] - 32

		if bins is None:
		for log_lambda, color, add_label in [
		(my_log_lambda0, None, None),
		(my_log_lambda1, "C0", " ledd=lsar-34"),
		(my_log_lambda2, "C1", " ledd=lsar-33"),
		(my_log_lambda3, "C2", " ledd=lsar-32"),
		]:
		my_plot(
		catalog, col, zmin, zmax, my_mbh, my_log_lambda, None, color="C0", fix_for_ctk=not add_ctk
		catalog,
		0,
		zmin,
		zmax,
		my_mbh,
		log_lambda,
		None,
		color=color,
		fix_for_ctk=not add_ctk,
		add_label=add_label,
		)
		continue

		# Plot Ananna+ 2022
		if zmin <= 0.25 < zmax:
		if reference == "ananna":
		phi = ananna2022.get_phi_bh(10*bins_bh, ananna2022.PARAMETERS["BHMF"]["All"][2][1:], h=0.70)
		axes[0, col].plot(bins_bh, phi, color="green", ls="dotted", label="Ananna+ 2022")
		axes[0, 0].plot(bins_bh, phi, color="green", ls="dotted", label="Ananna+ 2022")
		phis = []

		# mcmc the errors
		@@ -164,7 +199,7 @@ def plot_black_hole_mass_function_local(
		10*bins_bh, ananna2022.PARAMETERS["BHMF"]["All"][2][1:], h=0.70, sample=True
		)
		phis.append(phi)
		axes[0, col].fill_between(
		axes[0, 0].fill_between(
		bins_bh,
		np.quantile(phis, 0.16, axis=0),
		np.quantile(phis, 0.84, axis=0),
		@@ -172,39 +207,74 @@ def plot_black_hole_mass_function_local(
		color="green",
		)

		if reference == "schulze":
		from lsst_inaf_agile import schulze2015

		schulze1 = 10 ** schulze2015.get_log_Phi_bh(bins_bh, zmin)
		schulze2 = 10 ** schulze2015.get_log_Phi_bh(bins_bh, zmax)
		axes[0, 0].fill_between(bins_bh, schulze1, schulze2, alpha=0.20, color="C0", label="Schulze+15")

		# Set the text/limits for the BHMF plot
		axes[0, col].semilogy()
		# axes[0, col].legend(loc="lower left")
		axes[0, col].legend(loc="upper right")
		axes[0, col].text(
		axes[0, 0].semilogy()
		axes[0, 0].legend(loc="upper right")
		axes[0, 0].text(
		0.10,
		0.10,
		f"${zmin:.2f} \\leq z < {zmax:.2f}$",
		va="bottom",
		ha="left",
		transform=axes[0, col].transAxes,
		transform=axes[0, 0].transAxes,
		)
		axes[0, col].set_xlabel(r"$\log_{10} \left( M_{\rm BH} \,/\, M_\odot \right)$")
		axes[0, col].set_ylabel(r"$\Phi_{\rm BH}$ [$\mathrm{Mpc}^{-3}\,\mathrm{dex}^{-1}$]")
		axes[0, col].set_xlim(6.3, 10)
		axes[0, col].set_ylim(1e-6, 5e-3)
		axes[0, 0].set_xlabel(r"$\log_{10} \left( M_{\rm BH} \,/\, M_\odot \right)$")
		axes[0, 0].set_ylabel(r"$\Phi_{\rm BH}$ [$\mathrm{Mpc}^{-3}\,\mathrm{dex}^{-1}$]")

		if reference == "ananna":
		axes[0, 0].set_xlim(6.3, 10)
		axes[0, 0].set_ylim(1e-6, 5e-3)
		else:
		axes[0, 0].set_xlim(7.0, 10.5)
		axes[0, 0].set_ylim(1e-6, 1e-2)

		return fig, axes


		columns = (
		"Z",
		"M",
		"log_LX_2_10",
		"MBH",
		"log_L_bol",
		"is_optical_type2",
		"is_agn",
		"is_agn_ctn",
		"is_agn_ctk",
		"log_lambda_SAR",
		"log_lambda_Edd",
		)
		catalog = CatalogCombined(
		"data/catalog/dr1_new_new",
		sql_query="""
		SELECT Z, M, log_LX_2_10, MBH, log_L_bol,
		is_optical_type2, is_agn, is_agn_ctn, is_agn_ctk,
		log_lambda_SAR
		FROM Truth
		""",
		)
		fig, axes = plot_black_hole_mass_function_local([catalog], add_ctk=True, use_occupation_fraction=False)
		"data/catalog/dr1_new_new", columns=columns, sql_query=f"""SELECT {", ".join(columns)} FROM Truth"""
		)

		fig, axes = plot_black_hole_mass_function_local(
		catalog, add_ctk=True, use_occupation_fraction=False, reference="ananna"
		)
		# fig.savefig("fig/black_hole_mass_function_local_add_ctk_20251006.pdf")
		fig.savefig("fig/black_hole_mass_function_local_add_ctk_20251024.pdf")

		fig, axes = plot_black_hole_mass_function_local([catalog], add_ctk=True, use_occupation_fraction=True)
		fig, axes = plot_black_hole_mass_function_local(
		catalog, add_ctk=True, use_occupation_fraction=False, reference="ananna"
		)
		fig.savefig("fig/black_hole_mass_function_ananna_20260226.pdf")

		fig, axes = plot_black_hole_mass_function_local(
		catalog, add_ctk=True, use_occupation_fraction=False, reference="schulze"
		)
		fig.savefig("fig/black_hole_mass_function_schulze_20260226.pdf")
		plt.show()

		quit()


		fig, axes = plot_black_hole_mass_function_local(catalog, add_ctk=True, use_occupation_fraction=True)
		# fig.savefig("fig/black_hole_mass_function_local_add_ctk_focc_20251006.pdf")
		fig.savefig("fig/black_hole_mass_function_local_add_ctk_focc_20251024.pdf")