Maximum likelihood polarization fitting

cosipy/polarization/__init__.py

@@ -1,3 +1,3 @@
-from .polarization_asad import PolarizationASAD, calculate_uncertainties
+from .polarization_asad import PolarizationASAD
 from .conventions import PolarizationConvention, OrthographicConvention, StereographicConvention, IAUPolarizationConvention
 from .polarization_angle import PolarizationAngle

cosipy/polarization/polarization_angle.py

@@ -1,6 +1,7 @@
 import numpy as np
 from astropy.coordinates import SkyCoord, Angle
 import astropy.units as u
+from scoords import SpacecraftFrame
 
 from .conventions import PolarizationConvention
 
@@ -100,4 +101,70 @@ def transform_to(self, convention, *args, **kwargs):
                                  self.skycoord,
                                  convention = convention2)
 
+    @classmethod
+    def from_scattering_direction(cls, psichi, source_coord, convention):
+        """
+        Calculate the azimuthal scattering angle of a scattered photon.
+
+        Parameters
+        ----------
+        psichi : astropy.coordinates.SkyCoord
+            Scattered photon direction
+        source_coord : astropy.coordinates.SkyCoord
+            Source direction
+        convention :
+            cosipy.polarization.PolarizationConvention
+
+        Returns
+        -------
+        azimuthal_scattering_angle : cosipy.polarization.PolarizationAngle
+            Azimuthal scattering angle
+        """
+
+        if convention.frame.name == 'spacecraftframe':
+
+            source_coord = source_coord.transform_to(SpacecraftFrame(attitude=source_coord.attitude))
+            psichi = psichi.transform_to(SpacecraftFrame(attitude=source_coord.attitude))
+
+        elif convention.frame.name == 'icrs':
+
+            source_coord = source_coord.transform_to('icrs')
+            psichi = psichi.transform_to('icrs')
+
+        elif convention.frame.name == 'galactic':
+
+            source_coord = source_coord.transform_to('galactic')
+            psichi = psichi.transform_to('galactic')
+
+        else:
+
+            raise RuntimeError('Unsupported polarization convention frame "' + convention.frame.name + '"')
+
+        reference_coord = convention.get_basis(source_coord)[0]
+
+        source_vector_cartesian = source_coord.cartesian.xyz.value
+        reference_vector_cartesian = reference_coord.cartesian.xyz.value
+        scattered_photon_vector = psichi.cartesian.xyz.value
+
+        # Project scattered photon vector onto plane perpendicular to source direction
+        d = np.dot(scattered_photon_vector, source_vector_cartesian) / np.dot(source_vector_cartesian, source_vector_cartesian)
+        projection = [scattered_photon_vector[0] - (d * source_vector_cartesian[0]), 
+                      scattered_photon_vector[1] - (d * source_vector_cartesian[1]), 
+                      scattered_photon_vector[2] - (d * source_vector_cartesian[2])]
+
+        # Calculate angle between scattered photon vector & reference vector on plane perpendicular to source direction
+        cross_product = np.cross(projection, reference_vector_cartesian)
+        if np.dot(source_vector_cartesian, cross_product) < 0:
+            sign = -1
+        else:
+            sign = 1
+        normalization = np.sqrt(np.dot(projection, projection)) * np.sqrt(np.dot(reference_vector_cartesian, reference_vector_cartesian))
+
+        angle = Angle(sign * np.arccos(np.dot(projection, reference_vector_cartesian) / normalization), unit=u.rad)
+
+        azimuthal_scattering_angle = cls(angle, source_coord, convention=convention)
+
+        return azimuthal_scattering_angle
+
+
 

cosipy/response/FullDetectorResponse.py

