"""Part 2b: a search for total occultations in the public OGLE-II bulge DIA light curves. A. number of monitored stars per field and magnitude bin, and the fraction bluer than the red clump (foreground-disc proxy), from the OGLE-II VI maps B. metrics and cut flow for all 221,801 catalogue light curves; candidate list (baseline search: total flux below 50% of the bright level); epochs of the frames stacked into each reference image, recovered from the light curves C. characterisation of the step-like disappearances that all fall in the 1998/1999 seasonal gap D. looser variant (total flux below 70%: a totally occulted star that supplies only 30-50% of a blend) with objective vetting columns for its survivors E. regression checks of the load path (an injected occultation of an I = 16 star must be recovered through the same quality mask as real data; the 1998 decline of the R CrB star Terz V 2637 must be present in its loaded light curve; flags must not be set preferentially on faint-state points) F. injection-recovery efficiency on real light curves, as an effective exposure T_eff(duration), with the injection made into the raw arrays *before* the quality mask is built, an emulation of the variable finder, and the requirement that the star be present in the reference image G. figures 3 (search summary) and 5 (near misses) Outputs: data/ogle2_fields.csv, data/ogle2_cutflow.csv, data/ogle2_prefilter_metrics.csv, data/ogle2_candidates.csv, data/ogle2_c4_survivors.csv, data/ogle2_gap_steps.csv, data/ogle2_variant_cutflow.csv, data/ogle2_variant_survivors.csv, data/ogle2_efficiency.csv, data/ogle2_reference_epochs.csv, data/ogle2_reference_summary.csv, data/ogle2_search_summary.json, figures/fig3_ogle2_search.png, figures/fig5_ogle2_survivors.png Run with --quick for a faster, noisier efficiency estimate. """ import json import os import sys import warnings import numpy as np import pandas as pd import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from astropy.coordinates import SkyCoord import astropy.units as u from common import DATA, FIGS, RAW import ogle2lib as ol warnings.filterwarnings("ignore", category=RuntimeWarning) MAG_EDGES = [0, 15, 16, 17, 17.5, 18] DUR_GRID = [1, 3, 6, 10, 20, 40, 80, 160, 320, 640, 1280, 2560, 5120, 20000] RNG = np.random.default_rng(20261012) ID_STATS = np.zeros((len(MAG_EDGES) - 1, 2), int) # per magnitude bin: catalogue objects, and those with ID_RAD <= cut VARIANT_LOW = 0.7 BLUE_OFFSET = 0.4 # "bluer than the clump": V-I < (mode of V-I for I<16.5) - 0.4 mag REFW = {} # field -> weight of every epoch in the reference image (filled by search()) REF_ROWS = [] # per-field summary of the reference-epoch reconstruction FLAG_STATS = np.zeros((2, 3), int) # [inside sustained low runs, outside] x [epochs, error codes (flags 1-2), flags 3-11] GETVAR_STATS = np.zeros(2, int) # catalogue objects tested with the variable-finder emulation, and those passing I_REF_LIMITS = {"ref185": 18.5, "ref195": 19.5} # faintest reference-image magnitude at which the star is taken to be in the star list POSITIVE_CONTROL = (43, 2786) # Terz V 2637 (R CrB star): 1998 decline to zero flux # Notes from visual inspection of the survivors of the looser variant (see fig5, bottom rows). These # are descriptions, not cuts: the objective columns written next to them are computed by the code. VISUAL_NOTES = { (3, 1203): "bright semi-regular variable (OGLE BLG-LPV-104932) whose minimum grazes the 0.7 threshold", (28, 405): "bright semi-regular variable (OGLE BLG-LPV-35291) whose minimum grazes the 0.7 threshold", (37, 1526): "bright semi-regular variable (OGLE BLG-LPV-92040) whose minimum grazes the 0.7 threshold", (20, 2312): "V-shaped 8 d dip to 0.5-0.7 on the same nights (HJD 2450931-940) as sc21 #2369 in the adjacent field: common-mode systematic", (21, 2369): "ragged 13 d dip to 0.5-0.7 on the same nights (HJD 2450931-945) as sc20 #2312 in the adjacent field: common-mode systematic", (37, 5379): "step down to 0.46 across the 1998/99 gap, 14 arcsec from a G < 10.5 star (bright-star artefact family)", } def star_counts(): """Stars per field and magnitude bin in the OGLE-II maps, and the fraction of I<18 stars that are more than BLUE_OFFSET bluer than the red clump (cached as a small csv in RAW).""" cache = os.path.join(RAW, "ogle2_map_counts_v2.csv") if os.path.exists(cache): return pd.read_csv(cache) rows = [] for i in range(1, 50): d = pd.read_csv(os.path.join(ol.MAPS, f"bul_sc{i}.map.gz"), sep=r"\s+", header=None, usecols=[6, 7], names=["VI", "I"], dtype=np.float32, engine="c") I, VI = d.I.values, d.VI.values h = np.histogram(I, MAG_EDGES)[0] ok = (VI > -5) & (VI < 9) hh, e = np.histogram(VI[ok & (I < 16.5)], np.arange(0, 6, 0.05)) mode = float(e[np.argmax(hh)] + 0.025) b = ok & (I < 18) rows.append({"field": i, "n_map_total": len(d), **{f"n_I_{MAG_EDGES[k]}_{MAG_EDGES[k + 1]}": int(h[k]) for k in range(len(h))}, "n_I_18_18.5": int(((I >= 18) & (I < 18.5)).sum()), "frac_I18_with_V": float(b.sum() / max((I < 18).sum(), 1)), "VI_clump_mode": mode, "frac_I18_bluer_than_clump": float((b & (VI < mode - BLUE_OFFSET)).sum() / max(b.sum(), 1))}) print(f"[ogle2] map counts field {i}: {len(d)} stars, {int(h.sum())} with I<18", flush=True) out = pd.DataFrame(rows) out.to_csv(cache, index=False) return out def search(low_frac=ol.LOW_FRAC, keep="c1c3", collect_hosts=True): """Apply the cuts to every catalogue light curve. Returns (rows, cutflow, fields, hosts, lcs): lcs holds the normalised light curves of the objects that pass C1-C4.""" rows, cutflow = [], np.zeros(len(ol.CUT_NAMES) + 1, int) fields, hosts, lcs = [], [], {} n_total = 0 for i in range(1, 50): cat, t, F, E, flg, frames, jd = ol.load_field_raw(i) good = ol.quality_mask(F, E, flg) n = len(cat) n_total += n if collect_hosts: w, nrows, resid, resid_flat = ol.reference_weights(t, F, E, good, jd, np.array(cat["ID_RAD"], float)) REFW[i] = w sea = np.digitize(jd, [0, 800, 1160, 1560, 3000]) REF_ROWS.append({"field": i, "n_objects_used": nrows, "n_frames_weight_gt_0.01": int((w > 0.01).sum()), "weight_in_frames_gt_0.01": float(w[w > 0.01].sum()), "weight_1997": float(w[sea == 1].sum()), "weight_1998": float(w[sea == 2].sum()), "weight_1999": float(w[sea == 3].sum()), "rms_stack_mean_flux_over_amplitude": resid, "rms_flat_mean_flux_over_amplitude": resid_flat, "frames": frames, "jd": jd, "w": w}) ref = np.array(cat["REF_FLUX"], float) zp = np.array(cat["MAG_0"], float) idr = np.array(cat["ID_RAD"], float) c = SkyCoord(np.median(cat["RA"]) * u.deg, np.median(cat["DEC"]) * u.deg).galactic tt = t[0] fields.append({"field": i, "n_catalogue": n, "n_epochs": t.shape[1], "t_first": tt.min(), "t_last": tt.max(), "n_epochs_1997_1999": int((tt < ol.T_END_DETECTION).sum()), "l": float(c.l.wrap_at(180 * u.deg).deg), "b": float(c.b.deg)}) # quick vectorised prefilter: at least 3 low points Fm = np.where(good, F, np.nan) with np.errstate(all="ignore"): B = np.nanpercentile(Fm, 85, axis=1) tot = ref + B nlow = np.nansum(((ref[:, None] + Fm) / tot[:, None] < low_frac) & good, axis=1) mad = 1.4826 * np.nanmedian(np.abs(Fm - np.nanmedian(Fm, axis=1)[:, None]), axis=1) / tot imag = zp - 2.5 * np.log10(np.where(tot > 0, tot, np.nan)) if collect_hosts: # host stars for injection: bright state I<18, on a star, quiet, no low points at all hsel = np.where((imag < 18) & (idr <= ol.CUTS["idrad_max"]) & (nlow == 0) & (mad < 0.05) & (good.sum(1) > 60))[0] hosts += [(i, int(k), float(imag[k])) for k in hsel] okc = np.isfinite(imag) & (imag < 18) & (good.sum(1) > 60) bidx = np.digitize(imag[okc], MAG_EDGES) - 1 for b in range(len(MAG_EDGES) - 1): ID_STATS[b, 0] += int((bidx == b).sum()) ID_STATS[b, 1] += int(((bidx == b) & (idr[okc] <= ol.CUTS["idrad_max"])).sum()) # sanity check of the variable-finder emulation: every 10th catalogue object should pass it for k in range(0, n, 10): if good[k].sum() > 60: GETVAR_STATS[0] += 1 GETVAR_STATS[1] += int(ol.getvar_like(t[k], F[k], E[k], good[k])) cutflow[0] += n for k in np.where((nlow >= 3) & (tot > 0))[0]: m = ol.metrics(t[k], F[k], E[k], good[k], ref[k], zp[k], low_frac=low_frac) if m is None: continue fl = ol.cut_flags(m, idr[k]) flv = [fl[name] for name in ol.CUT_NAMES] ok = True for j in range(len(flv)): ok = ok and flv[j] if ok: cutflow[j + 1] += 1 if collect_hosts and flv[0] and flv[2]: # are quality flags set preferentially while a bright star is in a sustained faint state? ins = (t[k] >= m["t1"]) & (t[k] <= m["t2"]) err = (flg[k] & 3) != 0 oth = ((flg[k] & ol.QUALITY_BITS & ~3) != 0) & ~err for row, sel in ((0, ins), (1, ~ins)): FLAG_STATS[row] += [int(sel.sum()), int((sel & err).sum()), int((sel & oth).sum())] if (keep == "c1c3" and flv[0] and flv[2]) or (keep == "c1c5" and all(flv[:5])): rows.append({"field": i, "idx": int(k), "ra": float(cat["RA"][k]), "dec": float(cat["DEC"][k]), "x_tpl": float(cat["X_TPL"][k]), "y_tpl": float(cat["Y_TPL"][k]), "id_rad": float(idr[k]), "vtype": int(cat["VTYPE"][k]), "cat_flag": int(cat["FLAG"][k]), "ref_flux": float(ref[k]), **m, **{f"c{j + 1}": bool(v) for j, v in enumerate(flv)}}) if all(flv[:4]) or (keep == "c1c5" and all(flv[:5])): g = good[k] lcs[(i, int(k))] = (t[k][g], (ref[k] + F[k][g]) / m["tot"], E[k][g] / m["tot"]) print(f"[ogle2] searched field {i}/49 (low state < {low_frac:g}): {n} light curves", flush=True) assert n_total == 221801, n_total names = [("C1 at least 3 low points and bright state I<18" if j == 0 else nm) for j, nm in enumerate(ol.CUT_NAMES)] cf = pd.DataFrame({"cut": ["C0 all catalogue light curves"] + names, "n_remaining": cutflow}) return pd.DataFrame(rows), cf, pd.DataFrame(fields), hosts, lcs def bright_catalogue(): g = pd.read_csv(os.path.join(RAW, "gaia_bright_bulge.csv")) gb = g[g.phot_g_mean_mag < 10.5] return SkyCoord(gb.ra.values * u.deg, gb.dec.values * u.deg) def bright_sep(ra, dec, cb): c = SkyCoord(np.asarray(ra) * u.deg, np.asarray(dec) * u.deg) _, sep, _ = c.match_to_catalog_sky(cb) return sep.arcsec def classify_c4(pre, lcs): """Objective description of the objects that pass C1-C4 and fail C5 (no edge resolved in-season).""" from scipy.stats import spearmanr s = pre[pre.c1 & pre.c2 & pre.c3 & pre.c4].copy() rho = [] for f, k in zip(s.field, s.idx): t, r, e = lcs[(int(f), int(k))] rho.append(float(spearmanr(t, r)[0])) s["spearman_time_flux"] = rho in_gap = (s.t1 - s.tprev) > ol.CUTS["edge_gap"] step = in_gap & (s.pre_med > ol.CUTS["edge_level"]) & (s.r_run_med < 0.35) drift = ~step & (np.abs(s.spearman_time_flux) > 0.9) s["class"] = np.where(step, "step across a seasonal gap", np.where(drift, "monotonic drift through the whole light curve", "smooth variation")) return s def gap_steps(pre, cb): """Sustained low states whose edges fall in seasonal gaps (fail C5): when do they start, and are they associated with bright stars? (No centroid cut here: many are not on a star at all.)""" s = pre[pre.c1 & pre.c3 & pre.c4 & ~pre.c5 & (pre.r_run_med < 0.35)].copy() s["sep_bright_arcsec"] = bright_sep(s.ra.values, s.dec.values, cb) s["season_of_run_start"] = np.digitize(s.t1, [0, 800, 1160, 1560, 2000]) s["step_between_1998_and_1999"] = (s.tprev < 1160) & (s.t1 > 1200) & (s.t1 < 1300) return s def control_bright_sep(cb, n=20000): """Separation to the nearest G<10.5 star for random catalogue variables (control sample).""" from astropy.io import fits ra, dec = [], [] for i in range(1, 50): cat = fits.open(os.path.join(ol.DIA, f"bul_sc{i}_cat4.fts"))[1].data k = RNG.choice(len(cat), min(len(cat), n // 49), replace=False) ra.append(np.array(cat["RA"][k], float)) dec.append(np.array(cat["DEC"][k], float)) return bright_sep(np.concatenate(ra), np.concatenate(dec), cb) def cover_fraction(tk, ti, D, ramp=0.0): """Fraction of the star's light that is blocked at times tk by an occultation with first contact at ti and last contact at ti + D; ramp > 0 gives a linear ingress and egress of that length (a resolved stellar disc), ramp = 0 an instantaneous step (a point source).""" if ramp <= 0: return ((tk > ti) & (tk < ti + D)).astype(float) r = min(ramp, 0.5 * D) return np.clip(np.minimum((tk - ti) / r, (ti + D - tk) / r), 0.0, 1.0) def inject(Fraw, cover, depth): """Subtract depth x cover from every epoch that holds a measurement (error codes stay error codes).""" return np.where(Fraw != ol.ERR_CODE, Fraw - depth * cover, Fraw) def idfrac_at(imag): """Fraction of catalogue objects of this bright-state magnitude that pass the centroid cut C2 (objects fainter than 18 take the value of the faintest bin).""" b = int(np.clip(np.digitize(imag, MAG_EDGES) - 1, 0, len(MAG_EDGES) - 2)) return ID_STATS[b, 1] / max(ID_STATS[b, 0], 1) CONDS = ("any", "ref185", "ref195", "ingress") def injection(hosts, counts, n_host=1000, n_t0=6, blend=0.0, low_frac=ol.LOW_FRAC, seed=11, keep_example=False, ramp=0.0): """Inject total occultations into quiet catalogue stars and rerun the selection, end to end. The event is subtracted from the *raw* difference fluxes; the quality mask is then rebuilt with the same function that is applied to real data (ol.quality_mask), so that any flux-dependent loss of points in the load path shows up here. An event counts as recovered if it passes cuts C1-C7 and the emulated variable finder (ol.getvar_like). For a duration D the ingress time is uniform in [t_first - D, t_last]; the effective exposure is T_eff(D) = (t_last - t_first + D) * P(recovered). Hosts are weighted so that their bright-state magnitude distribution matches the OGLE-II maps (I < 18). blend = fraction of the catalogued flux that does not belong to the occulted star and therefore remains during the event. Four versions of T_eff are returned: any recovered, no further condition (as in the first version of this study) ref185 recovered, and the star is present in the reference image: its flux there is the stack-weighted mean over the reference frames (weights from ol.reference_weights); the corresponding magnitude I_ref must be brighter than 18.5, and the event is weighted by the ratio of the centroid-cut pass fractions at I_ref and at I ref195 same with I_ref < 19.5 ingress as ref185, and in addition the star was visible on the first frame (the event did not begin before the survey), i.e. only events whose ingress falls inside the survey""" rng = np.random.default_rng(seed) hosts = pd.DataFrame(hosts, columns=["field", "idx", "I"]) hosts["bin"] = np.digitize(hosts.I, MAG_EDGES) - 1 nb_ = len(MAG_EDGES) - 1 nmap = np.array([counts[f"n_I_{MAG_EDGES[k]}_{MAG_EDGES[k + 1]}"].sum() for k in range(nb_)], float) per_bin = max(1, n_host // nb_) pick = pd.concat([hosts[hosts.bin == b].sample(min(per_bin, int((hosts.bin == b).sum())), random_state=seed + b) for b in range(nb_)]) wbin = nmap / nmap.sum() nb = pick.bin.value_counts().reindex(range(nb_)).fillna(0).values pick["w"] = [wbin[b] / nb[b] for b in pick.bin] res = {(D, c): 0.0 for D in DUR_GRID for c in CONDS} # sum of w * T_eff resbin = {(D, b, c): 0.0 for D in DUR_GRID for b in range(nb_) for c in CONDS} nbin = np.zeros(nb_) wsum = 0.0 diag = dict(n_injected=0, n_pass_cuts=0, n_pass_cuts_fail_getvar=0, n_in_event_points_before=0, n_in_event_points_after=0) example = None for (fi, grp) in pick.groupby("field"): cat, t, F, E, flg, _, _ = ol.load_field_raw(int(fi)) ref = np.array(cat["REF_FLUX"], float) zp = np.array(cat["MAG_0"], float) idr = np.array(cat["ID_RAD"], float) wref = REFW[int(fi)] for h in grp.itertuples(): k = h.idx tk, Fk, Ek, fk = t[k], F[k], E[k], flg[k] gk = ol.quality_mask(Fk, Ek, fk) Bk = np.percentile(Fk[gk], 85) tot = ref[k] + Bk t0, t1 = tk.min(), tk.max() id0 = idfrac_at(h.I) wsum += h.w nbin[h.bin] += 1 for D in DUR_GRID: span = (t1 - t0) + D tin = rng.uniform(t0 - D, t1, n_t0) acc = dict.fromkeys(CONDS, 0.0) for ti in tin: cov = cover_fraction(tk, ti, D, ramp) diag["n_injected"] += 1 if (cov > 0).sum() < 3: continue F2 = inject(Fk, cov, tot * (1.0 - blend)) g2 = ol.quality_mask(F2, Ek, fk) # the mask is rebuilt after injection diag["n_in_event_points_before"] += int((gk & (cov > 0)).sum()) diag["n_in_event_points_after"] += int((g2 & (cov > 0)).sum()) m = ol.metrics(tk, F2, Ek, g2, ref[k], zp[k], low_frac=low_frac) if m is None or not ol.passes_all(m, idr[k]): continue diag["n_pass_cuts"] += 1 if not ol.getvar_like(tk, F2, Ek, g2): diag["n_pass_cuts_fail_getvar"] += 1 continue p_ref = blend + (1.0 - blend) * (1.0 - float(np.sum(wref * cov))) i_ref = h.I - 2.5 * np.log10(max(p_ref, 1e-9)) acc["any"] += 1.0 for c, lim in I_REF_LIMITS.items(): if i_ref < lim: acc[c] += min(1.0, idfrac_at(i_ref) / id0) if ti > t0 and i_ref < I_REF_LIMITS["ref185"]: acc["ingress"] += min(1.0, idfrac_at(i_ref) / id0) if (keep_example and example is None and D == 80 and 17.0 < h.I < 17.6 and ti > t0 + 30 and ti + D < t1 - 30 and np.max((ref[k] + Fk[gk]) / tot) < 1.2): example = dict(field=int(fi), idx=int(k), I=float(h.I), t=tk[g2], r=(ref[k] + F2[g2]) / m["tot"], e=Ek[g2] / m["tot"], t_in=float(ti), D=float(D)) for c in CONDS: teff = span * acc[c] / n_t0 / 365.25 res[(D, c)] += h.w * teff resbin[(D, h.bin, c)] += teff print(f"[ogle2] injection (low<{low_frac:g}, blend={blend:g}, ramp={ramp:g} d) field {fi} done ({len(grp)} hosts)", flush=True) rows = [] for D in DUR_GRID: row = {"low_frac": low_frac, "blend_residual": blend, "ramp_d": ramp, "duration_d": D, "n_hosts": int(len(pick)), "n_t0": n_t0} for c in CONDS: tag = "" if c == "any" else f"_{c}" row[f"T_eff_yr{tag}"] = res[(D, c)] / wsum for b in range(nb_): row[f"T_eff_yr{tag}_I_{MAG_EDGES[b]}_{MAG_EDGES[b + 1]}"] = resbin[(D, b, c)] / max(nbin[b], 1) rows.append(row) return pd.DataFrame(rows), example, diag def regression_checks(hosts): """Checks that would have caught the loader defect of the first version (see ol.quality_mask). 1. A total occultation injected into the raw fluxes of an I = 16 star, and passed through the same quality mask as real data, must be recovered, and no in-event point may be lost to the mask. The same event must NOT be recoverable with the defective mask (F > -90), which proves that the check is sensitive to the defect. 2. The loaded light curve of the R CrB star Terz V 2637 (sc43 #2786), which faded to zero flux in 1998, must contain that decline as a run of consecutive low points.""" out = {} h16 = [h for h in hosts if 15.9 < h[2] < 16.1][:20] n_ok = n_old = n_lost = n_pts = 0 loaded = None for fi, k, imag in h16: if loaded != fi: cat, t, F, E, flg, _, _ = ol.load_field_raw(fi) loaded = fi ref, zp, idr = float(cat["REF_FLUX"][k]), float(cat["MAG_0"][k]), float(cat["ID_RAD"][k]) gk = ol.quality_mask(F[k], E[k], flg[k]) tot = ref + np.percentile(F[k][gk], 85) cov = cover_fraction(t[k], 900.0, 160.0) # 160 d, inside the 1998 season F2 = inject(F[k], cov, tot) g2 = ol.quality_mask(F2, E[k], flg[k]) n_pts += int((gk & (cov > 0)).sum()) n_lost += int((gk & (cov > 0) & ~g2).sum()) m = ol.metrics(t[k], F2, E[k], g2, ref, zp) n_ok += int(m is not None and ol.passes_all(m, idr) and ol.getvar_like(t[k], F2, E[k], g2)) g_old = g2 & (F2 > -90) # the defective mask of the first version m_old = ol.metrics(t[k], F2, E[k], g_old, ref, zp) n_old += int(m_old is not None and ol.passes_all(m_old, idr)) out["injected_I16_160d"] = {"n_hosts": len(h16), "n_recovered": n_ok, "n_recovered_with_defective_mask": n_old, "in_event_points": n_pts, "in_event_points_lost_to_mask": n_lost} assert len(h16) >= 10 and n_ok >= 0.9 * len(h16), f"regression check 1 failed: {out}" assert n_lost == 0, f"quality mask removed injected points: {out}" fi, k = POSITIVE_CONTROL cat, t, F, E, good = ol.load_field(fi) ref, zp, idr = float(cat["REF_FLUX"][k]), float(cat["MAG_0"][k]), float(cat["ID_RAD"][k]) m = ol.metrics(t[k], F[k], E[k], good[k], ref, zp) assert m is not None and m["nrun"] >= 20 and m["r_run_med"] < 0.2 and 950 < m["t1"] < 1050, f"positive control missing: {m}" g_old = good[k] & (F[k] > -90) m_old = ol.metrics(t[k], F[k], E[k], g_old, ref, zp) fl = ol.cut_flags(m, idr) out["positive_control_TerzV2637"] = { "field": fi, "idx": k, "ra": float(cat["RA"][k]), "dec": float(cat["DEC"][k]), "I_bright": float(m["I_bright"]), "n_consecutive_low_points": int(m["nrun"]), "run_start_hjd": float(m["t1"]), "run_duration_d": float(m["dur"]), "median_flux_in_run_over_bright": float(m["r_run_med"]), "pre_med": float(m["pre_med"]), "post_med": float(m["post_med"]), "n_low_points_with_defective_mask": int(m_old["nlow"]) if m_old is not None else 0, "cuts": {name.split(" ")[0]: bool(v) for name, v in fl.items()}, "lc": (t[k][good[k]], (ref + F[k][good[k]]) / m["tot"], E[k][good[k]] / m["tot"])} return out def reference_tables(): """Write the recovered reference-image epochs and compare with the frame named in the header of each published reference image (its coordinate template, which is one of the stacked frames).""" from astropy.io import fits import re ep, summ = [], [] for r in REF_ROWS: w, frames, jd = r["w"], r["frames"], r["jd"] hdr_frame, hdr_w, hdr_rank = None, np.nan, np.nan path = os.path.join(RAW, "ogle2_ref", f"sc{r['field']}_ref_1_1_u.fts.gz") if os.path.exists(path): try: hist = " ".join(str(x) for x in fits.getheader(path).get("HISTORY", [])) mm = re.search(r"of (mr\d+)\.fts", hist) if mm: hdr_frame = mm.group(1) j = np.where(frames == hdr_frame)[0] if len(j): hdr_w = float(w[j[0]]) hdr_rank = int((w > w[j[0]]).sum() + 1) except Exception as e: # noqa: BLE001 print(f"[ogle2] could not read {path}: {e!r}") for j in np.where(w > 0.005)[0]: ep.append({"field": r["field"], "frame": frames[j], "hjd_minus_2450000": jd[j], "weight": w[j]}) summ.append({k: v for k, v in r.items() if k not in ("frames", "jd", "w")} | {"header_template_frame": hdr_frame, "header_frame_weight": hdr_w, "header_frame_rank": hdr_rank}) ep, summ = pd.DataFrame(ep), pd.DataFrame(summ) ep.to_csv(os.path.join(DATA, "ogle2_reference_epochs.csv"), index=False, float_format="%.5g") summ.to_csv(os.path.join(DATA, "ogle2_reference_summary.csv"), index=False, float_format="%.4g") return summ def figures(lcs_base, c4, lcs_var, var, eff, gs, ctrl, example, cf, pc): def lc(ax, key, title, ylim=(-0.35, 1.6)): t, r, e = key ax.errorbar(t, r, e, fmt="k.", ms=2.5, lw=0.3, ecolor="0.6") ax.axhline(0.0, color="b", lw=0.5) ax.axhline(0.5, color="r", lw=0.5, ls=":") ax.set_ylim(*ylim) ax.set_title(title, fontsize=8.5) fig, ax = plt.subplots(2, 3, figsize=(16.5, 8.6)) if example is not None: lc(ax[0, 0], (example["t"], example["r"], example["e"]), f"(a) Injected total occultation, {example['D']:.0f} d, in a real I = {example['I']:.1f} star\n" "(injected into the raw fluxes, masked like real data, recovered)") lc(ax[0, 1], pc["lc"], f"(b) Positive control: R CrB star Terz V 2637 (sc{pc['field']} #{pc['idx']}, I = {pc['I_bright']:.1f})\n" "its 1998 fade to zero is recorded; rejected for slow edges and variability", ylim=(-0.6, 1.6)) step = c4[c4["class"] == "step across a seasonal gap"].sort_values("dur") if len(step): r = step.iloc[len(step) // 2] lc(ax[0, 2], lcs_base[(int(r.field), int(r.idx))], f"(c) Real: step across the 1998/99 gap next to a bright star (artefact)\nsc{int(r.field)} #{int(r.idx)}, I = {r.I_bright:.1f}") for a in ax[0]: a.set_xlabel("HJD - 2450000") a.set_ylabel("total flux / bright level") a = ax[1, 0] s = eff[(eff.low_frac == 0.5) & (eff.blend_residual == 0.0) & (eff.ramp_d == 0.0)] for col, c, ls, lab in (("T_eff_yr", "0.5", ":", "recovered by the cuts (no reference-image condition)"), ("T_eff_yr_ref195", "#1f77b4", "-", "and star in the reference image at I < 19.5 (optimistic)"), ("T_eff_yr_ref185", "k", "-", "and star in the reference image at I < 18.5 (conservative)"), ("T_eff_yr_ingress", "#2ca02c", "--", "I_ref < 18.5 and ingress inside the survey")): a.semilogx(s.duration_d, s[col], ls, color=c, marker="o", ms=3, label=lab) for (lf, bl, rp), c, ls, lab in (((0.5, 0.4, 0.0), "#9467bd", "-.", "40% residual blend light (I_ref < 18.5)"), ((0.5, 0.0, 10.0), "#ff7f0e", "-.", "10-day ingress and egress (I_ref < 18.5)"), ((VARIANT_LOW, 0.5, 0.0), "#d62728", "-.", "looser variant, low < 0.7, 50% blend (I_ref < 18.5)")): v = eff[(eff.low_frac == lf) & (eff.blend_residual == bl) & (eff.ramp_d == rp)] if len(v): a.semilogx(v.duration_d, v.T_eff_yr_ref185, ls, color=c, lw=1, label=lab) a.set_xlabel("occultation duration [d]") a.set_ylabel("effective exposure per star T$_{eff}$ [yr]") a.set_title("(d) Injection-recovery: effective exposure per monitored star", fontsize=9) a.legend(fontsize=6.6, loc="upper left") a.set_ylim(0, 3.6) a.grid(alpha=0.3, which="both") a = ax[1, 1] g2 = gs[gs.step_between_1998_and_1999] for arr, lab, c in ((np.sort(g2.sep_bright_arcsec.values), f"{len(g2)} steps across the 1998/99 gap", "#d62728"), (np.sort(ctrl), "random catalogue variables", "0.4")): if len(arr): a.semilogx(arr, (np.arange(len(arr)) + 1) / len(arr), color=c, label=lab, drawstyle="steps-post") a.set_xlabel("separation from nearest Gaia G < 10.5 star [arcsec]") a.set_ylabel("cumulative fraction") a.set_title("(e) The step-like disappearances hug bright stars", fontsize=9) a.legend(fontsize=8) a.grid(alpha=0.3, which="both") a = ax[1, 2] a.axis("off") txt = "(f) Cut flow, baseline search\n\n" + "\n".join( f"{n:>8,d} {c.split('(')[0].strip()[:58]}" for c, n in zip(cf.cut, cf.n_remaining)) a.text(0.0, 1.0, txt, family="monospace", fontsize=8.2, va="top", transform=a.transAxes) fig.tight_layout() fig.savefig(os.path.join(FIGS, "fig3_ogle2_search.png"), dpi=140) plt.close(fig) # figure 5: the near misses. All objects passing C1-C5, all gap steps and monotonic drifts, the # deepest smooth variables, and every survivor of the looser variant (the full list of objects # passing C1-C4 is in data/ogle2_c4_survivors.csv). c4 = c4.copy() c4["group"] = np.where(c4.c5, "0 passes C5, fails C6", c4["class"]) smooth = c4[(c4.group == "smooth variation")].sort_values("r_run_med").head(14) c4s = pd.concat([c4[c4.group == "0 passes C5, fails C6"].sort_values("I_bright"), c4[c4.group == "step across a seasonal gap"].sort_values("dur"), c4[c4.group == "monotonic drift through the whole light curve"].sort_values("dur"), smooth]).reset_index(drop=True) n = len(c4s) + len(var) ncol = 6 nrow = int(np.ceil(n / ncol)) fig, ax = plt.subplots(nrow, ncol, figsize=(3.3 * ncol, 2.15 * nrow)) ax = ax.ravel() j = 0 short = {"0 passes C5, fails C6": "passes C5, fails C6", "step across a seasonal gap": "gap step", "monotonic drift through the whole light curve": "monotonic drift", "smooth variation": "deep smooth variable"} for _, r in c4s.iterrows(): key = (int(r["field"]), int(r["idx"])) lc(ax[j], lcs_base[key], f"sc{key[0]} #{key[1]} I={r['I_bright']:.1f}: {short[r['group']]}", ylim=(-0.5, 1.6)) j += 1 for _, r in var.iterrows(): key = (int(r["field"]), int(r["idx"])) lc(ax[j], lcs_var[key], f"VARIANT sc{key[0]} #{key[1]} I={r['I_bright']:.1f}") ax[j].axhline(VARIANT_LOW, color="r", lw=0.5, ls="--") j += 1 for a in ax[j:]: a.axis("off") for a in ax[:j]: a.tick_params(labelsize=6.5) fig.suptitle("OGLE-II near misses: objects passing C1-C5 (all fail the quiet-outside cut C6), steps across seasonal gaps, monotonic drifts, " "the deepest smooth variables passing C1-C4, and the survivors of the looser variant", fontsize=10) fig.tight_layout(rect=(0, 0, 1, 0.985)) fig.savefig(os.path.join(FIGS, "fig5_ogle2_survivors.png"), dpi=100) plt.close(fig) def main(): quick = "--quick" in sys.argv counts = star_counts() cb = bright_catalogue() # ---- B. baseline search pre, cf, fields, hosts, lcs = search() fields = fields.merge(counts, on="field") fields["n_I_lt_18"] = fields[[f"n_I_{MAG_EDGES[k]}_{MAG_EDGES[k + 1]}" for k in range(len(MAG_EDGES) - 1)]].sum(axis=1) fields.to_csv(os.path.join(DATA, "ogle2_fields.csv"), index=False, float_format="%.6g") cf.to_csv(os.path.join(DATA, "ogle2_cutflow.csv"), index=False) pre.to_csv(os.path.join(DATA, "ogle2_prefilter_metrics.csv"), index=False, float_format="%.6g") cand = pre[pre[[f"c{j}" for j in range(1, 8)]].all(axis=1)].copy() cand.to_csv(os.path.join(DATA, "ogle2_candidates.csv"), index=False, float_format="%.6g") print(cf.to_string(index=False)) print(f"[ogle2] candidates passing all cuts: {len(cand)}") c4 = classify_c4(pre, lcs) c4["sep_bright_arcsec"] = bright_sep(c4.ra.values, c4.dec.values, cb) c4.to_csv(os.path.join(DATA, "ogle2_c4_survivors.csv"), index=False, float_format="%.6g") print("[ogle2] objects passing C1-C4:", c4["class"].value_counts().to_dict()) c5pass = c4[c4.c5] refsum = reference_tables() print(refsum.drop(columns=["n_objects_used"]).to_string(index=False, float_format=lambda x: f"{x:.3g}")) # ---- C. step artefacts gs = gap_steps(pre, cb) gs.to_csv(os.path.join(DATA, "ogle2_gap_steps.csv"), index=False, float_format="%.6g") ctrl = control_bright_sep(cb) from scipy.stats import ks_2samp g9899 = gs[gs.step_between_1998_and_1999] ks = ks_2samp(g9899.sep_bright_arcsec, ctrl) if len(g9899) > 3 else None step_summary = { "n_sustained_faint_states_failing_C5": int(len(gs)), "n_step_down_in_1998_1999_gap": int(len(g9899)), "n_of_those_with_id_rad_gt_1.5px": int((g9899.id_rad > 1.5).sum()), "median_sep_to_G<10.5_star_arcsec": float(g9899.sep_bright_arcsec.median()) if len(g9899) else None, "max_sep_to_G<10.5_star_arcsec": float(g9899.sep_bright_arcsec.max()) if len(g9899) else None, "control_median_sep_arcsec": float(np.median(ctrl)), "frac_within_60arcsec_of_bright_star": float((g9899.sep_bright_arcsec < 60).mean()) if len(g9899) else None, "control_frac_within_60arcsec": float((ctrl < 60).mean()), "ks_p_value": float(ks.pvalue) if ks else None, "run_start_season_counts": {int(k): int(v) for k, v in gs.season_of_run_start.value_counts().sort_index().items()}, } print("[ogle2] gap steps:", json.dumps(step_summary, indent=1)) # ---- D. looser variant var, cfv, _, _, lcs_var = search(low_frac=VARIANT_LOW, keep="c1c5", collect_hosts=False) cfv.to_csv(os.path.join(DATA, "ogle2_variant_cutflow.csv"), index=False) c5 = var.copy() var = var[var.c6 & var.c7].copy() var["sep_bright_arcsec"] = bright_sep(var.ra.values, var.dec.values, cb) if len(var) else [] # how many C1-C5 objects (any field) start their low state on the same night (common-mode test) var["n_other_c5_objects_starting_within_1d"] = [int((np.abs(c5.t1 - t1) < 1.0).sum() - 1) for t1 in var.t1] var["run_ends_with_season"] = ~((var.tnext - var.t2) <= ol.CUTS["edge_gap"]) var["visual_note"] = [VISUAL_NOTES.get((int(f), int(k)), "not inspected") for f, k in zip(var.field, var.idx)] var.to_csv(os.path.join(DATA, "ogle2_variant_survivors.csv"), index=False, float_format="%.6g") print(cfv.to_string(index=False)) pd.set_option("display.width", 250) pd.set_option("display.max_colwidth", 120) if len(var): print(var[["field", "idx", "ra", "dec", "I_bright", "nrun", "dur", "t1", "r_run_med", "sep_bright_arcsec", "n_other_c5_objects_starting_within_1d", "run_ends_with_season", "visual_note"]].to_string(index=False)) # ---- E. regression checks of the load path reg = regression_checks(hosts) pc = reg["positive_control_TerzV2637"] print("[ogle2] regression checks:", json.dumps({k: {kk: vv for kk, vv in v.items() if kk != "lc"} for k, v in reg.items()}, indent=1)) fstat = {"epochs_inside_sustained_low_runs": int(FLAG_STATS[0, 0]), "epochs_outside": int(FLAG_STATS[1, 0]), "frac_error_code_inside": float(FLAG_STATS[0, 1] / max(FLAG_STATS[0, 0], 1)), "frac_error_code_outside": float(FLAG_STATS[1, 1] / max(FLAG_STATS[1, 0], 1)), "frac_flags_3_to_11_inside": float(FLAG_STATS[0, 2] / max(FLAG_STATS[0, 0], 1)), "frac_flags_3_to_11_outside": float(FLAG_STATS[1, 2] / max(FLAG_STATS[1, 0], 1)), "sample": "all epochs of the catalogue objects with a bright state I < 18 and a run of >= 5 consecutive low points (C1 and C3)"} print("[ogle2] quality flags inside and outside real faint states:", json.dumps(fstat, indent=1)) # ---- F. efficiency nh, nt = (300, 4) if quick else (1000, 6) effs, diags = [], {} e0, example, diags["baseline"] = injection(hosts, counts, n_host=nh, n_t0=nt, keep_example=True) effs.append(e0) nv = 150 if quick else 300 for lf, bl, rp in ((0.5, 0.3, 0.0), (0.5, 0.4, 0.0), (VARIANT_LOW, 0.0, 0.0), (VARIANT_LOW, 0.5, 0.0), (0.5, 0.0, 4.0), (0.5, 0.0, 10.0)): e, _, d = injection(hosts, counts, n_host=nv, n_t0=4, blend=bl, low_frac=lf, seed=23, ramp=rp) effs.append(e) diags[f"low<{lf:g} blend={bl:g} ramp={rp:g}d"] = d # same hosts and epochs as the variants, but instantaneous edges and no blend: reference for ratios e, _, d = injection(hosts, counts, n_host=nv, n_t0=4, blend=0.0, low_frac=0.5, seed=23, ramp=0.0) e["n_t0"] = 4 e["low_frac"] = 0.5 e["blend_residual"] = -1.0 # marker: comparison run for the ramp and blend variants (blend = 0) effs.append(e) eff = pd.concat(effs) eff.to_csv(os.path.join(DATA, "ogle2_efficiency.csv"), index=False, float_format="%.5g") print(eff[["low_frac", "blend_residual", "ramp_d", "duration_d", "T_eff_yr", "T_eff_yr_ref195", "T_eff_yr_ref185", "T_eff_yr_ingress", "n_hosts"]].to_string(index=False, float_format=lambda x: f"{x:.3g}")) print("[ogle2] injection diagnostics:", json.dumps(diags, indent=1)) ok_v = fields.frac_I18_with_V > 0.85 summary = { "n_light_curves": int(cf.n_remaining.iloc[0]), "n_candidates": int(len(cand)), "n_stars_I_lt_18": int(fields.n_I_lt_18.sum()), "n_stars_I_lt_17": int(fields[["n_I_0_15", "n_I_15_16", "n_I_16_17"]].sum().sum()), "n_stars_map_total": int(fields.n_map_total.sum()), "n_hosts_available": len(hosts), "n_pass_C1_to_C4": int(len(c4)), "c4_classes": c4["class"].value_counts().to_dict(), "n_pass_C1_to_C5": int(len(c5pass)), "c5_passers": {"I_bright_range": [float(c5pass.I_bright.min()), float(c5pass.I_bright.max())] if len(c5pass) else None, "r_run_med_range": [float(c5pass.r_run_med.min()), float(c5pass.r_run_med.max())] if len(c5pass) else None, "mad_out_range": [float(c5pass.mad_out.min()), float(c5pass.mad_out.max())] if len(c5pass) else None, "n_passing_C6": int(c5pass.c6.sum()) if len(c5pass) else 0}, "variant_low_frac": VARIANT_LOW, "variant_n_pass_C1_to_C5": int(len(c5)), "variant_n_pass_all": int(len(var)), "frac_I18_bluer_than_clump_median_fields_with_V": float(fields.frac_I18_bluer_than_clump[ok_v].median()), "frac_I18_bluer_than_clump_range": [float(fields.frac_I18_bluer_than_clump[ok_v].min()), float(fields.frac_I18_bluer_than_clump[ok_v].max())], "n_fields_with_V_for_85pct": int(ok_v.sum()), "frac_catalogue_objects_with_id_rad_le_cut_by_bin": {f"{MAG_EDGES[b]}_{MAG_EDGES[b + 1]}": float(ID_STATS[b, 1] / max(ID_STATS[b, 0], 1)) for b in range(len(MAG_EDGES) - 1)}, "frac_catalogue_objects_with_id_rad_le_cut_all": float(ID_STATS[:, 1].sum() / ID_STATS[:, 0].sum()), "wozniak2001_recovery_of_independent_microlensing_events": {"recovered": 174, "of": 214, "well_measured_missed": 6}, "gap_steps": step_summary, "regression_checks": {k: {kk: vv for kk, vv in v.items() if kk != "lc"} for k, v in reg.items()}, "quality_flags_in_real_faint_states": fstat, "variable_finder_emulation": {"catalogue_objects_tested": int(GETVAR_STATS[0]), "passing": int(GETVAR_STATS[1]), "fraction_passing": float(GETVAR_STATS[1] / max(GETVAR_STATS[0], 1))}, "injection_diagnostics": diags, "reference_image": { "n_fields": int(len(refsum)), "median_n_frames_weight_gt_0.01": float(refsum["n_frames_weight_gt_0.01"].median()), "mean_weight_by_season": {"1997": float(refsum.weight_1997.mean()), "1998": float(refsum.weight_1998.mean()), "1999": float(refsum.weight_1999.mean())}, "median_rms_stack_mean_flux_over_amplitude": float(refsum.rms_stack_mean_flux_over_amplitude.median()), "median_rms_flat_mean_flux_over_amplitude": float(refsum.rms_flat_mean_flux_over_amplitude.median()), "n_fields_with_header_frame": int(refsum.header_template_frame.notna().sum()), "n_fields_header_frame_recovered_weight_gt_0.01": int((refsum.header_frame_weight > 0.01).sum()), "i_ref_limits": I_REF_LIMITS}, "cuts": ol.CUTS, "quick": quick, } with open(os.path.join(DATA, "ogle2_search_summary.json"), "w") as f: json.dump(summary, f, indent=1) figures(lcs, c4, lcs_var, var, eff, gs, ctrl, example, cf, pc) if __name__ == "__main__": main()