<jupyter_start><jupyter_text>This notebook makes the main dark photon plots in our paper. The output arrays are obtained by running `cluster/firas_interface.py` on a cluster with several variations of parameters, which is done in `cluster/submit_firas.py`.<jupyter_code># Location of outputs of cluster runs
outputs_dir = "/Users/smsharma/chi2_arys/"

# Where to save plots
plots_dir = "../paper/draft-letter/plots/"
# Load a power spectrum
pspec_nonlin_matter = PowerSpectrumGridInterpolator("nonlin_matter_bj")

# Load an instance of FIRAS class
firas = FIRAS(pspec_nonlin_matter)
# Load existing constraints
eps_SRF, lim_SRF = np.transpose(pd.read_csv("../data/existing_constraints/SRF.cvs", header=None).to_numpy())
eps_SRF_dot, lim_SRF_dot = np.transpose(pd.read_csv("../data/existing_constraints/SRFdotted.cvs", header=None).to_numpy())
jup_m, jup_eps = 10 ** np.transpose(pd.read_csv("../data/existing_constraints/Jupiter.txt", sep='\t').to_numpy())<jupyter_output><empty_output><jupyter_text>Main constraints plot (Fig. 3, left)<jupyter_code>fig, ax = plt.subplots()

# Homogeneous curve

file_homo = np.load(outputs_dir + "firas_dp_scan_homo_pixie_0.npz")
m_A_ary = file_homo['m_A_ary'] / eV
limit_homo_ary = file_homo['limit']
limit_homo_ary[(limit_homo_ary > 1e-6) * (m_A_ary > 1e-13)] = savgol_filter(limit_homo_ary[(limit_homo_ary > 1e-6) * (m_A_ary > 1e-13)], 71, 3)

plt.plot(m_A_ary, limit_homo_ary, color='grey', ls='dotted', label="Homogeneous")

## Reset color cycle
plt.gca().set_prop_cycle(None)

# Lognormal with frankenspectrum curve

file_lower = np.load(outputs_dir + "firas_dp_scan_pdf_lognormal_franken_lower_cut_100.0_pixie_0_z_min_0.001_z_excise_max_20_b_1.npz")
m_A_ary = file_lower['m_A_ary'] / eV
limit_lower_ary = file_lower['limit']

file_upper = np.load(outputs_dir + "firas_dp_scan_pdf_lognormal_franken_upper_cut_100.0_pixie_0_z_min_0.001_z_excise_max_20_b_1.npz")
m_A_ary = file_upper['m_A_ary'] / eV
limit_upper_ary = file_upper['limit']

lim_LN_franken = np.max(np.transpose([limit_lower_ary, limit_upper_ary]), axis=1)

plt.fill_between(m_A_ary, lim_LN_franken, np.ones_like(limit_lower_ary), color=cols_default[0], alpha=0.05)
plt.plot(m_A_ary, lim_LN_franken, color=cols_default[0], label="Log-normal PDF")

# Analytic PDF curve

file = np.load(outputs_dir + "firas_dp_scan_pdf_log_analytic_pdf_interp_franken_upper_cut_100.0_pixie_0_z_min_0.001_z_excise_max_20.npz")
m_A_ary = file['m_A_ary'] / eV
limit_ary = file['limit']
plt.plot(m_A_ary, limit_ary, c=cols_default[1], label="Analytic PDF")

# Lognormal with frankenspectrum curve, PIXIE projection

file_lower = np.load(outputs_dir + "firas_dp_scan_pdf_lognormal_franken_lower_cut_100.0_pixie_1_z_min_0.001_z_excise_max_20_b_1.npz")
m_A_ary = file_lower['m_A_ary'] / eV
limit_lower_ary = file_lower['limit']

file_upper = np.load(outputs_dir + "firas_dp_scan_pdf_lognormal_franken_upper_cut_100.0_pixie_1_z_min_0.001_z_excise_max_20_b_1.npz")
m_A_ary = file_upper['m_A_ary'] / eV
limit_upper_ary = file_upper['limit']

# Homogeneous PIXIE curve, to stick at low masses

