<jupyter_start><jupyter_text>e-ASTROGAM's effective area<jupyter_code>from hazma.gamma_ray_limits.gamma_ray_limit_parameters import A_eff_e_ASTROGAM

e_gams = np.logspace(0, np.log10(500), 500)

plt.loglog(e_gams, [A_eff_e_ASTROGAM(e) for e in e_gams], label="Actual")
plt.loglog(e_gams[[0, -1]], [500, 500], label="Estimate")

plt.xlabel(r"$E_\gamma$ (MeV)")
plt.ylabel(r"$A_{\text{eff}}$ (cm$^2$)")
plt.legend()<jupyter_output><empty_output><jupyter_text>DM annihilation spectra<jupyter_code>e_gams = np.logspace(0, np.log10(600), 500)

# DM spectra
for mx in [110., 144., 285., 550., 770.]:
    sm.mx = mx
    spectra = sm.spectra(e_gams, 2.01*mx)
    plt.loglog(e_gams, e_gams * spectra["total"], label=r"$m_X = %.0f$ MeV" % mx)

# Background spectrum
plt.loglog(gamma_ray_limit_parameters.e_Bs, gamma_ray_limit_parameters.e_Bs * gamma_ray_limit_parameters.dPhi_dEdOmega_Bs, label="Background")

plt.xlim(e_gams[[0, -1]])
plt.legend()<jupyter_output><empty_output><jupyter_text>Unoptimized energy window<jupyter_code>mx_min = 100.
mx_max = 1000.
n_mxs = 100

mxs = np.logspace(np.log10(mx_min), np.log10(mx_max), n_mxs)

# Compute limits
sv_lims = sm.compute_limits(mxs)
plt.plot(mxs, sv_lims)

plt.grid()
plt.xlim(mxs[[0, -1]])
plt.ylim([1e-29, 1e-26])
plt.yscale("log")
plt.xlabel(r"$m_X$ (MeV)")
plt.ylabel(r"$\langle \sigma v \rangle$ (cm$^3$/s)")
plt.title("Projected e-ASTROGAM limits")<jupyter_output><empty_output><jupyter_text>Optimizing energy window<jupyter_code>from scipy.optimize import fsolve, minimize
from scipy.interpolate import interp1d
from scipy.integrate import quad

from hazma.gamma_ray_limits import gamma_ray_limit_parameters
from hazma.parameters import neutral_pion_mass as mpi0<jupyter_output><empty_output><jupyter_text>Optimizing lower bound only<jupyter_code>mx_min = 140.
mx_max = 1000.
n_mxs = 10

mxs = np.logspace(np.log10(mx_min), np.log10(mx_max), n_mxs)
optimal_limits = []
suboptimal_limits = []

for mx in mxs:
    sm.mx = mx
    spectrum = sm.spectra(e_gams, 2.01*mx)["total"]
    
    e_gams = np.logspace(0, np.log10(2000), 500)
    # Interpolate DM spectrum
    dnde_dm = interp1d(e_gams, spectrum)

    e_b = mx / 2. * (1. + np.sqrt(1. - mpi0**2 / mx**2))  # maximum photon energy from pi0's decay
    e_a_0 = mx / 2. * (1. - np.sqrt(1. - mpi0**2 / mx**2))  # minimum photon energy from pi0's decay

    def I_S(a, b=e_b):
        return quad(lambda e: gamma_ray_limit_parameters.A_eff_e_ASTROGAM(e) * dnde_dm(e), a, b)[0]

    def I_B(a, b=e_b):
        return quad(lambda e: gamma_ray_limit_parameters.A_eff_e_ASTROGAM(e) * gamma_ray_limit_parameters.dPhi_dEdOmega_B_default(e), a, b)[0]

    def to_optimize(a):
        return dnde_dm(a) - 0.5 * I_S(a) / I_B(a) * gamma_ray_limit_parameters.dPhi_dEdOmega_B_default(a)

    # Find optimal lower bound for energy window
    e_a_optimal = fsolve(to_optimize, e_a_0, xtol=1e-4)
    
    optimal_limits.append(compute_limits.compute_limit(dnde_dm, mx, e_a_optimal, e_b))
    suboptimal_limits.append(compute_limits.compute_limit(dnde_dm, mx, mpi0 / 2., mx))
    
optimal_limits = np.array(optimal_limits)
suboptimal_limits = np.array(suboptimal_limits)
plt.plot(mxs, optimal_limits, label="Optimal")
plt.plot(mxs, suboptimal_limits, label="Suboptimal")

plt.grid()
plt.xlim(mxs[[0, -1]])
plt.ylim([1e-29, 1e-26])
plt.yscale("log")
plt.xlabel(r"$m_X$ (MeV)")
plt.ylabel(r"$\langle \sigma v \rangle$ (cm$^3$/s)")
plt.legend()<jupyter_output><empty_output><jupyter_text>Optimizing upper and lower bounds<jupyter_code>mx_min = 144.
mx_max = 1000.
n_mxs = 10
mxs = np.logspace(np.log10(mx_min), np.log10(mx_max), n_mxs)

optimal_limits = []
suboptimal_limits = []

