Unverified Commit f4743b8d authored by Akke Viitanen's avatar Akke Viitanen
Browse files

add focc usage in the user guide and fix tests 

add focc usage in the user guide

add has_bh column

add test for has_bh column

fix bug in accessing star table length

remove silly debug logic from get_stars

increase test coverage to 100%

fix math formatting

re-write confusing maglim argument

Formerly, the maglim argument was used in the SQL query, and future
calls would return this catalog. This can lead to unexpected behavior if
maglim_new > maglim_old, since maglim_old would be the only one to take
effect. The fix is to always download the full catalog and extract a
magnitude-limited catalog from there.
parent a595f482

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