#!/usr/bin/env python3
"""
Voorspelling #17: Helix-Stabiliteitsdrempel bij Paarproductie
==============================================================
Spectrum Grid Model — Analyse Script

Hypothese:
    Het sub-quantum grid model voorspelt dat paarproductie (γ → e⁺e⁻)
    niet een gladde drempel heeft maar scherpe resonantiepieken op
    harmonischen van de fundamentele helix-frequentie:
        E_n = n × 2m_e·c²    voor n = 1, 2, 3, ...

    De pieken zouden een Breit-Wigner profiel hebben met onbekende
    breedte Γ (te bepalen uit data).

Methode:
    1. Bereken QED Bethe-Heitler paarproductie cross-section (standaardmodel)
    2. Voeg grid-model resonantiepieken toe op harmonischen
    3. Simuleer experimentele data met realistische fouten
    4. Analyseer residuelen (data - QED) op periodieke structuur
    5. Gebruik Lomb-Scargle periodogram voor signaaldetectie

Databron:
    - CERN Open Data Portal (opendata.cern.ch)
    - HEPData (hepdata.net) — pair production cross-sections
    - Specifiek: LEP/OPAL/DELPHI paarproductie-metingen

Auteur: Marald Bes / Spectrum van Alles
Datum: 2026-04-08
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks
from scipy.optimize import curve_fit
import os

# ─── Natuurconstanten ───────────────────────────────────────
m_e = 0.511  # MeV/c² — elektronmassa
c = 1.0      # natuurlijke eenheden
alpha = 1/137.036  # fijnstructuurconstante
r_e = 2.818e-13  # cm — klassieke elektronradius

THRESHOLD = 2 * m_e  # 1.022 MeV — minimale energie voor paarproductie

# ─── QED Bethe-Heitler Cross-Section (vereenvoudigd) ────────
def bethe_heitler_simplified(E_gamma, Z=1):
    """
    Vereenvoudigde Bethe-Heitler paarproductie cross-section.

    σ_BH ∝ α · r_e² · Z² · [28/9 · ln(2E/m_e) - 218/27]

    Geldig voor E_gamma >> 2m_e·c² (hoge-energie limiet).
    Nabij drempel gebruiken we de exacte Racah-formule benadering.

    Parameters:
        E_gamma: fotonenergie in MeV
        Z: atoomnummer van target (1 = waterstof)

    Returns:
        cross-section in cm² (of 0 als E < drempel)
    """
    if np.isscalar(E_gamma):
        E_gamma = np.array([E_gamma])

    sigma = np.zeros_like(E_gamma, dtype=float)
    mask = E_gamma > THRESHOLD

    E = E_gamma[mask]
    k = E / m_e  # gereduceerde energie

    # Nabij drempel: exacte benadering (Racah/Maximon)
    beta = np.sqrt(1 - (2*m_e/E)**2)  # snelheid van e+e- in CM frame
    beta = np.clip(beta, 1e-10, 1-1e-10)

    # Coulomb-correctie voor lage beta
    coulomb = (np.pi * alpha * Z / beta) / (1 - np.exp(-np.pi * alpha * Z / beta))

    # Hoge-energie limiet cross-section
    sigma_high = alpha * r_e**2 * Z**2 * (28.0/9.0 * np.log(2*k) - 218.0/27.0)
    sigma_high = np.maximum(sigma_high, 0)

    # Drempel-onderdrukking (fase-ruimte factor)
    phase_space = beta * (3 - beta**4) / 2  # benadering

    sigma[mask] = sigma_high * phase_space * coulomb

    return sigma


# ─── Grid-Model Resonantie ──────────────────────────────────
def breit_wigner(E, E_res, Gamma, A):
    """
    Breit-Wigner resonantieprofiel.

    Beschrijft een resonantiepiek bij energie E_res met breedte Γ
    en amplitude A. Dit is het verwachte profiel als een helixgolf
    in het grid een stabiele configuratie vormt bij een harmonische
    van de fundamentele massa-energie.

    Parameters:
        E: energie (MeV)
        E_res: resonantie-energie (MeV)
        Gamma: breedte (MeV) — onbekende parameter
        A: amplitude — onbekende parameter
    """
    return A * (Gamma/2)**2 / ((E - E_res)**2 + (Gamma/2)**2)


def grid_model_cross_section(E_gamma, Z=1, n_harmonics=6,
                              Gamma=0.05, amplitude_ratio=0.001):
    """
    Spectrum Grid Model paarproductie cross-section.

    QED (Bethe-Heitler) PLUS resonantiepieken op harmonischen
    van 2m_e·c² — de voorspelling van het helix-stabiliteitsmodel.

    Parameters:
        E_gamma: fotonenergie in MeV
        Z: atoomnummer target
        n_harmonics: aantal harmonischen om te berekenen
        Gamma: resonantiebreedte in MeV (vrije parameter)
        amplitude_ratio: piekamplitude als fractie van QED cross-section

    De amplitude_ratio is de cruciale vrije parameter:
        - 0.01  = 1%  afwijking (makkelijk detecteerbaar)
        - 0.001 = 0.1% afwijking (meetbaar met goede statistiek)
        - 0.0001 = 0.01% (extreem moeilijk, maar CERN-precisie)
    """
    sigma_qed = bethe_heitler_simplified(E_gamma, Z)
    sigma_grid = np.copy(sigma_qed)

    for n in range(1, n_harmonics + 1):
        E_res = n * THRESHOLD  # harmonische: n × 1.022 MeV
        # Amplitude schaalt met 1/n² (hogere harmonischen zwakker)
        A_n = amplitude_ratio * np.max(sigma_qed) / n**2
        sigma_grid += breit_wigner(E_gamma, E_res, Gamma, A_n)

    return sigma_grid


# ─── Simulatie van Experimentele Data ────────────────────────
def simulate_experiment(E_range, n_events=10000, grid_signal=True,
                        Gamma=0.05, amplitude_ratio=0.001):
    """
    Simuleer experimentele data met realistische Poisson-fouten.

    Parameters:
        E_range: tuple (E_min, E_max) in MeV
        n_events: totaal aantal events
        grid_signal: True = data bevat grid-resonanties
        Gamma: resonantiebreedte
        amplitude_ratio: signaalsterkte

    Returns:
        E_bins: energiebins (centers)
        counts: getelde events per bin
        errors: Poisson-fouten
        expected_qed: QED-verwachting per bin
    """
    n_bins = 200
    E_bins = np.linspace(E_range[0], E_range[1], n_bins)
    dE = E_bins[1] - E_bins[0]

    # Ware cross-section
    if grid_signal:
        sigma_true = grid_model_cross_section(E_bins, Gamma=Gamma,
                                               amplitude_ratio=amplitude_ratio)
    else:
        sigma_true = bethe_heitler_simplified(E_bins)

    # QED verwachting (altijd)
    sigma_qed = bethe_heitler_simplified(E_bins)

    # Normaliseer naar events
    if np.sum(sigma_true) > 0:
        rate = sigma_true / np.sum(sigma_true) * n_events
    else:
        rate = np.zeros_like(E_bins)

    rate_qed = sigma_qed / np.sum(sigma_qed) * n_events if np.sum(sigma_qed) > 0 else np.zeros_like(E_bins)

    # Poisson-sampling
    counts = np.random.poisson(np.maximum(rate, 0))
    errors = np.sqrt(np.maximum(counts, 1))

    return E_bins, counts, errors, rate_qed


# ─── Residueel-Analyse ──────────────────────────────────────
def analyze_residuals(E_bins, counts, expected_qed):
    """
    Analyseer residuelen (data - QED verwachting) op periodieke structuur.

    Als het grid-model correct is, zouden de residuelen pieken vertonen
    op regelmatige energieafstanden (harmonischen van 1.022 MeV).

    Returns:
        residuals: (counts - expected) / expected
        peaks: indices van gevonden pieken
        significance: significantie van elke piek (sigma)
    """
    # Vermijd deling door nul
    mask = expected_qed > 5  # minimaal 5 verwachte events per bin

    residuals = np.zeros_like(counts, dtype=float)
    residuals[mask] = (counts[mask] - expected_qed[mask]) / np.sqrt(expected_qed[mask])

    # Zoek pieken in residuelen (> 2σ)
    peaks, properties = find_peaks(residuals, height=2.0, distance=3)

    return residuals, peaks, properties.get('peak_heights', [])


# ─── Periodiciteitscanning ───────────────────────────────────
def scan_periodicity(E_bins, residuals, mask=None):
    """
    Scan residuelen op periodieke structuur met Lomb-Scargle-achtige methode.

    Als het grid model klopt, verwachten we een sterke piek bij
    periodiciteit = 1.022 MeV (= 2m_e·c²) en sub-harmonischen.

    Returns:
        test_periods: geteste periodiciteiten (MeV)
        power: power spectrum
        best_period: sterkste periodiciteit
    """
    if mask is None:
        mask = np.ones(len(E_bins), dtype=bool)

    E = E_bins[mask]
    R = residuals[mask]

    # Test periodiciteiten van 0.5 tot 5 MeV
    test_periods = np.linspace(0.3, 5.0, 500)
    power = np.zeros_like(test_periods)

    for i, P in enumerate(test_periods):
        # Phase-fold residuals op deze periodiciteit
        phase = (E % P) / P

        # Bereken power als gemiddelde van residuelen in fase-bins
        n_phase_bins = 10
        phase_bins = np.linspace(0, 1, n_phase_bins + 1)

        for j in range(n_phase_bins):
            in_bin = (phase >= phase_bins[j]) & (phase < phase_bins[j+1])
            if np.sum(in_bin) > 0:
                power[i] += np.mean(R[in_bin])**2

    best_idx = np.argmax(power)
    best_period = test_periods[best_idx]

    return test_periods, power, best_period


# ─── Hoofdvisualisatie ───────────────────────────────────────
def create_analysis_figure(output_dir):
    """
    Maak de volledige analyse-figuur met 4 panelen:

    A) QED vs Grid-model cross-sections
    B) Gesimuleerde data met grid-signaal
    C) Residuelen met gemarkeerde pieken
    D) Periodiciteitsscan
    """
    np.random.seed(42)

    fig, axes = plt.subplots(2, 2, figsize=(16, 12))
    fig.suptitle('Voorspelling #17: Helix-Stabiliteitsdrempel bij Paarproductie\n'
                 'Spectrum Sub-quantum Grid Model — Analyse',
                 fontsize=14, fontweight='bold', y=0.98)

    # ─── Panel A: Theorie ─────────────────────────────────────
    ax = axes[0, 0]
    E = np.linspace(0.8, 8.0, 2000)

    sigma_qed = bethe_heitler_simplified(E)
    sigma_grid = grid_model_cross_section(E, Gamma=0.03, amplitude_ratio=0.005)

    ax.plot(E, sigma_qed, 'b-', linewidth=2, label='QED (Bethe-Heitler)')
    ax.plot(E, sigma_grid, 'r-', linewidth=1.5, alpha=0.8, label='Grid-model (QED + harmonischen)')

    # Markeer harmonischen
    for n in range(1, 7):
        E_harm = n * THRESHOLD
        if E_harm < 8.0:
            ax.axvline(E_harm, color='red', alpha=0.3, linestyle='--', linewidth=0.8)
            ax.text(E_harm, ax.get_ylim()[1]*0.95 if n == 1 else ax.get_ylim()[1]*(0.95 - 0.07*(n-1)),
                    f'n={n}', fontsize=8, ha='center', color='red')

    ax.set_xlabel('Fotonenergie (MeV)')
    ax.set_ylabel('Cross-section (a.u.)')
    ax.set_title('A) QED vs Grid-model Cross-Section')
    ax.legend(fontsize=9)
    ax.set_xlim(0.8, 8.0)
    ax.grid(True, alpha=0.3)

    # ─── Panel B: Gesimuleerde data ───────────────────────────
    ax = axes[0, 1]

    E_bins, counts_signal, errors_s, expected_qed = simulate_experiment(
        (0.8, 8.0), n_events=50000, grid_signal=True,
        Gamma=0.03, amplitude_ratio=0.003
    )

    ax.errorbar(E_bins, counts_signal, yerr=errors_s, fmt='k.', markersize=2,
                alpha=0.5, linewidth=0.5, label='Gesimuleerde data (met grid-signaal)')
    ax.plot(E_bins, expected_qed, 'b-', linewidth=2, label='QED verwachting')

    ax.set_xlabel('Fotonenergie (MeV)')
    ax.set_ylabel('Events per bin')
    ax.set_title('B) Gesimuleerde Experimentele Data (N=50k)')
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)

    # ─── Panel C: Residuelen ──────────────────────────────────
    ax = axes[1, 0]

    residuals, peaks, peak_heights = analyze_residuals(E_bins, counts_signal, expected_qed)

    ax.plot(E_bins, residuals, 'k-', linewidth=0.8, alpha=0.7)
    ax.axhline(0, color='blue', linewidth=1)
    ax.axhline(2, color='orange', linewidth=0.8, linestyle='--', alpha=0.7, label='2σ drempel')
    ax.axhline(-2, color='orange', linewidth=0.8, linestyle='--', alpha=0.7)
    ax.axhline(3, color='red', linewidth=0.8, linestyle='--', alpha=0.7, label='3σ drempel')
    ax.axhline(-3, color='red', linewidth=0.8, linestyle='--', alpha=0.7)

    # Markeer pieken
    if len(peaks) > 0:
        ax.plot(E_bins[peaks], residuals[peaks], 'rv', markersize=8, label=f'Pieken > 2σ (n={len(peaks)})')

    # Markeer verwachte harmonischen
    for n in range(1, 7):
        E_harm = n * THRESHOLD
        if E_harm < 8.0:
            ax.axvline(E_harm, color='red', alpha=0.2, linestyle=':', linewidth=1)

    ax.set_xlabel('Fotonenergie (MeV)')
    ax.set_ylabel('Residueel (σ)')
    ax.set_title('C) Residuelen: (Data − QED) / √QED')
    ax.legend(fontsize=8)
    ax.set_ylim(-5, 5)
    ax.grid(True, alpha=0.3)

    # ─── Panel D: Periodiciteitsscan ──────────────────────────
    ax = axes[1, 1]

    mask = expected_qed > 5
    test_periods, power, best_period = scan_periodicity(E_bins, residuals, mask)

    ax.plot(test_periods, power, 'k-', linewidth=1)
    ax.axvline(THRESHOLD, color='red', linewidth=2, alpha=0.7,
               label=f'Verwacht: 2m_e·c² = {THRESHOLD:.3f} MeV')
    ax.axvline(best_period, color='green', linewidth=1.5, linestyle='--',
               label=f'Sterkste: {best_period:.3f} MeV')

    # Markeer sub-harmonischen
    for n in [2, 3]:
        ax.axvline(n * THRESHOLD, color='red', alpha=0.3, linestyle=':')

    ax.set_xlabel('Periodiciteit (MeV)')
    ax.set_ylabel('Power')
    ax.set_title('D) Periodiciteitsscan van Residuelen')
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)

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

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

    return filepath


# ─── Significantie-berekening ────────────────────────────────
def estimate_required_statistics(amplitude_ratio=0.001, target_sigma=3):
    """
    Schat het minimaal aantal events nodig om een grid-signaal
    van gegeven sterkte te detecteren met target significantie.

    Voor een fractional signal f in een bin met N events:
        significantie ≈ f × √N

    Dus voor 3σ detectie van 0.1% signaal:
        N = (3 / 0.001)² = 9,000,000 events per bin
    """
    N_per_bin = (target_sigma / amplitude_ratio)**2
    n_bins = 200  # typisch aantal bins
    N_total = N_per_bin * n_bins

    print(f"\n{'='*60}")
    print(f"Vereiste Statistiek voor Grid-Signaal Detectie")
    print(f"{'='*60}")
    print(f"Signalamplitude:    {amplitude_ratio*100:.2f}% van QED")
    print(f"Target significantie: {target_sigma}σ")
    print(f"Events per bin:     {N_per_bin:,.0f}")
    print(f"Totaal events:      {N_total:,.0f}")
    print(f"{'='*60}")

    if amplitude_ratio >= 0.01:
        print("→ Haalbaar met bestaande CERN data")
    elif amplitude_ratio >= 0.001:
        print("→ Vereist high-luminosity runs of gecombineerde datasets")
    else:
        print("→ Vereist toekomstige experimenten (FCC of soortgelijk)")

    return N_total


# ─── Main ────────────────────────────────────────────────────
if __name__ == "__main__":
    output_dir = os.path.dirname(os.path.abspath(__file__))

    print("Voorspelling #17: Helix-Stabiliteitsdrempel bij Paarproductie")
    print("=" * 65)

    # Genereer analyse-figuur
    filepath = create_analysis_figure(output_dir)

    # Statistiek-schatting
    for amp in [0.01, 0.001, 0.0001]:
        estimate_required_statistics(amplitude_ratio=amp)

    print(f"\n✓ Analyse voltooid. Figuur: {filepath}")