Create an example AGILE truth catalog

This notebooks demonstrates the basic usage of the truth catalog library, and generates an extermely small truth catalog for testing purposes.

In AGILE, each individual class of source (e.g. AGN, galaxy, star) is first contained in its own catalog. The final step in the truth catalog creation is to simply combine these individual source catalogs into one combined catalog that can be used for future image simulations.

The different classes of objects considered are:

Galaxies
AGNs
Stars
Binary stars

The following notebook may be accessed here

Initialize

[ ]:

import os

if not os.path.exists("src/lsst_inaf_agile"):
    os.chdir("../../")
    os.getcwd()
dirname = "data/tests/test_agile"

[ ]:

import logging
import matplotlib.pyplot as plt
import numpy as np
from lsst_inaf_agile.catalog_agn import CatalogAGN
from lsst_inaf_agile.catalog_combined import CatalogCombined
from lsst_inaf_agile.catalog_galaxy import CatalogGalaxy
from lsst_inaf_agile.catalog_star import CatalogStar
from lsst_inaf_agile.egg import Egg
from lsst_inaf_agile.image_simulator import ImageSimulator
from lsst_inaf_agile.merloni2014 import Merloni2014
from lsst_inaf_agile import util

Create the EGG galaxy catalog

[ ]:

# The following list of keyword arguments may be expanded on.
# Look into 'egg-gencat help'for all the available EGG arguments.
egg_kwargs = Egg.get_example_egg_kwargs(dirname + "/egg.fits")
egg_kwargs

[ ]:

egg = Egg(egg_kwargs)
egg.run()

[ ]:

catalog_egg = Egg.read(egg_kwargs["out"])

[ ]:

# Convert to the AGILE galaxy catalog class
catalog_galaxy = CatalogGalaxy(dirname, catalog_egg)

Create the AGN catalog

[ ]:

# Create the AGN catalog
kwargs_catalog_agn = dict(
    dirname=dirname,
    catalog_galaxy=catalog_galaxy,
    type_plambda="zou+2024",
    save_sed=1,
    seed=20251005,
    merloni2014=Merloni2014(1, 0, 0.05, 0.95),
)
catalog_agn = CatalogAGN(**kwargs_catalog_agn)

Examine an example AGN

[ ]:

# Select a type1 AGN
select = (catalog_agn["is_agn"] == 1) & (catalog_agn["is_optical_type2"] == 0)

# Select the first AGN by ID
ids = catalog_agn["ID"][select]
my_id = ids[0]

print(f"Found a type1 AGN with ID = {my_id}")

Examine an AGN SED

[ ]:

lam, flux = catalog_agn.get_sed(my_id)
plt.loglog(lam, flux)

Estimate a single AGN light curve

[ ]:

# Select a type1 AGN
select = (catalog_agn["is_agn"] == 1) & (catalog_agn["is_optical_type2"] == 0)

# Select the first AGN by ID
my_id = catalog_agn["ID"][select][0]

# Estimate its lightcurve.
# The lightcurve is estimated for 10-years with a cadence of one day.
# The returned flux is in the observer frame.
lc = catalog_agn.get_lightcurve(my_id, "lsst-r")
plt.plot(lc)
plt.xlabel("MJD")
plt.ylabel(r"flux $r$ [uJy]")

Estimate Damped Randow Walk parameters \(\tau\) and \(\mathrm{SF}_\infty\)

The following example shows how to estimate the Damped Random Walk parameters \(\tau\) and \(\mathrm{SF}_\infty\) for a single AGN using the LSST \(r\)-band. To estimate the parameters for the whole catalog, one can wrap

[ ]:

# Select a type1 AGN
select = (catalog_agn["is_agn"] == 1) & (catalog_agn["is_optical_type2"] == 0)

# Select the first AGN by ID
my_id = catalog_agn["ID"][select][0]

# LSST bandpass -- one of 'ugrizy'
b = "r"

tau, sf_inf = catalog_agn.get_lightcurve(my_id, f"lsst-{b}", return_tau_sf_inf=True)
print(f"{tau=}, {sf_inf=}")

Create the star and binary star catalogs

[ ]:

# Create the star catalogs
catalog_star = CatalogStar(dirname, catalog_galaxy, is_binary=False)
catalog_binary = CatalogStar(dirname, catalog_galaxy, is_binary=True)

