Commit e6464b43 authored by Akke Viitanen's avatar Akke Viitanen
Browse files

Fixes the reference catalog

parent b2108e1d

src/python/catalog_combined.py

+23 −13
Original line number Diff line number Diff line
@@ -9,6 +9,7 @@ Create Galaxy+AGN mocks for the LSST Italian AGN in-kind contribution 
"""



import logging

import os



import fitsio

import numpy as np
@@ -22,9 +23,9 @@ class CatalogCombined: 


    def __init__(

        self,

        dirname,

        catalog_agn,

        dirname: str,

        catalog_galaxy,

        catalog_agn,

        catalog_star,

        catalog_binary

    ):
@@ -32,8 +33,8 @@ class CatalogCombined: 
        self.dirname = dirname

        self.filename = f"{dirname}/catalog.fits"



        self.catalog_agn = catalog_agn

        self.catalog_galaxy = catalog_galaxy

        self.catalog_agn = catalog_agn

        self.catalog_star = catalog_star

        self.catalog_binary = catalog_binary
@@ -198,14 +199,20 @@ class CatalogCombined: 
        limit 'maglim' in the band 'selection_band'.

        """



        logger.info("Writing the reference catalog...")

        if os.path.exists(filename):

            logger.info("Reference catalog exists. Will not overwrite.")

            logger.info(f"{filename} exists, will not overwrite.")

            return



        # Cut at some magnitude limit e.g. r < 24 and non-variable sources

        select = util.flux_to_mag(self.get_flux_total(selection_band)) < float(maglim)

        select *= ~self.catalog["is_agn"]

        c = self.catalog[select]

        select *= ~self["is_agn"]

        c = self[select]



        # Solve the parallax

        parallax = np.zeros_like(c, dtype=float)

        select_star = ~np.isfinite(c["Z"])

        parallax[select_star] = util.distance_modulus_to_parallax(c["D"][select_star])



        # Initialize the data to be written

        data = {
@@ -217,7 +224,7 @@ class CatalogCombined: 
            "ra": c["RA"],

            "dec": c["DEC"],



            # Fluxes

            # magnitudes

            "u": util.flux_to_mag(c["lsst-u_total"]),

            "g": util.flux_to_mag(c["lsst-g_total"]),

            "r": util.flux_to_mag(c["lsst-r_total"]),
@@ -231,7 +238,7 @@ class CatalogCombined: 
            "pmra":           np.where(np.isfinite(c["pmracosd"]), c["pmracosd"], 0),

            "pmdec":          np.where(np.isfinite(c["pmdec"]),    c["pmdec"],    0),

            "epoch":          c["RA"].size * ["2000-01-01 12:00:00"],

            "parallax":       np.where(np.isfinite(c["Z"]), 0.0, util.distance_modulus_to_parallax(c["D"])),

            "parallax":       parallax,



            # Errors

            "ra_error":       np.zeros(c.size),
@@ -254,18 +261,21 @@ class CatalogCombined: 


        }



        logger.info("Writing the reference catalog...")

        # Write the reference catalog.csv

        with open(filename, 'w') as f:

            to_write = ','.join(data.keys())

            f.write(f"{to_write}\n")

            print(f"{to_write}", file=f)

            for i in range(c["RA"].size):

                to_write = ','.join([f"{v[i]}" for v in data.values()])

                f.write(f"{to_write}\n")

                print(f"{to_write}", file=f)



        # Run the LSST tools ontop

        dirname = f"{os.path.dirname(filename)}/reference_catalog"

        util.create_directory(dirname)

        os.system(f"convertReferenceCatalog {dirname} {ROOT}/src/python/config_reference_catalog.py {filename}")

        if os.path.exists(dirname):

            os.system(f"rm -rfv {dirname}")

        util.create_directory(dirname + "/")

        os.system(f"convertReferenceCatalog {dirname} src/python/config_reference_catalog.py {filename}")

        logger.info(f"Wrote reference catalog {filename}")



    def __getitem__(self, key):

        return self.catalog_combined[key]

src/python/catalog_galaxy.py

+1 −1
Original line number Diff line number Diff line
@@ -74,7 +74,7 @@ class CatalogGalaxy: 
            ("RA",            np.float64, "Right ascenscion", "deg"),

            ("DEC",           np.float64, "Declination", "deg"),

            ("Z",             np.float64, "Cosmological redshift", ""),

            ("D",             np.float64, "Luminosity distance", "Mpc"),

            ("D",             np.float64, "Luminosity distance OR distance modulus for stars", "Mpc OR mag"),

            ("M",             np.float64, "log10 of stellar mass", "Msun"),

            ("SFR",           np.float64, "Star-formation rate", "Msun/yr"),

            ("PASSIVE",       bool,       "Is passive (non-star-forming)", ""),