file_homo = np.load(outputs_dir + "firas_dp_scan_homo_pixie_1.npz")
m_A_ary = file_homo['m_A_ary'] / eV
limit_homo_ary = file_homo['limit']

lim_ary = np.max(np.transpose([limit_lower_ary, limit_upper_ary]), axis=1)

where = (limit_homo_ary < lim_ary) * (m_A_ary < 5e-15)
lim_ary[where] = limit_homo_ary[where]

plt.plot(m_A_ary, lim_ary, c=cols_default[0], ls='-.',lw=1)

# Existing and other projected bounds

## Jupiter bound
plt.fill_between(jup_m, jup_eps, np.ones_like(jup_eps), alpha=0.15, color=cols_default[6])

# Dark SRF projected bound
plt.plot([eps_SRF[0], eps_SRF[15]], [lim_SRF[0], lim_SRF[15]], ls='-.', c=cols_default[4], lw=1)

# Text labels

plt.text(1e-14, 4e-8, "PIXIE (projection)", color=cols_default[0], rotation=-20)
plt.text(1e-11, 2.5e-7, "Dark SRF (projection)", rotation=-34, color=cols_default[4])
plt.text(6e-15, 1e-1, "Jupiter", color=cols_default[6])
plt.text(2e-16, 2.5e-9, r"$10^{-2} < 1 + \delta < 10^2$", fontsize=14)

# Plot attributes

plt.xscale("log")
plt.yscale("log")

plt.xlabel(r"$m_{A^\prime}$\,[eV]") 
plt.ylabel(r"$\epsilon$")

plt.legend(fontsize=14, handlelength=1.3)

plt.title(r'\textbf{COBE/FIRAS $\gamma \rightarrow A^\prime$}')

plt.xlim(1e-16, 1e-9)
plt.ylim(1e-9, 0.3)

# Tick tweaks 

locmaj = matplotlib.ticker.LogLocator(base=10,numticks=12) 
ax.xaxis.set_major_locator(locmaj)

locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.2,0.4,0.6,0.8),numticks=12)
ax.xaxis.set_minor_locator(locmin)
ax.xaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())

locmaj = matplotlib.ticker.LogLocator(base=10,numticks=12) 
ax.yaxis.set_major_locator(locmaj)

locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.2,0.4,0.6,0.8),numticks=12)
ax.yaxis.set_minor_locator(locmin)
ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())

# Save

plt.tight_layout()
plt.savefig(plots_dir + "limit.pdf")<jupyter_output><empty_output><jupyter_text>Systematics plots  $1 + \delta$ systematics plots (Fig. 6. top left)<jupyter_code>def as_si(x, ndp, dolla=0):
    s = '{x:0.{ndp:d}e}'.format(x=x, ndp=ndp)
    m, e = s.split('e')
    if float(m) == 1.:
        s = r'10^{{{e:d}}}'.format(m=m, e=int(e))
    else:
        s = r'{m:s}\times 10^{{{e:d}}}'.format(m=m, e=int(e))
    if dolla:
        return '$' + s + '$'
    else:
        return s
colors_seq = sns.color_palette('Blues', n_colors=4) # sns.cubehelix_palette(4, start=1, rot=0)
colors_alt_seq = sns.color_palette('Reds', n_colors=4)
fig, ax=plt.subplots()

# 1 + \delta variations for log-normal with frankenspectrum (lower)
    
for i_d, one_plus_delta_bound in enumerate([10, 1e2, 1e3, 1e4]):

    file_lower = np.load(outputs_dir + "firas_dp_scan_pdf_lognormal_franken_lower_cut_" + str(one_plus_delta_bound) + "_pixie_0_z_min_0.001_z_excise_max_20_b_1.npz")
    m_A_ary = file_lower['m_A_ary'] / eV
    limit_lower_ary = file_lower['limit']

    plt.plot(m_A_ary, limit_lower_ary, label="$" + as_si(1 / one_plus_delta_bound, 1) + " < 1 + \delta < " + as_si(one_plus_delta_bound, 1) + "$", color=colors_seq[i_d])

