Specify ntransits

hawc_hal/HAL.py

@@ -50,16 +50,25 @@ class HAL(PluginPrototype):
     :param flat_sky_pixels_size: size of the pixel for the flat sky projection (Hammer Aitoff)
     """
 
-    def __init__(self, name, maptree, response_file, roi, flat_sky_pixels_size=0.17):
+    def __init__(self, name, maptree, response_file, roi, flat_sky_pixels_size=0.17, set_transits=-1):
 
         # Store ROI
         self._roi = roi
 
+        # optionally specify n_transits
+        if (set_transits > 0.0):
+            log.info("Setting Transits to {0}".format(set_transits))
+            n_transits=set_transits
+        else:
+            log.info("Using transits contained in the maptree")
+            n_transits=None
+
+
         # Set up the flat-sky projection
         self._flat_sky_projection = roi.get_flat_sky_projection(flat_sky_pixels_size)
 
         # Read map tree (data)
-        self._maptree = map_tree_factory(maptree, roi=roi)
+        self._maptree = map_tree_factory(maptree, roi, n_transits)
 
         # Read detector response_file
         self._response = hawc_response_factory(response_file)

hawc_hal/maptree/from_hdf5_file.py

@@ -12,13 +12,15 @@
 from ..healpix_handling import SparseHealpix, DenseHealpix
 from .data_analysis_bin import DataAnalysisBin
 
+import numpy as np
 
-def from_hdf5_file(map_tree_file, roi):
+def from_hdf5_file(map_tree_file, roi, n_transits):
     """
     Create a MapTree object from a HDF5 file and a ROI. Do not use this directly, use map_tree_factory instead.
 
     :param map_tree_file:
     :param roi:
+    :param n_transits:
     :return:
     """
 
@@ -62,10 +64,17 @@ def from_hdf5_file(map_tree_file, roi):
 
     bin_names = analysis_bins_df.index.levels[0]
 
-    # Loop over them and build a DataAnalysisBin instance for each one
-
     data_analysis_bins = collections.OrderedDict()
 
+    # Figure out the transits
+    # pick out the transits same as root file
+    if n_transits is None:
+        use_transits = meta_df["n_transits"].max()
+    else:
+        # specified n_transits
+        use_transits = n_transits
+
+
     for bin_name in bin_names:
 
         this_df = analysis_bins_df.loc[bin_name]
@@ -90,15 +99,15 @@ def from_hdf5_file(map_tree_file, roi):
 
             # This signals the DataAnalysisBin that we are dealing with a full sky map
             active_pixels_user = None
-
+            
         # Let's now build the instance
         this_bin = DataAnalysisBin(bin_name,
                                    observation_hpx_map=observation_hpx_map,
                                    background_hpx_map=background_hpx_map,
                                    active_pixels_ids=active_pixels_user,
-                                   n_transits=this_meta['n_transits'],
+                                   n_transits=use_transits,
                                    scheme='RING' if this_meta['scheme'] == 0 else 'NEST')
 
         data_analysis_bins[bin_name] = this_bin
 
-    return data_analysis_bins
+    return data_analysis_bins, use_transits

hawc_hal/maptree/from_root_file.py

@@ -44,12 +44,13 @@ def _get_bin_object(f, bin_name, suffix):
     return bin_tobject
 
 
-def from_root_file(map_tree_file, roi):
+def from_root_file(map_tree_file, roi, n_transits):
     """
     Create a MapTree object from a ROOT file and a ROI. Do not use this directly, use map_tree_factory instead.
 
     :param map_tree_file:
     :param roi:
+    :param n_transits:
     :return:
     """
 
@@ -94,16 +95,20 @@ def from_root_file(map_tree_file, roi):
             data_bins_labels = [ str(i) for i in data_bins_labels ]
 
 
-        # A transit is defined as 1 day, and totalDuration is in hours
-        # Get the number of transit from bin 0 (as LiFF does)
-
-        n_transits = root_numpy.tree2array(f.Get("BinInfo"), "totalDuration") / 24.0
 
         # The map-maker underestimate the livetime of bins with low statistic by removing time intervals with
         # zero events. Therefore, the best estimate of the livetime is the maximum of n_transits, which normally
         # happen in the bins with high statistic
-        n_transits = max(n_transits)
-
+        # Alternatively, specify n_transits
+
+        if n_transits is None:
+            # A transit is defined as 1 day, and totalDuration is in hours
+            # Get the number of transits from bin 0 (as LiFF does)
+            map_transits = root_numpy.tree2array(f.Get("BinInfo"), "totalDuration") / 24.0
+            use_transits = np.max(map_transits)
+        else:
+            use_transits = n_transits
+
         n_bins = len(data_bins_labels)
 
         # These are going to be Healpix maps, one for each data analysis bin_name
@@ -147,7 +152,7 @@ def from_root_file(map_tree_file, roi):
                                                          counts_hpx,
                                                          bkg_hpx,
                                                          active_pixels_ids=active_pixels,
-                                                         n_transits=n_transits,
+                                                         n_transits=use_transits,
                                                          scheme='RING')
 
             else:
@@ -161,12 +166,12 @@ def from_root_file(map_tree_file, roi):
                                                          DenseHealpix(counts),
                                                          DenseHealpix(bkg),
                                                          active_pixels_ids=None,
-                                                         n_transits=n_transits,
+                                                         n_transits=use_transits,
                                                          scheme='RING')
 
             data_analysis_bins[name] = this_data_analysis_bin
 
-    return data_analysis_bins
+    return data_analysis_bins, use_transits
 
 
 def _read_partial_tree(ttree_instance, elements_to_read):

hawc_hal/maptree/map_tree.py

@@ -19,36 +19,55 @@
 import astropy.units as u
 
 
-def map_tree_factory(map_tree_file, roi):
+def map_tree_factory(map_tree_file, roi, n_transits=None):
 
     # Sanitize files in input (expand variables and so on)
     map_tree_file = sanitize_filename(map_tree_file)
 
     if os.path.splitext(map_tree_file)[-1] == '.root':
 
-        return MapTree.from_root_file(map_tree_file, roi)
+        return MapTree.from_root_file(map_tree_file, roi, n_transits)
 
     else:
 
-        return MapTree.from_hdf5(map_tree_file, roi)
+        return MapTree.from_hdf5(map_tree_file, roi, n_transits)
 
 
 class MapTree(object):
 
-    def __init__(self, analysis_bins, roi):
+    def __init__(self, analysis_bins, roi, n_transits=None):
 
         self._analysis_bins = analysis_bins
         self._roi = roi
 
+        if n_transits != None:
+            self._n_transits = n_transits
+        else:
+            # we need to find the transits in the map
+            transits=[]
+            # loop through bins and get the n_transits in each 
+            for analysis_bin_key in self._analysis_bins:
+                transits.append(self._analysis_bins[analysis_bin_key].n_transits)
+
+            #Get the largest
+            self._n_transits = np.max(np.array(transits))
+
     @classmethod
-    def from_hdf5(cls, map_tree_file, roi):
+    def from_hdf5(cls, map_tree_file, roi, n_transits):
+        """
+        Create a MapTree object from a HDF5 file and a ROI. Do not use this directly, use map_tree_factory instead.
 
-        data_analysis_bins = from_hdf5_file(map_tree_file, roi)
+        :param map_tree_file:
+        :param roi:
+        :return:
+        """
 
-        return cls(data_analysis_bins, roi)
+        data_analysis_bins, transits = from_hdf5_file(map_tree_file, roi, n_transits)
+
+        return cls(data_analysis_bins, roi, transits)
 
     @classmethod
-    def from_root_file(cls, map_tree_file, roi):
+    def from_root_file(cls, map_tree_file, roi, n_transits):
         """
         Create a MapTree object from a ROOT file and a ROI. Do not use this directly, use map_tree_factory instead.
 
@@ -57,9 +76,9 @@ def from_root_file(cls, map_tree_file, roi):
         :return:
         """
 
-        data_analysis_bins = from_root_file(map_tree_file, roi)
+        data_analysis_bins, transits = from_root_file(map_tree_file, roi, n_transits)
 
-        return cls(data_analysis_bins, roi)
+        return cls(data_analysis_bins, roi, transits)
 
     def __iter__(self):
         """
@@ -98,6 +117,10 @@ def __len__(self):
 
         return len(self._analysis_bins)
 
+    @property
+    def n_transits(self):
+        return self._n_transits
+
     @property
     def analysis_bins_labels(self):
 

tests/conftest.py

@@ -92,6 +92,18 @@ def check_map_trees(m1, m2):
         assert p1.nside == p2.nside
         assert p1.n_transits == p2.n_transits
 
+def check_n_transits(maptree,ntransits):
+    # Should give maptree object and expected ntransits
+
+    #Check maptree is assigned proper transits
+    assert maptree.n_transits == ntransits
+
+    #Check each bin too
+    for bin_key in maptree:
+
+        the_bin = maptree[bin_key]
+        assert the_bin.n_transits == ntransits
+
 
 def check_responses(r1, r2):
 

tests/test_n_transits.py

@@ -0,0 +1,34 @@
+from hawc_hal.maptree.map_tree import map_tree_factory
+from hawc_hal import HealpixConeROI
+import os
+import pytest
+from conftest import check_n_transits
+
+
+def test_transits(maptree,
+                     roi):
+
+    roi_ = roi
+
+    #Case 1: not specifying transits
+
+    # Dont specify n_transits, use value in maptree
+    maptree_none = map_tree_factory(maptree, roi_)
+
+    n_transits = maptree_none.n_transits
+
+    #This is somewhat trivial, but something has definitely gone wrong if this fails
+    check_n_transits(maptree_none,n_transits)
+
+
+    #Case 2: specify transits
+
+    # specify transits, in this case 777.7
+    n_transits_specify = 777.7
+
+    # specify n_transits
+    maptree_specify = map_tree_factory(maptree, roi_, n_transits_specify)
+
+    # Does maptree return specified transits?
+    check_n_transits(maptree_specify,n_transits_specify)
+