<jupyter_start><jupyter_text>MeV RH Neutrino Dark Matter<jupyter_code>from scipy.integrate import quad
import numpy as np
import matplotlib.pyplot as plt
from collections import defaultdict
from hazma.gamma_ray_parameters import egret_diffuse, fermi_diffuse, comptel_diffuse, gecco_bg_model
from hazma.gamma_ray_parameters import A_eff_comptel, A_eff_egret, A_eff_fermi, energy_res_comptel, BackgroundModel
from hazma.gamma_ray_parameters import (
    gc_bg_model, default_bg_model
)
from hazma.gamma_ray_parameters import (
    gc_targets_optimistic, A_eff_e_astrogam, energy_res_e_astrogam, T_obs_e_astrogam, TargetParams, energy_res_adept, energy_res_amego, energy_res_mast, energy_res_pangu, energy_res_fermi, energy_res_egret
)
from hazma.gamma_ray_parameters import (
    A_eff_gecco, energy_res_gecco, energy_res_gecco_large, gc_targets, fornax_targets, m31_targets, draco_targets, A_eff_adept, A_eff_amego, A_eff_mast, A_eff_pangu
)
from hazma.rh_neutrino import RHNeutrino
from hazma.parameters import (
    neutral_pion_mass as mpi0,
    charged_pion_mass as mpi,
)<jupyter_output><empty_output><jupyter_text>Testings  Widths<jupyter_code>model = RHNeutrino(10.0, 1e-3, 'e')
mxs = np.logspace(np.log10(0.1), np.log10(1000.0), 300)
widths = {key: [] for key in model.decay_widths().keys()}

for mx in mxs:
    model.mx = mx
    for key, val in model.decay_widths().items():
        widths[key].append(val)

for key, val in model.decay_widths().items():
    widths[key] = np.array(widths[key])
plt.figure(dpi=150)

for key, val in widths.items():
    plt.plot(mxs, val, label=key)

plt.legend()
plt.yscale('log')
plt.xscale('log')
plt.ylabel(r"$\Gamma_{N} \ (\mathrm{MeV})$", fontsize=16)
plt.xlabel(r'$m_{N} \ (\mathrm{MeV})$', fontsize=16)

plt.xlim([np.min(mxs),np.max(mxs)])
plt.ylim([1e-40, 1e-13])
plt.savefig('figures/rh_neutrino_widths.pdf')<jupyter_output>findfont: Font family ['serif'] not found. Falling back to DejaVu Sans.<jupyter_text>Specta<jupyter_code>model = RHNeutrino(700.0, 1e-3, 'e', include_3body=True)
es = np.logspace(-3, np.log10(model.mx), 50)

plt.figure(dpi=150)

spectra = model.spectra(es)

for key, spec in spectra.items():
    plt.plot(es, es**2 * spec, label=key)

plt.ylabel(r"$E_{\gamma}^2 \frac{dN_{\gamma}}{dE_{\gamma}} \ (\mathrm{MeV})$", fontsize=16)
plt.xlabel(r"$E_{\gamma} \ (\mathrm{MeV})$", fontsize=16)
plt.xlim([np.min(es), np.max(es)])
plt.yscale('log')
plt.xscale('log')
plt.legend()

plt.savefig('figures/rh_neutrino_spectrum_electron.pdf')
model = RHNeutrino(100.0, 1e-3, 'e', include_3body=True)
es = np.logspace(-3, np.log10(model.mx), 50)

plt.figure(dpi=150)

spectra = model.spectra(es)

for key, spec in spectra.items():
    plt.plot(es, es**2 * spec, label=key)

plt.ylabel(r"$E_{\gamma}^2 \frac{dN_{\gamma}}{dE_{\gamma}} \ (\mathrm{MeV})$", fontsize=16)
plt.xlabel(r"$E_{\gamma} \ (\mathrm{MeV})$", fontsize=16)
plt.xlim([np.min(es), np.max(es)])
plt.ylim([1e-11, 1e-1])
plt.yscale('log')
plt.xscale('log')
plt.legend()

