Extended Source Response Generator

cosipy/response/ExtendedSourceResponse.py

@@ -28,6 +28,9 @@ def __init__(self, *args, **kwargs):
         """
         Initialize an ExtendedSourceResponse object.
         """
+        # Not to track the under/overflow bins
+        kwargs['track_overflow'] = False
+
         super().__init__(*args, **kwargs)
 
         if not np.all(self.axes.labels == ['NuLambda', 'Ei', 'Em', 'Phi', 'PsiChi']):

cosipy/response/FullDetectorResponse.py

@@ -1,5 +1,6 @@
 from .PointSourceResponse import PointSourceResponse
 from .DetectorResponse import DetectorResponse
+from .ExtendedSourceResponse import ExtendedSourceResponse
 from astromodels.core.model_parser import ModelParser
 import matplotlib.pyplot as plt
 from astropy.time import Time
@@ -11,6 +12,7 @@
 import numpy as np
 import mhealpy as hp
 from mhealpy import HealpixBase, HealpixMap
+import glob
 
 from scipy.special import erf
 
@@ -28,7 +30,8 @@
 
 from copy import copy, deepcopy
 import gzip
-from tqdm import tqdm
+#from tqdm import tqdm
+from tqdm.autonotebook import tqdm
 import subprocess
 import sys
 import pathlib
@@ -913,6 +916,175 @@
                 return psr[0]
             else:
                 return psr
