from herschelhelp_internal import git_version
print("This notebook was run with herschelhelp_internal version: \n{}".format(git_version()))
%matplotlib inline
#%config InlineBackend.figure_format = 'svg'
import matplotlib.pyplot as plt
plt.rc('figure', figsize=(10, 6))
plt.style.use('ggplot')
import locale
locale.setlocale(locale.LC_ALL, 'en_GB')
import os
import time
import itertools
from astropy.coordinates import SkyCoord
from astropy.table import Table
from astropy import units as u
from astropy import visualization as vis
import numpy as np
from matplotlib_venn import venn3
from herschelhelp_internal.masterlist import nb_compare_mags, nb_ccplots, nb_histograms, find_last_ml_suffix
OUT_DIR = os.environ.get('OUT_DIR', "./data")
SUFFIX = find_last_ml_suffix()
#SUFFIX = "20171016"
master_catalogue_filename = "master_catalogue_egs_{}.fits".format(SUFFIX)
master_catalogue = Table.read("{}/{}".format(OUT_DIR, master_catalogue_filename))
print("Diagnostics done using: {}".format(master_catalogue_filename))
flag_obs = master_catalogue['flag_optnir_obs']
flag_det = master_catalogue['flag_optnir_det']
venn3(
[
np.sum(flag_obs == 4),
np.sum(flag_obs == 2),
np.sum(flag_obs == 6),
np.sum(flag_obs == 1),
np.sum(flag_obs == 5),
np.sum(flag_obs == 3),
np.sum(flag_obs == 7)
],
set_labels=('mid-IR', 'near-IR', 'Optical' ),
subset_label_formatter=lambda x: "{}%".format(int(100*x/len(flag_obs)))
)
plt.title("Wavelength domain observations");
venn3(
[
np.sum(flag_det[flag_obs == 7] == 4),
np.sum(flag_det[flag_obs == 7] == 2),
np.sum(flag_det[flag_obs == 7] == 6),
np.sum(flag_det[flag_obs == 7] == 1),
np.sum(flag_det[flag_obs == 7] == 5),
np.sum(flag_det[flag_obs == 7] == 3),
np.sum(flag_det[flag_obs == 7] == 7)
],
set_labels=('mid-IR', 'near-IR', 'Optical'),
subset_label_formatter=lambda x: "{}%".format(int(100*x/np.sum(flag_det[flag_obs == 7] != 0)))
)
plt.title("Detection of the {} sources detected\n in any wavelength domains\n in regions observed in all domains "
"(among {} total sources)".format(
locale.format('%d', np.sum(flag_det != 0), grouping=True),
locale.format('%d', len(flag_det), grouping=True)));
The master list if composed of several catalogues containing magnitudes in similar filters on different instruments. We are comparing the magnitudes in these corresponding filters.
u_bands = ["Megacam u"]
g_bands = ["Megacam g", "GPC1 g"] #No Suprime g?
r_bands = ["Megacam r", "SUPRIME r", "GPC1 r"]
i_bands = ["Megacam i", "SUPRIME i", "GPC1 i"]
z_bands = ["Megacam z", "SUPRIME z", "GPC1 z"]
y_bands = [ "SUPRIME y", "GPC1 y"]
#TODO look at narrow hst bands. Also many other ugrizy bands
We compare the histograms of the total aperture magnitudes of similar bands.
for bands in [u_bands, g_bands, r_bands, i_bands, z_bands, y_bands]:
colnames = ["m_{}".format(band.replace(" ", "_").lower()) for band in bands]
nb_histograms(master_catalogue, colnames, bands)
We compare one to one each magnitude in similar bands.
for band_of_a_kind in [u_bands, g_bands, r_bands, i_bands, z_bands, y_bands]:
for band1, band2 in itertools.combinations(band_of_a_kind, 2):
basecol1, basecol2 = band1.replace(" ", "_").lower(), band2.replace(" ", "_").lower()
col1, col2 = "m_ap_{}".format(basecol1), "m_ap_{}".format(basecol2)
nb_compare_mags(master_catalogue[col1], master_catalogue[col2],
labels=("{} (aperture)".format(band1), "{} (aperture)".format(band2)))
col1, col2 = "m_{}".format(basecol1), "m_{}".format(basecol2)
nb_compare_mags(master_catalogue[col1], master_catalogue[col2],
labels=("{} (total)".format(band1), "{} (total)".format(band2)))
Cross-match the master list to SDSS and 2MASS to compare its magnitudes to SDSS and 2MASS ones.
master_catalogue_coords = SkyCoord(master_catalogue['ra'], master_catalogue['dec'])
The catalogue is cross-matched to SDSS-DR13 withing 0.2 arcsecond.
We compare the u, g, r, i, and z magnitudes to those from SDSS using fiberMag
for the aperture magnitude and petroMag
for the total magnitude.
sdss = Table.read("../../dmu0/dmu0_SDSS-DR13/data/SDSS-DR13_EGS.fits")
sdss_coords = SkyCoord(sdss['ra'] * u.deg, sdss['dec'] * u.deg)
idx, d2d, _ = sdss_coords.match_to_catalog_sky(master_catalogue_coords)
mask = (d2d < 0.2 * u.arcsec)
sdss = sdss[mask]
ml_sdss_idx = idx[mask]
for band_of_a_kind in [u_bands, g_bands, r_bands, i_bands, z_bands]:
for band in band_of_a_kind:
sdss_mag_ap = sdss["fiberMag_{}".format(band[-1])]
master_cat_mag_ap = master_catalogue["m_ap_{}".format(band.replace(" ", "_").lower())][ml_sdss_idx]
nb_compare_mags(sdss_mag_ap, master_cat_mag_ap,
labels=("SDSS {} (fiberMag)".format(band[-1]), "{} (aperture)".format(band)))
sdss_mag_tot = sdss["petroMag_{}".format(band[-1])]
master_cat_mag_tot = master_catalogue["m_ap_{}".format(band.replace(" ", "_").lower())][ml_sdss_idx]
nb_compare_mags(sdss_mag_ap, master_cat_mag_ap,
labels=("SDSS {} (petroMag)".format(band[-1]), "{} (total)".format(band)))
The catalogue is cross-matched to 2MASS-PSC withing 0.2 arcsecond. We compare the WIRCAM total J and K magnitudes to those from 2MASS.
The 2MASS magnitudes are “Vega-like” and we have to convert them to AB magnitudes using the zero points provided on this page:
Band | Fν - 0 mag (Jy) |
---|---|
J | 1594 |
H | 1024 |
Ks | 666.7 |
2MASS and WIRCAM both use Ks so no conversion is required.
# The AB zero point is 3631 Jy
j_2mass_to_ab = 2.5 * np.log10(3631/1595)
k_2mass_to_ab = 2.5 * np.log10(3631/666.7)
twomass = Table.read("../../dmu0/dmu0_2MASS-point-sources/data/2MASS-PSC_EGS.fits")
twomass_coords = SkyCoord(twomass['raj2000'], twomass['dej2000'])
idx, d2d, _ = twomass_coords.match_to_catalog_sky(master_catalogue_coords)
mask = (d2d < 0.2 * u.arcsec)
twomass = twomass[mask]
ml_twomass_idx = idx[mask]
nb_compare_mags(twomass['jmag'] + j_2mass_to_ab, master_catalogue['m_wircam_j'][ml_twomass_idx],
labels=("2MASS J", "WIRCAM J (total)"))
nb_compare_mags(twomass['kmag'] + k_2mass_to_ab, master_catalogue['m_wircam_k'][ml_twomass_idx],
labels=("2MASS Ks", "WIRCAM Ks (total)"))
From here, we are only comparing sources with a signal to noise ratio above 3, i.e. roughly we a magnitude error below 0.3.
To make it easier, we are setting to NaN in the catalogue the magnitudes associated with an error above 0.3 so we can't use these magnitudes after the next cell.
for error_column in [_ for _ in master_catalogue.colnames if _.startswith('merr_')]:
column = error_column.replace("merr", "m")
keep_mask = np.isfinite(master_catalogue[error_column])
keep_mask[keep_mask] &= master_catalogue[keep_mask][error_column] <= 0.3
master_catalogue[column][~keep_mask] = np.nan
nb_ccplots(
master_catalogue['m_megacam_r'],
master_catalogue['m_ap_megacam_r'] - master_catalogue['m_megacam_r'],
"r total magnitude (CFHT)", "r aperture mag - total mag (CFHT)",
master_catalogue["stellarity"],
invert_x=True
)
nb_ccplots(
master_catalogue['m_suprime_g'] - master_catalogue['m_suprime_i'],
master_catalogue['m_wircam_j'] - master_catalogue['m_wircam_k'],
"g - i (HSC)", "J - Ks (WIRCAM)",
master_catalogue["stellarity"]
)
nb_ccplots(
master_catalogue['m_suprime_i'] - master_catalogue['m_irac-egs_i1'],
master_catalogue['m_suprime_g'] - master_catalogue['m_suprime_i'],
"HSC i - IRAC1", "g - i (HSC)",
master_catalogue["stellarity"]
)
nb_ccplots(
master_catalogue['m_megacam_u'] - master_catalogue['m_megacam_g'],
master_catalogue['m_megacam_g'] - master_catalogue['m_megacam_r'],
"u - g (CFHT)", "g - r (CFHT)",
master_catalogue["stellarity"]
)
nb_ccplots(
master_catalogue['m_wircam_j'] - master_catalogue['m_wircam_k'],
master_catalogue['m_suprime_g'] - master_catalogue['m_wircam_j'],
"J - K (WIRCAM)", "g - J (HSC, WIRCAM)",
master_catalogue["stellarity"]
)
nb_ccplots(
master_catalogue['m_megacam_i'] - master_catalogue['m_megacam_z'],
master_catalogue['m_megacam_z'] - master_catalogue['m_wircam_j'],
"i - z (Megacam)", "z - J (Megacam, Wircam)",
master_catalogue["stellarity"]
)
nb_ccplots(
master_catalogue['m_irac-egs_i3'] - master_catalogue['m_irac-egs_i4'],
master_catalogue['m_irac-egs_i1'] - master_catalogue['m_irac-egs_i2'],
"IRAC3 - IRAC4", "IRAC1 - IRAC2",
master_catalogue["stellarity"]
)