#!/usr/bin/env python
"""Capstone figure: three independent methods converge — most (U)C HII ionizing
stars are NOT accreting. Left: per-source ALMA-IMF (continuum alpha vs hot-core
line S/N). Right: accreting fraction across the three methods."""
import numpy as np
import matplotlib; matplotlib.use('Agg')
import matplotlib.pyplot as plt
from astropy.table import Table

cb=Table.read('almaimf_alpha_commonbeam.fits')
ln=Table.read('phase6_line_signatures.fits')
# match by Name+field
lk={str(ln['Name'][i])+str(ln['field'][i]):i for i in range(len(ln))}
rows=[]
for i in range(len(cb)):
    k=str(cb['Name'][i])+str(cb['field'][i])
    if k in lk:
        j=lk[k]
        ch=max(float(ln['CH3CN'][j]),float(ln['CH3CNv8'][j]))
        rows.append((str(cb['Name'][i]),float(cb['alpha_cb'][i]),ch,float(ln['H41a'][j]),float(ln['SiO'][j])))
rows.sort(key=lambda r:r[1])
names=[r[0] for r in rows]; al=np.array([r[1] for r in rows]); ch=np.array([r[2] for r in rows])

fig,(axL,axR)=plt.subplots(1,2,figsize=(14,5.5),gridspec_kw={'width_ratios':[1.7,1]})

# ---- Left: per-source continuum alpha (color) + hot-core S/N ----
y=np.arange(len(names))
# accreting if dust (alpha>2.5) OR hot core (CH3CN S/N>5)
acc=(al>2.5)|(ch>5)
colors=['#c0392b' if a else '#2980b9' for a in acc]
axL.barh(y-0.18,al,height=0.34,color=['#c0392b' if al[i]>2.5 else '#85c1e9' for i in range(len(al))],label='continuum α (B3→B6)')
ax2=axL.twiny()
ax2.barh(y+0.18,np.clip(ch,0,80),height=0.34,color=['#c0392b' if ch[i]>5 else '#bdc3c7' for i in range(len(ch))],label='hot-core CH3CN S/N')
axL.axvline(2.5,ls='--',color='#c0392b',lw=1); axL.text(2.55,len(y)-0.5,'dust\nα>2.5',color='#c0392b',fontsize=7)
axL.axvline(-0.1,ls=':',color='#2471a3',lw=1)
ax2.axvline(5,ls='--',color='#7f8c8d',lw=1); ax2.text(5.5,-0.9,'hot-core\nthreshold 5σ',color='#7f8c8d',fontsize=7)
axL.set_yticks(y); axL.set_yticklabels(names,fontsize=8)
axL.set_xlabel('common-beam spectral index α  (blue bars)',fontsize=9)
ax2.set_xlabel('hot-core CH3CN line S/N  (clipped at 80, gray/red bars)',fontsize=9)
axL.set_xlim(-1.5,4.2)
for i,n in enumerate(names):
    if acc[i]: axL.text(4.15,y[i],'★ active',color='#c0392b',fontsize=8,va='center',ha='right',fontweight='bold')
axL.set_title('ALMA-IMF: dust + hot core at the ionizing source\n'
              '7/9 free-free & no hot core (not accreting); 2 active (★, cluster-confused)',fontsize=10)

# ---- Right: accreting fraction, three methods ----
meth=['Catalog\n(ALMAGAL dust\ncore at star,\nff-subtracted)\nN=30',
      'Continuum α\n(dust-dominated\nα>2.5)\nN=8',
      'Lines\n(hot-core CH3CN\n>5σ at star)\nN=9']
frac=[3,12,22]   # % accreting (1/30, 1/8, 2/9)
axR.bar(range(3),frac,color=['#16a085','#e67e22','#8e44ad'])
for i,f in enumerate(frac): axR.text(i,f+1,f'{f}%',ha='center',fontweight='bold')
axR.set_xticks(range(3)); axR.set_xticklabels(meth,fontsize=7.5)
axR.set_ylabel('% of (U)C HII with accretion signature\nat the ionizing star')
axR.set_ylim(0,100)
axR.axhspan(0,25,color='#2980b9',alpha=0.08)
axR.text(1,90,'Three independent methods\nagree: the MAJORITY are\nNOT accreting',ha='center',fontsize=9,
         bbox=dict(boxstyle='round',fc='#fdf2e9',ec='#e67e22'))
axR.set_title('Accreting fraction — triple convergence',fontsize=10)
plt.tight_layout(); plt.savefig('fig_convergence.png',dpi=135)
print('wrote fig_convergence.png')