+
+    def _setup_extended_source_response_params(self, coordsys, nside_image, nside_scatt_map):
+        """
+        Validate coordinate system and setup NSIDE parameters for extended source response generation.
+
+        Parameters
+        ----------
+        coordsys : str
+            Coordinate system to be used (currently only 'galactic' is supported)
+        nside_image : int or None
+            NSIDE parameter for the image reconstruction.
+            If None, uses the full detector response's NSIDE.
+        nside_scatt_map : int or None
+            NSIDE parameter for scatt map generation.
+            If None, uses the full detector response's NSIDE.
+
+        Returns
+        -------
+        tuple
+            (coordsys, nside_image, nside_scatt_map) : validated parameters
+        """
+        if coordsys != 'galactic':
+            raise ValueError(f'The coordsys {coordsys} not currently supported')
+
+        if nside_image is None:
+            nside_image = self.nside
+
+        if nside_scatt_map is None:
+            nside_scatt_map = self.nside
+
+        return coordsys, nside_image, nside_scatt_map
+
+    def get_point_source_response_per_image_pixel(self, ipix_image, orientation, coordsys = 'galactic', nside_image = None, nside_scatt_map = None, Earth_occ = True):
+        """
+        Generate point source response for a specific HEALPix pixel by convolving 
+        the all-sky detector response with exposure.
+
+        Parameters
+        ----------
+        ipix_image : int
+            HEALPix pixel index
+        orientation : cosipy.spacecraftfile.SpacecraftFile
+            Spacecraft attitude information
+        coordsys : str, default 'galactic'
+            Coordinate system (currently only 'galactic' is supported)
+        nside_image : int, optional
+            NSIDE parameter for image reconstruction.
+            If None, uses the detector response's NSIDE.
+        nside_scatt_map : int, optional
+            NSIDE parameter for scatt map generation.
+            If None, uses the detector response's NSIDE.
+        Earth_occ : bool, default True
+            Whether to include Earth occultation in the response
+
+        Returns
+        -------
+        :py:class:`PointSourceResponse`
+            Point source response for the specified pixel
+        """
+        coordsys, nside_image, nside_scatt_map = self._setup_extended_source_response_params(coordsys, nside_image, nside_scatt_map)
+
+        image_axes = HealpixAxis(nside = nside_image, coordsys = coordsys, scheme='ring', label = 'NuLambda') # The label should be 'lb' in the future
+
+        coord = image_axes.pix2skycoord(ipix_image)
+
+        scatt_map = orientation.get_scatt_map(target_coord = coord,
+                                              nside = nside_scatt_map,
+                                              scheme='ring',
+                                              coordsys=coordsys,
+                                              earth_occ=Earth_occ)
+
+        psr = self.get_point_source_response(coord = coord, scatt_map = scatt_map, Earth_occ = Earth_occ)
+
+        return psr
+
+    def get_extended_source_response(self, orientation, coordsys = 'galactic', nside_image = None, nside_scatt_map = None, Earth_occ = True):
+        """
+        Generate extended source response by convolving the all-sky detector 
+        response with exposure over the entire sky.
+
+        Parameters
+        ----------
+        orientation : cosipy.spacecraftfile.SpacecraftFile
+            Spacecraft attitude information
+        coordsys : str, default 'galactic'
+            Coordinate system (currently only 'galactic' is supported)
+        nside_image : int, optional
+            NSIDE parameter for image reconstruction.
+            If None, uses the detector response's NSIDE.
+        nside_scatt_map : int, optional
+            NSIDE parameter for scatt map generation.
+            If None, uses the detector response's NSIDE.
+        Earth_occ : bool, default True
+            Whether to include Earth occultation in the response
+
+        Returns
+        -------
+        :py:class:`ExtendedSourceResponse`
+            Extended source response covering the entire sky
+        """
+        coordsys, nside_image, nside_scatt_map = self._setup_extended_source_response_params(coordsys, nside_image, nside_scatt_map)
+
+        axes = [HealpixAxis(nside = nside_image, coordsys = coordsys, scheme='ring', label = 'NuLambda')] # The label should be 'lb' in the future
+        axes += list(self.axes[1:])
+        axes[-1].coordsys = coordsys
+
+        extended_source_response = ExtendedSourceResponse(axes, unit = u.Unit("cm2 s"))
+
+        for ipix in tqdm(range(hp.nside2npix(nside_image))):
+
+            psr = self.get_point_source_response_per_image_pixel(ipix, orientation, coordsys = coordsys, 
+                                                                 nside_image = nside_image, nside_scatt_map = nside_scatt_map, Earth_occ = Earth_occ)
+
+            extended_source_response[ipix] = psr.contents
+
+        return extended_source_response
+
+    def merge_psr_to_extended_source_response(self, basename, coordsys = 'galactic', nside_image = None):
+        """
+        Create extended source response by merging multiple point source responses.
+
+        Reads point source response files matching the pattern `basename` + index + file_extension. 
+        For example, with basename='histograms/hist_', filenames are expected to be like 'histograms/hist_00001.hdf5'.
+
+        Parameters
+        ----------
+        basename : str
+            Base filename pattern for point source response files
+        coordsys : str, default 'galactic'
+            Coordinate system (currently only 'galactic' is supported)
+        nside_image : int, optional
+            NSIDE parameter for image reconstruction.
+            If None, uses the detector response's NSIDE.
+
+        Returns
+        -------
+        :py:class:`ExtendedSourceResponse`
+            Combined extended source response
+        """
+        coordsys, nside_image, _ = self._setup_extended_source_response_params(coordsys, nside_image, None)
+
+        psr_files = glob.glob(basename + "*")
+
+        if not psr_files:
+            raise FileNotFoundError(f"No files found matching pattern {basename}*")
+
+        axes = [HealpixAxis(nside = nside_image, coordsys = coordsys, scheme='ring', label = 'NuLambda')] # The label should be 'lb' in the future
+        axes += list(self.axes[1:])
+        axes[-1].coordsys = coordsys
+
+        extended_source_response = ExtendedSourceResponse(axes, unit = u.Unit("cm2 s"))
+
+        filled_pixels = []
+
+        for filename in psr_files:
+
+            ipix = int(filename[len(basename):].split(".")[0])
+
+            psr = Histogram.open(filename)
+
+            extended_source_response[ipix] = psr.contents
+
+            filled_pixels.append(ipix)
+
+        expected_pixels = set(range(extended_source_response.axes[0].npix))
+        if set(filled_pixels) != expected_pixels:
+            raise ValueError(f"Missing pixels in the response files. Expected {extended_source_response.axes[0].npix} pixels, got {len(filled_pixels)} pixels")
+
+        return extended_source_response
 
     @staticmethod
     def _sum_rot_hist(h, h_new, exposure, axis = "PsiChi"):

