Point Source Fluxes and Multiple Sources

Once one has computed a spectral to a point source, getting the flux of that source is possible. In 3ML, we can obtain flux in a variety of units in a live analysis or from saved fits. There is no need to know exactly what you want to obtain at the time you do the fit.

Also, let’s explore how to deal with fitting multiple point sources and linking of parameters.

Let’s explore the possibilites.

[1]:
import warnings

warnings.simplefilter("ignore")
[2]:
%%capture
import matplotlib.pyplot as plt
import numpy as np

np.seterr(all="ignore")
import astropy.units as u
from threeML import *
from threeML.utils.OGIP.response import OGIPResponse
from threeML.io.package_data import get_path_of_data_file
[3]:
from jupyterthemes import jtplot

%matplotlib inline
jtplot.style(context="talk", fscale=1, ticks=True, grid=False)
silence_warnings()
set_threeML_style()

Generating some synthetic data

alt text

Let’s say we have two galactic x-ray sources, some accreting compact binaries perhaps? We observe them at two different times. These sources (imaginary) sources emit a blackbody which is theorized to always be at the same temperature, but perhaps at different flux levels.

Lets simulate one of these sources:

[4]:
np.random.seed(1234)

# we will use a demo response
response_1 = OGIPResponse(get_path_of_data_file("datasets/ogip_powerlaw.rsp"))


source_function_1 = Blackbody(K=5e-8, kT=500.0)
background_function_1 = Powerlaw(K=1, index=-1.5, piv=1.0e3)


spectrum_generator_1 = DispersionSpectrumLike.from_function(
    "s1",
    source_function=source_function_1,
    background_function=background_function_1,
    response=response_1,
)

fig = spectrum_generator_1.view_count_spectrum()
23:27:29 INFO      Auto-probed noise models:                                                    SpectrumLike.py:490
         INFO      - observation: poisson                                                       SpectrumLike.py:491
         INFO      - background: None                                                           SpectrumLike.py:492
23:27:32 INFO      Auto-probed noise models:                                                    SpectrumLike.py:490
         INFO      - observation: poisson                                                       SpectrumLike.py:491
         INFO      - background: None                                                           SpectrumLike.py:492
         INFO      Auto-probed noise models:                                                    SpectrumLike.py:490
         INFO      - observation: poisson                                                       SpectrumLike.py:491
         INFO      - background: poisson                                                        SpectrumLike.py:492
         INFO      Auto-probed noise models:                                                    SpectrumLike.py:490
         INFO      - observation: poisson                                                       SpectrumLike.py:491
         INFO      - background: poisson                                                        SpectrumLike.py:492
../_images/notebooks_flux_examples_5_12.png

Now let’s simulate the other source, but this one has an extra feature! There is a power law component in addition to the blackbody.

[5]:

response_2 = OGIPResponse(get_path_of_data_file("datasets/ogip_powerlaw.rsp")) source_function_2 = Blackbody(K=1e-7, kT=500.0) + Powerlaw_flux( F=2e2, index=-1.5, a=10, b=500 ) background_function_2 = Powerlaw(K=1, index=-1.5, piv=1.0e3) spectrum_generator_2 = DispersionSpectrumLike.from_function( "s2", source_function=source_function_2, background_function=background_function_2, response=response_2, ) fig = spectrum_generator_2.view_count_spectrum()
23:27:33 INFO      Auto-probed noise models:                                                    SpectrumLike.py:490
         INFO      - observation: poisson                                                       SpectrumLike.py:491
         INFO      - background: None                                                           SpectrumLike.py:492
         INFO      Auto-probed noise models:                                                    SpectrumLike.py:490
         INFO      - observation: poisson                                                       SpectrumLike.py:491
         INFO      - background: None                                                           SpectrumLike.py:492
         INFO      Auto-probed noise models:                                                    SpectrumLike.py:490
         INFO      - observation: poisson                                                       SpectrumLike.py:491
         INFO      - background: poisson                                                        SpectrumLike.py:492
         INFO      Auto-probed noise models:                                                    SpectrumLike.py:490
         INFO      - observation: poisson                                                       SpectrumLike.py:491
         INFO      - background: poisson                                                        SpectrumLike.py:492
../_images/notebooks_flux_examples_7_12.png

Make the model

Now let’s make the model we will use to fit the data. First, let’s make the spectral function for source_1 and set priors on the parameters.

