One of the key features of 3ML is the abil ity to fit multi-messenger data properly. A simple example of this is the joint fitting of two instruments whose data obey different likelihoods. Here, we have GBM data which obey a Poisson-Gaussian profile likelihoog ( PGSTAT in XSPEC lingo) and Swift BAT which data which are the result of a "fit" via a coded mask and hence obey a Gaussian ( $\chi^2$ ) likelihood.
In [1]:
%matplotlib notebook
from threeML import *
import os
We have data from the same time interval from Swift BAT and a GBM NAI and BGO detector. We have preprocessed GBM data to so that it is OGIP compliant. (Remember that we can handle the raw data with the TimeSeriesBuilder). Thus, we will use the OGIPLike plugin to read in each dataset, make energy selections and examine the raw count spectra.
In [3]:
gbm_dir = "../../examples/gbm"
bat_dir = "../../examples/bat"
bat = OGIPLike('BAT',
observation=os.path.join(bat_dir,'gbm_bat_joint_BAT.pha'),
response=os.path.join(bat_dir,'gbm_bat_joint_BAT.rsp'))
bat.set_active_measurements('15-150')
bat.view_count_spectrum()
In [4]:
nai6 = OGIPLike('n6',
os.path.join(gbm_dir,'gbm_bat_joint_NAI_06.pha'),
os.path.join(gbm_dir,'gbm_bat_joint_NAI_06.bak'),
os.path.join(gbm_dir,'gbm_bat_joint_NAI_06.rsp'),
spectrum_number=1)
nai6.set_active_measurements('8-900')
nai6.view_count_spectrum()
bgo0 = OGIPLike('b0',
os.path.join(gbm_dir,'gbm_bat_joint_BGO_00.pha'),
os.path.join(gbm_dir,'gbm_bat_joint_BGO_00.bak'),
os.path.join(gbm_dir,'gbm_bat_joint_BGO_00.rsp'),
spectrum_number=1)
bgo0.set_active_measurements('250-30000')
bgo0.view_count_spectrum()
In [7]:
band = Band()
model_no_eac = Model(PointSource('joint_fit_no_eac',0,0,spectral_shape=band))
In [8]:
data_list = DataList(bat,nai6,bgo0)
jl_no_eac = JointLikelihood(model_no_eac, data_list)
jl_no_eac.fit();
The fit has resulted in a very typical Band function fit. Let's look in count space at how good of a fit we have obtained.
In [10]:
threeML_config['ogip']['model plot cmap']='Set1'
In [12]:
display_spectrum_model_counts(jl_no_eac, step=False, min_rate=[.01,10.,10.],data_colors=['grey','k','k']);
It seems that the effective areas between GBM and BAT do not agree! We can look at the goodness of fit for the various data sets.
In [24]:
gof_object = GoodnessOfFit(jl_no_eac)
with parallel_computation(profile='default',start_cluster=False):
gof, res_frame, lh_frame = gof_object.by_mc(n_iterations=8000)
In [27]:
import pandas as pd
pd.Series(gof)
Out[27]:
Both the GBM NaI detector and Swift BAT exhibit poor GOF.
Now let's add an effective area correction between the detectors to see if this fixes the problem. The effective area is a nuissance parameter that attempts to model systematic problems in a instruments calibration. It simply scales the counts of an instrument by a multiplicative factor. It cannot handle more complicated energy dependent
In [35]:
# turn on the effective area correction and set it's bounds
nai6.use_effective_area_correction(.2,1.8)
bgo0.use_effective_area_correction(.2,1.8)
model_eac = Model(PointSource('joint_fit_eac',0,0,spectral_shape=band))
jl_eac = JointLikelihood(model_eac, data_list)
jl_eac.fit();
Now we have a much better fit to all data sets
In [36]:
display_spectrum_model_counts(jl_eac, step=False,min_rate=[.01,10.,10.],data_colors=['grey','k','k']);
In [32]:
gof_object = GoodnessOfFit(jl_eac)
with parallel_computation(profile='default',start_cluster=False):
gof, res_frame, lh_frame = gof_object.by_mc(n_iterations=8000,continue_on_failure=True)
In [33]:
import pandas as pd
pd.Series(gof)
Out[33]:
In [44]:
plot_point_source_spectra(jl_eac.results,
jl_no_eac.results,
fit_cmap = 'Set1',
contour_cmap = 'Set1',
flux_unit='erg2/(keV s cm2)',
equal_tailed=True);
We can easily see that the models are different