src/python/config_reference_catalog.py

+2 −15
Original line number Diff line number Diff line 
from lsst.meas.algorithms import convertRefcatManager

from convert_agile_manager import ConvertAgileManager



#config.manager.retarget(convertRefcatManager.ConvertRefcatManager)

config.manager.retarget(convertRefcatManager.ConvertAgileManager)

config.manager.retarget(ConvertAgileManager)

config.dataset_config.ref_dataset_name = "my-refcat"

config.n_processes = 10



#config.id_name = "source_id"

#config.ra_name = "ra"

#config.dec_name = "dec"

#config.mag_column_list = list("ugrizy")

#

#config.pm_ra_name = "pm_ra"

#config.pm_dec_name = "pm_dec"

#

#config.epoch_name = "epoch"

#config.epoch_format = "iso"

#config.epoch_scale = "utc"



config.full_position_information = True



config.id_name = "source_id"

src/python/util.py

+5 −47
Original line number Diff line number Diff line
@@ -176,46 +176,7 @@ def get_log_L_2500(log_L_2_keV, alpha=0.952, beta=2.138, scatter=True): 
        log_L_2500 += np.random.normal(loc=0, scale=0.4, size=log_L_2500.size)



    return log_L_2500

    #return (10 ** L_2500) * 2.998e18 / 2500





#import sed

#from astropy.cosmology import FlatLambdaCDM

#COSMO = FlatLambdaCDM(Om0=0.30, H0=70)

#def mock_lx_to_M_1450(lx, z, distance_cm, ebv, scatter, seed):

#    l_2_kev = get_log_L_2_keV(lx)

#    l_2500 = get_log_L_2500(l_2_kev, scatter=scatter)

#    wav, flux = sed.get_sed(

#        LogL2500=l_2500,

#        AGN_type=1,

#        ebv=ebv,

#        redshift=z,

#        distance_cm=distance_cm,

#        LogL2kev=l_2_kev,

#        flux_rf_1000_4000_gal=np.inf,

#        seed=seed,

#        wav_min=1449 * u.angstrom,

#        wav_max=1451 * u.angstrom,

#        dlog_wav=1e-4,

#    )[1]

#    flux_1450 = np.interp(1450, wav, flux)

#    M_1450 = flux_to_mag(flux_1450)

#    return M_1450

#

#def mock_lx_to_sed(lx, z, distance_cm, ebv, scatter, seed):

#    l_2_kev = get_log_L_2_keV(lx)

#    l_2500 = get_log_L_2500(l_2_kev, scatter=scatter)

#    wav, flux = sed.get_sed(

#        LogL2500=l_2500,

#        AGN_type=1,

#        ebv=ebv,

#        redshift=z,

#        distance_cm=distance_cm,

#        LogL2kev=l_2_kev,

#        flux_rf_1000_4000_gal=np.inf,

#        seed=seed,

#    )[1]

#    return wav, flux





def get_fraction_obscured(lx):
@@ -791,13 +752,10 @@ def get_galaxy_ab(reff, ratio): 




def create_directory(filename):



    if os.path.isdir(filename):

    if not os.path.isdir(filename):

        filename = os.path.dirname(filename)

    logger.info(f"Creating directory {filename}")

    os.makedirs(filename, exist_ok=True)

        return



    dirname = os.path.dirname(filename)

    os.makedirs(dirname, exist_ok=True)





def get_mjd_vec():

src/scripts/butler/ingest_reference_catalog.sh

+0 −4
Original line number Diff line number Diff line
@@ -7,9 +7,5 @@ if [ $# -ne 2 ] ; then 
	exit 1

fi



#butler register-dataset-type $REPO my-refcat SimpleCatalog htm7

#butler ingest-files -t direct $REPO gaia_dr2 refcats/DM-NNNNN $DIRNAME_REFERENCE_CATALOG/filename_to_htm.ecsv

#butler collection-chain REPO --mode extend refcats refcats/DM-NNNNN



butler register-dataset-type $REPO reference_catalog SimpleCatalog htm7

butler ingest-files -t direct $REPO reference_catalog refcats $DIRNAME_REFERENCE_CATALOG/filename_to_htm.ecsv