[ ]:

# NOTE: dal Tio+ simulated only 10% of the toal number of binary systems. They suggest a value of
# fbin=0.40 to account for this. To achieve this:
#   1) input star catalog is downsampled to 1 - fbin = 60% of the original size
#   2) the binary catalog is repeated 4 times
# The helper function below does these two steps automatically.
if "_once" not in globals():
    print(f"Old sizes are {len(catalog_star.stars)=}, {len(catalog_binary.stars)=}")
    catalog_star, catalog_binary = CatalogStar.get_star_binary_fbin(
        catalog_star, catalog_binary, fbin=0.40, nrepeat=4
    )
    _once = True
print(f"New sizes are {len(catalog_star.stars)=}, {len(catalog_binary.stars)=}")

Create the combined catalog of AGNs, galaxies, and stars

[ ]:

# NOTE: we set cache=False to prevent an existing catalog file to be loaded from the disk.
catalog_combined = CatalogCombined(
    dirname, catalog_galaxy, catalog_agn, catalog_star, catalog_binary, cache=False
)

Examine the truth catalog columns

Here we list all the available columns in the truth catalog. The meaning of these columns is described in the appendix of Viitanen+2026.

It is important to note that the truth catalog contains ONE row per object. That is, a single row corresponds to either an AGN, galaxy, or a star. Consequently, not all the different classes of objects have all the different columns available. For example, (host) galaxy stellar mass (‘M’) is only available for galaxies (and AGNs).

[ ]:

catalog_combined.get_dtype()

Examine the output files

At this stage, the static truth catalog is completed. To find out which files have been written to disk, we check the output directory.

[ ]:

os.listdir(dirname)

The explanation of the catalog files is as follows:

egg.fits – the EGG galaxy catalog
agn.fits – the AGN catalog
stars.fits – the full star catalog
binaries.fits – the full binary star catalog
catalog.fits – the combined truth catalog

In addition to:

egg-seds* – databse files needed to generate EGG SEDs (refer to EGG documentation)
seds/ – AGN SEDs in EGG format (refer to EGG documentation)
lightcurves/ – stored lightcurves (see the AGN catalog section above)

Note that by default, all galaxy SEDs only reside in the EGG database. These can be stored on disk on demand, but essentially doubles the disk space used which can be problematic on larger catalogs.

Also, light curves are only generated by request. In case images are simulated, then a subset of light curves will be estimated automatically for the sources within the region of interest.

Example plot: (host) galaxy stellar mass versus redshift

[ ]:

plt.figure(dpi=200)
select = catalog_combined["Z"] > 0
plt.plot(catalog_combined["Z"][select], catalog_combined["M"][select], ".")
plt.xlabel(r"$z$")
plt.ylabel(r"$\log (M\,/\,M_\mathrm{star})$");

Example plot: g-band flux distribution

[ ]:

# Example plot: distribution of the observed g-band flux.
plt.figure(dpi=200)
plt.hist(catalog_combined["lsst-g_total"], bins=np.logspace(-3, 3, 61))
plt.loglog()
plt.xlabel(r"lsst-g_total [uJy]")
plt.ylabel("frequency");

A note on the occupation fraction

