"""Step 8: account for every excess star -> data/candidates.csv. Evidence (the catalogue queries of step 6, the image centroids of step 7, and the manual image-inspection table data/visual_vetting.csv) is turned into one explanation class per star: young_star : independent evidence of youth: SIMBAD young-star type, or member of a recognised young association / moving group catalogued_excess_only : no independent youth evidence, but listed in a published infrared- excess / debris-disc catalogue. This is an identification, NOT an explanation: a swarm host would earn the same catalogue entry. The limits of step 9 are therefore given both with these stars accepted as natural discs and with them counted as unexplained. other_known_type : cataclysmic variable, white dwarf / hot subdwarf, evolved or extragalactic object according to SIMBAD background_catalogued : a catalogued galaxy / AGN / radio source / 2MASS extended source / extended optical source inside the WISE beam at the WISE epoch background_image_offset : the W3/W4 light centroid is displaced from the star by more than 99% of non-excess comparison stars of the same S/N (step 7) background_visual : (manual) contaminant, nebulosity or blend seen in the image cutouts comoving_companion : a co-moving Gaia source inside the beam bright enough to matter unexplained : none of the above Secondary indicators are recorded but do NOT change the class (so the 'unexplained' count is conservative): W1-W2 and J-Ks colours, the GCNS membership probability, Gaia photometric blending flags, a Gaia source within 4 arcsec that is not shown to be co-moving, and whether SIMBAD attaches the star to Sgro & Song (2021), who searched Gaia DR2 K/M dwarfs within 100 pc for W3/W4 excess (their list is a candidate list, not an explanation). """ import os import re import numpy as np import pandas as pd from common import RAW, DATA, C_LIGHT, LAM import importlib inj = importlib.import_module("04_injection") VET = os.path.join(RAW, "vet") VET_CMP = os.path.join(RAW, "vet_cmp") YOUNG_OTYPES = {"TT*", "TT?", "Y*O", "Y*?", "Or*", "pr*", "pr?", "Ae*", "Ae?", "HH"} OTHER_OTYPES = {"CV*", "CV?", "No*", "Sy*", "HS*", "Mi*", "LP*", "pA*", "C*", "S*", "AB*", "Be*"} # recognised young (< ~150 Myr) associations / moving groups / clusters (SIMBAD parent names) YOUNG_RE = re.compile(r"TW Hya|beta Pic|Mamajek|eta Cha|eps Cha|Columba|Carina(?!-near)|THA|Tuc|Argus|AB Dor|" r"Octans|32 Ori|Sco|Centaurus|Lupus|Fornax-Horologium|Alessi +13|SigMA|Pleiades|IC 2391|" r"Psc-Eri|MELANGE|Volans|Taurus|Cha |Oph|CrA|Corona", re.I) BG_SIMBAD = {"G", "GiG", "GiP", "GiC", "BiC", "AGN", "QSO", "Sy1", "Sy2", "SyG", "LIN", "SBG", "rG", "H2G", "EmG", "BLL", "Bla", "LSB", "IG", "PaG", "GrG", "ClG", "CGG", "bCG", "Rad", "mm", "smm", "cm", "G?", "AG?", "Q?", "MoC", "Cld", "DNe", "RNe", "HII", "ISM", "PoG", "FIR", "MIR", "IR", "PN", "GNe", "SFR", "Ir*"} G_FAINT = 14.0 # optical-neighbour evidence is used only for targets fainter than this # (survey catalogues are full of artefacts around brighter stars) R_BEAM = 6.0 # arcsec: inside the W3 beam R_OPT_EXT = 5.0 # extended optical source SELF_TRACK = 1.5 # arcsec: closer than this to the star's proper-motion track = the star itself R_GAIA_NEIGH = 4.0 # arcsec: a Gaia source this close that is not co-moving is recorded IPD_CLEAN = 2.0 # Gaia ipd_frac_multi_peak (%) above this hints at a resolved double NATURAL = ("young_star", "catalogued_excess_only", "other_known_type") def phot_clean(d): """Gaia photometry not obviously affected by blending: single-peaked image and a BP/RP flux excess below the usual colour-dependent limit 1.3 + 0.06 (BP-RP)^2.""" bprp = (d.BPmag - d.RPmag).values lim = 1.3 + 0.06 * np.where(np.isfinite(bprp), bprp, 0.0) ** 2 return (d.IPDfmp.fillna(0).values <= IPD_CLEAN) & (d.bprp_excess.fillna(9).values < lim) def track_distance(ra, dec, star): """Minimum angular distance (arcsec) of (ra, dec) from the star's track 1997-2021.""" ep = np.arange(1997.0, 2021.01, 0.25) - 2016.0 cosd = np.cos(np.radians(star.DE_ICRS)) sra = star.RA_ICRS + (star.pmRA if np.isfinite(star.pmRA) else 0) * ep / 3.6e6 / cosd sde = star.DE_ICRS + (star.pmDE if np.isfinite(star.pmDE) else 0) * ep / 3.6e6 dx = (np.asarray(ra)[:, None] - sra[None, :]) * cosd * 3600.0 dy = (np.asarray(dec)[:, None] - sde[None, :]) * 3600.0 return np.sqrt(dx ** 2 + dy ** 2).min(axis=1) def rd(name): p = os.path.join(VET, f"{name}.csv") if not os.path.exists(p) or os.path.getsize(p) < 5: return pd.DataFrame({"sid": []}) return pd.read_csv(p, low_memory=False) def evidence(c, vet_dir=VET): """Return a DataFrame of evidence columns indexed like c.""" global VET VET = vet_dir sb, so, sp = rd("simbad_basic"), rd("simbad_otypes"), rd("simbad_parents") xm = {n: rd("xm_" + n) for n in ("simbad", "gaia", "ps1", "sdss", "des", "catwise", "xsc", "nvss", "first", "vlass", "racs", "sumss", "cotten_prime", "cotten_reserved", "patel14", "chen14", "mcdonald17", "theissen14", "theissen17", "gagne18_members", "gagne18b_candidates", "gagne18b_bonafide")} sg = rd("simbad_sgro2021") sg_ids = set(sg.sid) if len(sg) else set() rows = [] for _, s in c.iterrows(): sid = s.sid e = dict(sid=sid) e["in_sgro_song_2021"] = sid in sg_ids b = sb[sb.sid == sid] e["simbad_id"] = b.main_id.iloc[0] if len(b) else "" e["simbad_otype"] = b.otype.iloc[0] if len(b) else "" e["sp_type"] = (b.sp_type.iloc[0] if len(b) and isinstance(b.sp_type.iloc[0], str) else "") ot = set(so[so.sid == sid].otype) if len(so) else set() e["otypes"] = "|".join(sorted(ot)) par = sp[(sp.sid == sid)] if len(sp) else sp pnames = sorted(set(par.parent)) if len(par) else [] young_par = [p for p in pnames if YOUNG_RE.search(p)] g1 = xm["gagne18_members"] gm = sorted(set(g1[g1.sid == sid].Assoc.astype(str))) if len(g1) else [] g2 = xm["gagne18b_candidates"] gc = sorted(set(g2[g2.sid == sid].Assoc.astype(str))) if len(g2) and "Assoc" in g2 else [] e["associations"] = "; ".join(pnames + ["BANYAN:" + x for x in gm + gc]) e["ev_young_type"] = bool(ot & YOUNG_OTYPES) e["ev_young_assoc"] = bool(young_par or gm or gc or len(xm["gagne18b_bonafide"].query("sid == @sid"))) cats = [n for n in ("cotten_prime", "cotten_reserved", "patel14", "chen14", "mcdonald17", "theissen14", "theissen17") if len(xm[n]) and (xm[n].sid == sid).any()] e["excess_catalogues"] = ";".join(cats) e["ev_known_excess"] = bool(cats) e["ev_other_type"] = bool(ot & OTHER_OTYPES) e["other_type"] = "|".join(sorted(ot & OTHER_OTYPES)) # ---- background evidence around the WISE-epoch position notes = [] t = xm["simbad"] t = t[t.sid == sid] if len(t) else t ev_bg_simbad = False if len(t): td = track_distance(t["ra.1"].values, t["dec.1"].values, s) for (_, r), dtrk in zip(t.iterrows(), td): if dtrk > SELF_TRACK and r.angDist <= R_BEAM and str(r.otype) in BG_SIMBAD: ev_bg_simbad = True notes.append(f"SIMBAD {r.main_id} ({r.otype}) at {r.angDist:.1f}\"") e["ev_bg_simbad"] = ev_bg_simbad ev_radio = False for n, rmax, fcol in (("first", 6.0, "Fint"), ("vlass", 6.0, "Ftot"), ("racs", 8.0, "Ftot"), ("nvss", 12.0, "S1.4"), ("sumss", 12.0, "St")): t = xm[n] t = t[(t.sid == sid) & (t.angDist <= rmax)] if len(t) else t if len(t): ev_radio = True notes.append(f"{n.upper()} {t[fcol].iloc[0]:.1f} mJy at {t.angDist.iloc[0]:.1f}\"") e["ev_bg_radio"] = ev_radio t = xm["xsc"] t = t[(t.sid == sid)] if len(t) else t ev_xsc = False if len(t): ok = t.angDist <= R_BEAM + t.r3sig.fillna(0) if ok.any(): ev_xsc = True notes.append(f"2MASX {t[ok]['2MASX'].iloc[0]} at {t[ok].angDist.iloc[0]:.1f}\"") e["ev_bg_xsc"] = ev_xsc # Gaia sources near this target that are stars (significant parallax or proper motion): # optical "galaxies" coinciding with them are stars (often the blended companion) gt = xm["gaia"] gt = gt[(gt.sid == sid) & (gt.Source != sid)] if len(gt) else gt if len(gt): pmsn = np.hypot(gt.pmRA / gt.e_pmRA, gt.pmDE / gt.e_pmDE) gstars = gt[(gt.RPlx > 3) | (pmsn > 4)] else: gstars = gt def not_star(ra, dec): if not len(gstars): return np.ones(len(ra), bool) cosd = np.cos(np.radians(s.DE_ICRS)) dx = (np.asarray(ra)[:, None] - gstars.RAdeg.values[None, :]) * cosd * 3600 dy = (np.asarray(dec)[:, None] - gstars.DEdeg.values[None, :]) * 3600 return np.sqrt(dx ** 2 + dy ** 2).min(axis=1) > 1.2 faint = s.Gmag >= G_FAINT ext, anyopt = [], [] t = xm["sdss"] t = t[(t.sid == sid)] if len(t) else t if len(t): td = track_distance(t.RA_ICRS.values, t.DE_ICRS.values, s) t = t.assign(trk=td) t = t[(t.trk > SELF_TRACK) & (t["mode"] == 1)] t = t[not_star(t.RA_ICRS.values, t.DE_ICRS.values)] if len(t) else t g = t[(t["class"] == 3) & (t.angDist <= R_OPT_EXT) & (t.rmag < 22.0)] if len(g) and faint: ext.append(f"SDSS galaxy r={g.rmag.min():.1f} at {g.angDist.min():.1f}\"") if (t.angDist <= 4.0).any(): anyopt.append("SDSS") t = xm["des"] t = t[(t.sid == sid)] if len(t) else t if len(t): td = track_distance(t.RAJ2000.values, t.DEJ2000.values, s) t = t.assign(trk=td) t = t[t.trk > SELF_TRACK] t = t[not_star(t.RAJ2000.values, t.DEJ2000.values)] if len(t) else t g = t[(t.ExtClsCoad >= 2) & (t.angDist <= R_OPT_EXT) & (t.rmag < 23.0)] if len(g) and faint: ext.append(f"DES extended r={g.rmag.min():.1f} at {g.angDist.min():.1f}\"") if (t.angDist <= 4.0).any(): anyopt.append("DES") t = xm["ps1"] t = t[(t.sid == sid)] if len(t) else t if len(t): td = track_distance(t.RAJ2000.values, t.DEJ2000.values, s) t = t.assign(trk=td) t = t[(t.trk > SELF_TRACK) & (t.Nd >= 5)] t = t[not_star(t.RAJ2000.values, t.DEJ2000.values)] if len(t) else t de = np.fmax((t.rmag - t.rKmag).fillna(-9), (t.imag - t.iKmag).fillna(-9)) g = t[(de > 0.05) & (t.angDist <= R_OPT_EXT) & (t.rmag.fillna(t.imag) < 21.5) & (t.rmag.fillna(t.imag) > 15.0)] if len(g) and faint: ext.append(f"PS1 extended r={g.rmag.min():.1f} at {g.angDist.min():.1f}\"") if (t.angDist <= 4.0).any(): anyopt.append("PS1") e["ev_bg_optical_ext"] = bool(ext) notes += ext # ---- Gaia neighbours: co-moving companions and other sources in the beam t = xm["gaia"] t = t[(t.sid == sid) & (t.Source != sid)] if len(t) else t comp, bgstar, ncsep = [], [], [] if len(t): for _, r in t.iterrows(): comoving = (np.isfinite(r.Plx) and abs(r.Plx - s.Plx) < max(0.2 * s.Plx, 3 * (r.e_Plx if np.isfinite(r.e_Plx) else 1))) if comoving and np.isfinite(r.pmRA): pm_s = np.hypot(s.pmRA, s.pmDE) comoving = np.hypot(r.pmRA - s.pmRA, r.pmDE - s.pmDE) < max(0.3 * pm_s, 5.0) if comoving: comp.append(f"Gaia companion G={r.Gmag:.1f} at {r.angDist:.1f}\"") elif r.angDist <= R_GAIA_NEIGH: bgstar.append(f"Gaia source G={r.Gmag:.1f} at {r.angDist:.1f}\"") ncsep.append(float(r.angDist)) anyopt.append("Gaia") e["ev_companion"] = bool(comp) e["ev_gaia_neighbour_4as"] = bool(bgstar) e["gaia_neighbour_minsep"] = min(ncsep) if ncsep else np.nan e["ev_opt_neighbour_4as"] = bool(anyopt) and bool(faint) notes += comp + bgstar t = xm["catwise"] t = t[(t.sid == sid)] if len(t) else t e["n_catwise_other_8as"] = 0 if len(t): td = track_distance(t.RAdeg.values, t.DEdeg.values, s) oth = t[(td > 2.0) & (t.angDist <= 8.0)] e["n_catwise_other_8as"] = int(len(oth)) if len(oth): notes.append(f"CatWISE source W1={oth.W1mproPM.min():.1f} at {oth.angDist.min():.1f}\"") e["notes"] = "; ".join(notes) rows.append(e) return pd.DataFrame(rows) def classify(ev, visual=None, offset=None): cls = np.full(len(ev), "unexplained", dtype=object) if offset is not None: cls[np.asarray(offset, bool)] = "background_image_offset" comp = ev.ev_companion.values bg = (ev.ev_bg_simbad | ev.ev_bg_radio | ev.ev_bg_xsc | ev.ev_bg_optical_ext).values other = ev.ev_other_type.values young = (ev.ev_young_type | ev.ev_young_assoc).values known = ev.ev_known_excess.values & ~young cls[comp] = "comoving_companion" cls[bg] = "background_catalogued" cls[other] = "other_known_type" cls[known] = "catalogued_excess_only" cls[young] = "young_star" vis = np.full(len(ev), "", dtype=object) if visual is not None and len(visual): v = visual.set_index("sid") for i, sid in enumerate(ev.sid.values): if sid in v.index: vis[i] = str(v.loc[sid, "visual_class"]) if cls[i] == "unexplained" and vis[i] in ("background", "nebulosity", "blend", "artefact"): cls[i] = "background_visual" return cls, vis def base_rates(d, out): """How often do the catalogue background criteria fire on matched NON-excess stars?""" pfile = os.path.join(VET_CMP, "positions.csv") if not os.path.exists(pfile): return None cp = pd.read_csv(pfile) cmp_ = d[d.sid.isin(cp.sid)].reset_index(drop=True) evc = evidence(cmp_, VET_CMP).merge(cmp_[["sid", "cls"]], on="sid") evidence(d.iloc[:0], VET) # restore the module-level directory flags = ["ev_bg_simbad", "ev_bg_radio", "ev_bg_xsc", "ev_bg_optical_ext"] evc["bg_any"] = evc[flags].any(axis=1) out = out.copy() out["bg_any"] = out[flags].any(axis=1) nonyoung = out[~out.explanation.isin(NATURAL)] unexpl = out[out.explanation == "unexplained"] rows = [] for cls in ("A", "FGK", "M"): a, b, c0 = out[out.cls == cls], nonyoung[nonyoung.cls == cls], evc[evc.cls == cls] u = unexpl[unexpl.cls == cls] # ev_gaia_neighbour_4as is NOT part of bg_any and is not used for classification: it is a # secondary indicator whose rate on the unexplained stars is compared with the base rate for name in flags + ["bg_any", "ev_gaia_neighbour_4as"]: p0 = c0[name].mean() if len(c0) else np.nan p1 = b[name].mean() if len(b) else np.nan pu = u[name].mean() if len(u) else np.nan rows.append(dict(cls=cls, criterion=name, n_comparison=len(c0), n_flag_comparison=int(c0[name].sum()), rate_comparison=p0, n_excess_all=len(a), n_flag_excess_all=int(a[name].sum()), n_excess_not_young=len(b), n_flag_excess_not_young=int(b[name].sum()), rate_excess_not_young=p1, min_background_fraction=(p1 - p0) / (1 - p0) if len(b) and p0 < 1 else np.nan, n_unexplained=len(u), n_flag_unexplained=int(u[name].sum()), n_unexplained_in_excess_of_base_rate=(len(u) * (pu - p0) / (1 - p0) if len(u) and p0 < 1 else np.nan))) br = pd.DataFrame(rows) br.to_csv(os.path.join(DATA, "background_base_rates.csv"), index=False, float_format="%.4g") print(br[br.criterion.isin(["bg_any", "ev_gaia_neighbour_4as"])].to_string()) return br def main(): d = inj.load_stars() c = d[d.det].copy().reset_index(drop=True) ev = evidence(c, VET) vfile = os.path.join(DATA, "visual_vetting.csv") visual = pd.read_csv(vfile) if os.path.exists(vfile) else None cfile = os.path.join(DATA, "image_centroids.csv") if os.path.exists(cfile): cen = pd.read_csv(cfile)[["sid", "offW3", "nW3", "offW4", "nW4", "flag_off_W3", "flag_off_W4", "flag_offset", "lost_W3", "lost_W4"]] ev = ev.merge(cen, on="sid", how="left") ev["flag_offset"] = ev.flag_offset.fillna(False).astype(bool) else: ev["flag_offset"] = False cls, vis = classify(ev, visual, ev.flag_offset.values) nu = {"3": C_LIGHT / (LAM["W3"] * 1e-6), "4": C_LIGHT / (LAM["W4"] * 1e-6)} out = c[["sid", "designation", "RA_ICRS", "DE_ICRS", "gl", "gb", "cls", "dist_pc", "Gmag", "M_G", "dM", "g_rp", "RUWE", "Lbol", "Ksmag", "w1mpro", "w3mpro", "w4mpro", "w3snr", "w4snr", "chi3", "chi4", "det3", "det4", "q3", "q4", "ebv", "w3rchi2", "w4rchi2"]].copy() out["E3_mJy"] = 1e3 * c.E3K out["E4_mJy"] = 1e3 * c.E4K out["ratio3"] = c.E3K / c.P3K out["ratio4"] = c.E4K / c.P4K # monochromatic nu F_nu / F_bol of the excess: a rough fractional-luminosity indicator out["nuFnu3_over_Fbol"] = c.E3K * 1e-26 * nu["3"] / c.fbol out["nuFnu4_over_Fbol"] = c.E4K * 1e-26 * nu["4"] / c.fbol # secondary indicators (not used for the class) out["W1mW2"] = (c.w1mpro - c.w2mpro).values out["JmKs"] = (c.Jmag - c.Ksmag).values out["GCNSprob"] = c.GCNSprob.values out["IPDfmp"] = c.IPDfmp.values out["bprp_excess"] = c.bprp_excess.values out["gaia_phot_clean"] = phot_clean(c) out = out.merge(ev, on="sid", how="left") out["visual_class"] = vis out["explanation"] = cls out = out.sort_values(["cls", "explanation", "chi3"], ascending=[True, True, False]) out.to_csv(os.path.join(DATA, "candidates.csv"), index=False, float_format="%.6g") print(pd.crosstab(out.explanation, out.cls, margins=True).to_string()) base_rates(d, out) return out if __name__ == "__main__": main()