add plot script for BHMF at higher redshifts

#!/usr/bin/env python3

# Author: Akke Viitanen

# Email: akke.viitanen@helsinki.fi

# Date: 2026-01-08 17:23:12

"""

Plots the (non-local) black hole mass function at higher redshifts.

"""

import numpy as np

from lsst_inaf_agile import kelly2013, util

from lsst_inaf_agile.catalog_combined import CatalogCombined

def plot_bhmf_catalog(fig, ax, catalog, redshift, dlogMBH=0.40, dz=0.30):

    def get_bins_cens(amin, amax, dbin=dlogMBH):

        bins = np.arange(amin, amax, dbin)

        cens = 0.50 * (bins[:-1] + bins[1:])

        return bins, cens

    bins_bh, cens_bh = get_bins_cens(5, 12)

    zmin = redshift - dz / 2

    zmax = redshift + dz / 2

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

        "MBH",

        bins=bins_bh,

        zmin=zmin,

        zmax=zmax,

        select=(

            (zmin <= catalog["Z"])

            & (catalog["Z"] < zmax)

            & (catalog["is_agn"] == 1)

            & (catalog["is_optical_type2"] == 0)

        ),

        use_occupation_fraction=False,

    )

    select = y > 0.0

    # convert to log

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

    ax.errorbar(

        x[select],

        y[select],

        (ylo[select], yhi[select]),

        label=r"Mock catalog, $24\,\mathrm{deg}^2$",

        marker=".",

        color="k",

        alpha=0.80,

        linestyle="none",

    )

    return fig, ax

def main():

    fig, axes = kelly2013._get_figax()

    kelly2013._plot_bhmf(fig, axes, kelly2013.REDSHIFTS)

    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

        """,

    )

    # Shen has deltaz of 0.20, 0.35, 0.55. Choose deltaz = 0.20 for now as the

    # binsize

    for redshift, ax in zip(kelly2013.REDSHIFTS, axes.flatten(), strict=False):

        plot_bhmf_catalog(fig, ax, catalog, redshift)

    fig.savefig("tmp.pdf")

if __name__ == "__main__":

    main()