plt.savefig('figures/rh_neutrino_spectrum_electron_mx=100.pdf')<jupyter_output><empty_output><jupyter_text>Compute Constraints<jupyter_code># This is approximately the same as `gc_bg_model`, but does not constrain
# the input energies.
gc_bg_model_approx = BackgroundModel(
    [0, 1e5], lambda e: 7 * default_bg_model.dPhi_dEdOmega(e)
)
# Other constants
T_obs_nt = 1e6  # s. From Alex.
v_mw = 1e-3
x_kd = 1e-6  # should compute this for actual models...

# Best selection of targets
nt_targets_ann = {
    "gc_ein_1_arcmin_cone_optimistic": (
        gc_targets_optimistic["ein"]["1 arcmin cone"], gc_bg_model_approx
    ),
    "gc_nfw_1_arcmin_cone": (gc_targets["nfw"]["1 arcmin cone"], gc_bg_model_approx),
    "m31_nfw_1_arcmin_cone": (m31_targets["nfw"]["1 arcmin cone"], gecco_bg_model),
    "draco_nfw_1_arcmin_cone": (
        draco_targets["nfw"]["1 arcmin cone"], gecco_bg_model
    )
}
nt_targets_dec = {
    "gc_ein_5_deg_optimistic": (
        gc_targets_optimistic["ein"]["5 deg cone"], gc_bg_model_approx
    ),
    "gc_nfw_5_deg": (gc_targets["nfw"]["5 deg cone"], gc_bg_model_approx),
    "m31_nfw_5_deg": (m31_targets["nfw"]["5 deg cone"], gecco_bg_model),
    "draco_nfw_5_deg": (
        draco_targets["nfw"]["5 deg cone"], gecco_bg_model
    )
}

def get_constraints(model, mxs, sts):
    constraints = defaultdict(lambda: np.zeros(len(mxs)))

    for i, mx in enumerate(mxs):
        model.mx = mx

        # Constraints from existing data
        constraints["comptel"][i] = model.binned_limit(comptel_diffuse)
        constraints["egret"][i] = model.binned_limit(egret_diffuse)
        constraints["fermi"][i] = model.binned_limit(fermi_diffuse)
        
        e_window_0=None

        # Projections
        nt_targets = nt_targets_ann if model.kind == "ann" else nt_targets_dec
        for key, (target, bg_model) in nt_targets.items():
            constraints["adept_" + key][i] = model.unbinned_limit(
                A_eff_adept, energy_res_adept, T_obs_nt, target, bg_model,
            )
            constraints["amego_" + key][i] = model.unbinned_limit(
                A_eff_amego, energy_res_amego, T_obs_nt, target, bg_model,
            )
            constraints["e-astrogam_" + key][i] = model.unbinned_limit(
                A_eff_e_astrogam, energy_res_e_astrogam, T_obs_nt, target, bg_model,
            )
            constraints["gecco_" + key][i] = model.unbinned_limit(
                A_eff_gecco, energy_res_gecco, T_obs_nt, target, bg_model,
            )
            constraints["mast_" + key][i] = model.unbinned_limit(
                A_eff_mast, energy_res_mast, T_obs_nt, target, bg_model,
            )
            constraints["pangu_" + key][i] = model.unbinned_limit(
                A_eff_pangu, energy_res_pangu, T_obs_nt, target, bg_model,
            )

        # CMB (s-wave)
        # TODO: not sure how this works for decaying DM...
        if model.kind == "ann":
            constraints["cmb"][i] = model.cmb_limit(x_kd=x_kd)

    return dict(constraints)

def get_constraint_label(key):
    label = key.replace("amego_", "AMEGO ")
    label = label.replace("adept_", "AdEPT ")
    label = label.replace("comptel", "COMPTEL ")
    label = label.replace("e-astrogam_", "e-ASTROGAM ")
    label = label.replace("egret", "EGRET ")
    label = label.replace("fermi", "Fermi ")
    label = label.replace("gecco_", "GECCO ")
    label = label.replace("mast_", "MAST ")
    label = label.replace("pangu_", "PANGU ")
    
    label = label.replace("gc_nfw_1_arcmin_cone", "(GC $1'$, NFW)")
    label = label.replace("gc_nfw_5_deg", "(GC $5^\circ$, NFW)")
    label = label.replace("gc_ein_1_arcmin_cone_optimistic", "(GC $1'$, Einasto)")
    label = label.replace("gc_ein_5_deg_optimistic", "(GC $5^\circ$, Einasto)")
    label = label.replace("m31_nfw_1_arcmin_cone", "(M31 $1'$)")
    label = label.replace("m31_nfw_5_deg", "(M31 $5^\circ$)")
    label = label.replace("draco_nfw_1_arcmin_cone", "(Draco $1'$)")
    label = label.replace("draco_nfw_5_deg", "(Draco $5^\circ$)")
    
    return label

