XMM-13hr: Validation Checks (FULL)¶
In [1]:
%matplotlib inline
#%config InlineBackend.figure_format = 'svg'
import matplotlib.pyplot as plt
plt.rc('figure', figsize=(10, 6))
plt.style.use('ggplot')
import numpy as np
from astropy.table import Table
import itertools
import time
In [2]:
t0 = time.time()
catname = "/data/help/master_catalogue_xmm-13hr_20180201.fits"
master_catalogue = Table.read(catname)
print('Elapsed time(secs): ', time.time() - t0)
print("Number of sources in master catalogue: ", len(master_catalogue))
In [3]:
field = master_catalogue["field"][0]
field = field.rstrip() # remove whitespaces at the end of the sting
print(field)
In [4]:
%%javascript
IPython.OutputArea.auto_scroll_threshold = 9999;
In [5]:
u_bands = []
g_bands = ["BASS g"]
r_bands = ["BASS r"]
i_bands = []
z_bands = ["BASS z"]
y_bands = []
J_bands = [ "WFCAM J"]
H_bands = []
K_bands = []
all_bands = []
irac_bands = []
opt_mags = u_bands + g_bands + r_bands + i_bands + z_bands + y_bands
ir_mags = J_bands + H_bands + K_bands + irac_bands
other_bands = []
all_mags = opt_mags + ir_mags + other_bands
1/ Magnitude errors¶
In [6]:
def mag_vs_err(x, y, fig, ax, labels=("x", "y"), savefig=False):
x_label, y_label = labels
print(x_label)
# Use only finite values
mask = np.isfinite(x) & np.isfinite(y) & (x!=99.) & (y!=99.)
x = np.copy(x[mask])
y = np.copy(y[mask])
if len(x) > 0:
print(" Error max: {:.0f}".format(np.max(y)))
err10 = y > 10
if len(x[err10]) > 0:
print(" magerr > 10: Number of objects = {:d}, min mag = {:.1f}".format(len(x[err10]), np.min(x[err10])))
else:
print(" magerr > 10: Number of objects = {:d}, min mag = {:.1f}".format(len(x[err10]), np.nan))
err100 = y > 100
if len(x[err100]) > 0:
print(" magerr > 100: Number of objects = {:d}, min mag = {:.1f}".format(len(x[err100]), np.min(x[err100])))
else:
print(" magerr > 100: Number of objects = {:d}, min mag = {:.1f}".format(len(x[err100]), np.nan))
err1000 = y > 1000
if len(x[err1000]) > 0:
print(" magerr > 1000: Number of objects = {:d}, min mag = {:.1f}".format(len(x[err1000]), np.min(x[err1000])))
else:
print(" magerr > 1000: Number of objects = {:d}, min mag = {:.1f}".format(len(x[err1000]), np.nan))
else:
print(" no data")
print("")
# Plot
ax.set_yscale('log') # to place before scatter to avoid issues
ax.scatter(x, y, marker='.', alpha=0.1, s=50)
ax.invert_xaxis()
#ax.get_yaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
#ax.get_xaxis().get_major_formatter().labelOnlyBase = False
ax.set_xlabel(labels[0])
ax.set_ylabel(labels[1])
# Save ex. fig
if savefig:
survey_label = ((x_label.replace(" ", "_")).replace("(", "")).replace(")", "")
figname = field + "_magVSmagerr_" + survey_label + ".png"
plt.savefig("/data/help/plots/" + figname, bbox_inches='tight')
#plt.show()
In [7]:
for mag in all_mags:
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(16, 6))
basecol = mag.replace(" ", "_").lower()
if basecol == "bass_g":
savefig = True
else:
savefig=False
col, ecol = "m_ap_{}".format(basecol), "merr_ap_{}".format(basecol)
mag_vs_err(master_catalogue[col], master_catalogue[ecol], fig, ax1,
labels=("{} mag (aperture)".format(mag), "{} magerr (aperture)".format(mag)), savefig=False)
col, ecol = "m_{}".format(basecol), "merr_{}".format(basecol)
mag_vs_err(master_catalogue[col], master_catalogue[ecol], fig, ax2,
labels=("{} mag (total)".format(mag), "{} magerr (total)".format(mag)), savefig=savefig)
display(fig)
plt.close()
2/ Magnitude and error issues¶
3/ Outliers¶
$$chi^2 = \frac{(mag_{1}-mag_{2})^2}{magerr_{1}^2 + magerr_{2}^2}$$
Interquartile range (IQR) and outliers:
We consider as outliers objects which have a high $chi^2$, about $5\sigma$ away from the mean.
$25th, 75th \;percentile = 0.6745\sigma$
$IQR = (75th \;percentile - 25th \;percentile) = 0.6745\sigma * 2 = 1.349\sigma$
$75th \;percentile + 3.2\times IQR = 0.6745\sigma + 3.2\times1.349\sigma = 5\sigma$
$$outliers == [chi^2 > (75th \;percentile + 3.2\times (75th \;percentile - 25th \;percentile))]$$
NB:
Bright sources tend to have their errors underestimated with values as low as $10^{-6}$, which is unrealistic. So to avoid high $chi^2$ due to unrealistic small errors, we clip the error to get a minimum value of 0.1% (i.e. all errors smaller then $10^{-3}$ are set to $10^{-3}$).
In [9]:
def outliers(x, y, xerr, yerr, labels=["x", "y"], savefig=False):
import matplotlib
import matplotlib.gridspec as gridspec
from astropy import visualization as vz
fig = plt.figure(figsize=(13, 6))
gs1 = gridspec.GridSpec(1, 1)
gs1.update(left=0.05, right=0.4, wspace=0.05)
ax1 = plt.subplot(gs1[:, :-1])
gs2 = gridspec.GridSpec(1, 3)
gs2.update(left=0.47, right=0.98, hspace=0.05, wspace=0.05)
ax2 = plt.subplot(gs2[:, :-1])
ax3 = plt.subplot(gs2[:, -1], sharey=ax2)
# Use only finite values
mask = np.isfinite(x) & np.isfinite(y) & np.isfinite(xerr) & np.isfinite(yerr)
x = np.copy(x[mask])
y = np.copy(y[mask])
xerr = np.copy(xerr[mask])
yerr = np.copy(yerr[mask])
# mag1 - mag2
diff = y - x
x_label, y_label = labels
# If the difference is all NaN there is nothing to compare.
if np.isnan(diff).all():
print("No sources have both {} and {} values.".format(
x_label, y_label))
print("")
return
# Set the minimum error to 10^-3
np.clip(xerr, 1e-3, np.max(xerr), out=xerr)
np.clip(yerr, 1e-3, np.max(yerr), out=yerr)
# Median, Median absolute deviation and 1% and 99% percentiles
diff_median = np.median(diff)
diff_mad = np.median(np.abs(diff - diff_median))
diff_1p, diff_99p = np.percentile(diff, [1., 99.])
diff_25p, diff_75p = np.percentile(diff, [25., 75.])
diff_label = "{} - {}".format(y_label, x_label)
print("{} ({} sources):".format(diff_label, len(x)))
print("- Median: {:.2f}".format(diff_median))
print("- Median Absolute Deviation: {:.2f}".format(diff_mad))
print("- 1% percentile: {}".format(diff_1p))
print("- 99% percentile: {}".format(diff_99p))
# Chi2 (Normalized difference)
ichi2 = np.power(diff, 2) / (np.power(xerr, 2) + np.power(yerr, 2))
# Use only non-null values of ichi2
mask2 = ichi2 != 0.0
diff, ichi2 = np.copy(diff[mask2]), np.copy(ichi2[mask2])
x, y, xerr, yerr = np.copy(x[mask2]), np.copy(y[mask2]), np.copy(xerr[mask2]), np.copy(yerr[mask2])
# Outliers (5sigma)
log_ichi2_25p, log_ichi2_75p = np.percentile(np.log10(ichi2), [25., 75.])
out_lim = log_ichi2_75p + 3.2*abs(log_ichi2_25p-log_ichi2_75p)
outliers = np.log10(ichi2) > out_lim
nb_outliers = len(x[outliers])
print("Outliers separation: log(chi2) = {:.2f}".format(out_lim))
print("Number of outliers: {}".format(nb_outliers))
print("")
# Comparing mag
ax1.scatter(x, y, marker='.', alpha=0.1, s=50)
ax1.scatter(x[outliers], y[outliers], marker='.', c='b', alpha=0.3, s=50, label='Outliers ({})'.format(nb_outliers))
min_val = np.min(np.r_[x, y])
max_val = np.max(np.r_[x, y])
ax1.autoscale(False)
ax1.plot([min_val, max_val], [min_val, max_val], "k:")
ax1.invert_xaxis()
ax1.invert_yaxis()
ax1.set_xlabel(x_label)
ax1.set_ylabel(y_label)
ax1.legend(loc='lower right', numpoints=1)
# Chi2 vs Diff
#ax1.set_yscale('log') # to place before scatter to avoid issues
ax2.scatter(diff, np.log10(ichi2), marker='.', alpha=0.1, s=50)
if nb_outliers != 0:
ax2.scatter(diff[outliers], np.log10(ichi2[outliers]), marker='.', alpha=0.3, s=50, color='b',\
label='Outliers ({})'.format(nb_outliers))
ax2.axhline(out_lim, color='grey', linestyle=':')
ax2.set_xlabel(diff_label)
ax2.set_ylabel('log(chi2)')
ax2.legend(loc='lower right', numpoints=1)
# Hist
n, bins, patches = vz.hist(np.log10(ichi2), ax=ax3, bins='knuth', facecolor='red', lw = 2, alpha=0.5,\
orientation="horizontal")
if nb_outliers > 3:
n, bins, patches = vz.hist(np.log10(ichi2[outliers]), ax=ax3, bins='knuth', facecolor='b', lw = 2, alpha=0.7,\
orientation="horizontal")
ax3.axhline(out_lim, color='grey', linestyle=':')
ax3.yaxis.set_tick_params(labelleft=False)
# Save ex. fig
if savefig:
survey_label = ((diff_label.replace(" ", "_")).replace("(", "")).replace(")", "")
figname = field + "_outliers_" + survey_label + ".png"
plt.savefig("/data/help/plots/" + figname, bbox_inches='tight')
display(fig)
plt.close()
In [10]:
for band_of_a_kind in all_bands:
for band1, band2 in itertools.combinations(band_of_a_kind, 2):
basecol1, basecol2 = band1.replace(" ", "_").lower(), band2.replace(" ", "_").lower()
if basecol1 == "gpc1_g" and basecol2 == "bass_g":
savefig = True
else:
savefig = False
# Aperture mag
col1, col2 = "m_ap_{}".format(basecol1), "m_ap_{}".format(basecol2)
ecol1, ecol2 = "merr_ap_{}".format(basecol1), "merr_ap_{}".format(basecol2)
outliers(master_catalogue[col1], master_catalogue[col2],
master_catalogue[ecol1], master_catalogue[ecol2],
labels=("{} (aperture)".format(band1), "{} (aperture)".format(band2)), savefig=savefig)
# Tot mag
col1, col2 = "m_{}".format(basecol1), "m_{}".format(basecol2)
ecol1, ecol2 = "merr_{}".format(basecol1), "merr_{}".format(basecol2)
outliers(master_catalogue[col1], master_catalogue[col2],
master_catalogue[ecol1], master_catalogue[ecol2],
labels=("{} (total)".format(band1), "{} (total)".format(band2)))
4. Aperture correction issues¶
In [11]:
def apcor_check(x, y, stellarity, labels=["x", "y"], savefig=False):
import matplotlib.gridspec as gridspec
from astropy import visualization as vz
#fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(16, 6)) #
fig = plt.figure(figsize=(13, 6))
gs1 = gridspec.GridSpec(1, 1)
gs1.update(left=0.05, right=0.4, wspace=0.05)
ax1 = plt.subplot(gs1[:, :-1])
gs2 = gridspec.GridSpec(1, 3)
gs2.update(left=0.47, right=0.98, hspace=0.05, wspace=0.05)
ax2 = plt.subplot(gs2[:, :-1])
ax3 = plt.subplot(gs2[:, -1], sharey=ax2)
# Use only finite values
mask = np.isfinite(x) & np.isfinite(y) & np.isfinite(stellarity)
x = np.copy(x[mask])
y = np.copy(y[mask])
stellarity = np.copy(stellarity[mask])
diff = y - x
x_label, y_label = labels
# If the difference is all NaN there is nothing to compare.
if np.isnan(diff).all():
print("No sources have both {} and {} values.".format(
x_label, y_label))
print("")
return
diff_label = "{} - {}".format(y_label, x_label)
print("{}:".format(diff_label))
# Subsample
zoom = (x > 16) & (x < 20)
# Comparing mag
ax1.scatter(x, diff, marker='.', alpha=0.1, s=50)
ax1.invert_xaxis()
ax1.set_ylabel(diff_label)
ax1.set_xlabel(x_label)
# Zoom Plot
y_min, y_max = np.percentile(diff[zoom], [1., 99.])
y_delta = .1 * (y_max - y_min)
y_min -= y_delta
y_max += y_delta
if len(x[zoom]) < 1000:
alpha = 0.4
else:
alpha = 0.1
print(len(x[zoom]))
pl = ax2.scatter(x[zoom], diff[zoom], marker='.', alpha=alpha, s=50, c=stellarity[zoom], cmap="jet")
ax2.invert_xaxis()
ax2.set_ylabel(diff_label)
ax2.set_xlabel(x_label)
ax2.set_ylim([y_min, y_max])
fig.colorbar(pl, label="stellarity (1=star)")
#ax2.legend(loc='lower right', numpoints=1)
# Hist
n, bins, patches = vz.hist(diff[zoom], ax=ax3, bins='knuth', facecolor='black', lw = 2, alpha=0.5,\
orientation="horizontal")
ax3.yaxis.set_tick_params(labelleft=False)
# Save ex. fig
if savefig:
survey_label = ((diff_label.replace(" ", "_")).replace("(", "")).replace(")", "")
figname = field + "_apcorrIssues_" + survey_label + ".png"
plt.savefig("/data/help/plots/" + figname, bbox_inches='tight')
display(fig)
plt.close()
In [12]:
for band_of_a_kind in all_bands:
for band1, band2 in itertools.combinations(band_of_a_kind, 2):
basecol1, basecol2 = band1.replace(" ", "_").lower(), band2.replace(" ", "_").lower()
if basecol1 == "gpc1_r" and basecol2 == "wfc_r":
savefig = True
else:
savefig = False
# Aperture mag
col1, col2 = "m_ap_{}".format(basecol1), "m_ap_{}".format(basecol2)
apcor_check(master_catalogue[col1], master_catalogue[col2], master_catalogue['stellarity'],
labels=("{} (aperture)".format(band1), "{} (aperture)".format(band2)), savefig=savefig)
for j in range(10):
plt.close()