Spectral fitting example (GRB)
To run this, you need the following files, which can be downloaded using the first few cells of this notebook: - orientation file (20280301_3_month.ori)
- binned data (grb_bkg_binned_data.hdf5, grb_binned_data.hdf5, & bkg_binned_data_1s_local.hdf5)
- detector response (SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.h5.zip)
The binned data are simulations of GRB090206620 and albedo photon background produced using the COSI SMEX mass model. The detector response needs to be unzipped before running the notebook.
This notebook fits the spectrum of a GRB simulated using MEGAlib and combined with background.
3ML is a high-level interface that allows multiple datasets from different instruments to be used coherently to fit the parameters of source model. A source model typically consists of a list of sources with parametrized spectral shapes, sky locations and, for extended sources, shape. Polarization is also possible. A “coherent” analysis, in this context, means that the source model parameters are fitted using all available datasets simultanously, rather than performing individual fits and finding a well-suited common model a posteriori.
In order for a dataset to be included in 3ML, each instrument needs to provide a “plugin”. Each plugin is responsible for reading the data, convolving the source model (provided by 3ML) with the instrument response, and returning a likelihood. In our case, we’ll compute a binned Poisson likelihood:
\[\log \mathcal{L}(\mathbf{x}) = \sum_i \log \frac{\lambda_i(\mathbf{x})^{d_i} \exp (-\lambda_i)}{d_i!}\]
where \(d_i\) are the counts on each bin and \(\lambda_i\) are the expected counts given a source model with parameters \(\mathbf{x}\).
In this example, we will fit a single point source with a known location. We’ll assume the background is known and fixed up to a scaling factor. Finally, we will fit a Band function:
\[\begin{split}f(x) = K \begin{cases} \left(\frac{x}{E_{piv}}\right)^{\alpha} \exp \left(-\frac{(2+\alpha)
* x}{x_{p}}\right) & x \leq (\alpha-\beta) \frac{x_{p}}{(\alpha+2)} \\ \left(\frac{x}{E_{piv}}\right)^{\beta}
* \exp (\beta-\alpha)\left[\frac{(\alpha-\beta) x_{p}}{E_{piv}(2+\alpha)}\right]^{\alpha-\beta}
* &x>(\alpha-\beta) \frac{x_{p}}{(\alpha+2)} \end{cases}\end{split}\]
where \(K\) (normalization), \(\alpha\) & \(\beta\) (spectral indeces), and \(x_p\) (peak energy) are the free parameters, while \(E_{piv}\) is the pivot energy which is fixed (and arbitrary).
Considering these assumptions:
\[\lambda_i(\mathbf{x}) = B*b_i + s_i(\mathbf{x})\]
where \(B*b_i\) are the estimated counts due to background in each bin of the Compton data space with \(B\) the amplitude and \(b_i\) the shape of the background, and \(s_i\) are the corresponding expected counts from the source, the goal is then to find the values of \(\mathbf{x} = [K, \alpha, \beta, x_p]\) and \(B\) that maximize \(\mathcal{L}\). These are the best estimations of the parameters.
The final module needs to also fit the time-dependent background, handle multiple point-like and extended sources, as well as all the spectral models supported by 3ML. Eventually, it will also fit the polarization angle. However, this simple example already contains all the necessary pieces to do a fit.
[1]:
from cosipy import COSILike, BinnedData
from cosipy.spacecraftfile import SpacecraftFile
from cosipy.response.FullDetectorResponse import FullDetectorResponse
from cosipy.util import fetch_wasabi_file
from scoords import SpacecraftFrame
from astropy.time import Time
import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.stats import poisson_conf_interval
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from threeML import Band, PointSource, Model, JointLikelihood, DataList
from cosipy import Band_Eflux
from astromodels import Parameter
from pathlib import Path
import os
12:04:24 WARNING The naima package is not available. Models that depend on it will not be functions.py:48 available
WARNING The GSL library or the pygsl wrapper cannot be loaded. Models that depend on it functions.py:69 will not be available.
12:04:25 WARNING We have set the min_value of K to 1e-99 because there was a postive transform parameter.py:704
WARNING We have set the min_value of K to 1e-99 because there was a postive transform parameter.py:704
WARNING We have set the min_value of K to 1e-99 because there was a postive transform parameter.py:704
WARNING We have set the min_value of K to 1e-99 because there was a postive transform parameter.py:704
WARNING We have set the min_value of F to 1e-99 because there was a postive transform parameter.py:704
WARNING We have set the min_value of K to 1e-99 because there was a postive transform parameter.py:704
12:04:25 INFO Starting 3ML! __init__.py:35
WARNING WARNINGs here are NOT errors __init__.py:36
WARNING but are inform you about optional packages that can be installed __init__.py:37
WARNING to disable these messages, turn off start_warning in your config file __init__.py:40
WARNING ROOT minimizer not available minimization.py:1345
WARNING Multinest minimizer not available minimization.py:1357
WARNING PyGMO is not available minimization.py:1369
WARNING The cthreeML package is not installed. You will not be able to use plugins which __init__.py:94 require the C/C++ interface (currently HAWC)
WARNING Could not import plugin HAWCLike.py. Do you have the relative instrument __init__.py:144 software installed and configured?
WARNING Could not import plugin FermiLATLike.py. Do you have the relative instrument __init__.py:144 software installed and configured?
WARNING No fermitools installed lat_transient_builder.py:44
WARNING Env. variable OMP_NUM_THREADS is not set. Please set it to 1 for optimal __init__.py:387 performances in 3ML
WARNING Env. variable MKL_NUM_THREADS is not set. Please set it to 1 for optimal __init__.py:387 performances in 3ML
WARNING Env. variable NUMEXPR_NUM_THREADS is not set. Please set it to 1 for optimal __init__.py:387 performances in 3ML
Download and read in binned data
Define the path to the directory containing the data, detector response, orientation file, and yaml files if they have already been downloaded, or the directory to download the files into
[2]:
data_path = Path("/path/to/files")
Download the orientation file (684.38 MB)
[6]:
fetch_wasabi_file('COSI-SMEX/DC2/Data/Orientation/20280301_3_month.ori', output=str(data_path / '20280301_3_month.ori'))
Download the binned GRB+background data (75.73 KB)
[8]:
fetch_wasabi_file('COSI-SMEX/cosipy_tutorials/grb_spectral_fit_local_frame/grb_bkg_binned_data.hdf5', output=str(data_path / 'grb_bkg_binned_data.hdf5'))
Download the binned GRB data (76.90 KB)
[19]:
fetch_wasabi_file('COSI-SMEX/cosipy_tutorials/grb_spectral_fit_local_frame/grb_binned_data.hdf5', output=str(data_path / 'grb_binned_data.hdf5'))
Download the binned background data (255.97 MB)
[20]:
fetch_wasabi_file('COSI-SMEX/cosipy_tutorials/grb_spectral_fit_local_frame/bkg_binned_data_1s_local.hdf5', output=str(data_path / 'bkg_binned_data_1s_local.hdf5'))
Download the response file (839.62 MB). This needs to be unzipped before running the rest of the notebook
[17]:
fetch_wasabi_file('COSI-SMEX/DC2/Responses/SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.h5.zip', output=str(data_path / 'SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.h5.zip'))
Read in the spacecraft orientation file & select the beginning and end times of the GRB
[4]:
ori = SpacecraftFile.parse_from_file(data_path / "20280301_3_month.ori")
tmin = Time(1842597410.0,format = 'unix')
tmax = Time(1842597450.0,format = 'unix')
sc_orientation = ori.source_interval(tmin, tmax)
Create BinnedData objects for the GRB only, GRB+background, and background only. The GRB only simulation is not used for the spectral fit, but can be used to compare the fitted spectrum to the source simulation
[5]:
grb = BinnedData(data_path / "grb.yaml")
grb_bkg = BinnedData(data_path / "grb.yaml")
bkg = BinnedData(data_path / "background.yaml")
Load binned .hdf5 files
[6]:
grb.load_binned_data_from_hdf5(binned_data=data_path / "grb_binned_data.hdf5")
grb_bkg.load_binned_data_from_hdf5(binned_data=data_path / "grb_bkg_binned_data.hdf5")
bkg.load_binned_data_from_hdf5(binned_data=data_path / "bkg_binned_data_1s_local.hdf5")
Define the path to the detector response
[7]:
dr = str(data_path / "SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.h5") # path to detector response
Perform spectral fit
Define time window of binned background simulation to use for background model
[8]:
bkg_tmin = 1842597310.0
bkg_tmax = 1842597550.0
bkg_min = np.where(bkg.binned_data.axes['Time'].edges.value == bkg_tmin)[0][0]
bkg_max = np.where(bkg.binned_data.axes['Time'].edges.value == bkg_tmax)[0][0]
Set background parameter, which is used to fit the amplitude of the background, and instantiate the COSI 3ML plugin
[9]:
bkg_par = Parameter("background_cosi", # background parameter
0.1, # initial value of parameter
min_value=0, # minimum value of parameter
max_value=5, # maximum value of parameter
delta=1e-3, # initial step used by fitting engine
desc="Background parameter for cosi")
cosi = COSILike("cosi", # COSI 3ML plugin
dr = dr, # detector response
data = grb_bkg.binned_data.project('Em', 'Phi', 'PsiChi'), # data (source+background)
bkg = bkg.binned_data.slice[{'Time':slice(bkg_min,bkg_max)}].project('Em', 'Phi', 'PsiChi'), # background model
sc_orientation = sc_orientation, # spacecraft orientation
nuisance_param = bkg_par) # background parameter
Define a point source at the known location with a Band function spectrum and add it to the model
[10]:
l = 93.
b = -53.
alpha = -1 # Setting parameters to something reasonable helps the fitting to converge\n",
beta = -3
xp = 450. * u.keV
piv = 500. * u.keV
K = 1 / u.cm / u.cm / u.s / u.keV
spectrum = Band()
spectrum.beta.min_value = -15.0
spectrum.alpha.value = alpha
spectrum.beta.value = beta
spectrum.xp.value = xp.value
spectrum.K.value = K.value
spectrum.piv.value = piv.value
spectrum.xp.unit = xp.unit
spectrum.K.unit = K.unit
spectrum.piv.unit = piv.unit
source = PointSource("source", # Name of source (arbitrary, but needs to be unique)
l = l, # Longitude (deg)
b = b, # Latitude (deg)
spectral_shape = spectrum) # Spectral model
# Optional: free the position parameters
#source.position.l.free = True
#source.position.b.free = True
model = Model(source) # Model with single source. If we had multiple sources, we would do Model(source1, source2, ...)
# Optional: if you want to call get_log_like manually, then you also need to set the model manually
# 3ML does this internally during the fit though
cosi.set_model(model)
... Calculating point source responses ...
Now converting to the Spacecraft frame...
Conversion completed!
--> done (source name : source)
--> all done
Gather all plugins and combine with the model in a JointLikelihood object, then perform maximum likelihood fit
[11]:
plugins = DataList(cosi) # If we had multiple instruments, we would do e.g. DataList(cosi, lat, hawc, ...)
like = JointLikelihood(model, plugins, verbose = False)
like.fit()
12:04:56 INFO set the minimizer to minuit joint_likelihood.py:1042
Adding 1e-12 to each bin of the expectation to avoid log-likelihood = -inf.
12:05:22 WARNING get_number_of_data_points not implemented, values for statistical plugin_prototype.py:128 measurements such as AIC or BIC are unreliable
Best fit values:
| result | unit | |
|---|---|---|
| parameter | ||
| source.spectrum.main.Band.K | (3.10 -0.20 +0.21) x 10^-2 | 1 / (cm2 keV s) |
| source.spectrum.main.Band.alpha | (-2.8 +/- 0.5) x 10^-1 | |
| source.spectrum.main.Band.xp | (4.75 +/- 0.05) x 10^2 | keV |
| source.spectrum.main.Band.beta | -6.8 +/- 1.2 | |
| background_cosi | (1.65 +/- 0.13) x 10^-1 |
Correlation matrix:
| 1.00 | 0.97 | -0.37 | 0.20 | -0.00 |
| 0.97 | 1.00 | -0.16 | 0.18 | -0.00 |
| -0.37 | -0.16 | 1.00 | -0.18 | -0.02 |
| 0.20 | 0.18 | -0.18 | 1.00 | 0.00 |
| -0.00 | -0.00 | -0.02 | 0.00 | 1.00 |
Values of -log(likelihood) at the minimum:
| -log(likelihood) | |
|---|---|
| cosi | 42920.049336 |
| total | 42920.049336 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 85838.098672 |
| BIC | 85840.098672 |
[11]:
( value negative_error positive_error \
source.spectrum.main.Band.K 0.030995 -0.001939 0.002116
source.spectrum.main.Band.alpha -0.276632 -0.050195 0.049825
source.spectrum.main.Band.xp 474.650732 -4.896060 4.873162
source.spectrum.main.Band.beta -6.756966 -1.215124 1.201150
background_cosi 0.164969 -0.012573 0.012434
error unit
source.spectrum.main.Band.K 0.002027 1 / (cm2 keV s)
source.spectrum.main.Band.alpha 0.050010
source.spectrum.main.Band.xp 4.884611 keV
source.spectrum.main.Band.beta 1.208137
background_cosi 0.012504 ,
-log(likelihood)
cosi 42920.049336
total 42920.049336)
Error propagation and plotting
Define Band function spectrum injected into MEGAlib
[12]:
alpha_inj = -0.360
beta_inj = -11.921
E0_inj = 288.016 * u.keV
xp_inj = E0_inj * (alpha_inj + 2)
piv_inj = 1. * u.keV
K_inj = 0.283 / u.cm / u.cm / u.s / u.keV
spectrum_inj = Band()
spectrum_inj.beta.min_value = -15.0
spectrum_inj.alpha.value = alpha_inj
spectrum_inj.beta.value = beta_inj
spectrum_inj.xp.value = xp_inj.value
spectrum_inj.K.value = K_inj.value
spectrum_inj.piv.value = piv_inj.value
spectrum_inj.xp.unit = xp_inj.unit
spectrum_inj.K.unit = K_inj.unit
spectrum_inj.piv.unit = piv_inj.unit
The summary of the results above tell you the optimal values of the parameters, as well as the errors. Propogate the errors to the “evaluate_at” method of the spectrum
[13]:
results = like.results
print(results.display())
parameters = {par.name:results.get_variates(par.path)
for par in results.optimized_model["source"].parameters.values()
if par.free}
results_err = results.propagate(results.optimized_model["source"].spectrum.main.shape.evaluate_at, **parameters)
print(results.optimized_model["source"])
Best fit values:
| result | unit | |
|---|---|---|
| parameter | ||
| source.spectrum.main.Band.K | (3.10 -0.20 +0.21) x 10^-2 | 1 / (cm2 keV s) |
| source.spectrum.main.Band.alpha | (-2.8 +/- 0.5) x 10^-1 | |
| source.spectrum.main.Band.xp | (4.75 +/- 0.05) x 10^2 | keV |
| source.spectrum.main.Band.beta | -6.8 +/- 1.2 | |
| background_cosi | (1.65 +/- 0.13) x 10^-1 |
Correlation matrix:
| 1.00 | 0.97 | -0.37 | 0.20 | -0.00 |
| 0.97 | 1.00 | -0.16 | 0.18 | -0.00 |
| -0.37 | -0.16 | 1.00 | -0.18 | -0.02 |
| 0.20 | 0.18 | -0.18 | 1.00 | 0.00 |
| -0.00 | -0.00 | -0.02 | 0.00 | 1.00 |
Values of -log(likelihood) at the minimum:
| -log(likelihood) | |
|---|---|
| cosi | 42920.049336 |
| total | 42920.049336 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 85838.098672 |
| BIC | 85840.098672 |
None
* source (point source):
* position:
* l:
* value: 93.0
* desc: Galactic longitude
* min_value: 0.0
* max_value: 360.0
* unit: deg
* is_normalization: false
* b:
* value: -53.0
* desc: Galactic latitude
* min_value: -90.0
* max_value: 90.0
* unit: deg
* is_normalization: false
* equinox: J2000
* spectrum:
* main:
* Band:
* K:
* value: 0.030994516909178687
* desc: Differential flux at the pivot energy
* min_value: 1.0e-99
* max_value: null
* unit: keV-1 s-1 cm-2
* is_normalization: true
* alpha:
* value: -0.27663221293105034
* desc: low-energy photon index
* min_value: -1.5
* max_value: 3.0
* unit: ''
* is_normalization: false
* xp:
* value: 474.6507320770641
* desc: peak in the x * x * N (nuFnu if x is a energy)
* min_value: 10.0
* max_value: null
* unit: keV
* is_normalization: false
* beta:
* value: -6.756965748051311
* desc: high-energy photon index
* min_value: -15.0
* max_value: -1.6
* unit: ''
* is_normalization: false
* piv:
* value: 500.0
* desc: pivot energy
* min_value: null
* max_value: null
* unit: keV
* is_normalization: false
* polarization: {}
Evaluate the flux and errors at a range of energies for the fitted and injected spectra, and the simulated source flux
[14]:
energy = np.geomspace(100*u.keV,10*u.MeV).to_value(u.keV)
flux_lo = np.zeros_like(energy)
flux_median = np.zeros_like(energy)
flux_hi = np.zeros_like(energy)
flux_inj = np.zeros_like(energy)
for i, e in enumerate(energy):
flux = results_err(e)
flux_median[i] = flux.median
flux_lo[i], flux_hi[i] = flux.equal_tail_interval(cl=0.68)
flux_inj[i] = spectrum_inj.evaluate_at(e)
binned_energy_edges = grb.binned_data.axes['Em'].edges.value
binned_energy = np.array([])
bin_sizes = np.array([])
for i in range(len(binned_energy_edges)-1):
binned_energy = np.append(binned_energy, (binned_energy_edges[i+1] + binned_energy_edges[i]) / 2)
bin_sizes = np.append(bin_sizes, binned_energy_edges[i+1] - binned_energy_edges[i])
expectation = cosi._expected_counts['source']
Plot the fitted and injected spectra
[15]:
fig,ax = plt.subplots()
ax.plot(energy, energy*energy*flux_median, label = "Best fit")
ax.fill_between(energy, energy*energy*flux_lo, energy*energy*flux_hi, alpha = .5, label = "Best fit (errors)")
ax.plot(energy, energy*energy*flux_inj, color = 'black', ls = ":", label = "Injected")
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel("Energy (keV)")
ax.set_ylabel(r"$E^2 \frac{dN}{dE}$ (keV cm$^{-2}$ s$^{-1}$)")
ax.legend()
[15]:
<matplotlib.legend.Legend at 0x29adcdf10>
Plot the fitted spectrum convolved with the response, as well as the simulated source counts
[16]:
fit_poisson_error = np.zeros((2,len(expectation.project('Em').todense().contents)))
fit_gaussian_error = np.zeros(len(expectation.project('Em').todense().contents))
inj_poisson_error = np.zeros((2,len(grb.binned_data.project('Em').todense().contents)))
inj_gaussian_error = np.zeros(len(grb.binned_data.project('Em').todense().contents))
for i, counts in enumerate(expectation.project('Em').todense().contents):
if counts > 5:
fit_gaussian_error[i] = np.sqrt(counts)
else:
poisson_error = poisson_conf_interval(counts, interval="frequentist-confidence", sigma=1)
fit_poisson_error[0][i] = poisson_error[0]
fit_poisson_error[1][i] = poisson_error[1]
for i, counts in enumerate(grb.binned_data.project('Em').todense().contents):
if counts > 5:
inj_gaussian_error[i] = np.sqrt(counts)
else:
poisson_error = poisson_conf_interval(counts, interval="frequentist-confidence", sigma=1)
inj_poisson_error[0][i] = poisson_error[0]
inj_poisson_error[1][i] = poisson_error[1]
fig,ax = plt.subplots()
ax.stairs(expectation.project('Em').todense().contents, binned_energy_edges, color='purple', label = "Best fit convolved with response")
ax.errorbar(binned_energy, expectation.project('Em').todense().contents, yerr=fit_poisson_error, color='purple', linewidth=0, elinewidth=1)
ax.errorbar(binned_energy, expectation.project('Em').todense().contents, yerr=fit_gaussian_error, color='purple', linewidth=0, elinewidth=1)
ax.stairs(grb.binned_data.project('Em').todense().contents, binned_energy_edges, color = 'black', ls = ":", label = "Source counts")
ax.errorbar(binned_energy, grb.binned_data.project('Em').todense().contents, yerr=inj_poisson_error, color='black', linewidth=0, elinewidth=1)
ax.errorbar(binned_energy, grb.binned_data.project('Em').todense().contents, yerr=inj_gaussian_error, color='black', linewidth=0, elinewidth=1)
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel("Energy (keV)")
ax.set_ylabel("Counts")
ax.legend()
[16]:
<matplotlib.legend.Legend at 0x29bdaa1f0>
Plot the fitted spectrum convolved with the response plus the fitted background, as well as the simulated source+background counts
[17]:
fit_bkg_poisson_error = np.zeros((2,len(expectation.project('Em').todense().contents+(bkg_par.value * bkg.binned_data.slice[{'Time':slice(bkg_min,bkg_max)}].project('Em').todense().contents))))
fit_bkg_gaussian_error = np.zeros(len(expectation.project('Em').todense().contents+(bkg_par.value * bkg.binned_data.slice[{'Time':slice(bkg_min,bkg_max)}].project('Em').todense().contents)))
inj_bkg_poisson_error = np.zeros((2,len(grb_bkg.binned_data.project('Em').todense().contents)))
inj_bkg_gaussian_error = np.zeros(len(grb_bkg.binned_data.project('Em').todense().contents))
for i, counts in enumerate(expectation.project('Em').todense().contents+(bkg_par.value * bkg.binned_data.slice[{'Time':slice(bkg_min,bkg_max)}].project('Em').todense().contents)):
if counts > 5:
fit_bkg_gaussian_error[i] = np.sqrt(counts)
else:
poisson_error = poisson_conf_interval(counts, interval="frequentist-confidence", sigma=1)
fit_bkg_poisson_error[0][i] = poisson_error[0]
fit_bkg_poisson_error[1][i] = poisson_error[1]
for i, counts in enumerate(grb_bkg.binned_data.project('Em').todense().contents):
if counts > 5:
inj_bkg_gaussian_error[i] = np.sqrt(counts)
else:
poisson_error = poisson_conf_interval(counts, interval="frequentist-confidence", sigma=1)
inj_bkg_poisson_error[0][i] = poisson_error[0]
inj_bkg_poisson_error[1][i] = poisson_error[1]
fig,ax = plt.subplots()
ax.stairs(expectation.project('Em').todense().contents+(bkg_par.value * bkg.binned_data.slice[{'Time':slice(bkg_min,bkg_max)}].project('Em').todense().contents), binned_energy_edges, color='purple', label = "Best fit convolved with response plus background")
ax.errorbar(binned_energy, expectation.project('Em').todense().contents+(bkg_par.value * bkg.binned_data.slice[{'Time':slice(bkg_min,bkg_max)}].project('Em').todense().contents), yerr=fit_bkg_poisson_error, color='purple', linewidth=0, elinewidth=1)
ax.errorbar(binned_energy, expectation.project('Em').todense().contents+(bkg_par.value * bkg.binned_data.slice[{'Time':slice(bkg_min,bkg_max)}].project('Em').todense().contents), yerr=fit_bkg_gaussian_error, color='purple', linewidth=0, elinewidth=1)
ax.stairs(grb_bkg.binned_data.project('Em').todense().contents, binned_energy_edges, color = 'black', ls = ":", label = "Total counts")
ax.errorbar(binned_energy, grb_bkg.binned_data.project('Em').todense().contents, yerr=inj_bkg_poisson_error, color='black', linewidth=0, elinewidth=1)
ax.errorbar(binned_energy, grb_bkg.binned_data.project('Em').todense().contents, yerr=inj_bkg_gaussian_error, color='black', linewidth=0, elinewidth=1)
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel("Energy (keV)")
ax.set_ylabel("Counts")
ax.legend()
[17]:
<matplotlib.legend.Legend at 0x2af24a160>
