fused_analysis.py

import numpy as np
import matplotlib.pyplot as plt
import scipy.odr as sodr
import scipy.stats as scs
from my_methods import *
import scipy.optimize
import scipy.constants as scc

r_detector_ring = 8.1/2
sigma_r_detector_ring = 0.01/np.sqrt(12) /2
r_detector_conv = 2.5/2
sigma_r_detector_conv = 0.01/np.sqrt(12) /2

epsilon_m = -1.37e-04
epsilon_b = 5.48e-01
sigma_epsilon_m = 0.18e-04
sigma_epsilon_b = 0.17e-01

rho_al = 2.699 #g/cm3
rho_air =  1.205e-03 #g/cm3
rho_fe = 7.874
N_perV_al = rho_al/13 * scc.N_A # n electrons per cm3
N_perV_fe = rho_fe/26 * scc.N_A # n electrons per cm3
F_d_ring = np.pi * (r_detector_ring)**2
F_d_conv = np.pi * (r_detector_conv)**2
sigma_F_d_ring = 2 * np.pi * r_detector_ring * sigma_r_detector_ring # syst
sigma_F_d_conv = 2 * np.pi * r_detector_conv * sigma_r_detector_conv # syst.
I = 0.851
sigma_I = 0.002
A = 12232e03
sigma_A = 18e03
T = 1200

E_gamma_0_theoretical = 661.657 # keV

eta_al_200keV = 1.223e-01 * rho_al
eta_al_800keV = 6.841e-02 * rho_al
eta_air_200keV = 1.233e-01 * rho_air
eta_air_800keV = 7.074e-02 * rho_air
eta_fe_200keV = 1.460e-01 * rho_fe
eta_fe_800keV = 6.699e-02 * rho_fe

'''
print(eta_al_200keV)
print(eta_al_800keV)
print(eta_fe_200keV)
print(eta_fe_800keV)
print(eta_air_200keV)
print(eta_air_800keV)
'''

sigma_x_mess = 0.1/np.sqrt(12)
sigma_R = np.sqrt(2) * 0.005/np.sqrt(12)
sigma_a_1 = 0.01/np.sqrt(12)

x_0 = 50
x_2_0 = 66.6
x_2_1 = 121.5
x_2_2 = x_2_1
x_4_0 = 113.6
x_4_1 = 168.5
x_4_2 = 211.5
a_1_0 = 3.36
a_1_1 = 3.5
a_1_2 = a_1_1
R_k = 6.745
R_m = 9.1975
R_g = 11.69

d_alfe_conv = 1.1 # cm
sigma_d_alfe_conv = 0.01 / np.sqrt(12)
height_alfe = 2.11 # cm
sigma_height_alfe = 0.01 / np.sqrt(12)
d_1_air_conv = 5.3  - d_alfe_conv/2 # cm
sigma_d_1_air_conv = 0.1 / np.sqrt(12)
d_2_air_conv = 37.5  - d_alfe_conv/2 # cm
sigma_d_2_air_conv = 0.1 / np.sqrt(12)

data_ring = parseCSV('./output/ring.log')
thetas_ring = data_ring[:, 0]
energies_ring = data_ring[:, 1]
sigma_thetas_ring = data_ring[:, 2]
sigma_energies_ring = data_ring[:, 3]
sigma_energies_syst_ring = data_ring[:, 4]
peak_count_ring = data_ring[:, 5]
sigma_peak_count_ring = data_ring[:, 6]
sigma_peak_count_syst_ring = data_ring[:, 7]

sigma_energies_ring[1] = 5 #keV
sigma_energies_ring[2] = 5 #keV

data_conv = parseCSV('./output/conv.log')
thetas_conv = data_conv[:, 0]
energies_conv = data_conv[:, 1]
sigma_thetas_conv = data_conv[:, 2]
sigma_energies_conv = data_conv[:, 3]
sigma_energies_syst_conv = data_conv[:, 4]
peak_count_conv = data_conv[:, 5]
sigma_peak_count_conv = data_conv[:, 6]
sigma_peak_count_syst_conv = data_conv[:, 7]

data_conv_fe = parseCSV('./output/conv_fe.log')
thetas_conv_fe = data_conv_fe[:, 0]
energies_conv_fe = data_conv_fe[:, 1]
sigma_thetas_conv_fe = data_conv_fe[:, 2]
sigma_energies_conv_fe = data_conv_fe[:, 3]
sigma_energies_syst_conv_fe = data_conv_fe[:, 4]
peak_count_conv_fe = data_conv_fe[:, 5]
sigma_peak_count_conv_fe = data_conv_fe[:, 6]
sigma_peak_count_syst_conv_fe = data_conv_fe[:, 7]

thetas = np.concatenate((thetas_ring, thetas_conv))
energies = np.concatenate((energies_ring, energies_conv))
sigma_thetas = np.concatenate((sigma_thetas_ring, sigma_thetas_conv))
sigma_energies = np.concatenate((sigma_energies_ring, sigma_energies_conv))
sigma_energies_syst = np.concatenate((sigma_energies_syst_ring, sigma_energies_syst_conv))
thetas_sorted = np.sort(thetas)

print('##### Compton-Energy')

def compton(E_gamma_0, theta):
    E_gamma_0 = E_gamma_0[0]
    a = E_gamma_0 / (511)
    return E_gamma_0 / (1 + a * (1 - np.cos(theta * np.pi / 180)))

model = sodr.Model(compton)
beta0 = [E_gamma_0_theoretical]
data = sodr.RealData(thetas, energies, sx = sigma_thetas, sy = sigma_energies)
odr = sodr.ODR(data, model, beta0 = beta0)
output = odr.run()
data = sodr.RealData(thetas, energies + sigma_energies_syst, sx = sigma_thetas, sy = sigma_energies)
odr = sodr.ODR(data, model, beta0 = beta0)
output_a = odr.run()
data = sodr.RealData(thetas, energies - sigma_energies_syst, sx = sigma_thetas, sy = sigma_energies)
odr = sodr.ODR(data, model, beta0 = beta0)
output_b = odr.run()
sigma_E_gamma_0_syst = np.abs(output_b.beta[0] - output_a.beta[0]) / 2

print(f'{round(output.beta[0], np.sqrt(output.cov_beta[0, 0]))} (stat.) +- {np.around(sigma_E_gamma_0_syst, 2)} (syst.) keV')

errorscale = 10
plt.errorbar(thetas_ring, energies_ring, xerr = sigma_thetas_ring * errorscale, yerr = sigma_energies_ring * errorscale , label = r'Ring-Geometry', linestyle = '', marker = '.')
plt.errorbar(thetas_conv, energies_conv, xerr = sigma_thetas_conv * errorscale, yerr = sigma_energies_conv * errorscale , label = r'Conventional-Geometry', linestyle = '', marker = '.')
plt.plot(thetas_sorted, compton([E_gamma_0_theoretical], thetas_sorted), label = 'Theoretical')
plt.plot(thetas_sorted, compton(output.beta, thetas_sorted), label = 'Fitted')
plt.ylabel(r'$Energy \, [keV]$')
plt.xlabel(r'$\theta$ [$\degree$]')
plt.title(r'Compton-Energy by $\theta$')
plt.grid()
plt.legend()
plt.tight_layout()
plt.savefig('output/compton_energy.pdf')
plt.clf()

print('##### Cross-Section')
print('### Ring')

print('# N_e')
N_e_ring = np.zeros((5), dtype = 'float64')
sigma_N_e_ring = np.zeros((5), dtype = 'float64')
V_k = (np.pi * R_k**2 * 1.5 - np.pi*(R_k - 1.5)**2 * 1.5)
V_m = (np.pi * R_m**2 * 1.5 - np.pi*(R_m - 1.5)**2 * 1.5)
V_g = (np.pi * R_g**2 * 1.5 - np.pi*(R_g - 1.5)**2 * 1.5)
sigma_V_k = np.pi * np.sqrt((2*np.pi*R_k*1.5 - 2*(R_k-1.5)*1.5)**2 * (0.01/np.sqrt(12))**2 + (np.pi*R_k**2 - np.pi *(R_k-1.5)**2)**2 * (0.01/np.sqrt(12))**2 + (2*np.pi*(R_k-1.5)*1.5)**2 * (0.01/np.sqrt(12))**2)
sigma_V_m = np.pi * np.sqrt((2*np.pi*R_m*1.5 - 2*(R_m-1.5)*1.5)**2 * (0.01/np.sqrt(12))**2 + (np.pi*R_m**2 - np.pi *(R_m-1.5)**2)**2 * (0.01/np.sqrt(12))**2 + (2*np.pi*(R_m-1.5)*1.5)**2 * (0.01/np.sqrt(12))**2)
sigma_V_g = np.pi * np.sqrt((2*np.pi*R_g*1.5 - 2*(R_g-1.5)*1.5)**2 * (0.01/np.sqrt(12))**2 + (np.pi*R_g**2 - np.pi *(R_g-1.5)**2)**2 * (0.01/np.sqrt(12))**2 + (2*np.pi*(R_g-1.5)*1.5)**2 * (0.01/np.sqrt(12))**2)
N_e_ring[0] = N_perV_al * V_g
N_e_ring[1] = N_perV_al * V_m
N_e_ring[2] = N_perV_al * V_k
N_e_ring[3] = N_perV_al * V_m
N_e_ring[4] = N_perV_al * V_k

