Merge branch '7-describe-focc-usage-in-the-documentation' into 'feature_documentation' (8c94172c) · Commits · Akke Esa Tapio Viitanen / lsst_inaf_agile

Makefile

+1 −1

Original line number	Diff line number	Diff line
		@@ -299,7 +299,7 @@ coverage:
		pytest --cov=src --cov-report=html tests/

		documentation:
		sphinx-build -M html ./docs ./build -j4
		cd docs && make html

		fig/black_hole_mass_function_redshift.pdf: src/scripts/plots/plot_black_hole_mass_function_redshift.py
		python3 $< $@

docs/userguide.rst

+27 −0

Original line number	Diff line number	Diff line
		@@ -4,6 +4,33 @@ User guide
		Creating truth catalogs of AGNs, galaxies, and stars
		----------------------------------------------------

		A note on the BH occupation fraction
		^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

		In the simulation, it is preliminarily assumed that each galaxy has a SMBH, and
		a corresponding value :math:`M_\mathrm{BH}`. Local observations suggest that
		this might not be the case, and especially at low values of Mstar, a BH might
		be missing altogether. To facilitate for this, in the truth catalog a flag
		``has_bh`` is provided. The formal definition of this flag is:

		.. math::
		\mathrm{has\_bh} = U < f_\mathrm{occ}(M_\mathrm{star}),

		where :math:`U` is a uniform random variable :math:`U \sim \mathrm{Unif}(0,
		1)`, and :math:`f_\mathrm{occ}` follows the observational results of Zou+2025.
		This column may be used as a weight for all calculations of statistical
		distributions. For example, the galaxy BH mass function can be calculated with
		and without the occupation fraction by simply weighting the mass function by
		unity (every galaxy is expected to host a BH), or ``has_bh`` (galaxies are
		expected to host BHs in accordance with the occupation fraction).

		For AGN, the logic remains the same. To weigh AGN by the occupation fraction,
		one would use the product :math:`\mathrm{is\_agn} \times \mathrm{has\_bh}`
		instead of simply :math:`\mathrm{is\_agn}`. Logically, this is equivalent to
		the condition "has BH AND BH is active". For :math:`L_\mathrm{X} >
		10^{42}\,\mathrm{erg}\,\mathrm{s}^{-1}`, this has a negligible effect on the
		AGN population at large, as was investigated in Viitanen+2026.

		Creating lightcurves
		--------------------

src/lsst_inaf_agile/catalog_agn.py

+35 −4

