From 9a776b51e80b483dd5147b5868def7d39f9c2c88 Mon Sep 17 00:00:00 2001 From: Eliza Neights Date: Mon, 17 Nov 2025 20:11:33 -0500 Subject: [PATCH 01/11] Fix mistake in point source response --- cosipy/response/PointSourceResponse.py | 1 + 1 file changed, 1 insertion(+) diff --git a/cosipy/response/PointSourceResponse.py b/cosipy/response/PointSourceResponse.py index 22dcbf5d..60fd1e83 100644 --- a/cosipy/response/PointSourceResponse.py +++ b/cosipy/response/PointSourceResponse.py @@ -89,6 +89,7 @@ def get_expectation(self, spectrum, polarization=None): polarized_weights[polarization_bin_index] = polarization_level weights = unpolarized_weights + polarized_weights + weights *= self.axes['Pol'].nbins contents = np.tensordot(weights, self.contents, axes=([0], [self.axes.label_to_index('Pol')])) From a9f1d19dcb4cb7df63302959e3380d17a782c576 Mon Sep 17 00:00:00 2001 From: Eliza Neights Date: Mon, 17 Nov 2025 20:11:54 -0500 Subject: [PATCH 02/11] Add polarization to COSILike --- cosipy/threeml/COSILike.py | 70 ++++++++++++++++++++++++++++++++------ 1 file changed, 60 insertions(+), 10 deletions(-) diff --git a/cosipy/threeml/COSILike.py b/cosipy/threeml/COSILike.py index 87740fff..dffd8710 100644 --- a/cosipy/threeml/COSILike.py +++ b/cosipy/threeml/COSILike.py @@ -6,6 +6,8 @@ from cosipy.response.FullDetectorResponse import FullDetectorResponse from cosipy.response.ExtendedSourceResponse import ExtendedSourceResponse +from cosipy.polarization.polarization_angle import PolarizationAngle +from cosipy.polarization.conventions import IAUPolarizationConvention, MEGAlibRelativeX, MEGAlibRelativeY, MEGAlibRelativeZ from scoords import SpacecraftFrame, Attitude @@ -59,9 +61,11 @@ class COSILike(PluginPrototype): Full path to precomputed point source response in Galactic coordinates earth_occ : bool, optional Option to include Earth occultation in fit (default is True). + response_pa_convention : str, optional + Polarization reference convention of response ('RelativeX', 'RelativeY', or 'RelativeZ'). Required if response contains polarization angle axis """ - def __init__(self, name, dr, data, bkg, sc_orientation, - nuisance_param=None, coordsys=None, precomputed_psr_file=None, earth_occ=True, **kwargs): + def __init__(self, name, dr, data, bkg, sc_orientation, nuisance_param=None, coordsys=None, + precomputed_psr_file=None, earth_occ=True, response_pa_convention=None, **kwargs): # create the hash for the nuisance parameters. We have none for now. self._nuisance_parameters = collections.OrderedDict() @@ -72,7 +76,7 @@ def __init__(self, name, dr, data, bkg, sc_orientation, # User inputs needed to compute the likelihood self._name = name self._rsp_path = dr - self._dr = FullDetectorResponse.open(dr) + self._dr = FullDetectorResponse.open(dr, pa_convention=response_pa_convention) self._data = data self._bkg = bkg self._sc_orientation = sc_orientation @@ -117,6 +121,22 @@ def __init__(self, name, dr, data, bkg, sc_orientation, logger.info("... loading the pre-computed image response ...") self.image_response = ExtendedSourceResponse.open(self.precomputed_psr_file) logger.info("--> done") + + if 'Pol' in self._dr.axes.labels: + self._response_pa_convention = response_pa_convention + if self._coordsys == 'spacecraftframe': + if self._response_pa_convention == 'RelativeX': + self._pa_convention = MEGAlibRelativeX(attitude=self._sc_orientation.get_attitude()[0]) + elif self._response_pa_convention == 'RelativeY': + self._pa_convention = MEGAlibRelativeY(attitude=self._sc_orientation.get_attitude()[0]) + elif self._response_pa_convention == 'RelativeZ': + self._pa_convention = MEGAlibRelativeZ(attitude=self._sc_orientation.get_attitude()[0]) + else: + raise RuntimeError("Response convention must be 'RelativeX', 'RelativeY', or 'RelativeZ'") + elif self._coordsys == 'galactic': + self._pa_convention = IAUPolarizationConvention() + else: + raise RuntimeError("Unknown coordinate system") def set_model(self, model): """ @@ -215,7 +235,7 @@ def set_model(self, model): dwell_time_map = self._get_dwell_time_map(coord) self._psr[name] = self._dr.get_point_source_response(exposure_map=dwell_time_map) elif self._coordsys == 'galactic': - scatt_map = self._get_scatt_map() + scatt_map = self._get_scatt_map(coord) self._psr[name] = self._dr.get_point_source_response(coord=coord, scatt_map=scatt_map) else: raise RuntimeError("Unknown coordinate system") @@ -227,9 +247,39 @@ def set_model(self, model): # Convolve with spectrum # See also the Detector Response and Source Injector tutorials - spectrum = source.spectrum.main.shape - total_expectation = self._psr[name].get_expectation(spectrum) + if hasattr(source.spectrum, 'main'): + + spectrum = source.spectrum.main.shape + total_expectation = self._psr[name].get_expectation(spectrum) + + else: + + component_counter = 0 + + for item in source.spectrum.to_dict(): + + spectrum = getattr(source.spectrum, item).shape + + if not 'Pol' in self._dr.axes.labels: + this_expectation = self._psr[name].get_expectation(spectrum) + else: + #polarization_level = source.components['grb'].polarization.degree.value / 100. + #polarization_angle = PolarizationAngle(coords.Angle(source.components['grb'].polarization.angle.value, unit=u.deg), source.position.sky_coord, convention=self._pa_convention) + if self._coordsys == 'spacecraftframe': + this_expectation = self._psr[name].get_expectation(spectrum, source.components['grb'].polarization) + elif self._coordsys == 'galactic': + scatt_map = self._get_scatt_map(source.position.sky_coord) + this_expectation = self._psr[name].get_expectation(spectrum, source.components['grb'].polarization) + else: + raise RuntimeError("Unknown coordinate system") + + if component_counter == 0: + total_expectation = this_expectation + else: + total_expectation += this_expectation + + component_counter += 1 # Save expected counts for source: self._expected_counts[name] = total_expectation.copy() @@ -319,8 +369,8 @@ def _get_dwell_time_map(self, coord): Dwell time map """ - self._sc_orientation.get_target_in_sc_frame(target_name = self._name, target_coord = coord) - dwell_time_map = self._sc_orientation.get_dwell_map(response = self._rsp_path) + self._sc_orientation.get_target_in_sc_frame(target_name=self._name, target_coord=coord) + dwell_time_map = self._sc_orientation.get_dwell_map(response=self._rsp_path, pa_convention=self._response_pa_convention) return dwell_time_map @@ -338,8 +388,8 @@ def _get_scatt_map(self, coord): scatt_map : cosipy.spacecraftfile.scatt_map.SpacecraftAttitudeMap """ - scatt_map = self._sc_orientation.get_scatt_map(nside = self._dr.nside * 2, target_coord = coord, - coordsys = 'galactic', earth_occ = self.earth_occ) + scatt_map = self._sc_orientation.get_scatt_map(nside=self._dr.nside * 2, target_coord=coord, + coordsys='galactic', earth_occ=self.earth_occ) return scatt_map From 982fa65b92e36783d150c37566f2f9f9b74c4dd8 Mon Sep 17 00:00:00 2001 From: Eliza Neights Date: Wed, 19 Nov 2025 16:03:19 -0500 Subject: [PATCH 03/11] Clean up COSILike --- cosipy/threeml/COSILike.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/cosipy/threeml/COSILike.py b/cosipy/threeml/COSILike.py index dffd8710..9f55e5a7 100644 --- a/cosipy/threeml/COSILike.py +++ b/cosipy/threeml/COSILike.py @@ -245,9 +245,6 @@ def set_model(self, model): # Get expectation for point sources: for name,source in point_sources.items(): - # Convolve with spectrum - # See also the Detector Response and Source Injector tutorials - if hasattr(source.spectrum, 'main'): spectrum = source.spectrum.main.shape @@ -264,8 +261,6 @@ def set_model(self, model): if not 'Pol' in self._dr.axes.labels: this_expectation = self._psr[name].get_expectation(spectrum) else: - #polarization_level = source.components['grb'].polarization.degree.value / 100. - #polarization_angle = PolarizationAngle(coords.Angle(source.components['grb'].polarization.angle.value, unit=u.deg), source.position.sky_coord, convention=self._pa_convention) if self._coordsys == 'spacecraftframe': this_expectation = self._psr[name].get_expectation(spectrum, source.components['grb'].polarization) elif self._coordsys == 'galactic': From 9d27cf7c782ab669b0e69ccf6c959559cbe7b65d Mon Sep 17 00:00:00 2001 From: Eliza Neights Date: Wed, 19 Nov 2025 16:16:06 -0500 Subject: [PATCH 04/11] Add example notebook --- .../maximum_likelihood_method.ipynb | 273 ++++++++++++++++++ 1 file changed, 273 insertions(+) create mode 100644 docs/tutorials/polarization/maximum_likelihood_method.ipynb diff --git a/docs/tutorials/polarization/maximum_likelihood_method.ipynb b/docs/tutorials/polarization/maximum_likelihood_method.ipynb new file mode 100644 index 00000000..641264cc --- /dev/null +++ b/docs/tutorials/polarization/maximum_likelihood_method.ipynb @@ -0,0 +1,273 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Polarization example - maximum likelihood method" + ] + }, + { + "cell_type": "markdown", + "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": null, + "metadata": {}, + "outputs": [], + "source": [ + "from cosipy import UnBinnedData, COSILike\n", + "from cosipy.spacecraftfile import SpacecraftFile\n", + "from cosipy.threeml.custom_functions import Band_Eflux\n", + "from astropy.time import Time\n", + "import numpy as np\n", + "from astropy.coordinates import SkyCoord\n", + "from astropy import units as u\n", + "from scoords import SpacecraftFrame\n", + "from cosipy.util import fetch_wasabi_file\n", + "from pathlib import Path\n", + "from threeML import LinearPolarization, SpectralComponent, PointSource, Model, JointLikelihood, DataList\n", + "from astromodels import Parameter" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Download and read in data" + ] + }, + { + "cell_type": "markdown", + "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", + "metadata": {}, + "source": [ + "Download the unbinned data (660.58 KB)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fetch_wasabi_file('COSI-SMEX/cosipy_tutorials/polarization_fit/grb_background.fits.gz', checksum = '21b1d75891edc6aaf1ff3fe46e91cb49')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Download the polarization response (217.47 MB)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "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", + "metadata": {}, + "source": [ + "Download the orientation file (1.10 GB)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fetch_wasabi_file('COSI-SMEX/DC3/Data/Orientation/DC3_final_530km_3_month_with_slew_1sbins_GalacticEarth_SAA.ori', checksum = 'b87fd41b6c28a5c0c51448ce2964e57c')" + ] + }, + { + "cell_type": "markdown", + "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": null, + "metadata": {}, + "outputs": [], + "source": [ + "data_path = Path(\"\") # Update to your path\n", + "\n", + "grb_background = UnBinnedData(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='grb_background_source_interval.fits.gz', output_name='grb_background_binned_galactic', psichi_binning='galactic')\n", + "grb_background.load_binned_data_from_hdf5('grb_background_binned_galactic.hdf5')\n", + "\n", + "background_before = UnBinnedData(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('background_before_binned_galactic.hdf5')\n", + "\n", + "background_after = UnBinnedData(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='background_after.fits.gz', output_name='background_after_binned_galactic', psichi_binning='galactic')\n", + "background_after.load_binned_data_from_hdf5('background_after_binned_galactic.hdf5')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Define the path to the detector response and read in the orientation file. The orientation is cut down to the time interval of the source." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "response_file = data_path/'ResponseContinuum.o3.pol.e200_10000.b4.p12.relx.s10396905069491.m420.filtered.binnedpolarization.11D.h5' # e.g. ResponseContinuum.o3.pol.e200_10000.b4.p12.s10396905069491.m441.filtered.binnedpolarization.11D.h5\n", + "\n", + "sc_orientation = SpacecraftFile.parse_from_file(data_path/'DC3_final_530km_3_month_with_slew_1sbins_GalacticEarth_SAA.ori') # e.g. DC3_final_530km_3_month_with_slew_1sbins_GalacticEarth_SAA.ori\n", + "sc_orientation = sc_orientation.source_interval(Time(1835493492.2, format = 'unix'), Time(1835493492.8, format = 'unix'))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Define the GRB position and spectrum." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "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": null, + "id": "e942447f", + "metadata": {}, + "outputs": [], + "source": [ + "polarization = LinearPolarization(75, 50) # 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", + "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": null, + "id": "7f853caa", + "metadata": {}, + "outputs": [], + "source": [ + "cosi = COSILike(\"cosi\", # COSI 3ML plugin\n", + " dr = response_file, # detector response\n", + " data = grb_background.binned_data.project('Em', 'Phi', 'PsiChi'), # data (source+background)\n", + " bkg = (background_before.binned_data.project('Em', 'Phi', 'PsiChi') + background_after.binned_data.project('Em', 'Phi', 'PsiChi')) * 0.0016, # background model \n", + " sc_orientation = sc_orientation, # spacecraft orientation\n", + " response_pa_convention = 'RelativeX') # polarization angle convention of response\n", + "\n", + "plugins = DataList(cosi)\n", + "\n", + "like = JointLikelihood(model, plugins, verbose=False)\n", + "\n", + "like.fit()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "cosipy-fork", + "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.10.0" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 3a6bc47cf9120a9938c31c8d94246a2c7ec0af3d Mon Sep 17 00:00:00 2001 From: Eliza Neights Date: Thu, 20 Nov 2025 12:39:23 -0500 Subject: [PATCH 05/11] Add unit test for MLM --- .../test_data/polarization_data_binned.hdf5 | Bin 0 -> 42392 bytes cosipy/test_data/polarization_data_mlm.yaml | 14 ++++ tests/polarization/test_polarization_mlm.py | 72 ++++++++++++++++++ 3 files changed, 86 insertions(+) create mode 100644 cosipy/test_data/polarization_data_binned.hdf5 create mode 100644 cosipy/test_data/polarization_data_mlm.yaml create mode 100644 tests/polarization/test_polarization_mlm.py diff --git a/cosipy/test_data/polarization_data_binned.hdf5 b/cosipy/test_data/polarization_data_binned.hdf5 new file mode 100644 index 0000000000000000000000000000000000000000..632dc42bf6cd19835137a7792f5ca326af4dcd79 GIT binary patch literal 42392 zcmeFY30RZYwl_{|TiQC%7A;4aY^|bILF>Sjp&qr=s-%b*5lB?3NErf#7(z<3RVm zjEagfBq}PC5D_H|8B`_}NMr~BLh?!&0tqA}Aqkm3DA0SZ_c`Z&&wc*iJ>UI*Zt-1v z`0cgVUVFc5?|1E;v~PEPv-ramAKJ`>_usd9$L6CyqVo^yT)1{;HfKIPUw(8h^41*v zZH^YcJsWw;X6B!b&0ps7Kc3_NI_}_s-Mek3k@NiL>0C9N_pi>z-k;BSL;s)Bz_+`1 zd_U9S!*ji!pQX>{=GnY5PDP>VA%`%>!-~YvC-c{Ak*+&DyrF{{X%_otnvid%-ND zg*GeyNWV49b~-*&{+7)TCr(BmOE@uYw`arqbC!0&+Y6>P&SuTjxVK=Q{hx2pT!sIS zNuWP$*IU!l%=X2`=IeL=c?0wB@g}c-j$1ib{JqaV$IZWoUwi*^+?&!j8u(Xe;PrKE z)>_!uT=?MMOS5S%xp|I$HAlD1(K$08|0brr=5R$*3SQFzu)?8@hl>4mdrgD*_`?#HY3Nv-OFS>%{B~)cE}X_q&9yX>H?xk+r)TPa^arQuJm=}ebm*QXe?DKg z_nT%I{*y32uU@+R6VPwFXrf)Il6X^e)+ZEPT0Ob=EvNC%+8hB@wb@{I}#(&GxYE0$xM1?e;tQ& zKQ>ovGdGQ=cR{luIx_aevD5#?XMSIDFf!)E3}SveFyG!EjZEjy*PqAFw?Ai=Z04uc z>+R3)Gv@2hw?A)A=jGoWGb`Whhpx?^az~vw^hd!w|JU1JG+XT-wm(G?fz>t@Sn{OMbfNXEB`g(y=nDdtbu=Yj&XR-Zp@3fX)f}A>KxJ!SYerTkAg$ME=NR^GR# z$&5d9K0VXLMK)hgKWENF-kbJ$zHfsLKNdaXkp73~C4Zo2{nGP%W~R*Re4>vXKN0=E z!RPPgnTh|Sdr|-0^K$;ZnbFxynD^_>$T$<`{q^&H<9UDlykB_U&pz*0pZ8zS`^^&( z;|-gqeZq?t+~{JD7I#-X-Z|aW#LgH&<@6uf^-Jp8AJ1UQKd{;KzRmq3R%_1_*LOAt ze?n_ev1dUsc6pUC#<6z~mb|~?59s0<=#4d9R^*w9cin#dHF0CKpmO}lEL^MXSntvI zY%%o@=*NFs`H{`0Y`>3n)U2P;?+$K@{N%)k(@mq|AH@i0gl_Fgs@68vUo*5?M~aJgcy26m@L+9kN${|j--8#`Umf9m_G45dep?PBH2??AKY|3&0=hU@ zdX=1Z8n`7lAmd0!EKDL1wwD437SkNw4ogI|^8HSCZ~Ew+uh-kIkJziSvPV;X*$a>P zEUxe?55k&q-vurvL(&rEsh0zLU}sbbZrtrjtV{=A{BkY!S;+;yFV0@IeJrfoc;D`r zYaQfQ!V>*gqfyW38J6HlXXtz8blYF8QFgt1D|+`{(f;0R|Dzn9ExR(Ty*FHuU7Mdx zyE1$dh7M@;-bVkjHUG=3Be7L%dHka+Do|NKN6EtkSwEg9B|-9Now-43u|5i8VWKt3 zJasf#QoyUm5h+;!U1#pA(1Lr7bX{CoE{!t7JMjn`bwszZ)R@!Gghl$_-iB(=4hg}) zI%;&0trs5-LqB>HpQPHyehU-rusw7~cH6TJqID-=UjgJR!v-0{X zy~n_uj6>))U+YN>#cH$h#inl(5s$j}97x%*C?(*pZ+~&@{hinPE&LyTV(N=S<+07y zVP_lq?KsAiHP&9yf>(d0e8;svQyKlUvAkr~m++Ad>u-_e!!lGi}f4P}J}aG@eA znAE@KBL37NIkmVC{AD$%__VHYP!xQ{xF!>S;*k8}xAJqveIHkxb~`xaUiiZK#w9cR zOFZNflli3)7o?3}W88KLy!^SS`3Gvo3zy~}kOA7IP;GPJ*p=X(RTYG%$ptNuFHFv@ z$A;IQAaaJ+s@*ejhg!D{`(K5&WPCgP*TIomiSHq11j~=*###@=l(+&xYk&_FkQ$Q6 zJ45bU8_ol(M4#hGZO-3dhjjMBpR9y@tKV0&lPAs>4LKJMg#`^c7yl4jGUOZ3j4imT zNDN|;uC&-aYH@wkvhhj__)$y1m6q^FeT#prhIXEADjbv-yhu;Iw1WyI((kFGas1<)^pgr2U^JZeM=$!ofSDL6&T^sK(9>;JyXq5SM~C zCUs3V-GzFBfyv_?$aKVZ<_Z;7j1jk9px23FlAwcicS+j`(FSZs!Cmfd>oUSR!;#62 zFDV>DjiC0THs)j(E!Y}_e}M?VuSU}&#!`y4U&AHPEjd_8U>^HhiF*z^x-8cY?E)V5 z-twZt`@AZ5sMpC4w>?Y9xSqMeqemYh!GqkmO}bHh7!O@2_1-V5hDAaBqudL&NjiFE z9W`!uF`R-&Zg(F>+ApB;t1R0kOu#S(LHU_ zcXhkf?;_Qzad8kzYLXPsCRuGf^X-mqP|#|Rsze~aH=!n~qMK|7CvB7E4_Ipt?0 zd|Y}$+2?s~V2nieYJdFK`6YB> z@ASGwFxU*&y|rc9LN`pW2h*?{dTYl**3n=7T6X-8rEhL+n~vKNJGt)$%^+tOzQ16l z%P{jp(YMVmw$i%dH#1U^Ynk)Z}?ySoy*&&@@6|e zvsL)JZNkz24s}fX@!p-iXFa!Hdi_55{61!${-?k9`49KNmv5Zy$jaCHbK$l0{~gKy zsE^qpnN!cSfAPP1A0BpPwugT=#D8R5`@f^(Z>0D?s)1>*7XNvyW57o{+rHo&^-bB* z%`sWo*LT;qzV%xYxX{yK(-A3r!=)bq<`hY|OoY6yfQeJrD+hKEYd&%H7fUkrLo3;o z7L=_%@e}G26>b&V&J=$EYHRtKz0evurF)v?!1|K0d@9Ke*q|EsXE`R;W|q>%{UW)0 z@b}sfm$duy!C$qBPCc|6ysYLqWqY8zR2Q8FvrVWa*(8qkSA3S7&`T0TXm44h_sXQp zzky+52%AUU;FzeI9}tq}P;H1WHsAc{J7d%>Lgjaa+V6C$gdMxk_pitaxAfNwNtNxU z@jKO5w5WV=--z*WJ`>!Apq>I~?LbOei6WSV>>+9MnEX?uhc26~rM4~OmjPaAwj276 z)lI5?t|%dP6@-{;r0gKq7|l<{GwhrATx}%#5%UpY7Z&q^G5i$t1FPkVBm$~uB_qN* z{j$@8ie+E6@-J6>qkD8^_+U0a%@BUea=nl}d3~#{-2{8+pw*^4bX_sN*jOjJ5a>Qh zhmc`8q)NkQ#t>mT&W>QmRTQg7NnCIJ3-E9?j#9I_z?)52Gk)hA>HuZMcuobinfzth{GayNrdfjZGL${0-J&4MLed zVhRI?C#@5GS)s!Q+i2G24-yeV-k#XHuNUoVYZ)*ra(35jw_AIEqu(Z{*cDfwam6oH z9d$HtsX8XF97IA$WIk-ZzRKVH6o-s3g>RIf$a}s@5}t>Y^C?i26FJ-T#|^%z6T+-d`3Y5S(sI~F1dDJ%!Mu>E>l#Nq0%stl69iQi`7{rZ zXGM6))067FGC+P2V6GF@49ewrX^f_!m%$*XrLzfV$~Ma`I#hUw$b<3J zr+aeEzYSFzUx9U$L};n6B~52M=^Rr=^y_6-vhAA`MU1;1XK1F1s`C8lM}&-Z6&mj` zlvk=_tH`GvX`^|GnVd!;4?hE`twA*`rZp(K2u1GIrTznM&Kng=9nz`sMSJ~Ad@+^y z9y#3#&p_o8Vl?|1Z1L8x4n_!aaG$cRL!8&VlxB$}o5~977|Okr!yT9eL?E(s{awyf zL+3IUNT#KDXI;~`c1)&r80a}9tVf9DD&I?8S(78ug7jml8aL!X{XGC56Dp9!^?6Q@ ziq1ENs_Mh(`*bl7P6N(}h#eS}IJMZb)S3Bg-OZtLuGF8Z?yLTh!zB*V@kQ>&DE7Ve z3ZbHlcFm|+aPZpFifM1%Yqpoi$bWT)dV$3|aCDLWW( z4o&nrPbcwgWosy)uKl^@=Wt~hRra(cK5rSV7=Kk_072C!y(z1%R$bxLGq1qoe{+); zCT$W&m`~~{enAl$$f&W@2L#a{Kw<6ohFDI4B{5i8Iw8oX$Dvxf?ymM2(S)R3f(+5~p7?t8s52HIk|W|SGs{<#1~Kt`z&5gb{P+2RDtH3xSP zj!}Jt9)O!P&Kl}!`HN&Mzqzm2A5E}^Lf>mrwAWl;Aemg%3e-HWV_9J=IFU+3UXz?; zF5?q$TUph3oY2R*9vW}&5r4maukL1hYEa*8U7J3vFIG5=jjzQ(9j$9)+iSBT@`jJl zs!Kp#d&mO^!p_&e-JG*DdU}w%n$J6JUbM4~5pl_6wVn*Xc-pQ6W!!dn>xg51blR=8 zJ%u`0mYKmv>I~xuHXsu)SwNJ9t=P8SpB)DUWLfwy+*f9{NJhpl=Bwk>3*itY66f#5 zJr8pj=+L!I1qz4TpzC9KbOQ=u3^w^7LHqs|EPVXH$=wwt*% zr*!HL`pmfy9p4+CI`%86XQB0h^<8Ulm4_&cIpKl-)bi9E!XPy$dNkAt*M@~|h=!r8 zyaWzFP*W0Re0t*O{RAFqG=h;-lxZ2>f?#IAY2u-%T!PmSnCDEv(7g27x_Uv?m8Pl& z58u`z+<3EC(wAgus~HH-)jVs$A>ySGX4edEZlG){P|U|3AKHeDZ?+zkoa;a$r1w?# zbl*8CiWDb+sKp+E3aJLqqw3|Qkb5ixqM?L( z0j}5pH3oq>RCEyu;G`x!IwQA8$vzR%wrX5mRBNYida6pe&3aL zigZ&8H7y(=!D)6oDSDoUXh17SNlo?kO*OeG81AQJ`}2&a#sz2TOEacVmb&)L%(qV zxz{%&D-jLV4f-!*6BoJ#YiW&|00yZQ)&Cdx=al%gd~=6atm%ja10gQMCsa(cR_J2CAlPe37mMJojOi2M&HVRs! zEv*`(Yk~=vQqv=d4s{=RD#J(=&6BFqe2wfOSzA%7|21pGBUJJ|<)oXUzQUDBtNyP3 zB1~RVU~tuBZtyjZB7&APo%}&v%NXh8fb58r74v%PiRo1`bb&@8PZ@ zOYDZxCH$vtRBhT+nwLCfmrKH31)~t2%-c*-1-oEt9A4BVH#CcwN5OFRIXJZcdaCg` zVFSBwXBuIAjG(N<8+f0UL(&xDGTm5ixUpRNVB$rfld;$;s0A)g)Y7KLQ3dm=4w zi!LZh+XQ-$=|@aeW&Yx!IKdznmum8V zkmZ@z(?ik^6Ca7>%hmhOAqlbrK(3!m%qph|g$iD4O&QB4%L2XNAj*j7ueGb+!&U4# zPH5O>k>1L~llt$WhlE1HXUIu6k34%_uIY+y^tq#K7_ki^>0^|l*}hO!buF>f8sp#M zsbAoSUFth?6M=c|V8BQ3y!v4A;zf%VojrDO$F8OOx3}v^D_nPUy8mK#jx|vBd(P4V zk~9IX%VX*i40@ee=21AM5}-hfrDQYc+SEpCLoTq0wc@f)bJV1IX8`7GNKr}DxBd$`3C^j*Qp?O63IEUTi z=#l}5v#03_MRgDQ3^4>z;cIJl%;pmnt%5j#ZZV3<2UeRGXXs4y#D~@b5nq*K7EQ6c zG<~P25V+hgodi8)(z*;-TNyOyq9l~qJtc$7Ag54m)vQ`r| zx!ks{#HFZN@yj4#Ib%^TX%f3cRp74K7dGlaWCNdAz5e11S={K#SacUDtF%m1qpE@G zFCfOtiJDDN({#PRb*IS7@l?gDTzkGD!P=_UAB3v&$cpB1q6VYbrAgqFiW+5o8Ca^l z!(|UlHCByvqK zN?WEcg@NrQxf%9DaK%(OU2&J1^cR}BdIHJ&7@}Fm8+Z{M)==pYV-=X6Y^OIyf%m&9 zu5hcmxhamMie`XQNj}Ky^f_QKyI{18cJt)i>FE18S(&p$0ncUF8sq zF)6J!yWH=lCL(y*cT=Tw^~Bf8 zjIVP}p$pl8B;l59R;^L1Q$T@!xqNi_C95Kl$Mt}5r|OPcLUIja(IU&{wE!$d@gB+P zGL$l5enm}*m>53_>KVwAsE=wIa+pc3w$>H*EQ#t0wwUN_BD}Of&RgnksQXm)R3>v< zs=@HpS_Rl(D5GpS95xK|AnC!jqt&0vEcDW&@u6||$t0xB)-=l#oNr$H7gju$n3ey|C zFCw+W7y7_hqLQ+MdtFkJtXzXfzJNP-0SwMjgIGpGAI;T$Rg1X-;MY+JKVXr!&{&G5 z4P;`)k-JL4L|v3C=LyqkQn`korEB)-%fbArOVQ#!CE)0gdwVZ8RV@a4^p^NB)u&N( zXQdQZR#X4gEbez+nK*=t^m>9%O{m^FWgS9|uJp4?&K({PbQv?|3h=EJU8e8tb zrwa2{Vy0)`6yjW?Z&pZ~cQjgL;>kRu1&<5k={dk%2Rcr8Za9FrE>@Xi`4G9@ zKbY|V*aQRzT%D^Tkcu81vNMBNuBJ^RirU1~s9k$lJ zTE3N*Oip2W3{qLZ`mPgDr2)cY2=@|s%vQi#8=GV9II6vGn7-+egVaBr>TT9S4hva= zSVTa`*s@xyo?A9I&=P1TERcz|(@aSO{=}aMP&8L#kmH}qGqDxj=#<2FoGdFr2 zOhh#BYtJ3rvGmveC4Vdl*(uLfe*Wc>>>h=6E9B=w{AcR^TgByh*9Rj2j{<+q^^Q>q zP|@s>%OV+D?H_cOG1^l-+)4cA!8X7|y>Q2pQIyI2f&+fuja-sr*1u#u?dG}e)P1yu z!wRv=cY0eIf93QHM@8cZB2GT3gf4&y$Y+M`!Vdvstz6b2bD5S`B?Gk8()E*l*HFBs zlIEQ~=;H)+h~IX+>Nd-gQDPV;bfNCRulDA_NYye2gH!BoV>eha?){KKZ}$u)43q(- z`hytS30a`|5#F=xc}6YK`l6}0Vz04Bh##Y+HctBp#(g*GqOVn{C+c1r!O2k2*-}=# z;;ikY=p5i|cxl{fm~#EzI;L4=7*t}ns&Z}ucjVh}6@F@8KQZeT$_PD+Jlkt8(KSde zWeI-IPRF6_$$SglR9Cv_sXMTt$BP zHOS{p0P~^Rp(m+nlX4>mFn4d^Z~h&O`6PlZave{e67`VdAwYd9rK5fbaoP@|n;h3S zyNr1}(3Pu~qev{f_4vsb6qNCIq`YJ*hnBZsd{~1}h2D1MB6C)qkzm>KV7#ibhSR3n zPDm%;LzIn?ck7q}xS)V|Q!)@+h{Pkh%8_c{RNN!H;xL?vThTh;y`OSMv;vM<0k_1s zn>*L%G);ZVAHi+nFb#MGzJzM{Mt3|kCj@tgUq;2`A}bTkKd37_B11FBo4s2~M6BDG z3{RLEwH4AwpoTEcm=6)8q>UlMZ_JnBhHrp@Fob{|FQpBYs!X)x+R(g;g1fqgbK^!& z)sLcX@=aU|)@6&z|?m>r?IPQxP#?zmM=dTiH5d)hPo*?4X9 zjH~ILpu77ne!ly|y7vzZ+AFQdOt*K=#4fbkv{JFME~9Tz)!?}Ey62xVV))$A({);- znj2H^k&?!~xy6;WYRRbe%cSH?yZj1B{L}||GBMUA*>Njy&-xMb+(4FQyTd9O`y;e# z`|xwU+SECkRozCFR9>_60ftwGO-8h{uG|ugb-4|DX1Pvhb-Q9H>y7t4LR|A8->5$; zA(X<}yBtd?+<`B9vNXpeuP}G1rPG%iTeip_)SlN})HHYvhF(MF)wW5CwyfER@`*Gb zgj6JXTBA7>RHY>#6x<~B$rr-za=%i!L?;+o~wQCpC@cW`IyIux9oy>Iil6vJ13@s4dFp$T8vh5K6`7pPtT_?^rQ(Kt^L9^r) zxI}WoLuJ@3N++xGPUxHjTH6AA%okhB)s=>v%;f4dSGnFiZfzQ`!=CyyNOTzJNaF-0 z)6zH{-Y;@g4RCcQ7AKBS@AK3M!bSU`L$HEj)Dqrjq!%7k9u4%7PL*Cqw&z8ftv=;7 z#cnPp-vCf$QZC(umWw}LGK`i>69;~(>6Qin{S(|>Mt77s8{St}4e1cN+Sk=4gWX&X zrV>|L8_YxXQ?HylsWot0p1qxCXs(}EQHRu9U!xu|C1Cd_l=WLB^@#{%f3Vq5PIVaKl7iuBT}`aE4f5cv_yc-GZZ@g`qFQg#{r*iKMS{=hR4rfOv?g zaf59q#2757MLsY}P4P4!Vh;4nvzq5ffCUcPIws~V+n*04vP z_%Kj)HzP{w2->X$ahIUte(<2hnVQt8PD6H;y2d6I$7(vm5x6nP3W$8S2rm5&z`6*W z2CFRt)!J3(4VQ=yN54Wwb6V+1NglC^iY1u)`U?CNg|kwTUFbeIu!qWF`V6wGs({gs zp~YD-5Ndk!U6|yN&?vrzI#{4cGf3|Gn)kUzK@3$2w;^FdmGZWC+D@~~Z`2Qwn^3{D z^>QVMv)3Z)^1=72tGFSlj5rt04yPCS#2j|jNS)zZf-e?xi+m7D~gEaK8 z)XQaP*pgIfNtraBruH4!Rt{=y(;d@Au(?I*xO*Z~xA8mCZ~2&KdzHoKGv&nv`WmPn zvYM#v(QEkm&=|`rP+ftA7a^&_{+2@=CAcfL1NUgdufX{T0{0QO)mrQ;j5b_Ai6Vx) z^VpjHYHhx?H!E4uh83b;cp&m5AgN+94?FM^@q0iM<*e!e3v`cF-BdLip*a@TTIXHJ z_O%yURercv78yPbQv)xl-LLDf;mOQ5t<|o;B%}?sy)U-}lbd)VjDzdQ^f?z%Kxk3i z2mx{s>USjt5P~%x5XfbQzT}(uCKCsvE=b(yFD^sWKhdqmjr4V-63aEC)x1%{pl+8W zQaIFMC+v`h3xFKp47o~C6Ojft>n(DVzj3S*`NBG(HW_+`#mu$tT^o^;8Oi&cR3BJH z3ejEe*}V`!wc#ewne@WV#Q-BrcEXZ-&_#8R9*NagR07m6m=;ns;X1IE*uAhJaY7V0^7fOM`3Fw&}@hbZc#Qi%YrF>;i~Ulp364^rba z#3%}92dmN6%8J{6Iko>YyZQl=eH>wDA|mg4)9J8jw`NCrM@~oh9f7~m%5MJLD`!i} z&De+i^~%hzKl-Im(!Acu1L)x=fi}og8OuRm|O1VeeQ6G`UweO(L*3(lb@_ldZixQOX(bdeJ;t@b?k~e;`VA4lA3}jbgMZ{8GM&iww@Z^ z7BO-ok4r>-@&tYHeTytjj^wuyzjki1H(C8By<! zzVfV^6E59CmDD^De3w8R2T8joe>wkq{^<4OCn@Y3C`m+KP3%|{;odIQZPd%=!3b@* zk%_xOYe|uzTn>G%_he779rPHgzZ?(4Y=Y?CP0sc+p7doaPde?l z+i$!6&%qZas7}4%Y6Ydr;kn1aq7;Zp|53-+22tbOX z_1C5%ScuC+lesf41*`I73&wmq`XU(_ipBy3H{WYC{7@Y~@u5o#&`En$PAz}DB~M_m zC-&-DDEKYhElR7RH}rn9E8z!I7gaWjy=uIz#3dl&sfSZbX-`HiXHi_u$PN8cFGqJ> znnZLrWjpkxYb|lyusy7%tg+&g>()@1`nXdoEELWiEi2J)L3U!Y;i16Yy=pK04arV* zx1k>Wp+p_~Hz+5@ph_Opt(KJS;uMro_U4eQSha+7Xbq<;9~WJdAPkFOQ-jM;b%3GA zYdn?Wci|vg;SRl3E)?F8Y*2LM}!&i$!f4BvLjXDZ3Wwg{S5%s!U&4Pz`YOX0AZM)Sc=pQ zeHZNow(Iq4@vS+PSh|qykPR6FUTMIw=ATijg}NN&b-3(tVq1!_m9V{)ZZYmt%d}># zd#vnjHblCr^CwUnU2;-aBx_5)Eu_9e`MgzNPVs&?+HLmKf77wQE|O|ffLR=SHNMxQ zWk4P(hCGDkAP2Cn2e1iCfDL^1sGv$fI7PuKLz)`7y$!y4QGjVYYyjU5?7Tm=zSNpc ze5}Bs98l_=*1QJGZ=}ERiQN+iY}K)O0>xe3@0h1Q17VOH(qX$8mzb(j^;JK$gN(x3 zMten-I|Ex-Tu0igpWwO>^-vo3Y(%?T@=wjA)SSsk(HUbAkDDCeDdkPRPpXBy%i4o_ zr7|uK9d3qn4sT(59LW>%hvQXuu8bI4rax#@jM!WjZA~uettDGykkjU~S8&N`_o82p zJ^tI~qxYW-AKet2Tr4PJq_C)@2Yt4vbgM9MfodZP&F$#!&kS*(}j}lYusj%eD z%GP@J8dlyEQQVY9f{hsx#KJTr%oHou*}Lvb1}O{Dps-{gB{|KlkLGsPm!Vj5wpNy& z)^9p_e9Fb&bWBXPZ!rEKR@m1qFUzt=?q5{EIIHH=bKTCWo72SRGsjJ_i36IfwEn9= zpm|7LSD)~`8Je?JPFr)9fZ{`gN$d zdd(4F^2_syW9<&xg|Wm7?n4{Hx~-~ZtYjqd9txCaVANQA+AQzGEY-@T9-~DV9#dAF z(#r4rNgs8`MSf9Z)30Cb!g*MPZ8)?Fq1=tQda2?qpMDcx*^d)Y{;>gtp z8D$~7V4^4Z?oPt)WR5GQTCEq)KIinH?468-DY+& zPx2S@mricbWn@P4&Z>LrE>A&d^_tjZ8Cg!x^UKAnK6hpDdrTYpR5x5Syttu={t5lY z&~b@Al;DN(Wbd^43kgQmYS(==sE@=A_)GP-(j_xp9slB#iE%o){?HUJ*t9BS=iwEu*rmFTlL1D*ga zR^23nkXQ@oe1-(NViMt{LX)E^1P1S*Gu8KvFG3ij>yeB# zNz|G`@>q`Cfdo`|bnd}yqspW547b%un5BL=jsC2Mi;lELrlnS>*RaO7PFBb5C<9M@1w{6v(f#oe#b80Ap}wnF_I$*NSmP#$273%DGh ztCS@K40+)%D!6qaqS39%6nwE?Qe5n?YZ!sllq6m&WL1Va;Zz6fWOOJk3KL}ngT(Pu1DxKV8CAA*OU7# zxeU%?LE~7%SVO4^Rn?9K)f*4NFiRn_VKu+2z=HHkdH@M4QzS|EyQ=y=ch&r)7NRL2 zBk4P>pGIm_URrX5mOSR5?v68p#QGw_U=i^&)cBx!nN?i5L{OvM&AMVObk%)KQdbft zuS6YZJftQ)a32m2%SFjUC}Gl*A^ap`?m*u|7sIzQ)5598Dc2JbcY2EtBtwFsT2}dXHWH?HqqHrC&hmDP|DP3 zy|2B(KY5d~B+bqZl07BdrH14jw$SaZ*QWY+#b4Hm>bYx(-fbqbTdWt6GWDbQY+a24 zwt=Q0fdiaiYlvF{K+t^;J^GC?;@ULQ{zVa>G&l<~ZFiSQ==P-SDPw)&S0r=IHY;hI z?_po}vFr4b@LNEF((xf`$u{Zxh?g4drE}=chkHD=k4;L}5ZEWITNO%>B1~ zA}jYVu5ABKV04_DraI_M-FI&ZhQ`k27vPFuV`%gK0^UtzJ!buR70SG!e?tg2EuNL+}uj6_du57h(+bCZ{oq5g~Dpjg((Z-M4oWUf_w zii@z3{ETb$mM;VF@A3+Pj>BTK)6cEdGmScc?m`~Co!G{kzxl?oSj zRSVNQf5q>14F(Ru4a_w!*k}YHx~dZF4$a{{%~EuvPi+V1##K59J*2A<0G2X+WFHL- zHlMG`Q;V}uy4#v)UETywl8AvMJYCKTL}&^<%F;BR>@2`|i{a%a%P~KQ)m*kKg)yTs z(w*vB67N@iDOMeT>$AZHuUW9YmvG@yW!dk$S3ErSaMDP7DH=yrNsI zHXy3>_#7_hgv_3;CdqW$0uaq`@F z<4^$9hu-yBtmFp1(q|Hr=nuL|lML_GkLs1Th`9?2bQM~%`e<{lQl9(5uWP-#qzH0V zH;LYVSySmLnv6F;LT-kIaqg zT5cUGQYcG2JYDh}hbw5j0n-CE2TmBmw#3z-D!m=;IG#8Mn4T^vl9kewOkw62nph%_ z65kzAHDkGOnra05g+i#|84Htx{cx4ltkFV-W=~iURqf>@Z;)Q4pGU{nx&;BlS<(&M zvm}9AKh$ikXvN9C2ehb;d`NT|kF7d`n$THDjMgaKam%AJ&KVe*9Aj)YX*|FO$p(R8 z(7Kybl8A=hu?JB11h_lr`y9k-qpHZ2U#}`A%dv86)GgyJRubE%b-yt(OukccP4_Xu zgUuT?QY~kQ`cTPzLi~~(x3Fk+fJyqK`Z!Ti!MaZpu%P4YIPU_OLTEUxPWQ?sIoc7) zYZcE6L<3_mP(mf{X6#T1B#*CcHr@rLb?8-gXBT4X5)oH=dLAqb`skNq0dY3(EL)y? z!+H7AWGh|$vmw^<+qp;MIp@4;@j-_Ub-g4&l2u9rb}S;C?PN?G^9*OZ8l%VT!bvX1 z(_`DTv~pS#Yq?R^tqa#%6UHRruD@?QuV{oUU&vLLLo+tICr&YawW4m1Mv<#AF4_!s zHV$+n8&&EbH;`N#uS#f*J`QZhSZ<8cyt&JaZ4}e=>qS6gtp14u!X;xP(exObac4pN zE2Cg~V+A3(3QSLSd$cw;W=>g_3)B;p4IRl=?=y9NeDM{mk3AIH>LwCRS0jdM;NUEY zssTDZ72GB}(G!8r;yN{xo&181u3zVF-UC9Cg9=&&>E))EzR;W{y{RKib6kl-6$mbx z*vJ{i_|gc4p!R|3X5pa;e8WrIK^UNA7%srfw{aVYxDt`_rG!>K@sM-G&U6_?EGwTn zLy`2iW6w3AP17gpw^*aidW%1)p)CnbNj&13xONGlhHzJ1qkKd%I9nS8=p%G_pqqwo zd{I)9(yVSpPk-Q;-Q}uI{XJH%+bb2`I7tg{NO2%}Am#a6NHj!mQ!!SeQv7XGT`6Lu*66QGtQkQZLio7(l6-UxLxjkVi&c{W@iAUW z@dzqV34HoX2bJs zpOVCnnqrhO5|lnZR`)HYEh!r0Nao@6nN8XVH~9*>(>} z8nAW6UZj%TVjUrYO zE0jV-Oh_L6d~rzPcBd&C(1kK@<4D=Yi|P>Mlt5JRB1-pAT5Q~-Yaga31)|1tgfJh6 z^a{qo)Kq5Ux(CCRZw7@BS=akIeU zQMez8+d-uroE{ZWUCu;32}F$iAj=sS3`Nx8GKAX8S~bTi{*n?B-^pfF`!fy}xG|%N z)0e!)(WNPJE=Pd?7+Fo)6m5a(g<`uy30u{hw6;ujZ%HrJpHQTDB3utMpTwgA6aq?R z>oKKT@uM&+!VrlF0WxUl@2Kmztt8#grHW{;N)YzE)poOJ#%^cy_?;GhKYfqTHkjib z@xmK#tyMqv$_)ASlGlzN>aJ5;4L@~`-m&sE<=94%Bc3#tU(Rh}xfoB53CfX;L}z38 zn4&zZ39%gSqwP8z?|Pm9()JuyA~Oh3ZCy7xoJ}_S9%2Z>>4r)#DY%i65fh3bIj$#D z!$}z%m$-%ZPaKVwNvO!GVzu;el~k5dE*;i29LBI57s4^41JrW&_}KHzwJeXOShI0! zW#V?Ss!$_xUS9^OgeW28jT{4li+@~-ztg#Y0 zv_tO)U11CKS~vP@Li)n*O@$!ChVm@^^hB`R5w~Pt$#_9YVu`Cr8cP_885mhlC?X}b z*G!+Qff{4I#N6#o6MYFKCJ#mP&iEsK`XyAAoqI-y`iN-U*BZB0WO~FN!pRCL2Njw- z7;zYvD@`&6f)%0LI9*oGHqds00#vfsiP)KfGzF1gq=*ZtVrj_~ZIhK6un zXU>cDBdhxV zBG8Xkzf^w%$5I`+l5^OLvFO0QUYO}=YK;6vzI^&(fu*GR&UulWLOqNXSkJ8@9xU*s@pey@ z(Bmhn^e9|WiMb(5zN~?+{){HCHOc6b8oC2*ILdCT-2|Fd!&Z*`!1Yic^$7CQ9WY(? zgo?ctEm6u%ZoGB4XLRF7$v&`#J4v5WMFqexJXc%6)zmsz3mxvM){K3Od%+g^PP`m(a?Xx2^=ZCtQSl!q?6XeDqXc?yG)Q z_y+Jj5^-?!61R&?aVWKg7!*I6SbmrNGy-%tFIi(J62b*A+;$gdtuISjkw`#|>OQf& z8VB_kBu=byMc?B{OO2K0v&gEYlHuoMlYe8ZF_dd(IOdAVqoAM$<^uUtH5gF%O~kP@ z-p1YWu|R8v`5lUM6wP~Ak5Xt%4_z>qJt`HX$g?E)TWgJT2w?6A|O(vkSafdMG^>vUjY)5 zkZ`$NvO97&ESuB%tb3l#eE+!jKKH$R?(Y5M?!xzu_qB&op3(?ZqxOS|Upl|x*uz}F z-Ywkq!^IM6Wnrc3IhAWkKx_+tW+VFo%O>DTI8a5~IU#gRmVMb+)W=AVcR?#UII=8e zt)+Uj+}RrHYUwH206ko(T@+X?_^rDy?XMHeZ)rlFX5`J>lI=D+TjN(&+#h!Ut@Q)3tP&j zM)(7l(m?)hyKnZ0P+Mz~`2!!wxYq@bR_JSZ6SYrPaIJAgLpkA8y!(NZDj}WAeCI~g zvDr7vc}9iw`qv<+jq%m@I1l+xKlI8o3xkg0-UuBIa6z>7!V>vv9^c$JXlSl#PP`7V z^R>aEmt93pDlhq1aZOTh{%219uQ^XTBp(VJkiIA1pg1m-(y z!+Ep?ETO4-n7RHdGV?tAVr^4?n?kbhh}QfEu~xs%(RW3)IbWmvh!Qjbt_F55r0LtC zb#h%LceE~;+i;l3oho3ZR*;-Y4Cdunx=2n=G0Bv`ZpGq!YnW6YZ8+-AD}Z0t&b0Q2 zDdkO)uH*Ku8;Z!MMB<>%Btsq%G`(8V=NIJwZcyq!5~PK>t@W)g-D|>;Kz*b;FbWK0 zeut;?auz8!OLY-J5`s`ycAgk>oq*362?(bK1mDnCa{IV1>vl|Ss|Gp>6r#vv7W8o? z6iuX9eEes+XujdVz7kSfbD-<2@r;W8;a@_tsgeNavi_&Zjgk`cVWl&TNa$P%<`Iat z=BAc)uG1{%)t-?2v1>eV#we#$i7Rb&?T*%UPWT><&jR#|<0+9jtdbT!-Pmkif73b; zEGqB8My9$V_EQg}v?_q6;|Q#Og+J%}U%L@Av|8hS0w&=vV~l50Yq^D=ejy&t#TQrhd|^>+?`$ z*OT$@hB(FRgmvo5w2@p#40zCP)|>~T$&}4XTLbX*YT|P=c~jZYi*&Q#0J%IAELs#S za%_2!FwAF7m#QvM&Js+lC56hFo5nozCsXg{Ypqp17gpLw4=^mF@+AAOGYT`a#}K$xSSzPk7)OSpA}zf zuH?M9*m|mo1Wr9qwGs*8#{Lu9(*snEjsJW}7Clw7xT~OX9a(s0Qj+B8;3;LfrN*B* z-VAS>Hm0@#v@%WG0!Mj^4>a_VMYml|i#SU@%N|%ov&x+}((i4dn5!2IKHZoLWm#esZa%WUe04S&152q9u@YGKXYV*N~Y?%~C!oE^6#ks77!l z{!#$%X%@AiS?tpI&h#p^)x^#d1(9enH6Ux&|D}?HAwKcg>&W%YkK3ZFXW&aZBO@&33<)*C{L)_m8r$0V2 zo)ai2rbDho7U`D(LY}sUy@l9Eo^(;aStTqfG^7rSAqR=C8x!O$F`O+Ee;#GvUSb2o zQu!$@uv5Eo0Ue|c_%d@h_fDz!=@Z%D=?m9`O&f(h)qD~?D3URl;~qa54xi|cO^&@m zHRl)G!sL5Nlxa7Gw>oHo(WJ=ai>sl@0=La)dWb@a#G42)V`6FQK4o1Ct)Z@_t1>-% zP|D4X1{Wwt3{*ez&oheQADJ&`s0P06qPRKOZ4H4wvEwo##Zx$PV=s_E z*IZKrV~@o)8Wl6$6WNb2-2R8ejsw{pOk=CnP8{46KZ9dO)cd43C!uay_N6l29=hA~ z$87FWlBwOUFLe%yxsf>|jW5{rK1xbMjs2~(H**q!n_$Sz)FAWm4@s2PDg#BZY_WCkaQD}K}rTvO@1R~n} z$P`Wq6Z|qf#>JyeXFif{RccwTr$ra%+ttl1{rCU1j>(@(QCZyEt&Q2$B?)i0hG+TX zXg7I3I`cmY^JOkDyE6>uWO5&zmD^2ENuIT3zLVT8A`^kh?IO~%9-CVn?)+`!H86Qv z4DVfDg8aUx6wAL?ksrfu5j#FSyd^{y1dAVn=jV94cbQ^DAY*=>yOk+M)ZDQBL8N`> z^8BzmA~J7S9z+`WPIV-~>mUuod#@t{t>cU*Av*mr>xd{q-daaQxp0oRd#@v+Ksi6p z-PZ9|{*iNXF3-K!5wR}BWWys+*}wVM0c4+%hT#vFyP0{?lVO(!(CL4=euy0+_WeYx z6_J`@F>=kHF{cta#`4`OQfo&h(-p=81W*5SfJLN<6 z5j!5lTJjgL^eECWTiN$eai977yOaI)-URs0ai{#Ze%}#`%K7#L&(3ZQhqM0~ok78K z95~$egnbFfBF?UaRD{=^$$R#&G7|QrBLUnnK9_sK`*!NJjrEk?Ya8p~0l7{f&1)O$ z;lVrgGT>#v%Yc^wF9Ti%ybO35@G{_Kz{|k@7X!&_W{+o$2t_f9QV$+#@E{IrkD~E% z6n}?e1d22ipGWaU6kkFy8pRk7?%#mMacKN18gE47*U|V5G~R^9o6$HPjhSf7LUF4H z*Ce3*+dOy)#r=t%c+GYX;&!0@DDK~hF2~`jp+7x!LSAsuedLBJaBDpw?g3c&8W}{w z9u#4Ii<#(v)K<6;Y4ZoJ9-m!`!<|A0$j}g%z}+!B_?E~{*|#HQ-#-T8uz{)wY{k1o z&x0{!nBLQNq1(gq%^#?{%e431YPa)sobO9|-x1?IU@rq+2D}X1rwn-f_I~ Date: Thu, 20 Nov 2025 12:43:54 -0500 Subject: [PATCH 06/11] Fix errors --- cosipy/threeml/COSILike.py | 5 ++--- .../polarization/maximum_likelihood_method.ipynb | 13 +++++-------- 2 files changed, 7 insertions(+), 11 deletions(-) diff --git a/cosipy/threeml/COSILike.py b/cosipy/threeml/COSILike.py index 9f55e5a7..b1969f9e 100644 --- a/cosipy/threeml/COSILike.py +++ b/cosipy/threeml/COSILike.py @@ -262,10 +262,9 @@ def set_model(self, model): this_expectation = self._psr[name].get_expectation(spectrum) else: if self._coordsys == 'spacecraftframe': - this_expectation = self._psr[name].get_expectation(spectrum, source.components['grb'].polarization) + this_expectation = self._psr[name].get_expectation(spectrum, source.components[item].polarization) elif self._coordsys == 'galactic': - scatt_map = self._get_scatt_map(source.position.sky_coord) - this_expectation = self._psr[name].get_expectation(spectrum, source.components['grb'].polarization) + this_expectation = self._psr[name].get_expectation(spectrum, source.components[item].polarization) else: raise RuntimeError("Unknown coordinate system") diff --git a/docs/tutorials/polarization/maximum_likelihood_method.ipynb b/docs/tutorials/polarization/maximum_likelihood_method.ipynb index 641264cc..7a5f5539 100644 --- a/docs/tutorials/polarization/maximum_likelihood_method.ipynb +++ b/docs/tutorials/polarization/maximum_likelihood_method.ipynb @@ -20,18 +20,15 @@ "metadata": {}, "outputs": [], "source": [ - "from cosipy import UnBinnedData, COSILike\n", + "from cosipy import BinnedData, COSILike\n", "from cosipy.spacecraftfile import SpacecraftFile\n", "from cosipy.threeml.custom_functions import Band_Eflux\n", "from astropy.time import Time\n", - "import numpy as np\n", "from astropy.coordinates import SkyCoord\n", "from astropy import units as u\n", - "from scoords import SpacecraftFrame\n", "from cosipy.util import fetch_wasabi_file\n", "from pathlib import Path\n", - "from threeML import LinearPolarization, SpectralComponent, PointSource, Model, JointLikelihood, DataList\n", - "from astromodels import Parameter" + "from threeML import LinearPolarization, SpectralComponent, PointSource, Model, JointLikelihood, DataList" ] }, { @@ -111,17 +108,17 @@ "source": [ "data_path = Path(\"\") # Update to your path\n", "\n", - "grb_background = UnBinnedData(data_path/'grb.yaml')\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='grb_background_source_interval.fits.gz', output_name='grb_background_binned_galactic', psichi_binning='galactic')\n", "grb_background.load_binned_data_from_hdf5('grb_background_binned_galactic.hdf5')\n", "\n", - "background_before = UnBinnedData(data_path/'background_before.yaml')\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('background_before_binned_galactic.hdf5')\n", "\n", - "background_after = UnBinnedData(data_path/'background_after.yaml') # e.g. background_after.yaml\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='background_after.fits.gz', output_name='background_after_binned_galactic', psichi_binning='galactic')\n", "background_after.load_binned_data_from_hdf5('background_after_binned_galactic.hdf5')" From 7a75c81e73e34e3d30ffd07e9dcc71bb7804fe52 Mon Sep 17 00:00:00 2001 From: Eliza Neights Date: Thu, 20 Nov 2025 12:51:49 -0500 Subject: [PATCH 07/11] Fix error in unit test --- tests/polarization/test_polarization_mlm.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/polarization/test_polarization_mlm.py b/tests/polarization/test_polarization_mlm.py index fdfda481..06356128 100644 --- a/tests/polarization/test_polarization_mlm.py +++ b/tests/polarization/test_polarization_mlm.py @@ -13,7 +13,6 @@ from cosipy import test_data analysis = BinnedData(test_data.path / 'polarization_data_mlm.yaml') -analysis.get_binned_data(unbinned_data='polarization_data.hdf5', output_name='polarization_data_binned', psichi_binning='galactic') analysis.load_binned_data_from_hdf5(test_data.path / 'polarization_data_binned.hdf5') response_path = test_data.path / 'test_polarization_response.h5' sc_orientation = SpacecraftFile.parse_from_file(test_data.path / 'polarization_ori.ori') From 7206b82568b2eaccf34ac6f13620b050c9ba7274 Mon Sep 17 00:00:00 2001 From: Eliza Neights Date: Thu, 20 Nov 2025 12:57:46 -0500 Subject: [PATCH 08/11] Fix error in COSILike that breaks it when there is no polarization --- cosipy/threeml/COSILike.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/cosipy/threeml/COSILike.py b/cosipy/threeml/COSILike.py index a0c32429..667e75d0 100644 --- a/cosipy/threeml/COSILike.py +++ b/cosipy/threeml/COSILike.py @@ -137,6 +137,8 @@ def __init__(self, name, dr, data, bkg, sc_orientation, nuisance_param=None, coo self._pa_convention = IAUPolarizationConvention() else: raise RuntimeError("Unknown coordinate system") + else: + self._response_pa_convention = None def set_model(self, model): """ From ed57498a8a7fe3fe460ef96ed46a848d00da101c Mon Sep 17 00:00:00 2001 From: Eliza Neights Date: Thu, 20 Nov 2025 13:17:15 -0500 Subject: [PATCH 09/11] Fix unit tests --- tests/polarization/test_polarization_asad.py | 16 ++++++++-------- tests/polarization/test_polarization_mlm.py | 4 ++-- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/tests/polarization/test_polarization_asad.py b/tests/polarization/test_polarization_asad.py index 22474fc5..6dc90b6f 100644 --- a/tests/polarization/test_polarization_asad.py +++ b/tests/polarization/test_polarization_asad.py @@ -66,9 +66,9 @@ def test_spacecraft_fit(): polarization_fit_spacecraft['fraction uncertainty'], polarization_fit_spacecraft['angle'].angle.rad, polarization_fit_spacecraft['angle uncertainty'].rad], - [0.8114804627334942, 0.8081587949263002, - 1.5713378840593466, 0.5340212799099183], - atol=[1.0, 0.5, 1.0, 0.1]) + [0.8248301194006322, 0.7889544546660269, + 1.5645814961205469, 0.5143016570410734], + atol=[0.2, 0.1, 0.2, 0.1]) # ASAD from binned data polarization_spacecraft = PolarizationASAD(source_direction, @@ -86,7 +86,7 @@ def test_spacecraft_fit(): polarization_fit_spacecraft['angle uncertainty'].rad], [0.9452187271167997, 0.9328483275998886, 1.993361180746714, 0.6416512077658346], - atol=[1.0, 0.5, 1.0, 0.1]) + atol=[0.2, 0.1, 0.2, 0.1]) def test_icrs_fit(): @@ -103,9 +103,9 @@ def test_icrs_fit(): polarization_fit_icrs['fraction uncertainty'], polarization_fit_icrs['angle'].angle.rad, polarization_fit_icrs['angle uncertainty'].rad], - [1.6268965437885632, 0.9763744515967512, - 1.8111679143685155, 0.40112053082203614], - atol=[1.0, 0.5, 1.0, 0.1]) + [1.496745801986812, 0.8751758133432178, + 1.84399578798795, 0.3812557749920544], + atol=[0.2, 0.1, 0.2, 0.1]) # ASAD from binned data polarization_icrs = PolarizationASAD(source_direction.transform_to('galactic'), @@ -122,4 +122,4 @@ def test_icrs_fit(): polarization_fit_icrs['angle uncertainty'].rad], [2.02118419504387, 0.7661035298627569, 1.6238519333293382, 0.22647693546905168], - atol=[1.0, 0.5, 1.0, 0.1]) + atol=[0.2, 0.1, 0.2, 0.1]) diff --git a/tests/polarization/test_polarization_mlm.py b/tests/polarization/test_polarization_mlm.py index 06356128..a3be0a54 100644 --- a/tests/polarization/test_polarization_mlm.py +++ b/tests/polarization/test_polarization_mlm.py @@ -37,7 +37,7 @@ source_direction = SkyCoord(0, 70, representation_type='spherical', frame=SpacecraftFrame(attitude=attitude), unit=u.deg).transform_to('galactic') -polarization = LinearPolarization(1, 100) +polarization = LinearPolarization(0.5, 100) spectral_component = SpectralComponent('test', spectrum, polarization) source = PointSource('test', @@ -68,4 +68,4 @@ def test_polarization_fit(): like.fit() assert np.allclose([source.spectrum.test.polarization.degree.value, source.spectrum.test.polarization.angle.value], - [100., 112.24], atol=[1.0, 1.0]) + [0.000012, 100.], atol=[0.1, 1.0]) From aff6d74fffc7ca6ee40a32ded7b1d8bd795174af Mon Sep 17 00:00:00 2001 From: Eliza Neights Date: Thu, 20 Nov 2025 13:30:40 -0500 Subject: [PATCH 10/11] Update point source response unit test --- tests/response/test_point_source_response.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/response/test_point_source_response.py b/tests/response/test_point_source_response.py index 799fa01b..71b43bdb 100644 --- a/tests/response/test_point_source_response.py +++ b/tests/response/test_point_source_response.py @@ -145,4 +145,4 @@ def test_get_expectation(): ## get expectation with polarization exp = psr_pol.get_expectation(const, polarization=LinearPolarization(angle=180, degree=100)) assert isinstance(exp, Histogram) - assert np.isclose(np.sum(exp.contents), 6.30823539e+11, rtol=1e-8) \ No newline at end of file + assert np.isclose(np.sum(exp.contents), 7.56988247e+12, rtol=1e-8) \ No newline at end of file From a059bb384d5154fa19461e5544554227e5fdd65d Mon Sep 17 00:00:00 2001 From: Eliza Neights Date: Thu, 20 Nov 2025 15:00:49 -0500 Subject: [PATCH 11/11] Update example notebook to include output --- .../maximum_likelihood_method.ipynb | 310 +++++++++++++++++- 1 file changed, 297 insertions(+), 13 deletions(-) diff --git a/docs/tutorials/polarization/maximum_likelihood_method.ipynb b/docs/tutorials/polarization/maximum_likelihood_method.ipynb index 7a5f5539..c61e2976 100644 --- a/docs/tutorials/polarization/maximum_likelihood_method.ipynb +++ b/docs/tutorials/polarization/maximum_likelihood_method.ipynb @@ -16,7 +16,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -110,18 +110,18 @@ "\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='grb_background_source_interval.fits.gz', output_name='grb_background_binned_galactic', psichi_binning='galactic')\n", - "grb_background.load_binned_data_from_hdf5('grb_background_binned_galactic.hdf5')\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('background_before_binned_galactic.hdf5')\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='background_after.fits.gz', output_name='background_after_binned_galactic', psichi_binning='galactic')\n", - "background_after.load_binned_data_from_hdf5('background_after_binned_galactic.hdf5')" + "#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')" ] }, { @@ -137,7 +137,7 @@ "metadata": {}, "outputs": [], "source": [ - "response_file = data_path/'ResponseContinuum.o3.pol.e200_10000.b4.p12.relx.s10396905069491.m420.filtered.binnedpolarization.11D.h5' # e.g. ResponseContinuum.o3.pol.e200_10000.b4.p12.s10396905069491.m441.filtered.binnedpolarization.11D.h5\n", + "response_file = data_path / 'ResponseContinuum.o3.pol.e200_10000.b4.p12.relx.s10396905069491.m420.filtered.binnedpolarization.11D.h5' # e.g. ResponseContinuum.o3.pol.e200_10000.b4.p12.s10396905069491.m441.filtered.binnedpolarization.11D.h5\n", "\n", "sc_orientation = SpacecraftFile.parse_from_file(data_path/'DC3_final_530km_3_month_with_slew_1sbins_GalacticEarth_SAA.ori') # e.g. DC3_final_530km_3_month_with_slew_1sbins_GalacticEarth_SAA.ori\n", "sc_orientation = sc_orientation.source_interval(Time(1835493492.2, format = 'unix'), Time(1835493492.8, format = 'unix'))" @@ -152,7 +152,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 9, "metadata": {}, "outputs": [], "source": [ @@ -188,12 +188,12 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 12, "id": "e942447f", "metadata": {}, "outputs": [], "source": [ - "polarization = LinearPolarization(75, 50) # polarization level (percentage out of 100), polarization angle (degrees)\n", + "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", @@ -226,10 +226,294 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 13, "id": "7f853caa", "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/html": [ + "
14:58:58 INFO      set the minimizer to minuit                                             joint_likelihood.py:1045\n",
+       "
\n" + ], + "text/plain": [ + "\u001b[38;5;46m14:58:58\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=596975;file:///Users/eneights/software/miniforge3/envs/cosipy-fork/lib/python3.10/site-packages/threeML/classicMLE/joint_likelihood.py\u001b\\\u001b[2mjoint_likelihood.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=647829;file:///Users/eneights/software/miniforge3/envs/cosipy-fork/lib/python3.10/site-packages/threeML/classicMLE/joint_likelihood.py#1045\u001b\\\u001b[2m1045\u001b[0m\u001b]8;;\u001b\\\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Adding 1e-12 to each bin of the expectation to avoid log-likelihood = -inf.\n" + ] + }, + { + "data": { + "text/html": [ + "
14:59:06 WARNING   get_number_of_data_points not implemented, values for statistical        plugin_prototype.py:130\n",
+       "                  measurements such as AIC or BIC are unreliable                                                   \n",
+       "
\n" + ], + "text/plain": [ + "\u001b[38;5;46m14:59:06\u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m get_number_of_data_points not implemented, values for statistical \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=818515;file:///Users/eneights/software/miniforge3/envs/cosipy-fork/lib/python3.10/site-packages/threeML/plugin_prototype.py\u001b\\\u001b[2mplugin_prototype.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=821737;file:///Users/eneights/software/miniforge3/envs/cosipy-fork/lib/python3.10/site-packages/threeML/plugin_prototype.py#130\u001b\\\u001b[2m130\u001b[0m\u001b]8;;\u001b\\\n", + "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mmeasurements such as AIC or BIC are unreliable \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\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.80 +/- 0.27) x 10
source.spectrum.grb.polarization.angle(9.00016 +/- 0.00006) x 10deg
\n", + "
" + ], + "text/plain": [ + " result unit\n", + "parameter \n", + "source.spectrum.grb.polarization.degree (2.80 +/- 0.27) x 10 \n", + "source.spectrum.grb.polarization.angle (9.00016 +/- 0.00006) 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.000.00
0.001.00
" + ], + "text/plain": [ + "1.00 0.00\n", + "0.00 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)
cosi19260.057111
total19260.057111
\n", + "
" + ], + "text/plain": [ + " -log(likelihood)\n", + "cosi 19260.057111\n", + "total 19260.057111" + ] + }, + "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
AIC38518.114221
BIC38520.114221
\n", + "
" + ], + "text/plain": [ + " statistical measures\n", + "AIC 38518.114221\n", + "BIC 38520.114221" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "( value negative_error \\\n", + " source.spectrum.grb.polarization.degree 27.984015 -2.817586 \n", + " source.spectrum.grb.polarization.angle 90.001559 -0.000649 \n", + " \n", + " positive_error error unit \n", + " source.spectrum.grb.polarization.degree 2.726043 2.771814 \n", + " source.spectrum.grb.polarization.angle 0.000627 0.000638 deg ,\n", + " -log(likelihood)\n", + " cosi 19260.057111\n", + " total 19260.057111)" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "cosi = COSILike(\"cosi\", # COSI 3ML plugin\n", " dr = response_file, # detector response\n",