docs/tutorials/response/extended_source_response_generator.py

@@ -0,0 +1,28 @@
+#!/usr/bin/env python
+# coding: UTF-8
+
+import sys
+import logging
+logger = logging.getLogger('cosipy')
+logger.setLevel(logging.INFO)
+logger.addHandler(logging.StreamHandler(sys.stdout))
+
+from cosipy.spacecraftfile import SpacecraftFile
+from cosipy.response import FullDetectorResponse, ExtendedSourceResponse
+
+# file path
+full_detector_response_path = "SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.h5"
+orientation_path = "20280301_3_month_with_orbital_info.ori"
+
+# load response and orientation
+full_detector_response = FullDetectorResponse.open(full_detector_response_path)
+orientation = SpacecraftFile.parse_from_file(orientation_path)
+
+# generate your extended source response
+extended_source_response = full_detector_response.get_extended_source_response(orientation,
+                                                                               coordsys='galactic',
+                                                                               nside_scatt_map=None,
+                                                                               Earth_occ=True)
+
+# save the extended source response
+extended_source_response.write("extended_source_response_continuum.h5", overwrite = True)

docs/tutorials/response/extended_source_response_generator_with_multiple_nodes.py

@@ -0,0 +1,39 @@
+#!/usr/bin/env python
+# coding: UTF-8
+
+import sys
+import logging
+logger = logging.getLogger('cosipy')
+logger.setLevel(logging.INFO)
+logger.addHandler(logging.StreamHandler(sys.stdout))
+
+from cosipy.spacecraftfile import SpacecraftFile
+from cosipy.response import FullDetectorResponse, ExtendedSourceResponse
+
+# file path
+full_detector_response_path = "SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.h5"
+orientation_path = "20280301_3_month_with_orbital_info.ori"
+
+# load response and orientation
+full_detector_response = FullDetectorResponse.open(full_detector_response_path)
+orientation = SpacecraftFile.parse_from_file(orientation_path)
+
+# set the healpix pixel index list
+ipix_image_list = [int(_) for _ in sys.argv[1:]]
+
+print(ipix_image_list)
+
+# generate a point source response at each pixel
+basename = "psr/psr_"
+
+for ipix_image in ipix_image_list:
+
+    psr = full_detector_response.get_point_source_response_per_image_pixel(ipix_image, orientation, 
+                                                                           coordsys='galactic',
+                                                                           nside_image=None,
+                                                                           nside_scatt_map=None,
+                                                                           Earth_occ=True)
+
+    psr.write(f"{basename}{ipix_image:08}.h5",overwrite = True)
+
+# see also merge_response_generated_with_mutiple_nodes.py to know how we can merge the above point source responses as a single extended source response.

docs/tutorials/response/merge_response_generated_with_mutiple_nodes.py

@@ -0,0 +1,24 @@
+#!/usr/bin/env python
+# coding: UTF-8
+
+import sys
+import logging
+logger = logging.getLogger('cosipy')
+logger.setLevel(logging.INFO)
+logger.addHandler(logging.StreamHandler(sys.stdout))
+
+from cosipy.spacecraftfile import SpacecraftFile
+from cosipy.response import FullDetectorResponse, ExtendedSourceResponse
+
+# load full detector response 
+full_detector_response_path = "SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.h5"
+full_detector_response = FullDetectorResponse.open(full_detector_response_path)
+
+# basename should be the same as one used before
+basename = "psr/psr_"
+
+# merge the point source responses
+extended_source_response = full_detector_response.merge_psr_to_extended_source_response(basename, coordsys = 'galactic', nside_image = None)
+
+# save the extended source response
+extended_source_response.write("extended_source_response_continuum_merged.h5", overwrite = True)