sigma_N_e_ring[0] = N_perV_al * sigma_V_g
sigma_N_e_ring[1] = N_perV_al * sigma_V_m
sigma_N_e_ring[2] = N_perV_al * sigma_V_k
sigma_N_e_ring[3] = N_perV_al * sigma_V_m
sigma_N_e_ring[4] = N_perV_al * sigma_V_k

'''
print(round(V_g, sigma_V_g))
print(round(V_m, sigma_V_m))
print(round(V_k, sigma_V_k))
print(round(N_e_ring[0], sigma_N_e_ring[0], scientific = True))
print(round(N_e_ring[1], sigma_N_e_ring[1], scientific = True))
print(round(N_e_ring[2], sigma_N_e_ring[2], scientific = True))
'''

print('# distances')
def distances(x_2, x_4, a_1, R):
    a = (x_4 - x_2) - a_1 - 4.33 + 2.5/2
    sigma_a = np.sqrt(2*sigma_x_mess**2 + sigma_a_1**2 + (0.01/np.sqrt(12))**2 + (0.01/2 /np.sqrt(12))**2)
    b = (x_2 - x_0) - 2.5 - 2.06/2 + 4.33
    sigma_b = np.sqrt(2*sigma_x_mess**2 + 3*(0.01/np.sqrt(12))**2 + 3*(0.01/2 /np.sqrt(12))**2)
    c = R
    sigma_c = sigma_R

    r_0 = np.sqrt(a**2 + c**2)
    sigma_r_0 = np.sqrt((a*sigma_a)**2 + (c*sigma_c)**2) / r_0

    r = np.sqrt(b**2 + c**2)
    sigma_r = np.sqrt((b*sigma_b)**2 + (c*sigma_c)**2) / r_0

    return (r_0, sigma_r_0, r, sigma_r)

r_0_ring = np.zeros((5), dtype = 'float64')
sigma_r_0_ring = np.zeros((5), dtype = 'float64')
r_ring = np.zeros((5), dtype = 'float64')
sigma_r_ring = np.zeros((5), dtype = 'float64')

r_0_ring[0], sigma_r_0_ring[0], r_ring[0], sigma_r_ring[0] = distances(x_2_0, x_4_0, a_1_0, R_g)
r_0_ring[1], sigma_r_0_ring[1], r_ring[1], sigma_r_ring[1] = distances(x_2_0, x_4_0, a_1_0, R_m)
r_0_ring[2], sigma_r_0_ring[2], r_ring[2], sigma_r_ring[2] = distances(x_2_0, x_4_0, a_1_0, R_k)
r_0_ring[3], sigma_r_0_ring[3], r_ring[3], sigma_r_ring[3] = distances(x_2_1, x_4_1, a_1_1, R_m)
r_0_ring[4], sigma_r_0_ring[4], r_ring[4], sigma_r_ring[4] = distances(x_2_2, x_4_2, a_1_2, R_k)

'''
for i in range(len(r_ring)):
    print(f'Ring & $ {round(r_0_ring[i], sigma_r_0_ring[i])} $ & $ {round(r_ring[i], sigma_r_ring[i])} $ \\\\')
'''

print('# absorption')
def eta_air(E, sigma_E):
    return ((eta_air_800keV - eta_air_200keV)/600 * (E - 200) + eta_air_200keV, (eta_air_800keV - eta_air_200keV)/600 * sigma_E)
def eta_al(E, sigma_E):
    return ((eta_al_800keV - eta_al_200keV)/600 * (E - 200) + eta_al_200keV, (eta_al_800keV - eta_al_200keV)/600 * sigma_E)