def get_fs_label(fs):
    if fs == "e e":
        return r"$e^+ e^-$"
    elif fs == "mu mu":
        return r"$\mu^+ \mu^-$"
    elif fs == "pi pi":
        return r"$\pi^+ \pi^-$"
    elif fs == "pi0 pi0":
        return r"$\pi^0 \pi^0$"
    elif fs == "pi0 g":
        return r"$\pi^0 \gamma$"
    elif fs == "g g":
        return r"$\gamma \gamma$"
    else:
        return fs

def get_formatted_fig(nrows, ncols, figsize, xlim, ylim):
    fig, axes = plt.subplots(nrows, ncols, figsize=figsize)

    for ax in axes[-1]:
        ax.set_xlabel(r"$m_\chi$ [MeV]")

    for ax in axes[:, 0]:
        ax.set_ylabel(r"$\langle \sigma v \rangle_{\bar{\chi}\chi, 0}$ [cm$^3$/s]")

    for ax in axes.flatten():
        ax.set_xlim(xlim)
        ax.set_ylim(ylim)
        ax.set_xscale("log")
        ax.set_yscale("log")
        ax.grid()

    return fig, axes
model = RHNeutrino(500.0, 1e-3, 'e')
mxs = np.logspace(np.log10(0.1), np.log10(500), 200)
constraints = get_constraints(model, mxs)
plt.figure(dpi=150)

plt.plot(mxs, constraints['comptel'],label='comptel')
plt.plot(mxs, constraints['egret'],label='egret')
plt.plot(mxs, constraints['fermi'],label='fermi')
plt.plot(mxs, constraints['gecco_draco_nfw_5_deg'],label=r'GECCO Draco')
plt.plot(mxs, constraints['gecco_gc_nfw_5_deg'], label='GECCO GC')
plt.plot(mxs, constraints['e-astrogam_draco_nfw_5_deg'],label=r'E-Astrogam Draco')
plt.plot(mxs, constraints['e-astrogam_gc_nfw_5_deg'], label='E-Astrogam GC')

plt.xlim([np.min(mxs),np.max(mxs)])
plt.ylim([1e-30,1e-14])

plt.xlabel(r'$m_{N} \ (\mathrm{MeV})$', fontsize=16)
plt.ylabel(r'$\Gamma_{N} \ (\mathrm{MeV})$', fontsize=16)

plt.legend(loc=1)
plt.yscale('log')
plt.xscale('log')

plt.savefig('figures/rh_neutrino_id_constraints_electron_width.pdf')
mev_to_sec_inv = 1.52e24 * 1e-3

plt.figure(dpi=150)

plt.plot(mxs, 1.0/(mev_to_sec_inv*constraints['comptel']),label='comptel')
plt.plot(mxs, 1.0/(mev_to_sec_inv*constraints['egret']),label='egret')
plt.plot(mxs, 1.0/(mev_to_sec_inv*constraints['fermi']),label='fermi')
plt.plot(mxs, 1.0/(mev_to_sec_inv*constraints['gecco_draco_nfw_5_deg']),label=r'GECCO Draco')
plt.plot(mxs, 1.0/(mev_to_sec_inv*constraints['gecco_gc_nfw_5_deg']), label='GECCO GC')
plt.plot(mxs, 1.0/(mev_to_sec_inv*constraints['e-astrogam_draco_nfw_5_deg']),label=r'E-Astrogam Draco')
plt.plot(mxs, 1.0/(mev_to_sec_inv*constraints['e-astrogam_gc_nfw_5_deg']), label='E-Astrogam GC')

plt.xlim([np.min(mxs),np.max(mxs)])
plt.ylim([1e-1,1e9])