# 1 + \delta variations for analytic PDF

for i_d, one_plus_delta_bound in enumerate([10, 1e2, 1e3, 1e4]):
        
    file_lower = np.load(outputs_dir + "firas_dp_scan_pdf_log_analytic_pdf_interp_franken_lower_cut_" + str(one_plus_delta_bound) + "_pixie_0_z_min_0.001_z_excise_max_20.npz")
    m_A_ary = file_lower['m_A_ary'] / eV
    limit_lower_ary = file_lower['limit']

    plt.plot(m_A_ary, limit_lower_ary, ls='--', color=colors_alt_seq[i_d])

# Jupyter bound

plt.fill_between(jup_m, jup_eps, np.ones_like(jup_eps), alpha=0.15, color=cols_default[6])

# Homogeneous curve

file_homo = np.load(outputs_dir + "firas_dp_scan_homo_pixie_0.npz")
m_A_ary = file_homo['m_A_ary'] / eV
limit_homo_ary = file_homo['limit']
limit_homo_ary[(limit_homo_ary > 1e-6) * (m_A_ary > 1e-13)] = savgol_filter(limit_homo_ary[(limit_homo_ary > 1e-6) * (m_A_ary > 1e-13)], 71, 3)

plt.plot(m_A_ary, limit_homo_ary, color='grey', ls='dotted', label="Homogeneous")

# Plot attributes

plt.plot([],[], label="Log-normal PDF", color=colors_seq[-1])
plt.plot([],[], ls='--', label="Analytic PDF", color=colors_alt_seq[-1])

plt.text(6e-15, 1e-1, "Jupiter", color=cols_default[6])

plt.xscale("log")
plt.yscale("log")

plt.xlabel(r"$m_{A^\prime}$\,[eV]") 
plt.ylabel(r"$\epsilon$")

lines = ax.get_lines()

## Create a legend for the first line.
first_legend = plt.legend(handles=lines[:-3], loc='upper right', fontsize=14, handlelength=1.3)

## Add the legend manually to the current Axes.
plt.gca().add_artist(first_legend)

## Create another legend for the second line.
plt.legend(handles=lines[-3:], loc='lower left', fontsize=14, handlelength=1.3)

plt.title(r'\textbf{$\gamma \rightarrow A^\prime$ PDF tails}')

plt.xlim(1e-16, 1e-9)
plt.ylim(1e-9, 0.3)

locmaj = matplotlib.ticker.LogLocator(base=10,numticks=12) 
ax.xaxis.set_major_locator(locmaj)

locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.2,0.4,0.6,0.8),numticks=12)
ax.xaxis.set_minor_locator(locmin)
ax.xaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())

locmaj = matplotlib.ticker.LogLocator(base=10,numticks=12) 
ax.yaxis.set_major_locator(locmaj)

locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.2,0.4,0.6,0.8),numticks=12)
ax.yaxis.set_minor_locator(locmin)
ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())

plt.tight_layout()
plt.savefig(plots_dir + "limit_dp_pdf_tails.pdf")<jupyter_output>/Users/smsharma/anaconda3/lib/python3.7/site-packages/matplotlib/legend.py:449: UserWarning: The handle <matplotlib.lines.Line2D object at 0x1c2f1ae5f8> has a label of '_line4' which cannot be automatically added to the legend.
  'legend.'.format(handle, label))
/Users/smsharma/anaconda3/lib/python3.7/site-packages/matplotlib/legend.py:449: UserWarning: The handle <matplotlib.lines.Line2D object at 0x1c2f1aed30> has a label of '_line5' which cannot be automatically added to the legend.
  'legend.'.format(handle, label))
/Users/smsharma/anaconda3/lib/python3.7/site-packages/matplotlib/legend.py:449: UserWarning: The handle <matplotlib.lines.Line2D object at 0x1c2f1aee80> has a label of '_line6' which cannot be automatically added to the legend.
  'legend.'.format(handle, label))
