$ python new_alpha_test03.py 

================================================================================

TESTING RAPIDITY SQUARED MODEL: Scale * [ln(γ)]^2

================================================================================

Best fit parameters:

  Scale Factor: 4.8510e-04 ± 1.1171e-06

  Power P:      FIXED at 2.0

Goodness of Fit:

  Reduced Chi-Squared: 3.8352

import numpy as np

import matplotlib.pyplot as plt

from scipy.optimize import curve_fit

# Physical Constants

c = 299792458

m_e = 9.10938e-31

e = 1.602176e-19

alpha_0 = 1/137.035999 

# Data (Same as before)

data_points = [

    (0.000511, 137.036, 0.001), 

    (5.0,      132.50,  0.60),   

    (10.0,     131.01,  0.50),   

    (20.0,     130.40,  0.40),   

    (40.0,     129.80,  0.30),   

    (58.0,     129.60,  0.20),   

    (91.1876,  127.94,  0.02),   

    (130.0,    127.60,  0.15),   

    (189.0,    127.10,  0.20),   

    (500.0,    125.70,  0.50),   

    (1000.0,   124.80,  0.90)    

]

energies_GeV = np.array([x[0] for x in data_points])

inv_alpha_val = np.array([x[1] for x in data_points])

inv_alpha_err = np.array([x[2] for x in data_points])

exp_alpha = 1 / inv_alpha_val

exp_alpha_err = inv_alpha_err / (inv_alpha_val**2)

def calculate_gamma(energy_GeV):

    energy_joules = energy_GeV * 1e9 * e

    gamma = energy_joules / (m_e * c**2)

    return gamma

# ==============================================================================

# RAPIDITY SQUARED MODEL

# ==============================================================================

def geometric_rapidity_squared(energy_GeV, scale_factor):

    """

    Hypothesis: α correction scales with Rapidity Squared

    Rapidity y ≈ ln(2*gamma) for high energy

    Model: α_eff = α_0 * (1 + scale_factor * [ln(gamma)]^2)

    """

    gamma = calculate_gamma(energy_GeV)

    gamma = np.maximum(gamma, 1.0)

    

    # FORCING POWER = 2

    geometric_contribution = scale_factor * (np.log(gamma)**2)

    

    return alpha_0 * (1 + geometric_contribution)

# ==============================================================================

# FITTING

# ==============================================================================

print("="*80)

print("TESTING RAPIDITY SQUARED MODEL: Scale * [ln(γ)]^2")

print("="*80)

# Only fitting Scale Factor now (1 parameter!)

popt, pcov = curve_fit(geometric_rapidity_squared, energies_GeV, exp_alpha, 

                       sigma=exp_alpha_err, absolute_sigma=True,

                       p0=[0.001],

                       maxfev=10000)

scale_fit = popt[0]

perr = np.sqrt(np.diag(pcov))[0]

print(f"\nBest fit parameters:")

print(f"  Scale Factor: {scale_fit:.4e} ± {perr:.4e}")

print(f"  Power P:      FIXED at 2.0")

# Calculate Goodness of Fit

residuals = (exp_alpha - geometric_rapidity_squared(energies_GeV, *popt)) / exp_alpha_err

chi_squared = np.sum(residuals**2)

dof = len(exp_alpha) - 1

reduced_chi_squared = chi_squared / dof

print(f"\nGoodness of Fit:")

print(f"  Reduced Chi-Squared: {reduced_chi_squared:.4f}")

# Plotting

plt.figure(figsize=(10, 7))

x_range = np.logspace(np.log10(0.0005), np.log10(2000), 500)

y_model = 1 / geometric_rapidity_squared(x_range, *popt)

plt.plot(x_range, y_model, 'b-', linewidth=2, label=f'Rapidity Squared (Fixed P=2)')

plt.errorbar(energies_GeV, inv_alpha_val, yerr=inv_alpha_err, 

             fmt='ko', capsize=5, label='Experimental Data')

plt.xscale('log')

plt.xlabel('Energy (GeV)')

plt.ylabel('1/α')

plt.title(f'Rapidity Squared Model (Fixed P=2)\nReduced χ² = {reduced_chi_squared:.3f}')

plt.legend()

plt.grid(True, alpha=0.3)

plt.savefig('alpha_rapidity_squared02.png', dpi=300)

plt.show()