Add schulze2015 plot script

fig/black_hole_mass_function_redshift.pdf: src/scripts/plots/plot_black_hole_mass_function_redshift.py

	python3 $< $@

fig/black_hole_mass_function_redshift_schulze2015_shen_kelly2013.pdf: src/scripts/plots/plot_black_hole_mass_function_schulze2015.py

	python3 $< $@

#!/usr/bin/env python3

# Author: Akke Viitanen

# Email: akke.viitanen@helsinki.fi

# Date: 2026-01-13 17:25:36

"""

Plot the Schulze black hole mass function

"""

import sys

import matplotlib.pyplot as plt

import numpy as np

from lsst_inaf_agile import kelly2013, schulze2015, util

from lsst_inaf_agile.catalog_combined import CatalogCombined

def main(savefig):

    columns = (

        "Z",

        "M",

        "log_lambda_SAR",

        "log_lambda_Edd",

        "MBH",

        "is_agn_ctn",

        "is_agn_ctk",

        "is_agn",

        "is_optical_type2",

    )

    catalog_combined = CatalogCombined(

        "data/catalog/dr1_new_new/", sql_query=f"SELECT {', '.join(columns)} FROM Truth"

    )

    log_m_bh = np.linspace(7, 10, 31)

    # Redshift bins from S+15 Fig. 8.

    redshifts = [

        (1.1, 1.3),

        (1.3, 1.5),

        (1.5, 1.7),

        (1.7, 1.9),

        (1.9, 2.1),

    ]

    fig, axes = plt.subplots(2, 3, sharex=True, sharey=True, figsize=(3 * 6.4, 2 * 4.8), dpi=150)

    axes[-1, -1].remove()

    for ax, (redshift1, redshift2) in zip(axes.flatten(), redshifts, strict=False):

        # Plot Schulze

        log_Phi_bh1 = schulze2015.get_log_Phi_bh(log_m_bh, redshift1)

        log_Phi_bh2 = schulze2015.get_log_Phi_bh(log_m_bh, redshift2)

        ax.fill_between(log_m_bh, log_Phi_bh1, log_Phi_bh2, label="S+15", color="k", alpha=0.4)

        # Plot Kelly

        x, y, _, _ = kelly2013.get_bhmf(0.5 * (redshift1 + redshift2))

        ax.plot(x, y, label="K&S13", color="r")

        def plot_one(ax, redshift1, redshift2, select, label, values=None):

            """

            Plot one selection of the mock catalog with a given label.

            """

            x, dx, y, dy = catalog_combined.get_luminosity_function(

                "MBH", log_m_bh, zmin=redshift1, zmax=redshift2, select=select, nmin=10, values=values

            )

            y, ylo, yhi = util.get_log_y_lo_hi(y, dy)

            ax.errorbar(x, y, xerr=dx, yerr=(ylo, yhi), label=label, linestyle="none", marker=".")

        # Select ALL AGN

        plot_one(ax, redshift1, redshift2, catalog_combined["is_agn"], "is_agn")

        # Select CTN AGN

        plot_one(ax, redshift1, redshift2, catalog_combined["is_agn_ctn"], "is_agn_ctn")

        # Select type1 AGN

        plot_one(

            ax,

            redshift1,

            redshift2,

            catalog_combined["is_agn_ctn"] & (catalog_combined["is_optical_type2"] == 0),

            "type1 AGN",

        )

        # Do stupid MBH test

        log_mbh = catalog_combined["M"] - np.log10(500)

        plot_one(

            ax,

            redshift1,

            redshift2,

            catalog_combined["is_agn_ctn"] & (catalog_combined["is_optical_type2"] == 0),

            "type1 AGN MBH'",

            values=log_mbh,

        )

        # Plot the redshift text

        ax.text(

            0.9, 0.9, f"${redshift1:.1f} < z < {redshift2:.1f}$", transform=ax.transAxes, ha="right", va="top"

        )

        # Set the limits and the legend

        ax.set_xlim(6.8, 10.2)

        ax.set_ylim(-9, -3)

        ax.legend()

    fig.supxlabel(r"$\log \left( M_\mathrm{BH} / M_\odot \right)$")

    fig.supylabel(r"$\log \left( \Phi_\mathrm{BH} / \mathrm{Mpc}^{-3}\,\mathrm{dex}^{-1} \right)$")

    fig.tight_layout()

    fig.savefig(savefig, bbox_inches="tight")

    return 0

if __name__ == "__main__":

    main(sys.argv[1])