/Users/smsharma/anaconda3/lib/python3.7/site-packages/matplotlib/legend.py:449: UserWarning: The handle <matplotlib.lines.Line2D object at 0x1c2f1bb0f0> has a label of '_line7' which cannot be automatically adde[...]<jupyter_text>PDF systematics (Fig. 6, top right)<jupyter_code>fig, ax=plt.subplots()

# Homogeneous curve

file_homo = np.load(outputs_dir + "firas_dp_scan_homo_pixie_0.npz")
m_A_ary = file_homo['m_A_ary'] / eV
limit_homo_ary = file_homo['limit']
limit_homo_ary[(limit_homo_ary > 1e-6) * (m_A_ary > 1e-13)] = savgol_filter(limit_homo_ary[(limit_homo_ary > 1e-6) * (m_A_ary > 1e-13)], 71, 3)

plt.plot(m_A_ary, limit_homo_ary, color='grey', ls='dotted', label="Homogeneous")

plt.gca().set_prop_cycle(None)

# Lognormal with frankenspectrum curve

file_lower = np.load(outputs_dir + "firas_dp_scan_pdf_lognormal_franken_lower_cut_100.0_pixie_0_z_min_0.001_z_excise_max_20_b_1.npz")
m_A_ary = file_lower['m_A_ary'] / eV
limit_lower_ary = file_lower['limit']

plt.plot(m_A_ary, limit_lower_ary, label=r"Log-normal PDF")
    
# b = 1.5 curve

b = 1.5
file_lower = np.load(outputs_dir + "firas_dp_scan_pdf_lognormal_franken_lower_cut_100.0_pixie_0_z_min_0.001_z_excise_max_20_b_" + str(b) + ".npz")
m_A_ary = file_lower['m_A_ary'] / eV
limit_lower_ary = file_lower['limit']

plt.plot(m_A_ary, limit_lower_ary, label="$b=" + str(b) + "$")

# Analytic PDF curve

file = np.load(outputs_dir + "firas_dp_scan_pdf_log_analytic_pdf_interp_franken_upper_cut_100.0_pixie_0_z_min_0.001_z_excise_max_20.npz")
m_A_ary = file['m_A_ary'] / eV
limit_ary = file['limit']
plt.plot(m_A_ary, limit_ary, label=r"Analytic PDF")

# Voids curve

file_lower = np.load(outputs_dir + "firas_dp_scan_pdf_voids_franken_lower_cut_-1_pixie_0_z_min_0.001_z_excise_max_20.npz")
m_A_ary = file_lower['m_A_ary'] / eV
limit_lower_ary = file_lower['limit']
plt.plot(m_A_ary, limit_lower_ary, label="Voids PDF")

# Gaussian curve

file_lower = np.load(outputs_dir + "firas_dp_scan_pdf_gaussian_franken_lower_cut_-1_pixie_0_z_min_0.001_z_excise_max_20.npz")
m_A_ary = file_lower['m_A_ary'] / eV
limit_lower_ary = file_lower['limit']
plt.plot(m_A_ary, limit_lower_ary, label="Gaussian PDF", ls='dotted')

# Jupyter bound
plt.fill_between(jup_m, jup_eps, np.ones_like(jup_eps), alpha=0.15, color=cols_default[6])
plt.text(6e-15, 1e-1, "Jupiter", color=cols_default[6])

# Plot attributes

plt.xscale("log")
plt.yscale("log")

plt.xlabel(r"$m_{A^\prime}$\,[eV]") 
plt.ylabel(r"$\epsilon$")

plt.legend(fontsize=14, handlelength=1.3)

plt.title(r'\textbf{$\gamma \rightarrow A^\prime$ PDF systematics}')

plt.xlim(1e-16, 1e-9)
plt.ylim(1e-9, 0.3)

locmaj = matplotlib.ticker.LogLocator(base=10,numticks=12) 
ax.xaxis.set_major_locator(locmaj)

locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.2,0.4,0.6,0.8),numticks=12)
ax.xaxis.set_minor_locator(locmin)
ax.xaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())

locmaj = matplotlib.ticker.LogLocator(base=10,numticks=12) 
ax.yaxis.set_major_locator(locmaj)

locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.2,0.4,0.6,0.8),numticks=12)
ax.yaxis.set_minor_locator(locmin)
ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())

plt.text(2e-16, 2.5e-9, r"$10^{-2} < 1 + \delta < 10^2$", fontsize=14)

plt.tight_layout()
plt.savefig(plots_dir + "limit_dp_pdfs.pdf")<jupyter_output><empty_output><jupyter_text>Redshift systematics<jupyter_code>fig, ax=plt.subplots(figsize=(8 / 1.3, 6.5 / 1.3))

# Homogeneous curve

file_homo = np.load(outputs_dir + "firas_dp_scan_homo_pixie_0.npz")
m_A_ary = file_homo['m_A_ary'] / eV
limit_homo_ary = file_homo['limit']
limit_homo_ary[(limit_homo_ary > 1e-6) * (m_A_ary > 1e-13)] = savgol_filter(limit_homo_ary[(limit_homo_ary > 1e-6) * (m_A_ary > 1e-13)], 71, 3)

plt.plot(m_A_ary, limit_homo_ary, color='grey', ls='dotted', label="Homogeneous")

plt.gca().set_prop_cycle(None)

# Lognormal with frankenspectrum curve

file_lower = np.load(outputs_dir + "firas_dp_scan_pdf_lognormal_franken_lower_cut_100.0_pixie_0_z_min_0.001_z_excise_max_20_b_1.npz")
m_A_ary = file_lower['m_A_ary'] / eV
limit_lower_ary = file_lower['limit']

plt.plot(m_A_ary, limit_lower_ary, label=r"Log-normal PDF")

# Lognormal with frankenspectrum curve, z_min = 20

file_lower = np.load(outputs_dir + "firas_dp_scan_pdf_lognormal_franken_lower_cut_100.0_pixie_0_z_min_20_z_excise_max_20_b_1.npz")
m_A_ary = file_lower['m_A_ary'] / eV
limit_lower_ary = file_lower['limit']
limit_lower_ary[(limit_lower_ary > 1e-6) * (m_A_ary > 1e-13)] = savgol_filter(limit_lower_ary[(limit_lower_ary > 1e-6) * (m_A_ary > 1e-13)], 71, 3)
plt.plot(m_A_ary, limit_lower_ary, label=r'Log-normal $z_\mathrm{min} = 20$')

# Lognormal with frankenspectrum curve, z_min = 0.1

file_lower = np.load(outputs_dir + "firas_dp_scan_pdf_lognormal_franken_lower_cut_100.0_pixie_0_z_min_0.1_z_excise_max_20_b_1.npz")
m_A_ary = file_lower['m_A_ary'] / eV
limit_lower_ary = file_lower['limit']
plt.plot(m_A_ary, limit_lower_ary, label=r'Log-normal $z_\mathrm{min} = 0.1$')

# Lognormal with frankenspectrum curve, z_min = 

file_lower = np.load(outputs_dir + "firas_dp_scan_pdf_lognormal_franken_lower_cut_100.0_pixie_0_z_min_1_z_excise_max_20_b_1.npz")
m_A_ary = file_lower['m_A_ary'] / eV
limit_lower_ary = file_lower['limit']
plt.plot(m_A_ary, limit_lower_ary, label=r'Log-normal $z_\mathrm{min} = 1$')

# Lognormal with frankenspectrum curve, excise 6 < z < 30

file_lower = np.load(outputs_dir + "firas_dp_scan_pdf_lognormal_franken_lower_cut_100.0_pixie_0_z_min_0.001_z_excise_max_30_b_1.npz")
m_A_ary = file_lower['m_A_ary'] / eV
limit_lower_ary = file_lower['limit']
plt.plot(m_A_ary, limit_lower_ary, label=r"Log-normal excise $6 < z < 30$", ls='--')

# Analytic PDF, z_min = 20.