[6]:
spectrum_1 = Blackbody()

spectrum_1.K.prior = Log_normal(mu=np.log(1e-7), sigma=1)
spectrum_1.kT.prior = Log_normal(mu=np.log(300), sigma=2)

ps1 = PointSource("src1", ra=1, dec=20, spectral_shape=spectrum_1)

We will do the same for the other source but also include the power law component

[7]:
spectrum_2 = Blackbody() + Powerlaw_flux(
    a=10, b=500
)  # a,b are the bounds for the flux for this model

spectrum_2.K_1.prior = Log_normal(mu=np.log(1e-6), sigma=1)
spectrum_2.kT_1.prior = Log_normal(mu=np.log(300), sigma=2)

spectrum_2.F_2.prior = Log_normal(mu=np.log(1e2), sigma=1)
spectrum_2.F_2.bounds = (None, None)

spectrum_2.index_2.prior = Gaussian(mu=-1.0, sigma=1)
spectrum_2.index_2.bounds = (None, None)

ps2 = PointSource("src2", ra=2, dec=-10, spectral_shape=spectrum_2)

Now we can combine these two sources into our model.

[8]:
model = Model(ps1, ps2)

Linking parameters

We hypothesized that both sources should have the a same blackbody temperature. We can impose this by linking the temperatures.

[9]:
model.link(
    model.src1.spectrum.main.Blackbody.kT, model.src2.spectrum.main.composite.kT_1
)

we could also link the parameters with an arbitrary function rather than directly. Check out the astromodels documentation for more details.

[10]:
model
[10]:
Model summary:

N
Point sources 2
Extended sources 0
Particle sources 0


Free parameters (5):

value min_value max_value unit
src1.spectrum.main.Blackbody.K 0.0001 0.0 None s-1 cm-2 keV-3
src2.spectrum.main.composite.K_1 0.0001 0.0 None s-1 cm-2 keV-3
src2.spectrum.main.composite.kT_1 30.0 0.0 None keV
src2.spectrum.main.composite.F_2 1.0 0.0 None s-1 cm-2
src2.spectrum.main.composite.index_2 -2.0 None None


Fixed parameters (8):
(abridged. Use complete=True to see all fixed parameters)


Properties (0):

(none)


Linked parameters (1):

src1.spectrum.main.Blackbody.kT
linked to src2.spectrum.main.composite.kT_1
function Line
current value 30.0
unit keV


Independent variables:

(none)

Linked functions (0):

(none)

Assigning sources to plugins

Now, if we simply passed out model to the BayesianAnalysis or JointLikelihood objects, it would sum the point source spectra together and apply both sources to all data.

This is not what we want. Many plugins have the ability to be assigned directly to a source. Let’s do that here:

[11]:
spectrum_generator_1.assign_to_source("src1")

spectrum_generator_2.assign_to_source("src2")

Now we simply make our our data list

[12]:
data = DataList(spectrum_generator_1, spectrum_generator_2)

Fitting the data

Now we fit the data as we normally would. We use Bayesian analysis here.

[13]:
ba = BayesianAnalysis(model, data)
ba.set_sampler("ultranest")
ba.sampler.setup(frac_remain=0.5)
_ = ba.sample()
23:27:34 INFO      sampler set to ultranest                                                bayesian_analysis.py:202
[ultranest] Sampling 400 live points from prior ...
[ultranest] Explored until L=-1e+03
[ultranest] Likelihood function evaluations: 66343
[ultranest]   logZ = -1339 +- 0.1967
[ultranest] Effective samples strategy satisfied (ESS = 955.4, need >400)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.46+-0.09 nat, need <0.50 nat)
[ultranest] Evidency uncertainty strategy is satisfied (dlogz=0.45, need <0.5)
[ultranest]   logZ error budget: single: 0.23 bs:0.20 tail:0.41 total:0.45 required:<0.50
[ultranest] done iterating.
23:28:17 INFO      fit restored to maximum of posterior                                         sampler_base.py:178
         INFO      fit restored to maximum of posterior                                         sampler_base.py:178
Maximum a posteriori probability (MAP) point:

result unit
parameter
src1.spectrum.main.Blackbody.K (4.64 -0.24 +0.33) x 10^-8 1 / (s cm2 keV3)
src2.spectrum.main.composite.K_1 (9.7 -0.7 +0.6) x 10^-8 1 / (s cm2 keV3)
src2.spectrum.main.composite.kT_1 (5.11 -0.11 +0.08) x 10^2 keV
src2.spectrum.main.composite.F_2 (2.06 +/- 0.08) x 10^2 1 / (s cm2)
src2.spectrum.main.composite.index_2 -1.512 -0.009 +0.012
Values of -log(posterior) at the minimum:

-log(posterior)
s1 -597.69557
s2 -692.11621
total -1289.81178
Values of statistical measures:

statistical measures
AIC 2589.863559
BIC 2607.349446
DIC 2614.605206
PDIC 4.903353
log(Z) -581.594127

Let’s examine the fits.

[14]:
fig = display_spectrum_model_counts(ba)
ax = fig.get_axes()[0]
ax.set_ylim(1e-6)
[14]:
(1e-06, 3395.2542994245655)
../_images/notebooks_flux_examples_25_1.png

Lets grab the result. Remember, we can save the results to disk, so all of the following operations can be run at a later time without having to redo all the above steps!

[15]:
result = ba.results
fig = result.corner_plot()
../_images/notebooks_flux_examples_27_0.png

Computing fluxes

Now we will compute fluxes. We can compute them an many different units, over any energy range also specified in any units.

The flux is computed by integrating the function over energy. By default, a fast trapezoid method is used. If you need more accuracy, you can change the method in the configuration.

[16]:
threeML_config.point_source.integrate_flux_method = "quad"

result.get_flux(ene_min=1 * u.keV, ene_max=1 * u.MeV, flux_unit="erg/cm2/s")
[16]:
flux low bound hi bound
src1: total 5.73311249625044e-06 erg / (s cm2) 5.394264964329162e-06 erg / (s cm2) 6.087393247337274e-06 erg / (s cm2)
src2: total 4.840910661251699e-05 erg / (s cm2) 4.770640184786882e-05 erg / (s cm2) 4.918210539119791e-05 erg / (s cm2)

We see that a pandas dataframe is returned with all the information. We could change the confidence region for the uncertainties if we desire. However, we could also sum the source fluxes! 3ML will take care of propagating the uncertainties (for any of these operations).

[17]:
threeML_config.point_source.integrate_flux_method = "trapz"

result.get_flux(
    ene_min=1 * u.keV,
    ene_max=1 * u.MeV,
    confidence_level=0.95,
    sum_sources=True,
    flux_unit="erg/cm2/s",
)
[17]:
flux low bound hi bound
total 5.4210003773202634e-05 erg / (s cm2) 5.264137588840751e-05 erg / (s cm2) 5.5800121533929216e-05 erg / (s cm2)

We can get the fluxes of individual components:

[18]:
result.get_flux(
    ene_min=10 * u.keV, ene_max=0.5 * u.MeV, use_components=True, flux_unit="1/(cm2 s)"
)
[18]:
flux low bound hi bound
src1: total 2.1310616431530063 1 / (s cm2) 2.005784949623805 1 / (s cm2) 2.2594472388463136 1 / (s cm2)
src2: Blackbody 4.3785337616985345 1 / (s cm2) 4.172855762768196 1 / (s cm2) 4.601569953725661 1 / (s cm2)
src2: Powerlaw_flux 219.70061771728942 1 / (s cm2) 211.41284850731898 1 / (s cm2) 228.40154210240914 1 / (s cm2)

As well as choose which source and component to compute

[19]:
result.get_flux(
    ene_min=10 * u.keV,
    ene_max=0.5 * u.MeV,
    sources=["src2"],
    use_components=True,
    components_to_use=["Blackbody"],
    flux_unit="1/(cm2 s)",
)
[19]:
flux low bound hi bound
src2: Blackbody 4.362547384558303 1 / (s cm2) 4.172998293955476 1 / (s cm2) 4.599261622012086 1 / (s cm2)

Finally, the returned flux object is a pandas table and can be manipulated as such:

[20]:
flux = result.get_flux(ene_min=1 * u.keV, ene_max=1 * u.MeV, flux_unit="erg/cm2/s")
[21]:
flux["flux"]
[21]:
src1: total    5.7375210409431574e-06 erg / (s cm2)
src2: total    4.8447074378816425e-05 erg / (s cm2)
Name: flux, dtype: object
[22]:
flux["flux"]["src1: total"]
[22]:
$5.737521 \times 10^{-6} \; \mathrm{\frac{erg}{s\,cm^{2}}}$