eta_temp, sigma_eta_temp = eta_air(E_gamma_0_theoretical, 0)
part_absorption_ring_0 = np.exp(- eta_temp * r_0_ring)
sigma_part_absorption_ring_0 = np.sqrt((eta_temp * sigma_r_0_ring * part_absorption_ring_0)**2 + (r_0_ring * sigma_eta_temp * part_absorption_ring_0)**2)

eta_temp, sigma_eta_temp = eta_al(E_gamma_0_theoretical, 0)
part_absorption_ring_1 = np.exp(- eta_temp * 1.5/2)
sigma_part_absorption_ring_1 = np.sqrt((eta_temp * sigma_a_1/2 * part_absorption_ring_1)**2 + (1.5/2 * sigma_eta_temp * part_absorption_ring_1)**2)

eta_temp, sigma_eta_temp = eta_air(energies_ring, sigma_energies_ring)
part_absorption_ring_2 = np.exp(- eta_temp * r_ring)
sigma_part_absorption_ring_2 = np.sqrt((eta_temp * sigma_r_ring * part_absorption_ring_2)**2 + (r_ring * sigma_eta_temp * part_absorption_ring_2)**2)

eta_temp, sigma_eta_temp = eta_al(energies_ring, sigma_energies_ring)
part_absorption_ring_3 = np.exp(- eta_temp * 1.5/2)
sigma_part_absorption_ring_3 = np.sqrt((eta_temp * sigma_a_1 * part_absorption_ring_3)**2 + (1.5/2 * sigma_eta_temp * part_absorption_ring_3)**2)

total_absorption_ring = part_absorption_ring_0 * part_absorption_ring_1 * part_absorption_ring_2 * part_absorption_ring_3
sigma_total_absorption_ring = np.sqrt((total_absorption_ring * sigma_part_absorption_ring_0 / part_absorption_ring_0)**2 + (total_absorption_ring * sigma_part_absorption_ring_1 / part_absorption_ring_1)**2 + (total_absorption_ring * sigma_part_absorption_ring_2 / part_absorption_ring_2)**2 + (total_absorption_ring * sigma_part_absorption_ring_3 / part_absorption_ring_3)**2)

'''
for i in range(len(total_absorption_ring)):
    print(f'Ring Al & abc & $ {round(total_absorption_ring[i], sigma_total_absorption_ring[i])} $ \\\\')
'''

print('# m')
m_ring = peak_count_ring / T
sigma_m_ring_stat = sigma_peak_count_ring / T
sigma_m_ring_syst = sigma_peak_count_syst_ring / T

'''
for i in range(len(m_ring)):
    print(f'Ring. Al & abc & $ {round(m_ring[i], sigma_m_ring_stat[i])} \\pm {np.around(sigma_m_ring_syst[i], 2)} $ \\\\')
'''

print('# epsilon')
epsilon_ring = epsilon_m * energies_ring + epsilon_b
sigma_epsilon_ring_stat = epsilon_m * sigma_energies_ring
sigma_epsilon_ring_syst = np.sqrt((epsilon_m * sigma_energies_syst_ring)**2 + (sigma_epsilon_m * energies_ring)**2 + (sigma_epsilon_b)**2)

print('# calculate cs')
cross_section_ring = m_ring * r_ring**2 * r_0_ring**2 * 4*np.pi /(A * I * total_absorption_ring * epsilon_ring * N_e_ring * F_d_ring)
sigma_cross_section_ring_stat = cross_section_ring * np.sqrt((sigma_m_ring_stat / m_ring)**2 + (sigma_r_ring / r_ring)**2 + (sigma_r_0_ring / r_0_ring)**2 + (sigma_total_absorption_ring / total_absorption_ring)**2 + (sigma_N_e_ring / N_e_ring)**2 + (sigma_epsilon_ring_stat / epsilon_ring)**2)
sigma_cross_section_ring_syst = cross_section_ring * np.sqrt((sigma_I / I)**2 + (sigma_A / A)**2 + (sigma_m_ring_syst / m_ring)**2 + (sigma_F_d_ring / F_d_ring)**2 + (sigma_epsilon_ring_syst / epsilon_ring)**2)

for i in range(len(cross_section_ring)):
    print(f'Ring / Al &  & $ {round(cross_section_ring[i] * 1e27, sigma_cross_section_ring_stat[i] * 1e27)} \\pm {np.around(sigma_cross_section_ring_syst[i] * 1e27, 2)} $ \\\\')
    
print('### Conventional')
print('# N_e')
V_conv = np.pi * d_alfe_conv**2 * height_alfe / 4
sigma_V_conv = V_conv * np.sqrt((sigma_d_alfe_conv / d_alfe_conv)**2 + (sigma_height_alfe / height_alfe)**2)

N_e_al_conv = N_perV_al * V_conv
N_e_fe_conv = N_perV_fe * V_conv
sigma_N_e_al_conv = N_perV_al * sigma_V_conv
sigma_N_e_fe_conv = N_perV_fe * sigma_V_conv

'''
print(round(V_conv, sigma_V_conv))
print(round(N_e_al_conv, sigma_N_e_al_conv, scientific = True))
print(round(N_e_fe_conv, sigma_N_e_fe_conv, scientific = True))
'''

print('# distances')
r_0_conv = d_1_air_conv + d_alfe_conv/2
sigma_r_0_conv = np.sqrt(sigma_d_1_air_conv**2 + (sigma_d_alfe_conv/2)**2)

r_conv = d_2_air_conv + d_alfe_conv/2
sigma_r_conv = np.sqrt(sigma_d_2_air_conv**2 + (sigma_d_alfe_conv/2)**2)

'''
print(f'Conv. & $ {round(r_0_conv, sigma_r_0_conv)} $ & $ {round(r_conv, sigma_r_conv)} $ \\\\')
'''

print('# absorption')
def eta_fe(E, sigma_E):
    return ((eta_fe_800keV - eta_fe_200keV)/600 * (E - 200) + eta_fe_200keV, (eta_fe_800keV - eta_fe_200keV)/600 * sigma_E)

# Al
eta_temp, sigma_eta_temp = eta_air(E_gamma_0_theoretical, 0)
part_absorption_conv_al_0 = np.exp(- eta_temp * r_0_conv)
sigma_part_absorption_conv_al_0 = np.sqrt((eta_temp * sigma_r_0_conv * part_absorption_conv_al_0)**2 + (r_0_conv * sigma_eta_temp * part_absorption_conv_al_0)**2)

eta_temp, sigma_eta_temp = eta_al(E_gamma_0_theoretical, 0)
part_absorption_conv_al_1 = np.exp(- eta_temp * d_alfe_conv/2)
sigma_part_absorption_conv_al_1 = np.sqrt((eta_temp * sigma_d_alfe_conv/2 * part_absorption_conv_al_1)**2 + (d_alfe_conv/2 * sigma_eta_temp * part_absorption_conv_al_1)**2)

eta_temp, sigma_eta_temp = eta_air(energies_conv, sigma_energies_conv)
part_absorption_conv_al_2 = np.exp(- eta_temp * r_conv)
sigma_part_absorption_conv_al_2 = np.sqrt((eta_temp * sigma_r_conv * part_absorption_conv_al_2)**2 + (r_conv * sigma_eta_temp * part_absorption_conv_al_2)**2)

eta_temp, sigma_eta_temp = eta_al(energies_conv, sigma_energies_conv)
part_absorption_conv_al_3 = np.exp(- eta_temp * d_alfe_conv/2)
sigma_part_absorption_conv_al_3 = np.sqrt((eta_temp * sigma_d_alfe_conv * part_absorption_conv_al_3)**2 + (d_alfe_conv/2 * sigma_eta_temp * part_absorption_conv_al_3)**2)

total_absorption_conv_al = part_absorption_conv_al_0 * part_absorption_conv_al_1 * part_absorption_conv_al_2 * part_absorption_conv_al_3
sigma_total_absorption_conv_al = np.sqrt((total_absorption_conv_al * sigma_part_absorption_conv_al_0 / part_absorption_conv_al_0)**2 + (total_absorption_conv_al * sigma_part_absorption_conv_al_1 / part_absorption_conv_al_1)**2 + (total_absorption_conv_al * sigma_part_absorption_conv_al_2 / part_absorption_conv_al_2)**2 + (total_absorption_conv_al * sigma_part_absorption_conv_al_3 / part_absorption_conv_al_3)**2)

# Fe
eta_temp, sigma_eta_temp = eta_air(E_gamma_0_theoretical, 0)
part_absorption_conv_fe_0 = np.exp(- eta_temp * r_0_conv)
sigma_part_absorption_conv_fe_0 = np.sqrt((eta_temp * sigma_r_0_conv * part_absorption_conv_fe_0)**2 + (r_0_conv * sigma_eta_temp * part_absorption_conv_fe_0)**2)

