diff --git a/hawc_hal/HAL.py b/hawc_hal/HAL.py index 22ed34f..ddbde26 100644 --- a/hawc_hal/HAL.py +++ b/hawc_hal/HAL.py @@ -12,6 +12,7 @@ import matplotlib.pyplot as plt import numpy as np import pandas as pd +from hawc_hal.dec_band_select import dec_index_search from astromodels import Parameter from astropy.convolution import Gaussian2DKernel from astropy.convolution import convolve_fft as convolve @@ -36,6 +37,10 @@ SparseHealpix, get_gnomonic_projection, ) +from hawc_hal.region_of_interest import ( + HealpixConeROI, + HealpixMapROI +) from hawc_hal.log_likelihood import log_likelihood from hawc_hal.maptree import map_tree_factory from hawc_hal.maptree.data_analysis_bin import DataAnalysisBin @@ -66,6 +71,7 @@ def __init__( roi, flat_sky_pixels_size=0.17, set_transits=None, + bin_list=None ): # Store ROI @@ -80,6 +86,11 @@ def __init__( n_transits = None log.info("Using transits contained in maptree") + if roi: + dec_arange=[roi.ra_dec_center[1]-roi.model_radius.to(u.deg).value, roi.ra_dec_center[1]+roi.model_radius.to(u.deg).value] + dec_list = dec_index_search(response_file, dec_arange, use_module=True) + print("Using declination list=", dec_list) + # Set up the flat-sky projection self.flat_sky_pixels_size = flat_sky_pixels_size self._flat_sky_projection = self._roi.get_flat_sky_projection(self.flat_sky_pixels_size) @@ -88,7 +99,7 @@ def __init__( self._maptree = map_tree_factory(maptree, roi=self._roi, n_transits=n_transits) # Read detector response_file - self._response = hawc_response_factory(response_file) + self._response = hawc_response_factory(response_file, bin_list, dec_list) # Use a renormalization of the background as nuisance parameter # NOTE: it is fixed to 1.0 unless the user explicitly sets it free (experimental) diff --git a/hawc_hal/convolved_source/convolved_extended_source.py b/hawc_hal/convolved_source/convolved_extended_source.py index f76cf58..089fa95 100644 --- a/hawc_hal/convolved_source/convolved_extended_source.py +++ b/hawc_hal/convolved_source/convolved_extended_source.py @@ -50,26 +50,31 @@ def __init__(self, source, response, flat_sky_projection): # Figure out maximum and minimum declination to be covered dec_min = max(min([dec1, dec2, dec3, dec4]), lat_start) dec_max = min(max([dec1, dec2, dec3, dec4]), lat_stop) - + log.info("Lat start = %.3f, Lat stop = %.3f" %(lat_start, lat_stop)) + log.info("Min = %.3f, Max = %.3f" %(dec_min, dec_max)) + # Get the defined dec bins lower edges # Get the defined dec bins lower edges lower_edges = np.array([x[0] for x in response.dec_bins]) upper_edges = np.array([x[-1] for x in response.dec_bins]) centers = np.array([x[1] for x in response.dec_bins]) dec_bins_to_consider_idx = np.flatnonzero((upper_edges >= dec_min) & (lower_edges <= dec_max)) - + log.info("Central bins = %s" %(dec_bins_to_consider_idx)) # Wrap the selection so we have always one bin before and one after. # NOTE: we assume that the ROI do not overlap with the very first or the very last dec bin # Add one dec bin to cover the last part dec_bins_to_consider_idx = np.append(dec_bins_to_consider_idx, [dec_bins_to_consider_idx[-1] + 1]) # Add one dec bin to cover the first part - dec_bins_to_consider_idx = np.insert(dec_bins_to_consider_idx, 0, [dec_bins_to_consider_idx[0] - 1]) - + if (dec_bins_to_consider_idx[0]!=0): + dec_bins_to_consider_idx = np.insert(dec_bins_to_consider_idx, 0, [dec_bins_to_consider_idx[0] - 1]) + #else (Rishi change code to remove dec band below 0) + log.info("The bins are %s" %(dec_bins_to_consider_idx)) self._dec_bins_to_consider = [response.response_bins[centers[x]] for x in dec_bins_to_consider_idx] log.info("Considering %i dec bins for extended source %s" % (len(self._dec_bins_to_consider), self._name)) - + log.info("The bins are %s" %(dec_bins_to_consider_idx)) + self.return_bins = dec_bins_to_consider_idx # Find central bin for the PSF dec_center = (lat_start + lat_stop) / 2.0 diff --git a/hawc_hal/dec_band_select.py b/hawc_hal/dec_band_select.py new file mode 100644 index 0000000..103b8b9 --- /dev/null +++ b/hawc_hal/dec_band_select.py @@ -0,0 +1,57 @@ +from __future__ import division +from builtins import zip +from past.utils import old_div +from builtins import object +import argparse as ap +import numpy as np +import scipy.stats as stats +import warnings +import os +import sys +import corner as cn +import healpy as hp +import scipy + +from pathlib import Path + +import uproot +from threeML.io.logging import setup_logger +from astropy.coordinates import SkyCoord + +with warnings.catch_warnings(): + import astromodels + +log = setup_logger(__name__) +log.propagate = False + +def dec_index_search(response_file, dec_var, use_module): + response = response_file + with uproot.open(response) as response_file: + dec_bins_lower_edge = response_file["DecBins/lowerEdge"].array().to_numpy() + dec_bins_upper_edge = response_file["DecBins/upperEdge"].array().to_numpy() + dec_bins_center = response_file["DecBins/simdec"].array().to_numpy() + if use_module==True: + dec_min = np.amin(dec_var) + dec_max = np.amax(dec_var) + else: + dec_min = -37.5 + dec_max = 77.5 + + log.info("Dec min= %.3f" %(dec_min)) + log.info("Dec max= %.3f" %(dec_max)) + #lower_edges = np.array([-37.5, -32.5, -27.5, -22.5, -17.5, -12.5, -7.5, -2.5, 2.5, 7.5, 12.5, 17.5, 22.5, 27.5, 32.5, 37.5, 42.5, 47.5, 52.5, 57.5, 62.5, 67.5, 72.5], dtype=np.float) + #upper_edges = np.array([-32.5, -27.5, -22.5, -17.5, -12.5, -7.5, -2.5, 2.5, 7.5, 12.5, 17.5, 22.5, 27.5, 32.5, 37.5, 42.5, 47.5, 52.5, 57.5, 62.5, 67.5, 72.5, 77.5], dtype=np.float) + #centers = np.array([-35., -30., -25., -20., -15., -10., -5., 0., 5., 10., 15., 20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75.], dtype=np.float) + + dec_bins_to_consider_idx = np.flatnonzero((dec_bins_upper_edge >= dec_min) & (dec_bins_lower_edge <= dec_max)) + #dec_bins_to_consider_idx = np.flatnonzero((upper_edges >= dec_min) & (lower_edges <= dec_max)) + #log.info("decbins before adding the extra bins : %s" %(dec_bins_to_consider_idx)) + #log.info("Decbins to from response file before adding the extra bins : %s" %(dec_bins_to_consider_idx2)) + + dec_bins_to_consider_idx = np.append(dec_bins_to_consider_idx, [dec_bins_to_consider_idx[-1] + 1]) + # Add one dec bin to cover the first part + dec_bins_to_consider_idx = np.insert(dec_bins_to_consider_idx, 0, [dec_bins_to_consider_idx[0] - 1]) + # Rescale bins to be remove dec bins less than 0 and greater than 22 + dec_bins_in_use=dec_bins_to_consider_idx[(dec_bins_to_consider_idx >=0) & (23>dec_bins_to_consider_idx)] + log.info("Declination bins to read response file: %s" %(dec_bins_in_use)) + return dec_bins_in_use \ No newline at end of file diff --git a/hawc_hal/response/response.py b/hawc_hal/response/response.py index 7f37bff..4027722 100644 --- a/hawc_hal/response/response.py +++ b/hawc_hal/response/response.py @@ -25,7 +25,7 @@ _instances = {} -def hawc_response_factory(response_file_name): +def hawc_response_factory(response_file_name, bin_list2=None, dec_list2=None): """ A factory function for the response which keeps a cache, so that the same response is not read over and over again. @@ -49,8 +49,8 @@ def hawc_response_factory(response_file_name): if extension == ".root": - new_instance = HAWCResponse.from_root_file(response_file_name) - + #new_instance = HAWCResponse.from_root_file(response_file_name) + new_instance = HAWCResponse.from_root_file(response_file_name, bin_list2, dec_list2) elif extension in [".hd5", ".hdf5", ".hdf"]: new_instance = HAWCResponse.from_hdf5(response_file_name) @@ -181,7 +181,7 @@ def from_hdf5(cls, response_file_name): return cls(response_file_name, dec_bins, response_bins) @classmethod - def from_root_file(cls, response_file_name: Path): + def from_root_file(cls, response_file_name: Path, bin_list2, dec_list2): """Build response from ROOT file. Do not use directly, use the hawc_response_factory instead. @@ -248,7 +248,8 @@ def from_root_file(cls, response_file_name: Path): # Now we create a dictionary of ResponseBin instances for each dec bin name response_bins = collections.OrderedDict() - for dec_id in range(len(dec_bins)): + #for dec_id in range(len(dec_bins)): + for dec_id in dec_list2: this_response_bins = collections.OrderedDict() min_dec, dec_center, max_dec = dec_bins[dec_id] @@ -261,7 +262,7 @@ def from_root_file(cls, response_file_name: Path): n_energy_bins = len(response_file[dec_id_label].keys(recursive=False)) response_bins_ids = list(range(n_energy_bins)) - + log.info("Dec ID= %s" %(dec_id)) for response_bin_id in response_bin_ids: this_response_bin = ResponseBin.from_ttree( @@ -272,7 +273,8 @@ def from_root_file(cls, response_file_name: Path): log_log_shape, min_dec, dec_center, - max_dec, + max_dec, + bin_list2 ) this_response_bins[response_bin_id] = this_response_bin diff --git a/hawc_hal/response/response_bin.py b/hawc_hal/response/response_bin.py index 06297f6..135a1d1 100644 --- a/hawc_hal/response/response_bin.py +++ b/hawc_hal/response/response_bin.py @@ -78,6 +78,7 @@ def from_ttree( min_dec: np.ndarray, dec_center: np.ndarray, max_dec: np.ndarray, + bin_list3 ): """ Obtain the information from Response ROOT file diff --git a/hawc_hal/tests/conftest.py b/hawc_hal/tests/conftest.py index f3706a0..f650f34 100644 --- a/hawc_hal/tests/conftest.py +++ b/hawc_hal/tests/conftest.py @@ -28,7 +28,7 @@ def fit_point_source(roi, maptree, response, point_source_model, liff=False): if not liff: # This is a 3ML plugin - hawc = HAL("HAWC", maptree, response, roi) + hawc = HAL("HAWC", maptree, response, roi, bin_list=None) else: diff --git a/hawc_hal/tests/test_copy.py b/hawc_hal/tests/test_copy.py index c42cb16..8016054 100644 --- a/hawc_hal/tests/test_copy.py +++ b/hawc_hal/tests/test_copy.py @@ -13,8 +13,8 @@ def deepcopy_hal(theMaptree, theResponse, extended=False): src_name = 'test_source' roi = HealpixConeROI(data_radius=5., model_radius=8., ra=src_ra, dec=src_dec) - - hawc = HAL('HAWC', theMaptree, theResponse, roi) + bins = ['1', '2', '3', '4', '5', '6', '7', '8', '9'] + hawc = HAL('HAWC', theMaptree, theResponse, roi, bin_list=bins) hawc.set_active_measurements(1, 9) data = DataList(hawc) diff --git a/hawc_hal/tests/test_declination_bin_selection.py b/hawc_hal/tests/test_declination_bin_selection.py new file mode 100644 index 0000000..3e36d13 --- /dev/null +++ b/hawc_hal/tests/test_declination_bin_selection.py @@ -0,0 +1,30 @@ +import os +from hawc_hal import HAL, HealpixConeROI +from threeML import * +import astromodels +import astropy.units as u +import pytest +from hawc_hal import HealpixConeROI +from hawc_hal.maptree.map_tree import map_tree_factory +from hawc_hal.region_of_interest import ( + HealpixConeROI, + HealpixMapROI +) +from hawc_hal.dec_band_select import dec_index_search +test_data_path = os.path.join(os.path.abspath(os.path.dirname(__file__)), "data") +response=os.path.join(test_data_path, "geminga_response.root") +ra_source, dec_source = 101.7, 16 +data_r = 3 +roi = HealpixConeROI(data_radius=data_r, + model_radius=data_r + 10.0, + ra=ra_source, dec=dec_source) + +dec_arange=[roi.ra_dec_center[1]-roi.model_radius.to(u.deg).value, roi.ra_dec_center[1]+roi.model_radius.to(u.deg).value] +response_file=response +#Use the declination band selection module +dec_list = dec_index_search(response_file, dec_arange, use_module=True) +print("Using declination list=", dec_list) + +#Using all declination bands +dec_list = dec_index_search(response_file, dec_arange, use_module=False) +print("Using declination list=", dec_list) \ No newline at end of file diff --git a/hawc_hal/tests/test_geminga_paper.py b/hawc_hal/tests/test_geminga_paper.py index ca6d44b..0bafb56 100644 --- a/hawc_hal/tests/test_geminga_paper.py +++ b/hawc_hal/tests/test_geminga_paper.py @@ -158,13 +158,14 @@ def go(args): roi = HealpixConeROI(data_radius=rad, model_radius=rad + 10.0, ra=ra_c, dec=dec_c) - + bins = ['1', '2', '3', '4', '5', '6', '7', '8', '9'] llh = HAL("HAWC", args.mtfile, args.rsfile, - roi) + roi, + bin_list=bins) - llh.set_active_measurements(args.startBin, args.stopBin) + llh.set_active_measurements(bin_list=bins) print(lm) diff --git a/hawc_hal/tests/test_plugin_creation.py b/hawc_hal/tests/test_plugin_creation.py index 1de128d..b26677f 100644 --- a/hawc_hal/tests/test_plugin_creation.py +++ b/hawc_hal/tests/test_plugin_creation.py @@ -13,12 +13,16 @@ import pytest -def test_plugin_from_root(geminga_maptree, geminga_response, geminga_roi): +@pytest.fixture +def bin_list(): + return bin_list; - hal = HAL("HAL", geminga_maptree, geminga_response, geminga_roi) +def test_plugin_from_root(geminga_maptree, geminga_response, geminga_roi, bin_list): + + hal = HAL("HAL", geminga_maptree, geminga_response, geminga_roi, bin_list=None) def test_plugin_from_hd5(maptree, response, roi): roi = HealpixConeROI(ra=82.628, dec=22.640, data_radius=5, model_radius=10) - hal = HAL("HAL", maptree, response, roi) + hal = HAL("HAL", maptree, response, roi, bin_list=None) diff --git a/hawc_hal/tests/test_radial_profiles.py b/hawc_hal/tests/test_radial_profiles.py index 73518da..6cad9e1 100644 --- a/hawc_hal/tests/test_radial_profiles.py +++ b/hawc_hal/tests/test_radial_profiles.py @@ -14,8 +14,8 @@ def test_fit(roi, maptree, response, point_source_model): pts_model = point_source_model - - hawc = HAL("HAWC", maptree, response, roi) + bins = ['1', '2', '3', '4', '5', '6', '7', '8', '9'] + hawc = HAL("HAWC", maptree, response, roi, bin_list=bins) hawc.set_active_measurements(1, 9)