@@ -810,7 +810,8 @@ def get_point_source_response(self,
                                   exposure_map = None,
                                   coord = None,
                                   scatt_map = None,
-                                  Earth_occ = True):
+                                  Earth_occ = True,
+                                  pol_convention = 'RelativeX'):
         """
         Convolve the all-sky detector response with exposure for a source at a given
         sky location.
@@ -830,6 +831,8 @@ def get_point_source_response(self,
             Option to include Earth occultation in the respeonce. 
             Default is True, in which case you can only pass one 
             coord, which must be the same as was used for the scatt map. 
+        pol_convention : str, optional
+            Polarization convention of response ('RelativeX', 'RelativeY', or 'RelativeZ') 
 
         Returns
         -------
@@ -901,7 +904,7 @@ def get_point_source_response(self,
 
                 dr_pix.axes['PsiChi'].coordsys = SpacecraftFrame(attitude = att)
 
-                self._sum_rot_hist(dr_pix, psr, exposure)
+                self._sum_rot_hist(dr_pix, psr, exposure, coord, pol_convention)
 
             # Convert to PSR
             psr = tuple([PointSourceResponse(psr.axes[1:],
@@ -1085,14 +1088,14 @@ def merge_psr_to_extended_source_response(self, basename, coordsys = 'galactic',
         return extended_source_response
 
     @staticmethod
-    def _sum_rot_hist(h, h_new, exposure, axis = "PsiChi"):
+    def _sum_rot_hist(h, h_new, exposure, coord, pol_convention, axis = "PsiChi"):
         """
         Rotate a histogram with HealpixAxis h into the grid of h_new, and sum
         it up with the weight of exposure.
 
         Meant to rotate the PsiChi of a CDS from local to galactic
         """
-        
+
         axis_id = h.axes.label_to_index(axis)
 
         old_axes = h.axes
@@ -1104,11 +1107,35 @@ def _sum_rot_hist(h, h_new, exposure, axis = "PsiChi"):
         # Convolve
         # TODO: Change this to interpolation (pixels + weights)
         old_pixels = old_axis.find_bin(new_axis.pix2skycoord(np.arange(new_axis.nbins)))
+
+        if 'Pol' in h.axes.labels and h_new.axes[axis].coordsys.name == 'galactic' or h_new.axes[axis].coordsys.name == 'icrs':
+
+            from cosipy.polarization.polarization_angle import PolarizationAngle
+            from cosipy.polarization.conventions import IAUPolarizationConvention
+
+            pol_axis_id = h.axes.label_to_index('Pol')
+
+            old_pol_axis = h.axes[pol_axis_id]
+            new_pol_axis = h_new.axes[pol_axis_id]
+
+            old_pol_indices = []
+            for i in range(h_new.axes['Pol'].nbins):
+
+                pa = PolarizationAngle(h_new.axes['Pol'].centers.to_value(u.deg)[i] * u.deg, coord.transform_to('icrs'), convention=IAUPolarizationConvention())
+                pa_old = pa.transform_to(pol_convention, attitude=coord.attitude)
+
+                if pa_old.angle.deg == 180.:
+                    pa_old = PolarizationAngle(0. * u.deg, coord, convention=IAUPolarizationConvention())
+
+                old_pol_indices.append(old_pol_axis.find_bin(pa_old.angle))
+
+            old_pol_indices = np.array(old_pol_indices)
+
         # NOTE: there are some pixels that are duplicated, since the center 2 pixels
         # of the original grid can land within the boundaries of a single pixel
         # of the target grid. The following commented code fixes this, but it's slow, and
         # the effect is very small, so maybe not worth it
-        # nulambda_npix = h.axes['NuLamnda'].nbins    
+        # nulambda_npix = h.axes['NuLamnda'].nbins
         # new_norm = np.zeros(shape = nulambda_npix)
         # for p in old_pixels:
         #     h_slice = h[{axis:p}]
@@ -1120,14 +1147,29 @@ def _sum_rot_hist(h, h_new, exposure, axis = "PsiChi"):
 
             #h_new[{axis:new_pix}] += exposure * h[{axis: old_pix}] # * norm_corr
             # The following code does the same than the code above, but is faster
-            # However, it uses some internal functionality in histpy, which is bad practice
-            # TODO: change this in a future version. We might need to modify histpy so that
-            # this is not needed
-
-            old_indices = tuple([slice(None)]*axis_id + [old_pix+1])
-            new_indices = tuple([slice(None)]*axis_id + [new_pix+1])
 
-            h_new._contents[new_indices] += exposure * h._contents[old_indices] # * norm_corr
+            if not 'Pol' in h.axes.labels:
+
+                old_index = (slice(None),)*axis_id + (old_pix,)
+                new_index = (slice(None),)*axis_id + (new_pix,)
+
+                h_new[new_index] += exposure * u.s * h[old_index] # * norm_corr
+
+            else:
+
+                for old_pol_bin,new_pol_bin in zip(old_pol_indices,range(new_pol_axis.nbins)):
+
+                    if pol_axis_id < axis_id:
+
+                        old_index = (slice(None),)*pol_axis_id + (old_pol_bin,) + (slice(None),)*(axis_id-pol_axis_id-1) + (old_pix,)
+                        new_index = (slice(None),)*pol_axis_id + (new_pol_bin,) + (slice(None),)*(axis_id-pol_axis_id-1) + (new_pix,)
+
+                    else:
+
+                        old_index = (slice(None),)*axis_id + (old_pix,) + (slice(None),)*(pol_axis_id-axis_id-1) + (old_pol_bin,)
+                        new_index = (slice(None),)*axis_id + (new_pix,) + (slice(None),)*(pol_axis_id-axis_id-1) + (new_pol_bin,)
+
+                    h_new[new_index] += exposure * u.s * h[old_index] # * norm_corr
 
 
     def __str__(self):

cosipy/response/PointSourceResponse.py

@@ -1,14 +1,14 @@
-from histpy import Histogram, Axes, Axis
+from histpy import Histogram#, Axes, Axis
 
 import numpy as np
 import astropy.units as u
-from astropy.units import Quantity
-from scipy import integrate
-
-from threeML import DiracDelta, Constant, Line, Quadratic, Cubic, Quartic, StepFunction, StepFunctionUpper, Cosine_Prior, Uniform_prior, PhAbs, Gaussian
+from scoords import SpacecraftFrame, Attitude
 
 from .functions import get_integrated_spectral_model
 
+import logging
+logger = logging.getLogger(__name__)
+
 class PointSourceResponse(Histogram):
     """
     Handles the multi-dimensional matrix that describes the expected