Original line number	Diff line number	Diff line
		@@ -157,6 +157,7 @@ class CatalogAGN:
		("MBH", np.float64, "Black hole mass", "Msun"),
		("log_L_bol", np.float64, "Bolometric AGN luminosity", "erg/s"),
		("occupation_fraction", np.float64, "SMBH occupation fraction", ""),
		("has_bh", bool, "Galaxy contains a SMBH according to the occupation fraction", ""),
		]
		+ [(b.strip() + "_point", np.float64, f"AGN {b} flux", "uJy") for b in bands]
		+ [
		@@ -572,19 +573,42 @@ class CatalogAGN:
		return flux_band

		def get_flux_agn(self, band, rest_frame, idxs=None):
		"""Return vector of AGN pass-band fluxes."""
		"""
		Return vector of AGN pass-band fluxes.

		Parameters
		----------
		band: str
		Band to estimate the flux in e.g. 'lsst-r'.
		rest_frame: bool
		Return rest-frame absolute magnitude instead of observed flux.
		idxs: List[int] or None
		If defined, return a vector containing the fluxes for the given
		indices. Otherwise returns the fluxes for the full catalog
		including zero fluxes for non-AGN galaxies.

		Returns
		-------
		flux_or_magabs: float
		Estimated AGN flux(es) or absolute magnitude(s) based on the given
		flags.
		"""
		flux = np.full(self.catalog.size, np.inf if rest_frame else 0.0)

		select_idx = None
		if idxs is None:
		idxs = self["ID"][self["is_agn"]]
		else:
		select_idx = np.isin(np.arange(flux.size), idxs)

		for idx in idxs:
		# NOTE: fluxes should have been initialized by now
		for idx in np.atleast_1d(idxs):
		if (idx, band, rest_frame) not in FLUXES:
		_ = self._get_sed(idx)
		flux[idx] = FLUXES[idx, band, rest_frame]

		return flux
		if idxs is not None:
		flux = flux[select_idx]
		return np.squeeze(flux)

		def get_log_L_bol(self):
		"""Return log10 of the AGN bolometric luminosity in erg/s."""
		@@ -687,6 +711,13 @@ class CatalogAGN:

		return get_occupation_fraction(self["M"])

		def get_has_bh(self):
		"""
		Return whether the galaxy contains a SMBH according to the occupation fraction.
		"""
		U = np.random.rand(self["M"].size)
		return self["occupation_fraction"] > U

		def get_lightcurve(self, i, band, args, *kwargs):
		"""
		Return an AGN lightcurve.

src/lsst_inaf_agile/catalog_star.py

+32 −32

Original line number	Diff line number	Diff line
		@@ -12,7 +12,6 @@ import numpy as np
		import pyvo as vo
		from astropy import constants
		from astropy import units as u
		from astropy.table import Table

		from . import lightcurve, util

		@@ -91,16 +90,24 @@ class CatalogStar:
		"""Return column corresponding to 'k'."""
		return self.catalog[k]

		def get_stars(self, selection="rmag", maglim=28, debug=False):
		"""Query NOIRlab to retrieve the stellar catalogs from LSST-SIM."""
		if os.path.exists(self.filename):
		return util.read_fits(self.filename)
		def get_stars(self, selection="rmag", maglim=28):
		"""
		Query NOIRlab to retrieve the stellar catalogs from LSST-SIM.

		Parameters
		----------
		selection: str
		Column name for selection in terms of magnitude.
		maglim: float
		Magnitude limit for the 'selection' column.
		"""

		table = "lsst_sim.simdr2"
		if self.is_binary:
		table = "lsst_sim.simdr2_binary"
		selection = "c3_" + selection

		if not os.path.exists(self.filename):
		ra, dec = (self.catalog_galaxy[k] for k in ("RA", "DEC"))
		ra_min = ra.min()
		ra_max = ra.max()
		@@ -109,8 +116,8 @@ class CatalogStar:

		query = f"""
		SELECT * FROM {table}
		WHERE {selection} < {maglim}
		AND {ra_min} < ra AND ra < {ra_max}
		WHERE
		{ra_min} < ra AND ra < {ra_max}
		AND {dec_min} < dec AND dec < {dec_max}
		"""
		logger.info(query)
		@@ -122,22 +129,15 @@ class CatalogStar:
		with warnings.catch_warnings():
		warnings.simplefilter("ignore", UnitsWarning)

		if debug:
		# Return dummy output
		if self.is_binary:
		table = Table.read("data/tests/test_catalog_star/bak/binaries.fits")
		else:
		table = Table.read("data/tests/test_catalog_star/bak/stars.fits")
		table = table[table[selection] < maglim]

		else:
		# Run the syncronous job
		tap_service = vo.dal.TAPService("https://datalab.noirlab.edu/tap")
		tap_results = tap_service.search(query)
		table = tap_results.to_table()

		table.write(self.filename)

		table = util.read_table(self.filename)
		if maglim is not None:
		return table[table[selection] < maglim]
		return table

		def is_cepheid(self, i):

src/lsst_inaf_agile/util.py

+13 −0

Original line number	Diff line number	Diff line
		@@ -410,3 +410,16 @@ def read_fits(filename, args, *kwargs):
		warnings.simplefilter("ignore", UnitsWarning)
		ret = fitsio.read(filename, args, *kwargs)
		return ret


		def read_table(filename):
		"""Read an astropy table but do so silently."""
		import warnings

		from astropy.table import Table
		from astropy.units import UnitsWarning

		with warnings.catch_warnings():
		warnings.simplefilter("ignore", UnitsWarning)
		table = Table.read(filename)
		return table