Skip to content
Snippets Groups Projects
Select Git revision
  • bbc9597e99e62a3eb710ef029f755bc9c3d6d153
  • main default protected
  • Coverage
  • gitlab-ci-docker
  • test-install-requirementsshell-executor
  • 1.0
6 results

VisPoissonBoltzmann.py

Blame
  • Forked from Jan Habscheid / fxdgm
    Source project has a limited visibility.
    Code owners
    Assign users and groups as approvers for specific file changes. Learn more.
    fused_analysis.py 23.31 KiB
    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()