"""Step 5: expected number of chance-alignment false positives from blank-sky control positions. For every GCNS star we queried AllWISE within 6 arcsec of a control position 90 arcsec north. A background source (galaxy, AGN, unrelated star) that fell inside a target's W3/W4 beam would add its flux to the target's profile-fit photometry. For every searchable star i and every control source j near *its own* control position (same sky region, same depth) we ask whether j, superposed on star i, would have produced a detection, using star i's actual noise and photosphere: excess vs Ks anchor = S_b,j excess vs W1 anchor = S_b,j - S_W1,j * (P_b / F_W1)_i (the source also raises W1) A source at offset r contributes S * exp(-r^2 / 4 sigma_psf^2) to a profile fit (Gaussian PSF), so it triggers a detection inside an area A = 4 pi sigma_psf^2 ln(S_eff / threshold), which we cap at the area of a circle of radius R_MAX (beyond that AllWISE detects and deblends it, and our nb = 1 / no-neighbour cuts remove the star). Expected false positives = sum_i sum_j A_ij / A_control. Two estimators are given, because the cap is an assumption: (a) analytic trigger area (above), caps 6 arcsec (W3) and 8 arcsec (W4); variant 6 / 6; (b) direct count: every control source is kept at its ACTUAL offset r from the control position, its flux is attenuated by exp(-r^2 / 4 sigma_psf^2), and it is counted if it would still trigger a detection and r <= cap, for caps of 3, 4 and 6 arcsec (the control cones are 6 arcsec, so larger caps cannot be counted directly). The spread between these numbers is the systematic uncertainty of the expectation. Neighbour census (empirical support for the cap): if AllWISE deblended a W3-bright, galaxy-coloured source lying 2-6 arcsec from a star, such sources would be catalogued next to the targets at the same surface density as around the control positions. The ratio of the two densities in rings of offset measures the fraction that is NOT merged into the star. Output: data/control_expectation.csv, data/control_sources.csv, data/control_direct_count.csv, data/neighbour_census.csv """ import glob import os import numpy as np import pandas as pd from common import RAW, DATA, F0, ZP_DN import importlib ex = importlib.import_module("03_excess") inj = importlib.import_module("04_injection") CTRL_DIR = os.path.join(RAW, "allwise_control_chunks") R_CTRL = 6.0 # arcsec, control search radius FWHM = {"3": 6.5, "4": 12.0} # arcsec (WISE W3, W4 beam FWHM) R_MAX = {"3": 6.0, "4": 8.0} # arcsec, radius beyond which the contaminant is separately catalogued R_MAX_ALT = {"3": 6.0, "4": 6.0} # robustness variant DIRECT_CAPS = (3.0, 4.0, 6.0) # arcsec, caps of the direct-count estimator CENSUS_RINGS = ((0.0, 2.0), (2.0, 4.0), (4.0, 6.0)) # arcsec CENSUS_W13 = 2.0 # W1-W3 (Vega) redder than this = galaxy-coloured (stars: 0 to 0.6) CENSUS_SNR3 = 5.0 DERIVED = os.path.join(RAW, "derived") def load_control(): files = sorted(glob.glob(os.path.join(CTRL_DIR, "aw_*.csv"))) c = pd.concat([pd.read_csv(f, dtype={"designation": str, "cc_flags": str, "ph_qual": str, "var_flg": str, "moon_lev": str}, low_memory=False) for f in files], ignore_index=True) c["cc_flags"] = c.cc_flags.fillna("0000").str.zfill(4) for b in ("W1", "W3", "W4"): k = 10 ** (-0.4 * ZP_DN[b]) * F0[b] c["S" + b] = np.clip(c[f"{b.lower()}flux"].fillna(0) * k, 0, None) # actual offset of the source from the control position (arcsec) cosd = np.cos(np.radians(c.dec_in.values)) c["r_ctrl"] = np.hypot((c.ra.values - c.ra_in.values) * cosd, c.dec.values - c.dec_in.values) * 3600.0 return c def neighbour_census(d, c): """Surface density of separately catalogued, W3-detected, galaxy-coloured AllWISE sources in rings around the M-dwarf targets relative to the same rings around their control positions.""" path = os.path.join(DERIVED, "neighbours.csv.gz") if not os.path.exists(path): return None nb = pd.read_csv(path) rows = [] for cls in ("FGK", "M"): tgt = d[(d.cls == cls) & d.has_aw] n_t, n_c = len(tgt), int((d.cls == cls).sum()) a = nb[nb.sid.isin(tgt.sid) & (nb.w3snr >= CENSUS_SNR3) & (nb.w1mpro - nb.w3mpro > CENSUS_W13)] cc = c[c.sid.isin(d.sid[d.cls == cls]) & (c.w3snr >= CENSUS_SNR3) & (c.w1mpro - c.w3mpro > CENSUS_W13)] for lo, hi in CENSUS_RINGS: k_t = int(((a.sep >= lo) & (a.sep < hi)).sum()) k_c = int(((cc.r_ctrl >= lo) & (cc.r_ctrl < hi)).sum()) exp_t = k_c * n_t / n_c # expected around the targets if nothing were merged rows.append(dict(cls=cls, ring_arcsec=f"{lo:g}-{hi:g}", n_targets=n_t, n_control_positions=n_c, n_sources_near_targets=k_t, n_sources_near_controls=k_c, expected_near_targets_if_deblended=exp_t, ratio=k_t / exp_t if exp_t > 0 else np.nan, ratio_err=np.sqrt(k_t + 1.0) / exp_t if exp_t > 0 else np.nan)) out = pd.DataFrame(rows) out.to_csv(os.path.join(DATA, "neighbour_census.csv"), index=False, float_format="%.4g") return out def main(thr=ex.THR, verbose=True, r_max=None, d=None, c=None): r_max = r_max or R_MAX d = inj.load_stars() if d is None else d c = load_control() if c is None else c n_ctrl_pos = len(d) # one control position per parent star m = c.merge(d[["sid", "cls", "q3", "q4", "FW1", "P3W", "P4W", "sig3K", "sig3W", "sig4K", "sig4W", "chi3K", "chi3W", "chi4K", "chi4W", "dist_pc", "w1mpro"]], on="sid", how="inner") rows, keep, drows, dkeep = [], [], [], [] a_ctrl = np.pi * R_CTRL ** 2 for b, i in (("3", 2), ("4", 3)): s = m[m["q" + b] & (m.cc_flags.str[i] == "0") & (m.ext_flg <= 1)].copy() S = s["SW" + b].values # significance the contaminant alone would add, for each anchor dK = S / s[f"sig{b}K"].values dW = (S - s.SW1.values * s[f"P{b}W"].values / s.FW1.values) / s[f"sig{b}W"].values # add to the star's own (noise) significance and take the minimum, as in the pipeline tot = np.minimum(s[f"chi{b}K"].values + dK, s[f"chi{b}W"].values + dW) base = np.minimum(s[f"chi{b}K"].values, s[f"chi{b}W"].values) ratio = np.minimum(dK, dW) / np.clip(thr - base, 0.5, None) # S_eff / (flux needed to reach threshold) sig_psf = FWHM[b] / 2.3548 area = 4 * np.pi * sig_psf ** 2 * np.log(np.clip(ratio, 1.0, None)) area = np.minimum(area, np.pi * r_max[b] ** 2) s["w_fp"] = np.where((ratio > 1) & (base < thr), area / a_ctrl, 0.0) s["band"] = "W" + b # direct count: the source at its actual offset, attenuated by the Gaussian beam s["trig_direct"] = (ratio * np.exp(-s.r_ctrl.values ** 2 / (4 * sig_psf ** 2)) > 1) & (base < thr) dkeep.append(s[s.trig_direct][["sid", "cls", "band", "designation", "r_ctrl"]]) for cap in DIRECT_CAPS: for cls in ("A", "FGK", "M", "all"): sel = ((s.cls == cls) if cls != "all" else np.ones(len(s), bool)) & s.trig_direct & (s.r_ctrl <= cap) drows.append(dict(band="W" + b, cls=cls, cap_arcsec=cap, n_direct=int(sel.sum()))) s["W1mW3"] = -2.5 * np.log10((s.SW1 / F0["W1"]) / np.clip(s.SW3 / F0["W3"], 1e-12, None)) keep.append(s[s.w_fp > 0][["sid", "cls", "band", "designation", "ra", "dec", "w1mpro_x", "w3mpro", "w4mpro", "SW1", "SW3", "SW4", "w_fp", "W1mW3", "ext_flg"]]) for cls in ("A", "FGK", "M", "all"): sel = (s.cls == cls) if cls != "all" else np.ones(len(s), bool) nstar = int((d["q" + b] & ((d.cls == cls) if cls != "all" else True)).sum()) rows.append(dict(band="W" + b, cls=cls, n_searchable=nstar, n_ctrl_sources_capable=int((s.w_fp[sel] > 0).sum()), expected_false_positives=float(s.w_fp[sel].sum()), expected_fp_if_full_6arcsec=float((s.w_fp[sel] > 0).sum()))) # either band: the two trigger areas are concentric, so the union is the larger one allk = pd.concat(keep) uni = allk.groupby(["sid", "designation", "cls"], as_index=False).w_fp.max() for cls in ("A", "FGK", "M", "all"): sel = (uni.cls == cls) if cls != "all" else np.ones(len(uni), bool) nstar = int(((d.q3 | d.q4) & ((d.cls == cls) if cls != "all" else True)).sum()) rows.append(dict(band="W3 or W4", cls=cls, n_searchable=nstar, n_ctrl_sources_capable=int(sel.sum()), expected_false_positives=float(uni.w_fp[sel].sum()), expected_fp_if_full_6arcsec=float(sel.sum()))) dall = pd.concat(dkeep).groupby(["sid", "designation", "cls"], as_index=False).r_ctrl.min() for cap in DIRECT_CAPS: for cls in ("A", "FGK", "M", "all"): sel = ((dall.cls == cls) if cls != "all" else np.ones(len(dall), bool)) & (dall.r_ctrl <= cap) drows.append(dict(band="W3 or W4", cls=cls, cap_arcsec=cap, n_direct=int(sel.sum()))) direct = pd.DataFrame(drows) direct["n_direct_err"] = np.sqrt(direct.n_direct.clip(lower=1)) out = pd.DataFrame(rows) out["rate_per_1e4_stars"] = 1e4 * out.expected_false_positives / out.n_searchable # Poisson uncertainty from the number of capable control sources out["expected_fp_err"] = out.expected_false_positives / np.sqrt(out.n_ctrl_sources_capable.clip(lower=1)) if verbose: out.to_csv(os.path.join(DATA, "control_expectation.csv"), index=False, float_format="%.4g") pd.concat(keep).to_csv(os.path.join(DATA, "control_sources.csv"), index=False, float_format="%.6g") print(out.to_string()) ks = pd.concat(keep) print("control sources capable of faking an excess: W1-W3 colour quantiles", ks[ks.band == "W3"].W1mW3.quantile([0.1, 0.5, 0.9]).round(2).to_dict()) direct.to_csv(os.path.join(DATA, "control_direct_count.csv"), index=False, float_format="%.4g") print(direct.pivot_table(index=["band", "cap_arcsec"], columns="cls", values="n_direct").to_string()) cen = neighbour_census(d, c) if cen is not None: print(cen.to_string()) return out if __name__ == "__main__": main()