@@ -43,32 +43,70 @@ def photon_energy_axis(self):
 
         return self.axes['Ei']
 
-    def get_expectation(self, spectrum):
+    def get_expectation(self, spectrum, polarization=None):
         """
-        Convolve the response with a spectral hypothesis to obtain the expected
+        Convolve the response with a spectral (and optionally, polarization) hypothesis to obtain the expected
         excess counts from the source.
 
         Parameters
         ----------
         spectrum : :py:class:`threeML.Model`
             Spectral hypothesis.
-
+        polarization : 'astromodels.core.polarization.LinearPolarization', optional
+            Polarization angle and degree. The angle is assumed to have same convention as point source response.
+
         Returns
         -------
         :py:class:`histpy.Histogram`
              Histogram with the expected counts on each analysis bin
         """
-
+
+        if polarization is None:
+
+            if 'Pol' in self.axes.labels:
+
+                raise RuntimeError("Must include polarization in point source response if using polarization response")
+
+            contents = self.contents
+            axes = self.axes[1:]
+
+        else:
+
+            if not 'Pol' in self.axes.labels:
+
+                raise RuntimeError("Response must have polarization angle axis to include polarization in point source response")
+
+            polarization_angle = polarization.angle.value
+            polarization_level = polarization.degree.value / 100.
+
+            if polarization_angle == 180.:
+                polarization_angle = 0.
+
+            unpolarized_weights = np.full(self.axes['Pol'].nbins, (1. - polarization_level) / self.axes['Pol'].nbins)
+            polarized_weights = np.zeros(self.axes['Pol'].nbins)
+
+            polarization_bin_index = self.axes['Pol'].find_bin(polarization_angle * u.deg)
+            polarized_weights[polarization_bin_index] = polarization_level
+
+            weights = unpolarized_weights + polarized_weights
+
+            contents = np.tensordot(weights, self.contents, axes=([0], [self.axes.label_to_index('Pol')]))
+
+            axes = self.axes['Em', 'Phi', 'PsiChi']
+
         energy_axis = self.photon_energy_axis
 
         flux = get_integrated_spectral_model(spectrum, energy_axis)
 
-        expectation = np.tensordot(self.contents, flux.contents, axes = ([0], [0]))
+        expectation = np.tensordot(contents, flux.contents, axes=([0], [0]))
 
         # Note: np.tensordot loses unit if we use a sparse matrix as it input.
         if self.is_sparse:
             expectation *= self.unit * flux.unit
 
-        hist = Histogram(self.axes[1:], contents = expectation)
-
+        hist = Histogram(axes, contents=expectation)
+
+        if not hist.unit == u.dimensionless_unscaled:
+            raise RuntimeError("Expectation should be dimensionless, but has units of " + str(hist.unit) + ".")
+
         return hist

cosipy/test_data/polarization_ori.ori

@@ -1,4 +1,4 @@
 Type OrientationsGalactic 
-OG 0.0 73.14907746670937 41.85821768724895 16.85092253329064 221.85821768724895 0.0 0.0 0.0 15.0 
-OG 15.0 73.14907746670937 41.85821768724895 16.85092253329064 221.85821768724895 0.0 0.0 0.0 0.0
+OG 0.0 73.14907746670937 41.85821768724895 16.85092253329064 221.85821768724895 523.6546162516867 29.730094567937144 229.94497777814732 15.0 
+OG 15.0 73.14907746670937 41.85821768724895 16.85092253329064 221.85821768724895 523.6546162516867 29.730094567937144 229.94497777814732 0.0
 EN