diff --git a/cosipy/response/ExtendedSourceResponse.py b/cosipy/response/ExtendedSourceResponse.py index 044760dd..433114ab 100644 --- a/cosipy/response/ExtendedSourceResponse.py +++ b/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']): diff --git a/cosipy/response/FullDetectorResponse.py b/cosipy/response/FullDetectorResponse.py index 6673e243..14a775e7 100644 --- a/cosipy/response/FullDetectorResponse.py +++ b/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 @@ def get_point_source_response(self, 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"): diff --git a/cosipy/test_data/test_full_detector_response_dense.h5 b/cosipy/test_data/test_full_detector_response_dense.h5 new file mode 100644 index 00000000..7e45a33c Binary files /dev/null and b/cosipy/test_data/test_full_detector_response_dense.h5 differ diff --git a/docs/tutorials/response/extended_source_response_generator.py b/docs/tutorials/response/extended_source_response_generator.py new file mode 100644 index 00000000..57143d91 --- /dev/null +++ b/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) diff --git a/docs/tutorials/response/extended_source_response_generator_with_multiple_nodes.py b/docs/tutorials/response/extended_source_response_generator_with_multiple_nodes.py new file mode 100644 index 00000000..5843079a --- /dev/null +++ b/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. diff --git a/docs/tutorials/response/merge_response_generated_with_mutiple_nodes.py b/docs/tutorials/response/merge_response_generated_with_mutiple_nodes.py new file mode 100644 index 00000000..8e279e5c --- /dev/null +++ b/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) diff --git a/tests/response/test_full_detector_response.py b/tests/response/test_full_detector_response.py index c4a663fc..a5136194 100644 --- a/tests/response/test_full_detector_response.py +++ b/tests/response/test_full_detector_response.py @@ -7,8 +7,10 @@ from cosipy import test_data from cosipy.response import FullDetectorResponse +from cosipy.spacecraftfile import SpacecraftFile -response_path = test_data.path / "test_full_detector_response.h5" +response_path = test_data.path / "test_full_detector_response_dense.h5" +orientation_path = test_data.path / "20280301_first_10sec.ori" def test_open(): @@ -16,10 +18,10 @@ def test_open(): assert response.filename == response_path - assert response.ndim == 7 + assert response.ndim == 5 assert arr_eq(response.axes.labels, - ['NuLambda', 'Ei', 'Em', 'Phi', 'PsiChi', 'SigmaTau', 'Dist']) + ['NuLambda', 'Ei', 'Em', 'Phi', 'PsiChi']) assert response.unit.is_equivalent('m2') @@ -30,10 +32,10 @@ def test_get_item(): drm = response[0] - assert drm.ndim == 6 + assert drm.ndim == 4 assert arr_eq(drm.axes.labels, - ['Ei', 'Em', 'Phi', 'PsiChi', 'SigmaTau', 'Dist']) + ['Ei', 'Em', 'Phi', 'PsiChi']) assert drm.unit.is_equivalent('m2') @@ -45,16 +47,56 @@ def test_get_interp_response(): lat = 0*u.deg, frame = SpacecraftFrame())) - assert drm.ndim == 6 + assert drm.ndim == 4 assert arr_eq(drm.axes.labels, - ['Ei', 'Em', 'Phi', 'PsiChi', 'SigmaTau', 'Dist']) + ['Ei', 'Em', 'Phi', 'PsiChi']) assert drm.unit.is_equivalent('m2') - +def test_get_extended_source_response(): + orientation = SpacecraftFile.parse_from_file(orientation_path) - - - + with FullDetectorResponse.open(response_path) as response: + + extended_source_response = response.get_extended_source_response(orientation, + coordsys = 'galactic', + nside_image = None, + nside_scatt_map = None, + Earth_occ = True) + + assert extended_source_response.ndim == 5 + + assert arr_eq(extended_source_response.axes.labels, + ['NuLambda', 'Ei', 'Em', 'Phi', 'PsiChi']) + + assert extended_source_response.unit.is_equivalent('cm2 s') + +def test_merge_psr_to_extended_source_response(tmp_path): + + orientation = SpacecraftFile.parse_from_file(orientation_path) + + with FullDetectorResponse.open(response_path) as response: + + for ipix_image in range(response.axes['NuLambda'].npix): + + psr = 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(tmp_path / f"psr_{ipix_image:08}.h5") + + + extended_source_response = response.merge_psr_to_extended_source_response(str(tmp_path / "psr_"), + coordsys = 'galactic', + nside_image = None) + + assert extended_source_response.ndim == 5 + + assert arr_eq(extended_source_response.axes.labels, + ['NuLambda', 'Ei', 'Em', 'Phi', 'PsiChi']) + + assert extended_source_response.unit.is_equivalent('cm2 s')