I need to display a scatter plot of the primary and secondary photon radiation of a black hole via the simulation software BlackHawk on windows under wsl. For my code I'm using matplot lib with jupyter as well as infrastructure provided by this github project https://github.com/vmmunoza/nuHawkHunter. I can graph the table into a line of best fit, but can't seem to graph it directly from a table to a scatter plot.
from Source.constants import *
# PBH masses
Mpbhs = [3.e12, 3.e15]
cols = ["r", "b", "g", "m", "y"]
# 1 for plotting only nu_e, 0 for the sum of the three flavours
only_nue = 1
for mm, Mpbh in enumerate(Mpbhs):
# Data folder
folder = folder_blackhawk + "{:.1e}".format(Mpbh) + "/"
data_primary = np.genfromtxt(folder+"instantaneous_primary_spectra.txt",skip_header = 2)
data_secondary = np.genfromtxt(folder+"instantaneous_secondary_spectra.txt",skip_header = 2)
if only_nue:
# Plot only photons
tot_sec = data_secondary[:,1]
tot_prim = data_primary[:,2]/3.
else:
# Plot sum of three neutrinos
tot_sec = 0.
# Sum over three neutrino species
for i in [3,4,5]:
tot_sec += data_secondary[:,i]
tot_prim = data_primary[:,6]
Esec = data_secondary[:,0]*1.e3 # in MeV
Eprim = data_primary[:,0]*1.e3
intsec = interp1d(Esec, tot_sec/1.e3) # in MeV
intprim = interp1d(Eprim, tot_prim/1.e3) # in MeV
#plt.loglog(Esec,intprim(Esec),linestyle=":",linewidth = 2, color=cols[mm])
#plt.loglog(Esec,intsec(Esec)-intprim(Esec), linestyle="--",linewidth = 2, color=cols[mm])
#plt.loglog(Esec,intsec(Esec), linestyle="-",linewidth = 2, color=cols[mm])
plt.scatter(Eprim, intprim(Eprim), marker='o', color=cols[mm])
# Figure legend
customlegend = []
for n, Mpbh in enumerate(Mpbhs):
customlegend.append( Line2D([0], [0], color=cols[n], lw=4, label = r"$M_{\rm PBH}=$"+scinot(Mpbh)+" g"))
#customlegend.append( Line2D([0], [0], color="black", linestyle=":", label="Primary"))
#customlegend.append( Line2D([0], [0], color="black", linestyle="--", label="Secondary"))
#customlegend.append( Line2D([0], [0], color="black", linestyle="-", label="Total"))
plt.ylim(1.e12,1.e21)
plt.xlim(1.e0,1.e3)
plt.xlabel(r'$E_\nu{\rm \,\, [MeV]}$')
plt.ylabel(r'$\frac{d^2N_\nu}{dE_\nu dt} \,\, [{\rm MeV}^{-1}{\rm s}^{-1}]$')
plt.tick_params(axis='both', which='both', top=True, right=True, direction="in")
plt.grid()
plt.legend(handles=customlegend)
plt.savefig("figures/photon_spectra_blackhawk.png", bbox_inches='tight', dpi=300)
plt.show()
We attempted to graph it, but it only did so one dimensionally across the x-axis.
