"""Step 7: merge the two samples into one list of unique disc galaxies with a single residual each. Priority: Cosmicflows-4 measurement where the galaxy is in the CF4 sample (by-eye inclination, re-measured photometry), otherwise the ALFALFA-SDSS gold measurement. Residual statistic D = least-underluminous band residual (dex) about the luminosity Tully-Fisher relation. Also computed for every galaxy: * automated quality flags (same-redshift SDSS DR12 spectroscopic neighbour within 2', a second a.100 detection within 4', HI centroid > 40" from the optical counterpart, line-width error > 20 km/s, a same-velocity SIMBAD galaxy within 5'), * a selection weight w_sel_f in [0, 1]: would the galaxy still lie inside the range of apparent magnitude AND mean surface brightness that its catalogue actually contains if it were dimmed by 0.75/1.5/2.5 mag? (limits = 98th percentile of each catalogue and band; 95th and 99.5th percentiles kept as variants; galaxies without a size measurement get the pass rate of comparable galaxies with one), * its baryonic Tully-Fisher residual if a fraction f of its starlight were removed, * the ratio nu L_nu(W1) / nu L_nu(i), needed to ask how much waste heat would fill in the W1 band. Outputs: data/union_sample.csv, data/selection_limits.csv """ import os import numpy as np import pandas as pd from scipy.spatial import cKDTree from common import RAW, DATA, C_KMS, radec_to_xyz, chord_for_arcsec, NULNU_SUN_I, NULNU_SUN_W1 from s05_simbad import simbad_neighbours, GALAXY_OTYPES FS = (0.5, 0.75, 0.9) PCTS = (98.0, 95.0, 99.5) # percentile defining the faint end of each catalogue; the first is the baseline E_W_MAX = 20.0 # km/s SIMBAD_RMIN, SIMBAD_RMAX = 0.4, 5.0 # arcmin; inner radius excludes the target's own SIMBAD entries BTF_CUT = -0.35 # dex; ~2 sigma below the baryonic relation counts as "baryon deficient" X0 = 2.2 def spec_neighbours(ra, dec, cz, w): tp = pd.read_csv(os.path.join(RAW, "tempel_dr12.csv")) tree = cKDTree(radec_to_xyz(tp.RAJ2000, tp.DEJ2000)) tcz = tp.zobs.values * C_KMS xyz = radec_to_xyz(ra, dec) nb = tree.query_ball_point(xyz, chord_for_arcsec(120.0)) rself = chord_for_arcsec(20.0) out = np.zeros(len(ra), int) for i, lst in enumerate(nb): if lst: lst = np.asarray(lst) sep = np.linalg.norm(tree.data[lst] - xyz[i], axis=1) out[i] = np.sum((np.abs(tcz[lst] - cz[i]) < w[i] / 2 + 200.0) & (sep >= rself)) return out def main(): cf = pd.read_csv(os.path.join(DATA, "cf4_resid.csv")) al = pd.read_csv(os.path.join(DATA, "alfalfa_resid.csv"), low_memory=False) al = al[al.gold].copy() bfit_cf = None # ---- CF4 part u1 = pd.DataFrame(dict( uid="P" + cf.PGC.astype(str), src="cf4", PGC=cf.PGC, AGC=cf.AGC, name=cf.Name.astype(str).str.strip(), ra=cf.ra, dec=cf.dec, cz=cf.HRV, W=cf.Wmx, e_W=cf.e_Wmx, inc=cf.Inc, V=cf.V, logV=cf.logV, Dist=cf.Dist, r_i=cf.r_i, r_W1=cf.r_W1, D=cf.D, r_bar=cf.r_bar, logMs=cf.logMs, logMHI=cf.logMHI, logMbar=cf.logMbar, fgas=cf.fgas, m_i=cf.i_star.where(cf.logLi.notna()), m_W1=cf.W1_star.where(cf.logLW1.notna()), in_cf4_table4=cf.in_table4, dmu=cf.dmu, mu_i=cf.mu_i.where(cf.logLi.notna()), mu_W1=cf.mu_W1.where(cf.logLW1.notna()), D_rm=cf.D_rm, # does the galaxy have an a.100 detection (3.5' Arecibo beam)? If not, the CF4 width necessarily comes # from another telescope (GBT ~9', Effelsberg ~9', Parkes ~14', Nancay 4'x22') or from older Arecibo # data, and the 5' SIMBAD neighbour flag does not cover the whole beam. in_a100=(cf.AGC > 0).astype(int), # nu L_nu(W1) / nu L_nu(i) from the galaxy's own photometry (NaN if one band is missing) w1_over_i=10 ** (cf.logLW1 - cf.logLi) * NULNU_SUN_W1 / NULNU_SUN_I)) u1["n_spec_nb"] = spec_neighbours(cf.ra.values, cf.dec.values, cf.HRV.values, cf.Wmx.values) alall = pd.read_csv(os.path.join(DATA, "alfalfa_all.csv"), low_memory=False, usecols=["AGC", "n_hi_nb", "hi_opt_off", "SNR"]).set_index("AGC") u1["n_hi_nb"] = u1.AGC.map(alall.n_hi_nb).fillna(0).astype(int) u1["hi_opt_off"] = u1.AGC.map(alall.hi_opt_off) u1["D_other"] = u1.AGC.map(al.set_index("AGC").D) # ALFALFA-SDSS value for the same galaxy, if gold # ---- ALFALFA-only part (gold galaxies that are not members of the CF4 sample) al2 = al[~al.AGC.isin(cf.AGC[cf.AGC > 0])].copy() u2 = pd.DataFrame(dict( uid="A" + al2.AGC.astype(str), src="alfalfa", PGC=-1, AGC=al2.AGC, name=al2.Name.astype(str).str.strip(), ra=al2.RAO, dec=al2.DEO, cz=al2.Vhel, W=al2.W50, e_W=al2.e_W50, inc=np.degrees(np.arcsin(al2.sini)), V=al2.V, logV=al2.logV, Dist=al2.Dist, r_i=al2.r_i, r_W1=al2.r_W1, D=al2.D, r_bar=al2.r_bar, logMs=al2.logMsM, logMHI=al2.logMHI, logMbar=al2.logMbar, fgas=al2.fgas, m_i=al2.imag, m_W1=np.nan, in_cf4_table4=False, n_spec_nb=al2.n_spec_nb, n_hi_nb=al2.n_hi_nb, hi_opt_off=al2.hi_opt_off, D_other=np.nan, dmu=np.nan, mu_i=al2.mu_i, mu_W1=np.nan, D_rm=al2.D_rm, in_a100=1, # logLW1 here is the McGaugh stellar mass = 0.45 x L_W1 [solar W1 units] w1_over_i=10 ** (al2.logLW1 - np.log10(0.45) - al2.logLi) * NULNU_SUN_W1 / NULNU_SUN_I)) u2["in_alfalfa_gold"] = True u1["in_alfalfa_gold"] = u1.AGC.isin(al.AGC) # uniformly HI-selected subsample (CF4 measurement used) u = pd.concat([u1, u2], ignore_index=True) # ---- same-velocity galaxies in SIMBAD (all-sky; matters for the large-beam widths used by CF4) nb = simbad_neighbours(u.uid, u.ra, u.dec, u.cz) nb = nb.merge(u[["uid", "cz", "W"]], on="uid") nb = nb[nb.otype.isin(GALAXY_OTYPES) & ((nb.rvz_radvel - nb.cz).abs() < nb.W / 2 + 200.0) & (nb.sep_arcmin >= SIMBAD_RMIN)] u["n_simbad_nb"] = u.uid.map(nb[nb.sep_arcmin < SIMBAD_RMAX].groupby("uid").size()).fillna(0).astype(int) u["n_simbad_nb10"] = u.uid.map(nb.groupby("uid").size()).fillna(0).astype(int) u["simbad_nb_sep"] = u.uid.map(nb.groupby("uid").sep_arcmin.min()) first = nb.sort_values("sep_arcmin").drop_duplicates("uid").set_index("uid") u["simbad_nb_name"] = u.uid.map(first.main_id) # annulus statistics: is there an excess of neighbours around under-luminous galaxies? rows = [] for lo, hi in ((0.4, 2), (2, 3.5), (3.5, 5), (5, 7), (7, 10)): has = u.uid.isin(nb.uid[(nb.sep_arcmin >= lo) & (nb.sep_arcmin < hi)]) for lab, sel in (("CF4 with a.100 detection", (u.src == "cf4") & (u.AGC > 0)), ("CF4 without a.100 detection", (u.src == "cf4") & ~(u.AGC > 0)), ("ALFALFA-only", u.src == "alfalfa")): o = sel & (u.D < np.log10(0.25)) rows.append(dict(subsample=lab, annulus_arcmin=f"{lo}-{hi}", N_all=int(sel.sum()), frac_all=has[sel].mean(), N_outliers=int(o.sum()), n_outliers_with_nb=int(has[o].sum()), frac_outliers=has[o].mean())) pd.DataFrame(rows).round(4).to_csv(os.path.join(DATA, "simbad_neighbour_excess.csv"), index=False) # width-error criterion: CF4 only accepts widths with errors <= 20 km/s; we apply the same to all u["auto_clean"] = ((u.n_spec_nb == 0) & (u.n_hi_nb == 0) & ~(u.hi_opt_off > 40.0) & (u.e_W <= E_W_MAX) & (u.n_simbad_nb == 0)) # ---- range of apparent magnitude and mean surface brightness actually populated by each catalogue. # A uniformly dimmed disc is fainter in BOTH by -2.5 log10(1-f). It counts as "still selectable" only if # every band it was measured in stays brighter than the catalogue's faint-end percentile in magnitude and # in surface brightness. (CF4: half-light surface brightness from the CF4 photometry; ALFALFA-SDSS: mean # surface brightness inside the SDSS Petrosian ellipse. The two definitions differ, so limits are per # catalogue.) This is a stand-in for the true selection functions, not a model of them. lim_rows = [] has = {"m_i": u.m_i.notna().values, "m_W1": u.m_W1.notna().values} for pct in PCTS: lim = {} for src in ("cf4", "alfalfa"): for col in ("m_i", "m_W1", "mu_i", "mu_W1"): v = u.loc[u.src == src, col].dropna() lim[(src, col)] = np.percentile(v, pct) if len(v) else np.nan lim_rows.append(dict(percentile=pct, catalogue=src, quantity=col, N=len(v), faint_limit=lim[(src, col)], median=v.median() if len(v) else np.nan)) if pct == PCTS[0]: print(" faint limits (98th percentile):", {f"{k[0]}:{k[1]}": round(float(v), 2) for k, v in lim.items()}) tag = "" if pct == PCTS[0] else f"_p{pct:g}" for f in FS: dm = -2.5 * np.log10(1 - f) magok = np.zeros(len(u), bool) sbok = np.zeros(len(u), bool) sbknown = np.zeros(len(u), bool) for src in ("cf4", "alfalfa"): sel = (u.src == src).values mi = (u.m_i + dm <= lim[(src, "m_i")]).values | ~has["m_i"] mw = (u.m_W1 + dm <= lim[(src, "m_W1")]).values | ~has["m_W1"] magok |= sel & mi & mw # surface brightness: every measured band with a known size must stay within range ki = has["m_i"] & u.mu_i.notna().values kw = has["m_W1"] & u.mu_W1.notna().values si = (u.mu_i + dm <= lim[(src, "mu_i")]).values | ~ki sw = (u.mu_W1 + dm <= lim[(src, "mu_W1")]).values | ~kw sbok |= sel & si & sw sbknown |= sel & (ki | kw) w = np.zeros(len(u)) for src in ("cf4", "alfalfa"): sel = (u.src == src).values ref = sel & sbknown & magok & (u.D < 0).values # comparable galaxies with a size p_unknown = sbok[ref].mean() if ref.any() else 1.0 w[sel] = np.where(sbknown[sel], (magok & sbok)[sel].astype(float), magok[sel] * p_unknown) if pct == PCTS[0]: print(f" f={f}: {src:7s} mag-range pass {magok[sel & (u.D < 0).values].mean():.3f}; " f"surface-brightness pass given mag {p_unknown:.3f}; size unknown for " f"{(~sbknown[sel]).mean():.3f} of galaxies") u[f"w_sel_{f}{tag}"] = w # baseline selection weight (magnitude AND surface brightness) if pct == PCTS[0]: u[f"inlim_{f}"] = magok # magnitude range only (first-pass definition; kept as a variant) u[f"sbknown_{f}"] = sbknown # baryonic residual after removing a fraction f of the starlight (gas unchanged) mb = np.log10((1 - f) * 10 ** u.logMs + 1.33 * 10 ** u.logMHI) u[f"rbar_dim_{f}"] = u.r_bar + (mb - u.logMbar) pd.DataFrame(lim_rows).round(3).to_csv(os.path.join(DATA, "selection_limits.csv"), index=False) u.to_csv(os.path.join(DATA, "union_sample.csv"), index=False) print(f" union sample: {len(u)} unique galaxies ({(u.src == 'cf4').sum()} CF4, {(u.src == 'alfalfa').sum()} ALFALFA-only);" f" auto-clean {u.auto_clean.sum()}") print(f" flag rates: spec neighbour {(u.n_spec_nb > 0).mean():.3f}, HI neighbour {(u.n_hi_nb > 0).mean():.3f}, " f"HI-optical offset>40\" {(u.hi_opt_off > 40).mean():.3f}, e_W>20 {(u.e_W > E_W_MAX).mean():.3f}, " f"SIMBAD neighbour<5' {(u.n_simbad_nb > 0).mean():.3f}") both = u1.D_other.notna() print(f" galaxies measured in both catalogues: {both.sum()}; D(CF4)-D(ALFALFA): median " f"{(u1.D - u1.D_other)[both].median():+.3f}, robust sigma " f"{1.4826 * (u1.D - u1.D_other)[both].sub((u1.D - u1.D_other)[both].median()).abs().median():.3f} dex") for t, f in zip((np.log10(0.5), np.log10(0.25), -1.0), FS): c = u[u.D < t] print(f" D<{t:+.3f}: {len(c)} ({(c.src == 'cf4').sum()} CF4 / {(c.src == 'alfalfa').sum()} ALFALFA-only); " f"auto-clean {c.auto_clean.sum()}; and baryon-deficient {(c.auto_clean & (c.r_bar < BTF_CUT)).sum()}") if __name__ == "__main__": main()