In the simulation, it is preliminarily assumed that each galaxy has a SMBH, and a corresponding value \(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:

:nbsphinx-math:`begin{equation}: mathrm{has_bh} = U < f_mathrm{occ}(M_mathrm{star}),

end{equation}`

where \(U\) is a uniform random variable \(U \sim \mathrm{Unif}(0, 1)\), and \(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). An example is provided below.

[ ]:

# Select all galaxies
is_galaxy = catalog_combined["Z"] > 0.0

# Select central BHs
has_bh = catalog_combined["has_bh"]

# Combine the two selections
is_galaxy_and_has_bh = is_galaxy & has_bh

# Do some statistical tests
print(f"Total number of galaxies: {is_galaxy.sum()}")
print(f"Total number of galaxies with BHs: {is_galaxy_and_has_bh.sum()}")

# Plot Mstar vs. occupation fraction
plt.title("BH occupation fraction vs. stellar mass")
plt.plot(catalog_combined["M"][is_galaxy], catalog_combined["occupation_fraction"][is_galaxy], ".")
plt.xlabel(r"$\log \left( M_\mathrm{star} / M_\odot \right)$")
plt.ylabel(r"$f_\mathrm{occ}$")

For AGN, the logic remains the same. To weigh AGN by the occupation fraction, one would use the product \(\mathrm{is\_agn} \times \mathrm{has\_bh}\) instead of simply \(\mathrm{is\_agn}\). Logically, this is equivalent to the condition “has BH AND BH is active”. For \(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.

[ ]:

# Select galaxies with BHs
has_bh = catalog_combined["has_bh"]

# Select AGNs
is_agn = catalog_combined["is_agn"]

# Select logLX > 42 objects for the sake of an example
is_loglx_gt_42 = catalog_combined["log_LX_2_10"] > 42

# Combine the selections
has_bh_and_is_agn = has_bh & is_agn
has_bh_and_is_agn = has_bh & is_agn & is_loglx_gt_42


print(f"Total number of AGNs: {is_agn.sum()}")
print(f"Total number of BHs with AGNs: {has_bh_and_is_agn.sum()}")

Simulate example images

This notebooks demonstrates how to use the AGILE truth catalog in order to simulate synthetic LSST images.

We use imSim (https://lsstdesc.org/imSim/index.html) to simulate the images. The purpose of this notebook is not to doucment the feature of imSim, but only to demonstrate how to pass our truth catalog as input to imSim, and how to run software using the tools available in AGILE. Note that imSim only simulates raw LSSTCam data (incl. e.g. sky and instrumental noise), which are not science ready products. The raw images are processed in subsequent steps by using the LSST Science pipelines.

Below, we provide helpful links and resources to further understand the LSST survey, LSSTCam, and the survey strategy.

VRO survey strategy: https://survey-strategy.lsst.io/
opSim: https://usdf-maf.slac.stanford.edu/
imSim: https://lsstdesc.org/imSim/index.html
LSSTCam: https://ctn-001.lsst.io/
LSSTCam detector layout: https://lsstdesc.org/imSim/lsst-camera.html

The following notebook may be accessed here

[ ]:

from lsst_inaf_agile.image_simulator import ImageSimulator

Initialize the ImageSimulator

[ ]:

# Select a baseline
# NOTE: other baselines are available here: https://usdf-maf.slac.stanford.edu/
filename_baseline = "data/baseline/baseline_v4.0_10yrs.db"
image_simulator = ImageSimulator(f"{dirname}/imsim", catalog_combined, filename_baseline)

Explore the visits database

[ ]:

visits = image_simulator.get_visit(limit=10)
visits.info()

[ ]:

visits[:5]

[ ]:

my_observation_id = 779
my_visit = image_simulator.get_visit(observation_id=my_observation_id)
my_visit

Write an instance catalog

Here, we write a single instance catalog. The instance catalog represents a snapshot of the truth catalog, where the magnitudes of the different objects (e.g. AGNs or binary stars) are modified according to the light curve.

[ ]:

# The following command writes a single instance catalog for the observation_id
image_simulator.write_instance_catalog(my_observation_id)

[ ]:

# The following command writes a single instance catalog for the observation_id
image_simulator.write_instance_catalog(my_observation_id)

[ ]:

# The following command writes a single instance catalog OVERRIDING
# the default ra, dec of the pointing -- which is useful for debugging small
# (<10deg2) fields
ra_dec = +150.11916667, +2.20583333
image_simulator.write_instance_catalog(my_observation_id, ra_dec=ra_dec)

Simulate a single visit and single detector

[ ]:

detector = 94
image_simulator.simulate_image(observation_id=my_observation_id, detector=detector)

Explore the output products

[ ]:

dirname_output = os.path.join(dirname, "imsim", str(observation_id), "output")
print("Dirname_output is {dirname_output}")

[ ]:

os.listdir(dirname_output)

On simulating multiple images and/or detectors

Simulating any large dataset is a computationally expensive task. To simulate any significantly large dataset (such as the AGILE DR1), one is adviced to wrap the image_simulator.simulate_image into a callable script which could be executed with multiple cores and/or threads. In AGILE DR1, this orchestration was done with slurm (https://slurm.schedmd.com/documentation.html).