{ "cells": [ { "cell_type": "markdown", "id": "2bcdf955", "metadata": {}, "source": [ "# Polarization example - maximum likelihood method" ] }, { "cell_type": "markdown", "id": "f6c9219c", "metadata": {}, "source": [ "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. " ] }, { "cell_type": "code", "execution_count": 1, "id": "9946365e", "metadata": {}, "outputs": [], "source": [ "%%capture\n", "from cosipy import BinnedData\n", "from cosipy.spacecraftfile import SpacecraftHistory\n", "from cosipy.statistics import PoissonLikelihood\n", "from cosipy.background_estimation import FreeNormBinnedBackground\n", "from cosipy.interfaces import ThreeMLPluginInterface\n", "from cosipy.response.FullDetectorResponse import FullDetectorResponse\n", "from cosipy.response import BinnedThreeMLModelFolding, BinnedInstrumentResponse, BinnedThreeMLPointSourceResponse\n", "from cosipy.data_io import EmCDSBinnedData\n", "from cosipy.threeml.custom_functions import Band_Eflux\n", "from cosipy.polarization import PolarizationAxis\n", "from astropy.time import Time\n", "from astropy.coordinates import SkyCoord\n", "from astropy import units as u\n", "from cosipy.util import fetch_wasabi_file\n", "from pathlib import Path\n", "import sys\n", "from threeML import LinearPolarization, SpectralComponent, PointSource, Model, JointLikelihood, DataList\n", "from astromodels import Parameter" ] }, { "cell_type": "markdown", "id": "6efdce51", "metadata": {}, "source": [ "### Download and read in data" ] }, { "cell_type": "markdown", "id": "87669a80", "metadata": {}, "source": [ "This will download the files needed to run this notebook. If you have already downloaded these files, you can skip this." ] }, { "cell_type": "markdown", "id": "03845b0e", "metadata": {}, "source": [ "Download the unbinned data (660.58 KB)" ] }, { "cell_type": "code", "execution_count": 2, "id": "4b12c604", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "A file named grb_background.fits.gz already exists with the specified checksum (21b1d75891edc6aaf1ff3fe46e91cb49). Skipping.\n" ] } ], "source": [ "fetch_wasabi_file('COSI-SMEX/cosipy_tutorials/polarization_fit/grb_background.fits.gz', checksum = '21b1d75891edc6aaf1ff3fe46e91cb49')" ] }, { "cell_type": "markdown", "id": "d4df5f37", "metadata": {}, "source": [ "Download the polarization response (217.47 MB)" ] }, { "cell_type": "code", "execution_count": 3, "id": "444d31c2", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "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.\n" ] } ], "source": [ "fetch_wasabi_file('COSI-SMEX/develop/Data/Responses/ResponseContinuum.o3.pol.e200_10000.b4.p12.relx.s10396905069491.m420.filtered.binnedpolarization.11D.h5', checksum = '46b006a6b397fd777dc561d3b028357f')" ] }, { "cell_type": "markdown", "id": "780f93a0", "metadata": {}, "source": [ "Download the orientation file (1.10 GB)" ] }, { "cell_type": "code", "execution_count": 4, "id": "9527501d", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "A file named DC3_final_530km_3_month_with_slew_1sbins_GalacticEarth_SAA.fits already exists with the specified checksum (1b851c042acf4c909798e2401e9d2e38). Skipping.\n" ] } ], "source": [ "fetch_wasabi_file('COSI-SMEX/develop/Data/Orientation/DC3_final_530km_3_month_with_slew_1sbins_GalacticEarth_SAA.fits', checksum = '1b851c042acf4c909798e2401e9d2e38')" ] }, { "cell_type": "markdown", "id": "ccbe09c0", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": 5, "id": "32f94da6", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "fill() discarded one or more values due to out-of-bounds coordinate in a dimension without under/overflow tracking\n", "fill() discarded one or more values due to out-of-bounds coordinate in a dimension without under/overflow tracking\n", "fill() discarded one or more values due to out-of-bounds coordinate in a dimension without under/overflow tracking\n" ] } ], "source": [ "data_path = Path(\"\") # Update to your path\n", "\n", "grb_background = BinnedData(data_path/'grb.yaml')\n", "grb_background.select_data_time(unbinned_data=data_path/'grb_background.fits.gz', output_name=data_path/'grb_background_source_interval') \n", "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')\n", "grb_background.load_binned_data_from_hdf5(data_path/'grb_background_binned_galactic.hdf5')\n", "\n", "background_before = BinnedData(data_path/'background_before.yaml')\n", "background_before.select_data_time(unbinned_data=data_path/'grb_background.fits.gz', output_name=data_path/'background_before')\n", "background_before.get_binned_data(unbinned_data='background_before.fits.gz', output_name='background_before_binned_galactic', psichi_binning='galactic')\n", "background_before.load_binned_data_from_hdf5(data_path/'background_before_binned_galactic.hdf5')\n", "\n", "background_after = BinnedData(data_path/'background_after.yaml') # e.g. background_after.yaml\n", "background_after.select_data_time(unbinned_data=data_path/'grb_background.fits.gz', output_name=data_path/'background_after')\n", "background_after.get_binned_data(unbinned_data=data_path/'background_after.fits.gz', output_name=data_path/'background_after_binned_galactic', psichi_binning='galactic')\n", "background_after.load_binned_data_from_hdf5(data_path/'background_after_binned_galactic.hdf5')" ] }, { "cell_type": "markdown", "id": "b188a60e", "metadata": {}, "source": [ "Read in the detector response and orientation file. The orientation is cut down to the time interval of the source." ] }, { "cell_type": "code", "execution_count": 6, "id": "fa953728", "metadata": {}, "outputs": [], "source": [ "response_file = data_path / 'ResponseContinuum.o3.pol.e200_10000.b4.p12.relx.s10396905069491.m420.filtered.binnedpolarization.11D.h5'\n", "dr = FullDetectorResponse.open(response_file, pa_convention='RelativeX')\n", "\n", "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'))" ] }, { "cell_type": "markdown", "id": "90aa0f59", "metadata": {}, "source": [ "Define the GRB position and spectrum." ] }, { "cell_type": "code", "execution_count": 7, "id": "9eba638f", "metadata": {}, "outputs": [], "source": [ "source_direction = SkyCoord(l=23.53, b=-53.44, frame='galactic', unit=u.deg)\n", "\n", "a = 100. * u.keV\n", "b = 10000. * u.keV\n", "alpha = -0.7368949\n", "beta = -2.095031\n", "ebreak = 622.389 * u.keV\n", "K = 300. / u.cm / u.cm / u.s\n", "\n", "spectrum = Band_Eflux(a = a.value,\n", " b = b.value,\n", " alpha = alpha,\n", " beta = beta,\n", " E0 = ebreak.value,\n", " K = K.value)\n", "\n", "spectrum.a.unit = a.unit\n", "spectrum.b.unit = b.unit\n", "spectrum.E0.unit = ebreak.unit\n", "spectrum.K.unit = K.unit" ] }, { "cell_type": "markdown", "id": "6497dabb", "metadata": {}, "source": [ "Define initial values of polarization level and angle, fix the spectral parameters to their true values, and create the source model." ] }, { "cell_type": "code", "execution_count": 8, "id": "e942447f", "metadata": {}, "outputs": [], "source": [ "polarization = LinearPolarization(80, 90) # polarization level (percentage out of 100), polarization angle (degrees)\n", "spectral_component = SpectralComponent('grb', spectrum, polarization)\n", "\n", "source = PointSource('source', # Name of source (arbitrary, but needs to be unique)\n", " l = source_direction.l.deg, # Longitude (deg)\n", " b = source_direction.b.deg, # Latitude (deg)\n", " components = [spectral_component]) # Spectral model\n", "\n", "source.components['grb'].shape.K.fix = True\n", "source.components['grb'].shape.E0.fix = True\n", "source.components['grb'].shape.alpha.fix = True\n", "source.components['grb'].shape.beta.fix = True\n", "\n", "model = Model(source)" ] }, { "cell_type": "markdown", "id": "90d6d112", "metadata": {}, "source": [ "### Polarization fit in ICRS frame" ] }, { "cell_type": "markdown", "id": "130e6df3", "metadata": {}, "source": [ "Instantiate the COSI 3ML plugin, combine with the model in a JointLikelihood object, then perform maximum likelihood fit." ] }, { "cell_type": "code", "execution_count": 9, "id": "7f853caa", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
15:21:42 INFO      set the minimizer to minuit                                             joint_likelihood.py:1017\n",
       "
\n" ], "text/plain": [ "\u001b[38;5;46m15:21:42\u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;49mINFO \u001b[0m \u001b[1;38;5;251m set the minimizer to minuit \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=507214;file:///Users/eneights/software/miniforge3/envs/cosipy-develop/lib/python3.12/site-packages/threeML/classicMLE/joint_likelihood.py\u001b\\\u001b[2mjoint_likelihood.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=857216;file:///Users/eneights/software/miniforge3/envs/cosipy-develop/lib/python3.12/site-packages/threeML/classicMLE/joint_likelihood.py#1017\u001b\\\u001b[2m1017\u001b[0m\u001b]8;;\u001b\\\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
Best fit values:\n",
       "\n",
       "
\n" ], "text/plain": [ "\u001b[1;4;38;5;49mBest fit values:\u001b[0m\n", "\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
resultunit
parameter
source.spectrum.grb.polarization.degree(2.65 +/- 0.32) x 10
source.spectrum.grb.polarization.angle(8.99980 +/- 0.00013) x 10deg
\n", "
" ], "text/plain": [ " result unit\n", "parameter \n", "source.spectrum.grb.polarization.degree (2.65 +/- 0.32) x 10 \n", "source.spectrum.grb.polarization.angle (8.99980 +/- 0.00013) x 10 deg" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n",
       "Correlation matrix:\n",
       "\n",
       "
\n" ], "text/plain": [ "\n", "\u001b[1;4;38;5;49mCorrelation matrix:\u001b[0m\n", "\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", "
1.00-0.53
-0.531.00
" ], "text/plain": [ " 1.00 -0.53\n", "-0.53 1.00" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n",
       "Values of -log(likelihood) at the minimum:\n",
       "\n",
       "
\n" ], "text/plain": [ "\n", "\u001b[1;4;38;5;49mValues of -\u001b[0m\u001b[1;4;38;5;49mlog\u001b[0m\u001b[1;4;38;5;49m(\u001b[0m\u001b[1;4;38;5;49mlikelihood\u001b[0m\u001b[1;4;38;5;49m)\u001b[0m\u001b[1;4;38;5;49m at the minimum:\u001b[0m\n", "\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
-log(likelihood)
cosi21743.138288324066
total21743.138288324066
\n", "
" ], "text/plain": [ " -log(likelihood)\n", "cosi 21743.138288324066\n", "total 21743.138288324066" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n",
       "Values of statistical measures:\n",
       "\n",
       "
\n" ], "text/plain": [ "\n", "\u001b[1;4;38;5;49mValues of statistical measures:\u001b[0m\n", "\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
statistical measures
AIC43490.276706860706
BIC43509.13913959999
\n", "
" ], "text/plain": [ " statistical measures\n", "AIC 43490.276706860706\n", "BIC 43509.13913959999" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "data = EmCDSBinnedData(grb_background.binned_data.project('Em', 'Phi', 'PsiChi'))\n", "\n", "total_bkg = background_before.binned_data.project('Em', 'Phi', 'PsiChi') + background_after.binned_data.project('Em', 'Phi', 'PsiChi')\n", "bkg_dist = {'total_bkg':total_bkg+sys.float_info.min}\n", "bkg = FreeNormBinnedBackground(bkg_dist, sc_history = sc_orientation, copy = False)\n", "\n", "instrument_response = BinnedInstrumentResponse(dr, data)\n", "\n", "psr = BinnedThreeMLPointSourceResponse(data = data,\n", " instrument_response = instrument_response,\n", " sc_history = sc_orientation,\n", " energy_axis = dr.axes['Ei'],\n", " polarization_axis = PolarizationAxis(dr.axes['Pol'], convention='RelativeX'),\n", " nside = 2*data.axes['PsiChi'].nside)\n", "\n", "response = BinnedThreeMLModelFolding(data = data, point_source_response = psr)\n", "\n", "like_fun = PoissonLikelihood(data, response, bkg)\n", "\n", "cosi = ThreeMLPluginInterface('cosi',\n", " like_fun,\n", " response,\n", " bkg)\n", "\n", "cosi.bkg_parameter['total_bkg'] = Parameter('total_bkg', # background parameter\n", " 0.0016, # initial value of parameter\n", " min_value=0, # minimum value of parameter\n", " max_value=100, # maximum value of parameter\n", " delta=0.05, # initial step used by fitting engine\n", " unit = u.Hz)\n", "cosi.bkg_parameter['total_bkg'].fix = True\n", "\n", "plugins = DataList(cosi)\n", "\n", "like = JointLikelihood(model, plugins, verbose=False)\n", "\n", "_ = like.fit()" ] }, { "cell_type": "code", "execution_count": null, "id": "e3b06940", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "cosipy-develop", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.12" } }, "nbformat": 4, "nbformat_minor": 5 }