In [1]:

import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline
%matplotlib notebook

from threeML import *
from astromodels.xspec import *

# The filter library takes a while to load so you must import it explicitly..

from threeML.plugins.photometry.filter_library import threeML_filter_library

Configuration read from /Users/jburgess/.threeML/threeML_config.yml
Loading xspec models...done
Loading optical filters

In [2]:

get_available_plugins()

Available plugins:

FermiGBMTTELike for Fermi GBM TTE (all detectors)
FermiLATLLELike for Fermi LAT LLE
OGIPLike for All OGIP-compliant instruments
FermipyLike for Fermi LAT (with fermipy)
XYLike for n.a.
EventListLike for Generic EventList data
PhotometryLike for Generic photometric data
SwiftXRTLike for Swift XRT
VERITASLike for VERITAS
SpectrumLike for General binned spectral data

Setup¶

We use speclite to handle optical filters: http://speclite.readthedocs.io/en/latest/ .

Therefore, you can easily build your own custom filters, use the built in speclite filters, or use the 3ML filter library that we have built thanks to Spanish Virtual Observatory (http://svo.cab.inta-csic.es/main/index.php).

If you use these filters, please be sure to cite the proper sources!

Simple example of building a filter¶

Let's say we have our own 1-m telescope with a Johnson filter and we happen to record the data. We also have simultaneous data at other wavelengths and we want to compare. Let's setup the optical plugin (we'll ignore the other data for now).

In [3]:

import speclite.filters as spec_filters

my_backyard_telescope_filter = spec_filters.load_filter('bessell-r')

# NOTE:
my_backyard_telescope_filter.name

Out[3]:

'bessell-R'

NOTE: the filter name is 'bessell-R'. The plugin will look for the name after the '-' i.e 'R'

Now let's build a 3ML plugin via PhotometryLike.

Our data are entered as keywords with the name of the filter as the keyword and the data in an magnitude,error tuple, i.e. R=(mag,mag_err):

In [4]:

my_backyard_telescope = PhotometryLike('backyard_astronomy',
                         filters=my_backyard_telescope_filter,
                         R=(20,.1) )

my_backyard_telescope.display_filters()

Using chi2 statistic with the provided errors.

3ML filter library¶

Explore the filter library. If you cannot find what you need, it is simple to add your own

In [5]:

threeML_filter_library.SLOAN

Out[5]:

SDSS:
- u
- g
- r
- i
- z

In [6]:

spec_filters.plot_filters(threeML_filter_library.SLOAN.SDSS)

In [7]:

spec_filters.plot_filters(threeML_filter_library.Herschel.SPIRE)

In [8]:

spec_filters.plot_filters(threeML_filter_library.Keck.NIRC2)

Build your own filters¶

Following the example from speclite, we can build our own filters and add them:

In [9]:

fangs_g = spec_filters.FilterResponse(
    wavelength = [3800, 4500, 5200] * u.Angstrom,
    response = [0, 0.5, 0], meta=dict(group_name='fangs', band_name='g'))
fangs_r = spec_filters.FilterResponse(
    wavelength = [4800, 5500, 6200] * u.Angstrom,
    response = [0, 0.5, 0], meta=dict(group_name='fangs', band_name='r'))

fangs = spec_filters.load_filters('fangs-g', 'fangs-r')

fangslike = PhotometryLike('fangs',filters=fangs,g=(20,.1),r=(18,.1))


fangslike.display_filters()

Using chi2 statistic with the provided errors.

GROND Example¶

Now we will look at GROND. We get the filter from the 3ML filter library.

(Just play with tab completion to see what is available!)

In [10]:

grond = PhotometryLike('GROND',
                       filters=threeML_filter_library.ESO.GROND,
                       #g=(21.5.93,.23),
                       #r=(22.,0.12),
                       i=(21.8,.01),
                       z=(21.2,.01),
                       J=(19.6,.01),
                       H=(18.6,.01),
                       K=(18.,.01))

Using chi2 statistic with the provided errors.

In [11]:

grond.display_filters()

Model specification¶

Here we use XSPEC's dust extinction models for the milky way and the host

In [12]:

spec =  Powerlaw() * XS_zdust() * XS_zdust()

data_list = DataList(grond)

model = Model(PointSource('grb',0,0,spectral_shape=spec))

spec.piv_1 = 1E-2
spec.index_1.fix=False
spec.redshift_2 = 0.347
spec.redshift_2.fix = True

spec.e_bmv_2 = 5./2.93
spec.e_bmv_2.fix = True
spec.rv_2 = 2.93
spec.rv_2.fix = True


spec.method_2 = 3
spec.method_2.fix=True



spec.e_bmv_3 = .002/3.08
spec.e_bmv_3.fix = True
spec.rv_3= 3.08
spec.rv_3.fix=True
spec.redshift_3 = 0
spec.redshift_3.fix=True
spec.method_3 = 1
spec.method_3.fix=True

jl = JointLikelihood(model,data_list)

We compute $m_{\rm AB}$ from astromodels photon fluxes. This is done by convolving the differential flux over the filter response:

$ F[R,f_\lambda] \equiv \int_0^\infty \frac{dg}{d\lambda}(\lambda)R(\lambda) \omega(\lambda) d\lambda$

where we have converted the astromodels functions to wavelength properly.

In [13]:

jl.set_minimizer('ROOT')
_ = jl.fit()

Best fit values:

	Value	Unit
grb.spectrum.main.composite.K_1	(4.61 +/- 0.12) x 10	1 / (cm2 keV s)
grb.spectrum.main.composite.index_1	-1.144 +/- 0.011

Correlation matrix:

1.00	0.99
0.99	1.00

Values of -log(likelihood) at the minimum:

	-log(likelihood)
GROND	954.636274
total	954.636274

Values of statistical measures:

	statistical measures
AIC	1919.272548
BIC	1912.491424

Examine the fit¶

In [14]:

_=display_photometry_model_magnitudes(jl)

In [15]:

new_grond = grond.get_simulated_dataset()

Using chi2 statistic with the provided errors.

In [16]:

_ = plot_point_source_spectra(jl.results,flux_unit='erg/(cm2 s keV)',
                              xscale='linear',
                              energy_unit='nm',ene_min=1E3, ene_max=1E5 )

In [ ]: