#!/usr/bin/env python3 """Turn the raw downloads into three standardised early-type-galaxy tables in data/: sample_6df.csv 6dFGSv J-band FP sample (Campbell+ 2014), WITHOUT the chi^2 / selection-probability clipping sample_sdss.csv our own SDSS DR17 Main-Galaxy-Sample early-type sample, no FP clipping of any kind sample_howlett.csv SDSS PV catalogue (Howlett+ 2022), which HAS been clipped of FP outliers (cross-check only) Every physical radius is computed here from the redshift distance (flat LCDM, H0=70, Om=0.3), never from an FP distance. Standard columns: id, ra, dec, zdist, logs, e_logs, logr [kpc], logi [Lsun/pc^2], msel (magnitude the survey selected on), mlim, snr, mu_obs (observed mean surface brightness). """ import glob import os import numpy as np import pandas as pd import fpk_common as F MSUN_R = 4.65 # AB, as used by Howlett et al. 2022 EVO_Q = 0.85 # mag per unit redshift luminosity evolution (Howlett et al. 2022) SDSS_MLIM = 17.77 # Main Galaxy Sample r-band Petrosian limit (extinction corrected) SIXDF_JLIM = 13.65 def build_6df(): d = pd.read_csv(os.path.join(F.RAW, "6df_FPsample.csv")) t2 = pd.read_csv(os.path.join(F.RAW, "6df_table2.csv")) # table2 has repeat spectra: keep the highest-S/N one per galaxy (Campbell+ adopt the best measurement) t2 = t2.sort_values("S/N", ascending=False).drop_duplicates("6dFGS") m = d.merge(t2[["6dFGS", "S/N", "R", "Vd", "e_Vd"]], on="6dFGS", how="left") code = m.Jcode.astype(int).astype(str).str.zfill(6) for k, nm in enumerate(["ok_sigma", "ok_cz", "ok_morph", "ok_mag", "ok_selprob", "ok_chi2"]): m[nm] = code.str[k].astype(int) out = pd.DataFrame({ "id": m["6dFGS"], "twomass": m["2MASX"], "ra": m.RAJ2000 * 15.0, "dec": m.DEJ2000, # VizieR column is in hours "cz_cmb": m.cz, "cz_group": m.cz_gr, }) out["zdist"] = np.where(m.cz_gr.notna() & (m.cz_gr > 0), m.cz_gr, m.cz) / F.C_KMS out["z"] = m.cz / F.C_KMS out["logs"] = m.logVd out["e_logs"] = m.e_logVd out["theta_e"] = 10 ** m.Jlogr # arcsec, circularised, PSF corrected out["logr"] = m.Jlogr + np.log10(F.angular_diameter_distance(out.zdist) * 1e3 / 206265.0) out["logr_cat_h"] = m.JlogRe # catalogue value (kpc/h) for a consistency check out["logi"] = m.JlogIe out["e_logi"] = m.e_JlogIe out["msel"] = m.Jtot out["mlim"] = SIXDF_JLIM out["Htot"] = m.Htot out["Ktot"] = m.Ktot out["HlogIe"] = m.HlogIe out["KlogIe"] = m.KlogIe out["Hlogr"] = m.Hlogr out["Klogr"] = m.Klogr out["mu_obs"] = m.Jtot + 2.5 * np.log10(2 * np.pi * (10 ** m.Jlogr) ** 2) out["snr"] = m["S/N"] out["Rtd"] = m["R"] out["sersic_n"] = m.n out["Mt"] = m.Mt out["group"] = m.Group.fillna(-1).astype(int) out["Nr"] = m.Nr.fillna(-1).astype(int) for nm in ["ok_sigma", "ok_cz", "ok_morph", "ok_mag", "ok_selprob", "ok_chi2"]: out[nm] = m[nm] out["in_6dFGSv_J"] = m.Js # --- sample definition: the 6dFGSv criteria 1-4 only (sigma, cz, morphology, magnitude); # criteria 5 (selection probability) and 6 (chi^2<12 FP-outlier clipping) are deliberately NOT applied out["in_sample"] = ((out.ok_sigma == 1) & (out.ok_cz == 1) & (out.ok_morph == 1) & (out.ok_mag == 1)).astype(int) out["clipped_by_catalogue"] = ((out.in_sample == 1) & (m.Js == 0)).astype(int) out.to_csv(os.path.join(F.DATA, "sample_6df.csv"), index=False) print("6dF: rows %d, in_sample %d, of which clipped by the published catalogue %d" % (len(out), out.in_sample.sum(), out.clipped_by_catalogue.sum())) return out def read_howlett(): fn = os.path.join(F.RAW, "sdss_pv", "SDSS_PV_public.dat") cols = open(fn).readline().lstrip("#").split() h = pd.read_csv(fn, sep=r"\s+", header=None, skiprows=1, names=cols, dtype={"objid": "int64"}) return h def kcorr_poly_terms(z, gr): gr = np.clip(gr, 0.55, 1.15) return np.column_stack([z, z * gr, z * gr ** 2, z ** 2, z ** 2 * gr, z ** 2 * gr ** 2, z ** 3, z ** 3 * gr, z * gr ** 3]) def fit_kcorr(h): """Empirical polynomial reproduction of the Chilingarian+ k_r(z, g-r) values tabulated by Howlett+ 2022.""" gr = (h.deVMag_g - h.extinction_g) - (h.deVMag_r - h.extinction_r) A = kcorr_poly_terms(h.zhelio.values, gr.values) c = np.linalg.lstsq(A, h.kcor_r.values, rcond=None)[0] rms = np.std(h.kcor_r.values - A @ c) print("k-correction polynomial: rms residual vs Howlett values %.4f mag" % rms) return c def build_howlett(h): out = pd.DataFrame({"id": h.objid, "ra": h.RA, "dec": h.Dec, "z": h.zcmb, "zdist": h.zcmb_group}) out["zdist"] = np.where(out.zdist > 0, out.zdist, out.z) theta = h.deVRad_r * np.sqrt(h.deVAB_r) out["theta_e"] = theta out["logs"] = h.s out["e_logs"] = h.es out["logr"] = np.log10(theta) + np.log10(F.angular_diameter_distance(out.zdist) * 1e3 / 206265.0) out["logi"] = h.i out["msel"] = h.deVMag_r out["mlim"] = 17.0 out["mu_obs"] = h.deVMag_r + 2.5 * np.log10(2 * np.pi * theta ** 2) out["logdist"] = h.logdist out["logdist_corr"] = h.logdist_corr out["group"] = h.IDgroupT17 out["Ngroup"] = h.NgroupT17 out["plate"] = h.plate out["in_sample"] = 1 out.to_csv(os.path.join(F.DATA, "sample_howlett.csv"), index=False) print("Howlett+2022: rows %d" % len(out)) return out SDSS_CUT_ORDER = ["mgs", "mag", "z", "dev", "conc", "ba", "colour", "sigma", "snr", "size", "d4000"] def sdss_colour_edges(Mr): """Blue and red edge of the g-r colour window of the sample definition (slides with absolute magnitude).""" return 0.63 - 0.02 * (Mr + 20), 1.03 - 0.02 * (Mr + 20) def sdss_cut_masks(out): """The cuts that define the SDSS sample, as boolean Series on the parent-derived table. Also used by revision_checks.py to build the 'red-extended' sample (colour window without its red edge).""" lo, hi = sdss_colour_edges(out.Mr) cut = {} cut["mgs"] = out.mgs == 1 cut["mag"] = (out.msel <= SDSS_MLIM) & (out.msel > 10) cut["z"] = (out.z >= 0.02) & (out.z <= 0.12) cut["dev"] = (out.fracdev >= 0.8) & (out.dlnL > 0) cut["conc"] = out.conc >= 2.6 cut["ba"] = out.ba >= 0.3 cut["colour"] = (out.gr >= lo) & (out.gr <= hi) cut["colour_blue_edge_only"] = out.gr >= lo cut["sigma"] = (out.veldisp >= 70) & (out.veldisp <= 420) & (out.e_logs > 0) cut["snr"] = out.snr >= 10 cut["size"] = (out.theta_e > 0.5) & np.isfinite(out.logi) & np.isfinite(out.logr) cut["d4000"] = out.d4000 >= 1.6 return cut def build_sdss(h, kc): fs = sorted(glob.glob(os.path.join(F.RAW, "sdss_own", "sdss_etg_ra*.csv"))) parts = [] for f in fs: try: x = pd.read_csv(f, comment="#") except Exception: continue if len(x) and "specObjID" in x.columns: parts.append(x) d = pd.concat(parts, ignore_index=True).drop_duplicates("specObjID") n_parent = len(d) zc = F.z_cmb_from_helio(d.z.values, d.ra.values, d["dec"].values) theta = d.deVRad_r * np.sqrt(d.deVAB_r) # circularised de Vaucouleurs radius, arcsec gr_dev = (d.deVMag_g - d.extinction_g) - (d.deVMag_r - d.extinction_r) kr = kcorr_poly_terms(d.z.values, gr_dev.values) @ kc m_r = d.deVMag_r - d.extinction_r - kr + EVO_Q * d.z with np.errstate(divide="ignore", invalid="ignore"): mu_e = m_r + 2.5 * np.log10(2 * np.pi * theta ** 2) - 10 * np.log10(1 + zc) logi = 0.4 * (MSUN_R - mu_e) + 8.629 logs = np.log10(d.velDisp) + 0.04 * np.log10(1.5 / (theta / 8.0)) logs_port = np.log10(d.sigmaStars.where(d.sigmaStars > 0)) + 0.04 * np.log10(1.5 / (theta / 8.0)) logr = np.log10(theta) + np.log10(F.angular_diameter_distance(zc) * 1e3 / 206265.0) petro = d.petroMag_r - d.extinction_r Mr = m_r - 5 * np.log10(F.luminosity_distance(zc) * 1e5) gr_mod = (d.modelMag_g - d.extinction_g) - (d.modelMag_r - d.extinction_r) out = pd.DataFrame({ "id": d.bestObjID, "specObjID": d.specObjID, "plate": d.plate, "mjd": d.mjd, "fiber": d.fiberID, "ra": d.ra, "dec": d["dec"], "z": zc, "zdist": zc, "zhelio": d.z, "logs": logs, "e_logs": d.velDispErr / d.velDisp / np.log(10), "logs_port": logs_port, "veldisp": d.velDisp, "sigma_port": d.sigmaStars, "theta_e": theta, "logr": logr, "logi": logi, "Mr": Mr, "msel": petro, "mlim": SDSS_MLIM, "devmag_r": d.deVMag_r, "kcor_r": kr, "mu50": petro + 2.5 * np.log10(2 * np.pi * d.petroR50_r ** 2), "mu_obs": d.deVMag_r + 2.5 * np.log10(2 * np.pi * theta ** 2), "snr": d.snMedian, "fibermag_r": d.fiberMag_r, "conc": d.petroR90_r / d.petroR50_r, "ba": d.deVAB_r, "fracdev": d.fracDeV_r, "dlnL": d.lnLDeV_r - d.lnLExp_r, "gr": gr_mod, "gr_dev": gr_dev, "petroR50": d.petroR50_r, "devrad": d.deVRad_r, "devraderr": d.deVRadErr_r, "dmag_dev_petro": d.deVMag_r - d.petroMag_r, "ew_ha_port": d.EW_Ha_6562, "ew_ha_mpa": -d.h_alpha_eqw, "d4000": d.d4000_n, "bpt": d.bpt, "clean": d.clean, "nchild": d.nChild, "has_parent": (d.parentID > 0).astype(int), "mgs": ((d.legacy_target1 & 64) > 0).astype(int), "ext_r": d.extinction_r, }) out["in_howlett"] = out.id.isin(set(h.objid.values)).astype(int) # ---------------- sample definition (see report.md, Data) ---------------- cut = sdss_cut_masks(out) keep = np.ones(len(out), bool) rows = [("parent download (GALAXY, zWarning=0, primary, legacy, 0.003350, fracDeV_r>=0.5)", n_parent)] for k in SDSS_CUT_ORDER: keep &= cut[k].values rows.append((k, int(keep.sum()))) out["in_sample"] = keep.astype(int) pd.DataFrame(rows, columns=["cut", "n_remaining"]).to_csv(os.path.join(F.DATA, "sdss_cut_ladder.csv"), index=False) for r in rows: print(" SDSS cut %-12s -> %d" % (r[0][:12], r[1])) # these tables are large (60-100 MB): keep them in the cache next to the raw downloads out.to_csv(os.path.join(F.CACHE, "sdss_own_parent_derived.csv"), index=False) out[out.in_sample == 1].to_csv(F.sample_path("sdss"), index=False, float_format="%.8g") return out def main(): build_6df() h = read_howlett() kc = fit_kcorr(h) build_howlett(h) build_sdss(h, kc) if __name__ == "__main__": main()