Example of a spectral analysis with 3ML
Giacomo Vianello (Stanford University) giacomov@stanford.edu
IPython Notebook setup (only if you are using the IPython Notebook)
This line will activate the support for inline display of matplotlib images:
%matplotlib inline
Import 3ML and see the available plug-ins
from threeML import *
getAvailablePlugins()
Define some general features for the source of interest
Here we define a name and the coordinates (Equatorial J2000). In this case, we are analyzing a GRB:
triggerName = 'bn090217206'
ra = 204.9
dec = -8.4
Setup the use of LAT data
Here we instanciate a plugin to handle Fermi/LAT data:
datadir = '/home/giacomov/science/hawc/jointLikelihood/examples/090217206/'
eventFile = os.path.join(datadir,'gll_ft1_tr_bn090217206_v00_filt.fit')
ft2File = '/home/giacomov/FermiData/bn090217206/gll_ft2_tr_bn090217206_v00.fit'
expomap = os.path.join(datadir,'gll_ft1_tr_bn090217206_v00_filt_expomap.fit')
ltcube = os.path.join(datadir,'gll_ft1_tr_bn090217206_v00_filt_ltcube.fit')
xmlModel = os.path.join(datadir,"bn090217206_LAT_xmlmodel.xml")
irf = "P7REP_TRANSIENT_V15"
fl = FermiLATUnbinnedLike("LAT_p7rep_transient",eventFile,ft2File,expomap,ltcube,xmlModel,irf)
Setup of GBM data
Here we create 3 instances of the plugin which handle GBM data, for 3 different GBM detectors (NaI6, NaI9, BGO1). We also define the channels (or energies) we want to use, by using setActiveMeasurements():
obsSpectrum = os.path.join(datadir,"bn090217206_n6_srcspectra.pha{1}")
bakSpectrum = os.path.join(datadir,"bn090217206_n6_bkgspectra.bak{1}")
rspFile = os.path.join(datadir,"bn090217206_n6_weightedrsp.rsp{1}")
#Plugin instance
NaI6 = FermiGBMLike("NaI6",obsSpectrum,bakSpectrum,rspFile)
#Choose chanels to use
NaI6.setActiveMeasurements("10.0-30.0 40.0-950.0")
obsSpectrum = os.path.join(datadir,"bn090217206_n9_srcspectra.pha{1}")
bakSpectrum = os.path.join(datadir,"bn090217206_n9_bkgspectra.bak{1}")
rspFile = os.path.join(datadir,"bn090217206_n9_weightedrsp.rsp{1}")
#Plugin instance
NaI9 = FermiGBMLike("NaI9",obsSpectrum,bakSpectrum,rspFile)
#Choose chanels to use
NaI9.setActiveMeasurements("10.0-30.0 40.0-950.0")
obsSpectrum = os.path.join(datadir,"bn090217206_b1_srcspectra.pha{1}")
bakSpectrum = os.path.join(datadir,"bn090217206_b1_bkgspectra.bak{1}")
rspFile = os.path.join(datadir,"bn090217206_b1_weightedrsp.rsp{1}")
#Plugin instance
BGO1 = FermiGBMLike("BGO1",obsSpectrum,bakSpectrum,rspFile)
#Choose chanels to use
BGO1.setActiveMeasurements("2-126")
Create a dataList object to collect all our data
datalist = dataList(fl,NaI6,NaI9,BGO1)
Define our model
Here we define our spectral model, and the init values for the parameters. Each value can be a straight number, or a full specification with a list of [value, initial delta, hard minimum, soft minimum, soft maximum, hard maximum]. For further information, see the Xspec manual for the newpar command.
mod = ModelManager(triggerName,ra,dec,"grbm",datalist)
mod.setParameterValue("grbm.alpha","-1.01482")
mod.setParameterValue("grbm.beta","-2.21990")
mod.setParameterValue("grbm.tem","526.247")
mod.setParameterValue("grbm.norm",[0.0173822,0.001,1e-12,1e-12,1e5,1e5])
print mod
Create a jointLikelihood object
jl = JointLikelihood(mod)
Profile likelihood maximization
If you want to use the profile likelihood, just perform the fit:
(p_mle, p_1SigmaErr, logLikeMin) = jl.fit()
Full Bayesian analysis
- Define your model and its initial values (as before):
mod = ModelManager(triggerName,ra,dec,"grbm",datalist)
mod.setParameterValue("grbm.alpha","-1.01482")
mod.setParameterValue("grbm.beta","-2.21990")
mod.setParameterValue("grbm.tem","526.247")
mod.setParameterValue("grbm.norm",[0.0173822,0.001,1e-12,1e-12,1e5,1e5])
- Define the priors for all your parameters:
#Define priors for model parameters
mod.setParameterPrior("grbm.alpha",Priors.UniformPrior(-10,10))
mod.setParameterPrior("grbm.beta",Priors.UniformPrior(-10,10))
mod.setParameterPrior("grbm.tem",Priors.UniformPrior(10,1e4))
#Use a log-uniform (Jeffrey's) prior for the normalization
mod.setParameterPrior("grbm.norm",Priors.LogUniformPrior(1e-3,1e6))
#Define priors for the nuisance parameters
mod.setNuisanceParameterPrior("IsotropicTemplate-Normalization",Priors.UniformPrior(0.85,1.15))
mod.setNuisanceParameterPrior("bn090217206-Normalization",Priors.UniformPrior(0.85,1.15))
mod.setNuisanceParameterValue("IsotropicTemplate-Normalization",1.0)
print mod
- Create a jointLikelihood object:
jl = jointLikelihood.JointLikelihood(mod)
- Run the Markov Chain MC and explore the parameter space: (this will take some time, depending on the number of samples)
nwalkers = 100
samplesPerWalker = 3000
jl.explore(nwalkers,samplesPerWalker)
Getting some useful information from the Bayesian analysis
Now that we have our MCMC sample of parameters, we can extract a lot of useful information. We can start by taking the usual median values, and the credibility intervals (the Bayesian equivalent of confidence intervals). This is made with the getPercentiles method, which uses 1-sigma intervals by default (i.e., 16 and 84 percentiles):
percentiles = jl.getPercentiles()
If we want instead the 90% interval, we can request the appropriate percentiles:
percentiles = jl.getPercentiles(levels=[50,5,95])
We can also easily produce a triangle plot which show all the monodimensional and bidimensional marginal distribution, with the latter being the equivalent of countour plots in frequentist analysis:
import triangle
fig = triangle.corner(jl.samples,
labels=[r'$\alpha$',r'$\beta$','tem',r'A$_{Band}$','N','N1'],
quantiles=[0.16, 0.5, 0.84],
show_titles=True, label_args={"fontsize": 16})
By using the MCMC samples we can also plot the best fit model along with its "1-sigma contour":
bestFitParam = map(lambda x:x[0],percentiles.values())
nufnu = mod.plotWithChain(bestFitParam,jl.samples,kind="nufnu",ymin=1e-19,emin=1e-3,emax=1e12,nbins=1000)
