#!/usr/bin/env python3 """How well do the stellar age scales agree? (added in revision) The fraction of stars older than the Sun is read off the cumulative age distribution where it is steep, so it inherits any error of the absolute age scale. This module compares, star by star, the ages of the catalogues used in the study: * Xiang & Rix 2022 (X&R; LAMOST spectroscopy + YY isochrones) * Nataf+ 2024 (broadband photometry + Gaia parallaxes, MIST isochrones) * StarHorse (Queiroz+ 2023; Gaia RVS sample, PARSEC isochrones) * Gaia DR3 FLAME * this study's own Gaia-only photometric ages (analysis_gaia_sg.py) Outputs data/age_scale_offsets.csv median difference and ratio for every pair, overall and in bins of the reference age data/age_swap_test.csv the LAMOST estimator re-run on the stars that also have Nataf+24 ages, with the two sets of ages swapped (same stars, same weights) """ import os import numpy as np import pandas as pd import common as C import fetch_data as F BINS = [(2, 3), (3, 4.567), (4.567, 6), (6, 8), (8, 10), (10, 14)] OUT = C.DATA + "/age_scale_offsets.csv" def _match(ra1, dec1, ra2, dec2, tol_arcsec): """Nearest neighbour of every (ra1, dec1) in catalogue 2; returns (index into 2, boolean 'matched').""" from astropy.coordinates import SkyCoord import astropy.units as u c1 = SkyCoord(np.asarray(ra1) * u.deg, np.asarray(dec1) * u.deg) c2 = SkyCoord(np.asarray(ra2) * u.deg, np.asarray(dec2) * u.deg) idx, sep, _ = c1.match_to_catalog_sky(c2) return idx, sep.arcsec < tol_arcsec def pair_rows(pair, ref_name, a_ref, a_other, note=""): """Rows: all stars, then bins of the reference age.""" a_ref, a_other = np.asarray(a_ref, float), np.asarray(a_other, float) rows = [] def row(subset, m): d = a_other[m] - a_ref[m] if m.sum() < 5: return rows.append(dict( pair=pair, reference=ref_name, subset=subset, n=int(m.sum()), median_diff_gyr=float(np.median(d)), median_ratio=float(np.median(a_other[m] / a_ref[m])), robust_scatter_gyr=float(1.4826 * np.median( np.abs(d - np.median(d)))), frac_older_than_sun_ref=float((a_ref[m] > C.T_SUN).mean()), frac_older_than_sun_other=float((a_other[m] > C.T_SUN).mean()), note=note)) ok = np.isfinite(a_ref) & np.isfinite(a_other) & (a_ref > 0) row("all", ok & (a_ref >= 2.0)) for lo, hi in BINS: row(f"{ref_name} age {lo:g}-{hi:g} Gyr", ok & (a_ref >= lo) & (a_ref < hi)) return rows def load_xr(): x = pd.read_csv(C.RAW + "xiangrix2022.csv").rename( columns={"[Fe/H]": "feh", "[alpha/Fe]": "afe"}) x["d"] = np.sqrt((x.X - 8.178) ** 2 + x.Y ** 2 + x.Z ** 2) return x def load_nataf(): from astropy.table import Table t = Table.read(F.NATAF_FITS).to_pandas() t["Sample"] = t.Sample.str.decode("ascii").str.strip() t["age_gyr"] = 10 ** (t.age - 9.0) t["quality"] = ((t.Blend == 0) & (t.VarExcess < 0.95) & ~(t.Chance_pvalue <= 0.001)) return t def xr_in_gaia_parent(x): """X&R stars within 230 pc matched by position to the Gaia 200 pc parent sample (positions propagated from epoch 2016 to J2000).""" xs = x[x.d < 0.23].reset_index(drop=True) p = pd.read_csv(C.RAW + "gaia_parent_200pc.csv") dt = -16.0 dec = p.dec + dt * p.pmdec.fillna(0) / 3.6e6 ra = p.ra + dt * p.pmra.fillna(0) / 3.6e6 / np.cos(np.radians(p.dec)) idx, ok = _match(xs.RAJ2000, xs.DEJ2000, ra, dec, 2.0) out = xs[ok].reset_index(drop=True) g = p.iloc[idx[ok]].reset_index(drop=True) out["source_id"] = g.source_id.values out["age_flame"] = g.age_flame.values return out def catalogue_pairs(): """Pairs that do not involve this study's own Gaia-only ages.""" rows = [] x = load_xr() have_nataf = os.path.exists(F.NATAF_FITS) if have_nataf: n = load_nataf() idx, ok = _match(x.RAJ2000, x.DEJ2000, n.ra2000, n.dec2000, 1.0) nm = n.iloc[idx].reset_index(drop=True) ok = ok & nm.quality.values rows += pair_rows("Nataf+24 vs X&R22", "X&R", x.Age.values[ok], nm.age_gyr.values[ok], "positional match < 1 arcsec, Nataf quality flags") lo_a = ok & (x.afe.values < 0.10) hi_a = ok & (x.afe.values > 0.16) for lab, m in (("low-alpha ([alpha/Fe] < 0.10)", lo_a), ("high-alpha ([alpha/Fe] > 0.16)", hi_a)): r = pair_rows("Nataf+24 vs X&R22", "X&R", x.Age.values[m], nm.age_gyr.values[m])[0] r["subset"] = lab rows.append(r) xg = xr_in_gaia_parent(x) sh = pd.read_csv(C.RAW + "starhorse_rvs_210pc.csv") ms = xg.merge(sh, left_on="source_id", right_on="GaiaEDR3", how="inner") ms = ms[ms.Age50.notna()] rows += pair_rows("StarHorse vs X&R22", "X&R", ms.Age, ms.Age50, "X&R stars within 230 pc") t = xg[xg.age_flame.notna()] rows += pair_rows("Gaia FLAME vs X&R22", "X&R", t.Age, t.age_flame, "X&R stars within 230 pc") if have_nataf: mn = sh.merge(n[n.quality], left_on="GaiaEDR3", right_on="dr3_source_id", how="inner") mn = mn[mn.Age50.notna()] rows += pair_rows("StarHorse vs Nataf+24", "Nataf", mn.age_gyr, mn.Age50, "StarHorse RVS stars within 210 pc") return rows def gaia_only_pairs(b, p_min=0.5): """Pairs involving this study's Gaia-only photometric ages. b: the Gaia subgiant sample (source_id, age, p_sg). Returns (rows, ratios), ratios = median(Gaia-only age / other catalogue age).""" rows, ratios = [], {} q = b[b.p_sg > p_min] x = load_xr() xg = xr_in_gaia_parent(x) m = xg.merge(q, on="source_id", how="inner", suffixes=("_xr", "")) r = pair_rows("this study's Gaia-only ages vs X&R22", "X&R", m.Age, m.age, f"p_sg > {p_min}") rows += r ratios["X&R"] = r[0]["median_ratio"] if r else np.nan if os.path.exists(F.NATAF_FITS): n = load_nataf() m = q.merge(n[n.quality], left_on="source_id", right_on="dr3_source_id", how="inner", suffixes=("", "_nataf")) r = pair_rows("this study's Gaia-only ages vs Nataf+24", "Nataf", m.age_gyr, m["age"], f"p_sg > {p_min}") rows += r ratios["Nataf"] = r[0]["median_ratio"] if r else np.nan sh = pd.read_csv(C.RAW + "starhorse_rvs_210pc.csv") m = q.merge(sh, left_on="source_id", right_on="GaiaEDR3", how="inner") m = m[m.Age50.notna()] r = pair_rows("this study's Gaia-only ages vs StarHorse", "StarHorse", m.Age50, m.age, f"p_sg > {p_min}") rows += r ratios["StarHorse"] = r[0]["median_ratio"] if r else np.nan return rows, ratios def measured_ratios(default=None): """Median age ratios relative to X&R, read from the offsets table.""" out = dict(default or {}) if os.path.exists(OUT): t = pd.read_csv(OUT) t = t[t.subset == "all"].set_index("pair") for key, pair in (("Nataf", "Nataf+24 vs X&R22"), ("StarHorse", "StarHorse vs X&R22"), ("FLAME", "Gaia FLAME vs X&R22")): if pair in t.index: out[key] = float(t.loc[pair, "median_ratio"]) return out def age_swap_test(): """Same stars, same weights, two age catalogues. The LAMOST shell sample is matched to Nataf+24. For the stars in common the LAMOST estimator (tracer weights x dynamical column weights, EM deconvolution, young fraction 0.10) is evaluated three times: with the X&R ages, with the Nataf+24 ages, and with the Nataf+24 ages multiplied by the inverse of the median Nataf/X&R ratio. The weights are those of the LAMOST analysis in all three cases, so the differences isolate the effect of the age catalogue.""" import sgweights as S import tracer as T import analysis_xr as A if not os.path.exists(F.NATAF_FITS): return None x = A.load() tw = S.TracerWeights("xr") n = load_nataf() s = A.select(x) idx, ok = _match(s.RAJ2000, s.DEJ2000, n.ra2000, n.dec2000, 1.0) nm = n.iloc[idx].reset_index(drop=True) ratio = measured_ratios({"Nataf": 0.95})["Nataf"] rows = [] for lab, sel in (("all stars in common", ok), ("in common, Nataf quality flags", ok & nm.quality.values)): ss = s[sel].reset_index(drop=True) nn = nm[sel].reset_index(drop=True) sm = A.Sample(ss, tw) w = sm.weights() aN = nn.age_gyr.values eN = np.maximum(0.5 * (10 ** (nn.age_84.values - 9) - 10 ** (nn.age_16.values - 9)), 0.05 * aN) for ages, a, e in (("X&R ages", sm.age, sm.err), ("Nataf+24 ages", aN, eN), (f"Nataf+24 ages / {ratio:.3f}", aN / ratio, eN / ratio)): k = a >= T.TAU_MIN st, _ = T.run_variant(a[k], e[k], e[k], w[k]) st2, _ = T.run_variant(a[k], e[k], e[k], (1 / sm.g)[k]) rows.append(dict( subset=lab, n=int(sel.sum()), ages=ages, f_older_column=st["f_gt_0"], f_gt_5_column=st["f_gt_5"], mean_headstart_older=st["mean_headstart_older"], f_older_midplane=st2["f_gt_0"], weight_share_xr_old_other_young=float(np.sum( w * ((sm.age > C.T_SUN) & (a <= C.T_SUN))) / w.sum()), weight_share_xr_young_other_old=float(np.sum( w * ((sm.age <= C.T_SUN) & (a > C.T_SUN))) / w.sum()))) df = pd.DataFrame(rows) df.to_csv(C.DATA + "/age_swap_test.csv", index=False) return df def main(): rows = catalogue_pairs() df = pd.DataFrame(rows) df.to_csv(OUT, index=False) pd.set_option("display.width", 250, "display.max_columns", 30, "display.max_colwidth", 60) print("[ages] catalogue age scales, star by star:") print(df[["pair", "subset", "n", "median_diff_gyr", "median_ratio", "robust_scatter_gyr"]].round(3).to_string(index=False)) return df def append_gaia_only(b): """Called by analysis_gaia_sg once its sample exists.""" rows, ratios = gaia_only_pairs(b) new = pd.DataFrame(rows) if os.path.exists(OUT): old = pd.read_csv(OUT) old = old[~old.pair.str.startswith("this study's Gaia-only")] new = pd.concat([old, new], ignore_index=True) new.to_csv(OUT, index=False) return ratios if __name__ == "__main__": main() print(age_swap_test())