eta_temp, sigma_eta_temp = eta_fe(E_gamma_0_theoretical, 0)
part_absorption_conv_fe_1 = np.exp(- eta_temp * d_alfe_conv/2)
sigma_part_absorption_conv_fe_1 = np.sqrt((eta_temp * sigma_d_alfe_conv/2 * part_absorption_conv_fe_1)**2 + (d_alfe_conv/2 * sigma_eta_temp * part_absorption_conv_fe_1)**2)

eta_temp, sigma_eta_temp = eta_air(energies_conv_fe, sigma_energies_conv_fe)
part_absorption_conv_fe_2 = np.exp(- eta_temp * r_conv)
sigma_part_absorption_conv_fe_2 = np.sqrt((eta_temp * sigma_r_conv * part_absorption_conv_fe_2)**2 + (r_conv * sigma_eta_temp * part_absorption_conv_fe_2)**2)

eta_temp, sigma_eta_temp = eta_fe(energies_conv_fe, sigma_energies_conv_fe)
part_absorption_conv_fe_3 = np.exp(- eta_temp * d_alfe_conv/2)
sigma_part_absorption_conv_fe_3 = np.sqrt((eta_temp * sigma_d_alfe_conv * part_absorption_conv_fe_3)**2 + (d_alfe_conv/2 * sigma_eta_temp * part_absorption_conv_fe_3)**2)

total_absorption_conv_fe = part_absorption_conv_fe_0 * part_absorption_conv_fe_1 * part_absorption_conv_fe_2 * part_absorption_conv_fe_3
sigma_total_absorption_conv_fe = np.sqrt((total_absorption_conv_fe * sigma_part_absorption_conv_fe_0 / part_absorption_conv_fe_0)**2 + (total_absorption_conv_fe * sigma_part_absorption_conv_fe_1 / part_absorption_conv_fe_1)**2 + (total_absorption_conv_fe * sigma_part_absorption_conv_fe_2 / part_absorption_conv_fe_2)**2 + (total_absorption_conv_fe * sigma_part_absorption_conv_fe_3 / part_absorption_conv_fe_3)**2)

'''
for i in range(len(total_absorption_conv_al)):
    print(f'Conv. Al & abc & $ {round(total_absorption_conv_al[i], sigma_total_absorption_conv_al[i])} $ \\\\')

for i in range(len(total_absorption_conv_fe)):
    print(f'Conv. Fe & abc & $ {round(total_absorption_conv_fe[i], sigma_total_absorption_conv_fe[i])} $ \\\\')
'''

print('# m')
m_conv_fe = peak_count_conv_fe / T
sigma_m_conv_fe_stat = sigma_peak_count_conv_fe / T
sigma_m_conv_fe_syst = sigma_peak_count_syst_conv_fe / T

m_conv_al = peak_count_conv / T
sigma_m_conv_al_stat = sigma_peak_count_conv / T
sigma_m_conv_al_syst = sigma_peak_count_syst_conv / T

'''
for i in range(len(m_conv_al)):
    print(f'Conv. Al & abc & $ {round(m_conv_al[i], sigma_m_conv_al_stat[i])} \\pm {np.around(sigma_m_conv_al_syst[i], 2)} $ \\\\')

for i in range(len(m_conv_fe)):
    print(f'Conv. Fe & abc & $ {round(m_conv_fe[i], sigma_m_conv_fe_stat[i])} \\pm {np.around(sigma_m_conv_fe_syst[i], 2)} $ \\\\')
'''

print('# epsilon')
epsilon_conv_fe = epsilon_m * energies_conv_fe + epsilon_b
sigma_epsilon_conv_fe_stat = epsilon_m * sigma_energies_conv_fe
sigma_epsilon_conv_fe_syst = np.sqrt((epsilon_m * sigma_energies_syst_conv_fe)**2 + (sigma_epsilon_m * energies_conv_fe)**2 + (sigma_epsilon_b)**2)

epsilon_conv_al = epsilon_m * energies_conv + epsilon_b
sigma_epsilon_conv_al_stat = epsilon_m * sigma_energies_conv
sigma_epsilon_conv_al_syst = np.sqrt((epsilon_m * sigma_energies_syst_conv)**2 + (sigma_epsilon_m * energies_conv)**2 + (sigma_epsilon_b)**2)

