#!/usr/bin/env python3 """Pillar 1c: a volume-complete, Gaia-only subgiant census within 200 pc. The parent sample is every Gaia DR3 source inside a colour - absolute- magnitude box within 200 pc (no spectroscopic targeting), ages are computed here from MIST isochrones using (M_G, BP-RP) and the Gaia GSP-Spec metallicity, and the vertical-orbit weights use Gaia radial velocities directly. STATUS AFTER REVISION. This is a cross-check, not an independent replication. For stars in common, the ages computed here run about 15% older than Xiang & Rix (2022) and about 18% older than Nataf+ (2024), who use the same MIST models with far more photometry (age_scales.py); 22% of the weight piles up at the oldest isochrone. The ages are therefore biased old, and the result is reported both as computed and with the age likelihoods rescaled onto the X&R and Nataf+24 scales. Sample quality: about 15% of the candidates have no Gaia radial velocity and are dropped, and 30-36% of the membership-weighted sample has RUWE > 1.4 (unresolved binaries); both fractions are tabulated by age in gaia_sg_quality.csv. Age estimate (per star) ----------------------- likelihood(tau) = sum over [M/H] grid, sum along the isochrone of N(M_G; model, s_M) N(BP-RP; model, s_c) xi(M) dM x N([M/H]; obs, s_Z) with s_M = quadrature of parallax error and 0.04 mag, s_c = 0.02 mag, s_Z = 0.12 dex for GSP-Spec metallicities. Stars without a GSP-Spec metallicity get the broad prior N(-0.15, 0.30) instead of being dropped, so that no spectroscopic selection is imposed. Extinction is neglected (d < 200 pc); a reddened-box test is part of the robustness table. """ import numpy as np import pandas as pd import common as C import sgweights as S import tracer as T # Parent window (Gaia DR3 photometry, no extinction correction) and the # definition of the tracer in model space: "subgiant" = MIST equivalent # evolutionary points EEP_SG[0]..EEP_SG[1] (454 = terminal-age main sequence; # ~472 = base of the red-giant branch). WINDOW = dict(c_lo=0.55, c_hi=1.25, M_hi=4.8) EEP_SG = (454, 472) DMAX = 200.0 S_COL, S_MAG0, S_MH, MH_PRIOR = 0.02, 0.04, 0.12, (-0.15, 0.30) LOGAGES = np.round(np.arange(9.0, 10.1501, 0.05), 2) # 1 - 14.1 Gyr EEPS = np.arange(380, 521) # neighbours of the SGB: late MS .. lower RGB MH_FINE = np.round(np.arange(-2.0, 0.5001, 0.05), 2) _TWIN = {} def twin_cut(c, margin=0.10): """Upper envelope of unresolved equal-mass main-sequence binaries: 0.753 mag above the main sequence of a metal-rich ([M/H] = +0.25, 1 Gyr) MIST isochrone, minus a margin. Stars (and model points) fainter than this line can be mimicked by binaries.""" if not _TWIN: arr = C.load_mist()[0.25] iso = arr[np.isclose(arr[:, 0], 9.0)] ms = iso[(iso[:, 1] >= 202) & (iso[:, 1] <= 454)] col = (ms[:, 9] - ms[:, 10]).astype(float) o = np.argsort(col) _TWIN["c"], _TWIN["M"] = col[o], ms[o, 8].astype(float) return np.interp(c, _TWIN["c"], _TWIN["M"]) - 0.753 - margin def in_window(c, M, w=None): w = WINDOW if w is None else w ok = (c > w["c_lo"]) & (c < w["c_hi"]) & (M < w["M_hi"]) if w.get("binary_safe"): ok = ok & (M < twin_cut(np.asarray(c, float))) return ok def sg_box_func(dc=0.0, dM=0.0, window=None): """Tracer definition for sgweights.TracerWeights: isochrone points with EEP in the subgiant range that fall inside the parent window. dc, dM shift the *model* colours/magnitudes (zero-point tests).""" def f(iso, feh): c = iso[:, 9].astype(float) - iso[:, 10].astype(float) + dc M = iso[:, 8].astype(float) + dM e = iso[:, 1] return (e >= EEP_SG[0]) & (e <= EEP_SG[1]) & in_window(c, M, window) return f def build_grid(): """Dense model arrays [n_mh_fine, n_age, n_eep] of M_G, colour and the IMF weight xi(M) dM, interpolated in [M/H] at fixed (age, EEP).""" mist = C.load_mist() fehs = np.array(sorted(mist.keys())) shape = (len(fehs), len(LOGAGES), len(EEPS)) MG = np.full(shape, np.nan); CO = np.full(shape, np.nan) MI = np.full(shape, np.nan) for a, feh in enumerate(fehs): arr = mist[feh] for j, la in enumerate(LOGAGES): iso = arr[np.isclose(arr[:, 0], la)] e = iso[:, 1].astype(int) ok = (e >= EEPS[0]) & (e <= EEPS[-1]) k = e[ok] - EEPS[0] MG[a, j, k] = iso[ok, 8] CO[a, j, k] = iso[ok, 9] - iso[ok, 10] MI[a, j, k] = iso[ok, 2] x = np.interp(MH_FINE, fehs, np.arange(len(fehs))) i0 = np.clip(np.floor(x).astype(int), 0, len(fehs) - 2) f = (x - i0)[:, None, None] MGf, COf, MIf = [(1 - f) * A[i0] + f * A[i0 + 1] for A in (MG, CO, MI)] dM = np.zeros_like(MIf) dM[:, :, 1:-1] = 0.5 * (MIf[:, :, 2:] - MIf[:, :, :-2]) dM[:, :, 0] = MIf[:, :, 1] - MIf[:, :, 0] dM[:, :, -1] = MIf[:, :, -1] - MIf[:, :, -2] wimf = np.where(np.isfinite(dM) & (dM > 0), MIf ** (-C.IMF_SLOPE) * dM, 0.0) wimf = np.nan_to_num(wimf) return (np.nan_to_num(MGf, nan=99.0).astype(np.float32), np.nan_to_num(COf, nan=99.0).astype(np.float32), wimf.astype(np.float32)) def prefilter(MGobs, cobs, grid, nsig=3.5, s_M=0.06): """Stars within nsig sigma of any model subgiant point (KD-tree).""" from scipy.spatial import cKDTree MGf, COf, wimf = grid sg = (EEPS >= EEP_SG[0]) & (EEPS <= EEP_SG[1]) pts = np.stack([COf[:, :, sg].ravel() / S_COL, MGf[:, :, sg].ravel() / s_M], axis=1) pts = pts[(wimf[:, :, sg].ravel() > 0) & (pts[:, 0] < 500)] tree = cKDTree(pts[::3]) dist, _ = tree.query(np.stack([cobs / S_COL, MGobs / s_M], axis=1)) return dist < nsig def age_likelihoods(MGobs, cobs, s_M, mh, s_mh, grid, chunk=300): """Per-star likelihoods on the LOGAGES nodes, marginalised over [M/H] and position along the isochrone. Returns L_sg [n, n_age] : likelihood restricted to subgiant EEPs p_sg [n] : posterior probability that the star is a subgiant (flat prior in linear age, IMF prior in mass) mh_post [n] : posterior-mean [M/H] given 'subgiant' """ MGf, COf, wimf = grid sg = (EEPS >= EEP_SG[0]) & (EEPS <= EEP_SG[1]) n = len(MGobs) L_sg = np.zeros((n, len(LOGAGES))) L_all = np.zeros((n, len(LOGAGES))) mh_num = np.zeros(n); mh_den = np.zeros(n) MGobs, cobs, s_M = (np.asarray(v, np.float32) for v in (MGobs, cobs, s_M)) for m, mhv in enumerate(MH_FINE): pz = np.exp(-0.5 * ((mhv - mh) / s_mh) ** 2) / s_mh act = np.where(pz > 1e-3 * np.exp(-0.5 * 0) / s_mh)[0] for a in range(0, len(act), chunk): ii = act[a:a + chunk] chi = (((MGobs[ii, None, None] - MGf[m][None]) / s_M[ii, None, None]) ** 2 + ((cobs[ii, None, None] - COf[m][None]) / np.float32(S_COL)) ** 2) like = np.exp(-0.5 * chi) * wimf[m][None] la = like.sum(axis=2).astype(float) * pz[ii, None] ls = like[:, :, sg].sum(axis=2).astype(float) * pz[ii, None] L_all[ii] += la L_sg[ii] += ls tot = ls.sum(axis=1) mh_num[ii] += tot * mhv mh_den[ii] += tot prior = 10.0 ** (LOGAGES - 9.0) # flat in linear age p_sg = (L_sg * prior).sum(axis=1) / np.maximum((L_all * prior).sum(axis=1), 1e-300) return L_sg, p_sg, mh_num / np.maximum(mh_den, 1e-300) def nodes_to_grid(Ln): """Interpolate node likelihoods (uniform in log age) onto the linear age grid; returns row-normalised L[n, n_bins].""" la = np.log10(np.maximum(C.AGE_MID, 1e-3)) + 9.0 x = np.interp(la, LOGAGES, np.arange(len(LOGAGES))) i0 = np.clip(np.floor(x).astype(int), 0, len(LOGAGES) - 2) f = x - i0 L = (1 - f)[None, :] * Ln[:, i0] + f[None, :] * Ln[:, i0 + 1] L[:, (la < LOGAGES[0]) | (la > LOGAGES[-1])] = 0.0 L = L + 1e-300 return L / L.sum(axis=1, keepdims=True) def percentiles_from_L(L, q=(0.16, 0.5, 0.84)): cs = np.cumsum(L, axis=1) out = [] for qq in q: k = np.clip((cs < qq).sum(axis=1), 0, L.shape[1] - 1) out.append(C.AGE_MID[k]) return out def load_parent(): p = pd.read_csv(C.RAW + "gaia_parent_200pc.csv") p["d"] = 1000.0 / p.parallax p["MG"] = p.phot_g_mean_mag + 5 * np.log10(p.parallax / 100.0) return p def prepare(dc=0.0, dM=0.0, mh_shift=0.0, grid=None, p_min=0.05, window=None): """Select candidate subgiants from the parent sample, compute ages, subgiant probabilities and kinematics. dc, dM are subtracted from the observed colour / absolute magnitude (de-reddening or zero-point tests). Returns (DataFrame, L on the common age grid, info).""" grid = grid or build_grid() p = load_parent() n_raw = len(p) p = p[in_window(p.bp_rp, p.MG, window) & (p.d < DMAX)].reset_index(drop=True) info = dict(n_parent_download=n_raw, n_parent=len(p)) cand = prefilter(p.MG.values - dM, p.bp_rp.values - dc, grid) b = p[cand].reset_index(drop=True) info["n_candidates"] = len(b) has = b.mh_gspspec.notna().values mh = np.where(has, b.mh_gspspec.values + mh_shift, MH_PRIOR[0]) s_mh = np.where(has, S_MH, MH_PRIOR[1]) s_M = np.sqrt(S_MAG0 ** 2 + (2.1715 * b.parallax_error / b.parallax) ** 2).values Ln, p_sg, mh_post = age_likelihoods(b.MG.values - dM, b.bp_rp.values - dc, s_M, mh, s_mh, grid) keep = p_sg > p_min b, Ln, p_sg, mh_post = b[keep].reset_index(drop=True), Ln[keep], p_sg[keep], mh_post[keep] L = nodes_to_grid(Ln) p16, p50, p84 = percentiles_from_L(L) b["age"], b["age_lo"], b["age_hi"] = p50, p16, p84 b["mh_fit"], b["p_sg"] = mh_post, p_sg rv = b.radial_velocity.notna().values W = np.full(len(b), np.nan); z = W.copy() _, _, w_, z_ = C.galactic_uvw(b.l[rv].values, b.b[rv].values, b.parallax[rv].values, b.pmra[rv].values, b.pmdec[rv].values, b.radial_velocity[rv].values, b.ra[rv].values, b.dec[rv].values) W[rv], z[rv] = w_ + C.W_SUN, z_ + C.Z_SUN b["W"], b["z"] = W, z info.update(n_psg_gt_min=len(b), n_subgiants_expected=float(p_sg.sum()), n_rv=int(rv.sum()), n_gspspec=int(b.mh_gspspec.notna().sum())) return b, L, info def starhorse_check(b, p_min=0.7): """Validation only: compare my isochrone ages with StarHorse (Queiroz+ 2023, Gaia RVS sample) for likely subgiants (p_sg > p_min).""" sh = pd.read_csv(C.RAW + "starhorse_rvs_210pc.csv") m = b.merge(sh, left_on="source_id", right_on="GaiaEDR3", how="inner") m = m[m.Age50.notna() & (m.p_sg > p_min)] d = m.age - m.Age50 rows = [dict(subset="all", n=len(m), median_diff_gyr=d.median(), robust_scatter_gyr=1.4826 * np.median(np.abs(d - d.median())), corr=np.corrcoef(m.age, m.Age50)[0, 1], median_this_work=m.age.median())] for lo, hi in [(2, 4), (4, 6), (6, 8), (8, 10), (10, 14)]: q = m[(m.Age50 >= lo) & (m.Age50 < hi)] dq = q.age - q.Age50 rows.append(dict(subset=f"StarHorse age {lo}-{hi} Gyr", n=len(q), median_diff_gyr=dq.median(), robust_scatter_gyr=1.4826 * np.median( np.abs(dq - dq.median())), corr=np.nan, median_this_work=q.age.median())) return pd.DataFrame(rows) def quality_table(b): """Membership-weighted fraction of candidates with a radial velocity and with RUWE > 1.4, by age (all candidates with p_sg > 0.05).""" rv = b.W.notna().values rows = [] for lo, hi in [(1, 2), (2, 3), (3, 4.567), (4.567, 8), (8, 11), (11, 14.1), (2, 14.1)]: m = ((b.age >= lo) & (b.age < hi)).values p = b.p_sg.values[m] rows.append(dict(age_lo=lo, age_hi=hi, n=int(m.sum()), sum_p_sg=float(p.sum()), rv_fraction=float(np.sum(p * rv[m]) / p.sum()), ruwe_gt_1p4_fraction=float( np.sum(p * (b.ruwe.values[m] > 1.4)) / p.sum()))) return pd.DataFrame(rows) class Rescaled: """A GaiaSample whose age likelihoods are mapped onto another age scale: age_new = factor * age_old. Tracer weights and living fractions are re-evaluated at the new ages; membership and column weights are kept.""" def __init__(self, sm, tw, factor, tau_min=T.TAU_MIN): told = C.AGE_MID / factor k = np.interp(told, C.AGE_MID, np.arange(len(C.AGE_MID))) i0 = np.clip(np.floor(k).astype(int), 0, len(C.AGE_MID) - 2) f = k - i0 L2 = (1 - f)[None, :] * sm.L[:, i0] + f[None, :] * sm.L[:, i0 + 1] L2[:, told > C.AGE_MID[-1]] = 0.0 L2 = L2 + 1e-300 L2 = L2 / L2.sum(axis=1, keepdims=True) # keep the stars that are still tracers (>= tau_min) on the new scale keep = (L2 * (C.AGE_MID[None, :] >= tau_min)).sum(axis=1) > 0.5 self.b = sm.b[keep].reset_index(drop=True) self.L = L2[keep] self.age = sm.age[keep] * factor mh = self.b.mh_fit.values I0 = S.imf_integral(0.8, 1.2) self.g = np.maximum(tw.g(self.age, mh) / I0, 0.01) self.wdyn = sm.wdyn[keep] self.alive = S.alive_fraction(tw, self.age, mh) self.psg = sm.psg[keep] weights = None # bound below (same formula as GaiaSample) class GaiaSample: def __init__(self, b, L, tw, pot_name="fiducial", tau_min=T.TAU_MIN, cut=None, dmax=DMAX): """cut: optional boolean Series selecting a sub-sample (quality or distance cuts; unlike a metallicity floor these do not change the ratio of old to young stars).""" use = (b.W.notna() & (b.age >= tau_min)).values if cut is not None: use = use & np.asarray(cut) self.b = b[use].reset_index(drop=True) self.L = L[use] self.pot = C.get_potential(pot_name) self.age = self.b.age.values I0 = S.imf_integral(0.8, 1.2) self.g = np.maximum(tw.g(self.age, self.b.mh_fit.values) / I0, 0.01) E = self.pot.energy(self.b.z.values, self.b.W.values) self.wdyn = 1.0 / self.pot.sphere_fraction(E, dmax) self.alive = S.alive_fraction(tw, self.age, self.b.mh_fit.values) self.psg = self.b.p_sg.values def weights(self, idx=None, living=False, column=True): """membership probability x column weight x tracer weight""" idx = np.arange(len(self.age)) if idx is None else idx w = self.psg[idx] * (self.wdyn[idx] if column else 1.0) / self.g[idx] return w * self.alive[idx] if living else w Rescaled.weights = GaiaSample.weights def variant(sm, name, keep=None, living=False, column=True, n_boot=0, phi_young=T.PHI_YOUNG[0], tau_min=T.TAU_MIN, seed=7): def stats(idx, py): w = sm.weights(idx, living, column) k = np.ones(len(idx), bool) if keep is None else keep[idx] Lk = sm.L[idx][k] * (C.AGE_MID[None, :] >= tau_min) Lk = Lk / np.maximum(Lk.sum(axis=1, keepdims=True), 1e-300) phi_old = C.em_deconvolve(Lk, w[k], n_iter=60) phi = T.combine_with_young(phi_old, py * tau_min / T.TAU_MIN, tau_min, old_keep=w[k].sum() / w.sum()) st = C.headstart_stats(phi) st["f_gt8Gyr"] = C.frac_older(phi, 8.0 - C.T_SUN) st["frac_pass"] = py + (1 - py) * w[k].sum() / w.sum() st["f_last_bin"] = float(phi[-4:].sum()) # pile-up at 13-14 Gyr return st, phi, w full = np.arange(len(sm.age)) st, phi, w = stats(full, phi_young) row = dict(variant=name, n=len(full) if keep is None else int(keep.sum()), n_eff=float(w.sum() ** 2 / (w ** 2).sum()), **st) if n_boot: rng = np.random.default_rng(seed) bs = [] for _ in range(n_boot): idx = rng.integers(0, len(full), len(full)) py = float(np.clip(rng.normal(*T.PHI_YOUNG), 0.02, 0.25)) bs.append(stats(idx, py)[0]) for k_ in st: lo, hi = np.percentile([x[k_] for x in bs], [16, 84]) row[k_ + "_lo"], row[k_ + "_hi"] = lo, hi return row, phi def main(n_boot=100): import age_scales grid = build_grid() tw = S.TracerWeights(sg_box_func()) b, L, info = prepare(grid=grid) print("[gaia-sg]", info) shc = starhorse_check(b) shc.to_csv(C.DATA + "/gaia_sg_vs_starhorse.csv", index=False) print("[gaia-sg] ages vs StarHorse (validation):") print(shc.round(2).to_string(index=False)) qt = quality_table(b) qt.to_csv(C.DATA + "/gaia_sg_quality.csv", index=False) print("[gaia-sg] radial-velocity and RUWE fractions by age:") print(qt.round(3).to_string(index=False)) rows, phis = [], {} def add(name, sm, **kw): r, phi = variant(sm, name, **kw) rows.append(r); phis[name] = phi print(f"[gaia-sg] {name:46s} N={r['n']:5d} f_older={r['f_gt_0']:.3f} " f"f>1={r['f_gt_1']:.3f} f>3={r['f_gt_3']:.3f} f>5={r['f_gt_5']:.3f}" f" ={r['mean_headstart_older']:.2f}") return r sm = GaiaSample(b, L, tw) pd.DataFrame([ dict(step="Gaia DR3 parent download (d < 200 pc, M_G < 4.8, " "0.55 < BP-RP < 1.25)", n=info["n_parent_download"], reason="volume-complete parent, no spectroscopic targeting"), dict(step="inside the analysis window", n=info["n_parent"], reason="colour-magnitude window used for the tracer definition"), dict(step="within 3.5 sigma of a model subgiant point", n=info["n_candidates"], reason="pre-filter for the age fit"), dict(step="subgiant membership probability p_sg > 0.05", n=info["n_psg_gt_min"], reason="membership-weighted candidates"), dict(step="with a Gaia radial velocity", n=info["n_rv"], reason="vertical velocity needed for the column weights " "(fraction flat with age, see gaia_sg_quality.csv)"), dict(step="median age >= 2 Gyr (sample used)", n=len(sm.age), reason="tracer limit"), dict(step="expected number of true subgiants (sum of p_sg)", n=int(round(float(sm.psg.sum()))), reason="information only"), dict(step="of which with GSP-Spec [M/H]", n=int(sm.b.mh_gspspec.notna().sum()), reason="others use a broad metallicity prior"), ]).to_csv(C.DATA + "/gaia_sg_cutflow.csv", index=False) # age scale of these home-made ages against published catalogues cols = ["source_id", "d", "MG", "bp_rp", "mh_gspspec", "mh_fit", "p_sg", "age", "age_lo", "age_hi", "W", "z", "ruwe", "age_flame"] sm.b[cols].round(4).to_csv(C.DATA + "/gaia_sg_sample.csv", index=False) ratios = age_scales.append_gaia_only(sm.b) print("[gaia-sg] median ratio of these ages to published ones:", {k: round(v, 3) for k, v in ratios.items()}) mh = sm.b.mh_fit.values add("formed: all [M/H]", sm, n_boot=n_boot) add("formed: [M/H]>-0.5", sm, keep=mh > -0.5, n_boot=n_boot) add("formed: [M/H]>-0.3", sm, keep=mh > -0.3, n_boot=n_boot) add("living MS: all [M/H]", sm, living=True, n_boot=n_boot) add("no column weights (mid-plane sphere)", sm, column=False) # the same sample with the age likelihoods mapped onto published scales for key, lab in (("X&R", "X&R22"), ("Nataf", "Nataf+24"), ("StarHorse", "StarHorse")): r_ = ratios.get(key, np.nan) if not np.isfinite(r_): continue s2 = Rescaled(sm, tw, 1.0 / r_) nb = n_boot if key == "X&R" else 0 add(f"ages on the {lab} scale: formed, all [M/H]", s2, n_boot=nb)["age_scale"] = 1.0 / r_ add(f"ages on the {lab} scale: [M/H]>-0.5", s2, keep=s2.b.mh_fit.values > -0.5)["age_scale"] = 1.0 / r_ add(f"ages on the {lab} scale: mid-plane", s2, column=False)["age_scale"] = 1.0 / r_ add("only stars with GSP-Spec [M/H]", GaiaSample(b, L, tw, cut=b.mh_gspspec.notna())) add("RUWE < 1.4 only", GaiaSample(b, L, tw, cut=b.ruwe < 1.4)) add("d < 150 pc", GaiaSample(b, L, tw, cut=b.d < 150, dmax=150.0)) add("d < 100 pc", GaiaSample(b, L, tw, cut=b.d < 100, dmax=100.0)) add("p_sg > 0.5 only", GaiaSample(b, L, tw, cut=b.p_sg > 0.5)) for pot in ("light", "heavy", "MW14"): add(f"{pot} potential", GaiaSample(b, L, tw, pot)) add("tracer ages > 3 Gyr only", GaiaSample(b, L, tw, tau_min=3.0), tau_min=3.0) # systematics that require recomputing ages and tracer weights for lab, kw in [("E(BP-RP)=0.02, A_G=0.04 de-reddening", dict(dc=0.02, dM=0.04)), ("[M/H] shifted by +0.1 dex", dict(mh_shift=0.1)), ("[M/H] shifted by -0.1 dex", dict(mh_shift=-0.1))]: b2, L2, _ = prepare(grid=grid, **kw) add(lab, GaiaSample(b2, L2, tw)) # isochrone colour zero-point: shift the models by +/-0.02 mag in colour for dc in (0.02, -0.02): tw2 = S.TracerWeights(sg_box_func(dc=dc)) b2, L2, _ = prepare(grid=grid, dc=-dc) # model + dc == data - dc add(f"isochrone colours {dc:+.2f} mag", GaiaSample(b2, L2, tw2)) # exclude the part of the colour-magnitude diagram that unresolved # equal-mass main-sequence binaries can reach (data and model alike) wb = dict(WINDOW, binary_safe=True) b2, L2, _ = prepare(grid=grid, window=wb) add("binary-safe window (above the twin sequence)", GaiaSample(b2, L2, S.TracerWeights(sg_box_func(window=wb)))) pd.DataFrame(rows).to_csv(C.DATA + "/gaia_sg_results.csv", index=False) out = pd.DataFrame({"age_lo": C.AGE_EDGES[:-1], "age_hi": C.AGE_EDGES[1:]}) for k in ["formed: all [M/H]", "formed: [M/H]>-0.5", "formed: [M/H]>-0.3", "living MS: all [M/H]", "no column weights (mid-plane sphere)", "ages on the X&R22 scale: formed, all [M/H]", "ages on the Nataf+24 scale: formed, all [M/H]"]: if k in phis: out[k] = phis[k] out.to_csv(C.DATA + "/gaia_sg_age_distribution.csv", index=False) return rows, phis, sm, info if __name__ == "__main__": main()