#!/usr/bin/env python3
"""
Voorspelling #17 v2: Model-Onafhankelijke Residueel-Analyse
=============================================================

Verbeterde aanpak: in plaats van een theoretische baseline (EPDL97),
fitten we een GLADDE SPLINE aan de experimentele data en zoeken we
naar systematische afwijkingen bij harmonische energieën.

Voordeel: model-onafhankelijk. We vragen niet "klopt de theorie?"
maar "zit er structuur in de data die een gladde curve niet vangt?"

Data: Basaglia et al. (2022), IEEE Trans. Nucl. Sci., 69(4), 858-873.
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import UnivariateSpline
from scipy.signal import find_peaks
import os

m_e = 0.511
THRESHOLD = 2 * m_e  # 1.022 MeV

# ═══════════════════════════════════════════════════════════════
# DATA — Geëxtraheerd uit Basaglia et al. figuren
# ═══════════════════════════════════════════════════════════════

ge_data = {
    'energy': np.array([
        1.12, 1.17, 1.22, 1.28, 1.33, 1.40, 1.46, 1.55,
        1.60, 1.69, 1.78, 1.88, 1.97,
        2.04, 2.15, 2.25, 2.35, 2.45, 2.55, 2.65, 2.75,
        2.85, 2.95,
        3.10, 3.30, 3.50, 3.70, 4.00, 4.20, 4.50, 4.80,
        5.00, 5.50, 6.00, 6.50
    ]),
    'sigma': np.array([
        18, 55, 85, 140, 190, 270, 340, 450,
        510, 620, 730, 830, 920,
        960, 1050, 1120, 1185, 1240, 1290, 1330, 1365,
        1395, 1420,
        1460, 1500, 1530, 1550, 1575, 1590, 1610, 1625,
        1635, 1650, 1660, 1665
    ]),
    'err': np.array([
        5, 10, 12, 18, 20, 25, 28, 32,
        34, 36, 40, 44, 46,
        48, 50, 50, 50, 50, 50, 50, 50,
        50, 50,
        50, 50, 50, 50, 50, 50, 50, 50,
        50, 50, 50, 50
    ]),
    'Z': 32, 'name': 'Germanium'
}

pb_data = {
    'energy': np.array([
        1.10, 1.12, 1.15, 1.19, 1.22, 1.28, 1.33, 1.40,
        1.50, 1.60, 1.70, 1.78, 1.88, 2.00, 2.10, 2.20,
        2.35, 2.50, 2.62, 2.75
    ]),
    'sigma': np.array([
        80, 130, 380, 680, 900, 1320, 1650, 2100,
        2680, 3100, 3450, 3650, 3900, 4100, 4250, 4350,
        4500, 4600, 4650, 4700
    ]),
    'err': np.array([
        30, 35, 55, 65, 70, 80, 90, 100,
        100, 100, 100, 100, 100, 100, 100, 100,
        100, 100, 100, 100
    ]),
    'Z': 82, 'name': 'Lead'
}

sn_data = {
    'energy': np.array([
        1.12, 1.17, 1.22, 1.28, 1.33, 1.40, 1.50, 1.60,
        1.70, 1.78, 1.88, 2.00, 2.10, 2.20, 2.35, 2.50, 2.75
    ]),
    'sigma': np.array([
        25, 70, 145, 250, 340, 450, 580, 700,
        800, 870, 950, 1040, 1090, 1130, 1185, 1220, 1270
    ]),
    'err': np.array([
        10, 18, 25, 30, 35, 40, 45, 50,
        50, 50, 50, 50, 50, 50, 50, 50, 50
    ]),
    'Z': 50, 'name': 'Tin'
}


def smooth_baseline(E, sigma, sigma_err, smoothing_factor=None):
    """
    Fit een gladde spline als baseline.

    De spline is opzettelijk GLAD — het vangt de globale trend maar
    niet lokale structuur. Als er resonantiepieken zijn, verschijnen
    ze als residuelen boven de spline.

    smoothing_factor: hogere waarde = gladdere fit (minder gevoelig
    voor lokale structuur). We kiezen bewust een gladde fit zodat
    eventuele pieken niet worden "weggefit."
    """
    weights = 1.0 / sigma_err

    if smoothing_factor is None:
        # Automatisch: gebruik len(data) als smoothing
        smoothing_factor = len(E) * 1.5

    spline = UnivariateSpline(E, sigma, w=weights, s=smoothing_factor)
    return spline


def harmonic_proximity_test(E, residuals, harmonics, window=0.08):
    """
    Test of residuelen systematisch HOGER zijn nabij harmonische
    energieën dan elders.

    Dit is de kerntest: als het grid-model klopt, verwachten we
    positieve residuelen rond n × 1.022 MeV.

    Returns:
        mean_at_harmonics: gemiddeld residueel nabij harmonischen
        mean_elsewhere: gemiddeld residueel elders
        p_value: significantie van het verschil (t-test benadering)
    """
    near_harmonic = np.zeros(len(E), dtype=bool)
    for E_h in harmonics:
        near_harmonic |= (np.abs(E - E_h) < window)

    res_at = residuals[near_harmonic] if np.any(near_harmonic) else np.array([0])
    res_away = residuals[~near_harmonic] if np.any(~near_harmonic) else np.array([0])

    mean_at = np.mean(res_at)
    mean_away = np.mean(res_away)

    # Vereenvoudigde significantie
    if len(res_at) > 1 and len(res_away) > 1:
        pooled_std = np.sqrt(
            (np.var(res_at) * len(res_at) + np.var(res_away) * len(res_away)) /
            (len(res_at) + len(res_away))
        )
        if pooled_std > 0:
            t_stat = (mean_at - mean_away) / (pooled_std * np.sqrt(1/len(res_at) + 1/len(res_away)))
        else:
            t_stat = 0
    else:
        t_stat = 0

    return mean_at, mean_away, t_stat, near_harmonic


def create_figure(output_dir):
    """Hoofdfiguur: 3×2 paneel analyse."""

    fig, axes = plt.subplots(3, 2, figsize=(18, 20))
    fig.suptitle(
        'Voorspelling #17: Helix-Resonanties — Model-Onafhankelijke Analyse\n'
        'Methode: gladde spline baseline → residuelen → test bij harmonischen van 1.022 MeV\n'
        'Data: Basaglia et al. (2022), 367 metingen, EPDL97 als referentie',
        fontsize=13, fontweight='bold', y=0.99
    )

    datasets = [ge_data, sn_data, pb_data]
    row_labels = ['Ge (Z=32)', 'Sn (Z=50)', 'Pb (Z=82)']

    print("\n" + "="*70)
    print("RESULTATEN: MODEL-ONAFHANKELIJKE RESIDUEEL-ANALYSE")
    print("="*70)

    for i, (data, label) in enumerate(zip(datasets, row_labels)):
        E = data['energy']
        sigma = data['sigma']
        err = data['err']
        Z = data['Z']

        # Harmonischen in het energiebereik
        harmonics = [n * THRESHOLD for n in range(1, 8)
                     if E.min() - 0.1 < n * THRESHOLD < E.max() + 0.1]

        # Fit gladde baseline
        spline = smooth_baseline(E, sigma, err)
        E_fine = np.linspace(E.min(), E.max(), 500)
        sigma_smooth = spline(E_fine)
        sigma_smooth_at_data = spline(E)

        # Residuelen
        residuals = (sigma - sigma_smooth_at_data) / err

        # Harmonische proximity test
        mean_at, mean_away, t_stat, near_harm = harmonic_proximity_test(
            E, residuals, harmonics, window=0.08
        )

        # ─── Links: Data + Spline ─────────────────────────────
        ax = axes[i, 0]
        ax.errorbar(E, sigma, yerr=err, fmt='ko', markersize=5,
                     capsize=3, linewidth=0.8, zorder=5,
                     label=f'Experiment ({data["name"]})')
        ax.plot(E_fine, sigma_smooth, 'b-', linewidth=2.5, alpha=0.8,
                label='Gladde spline baseline')

        # Markeer punten nabij harmonischen
        if np.any(near_harm):
            ax.errorbar(E[near_harm], sigma[near_harm], yerr=err[near_harm],
                         fmt='rs', markersize=8, capsize=3, linewidth=1,
                         zorder=6, label='Nabij harmonische')

        for n, E_h in enumerate(harmonics, 1):
            ax.axvline(E_h, color='red', alpha=0.35, linestyle='--', linewidth=1.5)
            ax.text(E_h, ax.get_ylim()[0] if i == 0 else sigma.min(),
                    f'  n={n}', fontsize=8, color='red', va='bottom',
                    rotation=90, alpha=0.7)

        ax.set_xlabel('Fotonenergie (MeV)', fontsize=11)
        ax.set_ylabel('Cross-section (mb)', fontsize=11)
        ax.set_title(f'{label} — Data + Gladde Baseline', fontsize=12)
        ax.legend(fontsize=9, loc='lower right')
        ax.grid(True, alpha=0.3)

        # ─── Rechts: Residuelen ───────────────────────────────
        ax = axes[i, 1]

        colors = np.where(near_harm, 'red', 'steelblue')
        ax.bar(E, residuals, width=(E.max()-E.min())/len(E)*0.6,
               color=colors, alpha=0.75, edgecolor='black', linewidth=0.3)

        ax.axhline(0, color='black', linewidth=1)
        ax.axhline(2, color='orange', linewidth=0.8, linestyle='--',
                    alpha=0.7, label='±2σ')
        ax.axhline(-2, color='orange', linewidth=0.8, linestyle='--', alpha=0.7)

        for n_h, E_h in enumerate(harmonics, 1):
            ax.axvline(E_h, color='red', alpha=0.5, linestyle='--', linewidth=2)
            ax.text(E_h, 3.8, f'n={n_h}', fontsize=9, color='red',
                    ha='center', fontweight='bold',
                    bbox=dict(boxstyle='round,pad=0.2', facecolor='lightyellow',
                              edgecolor='red', alpha=0.8))

        # Annoteer test-resultaat
        ax.text(0.98, 0.95,
                f'Nabij harmonischen: {mean_at:.2f}σ\n'
                f'Elders: {mean_away:.2f}σ\n'
                f't-statistiek: {t_stat:.2f}',
                transform=ax.transAxes, fontsize=10, va='top', ha='right',
                bbox=dict(boxstyle='round,pad=0.5', facecolor='lightyellow',
                          edgecolor='gray', alpha=0.9))

        ax.set_xlabel('Fotonenergie (MeV)', fontsize=11)
        ax.set_ylabel('Residueel (σ)', fontsize=11)
        ax.set_title(f'{label} — Residuelen bij Harmonischen', fontsize=12)
        ax.set_ylim(-4, 4)
        ax.legend(fontsize=9, loc='upper left')
        ax.grid(True, alpha=0.3)

        # Print resultaten
        print(f"\n{label}:")
        print(f"  Datapunten: {len(E)}, Harmonischen in bereik: {len(harmonics)}")
        print(f"  Residuelen nabij harmonischen: gemiddeld {mean_at:.2f}σ")
        print(f"  Residuelen elders:             gemiddeld {mean_away:.2f}σ")
        print(f"  t-statistiek (verschil):       {t_stat:.2f}")
        if abs(t_stat) > 2:
            print(f"  → SIGNIFICANT verschil (|t| > 2)")
        elif abs(t_stat) > 1:
            print(f"  → Suggestief verschil (|t| > 1)")
        else:
            print(f"  → Geen significant verschil")
        print(f"  Harmonischen: {', '.join([f'{h:.3f} MeV' for h in harmonics])}")

    plt.tight_layout(rect=[0, 0, 1, 0.95])

    filepath = os.path.join(output_dir, 'voorspelling_17_echte_data_v2.png')
    plt.savefig(filepath, dpi=150, bbox_inches='tight', facecolor='white')
    plt.close()
    print(f"\nSaved: {filepath}")

    # ─── Gecombineerde test over alle elementen ───────────────
    print(f"\n{'='*70}")
    print("GECOMBINEERDE TEST: ALLE ELEMENTEN")
    print(f"{'='*70}")

    all_res_at = []
    all_res_away = []

    for data in datasets:
        E = data['energy']
        sigma = data['sigma']
        err = data['err']
        spline = smooth_baseline(E, sigma, err)
        residuals = (sigma - spline(E)) / err

        harmonics = [n * THRESHOLD for n in range(1, 8)
                     if E.min() - 0.1 < n * THRESHOLD < E.max() + 0.1]

        _, _, _, near_harm = harmonic_proximity_test(E, residuals, harmonics)

        all_res_at.extend(residuals[near_harm])
        all_res_away.extend(residuals[~near_harm])

    all_res_at = np.array(all_res_at)
    all_res_away = np.array(all_res_away)

    if len(all_res_at) > 0 and len(all_res_away) > 0:
        combined_mean_at = np.mean(all_res_at)
        combined_mean_away = np.mean(all_res_away)
        pooled_std = np.sqrt(
            (np.var(all_res_at)*len(all_res_at) + np.var(all_res_away)*len(all_res_away)) /
            (len(all_res_at) + len(all_res_away))
        )
        combined_t = (combined_mean_at - combined_mean_away) / \
                     (pooled_std * np.sqrt(1/len(all_res_at) + 1/len(all_res_away))) \
                     if pooled_std > 0 else 0

        print(f"\nAlle datapunten nabij harmonischen: {len(all_res_at)}")
        print(f"Alle datapunten elders:            {len(all_res_away)}")
        print(f"Gemiddeld residueel nabij harm.:   {combined_mean_at:.3f}σ")
        print(f"Gemiddeld residueel elders:        {combined_mean_away:.3f}σ")
        print(f"Gecombineerde t-statistiek:        {combined_t:.3f}")

        if abs(combined_t) > 2.5:
            verdict = "SIGNIFICANT — data toont systematisch hogere cross-sections nabij harmonischen"
        elif abs(combined_t) > 1.5:
            verdict = "SUGGESTIEF — er is een trend maar niet significant met deze dataset"
        else:
            verdict = "NIET SIGNIFICANT — geen detecteerbaar verschil met deze dataset"

        print(f"\nVERDICT: {verdict}")

    print(f"""
{'─'*70}
EERLIJKE BEOORDELING:
{'─'*70}
Deze analyse is INDICATIEF, niet definitief, omdat:

1. De data is geëxtraheerd uit figuren, niet uit originele datasets
2. De foutschattingen zijn benaderingen
3. De gladde-spline methode is gevoelig voor de smoothing parameter
4. Het aantal datapunten nabij elke harmonische is klein (1-3 punten)

WAT WE NODIG HEBBEN:
→ Originele Ge data uit EXFOR (170 punten met exacte fouten)
→ Hogere energieresolutie scans nabij 1.022, 2.044, 3.066 MeV
→ Statistisch krachtigere test (bijv. Bayesiaanse model-vergelijking)

WAT DIT PAPER BEVESTIGT (onafhankelijk van onze analyse):
→ Basaglia et al. vinden dat XCOM, Geant4, PHOTX significant
  afwijken van experiment nabij de drempel
→ De afwijking is het sterkst bij lage energieën (< 1.5 MeV)
→ Dit is ONVERKLAARD in het paper
→ Het grid-model biedt een mogelijke verklaring: resonantiestructuur
""")

    return filepath


if __name__ == "__main__":
    output_dir = os.path.dirname(os.path.abspath(__file__))
    filepath = create_figure(output_dir)
    print(f"✓ Analyse voltooid: {filepath}")