#!/usr/bin/env python3 """Pillar 1a: column SFH and head-start distribution from the Xiang & Rix (2022) LAMOST subgiants (ages precise to ~7%). Selection handling ------------------ * distance shell 0.3 < d < 0.8 kpc: inside it the 1/G-weighted age mix does not depend on distance (bright young subgiants saturate LAMOST closer in, faint old ones drop out farther away; see the shell test written to data/xr_shell_test.csv); * |R - R0| < 0.5 kpc (solar annulus) and |Z| < z0 (default 200 pc); * tracer weights I0/G(tau, [M/H]) for the subgiant box (MIST isochrones); * vertical-orbit weights from the catalogue's vertical action J_z, converted to energies in the potential the actions were computed in (MWPotential2014). Revision (after the referee report) ----------------------------------- * cut-flow table with the count after every selection cut (xr_cutflow.csv), including the stars that end up with zero vertical weight; * alternative vertical weighting from STAR-COUNT scale heights (Mazzi+ 2024, Table 3) instead of the dynamical weights (modes "sc:exp", "sc:sech2", "sc:exp_thick"); * a direct test of both vertical models against the age mix LAMOST observes in layers above the slab (xr_layer_test.csv); * metallicity floors on total metallicity [M/H] as well as on [Fe/H]; * a systematics table (xr_systematics.csv): age scale, young fraction, vertical model, column cap. """ import numpy as np import pandas as pd import common as C import sgweights as S import tracer as T R0 = 8.178 # kpc, as adopted by Xiang & Rix SHELL = (0.3, 0.8) Z0 = 200.0 # pc DR = 0.5 # kpc G_FLOOR = 0.01 # floor on G/I0 (caps the tracer weight at 100) def load(): x = pd.read_csv(C.RAW + "xiangrix2022.csv").rename( columns={"[Fe/H]": "feh", "[alpha/Fe]": "afe"}) x["R"] = np.hypot(x.X, x.Y) x["d"] = np.sqrt((x.X - R0) ** 2 + x.Y ** 2 + x.Z ** 2) x["mh"] = C.mh_from_feh_alpha(x.feh, x.afe) return x def select(x, shell=SHELL, z0=Z0, dR=DR, tau_min=T.TAU_MIN, relerr=None): m = ((x.Age >= tau_min) & (x.d > shell[0]) & (x.d < shell[1]) & (np.abs(x.R - R0) < dR) & (np.abs(x.Z) * 1000 < z0) & (x.JZ > 0)) if relerr is not None: m &= (x.e_Age / x.Age) < relerr return x[m].reset_index(drop=True) class Sample: def __init__(self, s, tw, pot_name="MW14", z0=Z0, mode="vT", zcap=None): """zcap [pc]: if given, count only the part of each orbit below |z| = zcap (column integrated to zcap instead of infinity).""" self.s, self.tw, self.z0, self.mode = s, tw, z0, mode self.zcap = zcap self.pot = C.get_potential(pot_name) self.age = s.Age.values.astype(float) self.err = np.maximum(s.e_Age.values.astype(float), 0.05) self.zpc = s.Z.values * 1000.0 self.I0 = S.imf_integral(0.8, 1.2, tw_slope(tw)) self.g = np.maximum(tw.g(self.age, s.mh.values) / self.I0, G_FLOOR) # J_z was computed in MWPotential2014: recover (z, v_z) there, then # evaluate the energy in the potential used for the weights. mw = C.get_potential("MW14") E_mw = mw.E_of_J(s.JZ.values * 1000.0) vz2 = np.maximum(2 * (E_mw - mw.phi(self.zpc)), 0.0) self.E = 0.5 * vz2 + self.pot.phi(self.zpc) self.fcap = (1.0 if zcap is None else self.pot.slab_fraction(self.E, zcap)) self.alive = S.alive_fraction(tw, self.age, s.mh.values) # same for the narrower mass range of the GCS comparison sample self.alive09 = S.alive_fraction(tw, self.age, s.mh.values, m_lo=0.9) def weights(self, idx=None, living=False): """living: False (stars ever formed), True (0.8-1.2 Msun stars still on the main sequence) or "0.9" (0.9-1.2 Msun, as in the GCS census).""" idx = np.arange(len(self.age)) if idx is None else idx if self.mode.startswith("sc:"): # star-count scale heights wd = starcount_weights(self.age[idx], self.z0, self.mode[3:]) else: wd = T.vertical_weights(self.pot, self.E[idx], self.zpc[idx], self.z0, self.mode) w = wd / self.g[idx] if self.zcap is not None: w = w * self.fcap[idx] if living == "0.9": return w * self.alive09[idx] return w * self.alive[idx] if living else w def tw_slope(tw): return getattr(tw, "slope", C.IMF_SLOPE) def shell_test(x, tw): """1/G-weighted age mix versus distance (no vertical weights): tests that the LAMOST magnitude limits do not bias the age mix inside the shell.""" I0 = S.imf_integral(0.8, 1.2) rows = [] for lo, hi in [(0.15, 0.3), (0.3, 0.45), (0.45, 0.6), (0.6, 0.8), (0.8, 1.0), (1.0, 1.3), (1.3, 1.6), (1.6, 2.0)]: s = select(x, (lo, hi), 300.0) w = 1 / np.maximum(tw.g(s.Age.values, s.mh.values) / I0, G_FLOOR) a = s.Age.values rows.append(dict(d_lo=lo, d_hi=hi, n=len(s), f_gt_sun=np.sum(w * (a > C.T_SUN)) / w.sum(), f_gt_8=np.sum(w * (a > 8)) / w.sum())) return pd.DataFrame(rows) # Mazzi+ 2024 (MNRAS 527, 583), Table 3: exponential scale heights of # subgiant star counts between isochrones of the given log(age/yr). MAZZI24_HZ = [(9.3, 9.4, 211.4, 16.8), (9.4, 9.5, 264.2, 19.5), (9.5, 9.6, 296.3, 16.8), (9.6, 9.7, 318.7, 19.4), (9.7, 9.8, 388.1, 31.6), (9.8, 9.9, 445.3, 36.1), (9.9, 10.0, 542.8, 63.6), (10.0, 10.1, 367.7, 47.0)] # same table, sech^2 fits (h_z of sech^2(z / h_z)) MAZZI24_HZ_SECH2 = [278.6, 335.3, 376.8, 408.7, 463.3, 470.2, 576.5, 497.9] _HZ_EDGES = np.array([m[0] for m in MAZZI24_HZ] + [MAZZI24_HZ[-1][1]]) # canonical thick-disc scale height (Juric+ 2008; Bland-Hawthorn & Gerhard # 2016: 900 +/- 180 pc), used in the "exp_thick" variant for the two oldest # bins, where Mazzi+24's own single-exponential fit (368 pc above 10 Gyr) # is lower than for the 8-10 Gyr bin (543 pc) H_THICK = (600.0, 900.0) def starcount_height(age, kind="exp"): """Scale height [pc] by age from Mazzi+24 star counts. kind: "exp" | "sech2" | "exp_thick" (exp with 600 pc at 8-10 Gyr and 900 pc above 10 Gyr).""" la = np.log10(np.maximum(np.asarray(age, float), 0.1)) + 9.0 k = np.clip(np.digitize(la, _HZ_EDGES) - 1, 0, len(MAZZI24_HZ) - 1) if kind == "sech2": return np.asarray(MAZZI24_HZ_SECH2)[k] h = np.array([m[2] for m in MAZZI24_HZ])[k] if kind == "exp_thick": h = np.where(la >= 10.0, H_THICK[1], np.where(la >= 9.9, H_THICK[0], h)) return h def _profile_integral(h, a, b, kind): """int_a^b of the (un-normalised) vertical profile, a, b >= 0 [pc].""" if kind == "sech2": return h * (np.tanh(b / h) - np.tanh(a / h)) return h * (np.exp(-a / h) - np.exp(-b / h)) def starcount_weights(age, z0, kind="exp"): """Column / (count inside |z| < z0) for a population whose vertical profile is the star-count fit: no kinematics, and it assumes that the survey samples the slab uniformly in z.""" h = starcount_height(age, kind) return h / _profile_integral(h, 0.0, z0, kind) AGE_BINS = [(2, 3), (3, 4.567), (4.567, 6), (6, 8), (8, 10), (10, 20)] def cutflow(x, tw): """Number of stars left after every selection cut, with the reason for the cut, plus diagnostics of the stars that receive zero weight.""" steps = [ ("Xiang & Rix (2022) catalogue", np.ones(len(x), bool), "all rows of VizieR J/other/Nat/603.599"), ("age >= 2 Gyr", x.Age >= T.TAU_MIN, "no 0.8-1.2 Msun star has reached the subgiant branch earlier"), ("0.3 < d < 0.8 kpc", (x.d > SHELL[0]) & (x.d < SHELL[1]), "age mix independent of distance inside this shell (shell test)"), ("|R - R0| < 0.5 kpc", np.abs(x.R - R0) < DR, "solar annulus"), ("|Z| < 200 pc", np.abs(x.Z) * 1000 < Z0, "slab from which the column is reconstructed"), ("J_z > 0", x.JZ > 0, "vertical action needed for the orbit weights"), ] m = np.ones(len(x), bool) rows = [] for name, cut, why in steps: m &= np.asarray(cut) rows.append(dict(step=name, n=int(m.sum()), reason=why)) s = select(x) sm = Sample(s, tw) wd = T.vertical_weights(sm.pot, sm.E, sm.zpc, sm.z0, sm.mode) zero = wd == 0 rows.append(dict( step="of which: zero vertical weight", n=int(zero.sum()), reason="1-D energy implied by J_z lies below the potential at the " "observed height (v_z clipped to 0); " f"{100 * (sm.age[zero] > C.T_SUN).mean():.1f}% of them are " f"older than the Sun against " f"{100 * (sm.age > C.T_SUN).mean():.1f}% of all; median J_z " f"{np.median(s.JZ.values[zero]):.2f} against " f"{np.median(s.JZ.values):.2f} kpc km/s")) rows.append(dict(step="of which: tracer weight at the cap (G/I0 = 0.01)", n=int((sm.g <= G_FLOOR + 1e-12).sum()), reason="weight capped at 100")) rows.append(dict(step="of which: catalogue age > 13.8 Gyr", n=int((sm.age > C.T_MAX).sum()), reason="kept; counted at 13.8 Gyr in means")) return pd.DataFrame(rows) def bright_limit_test(x, tw): """Distance distribution of the tracers inside the shell, by age. If the brightest (youngest) subgiants are lost at the near side (LAMOST bright limit) they sit farther away on average than the older ones. Assuming the far part (d > 0.55 kpc) of every age bin is complete, the near part is sampled at a rate c relative to the tracers older than 3 Gyr and the completeness of the bin is 1 / ((1 - f) / c + f).""" s = select(x) a, d = s.Age.values, s.d.values ref = a >= 3.0 f_ref = float((d[ref] > 0.55).mean()) rows = [] for lo, hi in AGE_BINS: m = (a >= lo) & (a < hi) f = float((d[m] > 0.55).mean()) c = ((1 - f) / f) / ((1 - f_ref) / f_ref) rows.append(dict(age_lo=lo, age_hi=hi, n=int(m.sum()), median_d_kpc=float(np.median(d[m])), frac_beyond_0p55=f, frac_beyond_0p55_err=float( np.sqrt(f * (1 - f) / m.sum())), near_side_sampling=float(c), implied_completeness=float(1 / ((1 - f) / c + f)))) return pd.DataFrame(rows) def layer_test(x, tw, layers=((200.0, 350.0), (350.0, 500.0), (500.0, 800.0)), z0=Z0): """Which vertical model is right? A test that needs no star-count fit. The stars in the slab |Z| < z0, given their vertical orbits, predict the age mix at any other height. LAMOST also observes subgiants above the slab, inside the same distance shell (so with the same magnitude range). Its sampling efficiency varies with height but not with age, so the 1/G-weighted AGE MIX in a layer can be compared with the prediction of (i) the dynamical weights and (ii) the star-count scale heights. Returns one row per layer and model with the fraction of tracers older than the Sun and older than 8 Gyr.""" slab = select(x, SHELL, z0) sm = Sample(slab, tw) w_dyn = sm.weights() inv_g = 1.0 / sm.g rows = [] def frac(w, a, cut): return float(np.sum(w * (a > cut)) / w.sum()) def ferr(w, a, cut): f = frac(w, a, cut) return float(np.sqrt(np.sum(w ** 2 * ((a > cut) - f) ** 2)) / w.sum()) rows.append(dict(layer_lo=0.0, layer_hi=z0, model="observed", n=len(slab), f_older=frac(inv_g, sm.age, C.T_SUN), f_older_err=ferr(inv_g, sm.age, C.T_SUN), f_gt8=frac(inv_g, sm.age, 8.0), f_gt8_err=ferr(inv_g, sm.age, 8.0))) for a, b in layers: o = select(x, SHELL, b) o = o[np.abs(o.Z.values) * 1000 >= a].reset_index(drop=True) go = np.maximum(tw.g(o.Age.values, o.mh.values) / sm.I0, G_FLOOR) ao = o.Age.values.astype(float) rows.append(dict(layer_lo=a, layer_hi=b, model="observed", n=len(o), f_older=frac(1 / go, ao, C.T_SUN), f_older_err=ferr(1 / go, ao, C.T_SUN), f_gt8=frac(1 / go, ao, 8.0), f_gt8_err=ferr(1 / go, ao, 8.0))) fl = sm.pot.frac_weighted( sm.E, lambda z: ((np.abs(z) >= a) & (np.abs(z) < b)).astype(float)) preds = {"dynamical weights": w_dyn * fl} for kind in ("exp", "sech2", "exp_thick"): h = starcount_height(sm.age, kind) preds[f"star counts ({kind})"] = inv_g * ( _profile_integral(h, a, b, kind) / _profile_integral(h, 0.0, z0, kind)) for name, w in preds.items(): rows.append(dict(layer_lo=a, layer_hi=b, model=name, n=len(slab), f_older=frac(w, sm.age, C.T_SUN), f_older_err=ferr(w, sm.age, C.T_SUN), f_gt8=frac(w, sm.age, 8.0), f_gt8_err=ferr(w, sm.age, 8.0))) return pd.DataFrame(rows) def scale_height_check(x, tw, z0=100.0, pot="MW14"): """Validation of the dynamical column weights against star counts. For stars in the slab |z| < z0 the column density of a population is n0 * 2 z0 * <1/f_slab>, so its effective scale height (column / 2 n0, = h for an exponential profile) is z0 * <1/f_slab>. Compared with the scale heights Mazzi+24 fitted to Gaia subgiant counts versus |z|. Also returns the high-alpha fraction before/after column weighting.""" s = select(x, SHELL, z0) sm = Sample(s, tw, pot, z0, "slab") w = 1.0 / sm.pot.slab_fraction(sm.E, z0) la = np.log10(sm.age) + 9.0 rows = [] for lo, hi, h, eh in MAZZI24_HZ: m = (la >= lo) & (la < hi) if m.sum() < 20: continue rows.append(dict(logage_lo=lo, logage_hi=hi, n=int(m.sum()), h_eff_pc=z0 * w[m].mean(), h_eff_err_pc=z0 * w[m].std() / np.sqrt(m.sum()), h_mazzi24_exp_pc=h, h_mazzi24_err_pc=eh)) df = pd.DataFrame(rows) wt = w / sm.g hi_a = (s.afe.values > 0.16) df.attrs["high_alpha_midplane"] = float(np.sum(hi_a / sm.g) / np.sum(1 / sm.g)) df.attrs["high_alpha_column"] = float(np.sum(wt * hi_a) / wt.sum()) return df def one(x, tw, name, n_boot=0, pot="MW14", z0=Z0, mode="vT", shell=SHELL, dR=DR, tau_min=T.TAU_MIN, relerr=None, feh_floor=None, living=False, deconvolve=True, rmask=None, phi_young=T.PHI_YOUNG[0], zcap=None, age_scale=1.0, mh_floor=None, keep_func=None, wfunc=None): """feh_floor / mh_floor: metallicity floor on [Fe/H] or on the total metallicity [M/H] (Salaris+ 1993, from [Fe/H] and [alpha/Fe]); keep_func(s) -> bool mask (any other sub-population); wfunc(s) -> extra multiplicative weight (completeness corrections).""" if age_scale != 1.0: # systematic stretch of the stellar age scale # (stars pushed beyond the age grid are kept: they pile up in the # last bins and still count as older than the Sun) x = x.assign(Age=x.Age * age_scale, e_Age=x.e_Age * age_scale) s = select(x, shell, z0, dR, tau_min, relerr) if rmask is not None: s = s[rmask(s)].reset_index(drop=True) sm = Sample(s, tw, pot, z0, mode, zcap) keep = None if feh_floor is not None: keep = s.feh.values > feh_floor if mh_floor is not None: keep = s.mh.values > mh_floor if keep_func is not None: keep = np.asarray(keep_func(s)) extra = np.ones(len(s)) if wfunc is None else np.asarray(wfunc(s), float) w = sm.weights(living=living) * extra st, phi = T.run_variant(sm.age, sm.err, sm.err, w, keep, phi_young, deconvolve=deconvolve, tau_min=tau_min) row = dict(variant=name, n=len(s) if keep is None else int(keep.sum()), n_eff=float(w.sum() ** 2 / (w ** 2).sum()), **st) hi_a = s.afe.values > 0.16 row["high_alpha_share_tracers"] = float(np.sum(w * hi_a) / w.sum()) if n_boot: b = T.bootstrap_variant(sm.age, sm.err, sm.err, lambda idx: sm.weights(idx, living) * extra[idx], keep, n_boot=n_boot) for k, (p16, p50, p84) in b.items(): row[k + "_lo"], row[k + "_hi"] = p16, p84 return row, phi # Multiplicative age-scale factors applied to the Xiang & Rix ages in the # systematics table. The catalogue factors are the median age ratios # measured for stars in common (age_scales.py; the values below are only the # fall-back if that table has not been written yet). 0.87 and 1.146 are the # factors the referee derived from this study's own Gaia-only photometric # ages (X&R / Gaia-only and its inverse); +/-10% is the generic test of the # first pass. RATIO_DEFAULT = {"FLAME": 0.93, "Nataf": 0.95, "StarHorse": 1.04} SCALE_NAMES = {"FLAME": "ages on the Gaia FLAME scale", "Nataf": "ages on the Nataf+24 scale", "StarHorse": "ages on the StarHorse scale"} FIXED_SCALES = [(0.87, "ages x0.87 (referee: X&R / Gaia-only ages)"), (0.90, "ages x0.90 (generic -10%)"), (1.10, "ages x1.10 (generic +10%)"), (1.146, "ages x1.146 (referee: Gaia-only scale)")] SYS_KEYS = ["f_gt_0", "f_gt_1", "f_gt_3", "f_gt_5", "mean_headstart_older", "median_headstart_older", "mean_offset_all", "high_alpha_share_tracers"] def main(n_boot=100): x = load() tw = S.TracerWeights("xr") print(f"[xr] catalogue rows: {len(x)}") pd.set_option("display.width", 250, "display.max_columns", 30, "display.max_colwidth", 80) cf = cutflow(x, tw) cf.to_csv(C.DATA + "/xr_cutflow.csv", index=False) print(cf[["step", "n"]].to_string(index=False)) shell_test(x, tw).to_csv(C.DATA + "/xr_shell_test.csv", index=False) bl = bright_limit_test(x, tw) bl.to_csv(C.DATA + "/xr_bright_limit.csv", index=False) print("[xr] bright-limit test (distance distribution by age):") print(bl.round(3).to_string(index=False)) hz = scale_height_check(x, tw) hz["h_mazzi24_sech2_pc"] = MAZZI24_HZ_SECH2[:len(hz)] hz.to_csv(C.DATA + "/xr_scale_heights.csv", index=False) print("[xr] effective scale heights vs Mazzi+24 star counts:") print(hz.round(1).to_string(index=False)) print("[xr] high-alpha ([alpha/Fe]>0.16) fraction of formed Sun-like " f"stars: mid-plane {hz.attrs['high_alpha_midplane']:.3f}, " f"column {hz.attrs['high_alpha_column']:.3f}") lt = layer_test(x, tw) lt.to_csv(C.DATA + "/xr_layer_test.csv", index=False) print("[xr] age mix above the slab: observed vs predicted") print(lt.round(3).to_string(index=False)) rows, phis = [], {} def add(name, nb=0, tw_=None, **kw): r, phi = one(x, tw_ or tw, name, nb, **kw) rows.append(r) phis[name] = phi print(f"[xr] {name:46s} N={r['n']:6d} f_older={r['f_gt_0']:.3f} " f"f>1={r['f_gt_1']:.3f} f>3={r['f_gt_3']:.3f} " f"f>5={r['f_gt_5']:.3f} ={r['mean_headstart_older']:.2f}") return r # headline variants (with bootstrap) add("formed: all [Fe/H]", n_boot) add("formed: [Fe/H]>-0.5", n_boot, feh_floor=-0.5) add("formed: [Fe/H]>-0.3", n_boot, feh_floor=-0.3) add("formed: [M/H]>-0.5", n_boot, mh_floor=-0.5) add("formed: [M/H]>-0.3", n_boot, mh_floor=-0.3) add("living MS: all [Fe/H]", n_boot, living=True) add("living MS: [Fe/H]>-0.5", n_boot, living=True, feh_floor=-0.5) add("living MS: [Fe/H]>-0.3", n_boot, living=True, feh_floor=-0.3) # vertical structure: star-count scale heights instead of dynamics add("star-count heights (Mazzi+24 exp)", n_boot, mode="sc:exp") add("star-count heights (Mazzi+24 sech2)", mode="sc:sech2") add("star-count heights (exp, thick disc 600/900 pc)", n_boot, mode="sc:exp_thick") add("star-count heights (exp), [Fe/H]>-0.5", mode="sc:exp", feh_floor=-0.5) add("star-count heights (exp, thick), [Fe/H]>-0.5", mode="sc:exp_thick", feh_floor=-0.5) add("low-alpha stars only ([alpha/Fe] <= 0.16)", keep_func=lambda s: s.afe.values <= 0.16) add("living MS 0.9-1.2 Msun (GCS-like), column", living="0.9") add("living MS 0.9-1.2 Msun (GCS-like), mid-plane", living="0.9", mode="none") # robustness (point estimates) add("no column weights (mid-plane slab)", mode="none") add("no error deconvolution", deconvolve=False) add("slab-fraction weights (uniform S(z))", mode="slab") add("column to |z| < 0.8 kpc only (height verified by the layer test)", zcap=800.0) add("column to |z| < 1.3 kpc only", zcap=1300.0) add("column to |z| < 2 kpc only", zcap=2000.0) add("z0 = 150 pc", z0=150.0) add("z0 = 300 pc", z0=300.0) add("fiducial KG potential", pot="fiducial") add("light-disc potential", pot="light") add("heavy-disc potential", pot="heavy") add("shell 0.3-0.6 kpc", shell=(0.3, 0.6)) add("shell 0.45-0.8 kpc", shell=(0.45, 0.8)) add("shell 0.3-1.0 kpc", shell=(0.3, 1.0)) add("|R-R0| < 0.25 kpc", dR=0.25) add("inner half (R < R0)", rmask=lambda s: s.R < R0) add("outer half (R > R0)", rmask=lambda s: s.R >= R0) add("age error < 15%", relerr=0.15) add("tracer ages > 2.5 Gyr only", tau_min=2.5) import age_scales ratios = age_scales.measured_ratios(RATIO_DEFAULT) for key, lab in SCALE_NAMES.items(): add(lab, age_scale=ratios[key])["age_scale"] = ratios[key] for sc, lab in FIXED_SCALES: add(lab, age_scale=sc)["age_scale"] = sc add("phi_young = 0.06", phi_young=0.06) add("phi_young = 0.13 (tracers' own flat-rate value)", phi_young=0.13) add("phi_young = 0.15", phi_young=0.15) # bright-limit loss of the youngest tracers: up-weight them by the # inverse of the completeness implied by their distance distribution comp = {(r.age_lo, r.age_hi): r.implied_completeness for r in bl.itertuples()} c23 = comp[(2, 3)] add("2-3 Gyr tracers up-weighted for bright-limit loss", wfunc=lambda s: np.where(s.Age.values < 3.0, 1.0 / c23, 1.0)) for lab, kw in [("isochrones +50 K", dict(dteff=50.0)), ("isochrones -50 K", dict(dteff=-50.0)), ("isochrones +0.1 mag in M_K", dict(dmk=0.1)), ("isochrones -0.1 mag in M_K", dict(dmk=-0.1))]: add(lab, tw_=S.TracerWeights("xr", **kw)) for slope in (2.0, 2.7): tws = S.TracerWeights("xr", slope=slope); tws.slope = slope add(f"IMF slope {slope}", tw_=tws) res = pd.DataFrame(rows) res.to_csv(C.DATA + "/xr_results.csv", index=False) # ---- systematics table for the headline number r = res.set_index("variant") base = r.loc["formed: all [Fe/H]"] groups = [ ("baseline (dynamical weights, X&R ages, phi_young 0.10)", ["formed: all [Fe/H]"]), ("statistical + young-fraction prior (16-84% bootstrap)", None), ("age scale: other published catalogues (FLAME x%.2f, Nataf+24 " "x%.2f, StarHorse x%.2f)" % (ratios["FLAME"], ratios["Nataf"], ratios["StarHorse"]), list(SCALE_NAMES.values())), ("age scale: generic +/-10%", ["ages x0.90 (generic -10%)", "ages x1.10 (generic +10%)"]), ("age scale: this study's Gaia-only ages (x0.87, x1.146; referee)", ["ages x0.87 (referee: X&R / Gaia-only ages)", "ages x1.146 (referee: Gaia-only scale)"]), ("vertical model: star-count scale heights", ["star-count heights (Mazzi+24 exp)", "star-count heights (Mazzi+24 sech2)", "star-count heights (exp, thick disc 600/900 pc)"]), ("vertical model: column counted to 0.8 / 1.3 / 2 kpc", ["column to |z| < 0.8 kpc only (height verified by the layer test)", "column to |z| < 1.3 kpc only", "column to |z| < 2 kpc only"]), ("vertical model: potential (light / heavy / KG fiducial)", ["light-disc potential", "heavy-disc potential", "fiducial KG potential"]), ("young fraction 0.06 / 0.13 / 0.15", ["phi_young = 0.06", "phi_young = 0.13 (tracers' own flat-rate value)", "phi_young = 0.15"]), ("bright-limit loss of 2-3 Gyr tracers", ["2-3 Gyr tracers up-weighted for bright-limit loss"]), ("selection geometry (shell, slab, annulus)", ["shell 0.3-0.6 kpc", "shell 0.45-0.8 kpc", "shell 0.3-1.0 kpc", "z0 = 150 pc", "z0 = 300 pc", "|R-R0| < 0.25 kpc", "inner half (R < R0)", "outer half (R > R0)"]), ("tracer weights (isochrone shifts, IMF slope)", ["isochrones +50 K", "isochrones -50 K", "isochrones +0.1 mag in M_K", "isochrones -0.1 mag in M_K", "IMF slope 2.0", "IMF slope 2.7"]), ] srows = [] for label, names in groups: if names is None: srows.append(dict(source=label, f_min=base.f_gt_0_lo, f_max=base.f_gt_0_hi, mean_hs_min=base.mean_headstart_older_lo, mean_hs_max=base.mean_headstart_older_hi, variants="bootstrap")) continue f = r.loc[names, "f_gt_0"] h = r.loc[names, "mean_headstart_older"] srows.append(dict(source=label, f_min=f.min(), f_max=f.max(), mean_hs_min=h.min(), mean_hs_max=h.max(), variants="; ".join(f"{n}: {v:.3f}" for n, v in f.items()))) sysdf = pd.DataFrame(srows) sysdf["delta_min"] = sysdf.f_min - base.f_gt_0 sysdf["delta_max"] = sysdf.f_max - base.f_gt_0 sysdf.to_csv(C.DATA + "/xr_systematics.csv", index=False) print("[xr] systematics of the fraction older than the Sun:") print(sysdf[["source", "f_min", "f_max", "mean_hs_min", "mean_hs_max"]] .round(3).to_string(index=False)) out = pd.DataFrame({"age_lo": C.AGE_EDGES[:-1], "age_hi": C.AGE_EDGES[1:]}) for k in (["formed: all [Fe/H]", "formed: [Fe/H]>-0.5", "formed: [Fe/H]>-0.3", "formed: [M/H]>-0.5", "formed: [M/H]>-0.3", "living MS: all [Fe/H]", "living MS: [Fe/H]>-0.5", "living MS: [Fe/H]>-0.3", "no column weights (mid-plane slab)", "no error deconvolution", "star-count heights (Mazzi+24 exp)", "star-count heights (exp, thick disc 600/900 pc)", "star-count heights (exp), [Fe/H]>-0.5", "star-count heights (exp, thick), [Fe/H]>-0.5", "low-alpha stars only ([alpha/Fe] <= 0.16)"] + list(SCALE_NAMES.values())): out[k] = phis[k] out.to_csv(C.DATA + "/xr_age_distribution.csv", index=False) return rows, phis if __name__ == "__main__": main()