Update: Include declination bin and energy bin selection in IRF checks

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)

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

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

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

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

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:
 

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)
 

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)

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)
 

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)

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)