file = np.load(outputs_dir + "firas_dp_scan_pdf_log_analytic_pdf_interp_franken_upper_cut_100.0_pixie_0_z_min_20_z_excise_max_20.npz")
m_A_ary = file['m_A_ary'] / eV
limit_ary = file['limit']
limit_ary[(limit_ary > 1e-6) * (m_A_ary > 1e-13)] = savgol_filter(limit_ary[(limit_ary > 1e-6) * (m_A_ary > 1e-13)], 71, 3)
plt.plot(m_A_ary, limit_ary, label=r'Analytic $z_\mathrm{min} = 20$', ls='--', color=cols_default[9])

# Jupyter bound
plt.fill_between(jup_m, jup_eps, np.ones_like(jup_eps), alpha=0.15, color=cols_default[6])

plt.text(6e-15, 1e-1, "Jupiter", color=cols_default[6])

plt.xscale("log")
plt.yscale("log")

plt.xlabel(r"$m_{A^\prime}$\,[eV]") 
plt.ylabel(r"$\epsilon$")

plt.legend(fontsize=14, handlelength=1.3, loc='upper right')

plt.title(r'\textbf{$\gamma \rightarrow A^\prime$ redshift variations}')

plt.xlim(1e-16, 1e-9)
plt.ylim(1e-9, 0.3)

plt.text(2e-16, 2.5e-9, r"$10^{-2} < 1 + \delta < 10^2$", fontsize=14)

locmaj = matplotlib.ticker.LogLocator(base=10,numticks=12) 
ax.xaxis.set_major_locator(locmaj)

locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.2,0.4,0.6,0.8),numticks=12)
ax.xaxis.set_minor_locator(locmin)
ax.xaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())

locmaj = matplotlib.ticker.LogLocator(base=10,numticks=12) 
ax.yaxis.set_major_locator(locmaj)

locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.2,0.4,0.6,0.8),numticks=12)
ax.yaxis.set_minor_locator(locmin)
ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())

plt.tight_layout()
plt.savefig(plots_dir + "limit_dp_redshifts.pdf")<jupyter_output><empty_output><jupyter_text>$k_\mathrm{max}$ plot<jupyter_code>firas = FIRAS(PowerSpectrumGridInterpolator("franken_lower"))
m_ary = np.array([4e-15, 1e-13, 1e-12]) * eV
k_ary = np.logspace(-1, 4, 100)
lim_ary = np.zeros((len(m_ary), len(k_ary)))
for i_m, m in enumerate(m_ary):
    for i_k, k in enumerate(tqdm_notebook(k_ary)):
        lim_ary[i_m, i_k] = firas.get_lim([firas.chi2_FIRAS_scan(m, one_plus_delta_bound=1e2, z_int=[1e-3, 1e6], pdf="lognormal", k_max=k)])
# Homogeneous curve

file_homo = np.load(outputs_dir + "firas_dp_scan_homo_pixie_0.npz")
m_A_ary = file_homo['m_A_ary'] / eV
limit_homo_ary = file_homo['limit']
pylab.rcParams.update(params)

labels = [r"$4\times10^{-15}$\,eV","$10^{-13}$\,eV","$10^{-12}$\,eV"]
for i_m, m in enumerate(m_ary):
    plt.plot(k_ary, lim_ary[i_m], label=labels[i_m], color=cols_default[i_m])
    lim_homo = interp1d(m_A_ary, limit_homo_ary)(m / eV)
    if lim_homo < 1.:
        plt.axhline(lim_homo, color=cols_default[i_m], ls='--')
        
plt.plot([],[], c='k', ls='--', label='Homogeneous')

plt.xscale("log")
plt.yscale("log")

plt.xlabel("$k_\mathrm{max}$\,[$h$\,Mpc$^{-1}$]")
plt.ylabel(r"$\epsilon$")

plt.title(r"\textbf{$\gamma\to A'$ dependence on $k_\mathrm{max}$}")
plt.legend(fontsize=14, handlelength=1.3)
plt.xlim(min(k_ary), max(k_ary))
plt.ylim(1e-7, 0.03)

ax = plt.gca()

locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.2,0.4,0.6,0.8),numticks=12)
ax.xaxis.set_minor_locator(locmin)
ax.xaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())

locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.2,0.4,0.6,0.8),numticks=12)
ax.yaxis.set_minor_locator(locmin)
ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())


plt.tight_layout()
plt.savefig("../paper/draft-letter/plots/limit_dp_k_max.pdf")<jupyter_output>/Users/smsharma/anaconda3/lib/python3.7/site-packages/matplotlib/__init__.py:855: MatplotlibDeprecationWarning: 
examples.directory is deprecated; in the future, examples will be found relative to the 'datapath' directory.
  "found relative to the 'datapath' directory.".format(key))
/Users/smsharma/anaconda3/lib/python3.7/site-packages/matplotlib/__init__.py:846: MatplotlibDeprecationWarning: 
The text.latex.unicode rcparam was deprecated in Matplotlib 2.2 and will be removed in 3.1.
  "2.2", name=key, obj_type="rcparam", addendum=addendum)<jupyter_text>Plot for data points<jupyter_code>fig, ax=plt.subplots(figsize=(8 / 1.3, 6.5 / 1.3))

# Lognormal with frankenspectrum curve

file_lower = np.load(outputs_dir + "firas_dp_scan_widemass_1_pdf_lognormal_franken_lower_cut_100.0_pixie_0_z_min_0.001_z_excise_max_20_b_1.npz")
m_A_ary = file_lower['m_A_ary'] / eV
limit_lower_ary = file_lower['limit']

file_upper = np.load(outputs_dir + "firas_dp_scan_widemass_1_pdf_lognormal_franken_upper_cut_100.0_pixie_0_z_min_0.001_z_excise_max_20_b_1.npz")
m_A_ary = file_upper['m_A_ary'] / eV
limit_upper_ary = file_upper['limit']

lim_LN_franken = np.max(np.transpose([limit_lower_ary, limit_upper_ary]), axis=1)

plt.plot(m_A_ary, lim_LN_franken, label=r"Log-normal PDF")

# Jupyter bound
plt.fill_between(jup_m, jup_eps, np.ones_like(jup_eps), alpha=0.15, color=cols_default[6])

plt.text(6e-15, 1e-1, "Jupiter", color=cols_default[6])

plt.xscale("log")
plt.yscale("log")

plt.xlabel(r"$m_{A^\prime}$\,[eV]") 
plt.ylabel(r"$\epsilon$")

plt.legend(fontsize=14, handlelength=1.3, loc='upper right')

plt.title(r'\textbf{Fiducial limit}')

plt.xlim(1e-16, 1e-7)
plt.ylim(1e-9, 0.3)

plt.text(2e-16, 2.5e-9, r"$10^{-2} < 1 + \delta < 10^2$", fontsize=14)
df_fid = pd.DataFrame(np.transpose([m_A_ary[m_A_ary < 1.1e-9], lim_LN_franken[m_A_ary < 1.1e-9]]), columns=['m_Ap [eV]', '\epsilon'])
df_fid.to_csv("../data/constraints/fiducial_DP_FIRAS_one_plus_delta_1e2.csv", index = False)
file_homo = np.load(outputs_dir + "firas_dp_scan_homo_pixie_0.npz")
m_A_ary = file_homo['m_A_ary'] / eV
limit_homo_ary = file_homo['limit']

df_fid = pd.DataFrame(np.transpose([m_A_ary[m_A_ary < 1.1e-9], limit_homo_ary[m_A_ary < 1.1e-9]]), columns=['m_Ap [eV]', '\epsilon'])
df_fid.to_csv("../data/constraints/fiducial_DP_FIRAS_homo.csv", index = False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1571566265060241, 'fraction_numerical': 0.04308433734939759, 'mean_word_length': 6.637467795362532, 'pattern_counts': {'":': 0, '<': 52, '<?xml version=': 0, '>': 55, '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': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '654551', 'license': 'MIT', 'n_tokens_mistral': 9867, 'n_tokens_neox': 9218, 'n_words': 1230, 'path': 'notebooks/02_firas_dp_limits.ipynb', 'repo_name': 'andrea0292/Dark-Photon-Perturbations', 'size': '20750'}