for mx in mxs:
    sm.mx = mx
    spectrum = sm.spectra(e_gams, 2.01*mx)["total"]
    
    e_gams = np.logspace(0, np.log10(2000), 500)
    # Interpolate DM spectrum
    dnde_dm = interp1d(e_gams, spectrum)

    e_b_0 = mx / 2. * (1. + np.sqrt(1. - mpi0**2 / mx**2))  # maximum photon energy from pi0's decay
    e_a_0 = mx / 2. * (1. - np.sqrt(1. - mpi0**2 / mx**2))  # minimum photon energy from pi0's decay

    def I_S(a, b):
        return quad(lambda e: gamma_ray_limit_parameters.A_eff_e_ASTROGAM(e) * dnde_dm(e), a, b)[0]

    def I_B(a, b):
        return quad(lambda e: gamma_ray_limit_parameters.A_eff_e_ASTROGAM(e) * gamma_ray_limit_parameters.dPhi_dEdOmega_B_default(e), a, b)[0]

    def to_optimize(ab):
        return -I_S(ab[0], ab[1]) / np.sqrt(I_B(ab[0], ab[1]))

    # Find optimal energy window
    e_window = minimize(to_optimize, [e_a_0, e_b_0], tol=1e-4, bounds=[[e_gams[0], mpi0 / 2.0001], [mpi0 / 1.9999, mx]]).x
    
    if e_window[0] < e_window[1]:
        e_a_optimal = e_window[0]
        e_b_optimal = e_window[1]
    else:
        e_a_optimal = e_window[1]
        e_b_optimal = e_window[0]
        
    print("m_X: ", mx, " MeV")
    print("Optimal [E_a, E_b]: ", e_window)
    print("Suboptimal [E_a, E_b]: ", (e_a_0, e_b_0))
    print("Stefano's [E_a, E_b]: ", (mpi0 / 2., mx))
    print("")
    
    optimal_limits.append(compute_limits.compute_limit(dnde_dm, mx, e_a_optimal, e_b_optimal))
    suboptimal_limits.append(compute_limits.compute_limit(dnde_dm, mx, mpi0 / 2., mx))
    
optimal_limits = np.array(optimal_limits)
suboptimal_limits = np.array(suboptimal_limits)
plt.plot(mxs, optimal_limits, label="Optimal")
plt.plot(mxs, suboptimal_limits, label="Suboptimal")

plt.grid()
plt.xlim(mxs[[0, -1]])
plt.ylim([1e-29, 1e-26])
plt.yscale("log")
plt.xlabel(r"$m_X$ (MeV)")
plt.ylabel(r"$\langle \sigma v \rangle$ (cm$^3$/s)")
plt.legend()
mx_min = 144.
mx_max = 1000.
n_mxs = 5
mxs = np.logspace(np.log10(mx_min), np.log10(mx_max), n_mxs)

optimal_limits = []
suboptimal_limits = []

for mx in mxs:
    sm.mx = mx
    spectrum = sm.spectra(e_gams, 2.01*mx)["total"]
    
    e_gams = np.logspace(0, np.log10(2000), 500)
    # Interpolate DM spectrum
    dnde_dm = interp1d(e_gams, spectrum)

    e_b_0 = mx / 2. * (1. + np.sqrt(1. - mpi0**2 / mx**2))  # maximum photon energy from pi0's decay
    e_a_0 = mx / 2. * (1. - np.sqrt(1. - mpi0**2 / mx**2))  # minimum photon energy from pi0's decay

    def I_S(a, b):
        return quad(lambda e: gamma_ray_limit_parameters.A_eff_e_ASTROGAM(e) * dnde_dm(e), a, b)[0]

    def I_B(a, b):
        return quad(lambda e: gamma_ray_limit_parameters.A_eff_e_ASTROGAM(e) * gamma_ray_limit_parameters.dPhi_dEdOmega_B_default(e), a, b)[0]

    def to_optimize(ab):
#         return compute_limits.compute_limit(dnde_dm, sm.mx, ab[0], ab[1])
        return -I_S(ab[0], ab[1]) / np.sqrt(I_B(ab[0], ab[1]))

    # Find optimal energy window
    e_window = minimize(to_optimize, [e_a_0, e_b_0], bounds=[[e_gams[0], mpi0 / 2.0001], [mpi0 / 1.9999, mx]])
    e_a_optimal = e_window.x[0]
    e_b_optimal = e_window.x[1]
        
    print("m_X: ", mx, " MeV")
    print("Optimal [E_a, E_b]: ", e_window.x)
    print("Stefano's [E_a, E_b]: ", (mpi0 / 2., mx))
    print("")
    
#     optimal_limits.append(e_window.fun)
    optimal_limits.append(compute_limits.compute_limit(dnde_dm, mx, e_a_optimal, e_b_optimal))
    suboptimal_limits.append(compute_limits.compute_limit(dnde_dm, mx, mpi0 / 2., mx))
    
optimal_limits = np.array(optimal_limits)
suboptimal_limits = np.array(suboptimal_limits)
plt.plot(mxs, optimal_limits, label="Optimal")
plt.plot(mxs, suboptimal_limits, label="Suboptimal")

plt.grid()
plt.xlim(mxs[[0, -1]])
plt.ylim([1e-29, 1e-26])
plt.yscale("log")
plt.xlabel(r"$m_X$ (MeV)")
plt.ylabel(r"$\langle \sigma v \rangle$ (cm$^3$/s)")
plt.legend()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1607119139327932, 'fraction_numerical': 0.03984592907424625, 'mean_word_length': 4.291637385804638, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '662031', 'license': 'MIT', 'n_tokens_mistral': 3525, 'n_tokens_neox': 3270, 'n_words': 572, 'path': 'notebooks/dev/Gamma_ray_limits_testing_4_17_18.ipynb', 'repo_name': 'LoganAMorrison/Hazma', 'size': '7529'}