"""Column star-formation history from subgiant 'tracer' stars. Each subgiant i with age tau_i stands for w_i = w_dyn,i * I0 / G(tau_i, [M/H]_i) Sun-like (0.8-1.2 Msun) stars *formed* per unit area of the Galactic plane at look-back time tau_i: * I0/G converts one subgiant in the selection box into the number of Sun-like stars born at the same epoch (IMF + stellar-evolution speed; sgweights.py); * w_dyn converts a star seen near the plane into its contribution to the vertical column (orbit-time weighting; see vertical_weights()). Subgiants only trace ages above TAU_MIN (2 Gyr). The fraction of all stars formed in the last TAU_MIN Gyr, PHI_YOUNG, is taken from published CMD-fitting SFHs of the solar cylinder (Mor+19, Mazzi+24, del Alcazar-Julia+25): 0.10 +/- 0.04 (mean and scatter of six solutions; see published_sfh.py). """ import numpy as np import common as C TAU_MIN = 2.0 PHI_YOUNG = (0.10, 0.04) # mean, 1-sigma (truncated to 0.02-0.25) N_ITER = 100 def vertical_weights(pot, E, z_pc, z0, mode="vT", nbin=15): """Column weights for stars observed inside the slab |z| < z0. mode="vT" (default). A star of vertical energy E is found at height z with probability density p(z|E) = 2 / (T(E) |v_z(z)|). If the survey samples height z with an unknown efficiency S(z) that does not depend on age, then sum_i 1/p(z_i|E_i) = sum_i |v_z,i| T_i / 2 over the observed stars is proportional to the column density of the population times int_{|z| z0 and the weights of the cold orbits are corrected by Q = int_0^zmax S / int_0^z0 S. mode="slab". Assumes S(z) is uniform inside the slab: w = 1 / (fraction of the orbital period spent at |z| < z0). """ E = np.asarray(E, float) z = np.abs(np.asarray(z_pc, float)) if mode == "none": return np.ones(len(E)) if mode == "slab": return 1.0 / pot.slab_fraction(E, z0) zmax = pot.zmax_of_E(E) v = np.sqrt(np.maximum(2 * (E - pot.phi(z)), 0.0)) w0 = v * pot.T_of_E(E) / (4 * z0) hot = zmax > z0 edges = np.linspace(0, z0, nbin + 1) Sz, _ = np.histogram(z[hot], bins=edges, weights=w0[hot]) cum = np.concatenate([[0.0], np.cumsum(Sz)]) / max(Sz.sum(), 1e-30) Q = np.interp(np.minimum(zmax, z0), edges, cum) return w0 / np.maximum(Q, 0.05) def combine_with_young(phi_old, phi_young, tau_min=TAU_MIN, edges=C.AGE_EDGES, old_keep=1.0): """Full age distribution: tracer-based shape above tau_min, a uniform rate below it carrying a fraction phi_young of all stars formed. old_keep < 1 rescales the old part (metallicity floors remove part of the old population but none of the young one).""" lo, hi = edges[:-1], edges[1:] young_shape = np.clip(np.minimum(hi, tau_min) - lo, 0, None) young_shape = young_shape / young_shape.sum() old = np.where(lo >= tau_min - 1e-9, phi_old, 0.0) old = old / old.sum() phi = phi_young * young_shape + (1 - phi_young) * old_keep * old return phi / phi.sum() def tracer_sfh(age, s_lo, s_hi, w, n_iter=N_ITER, kind="splitnormal", tau_min=TAU_MIN): """Error-deconvolved, weighted age histogram of the tracers.""" L = C.likelihood_matrix(age, s_lo, s_hi, kind=kind) # tracers are only selected above tau_min: forbid younger true ages L = L * (C.AGE_MID[None, :] >= tau_min) L = L / np.maximum(L.sum(axis=1, keepdims=True), 1e-300) return C.em_deconvolve(L, w, n_iter=n_iter) def run_variant(age, s_lo, s_hi, w, keep=None, phi_young=PHI_YOUNG[0], n_iter=N_ITER, deconvolve=True, tau_min=TAU_MIN, rescale_young=True): """Head-start statistics for one weighting. keep : boolean mask (e.g. metallicity floor) applied to the tracers; the young (< tau_min) population is assumed to pass any floor. phi_young : fraction of all stars formed in the last TAU_MIN (2) Gyr; if tau_min differs from TAU_MIN it is rescaled assuming a constant recent star-formation rate, unless rescale_young is False, in which case phi_young is the fraction formed in the last tau_min Gyr. The young stars are spread uniformly over 0..tau_min. """ if rescale_young: phi_young = phi_young * tau_min / TAU_MIN keep = np.ones(len(age), bool) if keep is None else np.asarray(keep) old_keep = w[keep].sum() / w.sum() if deconvolve: phi_old = tracer_sfh(age[keep], s_lo[keep], s_hi[keep], w[keep], n_iter, tau_min=tau_min) else: phi_old = C.weighted_hist(age[keep], w[keep]) phi = combine_with_young(phi_old, phi_young, tau_min, old_keep=old_keep) st = C.headstart_stats(phi) st["f_gt8Gyr"] = C.frac_older(phi, 8.0 - C.T_SUN) old_only = np.where(C.AGE_EDGES[:-1] >= tau_min - 1e-9, phi_old, 0) st["f_old_given_tracer"] = C.frac_older(old_only, 0.0) # young fraction implied if the rate at tau < tau_min equals the tracer # rate at tau_min..2 tau_min (self-consistency check of PHI_YOUNG) lo = C.AGE_EDGES[:-1] q = old_only[(lo >= tau_min - 1e-9) & (lo < 2 * tau_min - 1e-9)].sum() \ / old_only.sum() * old_keep st["phi_young_if_flat"] = q / (1 + q) * TAU_MIN / tau_min st["phi_young_used"] = phi_young # fraction of all stars formed that pass the `keep` cut (metallicity # floor), under the assumption that every young star passes it st["frac_pass"] = phi_young + (1 - phi_young) * old_keep return st, phi def bootstrap_variant(age, s_lo, s_hi, weight_func, keep=None, n_boot=100, seed=11, phi_young=PHI_YOUNG, n_iter=60, tau_min=TAU_MIN, rescale_young=True): """Bootstrap over stars (weights recomputed by weight_func(idx)) and over the young-fraction prior. Returns dict key -> (p16, p50, p84).""" rng = np.random.default_rng(seed) n = len(age) keep = np.ones(n, bool) if keep is None else np.asarray(keep) rows = [] for _ in range(n_boot): idx = rng.integers(0, n, n) py = float(np.clip(rng.normal(*phi_young), 0.02, 0.25)) st, _ = run_variant(age[idx], s_lo[idx], s_hi[idx], weight_func(idx), keep[idx], py, n_iter=n_iter, tau_min=tau_min, rescale_young=rescale_young) rows.append(st) keys = rows[0].keys() return {k: tuple(np.percentile([r[k] for r in rows], [16, 50, 84])) for k in keys}