Specify ntransits

hawc_hal/HAL.py

@@ -47,16 +47,19 @@ 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, n_transits=None):
 
         # Store ROI
         self._roi = roi
 
+        # optionally specify n_transits
+        assert (n_transits==None or n_transits > 0.0), "You must specify n_transits >0"
+
         # 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

@@ -10,12 +10,13 @@
 from .data_analysis_bin import DataAnalysisBin
 
 
-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:
     """
 
@@ -88,12 +89,17 @@ 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
 
+        # Specify n_transits (or not), default value contained in map is this_meta['n_transits']
+        use_transits=this_meta['n_transits']
+        if n_transits!=None:
+            use_transits=n_transits
+
         # 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

hawc_hal/maptree/from_root_file.py

@@ -41,12 +41,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:
     """
 
@@ -93,15 +94,17 @@ def from_root_file(map_tree_file, roi):
 
 
         # 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
+        # Get the number of transits from bin 0 (as LiFF does)
+        map_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
+        use_transits = max(map_transits)
+        if n_transits != None:
+            use_transits = n_transits
+
         n_bins = len(data_bins_labels)
 
         # These are going to be Healpix maps, one for each data analysis bin_name
@@ -145,7 +148,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:
@@ -159,7 +162,7 @@ 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

hawc_hal/maptree/map_tree.py

@@ -17,18 +17,18 @@
 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):
@@ -39,14 +39,14 @@ def __init__(self, analysis_bins, roi):
         self._roi = roi
 
     @classmethod
-    def from_hdf5(cls, map_tree_file, roi):
+    def from_hdf5(cls, map_tree_file, roi, n_transits):
 
-        data_analysis_bins = from_hdf5_file(map_tree_file, roi)
+        data_analysis_bins = from_hdf5_file(map_tree_file, roi, n_transits)
 
         return cls(data_analysis_bins, roi)
 
     @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.
 
@@ -55,7 +55,7 @@ def from_root_file(cls, map_tree_file, roi):
         :return:
         """
 
-        data_analysis_bins = from_root_file(map_tree_file, roi)
+        data_analysis_bins = from_root_file(map_tree_file, roi, n_transits)
 
         return cls(data_analysis_bins, roi)