#!/usr/bin/env python3 """The natural mid-infrared baseline of the early-type sample, photometry checks, per-galaxy caps and stacks. Outputs (data/): phot_checks.json AllWISE profile-fit vs Legacy Surveys forced photometry; noise validation baseline_by_class.csv W3/W1 and W4/W1 flux ratios by emission-line class and for the 'clean' population feq_percentiles.csv natural emission expressed as an equivalent covering fraction (bright subsample) fmax_percentiles.csv per-galaxy upper caps f_max(T) (largest f not excluded), all galaxies fmax_cdf.csv fraction of galaxies whose cap exceeds f (not corrected for sensitivity) stack_limits.csv population-mean reprocessed fraction from summed fluxes """ import json import os import numpy as np import pandas as pd import mirq_common as M import s05_limits as L from s03_model import Model BRIGHT_F1 = 2.0 # mJy: W1 flux above which per-galaxy ratios are measured with small noise def vega_colour(ratio, b): return 2.5 * np.log10(ratio * M.WISE_F0[1] / M.WISE_F0[b]) def boot_ratio(num, den, nboot=300, seed=1): rng = np.random.default_rng(seed) n = len(num) out = np.empty(nboot) for k in range(nboot): i = rng.integers(0, n, n) out[k] = num[i].sum() / den[i].sum() return float(num.sum() / den.sum()), float(out.std()), float(np.percentile(out, 95)) def phot_checks(d): out = {} P, Q = L.phot(d, "psf"), L.phot(d, "ls") both = P["ok"] & Q["ok"] out["n_psf_ok"], out["n_ls_ok"], out["n_both_ok"] = int(P["ok"].sum()), int(Q["ok"].sum()), int(both.sum()) ext = d.ext_flg.values for name, m in (("point_like_ext0", both & (ext == 0)), ("xsc_ext5", both & (ext == 5)), ("theta_e_gt_5as", both & (d.theta_e.values > 5))): r = {"n": int(m.sum()), "L1_over_F1_median": float(np.median(Q["F1"][m] / P["F1"][m]))} for b in (3, 4): hi = m & (P["F%d" % b] / P["e%d" % b] > 10) & (Q["F%d" % b] / Q["e%d" % b] > 5) r["n_w%d_bright" % b] = int(hi.sum()) if hi.sum() > 10: r["L%d_over_F%d_median" % (b, b)] = float(np.median(Q["F%d" % b][hi] / P["F%d" % b][hi])) r["colour_ratio_w%d_median" % b] = float(np.median((Q["F%d" % b][hi] / Q["F1"][hi]) / (P["F%d" % b][hi] / P["F1"][hi]))) rr = np.log10((Q["F%d" % b][hi] / Q["F1"][hi]) / (P["F%d" % b][hi] / P["F1"][hi])) r["colour_ratio_w%d_scatter_dex" % b] = float(1.4826 * np.median(np.abs(rr - np.median(rr)))) out[name] = r # AllWISE elliptical-aperture (2MASS XSC shape) colours vs profile-fit colours, resolved galaxies for b in (3, 4): with np.errstate(invalid="ignore", divide="ignore"): hi = P["ok"] & (ext == 5) & (d["G%d" % b].values / d["ge%d" % b].values > 10) & (d.G1.values > 0) \ & (P["F%d" % b] / P["e%d" % b] > 10) rr = (d["G%d" % b].values / d.G1.values) / (P["F%d" % b] / P["F1"]) out["gmag_vs_psf_w%d" % b] = dict(n=int(hi.sum()), G1_over_F1_median=float(np.median((d.G1.values / P["F1"])[hi])), colour_ratio_median=float(np.median(rr[hi])), colour_ratio_p16=float(np.percentile(rr[hi], 16)), colour_ratio_p84=float(np.percentile(rr[hi], 84))) # noise validation: negative tail of F/sigma (a zero-flux Gaussian population gives 15.9%, 2.3%, 0.13%). # NOTE: the Legacy errors here already include the inflation factors L.LS_ERR_SCALE; the raw scale that # motivates them is re-measured below on line-less galaxies of the usable Legacy subsample (the adopted # factors are the larger of the two estimates per band, see s05_limits.py). ll = (d.line_class == "lineless").values for b in (3, 4): t = (Q["F%d" % b] / (Q["e%d" % b] / L.LS_ERR_SCALE[b]))[Q["ok"] & ll] q = np.percentile(t, [2.3, 16, 50]) t2 = (P["F%d" % b] / P["e%d" % b])[P["ok"] & ll] q2 = np.percentile(t2, [2.3, 16, 50]) out["empirical_error_scale_w%d" % b] = dict(legacy_from_p2_3=float((q[2] - q[0]) / 2), legacy_from_p16=float(q[2] - q[1]), allwise_from_p2_3=float((q2[2] - q2[0]) / 2), allwise_from_p16=float(q2[2] - q2[1]), legacy_scale_adopted=L.LS_ERR_SCALE[b]) for name, p in (("allwise", P), ("legacy_inflated", Q)): for b in (3, 4): t = (p["F%d" % b] / p["e%d" % b])[p["ok"]] out["%s_w%d_snr" % (name, b)] = dict(median=float(np.median(t)), frac_lt_m1=float((t < -1).mean()), frac_lt_m2=float((t < -2).mean()), frac_lt_m3=float((t < -3).mean()), frac_gt_3=float((t > 3).mean()), frac_gt_5=float((t > 5).mean())) # difference of the two photometries in units of the combined error (same WISE exposures -> correlated noise, # so this is expected to be narrower than 1 for sky-limited sources) for b in (3, 4): dz = ((Q["F%d" % b] - P["F%d" % b]) / np.hypot(Q["e%d" % b], P["e%d" % b]))[both] out["delta_w%d_ls_minus_psf" % b] = dict(median=float(np.median(dz)), mad_sigma=float(1.4826 * np.median(np.abs(dz - np.median(dz)))), frac_abs_gt3=float((np.abs(dz) > 3).mean())) # the usable Legacy subsample is NOT a random subset: fracflux_w3 < 0.5 preferentially keeps W3-bright galaxies with np.errstate(invalid="ignore", divide="ignore"): matched = (d.sep_ls.values <= 1.5) & (d.ls_type.values != "DUP") & np.isfinite(d.fracflux_w3.values) r31 = P["F3"] / P["F1"] out["legacy_subsample_selection"] = dict( n_matched_not_dup=int(matched.sum()), n_usable=int(Q["ok"].sum()), n_usable_and_allwise_ok=int(both.sum()), median_allwise_w3_over_w1_usable=float(np.nanmedian(r31[both])), median_allwise_w3_over_w1_all=float(np.nanmedian(r31[P["ok"]])), frac_allwise_w3_snr_gt5_usable=float(((P["F3"] / P["e3"])[both] > 5).mean()), frac_allwise_w3_snr_gt5_all=float(((P["F3"] / P["e3"])[P["ok"]] > 5).mean())) # how well do the Brown+14 template colours describe the observed profile-fit colours? (matters for the # photosphere-subtracted variant: subtracting the BLUEST template is conservative only if real galaxies # are not bluer than it) m = Model() z = d.zhelio.values fl = L.flags(d) clean = L.clean_mask(fl, True) tc = {} for b in (2, 3, 4): with np.errstate(invalid="ignore", divide="ignore"): r = P["F%d" % b] / P["F1"] e = P["e%d" % b] / P["F1"] cm, cb = m.C_of(b, z), m.C_of(b, z, lower=True) for name, sel in (("bright_all", P["ok"] & (P["F1"] > BRIGHT_F1)), ("bright_clean", P["ok"] & clean & (P["F1"] > BRIGHT_F1)), ("very_bright_clean_lineless", P["ok"] & clean & (P["F1"] > 5.0) & (d.line_class == "lineless").values)): tc["w%d_%s" % (b, name)] = dict( n=int(sel.sum()), median_obs_over_median_template=float(np.median((r / cm)[sel])), p05_obs_over_bluest_template=float(np.percentile((r / cb)[sel], 5)), median_obs_over_bluest_template=float(np.median((r / cb)[sel])), frac_below_bluest_template=float(((r / cb)[sel] < 1).mean()), frac_more_than_2sigma_below_bluest=float((((r - cb) / e)[sel] < -2).mean())) out["photosphere_template_check"] = tc json.dump(out, open(os.path.join(M.DATA, "phot_checks.json"), "w"), indent=1) print(json.dumps(out, indent=1)) return out def baseline(d, fl, m): z = d.zhelio.values clean = L.clean_mask(fl, True) pops = [("all", np.ones(len(d), bool)), ("clean", clean)] for c in ("lineless", "retired", "SF", "composite", "AGN"): pops.append(("line:" + c, (d.line_class == c).values)) rows = [] for kind in ("psf", "ls"): p = L.phot(d, kind) for name, sel in pops: s = sel & p["ok"] br = s & (p["F1"] > BRIGHT_F1) row = dict(phot=kind, pop=name, n=int(s.sum()), n_bright=int(br.sum())) for b in (2, 3, 4): Fb = p["F%d" % b] rs, se, _ = boot_ratio(Fb[s], p["F1"][s]) row["sum_F%d_over_F1" % b], row["sum_F%d_over_F1_err" % b] = rs, se q = np.percentile(Fb[br] / p["F1"][br], [5, 16, 50, 84, 95, 99]) for k, v in zip(("p05", "p16", "p50", "p84", "p95", "p99"), q): row["F%d_F1_%s" % (b, k)] = v row["W1-W%d_median_bright" % b] = vega_colour(q[2], b) row["template_F%d_F1" % b] = float(np.median(m.C_of(b, z[br]))) rows.append(row) out = pd.DataFrame(rows) out.to_csv(os.path.join(M.DATA, "baseline_by_class.csv"), index=False, float_format="%.5g") print(out[["phot", "pop", "n", "n_bright", "F3_F1_p16", "F3_F1_p50", "F3_F1_p84", "F3_F1_p99", "template_F3_F1", "F4_F1_p50", "F4_F1_p84", "template_F4_F1", "sum_F3_over_F1", "sum_F4_over_F1"]] .to_string(index=False, float_format=lambda v: "%.3g" % v)) return out def caps(d, fl, m): """Per-galaxy cap f_max(T): largest f for which the galaxy is 'compatible' (AllWISE profile-fit fluxes).""" z = d.zhelio.values P = L.phot(d, "psf") clean = L.clean_mask(fl, True) rows, cdf_rows, feq_rows = [], [], [] fgrid = np.concatenate([[0.001, 0.002], M.F_FINE, [0.7, 0.9]]) for T in M.T_GRID: R = {b: m.R_of(T, b, z) for b in (1, 2, 3, 4)} def gmax(p): g = np.full(len(d), np.inf) with np.errstate(invalid="ignore", divide="ignore"): for b in (2, 3, 4): X = (p["F%d" % b] + L.K_C * p["e%d" % b]) / ((1 - L.EPS) * p["F1"]) X = np.maximum(X, 0.0) den = R[b] - X * R[1] gb = np.where(den > 0, X / np.where(den > 0, den, 1.0), np.inf) g = np.minimum(g, gb) return np.where(p["ok"], g, np.nan) g = gmax(P) f = g / (1 + g) f = np.where(np.isinf(g), 1.0, f) ok = P["ok"] for name, sel in (("all", ok), ("clean", ok & clean)): q = np.nanpercentile(f[sel], [5, 16, 50, 84, 95, 99, 99.9]) rows.append(dict(T=T, pop=name, n=int(sel.sum()), **{k: v for k, v in zip( ("p05", "p16", "p50", "p84", "p95", "p99", "p99.9"), q)})) for ff in fgrid: cdf_rows.append(dict(T=T, pop=name, f=ff, frac_cap_ge_f=float(np.mean(f[sel] >= ff)), n_cap_ge_f=int(np.sum(f[sel] >= ff)))) # natural emission as an equivalent covering fraction, bright subsample, band of maximum leverage b = 4 if T <= 150 else 3 br = ok & (P["F1"] > BRIGHT_F1) with np.errstate(invalid="ignore", divide="ignore"): g_tot = P["F%d" % b] / (R[b] * P["F1"]) g_exc = (P["F%d" % b] - m.C_of(b, z) * P["F1"]) / (R[b] * P["F1"]) for name, sel in (("all", br), ("clean", br & clean)): for lab, gg in (("total", g_tot), ("excess_over_template", g_exc)): q = np.nanpercentile(gg[sel] / (1 + gg[sel]), [5, 16, 50, 84, 95, 99, 99.9]) feq_rows.append(dict(T=T, band="W%d" % b, pop=name, kind=lab, n=int(sel.sum()), **{k: v for k, v in zip( ("p05", "p16", "p50", "p84", "p95", "p99", "p99.9"), q)})) pd.DataFrame(rows).to_csv(os.path.join(M.DATA, "fmax_percentiles.csv"), index=False, float_format="%.4g") pd.DataFrame(cdf_rows).to_csv(os.path.join(M.DATA, "fmax_cdf.csv"), index=False, float_format="%.5g") feq = pd.DataFrame(feq_rows) feq.to_csv(os.path.join(M.DATA, "feq_percentiles.csv"), index=False, float_format="%.4g") print(pd.DataFrame(rows).to_string(index=False, float_format=lambda v: "%.3g" % v)) print(feq.to_string(index=False, float_format=lambda v: "%.3g" % v)) def stacks(d, fl, m): """Summed fluxes: if a fraction p of galaxies re-radiates f at T, then p * g(f) <= SUM F_b / SUM R_b F1 for each band (natural emission only adds). Quoted: the smaller of the W3 and W4 values.""" z = d.zhelio.values clean = L.clean_mask(fl, True) rows = [] for kind in ("psf", "ls"): p = L.phot(d, kind) for name, sel in (("all", p["ok"]), ("clean", p["ok"] & clean)): for T in M.T_GRID: best = None for b in (3, 4): Rb = m.R_of(T, b, z) for sub in (False, True): num = p["F%d" % b] - (m.C_of(b, z, lower=True) * p["F1"] if sub else 0.0) val, se, up = boot_ratio(num[sel], (Rb * p["F1"])[sel], nboot=200) rows.append(dict(phot=kind, pop=name, T=T, band="W%d" % b, photosphere_subtracted=sub, n=int(sel.sum()), mean_g=val, mean_g_err=se, mean_g_ul95=up)) out = pd.DataFrame(rows) out.to_csv(os.path.join(M.DATA, "stack_limits.csv"), index=False, float_format="%.5g") piv = out[(out.phot == "psf")].pivot_table(index=["pop", "photosphere_subtracted", "T"], columns="band", values="mean_g_ul95") print(piv.to_string(float_format=lambda v: "%.4f" % v)) def kbol_check(d): """Is L_bol/(nu L_nu)(W1) of the Brown+14 ellipticals appropriate for the sample? Compare the observed r-band / W1 flux ratio of the sample galaxies (Legacy Surveys Tractor model fluxes, Galactic-extinction corrected in r) with the same ratio synthesised from the templates at each galaxy's redshift (r band approximated by a top-hat 5600-7100 A, between the DECam and BASS r filters). The optical carries much of L_bol, so a sample that is redder in r-W1 than the templates has a lower L_bol per unit W1.""" import s03_model as S zs = np.linspace(0.0, 0.14, 15) ratios = [] for name in S.PASSIVE_E: lam, flam = S.load_template(name) fnu = flam * (lam * 1e4) ** 2 / 2.99792458e18 row = [] for z in zs: lr = np.linspace(0.56, 0.71, 60) / (1 + z) fr = np.mean(np.interp(lr, lam, fnu)) f1 = S.synth(lambda x: np.interp(x, lam, fnu), 1, z) row.append(fr / f1) ratios.append(row) tmpl = np.median(np.array(ratios), axis=0) # Legacy Surveys Tractor fluxes: r and W1 are measured with the SAME galaxy model, so the colour is a # total-light colour (AllWISE profile-fit W1 misses 20-40% of the light of these resolved galaxies) # (only r and W1 are needed here, so the W3 blending cut of L.phot(d, "ls") is not applied) L1 = d.L1.values with np.errstate(invalid="ignore"): sel = (d.sep_ls.values <= 1.5) & (d.ls_type.values != "DUP") & (L1 > 0) & (d.flux_r.values > 0) fr_jy = d.flux_r.values * M.NMGY_JY / 10 ** (-0.4 * d.extinction_r.values) with np.errstate(invalid="ignore", divide="ignore"): obs = fr_jy * 1e3 / L1 rel = obs[sel] / np.interp(d.zhelio.values[sel], zs, tmpl) out = dict(n=int(sel.sum()), median_obs_over_template=float(np.median(rel)), p16=float(np.percentile(rel, 16)), p84=float(np.percentile(rel, 84)), template_spread_rms=float(np.std(np.array(ratios)[:, 0] / tmpl[0]))) # by redshift for lo, hi in ((0.02, 0.05), (0.05, 0.08), (0.08, 0.125)): mz = sel & (d.zhelio.values >= lo) & (d.zhelio.values < hi) out["median_z%.2f_%.2f" % (lo, hi)] = float(np.median(obs[mz] / np.interp(d.zhelio.values[mz], zs, tmpl))) json.dump(out, open(os.path.join(M.DATA, "kbol_check.json"), "w"), indent=1) print("r/W1 flux ratio, sample over template:", json.dumps(out)) return out def main(): d = L.load_master() fl = L.flags(d) m = Model() phot_checks(d) try: kbol_check(d) except (FileNotFoundError, OSError) as e: # raw templates absent (run_all.py --analysis-only on a cold cache) print("kbol_check skipped, raw templates not in the cache (%s); data/kbol_check.json left as is" % e) baseline(d, fl, m) caps(d, fl, m) stacks(d, fl, m) if __name__ == "__main__": main()