Polarization example - maximum likelihood methodο
This notebook fits the polarization fraction and angle of a Data Challenge 3 GRB (GRB 080802386) simulated using MEGAlib and combined with albedo photon background. Itβs assumed that the start time, duration, localization, and spectrum of the GRB are already known. The GRB was simulated with 80% polarization at an angle of 90 degrees in the IAU convention, and was 20 degrees off-axis.
[1]:
%%capture
from cosipy import BinnedData
from cosipy.spacecraftfile import SpacecraftHistory
from cosipy.statistics import PoissonLikelihood
from cosipy.background_estimation import FreeNormBinnedBackground
from cosipy.interfaces import ThreeMLPluginInterface
from cosipy.response.FullDetectorResponse import FullDetectorResponse
from cosipy.response import BinnedThreeMLModelFolding, BinnedInstrumentResponse, BinnedThreeMLPointSourceResponse
from cosipy.data_io import EmCDSBinnedData
from cosipy.threeml.custom_functions import Band_Eflux
from cosipy.polarization import PolarizationAxis
from astropy.time import Time
from astropy.coordinates import SkyCoord
from astropy import units as u
from cosipy.util import fetch_wasabi_file
from pathlib import Path
import sys
from threeML import LinearPolarization, SpectralComponent, PointSource, Model, JointLikelihood, DataList
from astromodels import Parameter
Download and read in dataο
This will download the files needed to run this notebook. If you have already downloaded these files, you can skip this.
Download the unbinned data (660.58 KB)
[2]:
fetch_wasabi_file('COSI-SMEX/cosipy_tutorials/polarization_fit/grb_background.fits.gz', checksum = '21b1d75891edc6aaf1ff3fe46e91cb49')
A file named grb_background.fits.gz already exists with the specified checksum (21b1d75891edc6aaf1ff3fe46e91cb49). Skipping.
Download the polarization response (217.47 MB)
[3]:
fetch_wasabi_file('COSI-SMEX/develop/Data/Responses/ResponseContinuum.o3.pol.e200_10000.b4.p12.relx.s10396905069491.m420.filtered.binnedpolarization.11D.h5', checksum = '46b006a6b397fd777dc561d3b028357f')
A file named ResponseContinuum.o3.pol.e200_10000.b4.p12.relx.s10396905069491.m420.filtered.binnedpolarization.11D.h5 already exists with the specified checksum (46b006a6b397fd777dc561d3b028357f). Skipping.
Download the orientation file (1.10 GB)
[4]:
fetch_wasabi_file('COSI-SMEX/develop/Data/Orientation/DC3_final_530km_3_month_with_slew_1sbins_GalacticEarth_SAA.fits', checksum = '1b851c042acf4c909798e2401e9d2e38')
A file named DC3_final_530km_3_month_with_slew_1sbins_GalacticEarth_SAA.fits already exists with the specified checksum (1b851c042acf4c909798e2401e9d2e38). Skipping.
Read in and bin the data, which is a GRB placed within albedo photon background. A time cut is done for the duration of the GRB to produce the GRB+background data to fit. The time intervals before and after the GRB are used to produce a background model.
[5]:
data_path = Path("") # Update to your path
grb_background = BinnedData(data_path/'grb.yaml')
grb_background.select_data_time(unbinned_data=data_path/'grb_background.fits.gz', output_name=data_path/'grb_background_source_interval')
grb_background.get_binned_data(unbinned_data=data_path/'grb_background_source_interval.fits.gz', output_name=data_path/'grb_background_binned_galactic', psichi_binning='galactic')
grb_background.load_binned_data_from_hdf5(data_path/'grb_background_binned_galactic.hdf5')
background_before = BinnedData(data_path/'background_before.yaml')
background_before.select_data_time(unbinned_data=data_path/'grb_background.fits.gz', output_name=data_path/'background_before')
background_before.get_binned_data(unbinned_data='background_before.fits.gz', output_name='background_before_binned_galactic', psichi_binning='galactic')
background_before.load_binned_data_from_hdf5(data_path/'background_before_binned_galactic.hdf5')
background_after = BinnedData(data_path/'background_after.yaml') # e.g. background_after.yaml
background_after.select_data_time(unbinned_data=data_path/'grb_background.fits.gz', output_name=data_path/'background_after')
background_after.get_binned_data(unbinned_data=data_path/'background_after.fits.gz', output_name=data_path/'background_after_binned_galactic', psichi_binning='galactic')
background_after.load_binned_data_from_hdf5(data_path/'background_after_binned_galactic.hdf5')
fill() discarded one or more values due to out-of-bounds coordinate in a dimension without under/overflow tracking
fill() discarded one or more values due to out-of-bounds coordinate in a dimension without under/overflow tracking
fill() discarded one or more values due to out-of-bounds coordinate in a dimension without under/overflow tracking
Read in the detector response and orientation file. The orientation is cut down to the time interval of the source.
[6]:
response_file = data_path / 'ResponseContinuum.o3.pol.e200_10000.b4.p12.relx.s10396905069491.m420.filtered.binnedpolarization.11D.h5'
dr = FullDetectorResponse.open(response_file, pa_convention='RelativeX')
sc_orientation = SpacecraftHistory.open(data_path/'DC3_final_530km_3_month_with_slew_1sbins_GalacticEarth_SAA.fits', tstart=Time(1835493492.2, format = 'unix'), tstop=Time(1835493492.8, format = 'unix'))
Define the GRB position and spectrum.
[7]:
source_direction = SkyCoord(l=23.53, b=-53.44, frame='galactic', unit=u.deg)
a = 100. * u.keV
b = 10000. * u.keV
alpha = -0.7368949
beta = -2.095031
ebreak = 622.389 * u.keV
K = 300. / u.cm / u.cm / u.s
spectrum = Band_Eflux(a = a.value,
b = b.value,
alpha = alpha,
beta = beta,
E0 = ebreak.value,
K = K.value)
spectrum.a.unit = a.unit
spectrum.b.unit = b.unit
spectrum.E0.unit = ebreak.unit
spectrum.K.unit = K.unit
Define initial values of polarization level and angle, fix the spectral parameters to their true values, and create the source model.
[8]:
polarization = LinearPolarization(80, 90) # polarization level (percentage out of 100), polarization angle (degrees)
spectral_component = SpectralComponent('grb', spectrum, polarization)
source = PointSource('source', # Name of source (arbitrary, but needs to be unique)
l = source_direction.l.deg, # Longitude (deg)
b = source_direction.b.deg, # Latitude (deg)
components = [spectral_component]) # Spectral model
source.components['grb'].shape.K.fix = True
source.components['grb'].shape.E0.fix = True
source.components['grb'].shape.alpha.fix = True
source.components['grb'].shape.beta.fix = True
model = Model(source)
Polarization fit in ICRS frameο
Instantiate the COSI 3ML plugin, combine with the model in a JointLikelihood object, then perform maximum likelihood fit.
[9]:
data = EmCDSBinnedData(grb_background.binned_data.project('Em', 'Phi', 'PsiChi'))
total_bkg = background_before.binned_data.project('Em', 'Phi', 'PsiChi') + background_after.binned_data.project('Em', 'Phi', 'PsiChi')
bkg_dist = {'total_bkg':total_bkg+sys.float_info.min}
bkg = FreeNormBinnedBackground(bkg_dist, sc_history = sc_orientation, copy = False)
instrument_response = BinnedInstrumentResponse(dr, data)
psr = BinnedThreeMLPointSourceResponse(data = data,
instrument_response = instrument_response,
sc_history = sc_orientation,
energy_axis = dr.axes['Ei'],
polarization_axis = PolarizationAxis(dr.axes['Pol'], convention='RelativeX'),
nside = 2*data.axes['PsiChi'].nside)
response = BinnedThreeMLModelFolding(data = data, point_source_response = psr)
like_fun = PoissonLikelihood(data, response, bkg)
cosi = ThreeMLPluginInterface('cosi',
like_fun,
response,
bkg)
cosi.bkg_parameter['total_bkg'] = Parameter('total_bkg', # background parameter
0.0016, # initial value of parameter
min_value=0, # minimum value of parameter
max_value=100, # maximum value of parameter
delta=0.05, # initial step used by fitting engine
unit = u.Hz)
cosi.bkg_parameter['total_bkg'].fix = True
plugins = DataList(cosi)
like = JointLikelihood(model, plugins, verbose=False)
_ = like.fit()
15:21:42 INFO set the minimizer to minuit joint_likelihood.py:1017
Best fit values:
| result | unit | |
|---|---|---|
| parameter | ||
| source.spectrum.grb.polarization.degree | (2.65 +/- 0.32) x 10 | |
| source.spectrum.grb.polarization.angle | (8.99980 +/- 0.00013) x 10 | deg |
Correlation matrix:
| 1.00 | -0.53 |
| -0.53 | 1.00 |
Values of -log(likelihood) at the minimum:
| -log(likelihood) | |
|---|---|
| cosi | 21743.138288324066 |
| total | 21743.138288324066 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 43490.276706860706 |
| BIC | 43509.13913959999 |
[ ]: