"""Effective samples, vetted candidate counts, upper limits, robustness variants, survivor table and the conversion to numbers of beacons in the survey volume and in the Galaxy.""" import csv, json, os, time, urllib.request, urllib.parse import numpy as np import pandas as pd from scipy.stats import norm from common import * import s03_candidates as s3 import s04_vetting as s4 import s05_physics as s5 MAIN = s4.MAIN SALT = lambda T: T + s3.MARGIN # threshold on the observed excess -> salt level constrained MW_STELLAR_MASS = 5.0e10 # Msun (Bland-Hawthorn & Gerhard 2016, ARA&A 54, 529: 5 +/- 1 e10) RHO_STAR_LOCAL = 0.043 # Msun/pc^3 in stars and remnants (McKee, Parravano & Hollenbach 2015, ApJ 814, 13) SNR_HI = 50.0 # S/N cut of the high-S/N variant (CCD2, per pixel) def gcns_density(): out = RAW + 'gcns_phot.csv' if not (os.path.exists(out) and os.path.getsize(out) > 1000): q = 'SELECT Plx, "Gmag", BPmag, RPmag FROM "J/A+A/649/A6/table1c" WHERE Plx>=10' data = urllib.parse.urlencode({'REQUEST': 'doQuery', 'LANG': 'ADQL', 'FORMAT': 'csv', 'QUERY': q}).encode() with urllib.request.urlopen(urllib.request.Request('https://tapvizier.cds.unistra.fr/TAPVizieR/tap/sync', data=data), timeout=600) as r: open(out, 'wb').write(r.read()) g = pd.read_csv(out) MG = g.Gmag + 5 * np.log10(g.Plx / 100.0); col = g.BPmag - g.RPmag sel = (col >= 0.58) & (col <= 1.33) & (MG >= 2.3) & (MG <= 7.5) & ~((col > 0.95) & (MG < 3.5)) n = int(sel.sum()) return n, n / (4.0 / 3.0 * np.pi * 100.0 ** 3), int(len(g)) def simbad_cone(ra, dec, radius_arcsec=5.0): """Nearest SIMBAD object within the radius: main_id, otype, sp_type, nbref (cached on disk).""" cache_p = RAW + 'simbad_cache.json' cache = json.load(open(cache_p)) if os.path.exists(cache_p) else {} k = '%.5f_%.5f' % (ra, dec) if k in cache: return cache[k] q = ("SELECT TOP 1 main_id, otype, sp_type, nbref, DISTANCE(POINT('ICRS', ra, dec), POINT('ICRS', %.6f, %.6f)) AS d " "FROM basic WHERE CONTAINS(POINT('ICRS', ra, dec), CIRCLE('ICRS', %.6f, %.6f, %.6f))=1 ORDER BY d" % (ra, dec, ra, dec, radius_arcsec / 3600.0)) data = urllib.parse.urlencode({'REQUEST': 'doQuery', 'LANG': 'ADQL', 'FORMAT': 'csv', 'QUERY': q}).encode() res = ['query_failed', '', '', ''] for t in range(3): try: with urllib.request.urlopen(urllib.request.Request('https://simbad.cds.unistra.fr/simbad/sim-tap/sync', data=data), timeout=60) as r: rows = list(csv.reader(r.read().decode().splitlines())) res = [x.strip().strip('"') for x in rows[1][:4]] if len(rows) > 1 else ['not_in_simbad', '', '', ''] cache[k] = res json.dump(cache, open(cache_p, 'w')) break except Exception: time.sleep(2 + 2 * t) return res def twin_calibration(vet): """Efficiency of the twin test from the controls.""" out = {} nul = vet[(vet.kind == 'null_control') & (vet.status == 'ok')] out['null'] = {e: dict(n=int((nul.element == e).sum()), z_median=float(nul[nul.element == e].z_null.median()), z_robust_sigma=float(1.4826 * (nul[nul.element == e].z_null - nul[nul.element == e].z_null.median()).abs().median()), n_z_gt3=int((nul[nul.element == e].z_null > 3).sum()), frac_z_gt3=float((nul[nul.element == e].z_null > 3).mean()), frac_z_gt3_ul95=binom_upper(int((nul[nul.element == e].z_null > 3).sum()), int((nul.element == e).sum()))) for e in MAIN} pos = vet[(vet.kind == 'positive_control') & (vet.status == 'ok') & vet.rho.notna()].copy() def wrong(delta): # a real enhancement is wrongly refuted if the 2-sigma spectral upper bound falls more than # `delta` dex below the catalogue value (delta = 0: a star sitting exactly at the salt level) return pos.xs_hi < pos.x_cat - delta pos['wrong0'] = wrong(0.0); pos['wrongM'] = wrong(s3.MARGIN) pos['sig3'] = pos.rho >= 3 * pos.rho_err out['positive'] = {} for e in MAIN: p = pos[pos.element == e] if len(p): out['positive'][e] = dict(n=int(len(p)), rho_median=float(p.rho.median()), rho_p16=float(p.rho.quantile(0.16)), rho_p84=float(p.rho.quantile(0.84)), frac_rho_3sigma=float(p.sig3.mean()), n_wrongly_refuted_delta0=int(p.wrong0.sum()), n_wrongly_refuted_deltaM=int(p.wrongM.sum()), snr_mod_median=float(p.snr_mod.median())) n = len(pos); k0 = int(pos.wrong0.sum()); kM = int(pos.wrongM.sum()) out['pooled'] = dict(n=n, n_wrong_delta0=k0, n_wrong_deltaM=kM, eff_delta0=1 - k0 / n, eff_deltaM=1 - kM / n, eff_delta0_lower95=1 - binom_upper(k0, n), eff_deltaM_lower95=1 - binom_upper(kM, n)) # DR3-based completeness control: among dwarfs that the independent DR3 pipeline calls strongly # enhanced and that DR4 could detect, how many strong absorbers (z > 5 against twins) does the # DR4 catalogue fail to report? cnts = pd.read_csv(DATA + 'dr3_control_counts.csv').set_index('element') out['dr3_control'] = {} for e in s4.DR3_CTRL: o = {} for kind in ('dr3_only', 'dr3_both', 'null_control'): g = vet[(vet.kind == kind) & (vet.element == e) & (vet.status == 'ok')] z = g.z_null.values o[kind] = dict(n=int(len(g)), z_median=float(np.median(z)) if len(z) else None, frac_z_gt3=float((z > 3).mean()) if len(z) else None, frac_z_gt5=float((z > 5).mean()) if len(z) else None, dew_median_mA=float(g.dew_obs.median()) if len(z) else None) n_only = int(cnts.loc[e, 'n_dr4_no_detection']); n_both = int(cnts.loc[e, 'n_dr4_agrees']) k_only = int(round(o['dr3_only']['frac_z_gt5'] * o['dr3_only']['n'])); k_both = o['dr3_both']['frac_z_gt5'] * o['dr3_both']['n'] miss = n_only * k_only / max(o['dr3_only']['n'], 1); miss_ul = n_only * binom_upper(k_only, o['dr3_only']['n']) real = n_both * k_both / max(o['dr3_both']['n'], 1) o['n_dr3_enhanced_dr4_detectable'] = int(cnts.loc[e, 'n_dr3_enhanced_and_dr4_detectable']) o['n_dr4_no_detection'] = n_only; o['n_dr4_agrees'] = n_both o['strong_absorbers_missed_by_dr4'] = float(miss); o['strong_absorbers_reported_by_dr4'] = float(real) o['catalogue_completeness'] = float(real / (real + miss)) if real + miss > 0 else None o['catalogue_completeness_lower95'] = float(real / (real + miss_ul)) if real + miss_ul > 0 else None out['dr3_control'][e] = o return out def tally(d, c, vet, e, T, mask=None, margin=None, rule='clip'): """Classify the catalogue candidates of element e at threshold T. Returns counts and the list of non-refuted (idx, verdict). Candidates with excess in [T_SUB, T_VET) are only vetted for a random subsample; their surviving number is extrapolated (binomial upper limit).""" cm = s3.candidates(d, c, e, T) if mask is not None: cm = cm & mask idx = np.where(cm)[0] full = vet[(vet.kind == 'candidate') & (vet.element == e)].set_index('idx') sub = vet[(vet.kind == 'candidate_sub') & (vet.element == e)].set_index('idx') cnt = dict(n_catalogue=len(idx), refuted=0, supported=0, ambiguous=0, unverifiable=0, band_total=0, band_vetted=0, band_refuted=0, band_supported=0, band_ambiguous=0, band_unverifiable=0, not_vetted=0) surv = [] for i in idx: if i in full.index: k = s4.classify(full.loc[i], T, s3.MARGIN if margin is None else margin, rule); cnt[k] += 1 if k != 'refuted': surv.append((int(i), k, 'full')) else: cnt['band_total'] += 1 if i in sub.index: k = s4.classify(sub.loc[i], T, s3.MARGIN if margin is None else margin, rule) cnt['band_vetted'] += 1; cnt['band_' + k] += 1 if k != 'refuted': surv.append((int(i), k, 'subsample')) else: cnt['not_vetted'] += 1 return cnt, surv def upper_limit(cnt, neff, eff): """95% upper limit on the salted fraction. Fully vetted part: Poisson limit on the number of non-refuted candidates. Sub-sampled band: band size times the binomial upper limit on the non-refuted fraction (if nothing of the band was vetted, the whole band counts).""" k_full = cnt['supported'] + cnt['ambiguous'] + cnt['unverifiable'] n_ul = poisson_upper(k_full) k_exp = float(k_full) if cnt['band_total'] > 0: if cnt['band_vetted'] > 0: s = cnt['band_vetted'] - cnt['band_refuted'] n_ul += cnt['band_total'] * binom_upper(s, cnt['band_vetted']) k_exp += cnt['band_total'] * s / cnt['band_vetted'] else: n_ul += cnt['band_total']; k_exp += cnt['band_total'] return n_ul / (neff * eff) if neff > 0 else np.nan, k_exp T_NOEXTRAP = 0.75 # sensitivity is not credited beyond this observed-excess threshold in the 'no extrapolation' column def limits_table(d, c, vet, eff, label='baseline', mask=None, margin=None, x90=None, levels=None, rule='clip'): rows = []; surv_all = [] for e in MAIN: for T in (levels or s3.LEVELS): w = s3.eff_weight(d, c, e, T, margin=margin, x90=None if x90 is None else x90[e]) # conservative alternative: the detection model is an extrapolation above the abundances # that natural stars populate, so do not credit any star that would not already have # shown a 1.0 dex salt (T = 0.75) -- a stronger salt is certainly detectable in those w_ne = w if T <= T_NOEXTRAP else s3.eff_weight(d, c, e, T_NOEXTRAP, margin=margin, x90=None if x90 is None else x90[e]) if mask is not None: w = w * mask; w_ne = w_ne * mask neff = float(w.sum()); neff_ne = float(w_ne.sum()) cnt, surv = tally(d, c, vet, e, T, mask=mask, margin=margin, rule=rule) vetted = T >= s4.T_SUB if vetted: ul, k_exp = upper_limit(cnt, neff, eff) ul_ne = upper_limit(cnt, neff_ne, eff)[0] n_supported = cnt['supported'] else: ul = poisson_upper(cnt['n_catalogue']) / neff if neff > 0 else np.nan ul_ne = ul k_exp = cnt['n_catalogue']; n_supported = np.nan rows.append(dict(variant=label, element=e, T_excess=T, salt_level_dex=T + (s3.MARGIN if margin is None else margin), N_eff=neff, N_eff_no_extrapolation=neff_ne, f_ul95_no_extrapolation=ul_ne, efficiency=eff if vetted else 1.0, n_catalogue=cnt['n_catalogue'], n_refuted=cnt['refuted'], n_supported=cnt['supported'], n_ambiguous=cnt['ambiguous'], n_unverifiable=cnt['unverifiable'], band_total=cnt['band_total'], band_vetted=cnt['band_vetted'], band_refuted=cnt['band_refuted'], band_supported=cnt['band_supported'], band_ambiguous=cnt['band_ambiguous'] + cnt['band_unverifiable'], not_vetted=cnt['not_vetted'], n_not_refuted_est=k_exp, f_ul95=ul, f_ul95_catalogue_only=poisson_upper(cnt['n_catalogue']) / neff if neff > 0 else np.nan, vetted=bool(vetted))) if vetted: surv_all += [(e, T, i, k, how) for i, k, how in surv] return rows, surv_all def fission_rows(d, c, vet, eff): """Limit on stars salted with bulk fission waste, using the Nd probe. Fission waste also raises Ba (the anchor), so the observed Nd excess E is smaller than the Nd enrichment; the per-star efficiency is Phi((E - T) / sigma_Nd) instead of Phi(MARGIN / sigma). Revision (referee error 1): the first pass called tally() with the default margin, so a candidate counted as refuted whenever its spectrum excluded an excess of T + 0.25 = 1.0 dex, which a genuinely fission-salted star (E = 0.68-0.93 dex) would itself fail. The refutation threshold must not exceed what a salted star shows. Two choices are tabulated: refute_at = 'E' (referee's fix): refuted if the spectrum excludes an excess >= E; refute_at = 'T' : refuted only if it excludes an excess >= T (the catalogue threshold). The ADOPTED limit refutes at min(T, E): at T when E >= T (more conservative than the referee's fix, because a star scattered below E but above T is still kept), and at E when E < T (there a threshold of T would itself refute the salted stars, so only 'E' is self-consistent). Since a lower refutation threshold can only add non-refuted candidates, this is the larger of the two tabulated limits.""" fx = pd.read_csv(DATA + 'fission_salt_observed_excess.csv') out = [] for _, r in fx.iterrows(): if r['nd_enhancement'] > 2.0: continue for T in (0.5, 0.75): marg = float(r['ex_nd']) - T if marg < -0.10: # expected excess well below the threshold: efficiency under ~25%, not tabulated continue w = s3.eff_weight(d, c, 'nd', T, margin=marg) neff = float(w.sum()) row = dict(fuel=r['nuclide'], nd_enrichment_dex=float(r['nd_enhancement']), observed_ex_nd=float(r['ex_nd']), T_excess=T, N_eff=neff) for lab, mref in (('T', 0.0), ('E', marg)): cnt, surv = tally(d, c, vet, 'nd', T, margin=mref) ul, k_exp = upper_limit(cnt, neff, eff) if neff > 0 else (np.nan, np.nan) row['n_not_refuted_est_refute_at_' + lab] = k_exp; row['f_ul95_refute_at_' + lab] = ul row['refuted_at'] = 'T' if marg >= 0 else 'E' row['f_ul95'] = row['f_ul95_refute_at_T'] if marg >= 0 else row['f_ul95_refute_at_E'] out.append(row) return out M_REF = 0.10 # dex: adopted refutation margin of the multiplicative-definition limits M_CAT = s3.MARGIN def salt_added_rows(d, c, vet, eff, levels=(0.75, 1.0, 1.25, 1.5), m_refs=(0.0, 0.05, 0.10, 0.15, 0.25)): """Limits for the MULTIPLICATIVE definition of a salt (referee minor point (b)). The baseline table limits the fraction of stars whose excess over the natural prediction is >= L. Two readings of 'salted by S dex' differ in how that maps onto the excess: additive : an amount (10^S - 1) x [natural prediction] of the element was added. The excess is then log10(10^S - 1 + 10^-delta), within ~0.02 dex of S for S >= 1 whatever the star's own natural abundance (delta = prediction - natural), because the salt outweighs the natural content 9:1 or more. The baseline table applies directly, and the injection test (s12) measures its vetting efficiency at 0.93-0.97. multiplicative : the star's OWN natural abundance was multiplied by 10^S. The excess over the prediction is then S - delta, and with a good spectrum such a star is refuted at level L = S whenever delta > 0, i.e. about half the time (the referee's point). This function handles that case: the catalogue threshold stays at T = S - 0.25 but a candidate is refuted only if its spectrum excludes an excess >= S - m_ref; the salted star survives both steps if delta <= m_ref, so the per-star efficiency is Phi(m_ref / sigma_core) with the whole core scatter treated as intrinsic (worst case), or 1/0 if the star's own element is securely measured -- what eff_weight(margin=m_ref) computes. m_ref trades efficiency against confusion; M_REF is the adopted value and the sweep is tabulated to show the answer does not hinge on it.""" out = [] for e in MAIN: for S in levels: T = S - M_CAT for m in m_refs: w = s3.eff_weight(d, c, e, T, margin=m); neff = float(w.sum()) w_ne = w if T <= T_NOEXTRAP else s3.eff_weight(d, c, e, T_NOEXTRAP, margin=m) cnt, _ = tally(d, c, vet, e, T, margin=M_CAT - m) ul, k_exp = upper_limit(cnt, neff, eff) out.append(dict(element=e, salt_added_dex=S, T_excess=T, refutation_level_dex=S - m, m_ref=m, adopted=int(abs(m - M_REF) < 1e-9), N_eff_times_scatter_efficiency=neff, n_not_refuted_est=k_exp, f_ul95=ul, f_ul95_no_extrapolation=upper_limit(cnt, float(w_ne.sum()), eff)[0])) return out def consistency_table(vet): """What the twin test says about the catalogue candidates (excess >= T_VET), separating 'excluded as an excess of >= 1 dex' from 'inconsistent with the star's own catalogue value'.""" rows = [] for e in MAIN: g = vet[(vet.kind == 'candidate') & (vet.element == e)] ok = g[(g.status == 'ok') & g.rho.notna()] ver = pd.Series([s4.classify(r, s4.T_VET) for _, r in g.iterrows()], index=g.index) hi = ok.rho + 2 * ok.rho_err rows.append(dict(element=e, n_candidates=len(g), n_tested=len(ok), n_excluded_as_ge_1dex=int((ver == 'refuted').sum()), frac_excluded_as_ge_1dex=float((ver == 'refuted').mean()), n_inconsistent_with_catalogue_2sigma=int((hi < 1.0).sum()), frac_inconsistent_with_catalogue=float((hi < 1.0).mean()), n_less_than_half_the_catalogue_absorption=int((hi < 0.5).sum()), n_consistent_with_catalogue=int(((hi >= 1.0) & (ok.rho - 2 * ok.rho_err <= 1.0)).sum()), n_real_absorption_3sigma=int((ok.rho >= 3 * ok.rho_err).sum()), n_beyond_grid_not_refutable=int((ok.beyond_grid == 1).sum()))) return pd.DataFrame(rows) def main(d=None, c=None, verbose=True): if d is None: d = s3.load(); c = s3.compute(d) vet, nullstat = s4.evaluate() vet.drop(columns=['lines']).to_csv(DATA + 'vetting_twin_test.csv', index=False, float_format='%.4f') nullstat.to_csv(DATA + 'twin_test_null_line_noise.csv', index=False, float_format='%.3f') key = {} teff = d['teff'].astype(float) base = c['base']; clean = base & (d['flag_sp'] == 0) # --- cut flow cut = [('GALAH DR4 allstar rows', len(teff)), ('finite Teff, logg, [Fe/H]', int((np.isfinite(teff) & np.isfinite(d['logg']) & np.isfinite(d['fe_h'])).sum())), ('4500 <= Teff < 6500 K and logg >= 3.5', int(base.sum())), ('... and flag_sp == 0', int(clean.sum()))] for e in MAIN: cut.append(('... and anchor of %s securely detected' % e, int((clean & c['anchor_ok'][e]).sum()))) with open(DATA + 'cut_flow.csv', 'w', newline='') as fh: w = csv.writer(fh); w.writerow(['cut', 'n_stars']); w.writerows(cut) # --- twin-test calibration cal = twin_calibration(vet); key['twin_test'] = cal eff_twin = cal['pooled']['eff_delta0_lower95'] # conservative: 95% lower bound, star exactly at the salt level eff_cat = cal['dr3_control']['nd']['catalogue_completeness'] # measured for Nd, assumed for the other elements eff = eff_twin * eff_cat key['twin_efficiency_used'] = eff_twin; key['catalogue_completeness_used'] = eff_cat; key['efficiency_used'] = eff # --- baseline limits rows, surv_rows = limits_table(d, c, vet, eff) lim = pd.DataFrame(rows) lim.to_csv(DATA + 'limits_by_element.csv', index=False, float_format='%.5g') # --- fission-waste salt fis = pd.DataFrame(fission_rows(d, c, vet, eff)) fis.to_csv(DATA + 'limits_fission_salt.csv', index=False, float_format='%.5g') # --- limits for the multiplicative salt definition, and what the twin test says about the catalogue candidates sadd = pd.DataFrame(salt_added_rows(d, c, vet, eff)) sadd.to_csv(DATA + 'limits_salt_added.csv', index=False, float_format='%.5g') cons = consistency_table(vet) cons.to_csv(DATA + 'candidate_consistency.csv', index=False, float_format='%.4f') key['salt_added_m_ref'] = M_REF key['salt_added'] = sadd[sadd.adopted == 1].to_dict('records') key['candidate_consistency'] = cons.to_dict('records') key['scatter_efficiency'] = {e: float(norm.cdf(M_REF / c['sig'][e])) for e in MAIN} # --- robustness variants (levels T = 0.75 and 1.0 only) var = [] snr = d['snr_px_ccd2'].astype(float) lv = [0.75, 1.0] # T = 0.75 is the 1.0 dex level, T = 1.0 the 1.25 dex level var += limits_table(d, c, vet, eff, 'S/N(CCD2) >= %.0f' % SNR_HI, mask=(snr >= SNR_HI), levels=lv)[0] var += limits_table(d, c, vet, eff, 'Teff < 6250 K', mask=(teff < 6250), levels=lv)[0] var += limits_table(d, c, vet, eff, 'Teff >= 5250 K', mask=(teff >= 5250), levels=lv)[0] var += limits_table(d, c, vet, eff, 'RUWE < 1.4 and no second RV component', mask=(d['ruwe'] < 1.4) & ~np.isfinite(d['rv_comp_2'].astype(float)), levels=lv)[0] # G and K dwarfs proper: no subgiant dilution, slow settling (used for the time-integrated statement) gk = (d['logg'].astype(float) >= 4.0) & (teff < 6000) var += limits_table(d, c, vet, eff, 'G-K dwarfs only: logg >= 4.0 and Teff < 6000 K', mask=gk, levels=[0.75, 1.0, 1.25])[0] var += limits_table(d, c, vet, eff, 'first-pass rule (edge stars refuted by linear extrapolation)', levels=lv, rule='firstpass')[0] var += limits_table(d, c, vet, eff, 'saturation allowance for rho > 1 (N ~ W^2)', levels=lv, rule='sat')[0] w_all = s3.eff_weight(d, c, 'nd', 0.75) lg = d['logg'].astype(float) key['nd_neff_fractions'] = dict(logg_lt_4p0=float(w_all[lg < 4.0].sum() / w_all.sum()), logg_lt_3p8=float(w_all[lg < 3.8].sum() / w_all.sum()), teff_ge_6000=float(w_all[teff >= 6000].sum() / w_all.sum()), gk_dwarfs_logg_ge_4_teff_lt_6000=float(w_all[gk].sum() / w_all.sum())) c_alt = s3.compute(d, sfrac=S_FRAC_ALT) var += limits_table(d, c_alt, vet, eff, 's-fractions Travaglio+04 (T04+)', levels=lv)[0] d50 = dict(d) for e in NCAP: d50['x90_' + e] = d['x50_' + e] c50 = s3.compute(d50) var += limits_table(d50, c50, vet, eff, 'detection threshold x50 instead of x90', levels=lv)[0] vv = [] for e in MAIN: for T in lv: w = s3.eff_weight(d, c, e, T); neff = float(w.sum()) cnt, _ = tally(d, c, vet, e, T, rule='dropworst'); ul, k_exp = upper_limit(cnt, neff, eff) vv.append(dict(variant='refutation rule: no clipping, worst line dropped', element=e, T_excess=T, salt_level_dex=T + s3.MARGIN, N_eff=neff, n_catalogue=cnt['n_catalogue'], n_refuted=cnt['refuted'], n_supported=cnt['supported'], n_ambiguous=cnt['ambiguous'], n_unverifiable=cnt['unverifiable'], not_vetted=cnt['not_vetted'], n_not_refuted_est=k_exp, f_ul95=ul)) var += vv for mg in (0.15, 0.35): vv = [] for e in MAIN: for T in lv: w = s3.eff_weight(d, c, e, T, margin=mg); neff = float(w.sum()) cnt, _ = tally(d, c, vet, e, T, margin=mg); ul, k_exp = upper_limit(cnt, neff, eff) vv.append(dict(variant='margin %.2f dex: this row is the limit at salt level T + %.2f, not a variant of the same level' % (mg, mg), element=e, T_excess=T, salt_level_dex=T + mg, N_eff=neff, n_catalogue=cnt['n_catalogue'], n_refuted=cnt['refuted'], n_supported=cnt['supported'], n_ambiguous=cnt['ambiguous'], n_unverifiable=cnt['unverifiable'], not_vetted=cnt['not_vetted'], n_not_refuted_est=k_exp, f_ul95=ul)) var += vv base_rows = [dict(r, variant='baseline') for r in rows if r['T_excess'] in (0.75, 1.0)] vdf = pd.DataFrame(base_rows + var) vdf = vdf[['variant', 'element', 'T_excess', 'salt_level_dex', 'N_eff', 'n_catalogue', 'n_refuted', 'n_supported', 'n_ambiguous', 'n_unverifiable', 'not_vetted', 'n_not_refuted_est', 'f_ul95']] vdf.to_csv(DATA + 'robustness_variants.csv', index=False, float_format='%.5g') # --- N_eff by Teff bin tb = np.arange(TEFF_MIN, TEFF_MAX + 1, 250.0) with open(DATA + 'neff_by_teff.csv', 'w', newline='') as fh: wr = csv.writer(fh); wr.writerow(['element', 'T_excess', 'salt_level_dex', 'teff_lo', 'teff_hi', 'n_clean_stars', 'N_eff', 'median_x90', 'median_snr_ccd2']) for e in MAIN: for T in s3.LEVELS: w = s3.eff_weight(d, c, e, T) for t0 in tb[:-1]: m = clean & (teff >= t0) & (teff < t0 + 250) x90 = d['x90_' + e][m]; x90 = x90[np.isfinite(x90)] wr.writerow([e, T, T + s3.MARGIN, t0, t0 + 250, int(m.sum()), f'{w[m].sum():.0f}', f'{np.median(x90):.2f}' if len(x90) else '', f'{np.nanmedian(snr[m]):.0f}']) # --- detection limits vs Teff and S/N (median x90 of clean dwarfs) with open(DATA + 'detection_limit_grid.csv', 'w', newline='') as fh: wr = csv.writer(fh); wr.writerow(['element', 'teff_lo', 'teff_hi', 'snr_lo', 'snr_hi', 'n_stars', 'median_x90', 'median_x50']) sb = [10, 20, 30, 50, 80, 150, 1e9] for e in MAIN + ['ba', 'y']: for t0 in tb[:-1]: for a, b in zip(sb[:-1], sb[1:]): m = clean & (teff >= t0) & (teff < t0 + 250) & (snr >= a) & (snr < b) & (np.abs(d['fe_h']) < 0.3) if m.sum() >= 20: wr.writerow([e, t0, t0 + 250, a, b if b < 1e8 else '', int(m.sum()), f"{np.nanmedian(d['x90_' + e][m]):.2f}", f"{np.nanmedian(d['x50_' + e][m]):.2f}"]) # --- survivors: unique (element, star) at the highest level at which they are not refuted surv = {} for e, T, i, k, how in surv_rows: cur = surv.get((e, i)) if cur is None or T > cur[0]: surv[(e, i)] = (T, k, how) esc = set() if os.path.exists(RAW + 'escorza2026_srich.csv'): esc = set(pd.read_csv(RAW + 'escorza2026_srich.csv').GALAH.astype('int64').tolist()) d3 = dict(np.load(CACHE + 'dr3.npz')); o3 = np.argsort(d3['sobject_id']) vi = vet[vet.kind.isin(['candidate', 'candidate_sub'])].drop_duplicates(['element', 'idx']).set_index(['element', 'idx']) out = [] for (e, i), (T, k, how) in sorted(surv.items(), key=lambda kv: (kv[1][1] != 'supported', -kv[1][0], kv[0][0], kv[0][1])): r = vi.loc[(e, i)] p3 = np.searchsorted(d3['sobject_id'][o3], d['sobject_id'][i]); j = o3[min(p3, len(o3) - 1)] in3 = d3['sobject_id'][j] == d['sobject_id'][i] row = dict(element=e, max_level_T=T, verdict=k, vetting=how, sobject_id=int(d['sobject_id'][i]), gaiadr3_source_id=int(d['gaiadr3_source_id'][i]), ra=float(d['ra'][i]), dec=float(d['dec'][i]), teff=float(d['teff'][i]), logg=float(d['logg'][i]), fe_h=float(d['fe_h'][i]), vsini=float(d['vsini'][i]), snr_ccd2=float(d['snr_px_ccd2'][i]), ruwe=float(d['ruwe'][i]), G=float(d['phot_g_mean_mag'][i]), dist_pc=float(d['r_med'][i]), rv_comp_2=float(d['rv_comp_2'][i]), x_cat=float(r['x_cat']), pred=float(r['pred']), ex_cat=float(r['ex_cat']), rho=float(r['rho']), rho_err=float(r['rho_err']), n_lines_used=float(r['n_used']), snr_mod=float(r['snr_mod']), clipped_line=float(r['clipped_line']), chi2_red=float(r['chi2_red']), twin_mode=int(r['twin_mode']), ex_spec=float(r['ex_spec']), ex_spec_lo=float(r['ex_spec_lo']), ex_spec_hi=float(r['ex_spec_hi']), dr3_value=float(d3[e + '_fe'][j]) if in3 else np.nan, dr3_flag=int(d3['flag_' + e + '_fe'][j]) if in3 else -1, in_escorza2026_srich=int(int(d['sobject_id'][i]) in esc)) for el in ['y', 'zr', 'ba', 'la', 'ce', 'nd', 'sm', 'eu']: row[el + '_fe'] = float(d[el + '_fe'][i]); row['flag_' + el] = int(d['flag_' + el + '_fe'][i]) row['sec_' + el] = int(c['sec'][el][i]) for el in ['c_fe', 'ew_h_alpha']: row[el] = float(d[el][i]) out.append(row) sv = pd.DataFrame(out) sv.to_csv(DATA + 'surviving_candidates_all.csv', index=False, float_format='%.4f') # --- Galaxy conversion n_fgk100, n_loc, n_gcns = gcns_density() n_mw = n_loc * MW_STELLAR_MASS / RHO_STAR_LOCAL key['galaxy'] = dict(gcns_fgk_within_100pc=n_fgk100, n_local_per_pc3=n_loc, n_fgk_milky_way=n_mw, mw_stellar_mass=MW_STELLAR_MASS, rho_star_local=RHO_STAR_LOCAL) w = s3.eff_weight(d, c, 'nd', 0.75) > 0 r = d['r_med'][w].astype(float); r = r[np.isfinite(r)] key['nd_sample_distance_pc'] = dict(p10=float(np.percentile(r, 10)), p50=float(np.percentile(r, 50)), p90=float(np.percentile(r, 90))) # FGK dwarfs inside the sphere that contains 90% of the Nd sample (exponential disc, h_z = 300 pc) R = key['nd_sample_distance_pc']['p90'] zz = np.linspace(-R, R, 2001); hz = 300.0 yy = np.pi * (R ** 2 - zz ** 2) * np.exp(-np.abs(zz) / hz) vol_eff = float(np.sum(0.5 * (yy[1:] + yy[:-1]) * np.diff(zz))) key['n_fgk_in_survey_volume'] = float(n_loc * vol_eff); key['survey_volume_radius_pc'] = float(R) # numbers of beacons allowed: baseline (excess / additive definition) and multiplicative definition bk = {} for e in ('nd', 'la'): for r in rows: if r['element'] == e and r['vetted'] and r['salt_level_dex'] >= 1.0: m = sadd[(sadd.element == e) & (sadd.adopted == 1) & (np.abs(sadd.salt_added_dex - r['salt_level_dex']) < 1e-6)] f_lo = r['f_ul95']; f_hi = float(max(m.f_ul95_no_extrapolation.max(), r['f_ul95_no_extrapolation'])) if len(m) else r['f_ul95_no_extrapolation'] bk['%s_%.2f' % (e, r['salt_level_dex'])] = dict( f_baseline=f_lo, f_most_conservative=f_hi, n_survey_volume_baseline=f_lo * key['n_fgk_in_survey_volume'], n_survey_volume_conservative=f_hi * key['n_fgk_in_survey_volume'], n_galaxy_baseline=f_lo * n_mw, n_galaxy_conservative=f_hi * n_mw) key['beacons'] = bk key['cut_flow'] = cut key['limits'] = rows key['fission'] = fis.to_dict('records') key['sigma_core'] = {e: c['sig'][e] for e in NCAP} key['zero_points'] = {e: c['zp'][e] for e in NCAP} key['n_survivor_pairs'] = int(len(sv)); key['n_supported_pairs'] = int((sv.verdict == 'supported').sum()) if len(sv) else 0 json.dump(key, open(DATA + 'key_numbers.json', 'w'), indent=1, default=float) if verbose: pd.set_option('display.width', 260); pd.set_option('display.max_columns', 40); pd.set_option('display.max_rows', 500) print(lim[['element', 'T_excess', 'salt_level_dex', 'N_eff', 'n_catalogue', 'n_refuted', 'n_supported', 'n_ambiguous', 'n_unverifiable', 'band_total', 'band_vetted', 'band_refuted', 'band_supported', 'n_not_refuted_est', 'f_ul95', 'f_ul95_catalogue_only']].to_string(index=False)) print(json.dumps(cal, indent=1)) print(fis.to_string(index=False)) print(vdf.to_string(index=False)) print('survivor (element, star) pairs:', len(sv)) return key, lim, sv if __name__ == '__main__': main()