implement occupation_fraction in XLF plot (1ee099c9) · Commits · Akke Esa Tapio Viitanen / lsst_inaf_agile

src/lsst_inaf_agile/catalog_combined.py

+39 −9

Original line number	Diff line number	Diff line
		@@ -5,6 +5,7 @@

		"""Module for the combined catalog of AGN, galaxies, and stars."""

		import glob
		import logging
		import os
		import re
		@@ -17,7 +18,12 @@ logger = logging.getLogger(__name__)


		class CatalogCombined:
		"""Combined catalog class of AGN, galaxies, and stars.
		"""
		Combined catalog class of AGN, galaxies, and stars.

		If the catalog does not exist, then create it. If the catalog exists as a
		FITS file, try to load it. Finally, if there is a database containing the
		catalog, return that one instead.

		Parameters
		----------
		@@ -47,7 +53,13 @@ class CatalogCombined:
		"""

		def __init__(
		self, dirname: str, catalog_galaxy=None, catalog_agn=None, catalog_star=None, catalog_binary=None
		self,
		dirname: str,
		catalog_galaxy=None,
		catalog_agn=None,
		catalog_star=None,
		catalog_binary=None,
		sql_query=None,
		):
		self.dirname = dirname
		self.catalog_galaxy = catalog_galaxy
		@@ -57,9 +69,25 @@ class CatalogCombined:

		# NOTE: short-circuit for an existing truth catalog
		if os.path.exists(f := self.get_filename()):
		logger.info(f"Found catalog FITS file {f}")
		self.catalog_combined = fitsio.read(f)
		return

		# NOTE: short-circuit for an existing database truth catalog
		if filenames := glob.glob(f"{self.dirname}/db/**/master.db", recursive=True):
		import sqlite3

		import pandas as pd
		from astropy.table import Table

		filename = sorted(filenames)[-1]
		logger.info(f"Found database file {filename}")
		with sqlite3.connect(filename) as con:
		df = pd.read_sql_query(sql_query, con)
		self.catalog_combined = Table.from_pandas(df)
		return

		# Create the catalog
		self.catalog_combined = self.get_catalog_combined()
		self.catalog_combined = self.postprocess()
		self.write()
		@@ -344,15 +372,17 @@ class CatalogCombined:
		logger.info(f"Wrote reference catalog {filename}")

		def get_area(self):
		"""
		Return the EGG catalog area in deg2
		"""
		"""Return the EGG catalog area in deg2."""
		try:
		cmd = fitsio.read(f"{self.dirname}/egg.fits", columns="CMD")[0]
		area_string = re.findall("area=([0-9.]+)", cmd)[0]
		return float(area_string)
		except OSError:
		logger.warning("Could not find egg.fits, defaulting to DR1 area of 24deg2")
		return 24.0

		def __getitem__(self, key):
		"""
		Return combined catalog value corresponding to 'key'
		Return combined catalog value corresponding to 'key'.
		"""
		return self.catalog_combined[key]

src/lsst_inaf_agile/mbh.py

+10 −5

Original line number	Diff line number	Diff line
		@@ -302,15 +302,20 @@ def get_log_mbh(log_mstar, reference="Shankar+16", z=0.0):


		def get_occupation_fraction(logMstar):
		"""
		Returns the black hole occupation fraction 'f_occ' as a function of
		logMstar. The occupation fraction is the percentage of host galaxies
		that are expected to host SMBH, which could be <1 for dwarf galaxies
		(Miller+15, Burke+25, Zou+).
		"""Returns the black hole occupation fraction.

		The occupation fraction is the percentage of host galaxies that are
		expected to host SMBH, which could be <1 for dwarf galaxies (Miller+15,
		Burke+25, Zou+).

		The updated AGN duty cycle can be estimated using:

		is_agn_new = occupation_fraction(logMstar) * is_agn_old

		Parameters
		----------
		logMstar: float or array_like
		Log10 of host galaxy stellar mass in Msun.
		"""

		table = Table.read("data/from_zou/table2.csv")

src/scripts/plots/plot_xray_luminosity_function.py

+13 −3

Original line number	Diff line number	Diff line
		@@ -19,6 +19,9 @@ import zou2024


		def main(dirname, filename_savefig):
		global catalog
		catalog = None

		plot_catalog = True
		plot_catalog_focc = True

		@@ -43,6 +46,8 @@ def main(dirname, filename_savefig):
		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):
		global catalog

		print(ax, zmin, zmax, zcen)

		if plot_miyaji15:
		@@ -125,7 +130,8 @@ def main(dirname, filename_savefig):
		from catalog_combined import CatalogCombined
		from util import get_key_function

		catalog = CatalogCombined(dirname)
		if catalog is None:
		catalog = CatalogCombined(dirname, sql_query="SELECT Z, M, log_LX_2_10 FROM Truth")
		values = (catalog["log_LX_2_10"] > 42) & (zmin <= catalog["Z"]) & (catalog["Z"] < zmax)
		x, dx, y, dy = get_key_function(
		bins,
		@@ -134,6 +140,7 @@ def main(dirname, filename_savefig):
		area_deg2=catalog.get_area(),
		zmin=zmin,
		zmax=zmax,
		nmin=1,
		)
		select = y - dy > 0.0
		x, dx, y, dy = map(lambda a: a[select], [x, dx, y, dy])
		@@ -153,11 +160,13 @@ def main(dirname, filename_savefig):
		# Plot the Xray luminosity function weighted by the black hole
		# occupation fraction
		from catalog_combined import CatalogCombined
		from mbh import get_occupation_fraction
		from util import get_key_function

		catalog = CatalogCombined(dirname)
		if catalog is None:
		catalog = CatalogCombined(dirname, sql_query="SELECT Z, M, log_LX_2_10 FROM Truth")
		values = (catalog["log_LX_2_10"] > 42) & (zmin <= catalog["Z"]) & (catalog["Z"] < zmax)
		values = values.astype(np.float64) * catalog["occupation_fraction"]
		values = values.astype(np.float64) * get_occupation_fraction(catalog["M"])
		x, dx, y, dy = get_key_function(
		bins,
		catalog["log_LX_2_10"],
		@@ -165,6 +174,7 @@ def main(dirname, filename_savefig):
		area_deg2=catalog.get_area(),
		zmin=zmin,
		zmax=zmax,
		nmin=1,
		)
		select = y - dy > 0.0
		x, dx, y, dy = map(lambda a: a[select], [x, dx, y, dy])

tests/test_mbh.py

0 → 100644

+21 −0

Original line number	Diff line number	Diff line
		#!/usr/bin/env python3

		"""Test MBH module."""

		import unittest

		import mbh
		import numpy as np


		class TestMBH(unittest.TestCase):
		def test_get_occupation_fraction(self):
		logMstar = np.array([8.0, 9.0, 10.0])
		test = mbh.get_occupation_fraction(logMstar)
		expected = np.array([0.20, 0.45, 0.80])
		diff = test - expected
		self.assertTrue(np.all(diff < 0.10))


		if __name__ == "__main__":
		unittest.main()