print('# calculate cs')
cross_section_conv_fe = m_conv_fe * r_conv**2 * r_0_conv**2 * 4*np.pi /(A * I * total_absorption_conv_fe * epsilon_conv_fe * N_e_fe_conv * F_d_conv)
sigma_cross_section_conv_fe_stat = cross_section_conv_fe * np.sqrt((sigma_m_conv_fe_stat / m_conv_fe)**2 + (sigma_r_conv / r_conv)**2 + (sigma_r_0_conv / r_0_conv)**2 + (sigma_total_absorption_conv_fe / total_absorption_conv_fe)**2 + (sigma_N_e_fe_conv / N_e_fe_conv)**2 + (sigma_epsilon_conv_fe_stat / epsilon_conv_fe)**2)
sigma_cross_section_conv_fe_syst = cross_section_conv_fe * np.sqrt((sigma_I / I)**2 + (sigma_A / A)**2 + (sigma_m_conv_fe_syst / m_conv_fe)**2 + (sigma_F_d_conv / F_d_conv)**2 + (sigma_epsilon_conv_fe_syst / epsilon_conv_fe)**2)

cross_section_conv_al = m_conv_al * r_conv**2 * r_0_conv**2 * 4*np.pi /(A * I * total_absorption_conv_al * epsilon_conv_al * N_e_al_conv * F_d_conv)
sigma_cross_section_conv_al_stat = cross_section_conv_al * np.sqrt((sigma_m_conv_al_stat / m_conv_al)**2 + (sigma_r_conv / r_conv)**2 + (sigma_r_0_conv / r_0_conv)**2 + (sigma_total_absorption_conv_al / total_absorption_conv_al)**2 + (sigma_N_e_al_conv / N_e_al_conv)**2 + (sigma_epsilon_conv_al_stat / epsilon_conv_al)**2)
sigma_cross_section_conv_al_syst = cross_section_conv_al * np.sqrt((sigma_I / I)**2 + (sigma_A / A)**2 + (sigma_m_conv_al_syst / m_conv_al)**2 + (sigma_F_d_conv / F_d_conv)**2 + (sigma_epsilon_conv_al_syst / epsilon_conv_al)**2)

for i in range(len(cross_section_conv_al)):
    print(f'Conv. / Al &  & $ {round(cross_section_conv_al[i] * 1e27, sigma_cross_section_conv_al_stat[i] * 1e27)} \\pm {np.around(sigma_cross_section_conv_al_syst[i] * 1e27, 2)} $ \\\\')

for i in range(len(cross_section_conv_fe)):
    print(f'Conv. / Fe &  & $ {round(cross_section_conv_fe[i] * 1e27, sigma_cross_section_conv_fe_stat[i] * 1e27)} \\pm {np.around(sigma_cross_section_conv_fe_syst[i] * 1e27, 2)} $ \\\\')

print('### plot cs')
def cs_theory(theta):
    theta = theta * np.pi/180
    a = E_gamma_0_theoretical / 511.
    rho = 1 + a * (1-np.cos(theta))
    lambda_e = 2.4e-10 # cm
    cs = 1/137**2 * lambda_e**2 / (8 * np.pi**2) * 1/rho**2 * (rho + 1/rho - np.sin(theta)**2) # cm2
    cs = cs * 1e24 * 1e3 # mbarn
    return cs

plt.errorbar(thetas_ring, cross_section_ring*1e24 /13 *1e03, yerr = np.sqrt(sigma_cross_section_ring_stat**2 + sigma_cross_section_ring_syst**2)*1e24 /13 *1e03, xerr = sigma_thetas_ring, linestyle = '', label = 'differential cs data, RG, Al')
plt.errorbar(thetas_conv, cross_section_conv_al*1e24 /13 *1e03, yerr = np.sqrt(sigma_cross_section_conv_al_stat**2 + sigma_cross_section_conv_al_syst**2)*1e24 /13 *1e03, xerr = sigma_thetas_conv, linestyle = '', label = 'differential cs data, CG, Al')
plt.errorbar(thetas_conv_fe, cross_section_conv_fe*1e24 /26 *1e03, yerr = np.sqrt(sigma_cross_section_conv_fe_stat**2 + sigma_cross_section_conv_fe_syst**2)*1e24 /26 *1e03, xerr = sigma_thetas_conv, linestyle = '', label = 'differential cs data, CG, Fe')
plt.plot(np.linspace(0, 180, 100), cs_theory(np.linspace(0, 180, 100)), label = 'differential cs theory')
plt.ylabel(r'$\frac{d\sigma}{d\Omega} \, [mbarn]$')
plt.xlabel(r'$\theta$ [$\degree$]')
plt.title(r'Differential cross section by $\theta$, ring-geometry, Al')
plt.grid()
plt.legend()
plt.tight_layout()
plt.savefig('output/cs_ring.pdf')
plt.clf()

