Unverified Commit 769b6f6c authored by Akke Viitanen's avatar Akke Viitanen
Browse files

update plots

parent 48b657df

src/scripts/plots/plot_black_hole_mass_function.py

+22 −16
Original line number Diff line number Diff line
@@ -7,10 +7,13 @@ 
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")





def plot_black_hole_mass_function_local(

    catalogs,
@@ -30,7 +33,7 @@ def plot_black_hole_mass_function_local( 
        return bins, cens



    bins_bh, cens_bh = get_bins_cens(5, 12)

    fig, axes = plt.subplots(1, 1, figsize=(1 * 6.4, 1 * 4.8))

    fig, axes = plt.subplots(1, 1, figsize=(8.3 / 2, 5.85 / 2))

    axes = np.atleast_2d(axes)



    # Define the selection
@@ -101,7 +104,7 @@ def plot_black_hole_mass_function_local( 
                        x[my_select],

                        y[my_select],

                        dy[my_select],

                        label="Mock catalog, 24 deg2" + add_label,

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

                        marker=marker,

                        color="k",

                        alpha=alpha,
@@ -125,7 +128,7 @@ def plot_black_hole_mass_function_local( 
                            x[my_select],

                            y[my_select],

                            dy[my_select],

                            label=r"Mock catalog, 24 deg2 $f_\mathrm{occ}$ weighted",

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

                            marker=marker,

                            color="r",

                            alpha=alpha,
@@ -151,14 +154,14 @@ def plot_black_hole_mass_function_local( 


        # Plot Ananna+ 2022

        if zmin <= 0.25 < zmax:

            phi = ananna2022.get_phi_bh(10**bins_bh, *ananna2022.parameters["BHMF"]["All"][2][1:], h=0.70)

            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")

            phis = []



            # mcmc the errors

            for _ in range(1000):

                phi = ananna2022.get_phi_bh(

                    10**bins_bh, *ananna2022.parameters["BHMF"]["All"][2][1:], h=0.70, sample=True

                    10**bins_bh, *ananna2022.PARAMETERS["BHMF"]["All"][2][1:], h=0.70, sample=True

                )

                phis.append(phi)

            axes[0, col].fill_between(
@@ -171,19 +174,20 @@ def plot_black_hole_mass_function_local( 


        # Set the text/limits for the BHMF plot

        axes[0, col].semilogy()

        axes[0, col].legend(loc="lower left")

        # axes[0, col].legend(loc="lower left")

        axes[0, col].legend(loc="upper right")

        axes[0, col].text(

            0.90,

            0.90,

            0.10,

            0.10,

            f"${zmin:.2f} \\leq z < {zmax:.2f}$",

            va="top",

            ha="right",

            va="bottom",

            ha="left",

            transform=axes[0, col].transAxes,

        )

        axes[0, col].set_xlabel(r"$\log M_{\rm BH}$ [Msun]")

        axes[0, col].set_ylabel(r"$\phi_{\rm BH}$ [1/Mpc3/dex]")

        axes[0, col].set_xlim(6.5, 9.5)

        axes[0, col].set_ylim(1e-6, 1e-2)

        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)



    return fig, axes
@@ -198,7 +202,9 @@ FROM Truth 
""",

)

fig, axes = plot_black_hole_mass_function_local([catalog], add_ctk=True, use_occupation_fraction=False)

fig.savefig("fig/black_hole_mass_function_local_add_ctk_20250909.pdf")

# 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.savefig("fig/black_hole_mass_function_local_add_ctk_focc_20250909.pdf")
 
# 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")

src/scripts/plots/plot_flux_sigma_eta_ratio_magnitudes.py

+1 −1
Original line number Diff line number Diff line
@@ -288,7 +288,7 @@ def do_one(fig, axes, row, col, name_truth, name_flux, is_extended, plot=True): 
    fig.supxlabel(r"$m_\mathrm{truth}$ [ABmag]", va="bottom")

    if col == 0:

        # axes[row, col].set_ylabel(f"{name_flux}Flux / truth")

        mag_str = f"$\\Delta m_{{{name_flux}}}$"

        mag_str = f"$\\Delta m_\\mathrm{{{name_flux}}}$"

        axes[row, col].set_ylabel(f"{mag_str}")

src/scripts/plots/plot_number_counts_gband.py

+21 −16
Original line number Diff line number Diff line 
import sqlite3

import sys



import matplotlib as mpl

import matplotlib.pyplot as plt

import numpy as np

import pandas as pd

from astropy.table import Table

from lsst_inaf_agile import util



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



###############################################################################

# Settings

DPI = 300  # dots per inch
@@ -17,7 +21,9 @@ NCOLS = 1 # number of cols for subplot 
# SAVEFIG = "fig/number_counts_lsst_gband_20250707.pdf"  # output filename

# SAVEFIG = "fig/number_counts_lsst_gband_20250811.pdf"  # output filename

# SAVEFIG = "fig/number_counts_lsst_gband_20250930.pdf"  # output filename

SAVEFIG = "fig/number_counts_lsst_gband_20251017.pdf"  # output filename

# SAVEFIG = "fig/number_counts_lsst_gband_20251017.pdf"  # output filename

# SAVEFIG = "fig/number_counts_lsst_gband_20251024.pdf"  # output filename

SAVEFIG = sys.argv[1]  # output filename

###############################################################################
@@ -65,21 +71,19 @@ class SinglePlotter: 
        )



    def plot(self):

        """

        Plots something

        """

        """Plots something."""



        dmag = 0.25

        bins = np.arange(5, 50, dmag)



        def plot_one(ax, b):

        def plot_one(ax, b, plot_focc=False):

            print(b)



            # Plots the mock

            for label, select, ls in (

                ("AGN1", (self.catalog["is_agn"] == 1) & (self.catalog["is_optical_type2"] == 0), "-"),

                ("AGN2", (self.catalog["is_agn"] == 1) & (self.catalog["is_optical_type2"] == 1), "--"),

                ("Gal", (self.catalog["is_agn"] == 0) & (self.catalog["Z"] > 0), "-."),

                ("Galaxy", (self.catalog["is_agn"] == 0) & (self.catalog["Z"] > 0), "-."),

                ("Star", self.catalog["Z"].mask, ":"),

            ):

                from lsst_inaf_agile import mbh
@@ -98,9 +102,9 @@ class SinglePlotter: 
                    select_m = mag < my_bin

                    Ns_per_deg2.append(np.sum(select * select_m) / 24.0)

                print(label, select.sum(), np.interp(28.56, bins, Ns_per_deg2))

                ax.plot(bins, np.ma.log10(Ns_per_deg2), label=label, ls=ls, lw=3.0)

                ax.plot(bins, np.ma.log10(Ns_per_deg2), label=label, ls=ls, lw=2.0)



                if "AGN1" in label:

                if plot_focc and "AGN1" in label:

                    Ns_per_deg2 = []

                    for my_bin in bins:

                        select_m = mag < my_bin
@@ -127,17 +131,18 @@ class SinglePlotter: 
            x = limiting_magnitude_cosmos[b]

            # y = qso_number_counts_cosmos[b]

            y = qso_number_counts_cosmos_baseline_v4p0[b]

            lsst_fov = 9.6  # in deg2

            ax.plot(x, np.log10(y / lsst_fov), marker="o", fillstyle="none", color="C0", markersize=10)

            if False:

                ax.plot(x, np.log10(y / 9.6), marker="o", fillstyle="none", color="C0", markersize=10)



            # Plot another example with using 18.55 as the FoV from the COSMOS dither pattern

            # (FoV estimated with topcat on 20251017 from baseline v4.0)

            ax.plot(x, np.log10(y / 18.55), marker="^", fillstyle="none", color="C0", markersize=10)

            ax.plot(x, np.log10(y / 18.55), marker="^", fillstyle="none", color="k", markersize=10)

            print("baseline v4.0:", x, y / 18.55)



            # Plot Li+2025 (submitted) values, see their Table 1 for the g-band

            x = 26.7

            y = 349.0

            ax.plot(x, np.log10(y), marker="s", fillstyle="none", color="C0", markersize=10)

            ax.plot(x, np.log10(y), marker="s", fillstyle="none", color="k", markersize=10)



            # Plot literature values

            if True:
@@ -174,7 +179,7 @@ class SinglePlotter: 
                            marker=marker,

                            label=None,

                            color=color,

                            markersize=3,

                            markersize=5,

                            alpha=alpha,

                        )
@@ -197,7 +202,7 @@ class SinglePlotter: 
                        x, y = np.loadtxt(f).T

                        # ax.plot(x, y, '.', label=label, color=color)

                        ax.plot(

                            x, y, label=None, marker=marker, color=color, ls="none", markersize=3, alpha=alpha

                            x, y, label=None, marker=marker, color=color, ls="none", markersize=5, alpha=alpha

                        )



            # ax.text(0.10, 0.90, f"lsst-${b}$", transform=ax.transAxes, ha="left", va="top")
@@ -206,7 +211,7 @@ class SinglePlotter: 
            ax.set_xlim(9, 30)

            # ax.set_xticks(np.arange(10, 40))



            ax.set_ylim(-1.5, None)

            ax.set_ylim(-1.5, 5.5)

            # ax.set_yticks(10 ** np.arange(-10.0, 10.0))



        # plot_one(self.axes[0, 0], "u")
@@ -219,7 +224,7 @@ class SinglePlotter: 
        # self.axes.legend(loc="lower right", fontsize="small")

        self.axes.legend(loc="upper left", fontsize="small")

        self.axes.set_xlabel(r"$g$ [ABmag]")

        self.axes.set_ylabel(r"$\log N(<g)$ [1/deg$^2$]")

        self.axes.set_ylabel(r"$\log_{10} \left[ N(<g) \,/\, \mathrm{deg}^{-2} \right]$")

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



        return self.fig, self.axes

src/scripts/plots/plot_plambda_qu_sf_ctn_ctk.py

+12 −6
Original line number Diff line number Diff line
@@ -4,11 +4,12 @@ 


import sys



import matplotlib as mpl

import matplotlib.pyplot as plt

import numpy as np

from lsst_inaf_agile import zou2024



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

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





def plot_plambda_ctn_ctk(fig, ax, log_mstar, z):
@@ -44,7 +45,7 @@ def plot_plambda_ctn_ctk(fig, ax, log_mstar, z): 


    # Set the text

    # NOTE: log(Mstar/Modot) in LaTeX

    mstar_text = r"\log (M_\mathrm{star} / M_\odot)"

    mstar_text = r"\log_{10} \left( M_\mathrm{star} \,/\, M_\odot \right)"

    text = "\n".join(

        [

            f"${mstar_text} = {log_mstar:.1f}$",
@@ -61,8 +62,10 @@ def plot_plambda_ctn_ctk(fig, ax, log_mstar, z): 
if __name__ == "__main__":

    # Init

    savefig = sys.argv[1]



    log_mstars = 9.5, 10.5, 11.5

    zs = 0.5, 1.5, 2.5, 3.5

    # zs = 0.5, 1.5, 2.5, 3.5

    zs = 1.0, 3.0



    fig, axes = plt.subplots(

        len(zs),
@@ -82,11 +85,14 @@ if __name__ == "__main__": 
    axes[0, 0].set_ylim(-6, 0.5)



    # set the labels

    for ax in axes[3, :]:

    for ax in axes[len(zs) - 1, :]:

        # NOTE: xlabel was Msun instead of M_odot

        ax.set_xlabel(r"$\log(\lambda_{\rm SAR}/\mathrm{erg\,s^{-1}\,M_\odot^{-1}})$", fontsize="large")

        ax.set_xlabel(

            r"$\log_{10} \left( \lambda_{\rm SAR} \,/\, \mathrm{erg}\,\mathrm{s}^{-1}\,M_\odot^{-1} \right)$"

        )



    for ax in axes[:, 0]:

        ax.set_ylabel(r"$\log p(\lambda_{\rm SAR})$ [1/dex]", fontsize="large")

        ax.set_ylabel(r"$\log_{10} \left[ p(\lambda_{\rm SAR}) \,/\, \mathrm{dex}^{-1} \right]$")



        # fix the yticks to once per dex

        ax.set_yticks(np.arange(-6, 1))

src/scripts/plots/plot_quasar_luminosity_function.py

+15 −12
Original line number Diff line number Diff line
@@ -2,9 +2,12 @@ import sys 
from copy import deepcopy



import astropy.units as u

import matplotlib as mpl

import matplotlib.pyplot as plt

import numpy as np



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



catalog = None
@@ -27,7 +30,7 @@ def _plot_quasar_luminosity_function( 
            "magabs_mock-4400_point",

            np.arange(-28.0, -18 + 1e-6, dbin),

            r"$M_{\rm B}$ [ABmag]",

            r"$\log \Phi_{\rm B}$  [1/Mpc$^3$/mag]",

            r"$\log_{10} \left( \Phi_{\rm B} \,/\, \mathrm{Mpc}^{-3}\,\mathrm{mag}^{-1} \right)$",

            4400 * u.angstrom,

            1,

        ),
@@ -37,7 +40,7 @@ def _plot_quasar_luminosity_function( 
    from lsst_inaf_agile.shen2020 import get_phi_shen2020



    x, dx, y, dy = get_phi_shen2020(z)[idx]

    ax.errorbar(x, y, yerr=dy, marker=".", linestyle="none", label="Observed (Shen+ 2020)", color=color)

    ax.errorbar(x, y, yerr=dy, marker=".", linestyle="none", label="Shen+ 2020", color=color)



    # Plot the mock

    from lsst_inaf_agile.catalog_combined import CatalogCombined
@@ -54,7 +57,7 @@ def _plot_quasar_luminosity_function( 
            """,

        )

    select = (catalog["is_agn"] == 1) & (catalog["is_optical_type2"] == 0)

    label = r"Mock catalog, $24\,{\rm deg}^2$, AGN type1"

    label = r"Mock catalog, $24\,{\rm deg}^2$"  # , AGN type1"

    deredden = True

    color = "C0" if "B band" in band else "C4"

    ls = "dotted"
@@ -127,31 +130,31 @@ def _plot_quasar_luminosity_function( 
        ax.plot(x, np.ma.log10(y), label=label, color=color, ls=ls)

        ax.fill_between(x, np.ma.log10(y - dy), np.ma.log10(y + dy), color=color, alpha=0.20)



    # Set the text exc...

    ax.set_xlim(np.min(bins), np.max(bins))



    if show_xlabel:

        ax.set_xlabel(xlabel, fontsize="large")

    if show_ylabel:

        ax.set_ylabel(ylabel, fontsize="large")

    if show_legend:

        ax.legend(frameon=False, loc="lower right")

        ax.legend(frameon=False, loc="upper left")



    ax.text(

        0.10,

        0.90,

        f"$z \\sim {z}$",

        0.10,

        f"$z \\approx {z}$",

        fontsize="large",

        transform=ax.transAxes,

        horizontalalignment="left",

        horizontalalignment="right",

    )



    ax.set_ylim(-8.0, None)

    # ax.set_xlim(np.min(bins), np.max(bins))

    ax.set_xlim(-28.5, -21.2)

    ax.set_ylim(-8.6, -4.9)



    return fig, ax





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

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

fig, axes = plt.subplots(2, 2, sharex=True, sharey=True)

zs = 0.5, 1.5, 2.5, 3.5

savefig = sys.argv[1]

show_occupation_fraction = int(sys.argv[2])