diff --git a/cosipy/response/PointSourceResponse.py b/cosipy/response/PointSourceResponse.py index d139a5fb..6da6a3d9 100644 --- a/cosipy/response/PointSourceResponse.py +++ b/cosipy/response/PointSourceResponse.py @@ -88,10 +88,11 @@ def get_expectation(self, spectrum, polarization=None, flux=None): # unpolarized weights weights = np.full(pol_axis.nbins, (1. - polarization_level) / pol_axis.nbins) - + # add polarized weights polarization_bin_index = pol_axis.find_bin(polarization_angle * u.deg) weights[polarization_bin_index] += polarization_level + weights *= self.axes['Pol'].nbins contents = np.tensordot(weights, self.contents, axes=(0, self.axes.label_to_index('Pol'))) diff --git a/cosipy/test_data/polarization_data_binned.hdf5 b/cosipy/test_data/polarization_data_binned.hdf5 new file mode 100644 index 00000000..632dc42b Binary files /dev/null and b/cosipy/test_data/polarization_data_binned.hdf5 differ diff --git a/cosipy/test_data/polarization_data_mlm.yaml b/cosipy/test_data/polarization_data_mlm.yaml new file mode 100644 index 00000000..7feeccc3 --- /dev/null +++ b/cosipy/test_data/polarization_data_mlm.yaml @@ -0,0 +1,14 @@ +#----------# +# Data I/O: + +data_file: "abc" # full path +ori_file: 'data3/ori.ori' +unbinned_output: 'hdf5' # 'fits' or 'hdf5' +time_bins: 1 # time bin size in seconds. Takes int or list of bin edges. +energy_bins: [200., 500., 1000., 2000., 10000.] # Takes list. Needs to match response. +phi_pix_size: 6 # binning of Compton scattering angle [deg] +nside: 1 # healpix binning of psi chi local +scheme: 'ring' # healpix binning of psi chi local +tmin: 0. # Min time cut in seconds. +tmax: 13. # Max time cut in seconds. +#----------# diff --git a/cosipy/threeml/COSILike.py b/cosipy/threeml/COSILike.py index f7cf9f09..667e75d0 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,24 @@ 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") + else: + self._response_pa_convention = None def set_model(self, model): """ @@ -215,7 +237,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") @@ -225,11 +247,35 @@ 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 - spectrum = source.spectrum.main.shape + 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(): - total_expectation = self._psr[name].get_expectation(spectrum) + spectrum = getattr(source.spectrum, item).shape + + if not 'Pol' in self._dr.axes.labels: + this_expectation = self._psr[name].get_expectation(spectrum) + else: + if self._coordsys == 'spacecraftframe': + this_expectation = self._psr[name].get_expectation(spectrum, source.components[item].polarization) + elif self._coordsys == 'galactic': + this_expectation = self._psr[name].get_expectation(spectrum, source.components[item].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() @@ -321,7 +367,8 @@ def _get_dwell_time_map(self, coord): src_path = self._sc_orientation.get_target_in_sc_frame(coord) dwell_time_map = self._sc_orientation.get_dwell_map(response = self._rsp_path, - src_path = src_path) + src_path = src_path, + pa_convention = self._response_pa_convention) return dwell_time_map diff --git a/docs/tutorials/polarization/maximum_likelihood_method.ipynb b/docs/tutorials/polarization/maximum_likelihood_method.ipynb new file mode 100644 index 00000000..c61e2976 --- /dev/null +++ b/docs/tutorials/polarization/maximum_likelihood_method.ipynb @@ -0,0 +1,554 @@ +{ + "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": 2, + "metadata": {}, + "outputs": [], + "source": [ + "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", + "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", + "from threeML import LinearPolarization, SpectralComponent, PointSource, Model, JointLikelihood, DataList" + ] + }, + { + "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 = 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", + "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": 9, + "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": 12, + "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", + "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": 13, + "id": "7f853caa", + "metadata": {}, + "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", + " 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 +} 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 new file mode 100644 index 00000000..a3be0a54 --- /dev/null +++ b/tests/polarization/test_polarization_mlm.py @@ -0,0 +1,71 @@ +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 +from threeML import LinearPolarization, SpectralComponent, PointSource, Model, JointLikelihood, DataList +from scoords import SpacecraftFrame +import numpy as np + +from cosipy import BinnedData, COSILike +from cosipy.spacecraftfile import SpacecraftFile +from cosipy.threeml.custom_functions import Band_Eflux +from cosipy import test_data + +analysis = BinnedData(test_data.path / 'polarization_data_mlm.yaml') +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') +attitude = sc_orientation.get_attitude()[0] + +a = 10. * u.keV +b = 10000. * u.keV +alpha = -1. +beta = -2. +ebreak = 350. * u.keV +K = 50. / 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 + +source_direction = SkyCoord(0, 70, representation_type='spherical', frame=SpacecraftFrame(attitude=attitude), unit=u.deg).transform_to('galactic') + +polarization = LinearPolarization(0.5, 100) +spectral_component = SpectralComponent('test', spectrum, polarization) + +source = PointSource('test', + l = source_direction.l.deg, + b = source_direction.b.deg, + components = [spectral_component]) + +source.components['test'].shape.K.fix = True +source.components['test'].shape.E0.fix = True +source.components['test'].shape.alpha.fix = True +source.components['test'].shape.beta.fix = True + +model = Model(source) + +def test_polarization_fit(): + + cosi = COSILike("cosi", + dr = response_path, + data = analysis.binned_data.project('Em', 'Phi', 'PsiChi'), + bkg = analysis.binned_data.project('Em', 'Phi', 'PsiChi') * 0., + sc_orientation = sc_orientation, + response_pa_convention = 'RelativeZ') + + plugins = DataList(cosi) + + like = JointLikelihood(model, plugins, verbose=False) + + like.fit() + + assert np.allclose([source.spectrum.test.polarization.degree.value, source.spectrum.test.polarization.angle.value], + [0.000012, 100.], atol=[0.1, 1.0]) 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