print('##### Electron mass')
energies = np.concatenate((energies, energies_conv_fe))
sigma_energies = np.concatenate((sigma_energies, sigma_energies_conv_fe))
sigma_energies_syst = np.concatenate((sigma_energies_syst, sigma_energies_syst_conv_fe))
thetas = np.concatenate((thetas, thetas_conv_fe))
sigma_thetas = np.concatenate((sigma_thetas, sigma_thetas_conv_fe))

scale = 1e-03
# y = 1/E' - 1/E
y = 1/(energies * scale) - 1/(E_gamma_0_theoretical * scale)
sigma_y = np.abs((sigma_energies * scale)/(energies * scale)**2)
sigma_y_syst = np.abs((sigma_energies_syst * scale)/(energies * scale)**2)
# x = 1 - cos(theta)
x = 1 - np.cos(thetas * np.pi/180)
sigma_x = np.abs(np.sin(thetas * np.pi/180) * sigma_thetas)

'''
start = 0
end = 20
x = x[start:end:]
sigma_x = sigma_x[start:end:]
y = y[start:end:]
sigma_y = sigma_y[start:end:]
sigma_y_syst = sigma_y_syst[start:end:]
'''

output_a = linear_fit(x, y + sigma_y_syst, xerr = sigma_x, yerr = sigma_y)
output_b = linear_fit(x, y - sigma_y_syst, xerr = sigma_x, yerr = sigma_y)
sigma_m_syst = np.abs(output_a.beta[0] - output_b.beta[0])/2

output = linear_fit(x, y, xerr = sigma_x, yerr = sigma_y)
chiq_dof, p_value = chi_sq_dof_p(x, y, xerr = sigma_x, yerr = sigma_y, beta = output.beta, model = linear)
error_2d = calc_2d_error(x, y, xerr = sigma_x, yerr = sigma_y, beta = output.beta, model = linear)

# m = 1/(m_e*c**2)
m_e_compton = 1 / output.beta[0] / scale # keV
sigma_m_e_compton = np.sqrt(output.cov_beta[0, 0]) / output.beta[0]**2 / scale
sigma_m_e_compton_syst = sigma_m_syst / output.beta[0]**2 / scale

print(f'\t m_e: {round(m_e_compton, sigma_m_e_compton)} (stat.) +- {np.around(sigma_m_e_compton_syst, 3)} (syst.)')

errorscale = 10
fig, ax = plt.subplots(nrows = 2, sharex = True)
ax[0].plot(x, linear(output.beta, x), label = r'fit $y = m \cdot x + b$', color='C1')
ax[0].errorbar(x, y, xerr = sigma_x / errorscale, yerr = sigma_y * errorscale, label='data', linestyle = '', marker = '.')
ax[1].errorbar(x , y - linear(output.beta, x), yerr = error_2d, label = 'residuals', linestyle = '', marker = '.')
ax[0].set(ylabel = r"$\frac{1}{E_{\gamma}'} - \frac{1}{E_{\gamma}} \, [\frac{1}{MeV}]$", title = r'Fit for $m_e$')
ax[1].set(xlabel = r'$1 - cos \, \theta$', ylabel = r'Data - Fit [$\frac{1}{MeV}$]')
ax[0].text(0.985, 0.045, r'$m$'+f': {round(output.beta[0], np.sqrt(output.cov_beta[0, 0]))}(stat.)'+r' $\pm$ '+f'{np.around(sigma_m_syst, 3)}(syst.)\n'+r'$b$'+f': {round(output.beta[1], np.sqrt(output.cov_beta[1, 1]))}\n'+r'$\frac{\chi^2}{dof}$'+f': {np.round(chiq_dof, 1)}\tp: {np.around(p_value, 3)}' , horizontalalignment='right', verticalalignment='bottom', transform=ax[0].transAxes, bbox=dict(facecolor='white', edgecolor='gray'))
ax[1].set_xlim(-np.max(x)*0.02, np.max(x)*1.15)
ax[0].grid()
ax[1].grid()
ax[0].legend(loc = 'upper left')
ax[1].legend(loc = 'upper left')
fig.tight_layout()
plt.savefig('output/m_e_compton.pdf')
plt.clf()