Loading src/lsst_inaf_agile/catalog_agn.py +202 −57 Original line number Diff line number Diff line Loading @@ -3,9 +3,7 @@ # Email: akke.viitanen@helsinki.fi # Date: 2023-02-13 16:20:51 """ Create Galaxy+AGN mocks for the LSST Italian AGN in-kind contribution """ """Create Galaxy+AGN mocks for the LSST Italian AGN in-kind contribution.""" import logging import multiprocessing Loading Loading @@ -49,21 +47,7 @@ except FileNotFoundError: class CatalogAGN: """ AGN catalog class """ def __init__( self, dirname: str, catalog_galaxy: CatalogGalaxy, type_plambda: str, save_sed: bool, seed: int, merloni2014: Merloni2014, filter_db: str = "data/egg/share/filter-db/db.dat", ): """ Initialize an AGN catalog. AGN catalog class. Parameters ---------- Loading @@ -81,7 +65,38 @@ class CatalogAGN: assumed Merloni2014 obscuration model filter_db: str, optional filename of the EGG filter database (db.dat) Attributes ---------- dirname: str directory name to write the catalog to catalog_galaxy: CatalogGalaxy underlying galaxy catalog type_plambda: str assumed p(lambda) save_sed: bool save AGN SEDs to file seed: int random number seed merloni2014: Merloni2014 assumed Merloni2014 obscuration model filter_db: str, optional filename of the EGG filter database (db.dat) catalog: np.ndarray catalog data """ def __init__( self, dirname: str, catalog_galaxy: CatalogGalaxy, type_plambda: str, save_sed: bool, seed: int, merloni2014: Merloni2014, filter_db: str = "data/egg/share/filter-db/db.dat", ): """Initialize the AGN catalog.""" self.dirname = dirname self.catalog_galaxy = catalog_galaxy self.type_plambda = type_plambda Loading @@ -105,11 +120,11 @@ class CatalogAGN: Parameters ---------- key: str name of column Name of column. Returns ------- value of the column Value of the column. """ if key not in self.catalog.dtype.names: return self.catalog_galaxy[key] Loading @@ -118,7 +133,7 @@ class CatalogAGN: @staticmethod def get_columns(bands, rfbands, is_public=False): """ Return the names, types, descriptions, and units of each column Return the names, types, descriptions, and units of each column. Parameters ---------- Loading @@ -131,7 +146,7 @@ class CatalogAGN: Returns ------- List of 4-tuples of column names, types, descriptions, and units List of 4-tuples of column names, types, descriptions, and units. """ return ( [ Loading Loading @@ -164,12 +179,14 @@ class CatalogAGN: ) def get_dtype(self): """Return the numpy dtype corresponding to the column names and types.""" i = [] for n, t, _, _ in self.get_columns(self.catalog_galaxy.bands, self.catalog_galaxy.rfbands): i += [(n, t)] return np.dtype(i) def get_catalog(self): """Build the catalog column-by-column and write to FITS.""" # Short-circuit for existing catalog filename = f"{self.dirname}/agn.fits" if os.path.exists(filename): Loading @@ -189,6 +206,12 @@ class CatalogAGN: @staticmethod def _get_seed(key): """ Set the random number seed. The random number seed is generated programmatically from the name of the column. """ # NOTE: numpy random seed must be between 0 and 2 ** 32 - 1 # NOTE: hash is not persistent between python runs. Use a custom hash # function instead based on the name of the key Loading @@ -196,6 +219,7 @@ class CatalogAGN: return seed % (2**32 - 1) def get_agn(self, key): """Get AGN attribute corresponding to the column 'key'.""" logger.info(f"Getting AGN attribute {key}") ret = None Loading @@ -217,6 +241,7 @@ class CatalogAGN: return ret def get_log_lambda_SAR(self, add_ctk=True, idxs=None): """Return log10 of the specific black hole accretion rate in erg/s/Msun.""" if self.type_plambda != "zou+2024": raise ValueError Loading Loading @@ -257,11 +282,13 @@ class CatalogAGN: return log_lambda_SAR def get_log_LX_2_10(self): """Return log10 of the intrinsic 2-10 keV X-ray luminosity in erg/s.""" # NOTE: convention of assigning negative lambda to CTK AGN x = np.where(np.isfinite(self["log_lambda_SAR"]), np.abs(self["log_lambda_SAR"]), -np.inf) return x + self["M"] def get_log_FX_2_10(self, Gamma=1.9): """Return log10 of the intrinsic 2-10 keV X-ray flux in erg/s/cm2 in the observer frame.""" z = self["Z"].astype(np.float64) d = self["D"].astype(np.float64) * u.Mpc.to(u.cm) Loading @@ -269,12 +296,24 @@ class CatalogAGN: return self["log_LX_2_10"] - np.log10(4 * np.pi * d**2) + (2 - Gamma) * np.log10(1 + z) def get_log_FX_2_7(self, Gamma=1.9): """ Return log10 of the intrinsic 2-7 keV X-ray flux in erg/s/cm2 in the observer frame. Band conversion is done from the 2-10 keV assuming power-law index 'Gamma'. Parameters ---------- Gamma: float Assumed power-law index for the band conversion. """ return np.log10(util.convert_flux(10 ** self["log_FX_2_10"], 2.0, 10.0, 2.0, 7.0, Gamma=Gamma)) def get_log_L_2_keV(self, Gamma=1.9, wavelength=6.2): """ Returns monochromatic X-ray luminosity at lambda = wavelength in erg/s Hz^-1 to be used for the alpha_ox Lx = restframe 2-10 kev luminosity. Return monochromatic X-ray luminosity at lambda = wavelength in erg/s/Hz. To be used for the alpha_ox Lx = restframe 2-10 kev luminosity. maybe we can work directly with frequencies and have it in one line """ Lx = 10 ** self["log_LX_2_10"] Loading @@ -282,20 +321,21 @@ class CatalogAGN: return np.log10((K * wavelength ** (Gamma - 1)) / 2.998e18) def get_log_L_2500(self, alpha=0.952, beta=2.138, scatter=0.40): # def get_log_L_2500(self, alpha=0.952, beta=2.138, scatter=0.20): """ Returns the 2500 ang° monochromatic luminosity (in erg/s). It uses Lusso+10 eq. 5 (inverted) Lx = alpha L_opt - beta """ Return the 2500 angstrom° monochromatic luminosity (in erg/s). Uses Lusso+10 Eq. 5 (inverted) assuming Lx = alpha L_opt - beta """ return lusso2010.get_log_l_2500(self["log_L_2_keV"], alpha, beta, scatter) def get_is_optical_type2(self, func_get_f_obs=None, use_f_obs_for_ctk_agn=False): # Get obscured AGN fraction from Merloni+2014 as a function of z, LX # select_ctn = self["is_agn_ctn"] select_ctk = self["is_agn_ctk"] """Return boolean vector selecting optical type2 AGN. The optical type2 AGN fraction is evaluated based on Merloni+2014 assuming obscured fraction as a function of intrinsic X-ray luminosity and redshift. """ select_ctk = self["is_agn_ctk"] if func_get_f_obs is None: func_get_f_obs = self.merloni2014.get_f_obs Loading @@ -319,6 +359,11 @@ class CatalogAGN: mu_type_2=0.3, ebv_ctk=9.0, ): """Return AGN reddening E(B-V) in ABmag. The E(B-V) is assumed to be different for type1 and type2 AGN, while CTK AGN are assumed to be all type2. """ type_1_ebv = (np.linspace(0, 1, 101),) type_2_ebv = (np.linspace(0, 3, 301),) Loading @@ -328,7 +373,6 @@ class CatalogAGN: return np.interp(np.random.rand(N_AGN), cumulative, ebv_range) def hopkins04(x, alpha, n): """p(E_BV)""" y = 1 / (1 + (x * alpha) ** n) return y / np.trapezoid(y, x) Loading @@ -355,8 +399,9 @@ class CatalogAGN: @staticmethod def get_plambda_ctk(loglambda, mstar, z, t): r""" Return the accretion rate distribution of the CTK AGN population by combining the accretion rate distribution of Zou+2024 (CTN AGN) and the Return the accretion rate distribution of the CTK AGN population. Combine the accretion rate distribution of Zou+2024 (CTN AGN) and the CTK AGN fraction of Ueda+2014. The p_CTK is defined as: p_CTK = f_CTK_AGN / (f_CTK_AGN - 1) * p_CTN Loading @@ -367,7 +412,6 @@ class CatalogAGN: and p_CTN is the accretion rate distribution of CTN AGN. """ p_ctn = 10 ** zou2024.get_log_plambda(loglambda, mstar, z, t) lx = loglambda + mstar Loading @@ -385,27 +429,28 @@ class CatalogAGN: @staticmethod def get_plambda_ctk_prime(loglambda, mstar, z, t, dloglambda): """ Return the CTK AGN accretion rate distribution for a REDUCED galaxy population characterized by Ngal' = Ngal - Nctn. Return the reduced CTK AGN accretion rate distribution. The reduced galaxy population refers to a population of objects reduced already by the number of CTN AGN i.e. Ngal' = Ngal - Nctn. So that the total AGN number remains to be given by CTN + CTK. The resulting p(lambda) is given by p_ctk_prime = p_ctk / (1 - p_ctn * dloglambda) """ p_ctk = CatalogAGN.get_plambda_ctk(loglambda, mstar, z, t) p_ctn = 10 ** zou2024.get_log_plambda(loglambda, mstar, z, t) p_ctk_prime = p_ctk / (1 - p_ctn * dloglambda) return p_ctk_prime def get_is_agn_ctn(self): """Render a random sample of AGN "active" according to the duty cycle""" """Render a random sample of AGN "active" according to the duty cycle.""" # NOTE: avoid non-AGN and CTK AGN when assigning this flag return np.isfinite(self["log_lambda_SAR"]) * self["log_lambda_SAR"] >= 32.00 def get_is_agn_ctk(self): """Assign CTK AGN fraction randomly""" """Assign CTK AGN fraction randomly.""" # NOTE: convention of assigning negative lambda to CTK AGN x = self["log_lambda_SAR"] is_agn_ctk = np.isfinite(x) * (x <= -32) Loading @@ -413,9 +458,24 @@ class CatalogAGN: return is_agn_ctk def get_is_agn(self): """Return boolean vector selecting CTN or CTK AGN.""" return self["is_agn_ctn"] + self["is_agn_ctk"] def get_sed(self, i): """Return an AGN SED. Parameters ---------- i: int ID of AGN for which SED is returned. Returns ------- lambda: array_like Wavelength in microns. flux: array_like observer frame flux in microjanskies. """ filename = f"{self.dirname}/seds/agn-seds-{i}.fits" if not os.path.exists(filename): return None Loading @@ -424,6 +484,17 @@ class CatalogAGN: return fits["LAMBDA"][0], fits["FLUX"][0] def _get_sed(self, i, ratio_max=0.90, dlog_wav=7.65e-4): """Generate an AGN SED. Parameters ---------- i: int ID of AGN for which SED is returned. ratio_max: float Maximum allowed ratio flux_agn / flux_total for type2 AGN. dlog_wav: float SED wavelength resolution in dex. """ logger.debug(f"Getting AGN SED {i} {self.catalog.size}") # Get the wavlen in angstrom Loading Loading @@ -497,6 +568,25 @@ class CatalogAGN: return agn_sed def _get_flux_single(self, my_sed, band, rest_frame, redshift, distance): """Return a single bandbass flux by integrating the SED. Take a rest-frame intrinsic SED and integrate it along a filter passband. The SED is shifted to either rest-frame (10pc) or observer frame in accordance with the parameters. Parameters ---------- my_sed: QSOGEN SED object Input QSOGEN SED. band: str EGG-style passband name. rest_frame: bool Is rest_frame (else observer frame). redshift: float Redshift of the SED distance: float Distance at which to evaluate the flux. """ # Get the flux in observed frame or rest frame if rest_frame: redshift = 0 Loading @@ -519,6 +609,7 @@ class CatalogAGN: return flux_band def get_flux_agn(self, band, rest_frame, idxs=None): """Return vector of AGN pass-band fluxes.""" flux = np.full(self.catalog.size, np.inf if rest_frame else 0.0) if idxs is None: Loading @@ -533,6 +624,7 @@ class CatalogAGN: return flux def get_log_L_bol(self): """Return log10 of the AGN bolometric luminosity in erg/s.""" # NOTE: this function initializes the SEDs... log_L_bol = np.full(self.catalog.size, -np.inf) for i in self["ID"][self["is_agn"]]: Loading @@ -544,6 +636,16 @@ class CatalogAGN: def get_MBH( self, sigma_total=0.50, offset=0.00, sigma_intrinsic=False, sigma_observational=None, select=None ): """Return log10 of the supermassive black hole mass in Msun. The MBH is assigned in accordance with the continuity equation after which log-normal scatter is added to the final value of the MBH. Parameters ---------- sigma_total: float Magnitude of log-normal scatter in dex. """ # Get MBH log_mbh = get_log_mbh_continuity_new2(self["M"], self["Z"]) Loading @@ -554,10 +656,7 @@ class CatalogAGN: return log_mbh + sigma_tot + offset def get_log_L_AGN_5_GHz(self, scatter=True): """ AGN Radio luminosity according to the fundamental plane Gultekin+19 """ """Return AGN radio luminosity according to the fundamental plane Gultekin+19.""" # Eq. 5 # mu = mu0 + xi_mu_R * R + xi_mu_X * X # mu = log10(M / 1e8) Loading @@ -581,15 +680,15 @@ class CatalogAGN: return R + 38 def get_log_lambda_Edd(self, log_l_edd=LOG_LAMBDA_EDD): """ Get the accretion ratio with respect to Eddington """ """Get the accretion ratio with respect to Eddington.""" return self["log_L_bol"] - log_l_edd - self["MBH"] def get_magnitude_agn_dereddened(self, key): """ Return a dereddened optical/UV magnitude. If a different key is given, simply return the value corresponding to the key. Return a dereddened optical/UV magnitude. If a different key is given, simply return the value corresponding to the key. """ # NOTE: discussion with Ivano on 20250227 # updated with 12.6236 for UV from Ivano 20250303 Loading @@ -601,11 +700,35 @@ class CatalogAGN: return self[key] - a * self["E_BV"] def get_occupation_fraction(self): """ Return black hole occupation fraction. The black hole occupation fraction corresponds to the probability that a host galaxy with mass Mstar hosts a (super)massive black hole. This could especially prevalent in the dwarf galaxy (<1e9 Msun) regime. The black hole occupation fraction applied here is based on Zou+2025 (accepted). """ from mbh import get_occupation_fraction return get_occupation_fraction(self["M"]) def get_lightcurve(self, i, band, mjd=None, *args, **kwargs): """ Return an AGN lightcurve. The lightcurve for an AGN is returned for the given band and MJD. Parameters ---------- i: int ID of AGN for which SED is returned. band: str EGG-style passband name. mjd: float or array_like or None MJD of observation. """ if mjd is None: mjd = util.get_mjd_vec() Loading Loading @@ -638,5 +761,27 @@ class CatalogAGN: return lc def get_lightcurve_mjd(self, i, band, mjd, mjd0): """ Return an AGN lightcurve at given MJD. The lightcurve is first estimated assuming the timespan starting from 0. The lightcurve is then evaluated at MJD - MJD0 to return the final value of the lightcurve. Parameters ---------- i: int ID of AGN for which SED is returned. band: str EGG-style passband name. mjd: float or array_like MJD of observation. mjd0: float MJD offset corresponding to the first light of the survey. Returns ------- The value of the lightcurve at given MJD in microjanskies. """ lc = self.get_lightcurve(i, band) return np.interp(mjd - mjd0, util.get_mjd_vec(), lc) src/lsst_inaf_agile/catalog_combined.py +114 −22 Original line number Diff line number Diff line Loading @@ -3,9 +3,7 @@ # Email: akke.viitanen@helsinki.fi # Date: 2023-02-13 16:20:51 """ Create Galaxy+AGN mocks for the LSST Italian AGN in-kind contribution """ """Module for the combined catalog of AGN, galaxies, and stars.""" import logging import os Loading @@ -19,27 +17,65 @@ logger = logging.getLogger(__name__) class CatalogCombined: """Combined catalog class of AGN, galaxies, and stars. Parameters ---------- dirname: str Directory name to contain the catalog. catalog_galaxy: catalog_galaxy.CatalogGalaxy or None Galaxy catalog. catalog_agn: catalog_agn.CatalogAGN or None AGN catalog. catalog_star: catalog_star.CatalogStar or None Single star catalog. catalog_binary: catalog_star.CatalogStar or None Binary star catalog. Attributes ---------- dirname: str Directory name to contain the catalog. catalog_galaxy: catalog_galaxy.CatalogGalaxy or None Galaxy catalog. catalog_agn: catalog_agn.CatalogAGN or None AGN catalog. catalog_star: catalog_star.CatalogStar or None Single star catalog. catalog_binary: catalog_star.CatalogStar or None Binary star catalog. """ def __init__( self, dirname: str, catalog_galaxy=None, catalog_agn=None, catalog_star=None, catalog_binary=None ): self.dirname = dirname self.filename = f"{dirname}/catalog.fits" self.catalog_galaxy = catalog_galaxy self.catalog_agn = catalog_agn self.catalog_star = catalog_star self.catalog_binary = catalog_binary # NOTE: short-circuit for an existing truth catalog if os.path.exists(self.filename): self.catalog_combined = fitsio.read(self.filename) if os.path.exists(f := self.get_filename()): self.catalog_combined = fitsio.read(f) return self.catalog_combined = self.get_catalog_combined() self.catalog_combined = self.postprocess() self.write() def get_filename(self): """Return the FITS filename of the combined catalog""" return f"{self.dirname}/catalog.fits" def get_dtype(self): """ Return the combined catalog dtype constructed from the individual input catalogs. Returns ------- List of (name, type) corresponding to the names and types of the columns. """ dtype = {} for catalog in (self.catalog_galaxy, self.catalog_agn, self.catalog_star, self.catalog_binary): for name in catalog.catalog.dtype.names: Loading @@ -55,7 +91,8 @@ class CatalogCombined: def get_number_galaxy_star_binary(self): """ Return the number of galaxies, stars, and binaries Return the tuple (N_galaxy, N_star, N_binary) where N is the number of objects in the combined catalog. """ return ( len(self.catalog_galaxy.catalog), Loading @@ -64,7 +101,7 @@ class CatalogCombined: ) def get_catalog_combined(self): """ r""" Constructs the combined catalog based on the dtype. For each column in the catalog, the column is looked for in the Loading @@ -72,14 +109,13 @@ class CatalogCombined: are assigned the values from the reference catalog. Note that some columns may appear in multiple catalog. The indices correspond to the following type idx_first idx_last The indices correspond to the following object types: AGN \in [0, n_galaxy[, galaxy \in [0, n_galaxy[, star \in [n_galaxy, n_galaxy + n_star], and binary \in [n_galaxy + n_star, n_total[. AGN 0 n_galaxy galaxy 0 n_galaxy star n_galaxy n_galaxy + n_star binary n_galaxy + n_star n_total Returns ------- The combined catalog as a numpy ndarray. """ n_galaxy, n_star, n_binary = self.get_number_galaxy_star_binary() Loading @@ -104,12 +140,14 @@ class CatalogCombined: # Handles star if name in self.catalog_star.catalog.dtype.names: catalog_combined[name][n_galaxy : n_galaxy + n_star] = self.catalog_star[name] # NOTE: set the flags to False for stars elif catalog_combined.dtype[name] == np.dtype("bool"): catalog_combined[name][n_galaxy : n_galaxy + n_star] = False # Handles binary star if name in self.catalog_binary.catalog.dtype.names: catalog_combined[name][n_galaxy + n_star :] = self.catalog_binary[name] # NOTE: set the flags to False for stars elif catalog_combined.dtype[name] == np.dtype("bool"): catalog_combined[name][n_galaxy + n_star :] = False Loading @@ -131,7 +169,18 @@ class CatalogCombined: return self.catalog_combined def get_flux_total(self, band, rest_frame=False): # NOTE: override rest-frame behavior as EGG does not have absmags for disk/bulge separately """ Return the total flux. The total flux is defined for each type of object as follows: - AGN: sum of galaxy bulge, galaxy disk, and point fluxes - galaxy: sum of galaxy bulge and disk fluxes - star: point fluxes - binary star: sum of the component point fluxes """ # NOTE: override rest-frame behavior as EGG does not have absmags for # disk/bulge separately keys = ["disk", "bulge", "point"] if rest_frame: keys = ["", "point"] Loading @@ -151,19 +200,56 @@ class CatalogCombined: return ret def write(self): fitsio.write(self.filename, self.catalog_combined, clobber=True) """ Writes the combined catalog to a FITS file """ fitsio.write(self.get_filename(), self.catalog_combined, clobber=True) return self.catalog_combined def ingest(self, filename_database, name_table, if_exists="replace"): """ Ingest the truth catalog into a sqlite3 database. Parameters ---------- filename_database: str filename of the sqlite3 database name_table: str name of the table in the sqlite3 database """ import sqlite3 from astropy.table import Table table = Table(self.catalog) df = table.to_pandas() df.reset_index() with sqlite3.connect(filename_database) as con: df.to_sql(name_table, con, index=False, if_exists=if_exists, chunksize=1024, method="multi") def get_is_galaxy(self): """ Return boolean vector selecting galaxies """ return np.isfinite(self.catalog_combined["Z"]) def get_is_star(self): """ Return boolean vector selecting stars """ return ~self.get_is_galaxy() def get_index_star(self, i): """ Return the index in the stellar catalog corresponding to 'i' in the combined catalog combined catalog. Parameters ---------- i: int ID of the truth catalog object """ n_galaxy, n_star, n_binary = self.get_number_galaxy_star_binary() Loading @@ -181,7 +267,7 @@ class CatalogCombined: def write_reference_catalog(self, filename: str, maglim: float, selection_band: str): """ Write a mock LSST reference catalog in csv format to the given Write a mock LSST reference catalog in CSV format to the given 'filename'. The reference catalog is cut to brighter than the magnitude limit 'maglim' in the band 'selection_band'. """ Loading @@ -198,7 +284,7 @@ class CatalogCombined: # Solve the parallax parallax = np.zeros_like(c, dtype=float) select_star = ~np.isfinite(c["Z"]) select_star = self.get_is_star() parallax[select_star] = util.distance_modulus_to_parallax(c["D"][select_star]) # Initialize the data to be written Loading Loading @@ -241,7 +327,7 @@ class CatalogCombined: "pmra_pmdec_corr": np.zeros(c.size), } # Write the reference catalog.csv # Write reference catalog.csv with open(filename, "w") as f: to_write = ",".join(data.keys()) print(f"{to_write}", file=f) Loading @@ -258,9 +344,15 @@ class CatalogCombined: logger.info(f"Wrote reference catalog {filename}") def get_area(self): """ Return the EGG catalog area in deg2 """ cmd = fitsio.read(f"{self.dirname}/egg.fits", columns="CMD")[0] area_string = re.findall("area=([0-9.]+)", cmd)[0] return float(area_string) def __getitem__(self, key): """ Return combined catalog value corresponding to 'key' """ return self.catalog_combined[key] Loading
src/lsst_inaf_agile/catalog_agn.py +202 −57 Original line number Diff line number Diff line Loading @@ -3,9 +3,7 @@ # Email: akke.viitanen@helsinki.fi # Date: 2023-02-13 16:20:51 """ Create Galaxy+AGN mocks for the LSST Italian AGN in-kind contribution """ """Create Galaxy+AGN mocks for the LSST Italian AGN in-kind contribution.""" import logging import multiprocessing Loading Loading @@ -49,21 +47,7 @@ except FileNotFoundError: class CatalogAGN: """ AGN catalog class """ def __init__( self, dirname: str, catalog_galaxy: CatalogGalaxy, type_plambda: str, save_sed: bool, seed: int, merloni2014: Merloni2014, filter_db: str = "data/egg/share/filter-db/db.dat", ): """ Initialize an AGN catalog. AGN catalog class. Parameters ---------- Loading @@ -81,7 +65,38 @@ class CatalogAGN: assumed Merloni2014 obscuration model filter_db: str, optional filename of the EGG filter database (db.dat) Attributes ---------- dirname: str directory name to write the catalog to catalog_galaxy: CatalogGalaxy underlying galaxy catalog type_plambda: str assumed p(lambda) save_sed: bool save AGN SEDs to file seed: int random number seed merloni2014: Merloni2014 assumed Merloni2014 obscuration model filter_db: str, optional filename of the EGG filter database (db.dat) catalog: np.ndarray catalog data """ def __init__( self, dirname: str, catalog_galaxy: CatalogGalaxy, type_plambda: str, save_sed: bool, seed: int, merloni2014: Merloni2014, filter_db: str = "data/egg/share/filter-db/db.dat", ): """Initialize the AGN catalog.""" self.dirname = dirname self.catalog_galaxy = catalog_galaxy self.type_plambda = type_plambda Loading @@ -105,11 +120,11 @@ class CatalogAGN: Parameters ---------- key: str name of column Name of column. Returns ------- value of the column Value of the column. """ if key not in self.catalog.dtype.names: return self.catalog_galaxy[key] Loading @@ -118,7 +133,7 @@ class CatalogAGN: @staticmethod def get_columns(bands, rfbands, is_public=False): """ Return the names, types, descriptions, and units of each column Return the names, types, descriptions, and units of each column. Parameters ---------- Loading @@ -131,7 +146,7 @@ class CatalogAGN: Returns ------- List of 4-tuples of column names, types, descriptions, and units List of 4-tuples of column names, types, descriptions, and units. """ return ( [ Loading Loading @@ -164,12 +179,14 @@ class CatalogAGN: ) def get_dtype(self): """Return the numpy dtype corresponding to the column names and types.""" i = [] for n, t, _, _ in self.get_columns(self.catalog_galaxy.bands, self.catalog_galaxy.rfbands): i += [(n, t)] return np.dtype(i) def get_catalog(self): """Build the catalog column-by-column and write to FITS.""" # Short-circuit for existing catalog filename = f"{self.dirname}/agn.fits" if os.path.exists(filename): Loading @@ -189,6 +206,12 @@ class CatalogAGN: @staticmethod def _get_seed(key): """ Set the random number seed. The random number seed is generated programmatically from the name of the column. """ # NOTE: numpy random seed must be between 0 and 2 ** 32 - 1 # NOTE: hash is not persistent between python runs. Use a custom hash # function instead based on the name of the key Loading @@ -196,6 +219,7 @@ class CatalogAGN: return seed % (2**32 - 1) def get_agn(self, key): """Get AGN attribute corresponding to the column 'key'.""" logger.info(f"Getting AGN attribute {key}") ret = None Loading @@ -217,6 +241,7 @@ class CatalogAGN: return ret def get_log_lambda_SAR(self, add_ctk=True, idxs=None): """Return log10 of the specific black hole accretion rate in erg/s/Msun.""" if self.type_plambda != "zou+2024": raise ValueError Loading Loading @@ -257,11 +282,13 @@ class CatalogAGN: return log_lambda_SAR def get_log_LX_2_10(self): """Return log10 of the intrinsic 2-10 keV X-ray luminosity in erg/s.""" # NOTE: convention of assigning negative lambda to CTK AGN x = np.where(np.isfinite(self["log_lambda_SAR"]), np.abs(self["log_lambda_SAR"]), -np.inf) return x + self["M"] def get_log_FX_2_10(self, Gamma=1.9): """Return log10 of the intrinsic 2-10 keV X-ray flux in erg/s/cm2 in the observer frame.""" z = self["Z"].astype(np.float64) d = self["D"].astype(np.float64) * u.Mpc.to(u.cm) Loading @@ -269,12 +296,24 @@ class CatalogAGN: return self["log_LX_2_10"] - np.log10(4 * np.pi * d**2) + (2 - Gamma) * np.log10(1 + z) def get_log_FX_2_7(self, Gamma=1.9): """ Return log10 of the intrinsic 2-7 keV X-ray flux in erg/s/cm2 in the observer frame. Band conversion is done from the 2-10 keV assuming power-law index 'Gamma'. Parameters ---------- Gamma: float Assumed power-law index for the band conversion. """ return np.log10(util.convert_flux(10 ** self["log_FX_2_10"], 2.0, 10.0, 2.0, 7.0, Gamma=Gamma)) def get_log_L_2_keV(self, Gamma=1.9, wavelength=6.2): """ Returns monochromatic X-ray luminosity at lambda = wavelength in erg/s Hz^-1 to be used for the alpha_ox Lx = restframe 2-10 kev luminosity. Return monochromatic X-ray luminosity at lambda = wavelength in erg/s/Hz. To be used for the alpha_ox Lx = restframe 2-10 kev luminosity. maybe we can work directly with frequencies and have it in one line """ Lx = 10 ** self["log_LX_2_10"] Loading @@ -282,20 +321,21 @@ class CatalogAGN: return np.log10((K * wavelength ** (Gamma - 1)) / 2.998e18) def get_log_L_2500(self, alpha=0.952, beta=2.138, scatter=0.40): # def get_log_L_2500(self, alpha=0.952, beta=2.138, scatter=0.20): """ Returns the 2500 ang° monochromatic luminosity (in erg/s). It uses Lusso+10 eq. 5 (inverted) Lx = alpha L_opt - beta """ Return the 2500 angstrom° monochromatic luminosity (in erg/s). Uses Lusso+10 Eq. 5 (inverted) assuming Lx = alpha L_opt - beta """ return lusso2010.get_log_l_2500(self["log_L_2_keV"], alpha, beta, scatter) def get_is_optical_type2(self, func_get_f_obs=None, use_f_obs_for_ctk_agn=False): # Get obscured AGN fraction from Merloni+2014 as a function of z, LX # select_ctn = self["is_agn_ctn"] select_ctk = self["is_agn_ctk"] """Return boolean vector selecting optical type2 AGN. The optical type2 AGN fraction is evaluated based on Merloni+2014 assuming obscured fraction as a function of intrinsic X-ray luminosity and redshift. """ select_ctk = self["is_agn_ctk"] if func_get_f_obs is None: func_get_f_obs = self.merloni2014.get_f_obs Loading @@ -319,6 +359,11 @@ class CatalogAGN: mu_type_2=0.3, ebv_ctk=9.0, ): """Return AGN reddening E(B-V) in ABmag. The E(B-V) is assumed to be different for type1 and type2 AGN, while CTK AGN are assumed to be all type2. """ type_1_ebv = (np.linspace(0, 1, 101),) type_2_ebv = (np.linspace(0, 3, 301),) Loading @@ -328,7 +373,6 @@ class CatalogAGN: return np.interp(np.random.rand(N_AGN), cumulative, ebv_range) def hopkins04(x, alpha, n): """p(E_BV)""" y = 1 / (1 + (x * alpha) ** n) return y / np.trapezoid(y, x) Loading @@ -355,8 +399,9 @@ class CatalogAGN: @staticmethod def get_plambda_ctk(loglambda, mstar, z, t): r""" Return the accretion rate distribution of the CTK AGN population by combining the accretion rate distribution of Zou+2024 (CTN AGN) and the Return the accretion rate distribution of the CTK AGN population. Combine the accretion rate distribution of Zou+2024 (CTN AGN) and the CTK AGN fraction of Ueda+2014. The p_CTK is defined as: p_CTK = f_CTK_AGN / (f_CTK_AGN - 1) * p_CTN Loading @@ -367,7 +412,6 @@ class CatalogAGN: and p_CTN is the accretion rate distribution of CTN AGN. """ p_ctn = 10 ** zou2024.get_log_plambda(loglambda, mstar, z, t) lx = loglambda + mstar Loading @@ -385,27 +429,28 @@ class CatalogAGN: @staticmethod def get_plambda_ctk_prime(loglambda, mstar, z, t, dloglambda): """ Return the CTK AGN accretion rate distribution for a REDUCED galaxy population characterized by Ngal' = Ngal - Nctn. Return the reduced CTK AGN accretion rate distribution. The reduced galaxy population refers to a population of objects reduced already by the number of CTN AGN i.e. Ngal' = Ngal - Nctn. So that the total AGN number remains to be given by CTN + CTK. The resulting p(lambda) is given by p_ctk_prime = p_ctk / (1 - p_ctn * dloglambda) """ p_ctk = CatalogAGN.get_plambda_ctk(loglambda, mstar, z, t) p_ctn = 10 ** zou2024.get_log_plambda(loglambda, mstar, z, t) p_ctk_prime = p_ctk / (1 - p_ctn * dloglambda) return p_ctk_prime def get_is_agn_ctn(self): """Render a random sample of AGN "active" according to the duty cycle""" """Render a random sample of AGN "active" according to the duty cycle.""" # NOTE: avoid non-AGN and CTK AGN when assigning this flag return np.isfinite(self["log_lambda_SAR"]) * self["log_lambda_SAR"] >= 32.00 def get_is_agn_ctk(self): """Assign CTK AGN fraction randomly""" """Assign CTK AGN fraction randomly.""" # NOTE: convention of assigning negative lambda to CTK AGN x = self["log_lambda_SAR"] is_agn_ctk = np.isfinite(x) * (x <= -32) Loading @@ -413,9 +458,24 @@ class CatalogAGN: return is_agn_ctk def get_is_agn(self): """Return boolean vector selecting CTN or CTK AGN.""" return self["is_agn_ctn"] + self["is_agn_ctk"] def get_sed(self, i): """Return an AGN SED. Parameters ---------- i: int ID of AGN for which SED is returned. Returns ------- lambda: array_like Wavelength in microns. flux: array_like observer frame flux in microjanskies. """ filename = f"{self.dirname}/seds/agn-seds-{i}.fits" if not os.path.exists(filename): return None Loading @@ -424,6 +484,17 @@ class CatalogAGN: return fits["LAMBDA"][0], fits["FLUX"][0] def _get_sed(self, i, ratio_max=0.90, dlog_wav=7.65e-4): """Generate an AGN SED. Parameters ---------- i: int ID of AGN for which SED is returned. ratio_max: float Maximum allowed ratio flux_agn / flux_total for type2 AGN. dlog_wav: float SED wavelength resolution in dex. """ logger.debug(f"Getting AGN SED {i} {self.catalog.size}") # Get the wavlen in angstrom Loading Loading @@ -497,6 +568,25 @@ class CatalogAGN: return agn_sed def _get_flux_single(self, my_sed, band, rest_frame, redshift, distance): """Return a single bandbass flux by integrating the SED. Take a rest-frame intrinsic SED and integrate it along a filter passband. The SED is shifted to either rest-frame (10pc) or observer frame in accordance with the parameters. Parameters ---------- my_sed: QSOGEN SED object Input QSOGEN SED. band: str EGG-style passband name. rest_frame: bool Is rest_frame (else observer frame). redshift: float Redshift of the SED distance: float Distance at which to evaluate the flux. """ # Get the flux in observed frame or rest frame if rest_frame: redshift = 0 Loading @@ -519,6 +609,7 @@ class CatalogAGN: return flux_band def get_flux_agn(self, band, rest_frame, idxs=None): """Return vector of AGN pass-band fluxes.""" flux = np.full(self.catalog.size, np.inf if rest_frame else 0.0) if idxs is None: Loading @@ -533,6 +624,7 @@ class CatalogAGN: return flux def get_log_L_bol(self): """Return log10 of the AGN bolometric luminosity in erg/s.""" # NOTE: this function initializes the SEDs... log_L_bol = np.full(self.catalog.size, -np.inf) for i in self["ID"][self["is_agn"]]: Loading @@ -544,6 +636,16 @@ class CatalogAGN: def get_MBH( self, sigma_total=0.50, offset=0.00, sigma_intrinsic=False, sigma_observational=None, select=None ): """Return log10 of the supermassive black hole mass in Msun. The MBH is assigned in accordance with the continuity equation after which log-normal scatter is added to the final value of the MBH. Parameters ---------- sigma_total: float Magnitude of log-normal scatter in dex. """ # Get MBH log_mbh = get_log_mbh_continuity_new2(self["M"], self["Z"]) Loading @@ -554,10 +656,7 @@ class CatalogAGN: return log_mbh + sigma_tot + offset def get_log_L_AGN_5_GHz(self, scatter=True): """ AGN Radio luminosity according to the fundamental plane Gultekin+19 """ """Return AGN radio luminosity according to the fundamental plane Gultekin+19.""" # Eq. 5 # mu = mu0 + xi_mu_R * R + xi_mu_X * X # mu = log10(M / 1e8) Loading @@ -581,15 +680,15 @@ class CatalogAGN: return R + 38 def get_log_lambda_Edd(self, log_l_edd=LOG_LAMBDA_EDD): """ Get the accretion ratio with respect to Eddington """ """Get the accretion ratio with respect to Eddington.""" return self["log_L_bol"] - log_l_edd - self["MBH"] def get_magnitude_agn_dereddened(self, key): """ Return a dereddened optical/UV magnitude. If a different key is given, simply return the value corresponding to the key. Return a dereddened optical/UV magnitude. If a different key is given, simply return the value corresponding to the key. """ # NOTE: discussion with Ivano on 20250227 # updated with 12.6236 for UV from Ivano 20250303 Loading @@ -601,11 +700,35 @@ class CatalogAGN: return self[key] - a * self["E_BV"] def get_occupation_fraction(self): """ Return black hole occupation fraction. The black hole occupation fraction corresponds to the probability that a host galaxy with mass Mstar hosts a (super)massive black hole. This could especially prevalent in the dwarf galaxy (<1e9 Msun) regime. The black hole occupation fraction applied here is based on Zou+2025 (accepted). """ from mbh import get_occupation_fraction return get_occupation_fraction(self["M"]) def get_lightcurve(self, i, band, mjd=None, *args, **kwargs): """ Return an AGN lightcurve. The lightcurve for an AGN is returned for the given band and MJD. Parameters ---------- i: int ID of AGN for which SED is returned. band: str EGG-style passband name. mjd: float or array_like or None MJD of observation. """ if mjd is None: mjd = util.get_mjd_vec() Loading Loading @@ -638,5 +761,27 @@ class CatalogAGN: return lc def get_lightcurve_mjd(self, i, band, mjd, mjd0): """ Return an AGN lightcurve at given MJD. The lightcurve is first estimated assuming the timespan starting from 0. The lightcurve is then evaluated at MJD - MJD0 to return the final value of the lightcurve. Parameters ---------- i: int ID of AGN for which SED is returned. band: str EGG-style passband name. mjd: float or array_like MJD of observation. mjd0: float MJD offset corresponding to the first light of the survey. Returns ------- The value of the lightcurve at given MJD in microjanskies. """ lc = self.get_lightcurve(i, band) return np.interp(mjd - mjd0, util.get_mjd_vec(), lc)
src/lsst_inaf_agile/catalog_combined.py +114 −22 Original line number Diff line number Diff line Loading @@ -3,9 +3,7 @@ # Email: akke.viitanen@helsinki.fi # Date: 2023-02-13 16:20:51 """ Create Galaxy+AGN mocks for the LSST Italian AGN in-kind contribution """ """Module for the combined catalog of AGN, galaxies, and stars.""" import logging import os Loading @@ -19,27 +17,65 @@ logger = logging.getLogger(__name__) class CatalogCombined: """Combined catalog class of AGN, galaxies, and stars. Parameters ---------- dirname: str Directory name to contain the catalog. catalog_galaxy: catalog_galaxy.CatalogGalaxy or None Galaxy catalog. catalog_agn: catalog_agn.CatalogAGN or None AGN catalog. catalog_star: catalog_star.CatalogStar or None Single star catalog. catalog_binary: catalog_star.CatalogStar or None Binary star catalog. Attributes ---------- dirname: str Directory name to contain the catalog. catalog_galaxy: catalog_galaxy.CatalogGalaxy or None Galaxy catalog. catalog_agn: catalog_agn.CatalogAGN or None AGN catalog. catalog_star: catalog_star.CatalogStar or None Single star catalog. catalog_binary: catalog_star.CatalogStar or None Binary star catalog. """ def __init__( self, dirname: str, catalog_galaxy=None, catalog_agn=None, catalog_star=None, catalog_binary=None ): self.dirname = dirname self.filename = f"{dirname}/catalog.fits" self.catalog_galaxy = catalog_galaxy self.catalog_agn = catalog_agn self.catalog_star = catalog_star self.catalog_binary = catalog_binary # NOTE: short-circuit for an existing truth catalog if os.path.exists(self.filename): self.catalog_combined = fitsio.read(self.filename) if os.path.exists(f := self.get_filename()): self.catalog_combined = fitsio.read(f) return self.catalog_combined = self.get_catalog_combined() self.catalog_combined = self.postprocess() self.write() def get_filename(self): """Return the FITS filename of the combined catalog""" return f"{self.dirname}/catalog.fits" def get_dtype(self): """ Return the combined catalog dtype constructed from the individual input catalogs. Returns ------- List of (name, type) corresponding to the names and types of the columns. """ dtype = {} for catalog in (self.catalog_galaxy, self.catalog_agn, self.catalog_star, self.catalog_binary): for name in catalog.catalog.dtype.names: Loading @@ -55,7 +91,8 @@ class CatalogCombined: def get_number_galaxy_star_binary(self): """ Return the number of galaxies, stars, and binaries Return the tuple (N_galaxy, N_star, N_binary) where N is the number of objects in the combined catalog. """ return ( len(self.catalog_galaxy.catalog), Loading @@ -64,7 +101,7 @@ class CatalogCombined: ) def get_catalog_combined(self): """ r""" Constructs the combined catalog based on the dtype. For each column in the catalog, the column is looked for in the Loading @@ -72,14 +109,13 @@ class CatalogCombined: are assigned the values from the reference catalog. Note that some columns may appear in multiple catalog. The indices correspond to the following type idx_first idx_last The indices correspond to the following object types: AGN \in [0, n_galaxy[, galaxy \in [0, n_galaxy[, star \in [n_galaxy, n_galaxy + n_star], and binary \in [n_galaxy + n_star, n_total[. AGN 0 n_galaxy galaxy 0 n_galaxy star n_galaxy n_galaxy + n_star binary n_galaxy + n_star n_total Returns ------- The combined catalog as a numpy ndarray. """ n_galaxy, n_star, n_binary = self.get_number_galaxy_star_binary() Loading @@ -104,12 +140,14 @@ class CatalogCombined: # Handles star if name in self.catalog_star.catalog.dtype.names: catalog_combined[name][n_galaxy : n_galaxy + n_star] = self.catalog_star[name] # NOTE: set the flags to False for stars elif catalog_combined.dtype[name] == np.dtype("bool"): catalog_combined[name][n_galaxy : n_galaxy + n_star] = False # Handles binary star if name in self.catalog_binary.catalog.dtype.names: catalog_combined[name][n_galaxy + n_star :] = self.catalog_binary[name] # NOTE: set the flags to False for stars elif catalog_combined.dtype[name] == np.dtype("bool"): catalog_combined[name][n_galaxy + n_star :] = False Loading @@ -131,7 +169,18 @@ class CatalogCombined: return self.catalog_combined def get_flux_total(self, band, rest_frame=False): # NOTE: override rest-frame behavior as EGG does not have absmags for disk/bulge separately """ Return the total flux. The total flux is defined for each type of object as follows: - AGN: sum of galaxy bulge, galaxy disk, and point fluxes - galaxy: sum of galaxy bulge and disk fluxes - star: point fluxes - binary star: sum of the component point fluxes """ # NOTE: override rest-frame behavior as EGG does not have absmags for # disk/bulge separately keys = ["disk", "bulge", "point"] if rest_frame: keys = ["", "point"] Loading @@ -151,19 +200,56 @@ class CatalogCombined: return ret def write(self): fitsio.write(self.filename, self.catalog_combined, clobber=True) """ Writes the combined catalog to a FITS file """ fitsio.write(self.get_filename(), self.catalog_combined, clobber=True) return self.catalog_combined def ingest(self, filename_database, name_table, if_exists="replace"): """ Ingest the truth catalog into a sqlite3 database. Parameters ---------- filename_database: str filename of the sqlite3 database name_table: str name of the table in the sqlite3 database """ import sqlite3 from astropy.table import Table table = Table(self.catalog) df = table.to_pandas() df.reset_index() with sqlite3.connect(filename_database) as con: df.to_sql(name_table, con, index=False, if_exists=if_exists, chunksize=1024, method="multi") def get_is_galaxy(self): """ Return boolean vector selecting galaxies """ return np.isfinite(self.catalog_combined["Z"]) def get_is_star(self): """ Return boolean vector selecting stars """ return ~self.get_is_galaxy() def get_index_star(self, i): """ Return the index in the stellar catalog corresponding to 'i' in the combined catalog combined catalog. Parameters ---------- i: int ID of the truth catalog object """ n_galaxy, n_star, n_binary = self.get_number_galaxy_star_binary() Loading @@ -181,7 +267,7 @@ class CatalogCombined: def write_reference_catalog(self, filename: str, maglim: float, selection_band: str): """ Write a mock LSST reference catalog in csv format to the given Write a mock LSST reference catalog in CSV format to the given 'filename'. The reference catalog is cut to brighter than the magnitude limit 'maglim' in the band 'selection_band'. """ Loading @@ -198,7 +284,7 @@ class CatalogCombined: # Solve the parallax parallax = np.zeros_like(c, dtype=float) select_star = ~np.isfinite(c["Z"]) select_star = self.get_is_star() parallax[select_star] = util.distance_modulus_to_parallax(c["D"][select_star]) # Initialize the data to be written Loading Loading @@ -241,7 +327,7 @@ class CatalogCombined: "pmra_pmdec_corr": np.zeros(c.size), } # Write the reference catalog.csv # Write reference catalog.csv with open(filename, "w") as f: to_write = ",".join(data.keys()) print(f"{to_write}", file=f) Loading @@ -258,9 +344,15 @@ class CatalogCombined: logger.info(f"Wrote reference catalog {filename}") def get_area(self): """ Return the EGG catalog area in deg2 """ cmd = fitsio.read(f"{self.dirname}/egg.fits", columns="CMD")[0] area_string = re.findall("area=([0-9.]+)", cmd)[0] return float(area_string) def __getitem__(self, key): """ Return combined catalog value corresponding to 'key' """ return self.catalog_combined[key]
