"""Step 3: photosphere calibration, excess significance and quality flags for every star. Outputs RAW/derived/stars.csv.gz : parent sample + bolometric flux, predictions, chi, quality flags data/calibration.json : colour loci, saturation corrections, error floors data/chi_diagnostics.csv : robust width / tail counts of the significance distributions """ import json import os import numpy as np import pandas as pd import photosphere as ph from common import RAW, DATA, mag2jy, jy2mag, MBOL_SUN import importlib bs = importlib.import_module("02_build_sample") DERIVED = os.path.join(RAW, "derived") THR = 5.0 # detection threshold on min(chi_Ks-anchored, chi_W1-anchored) EBV_MAX = 0.15 # SFD E(B-V) limit (cirrus / crowded regions); chosen from the NEGATIVE tail only SEP0, SEP_PM = 0.4, 0.35 # max Gaia(2010.5)-AllWISE separation: SEP0 + SEP_PM * proper motion ["/yr] def bolometric(d): """Bolometric flux (W m^-2) and luminosity from Ks + BC_Ks(M_Ks) (Mamajek dwarf table); stars without Ks use G + BC_G(M_G).""" m = bs.read_mamajek().dropna(subset=["Mbol", "M_Ks", "M_G"]) mk = m.sort_values("M_Ks") mg = m.sort_values("M_G") M_Ks = d.Ksmag - 5 * np.log10(d.dist_pc / 10.0) bcK = np.interp(M_Ks, mk.M_Ks.values, (mk.Mbol - mk.M_Ks).values) bcG = np.interp(d.M_G, mg.M_G.values, (mg.Mbol - mg.M_G).values) mbol = np.where(np.isfinite(d.Ksmag), d.Ksmag + bcK, d.Gmag + bcG) fbol = 10 ** (-0.4 * (mbol + 18.997351)) Mbol = mbol - 5 * np.log10(d.dist_pc / 10.0) return fbol, 10 ** (-0.4 * (Mbol - MBOL_SUN)) def quality(d): """Per-band 'searchable' flags. A star is searchable in a band when its AllWISE measurement in that band is free of flagged artefacts, blending and confusion. Every cut was chosen from a-priori contamination arguments and from the NEGATIVE tail of the significance distribution (never from the positive tail): * unique AllWISE match within SEP0 + SEP_PM*mu of the Gaia position at epoch 2010.5; * 2MASS Ks with a quoted uncertainty and a positive W1 profile-fit flux; * SFD E(B-V) < EBV_MAX (cirrus, molecular clouds, the Galactic plane, the Magellanic Clouds); * W1 contamination flag clear (lower-case flags tolerated for W1 < 8, where bright stars flag themselves), not blended (nb = 1, na = 0), not inside a 2MASS extended source; * band flag clear, profile-fit reduced chi2 < 3 + (0.04 S/N)^2 (the second term allows for PSF systematics of bright point sources), and no other AllWISE source within 12 arcsec whose flux in that band exceeds 3 sigma of the target's noise. """ cc = d.cc_flags.fillna("9999") cc1ok = (cc.str[0] == "0") | (cc.str[0].isin(list("hdpo")) & (d.w1mpro < 8.0)) base = (d.has_aw & np.isfinite(d.e_Ksmag) & np.isfinite(d.FW1) & (d.FW1 > 0) & np.isfinite(d.eW1) & (d.n_share == 1) & cc1ok & (d.nb == 1) & (d.na == 0) & (d.ext_flg <= 1) & (d.sep <= SEP0 + SEP_PM * d.pm_asyr) & (d.ebv < EBV_MAX)) q = {} for b, i in (("3", 2), ("4", 3)): snr = d[f"w{b}snr"].clip(lower=0).fillna(0) q[b] = (base & (cc.str[i] == "0") & np.isfinite(d[f"FW{b}"]) & np.isfinite(d[f"eW{b}"]) & (d[f"w{b}rchi2"] < 3.0 + (0.04 * snr) ** 2) & (d[f"neighFW{b}"].fillna(0) < 3 * d[f"eW{b}"])) return base.values, q["3"].values, q["4"].values def main(): d = pd.read_csv(os.path.join(DERIVED, "parent_sample.csv.gz"), dtype={"cc_flags": str, "moon_lev": str, "ph_qual": str, "designation": str}, low_memory=False) d["has_aw"] = d.has_aw.astype(bool) d["fbol"], d["Lbol"] = bolometric(d) d["FKs"] = mag2jy(d.Ksmag, "Ks") use = d.has_aw & np.isfinite(d.e_Ksmag) & np.isfinite(d.FW1) & (d.FW1 > 0) & np.isfinite(d.eW1) m = d[use].copy() cal = ph.calibrate(m) pred = ph.predict(cal, m.M_G.values, m.g_rp.values, m.FKs.values, m.FW1.values) cal = ph.fit_error_floor(cal, m, pred) sig = ph.significance(cal, m.e_Ksmag.values, m.FW1.values, m.eW1.values, m.FW3.values, m.eW3.values, m.FW4.values, m.eW4.values, pred) new = pd.DataFrame({k: v for k, v in {**pred, **sig}.items()}, index=m.index) d = pd.concat([d, new], axis=1) base, q3, q4 = quality(d) det3 = q3 & (d.chi3.values >= THR) det4 = q4 & (d.chi4.values >= THR) d = pd.concat([d, pd.DataFrame(dict(q_base=base, q3=q3, q4=q4, det3=det3, det4=det4, det=det3 | det4), index=d.index)], axis=1) d.to_csv(os.path.join(DERIVED, "stars.csv.gz"), index=False, float_format="%.9g") with open(os.path.join(DATA, "calibration.json"), "w") as f: json.dump({k: [list(map(float, v[0])), list(map(float, v[1]))] for k, v in cal.items()}, f, indent=1) # ---------------------------------------------------------------- diagnostics rows = [] for band, q in (("3", d.q3), ("4", d.q4)): for cls in ("A", "FGK", "M", "all"): s = q & ((d.cls == cls) if cls != "all" else True) chi = d.loc[s, f"chi{band}"].values chin = d.loc[s, f"chi{band}_neg"].values for anc in ("K", "W"): c = d.loc[s, f"chi{band}{anc}"].values rows.append(dict(band="W" + band, cls=cls, stat="chi_" + anc, n=int(s.sum()), median=np.median(c), robust_sigma=ph.robust_sigma(c), n_pos5=int((c >= 5).sum()), n_neg5=int((c <= -5).sum()), n_pos4=int((c >= 4).sum()), n_neg4=int((c <= -4).sum()), n_pos3=int((c >= 3).sum()), n_neg3=int((c <= -3).sum()))) rows.append(dict(band="W" + band, cls=cls, stat="min/max", n=int(s.sum()), median=np.median(chi), robust_sigma=ph.robust_sigma(chi), n_pos5=int((chi >= 5).sum()), n_neg5=int((chin <= -5).sum()), n_pos4=int((chi >= 4).sum()), n_neg4=int((chin <= -4).sum()), n_pos3=int((chi >= 3).sum()), n_neg3=int((chin <= -3).sum()))) diag = pd.DataFrame(rows) diag.to_csv(os.path.join(DATA, "chi_diagnostics.csv"), index=False, float_format="%.4g") print(diag.to_string()) print("parent stars:", len(d), "| base-searchable:", int(base.sum()), "| W3 searchable:", int(q3.sum()), "| W4 searchable:", int(q4.sum())) print("detections: W3", int(d.det3.sum()), "W4", int(d.det4.sum()), "either", int(d.det.sum())) print(d[d.det].groupby("cls").size()) if __name__ == "__main__": main()