diff --git a/hawc_hal/HAL.py b/hawc_hal/HAL.py index 0e3b51d..5a56609 100644 --- a/hawc_hal/HAL.py +++ b/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) diff --git a/hawc_hal/maptree/from_hdf5_file.py b/hawc_hal/maptree/from_hdf5_file.py index 0e23979..925cd7b 100644 --- a/hawc_hal/maptree/from_hdf5_file.py +++ b/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 diff --git a/hawc_hal/maptree/from_root_file.py b/hawc_hal/maptree/from_root_file.py index 00c705e..ab8f0cd 100644 --- a/hawc_hal/maptree/from_root_file.py +++ b/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): diff --git a/hawc_hal/maptree/map_tree.py b/hawc_hal/maptree/map_tree.py index 9e3e888..f71b35a 100644 --- a/hawc_hal/maptree/map_tree.py +++ b/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): diff --git a/tests/conftest.py b/tests/conftest.py index f136176..e3b96d3 100644 --- a/tests/conftest.py +++ b/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): diff --git a/tests/test_n_transits.py b/tests/test_n_transits.py new file mode 100644 index 0000000..be9a130 --- /dev/null +++ b/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) +