The SpectrumLike plugin is designed to handle binned photon/particle spectra. It comes in three basic classes:
The functionality of all three plugins is the same.
The most basic spectrum plugin is SpectrumLike which handles spectra with and without backgrounds. There are six basic features of a spectrum:
Let's start by examining an observation where the total counts are Poisson distributed and the measured background ground has been observed by viewing an off-source region and hence is also Poisson.
In [1]:
from threeML import *
%matplotlib notebook
import matplotlib.pyplot as plt
import numpy as np
We will construct a simulated spectrum over the energy range 10-1000 keV. The spectrum will have logrithmic energy boundaries.
We will simulate a blackbody source spectrum on top of powerlaw background.
In [2]:
energies = np.logspace(1,3,51)
low_edge = energies[:-1]
high_edge = energies[1:]
# get a blackbody source function
source_function = Blackbody(K=9E-2,kT=20)
# power law background function
background_function = Powerlaw(K=1,index=-1.5, piv=100.)
spectrum_generator = SpectrumLike.from_function('fake',
source_function=source_function,
background_function=background_function,
energy_min=low_edge,
energy_max=high_edge)
In [3]:
spectrum_generator.display()
These properties are accessible from the object. For example:
In [4]:
print(spectrum_generator.exposure)
print(spectrum_generator.significance)
print(spectrum_generator.observed_counts)
To view the count spectrum, we call the view_count_spectrum method:
In [5]:
fig = spectrum_generator.view_count_spectrum()
It is even possible see which channels are above a given significance threshold. Red regions are below the supplied significance regions.
In [6]:
fig = spectrum_generator.view_count_spectrum(significance_level=5)
Note: In 3ML, the Significance module is used to compute significnaces. When total counts ($N_{\rm on}$) are Poisson distributed and the background or off-source counts ($N_{\rm off}$) are also Poisson distributed, the significance in $\sigma$ is calculated via the likelihood ratio derived in Li & Ma (1980):
$$ \sigma = \sqrt{-2 \log \lambda} = \sqrt{2} \left( N_{\rm on} \log \left[ \frac{1+\alpha}{\alpha} \frac{N_{\rm on}}{N_{\rm on}+N_{\rm off}} \right] + N_{\rm off} \log \left[ (1 + \alpha)\frac{N_{\rm off}}{N_{\rm on}+N_{\rm off}} \right] \right)$$In the case that the background is Gaussian distributed, an equivalent likelihood ratio is used (see Vianello in prep).
Many times, there are channels that we are not valid for analysis due to poor instrument characteristics, overflow, or systematics. We then would like to mask or exclude these channels before fitting the spectrum. We provide several ways to do this and it is useful to consult the docstring. However, we review the process here.
NOTE to Xspec users: while XSpec uses integers and floats to distinguish between energies and channels specifications, 3ML does not, as it would be error-prone when writing scripts. Read the following documentation to know how to achieve the same functionality.
They are specified as 'emin-emax'. Energies are in keV.
In [7]:
spectrum_generator.set_active_measurements('10-12.5','56.0-100.0')
fig = spectrum_generator.view_count_spectrum()
In [8]:
spectrum_generator.set_active_measurements('c10-c12','c20-c50')
fig = spectrum_generator.view_count_spectrum()
In [9]:
spectrum_generator.set_active_measurements('0.2-c10','c20-1000')
fig = spectrum_generator.view_count_spectrum()
In [10]:
spectrum_generator.set_active_measurements('all')
fig = spectrum_generator.view_count_spectrum()
In [11]:
spectrum_generator.set_active_measurements(exclude=["c2-20", "50-c40"])
fig = spectrum_generator.view_count_spectrum()
In [12]:
spectrum_generator.set_active_measurements("0.2-c10",exclude=["c30-c50"])
fig = spectrum_generator.view_count_spectrum()
We can rebin the spectra based off a minimum total or background rate requried. This is useful when using profile likelihoods, however, we do not change the underlying likelihood by binning up the data. For more information, consult the statistics section.
To rebin a spectrum based off the total counts, we specify the minimum counts per bin we would like, say 100:
In [13]:
spectrum_generator.set_active_measurements("all")
spectrum_generator.rebin_on_source(100)
fig = spectrum_generator.view_count_spectrum()
We can remove the rebinning this way:
In [14]:
spectrum_generator.remove_rebinning()
Instead, when using a profile likelihood which requires at least one background count per bin to be valid, we would call:
In [15]:
spectrum_generator.rebin_on_background(10)
fig = spectrum_generator.view_count_spectrum()
In [16]:
spectrum_generator.remove_rebinning()
In [17]:
bb = Blackbody()
pts = PointSource('mysource',0,0,spectral_shape=bb)
model = Model(pts)
# MLE fitting
jl = JointLikelihood(model,DataList(spectrum_generator))
result = jl.fit()
In [18]:
count_fig1 = spectrum_generator.display_model(min_rate=10)
In [19]:
_ = plot_spectra(jl.results, flux_unit='erg/(cm2 s keV)')
_ = plt.ylim(1E-20)
Perhaps we want to fit a different model and compare the results. We change the spectral model and will overplot the fit's expected counts with the fit to the blackbody.
In [20]:
import warnings
pl = Powerlaw()
pts = PointSource('mysource',0,0,spectral_shape=pl)
model = Model(pts)
# MLE fitting
jl = JointLikelihood(model,DataList(spectrum_generator))
with warnings.catch_warnings():
warnings.simplefilter('ignore')
result = jl.fit()
In [21]:
spectrum_generator.display_model(min_rate=10,
show_data=False,
show_residuals=True,
data_color='g',
model_color='b',
model_label='powerlaw',
model_subplot=count_fig1.axes)
Out[21]:
Examining the fit in count space lets us easily see that the fit with the powerlaw model is very poor. We can of course deterimine the fit quality numerically, but this is saved for another section.
Instruments that exhibit energy dispersion must have their spectra fit through a process called forward folding. Let $R(\varepsilon,E)$ be our response converting between true (monte carlo) energy ($E$) and detector channel/energy ($\varepsilon$), $f(E, \vec{\phi}_{\rm s})$ be our photon model which is a function of $E$ and source model parameters $\vec{\phi}_{\rm s}$. Then, the source counts ($S_{c} (\vec{\phi}_{\rm s})$) registered in the detector between channel (c) with energy boundaries $E_{{\rm min}, c}$ and $E_{{\rm max}, c}$ (in the absence of background) are given by the convolution of the photon model with the response:
$$S_{c} (\vec{\phi}_{\rm s}) = \int_{0}^\infty {\rm d} E \, f(E, \vec{\phi}_{\rm s}) \int_{E_{{\rm min}, c}}^{E_{{\rm max}, c}} {\rm d} \varepsilon \, R(\varepsilon, E) $$Therefore, to fit the data in count space, we assume a photon model, fold it through the response, and calculate the predicted counts. This process is iterated on the source model parameters via likelihood minimization or posterior sampling until an optimal set of parameters is found.
To handle dispersed spectra, 3ML provides the DispersionSpectrumLike plugin.
In [22]:
from threeML.plugins.DispersionSpectrumLike import DispersionSpectrumLike
from threeML.utils.OGIP.response import OGIPResponse
from threeML.io.package_data import get_path_of_data_file
# we will use a demo response
response = OGIPResponse(get_path_of_data_file('datasets/ogip_powerlaw.rsp'))
source_function = Broken_powerlaw(K=1E-2,
alpha=0,
beta=-2,
xb=2000,
piv=200)
background_function = Powerlaw(K=10,
index=-1.5,
piv=100.)
dispersion_spectrum_generator = DispersionSpectrumLike.from_function('test',
source_function=source_function,
response=response,
background_function=background_function)
We can view the response and the count spectrum created.
In [23]:
_ = dispersion_spectrum_generator.display_rsp()
In [24]:
fig = dispersion_spectrum_generator.view_count_spectrum()
All the functionality of SpectrumLike is inherited in DispersionSpectrumLike. Therefore, fitting, and examination of the data is the same.
Finally, many x-ray mission provide data in the form of fits files known and pulse-height analysis (PHA) data. The keywords for the information in the data are known as the Office of Guest Investigators Program (OGIP) standard. While these data are always a form of binned spectra, 3ML provide a convience plugin for reading OGIP standard PHA Type I (single spectrum) and Type II (multiple spectrum) files.
The OGIPLike plugin inherits from DispersionSpectrumLike and thus needs either a full response or a redistribution matrix (RMF) and ancillary response (ARF) file. The plugin will prove the keywords in the data files to automatically figure out the correct likelihood for the observation.
In [25]:
ogip_data = OGIPLike('ogip',
observation=get_path_of_data_file('datasets/ogip_powerlaw.pha'),
background = get_path_of_data_file('datasets/ogip_powerlaw.bak'),
response=get_path_of_data_file('datasets/ogip_powerlaw.rsp'))
In [26]:
ogip_data.display()
In [27]:
fig = ogip_data.view_count_spectrum()
When using PHA Type II data, a spectrum number must be supplied to indicate which spectrum from the file to use. Users can also follow the XSPEC convention of specifying the spectrum number in the filename (e.g. 'my_pha.fits{1}')