re-implements xray luminosity function plot (3c1c4831) · Commits · Akke Esa Tapio Viitanen / lsst_inaf_agile

src/lsst_inaf_agile/catalog_agn.py

+28 −29

Original line number	Diff line number	Diff line
		@@ -48,7 +48,6 @@ except FileNotFoundError:


		class CatalogAGN:

		"""
		AGN catalog class
		"""
		@@ -100,40 +99,40 @@ class CatalogAGN:
		self.catalog = self.get_catalog()

		def __getitem__(self, key):
		"""
		Return AGN catalog column.

		Parameters
		----------
		key: str
		name of column

		Returns
		-------
		value of the column
		"""
		if key not in self.catalog.dtype.names:
		return self.catalog_galaxy[key]
		return self.catalog[key]

		def _get_select(self, redshift=None, mstar=None, luminosity=None, t="all"):
		select = np.ones_like(self["Z"], dtype=bool)

		# Perform the selection in redshift and Mstar
		if redshift is not None:
		if redshift not in self.SZ:
		self.SZ[redshift] = (redshift <= self["Z"]) * (self["Z"] < redshift + self.dz)
		select *= self.SZ[redshift]

		if mstar is not None:
		if mstar not in self.SM:
		self.SM[mstar] = (mstar <= self["M"]) * (self["M"] < mstar + self.dm)
		select *= self.SM[mstar]

		# Perform the selection in lx if requested
		if luminosity is not None:
		if luminosity not in self.SL:
		self.SL[luminosity] = (luminosity <= self.catalog["log_LX_2_10"]) * (
		self.catalog["log_LX_2_10"] < luminosity + self.dl
		)
		select *= self.SL[luminosity]

		# Perform the selection in type
		if t != "all":
		select *= self.catalog["PASSIVE"] == (t == "quiescent")

		return select

		@staticmethod
		def get_columns(bands, rfbands, is_public=False):
		"""
		Return the names, types, descriptions, and units of each column

		Parameters
		----------
		bands: list
		egg-style list of apparent magnitude bands
		rfbands: list
		egg-style list of absolute magnitude bands
		is_public: bool
		return a public column listing instead of an internal one

		Returns
		-------
		List of 4-tuples of column names, types, descriptions, and units
		"""
		return (
		[
		("ID", np.int64, "unique ID", ""),

src/lsst_inaf_agile/catalog_combined.py

+15 −2

Original line number	Diff line number	Diff line
		@@ -9,6 +9,7 @@ Create Galaxy+AGN mocks for the LSST Italian AGN in-kind contribution

		import logging
		import os
		import re

		import fitsio
		import numpy as np
		@@ -18,7 +19,9 @@ logger = logging.getLogger(__name__)


		class CatalogCombined:
		def __init__(self, dirname: str, catalog_galaxy, catalog_agn, catalog_star, catalog_binary):
		def __init__(
		self, dirname: str, catalog_galaxy=None, catalog_agn=None, catalog_star=None, catalog_binary=None
		):
		self.dirname = dirname
		self.filename = f"{dirname}/catalog.fits"

		@@ -27,6 +30,11 @@ class CatalogCombined:
		self.catalog_star = catalog_star
		self.catalog_binary = catalog_binary

		# NOTE: short-circuit for an existing truth catalog
		if os.path.exists(self.filename):
		self.catalog_combined = fitsio.read(self.filename)
		return

		self.catalog_combined = self.get_catalog_combined()
		self.catalog_combined = self.postprocess()
		self.write()
		@@ -249,5 +257,10 @@ class CatalogCombined:
		os.system(f"convertReferenceCatalog {dirname} etc/config_reference_catalog.py {filename}")
		logger.info(f"Wrote reference catalog {filename}")

		def get_area(self):
		cmd = fitsio.read(f"{self.dirname}/egg.fits", columns="CMD")[0]
		area_string = re.findall("area=([0-9.]+)", cmd)[0]
		return float(area_string)

		def __getitem__(self, key):
		return self.catalog_combined[key]

src/lsst_inaf_agile/util.py

+1 −1

Original line number	Diff line number	Diff line
		@@ -20,7 +20,7 @@ from astropy.time import Time
		from igm_absorption import my_get_IGM_absorption
		from scipy.stats import binned_statistic

		logging.basicConfig(format="%(asctime)s - %(name)s - %(levelname)s - %(message)s", level=logging.DEBUG)
		logging.basicConfig(format="%(asctime)s - %(name)s - %(levelname)s - %(message)s", level=logging.INFO)
		logger = logging.getLogger(__name__)

		ROOT = "/staff2/viitanen/agile"

src/scripts/plots/plot_xray_luminosity_function.py

0 → 100644

+232 −0

Original line number	Diff line number	Diff line
		#!/usr/bin/env python3
		# Author: Akke Viitanen
		# Email: akke.viitanen@helsinki.fi

		"""
		Plots the X-ray luminosity function as a function of redshift. The X-ray
		luminosity function of the input catalog can be plotted along with relevant
		literature references.
		"""

		import glob
		import re
		from argparse import ArgumentParser

		import matplotlib.pyplot as plt
		import numpy as np
		import util
		import zou2024


		def main(dirname, filename_savefig):
		plot_catalog = True
		plot_catalog_focc = True

		plot_miyaji15 = False
		plot_ueda14 = False
		plot_shen20 = True
		plot_zou24 = False

		dbin = 0.20
		bins = np.arange(42, 46 + dbin / 2, dbin)
		cens = 0.5 * (bins[:-1] + bins[1:])
		log_lx = cens

		miy15 = np.loadtxt("data/miyaji2015/xlf_miyaji2015_table4.dat")
		zcen = miy15[:, 2]

		dz = 0.20
		zcens = 0.5, 1.5, 2.5, 3.5
		zset = [(z - dz, z + dz, z) for z in zcens]

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

		def plot_one(ax, zmin, zmax, zcen):
		print(ax, zmin, zmax, zcen)

		if plot_miyaji15:
		# Plot Miy+15 data points
		# Find closest miyaji
		zbest = miy15[np.argmin(np.abs(miy15[:, 2] - zcen)), 2]
		select = miy15[:, 2] == zbest
		y, dy1, dy2 = util.get_log_y_lo_hi(miy15[select, -2], miy15[select, -1])
		ax.errorbar(
		miy15[select, 5],
		y,
		(dy1, dy2),
		linestyle="none",
		marker=".",
		label="Miyaji+ 2015",
		color="C2",
		fillstyle="none",
		markersize=12,
		)

		if plot_ueda14:
		# Plot Ueda+2014 data points
		filename_ueda = sorted(glob.glob("opt/quasarlf/data/ueda2014/hard*.dat"))
		zmin_ueda = np.array([float(re.findall("z(.)-(.).dat", f)[0][0]) for f in filename_ueda])
		zmax_ueda = np.array([float(re.findall("z(.)-(.).dat", f)[0][1]) for f in filename_ueda])
		zcen_ueda = (zmin_ueda * zmax_ueda) ** 0.5
		idx = np.argmin(np.abs(zcen - zcen_ueda))
		x, _y, _y1, _y2 = np.loadtxt(filename_ueda[idx]).T
		y = np.log10(_y)
		y_lo = np.log10(_y1)
		y_hi = np.log10(_y2)
		ax.errorbar(
		x,
		y,
		(y - y_lo, y_hi - y),
		label="Ueda+ 2014",
		marker=".",
		linestyle="none",
		color="C3",
		fillstyle="none",
		markersize=12,
		)

		if plot_zou24:
		raise NotImplementedError("Re-implement the SMF module to enable Zou+2024")

		# Plot Zou+2024 model
		# Initialize the SMF
		import smf

		my_smf = smf.StellarMassFunctionCOSMOS2020AGN(type_plambda="zou+2024")
		kwargs = {
		"log_lx": log_lx,
		"fun_phi_star": lambda m, z, t: my_smf.get_stellar_mass_function(
		10**m, z, t, loglambda_min=32.0
		),
		"log_mstar_min": 9.0,
		"log_mstar_max": 12.5,
		}

		xlf1 = np.zeros_like(log_lx)
		xlf1 += zou2024.get_xray_luminosity_function(z=zmin, t="star-forming", **kwargs)
		xlf1 += zou2024.get_xray_luminosity_function(z=zmin, t="quiescent", **kwargs)

		xlf2 = np.zeros_like(log_lx)
		xlf2 += zou2024.get_xray_luminosity_function(z=zmax, t="star-forming", **kwargs)
		xlf2 += zou2024.get_xray_luminosity_function(z=zmax, t="quiescent", **kwargs)

		ax.fill_between(
		log_lx,
		np.log10(np.minimum(xlf1, xlf2)),
		np.log10(np.maximum(xlf1, xlf2)),
		alpha=0.20,
		color="k",
		label="Zou+ 2024",
		)

		if plot_catalog:
		# Plot the Xray luminosity function
		from catalog_combined import CatalogCombined
		from util import get_key_function

		catalog = CatalogCombined(dirname)
		values = (catalog["log_LX_2_10"] > 42) & (zmin <= catalog["Z"]) & (catalog["Z"] < zmax)
		x, dx, y, dy = get_key_function(
		bins,
		catalog["log_LX_2_10"],
		values=values,
		area_deg2=catalog.get_area(),
		zmin=zmin,
		zmax=zmax,
		)
		select = y - dy > 0.0
		x, dx, y, dy = map(lambda a: a[select], [x, dx, y, dy])
		log_y_lo = np.log10(y - dy)
		log_y_hi = np.log10(y + dy)

		ax.fill_between(
		x,
		log_y_lo,
		log_y_hi,
		label=r"Mock catalog, $24\,{\rm deg}^2$",
		alpha=0.20,
		color="k",
		)

		if plot_catalog_focc:
		# Plot the Xray luminosity function weighted by the black hole
		# occupation fraction
		from catalog_combined import CatalogCombined
		from util import get_key_function

		catalog = CatalogCombined(dirname)
		values = (catalog["log_LX_2_10"] > 42) & (zmin <= catalog["Z"]) & (catalog["Z"] < zmax)
		values = values.astype(np.float64) * catalog["occupation_fraction"]
		x, dx, y, dy = get_key_function(
		bins,
		catalog["log_LX_2_10"],
		values=values,
		area_deg2=catalog.get_area(),
		zmin=zmin,
		zmax=zmax,
		)
		select = y - dy > 0.0
		x, dx, y, dy = map(lambda a: a[select], [x, dx, y, dy])
		log_y_lo = np.log10(y - dy)
		log_y_hi = np.log10(y + dy)

		ax.fill_between(
		x,
		log_y_lo,
		log_y_hi,
		label=r"Mock catalog, $24\,{\rm deg}^2$, $f_{\rm occ}$-weighted",
		alpha=0.20,
		color="r",
		)

		if plot_shen20:
		# Plot Shen+ 2020
		from shen2020 import get_phi_shen2020

		# NOTE: 3 is for 2-10 keV xray
		x, dx, y, dy = get_phi_shen2020(zcen)[3]
		ax.errorbar(
		x,
		y,
		yerr=dy,
		marker=".",
		linestyle="none",
		label="Observed (Shen+ 2020)",
		color="green",
		)

		# Set the text etc..
		ax.set_xlim(42.0 - 0.20, 46.0 + 0.20)
		ax.set_ylim(-8.5 + 0.10, -3.0 - 0.10)
		ax.text(
		0.9,
		0.9,
		f"$z \\sim {zcen:.1f}$",
		transform=ax.transAxes,
		horizontalalignment="right",
		fontsize="large",
		)

		if ax == axes[0, 0]:
		ax.legend(loc="lower left")

		if ax in axes[-1, :]:
		ax.set_xlabel(r"$\log L_{\rm X}$ [erg/s]", fontsize="large")

		if ax in axes[:, 0]:
		ax.set_ylabel(r"$\log \Phi_{\rm X}$ [1/Mpc$^3$/dex]", fontsize="large")

		for ax, (zmin, zmax, zcen) in zip(axes.flatten(), zset):
		plot_one(ax, zmin, zmax, zcen)

		fig.tight_layout(pad=0.20)
		fig.savefig(filename_savefig)


		if __name__ == "__main__":
		parser = ArgumentParser()
		parser.add_argument("dirname", help="Directory name containing the truth catalog")
		parser.add_argument("filename_savefig")
		args = parser.parse_args()
		main(args.dirname, args.filename_savefig)

tests/lsst_inaf_agile/test_catalog_combined.py

+19 −3

Original line number	Diff line number	Diff line
		@@ -4,12 +4,13 @@
		# Date: 2025-07-03 19:33:15

		"""
		Tests the galaxy catalog
		Tests the combined catalog
		"""

		import logging
		import os

		import numpy as np
		from catalog_agn import CatalogAGN
		from catalog_combined import CatalogCombined
		from catalog_galaxy import CatalogGalaxy
		@@ -21,8 +22,11 @@ from merloni2014 import Merloni2014
		logger = logging.getLogger(__name__)


		def main():
		dirname = "data/catalog/test"
		def get_catalog_combined(dirname):
		"""
		Returns a test CatalogCombined object
		"""

		os.makedirs(dirname, exist_ok=True)

		# Initialize EGG
		@@ -45,6 +49,18 @@ def main():
		# Create the combined catalog
		catalog_combined = CatalogCombined(dirname, catalog_galaxy, catalog_agn, catalog_star, catalog_binary)

		assert np.isclose(catalog_combined.get_area(), 0.10)

		return catalog_combined


		def main():
		# Setup the dirname
		dirname = "data/catalog/test"

		# Get the combined catalog
		catalog_combined = get_catalog_combined(dirname)

		# Create the reference catalog
		maglim = 24
		selection_band = "lsst-r"