"""Astrophysical vetting of the candidates whose excess is confirmed in the spectrum. For every (star, element) with verdict 'supported' this collects, in one table: * SIMBAD identification / object type within 5 arcsec; * binarity indicators (Gaia RUWE, second RV component found by the GALAH pipeline); * GALAH DR3 (independent pipeline) value; * membership of the Escorza et al. (2026) s-process-rich list; * a pattern test: the best non-negative s + r mixture fitted to ALL securely detected heavy neutron-capture elements (Ba, La, Ce, Nd, Sm, Eu) with f_s and f_r free, i.e. not tied to the Ba anchor or to the disc Eu trend. rms and the largest residual say whether the star 'follows an s/r pattern' once its other lanthanides are used; the [Eu/Fe] that the fit implies is compared with the star's own Eu detection threshold (an r-process explanation that needs more Eu than the spectrum allows is ruled out); * the same fit for 'solar mix + U-235 fission waste'; * an automatic list of natural / instrumental explanations that apply. A star is called 'unexplained' only if none applies. """ import csv, json, os import numpy as np import pandas as pd from common import * import s03_candidates as s3 import s04_vetting as s4 import s05_physics as s5 import s06_limits as s6 HEAVY = ['ba', 'la', 'ce', 'nd', 'sm', 'eu'] DR3_OFF = {} FIELD_COUNT = {} # plate/night (sobject_id // 1000) -> number of catalogue candidates of any element, clean or flagged FOLLOWUP_SIGMA = 4.5 # significance from which an unexplained candidate without a star-specific defect is listed for follow-up def natural_fit(obs, feh): """obs: dict element -> zero-point-corrected [X/Fe] of securely detected heavy elements. Grid search for the non-negative s + r mixture; returns (log f_s, log f_r, rms, worst element, worst residual).""" els = list(obs.keys()) y = np.array([obs[e] for e in els]) S = np.array([S_FRAC[e] for e in els]); R = 1 - S lfs = np.arange(-1.5, 2.81, 0.05); lfr = np.arange(-1.0, 2.81, 0.05) FS, FR = np.meshgrid(10 ** lfs, 10 ** lfr, indexing='ij') mod = np.log10(S[None, None, :] * FS[:, :, None] + R[None, None, :] * FR[:, :, None]) chi = ((mod - y[None, None, :]) ** 2).sum(-1) i, j = np.unravel_index(np.argmin(chi), chi.shape) res = y - mod[i, j] k = int(np.argmax(np.abs(res))) return float(lfs[i]), float(lfr[j]), float(np.sqrt(np.mean(res ** 2))), els[k], float(res[k]) def fission_fit(obs, ey): """Best 'scaled-solar + fission waste' model: one free parameter (amount of waste).""" els = list(obs.keys()); y = np.array([obs[e] for e in els]) best = (np.inf, 0.0) for dx in np.arange(0.0, 3.01, 0.05): p = s5.fission_pattern(ey, 'nd', max(dx, 1e-3))[0] m = np.array([p[e] for e in els]) chi = float(np.sqrt(np.mean((y - m) ** 2))) if chi < best[0]: best = (chi, dx) return best def anchor_check(d, c, idx_list, elements): """Twin test on the ANCHOR element (Ba for lanthanides, Y for Zr) of each supported candidate: a candidate whose excess exists only because its anchor abundance is low must have anchor lines that are really weaker than those of its twins.""" import spectra masks = s4.line_masks() tw = s4.pick_twins(d, c, np.array(sorted(set(idx_list)))) out = {} for i, e in zip(idx_list, elements): a = s3.ANCHOR[e] S = spectra.read(int(d['sobject_id'][i])) T = [spectra.read(int(d['sobject_id'][j])) for j in tw[i]] xt = np.array([d[a + '_fe'][j] for j, t in zip(tw[i], T) if t is not None], float) T = [t for t in T if t is not None] xa = float(d[a + '_fe'][i]); x_tw = float(np.nanmedian(xt)) if np.isfinite(xt).any() else 0.0 res = dict(anchor=a, anchor_cat=xa, anchor_twin=x_tw, anchor_rho=np.nan, anchor_rho_err=np.nan, anchor_spec=np.nan) if S is not None and len(T) >= 2: r = s4.twin_test(S, T, masks.get(a, [])) L = np.array(r['lines']) if r['lines'] else np.zeros((0, 4)) if len(L) and np.sum(L[:, 3] ** 2 / L[:, 2] ** 2) > 0: den = np.sum(L[:, 3] ** 2 / L[:, 2] ** 2) rho = float(np.sum(L[:, 1] * L[:, 3] / L[:, 2] ** 2) / den) res.update(anchor_rho=rho, anchor_rho_err=float(den ** -0.5), anchor_snr_mod=float(np.sqrt(den)), anchor_spec=float(s4.x_spec(rho, xa, x_tw))) out[(int(i), e)] = res return out def main(verbose=True): d = s3.load(); c = s3.compute(d) sv = pd.read_csv(DATA + 'surviving_candidates_all.csv') try: o = json.load(open(DATA + 'other_surveys_summary.json'))['dr3'] for e in s4.MAIN: if isinstance(o.get(e), dict) and o[e].get('median_dr4_minus_dr3') is not None: DR3_OFF[e] = float(o[e]['median_dr4_minus_dr3']) except Exception: pass key = json.load(open(DATA + 'key_numbers.json')) # candidates per plate and night (clean and flagged stars, any probe element, excess >= 0.75 dex) anyc = np.zeros(len(d['teff']), bool) for e_ in s4.MAIN: ok_ = c['base'] & c['sec'][e_] & c['anchor_ok'][e_] anyc |= ok_ & ((c['ex'][e_] >= 0.75) | (s3.at_edge(d, e_) & (c['ex'][e_] >= 0.5))) uf, cnt = np.unique(d['sobject_id'][anyc] // 1000, return_counts=True) FIELD_COUNT.update(dict(zip(uf.tolist(), cnt.tolist()))) sup = sv[sv.verdict == 'supported'].copy().reset_index(drop=True) idx_of = {int(s): i for i, s in enumerate(d['sobject_id'])} E, chain = s5.chain_yields('U_235'); ey = s5.elemental_yields(chain) tb = c['tb'] anc = anchor_check(d, c, [idx_of[int(x)] for x in sup.sobject_id], list(sup.element)) rows = [] for _, r in sup.iterrows(): i = idx_of[int(r.sobject_id)]; e = r.element t = float(d['teff'][i]); feh = float(d['fe_h'][i]) kb = int(np.clip(np.digitize(t, tb) - 1, 0, len(tb) - 2)) obs = {} for h in HEAVY: if c['sec'][h][i]: zp = c['zp'][h][kb] if h not in s3.BOUND_PROBES else 0.0 obs[h] = float(d[h + '_fe'][i]) - zp if e not in obs and e in HEAVY: obs[e] = float(d[e + '_fe'][i]) out = dict(r) heavy_probe = e in HEAVY if heavy_probe and len(obs) >= 2: lfs, lfr, rms, worst, wres = natural_fit(obs, feh) eu_implied = float(np.log10(S_FRAC['eu'] * 10 ** lfs + (1 - S_FRAC['eu']) * 10 ** lfr)) frms, fdx = fission_fit(obs, ey) out.update(n_secure_heavy=len(obs), secure_heavy='+'.join(obs.keys()), nat_log_fs=lfs, nat_log_fr=lfr, nat_rms=rms, nat_worst=worst, nat_worst_res=wres, nat_eu_implied=eu_implied, eu_x90=float(d['x90_eu'][i]), eu_secure=int(c['sec']['eu'][i]), fission_rms=frms, fission_nd_dex=fdx) else: out.update(n_secure_heavy=len(obs), secure_heavy='+'.join(obs.keys()), nat_log_fs=np.nan, nat_log_fr=np.nan, nat_rms=np.nan, nat_worst='', nat_worst_res=np.nan, nat_eu_implied=np.nan, eu_x90=float(d['x90_eu'][i]), eu_secure=int(c['sec']['eu'][i]), fission_rms=np.nan, fission_nd_dex=np.nan) sim = s6.simbad_cone(float(r.ra), float(r.dec)) out.update(simbad_id=sim[0], simbad_otype=sim[1], simbad_sptype=sim[2], simbad_nbref=sim[3]) # --- natural / instrumental explanations (each one is sufficient) why = [] A = anc[(int(i), e)] out.update(A) # how much of the excess would vanish if the anchor were at the level of the twins? if np.isfinite(A['anchor_spec']) and A['anchor_cat'] < A['anchor_twin'] - 0.15 and A.get('anchor_snr_mod', 0) >= 3: if A['anchor_rho'] + 2 * A['anchor_rho_err'] < 0.5: why.append('anchor [%s/Fe] = %.2f is not confirmed by its lines (rho = %.2f +- %.2f): prediction too low' % (A['anchor'].capitalize(), A['anchor_cat'], A['anchor_rho'], A['anchor_rho_err'])) if np.isfinite(r.rv_comp_2): why.append('SB2 (second RV component in GALAH fit)') if r.ruwe >= 1.4: why.append('Gaia RUWE >= 1.4 (unresolved companion: composite spectrum or mass-transfer binary)') if 'SB' in str(sim[1]) or '**' in str(sim[1]): why.append('SIMBAD: binary (%s)' % sim[1]) if feh < -1.0: why.append('[Fe/H] < -1: halo star, r-process-rich patterns are natural') if heavy_probe and len(obs) >= 3 and out['nat_rms'] <= 0.20 and (out['nat_eu_implied'] <= max(0.6, -0.4 * feh + 0.3) or (out['eu_secure'] == 0 and out['nat_eu_implied'] <= out['eu_x90'])): why.append('fits a free s+r mixture (rms %.2f dex) with an [Eu/Fe] the data allow' % out['nat_rms']) if float(d['ba_fe'][i]) >= 0.6 and c['sec']['ba'][i]: why.append('s-process-rich star ([Ba/Fe] >= +0.6): the probe exceeds the solar-s-pattern prediction, not the s-process') if int(r.in_escorza2026_srich): why.append('in the Escorza+2026 s-process-rich list') if r.vsini >= 20: why.append('vsini >= 20 km/s (blended, shallow lines)') # --- caveats (weaken the candidate without explaining it). Revision (referee minor point): # star-specific caveats (something seen in THIS star's data) are kept apart from generic # ones that are assigned from Teff/logg alone, and marginal significance is recorded # separately, so that 'every candidate has a caveat' is not true by construction. cav = []; gen = [] sigma = float(r.rho / r.rho_err) if r.rho_err > 0 else np.nan nfield = int(FIELD_COUNT.get(int(r.sobject_id) // 1000, 0)) if nfield >= 5: cav.append('%d catalogue candidates on the same plate and night (possible reduction problem)' % nfield) if r.n_lines_used < 2: cav.append('single line') if r.rho > 2.0 and not bool(s3.at_edge(d, e)[i]): cav.append('absorption %.1fx stronger than the catalogue abundance implies (blend or defect)' % r.rho) if r.n_lines_used >= 3 and np.isfinite(r.chi2_red) and r.chi2_red > 3.0: cav.append('lines mutually inconsistent (chi2/dof = %.1f)' % r.chi2_red) if int(r.twin_mode) == 2: cav.append('no close twins (nearest-neighbour comparison stars)') if t < 4900 and float(d['logg'][i]) > 4.2: gen.append('K dwarf: line windows sit in MgH/CN-affected regions (generic, assigned from Teff and logg)') dv = r.dr3_value; df_ = int(r.dr3_flag) off = DR3_OFF.get(e, 0.0) if df_ == 0 and np.isfinite(dv): dr3 = 'agrees' if dv + off >= r.x_cat - 0.3 else 'lower by %.2f' % (r.x_cat - dv - off) elif df_ == 1: dr3 = 'upper limit' elif df_ == -1: dr3 = 'not in DR3' else: dr3 = 'not measured (flag %d)' % df_ if dr3.startswith('lower') or dr3 == 'upper limit': cav.append('GALAH DR3 (independent pipeline): ' + dr3) out.update(significance=sigma, dr3_status=dr3, explanations='; '.join(why), caveats='; '.join(cav + gen), specific_caveats='; '.join(cav), generic_caveats='; '.join(gen), marginal_lt_5sigma=int(sigma < 5), n_candidates_same_field=nfield, unexplained=int(len(why) == 0), unexplained_no_caveat=int(len(why) == 0 and len(cav) == 0 and len(gen) == 0 and sigma >= 5), unexplained_no_specific_caveat=int(len(why) == 0 and len(cav) == 0), followup=int(len(why) == 0 and len(cav) == 0 and sigma >= FOLLOWUP_SIGMA)) rows.append(out) res = pd.DataFrame(rows) res.to_csv(DATA + 'supported_candidates_vetted.csv', index=False, float_format='%.4f') # follow-up list: unexplained, no star-specific defect, significance >= FOLLOWUP_SIGMA, salt level >= 1 dex fu = res[(res.followup == 1) & (res.max_level_T >= 0.75) & (res.vetting == 'full')].sort_values('significance', ascending=False) fu[['element', 'sobject_id', 'gaiadr3_source_id', 'ra', 'dec', 'teff', 'logg', 'fe_h', 'G', 'snr_ccd2', 'x_cat', 'ex_cat', 'rho', 'rho_err', 'significance', 'ex_spec', 'dr3_status', 'generic_caveats']].to_csv(DATA + 'followup_targets.csv', index=False, float_format='%.4f') summ = {} vet = pd.read_csv(DATA + 'vetting_twin_test.csv') if len(res): for T in (0.5, 0.75, 1.0): g = res[(res.max_level_T >= T) & (res.vetting == 'full')] if T >= 0.75 else res[res.max_level_T >= T] by = {} for e in s4.MAIN: ge = g[g.element == e] n_tested = int(((vet.kind == 'candidate') & (vet.element == e) & (vet.status == 'ok') & (vet.ex_cat >= T)).sum()) fp = key['twin_test']['null'][e] by[e] = dict(n_supported=int(len(ge)), n_explained=int((ge.unexplained == 0).sum()), n_unexplained=int(ge.unexplained.sum()), n_unexplained_no_caveat=int(ge.unexplained_no_caveat.sum()), n_ge_5sigma=int((ge.significance >= 5).sum()), n_candidates_tested=n_tested, null_false_positive_rate=fp['frac_z_gt3'], null_false_positive_rate_ul95=fp['frac_z_gt3_ul95'], expected_chance_supports=n_tested * fp['frac_z_gt3'], expected_chance_supports_ul95=n_tested * fp['frac_z_gt3_ul95']) summ['T>=%.2f' % T] = dict(n_supported=int(len(g)), n_unexplained=int(g.unexplained.sum()), n_unexplained_no_caveat=int(g.unexplained_no_caveat.sum()), n_unexplained_no_specific_caveat=int(g.unexplained_no_specific_caveat.sum()), n_unexplained_only_generic_or_marginal=int(((g.unexplained == 1) & (g.specific_caveats.fillna('') == '')).sum()), n_followup=int(g.followup.sum()), followup=[dict(element=r_.element, sobject_id=int(r_.sobject_id), significance=float(r_.significance)) for _, r_ in g[g.followup == 1].iterrows()], by_element=by) json.dump(summ, open(DATA + 'supported_candidates_summary.json', 'w'), indent=1) if verbose: pd.set_option('display.width', 280); pd.set_option('display.max_columns', 60); pd.set_option('display.max_rows', 400) print(json.dumps(summ, indent=1)) cols = ['element', 'max_level_T', 'vetting', 'sobject_id', 'teff', 'logg', 'fe_h', 'snr_ccd2', 'x_cat', 'ex_cat', 'rho', 'rho_err', 'significance', 'ex_spec', 'ba_fe', 'nat_rms', 'dr3_status', 'simbad_otype', 'n_candidates_same_field', 'followup', 'explanations', 'specific_caveats', 'generic_caveats'] if len(res): print(res[cols].to_string(index=False)) return res, summ if __name__ == '__main__': main()