diff --git a/Makefile b/Makefile index 4fc43b7fbd2928e523d2860797a31e4e48ec5d35..d0104ce6640062046c91a2fda5162f757506b05b 100644 --- a/Makefile +++ b/Makefile @@ -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 $< $@ diff --git a/docs/userguide.rst b/docs/userguide.rst index cec333a8b67042f083e529edb32344f2c88ac6fa..faea5ac651f8e008fbc50bb1c41204b540f42b46 100644 --- a/docs/userguide.rst +++ b/docs/userguide.rst @@ -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 -------------------- diff --git a/src/lsst_inaf_agile/catalog_agn.py b/src/lsst_inaf_agile/catalog_agn.py index 78c7e548312e5ba1ab6d25180b7c301ea8421f23..ea164b514557edc89c66c86aaee41b47e69038d8 100644 --- a/src/lsst_inaf_agile/catalog_agn.py +++ b/src/lsst_inaf_agile/catalog_agn.py @@ -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. diff --git a/src/lsst_inaf_agile/catalog_star.py b/src/lsst_inaf_agile/catalog_star.py index bb982dfc97ddd12ebc2e0263646e4ab8eb5f0091..f83dc020e1ec90710fd5a55ba0545c5d6f44356c 100644 --- a/src/lsst_inaf_agile/catalog_star.py +++ b/src/lsst_inaf_agile/catalog_star.py @@ -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,53 +90,54 @@ 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 - ra, dec = (self.catalog_galaxy[k] for k in ("RA", "DEC")) - ra_min = ra.min() - ra_max = ra.max() - dec_min = dec.min() - dec_max = dec.max() + 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() + dec_min = dec.min() + dec_max = dec.max() - query = f""" - SELECT * FROM {table} - WHERE {selection} < {maglim} - AND {ra_min} < ra AND ra < {ra_max} - AND {dec_min} < dec AND dec < {dec_max} - """ - logger.info(query) + query = f""" + SELECT * FROM {table} + WHERE + {ra_min} < ra AND ra < {ra_max} + AND {dec_min} < dec AND dec < {dec_max} + """ + logger.info(query) - import warnings + import warnings - from astropy.units import UnitsWarning + from astropy.units import UnitsWarning - with warnings.catch_warnings(): - warnings.simplefilter("ignore", UnitsWarning) + 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.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): diff --git a/src/lsst_inaf_agile/util.py b/src/lsst_inaf_agile/util.py index c80767d500995c5bef51dcdf4e6d3d969cbbaac0..8f23726ef2d4d4561a3d8cbcea79f7276077e568 100644 --- a/src/lsst_inaf_agile/util.py +++ b/src/lsst_inaf_agile/util.py @@ -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 diff --git a/tests/lsst_inaf_agile/test_catalog_agn.py b/tests/lsst_inaf_agile/test_catalog_agn.py index 82dd33393c7b063e3f34b4811e52ca1f4655bfa6..eabfd25f21d4d513ca058725c32890af89dc251f 100644 --- a/tests/lsst_inaf_agile/test_catalog_agn.py +++ b/tests/lsst_inaf_agile/test_catalog_agn.py @@ -180,6 +180,17 @@ class TestCatalogAGN(TestCase): ebv2 = self.catalog_agn["E_BV"][my_id2] self.assertGreater(ebv2, ebv1) + # test key not in FLUXES + from lsst_inaf_agile.catalog_agn import FLUXES + + id0 = self.catalog_agn["ID"][self.catalog_agn["is_agn"]][0] + key = id0, "lsst-r", False + flux1 = FLUXES[key] + FLUXES.pop(key) + flux2 = self.catalog_agn.get_flux_agn("lsst-r", False, idxs=id0) + self.assertIn(key, FLUXES) + self.assertTrue(flux1 == flux2) + def test_get_flux_agn(self): # range check on flux gflux = self.catalog_agn["lsst-g_point"] @@ -266,3 +277,16 @@ class TestCatalogAGN(TestCase): mean = np.mean(lc_all) flux = self.catalog_agn["lsst-g_point"][id_agn] self.assertTrue(np.isclose(mean, flux)) + + def test_get_occupation_fraction(self): + ret = self.catalog_agn.get_occupation_fraction() + self.assertTrue(np.all((ret >= 0.0) & (ret <= 1.0))) + select1 = self.catalog_agn["M"] < 10.0 + select2 = self.catalog_agn["M"] > 10.0 + self.assertLess(np.mean(ret[select1]), np.mean(ret[select2])) + + def test_get_has_bh(self): + np.random.seed(20260112) + Ngal = self.catalog_agn["M"].size + has_bh = self.catalog_agn.get_has_bh() + self.assertLess(has_bh.sum(), Ngal) diff --git a/tests/lsst_inaf_agile/test_catalog_star.py b/tests/lsst_inaf_agile/test_catalog_star.py index cf20f0a74523fc75d768b34d20fed3f742069c62..bb3d54025e28eb2f1be5d895ea13ae25beb17a35 100644 --- a/tests/lsst_inaf_agile/test_catalog_star.py +++ b/tests/lsst_inaf_agile/test_catalog_star.py @@ -60,16 +60,34 @@ class TestCatalogStar(TestCase): def test_get_stars(self): # test existing filename and single/binary + cs1 = [] for c in self.catalog_star, self.catalog_binary: os.remove(c.filename) - c.get_stars(debug=True) + cs1.append(c.get_stars()) self.assertTrue(os.path.exists(c.filename)) + # call a second time + cs2 = [] + for c in self.catalog_star, self.catalog_binary: + cs2.append(c.get_stars()) + self.assertTrue(os.path.exists(c.filename)) + self.assertTrue(np.all(cs1[0] == cs2[0])) + self.assertTrue(np.all(cs1[1] == cs2[1])) + + # call a third time without a maglim + cs3 = [] + for c in self.catalog_star, self.catalog_binary: + os.remove(c.filename) + cs3.append(c.get_stars(maglim=None)) + self.assertTrue(os.path.exists(c.filename)) + self.assertTrue(np.all(cs1[0] == cs3[0])) + self.assertTrue(np.all(cs1[1] == cs3[1])) + # test selection and maglim maglim = 28 for b in "ugrizy": os.remove(self.catalog_star.filename) - c = self.catalog_star.get_stars(selection=f"{b}mag", maglim=maglim, debug=True) + c = self.catalog_star.get_stars(selection=f"{b}mag", maglim=maglim) self.assertTrue(np.all(c[f"{b}mag"] <= maglim)) def test_is_cepheid(self): @@ -87,7 +105,7 @@ class TestCatalogStar(TestCase): val2 = self.catalog_star.is_c_rich(None) val3 = self.catalog_star.is_c_rich() - self.assertEqual(np.array(val2).shape, self.catalog_star.stars.shape) + self.assertEqual(np.array(val2).shape[0], len(self.catalog_star.stars)) self.assertTrue(np.all(val2 == val3)) def test_is_lpv(self):