plt.xlabel(r'$m_{N} \ (\mathrm{MeV})$', fontsize=16)
plt.ylabel(r'$\tau_{N} \ (\mathrm{s})$', fontsize=16)

plt.legend()
plt.yscale('log')
plt.xscale('log')

plt.savefig('figures/rh_neutrino_id_constraints_electron_lifetime.pdf')
model = RHNeutrino(500.0, 1e-3, 'mu')
mxs = np.logspace(np.log10(0.1), np.log10(500), 200)
constraints_mu = get_constraints(model, mxs)
plt.figure(dpi=150)

plt.plot(mxs, constraints_mu['comptel'],label='comptel')
plt.plot(mxs, constraints_mu['egret'],label='egret')
plt.plot(mxs, constraints_mu['fermi'],label='fermi')
plt.plot(mxs, constraints_mu['gecco_draco_nfw_5_deg'],label=r'GECCO Draco')
plt.plot(mxs, constraints_mu['gecco_gc_nfw_5_deg'], label='GECCO GC')
plt.plot(mxs, constraints_mu['e-astrogam_draco_nfw_5_deg'],label=r'E-Astrogam Draco')
plt.plot(mxs, constraints_mu['e-astrogam_gc_nfw_5_deg'], label='E-Astrogam GC')

plt.xlim([np.min(mxs),np.max(mxs)])
plt.ylim([1e-30,1e-14])

plt.xlabel(r'$m_{N} \ (\mathrm{MeV})$', fontsize=16)
plt.ylabel(r'$\Gamma_{N} \ (\mathrm{MeV})$', fontsize=16)

plt.legend(loc=1)
plt.yscale('log')
plt.xscale('log')

plt.savefig('figures/rh_neutrino_id_constraints_muon.pdf')
mev_to_sec_inv = 1.52e24 * 1e-3

plt.figure(dpi=150)

plt.plot(mxs, 1.0/(mev_to_sec_inv*constraints_mu['comptel']),label='comptel')
plt.plot(mxs, 1.0/(mev_to_sec_inv*constraints_mu['egret']),label='egret')
plt.plot(mxs, 1.0/(mev_to_sec_inv*constraints_mu['fermi']),label='fermi')
plt.plot(mxs, 1.0/(mev_to_sec_inv*constraints_mu['gecco_draco_nfw_5_deg']),label=r'GECCO Draco')
plt.plot(mxs, 1.0/(mev_to_sec_inv*constraints_mu['gecco_gc_nfw_5_deg']), label='GECCO GC')
plt.plot(mxs, 1.0/(mev_to_sec_inv*constraints_mu['e-astrogam_draco_nfw_5_deg']),label=r'E-Astrogam Draco')
plt.plot(mxs, 1.0/(mev_to_sec_inv*constraints_mu['e-astrogam_gc_nfw_5_deg']), label='E-Astrogam GC')

plt.xlim([np.min(mxs),np.max(mxs)])
plt.ylim([1e-1,1e9])

plt.xlabel(r'$m_{N} \ (\mathrm{MeV})$', fontsize=16)
plt.ylabel(r'$\tau_{N} \ (\mathrm{s})$', fontsize=16)

plt.legend()
plt.yscale('log')
plt.xscale('log')

plt.savefig('figures/rh_neutrino_id_constraints_muon_lifetime.pdf')<jupyter_output><empty_output><jupyter_text>Spectra<jupyter_code>model=RHNeutrino(500.0, 1e-3, 'e')
es = np.logspace(-3, np.log10(model.mx), 500)
spec = model.total_spectrum(es)

mxs = [150, 300, 500]

for mx in mxs:
    model.mx = mx
    spec = model.total_spectrum(es)
    plt.plot(es, spec)
plt.yscale('log')
plt.xscale('log')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17297432477492497, 'fraction_numerical': 0.026592197399133044, 'mean_word_length': 4.790057915057915, 'pattern_counts': {'":': 15, '<': 20, '<?xml version=': 0, '>': 20, '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': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '544156', 'license': 'MIT', 'n_tokens_mistral': 5375, 'n_tokens_neox': 5048, 'n_words': 745, 'path': 'notebooks/dev/RHNeutrino.ipynb', 'repo_name': 'LoganAMorrison/Hazma', 'size': '11996'}